Superconducting Phases in Lithium Decorated Graphene LiC 6
Rouhollah Gholami, Rostam Moradian, Sina Moradian, Warren E. Pickett
aa r X i v : . [ c ond - m a t . s up r- c on ] S e p Superconducting Phases in Lithium DecoratedGraphene LiC Rouhollah Gholami , Rostam Moradian , Sina Moradian , and Warren E. Pickett Physics Department, Faculty of Science Razi University, Kermanshah, Iran Nano science and nano technology research center, Razi University, Kermanshah, Iran * Correspondence and requested materials should be addressed to R. M. [[email protected]] Department of Electrical and Computer Engineering, University of Central Florida, Orlando, Florida, USA. Department of Physics UC Davis, One Shield Avenue, Davis, CA 95616, USA.
ABSTRACT
A study of possible superconducting phases of graphene has been constructed in detail. A realistic tight binding model,fit to ab initio calculations, accounts for the Li-decoration of graphene with broken lattice symmetry, and includes s and d symmetry Bloch character that influences the gap symmetries that can arise. The resulting seven hybridized Li-C orbitalsthat support nine possible bond pairing amplitudes. The gap equation is solved for all possible gap symmetries. One bandis weakly dispersive near the Fermi energy along Γ → M where its Bloch wave function has linear combination of d x − y and d xy character, and is responsible for d x − y and d xy pairing with lowest pairing energy in our model. These symmetries almostpreserve properties from a two band model of pristine graphene. Another part of this band, along K → Γ , is nearly degeneratewith upper s band that favors extended s wave pairing which is not found in two band model. Upon electron doping to a criticalchemical potential µ = . eV the pairing potential decreases, then increases until a second critical value µ =1.3 eV at whicha phase transition to a distorted s -wave occurs. The distortion of d - or s-wave phases are a consequence of decoration whichis not appear in two band pristine model. In the pristine graphene these phases convert to usual d-wave or extended s -wavepairing. Two dimensional superconducting phases have become of great interest since the discovery of the high temperature supercon-ducting (HTS) cuprates and subsequent finding of Fe-pnictide and -chalcogenide HTSs. Interest was re-invigorated by thediscovery of superconductivity onsets up to 75K in single layer FeSe grown on
SrTiO and related substrates. With theenormous research activity focused on graphene in recent years, it is not surprising that graphene-based superconductivityhas become an active area of research. Very recently superconductivity up to 1.7K has been reported in magic angle bilayergraphene, which will buttress activity on two dimension superconductors and especially the related type that we discuss here.Superconductivity has been known for some time in intercalated graphite compounds such as C Ca and C Y b . With themany remarkable properties of graphene, it has been anticipated that doping by gating or by decorating with electro-positiveelements, thereby moving the chemical potential away from the Dirac points, might induce superconductivity. However,graphene decorated with alkali metals has three valence bands with one weakly dispersive band near Fermi energy. Due tothis flat band, there are additional available states around the Fermi level and the required pairing potential is reduced.Discussion of superconductivity in doped graphene has been primarily within theoretical models, as we review below,but some encouraging data have been reported. Experimental evidence for a superconducting gap in Li-decorated monolayergraphene around 6 K has been reported by Ludbrook et al. based on angle-resolved photoemission spectroscopy (ARPES).Scanning tunneling spectroscopy (STS) was applied by Palinkas et al. to graphene suspended on tin nanoparticles, whoconcluded that superconductivity is induced in the graphene layer. Evidence of superconductivity in Li-decorated few layergraphene at 7.4 K has been reported by Tiwari and collaborators . Low temperature mobility of K and Li atoms on graphenewas observed by Woo et al. , and suggest that mobility may persist at lower temperatures, which would provide new challengesfor theory.Various mechanisms of pairing have been proposed. Uchoa and Castro-Neto modeled pristine and doped graphene withelectron-phonon coupling or plasmon mediated in mind . Repulsive electron-electron interactions were modeled by Nand-kishore and collaborators. Beginning from pristine graphene, varying the chemical potential leads to dominant chiralsinglet d x − y + id xy pairing for nearest neighbor interaction, according to Black-Schaffer et al. a triplet f -wave state hasbeen proposed to arise from next-nearest neighbor interaction with chemical potential near van Hove peak. Both chiral andconventional p -wave states in graphene have been discussed , with the many pictures raising various possibilities but little of certain nature.More specific predictions have begun to appear. Profeta et al. predicted based on Eliashberg theory that decorationby electron donating atoms such as Ca and Li would make single layer graphene superconducting, with modest criticaltemperatures in the 1-8 K range. In somewhat related work, Wong et al. have predicted from an ab initio treatment a criticaltemperature around T c =14K for carbon nanotubes, which was increased to above 100 K for a certain type of carbon ring.Expectations of adjusting the chemical potential include gating, but the main focus has been on decoration of grapheneby electropositive atoms, viz. alkalis or alkaline earths. Charge migration from such decorating atoms to the graphene layerwill affect the C-C bonding, leading to contraction or expansion of the graphene hexagons that are centered by the decoratingatoms, thus breaking the symmetry of C-C hopping integrals around the honeycomb loop. This asymmetric graphene layerwill be referred to in this paper as “shrunken graphene”. Taking LiC for illustration, each cell site has six C atoms in ahexagon with an alkaline atom lying above the center of the hexagon. The C π orbitals and alkaline atom’s s -orbital hybridizeto give seven “molecular” orbitals. For two dimensional graphene-like structures effects, differences in nearest neighbourhopping integrals affect the band structure near the important Dirac point, which is folded back to the Γ point of the shrunkengraphene superlattice investigated by Hou et al. and Long-Hua et al. . For such systems not even the full analytic tight-binding band structure has yet been reported. The intent here is to extend study of this system, with representative LiC , fromthe underlying electronic structure to investigation of the possible superconducting phases.The organization of the paper is as follows. In Sec. II the interacting seven orbital model Hamiltonian is presented. Theexact band structure of the normal state of this shrunken graphene system is described in Sec. III. Perturbation theory is appliedto obtain the band structures in analytic form. Applying the Hubbard model and minimizing free energy of the superconductorstate, we obtain in Sec. IV the gap equations and approximate critical temperature. These equations are solved analytically toestablish the possible pairing symmetries and other properties of the superconducting states. A summary is provided in Sec.V. Because the unit cell contains several atoms with important specific aspects, we provide many of the details of the expressionsthat can be obtained analytically. LiC , as illustrated in Fig. 1, consists of a graphene layer decorated by a lithium layer inwhich Li atoms are located at the center of a carbon hexagon surrounded by six empty center hexagons. The height of Liabove the carbon layer is calculated to be h z = .
85, somewhat smaller than the value 1.93Å obtained by Profeta et al. . Thenearest Li-C distances are h = .
40. Since the Li 2 s orbital energy is higher than the C 2 p z orbital, charge transfer occurs. It iscalculated that 0 . e from Li transfers to the six C atoms equally. The positive Li ion and negative C ion provide a relativeCoulomb (Madelung) shift in site potentials of the two atoms.The attractive interaction between Li and C ions after charge transfer contracts the Li-C distance and reduces the C-C bondlengths in the Li-centered hexagon to a = . a = . t , with that between stretched carbon sites is denoted t ′ . We refer to thisbroken symmetry situation as “shrunk graphene”. The difference in hopping amplitudes indicates that the new Li-C hoppingparameter is the central new feature in LiC compared to graphene. Symmetry breakdown leads to the opening of a smallenergy gap at the Γ point.The lattice then becomes a two dimensional hexagonal Bravais lattice with seven atomic sites. These will be labeled as A , A , A , B , B , B and Li , as illustrated in Fig. 1. The Hamiltonian of this system isˆ H = − ∑ i α ∑ j β , σ t σ , σ i α , j β ˆ c † i ασ ˆ c j βσ + ∑ i α , σ ( ε i α − µ o ) ˆ n i ασ + ∑ i α , σ ∑ j β , σ ′ U σ , σ ′ i α , j β ˆ n i ασ ˆ n j βσ ′ = ˆ H N + ˆ H P . (1)Here H N and H P denote the non-interacting and interaction Hamiltonians respectively. In these expressions α and β run over A i , B i and Li . Here ˆ c † i ασ , ˆ c i ασ are creation and annihilation operators of an electron with spin σ on subsite α of i th lattice site,and ˆ n i ασ = ˆ c † i ασ ˆ c i ασ is the electron number operator. The noninteracting chemical potential is µ and t i α , j β is the hoppingintegral from the α site of i th cell to the β site of j th cell. We denote the on-site energy by ε α .The interaction stated above corresponds to an extended (negative U ) Hubbard model, which allows a variety of phe-nomenological values to be chosen and studied. It is largely for this reason that we provide substantial detail of the underlying,non-interacting C-Li lattice and electronic structure. The interactions that we study are introduced in Sec. IV. Many studies of graphene rely on tight binding parametrization of the band structure. The early parametrization of Wallace already employed both first and second neighbors. Extensions in various ways have followed, culminating in the applica- ion of Wannier functions by Jung and MacDonald to provide simple but realistic five parameter model and a more accuratebut more involved 15 parameter model. Our aim in this section is to construct a realistic seven band model for distorted LiC ,while also developing the formalism to allow exploration of superconducting phases once the interaction has been included.The distortion of the graphene layer to shrunken graphene and the coupling to Li requires a considerable generalizationof the underlying tight binding model Hamiltonian, and many of the details are relegated to appendices. The Hamiltonian ofnon-interacting LiC isˆ H N = − ∑ i α ∑ j β , σ t σ , σ i α , j β c † i ασ c j βσ + ∑ i α , σ ( ε i α − µ o ) ˆ n i ασ . (2)Eq. 2 incorporates broken symmetries in the on-site energies, hopping integrals, and bond lengths. Here, it has been assumedthat on site energies ε A i = ε A and ε B i = ε B . It is diagonalized in terms of Bloch eigenfunction of the form Eq. A.2 . In matrixrepresentation, the equation for the coefficients becomes ε ( ~ k ) d c ( ~ k ) d c ( ~ k ) d c ( ~ k ) d ∗ c ( ~ k ) d ∗ c ( ~ k ) d ∗ c ( ~ k ) d ∗ c ( ~ k ) ε ( ~ k ) β ( ~ k ) γ ( ~ k ) τ ( ~ k ) d ( ~ k ) d ( ~ k ) d ∗ c ( ~ k ) β ∗ ( ~ k ) ε ( ~ k ) θ ( ~ k ) d ( ~ k ) τ ( ~ k ) d ( ~ k ) d ∗ c ( ~ k ) γ ∗ ( ~ k ) θ ∗ ( ~ k ) ε ( ~ k ) d ( ~ k ) d ( ~ k ) τ ( ~ k ) d c ( ~ k ) τ ∗ ( ~ k ) d ∗ ( ~ k ) d ∗ ( ~ k ) ε ( ~ k ) β ∗ ( ~ k ) γ ∗ ( ~ k ) d c ( ~ k ) d ∗ ( ~ k ) τ ∗ ( ~ k ) d ∗ ( ~ k ) β ( ~ k ) ε ( ~ k ) θ ∗ ( ~ k ) d c ( ~ k ) d ∗ ( ~ k ) d ∗ ( ~ k ) τ ∗ ( ~ k ) γ ( ~ k ) θ ( ~ k ) ε ( ~ k ) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) = E i ( ~ k ) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) C ( E i ( ~ k )) (3)where d ci ( ~ k ) , ε i ( ~ k ) , β ( ~ k ) , θ ( ~ k ) , γ ( ~ k ) , d i ( ~ k ) and τ i ( ~ k ) functions are defined in supplementary materials Eqs. A.7, A.8, A.9, andA.10 respectively. For general ~ k vectors, it is challenging to obtain an exact analytical expression for the full Hamiltonianin Eq. 3 and it would not be transparent anyway. However, analytical expression for Eq. 3 can be achieved in two steps.Since hopping from Li atoms to nearest neighbor carbon sites t LiC is small with respect to C-C nearest neighbor hopping t , by first neglecting the lithium-carbon hopping t LiC →
0, first column and row of the Hamiltonian matrix in Eq. 3 areomitted, the remaining part given by Eq. B.1 is uncoupled shrunken graphene Hamiltonian which can be diagonalized exactlyto obtain E sh , n . Finally, Li-C coupling is taken into account by perturbation theory to obtain eigenvalues E n , as presented inthe appendices. C Dispersion Relations
By first neglecting the lithium-carbon hopping, t LiC →
0, the uncoupled shrunken graphene Hamiltonian given by Eq. B.1 canbe diagonalized exactly. Even though Li-C hopping has been neglected but still remaining part of shrunken Hamiltonian inthe most general case, include broken symmetries in the hopping integrals, bond lengths and on-site energies. The non trivialeigenvalues of uncoupled shrunken graphene Hamiltonian in general form are given by E sh , ml ( t i ,~ ξ i ,~ k ) = − µ o + α ( ~ k ) + u m Π ( ~ k ) + u ∗ m Π ∗ ( ~ k ) + (cid:20) ε A + ε B + ( − ) l q ( ε A − ε B ) + w m ( ~ k ) (cid:21) (4)with details presented in supplementary materials Appendix B. However, the obtained equations are often complicated. Toprovide insight into the method, uncoupled shrunken graphene Hamiltonian can diagonalized in some particular cases. TheBrillouin zone (BZ) of C is one third of that of graphene, with the Dirac points folded back to the Γ point. In this mini-BZ, thetwo π bands of pristine graphene i.e. E ± = ± t | η | folds to six branches as illustrated inFig. 2. These branches are solutionsof Eq. B.1 in the limited case of pristine which in the nearest neighbor approximation they are given by, E ± γ ( ~ k ) = ± t | η ( ~ k ) | , η ( ~ k ) = e i ~ k .~ δ + e i ~ k .~ δ + e i ~ k .~ δ E ± β ( ~ k ) = ± t | η ( ~ k ) | , η ( ~ k ) = e i ~ k .~ δ + e i π e i ~ k .~ δ + e i π e i ~ k .~ δ E ± α ( ~ k ) = ± t | η ( ~ k ) | , η ( ~ k ) = e i ~ k .~ δ + e i π e i ~ k .~ δ + e i π e i ~ k .~ δ . (5)Exact analytical solutions for pristine graphene wherein next neighbor hopping integrals are taken into account are presentedin supplementary materials Eqs. B.7 and B.8. As shown in Fig. 2 one sees that E ± β ( ~ k ) is weakly dispersive near the vanHove singularity at the saddle points M at 3/8 or 5/8 filling (0.25 electron per carbon doping), this band plays a major rolein the formation of superconductivity in graphene. Also, one can observe that the band structure is four-fold degenerateat the charge neutral Dirac points. Solution of the Schrödinger equation for pristine graphene in the mini-BZ has anotheradvantage: the Bloch-wave symmetry character of each branch can be distinguished. The Bloch coefficients of the branch abeled by E γ are of s -wave character, C A i = ( , , ) while for those labeled as E α and E β are of the form d ± id -wavei.e. C A i = ( , e ± i i π , e ± i i π ) as illustrated inFig. 2 and demonstrated in more detail in Appendix B , Eqs. B.4 and B.6. Thisbecomes important when it is shown that different superconducting phases of graphene in a variety of doping regimes are dueto electron pairing in each of these branches.Decoration of graphene with metals reduces symmetries that lead to removal of bands degeneracy in some regions. Whiledecoration causes expansion and contraction of bonds length in three inequivalent directions in the honeycomb lattice i.e. | ~ τ i | 6 = | ~ δ i | , eigenenergies E sh , ml ( t i ,~ ξ i ,~ k ) in Eq. 4 do not depend on the bond lengths ~ τ i and ~ δ i separately but are functionsof LiC lattice bases length | ~ ξ i = ~ τ i + ~ δ i | , so symmetry breakdown of bond lengths does not break symmetries of bands.Symmetry reduction of hopping integrals removes degeneracies occurring in pristine graphene band structure, with the mostimportant effect being to open a gap E g = | t ′ − t | at the Dirac point which has been folded back to the Γ point. This gaparises from symmetry breaking of the nearest neighbor hopping and dose not affected by the other next neighbors hopping norby the Li-C hopping integral. Comparison with DFT band structures gives E g = . eV . Another gap can arise at the Γ pointbecause of symmetry breaking of on-site energies ε A = ε B , seen from Eq. 4. For the case t = t ′ the gap becomes 2 | ε A − ε B | .In Li decorated graphene that we consider here, all carbon on-site energies are equal so this type gap does not arise.While for folded but pristine graphene Bloch wave solutions are pure s -wave or chiral d ± id -wave and there are nomixed states, when symmetries in hopping integrals are broken by decoration, Bloch functions are linear combinations ofall these phases, Eq. B.6. Equation B.7 demonstrates that for a general ~ k all probabilities are equal in pristine graphene i.e. | C A i ( E m ) | = | C B i ( E m ) | = . In shrunken graphene these probabilities are ~ k dependent and unequal in general. It will be seenthat these small deviations influence the superconducting gap equation symmetries. Dispersion Relations
Li-C hopping adds a perturbation term to the shrunken graphene Hamiltonian. Obtaining exact dispersions from Eq. 3 is verychallenging, so perturbation theory is applied to obtain approximate solutions, as presented in Appendix C. However, to getsome insight into effects of the coupling, Eq. 3 can be solved exactly at the Γ point. At ~ k =0 only the isolated Li band, E Li , ( ) and the lowest valance band, E sh , ( ) , are mutually affected. The energies of these bands are, with E ( ) ≡ E + , E ( ) ≡ E − , E ± ( ) = (cid:0) E Li , ( ) + E sh , ( ) (cid:1) ± s (cid:2) E Li , ( ) − E sh , ( ) (cid:3) + ( t LiC ) (6)and other shrunk graphene bands given by (supplementary) Eq. B.5 remain unchanged. Comparing the fit results from DFTto these equations suggests that t Li − C is in the 0.3-0.5 eV range, and other next neighbor hopping from Li atoms to C sites arenegligible.There are two critical points in the pure graphene band structure which are affected by decoration and become important:the charge neutrality Dirac points folded at the Γ point, and the van Hove singularity at the M point. We define a hoppingintegral symmetry breaking index, w t = t ′ t = w t open a small gap E g at Γ , whichdoes not depend on t LiC . Depending on doping level, Li-C hopping affects the band structure near the points that the isolatedLi band E Li , ( ~ k ) and uncoupled shrunken graphene bands intersect. These impurity effects causes not only changes in energylevel but alter the density of states. Superconductivity emerges from pairing of electrons near the Fermi energy and it isimportant to know how the density of states at the Fermi energy N ( ) changes with decoration. The seven band tight binding model of LiC was fit to the DFT band structure, with results illustrated in Fig. 3. In the graphenelayer shown in Fig. 1(a) and (c), A subsite chosen as central site labeled by 0 and B subsite in adjacent hexagon consideredas second neighbor while just slightly longer than the first neighbors atoms B and B in same hexagon, this neighbor labeledby n = n th neighbor has been shown by t CC n .In-plane Li-Li hopping, t LiLi m obtained up to m = ]. The fitted hopping amplitudes and on-site energies are presentedin Tables 1. Note that by comparing band structure of LiC with pristine graphene in ref. , it is observed that Li decorationonly slightly changes the pristine graphene band structure. These changes are due to electron transfer from Li to graphene,which changes the pristine on site ε pristine = ε A = ε B = ε c . Superconducting Pairing and States
LiC presents a multiband system in which three bands cross the Fermi level. We presume singlet pairing that can be bothintraband and interband in nature. We adopt a local viewpoint in which pairing occurs between electrons on carbon atoms.Seven hybridized Li-C orbitals, support nine possible bond pairing amplitudes in real space. Figure 4 (a) illustrates allthe nearest neighbour order parameters possibilities. Leaving the analytical derivation details to supplementary materialsAppendices D and E, the quasiparticle energies are obtained by Bogoliubov-de Gennes unitary transformation in the sevenband space, E Qm , s ( ~ k ) = s E m ( ~ k ) + ∑ i = (cid:12)(cid:12)(cid:12) ∆ mi ( ~ k ) (cid:12)(cid:12)(cid:12) E m ( ~ k ) + E i ( ~ k ) s = ± s = s = − E m are the normal state eigenvalues. The ~ k -dependent gap | ∆ mi ( ~ k ) | in the spectrum are expressed as ∆ mn ( ~ k ) = ∑ α = Ω α mn ( ~ k ) ∆ α (8)in which m and n are band indexes. The band pair order parameter ∆ mn ( ~ k ) denotes pairing between electrons in the m -th and n -th bands in LiC . Also, ( ∆ , ∆ , ∆ ) = ( ∆ ′′ , ∆ ′′ , ∆ ′′ ) ; ( ∆ , ∆ , ∆ ) = ( ∆ , ∆ , ∆ ) ; ( ∆ , ∆ , ∆ ) = ( ∆ ′ , ∆ ′ , ∆ ′ ) are shownin Fig. 4(a), and Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ τ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ τ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ τ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ τ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ τ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ τ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ δ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ δ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ δ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ δ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ δ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ δ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ δ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ δ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ δ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ δ Ω i j ( ~ k ) = C ∗ ( E i ) C ( E j ) e i ~ k .~ δ + C ∗ ( E i ) C ( E j ) e − i ~ k .~ δ . (9)where C i ( E j ) are Bloch wave coefficients of the j -th band. Possible order parameter symmetries in Eq. 8 are related tosymmetries of Bloch wave functions, through Ω i j ( ~ k ) functions in Eq. 9. In the limiting case of (folded) six band pristinegraphene, the symmetry character of different conduction bands along high symmetry lines were provided in Fig. 2. Blochsymmetry character of non-interacting bands specifies the symmetry of the band order parameter. The linearized gap equation, obtained by minimizing the quasiparticle free energy with respect to nearest neighbor orderparameters, is J β ∆ β = − N ∑ α = ∑ ~ k ∑ n = ∑ i = tanh ( E Qn k B T ) E n ( ~ k ) + E i ( ~ k ) (cid:16) Ω α ni ( ~ k ) Ω ∗ β ni ( ~ k ) + Ω β ni ( ~ k ) Ω ∗ α ni ( ~ k ) (cid:17) ∆ α ≡ − ∑ α = Γ βα ∆ α . (10)This equation can be written in matrix form as A × B × B × B × C × D × B × D × C × g V g V g V = − V V V (11) here A × = Γ Γ Γ Γ Γ Γ Γ Γ Γ , C × = Γ Γ Γ Γ Γ Γ Γ Γ Γ , B × = Γ Γ Γ Γ Γ Γ Γ Γ Γ , D × = Γ Γ Γ Γ Γ Γ Γ Γ Γ (12)and g V = (cid:16) ∆ ′′ ∆ ′′ ∆ ′′ (cid:17) T , g V = ( ∆ ∆ ∆ ) T and g V = (cid:16) ∆ ′ ∆ ′ ∆ ′ (cid:17) T ; the subscripts < i j > has been dropped for brevity.The A × , B × , C × , and D × matrices, given by Eq. 12, have identical structures, hence they share eigenvectors: V Ts =( ) , V Td xy = ( − ) , and V Td x − y = ( − ) , where the latter two are degenerate. Their eigenvalues, in obviousnotation, are a s = Γ + Γ , b s = Γ + Γ , c s = Γ + Γ , d s = Γ + Γ a d = Γ − Γ , b d = Γ − Γ , c d = Γ − Γ , d d = Γ − Γ . (13)For folded six band pure graphene g = g , the Bloch wave coefficients appearing in Eq. 9 can be replaced from Eq. B.7 toshow that Ω i j ( ~ k ) = Ω ji ( ~ k ) = Ω i j ( ~ k ) and similarly relations for other elements, hence C × = A × and D × = B × . Eq. 11takes the more symmetric form A × B × B × B × A × B × B × B × A × V V V = − g V V V (14)For the case V = V = V = V sy where sy subscripts indicates each of the s , d xy or d x − y symmetry, the six band gap Eq. 14 re-duces to ( A + B ) V sy = − g V sy , i.e. the linearized gap equation of the two band model of pristine graphene in Ref. [ ]. Thesethree solutions preserve symmetry of the two band unit cell as illustrated in Fig. 4 (b), (c). In addition to these three states,there are six more non-orthogonal solutions Φ n = ( V sy − V sy ) and Φ n = ( V sy − V sy ) that break symmetries of pristinegraphene two band model. Inserting these solutions into Eq. 14 leads to a new two band gap equation, ( A − B ) V sy = − g V sy ,which is unphysical because of an unreachably high energy pairing potential g . In the following section the superconductinggap equation has been solved for LiC and it is demonstrated how Li-C coupling influences superconducting phases. Self-consistent solutions of the gap equation Eq. 11 can be obtained analytically. There are three superconducting states withisland character (discussed in more detail below) that can be expressed in compact form as [ Φ n ] T = [ V sy − V sy ] , J sy = c sy − d sy (15)where V sy refers to one of the V s , V d xy or V d x − y -wave symmetries. Pairing in these phases cannot propagate, as may be picturedin Fig. 5. The other six superconducting states of Eq. 11 have the explicit form [ Φ n ] T = [ α ± sy V sy V sy V sy ] (16)corresponding to the interaction potential is g = J sy wherein α ± sy = J ± sy − c sy − d sy b sy , J ± sy = (cid:18) κ a sy + c sy + d sy ± q κ b sy + [ c sy + d sy − κ a sy ] (cid:19) (17)In these expressions J ± d ( s ) and α ± d ( s ) are obtained from Eq. 17 by substituting a sy , b sy , c sy , and d sy by a d ( s ) , b d ( s ) , c d ( s ) , and d d ( s ) respectively.By comparing the gap equations introduced in Eq. 11 and Eq. 14 the gap equation symmetry reduction of decoratedgraphene with respect to folded but pristine graphene becomes clear. This symmetry reduction results in an α sy coefficientappearing in the pairing amplitudes of stretched bonds as shown in Eq. 16 and Fig. 4. We refer to these symmetry reductionphases as “distorted phases”.The six bands of pristine graphene support nine pairing amplitudes while in the two band model there are three possiblepairing amplitudes along three different bonds. These two notions can be mapped onto each other only if α sy = κ = a sy = c cy and b sy = d cy from Eq. 17, and it follows that if b sy > α + sy = α − sy = −
2. Also g + < g − so in this case (+) signpreserves the two band model while the ( − ) sign phases are unphysical. Numerical calculation shows that b + sy >
0. Thesesuperconducting states can be categorized into three groups according to their corresponding pairing potential. electron pairing states with island character and very high pairing potential; (unphysical solutions) Φ p x = [( ) ( − )( − )] T , g = c d − d d , Φ p y = √ [( ) ( − )( − − )] T , g = c d − d d , Φ f = √ [( ) ( )( − − − )] T , g = c s − d s . (18) Φ − d x − y = (cid:0) ( α − d ) + (cid:1) − (cid:2) α − d ( − ) ( − ) ( − ) (cid:3) T , g = J − d Φ − d xy = (cid:0) ( α − d ) + (cid:1) − (cid:2) α − d ( − ) ( − ) ( − ) (cid:3) T , g = J − d Φ − s = (cid:0) ( α − s ) + (cid:1) − (cid:2) α − s ( ) ( ) ( ) (cid:3) T , g = J − s (19) Φ + d x − y = (cid:0) ( α + d ) + (cid:1) − (cid:2) α + d ( − ) ( − ) ( − ) (cid:3) T , g = J + d Φ + d xy = (cid:0) ( α + d ) + (cid:1) − (cid:2) α + d ( − ) ( − ) ( − ) (cid:3) T , g = J + d Φ + s = (cid:0) ( α + s ) + (cid:1) − (cid:2) α + s ( ) ( ) ( ) (cid:3) T , g = J + s . (20)All states are orthogonal except those with same subscript, viz. Φ − s and Φ + s . Such solutions are orthogonal if κ = i.e.g = g . Only for this case the matrix gap equation becomes Hermitian, then band order parameters takes the following formin terms of the band Green function and g , ∆ i j ( ~ k ) = g h d ↑ i ( ~ k ) d ↓ j ( ~ k ) i . (21)Here ˆ d σ i ( ~ k ) = ∑ m = C ∗ m ( E i ( ~ k )) ˆ c σ m ( ~ k ) annihilates an electron with spin σ in the i th band with energy E i ( ~ k ) . Although it isassumed that g = g but deviation from pristine leads to distortion of Green’s functions < ˆ c σ i α ˆ c σ j β > along different bonds.Phases Φ − d x − y and Φ − d xy are degenerate with eigenvalue J − d x − y = J − d xy = J − d , and similarly Φ + d x − y and Φ + d xy with eigenvalue J + d . For Li decorated graphene, numerical calculation shows J + sy > J − sy , so g in the (+) states is lower than g in the ( − ) stateshence pairing in this modes are dominant. From Eqs. 19 and 20 we observe that probability amplitudes for pairing on differentbonds in real space differ for the various states. For the long C-C bonds the probability is proportional to ( α ± sy ) while for theothers is unity. Numerical results are shown inFig. 6. The possibility of a superconductivity state in metal decorated graphene has been suggested theoretically by a few groups.
Some have suggested phonon-mediated superconductivity in single layer graphene. Most prominently, Profeta et al. calcu-lated on the basis of density functional theory for superconductors that decoration by electron donating atoms such as Ca andLi will make single layer graphene superconducting, up to 8K for the case of Li. The ab initio anisotropic Migdal-Eliashbergformalism was used by Zheng and Margine , who predicted a single anisotropic superconducting gap with critical tempera-ture T c = . − . K , in surprisingly good agreement with experimental reported superconductivity around 6K in LiC . Using a phenomenological microscopic Hamiltonian in a nearest-neighbor tight-binding approximation, possible super-conducting phases of pristine graphene have been discussed by Uchoa and Castro-Neto and also by Black-Schaffer andDoniach . The possibility of a singlet p + ip phase pairing near the Dirac points between nearest neighbors subsites weresuggested by Uchoa and Castro-Neto . They worked in terms of a plasmon mediated mechanism for metal coated graphene,and discussed the conditions under which attractive electron-electron interaction can be mediated by plasmons. inglet superconducting gap phases of pristine graphene have been proposed and discussed by Black-Schaffer and Do-niach. For the nearest neighbors pairing amplitudes ∆ < iA jB > = ∆ iA , iA + ~ δ j where ~ δ j are the vectors that connects the iA siteto its three nearest neighbors, it was observed that there are three states that minimize the free energy in various regimes ofthe parameters, which here have been denoted by V s = ( , , ) T , V d x − y = ( , − , − ) T , and V d xy = ( , − , ) T . Pairing sym-metries d xy and d x − y are degenerate, and only the linear combination of d x − y + id xy ≡ d + id preserves the graphene bandsymmetry. Depending on the position of the Fermi energy with respect to Dirac points, d + id or s states tend to dominate.Their numerical calculation showed that d -wave solutions will always be favored for electron or hole doping in the regime0 < ¯ n c < . n α = < ˆ c † i α ˆ c i α > −
1. In this regime, superconductivity can emerge from electroniccorrelation effects. Near the van Hove singularity at the saddle point M corresponding to 3 / / n c = . d + id superconductivity, which breaks time-reversal symmetry, can be stabilized. In this regime d -wave superconductivity may arise from repulsive electron-electron interaction .Although doping by a gate voltage is normally considered to change only the chemical potential but not the band structure,gating cannot be expected to push the Fermi energy to the van Hove singularity without altering the band dispersion. The mostlikely way to do this is by decoration with electropositive atoms, which has been our focus. We note that doping is essential,when graphene decorated, in addition to the expected charge migration from the decorating atoms to the graphene sheet, itis then necessary the interlayer state is partially occupied to induce superconductivity as happens in GICs. Hybridization ofinterlayer s -band and graphene π bands changes the graphene band structure. The s orbitals of Ca have more overlap with Corbitals than Li and lead to stronger and longer range interactions as well as increasing the doping level, effects that becomedetrimental to superconductivity. For this reason our emphasis here is on the Li decorated graphene.We review some of our main points. When graphene is decorated by Li, electron transfer from Li atoms to C contracts theLi-C distance and reduces the C-C bond lengths in the Li-centered hexagon. In this kekulé -type structure, hopping amplitudesymmetries of all C-C neighbors are broken (our “ shrunken graphene”). This model allows study of multiband effects onthe superconducting phase diagram. To gain insight into our model, solutions of superconducting gap equation in both casesof folded bands otherwise pristine C and the usual two band model of C were compared. These two viewpoints coincideif the same pairing paradigms are considered. For pristine graphene with its two site cell, in real space picture electronscan pair with near neighbors in three inequivalent directions, ∆ i , i + ~ δ = V sy = ( ∆ ∆ ∆ ) T which must respect honeycombsymmetries. The V sy quantities are the three vectors that belong to the irreducible representation of crystal point group D h i.e. V Tsy = (1,1,1), (-1,1,0) and (2,-1,-1) for which the sy subscript stands for symmetries s , d xy and d x − y . Permutation of s -wavesolution (1,1,1) along three different bonds constructs just one state while permutation of d xy solution (-1,1,0) up to a minussign constructs two nonorthogonal linear independent states viz. (-1,1,0) and (-1,0,1) which orthogonal linear combination ofthem are d Txy = (-1, 1, 0) and d Tx − y = (2 -1 -1).A similar procedure again can be applied to pristine graphene but now in enlarged six site unit cell. Unit cell of C includessix carbon subsites and nine different bonds that support nine possible nearest neighbor bond pairing amplitudes as illustratedin Fig. 4 and denoted them by Φ Tsy = [( ∆ ′′ , ∆ ′′ ∆ ′′ )( ∆ , ∆ , ∆ )( ∆ ′ , ∆ ′ ∆ ′ )] . The gap equation is a 9 × ( , , , , , , , , ) the other eight solutions are constructed by all possible permutations of (-1,1,0) along these bonds thatpreserve our supercell symmetry. There are only three solutions which can preserve symmetry of both two and six atoms cellssimultaneously which they are of the form Φ + sy = ( V sy V sy V sy ) T as illustrated in Fig. 4. For these solutions, the folded 9 × × Φ f , Φ p x and Φ p y , of the form Φ sy = ( V sy − V sy ) T . These phases have been designated as islandphases, as illustrated in Fig. 5(b) for Φ f , within which a pairing amplitude is localized within island hexagons and cannotpropagate. For these island phases, numerical calculation of the electron pair potential energy g shows that g is large. Thiskind of solutions is a consequence of the six atom basis and does not appear for the two atom basis. Also, there are threesolutions of the form Φ − sy = ( − V sy V sy V sy ) which also break symmetry of two atom cell. For these reasons, in associationwith the normal state band structure of graphene, we concentrate on superconductivity in the three Φ + sy symmetry phases.For pristine graphene C , two normal bands are E ± = ± t | η | which fold to six branches in mini-BZ of C i.e. E ± γ = ± t | η | , E ± β = ± t | η | and E ± α = ± t | η | as shown in figure:eta-k, also Bloch-wave symmetry character of each branch hasbeen distinguished. The Bloch coefficients of the branch labeled by γ are of s -wave character, C A i = ( , , ) and for thoselabeled as α and β are of the form d ± id type, i.e. C A i = ( , e ± i i π , e ± i i π ) . Based on Bloch wave character of these branchesone can obtain the dominant superconducting phases of pristine graphene in various doping regimes. d -wave pairing emergesfrom the d -wave branches of the folded band structure E ± α and E ± β , while s -wave pairing arises from the s -wave branch E ± γ . or folded but otherwise pristine graphene, Fig. 2 illustrates that the lowest conduction band, weakly dispersive along Γ → M ,is responsible for dominant singlet superconductivity in chiral d ± id symmetry. Upon electron doping to the critical vHs at¯ n c = .
25, the pairing potential g in the d ± id phase decreases, beyond which density of states decreases. g increases until asecond critical value of doping ¯ n c = . s -wave pairing occurs. Bloch states in higher conductionbands include combinations of s and f symmetries that favor extended s wave pairing. The multiband character is responsiblefor stabilizing singlet s superconductivity at high electron or hole doping.To understand how superconducting phases of graphene can be affected by decoration by Li, one can compare the LiC gapsolutions with those of folded bands C at the same doping. Numerical results for pristine graphene gap equation performedin the nearest neighbor approximation in Ref. have been extended by applying a more accurate tight binding model fit tothe DFT band structure of pristine graphene . Although a quantum critical point for zero doping reported by Black-Schafferand Doniach at dimensionless coupling g t = .
91 which d - and s -wave solutions are degenerate. In the more realistic tightbinding model we applied, this degeneracy is not observed at the Γ point, and the d -wave solution is dominant. This differencemay be consequence of particle-hole symmetry breaking of valence and conduction bands. Also the van Hove singularity at theM point is moved from 0.25 doping for nearest neighbor hopping to 0.16 doping in the accurate model. The phase transitionfrom d -wave to s -wave is shifted to 0.35 doping instead of the 0.4 doping reported for nearest neighbor hopping. . Numericalcalculations for this more detailed model are illustrated in Fig. 7.When graphene is decorated by Li, around 0.68 electron per lithium atom transfers to neighboring C sites, viz. ¯ n c = . Γ move to -1.52 eV. Symmetry breaking of the hopping partially removes degeneracies of bandstructure of pristine graphene, which leads to creation of the small gap at Γ , with energy E g = | t − t ′ | = . eV . Also twoof four-fold degeneracies between valence and conduction bands at the Dirac points are removed. Compression between bandstructure of decorated graphene and folded pristine graphene at the same doping shows that hybridization of the Li s bandand C π band is small. This means nearest neighbor Li-C hopping is in the range t LiC ∼ . − .
5, and further hoppings arenegligible.Li decoration of graphene changes not only the band structure but also the Bloch wave coefficients from those of pristinegraphene. While pristine graphene Bloch wave coefficients have pure s - or d -wave character and their magnitudes are ~ k -independent. In the case of LiC they become mixed and vary with ~ k , hence gap equation symmetry is reduced. Becauseof this symmetry reduction, for the longer C-C bonds, a new coefficient α sy appears in the pairing amplitudes. In terms ofthis coefficient we have classified superconducting phase symmetries into three groups. Eqs. 18, 19, and 20 present all ninepossible pairing phases of LiC . There are three categories of solutions which have not appeared in complete form in theliterature. The total of nine phases arise from spatial, and therefore hopping parameter, symmetry breaking.In the first category Φ f , Φ p x and Φ p y , there is α sy = . For the second category, α sy (denoted by α − ) is negative, in the case of pristine α − = − g must be large to realize these phases. For the lastcategory α + is positive. Three phases which correspond to α + > Φ + d x − y , Φ + d xy , and Φ + s , and these have the lowestpairing potentials with respect to the other six phases.In the limiting case of folded six band pristine graphene α + d x − y , α + d xy , and α + s are all equal to unity, which maps theresults to the two-band symmetries as it should. But when Li decorated, depending on doping strength viz. w t and t LiC thesecoefficients α + sy no longer remain unity. The pairing amplitude distortion along longer C-C bonds α + , for s-wave phase issignificant due to its spatial isotropic symmetry. In spite of the pristine nature this phase no longer preserves two band modelsymmetry. On the other hand, d -wave phases are hardly affected by doping and their superconductivity is more persistentagainst perturbation. The chirality or non-chirality of Cooper pairs in these phases is undetermined, however. As shown inFig. 6(b), at low temperature α + s ≈ . Φ + s , and α + d x − y = α + d xy ≡ α + d is approximately equal to unity and varies little withtemperature.At a given critical temperature T c and chemical potential µ , for each of nine possible superconducting phases, Eqs. 10, 13,and 17 were evaluated numerically over the BZ of Li C to find the corresponding pairing potential g = J sy and α sy coefficient.Smaller g means less Cooper pair formation energy is required. Fig. 6(a) provides the phase boundaries for T c in terms of thepairing potential g for Li C in which µ =
0. For a given transition temperature T c , by changing the chemical potential µ ofLi C via gating, one can engineer the pairing potential g . Figure 8 gives a g - µ phase boundary diagram at T c = . d to “ distorted s -wave.” Changing µ o up to µ o − v ≈ . eV sothat the distance between the Fermi energy and the saddle points decreases, leads to a decrease in g . Continuously increasing µ o up to 0.5 eV causes g to increase for both d -wave and “distorted s -wave” pairing, and after that a smooth decrease proceeds.For both symmetries at critical µ o − c = . µ o − c = . unction of this band consists of d and p character, therefore Γ , Γ , Γ and Γ in Eq. 13 carry minus signs. This makes itevident from Eq. 17 that d − wave pairing is dominant. Beyond that, the uneven part of the “flat band” and also upper bandsassume a major role. These bands consist of d , p , s , and f character Bloch wave functions (as defined in earlier sections) witha significantly low density of states. In this case Γ , Γ , Γ and Γ change their sign, hence s -wave pairing is favored.Numerically we have demonstrated that electron pairing g in the limit of pristine graphene is minimal for all dopings. Ourcalculations indicate that any perturbation of the flat band reduces T c . The flat band can be perturbed through electron hoppingfrom decorating atoms to carbon sites ( t LiC ) or by hopping symmetry breaking index w t . For fixed doping at ¯ n = .
11 electronper carbon site and for fixed w t = .
94 as obtained for lithium decorated, in a variety of Li-C hopping between 0.3-0.4eV,numerical calculation doesn’t show significant altering of pair interaction potential g in s - and d -wave phases. But , as onecould see there is not an explicit behavior in a general coupling strength. A result is that a general aspect of superconductingpairing in LiC and pristine graphene is almost the same in the d x − y and d xy phases due to robustness of the flat band againstperturbation.To summarize, our calculations indicate that d -wave phases exist and are dominant symmetry of pairing in both pristineand Li decorated graphene. Pure s -wave phase does not appear in LiC , and s -wave superconductivity in metal decoratedgraphene is disfavored because of spatially increased overlap for s -wave symmetry. These results show that while degree ofdoping plays a major role in the graphene superconductivity, perturbation effects of decorating atoms finally determine thephase diagram. Our work also provides a new type of classification of superconducting phases in Li C -like nanostructures,and certain aspects of the formalism may be useful in modeling the recently observed superconductivity in magic angle bilayergraphene. References Wang, Q. Y. et al.
Interface-Induced High-Temperature Superconductivity in Single Unit-Cell FeSe Films on SrTiO . Chin.Phys. Lett. , 037402 (2012). He, S. L. et al.
Phase diagram and electronic indication of high-temperature superconductivity at 65 K in single-layer FeSefilms. Nat. Mater. , 605 (2013). Cao, Y. et al.
Unconventional superconductivity in magic-angle graphene superlattices. Nature , 43-50 (2018). Weller, T. E. et al.
Superconductivity in the intercalated graphite compounds C Yb and C Ca. Nature Phys. , 39 (2005). Ludbrook, B.M. et al.
Evidence for superconductivity in Li-decorated monolayer graphene. Proc. Natl. Acad. Sci. USA , 11795-11799 (2015). Palinkas, A. et al.
Novel graphene/Sn and graphene/SnO x hybrid nanostructures: induced superconductivity and band gapsrevealed by scanning probe measurements. Carbon , 611 (2017). Tiwari, A. P. et al.
Superconductivity at 7.4K in few layer graphene by Li intercalation. J. Phys.: Condens. Matt. ,445701 (2017). Woo, S. et al.
Temperature-dependent transport properties of graphene deorated by alkali metal adatoms (Li,K). Appl.Phys. Lett. , 263502 (2017). Uchoa, B. & Castro Neto, A. Superconducting States of Pure and Doped Graphene. Phy. Rev. Lett. , 146801 (2007). Nandkishore, R., Levitov, L. S. & Chubukov, A. V. Chiral superconductivity from repulsive interactions in doped grapheneNat. Phys. , 158–163 (2012). Nandkishore, R., Thomale R. & Chubukov A. V. Superconductivity from weak repulsions in hexagonal lattice systems.Phys. Rev. B 89, 144501 (2014).
Black-Schaffer, A. M. & Doniach, S. Resonating valence bonds and mean field d-wave superconductivity in graphene.Phys. Rev. B , 134512 (2007). Kiesel M. L. , Platt C., Hanke W. , Abanin D. A., & Thomale R. Competing many-body instabilities and unconventionalsuperconductivity in graphene. Phys. Rev. B , 020507R (2012). Ma, T., Yang, F., Yao, H. & Lin H. Q. Possible triplet p + ip superconductivity in graphene at low filling. Phys. Rev. B ,245114 (2014). Profeta, G., Calandra, M. & Mauri, F. Phonon-mediated superconductivity in graphene by Lithium deposition. Nat. Phys. , 131-134 (2012). Wong, C. H., Lortz, R., Buntov, E. A., Kasimova, R. E. & Zatsepin, A. F. A theoretical quest for high temperaturesuperconductivity on the example of low-dimensional carbon structures. Sci. Rep. , 15805 (2017). Hou C. Y., Chamon, C., & Mudry C. Electron Fractionalization in Two-Dimensional Graphenelike Structures. Phys. Rev.Lett. , 186809 (2007). Wu, L. -H. & Hu, X. Topological Properties of Electrons in Honeycomb Lattice with Detuned Hopping Energy. Sci. Rep. , 24347 (2016). Guzman, D. M., Alyahyaei H. M. & Jishi, R. A. Superconductivity in graphene-lithium. 2D Materials , (2014) 021005. Wallace, P. R. The Band Theory of Graphite. Phys. Rev. , 622 (1947). Reich, S., Maultzsch, J., Thomsen, C. & Ordejón, P. Tight-binding description of graphene. Phys. Rev. B , 035412(2002). Kundu, R. Tight Binding Parameters for Graphene. Mod. Phys. Lett. B , 163 (2011). Jung J. & MacDonald, A. H. Tight-binding model for graphene π -bands from maximally localized Wannier functions.Phys.Rev. B , 195450 (2013). Zheng, J. -J. & Margine, E. R. First-principles calculations of the superconducting properties in Li-decorated monolayergraphene within the anisotropic Migdal-Eliashberg formalism. Phys.Rev. B , 064509 (2016). Acknowledgments
R. Gholami acknowledges support that allowed an extended visit to the University of California Davis during part of this work.W.E.P. was supported by NSF grant DMR-1207622.
Author contributions
R. M. proposed the idea. R. G. and R. M. did the analytical derivation. R. G. carried out the numerical calculations undersupervision of R. M. and W. E. P. and S. M. performed the DFT calculations. All authors analyzed the results and wrote thearticle.
Additional information
Competing interests:
The authors declare no competing interests.
AppendicesA Accurate Tight Binding Model for Lithium decorated Graphene
The Hamiltonian of non-interacting LiC isˆ H N = − ∑ i α ∑ j β , σ t σ , σ i α , j β c † i ασ c j βσ + ∑ i α , σ ( ε i α − µ o ) ˆ n i ασ . (22)Eq 22 can be diagonalized in terms of Bloch eigenfunctions of the form (cid:12)(cid:12) Ψ ~ k (cid:11) = √ N N ∑ n = ∑ α = C α e i ~ k .~ r n α | φ n α i (23)in which ~ r n α = ~ r n + ~ d α and ~ r n is n th Bravais lattice site vector position and ~ d α is vector position of the α -th subsite withrespect to unit cell n . | φ n α i = (cid:12)(cid:12)(cid:12) φ n α ( ~ r − ~ r n − ~ d α ) E is the atomic π electron ket of atom α of cell n . By inserting Eqs. 22 and 23into ˆ H N ( ~ k ) (cid:12)(cid:12) Ψ ~ k (cid:11) = E ( ~ k ) (cid:12)(cid:12) Ψ ~ k (cid:11) , the equation for the coefficients becomes ∑ β = ε αβ ( ~ k ) C β + ( ε α − µ o ) C α = E ( ~ k ) C α where ε αβ ( ~ k ) = − N ∑ i j t σσ i α j β e i ~ k . ( ~ r i α − ~ r j β ) . (24)Eq. 24 can be expressed in following matrix form ε LiLi ( ~ k ) + ε Li − µ o h LiA ( ~ k ) h LiB ( ~ k ) h † LiA ( ~ k ) h AA ( ~ k ) + ε A − µ o h AB ( ~ k ) h † LiB ( ~ k ) h BA ( ~ k ) h BB ( ~ k ) + ε B − µ o C Li ( E i ( ~ k )) ξ A ( E i ( ~ k )) ξ B ( E i ( ~ k )) = E i ( ~ k ) C Li ( E i ( ~ k )) ξ A ( E i ( ~ k )) ξ B ( E i ( ~ k )) (25) here ε LiLi ( ~ k ) = t LiLi (cid:16) cos ~ k .~ ξ + cos ~ k .~ ξ + cos ~ k .~ ξ (cid:17) + t LiLi (cid:16) cos ~ k . ( ~ ξ − ~ ξ ) + cos ~ k . ( ~ ξ − ~ ξ ) + cos ~ k . ( ~ ξ − ~ ξ ) (cid:17) + t LiLi (cid:16) cos 2 ~ k .~ ξ + cos2 ~ k .~ ξ + cos2 ~ k .~ ξ (cid:17) + . . . (26)In Eq. 25, h -sub-block matrices are C-C or Li-C dispersion matrices and we have h AA = h ∗ BB , h AB = h † BA . The carbon-carbondispersion matrices i.e. ε ( A , B ) i ( A , B ) j ( ~ k ) are h AA ( ~ k ) = α ( ~ k ) β ( ~ k ) γ ( ~ k ) β ∗ ( ~ k ) α ( ~ k ) θ ( ~ k ) γ ∗ ( ~ k ) θ ∗ ( ~ k ) α ( ~ k ) , h AB ( ~ k ) = τ ( ~ k ) d ( ~ k ) d ( ~ k ) d ( ~ k ) τ ( ~ k ) d ( ~ k ) d ( ~ k ) d ( ~ k ) τ ( ~ k ) . (27)The Li-C dispersion row matrices i.e. ε LiA i ( ~ k ) and ε LiB i ( ~ k ) are h LiA ( ~ k ) = [ d c ( ~ k ) d c ( ~ k ) d c ( ~ k )] = − t LiC e ik z h (cid:16) e i ~ k .~ δ e i ~ k .~ δ e i ~ k .~ δ (cid:17) , h LiB ( ~ k ) = e ik z h h LiA ∗ ( ~ k ) (28)where e k z h factor takes 1 by confinement. New variables, ~ k dependent on-site energy and chemical potential has been definedas ε ( ~ k ) = ε Li − µ o + ε LiLi ( ~ k ) ε ( ~ k ) = ε A − µ o + ε A A ( ~ k ) ε ( ~ k ) = ε B − µ o + ε B B ( ~ k ) (29)Shorthand notation has been introduced as follows, α ( ~ k ) ≡ ε A i A i ( ~ k ) = ε B i B i ( ~ k ) = − t − t (cid:16) cos ~ k .~ ξ + cos ~ k .~ ξ + cos ~ k .~ ξ (cid:17) β ( ~ k ) ≡ ε A A ( ~ k ) = ε ∗ A A ( ~ k ) = − t e i ~ k . ( ~ δ − ~ δ ) h + w t (cid:16) e − i ~ k .~ ξ + e i ~ k .~ ξ (cid:17)i γ ( ~ k ) ≡ ε A A ( ~ k ) = ε ∗ A A ( ~ k ) = − t e i ~ k . ( ~ δ − ~ δ ) h + w t (cid:16) e − i ~ k .~ ξ + e i ~ k .~ ξ (cid:17)i θ ( ~ k ) ≡ ε A A ( ~ k ) = ε ∗ A A ( ~ k ) = − t e i ~ k . ( ~ δ − ~ δ ) h + w t (cid:16) e − i ~ k .~ ξ + e i ~ k .~ ξ (cid:17)i (30)in which w t = t ′ t = t ′ t and ~ ξ i = ~ τ i + ~ δ i the d and τ ,functions are given by τ ( ~ k ) = − t ′ e i ~ k .~ τ " + t t ′ e − i ~ k .~ ξ + t t ′ (cid:16) e i ~ k .~ ξ + e i ~ k .~ ξ (cid:17) , d ( ~ k ) = − t e i ~ k .~ δ (cid:20) + t t e − i ~ k .~ ξ + t t (cid:16) e i ~ k .~ ξ + e i ~ k .~ ξ (cid:17)(cid:21) τ ( ~ k ) = − t ′ e i ~ k .~ τ " + t t ′ e − i ~ k .~ ξ + t t ′ (cid:16) e i ~ k .~ ξ + e i ~ k .~ ξ (cid:17) , d ( ~ k ) = − t e i ~ k .~ δ (cid:20) + t t e − i ~ k .~ ξ + t t (cid:16) e i ~ k .~ ξ + e i ~ k .~ ξ (cid:17)(cid:21) τ ( ~ k ) = − t ′ e i ~ k .~ τ " + t t ′ e − i ~ k .~ ξ + t t ′ (cid:16) e i ~ k .~ ξ + e i ~ k .~ ξ (cid:17) , d ( ~ k ) = − t e i ~ k .~ δ (cid:20) + t t e − i ~ k .~ ξ + t t (cid:16) e i ~ k .~ ξ + e i ~ k .~ ξ (cid:17)(cid:21) . (31)Also ξ A ( ~ k ) = C A ( E i ( ~ k )) C A ( E i ( ~ k )) C A ( E i ( ~ k )) ; ξ B ( ~ k ) = C B ( E i ( ~ k )) C B ( E i ( ~ k )) C B ( E i ( ~ k )) . (32)Since the Li is an inversion center of the two-dimension tight binding model, the Bloch wave function should respect inversionsymmetry. Using Ψ ~ k ( ~ r ) = h ~ r | Ψ ~ k i , the condition is Ψ − ~ k ( − ~ r ) = Ce i φ Ψ ~ k ( ~ r ) , (33)where C is ± ε A i = ε B i . This condition is satisfied if C A i is proportional to C ∗ B i C A m ( E i ( ~ k )) = f i ( ~ k ) C ∗ B m ( E i ( ~ k )) , and C Li ( E i ( ~ k )) = f i ( ~ k ) C ∗ Li ( E i ( ~ k )) , i = , , ..., here f i ( ~ k ) is a coefficient to be determined. By inserting C A i = | C A i | e i φ Ai and C B i = | C B i | e i φ Bi into Eq. 34 we have f i ( ~ k ) = | C A m ( E i ( ~ k )) || C B m ( E i ( ~ k )) | e i ( φ Ai + φ Bi ) = C i e i φ i , i = , , ..., . (35)In the next subsection we use Eqs. 34 and 35 to reduce the eigenvalue Eq. 24 in matrix form to 3 × B Uncoupled C Dispersion Relations
By first neglecting the lithium-carbon hopping t LiC →
0, the shrunken graphene Hamiltonian Eq. 25 can be diagonalizedexactly. Our notation is ε ( ~ k ) ε ( ~ k ) β ( ~ k ) γ ( ~ k ) τ ( ~ k ) d ( ~ k ) d ( ~ k ) β ∗ ( ~ k ) ε ( ~ k ) θ ( ~ k ) d ( ~ k ) τ ( ~ k ) d ( ~ k ) γ ∗ ( ~ k ) θ ∗ ( ~ k ) ε ( ~ k ) d ( ~ k ) d ( ~ k ) τ ( ~ k ) τ ∗ ( ~ k ) d ∗ ( ~ k ) d ∗ ( ~ k ) ε ( ~ k ) β ∗ ( ~ k ) γ ∗ ( ~ k ) d ∗ ( ~ k ) τ ∗ ( ~ k ) d ∗ ( ~ k ) β ( ~ k ) ε ( ~ k ) θ ∗ ( ~ k ) d ∗ ( ~ k ) d ∗ ( ~ k ) τ ∗ ( ~ k ) γ ( ~ k ) θ ( ~ k ) ε ( ~ k ) C Li ( E i ( ~ k )) C A ( E i ( ~ k )) C A ( E i ( ~ k )) C A ( E i ( ~ k )) C B ( E i ( ~ k )) C B ( E i ( ~ k )) C B ( E i ( ~ k )) = E i ( ~ k ) C Li ( E i ( ~ k )) C A ( E i ( ~ k )) C A ( E i ( ~ k )) C A ( E i ( ~ k )) C B ( E i ( ~ k )) C B ( E i ( ~ k )) C B ( E i ( ~ k )) (36)The non trivial eigenvalues of uncoupled Hamiltonian Eq. 36 are given by E sh ; n ( ~ k ) = E sh , ml ( t i ,~ ξ i ,~ k ) = − µ o + α ( ~ k ) + u m Π ( ~ k ) + u ∗ m Π ∗ ( ~ k ) + (cid:20) ε A + ε B + ( − ) l q ( ε A − ε B ) + w m ( ~ k ) (cid:21) (37)where n is band index defined as n = ( l + ) + ( − ) l m ; m = (cid:26) , , − , − , − E sh , , E sh , and E sh , are conduction bands which corresponds to l = m = , , E sh , , E sh , and E sh , arevalence bands which correspond to l = , m = − , − , − Γ point, when ε A = ε B = ε c , the shrunk graphene eigenstates | φ n ( ) i = ( C A C A C A C B C B C B ) T over thehexagonal subsites take the forms similar to conventional s , d and p orbitals, | f i = ( − − − ) T , | p x i = ( − − ) T , (cid:12)(cid:12) p y (cid:11) = ( − − − ) T , (cid:12)(cid:12)(cid:12) d x − y E = ( − − ) T , (cid:12)(cid:12) d xy (cid:11) = ( − − ) T , | S i = ( ) T , (39)with energies, E f = E + γ ( ) = µ s + [( t ′ + t ) + t + t ] E s = E − γ ( ) = µ s − [( t ′ + t ) + t + t ] E p = E + α ( ) = E + β ( ) = µ d + ( t ′ − t ) E d = E − α ( ) = E − β ( ) = µ d − ( t ′ − t ) (40)where µ s = ε c − µ − ( t + t ′ ) − t and µ d = ε c − µ + ( t + t ′ ) − t . For general ~ k the uncoupled shrunken grapheneeigenfunction Eq. 59 can be written in terms of | S i , | f i , | p x i , (cid:12)(cid:12) p y (cid:11) , (cid:12)(cid:12) d xy (cid:11) and (cid:12)(cid:12)(cid:12) d x − y E as (cid:12)(cid:12)(cid:12) φ n ( ~ k ) E = (cid:16) f ns ( ~ k ) | S i + f nf ( ~ k ) | f i (cid:17) + (cid:16) f np y ( ~ k ) (cid:12)(cid:12) p y (cid:11) + i f nd xy ( ~ k ) (cid:12)(cid:12) d xy (cid:11)(cid:17) + (cid:16) f nd x − y ( ~ k ) (cid:12)(cid:12)(cid:12) d x − y E + i f np x ( ~ k ) | p x i (cid:17) . (41)In the particular case of pristine graphene in which w t = ~ τ = ~ δ , ~ τ = ~ δ , ~ τ = ~ δ , hence β = θ = γ ∗ and also ε A = ε B so C = ±
1. Therefore eigenvectors take following form (cid:12)(cid:12)(cid:12) ϕ m , l ( ~ k ) E = √ (cid:18) ω m ω m ∗ ( − ) l η ∗ m | η m | ω m ∗ ( − ) l η ∗ m | η m | ω m ( − ) l η ∗ m | η m | (cid:19) T , ω m = e i π m / , (42) here m = , , l = , η m ( ~ k ) = d ( ~ k ) + ω m d ( ~ k ) + ω ∗ m d ( ~ k ) . The eigenvalues are E m , l = ε A A ( ~ k ) + ω m β ( ~ k ) + ω ∗ m β ∗ ( ~ k ) + ( − ) l t | η m ( ~ k ) | (43)By comparing pristine eigenvectors Eq. 42 with general shrunken graphene eigenvectors Eq. 41 it is found that for m = f d x − y ( ~ k ) = − f p x ( ~ k ) = ( − e i φ ) √ which corresponds to d x − y − ip x , also f p y ( ~ k ) = − f d xy ( ~ k ) = ( + e i φ ) √ which corresponds to p y − id xy , and the coefficients of s and f are zero. For m = f d x − y ( ~ k ) = f p x ( ~ k ) = ( − e i φ ) √ which corresponds to d x − y + ip x also f p y ( ~ k ) = f d xy ( ~ k ) = ( + e i φ ) √ which corresponds to p y + id xy while s and f coefficients are zero. For m = f s ( ~ k ) = ( − e i φ ) and f f ( ~ k ) = ( + e i φ ) and other coefficients are zero. Here, e i φ m = η ∗ m | η m | .The Hamiltonian ˆ H shrN for the broken symmetry (shrunk graphene) is 6 × × (cid:18) h AA ( ~ k ) + ε A − µ o h AB ( ~ k ) h BA ( ~ k ) h ∗ AA ( ~ k ) + ε B − µ o (cid:19) (cid:18) ξ A ( E sh ; i ( ~ k )) ξ B ( E sh ; i ( ~ k )) (cid:19) = E sh ; i ( ~ k ) (cid:18) ξ A ( E sh ; i ( ~ k )) ξ B ( E sh ; i ( ~ k )) (cid:19) (44)To solve Schrödinger Eq. 44 we first separate the left hand side of Eq. 44 into two terms H sh ( ~ k ) = H sh ( ~ k ) + H sh ( ~ k ) where H sh ( ~ k ) = (cid:18) ε ( ~ k ) I × h AB ( ~ k ) h BA ( ~ k ) ε ( ~ k ) I × (cid:19) , H sh ( ~ k ) = (cid:18) h ′ aa ( ~ k ) × × h ′ bb ( ~ k ) (cid:19) (45)where h ′ aa ( ~ k ) = ( h AA ( ~ k ) − α ( ~ k ) I × ) . In the nearest neighbor approximation it is straightforward to show that H sh ( ~ k ) and H sh ( ~ k ) commute with each other, hence they have the same eigenvectors. In more generality one can consider approximatelythe eigenvectors of H sh to be the same as the eigenvectors of H sh , therefore one obtains ξ B ( ~ k ) = Ce i φ ξ ∗ A : E sh ( ~ k ) ≈ E sh ( ~ k ) + E sh ( ~ k ) . (46)with the first equation arising from similar conditions as for Eq. 34. To find E sh ( ~ k ) we solve the following eigenvalue problems (cid:18) h ′ aa ( ~ k ) × × h ′ ∗ aa ( ~ k ) (cid:19) (cid:18) ξ A ( ~ k ) ξ B ( ~ k ) (cid:19) = E sh ( ~ k ) (cid:18) ξ A ( ~ k ) ξ B ( ~ k ) (cid:19) (47)and (cid:18) ε ( ~ k ) I × h AB ( ~ k ) h BA ( ~ k ) ε ( ~ k ) I × (cid:19) (cid:18) ξ A ( ~ k ) ξ B ( ~ k ) (cid:19) = E sh ( ~ k ) (cid:18) ξ A ( ~ k ) ξ B ( ~ k ) (cid:19) (48)Eq. 47 converts to following eigenvalue problem h ′ aa ( ~ k ) ξ A ( ~ k ) = E sh ( ~ k ) ξ A ( ~ k ) and its complex conjugate with A ↔ B , defining c ( t ,~ ξ i ,~ k ) = β ( ~ k ) θ ( ~ k ) γ ∗ ( ~ k ) + γ ( ~ k ) β ∗ ( ~ k ) θ ∗ ( ~ k ) , c ( t ,~ ξ i ,~ k ) = (cid:12)(cid:12)(cid:12) β ( ~ k ) (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) θ ( ~ k ) (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) γ ( ~ k ) (cid:12)(cid:12)(cid:12) Π ( t ,~ ξ i ,~ k ) = c ( t ,~ ξ i ,~ k ) + i s(cid:18) c ( t ,~ ξ i ,~ k ) (cid:19) − (cid:18) c ( t ,~ ξ i ,~ k ) (cid:19) / . (49)eigenvalue equation Eq. 47 has the three different eigenvalues E sh , m ( t ,~ ξ i ,~ k ) = u m Π ( t ,~ ξ i ,~ k ) + u ∗ m Π ∗ ( t ,~ ξ i ,~ k ) ; u m = √ = e im π / ; m = , , E sh , m ( t ,~ ξ i ,~ k ) is function of second neighbor hopping t < iA i , jA j > = t and ~ ξ i = ~ τ i + ~ δ i i.e. lattice bases vector and itdoes not depend on ~ τ i and ~ δ i separately. Now we calculate E sh ( ~ k ) . Eq. 48 can be separated into following equations h AB ( ~ k ) ξ B ( ~ k ) = ( E sh ( ~ k ) − ε ( ~ k )) ξ A ( ~ k ) , h BA ( ~ k ) ξ A ( ~ k ) = ( E sh ( ~ k ) − ε ( ~ k )) ξ B ( ~ k ) . (51)By multiplying first equation of Eq. 51 by h BA and second by h AB we have h AB ( ~ k ) h BA ( ~ k ) ξ A ( ~ k ) = ( E sh ( ~ k ) − ε ( ~ k ))( E sh ( ~ k ) − ε ( ~ k )) ξ A ( ~ k ) h BA ( ~ k ) h AB ( ~ k ) ξ B ( ~ k ) = ( E sh ( ~ k ) − ε ( ~ k ))( E sh ( ~ k ) − ε ( ~ k )) ξ B ( ~ k ) . (52) q. 52 is an eigenvalue problem where second equation is just complex conjugated of first one. By defining new matrix G ( ~ k ) = h AB ( ~ k ) h BA ( ~ k ) and w i ( ~ k ) = [ E sh , i ( ~ k )] − [ ε ( ~ k ) + ε ( ~ k )] E sh , i ( ~ k ) + ε ( ~ k ) ε ( ~ k ) , Eq. 52 takes the form G ( ~ k ) ξ A ( E i ( ~ k )) = w i ( ~ k ) ξ A ( E i ( ~ k )) (53)Schrödinger Eq. 53 can be solved to find eigenvalues of Eq. 48 i.e. E sh ( ~ k ) . Defining C ( t i ,~ ξ i ,~ k ) = G + G + G C ( t i ,~ ξ i ,~ k ) = | G | + | G | + | G | − ( G G + G G + G G ) C ( t i ,~ ξ i ,~ k ) = G ( G G ) ∗ + G ∗ ( G G ) − G | G | − G | G | − G | G | + G G G (54)where G i j = ∑ m = ε A i B m ( ~ k ) ε B m A j ( ~ k ) . Also introducing Π ( t i ,~ ξ i ,~ k ) = (cid:18) Q ( t i ,~ ξ i ,~ k ) + i q P ( t i ,~ ξ i ,~ k ) − Q ( t i ,~ ξ i ,~ k ) (cid:19) Q ( t i ,~ ξ i ,~ k ) = C ( t i ,~ ξ i ,~ k ) + C ( t i ,~ ξ i ,~ k ) C ( t i ,~ ξ i ,~ k ) + C ( t i ,~ ξ i ,~ k ) P ( t i ,~ ξ i ,~ k ) = C ( t i ,~ ξ i ,~ k ) + C ( t i ,~ ξ i ,~ k ) . (55)where t i = t iA jB are first, 3rd and 4th neighbor hopping integrals. Hence eigenvalues of Eq. 48 can be obtained, they are E sh ; m , l ( t i ,~ ξ i ,~ k ) = ε A A ( ~ k ) − µ o + (cid:20) ε A + ε B + ( − ) l q ( ε A − ε B ) + w m ( t i ,~ ξ i ,~ k ) (cid:21) (56)where l = , w m ( t i ,~ ξ i ,~ k ) i.e. solutions of Schrödinger Eq. 53 are w m ( t i ,~ ξ i ,~ k ) = C ( t i ,~ ξ i ,~ k ) + u m Π ( t i ,~ ξ i ,~ k ) + u ∗ m Π ∗ ( t i ,~ ξ i ,~ k ) ; u m = √ = e im π / ; m = , , . (57)Hence from Eqs. 46, 50 and 56 eigenvalues of Schrödinger Eq. 44 can be obtained, E sh ; m , l ( t i ,~ ξ i ,~ k ) = E sh , m ( t i ,~ ξ i ,~ k ) + ε A A ( ~ k ) − µ o + (cid:20) ε A + ε B + ( − ) l q ( ε A − ε B ) + w m ( t i ,~ ξ i ,~ k ) (cid:21) (58)The corresponding orthogonal eigenvectors are (cid:12)(cid:12)(cid:12) φ ( E sh , n ( ~ k )) E = C A ( E sh ; n ( ~ k )) (cid:20) (cid:18) C A ( E sh ; n ( ~ k )) C A ( E sh ; n ( ~ k )) C A ( E sh ; n ( ~ k )) C A ( E sh ; n ( ~ k )) (cid:19) C η ∗ ( E sh ; n ( ~ k )) | η ( E sh ; n ( ~ k )) | (cid:18) C ∗ A ( E sh ; n ( ~ k )) C ∗ A ( E sh ; n ( ~ k )) C ∗ A ( E sh ; n ( ~ k )) C ∗ A ( E sh ; n ( ~ k )) (cid:19) (cid:21) T (59)where C A ( E sh , n ( ~ k )) can be found from orthogonality condition. Also, C A ( E sh , i ( ~ k )) = − (cid:16) G − w i ( ~ k ) (cid:17) G + G G (cid:16) G − w i ( ~ k ) (cid:17) (cid:16) G − w i ( ~ k ) (cid:17) − | G | C A ( E sh , i ( ~ k )) C A ( E sh , i ( ~ k )) = − (cid:16) G − w i ( ~ k ) (cid:17) G + G G (cid:16) G − w i ( ~ k ) (cid:17) (cid:16) G − w i ( ~ k ) (cid:17) − | G | C A ( E sh , i ( ~ k )) . (60)To find ξ B ( E n ( ~ k )) it has been used symmetry condition in Eq. 46 C B m ( E sh ; i ( ~ k )) = Ce i ϕ ( E sh ; i ( ~ k )) C ∗ A m ( E sh ; i ( ~ k )) . (61)Replacing Eq. 61 into second equation of Eq. 51 we get e i ϕ ( E sh ; i ( ~ k )) = η ∗ ( E sh ; i ( ~ k )) (cid:12)(cid:12)(cid:12) η ( E sh ; i ( ~ k )) (cid:12)(cid:12)(cid:12) C A ( E sh ; i ( ~ k )) C ∗ A ( E sh ; i ( ~ k )) ; C = E sh , i ( ~ k ) − ε ( ~ k ) E sh , i ( ~ k ) − ε ( ~ k ) η ( E sh ; i ( ~ k )) = − t d C ∗ A ( E sh ; i ( ~ k )) C ∗ A ( E sh ; i ( ~ k )) + d C ∗ A ( E sh ; i ( ~ k )) C ∗ A ( E sh ; i ( ~ k )) + τ ! . (62)It is easy to show that (cid:12)(cid:12)(cid:12) η ( E sh ; i ( ~ k )) (cid:12)(cid:12)(cid:12) = w i ( ~ k ) Coupled Li-C dispersion relations By applying the following unitary transformation, P †0 N ˆ H N P N (cid:16) P †0 N (cid:12)(cid:12)(cid:12) Ψ ~ k , n ( ~ r ) E(cid:17) = E n ( ~ k ) P †0 N (cid:12)(cid:12)(cid:12) Ψ ~ k , n ( ~ r ) E , where ˆ P N is the operatorthat diagonalize Eq. 36, the Schrödinger Eq. 25 is written in a new matrix representation as E Li , ( ~ k ) γ ( ~ k ) γ ( ~ k ) γ ( ~ k ) γ ( ~ k ) γ ( ~ k ) γ ( ~ k ) γ ∗ ( ~ k ) E sh ;1 ( ~ k ) γ ∗ ( ~ k ) E sh ;2 ( ~ k ) γ ∗ ( ~ k ) E sh ;3 ( ~ k ) γ ∗ ( ~ k ) E sh ;4 ( ~ k ) γ ∗ ( ~ k ) E sh ;5 ( ~ k ) γ ∗ ( ~ k ) E sh ;6 ( ~ k ) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) = E i ( ~ k ) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) A ( E i ( ~ k )) (63)where the relation between the column matrix eigenstate of Eq. 63, A ( E i ( ~ k )) , and the eigenstates of Eq. 25, C is A = P †0 N C , and A j ( E i ( ~ k )) = γ ∗ j E i ( ~ k ) − E sh , j ( ~ k ) A ( E i ( ~ k )) (64)in which A ( E i ( ~ k )) is determined from the normalization condition, and also γ i ( ~ k ) = − t LiC ∑ m = ( C A m ( E sh , i ( ~ k )) e i ~ k .~ ζ Am + C B m ( E sh , i ( ~ k )) e i ~ k .~ ζ Bm ) (65)where ~ ζ A m is a vector that connects Li to the A m carbon atom and ~ ζ B m is a vector which connects Li to the B m carbon atom.At the Γ point, γ ( ) = −√ t LiC and γ i ( ) = i = , ...,
5. These results show that just the isolated intercalant band, E Li , ( ) and the lowest valance band, E sh , ( ) , are mutually affected. The energies of these bands are, with E ( ) ≡ E + , E ( ) ≡ E − , E ± ( ) = (cid:0) E Li , ( ) + E sh , ( ) (cid:1) ± s (cid:2) E Li , ( ) − E sh , ( ) (cid:3) + ( t LiC ) (66)and other shrunk graphene bands Eq. 39 remain unchanged. This means the energy gap at Γ , E g ( ) , depends only on thenearest neighbor hopping difference rather than on t LiC . That is because the overlap between the Li s band and the valanceband of uncoupled shrunken graphene (which is linear combination of s and f , | φ ( ) i ) is significant while others are zero.For general ~ k vectors it is challenging to obtain an exact analytical expression for the full Hamiltonian in Eq. 25 and itwould not be transparent anyway. One can use perturbation theory to obtain useful results. For H ND = H D + H D we can usenon degenerate perturbation theory, obtaining H ND ( ~ k ) | ψ N ( E i ) i = E i ( ~ k ) | ψ N ( E i ) i (67)where expansion of E i ( ~ k ) in terms of perturbation parameter is E i ( ~ k ) = E i ( ~ k ) + E i ( ~ k ) + E i ( ~ k ) + .... From non degenerateperturbation we have E i ( ~ k ) = h φ i | H D | φ i i , E i ( ~ k ) = ∑ j = i (cid:12)(cid:12)(cid:10) φ i | H D | φ j (cid:11)(cid:12)(cid:12) E i ( ~ k ) − E j ( ~ k ) (68)where | φ i i is i th eigenstate of diagonal H D . Perturbed system eigenstates up to first order are | ψ N ( E i ( ~ k )) > = (cid:12)(cid:12)(cid:12) φ i ( E i ( ~ k ) E + ∑ j = i E i − E j (cid:12)(cid:12)(cid:12) φ j ( E j ( ~ k ) E . (69)Non degenerate perturbation theory can be used in Eq. 63 for completely filled or empty bands that are far from lithium band, E Li , ( ~ k ) , and also without overlap. Therefore, except E sh , and E sh , which are nearly degenerate with lithium band in someregions, non degenerate approximation can be used for other bands of H N . We denote Hamiltonian in Eq. 63 as H DN . This amiltonian can be separated to H ND = H D + H D where H D = E Li , ( ~ k ) γ ( ~ k ) γ ( ~ k ) E sh , ( ~ k ) γ ∗ ( ~ k ) E sh , ( ~ k ) γ ∗ ( ~ k ) E sh , ( ~ k ) E sh , ( ~ k ) E sh , ( ~ k )
00 0 0 0 0 0 E sh , ( ~ k ) . (70) H D = γ ( ~ k ) γ ( ~ k ) γ ( k ) γ ( ~ k ) γ ∗ ( ~ k ) γ ∗ ( ~ k ) γ ∗ ( ~ k ) γ ∗ ( ~ k ) . (71)introducing below coefficients c = E Li , ( ~ k ) + E sh , ( ~ k ) + E sh , ( ~ k ) c = − ( E Li , ( ~ k ) E sh , ( ~ k ) + E Li , ( ~ k ) E sh , ( ~ k ) + E sh , ( ~ k ) E sh , ( ~ k ) − | γ ( ~ k ) | − | γ ( ~ k ) | ) c = E Li , ( ~ k ) E sh , ( ~ k ) E sh , ( ~ k ) − E sh , ( ~ k ) | γ ( ~ k ) | − E sh , ( ~ k ) | γ ( ~ k ) | Π = ( q + i p p − q ) , q = c + c c + c , p = c + c H D are E ( ~ k ) = c + Π + Π ∗ , E ( ~ k ) = c + e i π / Π + e − i π / Π ∗ , E ( ~ k ) = c + e − i π / Π + e i π / Π ∗ (73)and corresponding eigenstates are (cid:12)(cid:12)(cid:12) φ ( E ( ~ k )) E = ( C p ( E ( ~ k )) C p ( E ( ~ k )) C p ( E ( ~ k )) ) T (cid:12)(cid:12)(cid:12) φ ( E ( ~ k )) E = ( C p ( E ( ~ k )) C p ( E ( ~ k )) C p ( E ( ~ k )) ) T (cid:12)(cid:12)(cid:12) φ ( E ( ~ k )) E = ( C p ( E ( ~ k )) C p ( E ( ~ k )) C p ( E ( ~ k )) ) T (74) t is easy to show that E =
0. So up to second order perturbation parameter the eigenenergies are E ( ~ k ) = E ( ~ k ) − | C p ( E ) | | γ | E sh , − E + | γ | E sh , − E + | γ | E sh , − E + | γ | E sh , − E ! E ( ~ k ) = E sh , ( ~ k ) + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E E ( ~ k ) = E ( ~ k ) − | C p ( E ) | | γ | E sh , − E + | γ | E sh , − E + | γ | E sh , − E + | γ | E sh , − E ! E ( ~ k ) = E ( ~ k ) − | C p ( E ) | | γ | E sh , − E + | γ | E sh , − E + | γ | E sh , − E + | γ | E sh , − E ! E ( ~ k ) = E sh , ( ~ k ) + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E E ( ~ k ) = E sh , ( ~ k ) + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E E ( ~ k ) = E sh , ( ~ k ) + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E + | C p ( E ) | | γ | E sh , − E (75) D Bogoliubov-de Gennes Transformation
The electron-electron interaction part of Hamiltonian, H P , in the mean field approximation becomesˆ H MFP = ∑ i ασ ∑ j βσ ′ ∆ σσ ′ i α j β ˆ c † i ασ ˆ c † j βσ ′ + h . c . + F = ∑ ~ k ασβσ ′ ∆ σσ ′ αβ ( ~ k ) ˆ c † ασ ( ~ k ) ˆ c † βσ ′ ( − ~ k ) + h . c . + F (76)in which ∆ σσ ′ i α j β = U σσ ′ i α , j β (cid:10) ˆ c i ασ ˆ c j βσ ′ (cid:11) is the matrix of order parameters in real space. Fourier transformation of real space orderparameters are given by ∆ σσ ′ αβ ( ~ k ) = N ∑ i j ∆ σσ ′ i α j β e i ~ k . ( ~ r i α − ~ r j β ) . (77)here Latin subscripts, α and β refers to A i or B i subsites. The interacting Hamiltonian in Nambu space isˆ H SU = ∑ ~ k ˆ Ψ † ( ~ k ) H SU ( ~ k ) ˆ Ψ ( ~ k ) (78)where ˆ Ψ † ( ~ k ) = ( c †0 ↑ ( ~ k ) c †1 ↑ ( ~ k ) , ..., c †6 ↑ ( ~ k ) c ↓ ( − ~ k ) c ↓ ( − ~ k ) , ..., c ↓ ( − ~ k )) where Li , A , A , A , B , B and B are labeled by0 , , , , , , H SU ( ~ k ) = (cid:18) H N ( ~ k ) H P ( ~ k ) H † P ( ~ k ) − H ∗ N ( − ~ k ) (cid:19) . (79)The coupling is given by H P ( ~ k ) = ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) ∆ ↑↓∗ A B ( ~ k ) (80) or singlet ∆ ↑↓ βα ( ~ k ) = ∆ ↑↓ αβ ( − ~ k ) = − ∆ ↓↑ αβ ( − ~ k ) = ∆ ↑↓ αβ ∗ ( ~ k ) . The order parameters according to (main text) Fig. 4 are ∆ ↑↓ A B ( ~ k ) = ∆ ′′ e i ~ k .~ τ , ∆ ↑↓ A B ( ~ k ) = ∆ e i ~ k .~ δ , ∆ ↑↓ A B ( ~ k ) = ∆ ′ e i ~ k .~ δ ; ∆ ↑↓ A B ( ~ k ) = ∆ ′ e i ~ k .~ δ , ∆ ↑↓ A B ( ~ k ) = ∆ ′′ e i ~ k .~ τ , ∆ ↑↓ A B ( ~ k ) = ∆ e i ~ k .~ δ ; ∆ ↑↓ A B ( ~ k ) = ∆ e i ~ k .~ δ , ∆ ↑↓ A B ( ~ k ) = ∆ ′ e i ~ k .~ δ , ∆ ↑↓ A B ( ~ k ) = ∆ ′′ e i ~ k .~ τ . (81) U < iA ↑ jB ↓ > = U < iA ↑ jB ↓ > = U < iA ↑ jB ↓ > g Quasiparticle energies are obtained by unitary transformation in the seven bandspace ˆ H SU = ∑ ~ k ˆ Ψ † ( ~ k ) Q h Q † H SU ( ~ k ) Q i Q † ˆ Ψ ( ~ k ) = ∑ ~ k Λ † ( ~ k ) H NSU ( ~ k ) Λ ( ~ k ) (82)where in matrix notation, H NSU ( ~ k ) = (cid:18) H ND ( k ) H PD ( ~ k ) H † PD ( ~ k ) − H ∗ ND ( − ~ k ) (cid:19) , Q = ˆ P N ( ~ k ) ⌢ ⌢ P ∗ N ( − ~ k ) ! . (83)ˆ P N ( ~ k ) is a 7 × [ ˆ P † N ] i , j = C i ( E j ( ~ k )) . H ND ( ~ k ) is the corresponding diagonal seven band Hamiltonian. Here ˆ P ∗ N ( − ~ k ) = ˆ P N ( ~ k ) .In the normal band space the matrix elements of the off-diagonal array are given by [ H DP ( ~ k )] i , j = ∆ i j ( ~ k ) = ∑ α = Ω α i j ( ~ k ) ∆ α (84)Using the fact that the gap is small, applying perturbation up to second order in the order parameter gives quasiparticle energies E Qm , s ( ~ k ) = s E m ( ~ k ) + ∑ i = (cid:12)(cid:12)(cid:12) ∆ mi ( ~ k ) (cid:12)(cid:12)(cid:12) E m ( ~ k ) + E i ( ~ k ) s = ± s = s = − E Superconducting States
By minimizing the quasiparticle free energy with respect to nearest neighbor order parameters the gap equation is obtained.The free energy is F = − β ∑ ~ k ∑ n = ln " ( E Qn k B T ) + F . (86)where F is the system condensation energy, F = − ∑ i ασ ∑ j βσ ′ ∆ σσ ′ i α j β G σσ ′ i α j β † = − N ∑ α = J α ( ∆ α ) (87)where J = J = J = g and J = J = J = J = J = J = g . The linearized gap equation, obtained by minimizing freeenergy of the system, is J β ∆ β = − N ∑ α = ∑ ~ k ∑ n = ∑ i = tanh ( E Qn k B T ) E n ( ~ k ) + E i ( ~ k ) (cid:16) Ω α ni ( ~ k ) Ω ∗ β ni ( ~ k ) + Ω β ni ( ~ k ) Ω ∗ α ni ( ~ k ) (cid:17) ∆ α ≡ − ∑ α = Γ βα ∆ α . (88)We have used that at T c , where | ∆ i j | can be neglected, E Qn → E n .In the general case, ( see main text Fig. 4) , the system is invariant under interchange 2 ⇋
5, 3 ⇋ ⇋
7, whichmeans ∆ ′ i ⇋ ∆ i . These relations correspond to ∆ ⇋ ∆ , ∆ ⇋ ∆ and ∆ ⇋ ∆ in Eq. 88. This symmetry and that of thesymmetric Γ matrix Γ βα = Γ αβ allows Eq. 88 to be written in matrix form as A × B × B × B × C × D × B × D × C × g V g V g V = − V V V (89) here A × = Γ Γ Γ Γ Γ Γ Γ Γ Γ , C × = Γ Γ Γ Γ Γ Γ Γ Γ Γ , B × = Γ Γ Γ Γ Γ Γ Γ Γ Γ , D × = Γ Γ Γ Γ Γ Γ Γ Γ Γ (90)Eq. 89 can be written as a non-Hermitian eigenvalue problem κ A κ B κ BB C DB D C g V g V g V = − g g V g V g V (91)where κ = g g . Substates V , V and V in gap equation Eq. 91, cannot be have different symmetries, so each of the eigenvectorsof Eq. 91 can be expressed in the compact form [ Φ n ] T = [ α sy V sy β sy V sy γ sy V sy ] (92)where subscript sy refers to one of the s -wave, d xy -wave or d x − y -wave symmetries, and the coefficients α sy , β sy and γ sy areto be determined. By inserting eigenvectors from Eq. 92 into the gap equation Eq. 91 one finds α sy β sy = α sy γ sy , γ sy β sy = β sy γ sy , J sy = c sy + β sy γ sy d sy + α sy γ sy b sy (93)where c sy = c s , c d , b sy = b s , b d , d sy = d s , d d and J sy = − g for each symmetry. Eq. 93 has two classes of solutions β sy = − γ sy ≡ ⇒ α sy = , J sy = c sy − d sy , β sy = + γ sy ≡ ⇒ b sy α sy + ( c sy + d sy − κ a sy ) α sy − κ b sy = α sy = α sy = −
2. The last solution in addition to island states α sy = Φ n and Φ n . But in the general case of symmetry breaking LiC characterized by gapequation Eq. 89, the quadratic equation in Eq. 94, has two temperature dependent solutions, α ± sy = J ± sy − c sy − d sy b sy , J ± sy = (cid:18) κ a sy + c sy + d sy ± q κ b sy + [ c sy + d sy − κ a sy ] (cid:19) (95) n t CC n ε c = − . t = . t ′ = . t t = − . t ′ = . t t = . t ′ ≈ t t = − . t ′ ≈ t t = − . m t LiLi m ε Li = . − .
30 0 .
09 0 . − . t LiC m − . Table 1.
The C-C hopping parameters (eV) for LiC are denoted by t CC n where the index n indicates n -th neighbours. In theLi plane, Li-Li hopping parameters are denoted by t LiLi m where m is m -th Li neighbor of central Li. The Li-C hoppingparameter is t LiC m . (cid:0)(cid:1)(cid:1)(cid:0)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:1) ξξ ξ123
01 1 23 3 444 45 66 778 88 89999 9 9 δδ δ1 23 ττ τ
LiC t LiC ξ ξ3 21 Β2Α2Α3Β1ΑΒ3Β2 3Β Α2Α1 h Α3 Β1ξ1 (a) (b) (c)
Figure 1. (Color online) (a) Schematic diagram of the shrunk graphene lattice, with the distortion emphasized. (b) Thehexagonal Li sheet, indicating the circles that Li neighbors lie on. (c) Diagram of the graphene decorated by lithium. The redLi atoms lie above the centers of the C hexagons. |η | m Γ Γ KM m=2m=1m=3 E β E αγ (k)(k)(k) E Figure 2. (Color online) A plot of the dispersion expressions | η m ( ~ k ) | , the three folded branches of pristine π ∗ band structurein the mini-Brillouin zone of graphene C . The Bloch wave character of E β is f d | d − ip > + f p | p − id > , of the E α is f d | d + ip > + f p | p + id > and for E γ is f s | s > + f f | f > . Here we use abbreviated notation f ( ) d | d ± ip > = f ( ) d x − y ( | d x − y > ± i | p x > ) and f ( ) p | p ± id > = f ( ) p y ( | p y > ± i | d xy > ) (a) (b) Γ Γ KM E ( e V ) k −1.5012.5 eV kx E ( e V ) Γ M K k k y Figure 3. (Color online) The left panel provides the band structure of lithium decorated graphene. The dashed lines indicatethe DFT bands, while the fitted bands are shown in color. The Fermi energy set to zero at µ = . eV . A small gap, E g = . eV is opened at the Γ point around -1.12 eV. The right panel provides a surface plot of the relatively flat band ofLiC . d -wave pairing dominates due to electrons in the valleys around saddle points at M . ααα α α α11 1111
23 467 5 αα α11α0 000 11 (c)(b) ∆∆∆ ∆∆∆∆ " """ ""
23 467
12 3123 1 1 ∆ ∆3 ∆ ∆ ∆ (a) Figure 4. (Color online) (a) Designation of the pairing amplitudes considered in this study, which cover all nearest neighborpairing possibilities denoted by ∆ n < i j > , ∆ ′ n < i j > and ∆ ′′ n < i j > where subscript < i j > has been dropped for brevity. (b) showsthe pairing amplitude for Φ + S phase with α ≈ . Φ − S phase with α ≈ − .
4. Both phases broken two band graphenesymmetry as can be seen by comparing symmetries along different bonds in seven atoms unit cell and two bands unit cellwhere its Bravais lattice points are labeled by 5, 6 and 7. (c) shows the pairing amplitude Φ + d xy where α ≈ Φ − d xy where α ≈ −
2. The first phase approximately preserves two band graphene symmetry while the others arise from broken symmetry. α ∆ α ∆−α ∆−α ∆ 0 0−∆−∆∆∆0 0 000 0 0 0∆−∆∆ ∆−∆−∆ (a) (b)
Figure 5. (Color online) Schematic diagram illustrating Cooper pair propagation for (a) Φ Jxy along a chain as shown by thegreen dashed line and arrow, and (b) localized “island” pairs for Φ f . (b)(a) −63*103*103*10 −4−2 210−1−2−3−4 B k T ( e V ) c α Φ S Φ S ΦΦ , xy Φ x − y θ s , xy Φ x − y ++ fpx , py Φ Φ Φ ΦΦΦ k T B c ( e V ) g (eV) Φ ++ + p y p xS fS dd x−y2 2xy ΦΦ dd x − y22xy Figure 6. (Color on line) (a) This phase diagram illustrates the relation between T c and the pairing potential g for Li C inwhich µ =
0. Panel (b) shows T c in terms of α sy s−wave−5d−wave−5s−wave−1d−wave−1 g (eV) Figure 7. (Color online) shows cooper pair interaction g in terms of doping ¯ n for d and s-wave phases for pristine grapheneat T= 0.1K. The solid (dashed) red line indicates d- wave (s- wave) pairing interaction in first nearest neighbor hopping t = . eV and similarly green line for accurate tight binding model can fit on DFT. For red line at the charge neutrality s-and d- wave are degenerate with g = .
76 while for full approximation they are not degenerate.
Φ ΦΦΦΦΦ Φ f , Φ , Φ , g (eV)0 + + + pypx (eV) µ d x − y2 2 d xy sd x − y2 2 d xy s Figure 8. (Color online) This diagram illustrates interaction potential g in terms of chemical potential µ at T c = . µ o − v = . eV (van Hove singularity) for symmetries Φ + d x − y , Φ + d xy ,and Φ + s the pairing potential decreases, then increases until a second critical value µ o − c =1.3 eV at which a phase transition to Φ + s occurs.occurs.