Signatures of unconventional pairing in near-vortex electronic structure of LiFeAs

For the non-self-consistent calculation, we impose a gap function in the form given by equation (3) and find the low lying eigenvalues and eigenstates of $\mathcal {H}^{\mathrm {BdG}}$ using the Lanczos algorithm². The LDOS can be calculated from the eigenenergies Eⁿ and eigenstates [uⁿ(r),vⁿ(r)] as

$$\begin{equation} N({\mathbf{r}}, E) =\! \sum_{n} | u^n ({\mathbf{r}}) |^2 \delta( E - E^n ) + | v^n ({\mathbf{r}}) |^2 \delta( E + E^n ). \end{equation} \tag{ 7 }$$

Since we are not interested in the absolute value of the LDOS but rather in the spatial profile at a given energy, we normalize the LDOS such that for a given energy E, the maximum value of N(r,E) is unity.

For the self-consistent calculation, we assume initial gap functions and use the eigenvalues and eigenvectors of equation (1) to calculate the gap functions given by equation (2). We proceed iteratively until self-consistency is achieved. In diagonalizing $\mathcal {H}^{\mathrm {BdG}}$, we can no longer make use of the crystal momentum basis to simplify the problem since the Peierls phase factor prevents the kinetic part of the Hamiltonian from commuting with the ordinary lattice translation operator T_R. However, the kinetic part commutes with the magnetic translation operator

$$\begin{equation} \hat{T}_{{\mathbf{R}}}\equiv \,\mathrm{e}^{ -\mathrm{i} \frac{\pi}{\Phi_{0}} {\mathbf{A}} ({\mathbf{R}}) \cdot {\mathbf{r}} } T_{{\mathbf{R}}} \end{equation} \tag{ 8 }$$

for a magnetic lattice vector R whose unit cell contains two magnetic fluxoids.

The pairing term in general does not commute with $\hat {T}_{\mathbf {R}}$. Nevertheless, when vortices form a lattice, $\hat {T}_{\mathbf {R}}$ commutes with the pairing term when R is a vector of a vortex sublattice containing every other vortex. Since we focus on the electronic structure near a single vortex, we expect the shape of the vortex lattice to have little influence on our results. Therefore, we make an arbitrary choice for its primitive vectors to be $L_x \hat {{\mathbf {x}}}$ and $L_y \hat {{\mathbf {y}}}$, such that R forms a rectangular lattice R = (m_xL_x,m_yL_y), where m_α = 0···M_α − 1 and M_α ≡ N_α/L_α.³ Note that periodic boundary conditions in the Landau gauge ${\mathbf {A}}({\mathbf {r}}) = -H y \hat {{\mathbf {x}}}$ require the total magnetic flux through the system to be an integer multiple of 2Φ₀N_x. In addition, one magnetic unit cell contains a magnetic flux of 2Φ₀, i.e. H = 2Φ₀/L_xL_y. We satisfy these two requirements by choosing M_x = L_y, M_y = L_x.

Working with the magnetic Bloch states

$$\begin{equation} \Psi_{\mathbf{k}} ({\mathbf{r}}) = \sum_{\mathbf{R}} \mathrm{e}^{-\mathrm{i} \mathbf{k} \cdot \mathbf{R}} \; \hat{T}_{\mathbf{R}} \Psi ({\mathbf{r}}) \hat{T}_{\mathbf{R}}^{-1} \end{equation} \tag{ 9 }$$

allows us to block diagonalize the Hamiltonian

$$\begin{equation} \mathcal{H}^{\rm BdG} = \frac{1}{M_{x} M_{y}} \sum_{\mathbf{k}}\sum_{\mathbf{r}, \mathbf{r}' } \Psi_{\mathbf{k}}^\dagger ( \mathbf{r} ) H_{\mathbf{k}} (\mathbf{r}, \mathbf{r}') \Psi_{\mathbf{k}} ( \mathbf{r}' ). \end{equation} \tag{ 10 }$$

The indices k and r from here on are defined in the magnetic Brillouin zone and magnetic unit cell, respectively, that is

$$\begin{equation} \mathbf{k} = \left( 2\pi \frac{m_{x}}{L_{x}M_{x}},2\pi \frac{m_{y}}{L_{y}M_{y}} \right),\quad m_{\alpha} = 0 \cdots M_{\alpha}-1, \end{equation} \tag{ 11a }$$

$$\begin{equation} \mathbf{r}= \left( \ell_{x}, \ell_{y} \right), \quad \ell_{\alpha} = 0 \cdots L_{\alpha}-1. \end{equation} \tag{ 11b }$$

$\mathcal {H}^{\mathrm { BdG}}$

Figure 2 shows the near-vortex LDOS calculated by diagonalizing $\mathcal {H}^{\mathrm {BdG}}$ of equation (1) with fixed gap functions as given by equation (3) on a system of dimension (N_x,N_y) = (301,301). We choose realistic values of the parameters for the coherence length ξ = 16.4a₀ [31, 32], as well as gap values Δ⁰_s = 3 meV for on-site pairing and Δ⁰_s^± = 1.5 meV for NNN pairing [21, 22, 30].

Zoom In Zoom Out Reset image size

We can interpret the vortex bound states in this non-self-consistent BdG calculation as bound states in a potential well given by equation (3), where only states around the normal-state Fermi surface constitute the bound states. There are then two sources of anisotropy: anisotropic, quasi-one-dimensional low-energy properties of the normal state, and an anisotropic gap, both defined in the momentum space. The geometric distribution of LDOS will be dominated by one or the other source of anisotropy at different energies.

At low energies, the normal state properties dominate the distribution of LDOS (figures 2(a) and (b)). Hence irrespective of pairing structure, the bound state is located at the center of the potential well. Since the Bloch states making up this bound state have two main velocities due to the quasi-one-dimensional parts of the Fermi surface, the bound state mainly extends in these two directions out of the well, resulting in the rays in figures 2(a) and (b). The gap is suppressed near the vortex center, and its anisotropy is of little importance. Hence the flat (quasi-one-dimensional) parts of the electronic structure in figure 1(a) (solid line) dominate over the small anisotropy of the s^± gap (see figure 1(b)). For a better comparison with experiment, we present results of reduced spatial resolution by Gaussian filtering (σ = 3a₀) in the insets. The low resolution result is consistent with results of the quasi-classical analysis by Wang et al [14] and in good agreement with experiment shown in figure 2(c).

At higher energies on the other hand, the bound state is located away from the vortex core. The quasi-one-dimensionality of the Fermi surface allows for localization in one direction and extension in the other. This leads to a square-like inner ring in the LDOS for both pairings (figures 2(e) and (f)). The difference, however, results from the anisotropy of the gap function. While the isotropic s-wave gap is analogous to a potential that is independent of momentum, the anisotropic gap is one for which different states around the Fermi surface experience different potentials depending on their momenta. With the gap function of s^± form, the quasi-one-dimensional portion of the Fermi surface experiences a stronger trap potential, leading to a suppression of its contribution to the bound-state wave function. As a result, the bound state exhibits pronounced isolated segments, 'hot spots', within the inner ring that point in the Fe–Fe direction, as shown in figure 2(f). We again Gaussian filter the images and show them in the insets. Note the 'hot spot' are robust and even more pronounced in the low resolution insets in figure 2(f) in good agreement with the experimental data figure 2(g).

We now turn to the LDOS at the core of the vortex and its particle–hole asymmetry. This turns out to be largely insensitive to anisotropy of pairing. The LDOS at the core of the vortex for the on-site pairing shown in figure 2(d) exhibits particle–hole asymmetry with the highest peak at negative energy. Such asymmetry appears in the so-called 'quantum-limit' vortex bound state [33], whose highest LDOS peak is at energy Δ²/2E_F above(below) the Fermi energy for an electron(hole)-like band, where E_F is the energy difference between the Fermi energy and the bottom(top) of the band. The energy of the LDOS peak being 0.05 meV below the Fermi energy is expected given E_F = 98 meV and Δ = 3 meV within our input bandstructure. Though similar particle–hole asymmetry has been observed in [8] (see figure 2(h)) the energy at which the peak was observed suggests that other hole pockets with larger gap values may be responsible.

Now, we check whether the single-band model is sufficient to describe vortex bound states within the energy range of interest. A simple insight can be gained by treating each band independently and estimating the energy of its lowest bound state to be Δ²/2E_F following Caroli et al [33] for the gap size Δ and the Fermi energy E_F specific to each band. Using measured Fermi energies and gap parameters [10, 21, 22, 30], we estimate the energies of the lowest bound states of the γ pocket and the electron pockets to be of the same order. However, the lowest bound state energies of the two smaller hole pockets are an order of magnitude larger. This rough estimate implies that the LDOS within the energy below 1 meV should be dominated by bound states coming from the γ band and those coming from the two electron bands. If indeed each bound state comes from a single band, we expect to find bound states with LDOS distribution resembling what we predicted in section 4.1.

For concreteness, we carry out a non-self-consistent BdG calculation using the band structure given by [20] with five bands. Unfortunately, the γ-pocket Fermi surface of this band structure (dashed line in figure 1(a)) is far more isotropic compared to what has been measured in [30] and guided the band structure we use in the rest of this paper. Hence we do not expect as pronounced 'ray' features at low energies compared to what is shown in figure 2 from our (single-band) calculations and experiment. Another issue we face with a five-band calculation is the limitation on the accessible system size. For a system of size (101,101), we impose NNN pairing that is trivial in the orbital space having magnitude Δ_s^± = 15 meV in order to fit the vortex bound states within the system and minimize the boundary effect. As in the single-band calculation, we create a vortex at the center of the form given in equation (3), however with ξ = 10a₀. Figure 3 shows the resulting LDOS at different bound state energies. At the lowest energy there is no clear sign of 'rays' though a small amount of anisotropy is still present, as expected from the smaller γ-band anisotropy (see figure 3(a)). Figures 3(b) and (c) show typical LDOS images of vortex bound states at higher energies. Figure 3(b) looks very different from the LDOS distribution obtained in section 4.1 and we hence assign the corresponding bound state to the electron pockets. However, the LDOS shown in figure 3(c) shows the same 'hot spots' as obtained within our single-band calculation and shown in figure 2(f). Focusing on the γ band should thus suffice to capture the features observed in [8].

Zoom In Zoom Out Reset image size

Figure 4 shows the results from the (single-band) self-consistent calculation. We compare two pairing interactions—on-site attraction U = −0.35 eV, and NNN attraction V ' = −0.3 eV—for a system with magnetic unit cell of dimensions (L_x,L_y) = (19,38). This corresponds to a full lattice size of (N_x,N_y) = (38 × 19,19 × 38). In zero field, the two cases lead to a uniform superconducting gap of Δ⁰_s = 27 meV and Δ⁰_s^± = 10 meV, respectively. We have chosen U and V ' such that the coherence length ξ∝Δ⁻¹ is small compared to the inter-vortex spacing. This allows us to focus on a nearly isolated vortex within the computationally feasible size of the magnetic unit cell. Although the resulting gap values are an order of magnitude larger than what is known experimentally, this should not affect the validity of the results in a qualitative manner. Both at low energy and at higher energy close to the gap value, we observe features that qualitatively agree with the results obtained in the previous section.

Zoom In Zoom Out Reset image size

The self-consistent calculation also allows us to study the local order parameters of a given structure near a vortex. Unlike for the on-site attraction, where Δ_s(r) is the only allowed gap function, order parameters of different symmetries can mix near a vortex for NNN attraction. A near-vortex map of Δ_s(r) for on-site pairing shown in figure 5(a) indeed shows almost isotropic healing of the order parameter away from the vortex core. However, for the NNN attraction which leads to uniform s^±-wave pairing in zero-field, the secondary order parameter Δ_dxy(r) is induced near the vortex. Coupling between this secondary order parameter and the primary Δ_s^±(r) leads to a strong angular variation of both components as can be seen in figures 5(b) and (c).

Zoom In Zoom Out Reset image size

To gain further insight into the admixing of a secondary order parameter near a vortex for the anisotropic pairing, we analyze the Ginzburg–Landau free-energy density. The free-energy density for s-wave and d-wave order parameters reads

$$\begin{eqnarray} &&\fl f = \alpha_{s} |s|^{2} + \alpha_{d} |d|^{2} + \beta_{1} |s|^{4} + \beta_{2} |d|^{4} + \beta_{3} |s|^{2} |d|^{2} + \beta_{4} ( s^{*2} d^{2} + \mathrm{c.c.}) + \gamma_{s} | \vec{D} s |^{2} + \gamma_{d} | \vec{D} d |^{2} \nonumber\\ &&+ \gamma_{v} ( D_{x} s D_{y} d^{*} + D_{y} s D_{x} d^{*} + \mathrm{c.c.}), \end{eqnarray} \tag{ 12 }$$

where s and d are shorthands for s(r) and d(r), the order-parameter fields for the s^± and d_xy gaps, respectively, and D_i = ∂_i − ieA_i is the covariant derivative. The fields s(r) and d(r) can be thought of as Δ_s/s^±(r) and Δ_dxy(r) after coarse graining. This type of admixing near a vortex has previously been studied in the context of cuprates, leading to the prediction of a fourfold-anisotropic order parameter around a vortex [15–18].⁴ As the large halo around vortices in cuprates [34] hindered the observation of this admixing, LiFeAs presents an opportunity for this observation.

The spatial variation of the secondary component d_xy in figure 5(c) is largely due to the derivative coupling, the term proportional to γ_v in equation (12). This intermixing term is expected to be large when the s-wave order parameter is of s^± type, since the same NNN pairing interaction is responsible for both s-wave and d-wave order parameter. For |s| ≫ |d| and $|\vec {D} s| \gg |\vec {D} d|$ the spatial structure of the d_xy component is determined largely by the structure of the s-wave component. Minimizing equation (12) with respect to d(r) and keeping only terms up to linear order in d(r), we find

$$\begin{equation} -\gamma_{d} \vec{D}^{2} d+\alpha_{d} d + \beta_{3} |s|^2 d + \beta_{4} s^2 d^{*}= \gamma_{v} ( D_{x} D_{y} + D_{y} D_{x} ) s. \end{equation} \noindent \tag{ 13 }$$

Hence, the curvature in the leading s-wave component will induce the secondary (d_xy) component. Now, consider a single isolated vortex. As s(r) is recovered at the length scale of the coherence length ξ away from the core of the vortex, we expect a large d(r) due to coupling to the large curvature of s(r) at this distance. Since ξ = ℏv_F/πΔ ∼ 3.0a₀ for the uniform gap value with V ' = −0.3eV, this is in agreement with the positions of the maxima of d(r) in figure 5(c) as a function of |r| setting the vortex core at the origin. We can also explain the angular variation and the form of the anisotropy of d(r) in this framework. If we assume s(r) = f(r)e^iθ with a slowly changing f(r) and the azimuthal angle θ measured from the Fe–Fe direction, we find from equation (13)

$$\begin{equation} d({\mathbf{r}}) \sim \partial_{x} \partial_{y} s({\mathbf{r}})\sim \,\mathrm{e}^{-\mathrm{i} \theta} (1 + 3 \,\mathrm{e}^{4 \mathrm{i} \theta}), \end{equation} \noindent \tag{ 14 }$$

ignoring the phase due to the magnetic field. The structure of the derivative hence gives rise to a four-fold anisotropy, which explains the fact that |d(r)| is maximum in the Fe–Fe direction, while it is suppressed along the 45° direction. Coupling to d(r) gives then in turn cause for the four-fold anisotropy in s(r).