Spin excitations and thermodynamics of the t-J model on the honeycomb lattice

We present a spin-rotation-invariant Green-function theory for the dynamic spin susceptibility in the spin-1/2 antiferromagnetic t-J Heisenberg model on the honeycomb lattice. Employing a generalized mean-field approximation for arbitrary temperatures and hole dopings, the electronic spectrum of excitations, the spin-excitation spectrum and thermodynamic quantities (two-spin correlation functions, staggered magnetization, magnetic susceptibility, correlation length) are calculated by solving a coupled system of self-consistency equations for the correlation functions. The temperature and doping dependence of the magnetic (uniform static) susceptibility is ascribed to antiferromagnetic short-range order. Our results on the doping dependencies of the magnetization and susceptibility are analyzed in comparison with previous results for the t-J model on the square lattice.


Introduction
In recent years the two-dimensional carbon honeycomb lattice, the graphene, has been extensively studied due to its peculiar electronic properties (for a review see [1,2]). Studies of the graphene beyond the simple model of noninteracting electrons by taking into account the Coulomb interaction (CI) reveal a rich phase diagram with phase transitions to the antiferromagnetic (AF) state, spindensity wave (SDW), charge-density wave (CDW), and nonconventional superconductivity (SC).
In many papers the electronic properties of the Hubbard model on the honeycomb lattice were investigated. It was found that at a sufficiently large single-site Coulomb repulsion U > U c ≈ 4t, the AF long-range order (LRO) emerges for a single layer close to half-filling [3][4][5]. The phase diagram and spin excitations of the Hubbard model for graphene layers using the mean-field approximation (MFA) and the random-phase approximation were considered in reference [6]. Depending on the value of U and electronic density n, various phases were observed: at large U > U c ≈ 3.8t, the AF phase for n 1 was found, while at larger doping the ferromagnetic and spiral phases were obtained. The quantum phase transition in the half-filled Hubbard model on the honeycomb lattice at U > U c with U c /t ≈ 4 − 5 was found in [7] using the quantum Monte Carlo (QMC) and series expansion techniques. The temperature dependence of the specific heat also points to the AF phase transition at U > U c . In reference [8] phase a e-mail: plakida@theor.jinr.ru transitions in the Hubbard model of N -flavor electrons on the honeycomb lattice have been discussed in the limit of large N . There, a semimetal to AF insulator phase transition at the quantum critical point in the universality class of the Gross-Neveu model was found. A general lowenergy theory of electrons with repulsive short-range CI on the honeycomb lattice at half-filling is presented in [9].
The phase diagram of extended Hubbard models with nearest-neighbor (nn) and next nearest-neighbor (nnn) repulsive interactions V 1 and V 2 , respectively, on the honeycomb lattice in MFA was obtained in reference [10]. A phase transition from the semimetal to Mott insulating phases at half-filling was found at large U > U c ≈ 3.8t. For small V 1 and V 2 the AF phase appears, while for larger V 1 and V 2 the renormalization group (RG) analysis shows transitions to the SDW or CDW.
In a more recent QMC calculation [11] a gapped AF phase at half-filling for U/t > 4.3 was found, and for the intermediate coupling 3.5 < U/t < 4.3, an insulating gapped spin-liquid state formed by short-range resonating valence bonds was predicted. But later QMC studies of larger clusters have not confirmed this transition to the spin-liquid state [12]. The two-particle self-consistent approach for the Hubbard model on the honeycomb lattice in reference [13] shows the semimetal to spinliquid transition before the transition to the AF state. In reference [14], effective spin models for the Hubbard model on the honeycomb lattice at half-filling were derived. It was observed that the six-spin interactions frustrate the AF order and may lead the spin-liquid state behavior. But the spin-liquid state has not been found in other publications. The transition from the weak-coupling semimetal to the strong-coupling insulating phase was studied in [15] using QMC simulations for the SU(N )symmetric Heisenberg model with the nn flavor exchange interaction on the honeycomb lattice at half filling. In the SU(2) case a direct transition between the semimetal and an AF insulator was obtained. In reference [16] a continuous quantum phase transition between the semimetallic and the insulating AF states was found at U c /t = 3.78 by considering a staggered magnetic moment in the local magnetic field. A direct transition from a Dirac semimetal to an AF Mott insulator was confirmed in reference [17] by using the projective auxiliary-field QMC simulations and a finite-size scaling analysis. Although the existence of the spin-liquid state in the Hubbard model on the honeycomb lattice is still under discussion, the transition from the semimetal to the AF LRO phase is proved for a large enough single-site CI, U > U c ≈ 4t.
Superconducting phase transitions in the Hubbard model on the honeycomb lattice have been considered in several publications. The RG approach was used in [18] to study phase transitions in the extended Hubbard model with the on-site interaction U , the nn intersite repulsion V , and the spin-exchange interaction J. Close to half-filling, the SDW or CDW orders occur for large U and V , while for a large doping f -wave triplet-pairing and d + id-wave singlet-pairing emerge. Chiral triplet superconductivity on the graphene lattice was considered in [19]. 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 [20]. More references on studies of superconduction phase transitions in the Hubbard model on the honeycomb lattice may be found in reference [20].
In the limit of strong correlations, U t, the conduction band of electrons on the honeycomb lattice splits into the singly-and doubly-occupied Hubbard subbands. In this limit the Hubbard model can be reduced to the t-J model for the projected electron operators in one subband. This model was investigated by several authors. In [21] a single-hole excitation was considered within the t-J spinpolaron model. The results obtained for the honeycomb lattice are qualitatively similar to those for the square lattice. A detailed study of the t-J model on the honeycomb lattice was presented in [22]. 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 (see also [23]).
In the papers cited above mostly the phase diagram of the correlation models on the honeycomb lattice at zero temperature was studied. Less attention has been paid to the investigation of electron-and spin-excitation spectra and of thermodynamic properties as functions of temperature and electron concentration. Motivated by this shortcoming, in the present paper we report results of investigations of these spectra and of the thermodynamics in the limit of strong correlations within the t-J model. In our previous paper [24] we have studied the honeycomb Heisenberg model at half-filling over the whole temperature region both in the AF and paramagnetic phases. Thereby, we have calculated the dynamic spin susceptibility (DSS) within the spin-rotation-invariant relaxationfunction theory [25][26][27] using the generalized mean-field approximation (GMFA). Let us point out that the GMFA has been successfully applied to several quantum spin systems (see, e.g., [24] and references therein). In the present paper we consider the effects of doping on AF order within the GMFA for the DSS. Similar studies have been done for the t-J model on the square lattice in our paper [26].
In Section 2 we formulate the t-J model in terms of Hubbard operators. The electronic excitation spectrum is calculated in Section 3. The spin-excitation spectrum and thermodynamic quantities are considered in Section 4. The numerical results and discussion are given in Section 5. The conclusion can be found in Section 6.

The t-J model
We study the Hubbard model on the honeycomb lattice in the limit of strong electron correlations U >> t, when it can be reduced to the one-subband t-J model: are projected creation and annihilation electron operators with spin σ/2 (σ = ±1,σ = −σ) in the singly occupied Hubbard subband, n i,σ =ã + i,σã i,σ . Here, t is the nn electron hopping energy.
The Heisenberg Hamiltonian in (1) is given by where J = 4t 2 /U is the nn AF exchange interaction and To take into account on a rigorous basis the projected character of electron operatorsã + i,σ , we employ the Hubbard operator (HO) technique [28]. The HOs are defined as for three possible states at 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 no double occupancy of the quantum state |i, n on any lattice site i. From the multiplication rule X nm i X kl i = δ mk X nl i 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 (4) or the spin operators (5).
Using the Hubbard operator representation, j and equations (4) and (5), we write the Hamiltonian of the t − J model (1) in the form: The Hamiltonian has the conventional form of the t-J model in terms of Hubbard operators (see, e.g., [29]). We consider the honeycomb lattice shown in Figure 1. The lattice is bipartite with two triangular sublattices A and B. Each of the N sites on the A sublattice is connected to three nn sites belonging to the B sublattice by vectors δ j , and N sites on B are connected to A by vectors −δ j : The basis vectors are a 1 = δ 3 − δ 2 = (a 0 /2)( √ 3, 3) and where a 0 is the nn distance (see Fig. 1); 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 the two-sublattice representation it is convenient to split the site indices into the unit cell and sublattice indices, The chemical potential µ depends on the average electron occupation number where N is the number of unit cells and . . . denotes the statistical average with the Hamiltonian (7).

Electronic excitation spectrum
To calculate the electron excitation spectrum within the model (7), we consider the anticommutator two-time matrix Green function (GF) [30]  where is the Heaviside function. Here we introduce the Hubbard operators in the two-sublattice representation: The Fourier representation in (k, ω)-space is defined bŷ Differentiating the GF (10) with respect to the time t we get , τ 0 is the unity matrix, and Q = X 00 iα + X σσ iα = 1 − n α /2. Now, we project the many-particle GF in (13) on the single-electron GF by introducing the irreducible part of theẐ 0σ i operator, which is orthogonal to the right-hand side operator: {Ẑ 0σ(irr) i ,X σ0 j } = 0. This results in the equation for the frequency matrix, Using the Fourier transformation of the GF (12) we obtain the equation for the GF in the GMFA neglecting the last term in (14) which describes inelastic scattering: Here the electronic excitation spectrum in GMFA is determined by the matrix of correlation functions: where kA , H], X σ0 kB } . The solution of the matrix equation for the GF (16) reads: The electronic spectrum is defined from the equation and is given by The calculation of the matrix elements in (17) gives the following result for ε(k): where we introduce the nn correlation functions for electrons and spins, For the off-diagonal energy we have: where Note that equations (21), (23), (24) are similar to the results obtained for the spectrum of the t-J model on the square lattice in [29].
Thus, the electronic spectrum has two branches: It agrees with the spectrum of graphene (see, e.g., [1,31]), except for the renormalization of the chemical potential µ and the hopping parameter t due to strong correlations. Therefore, in the strong-correlation limit the cone-type dispersion is conserved, i.e., the spectrum reveals Dirac cones at the corners (K points) of the Brillouin zone (BZ). The BZ and the Fermi surfaces (FS) for holes at the electronic occupation numbers n = 0.95, 0.76, and 0.7 are shown in Figure 2. At n 1 the hole FS is small and centered at the Γ point. With decreasing n the FS becomes larger, and at some characteristic value n 0 = 0.76 the FS touches the BZ at M-points. At 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 off-diagonal GF in (18) we obtain and for the corresponding correlation function we get

Spin-excitation spectrum and thermodynamics
To calculate the spin-excitation spectrum and to evaluate the thermodynamic quantities in the model (7), we consider the two-time matrix commutator GF [30]: where [S, Y ] = SY − Y S. The Fourier representation of the spin GF, S + q | S − −q ω , is defined by the same relation as for the electronic GF (12).
In the relaxation-function theory developed on the basis of the equation of motion method in [25,26] we obtain the following representation of the DSSχ(q, ω) = −Ĝ +− (q, ω): Here,F (q) is the frequency matrix of spin excitations in the GMFA, where the approximation The self-energyΣ(q, ω) can be expressed exactly by a multispin GF (see [25,26]).
We consider the GMFA for the DSS neglecting the selfenergy in (32). Then, for a lattice with basis the zero-order DSS is given by (cf. [24]) For the static spin susceptibility we obtain The direct calculation of the matrix elements m αβ (q) yieldsm where To calculate the frequency matrixF (q) in equation (33) (7), we obtain where We do not consider the off-diagonal terms F tJ i , F Jt i as discussed in [26].
We perform the following decoupling procedure preserving the local correlations. Decoupling the operators in H σ ijn we introduce the parameter λ: where for n = i, X σ0 i X 0σ i = X σσ i = n/2 , and the second term of H σ ijn with n = i is neglected (cf. [32]). Analogously, the operators in Π ijn=i and P ij,n=i are decoupled as where δ = X 00 i = 1 − n is the hole concentration, and we used the equations: ) . The parameter λ describes the renormalization of the vertex for spin scattering on charge fluctuations.
The contribution F JJ i is proportional to products of three spin operators on different lattice sites along nn sequences, e.g., iA, jB, kA . We perform the decoupling of them as follows Here, the vertex renormalization parameters α 1 and α 2 are attached to the nn and the nnn correlation functions C 1 and C 2 , respectively, and the equality S z i S z j = 1 2 S + i S − j due to spin-rotation invariance is taken into account. Using these decouplings we obtain the frequency matrix F (q):F where F AA (q) = J 2 (3(1 − λδ) + 12α 2 C 2 + 2γ 2 (q)α 1 C 1 ) with and the nnn correlation function D 2 = X σ0 iα X 0σ i+a1,α . Since the matricesm(q) andF (q) commute, it is convenient to solve equation (33) by introducing the eigenvalues m ± (q) and the normalized eigenvectors |E ± (q) of the matrix (35): which are given by For the same eigenvectors |E ± (q) the spin-excitation frequencies are obtained as In this notation the DSS reads where and α|E ± E ± |β = 1/2 for α = β, otherwise α|E ± E ± |β = ∓γ 1 (q)/(2|γ 1 (q)|).

Results
To evaluate the spin-excitation spectrum and the thermodynamic properties, the correlation functions C 1 , C 2 , the transfer amplitudes D 1 , D 2 , and the vertex parameters α 1 , α 2 , and λ, appearing in the spectrum ω ± (q) as well as the condensation term C in the LRO phase have to be determined as solution of a coupled system of selfconsistency equations. Besides equations (59) and (30) for calculating the correlation functions and the transfer amplitudes, respectively, we have the sum rule C r=0,αα = S + iα S − iα = (1 − δ)/2 and the LRO condition ω + (Q) = 0. That is, we have more parameters than equations. To obtain a closed system of equations, we take the following choice of the vertex parameters.

Spin excitations
Let us first present the spectrum of spin excitations ω ± (q) at T = 0 calculated from equation (56). The spectrum is plotted in Figure 3 along the symmetry direc- Figure 2. In the LRO phase, i.e., at δ < δ c , the spin excitations are spin waves with gapless branches depicted in Figure 3a. In the paramagnetic phase, spin waves propagating in AF SRO can exist, if their wavelength is smaller than the correlation length, i.e., if q > q c = 2πξ −1 . For dopings slightly above δ c , where the correlation length is large enough as discussed below, this condition can be fulfilled. With increasing doping, we may have q < q c so that the spin-wave picture breaks down, and "paramagnon" excitations with the energies ω ± (q) appear. Thus, our spin-excitation spectra may reveal a smooth crossover from spin-wave to paramagnon behavior depending on the wavenumber and doping. In the upper (optical) branch, at δ c a gap is opening at the Γ point, i.e., at the AF wave vector Q = (0, 0) characterizing the LRO phase in the two-sublattice model. As can be seen in Figure 3, the spin-excitation energies are decreasing with increasing doping. Interestingly, at the K points the spin-excitation spectrum has a maximum.

Correlation functions and magnetization
In Figure 4 our results for the doping dependence of the spin correlation functions for the first (C 1 ), the second (C 2 ), and the third (C 3 , C 4 ) neighbors at T = 0 are presented. The different sign of C 1 and C 2 reflects the AF order which gradually decreases with increasing doping, where at δ 0.6 only nn spin correlations survive. The third neighbor correlation functions C 3 and C 4 have similar behavior because they connect spins on the same sublattices A, B and have close distances r between the spins: C 3 (r = 2 a 0 ), C 4 (r = a 0 √ 7 ≈ 2.65 a 0 ). Considering the staggered magnetization m st (δ) at T = 0 depicted in Figure 5, the AF LRO is suppressed with increasing doping due to the spin-hole interaction.  At the critical doping δ c (J/t) we obtain a smooth phase transition from the LRO phase to a paramagnetic phase with AF SRO, where δ c (J/t = 1/3) = 0.069 and δ c (J/t = 1/2) = 0.117. For J/t = 1/3 the critical doping lies near the value δ c 0.1 found in [22] by density-matrix renormalization group calculations and by a variational method. For comparison, the staggered magnetization for the square lattice for J/t = 1/3 is also shown in Figure 5 by a dotted line, where δ c = 0.042 calculated as described in [26].
Let us compare the results for δ c at J/t = 1/3 in the honeycomb lattice (δ hl c ) with those found in the square lattice (δ sl c ). Taking δ hl c from the approach of [22], δ hl c 0.1, and δ sl c from the cumulant approach of [37], δ sl c 0.045, we have δ hl c /δ sl c 2.2. In the GMFA calculations at  J/t = 1/3 we get δ hl c = 0.069 and δ sl c = 0.042, i.e., the ratio δ hl c /δ sl c 1.64 nearly agrees with that given by the approaches of [22] and [37]. That means, in the honeycomb lattice the LRO is favored as compared with the square lattice, although the coordination number z hl = 3 is smaller than z sl = 4. This behavior may be due to different geometries of the hopping paths in the two lattices.

Static susceptibility and correlation length
In Figure 6 the uniform static spin susceptibility χ is plotted as a function of doping at various temperatures. The increase in χ upon doping is caused by the decrease of AF SRO (cf. Fig. 4), i.e., of the spin stiffness against the orientation along a homogeneous external magnetic field. At large doping, δ 0.6, χ decreases with increasing δ due to the decreasing number of spins. Note that the position of the SRO-induced maximum of χ at δ max (T ) nearly agrees with the doping at which the further-distance correlation functions C n (n = 2, 3, 4) become vanishingly small (see Fig. 4). As compared with the results for the t − J model on the square lattice, the values of δ max (T ) are much higher. That means, similar to the behavior of LRO, in the honeycomb lattice the SRO is favored as compared with the square lattice.
Concerning the temperature dependence of χ at fixed doping, from Figure 6 it can be seen that there appears a maximum at T max (δ). This maximum can be understood as a SRO effect in analogy to the explanation of the doping dependence of χ. Figure 7 shows the inverse correlation length ξ −1 . Considering the doping dependence (Fig. 7a) at T = 0, in the limit δ → δ c +, AF LRO emerges which is connected with the closing of the AF gap, ω + (Q) → 0, and by (63), with the divergence of ξ. At T > 0, there is no LRO so that ξ remains finite, where ξ −1 almost linearly increases with δ. This corresponds to the weakening of AF correlations (see Fig. 4). Let us consider the temperature dependence of ξ (Fig. 7b). At δ < δ c , ξ diverges in the limit T → 0, where at δ = 0, ξ exhibits the known exponential increase (see [24]).

Conclusion
In the paper we have calculated the electronic and spinexcitation spectra and evaluated thermodynamic quantities in the strong-correlation limit within the t-J model on the honeycomb lattice. The electronic spectrum shows the graphene-type dispersion with the renormalized chemical potential and the hopping parameter. The Fermi surface has 6 cones at K-points in the BZ at half-filling for the electronic occupation number n = 2/3. The dynamic spin susceptibility is calculated within a spin-rotation-invariant generalized mean-field approach for arbitrary temperatures and hole dopings. Our main focus was the analysis of the doping dependence of the zero-temperature magnetization, of the uniform static spin susceptibility and of the AF correlation length which we have explained in terms of SRO. The obtained results agree with our previous results for the t-J model on the square lattice [26], but both the AF LRO and SRO are found to exist in a larger doping region.
Using the obtained data for the electronic and spinexcitation spectra we can consider superconductivity as has been done for the square lattice in [29]. In particular, we expect a high superconduction temperature T c at n = 0.76 when the chemical potential occurs close to the van Hove singularity. The spin-excitation spectrum found in the paper can be used to consider the spin-fluctuation mechanism of pairing induced by strong kinematical interaction of electrons beyond the commonly used MFA for the exchange interaction pairing. As shown for the Hubbard model for U t on the square lattice [38,39], the spin-fluctuation mechanism of pairing results in high-T c . Of particular interest is the study of a coexistence of SC and AF ground-state order at low doping found in [22]. This issue may be further elaborated for finite temperatures within the layered honeycomb t-J model. So our investigation forms a good basis for further analytical studies of superconductivity in the honeycomb t-J model.