Single vortex structure in two models of iron pnictide $s^\pm$ superconductivity

The structure of a single vortex in a FeAs superconductor is studied in the framework of two formulations of superconductivity for the recently proposed sign-reversed $s$ wave ($s^\pm$) scenario: {\it (i)} a continuum model taking into account the existence of an electron and a hole band with a repulsive local interaction between the two; {\it (ii)} a lattice tight-binding model with two orbitals per unit cell and a next-nearest-neighbour attractive interaction. In the first model, the local density of states (LDOS) at the vortex centre, as a function of energy, exhibits a peak at the Fermi level, while in the second model such LDOS peak is deviated from the Fermi level and its energy depends on band filling. An impurity located outside the vortex core has little effect on the LDOS peak, but an impurity close to the vortex core can almost suppress it and modify its position.


Introduction
The discovery of superconductivity in iron based materials at relatively high temperatures quickly spurred intensive theoretical and experimental work [1]. The electron-phonon interaction has early been ruled out as the pairing mechanism [2] and spin fluctuation mediated interactions have been proposed [3,4,5]. The most studied materials can be divided in three families: the "lanthanum" family LaFeAsO 1−x F x , the "rare-earth" family (Sm, Nd, Ce, Pr)FeAsO 1−x F x , and the (Ba, Ca, Sr)Fe 2 As 2 family (also known as "122" compounds).
Iron pnictide superconductors have a complex band structure, with a Fermi surface (FS) consisting of four sheets, two of them hole-like, and the other two electron-like [3,6,7,8]. S-wave, d-wave and p-wave pairing scenarios have early been proposed to describe the superconducting state [9,10,11,12], as well as extended s-wave [13]. More recently, two possible pairing scenarios have been in dispute: the "sign-reversed s-wave state" (s ± -state), where the gap function has no nodes on any of the FS sheets but takes on opposite signs on the hole-like and electron-like sheets of the FS [3,4]; the extended s-wave symmetry, in which the gap function has four nodes on the electron FS sheets, no nodes on the hole sheets, and has the same sign on both hole FS's [5]. The existing experimental evidence for lanthanum based superconductors seems to be mostly consistent with unconventional gap functions with nodes on the FS, while the experiments with superconductors containing rare-earth atoms as well as with the 122 compounds tend to suggest a nodeless nearly isotropic gap across the FS. However, Andreev reflection studies have lead to different conclusions regarding the pairing symmetry for the same (rare-earth based) material [14,15].
The mixed phase of the iron pnictide superconductors has been less studied. Specific heat measurements in LaFeAsO 1−x F x under a magnetic field revealed a √ H dependence of the specific heat coefficient, which has been interpreted as due to the Doppler shift of low energy excitations (close to the nodes of the gap function) caused by the superfluid flow far from the vortex cores [16,17,18]. More recently, a scanning tunnelling microscope (STM) observation of vortices in BaFe 1.8 Co 0.2 As 2 has been made [19], where the differential conductance versus bias voltage was measured inside and outside the vortex cores. On the sample surface scanned by the STM, the gap value (6 meV) was almost uniform across the surface, with a fractional variance of 12%, and the vortex locations were not coincident with the impurities present in the material. The observed differential conductance inside the vortex cores displayed no peaks signalling the presence of core bound states, only the V-shaped background.
The s ± superconducting state, because it is novel, has been the subject of intense theoretical investigation. Such studies often predict unconventional behaviour of various physical observables, qualitatively different from the usual s-wave state. This is because the gap function, despite having no nodes on the FS sheets, changes sign between them, effectively leading to a Josephson-like coupling between two s-wave gap functions with a phase difference of π, and also because of quantum interference effects between Bloch waves belonging to different FS pockets.
In this work we study the single vortex structure in the s ± state, using two different formulations. In the continuum formulation, only two FS sheets are considered, having parabolic dispersions with opposite masses. In the discrete formulation, a tight-binding model that reproduces the band structure close to the Fermi level is used. In this model, we find the lowest energy (subgap) excitations to be spatially extended, in contrast to the case with conventional s-wave vortices and the results of the continuum formulation, where the excitations are localized near the vortex core. Only higher energy excitations (but still subgap) are localized. This effect is a manifestation of interband interference effects which have earlier been found to produce Andreev bound states at finite energy [20]. We also study the effect of impurities, close to the vortex core, on the local density of states and find that they have a strong effect on local density of states at the vortex core.

Model
We begin by considering a simplified model for the FeAs as described in ref. [21]. The Hamiltonian is of the form The model describes two bands, one of positive mass particles (electrons),ê k,σ , with dispersion ǫ e (k) and another of negative mass particles (holes),ĥ k,σ , with dispersion ǫ h (k) in two dimensions. The hole band is centred at the origin of the Brillouin zone (Γ point), while the electron band is centred at the M point, (π, π). It is then assumed that the two types of particles interact via pair intraband couplings (V hh , V ee ) and pair interband (or Josephson) couplings (V he , V eh ). These interactions are supposed to have a magnetic origin and are taken as repulsive [21,22]. In order to do a continuum formulation of the vortex problem in real space, we take the h band as parabolic with negative mass and the e band as a parabola of positive mass. The bottom of the e band is located at E e 0 = −0.58, the top of the h band is located at E h 0 = 0.24 (the energies are in eV units). The chemical potential is located at E F = 0. The masses of the two bands are taken as m e = 1, m h = −1.8. Since the h particles have negative mass, it is convenient to make a particle-hole transformation asĥ σ (r) →h † −σ (r). We also perform the gauge transformationē σ (r) → e iQ·rê σ (r), with Q = (π, π), on the electrons. A standard mean field decoupling of the interaction terms leads to the Bogolubov-de Gennes (BdG) equations for the wave functions of the quasiparticles [23] wherep = −i ∇ in two dimensions. Here,ū e (r),v e (r) are the usual BdG components of the wave function for the electronsê andũ h (r),ṽ h (r) are the wave function components for the holesh (which now have positive mass). The self-consistency of the mean-field approach implies that where f (E) is the Fermi function. We have approximated the pair couplings by two constants V 1 = V hh = V ee and V 2 = V he = V eh . We may consider the insertion of a magnetic field in the usual minimal coupling way taking p → p − q c A, where A is the vector potential. The charge q is negative forê electrons (q = e < 0) and positive for theh holes (q = −e). The vector potential is given by Maxwell's equations which, in the Coulomb gauge ( ∇.A = 0 ), is given by The current density originated in the supercurrents is obtained self-consistently by We write the order parameters using polar coordinates in the form where ρ and ϕ and the radial and angular polar coordinates. Such a two-component vortex has magnetic flux equal to one flux quantum ( Φ = Φ 0 = hc 2|e| ): this is because the diamagnetic currents from electrons and holes add up, screening the external magnetic field to one flux quantum [24]. The system is placed in a cylinder of radius R. Given the azimuthal symmetry of the system, neither ∆(ρ) nor A depend on ϕ. Therefore the Hamiltonian may be diagonalized separately for each value of the angular momentum, µ, which becomes a good quantum number. Following Ref. [25], the wave functionsū n andṽ n are written in the form where the radial functions are expanded in a way similar to ref. [25,26] using as basis functions The functions J m are the cylindrical Bessel functions and α jm is the j th zero of the Bessel function of order m. The set of values of the angular momentum is given by To illustrate the procedure we consider the standard case of the e bands. For each eigenvalue E n , we have a single value of µ and it is enough to diagonalize the matrix, defined in the subspace of the zeros of the Bessel function, where with and Also we have In the case of the h band the angular momenta are interchanged. It is important to note that the symmetry of the BdG equations u n (r) → v * n (r), v n (r) → −u * n (r), E n → −E n allows to reduce the solution to the positive values of µ. We obtain the eigenvectors and eigenvalues for negative values of µ using the above symmetry. The vector potential is also solved self-consistently, using Poisson's equation, as explained in [26]. From the selfconsistent solution of the BdG equations we obtain the energy spectrum, the magnetic flux, the current, the magnetic field and the order parameter profiles as a function of distance from the vortex core. In general the contribution of the vector potential in the BdG equations is negligible but not in the calculation of the supercurrents. Also, to obtain the exponential decay of the magnetic field we have to take into account its effect carefully [25,26]. The results are shown below.
We also compute the local density of states (LDOS) at site r and energy E, defined as for the e electrons, and for the holes. The local density of states is observable through STM (for a recent review see ref. [27]).

Results
In Figure 1 we show the energy spectrum for the two bands for a characteristic set of parameters. We take V 1 = 0.6 and V 2 = 0.07. In the case of the electron (e) band there is a bound state at positive energies for each value of the angular momentum. These are the well known Caroli-de Gennes-Matricon states [23]. The other states are extended.   In the case of the hole band (h) band the bound states appear at negative energies (or at positive energies for negative values of the angular momentum). The spectrum is not symmetric around the Fermi level, as expected [25]. The bound states branch is symmetric with respect to the sign of the angular momentum.
In Figs. 2, 3 we show the spatial dependence of the radial part of the gap functions    for the electron and the hole bands. As shown before, the two gap functions have opposite signs, characteristic of the sign-reversed s-wave scenario. At short distances and low temperatures the gap functions oscillate due to the states inside the core [28]. As the temperature grows the oscillations decrease in size, and the size of the core increases (Kramer-Pesch effect [29]). The highest temperature considered is of the order of half the critical temperature.
In Fig. 4 we show the magnetic flux enclosed by the vortex as a function of the distance to the vortex core, saturating to a quantum of flux, as determined by the choice of the vorticity of the gap functions [26,30,31]. Fig. 5 shows the exponential decay of the absolute value of the magnetic field.
The current is the result of the electron and hole contributions. This is shown in Fig. 6. The superfluid current shields the external field that penetrates through the vortex. Both contributions to the supercurrent have the same (positive) sign.
Since the gap function of the holes is smaller, the oscillations at small temperatures are larger. Also, since the energy scale where the gap function tends to the bulk value is smaller, the spatial length scale is larger in the sense that the fluctuations extend to larger distances. This also affects the oscillations of the current, which occur, at low temperatures, at two length scales, one associated with the e band and the other with the h band.
In Fig. 7 we show the LDOS at a set of points whose distance from the vortex core increases. First we consider a point at the vortex core (simulating the location of the tip of the STM on top of the vortex core location). We see that the contributions from the electrons and the holes are not symmetric around the zero bias, reflecting the existence of bound states at positive and negative energies, respectively. However, the total LDOS resulting from the sum of the two contributions is quite symmetric. We find, therefore, that in this simple model there is a zero bias peak at the vortex core. As we move further from the vortex core position, the central peak splits, leading to two contributions that clearly separate. This is due to the "coherence" peaks which reflect the size of the gap function as we move away from the vortex core and approximate the LDOS of the bulk. The LDOS has a gap which is due to the smaller of the gaps (in this case originates in the h band) and shows a double peak structure due to the different h and e band gaps.

Model
A recent tight-binding model [32] for the band structure of iron pnictide superconductors assumes two orbitals per unit cell, d xz and d yz , in a square lattice (we take the lattice constant equal to unity). The tight-binding Hamiltonian in real space can be written aŝ where σ =↑, ↓ denotes spin projections, the indices i, j denote lattice sites, ν = x, y denotes the atomic orbitals d xz and d yz at a lattice site, and {ν,ν} = {x, y}. The first neighbour hoppings along the x direction are t ijσ,x = t 1 , t ijσ,y = t 2 ; first neighbour hoppings along the y direction are t ijσ,x = t 2 , t ijσ,y = t 1 . The second neighbour hoppings It has been pointed out that the Fermi velocity in the electron Fermi surface sheets becomes underestimated in model (23), however [5]. The above Hamiltonian contains   in momentum space ∆(k) = ∆ cos(k x ) cos(k y ), corresponding to the proposed signreversed s wave (s ± ) scenario of superconductivity [3,4]. The superconducting term in the model (23) is derived from the interaction term and the gap equation is We here assume V ij to take on a negative value, V , for second neighbours and be equal to zero otherwise. The field operators are related to the excitation operators,γ n , with energy E n , through the equations [23] d i↑,ν = n u n,ν (i)γ n − v * n,ν (i)γ † n (26) The Bogolubov-de Gennes equations for the amplitudes u and v, using the model (23), are: Equations (28)- (29) and (25) have been solved self-consistently, by iterations until convergence was achieved, in a 30×30 lattice, with open boundary conditions. The vortex core is located at the centre of a plaquette (see Figure 8). We take V = −4 and the chemical potential E F in the range 1.35-1.8, below. The gap function describing a vortex can be written in the form [33] where θ i is the azimuthal angle of site i. We define the s wave component of ∆ ij as and the d wave component of ∆ ij as

Results
For this choice of V ij both components are always present, even in the uniform (no vortex) case. The d wave component is, nevertheless, rather small, as can be seen from Figures 9 and 10, which show ∆ s (i) and ∆ d (i) in a region including the vortex core. The d wave component displays the more complex behaviour: ∆ d (i) has four minima, symmetric with respect to rotations through the angle π (not π/2) around the vortex centre. The s wave component, on the other hand, has a single minimum.   . At the latter location, where ∆ ij is almost uniform, the LDOS displays the typical behaviour: it is small at low energy and has two maxima at energies close to ±∆. Due to the band structure, the curve is not quite symmetrical around E = 0. The LDOS at the vortex centre is not typical: contrary to that obtained in the continuum model above, it has a fairly sharp maximum at finite energy above the Fermi level. The peak energy is smaller than but of the order of ∆ and is also dependent on the chemical potential. In addition to this peak, the LDOS also displays a V-shaped background with weak oscillations which are due to finite size effects.
The LDOS peak reveals the existence of localized excitations, at the vortex core, at finite energy and Figure 13 shows one such localized excitation. The low-lying excitations are spatially delocalized. Such behaviour can be qualitatively understood if we see the core excitations as Andreev bound states that form in the region of small ∆ that is in contact with the region of large ∆. The Andreev reflection problem has recently been investigated in the context of the s ± superconductivity [20]. It was found that Andreev bound states form at finite energy because of interference effects between waves in different bands of the superconductor. In the continuum formulation used in the previous section, the bands were treated separately and such interference effects were lost. We note that the familiar Caroli-de Gennes-Matricon expression ∆ 2 /E F for the bound states energy cannot be straightforwardly used here because it applies to a single electron parabolic band with the Fermi energy measured from the bottom of the band, assumed to be zero. The Fermi energy values, E F , above, are not measured from the band bottom. The band structure has half-metal character, with hole and electron bands and quantum interference effects between them.  In the STM experiment of Ref [19], the vortex positions were not correlated with the positions of impurities on the sample surface. We find that the presence of an impurity far from the vortex core only affects the LDOS in the impurity's vicinity. This is in agreement with results in other systems where the effect of impurities is rather local [27]. However, we may consider a situation where the vortices are pinned by impurities even though there seems to be bulk pinning in the experimental system studied in Ref. [19]. The effect of an impurity located at the vortex core is shown in Figures 14 and 15. While an impurity located outside the vortex has little effect on the core LDOS peak, an attractive impurity close to the core shifts the peak to negative energy, as shown in Figure 14, whereas the LDOS outside the core remains unchanged.
A repulsive impurity close to the vortex core has a stronger effect: the former LDOS peak is almost suppressed and shifted to higher energy, as shown in Figure 15, while the LDOS outside the core remains unchanged. The location of the low energy peak is therefore shifted depending on the type and strength of the impurity and on the band filling.

Summary and discussion
Both the continuum and the lattice formulations presented here lead to localized excitations in the vortex core. In the first case there is a zero energy peak in the LDOS at the vortex location resulting from the addition of two off-centre peaks, one from the electron band and one from the hole band. In the lattice formulation, interference effects from the two bands are taken into account, and only one gap function is introduced containing both bands. The low energy states are extended and localized modes at the vortex core appear at higher energy producing a LDOS peak deviated from the Fermi level. This may be understood from the recent prediction that Andreev bound states have finite energy in the s ± scenario. It also implies that the density of (low lying extended) states should be proportional to the density of vortices, hence proportional to the applied magnetic field, H. Therefore, if the quasiparticle mean-free-path is longer than intervortex spacing [34] then linear behavior of zero temperature heat transport, κ 0 /T (H) ∝ H, may be expected [35].
Our calculations shows that the effect of impurities at the vortex location has a significant effect on the peak in the LDOS but does not remove it.
The prediction of localized excitations in the vortex core is in conflict with the experimentally observed LDOS, with a STM, in BaFe 2 Co 0.2 As 2 , where no peak in ρ(E, r) was observed at the vortex centre [19]. Indeed, the measured differential conductance displays only the V-shaped background line, similar to that visible in figures 11 and 12 if one disregards the central peak and weak oscillations (due to finite size effects). The observed LDOS away from the vortex centre is similar to that in figures 11 and 12, showing a depression below the gap energy.
In other superconductors either conventional or unconventional there are low energy peaks in the LDOS at the vortex location. In the case of s-wave conventional superconductors there are peaks corresponding to the Caroli-de Gennes-Matricon states [25,28] and in d-wave unconventional superconductors there are two low-energy peaks that replace the coherence peaks at the vortex location and immediate vicinity [27]. Most superconductors have some sort of disorder leading in general to positional disorder of the vortices. This leads to an increase of the low energy LDOS and possible closing of the gap, in the s-wave superconductors, and to a finite density of states in the d-wave case [18,36]. The effect of a moderate concentration of impurities was considered in the d-wave case showing that the effect of the vortex disorder is dominant [37,38]. However, in the very dirty superconductor Nb 1−x Ta x Se 2 , where the quasiparticle mean path is smaller than the coherence length, the low energy peaks in the LDOS have been observed to give way to an almost flat conductance similar to the normal state one [39]. Similar results were obtained in the case of MgB 2 , even though in this system the mean path is expected to be large [40].
In the course of completion of this work, two preprints appeared where the LDOS in a vortex lattice was theoretically studied [41,42]. Our results for a single vortex are in agreement with those studies.