Local spin-density-wave order inside vortex cores in multiband superconductors
LLocal spin-density-wave order inside vortex cores in multiband superconductors
Vivek Mishra and Alexei E. Koshelev
Materials Science Division, Argonne National Laboratory, Lemont, IL 60439, USA (Dated: October 21, 2018)Coexistence of antiferromagnetic order with superconductivity in many families of newly dis-covered iron-based superconductors has renewed interest to this old problem. Due to competitionbetween the two types of order, one can expect appearance of the antiferromagnetism inside thecores of the vortices generated by the external magnetic field. The structure of a vortex in typeII superconductors holds significant importance from the theoretical and the application points ofview. Here we consider the internal vortex structure in a two-band s ± superconductor near a spin-density-wave instability. We treat the problem in a completely self-consistent manner within thequasiclassical Eilenberger formalism. We study the structure of the s ± superconducting order andmagnetic field-induced spin-density-wave order near an isolated vortex. We examine the effect ofthis spin-density-wave state inside the vortex cores on the local density of states. PACS numbers: 74.20.De, 74.25.Op, 74.25.Ha
I. INTRODUCTION
The emergence of superconductivity at the onset ofmagnetism is a hallmark of many families of unconven-tional superconductors. Recently discovered iron basedsuperconductors (FeSCs) provided a new addition to thislist. Parent compounds for many of the FeSCs havethe spin-density-wave (SDW) order, and superconduc-tivity (SC) appears upon doping or under pressure.
The doping-temperature phase diagrams describing lo-cation of these two phases vary from material to ma-terial. In some systems, for example in the 1111 fam-ily RFeAsO − x F x (R is a rare-earth element), the SDWphase abruptly disappears once the SC phase devel-ops. In other systems, such as 122 compounds basedon BaFe As , the SDW order coexists with SC withinsome range of parameters.To understand how the SDW and SC phases interactwith each other and what triggers superconductivity withsuch high transition temperatures, knowledge of the pre-cise structure of order parameters is essential. The localelectronic structure near defects can be extracted fromscanning tunneling spectroscopy (STM) measurementsand it provides vital information about the order parame-ter. A Defect could be either an impurity or a topologicalsingularity like a vortex induced by magnetic field. Herewe focus on structure of an isolated vortex.It is well known that the shape of the vortex andthe electronic structure close to it are very sensi-tive to the gap structure. Copper oxide based high-temperature superconductors have been subjected to ex-tensive research for various kinds of competing ordersinside the vortex cores. The signatures of such vortex-core orders have been reported in Bi Sr CaCu O δ , La − δ Sr δ CuO , YBa Cu O − δ , YBa Cu O , and Tl Ba CuO δ . From the standpoint of the the-ory, several different approaches have been adopted toexplain these experimental observations. Arovas et al. and Sachdev et al. studied antiferromagnetism in thevortex cores within phenomenological Ginzburg-Landau free-energy functional method. Ghosal et al. used mi-croscopic Bogoliubov-de Gennes (BdG) technique to in-vestigate this problem for the superconductors with ad-wave symmetric order parameter. The BdG techniquewas used heavily by many researchers to understand var-ious aspects of the competing orders inside the vortexcores. We also mention the work of Garkusha et al. where the Usadel equation formalism has been used toexplore the problem of antiferromagnetic vortex cores.The Usadel equations, however, are only applicable inthe dirty limit, when the electronic mean free path isshorter than the coherence length, hence only appropri-ate for s-wave superconductors.The vortex state in FeSCs has been studied extensivelyby the STM and several novel features near vortexcores have been revealed. In the first study of the vor-tex structure in the optimally-doped BaFe . Co . As byYin et al. no subgap states have been found. This ismost probably due to large quasiparticle scattering ratein this material. On the other hand, optimally-dopedBa . K . Fe As does show peak at the vortex center,which is shifted from the Fermi level to lower energy. This shift was attributed to the quantum effect, whichis realized in the materials with moderate values of theproduct of the Fermi momentum k F and the coherencelength ξ at temperatures lower than T c / ( k F ξ ). Alter-natively, such energy shift of the localized state whichbreaks the particle-hole symmetry can be caused by mag-netic field-induced order in the vortex core. This sce-nario is very likely when a superconductor is close to aSDW instability. This possibility, however, has not beenconsidered in Ref. 34. Similar downshift was found inLiFeAs by Hanaguri et al. , even though this materialdoes not have obvious proximity to magnetism. Song et al. studied the vortex state in FeSe and found en-hanced C symmetry breaking in the vortex core whichis probably related to orbital order in this material.These compelling features have motivated many the-oretical works. One class of theories has associatedthe particle-hole asymmetric finite-energy peaks to the a r X i v : . [ c ond - m a t . s up r- c on ] A ug normal-state band structure of the materials. In thiscase the mechanism of particle-hole symmetry breaking isdue to the quantum effect discussed in Ref. 39. Contraryto this proposal, several other authors have consideredorbital , nematic or SDW order . Hung et al. have included orbital ordering within a self-consistentBdG approach, and explained the enhanced C symme-try breaking observed in FeSe by Song et al. in Ref. 36.Similar results were reported by Jiang et al. and Hu etal. for the SDW order also using the BdG method.In this paper, we consider the emergence of the SDWorder in the vortex cores and its spectroscopic conse-quences. We use the quasiclassical Eilenberger approachto study the problem of the field-induced SDW order in-side an isolated vortex. Both the BdG and Eilenbergerapproaches have their own advantages and complementeach other. The BdG method is more microscopic. Onthe other hand, the Eilenberger approach relies on fewmost essential physical parameters. It is numericallyless expensive and allows to study more complex prob-lems. We compute distribution of the superconductingand SDW order parameters inside the core and typicallength scales for both order parameters. We also investi-gate influence of emerging SDW order on the density ofstates (DOS) near the vortex. This paper is organized inthe following manner. In the next Sec. II, we describethe details of the model and the method. In Sec. III wediscuss the results and conclude in Sec. IV. II. MODEL & METHODA. Quasiclassical equations for a two-bandsuperconductor with spin-density wave
We consider a simple minimal model with two cylindri-cal Fermi surfaces, which allows us to capture qualitativeunderstanding of the problem. For the dispersion of theholelike Fermi surface, we take ξ h ( k ) ≡ ξ ( k ) = µ h − k m h , (1)and for the electronlike Fermi surface we consider follow-ing dispersion, ξ e ( k ) ≡ ξ ( k − Q ) = ( k x − Q x ) m e (1 − (cid:15) ) + ( k y − Q y ) m e (1+ (cid:15) ) − µ e , (2)where ( Q x , Q y ) is the SDW ordering vector and µ h , µ e are the energy offset for the hole and the electron bandrespectively. Fig. 1 shows a schematic picture of the twoFermi surfaces. It is useful to write these dispersions as, ξ h ( k ) = − ξ, (3) ξ e ( k + Q ) = ξ + 2 δ, (4)where δ is the energy scale, which measures the deviationfrom perfect nesting. In general, δ is a function of theangle on the Fermi surface φ and goes to zero at the hot (cid:48)(cid:42) Q FIG. 1. (Color online) Schematic representation of the hole-like and electronlike Fermi surfaces (solid lines) centeredaround the Γ point and M point respectively. A shifted hole-like Fermi surface is shown with dashed line. Filled circles arethe hot spots, where the nesting is perfect. spots (shown in Fig 1). For the dispersions consideredhere, δ ( φ ) = δ iso + δ ani cos 2 φ, (5)with δ iso = 12 (cid:18) m h µ h m e (1 − (cid:15) ) − µ e (cid:19) ,δ ani = m h µ h m e (cid:15) − (cid:15) , and we treat δ iso and δ ani as tuning parameters.The model Hamiltonian is same as used by severalother groups , H = H kin + H sc + H sdw , (6) H kin = (cid:88) i, k ,α ξ i ( k ) c † i, k ,α c i, k ,α , (7) H sc = (cid:88) i, k ,α,β (cid:104) ∆ i ( ıσ y ) αβ c † i, k ,α c † i, − k ,β + h.c. (cid:105) , (8) H sdw = (cid:88) k ,α,β (cid:104) M ∗ ( σ z ) αβ c † , k ,α c , k ,β + h.c. (cid:105) , (9)where ∆ i are the SC order parameters for two bandswith i = 1 , M is the SDWorder parameter. We only consider singlet superconduc-tivity. For incommensurate SDW order M is a complexquantity. Here we consider only the commensurate SDWorder, which makes M a real quantity. We will brieflydiscuss the consequences of incommensurability in theSDW order. The indices α , β denote the spin states,and c † i, k ,α ( c i, k ,α ) is the fermionic creation (annihilation)operator for a fermion in the i th band with spin α . Weconsider s ± state for superconductivity with equal gapmagnitudes in two bands with a relative sign change.The self-consistency conditions read,∆ i = (cid:88) j, k ,α,β V sc ij ( − ıσ y ) αβ (cid:104) c j, − k ,α c j, k ,β (cid:105) , (10) M = (cid:88) k,α,β V sdw ( σ z ) αβ (cid:68) c † , k ,α c , k ,β (cid:69) . (11)Here V sc and V sdw are the pairing interactions for theSC and SDW phases respectively and assumed to be mo-mentum independent.In the extended particle-hole basis, Ψ † = (cid:16) c † , k , ↑ , c , − k , ↓ , c † , k , ↑ , c , − k , ↓ (cid:17) , the Hamiltonian reads, H = (cid:88) k Ψ † · ˆ H · Ψ , (12)ˆ H = ξ ∆ M ∗ − ξ MM ξ ∆ M ∆ ∗ − ξ . (13)The 4 × G = (cid:16) iω ˆ − ˆ H (cid:17) − , (14)where ω = 2 πT ( n + 1 /
2) is the fermionic Matsubara fre-quency and ˆ is 4 × . These are transportlike equationsfor the kinetic energy integrated Green’s functions,ˆ g = iπ (cid:90) dξ ˆ γ · ˆ G , (15)where ˆ γ is a 4 × , − , − , (cid:104)(cid:16) ω + eic v F · A (cid:17) ˆ γ, ˆ g (cid:105) + v F · ∇ ˆ g + i (cid:104) ( ˆ H δ + ˆ H sc + ˆ H sdw )ˆ γ, ˆ g (cid:105) = 0 , (16)where v F is the Fermi velocity and A is the vector po-tential. ˆ H sc , ˆ H sdw are the SC and the SDW componentsof the mean-field Hamiltonian in the basis spanned by Ψ.ˆ H δ is a 4 × δ ,-2 δ )containing information about nesting between the Fermisurfaces. Its contribution drops out in the absence ofthe SDW order. This equation agrees with one derivedby Moor et al. . In the pure superconducting limit Eq.(16) reduces to the well-known Eilenberger equation, v F · ∇ ˆ g + (cid:104)(cid:16) ω + eic v F · A (cid:17) ˆ γ, ˆ g (cid:105) + i (cid:104) ˆ H sc ˆ γ, ˆ g (cid:105) = 0 . (17)Equation (16) has to be supplemented with the normal-ization condition, ˆ g = ˆ g bulk . (18)In particular, ˆ g bulk = ˆ for a uniform superconductorwithout SDW order. The self-consistency conditions for FIG. 2. (Color online) Coordinate systems used to solve theEilenberger equations. The real space lab frame is shown withsolid lines, while the dashed lines represent the new coordi-nate system. Filled circle is the location of the vortex coreand chosen as the origin. Open circle is the point where thesolution is required. A point in the lab frame r = x p ˆ x + y p ˆ y maps to λ p ˆ v + µ p ˆ u in the new coordinate system. the order parameters can be expressed in terms of theEilenberger functions as,∆ = iπT (cid:88) ω (cid:104) V sc11 g + V sc12 g (cid:105) F.S. , (19)∆ = iπT (cid:88) ω (cid:104) V sc12 g + V sc22 g (cid:105) F.S. , (20) M = iπT (cid:88) ω V sdw (cid:104) g + g (cid:105) F.S. . (21)Here g ij are components of 4 × i, j = (1 , ,
4) correspond to the hole/electronband). (cid:104) g ij (cid:105) F.S. means angular average over the Fermisurface for the respective bands weighted with the den-sity of states, which is approximately the same for theboth bands for the Fermi surfaces considered here. Asfor the pairing-interaction matrix, to get the s ± super-conducting state we take into account only the interbandrepulsive interaction and neglect intraband terms, i.e., V sc11 = V sc22 = 0. For pure SC state, Eq. (17) equationscan be transformed to a set of Riccati equations, whichmakes the numerical solution much easier. However thistransformation is not useful for the problem consideredhere. In the following section, we discuss the strategy tosolve these equations for the present case. δ ani /2 π T s0 ( T s , T c ) / T c SDW SDW +SC SC δ iso =0.16(2 π T s0 ) FIG. 3. (Color online) The phase diagram in the (T, δ ani )plane for the two phases at zero magnetic field. Thickblack/red line indicates transition to the SDW/SC state. Thethin dashed line shows the T s in the absence of the SC corre-lations. Filled squares are considered as representative casesin this paper. B. Numerical solution
The Eilenberger equations are first-order partial differ-ential equations. A standard tool for solution of this typeof equations is the method of characteristics. The basicidea of this method is to introduce the new coordinatesystem, in which the partial differential equation reducesto an ordinary differential equation. Hence, it is usefulto introduce a coordinate system spanned by the two or-thogonal vectors, a unit vector along the direction of theFermi velocity ˆ v and a unit vector ˆ u orthogonal to ˆ v , seeFig. 2. The unit vectors spanning the new coordinatesystem read, ˆ v = cos θ ˆ x + sin θ ˆ y, (22)ˆ u = − sin θ ˆ x + cos θ ˆ y. (23)Here θ is the angle between the Fermi velocity and the x axis in the lab frame. A point in the lab frame r = ( x, y )maps to ( λ, µ ) in the ˆ v -ˆ u frame. A point ( x p , y p ) at whichsolution is desired transforms to ( λ p , µ p ) as, λ p = x p cos θ + y p sin θ, (24) µ p = − x p sin θ + y p cos θ, (25)where the parameter µ p has the meaning of an impactparameter. For a fixed trajectory, this impact parameteris uniquely determined by ( x p , y p ) and it does not changewith change of λ . The quasiclassical equations are solvedalong these classical trajectories (ˆ v ) in the real space.Along such trajectories quasiclassical equations reduce tosystem of ordinary differential equations, which are mucheasier to handle than solving a set of partial differentialequations. Far away from the bulk, the system is ho-mogeneous. The homogeneous values are used as initial values. Now on a given trajectory, there are two possi-bilities. One can integrate towards the defect (vortex inthis case) from the two extreme ends ( λ = ±∞ ) of thetrajectory. Due to the first-order nature of the equations,the numerical solution readily converges to exponentiallygrowing function. Of course, these exponentially growingsolutions are unphysical. However, it is possible to con-struct the physically bounded solution at any point usingthe exponentially growing solution using the explosionmethod. The explosion method exploits exponentiallygrowing solutions to obtain the physical solution. (See Appendix A for details) For each point in the realspace, one has to solve the Eilenberger equations for allthe trajectories and for each Matsubara frequency. Toobtain a physical solution, we solve the Eilenberger equa-tions from two opposite directions λ = ±∞ towards thepoint, where the solution is desired. As shown in the Ap-pendix, the two exploding solutions ˆ g ± diverging in the ±∞ limits, provide the physical solution ˆ g p asˆ g p = ˆ g − ˆ g + − ˆ g + ˆ g − ˆ g − ˆ g + + ˆ g + ˆ g − . (26)Once all the Eilenberger functions are computed for aninitial guess for the order parameters, an updated set oforder parameters is recalculated, and this process con-tinues till it converges to a solution. It is convenient tonormalize all the energy scales to T c and all the lengthsare measured in the unit of superconducting coherencelength ξ = v F / (2 πT c ). Here v F is the average Fermi ve-locity of the two bands. We consider weak ellipticity forthe electronlike Fermi surface, and the Fermi velocitiesof the two bands are roughly equal. All our results arepresented in these units. III. RESULTS & DISCUSSION
Coexistence of the SDW state and the superconduc-tivity is very sensitive to the underlying electronic struc-ture. For the two-band model we consider, the nestingfunction δ ( φ ) in Eq. (5) can be tuned to get a co-existingphase. Here we are interested in a situation wherethere is no long-range SDW order in the absence of themagnetic field. Fig. 3 shows the phase diagram as a func-tion of the anisotropic nesting parameter δ ani for a fixedvalue of δ iso = 0 . πT s ). The presence of the super-conductivity strongly modifies the SDW phase bound-ary (thin dashed line in Fig. 3). The region betweenthe original and SC-renormalized phase boundaries pro-vides a possibility of the SDW order in the regions, wherethe SC phase get suppressed locally. The phase diagramshown in Fig. 3 is only includes the commensurate SDWphase. Vorontsov et al. have shown that the incom-mensurate SDW phase may co-exist with the SC in alarger area of the phase diagram. We consider few rep-resentative cases with δ ani / πT s = 0 . , .
26 and 0 . T s is the SDW transition temperature for a sys-tem with perfect nesting and with the same interaction c ( M , ∆ ) / T c ˜ δ ani = 0 . δ ani = 0 . δ ani = 0 . FIG. 4. (Color online) The temperature dependence of thebulk superconducting order parameter and the SDW ordersat the center of the vortex for three values of the parame-ter δ ani / πT s = 0 . .
26 and 0 .
27. All the energy scales arenormalized to T c . The thin line shows the temperature de-pendence of the magnitude of the SC order in the bulk, whichis almost identical for two bands. This signs of SC order pa-rameters are opposite for two bands. The mean-field SDWtransition temperature in the absence of superconductivity isindicated with filled square for each case. T=0.3T c x/ ξ y / ξ −10 −5 0 5 10−10−50510 −10 0 10−10123 ˆ x ˆ y FIG. 5. (Color online) The spatial variation of the magni-tude of the SDW order at T = 0 . T c for δ ani = 0 . πT s ),where T s is the SDW transition temperature for a systemwith perfect nesting and with the same interaction strength.The magnitude of the SDW order is normalized to T c . Thehot spots are located near the y-axis. strength. We have set T s = 2 T c , where T c is the su-perconducting temperature. As we are mostly interestedin low-temperature behavior, we restrict our calculationsbelow T = 0 . T c . A. Order parameter profiles
Fig. 4 shows the temperature dependence of the bulksuperconducting gap and the SDW order in the core. Itis evident that the SDW order appears roughly belowthe temperature one would expect from the phase di-agram shown in Fig. 3. The temperature dependenceof the SDW order parameter deviates strongly from themean-field behavior ( ∝ (cid:112) − T /T c ). We see two differ-ent temperature regimes. At lower temperature the SDWorder grows strongly, but slightly below the mean-fieldSDW transition T s it develops a tail, which survives evenabove T s . It should be noted here that the phase diagramis based on the commensurate SDW phase for the normalstate electronic structure. However, the electronic struc-ture of the vortex core states is not the same as in thenormal state. The onset of the SDW order is mostly de-termined by the core bound states. The SDW transitiontemperature T s ∝ exp[ − /N V sdw f ( δ iso , δ ani )], wherethe function f ( δ iso , δ ani ) depends on nesting parameters, V sdw is the SDW interaction and N is the density ofstates. As in the vortex core the density of states ishigher than the normal-state value, the SDW onset tem-perature may exceed the mean-field value.We restrict ourselves to the commensurate case only,but for the incommensurate case the phase boundaryshifts towards slightly higher temperatures. As reportedby Vorontsov et al. in Ref. 46, an incommensurate or-der may exist in a larger portion of the phase diagram.We have also performed calculations, where we allow in-commensurability in the SDW order. With incommen-surability the SDW order parameter becomes complexand acquires a finite phase. We found that the phaseof the SDW order parameter is temperature dependentand varies very weakly in the real space. The incommen-surate order persists above the phase boundary shownin Fig. 3. Since we did not find anything qualitativelydifferent, we will focus on the commensurate case only.Furthermore, there is no qualitative difference betweenthe cases considered here, except for the temperature de-pendence. Therefore, we continue our discussion with δ ani / πT s = 0 .
25. For this δ ani and δ iso / πT s = 0 . φ ≈ . ◦ is close to the y di-rection which strongly influences anisotropic propertiesof the vortex. Figure 5 shows the magnetic field-inducedSDW order parameter in the real space at T = 0 . T c .Spatial coordinates have been normalized to the super-conducting coherence length ξ and the SDW order isnormalized to T c . An important feature is the oscilla-tions of the SDW order along the x direction which ismost clearly seen in the inset. At lower temperatureswhen the SC vortex is very small, the SDW state is lo-calized very close to the vortex core. As the tempera-ture grows, vortex becomes larger and the region withthe SDW order also increases, due to a larger region ofthe suppressed superconductivity. The size of the vortexis larger in the presence of the field-induced SDW order.This property can be seen in Fig. 6, in which we compare x/ ξ y / ξ T=0.3T c −10 −5 0 5 10−10−50510 00.20.40.60.81 x/ ξ y / ξ T=0.3T c −10 −5 0 5 10−10−50510 00.20.40.60.81 FIG. 6. (Color online) The left panel shows the gap magnitude normalized to its bulk value at T = 0 . T c in the presence ofthe SDW order and the right panel shows the gap structure at the same temperature with no SDW order for the hole band. a SC vortex with and without the magnetic field-inducedSDW order. The SDW order makes vortex larger andanisotropic. The intrinsic anisotropy of the underlyingband structure is weak. Hence the large anisotropy inthe real space is mainly due to the field-induced SDWorder. Strong enhancement of the anisotropy is reflectedin the characteristic length scales associated with the SCorder along the two principal directions. Fig. 7 shows thelength scales associated with the SDW ( ξ sdw x/y ) and the SCorder ( ξ sc x/y ) along the x and y axis. We define the super-conducting coherence length ξ sc as a distance from thevortex core, where the order parameters reaches half ofits bulk value. Similarly, the magnitude of the SDW or-der drops to half of its value at the core at a distance ξ sdw from the vortex core. As illustrated in Fig. 7, theSC length scales along x and y directions become differ-ent in the presence of the SDW order and this is mainlydue to anisotropy in the field-induced SDW order. Notethat the SDW correlations are stronger along the y di-rection which is closer to the nesting hot spots. Thiscauses stronger suppression of the SC order and reducesthe SC characteristic length along this direction. Next,we discuss the density of states near the vortex core. B. Density of states
All the previous calculations were done using the Mat-subara frequencies. For the DOS calculation it is nec-essary to go to the real frequencies. We used the orderparameter profiles calculated in the Matsubara represen-tation. The analytic continuation, iω → E + iη , is donewith an artificial broadening η = 0 . T c . This processalso requires solution of the Eilenberger equations forthe real frequencies using the same approach as in thecalculation of the order parameters in the previous sec-tion. The total DOS is the sum of the partial DOSs for c ξ / ξ ξ ξ xsc ξ ysc ξ xsdw ξ ysdw FIG. 7. (Color online) The temperature dependence of char-acteristic length scale for SDW order ( ξ sdw x/y ) and the charac-teristic length scale of the SC order ( ξ sc x/y ) on the hole bandfor δ ani / πT s = 0 .
25. Subscripts denotes x and y directionsin the real space. The same characteristic length scale for apure superconductor is plotted with a dotted dashed line forcomparison. each band and, in terms of the Eilenberger functions, itis given by N ( E ) = Re (cid:104) [ g ( E + iη ) + g ( E + iη )] (cid:105) F S . (27)Fig. 8 presents evolutions of the density of states alongthe two principal directions for two temperatures. Thefirst row of Fig. 8 shows DOS along the x axis, which isaway from the hot spots and the second row shows DOSalong the y axis which is closer to the hot spots. Thefield-induced SDW order enhances the violation of theC rotational symmetry the vortex center. For conven-tional superconductors, the DOS is always particle-hole ω /T c x / ξ sc T=0.3T c (a)−2 −1 0 1 2−10−50510 00.511.522.53 ω /T c x / ξ sc T=0.5T c (c)−2 −1 0 1 2−10−50510 123456 ω /T c y / ξ sc T=0.3T c (b)−2 −1 0 1 2−10−50510 00.511.522.53 ω /T c y / ξ sc T=0.5T c (d)−2 −1 0 1 2−10−50510 123456 FIG. 8. (Color online) DOS along the x direction (first row) and along y direction (second row) at T=0.3 T c in panel (a), (b)and at T=0.5 T c in panel (c) and (d) for δ ani / πT c =0.26. symmetric, unless the superconductor is in the quantumregime when k F ξ is not too large and T (cid:28) T c /k F ξ . Another key feature of the classical clean-limit DOS inthe vortex core is the sharp peak at zero energy corre-sponding to the localized state. The quantum effects shiftthis zero-bias peak to the finite energy but, do not breakthe rotational symmetry. The emergence of the SDWorder leads to particle-hole asymmetry in the DOS andalso strongly violates the C symmetry. The particle-holeasymmetry is bigger along x direction because of strongerdeviation from nesting in this direction. Another impor-tant feature which is visible in the Fig. 8 is suppressedspectral weight at the core indicating that the energy oflocalized state is shifted from zero to a finite value cor-responding to opening of a minigap in the core. This isshown more clearly in Fig. 9 in which we show the DOSplots at several representative points. This small gap inthe DOS near the vortex core is particle-hole asymmetric, which is a hallmark of energy gap due to the SDW order.This gap vanishes away from the vortex core indicatingpresence of a state with energy close to zero localizedoutside the core. As the temperature increases and theSDW order weakens, the apparent gap in the core dis-appears as shown in panel (c) and (d) of Fig. 8 at T =0.5 T c . The described features are the keys to distinguishbetween the quantum effect and the field-induced SDWorder. Figure 8 also shows the DOS for two differenttemperatures. The temperature dependence of the DOSis easy to understand. As the temperature increases, theSDW order weakens, which reduces the degree of C sym-metry breaking and the particle-hole asymmetry in theDOS. −2 −1 0 1 201234 ω /T c N ( ω , x / ξ sc ) (a)T=0.3T c x/ ξ sc =0.00x/ ξ sc =0.25x/ ξ sc =0.50x/ ξ sc =0.75x/ ξ sc =1.00 −2 −1 0 1 201234 ω /T c N ( ω , y / ξ sc ) (b)T=0.3T c y/ ξ sc =0.00y/ ξ sc =0.25y/ ξ sc =0.50y/ ξ sc =0.75y/ ξ sc =1.00 FIG. 9. (Color online) Panels (a) and (b) show DOS near the vortex core along x and y directions respectively at T=0.3 T c . IV. SUMMARY AND CONCLUSION
We study the structure of an isolated superconductingvortex near a SDW instability inside the superconduct-ing dome. We show that the SDW order develops insidethe vortex below a critical temperature determined bythe strength of the SDW instability. This leads to C symmetry breaking near the vortex core. If there is al-ready C breaking in underlying band structure, then itgets strongly enhanced due to the SDW order near thevortex core. The corresponding deformation of the vor-tex shape can be imaged by the STM technique. Wefind that the field-induced SDW order persists beyondthe superconducting vortex region. The tunneling DOScarries very strong signatures of this field-induced order.A small energy gap develops inside the core and gives riseto strong particle hole asymmetry, which is pronouncedalong the directions away from the hot spots. Our re-sults are in qualitative agreement with the STM data onon Ba . K . Fe As which may indicate the presence ofthe vortex-core SDW order in this material. Our findingsalso agree with BdG-based works by other groups. ACKNOWLEDGMENTS
This work was supported by the Center for EmergentSuperconductivity, an Energy Frontier Research Centerfunded by the US DOE, Office of Science, under AwardNo. DE-AC0298CH1088.
Appendix A: Explosion method
The general structure of the Eilenberger equations is, d ˆ gdλ = [ T , ˆ g ] , (A1)which is a first-order ordinary differential equation. Itis straightforward to show that if ˆ g is a solution of thisequation then ˆ g is also a solution, which impliesˆ g = ˆ a + a ˆ g, (A2)where ˆ a is a constant matrix and a is a complex num-ber. It can be further shown that the product of any twosolutions of the Eilenberger equations is also a solutionof these equations. This gives a very powerful relation,ˆ g = ˆ a = ˆ g bulk , (A3)which is very useful in obtaining the numerical solu-tions of these equations. There are multiple solutionsto this system of equations. In pure superconductingstate, there are three independent solutions. There aretwo divergent solutions along with a bounded physicalsolution. All the unphysical solutions decay to zero inthe bulk. Let’s consider two such unphysical solutions,ˆ g ± ∝ e ± νλ . It can be shown that commutator of thesetwo unphysical solutions gives the physically bounded so-lution, ˆ X ± = ˆ g − ˆ g + ± ˆ g + ˆ g − , (A4)ˆ˙ X ± = [ T , ˆ X ± ] . (A5)The bounded physical solution is,ˆ g p = ˆ c p ˆ X − , (A6)the constant ˆ c p can be determined using Eq. (A3) and(A6) and it reads, ˆ c p = 1ˆ X + . (A7)We use these unphysical solutions ˆ g + , − in the bulk andintegrate towards the vortex core starting from the bulk.Since these solutions grow exponentially, they can beeasily computed numerically with appropriate boundaryconditions. Far away from defects, we can ignore the spa-tial dependence of the order parameters. Since there is nolong-range SDW order, we have the standard Eilenbergerequations for the pure superconducting state in the bulkand in the basis we consider here, the Eilenberger Green’sfunction is a block diagonal matrix. The two bands arecoupled through the self-consistency condition. There-fore it is sufficient to illustrate the idea for one band, forwhich we write down the equations explicitly,˙ g = i (cid:0) ∆ ∗ f + ∆ f † (cid:1) , (A8)˙ f = − ωf − i ∆ g, (A9)˙ f † = 2 ωf † − i ∆ ∗ g. (A10)We first find two unphysical solutions, which can be de-termined easily. Let’s consider, g = c e ζλ , (A11) f = c e ζλ , (A12) f † = c e ζλ . (A13)where ζ = ± ν . Normalization condition requires,ˆ g ± = 0 . (A14) This ensures that all the unphysical solution decay tozero in the bulk. Which gives, c + c c = 0 , (A15) c = i √ c c . (A16)Using these conditions,( ζ + 2 ω ) c = 2∆ √ c c , (A17)( ζ − ω ) c = 2∆ ∗ √ c c . (A18)These two equations give the value of ζ = ± (cid:112) ω + | ∆ | /v F and c = − i ∆ ζv F +2 ω c , (A19) c = − i ∆ ∗ ζv F − ω c . (A20)Here we fix c = 1 and write the exploding solutions,ˆ g + = exp (cid:20) + 2 Qλv F (cid:21) (cid:20) − i ∆ p + i ∆ ∗ p − − (cid:21) , (A21)ˆ g − = exp (cid:20) − Qλv F (cid:21) (cid:20) − i ∆ p − i ∆ ∗ p + − (cid:21) , (A22) p ± = 1 ω ± Q , (A23) Q = (cid:112) ω + | ∆ | . (A24)Now we can write down the physical solution,ˆ g p = ˆ X − ˆ X + = 1 (cid:112) ω + | ∆ | (cid:20) ω − i ∆ i ∆ ∗ − ω (cid:21) (A25)Once we get the values of these two diverging solutions inbulk then we can use the bulk values to integrate towardsthe vortex core and find the physical solution using twounphysical solutions. D. C. Johnston, Advances in Physics , 803 (2010). P. C. Canfield and S. L. Bud’ko, Annual Review of Con-densed Matter Physics , 27 (2010). P. J. Hirschfeld, M. M. Korshunov, and I. I. Mazin, Re-ports on Progress in Physics , 124508 (2011). G. R. Stewart, Rev. Mod. Phys. , 1589 (2011). H.-H. Wen and S. Li, Annual Review of Condensed MatterPhysics , 121 (2011). A. Chubukov, Ann. Rev. Cond. Mat. Phys. , 57 (2012). Ø. Fischer, M. Kugler, I. Maggio-Aprile, C. Berthod, andC. Renner, Rev. Mod. Phys. , 353 (2007). J. E. Hoffman, E. W. Hudson, K. M. Lang, V. Madhavan,H. Eisaki, S. Uchida, and J. C. Davis, Science , 466(2002). B. Lake, G. Aeppli, K. N. Clausen, D. F. McMorrow,K. Lefmann, N. E. Hussey, N. Mangkorntong, M. Nohara,H. Takagi, T. E. Mason, and A. Schr¨oder, Science ,1759 (2001). B. Khaykovich, Y. S. Lee, R. W. Erwin, S.-H. Lee, S. Waki-moto, K. J. Thomas, M. A. Kastner, and R. J. Birgeneau,Phys. Rev. B , 014528 (2002). B. Lake, H. M. Rønnow, N. B. Christensen, G. Aep-pli, K. Lefmann, D. F. McMorrow, P. Vorderwisch,P. Smeibidl, N. Mangkorntong, T. Sasagawa, M. Nohara,H. Takagi, and T. E. Mason, Nature (London) , 299(2002). J. M. Tranquada, C. H. Lee, K. Yamada, Y. S. Lee, L. P.Regnault, and H. M. Rønnow, Phys. Rev. B , 174507(2004). B. Khaykovich, S. Wakimoto, R. J. Birgeneau, M. A. Kast-ner, Y. S. Lee, P. Smeibidl, P. Vorderwisch, and K. Ya-mada, Phys. Rev. B , 220508 (2005). J. Chang, A. P. Schnyder, R. Gilardi, H. M. Rønnow,S. Pailhes, N. B. Christensen, C. Niedermayer, D. F. Mc-Morrow, A. Hiess, A. Stunault, M. Enderle, B. Lake,O. Sobolev, N. Momono, M. Oda, M. Ido, C. Mudry, andJ. Mesot, Phys. Rev. Lett. , 077004 (2007). J. Chang, C. Niedermayer, R. Gilardi, N. B. Christensen,H. M. Rønnow, D. F. McMorrow, M. Ay, J. Stahn,O. Sobolev, A. Hiess, S. Pailhes, C. Baines, N. Momono,M. Oda, M. Ido, and J. Mesot, Phys. Rev. B , 104525(2008). V. F. Mitrovic, E. E. Sigmund, W. P. Halperin, A. P.Reyes, P. Kuhns, and W. G. Moulton, Phys. Rev. B ,220503 (2003). D. Haug, V. Hinkov, A. Suchaneck, D. S. Inosov, N. B.Christensen, C. Niedermayer, P. Bourges, Y. Sidis, J. T.Park, A. Ivanov, C. T. Lin, J. Mesot, and B. Keimer,Phys. Rev. Lett. , 017001 (2009). T. Wu, H. Mayaffre, S. Kr¨amer, M. Horvati´c, C. Berthier,W. N. Hardy, R. Liang, D. A. Bonn, and M.-H. Julien,Nature (London) , 191 (2011). J. Chang, E. Blackburn, A. T. Holmes, N. B. Christensen,J. Larsen, J. Mesot, R. Liang, D. A. Bonn, W. N. Hardy,A. Watenphul, M. V. Zimmermann, E. M. Forgan, andS. M. Hayden, Nature Physics , 871 (2012). D. Leboeuf, S. Kr¨amer, W. N. Hardy, R. Liang, D. A.Bonn, and C. Proust, Nature Physics , 79 (2012). T. Wu, H. Mayaffre, S. Kr¨amer, M. Horvati´c, C. Berthier,P. L. Kuhns, A. P. Reyes, R. Liang, W. N. Hardy, D. A.Bonn, and M.-H. Julien, Nature Communications , 2113(2013). K. Kakuyanagi, K.-i. Kumagai, and Y. Matsuda, Phys.Rev. B , 060503 (2002). K. Kakuyanagi, K. Kumagai, Y. Matsuda, andM. Hasegawa, Phys. Rev. Lett. , 197003 (2003). D. P. Arovas, A. J. Berlinsky, C. Kallin, and S.-C. Zhang,Phys. Rev. Lett. , 2871 (1997). S. Sachdev and E. Demler, Phys. Rev. B , 144504 (2004). A. Ghosal, C. Kallin, and A. J. Berlinsky, Phys. Rev. B , 214502 (2002). J.-X. Zhu and C. S. Ting, Physical Review Letters ,147002 (2001). H.-Y. Chen and C. S. Ting, Phys. Rev. B , 220510(2005). W. A. Atkinson and J. E. Sonier, Phys. Rev. B , 024514(2008). M. Schmid, B. Andersen, A. Kampf, and P. Hirschfeld,New Journal of Physics , 053043 (2010). V. V. Garkusha and V. N. Krivoruchko, Journal of LowTemperature Physics , 37 (2005). J. E. Hoffman, Rep. Prog. Phys. , 124513 (2011). Y. Yin, M. Zech, T. L. Williams, X. F. Wang, G. Wu, X. H.Chen, and J. E. Hoffman, Phys. Rev. Lett. , 097002(2009). L. Shan, Y.-L. Wang, B. Shen, B. Zeng, Y. Huang, A. Li,D. Wang, H. Yang, C. Ren, Q.-H. Wang, S. Pan, and H.-H. Wen, Nat. Phys. , 325 (2011). T. Hanaguri, K. Kitagawa, K. Matsubayashi, Y. Mazaki,Y. Uwatoko, and H. Takagi, Phys. Rev. B , 214505(2012). C.-L. Song, Y.-L. Wang, P. Cheng, Y.-P. Jiang, W. Li,T. Zhang, Z. Li, K. He, L. Wang, J.-F. Jia, H.-H. Hung,C. Wu, X. Ma, X. Chen, and Q.-K. Xue, Science ,1410 (2011). M. A. N. Ara´ujo, M. Cardoso, and P. D. Sacramento, NewJ. Phys. , 113008 (2009). D. Wang, J. Xu, Y.-Y. Xiang, and Q.-H. Wang, Phys.Rev. B , 184519 (2010). N. Hayashi, T. Isoshima, M. Ichioka, and K. Machida,Phys. Rev. Lett. , 2921 (1998). H.-H. Hung, C.-L. Song, X. Chen, X. Ma, Q.-k. Xue, andC. Wu, Phys. Rev. B , 104510 (2012). D. Chowdhury, E. Berg, and S. Sachdev, Phys. Rev. B ,205113 (2011). H.-M. Jiang, J.-X. Li, and Z. D. Wang, Phys. Rev. B ,134505 (2009). X. Hu, C. S. Ting, and J.-X. Zhu, Phys. Rev. B , 014523(2009). Y. Gao, H.-X. Huang, C. Chen, C. S. Ting, and W.-P. Su,Phys. Rev. Lett. , 027004 (2011). A. B. Vorontsov, M. G. Vavilov, and A. V. Chubukov,Phys. Rev. B , 060508 (2009). A. B. Vorontsov, M. G. Vavilov, and A. V. Chubukov,Phys. Rev. B , 174538 (2010). R. M. Fernandes and J. Schmalian, Phys. Rev. B ,014521 (2010). G. Eilenberger, Z. Phys. , 195 (1968). A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP ,1200 (1969). A. Moor, A. F. Volkov, and K. B. Efetov, Phys. Rev. B , 134524 (2011). N. Schopohl, arXiv:cond-mat/9804064 , (1998). E. V. Thuneberg, J. Kurkij¨arvi, and D. Rainer, Phys. Rev.B , 3913 (1984). W. Zhang, J. Kurkij¨arvi, and E. V. Thuneberg, Phys. Rev.B , 1987 (1987). U. Klein, J. of Low Temp. Phys.69