Electronic spectrum and superconductivity in the $t$-$J$ model on the honeycomb lattice

A microscopic theory of electronic spectrum and superconductivity within the $t$-$J$ model on the honeycomb lattice is formulated. The Dyson equation for the normal and anomalous Green functions for the two-band model in terms of the Hubbard operators is derived by applying the Mori-type projection technique. The self-energy is evaluated in the self-consistent Born approximation for electron scattering on spin and charge fluctuations induced by the kinematical interaction for the Hubbard operators. Superconducting pairing mediated by the antiferromagnetic exchange and spin fluctuations is discussed.


Introduction
One of the crucial issues in the current theory of condensed matter is the study of superconductivity in strongly correlated electronic systems, as in the cuprate high-temperature superconductors (see, e.g., [1]). The basic model commonly used in these studies is the Hubbard model which is specified by two parameters: the single-electron hopping matrix element t and the single-site Coulomb energy U. A rigorous analytical treatment of the model is based on the Hubbard operator (HO) technique [2] since in this representation the local constraint of no double occupancy of any lattice site is rigorously implemented by the Hubbard operator algebra. A superconduction pairing due to the kinematical interaction in the Hubbard model in the limit of strong electron correlations (U → ∞) was first considered by Zaitsev and Ivanov [3] who used a diagram technique for the HOs. Subsequently, superconducting pairing in the Hubbard model was studied by Plakida and Stasyuk [6] by applying the equation of motion method for the thermodynamic Green functions (GFs) [7,8]. In the later studies [9,10] of the t-J model within the Mori-type projection technique [11,12], a formula for the superconduction temperature similar to [3] was obtained in the mean-field approximation (MFA). A self-consistent solution of Dyson equations beyond the MFA for normal and anomalous GFs and respective self-energies within the t-J model was performed in [13] that confirmed the results of the previous studies in MFA [9,10]. Detailed studies of the extended Hubbard model on the square lattice were performed in [14,15]. Taking into account an intersite Coulomb repulsion and electron-phonon interaction, a microscopic theory of superconductivity in cuprates was formulated. It was proved that the spin-fluctuation pairing mechanism induced by the strong kinematical interaction is the leading interaction which results in high-temperature superconductivity.
In recent years, the two-dimensional carbon honeycomb lattice, the graphene, has been extensively studied due its peculiar electronic properties (for a review see [16,17]). Studies of the graphene beyond the simple model of noninteracting electrons by taking into account the Coulomb interaction reveal a rich phase diagram with phase transitions to the antiferromagnetic (AF) state, spin-density wave (SDW), charge-density wave (CDW), and nonconventional superconductivity (SC). Superconducting phase transitions in the Hubbard model on the honeycomb lattice have been considered in several publications. The renormalization group approach was used in [18] to study phase transitions in the extended Hubbard model with a moderate on-site interaction U. Close to half-filling, the SDW or CDW orders occur, while for a large doping, f -wave triplet-pairing and d + id-wave singlet-pairing appear. Using the dynamic cluster approximation for the Hubbard model with U/t = 2 − 6, a transition from the d + id-wave singlet pairing at weak coupling to the p-wave triplet pairing at larger coupling was observed in [19].
In the limit of strong correlations, U ≫ t, two conduction bands of electrons in the honeycomb lattice split into the singly-and doubly-occupied Hubbard subbands. In this limit, the Hubbard model can be reduced to the two-band t-J model for the projected electron operators. A detailed study of the t-J model on the honeycomb lattice was presented in [20]. The ground-state energy and the staggered magnetization in the AF phase as function of doping δ have been calculated using the Grassmann tensor product state approach, exact diagonalization and density-matrix renormalization methods. The occurrence of the time-reversal symmetry breaking d + id-wave SC at large doping was found. Moreover, a coexisting of the SC and AF order was observed for low doping, 0 < δ < 0.1, where the triplet pairing is induced.
In the above cited papers, mostly the phase diagram of the models with the Coulomb interaction on the honeycomb lattice was studied at zero temperature by numerical methods. To investigate temperature and doping dependence of the phase diagram and, in particular, superconduction phase transition analytical methods should be used. The honeycomb Heisenberg model at half-filling over the whole temperature region both in the AF and paramagnetic phases was studied in [21]. The electronic and spin-fluctuation spectra in this model were studied in [22] using the generalized mean-field approximation (GMFA). In the present paper, we formulate a microscopic theory of superconductivity within the two-band t-J model on the honeycomb lattice. We derive Dyson equation for the normal and anomalous matrix GFs with the self-energy evaluated in the self-consistent Born approximation (SCBA) as for the one-band t-J model in [13] and the Hubbard model in [14,15] on the square lattice. Electronic spectrum and superconductivity is considered in the MFA.
In section 2, we formulate the t-J model in terms of Hubbard operators. The Dyson equation is derived in section 3. The results of calculation of the electronic spectrum and the superconduction gap function are presented in section 4. The conclusion can be found in section 5.

The t-J model
We consider the t-J on the honeycomb lattice. The honeycomb lattice is bipartite with two triangular sublattices A and B. Each of the N sites on the A sublattice is connected to three nearest-neighbor (nn) sites belonging to the B sublattice by vectors δ α , and N sites on B are connected to A by vectors −δ α : The basis vectors are a 1 = δ 3 − δ 2 = (a 0 /2)( √ 3, 3) and a 2 = δ 3 − δ 1 = (a 0 /2)(− √ 3, 3), the lattice constant is a = |a 1 | = |a 2 | = √ 3a 0 , where a 0 is the nn distance; hereafter we put a 0 = 1. The reciprocal lattice vectors are k 1 = (2π/3)( √ 3, 1) and k 2 = (2π/3)(− √ 3, 1). In conventional notation, the t-J model reads: where for the two-sublattice representation A, B of the site indices on the honeycomb lattice, a shorter notation iα → i is used. The projected electron operators in the singly occupied Hubbard subband in the model (2.2) are defined asã + i,σ = a + i,σ (1 − n i,σ ) andã iσ = a iσ (1 − n i,σ ) for creation and annihilation of an electron with spin σ/2 (σ = ±1,σ = −σ) , n i,σ =ã + i,σã i,σ , and n i = σ n i,σ . Here, t is the nn electron hopping energy and J = 4t 2 /U is the nn AF exchange interaction for electronic spins S i .
To take into account, on a rigorous basis, the no-double-occupancy constraint for the projected electron operatorsã + i,σ , we employ the Hubbard operator (HO) technique [2]. The HOs X nm i = |i, n i, m| 33706-2 determine transitions between three possible states on a lattice site i: |i, n = |i, 0 and |i, σ for an empty site and for a singly occupied site by an electron with spin σ/2, respectively. The electron number operator and the spin operators in terms of HOs are defined as The completeness relation for the HOs, X 00 i + σ X σσ i = 1, rigorously preserves the constraint of nodouble-occupancy for a quantum state |i, n on any lattice site i. From the multiplication rule X nm i X kl i = δ mk X nl i for HOs, there follow the commutation relations: The upper sign refers to Fermi-type operators such as X 0σ i , while the lower sign refers to Bose-type operators such as n i (2.3) or the spin operators (2.4).
Using the Hubbard operator representation forã + iσ = X σ0 i ,ã jσ = X 0σ j and equations (2.3) and (2.4), we write the Hamiltonian of the t-J model (2.2) in the form: The chemical potential µ depends on the average electron occupation number n = n A = n B , where N is the number of unit cells and ... denotes the statistical average with the Hamiltonian (2.6).

Dyson equation
To consider superconductivity within the model (2.6), we introduce the anticommutator two-time matrix Green function (GF) [7,8] where X(t) = e iHt Xe −iHt and θ(x) is the Heaviside function. Here, we use Nambu notation and introduce the Hubbard operators in the two-sublattice representationX iσ andX † iσ wherê The Fourier representation in (k, ω)-space is defined by The 4 × 4 matrix GF (3.1) can be written as

33706-3
where the normalĜ i jσ and anomalousF i jσ components of the GF are 2 × 2 matrices, which are coupled by the symmetry relations for the anticommutator retarded GF [7,8].
To calculate the GF (3.1), we use the equation of motion method. Differentiating the GF with respect to time t, its Fourier representation leads to the equation Here,τ 0 is the 4 × 4 unit matrix and in a paramagnetic state, the coefficient Q = X 00 iα + X σσ iα = 1 − n α /2 depends on the occupation number of electrons (2.7) only. For a system of strongly correlated electrons as in the t-J model, there is no well-defined quasiparticle (QP) excitations specified by zeroth-order kinetic energy. Therefore, it is convenient to choose the mean-field contribution to the energy of QPs in the equations of motion (3.5) as the zeroth-order QP energy. To identify this contribution, we use the projection operator method developed for the GF [12]. To this end, we write the operatorẐ iσ = [X iσ , H] in (3.5) as a sum of the linear part, proportional to the original operatorX iσ , and the irreducible partẐ The orthogonality condition iσ = 0 determines the linear part, the frequency matrix: Frequency matrix (3.7) determines QP spectrumÊ i j and the gap function∆ i jσ in the GMFA. The corresponding zeroth-order GF in the Fourier representation reads .
Differentiating the multiparticle GF Ẑ iσ (t),X † jσ (t ′ ) in (3.5) with respect to the second time t ′ and using the same projection procedure as in (3.6) leads to the Dyson equation for the GF (3.1). In the (k, ω)-representation, the Dyson equation reads The self-energy operator Σ σ (k, ω) is defined by the proper part of the scattering matrix which has no parts connected by the zeroth order GF (3.8): The self-energy operator (3.10) can be written in the same form as GF (3.4): The normalM and anomalous (pair)Φ components of the self-energy operator (3.11) are given by the 2 × 2 matrices [see (3.22), (3.23)].
The system of equations for the (4 × 4) matrix GF (3.4) and the self-energy (3.11) can be reduced to a system of equations for the normalĜ σ (k, ω) and anomalousF σ (k, ω) (2 × 2) matrix components. By using representations for the frequency matrix (3.7) and the self-energy (3.11), we derive the following system of matrix equations for them: In (3.12), we introduced the normal state matrix GF and the matrix superconduction gap function: To obtain a closed system of equations, we should evaluate the multiparticle GF in the self-energy operator (3.11) which describes the processes of inelastic scattering of electrons on charge and spin fluctuations due to the kinematic interaction.

Mean-field approximation
The superconducting pairing in the Hubbard model already occurs in the MFA and is caused by the kinetic exchange interaction as proposed by Anderson [24,25]. It is, therefore, reasonable to consider the MFA described by zeroth-order GF (3.8) separately. Using commutation relations for Hubbard operators (2.5), we evaluate the frequency matrix (3.7). The matrixÊ(k) determines the QP spectrum in the normal phase (see [22]): The solution of the matrix equation for the zero-order GF (3.8) in the normal state reads: The electronic spectrum has two branches: where we take into account the fact that in the paramagnetic state, the sublattices are equivalent and, therefore, ε A (k) = ε B (k) = ε(k). We calculate the spectrum and consider the Fermi surface in the normal state in the next section. Now, we evaluate the anomalous component∆ i jσ of the matrix (3.7), which determines the superconducting gap. Similar to the normal-state frequency matrix (3.16), we obtain the representation Diagonalization of the gap matrix (3.20) shows a two-gap structure: where we take into account the fact that in the paramagnetic state, the sublattices are equivalent and, therefore, ∆ Aσ (k) = ∆ Bσ (k) = ∆ σ (k). We consider the gap components (3.20) in the next section.

Self-energy operator
The normalM and anomalous (pair)Φ components of the self-energy operator (3.11) are given by the 2 × 2 matrices:M where [ X n,m iα , H] are the irreducible parts of the commutators determined by equation (3.6). Taking into account equations for the [X 0σ iα , H] operators, we obtain the following representation for the matrix elements of the normal self-energy To calculate the self-energy matrix elements in (3.24), (3.25), and (3.27), (3.28), we use the SCBA for the multiparticle GFs. For this, it is convenient to write the self-energy in terms of time-dependent correlation functions. In particular, for the normal component (3.24), we have the spectral representation: In the SCBA, a propagation of excitations described by the Fermi-like operatorsX σ0 i and the Bose-like operators B σ ′ σ j for i j is assumed to be independent. Therefore, the corresponding time-dependent multiparticle correlation functions can be written as a product of fermionic and bosonic correlation functions, (3.30)

33706-6
The time-dependent correlation functions in (3.30) are calculated self-consistently using the corresponding GFs, as e.g., where n(ω) and N(ω) are the Fermi-and Bose-function, respectively. Integration over time t in (3.29) yields Taking into account the definition of the bosonic operator (3.26), the bosonic GFs in these equations can be written as After summation over σ ′ in equation (3.33) for the bosonic GF (3.34) and taking into account that the normal GF in the paramagnetic state is independent of spin, the spin-fluctuation contribution can be written in the form: Introducing the k-representation for the GFs in the self-energy (3.33), we obtain the expression: The SCBA for M i jAB (ω) (3.25) and for the anomalous self-energies (3.27), (3.28) gives similar results.

Results
To find out the normal state GF (3.14) beyond the MFA, the self-energy matrix elements (3.24), (3.25) should be calculated. Similar calculations of the anomalous self-energy (3.27), (3.28) give an equation for the frequency dependent gap (3.15). However, implementation of this program demands complicated self-consistent numerical computations for the two-band t-J model which are much more complicated than similar calculations performed for the one-band t-J model on the square lattice in [13]. Therefore, in the present paper, we restrict our consideration only to the results in MFA and study the complete system of equations later.

Electronic spectrum
Calculation of the matrix elements (3.16) gives the following results: where we introduce the nn correlation functions for electrons and spins, For the off-diagonal energy, we have: where γ 1 (k) = α exp(ikδ α ) and |γ 1 (k)| 2 = 1 + 4 cos( √ 3k x /2)[cos( √ 3k x /2) + cos(3k y /2)]. Thus, the electronic spectrum has two branches ε ± (k) = −µ ± t |γ 1 (k)| , (4.4) which is similar to the spectrum of the graphene (see, e.g., reference [16]) but with the renormalized chemical potential µ due to strong correlations, equation (4.1), and the hopping parameter t , equation (4.3). Therefore, the spectrum shows a cone-type behaviour at Dirac points K and K ′ at the corners of the graphene Brillouin zone (BZ) as shown in figure 1. In figure 2, the BZ and the Fermi surfaces (FS) are shown for holes at the electronic occupation numbers n = 0.95, 0.76, and 0.7. At n 1, the hole FS is small and centered at the Γ point. With a decreasing n, the FS becomes larger, and at some characteristic value n 0 = 0.76, the FS touches the BZ at M-points. At a larger hole doping, six pockets centered at the K-points emerge which shrink to points for the half-filled band at n = 2/3. For the diagonal GF, we have The average occupation number of electrons in one sublattice is equal to which yields in MFA where N(ε ± (k)) = {exp[ε ± (k)/T] + 1} −1 , and the average number of electrons in one sublattice is given by with n α 1. For the off-diagonal GF in equation (3.10), we obtain For the nn correlation function, we get (4.10) These results can be used in calculating the self-energy in the normal state (3.36).

Gap equation
The matrix elements of the gap function (3.20) are given by Using the commutation relations for the HOs, for these functions we obtain the following representation: where the symmetry relations were used: X 0σ jB X 0σ iA = − X 0σ iA X 0σ jB and X 0σ jB X 0σ iA = − X 0σ jB X 0σ iA .

33706-9
Fourier transformation results in the gap function components: where δ α = δ 1 , δ 2 , δ 3 . Introducing the bond-dependent anomalous correlation functions the gaps can be written as Here, the nn correlation function (4.16) is calculated using the spectral representation: For the bond-independent correlation functions (s-wave pairing), F α ABσ = F ABσ , we have the results similar to the normal state energy: Using the GMFA for the anomalous GF in (4.18), we can derive the Bardeen-Cooper-Schrieffer (BCS) type equations as was proposed in several publications for electrons with the Dirac spectrum. In particular, in [23], the s-wave BCS model was considered with a phenomenological coupling constant. In the case of the two-subband honeycomb lattice, we expect to obtain a more complicated pairing symmetry, e.g., d + id singlet pairing or p triplet pairing as proposed in [20]. In the case of the t-J on the square lattice in [13], we have also found a more complicated d-wave pairing induced by spin-fluctuations. A self-consistent solution of the gap equations (4.17) with microscopic coupling constants and anomalous GF (3.13) and a resulting equation for determining the superconduction transition temperature will be considered in a subsequent publication.

Conclusions
In the present paper, a microscopic theory of the electron spectrum and superconductivity within the two-band t-J model for strongly correlated electrons on the honeycomb lattice is presented. Using the projection operator method for thermodynamic GFs, we derived a self-consistent system of equations for the matrix GF and the self-energy in the SCBA. The latter is similar to the Migdal-Eliasberg strongcoupling theory for the electron-phonon system. In the present paper, we consider only the results obtained in the MFA. However, as it was pointed out in the microscopic theory of superconductivity formulated for the Hubbard model on the square lattice in [14,15], consideration of the self-energy effects is very important. The normal self-energy component determines a reduction of the QP weight which can be quite large. The anomalous (pair) self-energy component strongly enhances the superconduction pairing induced by the kinematical interaction of electrons with spin-fluctuations. This interaction being proportional to the hopping parameter t [see equation (3.36)] is an order of magnitude larger than the exchange interaction J considered in the Anderson's theory of the resonating valence bonds [24,25]. We are planning to take into account the self-energy effects in a subsequent publication.