Superconductivity of strongly correlated electrons on the honeycomb lattice
aa r X i v : . [ c ond - m a t . s up r- c on ] A p r Superconductivity of strongly correlated electrons on the honeycomb lattice
A.A. Vladimirov a , D. Ihle b and N. M. Plakida a,c a Joint Institute for Nuclear Research, 141980 Dubna, Russia b Institut f¨ur Theoretische Physik, Universit¨at Leipzig, D-04109, Leipzig, Germany and c Max-Planck-Institut f¨ur Physik Komplexer Systeme D-01187, Dresden, Germany ∗ (Dated: May 1, 2019)A microscopic theory of the electronic spectrum and of superconductivity within the t - J model onthe honeycomb lattice is developed. We derive the equations for the normal and anomalous Greenfunctions in terms of the Hubbard operators by applying the projection technique. Superconductingpairing of d + id ′ -type mediated by the antiferromagnetic exchange is found. The superconducting T c as a function of hole doping exhibits a two-peak structure related to the van Hove singularitiesof the density of states for the two-band t - J model. At half-filling and for large enough values ofthe exchange coupling, gapless superconductivity may occur. For small doping the coexistence ofantiferromagnetic order and superconductivity is suggested. It is shown that the s -wave pairing isprohibited, since it violates the constraint of no-double-occupancy. PACS numbers: 71.27.+a,71.10.Fd,74.20.-z,74.20.Mn
I. INTRODUCTION
Graphene, the two-dimensional carbon honeycomb lat-tice, has been recently extensively studied due its pecu-liar electronic properties caused by the low energy cone-type electronic spectrum at the Dirac points K and K ′ in the Brillouin zone (BZ) (for a review see [1]). The ef-fects of electronic interactions play a minor role close tohalf-filling at low density of electronic states at the Diracpoints, but for large doping interactions appear to be im-portant. Studies of graphene beyond the simple modelof noninteracting electrons by taking into account theCoulomb interaction reveal a rich phase diagram withphase transitions to the antiferromagnetic (AF) state,spin-density wave (SDW), charge-density wave (CDW),and unconventional superconductivity (SC) (for a reviewsee [2]).The superconducting order parameters in the two-dimensional honeycomb lattice described by the hexag-onal symmetry group D h have a complex character. Ageneral symmetry analysis of available irreducible rep-resentations (IR) and superconducting order parametersis given in Ref. [3]. In the case of the singlet pairing,the extended s -wave s x + y order parameter ( E g IR), d x − y and d xy order parameters ( E g IR) preserving thetime-reversal symmetry and its time-reversal symmetrybreaking complex combination d x − y ± id xy ( d ± id ′ )are commonly discussed. The triplet pairing with the p x + ip y ( p + ip ′ ) order parameter ( E u IR) is also oftendiscussed in the literature. This symmetry considerationis very important in discussing SC in graphene. But itcan be also applied in studies of other electronic systemswith the two-dimensional honeycomb lattice, such as thetransition metal dichalcogenides [4].Several models of the Bardeen-Cooper-Schrieffer-type ∗ E-mail: [email protected] (BCS) were discussed. In Ref. [5] a model with the on-siteand the nearest-neighbor (nn) electron interactions of theBCS-type was considered. Assuming bond-independentanomalous correlation functions on the nn lattice sites,the s -wave pairing with k -independent gap induced bythe on-site interaction and the extended s -wave pairinginduced by the nn interaction were found. At the Diracpoints close to half filling the latter can be describedas p + ip phase. At large value of coupling constants agapless SC emerges at half-filling.The superconducting phase transition caused by theCoulomb interaction was studied within the Hubbardmodel[6] on the honeycomb lattice in a wide range of theon-site Coulomb repulsion U/t from weak to strong cou-pling (in graphene
U/t = 4 − t ∼ . U , the nn repulsion V , andthe spin-exchange interaction J . Close to half-filling, theSDW or CDW orders occur for large U and V , while fora large doping f -wave triplet-pairing and d + id ′ -wavesinglet-pairing emerge. In Ref. [8] the extended Hubbardmodel for graphene with the nn repulsion V and on-siteinteraction U was considered. Using the variational clus-ter approximation and the cellular dynamical mean-fieldtheory the SC of different symmetries was studied. De-pending on the values of V and U , the triplet p -wavesymmetry and the chiral combination p + ip ′ were found,while the singlet SC (extended s - or d -wave) was notclearly detected.Superconductivity on the honeycomb lattice is com-monly studied also within the phenomenological t – J model, where the exchange interaction ( J ∼ t ) is con-sidered as a fitting parameter. In Ref. [9] SC was studiedin a graphite layer within the resonating valence bonds(RVB) approach [10] for the t – J model. Both the ex-tended s -wave and d -wave pairing with the order param-eter determined by bond-dependent anomalous correla-tion functions were considered. The superconducting T c for the d -wave pairing appears to be much larger than forthe extended s -wave pairing and has a high value witha maximum at doping δ = 0 .
25. A similar model wasconsidered in Ref. [11] for the d + id ′ -pairing with thebond-dependent anomalous correlation functions. Thespectrum of electronic excitations in the superconductingstate is determined by two gaps g ± = | ∆( ± k ) | with spinup in one valley and spin down in the other valley withzeros at K and K ′ points of the BZ, respectively. Theexcitations are gapless at half filling for any value of thecoupling constant. The t – J model with the on-site inter-action of intermediate strength U/t = 2 . d + id ′ pairing fordoping 0 < δ < . T c can reach room temperatures atan optimal doping around δ = 0 .
15. The extended t – J model with the nn and the next-nearest-neighbor (nnn)exchange interaction J and J , respectively, was consid-ered in Ref. [13]. In the heavily doped case (around 3/8and 5/8 filling), a chiral d + id ′ symmetry was obtained.The competition between antiferromagnetism and SC inthe vicinity of half filling was considered by applying thefunctional RG.The quantum phase diagrams of both the Hubbardmodel and the t – J model on the honeycomb lattice at1/4 doping were studied in Ref. [14]. At this doping,in the nn tight-binding model the nested Fermi surfaceemerges which is unstable in the presence of a weakinteraction. Using a combination of exact diagonaliza-tion, density matrix renormalization group, the varia-tional Monte Carlo method, and quantum field theories,it was shown that in a wide range of the Hubbard re-pulsion, 1 < U/t <
40, or the exchange interaction,0 . < J/t < .
8, the quantum ground state is either achiral SDW state or a spin-charge-Chern liquid, but nota d + id ′ superconductor. For the t – J model at larger J/t > . d + id ′ superconductor occurs.A detailed study of the t – J model on the honeycomblattice was presented in Ref. [15]. Using the Grass-mann tensor product state approach, exact diagonaliza-tion and density-matrix renormalization methods, theground-state energy, the staggered magnetization in theAF phase, and the SC order parameter as a function ofdoping δ have been calculated. The occurrence of thetime-reversal symmetry breaking d + id ′ -wave SC up to δ = 0 . < δ < .
1, wherethe triplet pairing is induced (see also [16]). In Ref. [3] SCon the honeycomb lattice close to the Mott state at halffilling was studied within the t – J model using the renor-malized mean-field theory and in the Hubbard model byquantum Monte Carlo calculations. It was shown thatthe chiral d + id ′ -wave state is the most favorable statefor a wide range of the on-site repulsion U . At the sametime, a mixed chirality d -wave state, such as a state with d + id ′ -wave symmetry in one valley but d − id ′ -wave sym- metry in the other valley, is not possible in the t – J modelwithout reducing the translational symmetry. No ener-getically favorable d -wave solution with an overall zerochirality was found within the t – J model.The van Hove singularity (VHS) scenario of SC beingdeveloped for cuprates (see, e.g., [17]) was discussed inseveral publications. In Ref. [18] the extended VHS indoped graphene was found using the angle-resolved pho-toemission spectroscopy. Considering the conventionalfluctuation exchange approximation [19] with the weakHubbard interaction U/t .
4, the competition betweenmagnetic instability and SC was analyzed. It was foundthat SC plays a dominant role when the Fermi level isplaced close enough to the extended VHS, where the tran-sition temperature T c can be quite high. In Ref. [20] itwas shown that, due to the strong anisotropy of the elec-tron scattering at the VHS, attractive coupling channelsappear from the originally repulsive interaction that re-sults in the superconducting pairing with T c ∼
10 K.In Ref. [21] studies of the Hubbard model on the hon-eycomb lattice with nn and nnn interactions show theappearance of the extended VHS, where the density ofstates diverges in a power law. Using the random-phase-approximation and determinant quantum Monte Carloapproaches a possible triplet p + ip ′ SC with relativelyhigh T c was found at low filling. The interplay betweenSC and SDW order in graphene close to the VHS wasconsidered in Ref. [22]. The instabilities to both the chi-ral d + id ′ SC and the uniaxial SDW were found in amodel with four different interactions between fermionsnear saddle points. The SC is strongest at the VHS,but slightly away from it SDW appears first. To inves-tigate the possibility of coexistence of SC and SDW, theLandau-Ginzburg functional was derived. It was shownthat SDW does not coexist with SC, because both phasesare separated by first-order transitions. The dynamiccluster approximation was used in Ref [23] to study SCin the Hubbard model with
U/t = 2 −
6. A transitionfrom the d + id ′ -wave singlet pairing, which dominatesclose to the VHS filling, to the p -wave triplet pairing atlarger coupling was found.In several studies the renormalized mean field theoryfor the t – J model was employed. To take into accountstrong correlations of electrons in the singly occupiedband, the hopping parameter t and exchange interaction J were renormalized by the statistical weighting factors g t = 2 δ/ (1 + δ ) and g J = 4 / (1 + δ ) , as in the RVBtheory for cuprates [24]. We point out that this renor-malization is rather crude, e.g., for δ = 0 it results in thezero band-width ∼ δt though at low doping spin-polaronquasiparticles appear with a finite bandwidth of the order J (see, e.g. [25]). Moreover, in the undoped case the fourtimes stronger exchange interaction 4 J results while theHeisenberg model with the original exchange interaction J should be found. The slave-boson approximation alsostrongly violates the statistics of the projected electronsin the original t – J model [26]. To take into account therestriction of no-double-occupancy in the t – J model, atechnique for the projected electron operators, the Hub-bard operators [27], should be used.In the present paper we develop a microscopic the-ory of SC of strongly correlated electrons on the hon-eycomb lattice employing the Hubbard operator tech-nique. This technique was used in our previous paperfor the calculation of the electronic spectrum, the spin-excitation spectrum, and of thermodynamic quantitieswithin the t – J model [28]. Using the projection operatortechnique [29] developed for the thermodynamic Greenfunctions (GFs) [30] in Ref. [31], we derive equations forthe normal and anomalous (pair) GFs. In the general-ized mean-field approximation (GMFA) a self-consistentsystem of equations for the singlet order parameters isobtained and the superconducting T c as a function ofdoping δ is calculated. It is shown that the condition ofthe no-double-occupancy of quantum states in the t – J model is violated for the s -wave pairing, while the d + id ′ pairing obeys this restriction.The paper is organized as following. In Section II the t – J model for the honeycomb lattice is formulated. Equa-tions for the GFs are derived in Section III. Gap equa-tions and the calculation of T c are given in Section IV.The conclusion is presented in Section V. II. THE t - J MODEL
The Hubbard model [6] is commonly used in studiesof correlated electronic systems. In the limit of strongcorrelations the model is reduced to the one-subband t - J model [26]. In the lattice site representation the modelreads: H = − t X h i,j i σ ˜ a + i,σ ˜ a j,σ − µ X i,σ n i,σ + J X h i,j i (cid:18) S i S j − n i n j (cid:19) , (1)where ˜ a + i,σ = a + i,σ (1 − n i, ¯ σ ) and ˜ a iσ = a iσ (1 − n i, ¯ σ ) areprojected creation and annihilation electron operators onthe site i with spin σ/ σ = ± , ¯ σ = − σ ), and thenumber operator n i = P σ ˜ a + i,σ ˜ a i,σ . Here h i, j i denotethe nearest neighbors for electron hopping with energy t and for spins S i with AF exchange interaction J .It is convenient to employ the Hubbard operator (HO)technique [27] where the projected electron operatorsare written as: ˜ a + iσ = X σ i , ˜ a jσ = X σj . The HOs X nmi = | i, n ih i, m | describe transitions between threepossible states at a lattice site i : | i, n i = | i, i and | i, σ i for an empty site and for a singly occupied site by anelectron with spin σ/
2, respectively.The electron number operator and the spin operatorsare defined as n i = X σ X σσi = X ++ i + X −− i , (2) S σi = X σ ¯ σi , S zi = ( σ/
2) [ X σσi − X ¯ σ ¯ σi ] . (3) The commutation relations for the HOs read: (cid:2) X nmi , X klj (cid:3) ± = δ ij (cid:0) δ mk X nli ± δ nl X kmi (cid:1) . (4)The upper sign refers to Fermi-type operators such as X σi , while the lower sign refers to Bose-type operatorssuch as n i (2) or the spin operators (3). The completenessrelation for the HOs, X i + P σ X σσi = 1, rigorouslypreserves the constraint of no-double-occupancy of thequantum state | i, n i on any lattice site i .In terms of HOs the t - J model (1) takes the form: H = − t X h i,j i σ X σ i X σj − µ X iσ X σσi + J X h i,j i σ (cid:0) X σ ¯ σi X ¯ σσj − X σσi X ¯ σ ¯ σj (cid:1) . (5)We consider the t - J model on the honeycomb latticewhich is bipartite with two triangular sublattices A and B . Each of the N sites on the A sublattice is connectedto three nn sites α = 1 , , B sublatticeby vectors δ α , and N sites on B are connected to A byvectors − δ α : δ = a √ , − , δ = − a √ , , δ = a (0 , . (6)The basis vectors are a = δ − δ = ( a / √ ,
3) and a = δ − δ = ( a / −√ , a = | a | = | a | = √ a , where a is the nn distance;hereafter we put a = 1. The reciprocal lattice vectorsare k = (2 π/ √ ,
1) and k = (2 π/ −√ , i → iA, iB .The chemical potential µ depends on the average elec-tron occupation number n = n A = n B = 1 N X i,σ h X σσiA i , (7)where N is the number of unit cells and h ... i denotes thestatistical average with the Hamiltonian (5). III. GREEN FUNCTIONSA. Equations for the Green functions
To consider SC within the model (5), we introduce theanticommutator retarded matrix GF [30] G ijσ ( t − t ′ ) = − iθ ( t − t ′ ) h{ ˆ X iσ ( t ) , ˆ X † jσ ( t ′ ) }i≡ hh ˆ X iσ ( t ) , ˆ X † jσ ( t ′ ) ii , (8)where { X, Y } = XY + Y X , X ( t ) = e iHt X e − iHt ( ~ = 1),and θ ( x ) is the Heaviside function. Here we use Nambunotation and introduce the vector Hubbard operatorsˆ X iσ and ˆ X † jσ with 4 components:ˆ X iσ = X σiA X σiB X ¯ σ iA X ¯ σ iB , ˆ X † jσ = (cid:0) X σ jA X σ jB X σjA X σjB (cid:1) . (9)The Fourier representation in ( k , ω )-space is defined by G ijσ ( t − t ′ ) = Z ∞−∞ dω π e − iω ( t − t ′ ) G ijσ ( ω ) , G ijσ ( ω ) = 1 N X k e i k ( r i − r j ) G σ ( k , ω ) . (10)The 4 × G σ ( k , ω ) = (cid:18) ˆ G ( k , ω ) ˆ F σ ( k , ω )ˆ F † σ ( − k , ω ) − ˆ G ( k , − ω ) (cid:19) , (11)where the components of the normal GF read asˆ G ( k , ω ) = (cid:18) G AA ( k , ω ) G AB ( k , ω ) G BA ( k , ω ) G BB ( k , ω ) (cid:19) , (12)and the components of the anomalous GF are given byˆ F σ ( k , ω ) = (cid:18) F σAA ( k ω ) F σAB ( k , ω ) F σBA ( k , ω ) F σBB ( k , ω ) (cid:19) . (13)To calculate the GF (8), we use the equation of motionmethod. Differentiating the GF with respect to time t weobtain ω G ijσ ( ω ) = δ ij Q + hh [ ˆ X iσ , H ] , ˆ X † jσ ii ω , (14)where Q = h{ ˆ X iσ , ˆ X † iσ }i = ˜ τ Q . Here, ˜ τ is the 4 × Q = h X iβ + X σσiβ i = 1 − n/ t – J model it is convenient to choose the mean-field con-tribution in the equations of motion (14) as the zeroth-order quasiparticle (QP) energy. We calculate it in theGMFA using the projection operator method [31]. In thisapproach we write the operator [ ˆ X iσ , H ] in (14) as a sumof the linear part, proportional to the single-particle op-erator ˆ X iσ , and the irreducible part ˆ Z ( ir ) iσ orthogonal toˆ X iσ : ˆ Z iσ = [ ˆ X iσ , H ] = X l E ilσ ˆ X lσ + ˆ Z ( ir ) iσ . (15)The orthogonality condition h{ ˆ Z ( ir ) iσ , ˆ X † jσ }i = 0 deter-mines the linear part, the zeroth-order QP energy: E ijσ = h{ [ ˆ X iσ , H ] , ˆ X † jσ }i Q − = (cid:18) ˆ E ij ˆ∆ ijσ ˆ∆ ∗ jiσ − ˆ E ji (cid:19) , (16) where ˆ E ij and ˆ∆ ijσ are the normal and anomalous com-ponents of the matrix. The corresponding zeroth-orderGF in (14) in the Fourier representation (10) is given by G σ ( k , ω ) = (cid:16) ω ˜ τ − E σ ( k ) (cid:17) − Q , (17) E σ ( k ) = (cid:18) ˆ E ( k ) ˆ∆ σ ( k )ˆ∆ ∗ σ ( k ) − ˆ E ( k ) (cid:19) . (18)It is possible to calculate the self-energy operator givenby the many-particle GF hh ˆ Z ( ir ) iσ | ˆ X † jσ ii ω in (14) and toderive the Dyson equation for the GF (8), as has beendone in our previous publications (see [32, 33]). In thepresent paper we consider the theory in GMFA withinthe zeroth-order approximation for the GF (17).The components of the energy matrix (16) are deter-mined by the commutators:ˆ E ij = h{ [ (cid:18) X σiA X σiB (cid:19) , H ] , (cid:0) X σ jA X σ jB (cid:1) }i Q − , (19)ˆ∆ ij,σ = h{ [ (cid:18) X σiA X σiB (cid:19) , H ] , (cid:0) X σjA X σjB (cid:1) }i Q − . (20)Performing commutations and introducing the Fourierrepresentation, X σiA = (1 / √ N ) P k e i kr i X σ k A , we obtain:ˆ E ( k ) = (cid:18) ε A ε AB ( k ) ε BA ( k ) ε B (cid:19) , (21) ε A = h{ [ X σ k A , H ] , X σ k A }i Q − = − e µ,ε B = h{ [ X σ k B , H ] , X σ k B }i Q − = − e µ,ε AB ( k ) = h{ [ X σ k A , H ] , X σ k B }i Q − = − e t γ ( k ) ,ε BA ( k ) = ε AB ( k ) ∗ = − e t γ ∗ ( k ) , (22)where γ ( k ) = P α exp( i k −→ δ α ) and | γ ( k ) | =1 + 4 cos( √ k x / √ k x /
2) + cos(3 k y / e µ and hopping param-eter e t were calculated in Ref. [28] and are given by therelations: e µ = µ − tQ D + 3 J n − J Q C , (23) e t = t Q (cid:18) C Q (cid:19) + J D Q . (24)Here the nn correlation functions for electrons and spinsare: D = h X σ iA X σi + δ ,B i , C = h S ziA S zi + δ ,B i . (25)The anomalous energy matrix for the gaps reads:ˆ∆ σ ( k ) = (cid:18) ∆ Aσ ∆ ABσ ( k )∆ BAσ ( k ) ∆ Bσ (cid:19) , (26)∆ Aσ = ∆ Bσ ≡ ∆ σ = h{ [ X σ k A , H ] , X σ − k A }i Q − = 2 t X l h X σlB X σiA i Q − , (27)∆ ABσ ( k ) = h{ [ X σ k A , H ] , X σ − k B }i Q − = − JQ X r j = r i + δ α exp[ i k ( r j − r i )] h X σjB X σiA i , (28)∆ BAσ ( k ) = h{ [ X σ k B , H ] , X σ − k A }i Q − = − ∆ AB ¯ σ ( − k ) = ∆ ABσ ( − k ) . (29)Note that both gap functions ∆ σ and ∆ ABσ ( k ) are deter-mined by the nn correlation functions h X σjB X σiA i , sincewe have no pairing on one lattice site contrary to Ref. [5].Using Eqs. (21) and (26) we obtain the energy matrix: E σ ( k ) = − e µ − e t γ ( k ) ∆ σ ∆ ABσ ( k ) − e t γ ∗ ( k ) − e µ ∆ ABσ ( − k ) ∆ σ ∆ ∗ σ ∆ ∗ ABσ ( − k ) e µ e t γ ( k )∆ ∗ ABσ ( k ) ∆ ∗ σ e t γ ∗ ( k ) e µ . (30)We point out that the matrix (30) is similar to the ma-trix in the superconducting state in MFA for grapheneobtained in Ref. [5] and in Ref. [11] but with the renor-malized chemical potential e µ and hopping parameter e t .The GF in Eq. (17) is defined by the inverse matrix (cid:16) ω ˜ τ − E σ ( k ) (cid:17) − . Its calculation results in the GF: G σ ( k , ω ) = Q e D ( k , ω ) A T σ ( k , ω ) , (31)where A T σ ( k , ω ) is the transposed matrix of the cofactorsof the matrix (cid:16) ω ˜ τ − E σ ( k ) (cid:17) .The diagonal and off-diagonal normal G αβ ( k , ω ) andanomalous F αβσ ( k , ω ) GFs components in the matrix(11) read: G AA ( k , ω ) = hh X σ k A , X σ k A ii ω = A Q e D ( k , ω ) , (32) G AB ( k , ω ) = hh X σ k A , X σ k B ii ω = A Q e D ( k , ω ) , (33) F σAA ( k , ω ) = hh X σ k A , X σ − k A ii ω = A σ Q e D ( k , ω ) , (34) F σAB ( k , ω ) = hh X σ k A , X σ − k B ii ω = A σ Q e D ( k , ω ) . (35)The coefficients A nm are given by the equations: A ( k ) = ( ω − e µ )( ω − e µ ) − ∆ ∗ ABσ ( − k ) ∆ ∗ σ e tγ ( k ) − ∆ σ ∆ ∗ ABσ ( − k ) e t γ ∗ ( k ) − | ∆ σ ) | ( ω − e µ ) − ( ω + e µ ) e t | γ ( k ) | − | ∆ ABσ ( − k ) | ( ω − e µ ) , (36) A ( k ) = − e t γ ( k ) ( ω − e µ ) + | ∆ σ | e t γ ( k )+∆ ABσ ( k ) e t γ ∗ ( k ) + ∆ ABσ ( k ) ∆ ∗ σ ( ω − e µ )+ e t γ ( k ) | γ ( k )) | + ∆ σ ∆ ABσ ( k ) ( ω − e µ ) , (37) A σ ( k ) = − ∆ ABσ ( − k ) e t γ ( k ) ( ω − e µ ) − ∆ σ | ∆ Bσ | +∆ ABσ ( k ) e t γ ∗ ( k ) ( ω + e µ ) + ∆ ABσ ( k ) ∆ ABσ ( − k ) ∆ ∗ σ +∆ σ ( ω − e µ − e t | γ ( k ) | ) , (38) A σ ( k ) = − ∆ ABσ ( − k ) e t γ ( k ) + ∆ σ ∆ ABσ ( k )+∆ ABσ ( k )( ω − e µ ) − ∆ ABσ ( k ) ∆ ABσ ( − k )+∆ σ e µ e t γ ( k ) . (39)The determinant e D ( k , ω ) of the matrix, e D ( k , ω ) = k ω ˜ τ − E σ ( k ) k , (40)gives the equation for the spectrum in the superconduct-ing state e D ( k , ω ) = 0 which can be written in a generalcase (with notations ∆ σ ± ≡ | ∆ ABσ ( ± k ) | ) as: ω − ω (∆ σ + 12 (∆ σ + + ∆ σ − ) + e t | γ ( k ) | + e µ )+∆ σ + ∆ σ − + e t | γ ( k ) | + e µ + ∆ σ +2 e t Re[ γ ( k ) ∆ ∗ ABσ ( k )∆ ABσ ( − k )] − e t | γ ( k ) | e µ +(∆ σ + + ∆ σ − ) e µ + 2∆ σ ( e µ + e t | γ ( k ) | ) − e t e µ ∆ σ Re[∆
ABσ ( − k ) γ ( k ) + ∆ ABσ ( k ) γ ∗ ( k )] − σ Re[∆
ABσ ( k )∆ ABσ ( − k )] = 0 . (41)The solution of this equation for the extended s -pairing(see Eq. (54)) gives the spectrum of excitations:Ω ( k ) = ( e µ ± e t | γ ( k ) | ) + (∆ σ ± | ∆ ABσ ( k ) | ) , (42)which coincides with the spectrum for graphene in MFAfound in Ref. [5]. Note that for the gap ∆ σ = 0 thespectrum is gapless at e µ = 0 at six corners of the BZat K, K ′ points: Ω s ( k ) = e t | γ ( k ) | q ABσ / e t , as waspointed out in [5].The spectrum for the d -pairing (see Eq. (55)) result-ing from Eq. (41) with ∆ σ = 0 has a more complicatedstructure with two gaps ∆ σ ± :Ω ( k ) = 12 (∆ σ + + ∆ σ − ) + e t | γ ( k ) | + e µ ± ne t | γ ( k ) | (4 e µ + ∆ σ + + ∆ σ − ) + 14 (∆ σ + − ∆ σ − ) − e t Re[ γ ( k ) ∆ ∗ ABσ ( k )∆ ABσ ( − k )] o / . (43)A similar spectrum was found for the d -wave symmetryin Ref. [11]. Since the d -wave gap (55) is zero at three cor-ners of the BZ at K points and has a maximum value atanother three corners of the BZ at K ′ points, the spec-trum (43) is gapless at e µ = 0 at K points and has amaximum value at K ′ points (see also Refs. [2, 13]). B. Normal state Green function
The determinant for the normal state GFs (32) and(33) reads e D ( k , ω ) = D ( k , ω ) D ( k , − ω ) ,D ( k , ω ) = [ ε + ( k ) − ω ][ ε − ( k ) − ω ] . (44) . . . . . d t N ( m ~ ) d − − − m ~ /t FIG. 1: (Color online) Density of electronic states at theFermi energy t N ( e µ ) (solid blue line, left axis) and the chem-ical potential e µ/t (dashed red line, right axis) vs doping δ . Here the electronic spectrum is given by the matrix (30)in the normal state (cf. Ref. [28]): ε ± ( k ) = − e µ ± e t | γ ( k ) | . (45)The spectrum has the Dirac cone-type behaviour at K and K ′ points at the corners of the BZ as in graphene.The cones touch the Fermi surface (FS) at e µ = 0 for halffilling at the electron occupation number n = 2 / t - J model. The detailed doping dependence of the FSwas considered in Ref. [28].For the normal state GFs (32) and (33) with the coef-ficients A (36) and A (37) we obtain: G AA ( k , ω ) = Q h ω − ε + ( k ) + 1 ω − ε − ( k ) i , (46) G AB ( k , ω ) = − Qγ ( k )2 | γ ( k ) | h ω − ε + ( k ) − ω − ε − ( k ) i . (47)The density of electronic states (DOS) at the Fermi en-ergy e µ in the normal state is determined by the GF (46),and for the dispersion ε ± ( k ) (45), is given by the equa-tion: N ( e µ ) = 1 N X k ,σ [ − π ] Im G AA ( k , iǫ )= QN X k [ δ ( e µ − e t | γ ( k ) | ) + δ ( e µ + e t | γ ( k ) | )] . (48)The DOS at the Fermi energy N ( e µ ) and the chemicalpotential e µ for − ≤ e µ/t ≤ . ≤ δ ≤
1. The two peaks ofthe DOS correspond to the VHS in two bands. The DOS is nonsymmetric with respect to e µ = 0 in comparisonwith the DOS of graphene due to the renormalization ofthe hopping parameter in (48), where we use e t = t Q = t (1 + δ ) / n Aσ ( k ) = h X σ k A X σ k A i and h X σ k B X σ k A i in Ref. [28] are correct. C. Anomalous Green functions
To calculate the superconducting T c we use the lin-earized approximation for the anomalous GFs (34) and(35). They are determined by the relations: F σAA ( k , ω ) = Q A σ ( k , ω )[ ω − ε ( k )][ ω − ε − ( k )] , (49) A σ ( k , ω ) = − ∆ ABσ ( − k ) ( ω − e µ ) e t γ ( k )+∆ ABσ ( k )( ω + e µ ) e t γ ∗ ( k )+∆ σ ( ω − e µ − e t | γ ( k ) | ) , (50) F σAB ( k , ω ) = Q A σ ( k , ω )[ ω − ε ( k )][ ω − ε − ( k )] , (51) A σ ( k , ω ) = ∆ ABσ ( k )( ω − e µ ) − ∆ ABσ ( − k ) e t γ ( k ) + ∆ σ e µ e t γ ( k ) . (52)The corresponding anomalous correlation function F σBA ( k ) = h X σ − k B X σ k A i in the gap equations (27), (28)determined by the GF (51) is given by F σBA ( k ) = − π Z ∞−∞ dω exp( ω/T ) + 1 Im F σAB ( k , ω )= Qε − ( k ) − ε ( k ) n A σ ( ω = ε + ( k ))2 ε + ( k ) tanh ε + ( k )2 T − A σ ( ω = ε − ( k ))2 ε − ( k ) tanh ε − ( k )2 T o . (53) IV. GAP EQUATIONS AND T c Let us consider solutions for the gaps of different sym-metries. In particular, the extended s -wave gap (28) de-termined by the bond-independent anomalous correlationfunction (53) can be written as:∆ ABσ ( k ) = ∆ ABσ X α exp[ i ( k δ α ] = ∆ ABσ γ ( k ) . (54)The d + id ′ -wave pairing is determined by the gap whichhas different phases on bonds (Refs. [9, 11, 15]):∆ ABσ ( k ) = X α ∆ ασ exp[ i k δ α ] , ∆ ασ = ∆ σ e i (2 π/ α . (55)The numerical solution of the gap equations yields T c as a function of doping given below, where we use therenormalized hopping parameter e t = t (1 + δ ) / T c is the temperature of Cooper-pair forma-tion without superconducting long-range order (see, e.g.,Refs. [2, 5]).It should be pointed out that for the conventionalFermi-liquid we can have the s -wave pairing on a singlesite given by the anomalous correlation function h a i ¯ σ a iσ i .In terms of HOs this pairing can be described by theequation: h a i ¯ σ a iσ i = h X i i = h X σi X ¯ σ i i 6 = 0 , (56)which shows that the double occupancy of one lattice siteis permitted but in the two Hubbard subbands. For theHubbard model in the limit of strong correlations, i.e.,for the t - J model, this pairing is prohibited due to theno-double-occupancy restriction which in terms of HOs,as was proposed in Ref. [34], can be written as: h X σi X σi i = h a i ¯ σ (1 − n i,σ ) a iσ (1 − n i, ¯ σ ) i ≡ . (57)As shown in Sect. IV D, both the k -independent s -wavepairing with the gap (27) and the extended s -wave pair-ing with the gap (54) violate this condition and must beexcluded from a rigorous point of view. For the d -pairingwith the gap (55) the condition (57) is fulfilled.However, in several publications this constraint hasnot been rigorously taken into account, for example us-ing the phenomenological t - J model (see Ref. [9] wherethe double-site occupancy is included in the Hilbertspace), or introducing the statistical weighting factors g t = 2 δ/ (1 + δ ) and g J = 4 / (1 + δ ) in the t - J model(see, Refs. [3, 13]), as discussed in Sect. I, or using theGutzwiller projector P g = Q i (1 − g n i ↑ n i ↓ ) with g < k -independent s -wave pairing in Sect. IV A andthe extended s -wave pairing in Sect. IV B. A. k-independent s -wave gap equation For the k -independent s -wave gap (27) we have∆ σ = 2 tQ X l h X σlB X σiA i = 2 tQN X α X q exp[ − i q δ α ] F σBA ( q ) . (58)Assuming that t ≫ J we can solve the equation for ∆ σ considering the q -dependent part of the gap ∆ ABσ ( q ) in A σ ( ω, q ) (52) as a small perturbation. We obtain thegap equation:∆ σ = ∆ σ tN X q | γ ( q ) | h ε + ( q ) tanh ε + ( q )2 T c − ε − ( q ) tanh ε − ( q )2 T c io . (59)We find a solution for T c only for hole doping δ ≤ . µ >
0. It has a high maximum value, T c ∼ . t ,due to the strong pairing interaction t ≫ J . Note thatthe coupling proportional to the hopping parameter t isdue to the kinematical interaction induced by the com-mutation relations (4) for the HOs. The same type ofpairing occurs in the t – J model on the square lattice(see Ref. [34]) which has been disregarded, since it vi-olates the constraint of no-double-occupancy (57). Asshown in Sect. IV D, the s -wave pairing described by the k -independent gap ∆ σ (58) on the honeycomb lattice alsoviolates the constraint (57), and in what follows we takethe solution ∆ σ = 0. B. Extended s -wave gap equation Using Eq. (53) the gap equation (28) reads∆
ABσ ( k ) = J N X q X α exp[ i ( k − q ) δ α ]∆ ABσ ( q ) × h ε + ( q ) tanh ε + ( q )2 T c + 12 ε − ( q ) tanh ε − ( q )2 T c i − JN X q X α exp[ i ( k − q ) δ α ] e t µ ˜ t | γ ( q ) |× h ∆ ABσ ( q ) | γ ( q ) | − ∆ ABσ ( − q ) γ ( q ) i × h ε + ( q ) tanh ε + ( q )2 T c − ε − ( q ) tanh ε − ( q )2 T c i . (60)Taking into account that the second term vanishes for thebond-independent gap ∆ ABσ ( k ) = ∆ ABσ γ ( k ) becauseof γ ( q ) | γ ( q ) | − γ ∗ ( q ) γ ( q ) = 0 , we obtain the gapequation γ ( k ) = J N X q X α exp[ i ( k − q ) δ α ] γ ( q ) × h ε + ( q ) tanh ε + ( q )2 T c + 12 ε − ( q ) tanh ε − ( q )2 T c i . (61)The solution of this equation shows that T c ( δ ) linearly in-creases with doping and vanishes for large δ dependingon J , e.g., T c ( δ ) = 0 at δ > . J = t/
2. The maxi-mum value of T c rapidly increases with the interaction J ,e.g., T max c = 0 . t (0 . t ) for J = t/ t/ s -wave pairing also violates the constraintof no-double-occupancy (57), we consider it to compareour results with calculations for the t - J model, where thisrestriction has not been rigorously taken into account. C. d -wave gap equation In the case of d -wave symmetry ( d + id ′ ) for the bond-dependent gap (55) we have the same equation (60). Thedirect numerical calculation for this gap reveals that thesecond term in the integral gives no contribution. It isconvenient to consider the equation for a particular valueof α for the phase sensitive contribution:∆ ασ exp[ i ( k δ α ] = ∆ σ e i (2 π/ α exp[ i k δ α ]= J N X q exp[ i ( k − q ) δ α ] X β ∆ σ e i (2 π/ β exp[ i q δ β ] × h ε + ( q ) tanh ε + ( q )2 T c + 12 ε − ( q ) tanh ε − ( q )2 T c i . (62)Cancelling out the terms ∆ σ exp[ i k δ α ] we write thisequation as1 = J N X q X β exp[ i q ( δ β − δ α )] e i (2 π/ β − α ) × h ε + ( q ) tanh ε + ( q )2 T c + 12 ε − ( q ) tanh ε − ( q )2 T c i . (63)Considering the vectors a βα ≡ δ β − δ α we can see thatdifferent values of β − α correspond to a set of threevectors, one is equal to zero, and the other two are thelattice vectors. The change of α rotates the whole setby 2 π/ C symmetry of the lattice so thatthe above equation does not depend on α . Therefore, allthree equations for various α are the same. The equationfor α = 0 reads:1 = JN X q X β exp[ i ( q ( δ β − δ )] e i (2 π/ β × h ε + ( q ) tanh ε + ( q )2 T c + 12 ε − ( q ) tanh ε − ( q )2 T c i . (64)This equation was solved numerically for various valuesof J/t . The results for T c as a function of the hole doping δ are depicted in Fig. 2 and Fig. 3. T c ( δ ) has a two-peakstructure with the maxima corresponding to the VHS inthe DOS shown in Fig. 1.The transition temperature T c rapidly increases withthe interaction J/t in Eq. (64), as was also found inRef. [9]. For small
J/t = 1 / T c with δ and T c = 0 (or exponentially small) for δ < .
05 similar to the d -wave pairing in the one-band t – J model on the square lattice (see Ref. [32]). In Ref. [15]for J/t = 1 / T = 0 wasfound for 0 < δ < .
4, but in Ref. [9] T c is exponentiallysmall for δ < .
05 when J = 0 . t . For J = 3 t/ T c ≈ . t is comparable withthe maximum value of T c ≈ . t in Ref. [3] for J/t = 0 . T c ≈ . t in Ref. [9] for J/t = 0 .
8. Aswas claimed in these publications, such values of T c withthe hopping parameter t ≈ . T c SC of electrons on the honeycomblattice at optimal doping. The variational Monte Carlo . . . . . d T c /t J = t/2J = t/3
FIG. 2: (Color online) T c ( δ ) for the d -wave pairing with J = t/ J = t/ . . . . d T c /t J = tJ = 3t/4
FIG. 3: (Color online) T c ( δ ) for the d -wave pairing with J =3 t/ J = t (bold line). calculation for the RVB theory in Ref. [12] shows that T c can be estimated as about twice of the room temper-ature. For small values of J < J c = 0 . t there is nopairing at half-filling for ˜ µ = 0. For large values of J/t a sharp increase of T c with δ is found in Fig. 3, and for J > J c T c is nonzero at ˜ µ = 0. In this region we canobserve the gapless SC, as was also found in Ref. [5] forthe s -wave pairing in graphene for the nn BCS couplingparameter g > g c .Taking into account the results of Ref. [28] for the t – J model on the honeycomb lattice, where AF long-range or-der at T = 0 was observed for δ . .
1, we argue that theAF and superconducting ground-state order may coexistin the range δ . .
1, in agreement with the numericalcalculations in Ref. [15].Let us compare the results for the extended s -wavepairing and d -wave pairing. As was found in Sect. IV B,the maximum value of T max c ≈ . t for J = t/ d -wave pairing in Fig. 2.It contradicts to the results of Ref. [9], where T c for the d -wave pairing is much larger than for the extended s -wavepairing. Note that in Ref. [9] the mean-field RVB theorywas used, where the restriction of no-double-occupancyhas not been taken into account. On the other hand, inRefs. [3, 13], where the renormalized mean-field approx-imation determined by the parameters g t and g J for the t – J model was employed, the maximum value of the orderparameter (and T c ) for the extended s -wave pairing andthe d -wave pairing are quite close. So we can say that ac-counting for strong correlations both in our approach andin the renormalized mean-field approximation results inclose maximum values of T c for the extended s -wave and d -wave pairing. However, it should be stressed that theconstraint (57) excludes the s -wave pairing. This conclu-sion is supported by the calculations for the t – J modelin Ref. [15], where the projected character of the electronoperators has been implemented and only d -wave pairinghas been found. D. Constraint for the s-wave pairing
Let us consider the restriction of no-double-occupancy(57) for the Hubbard operators: h X σiA X σiA i = 0 . (65)Using the GF F σAA ( k , ω ) (49) this condition reads h X σiA X σiA i = 1 N X k Z ∞−∞ dω exp( ω/T ) + 1 × [ − π ] Im Q A ( k , ω )[ ω − ε ( k )][ ω − ε − ( k )] = 0 . (66)where the function A σ ( k , ω ) is given by (50). This con-dition for the s -wave pairing determined by the gap func-tion ∆ σ (27), where A σ ( ω = ε ± ( k )) = ∓ σ e µ e t | γ ( k ) | ,results in the relation: h X σiA X σiA i = − Q ∆ σ N X k h ε + ( k ) tanh ε + ( k )2 T + 12 ε − ( k ) tanh ε − ( k )2 T i . (67)The summation over k does not vanish for a positivelydefined integrand. Therefore, the s -wave pairing with the k -independent gap ∆ σ violates the kinematic restriction(66) and is ruled out.For the extended s -wave pairing with the gap (54) wehave A σ ( ω = ε ± ( k )) = 2∆ ABσ e µ e t | γ ( k ) | , and the re-lation (66) reads: h X σiA X σiA i = − Q ∆ ABσ N X k | γ ( k ) | h ε + ( k ) tanh ε + ( k )2 T − ε − ( k ) tanh ε − ( k )2 T i . (68)The summation over k does not vanish except for e µ = 0when ε + ( k ) = − ε − ( k ). For other doping the correlationfunction (68) is non-zero that violates the kinematic re-striction (66), and the extended s -wave pairing cannotbe realized.Finally, let us consider the d -wave pairing with thebond-dependent gap (55). In this case, for the function(50) we have A ( ω = ε ± ( k )) = 2 e µ e t γ ( k ) X α ∆ ασ exp[ − i k δ α ] ± e t | γ ( k ) | X α ∆ ασ i Im[ γ ∗ ( k ) exp( i k δ α )] , (69)where the relation γ ∗ ( k ) exp( i k δ α ) − γ ( k ) exp( − i k δ α ) =2 i Im( γ ∗ ( k ) exp( i k δ α )) was used. For the correlationfunction (66) we obtain: h X σiA X σiA i = Q N X α ∆ ασ X k γ ( k ) | γ ( k ) | exp( − i k δ α ) × h ε + ( k ) tanh ε + ( k )2 T − ε − ( k ) tanh ε − ( k )2 T i + Q e t e µ N X k X α ∆ ασ i Im[ γ ∗ ( k ) exp( i k δ α )] × h ε + ( k ) tanh ε + ( k )2 T + 12 ε − ( k ) tanh ε − ( k )2 T i . (70)In the first term the summation over k does not dependon α due to the C symmetry of γ ( k ) and ε ± ( k ). There-fore, the summation over α of the gap function ∆ ασ can bedone independently that results in the vanishing of thefirst term: P α ∆ ασ = ∆ σ P α exp( i (2 π/ α ) = 0. Thesecond term also gives no contribution due to summationover k of the odd in k function Im[ γ ∗ ( k ) exp( i k δ α )] = i P β sin[ k ( δ α − δ β )]. So, the d -wave pairing with thegap (55) does not violate the condition of no-double-occupancy (66). This conclusion was checked by the di-rect integration over k in Eq. (70). V. CONCLUSION
In this paper a microscopic theory of superconductivityin electronic systems with strong correlations is presentedwithin the t – J model on the honeycomb lattice. Theconstraint of no-double-occupancy in the two-band t – J model is rigorously taken into account by employing theHO technique. The superconducting T c as a functionof doping is calculated for the d + id ′ gap function. Itreveals a two-peak structure related to the two VHS inthe two-band electronic spectrum. For large values of J > J c = 0 . t a gapless superconductivity is found at˜ µ = 0. It is suggested that for small doping, δ . .
1, the0AF long-range order found in [28] may coexist with the d + id ′ superconductivity.We have also calculated T c for the extended s -wavepairing to compare it with the results obtained for the t – J model, where the constraint has been neglected orconsidered approximately. In the latter case the resultsfor the maximum value of T c are comparable. However,we have shown that the s -wave pairing violates the con-straint and should be ruled out.In graphene with the large hopping parameter t ≈ . U/t = 4 − t – J model to graphene is questionable. In our the-ory we have the renormalized hopping parameter ˜ t (24)which is small in the region δ ≤ .
4, where the nn cor-relation function C ≤ − . U/ ˜ t ≫
1, and the application of the t – J model can bejustified. Moreover, complicated structures like sulfur-graphite composites [9] may result in a larger value of U/t , and high- T c superconductivity can be achieved. Acknowledgments
The financial support by the Heisenberg-Landauprogram of JINR is acknowledged. One of the authors(N. P.) thanks the Directorate of the MPIPKS forthe hospitality extended to him during his stay at theInstitute, where a part of the present work has been done. [1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, A. K. Geim, Rev. Mod. Phys. , 109 (2009)[2] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, A. H.Castro Neto, Rev. Mod. Phys. , 1067 (2012)[3] A. M. Black-Schaffer, Wei Wu, K. Le Hur, Phys. Rev. B , 054521 (2014)[4] B. Uchoa, G. G. Cabrera, A. H. Castro Neto, Phys. Rev.Lett. , 184509 (2005)[5] B. Uchoa, A. H. Castro Neto, Phys. Rev. Lett. , 146801(2007)[6] J. Hubbard, Proc. Roy. Soc. (London) A , 238 (1963)[7] C. Honerkamp, Phys. Rev. Lett. , 146404 (2008)[8] J. P. L. Faye, P. Sahebsara, D. S`en`echal, Phys. Rev. B , 085121 (2015)[9] A. M. Black-Schaffer, S. Doniach, Phys. Rev. B ,134512 (2007)[10] G. Baskaran, Phys. Rev. B , 212505 (2002)[11] Y. Jiang, Dao-Xin Yao, E. W. Carlson, Han-Dong Chen,Jiang Ping Hu, Phys. Rev. B , 235420 (2008)[12] S. Pathak, V. B. Shenoy, G. Baskaran, Phys. Rev. B ,085431 (2010)[13] Wei Wu, M. M. Scherer, C. Honerkamp, K. Le Hur, Phys.Rev. B , 094521 (2013)[14] S. Jiang, A. Mesaros, Y. Ran, Phys. Rev. X , 031040(2014)[15] Z.-C. Gu, H.-C. Jiang, D. N. Sheng, H. Yao, L. Balents,X.-G. Wen, Phys. Rev. B , 155112 (2013)[16] Z.-C. Gu, H.-C. Jiang, G. Baskaran, arXiv:1408.6820[17] R.S. Markiewicz, J. Phys. Chem. Solids , 1179 (1997)[18] J. L. McChesney, A. Bostwick, T. Ohta, T. Seyller, K.Horn, J. Gonzalez, E.Rotenberg, Phys. Rev. Lett. ,136803 (2010) [19] D. J. Scalapino, E. Loh, J. E. Hirsch, Phys. Rev. B ,6694 (1987)[20] J. Gonz´alez, Phys. Rev. B , 205431 (2008)[21] T. Ma, F. Yang, H. Yao, Hai-Qing Lin, Phys. Rev. B ,245114 (2014)[22] R. Nandkishore, A. V. Chubukov, Phys. Rev. B ,115426 (2012)[23] X.-Y. Xu, S. Wessel, Z.-Y. Meng, Phys. Rev. B ,115105 (2016)[24] F.-C. Zhang, C. Gros, T. M. Rice, H. Shiba, Supercond.Sci. Technol. , 36 (1988)[25] G. Martinez, P. Horsch Phys. Rev. B , 317 (1991)[26] Anderson P.W., Science , 1196 (1987)[27] J. Hubbard, Proc. Roy. Soc. A (London) , 542 (1965)[28] A. A. Vladimirov, D. Ihle, N. M. Plakida, Eur. Phys. J.B , 195 (2018)[29] H. Mori, Prog. Theor. Phys. , 399 (1965)[30] D. N. Zubarev, Usp. Phys. Nauk , 71 (1960) [trans-lated in Sov. Phys. Usp. , 320 (1960)], NonequilibriumStatistical Thermodynamics, Consultants Bureau, New-York, 1974.[31] Plakida N. M., In: Theoretical Methods for Strongly Cor-related Systems, Chap. 6, pp. 173–202. Eds. A. Avellaand F. Mancini, Springer Series in Solid–State SciencesVol. , Springer, Heidelberg, 2011.[32] N.M. Plakida, V.S. Oudovenko, Phys. Rev. B , 11949(1999)[33] N.M. Plakida, V.S. Oudovenko, Eur. Phys. J. B , 115(2013)[34] N.M. Plakida, V. Yushankhai, I.V. Stasyuk, Physica C160