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

A major question in Fe-based superconductors remains the structure of the pairing, in particular whether it is of unconventional nature. The electronic structure near vortices can serve as a platform for phase-sensitive measurements to answer this question. By solving Bogoliubov-de Gennes equations for LiFeAs, we calculate the energy-dependent local electronic structure near a vortex for different nodeless gap-structure possibilities. At low energies, the local density of states (LDOS) around a vortex is determined by the normal-state electronic structure. However, at energies closer to the gap value, the LDOS can distinguish an anisotropic from a conventional isotropic s-wave gap. We show within our self-consistent calculation that in addition, the local gap profile differs between a conventional and an unconventional pairing. We explain this through admixing of a secondary order parameter within Ginzburg-Landau theory. In-field scanning tunneling spectroscopy near vortices can therefore be used as a real-space probe of the gap structure.


Introduction
The gap structure in the Fe-based superconductors and its possible unconventional nature is still a key issue in the field four years after their discovery. In most compounds, the pairing is believed to be of the so-called s ± type, for which the order parameter changes sign between the electron-like and the hole-like Fermi surfaces [1,2]. Some experimental results are consistent with this prediction [3][4][5][6][7]. However, a major difficulty in distinguishing such an unconventional pairing state from a trivial s-wave gap is that both states are nodeless and transform trivially under all the symmetry operations of the material's point group. As the experimental probes that are usually used to distinguish various gap structures, such as phase-sensitive probes, are not Fermi pocket specific, an unambiguous evidence of the unconventional s ± pairing remains evasive.
One route to accessing phase information using a phase-insensitive probe would be through vortex bound states, as a vortex introduces a spatial texture to the superconducting order parameter.
Advancements in in-field scanning tunneling spectroscopy (STS) have enabled the study of vortex bound states. Indeed, a recent STS experiment on LiFeAs under a magnetic field has shown an intriguing energy dependence in the spatial distribution of the local density of states (LDOS) near a vortex [8]. The remaining question is whether the observed LDOS distribution near vortex can be instrumental in selecting one of the proposed gap structures: s ± -wave [1], s ++ -wave [9,10], and (spin-triplet) p-wave [11,12]. At zero bias, the LDOS shows a fourfold star shape with high-intensity 'rays' along the Fe-As direction. Similar features in NbSe 2 [13] were interpreted as a sign of gap minima along this direction. However, a quasi-classical analysis by Wang et al. [14] pointed out that the normal-state band structure of LiFeAs -namely a highly anisotropic hole pocket around the Γ pointcould be producing these rays irrespective of gap structure. By contrast, little attention has been given to the high energy LDOS distribution observed in Ref. [8]: hot spots appearing at the intersection of split rays.
Motivated by these observations, we present a study of the near-vortex electronic structure and signatures of unconventional pairing therein within the Bogoliubov-de Gennes (BdG) framework. By (non-self-consistently) imposing a gap structure and solving the BdG Hamiltonian, we first show that the isotropic s-wave and s ± -wave pairing result in different spatial distributions of the LDOS at energies approaching the gap value. In particular, we find s ± -wave pairing to yield the observed hot spots. Then we solve the BdG equations self-consistently, and based on our results propose detecting the spatial distribution of the gap around a vortex for a more direct evidence of unconventional s ± -wave pairing. A vortex not only suppresses the order-parameter amplitude at its core and introduces a singular point in space around which the phase of the order parameter winds, but it also induces a secondary order parameter in its vicinity [15][16][17][18][19]. Due to the induced secondary order parameter near the vortex, the gap recovery should show a strong angular dependence. Detection of such anisotropy will be an unambiguous evidence of unconventional pairing.  [20]. For the most part of this work, we focus on the γ band that is around the Γ point, whose Fermi surface is shown as a solid line. (b) Sketch of the three gap functions with s-, s ± -, and d xy -wave momentum structure around the γ-band Fermi surface.
The remainder of this paper is organized as follows: In sections 2 and 3, we introduce the microscopic model and describe the Bogoliubov-de Gennes calculations, respectively. In section 4, we present the results of the BdG calculations and discuss them within Ginzburg-Landau theory. In section 5, we summarize our findings and remark on future directions. Throughout the paper we focus on the large hole pocket and study the single band problem. However, we also present results from non-self-consistent BdG calculations on a five-band model in section 4, which show good agreement with observations from single-band model calculations in the energy range of our interest.

Model
We describe LiFeAs in the superconducting state with the (mean-field) BdG Hamiltonian Here, Ψ i ≡ (c i↑ , c † i↓ ) T is a Nambu spinor, and c is (c † is ) annihilates (creates) an electron at lattice site i with spin s within a single-band model for the large hole pocket around the Γ point: the so-called γ band. However, Eq. (1) can easily be generalized for a multiband model. In this paper, we focus on the single-band model for the most part since the superconducting gap is the smallest on the γ band [21,22] and hence we expect low energy physics to be dominated by this band. Moreover, this band mainly stems from the (in-plane) d xy orbitals, and thus shows little k z dependence [23]. It is therefore a natural choice for LiFeAs. Note that previous BdG calculations on different Fe-pnictides focused on two-band models for the d xz / d yz orbitals [24][25][26][27][28][29]. Our choice of the hopping matrix t ij is guided by the experimental observations on the γ pocket [21,22,30] to be t = −0.25eV for nearest-neighbor hopping, t ′ = 0.082eV for next-nearest-neighbor hopping, and t ii = µ = 0.57eV for the chemical potential. Figure 1(a) shows the resulting Fermi surface in solid red line. Though we stay within this single-band model for the self-consistent BdG studies, we have also used a five-band model for the non-selfconsistent calculation with tight-binding parameters from Ref. [20] to test the validity of focussing on the γ band for the energy range of our interest (see section 4.2). Figure  1 The ∆ ij are the (bond) gap functions. For a self-consistent solution of H BdG , we require the gap functions to satisfy where V ij < 0 is the attractive interaction between sites i and j in the singlet channel, and · denotes the thermal expectation value. Restricting the interaction V ij to a specific form constrains the momentum structure of the gap function, since ∆ ij = 0 only if V ij = 0. In the uniform case, an on-site attraction V ij = Uδ ij leads to a spin-singlet s-wave gap ∆(k) = ∆ 0 s , while a next-nearest-neighbor (NNN) attraction V ij = V ′ δ i,j allows for the singlet gap functions of s ± form, ∆(k) = 4∆ 0 s ± cos k x cos k y , and d xy form, ∆(k) = 4∆ 0 dxy sin k x sin k y . Figure 1(b) shows sketches of these gap functions. We restrict our calculations in the following to these 'pure' gap structures. Even though the true gap function is a (symmetry-allowed) mixture of such gap functions, the dominant channel (on-site or NNN interactions) will determine whether an s ± -or an s ++ -wave gap is realized in the presence of the electron pockets.
For the non-self-consistent BdG study, the vortex will be imposed through the gap-function configuration of where the vector r ij points to the midpoint of sites i and j, and θ ij is its azimuthal angle measured from the Fe-Fe direction. This corresponds to a single vortex located at the origin suppressing locally the order-parameter amplitude. In addition, the orderparameter phase winds around the vortex core. For the self-consistent BdG study, we induce the vortices by applying a magnetic field Hẑ. Assuming minimal coupling between an electron and the field, the hopping between sites i and j acquires a Peierls phase where Φ 0 = h/2e is the magnetic fluxoid and r i is the vector pointing to the site i. We assume a uniform magnetic field H and write the vector potential in the Landau gauge A(r) = −Hyx. From the self-consistent solution ∆ ij , we can define local gap order parameters of different symmetries. For an on-site interaction, the local s-wave order parameter is defined as ∆ s (r) = ∆ r,r . Note that from here on, we use r without any site index to denote both a lattice site and the vector pointing to it in units of the lattice constant a 0 . With NNN interaction, we define local order parameters of s ± form and d xy form where ∆ rr ′ ≡ ∆ rr ′ exp[−iϕ(r, r ′ )] ensures that order parameters of different symmetries do not mix under magnetic translations. Note that for the uniform case, ∆ s (r) = ∆ 0 s , ∆ s ± (r) = ∆ 0 s ± , and ∆ dxy (r) = ∆ 0 dxy as defined above.

Method
In this section, we elaborate on our two approaches to solve the BdG equations and obtain the LDOS near a vortex. For both, diagonalizing the Hamiltonian H BdG in Eq. (1) for a system of size (N x , N y ) is computationally the most expensive part.

Non-Self-Consistent Approach
For the non-self-consistent calculation, we impose a gap function in the form given by Eq.
(3) and find the low lying eigenvalues and eigenstates of H BdG using the Lanczos algorithm ‡. The LDOS can be calculated from the eigenenergies E n and eigenstates [u n (r), v n (r)] as 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.

Self-Consistent Approach
For the self-consistent calculation, we assume initial gap functions and use the eigenvalues and eigenvectors of Eq. (1) to calculate the gap functions given by Eq. (2). We proceed iteratively until self-consistency is achieved. In diagonalizing H 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 ‡ We suppress low energy states from forming at the boundary by imposing an on-site potential of 10 eV to the sites at the boundary.
the ordinary lattice translation operator T R . However, the kinetic part commutes with the magnetic translation operatorT for a magnetic lattice vector R whose unit cell contains two magnetic fluxoids. The pairing term in general does not commute withT R . Nevertheless, when vortices form a lattice,T 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 xx and L yŷ , such that R forms a rectangular lattice R = (m x L x , m y L y ), where m α = 0 · · · M α − 1 and M α ≡ N α /L α §. Note that periodic boundary conditions in the Landau gauge A(r) = −Hyx require the total magnetic flux through the system to be an integer multiple of 2Φ 0 N x . In addition, one magnetic unit cell contains a magnetic flux of 2Φ 0 , i.e. H = 2Φ 0 /L x L y . We satisfy these two requirements by Working with the magnetic Bloch states allows us to block diagonalize the Hamiltonian The indices k and r from here on are defined in the magnetic Brillouin zone and magnetic unit cell, respectively, that is (10), we can compute the eigenstates and eigenenergies of H BdG . These are then used to calculate ∆ ij with Eq. (2) closing the self-consistency cycle. Finally, we use the self-consistent solution ∆ ij to calculate the local order parameters of s-, s ± -and d xy -wave symmetry and also the LDOS of the electronic degrees of freedom, as defined in Eq. (7).  We choose realistic values of the parameters for the coherence length ξ = 16.4a 0 [31,32], as well as gap values ∆ 0 s = 3meV for on-site pairing and ∆ 0 s ± = 1.5meV for NNN pairing [21,22,30].

Non-Self-Consistent Approach on Single Band Model
We can interpret the vortex bound states in this non-self-consistent BdG calculation as bound states in a potential well given by Eq. (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 [ Figs. 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 Figs. 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 Fig. 1(a) (solid line) dominate over the small anisotropy of the s ± gap [see Fig. 1(b)]. For a better comparison with experiment, we present results of reduced spatial resolution by gaussian filtering (σ = 3a 0 ) 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 Fig. 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 [Figs. 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 Fig. 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 Fig. 2(f) in good agreement with the experimental data Fig. 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 Fig 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 ∆ 2 /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.05meV below the Fermi energy is expected given E F = 98meV and ∆ = 3meV within our input bandstructure. Though similar particle-hole asymmetry has been observed in Ref. [8] [see Fig. 2(h)] the energy at which the peak was observed suggests that other hole pockets with larger gap values may be responsible.

Non-Self-Consistent Approach on Five Band Model
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 ∆ 2 /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 1meV 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 Ref. [20] with five bands. Unfortunately, the γ-pocket Fermi surface of this band structure [dashed line in Fig 1(a)] is far more isotropic compared to what has been measured in Ref. [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 Fig. 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 ± = 15meV 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 Eq. (3), however with ξ = 10a 0 . 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 Fig. 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 Fig. 3(c) shows the same 'hot spots' as obtained within our single-band calculation and shown in Fig. 2(f).
Focussing on the γ band should thus suffice to capture the features observed in Ref. [8].  , (c). The equal-amplitude contours in red go from 0.825 for the innermost to 0.925 for the outermost contours (after normalization) with equal intervals between the contours in between. The insets again indicate the structure of the local order parameter. Note that the color-scale for ∆ dxy (r) is much smaller than for ∆ s (r).

Self-Consistent Approach
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 ∆ 0 s = 27meV and ∆ 0 s ± = 10meV, respectively. We have chosen U and V ′ such that the coherence length ξ ∝ ∆ −1 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.
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 Fig. 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 Figs. 5(b) and (c).
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 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][16][17][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 Fig. 5(c) is largely due to the derivative coupling, the term proportional to γ v in Eq. (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 reponsible for both s-wave and d-wave order parameter. For |s| ≫ |d| and | Ds| ≫ | Dd| the spatial structure of the d xy component is determined largely by the structure of the s-wave component. Minimizing Eq. (12) with respect to d(r) and keeping only terms up to linear order in d(r), we find 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 ξ =hv F /π∆ ∼ 3.0a 0 for the uniform gap value with V ′ = −0.3eV, this is in agreement with the positions of the maxima of d(r) in Fig. 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 The microscopic model we consider is related to the single-band model of cuprates through rotation by 45 • , the roles played by s-wave and d-wave order parameters are reversed and our d-wave order parameter is of d xy form rather than d x 2 −y 2 . angle θ measured from the Fe-Fe direction, we find from Eq. (13) 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).

Conclusion
We have contrasted the effects of anisotropic s ± -wave (NNN) pairing and isotropic swave (on-site) pairing on the near-vortex local density of states in LiFeAs by solving Bogoliubov-de Gennes equations both non-self-consistently and self-consistently. We have found qualitative changes in the geometric distribution of the density of states as a function of energy. At low energies, the anisotropy of the vortex bound state, and hence the LDOS, is determined by the normal state low energy electronic structure, independent of the gap structure. Different pairing structures, however, lead to qualitatively different LDOS distributions at higher energies: While the isotropic swave shows a square-like feature of roughly equal intensity, four 'hot spots' develop in the case of an (anisotropic) s ± -wave gap. Indeed, our results for the latter case qualitatively agree with recent experiments [8].
From the self-consistent treatment we have further found a difference in the recovery of the order parameter away from the vortex core: a pronounced angular dependence of the s ± -wave gap compared to isotropic behavior for the s-wave gap. Employing a Ginzburg-Landau analysis, we have explained this difference through admixing of a secondary order parameter supported by the NNN interaction. Note that such intermixing is negligible for an s-wave pairing with a dominant on-site pairing interaction, as no other pairing instabilities are nearby. For the NNN interaction, however, s ± -and d xy -wave instabilities have comparable transition temperatures. Detection of the anisotropy or even the secondary order parameter would be a strong proof of the unconventional nature of the pairing.
In this work, we focused on the γ band with interest in low energy properties, as this is the band with the smallest gap [21,22]. Hence, for features at energies less than the gap scale, we expect our calculation to capture salient features of in-field STS experiments. The comparison between the calculated LDOS for the single-and the five-band models and the results in Ref. [8] supports this conjecture.
In closing we note that our calculation captures Friedel-like oscillations, frequently referred to as quasi-particle interference (QPI), due to vortices. QPI in the presence of vortices was successfully used to access phase information with STS in cuprates [35]. Recent in-field QPI experiments on FeSe have been interpreted to be consistent with an s ± scenario when a vortex is treated as a magnetic scatterer for BdG quasiparticles [4]. However, a vortex is at once a point of gap suppression, a point with magnetic flux, and a point around which the order-parameter phase winds. While we treated vortices faithfully in the self-consistent calculation, we could not investigate effects of interpocket sign change as we only considered one pocket. An extension of the present work with the full band structure would be necessary to work out what to expect for different order-parameter possibilities, especially how the phase difference between different pockets affects in-field QPI.