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.


I. INTRODUCTION
In recent years the two-dimensional carbon honeycomb lattice, the graphene, has been extensively studied due 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) states, 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 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 crit- * E-mail: plakida@theor.jinr.ru ical point in the universality class of the Gross-Neveu model was found. A general low-energy 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 selfconsistent approach for the Hubbard model on the honeycomb lattice in reference [13] shows the semimetal to spin-liquid 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 spinliquid 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 singlesite 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 superconducting 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 spin-polaron 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 groundstate 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-rotationinvariant relaxation-function 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 done for the DSS. Similar studies have been done for the t-J model on the square lattice in our paper [26].
In Section II we formulate the t-J model in terms of Hubbard operators. The electronic excitation spectrum is calculated in Section III. The spin-excitation spectrum and thermodynamic quantities are considered in Section IV. The numerical results and discussion are given in Section V. The conclusion can be found in Section VI.

II. 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 n i = σ n i,σ . 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, equa- 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 where a 0 is the nn distance (see Figure 1); hereafter we put a 0 = 1. The reciprocal lattice vectors are k 1 = (2π/3)( √ 3, 1) and . In the two-sublattice representation it is convenient to split the site indices into the unit cell and sublattice indices, i → iα, α = A, B.
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).

III. 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 {X, Y } = XY + Y X, X(t) = e iHt Xe −iHt , and θ(x) 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: {Ẑ 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: where 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). For the diagonal GF we have The mean occupation number of electrons is equal to with n ασ (k) = X σσ kα = X σ0 kα X 0σ where N (ε ± (k)) = [exp[ε ± (k)/T ] + 1] −1 , and n α ≤ 1.
For the off-diagonal GF in (18) we obtain and for the corresponding correlation function we get

IV. 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-energŷ Σ(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 zeroorder 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), we start from the second derivative −S + i . Taking into account only the diagonal contributions is the hopping (exchange) part of the model (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 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 conribution 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 matrixF (q):F where 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)|).

V. 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 self-consistency 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.
In Figure 3 our results for the doping dependence of the spin correlation functions C n up to n = 4 at T = 0 are presented. The different sign of C n reflects the AF order which gradually decreases with increasing doping, where at δ 0.6 only nn spin correlations survive. Considering the staggered magnetization m st (δ) at T = 0 depicted in Figure 4, 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.
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 we get δ hl c = 0.069 and δ sl c = 0.038 (see [26]), i.e., the ratio δ hl c /δ sl c ≃ 1.8 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.
In Figure 5 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. Figure 3), 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 Figure 3). 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 5 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 6 show the inverse correlation length ξ −1 . Considering the doping dependence ( Figure 6(a) 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 Figure 3). Let us consider the temperature dependence of ξ ( Figure 6(b). At δ < δ c , ξ diverges in the limit T → 0, where at δ = 0, ξ exhibits the known exponential increase (see [24]).
The spectrum of spin excitations ω ± (q) at T = 0 is shown in Figure 7 along the symmetry directions X(−1, 0) → K(−2/3, 0) → Γ(0, 0) → Y (0, 1) → K(1/3, 1) → Γ ′ (1, 1) → M (1/2, 1/2) → Γ of the BZ shown in Figure 2. In the LRO phase, i.e., at δ < δ c , the spin excitations are spin waves with gapless branches depicted in Figure 7(a). 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, 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 7, the spin-excitation energies are decreasing with increasing doping. Interestingly, at the K points the spin-excitation spectrum has a maximum.

VI. CONCLUSION
In this paper we have evaluated thermodynamic quantities as well as the electronic and spin-excitation spectra in the strong-correlation limit within the t − J model on the honeycomb lattice. 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 and of the uniform static spin susceptibility which we have explained in terms of AF SRO. As compared with previous results for the t − J model on the square lattice, both the AF LRO and SRO are found to exist in a larger doping region. We conclude that our investigation forms a good basis for forthcoming studies of superconductivity in the honeycomb t − J model.