Spin excitations and thermodynamics of the antiferromagnetic Heisenberg model on the layered honeycomb lattice

We present a spin-rotation-invariant Green-function theory for the dynamic spin susceptibility in the spin-1/2 antiferromagnetic Heisenberg model on a stacked honeycomb lattice. Employing a generalized mean-field approximation for arbitrary temperatures, the thermodynamic quantities (two-spin correlation functions, internal energy, magnetic susceptibility, staggered magnetization, Néel temperature, correlation length) and the spin-excitation spectrum are calculated by solving a coupled system of self-consistency equations for the correlation functions. The temperature dependence of the magnetic (uniform static) susceptibility is ascribed to antiferromagnetic short-range order. The Néel temperature is calculated for arbitrary interlayer couplings. Our results are in a good agreement with numerical computations for finite clusters and with available experimental data on the β-Cu2V2O2 compound.


Introduction
In recent years the low-dimensional quantum Heisenberg antiferromagnets have been extensively studied motivated by experimental research of such systems. In particular, there is a vast literature devoted to the study of the two-dimensional (2D) antiferromagnetic Heisenberg model (AFHM) on the square lattice initiated by the discovery of AF long-range order (LRO) in high-temperature superconductors (for a review see [1]). According to the Mermin-Wagner theorem [2], in 2D isotropic Heisenberg magnets quantum spin fluctuations destroy the magnetic LRO at finite temperature. However, it has been shown that for the spin-1/2 2D AFHM with nearest-neighbor (nn) interaction on the square lattice, LRO can occur at zero temperature, while at finite temperature only the exponential increase of the AF correlation length with decreasing temperature is observed [3]. For the 2D honeycomb lattice with the lower coordination number z = 3, quantum spin fluctuations are more intensive and detrimental for the occurrence of LRO. Studies of the frustrated Heisenberg model with the AF interaction between the second-nearest (J 2 ) and the third-nearest (J 3 ) neighbors on the honeycomb lattice have revealed several phases with LRO, and a more complicated phase diagram occurs. For instance, using the coupled cluster method [4,5], in the J 1 -J 2 -J 3 Heisenberg model on the 2D honeycomb lattice four competing magnetic phases a e-mail: plakida@theor.jinr.ru (Néel, stripe, Néel-II collinear AF, and spiral phases) were found depending on the model parameters. Much attention has been paid to studies of the 2D Kitaev-Heisenberg model with the isotropic nn interaction J and the bonddepending Kitaev interaction K α [6]. Besides the spinliquid Kitaev phase at J = 0, a rich phase diagram at zero temperature with competing LRO (ferromagnetic, AF, stripe and zigzag phases) emerges (see, e.g., [7,8]). In a model with anisotropic K α interactions, e.g., K z > K x = K y , the LRO survives at finite temperature as shown in reference [9].
For the isotropic 2D honeycomb Heisenberg model with nn AF interaction, the LRO at zero temperature, similar to the square lattice, was confirmed in a number of studies mostly performed for finite-lattice systems. The ground-state energy E 0 , staggered magnetization m st , uniform static susceptibility χ(0) at T = 0, and other properties of the 2D honeycomb AFHM were calculated using various methods, such as extrapolations of finite-lattice diagonalizations and quantum Monte Carlo (QMC) simulations [10][11][12][13], spin-wave theory [13,14], series expansions around the Ising limit [15,16], Schwinger boson method [17], and variational RVB wave function approach [18]. In reference [19] it was suggested to consider the β-Cu 2 V 2 O 2 compound as the best available experimental realization of the spin-1/2 AFHM on the honeycomb-like lattice. The system can be characterized by the nearly isotropic nn exchange interaction J = (60-66) K and the interlayer coupling J z = 0.2J which results in the Néel temperature T N 26 K.
Besides the honeycomb AFHM, the honeycomb Hubbard model has been studied too. In reference [20] the AF LRO was found for a single layer close to half-filling at large Coulomb repulsion U , where the staggered magnetization m st = 0.335 was obtained. In reference [21], using the two-particle self-consistent approach for the Hubbard model on the honeycomb lattice, the semimetal to spinliquid transition was found before the transition to the AF state.
The AFHM on the honeycomb lattice has been less well studied as compared with the same model on the square lattice. Most of numerical computations have been performed for 2D finite-lattice systems at zero temperature, where such computations for 3D systems are difficult to realize. Moreover, the thermodynamics at arbitrary interlayer coupling, e.g., the dependence of T N on J z , is not yet developed. Therefore, analytical approaches which are capable to evaluate the thermodynamics of the AFHM on the layered honeycomb lattice both in the AF phase and in the paramagnetic phase with a temperature-dependent AF short-range order (SRO) are desirable.
To this end, in this paper we present a theory of magnetic order in the honeycomb AFHM over the whole temperature region that is based on the calculation of the dynamic spin susceptibility (DSS) within the spin-rotation-invariant (SRI) relaxation-function theory [22][23][24] using the generalized mean-field approximation (GMFA), as has been done in our study of the compass-Heisenberg model on the square lattice [25]. Using the result for the DSS, the staggered magnetization, the static spin susceptibility, the Néel temperature as a function of the interlayer coupling, the AF correlation length, and the spin-excitation dispersion both in the AF phase and in the paramagnetic phase are calculated selfconsistently, where we pay particular attention to a proper description of AF SRO. It should be pointed out that the GMFA has been successfully applied to several quantum spin systems (see, e.g., Refs. ). In particular, let us mention an application of the GMFA to study the spin-1/2 AFHM on the stacked kagomé lattice [46], where no LRO was found even for a strong interlayer coupling due to the frustration character of the AF interaction on this lattice.
In Section 2 we formulate the model and give equations for the DSS. The solution of the self-consistency equations for the spin correlation functions and the spin-excitation spectrum in the GMFA is presented in Section 3. The numerical results and discussion are given in Section 4. The conclusion can be found in Section 5.

Model and dynamic spin susceptibility
We consider congruently stacked honeycomb layers shown in Figure 1. The lattice is bipartite with two triangular sublattices A and B. Each site on the A sublattice is connected to three nn sites belonging to the B sublattice by vectors δ j , and sites on B are connected to A by vectors −δ j : The basis vectors in the layer are where a 0 is the nn distance (see Fig. 1); hereafter we put a 0 = 1. The reciprocal lattice in the layer is defined by the vectors ). The Heisenberg model on this layered honeycomb lattice can be written as where i, j count the unit cells, α and β are the sublattice indexes, iα, jβ xy and iα, jβ z denote nn sites in the xy plane and along the z direction, respectively, J > 0 is the AF intralayer exchange interaction, and J z is the coupling between the layers. To calculate the spin-excitation spectrum and to evaluate the thermodynamic quantities in the model (2), we consider the retarded two-time commutator Green function (GF) [48]: , and θ(x) is the Heaviside function. In the SRI theory all GF components ν = x, y, z are equivalent; therefore, we consider only one of them, e.g., ν = x. The Fourier representation in (q, ω)-space is defined by the relation: where N and R i is the number and position of unit cells, respectively, and the GF matrix reads In the relaxation-function theory developed on the basis of the equation of motion method in references [22][23][24] we obtain the following representation of the DSS χ(q, ω) = −G ν (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 Refs. [22][23][24]).

Generalized mean-field approximation
We consider the GMFA for the DSS neglecting the selfenergy Σ(q, ω) in equation (6). Then, for a lattice with basis the zero-order DSS is given by [46]: For the static spin susceptibility we obtain The direct calculation of the matrix elements m αβ (q) yields where From the symmetry of our model it is obvious that we have only two different nn correlation functions: C 1 within the plane and C z between neighboring planes. To calculate the frequency matrix F (q) in equation (7), we start from the second derivative −S ν iα that 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 where the index i + i z denotes the unit cell with position R i+iz = R i + a 3 . Here the vertex renormalization parameters α 1 and α 2 are attached to the nn and the next-nn correlation functions C 1 and C 2 respectively. In the 3D case we introduce α z associated with the nn interlayer correlation function C z . For the next-nn correlation functions between the layers C zz = S y iA S y i+2iz A and C 1z = S y iA S y j+iz B we attach the same vertex parameter α 2 as for the next-nn within the layer. Using these decouplings we obtain the frequency matrix F (q): where with Since the matrices m(q) and F (q) commute, it is convenient to solve equation (7) by introducing the eigenvalues m ± (q) and the normalized eigenvectors |E ± (q) of the matrix (9): 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)|). Introducing the Fourier representation for the correlation function as in equation (4), and using the spectral representation for the GF, we calculate the spin correlation function: The correlation functions C R are written as: where C αα = −C α =β = C, and the wave-vector Q characterizes the LRO. The condensation part C appears in the ordered phase when ω + (q) condensates at Q which determines the LRO, ω + (Q) = 0. In the case of AF order in the two-sublattice model we have Q = (0, 0, π), and the staggered magnetization m st is determined by Let us consider the uniform static susceptibility χ = χ AA (0) + χ AB (0) = χ − (0) and the staggered susceptibility χ st = χ AA (Q) − χ AB (Q) = χ + (Q). Expanding χ − (q) around q = 0 we get Considering q ⇒ 0 from different directions in momentum space we can write the isotropy condition for the uniform static spin susceptibility which should not depend on the direction of q ⇒ 0: lim qx,y →0 χ − (q)| qz =0 = lim qz→0 χ − (q)| qx,y =0 (see Refs. [32][33][34]38,42,46]). This condition gives us the relation between the intralayer and interlayer correlation functions: We expand χ + (q) in the neighborhood of the AF vector Q and obtain χ + (Q+k) = χ + (Q)[1+ξ 2 xy (k 2 x +k 2 y )+ξ 2 z k 2 z ] −1 , where for the intralayer correlation length ξ xy we get: At zero temperature in the 2D case or at T = T N in the 3D case, the LRO occurs when both the correlation length (27) and χ + (Q) diverge. Table 1. Results for the 2D honeycomb AF Heisenberg model at T = 0: ground-state energy per site E0/Ns, staggered magnetization mst and uniform static susceptibility χ(0) at T = 0, obtained by QMC [10][11][12], spin-wave theory [14], series expansion [15,16], slave boson method [17], variational RVB wave function [18], coupled cluster method [4], and by our theory using mst as input. References

Results
To evaluate the spin-excitation spectrum and the thermodynamic properties, the correlation functions C R,αβ and the vertex parameters α 1 , α 2 , and α z 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 equation (23) for calculating the correlation functions, we have the sum rule C R=0,αα = S ν iα S ν iα = 1/4, the isotropy condition (26), and the LRO condition ω + (Q) = 0. That is, we have more parameters than equations. To obtain a closed system of equations, we adjust the staggered magnetization at T = 0 to the QMC values m st (0) = 0.2681 for J z = 0 and m st (0) = 0.41 for J z = J [11], where we take the latter value in the region 0 < J z ≤ 1. For finite temperatures we follow references [27,31,33,34,42] and use the ansatz The results of our self-consistent calculations are presented in Table 1 and in Figures 2-8. As can be seen from Table 1, our results for the ground-state energy E 0 /JN s per site and for the zero-temperature uniform static susceptibility χ(0) rather well agree with the calculations by various methods. Note that the ground-state energy per site taken from reference [11], E 0 /JN s = −0.5445, is related to the ground-state energy e 0 /J per bond, e 0 /J = −0.36303 obtained in reference [11], by E 0 /N s = 3 2 e 0 . In Figure 2 the nn and next-nn correlation functions in the 2D model are plotted. With increasing temperature the AF SRO, reflected in alternating signs of C 1 and C 2 , is weakened which corresponds to a decrease of the correlation functions and vertex parameters (see inset of Fig. 2).
The uniform static spin susceptibility χ is shown in Figure 3 for the 2D and 3D (J z = J) cases. The increase of χ with temperature is due to the reduction of AF short-range correlations which is connected with a weakening of the spin stiffness against the orientation of spins along a homogeneous magnetic field. The further decrease of AF SRO results in a crossover to the high-temperature Curie-Weiss behavior, so that χ reveals   Fig. 3. Uniform static spin susceptibility Jχ versus temperature for Jz = 0 (solid) and Jz = J (dashed), compared with the QMC data [11] for Jz = 0 (diamonds) and Jz = J (circles). a maximum at a temperature of the order of the exchange interaction J. For J z = 0 we obtain the maximum value χ max = 0.122/J at T max = 0.96J. This is close to the QMC result of reference [11]: χ max = 0.117/J at T max = 0.72J. In the 3D case with J z = J our results for the temperature dependence of χ in Figure 3 are in a very good agreement with the QMC calculations in reference [11]. The interlayer coupling J z > 0 lowers the susceptibility and slightly shifts the maximum of χ to higher temperatures in comparison with the 2D case, which is due to the J z -induced enhancement of AF SRO.  Let us compare the maximum temperature T max with the value obtained by susceptibility experiments on the the β-Cu 2 V 2 O 2 compound, T max exp = 50 K [49]. Taking J = 60 K [19] we get T max = 58 K which is near to the experimental value.

7-
In Figure 4 the staggered magnetization m st is depicted. It reveals a second-order transition to the AF phase below the Néel temperature T N . For J z = 0.2J we obtain T N = 0.52J, whereas the QMC simulations of reference [19] yield the lower value T QMC N = 0.41J. The inequality T GMF A N > T QMC N was also found for the AFHM on the stacked square lattice [42]. Considering β-Cu 2 V 2 O 2 with J z = 0.2J and J = (60-66) K [19] we have T N (30)(31)(32)(33)(34)(35) K that is close to the experimental value T exp N = 26 K [49]. For J z = J we find T N = 0.78J and a good agreement of the temperature dependence of m st with the QMC data for the largest system size obtained in reference [11].
The Néel temperature as a function of the interlayer coupling J z is shown in Figure 5, where T N reveals the strongest decrease with decreasing J z for J z /J 1 and tends to zero logarithmically for J z → 0. The J z dependence of T N in the region 0 ≤ J z ≤ 0.2J can be approximated with the accuracy of 1% by the formula where A = 1.48, B = 1.25. This behavior resembles the empirical formula proposed in reference [50] and the result of references [42,47]. Note that qualitatively the same law (28) with lim Jz→0 ∂T N /∂J z = ∞ was also found in the random phase approximation (see Ref. [42]). In Figure 6 the influence of the interlayer coupling on the temperature dependence of the intralayer correlation length ξ xy is illustrated. In the 2D case, ξ xy diverges exponentially at T = 0. For J z > 0, ξ xy diverges at T N , since the gap ω + (Q) closes as T approaches T N from above. In the vicinity of T N , ξ −1 xy behaves as T − T N (corresponding to the critical index ν = 1), as it was also found by previous mean-field approaches (see Refs. [33,34,42]).
The spectrum of spin excitations ω ± (q) is shown in Figures 7 and 8 along the symmetry directions X(−1, 0) → Γ (0, 0) → Y (0, 1) → Γ (1, 1) → M (1/2, 1/2) → Γ of the BZ using the same notation as in reference [8]. In the LRO phase, i.e., for T = 0 at J z = 0 and for T < T N = 0.52J at J z = 0.2J, the spin excitations are spin waves with gapless branches depicted in  Figures 7a and 8a. 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 temperatures slightly above the transition temperature, where the correlation length is large enough, this condition can be fulfilled. The validity region of the spin-wave picture shrinks with increasing temperature, where predominantly high-energy magnons may be observed. Such a situation was encountered in the study of dibromo Ni complexes given in reference [51]. For q < q c the spin-wave picture breaks down. Considering T = 0.6J in Figures 7b  and 8b, for J z = 0(0.2J) we get ξ −1 = 0.86(0.41) (see Fig. 6). Correspondingly, in the whole q range of Figure 7b we have q < q c , and the spin excitations may be named "paramagnon" excitations with the energies ω ± (q). Note that such excitations have been measured, for example, by resonant inelastic X-ray scattering on high-T c superconductors [52,53]. In Figure 8b, due to the lower value of q c as compared with Figure 7b, we have q region with q > q c , where the spin excitations are high-energy magnons. Thus, our spin-excitation spectra reveal a smooth crossover from spin-wave to paramagnon behavior depending on the wave number and temperature. In the upper (optical) branch, at T N a gap is opening at the Γ point, i.e., at the AF wave vector Q characterizing the LRO phase in the twosublattice model (Q = (0, 0) at J z = 0 and Q = (0, 0, π) at J z > 0). As can be seen in Figures 7 and 8, the spinexcitation energies are decreasing with increasing temperature. Let us point out that the spectrum of spin excitations is calculated in the GMFA which neglects finite-time effects. To take them into account, one has to determine the self-energy, as has been done in references [23,24]. As it turns out, in the square-lattice Heisenberg model the damping of spin excitations is quite small. Therefore, we believe that in our model the inclusion of damping will not qualitatively change our results for the spin-excitation spectrum. Provided that a spin-1/2 antiferromagnet with a real honeycomb-lattice structure may be found, it is desirable to confirm our findings on the optical branch ω + (q) by scattering experiments, where the intensity of the scattering is determined by the staggered susceptibility χ + (Q) given by equation (20), so that the gap ω + (Q) may be measured.

Conclusion
In this paper we have evaluated thermodynamic quantities and spin-excitation spectra of the AF Heisenberg model on the stacked honeycomb lattice by calculating the dynamic spin susceptibility within a spin-rotation-invariant generalized mean-field approach for arbitrary temperatures. Our main focus was the analysis of the temperature dependence of the uniform static susceptibility in the paramagnetic phase, which we have explained in terms of AF SRO, and the calculation of the Néel temperature in dependence on the interlayer coupling. Our results are in a good agreement with available QMC and experimental data. From this we conclude that our investigation forms a good basis for forthcoming studies of extended layered honeycomb Heisenberg models (e.g., hole hopping).
The authors would like to thank J. Richter for valuable discussions. The financial support by the Heisenberg-Landau program of JINR is acknowledged. One of the authors (N.P.) thanks the Directorate of the MPIPKS for the hospitality extended to him during his stay at the Institute. Open access funding provided by Max Planck Society.