Orbital-resolved vortex core states in FeSe Superconductors: calculation based on a three-orbital model

2.8 MB·

by  Q. E. Wang

Orbital-resolved vortex core states in FeSe Superconductors: calculation based on a
three-orbital model
Q. E. Wang1 and F. C. Zhang1,2,3,∗
1Department of Physics and Center of Theoretical and Computational Physics,
The University of Hong Kong, Hong Kong, China
2Department of Physics, Zhejiang University, Hangzhou, China
3Collaborative Innovation Center of Advanced Microstructures, Nanjing 210093, China
(Dated: June17, 2015)
5
1
Westudyelectronic structureofvortexcorestatesofFeSesuperconductorsbased onat2g three-
0
orbital model by solving the Bogoliubov-de Gennes(BdG) equation self-consistently. The orbital-
2
resolved vortex core states of different pairing symmetries manifest themselves as distinguishable

n structuresduetodifferentquasi-particlewavefunctions. Theobtainedvorticesareclassifiedinterms
u of the invariant subgroups of the symmetry group of the mean-field Hamiltonian in the presence
J of magnetic field. Isotropic s and anisotropic s wave vortices have G5 symmetry for each orbital,
6 whereas dx2−y2 wave vortices show G∗6 symmetry for dxz/yz orbitals and G∗5 symmetry for dxy
1 orbital. In the case of dx2−y2 wave vortices, hybridized-pairing between dxz and dyz orbitals gives
rise to a relative phase difference in terms of gauge transformed pairing order parameters between

] d and d orbitals, which is essentially caused by a transformation of co-representation of
xz/yz xy
on qGu∗5alaintadtiGve∗6lysusibmgirloaurpp.atTtherencwalictuhlaetxepderliomcaenltdreenssuitlyts.ofTshteatpehs(aLseDdOiffSe)roefncdex2o−fyπ2 bweatvweeevnordtices shaonwd
c 4 xz/yz
- dxy orbital-resolved dx2−y2 wave vortices can be verified byfurther experiment observation.
r
p
PACSnumbers: 74.70.Xa,74.25.Wx,74.20.-z
u
s
.
t I. INTRODUCTION outholepocketisa greatchallengeto the weak-coupling
a
m theory where the superconductivity is proposed to be
drivenby the nesting of the electronandhole Fermi sur-
- Quantized vortices, as stable topological defects, ob-
d faces. The absence of the hole pocket in FeSe makes the
served in a variety of quantum systems such as super-
n argument difficult.
conductor and superfluid, are characterized by their na-
o
c ture of soliton solutions of dynamical systems1. Elec- Vortex structures in iron-based superconductors have
[ tronic structures of vortices in cuprate superconduc- beenstudied bya number ofauthors25–29. These studies
tors exhibit charging effects2–4 and are anisotropic due are mainly based on band structures having hole pocket
4
v to dx2−y2 wave pairing symmetry5. Earlier theoretical at Γ. We expect the vortex structure be affected by the
9 works have investigated the vortex line states based on Fermi surface topology. In this work we use a three-
5 microscopic models6–9. In iron-based superconductors, orbitalmodeltostudyvortexstructureofFeSeSCstate.
1 band structure and multi-orbital pairings play an im- The three-orbital microscopic model reproduces qualita-
7 portant role and the vortex structures may be richer tively the correctFermi surface with only electron pock-
1. due to multi-orbital dependency. Vortex core states ets. WesolvetheBdGequationself-consistentlytostudy
0 of two-fold rotational symmetry, which is proposed to orbital-resolved vortex core states for various SC pair-
4 be attributed to the orbital-dependent reconstruction in ing symmetries: isotropic s(on-site pairing), anisotropic
v:1 SFoeSneg esutpaelr.cofnrodmuctsocrasn1n0,inghatvuennbeeleinngrempiocrrtoesdcopbyy(SCT.ML). sn(enigexhtbonreasriteestpnaeiriginhgb)owrasivtees.paWireincgo)m,apnadredrxe2s−uyl2t(snoefarceaslt-
i experiment11. culations with that observed from the recent STM ex-
X
In most of iron-based superconductors, the Fermi sur- periment on FeSe vortex and suggest that the pairing
r
a faces consist of both electron pockets aroundM point at symmetry to be dx2−y2 wave. We predict that there is a
thecornersandholepocketsaroundΓpointofthefolded relative phase difference about π between pairing order
4
Brilliouin zone (BZ). An s wave superconducting (SC) parametersdefinedond andd orbitalsinthecase
± xz/yz xy
pairing symmetry has been proposed, where the pair- ofdx2−y2 orbital-resolvedvortices,whilesuchaphasedif-
ing order parameters at the electron and hole pockets ferenceistrivialinthecaseofisotropicsandanisotropic
have opposite signs13–22. FeSe superconductor is inter- s wave vortices. The paper is organized as follows. An
esting for its unique electronic structure in which only introduction of magnetic translation group and classi-
electron pockets are found and the hole pockets are well fication of vortex solutions are given in Section II. In
below the Fermi level, as angle-resolved photoemission section III we present the three-orbital model and the
spectroscopy(ARPES) shows12. When the hole pocket self-consistent BdG approach. In Section IV, we discuss
at Γ vanishes, we don’t have s wave pairing state any- properties of the vortex core states for different pairing
±
moreandthe s anddx2−y2 wavepairingstatesshouldbe symmetries and compare our results with experimental
considered23,24. The iron-based superconductivity with- observations. Finally, a summary is given in Section V.
2
II. MAGNETIC TRANSLATIONAL paid to the difference of vortex states defined between
SYMMETRY AND WINDING STRUCTURES OF A1g(isotropic s and anisotropic s wave) and B1g(dx2−y2
SINGLE VORTEX
wave) irreducible unitary representations of D group.
4
TheHamiltonianoftheSCsysteminthepresenceofa
homogeneousmagneticfieldalongzˆdirectionisobtained
From a theoretical point of view, vortex lattice in
from its zero-field form by modifying the hopping and
mixed states of type II superconductors is ground state
pairingterms with Peierlsphase35, respectively,whichis
offermionic systemwhichis characterizedbyinteraction
of the following form
betweenahomogenousmagneticfieldwithC symmetry
∞
andCooperpairswithadefiniteSCpairingsymmetry30.
H =H +H
0 pair
In iron-based superconductors, situation becomes com-
plicated because of orbital degrees of freedom. Conse- H0 = [t˜σσ(iα,jβ)−µδijδαβ]a†iασajβσ
quently, the crystal symmetry, band structure, and SC i,j,α,β,σ (1)
X
pairing symmetry, determine the electronic structure of H = [∆˜ (iα,jβ)a† a† +h.c.]
vortices. Among these constraints of symmetry, the vor- pair ↑↓ iα↑ jβ↓
i,j,α,β
tex structures are mainly dominated by magnetic trans- X
lation invariance, whose generator are crystal momen- in which
tum and vector potential of magnetic field31. However,
ie i
such conventional magnetic translation group defines a t˜ (iα,jβ)=t (iα,jβ)exp[ A~(~r)·d~r]
magnetic unit cell containing two vortices. It is not the σσ σσ ~cZj (2)
symmetrygroupofAbrikosovlatticeinwhichonlysingle ∆˜ (iα,jβ)=∆ (iα,jβ)exp[iφ(i,j)]
↑↓ ↑↓
vortexisstabilizedwithinonemagneticunitcell. Break-
through of this difficulty was presented by M. Ozaki et where a† (a ) denotes the creation(annihilation) op-
al.32. In their work the magnetic translation group de- erator oifασelecitαrσons with spin σ =↑,↓ and orbital α at
scribing single vortex was discovered to be a subgroup site i. t (iα,jβ) are hopping integrals and µ is the
σσ
of direct product of conventional magnetic translation chemical potential. We assume that the screening mag-
group and gauge transformation group U(1). Therefore, netic field inside the superconductor can be neglected
stable vortex structure can be solved numerically in one except for that the magnetic field is close to the up-
magnetic unit cell taking advantages of nontrivial wind- per critical field. The SC pairing mechanism has been
ing boundaryconditionsderivedfrompropertiesofmag- proposed to be of magnetic origin. In this paper, how-
netic translation group33. ever, we shall focus on the vortex core state and start
Instead of doing calculations of two vortices in one from an extended attractive Hubbard model for sim-
magnetic unit cell, we follow the method given by M. plicity. The SC order parameter stemming from the
Ozaki et al.32,33, in which only single vortex structures mean-field decoupling of the paired scattering term is
are calculated in one magnetic unit cell, so that the cal- expressedas ∆↑↓(iα,jβ)=V↑↓(iα,jβ)hajβ↓aiα↑i for sin-
culated results can be classified by irreducible represen- glet pairing channel. The Peierls phase35 in hopping
tations of magnetic translation group. The numerical terms comes from the fact that the Lagrangian of elec-
calculations in previous works,as mentioned above25–27, troninamagneticfieldcontainsadynamicalterm e~v·A~,
c
are mostly carried out for two vortices in one magnetic which gives rise to the phase accumulationin the propa-
unitcell. Thesevortexstates,however,cannotbeidenti- gatorofelectrondescribingthe hopping processbetween
fiedbyinvariantsubgroupsofmagnetictranslationgroup two lattice sites. The modification of pairing order pa-
because they belong to the irreducible representationsof rametersaccounts for eliminating the mixing of different
conventional magnetic translation group. Furthermore, pairing states under the action of magnetic translation
two vortices in one magnetic unit cell are not indepen- group. The mathematical interpretation of doing this is
dent because the induction of interaction between them. essentiallysearchingforgaugetransformedorderparam-
It is well-known that the topological defects in uncon- eters, which span a representation of magnetic transla-
ventional superconductors and superfluids with certain tion group32,33. The gauge transformation, carried out
symmetrybreakingbehavedistinguishablyfromthecon- by phase φ(i,j), has different definition with respect to
ventional singular(hard core) vortices34. For instance, a anisotropic s and dx2−y2 wave pairing states, whereas in
vortex in 3He has a finite amplitude of order parame- thecaseofisotropicswavepairingitistrivial. Thegauge
ters in the soft core region whose size is larger than the transformedorderparameterfordx2−y2 wavepairinghas
coherent length, whereas the winding structure is non- beenderivedbymeans ofgrouptheoreticalanalysis32,33.
trivial. Therefore in our numerical calculation, we con- The magnetic translation operator takes following form
centrate on winding structures of vortices for each or- in symmetric gauge A~ = −1~r×B~, when it acts on cre-
2
bitals, although the vortices in iron-based superconduc- ation operators31,32
torsaremostlyofhardcorefeature,andclassifythevor-
tex structures of isotropic s, anisotropic s, and dx2−y2 L(R~λ)a†iασ =eiπ2(Nvλxλy)T(R~λ)a†iασ
wofamveapganiertiincgtsryamnsmlaettiorinesgirnoutper.mSspoefciianlvaarttiaennttisounbwgriollubpes =eiπ2Nv[λxλy+N1(λxiy−λyix)]a†i+λ,ασ (3)
3
and the resultant transformation of order parameter is Here we follow method given by M. Ozaki et al.32,33
to derive the gauge transformed order parameters for
ha a i (4) anisotropic s ∼ cos(k ) · cos(k ) wave pairing symme-
j+λ,β↓ i+λ,α↑ x y
=eiπNv[λxλy+21Nλx(iy+jy)−21Nλy(ix+jx)]hajβ↓aiα↑i torny.pTaihriengrebsuolntsdoafloancgtioxˆn+ofyˆcdoinrjeucgtiaotne raorteationsubgroup
where R~λ = λxNxˆ+λxNyˆ is the basis vector of mag- C4z(ix,iy)haiα↓ai+xˆ+yˆ,β↑i
netic unit cell containing N lattice sites and Nv is the =e−2iKiyha a i
iα↓ i−xˆ+yˆ,β↑
number of vortices within one magnetic unit cell. We
C (i ,i )ha a i
have restricted ourselves to the cases of square vortex 2z x y iα↓ i+xˆ+yˆ,β↑
(7)
lattice with lattice constant set to unity. Eq. (3) de- =e2iK(ix−iy)ha a i
iα↓ i−xˆ−yˆ,β↑
fines actions of magnetic translation group {L(R~λ)} on C43z(ix,iy)haiα↓ai+xˆ+yˆ,β↑i
field operators, and all of the operations form a group
=e2iKixha a i
in representation space spanned by gauge transformed iα↓ i+xˆ−yˆ,β↑
order parameters, provided that certain group condi-
then a symmetric phase rearrangement can be made by
tion is satisfied. Note that the gauge transformation, multiplying a PeierlsphaseeiK(iy−ix) to regainthe mag-
as an internal symmetry transformation, takes its com-
netic translational symmetry as following
plexconjugateformwhenactsonannihilationoperators.
Different from the situation for conventional magnetic
translation group31 {T(R~λ)}: Nv = 2, the group condi- ∆˜a↑↓nis. s(iα,jβ)
tion of magnetic translation group, which is the symme- V (jβ,iα)
= ↑↓ ha a i[eiK(iy−ix)δ
try group of Abrikosov lattice, is that only single mag- iα↓ jβ↑ i+xˆ+yˆ,j
4 (8)
netic flux ϕ = hc is contained in one magnetic unit
cell32,33, i.e.0, N 2=e 1. It has been pointed out that +e−iK(ix+iy)δi−xˆ+yˆ,j
v
dx2−y2 ∼ cos(kx)−cos(ky) wave order parameters will +e−iK(iy−ix)δi−xˆ−yˆ,j +eiK(ix+iy)δi+xˆ−yˆ,j]
mixwithextendeds∗ ∼cos(k )+cos(k ), p ∼isin(k ),
x y x x The magnetic translationproperty of gaugetransformed
and p ∼ isin(k ) wave order parameters under opera-
y y order parameters for anisotropic s wave pairing state,
tion of magnetic translation group32,33. Such a mixing
originates from the fact the symmetry group of normal which is consistent with dx2−y2 wave, is
state Hamiltonian contains a local gauge transformation ∆˜anis. s(i+λ,α,j+λ,β)
↑↓
generatedbythevectorpotentialofamagneticfield. The (9)
re-definedSCgaugetransformedorderparameterstrans- =eiπNv[λxλy+N1(λxiy−λyix)]∆˜a↑↓nis. s(iα,jβ)
forming according to invariant subgroups of D group
4 where j is always related to i as next nearest neighbor
without any gauge component, as order parameters do
site pairing. Compare this expression with Eq. (4),
in the absence of magnetic field, are obtained by gener-
it is obvious that the gauge transformed order param-
ating all of them with the action of a conjugate rotation
eters(referring to order parameters thereafter) now form
subgroup {Ck (i ,j ),k = 1,2,3,4} on one of the pair-
4z x y a basis of representation of magnetic translation group
ing bonds of every local order parameters accompanied
andthemixingbetweenanisotropicsandd wavepair-
by a Peierls phase factor35. The generator of conjugate xy
ingstatesunderactionofmagnetictranslationgrouphas
rotation subgroup is defined as
been eliminated.
TheSCgroundstates,intheabsenceofmagneticfield,
C (i ,j )=T(i ,j )C T−1(i ,j ) (5)
4z x y x y 4z x y canbe classifiedbyfinding allthe invariantsubgroupsof
the symmetrygroupD ⊗U(1),whichhaveaone-to-one
whereC is4-foldrotationaroundtheoriginofthecoor- 4
4z
correspondencetotheirreducibleunitaryrepresentations
dinatesystem. Thereforethemixingoforderparameters
ofthesymmetrygroupofnormalstateHamiltonian36,37.
under magnetic translation is eliminated by re-defining
In the case of D point group symmetry, such a classifi-
rotations of all local order parameters at different sites 4
cation is obtained by the fact that D has three invari-
backtoorigin. Thedx2−y2 wavegaugetransformedorder antsubgroupsofindex 2,andthe two4dimensionalcyclic
parameter is consequently re-defined as
group, as a subgroup of U(1), compensate the phase
∆˜dx2−y2(iα,jβ) (6) change of order parameters by eiπ when the elements of
↑↓ coset representative acts on them. In the same manner,
V (jβ,iα) the ground state of a vortex structure can also be clas-
= ↑↓ ha a i(e±iKiyδ −e∓iKixδ )
2 iα↓ jβ↑ i±xˆ,j i±yˆ,j sifiedby finding allthe invariantsubgroupsof symmetry
groupof the Hamiltonian in a magnetic field32, and con-
where K = πNv and xˆ(yˆ) denote the unit vectors of sequently the winding structure of the vortexcore states
2N2
two-dimensionallattice. Notethatforsingletpairingthe havesymmetryconstraintsofdifferentclasses. Thetopo-
order parameters are symmetric under exchange of site- logicalcharacteristicsofvortexstatesarelocationofpin-
orbital quantum number. ning center, phase distribution of order parameters, and
4
TABLE I. Winding number of order parameters for differ-
ent pairing states. G ,i=1,2,3,4,5,6 are six maximal little
i
groups. G∗5,6 differs from G5,6 by taking the complex conju-
gate of gauge transformation. The order parameters trans-
form according to basis functions of D4 group as s wave:
∼ const., anisotropic s wave: ∼ cos(k )cos(k ), extended s∗
x y
wave: ∼cos(kx)+cos(ky),dx2−y2 wave: ∼cos(kx)−cos(ky),
and d wave: ∼sin(k )sin(k ), respectively. The index l is
xy x y
defined in Eq. (10).
(a) (b)
l W(s, anis. s, s∗) W(dx2−y2, dxy)
G1 ∼G2 0 0 2 or −2 FIG. 1. (color online) Schematic pictures showing the phase
G5 ∼G6 1 1 3 or −1 difference between G∗5 (a) and G∗6 (b) winding structures in
G∗5 ∼G∗6 −1 −1 1 or −3 thevicinityofthevortexcorecenter32. Thepurplecirclesde-
G3 ∼G4 2 2 4 or 0 pictlatticesitesonwhichtheSCorderparametersaredefined
and the arrows show the phase distribution of ∆↑↓(iα,jβ).
The blue arrows in (b) denote phase difference of −3π and
4
red arrows π4 between G∗6 and G∗5 winding structures.
windingnumber. Itturnsoutthatthewindingnumberof
vorticesofdifferentsymmetryproperties,havingastruc-
tural vanishing region, can be calculated from the sym-
the vortex core states based on an effective three-orbital
metryconstraintsofcorrespondingmaximallittlegroups. model39. Taking advantage of the 4-fold rotational sym-
In work of M. Ozaki et al.32,33, winding numbers W of
metry,theBlo¨chHamiltoniancanbewrittenasfollowing
s∗ and dx2−y2 wave vortices have been calculated. Here
wecalculateW foranisotropicsandd wavestatesand H = ψ†(k)M(k)ψ(k)
xy 0
list all the results in Table I, in which
k
X
M(k)=K +K eikx +C K C3 eiky
Gl =(e+tC )C˜l∧L 0 1 4z 1 4z
2x
(10) +C K C e−ikx +C3 K C e−iky
C˜l ={e−π2lkCk ,k =1,2,3,4} 2z 1 2z 4z 1 4z
4z +K ei(kx+ky)+C K C3 ei(−kx+ky)
2 4z 2 4z
NotethatC˜l alwaysactsonpairedfieldoperatorsrather +C K C ei(−kx−ky)+C3 K C ei(kx−ky)
2z 2 2z 4z 2 4z
than single particle operator. The derivation is based on
(12)
the fact that the generator of C˜l, as a symmetry trans-
formation of order parameters, leaves them invariant33. where ψ†(k) = [a† (k),a† (k),a† (k)] and the 4-fold ro-
xz yz xy
The winding structures of G∗ and G∗ vortices, which tationiscarriedoutbyoneofthegeneratorsofD group
5 6 4
have been obtained from our numerical calculations, are
0 −1 0
shown in Fig. 1, respectively, where they differ by a co-
C = 1 0 0 (13)
representation transformation as 4z  
0 0 1
−1  
3π 3π The irreducible hopping subsets40(in unit: eV) corre-
G∗ = G∗
6 4 5 4 sponding to on-site atomic energies, hopping along xˆ,
!
(11)
c c and xˆ+yˆdirections are
−1
π π
G∗6 = 4G∗5 4 K0 =diag(−µ,−µ,0.4−µ)
(cid:18) (cid:19)
b b 0.05 0.00 −0.20
The gauge transformation of field operator is defined as K = 0.00 0.01 0.00
1
φ·aiασ =e−iφ2aiασ32. Note that the global gauge trans- 0.20 0.00 0.20  (14)
formationof−3π or π arebothallowedbygrouptheory.  
4 4 0.02 0.01 0.10
Bbut it turns out from our numerical calculationthat the
K = 0.01 0.02 0.10
phase difference of π is more energetically favorable. 2  
4 −0.10 −0.10 0.20
 
For simplicity, the spin indices have been dropped. In-
III. METHODOLOGY AND BAND MODEL
stead of going along the boundary of the irreducible BZ,
an alternative path has been used to show the band
It has been reported that the electronic structure of structure with dominating orbital weights in Fig. 2 (a).
iron-based superconductors in the vicinity of the Fermi The projected density of states(PDOS) reveals strongly-
levelisdominatedbyd ,d ,andd orbitalsfromfirst- hybridized bands which are composedof d and d or-
xz yz xy xz yz
principlecalculation38,thereforeitisfeasibletocalculate bitals along the off-diagonalline of the extend BZ below
5
PDOS (a rb.units) where h˜σσ(iα,jβ)=t˜σσ(iα,jβ)−µδijδαβ and the order
2.0
parameters defined on different orbitals are
1.5
V (iα,jβ) ǫ
1.0 ∆˜↑↓(iα,jβ)=− ↑↓ 2 uniα↑↑vjnβ∗↓↑tanh(2kn↑T)
B
V) 0.5 ǫn↑X>0,<0 (17)
e
(
E 0.0
Eq. (3) and (15) give a nontrivial winding boundary
-0.5 condition to quasi-particle amplitudes as
-1.0
-1(.5 ) ( ) ( ) ( ) ( ) ( ) (cid:20) uvinni++λλ,,αα↓↑↑↑ (cid:21)="ee−iiπ2π2NNvv[λ[λxxλλyy++N1N1(λ(λxxiyiy−−λλyyixix)])u]vniiαnα↑↓↑↑ #
(a)
(18)
π
The order parameters are calculated by BdG equa-
tionself-consistentlywiththeaboveboundarycondition,
which is assigned to the matrix element h˜ (iα,jβ) and
σσ
∆˜ (iα,jβ) for N = 1. The self-consistent calculation
↑↓ v
starts with arbitrarily distributed order parameters and
ky 0 the iteration is performed with a convergence criterion
that the order parameters have relative difference less
that 10−3 between two consecutive steps. The particle
densityarecalculatedviaquasi-particlewavefunctionsas
−π 1 ǫ
−π 0 π hn i= |un |2[1−tanh( n↑ )]
(b)kx iα↑ 2ǫn↑X>,<0 iα↑↑ 2kBT (19)
1 ǫ
hn i= |vn |2[1+tanh( n↑ )]
iα↓ 2 iα↓↑ 2k T
B
FIG.2. (coloronline)Orbital-resolvedbandstructure,PDOS ǫn↑X>,<0
(a) and Fermi Surface (b). The red (d ), green (d ), and
xz yz
Theenergyspectrumofthequasi-particle,i.e.,theLDOS
blue (d ) curves represent wight-dominating orbitals. The
xy
Fermi level has been set to zero. at site i for orbital α is calculated via
1
ρ (ǫ)= |un |2δ[ǫ−ǫ (~k)] (20)
iα M M iα↑↑ n↑
x y
theFermilevel. TheFermisurface(Fig. 2(b)),obtained ~k∈XFBZǫn↑X>,<0
with a chemical potential µ=0.312 eV corresponding to +|vinα↓↑|2δ[ǫ+ǫn↑(~k)]
a filling factor n=4.23, has four electron pockets which
where the supercell method has been used43 for M =
do not have any SC gap node in cases of anisotropic s x
M = 10. The Lorentzian smearing method is used to
and dx2−y2 wave pairing sates. The absence of electron y
visualize the LDOS with a broadening width σ = 0.001.
orholepocketsatΓpointisconsistentwithexperimental
observation12. All the self-consistent calculations are performed on a
28×28 lattice at temperature T=0.1K.
The Hamiltonian in Eq. (1) can be diagonalized by
conducting the Bogoliubov-Valatin transformation41,42 Calculationofmagneticexchangecouplingsshowsthat
the leading pairing instability comes from the intra-
containing t orbital degrees of freedom as
2g
orbital pairing contribution, whereas the inter-orbital
a = un γ +σ¯vn∗ γ† (15) componentsarefoundtobesignificantlysmall39. Conse-
iασ iασσ nσ iασσ¯ nσ¯ quently,onlyintra-orbitalpairingpotentialisconsidered
ǫnX↑>0 in our numerical calculation. The SC gap function for a
multi-orbital superconductor is generally defined in mo-
where the quasiparticle creation operator γ† is the lad-
nσ mentum space as
der operator of the eigen-spectrum of the Hamiltonian
which satisfies [H,γn†σ]− = ǫnσγn†σ. The diagonal condi- ∆i (~k)=gi(~k)Γ (iσ ) (21)
tion of the Hamiltonian is the BdG equation αβ αβ 2
h˜ (iα,jβ) ∆˜ (iα,jβ) un un wheregi(~k)isbasisoftheirreducibleunitaryrepresenta-
∆˜↑∗↑(iα,jβ) −h˜↑∗↓ (iα,jβ) vnjβ↑↑ =ǫn↑ vniα↑↑ tionsofD4 pointgroup,iσ2 definesatensorstateforsin-
Xj,β (cid:20) ↑↓ ↓↓ (cid:21)(cid:20) jβ↓↑ (cid:21) (cid:20) iα↓↑ (cid:21) gletpairing,and Γαβ is the orbitalbasisfor D4 transfor-
(16) mation. Thetransformationpropertiesofbandstructure
6
π 1.1 28
28 0.056
1.01
0.048 21
21 0.04
0.92
ky 0 14 0.032 14
0.024
0.83
0.016 7
7
0.74 0.008
7 14 21 28 0 7 14 21 28
−π 0.65 (a) (b)
−π 0 π
kx 28
28 0.056
FIG.3. (coloronline)ColormappingofFermivelocity~v (in
F 0.048 21
unit eV·m).
21 0.04
0.032 14
14
0.024
determine all the symmetry transformation of SC order
parameters39,45. Another reason that the inter-orbital 7 0.016 7
pairing has been omitted in our calculation is that only 0.008
if Γαβ transform according to A1g representation, then 7 14 21 28 0 7 14 21 28
symmetry of pairing state canbe exclusively determined (c) (d)
by its spatial component gi(~k), such that the calculated
28
vortex sate has a classificationof Table I. For isotropics
28 0.032
wave pairing,
21
V↑↓(iα,jα)=−g0δij (22) 21 0.024
for anisotropic s wave pairing, 14
14 0.016
g
1
V (iα,jα)=− (δ +δ (23)
↑↓ i+xˆ+yˆ,j i−xˆ+yˆ,j
4 7 0.008 7
+δ +δ ) (24)
i−xˆ−yˆ,j i+xˆ−yˆ,j
and for dx2−y2 wave pairing, 7 14(e) 21 28 0 7 ( 1f4) 21 28
g
2
V (iα,jα)=− (δ +δ +δ +δ )
↑↓ i+xˆ,j i+yˆ,j i−xˆ,j i−yˆ,j
2 FIG. 4. (color online) Amplitudes(color mapping) and phase
(25)
distribution of order parameters for isotropic s wave pairing
state for d orbital (a) and (b), d orbital (c) and (d), and
whereg arepairingamplitudes foreachpairingsym- xz yz
0,1,2 d orbital(e)and(f),respectively. Thephasedistributionof
metry. Fig. 3showstheFermivelovity~~v (~k)=∇ ǫ (~k) xy
n ~k n orderparametershavebeenmappedtoavectorfield. Length
whichisusedtodeterminethepairingpotential. Inorder of arrows represent theamplitude of order parameters.
to mimic the intermediate coupling cases for FeSe10 and
A Fe Se (A=K, Rb, orCs)44 superconductorswhose
y 2−x 2
coherentlengthξ = ~vF rangesfrom4ato12a,wherea
π∆(0)
is lattice constant, the maximum pairing amplitudes are IV. RESULTS AND DISCUSSION
taken to be g = 0.62, g = 2.60, and g = 1.28, respec-
0 1 2
tively, which result in two SC order parameters(eV) due
Thevortexstructuresforisotropicswavepairingstate
to orbital anisotropy in zero-field case as for isotropic s
are shown in Fig. 4 for different orbitals, respectively.
wave
Thevortexstatesexhibitorbitalanisotropy. Ford and
xz
|∆s (0)|=0.047; |∆s (0)|=0.026 (26) d orbitalsthe amplitudes havetwoplateauswithadif-
xz,yz xy yz
ference about 0.005eV along yˆ and xˆ directions on both
for anisotropic s wave
sides of the core region and the pinning center deviates
|∆anis. s(0)|=0.048; |∆anis. s(0)|=0.023 (27) slightly from the center of magnetic unit cell. The phase
xz,yz xy
distribution shows a winding number W = 1, such that
and for dx2−y2 wave the symmetry subgroup of the vortex structure is G532.
The winding structure of the s wave vortex, as mapped
|∆dx2−y2(0)|=0.048; |∆dx2−y2(0)|=0.025 (28)
xz,yz xy to a vector field, has a sink-type core center. Fig. 5
7
0.06
7 7
6 6
5 5
0.04 000...000000000 0246 01234 0000....000000000000 02468 01234
0.02 28 28
E(eV) 0 vorzteerxo s-ftiaetlde 28 24 20 16 12 8 ( 4a) 0 0 4 8 12 16 20 24 28 24 20 16 12 8 ( 4b) 0 0 4 8 12 16 20 24
-0.02
-0.04 67 67
-0.06 000000 ......0000000.00000000011220 04826482 012345 0000....00000000 01234 012345
2324 2328 2332 2336 2340 2344 2348 2352 2356 2360 2364 2368 2372 2376 2380 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28
FIG.5. (coloronline)EigenvaluesofBdGequationataround
Fermi level in the cases of isotropic s wave pairing state for (c) (d)
zero-field states, shown in green circles, and vortex states,
shown in red squares, respectively. The eigenvalues are plot-
ted in an ascending sequencein horizontal axis. 7 7
6 6
5 0.004 5
00..000012 234 00..000023 234
1 0.001 1
0 0 0 0
sshtaotwess,twheheerigeeintvhaalusebseoebntfaoiunneddftrhoemrevaorrete1x6ainnd-gzaeproe-igfieenld- 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28
states for both positive and negative eigenvalues. We (e) (f)
examine the behavior of the quasi-particle wavefunction
un and vn and it turns out that all the 32 in-gap FIG.6. (coloronline)Amplitudesandphases(colormapping)
iα↑↑ iα↓↑
states are extended to the entire magnetic unit cell(Fig. of quasi-particle wavefunctions uniα↑↑ and vinα↓↑ for isotropic
6, eigenstate |ǫ2353↑i). The orbital anisotropy again ap- s wave pairing symmetry of index n=2353 for dxz orbital (a)
and (b), d orbital (c) and (d), and d orbital (e) and (f),
pears as for d and d orbitals, the wavefunction ex- yz xy
xz yz
respectively.
tends to xˆ and yˆ direction because the spatial orienta-
tion of d-orbital harmonics, whereas for d orbital, the
xy
spreading of wavefunction is symmetric in xˆ and yˆ di-
rections. These extended wavefunctions amount to large
scalevariationoforderparameterswithintheentiremag- topy group of order parameter space of a vortex state is
netic unit cell and consequently a relatively large vortex π [U(1)] = Z and that the winding number N = 1 has
1 v
core region. been fixed when the self-consistent calculation is carried
Thestructuresofanisotropicswavevorticesareshown out, we expect that the variation along the loop around
in Fig. 7. The core regions of d orbital vortices are the corner is definitely not homotopic equivalent to that
xz/yz
notageometricpointanymore. Instead,theyhavebeen aroundvortexatcenter. Fig. 8(a)showsthephasevari-
stretched along xˆ and yˆ directions due to the fact that ation around the vortex core, where the phases change
although the pairing bonds are defined on next nearest slowly on a number of lattice sites at the very begin-
neighbor sites, the electrons forming Cooper pairs come ning of the loop as shown in Fig. 7 (a) in the vicinity
from distinguishable oriented orbitals. The symmetry of site (25,3). We have deliberately chosen a loop far
subgroup of anisotropic s wave vortices is still G , but away from the core region, since a stable topological de-
5
orbital asymmetry results in a line-type topological de- fectalwaysleavesitssignatureanywherearbitrarilyaway
fect for d orbital vortices, whereas d orbital vor- from it1. However, the phase variation of order parame-
xz/yz xy
tex is still of sink-type. One special fact worth noting is tersaroundthecornerofmagneticunitcellexhibitssome
that there is a suppression of order parameters at cor- turning-backpoints,fromwhichtheclockwiseincrements
ners of magnetic unit cell, which also exists for pairing contribute negative phase winding. Therefore the total
bond along −xˆ ± yˆ and xˆ − yˆ directions. In order to winding aroundthe corneris zero,whichprovesthat the
understand the physical origin of this phenomena, we suppression of order parameters at corners of magnetic
examine the phase variation along two loops around the unit cell is not a vortex. Detailed analysis about the
centerandcornerofmagneticunitcell,respectively. The phasedifferenceoneachlatticesitesshowsthatsuchsin-
looparoundthecorneriswell-definedinorderparameter gularitiesatcornersisactuallycausedbythediscontinu-
space because the nontrivial winding periodic boundary ityofboundaryconditionofwavefunctionofeachorbitals
condition Eq. (18) has been applied. Since the homo- when the calculationis carriedout ona N ×N lattice.
x y
8
0.05 90 0.05 30
28
80
28 0.07 0.03 70 0.03 25
21 000...000456 1241 Im- 00..0011 34560000 loop direction Im- 00..0011 112050 loop direction
14 -0.03 20 -0.03 5
0.03 10
-0.05 0 -0.05 0
7 0.02 7 -0.05 -0.03 -0.01Re 0.01 0.03 0.05 -0.05 -0.03 -0.01Re 0.01 0.03 0.05
0.01 (a) (b)
7 14 21 28 0 7 14 21 28
(a) (b) FIG. 8. (color online) Phase mapping onto complex plane of
anisotropic s wave pairing bond along xˆ+yˆdirection. Loop
28 aroundcenterofmagneticunitcell(a)is(25,3)→(25,25)→
28 0.07 (3,25)→(3,3)→(25,3) andaroundcornerofmagnetic unit
0.06 21 cell (b) is (3,1) → (3,3) → (1,3) → (28,3) → (25,3) →
(25,1) → (25,28) → (25,25) → (28,25) → (1,25) →
21 0.05 (3,25)→(3,28)→(3,1). Theloopdirectionhasbeenshown
0.04 14 bycolor mapping of each steps.
14
0.03
7 0.02 7
0.01
0.06
7 14 21 28 0 7 14 21 28
(c) (d)
0.04
28
28 0.03 0.02
21 00..00225 21 E(eV) 0 vorzteerxo s-ftiaetlde
14 -0.02
14 0.015
0.01
-0.04
7 7
0.005
-0.06
0 0 7 14 21 28 0 7 14 21 28
(e) (f) 2324 2328 2332 2336 2340 2344 2348 2352 2356 2360 2364 2368 2372 2376 2380
FIG.9. (coloronline)EigenvaluesofBdGequationataround
FIG. 7. (color online) Amplitudes(color mapping) and phase
Fermilevelinthecasesofanisotropicswavepairingstatefor
distribution of anisotropic s wave pairing bonds along xˆ+yˆ
zero-field states, shown in green circles, and vortex states,
direction for d orbital (a) and (b), d orbital (c) and (d),
xz yz shown in red squares, respectively.
and d orbital (e) and (f), respectively. Results of pairing
xy
bondsalong−xˆ+yˆ,−xˆ−yˆ,andxˆ−yˆdirectionsarethesame
with these results.
variant subgroups of magnetic translation group, which
isoriginallyaimedatdescribingtheAbrikosovlatticefor
From Eq. (18), we know that the variation of boundary N = 1. There are 12 in-gap eigenstates, as shown in
v
condition along xˆ direction for adjacent (λ =1,λ =0) Fig. 9, which locate symmetrically on both sides of the
x y
magnetic unit cell is eiKNxiy, and it will come back to Fermi level. The wavefunctions of these states are typ-
ei(KNx+2π) when the condition iy =4Ny+1 is satisfied. ically localized for dxz/yz orbitals and extended for dxy
Itisobviouslythatsuchaconditioncannotberealizedin orbital, as shown in Fig. 10 for eigenstate |ǫ i. It has
2353
numericalcalculationforanygivenN ,thereforethedis- been observed that the wavefunctions for each orbitals
y
continuity,whichcanberegardedasanimpurityinduced show particle-hole asymmetry. Although the difference
by winding boundary condition, cannotbe avoided. The of vortices between isotropic s and anisotropic s wave
impurity nature of these singularities can also be recog- pairing states has been observed from the hitherto re-
nized as the suppression of order parameters occurs on sults, such a difference may rely on the limitation of our
single site at corners, which is different from a genuine modelcalculationinthatsincetheHamiltonianisdefined
vortex having an effective core region. We also noted on site-orbital representation, there is no well-defined k-
that such a singularity does not exist for N = 4, but space energy cut-off in the vicinity of the Fermi level for
v
in this case the vortex states cannot be classified by in- the attractive pairing potential. Therefore, pairing elec-
9
28
28 ’./op_r_d_orb1.dat’ u 1:2:3 0.08
7 7
000...000000123 123456 0000....000000001234 123456 21 00..0067 21
0 0 0 0
0.05
28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 14 00..0034 14
(a) (b) 7 0.02 7
0.01
7 14 21 28 0 0 0 7 14 21 28
7 7 (a) (b)
0.001 56 56
4 0.001 4
0.0005 123 0.0005 123 28
0 0 0 0 28 ’./op_r_d_orb1.dat’ u 1:2:4 0.05
28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 21 0.04 21
(c) (d) 0.03
14
14
0.02
7
67 67 7 0.01
00000.....0000000000 012345 012345 0000....00000000 02468 012345 7 14 21 28 0 0 0 7 14 21 28
28 28 (c) (d)
28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28 24 20 16 12 8 4 0 0 4 8 12 16 20 24 28
(e) (f) 28 ’./op_r_d_orb2.dat’ u 1:2:3 0.06
0.05 21
21
FIG. 10. (color online) Amplitudes and phases(color map-
0.04
ping) of quasi-particle wavefunctions un and vn for
iα↑↑ iα↓↑ 14
anisotropic s wave pairing state of index n=2353 for dxz or- 14 0.03
bital (a) and (b), d orbital (c) and (d), and d orbital (e)
yz xy 0.02
and (f), respectively. 7
7
0.01
7 14 21 28 0 0 0 7 14 21 28
(e) (f)
trons may come from the region which is far away from
the four electron pockets. Consequently, the absence of 28
pocket at Γ point may induce ambiguity for anisotropic 28 ’./op_r_d_orb2.dat’ u 1:2:4 0.08
s wavepairing state in a frameworkof BCS-type pairing 0.07 21
scheme. 21 0.06
Results of dx2−y2 wave vortices are different from A1g 0.05 14
vortices discussed above in many aspects. The orbital 14 0.04
anisotropy dominates the vortex structures. Fig. 11 and 0.03
7 0.02 7
0.01
TABLEII.Values(inunit: 10−1eV)oforbital-resolveddx2−y2 7 14(g) 21 28 0 0 0 7 ( h14) 21 28
wave pairing order parameters(pairing bonds) π and σ
x,y x,y
asdefinedinFig. 13forsite(3,3)forzero-fieldSCandvortex
FIG.11. (coloronline)Amplitudes(colormapping)andphase
states. The spin and site indices havebeen omitted.
distributionofdx2−y2 wavepairingbondsfordxzorbitalalong
xˆ direction (a) and (b), yˆ direction (c) and (d), and for d
yz
Zero-field SC state Vortex state orbital along xˆ direction (e) and (f), yˆdirection (g) and (h),
∆xz(σx) (0.43, 0.43) (-0.12, 0.58) respectively. Results of pairing bonds along the other two
∆xz(πy) (0.032, 0.032) (-0.34, 0.11) directions of next nearest site pairing are same with these
∆yz(πx) (-0.032, -0.032) (-0.12, 0.34) results.
∆ (σ ) (-0.43, -0.43) (-0.58, 0.11)
yz y
∆ (xˆ) (0.17, 0.17) (0.091, 0.24)
xy
∆ (yˆ) (-0.17, -0.17) (-0.24, 0.096)
xy
12showthe amplitudes andphase distributionof dx2−y2
10
28
28 ’./op_r_d_orb3.dat’ u 1:2:3 0.04 d(cid:13)yz(cid:13)(cid:13) d(cid:13)yz(cid:13)(cid:13)
0.035
21 d(cid:13) d(cid:13)
xz(cid:13)(cid:13) xz(cid:13)(cid:13)
21 0.03
0.025
14
14 0.02 π(cid:13) σ(cid:13) y(cid:13)
0.015 y(cid:13) y(cid:13)
7 π(cid:13) x(cid:13)
7 0.01 x(cid:13)
0.005
7 14 21 28 0 0 0 7 14 21 28 d(cid:13)xz(cid:13)(cid:13) d(cid:13)xz(cid:13)(cid:13)
(a) (b) σ(cid:13)
d(cid:13)yz(cid:13)(cid:13) x(cid:13) d(cid:13)yz(cid:13)(cid:13)
28
28 ’./op_r_d_orb3.dat’ u 1:2:4 0.04
0.035
21 FIG. 13. (color online) A schematic picture illustrates that
21 0.03 dx2−y2 wave pairing state is re-defined between dxz and dyz
0.025 orbitals dueto4-fold rotational symmetry. Thered andblue
14
14 0.02 colorindicatepositiveandnegativesignsoforbitalwavefunc-
0.015 tions. The long and short double-headed arrows correspond-
7 ing π and σ pairing bondsalong xˆ and yˆdirections show the
7 0.01
0.005 exchangeof orbital states underC4z rotation.
7 14 21 28 0 0 0 7 14 21 28
(c) (d)
0.1
FIG.12. (coloronline)Amplitudes(colormapping)andphase
distributionofdx2−y2 wavepairingbondfordxy orbitalalong
xˆ direction (a) and (b), and yˆ direction (c) and (d), respec-
0.05
tively. Resultsofpairingbondsalongtheothertwodirections
of next nearest site pairing are same with theseresults.
0
V)
E(e vorzteerxo s-ftiaetlde
wavepairingbondsforeachorbitals. Ithasbeenpointed -0.05
outinprevioussectionthatthesymmetryofbandstruc-
ture gives constraints to symmetry of pairing states. A
-0.1
strong hybridizationof d and d orbitals, as shownin
xz yz
PDOSinFig. 2,resultsinare-defineddx2−y2 wavepair-
ing state, as shown in Fig. 13, since the wavefunctions
2333 2338 2343 2348 2353 2358 2363 2368 2373
ofthese twoorbitalstransformunderactionofgenerator
C4z as FIG. 14. (color online) Eigenvalues of BdG equation at
around Fermi level in the cases of dx2−y2 wave pairing state
C4z|dxzi=|dyzi (29) forzero-fieldstates,showningreencircles,andvortexstates,
C |d i=−|d i shown in red squares, respectively.
4z yz xz
while d orbital does not mix with them under such a
xy
transformation. Here we give an example of numerical
resultsoforderparametersforeachorbitalsonsite(3,3), of magnetic field, the signchange of dx2−y2 wavepairing
as shown in Table II. In zero-field case, phase difference symmetry, along with the orbital-hybridized order pa-
of eiπ is observed between π and π , σ and σ bonds, rameters together give rise to a G∗ winding structure
x y x y 6
whicharedefinedondifferentorbitals,whereasinvortex for d orbitals, which seems like a solenoidal vector
xz/yz
states, such a phase will undergo a gauge modification field. Such phase difference has been observed between
which is induced by magnetic field. The winding struc- ∆ (σ ) as shown in Fig. 11 (b) and ∆ (σ ) as shown
xz x yz y
tures showninFig. 11 and12for differentorbitals share in Fig. 11 (h), and also between ∆ (π ) as shown in
xz y
this common feature for all order parameters defined on Fig. 11(d)and∆ (π )asshowninFig. 11(f). Among
yz x
entire magnetic unit cell. For d orbital, G∗ vortices 17(positive)in-gapstatesassociatedwithorbital-resolved
xy 5
which are defined on pairing bonds ∆xy(xˆ) and ∆xy(yˆ) dx2−y2 wave vortices as shown in Fig. 14, the wavefunc-
areofsink-andsource-type,respectively,becausetheor- tions of eigenvalue |ǫ i for d and d orbitals, and
2353↑ xz yz
der parameters change sign as they transform according |ǫ iford orbitalareshowninFig. 15. Theparticle-
2354↑ xy
toB irreducibleunitaryrepresentation. Inthepresence hole asymmetry is evidently for d and d orbitals in
1g xz yz