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

TL;DR

This paper develops a Green-function theoretical approach to study spin excitations and thermodynamics in the antiferromagnetic Heisenberg model on a layered honeycomb lattice, providing insights into magnetic properties and phase transitions.

Contribution

It introduces a spin-rotation-invariant Green-function method with a generalized mean-field approximation for arbitrary temperatures, enabling calculation of thermodynamic quantities and spin spectra.

Findings

Good agreement with numerical cluster computations

Accurate prediction of N'eel temperature for various interlayer couplings

Reproduction of experimental data on -Cu2V2O2

Abstract

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'eel 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\'{e}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 \beta-Cu2V2O2 compound.

Tables

Table 1. Table 1: Results for the 2D honeycomb AF Heisenberg model at T = 0 𝑇 0 T=0 : ground-state energy per site E 0 / N s subscript 𝐸 0 subscript 𝑁 𝑠 E_{0}/N_{s} , staggered magnetization m s ​ t subscript 𝑚 𝑠 𝑡 m_{st} and uniform static susceptibility χ ​ ( 0 ) 𝜒 0 \chi(0) at T = 0 𝑇 0 T=0 , obtained by QMC  Reger89 ; Low09 ; Jiang12 , spin-wave theory  Weihong91 , series expansion Oitmaa92 ; Oitmaa12 , slave boson method  Mattsson94 , variational RVB wave function  Noorbakhsh09 , coupled cluster method  Bishop12 , and by our theory using m s ​ t subscript 𝑚 𝑠 𝑡 m_{st} as input.

Refs	$E_0/J{N}_s$	$m_{st}$	$\chi \left(0\right)J$
Reger89	- 0.5445	0.22(3)	-
Low09	- 0.5445	0.2681	0.05188
Jiang12	-	0.2688	-
Weihong91	-0.5485	0.2418	0.0457
Oitmaa92	-0.5443	0.266(9)	0.0756
Oitmaa12	-0.54	0.265	-
Mattsson94	-	-	0.0683
Noorbakhsh09	-0.5440	0.26	-
Bishop12	-0.53	0.272	-
Our results	- 0.59	0.2681	0.056

Equations

$${\delta }_1=\frac{a_0}{2}(\sqrt{3},-1),{\delta }_2=-\frac{a_0}{2}(\sqrt{3},1),{\delta }_3=a_0(0,1).$$

$$H=\frac{J}{2}\sum _{⟨i\alpha ,j\beta {⟩}_{xy}}{\mathbf{S}}_{i\alpha }{\mathbf{S}}_{j\beta }+\frac{J_z}{2}\sum _{⟨i\alpha ,j\alpha {⟩}_z}{\mathbf{S}}_{i\alpha }{\mathbf{S}}_{j\alpha },$$

$$G_{ij,\alpha \beta }^{\nu }(t-t^{\prime })$$

$$\displaystyle{G}^{\nu}({\bf q},\omega)=\left(\begin{array}[]{cccc}\langle\!\langle S^{\nu}_{{\bf q}A}|S^{\nu}_{{\bf q}A}\rangle\!\rangle_{\omega}\;\langle\!\langle S^{\nu}_{{\bf q}A}|S^{\nu}_{{\bf q}B}\rangle\!\rangle_{\omega}\\ \langle\!\langle S^{\nu}_{{\bf q}B}|S^{\nu}_{{\bf q}A}\rangle\!\rangle_{\omega}\;\langle\!\langle S^{\nu}_{{\bf q}B}|S^{\nu}_{{\bf q}B}\rangle\!\rangle_{\omega}\\ \end{array}\right).$$

$$\chi (\mathbf{q},\omega )=[F(\mathbf{q})+\Sigma (\mathbf{q},\omega )-{\omega }^2{\tau }_0{]}^{-1}\times m(\mathbf{q}).$$

$$[F(\mathbf{q})-{\omega }^2{\tau }_0]\times \chi (\mathbf{q},\omega )=m(\mathbf{q}).$$

$${\chi }_{(}\mathbf{q},0)\equiv {\chi }_{(}\mathbf{q})=F^{-1}(\mathbf{q})\times m(\mathbf{q}).$$

$$m({\bf q})=\left(\begin{array}[]{cc}m_{AA}({\bf q})&m_{AB}({\bf q})\\ m_{AB}^{*}({\bf q})&m_{AA}({\bf q})\end{array}\right),$$

$$m_{AA}(\mathbf{q})=-6J{C}_1-4J_z{C}_z(1-\cos q_z),$$

$$m_{AB}(\mathbf{q})=2J{C}_1{\gamma }_1(\mathbf{q}),{\gamma }_1(\mathbf{q})=\sum _i\exp (i\mathbf{q}{\delta }_i),$$

$$|{\gamma }_1(q){|}^2=1+4\cos (\frac{\sqrt{3}}{2}{q}_x)[\cos (\frac{\sqrt{3}}{2}{q}_x)+\cos (\frac{3}{2}{q}_y)].$$

$$S_{iA}^x{S}_{jB}^y{S}_{kA}^y$$

$$S_{jB}^x{S}_{iA}^y{S}_{kA}^y$$

$$S_{jB}^x{S}_{iA}^y{S}_{i+i_{z}A}^y$$

$$F({\bf q})=\frac{1}{2}\left(\begin{array}[]{cc}F_{AA}({\bf q})&F_{AB}({\bf q})\\ F^{*}_{AB}({\bf q})&F_{AA}({\bf q})\end{array}\right),$$

$$F_{AA}(\mathbf{q})$$

$$F_{AB}(\mathbf{q})$$

$$f_{AB}(\mathbf{q})$$

$${\gamma }_2(\mathbf{q})$$

$$[m(\mathbf{q})-{\tau }_0{m}_{\pm }(\mathbf{q})]|E_{\pm }(\mathbf{q})⟩=0,$$

$$m_{\pm }(\mathbf{q})$$

$$|E_{\pm }(\mathbf{q})⟩$$

$${\omega }_{\pm }^2(\mathbf{q})=\frac{1}{2}(F_{AA}(\mathbf{q})\pm |{\gamma }_1(\mathbf{q})|f_{AB}(\mathbf{q})).$$

$${\chi }_{\alpha \beta }(\mathbf{q},\omega )=\sum _{j=\pm }{\chi }_j(\mathbf{q},\omega )⟨\alpha |E_j(\mathbf{q})⟩⟨E_j(\mathbf{q})|\beta ⟩,$$

$${\chi }_{\pm }(\mathbf{q},\omega )=\frac{m_{\pm }(\mathbf{q})}{{\omega }_{\pm }^2(\mathbf{q})-{\omega }^2},$$

$$C_{ij,\alpha \beta }=⟨S_{i\alpha }^{\nu }{S}_{j\beta }^{\nu }⟩=\frac{1}{N}\sum _{\mathbf{q}}{\mathrm{e}}^{i\mathbf{q}({\mathbf{R}}_{\mathbf{i}}-{\mathbf{R}}_{\mathbf{j}})}{C}_{\alpha \beta }(\mathbf{q}),$$

$$C_{\alpha \beta }(\mathbf{q})=⟨S_{\mathbf{q}\alpha }^{\nu }{S}_{\mathbf{q}\beta }^{\nu }⟩$$

$$=\sum _{j=\pm }\frac{m_j(\mathbf{q})}{2{\omega }_j(\mathbf{q})}\coth \frac{{\omega }_j(\mathbf{q})}{2T}⟨\alpha |E_j(\mathbf{q})⟩⟨E_j(\mathbf{q})|\beta ⟩.$$

$$C_{\mathbf{R},\alpha \beta }=⟨S_{0\alpha }^{\nu }{S}_{\mathbf{R},\beta }^{\nu }⟩$$

$$=\frac{1}{N}\sum _{\mathbf{q}\ne \mathbf{Q}}{C}_{\alpha \beta }(\mathbf{q})e^{i\mathbf{qR}}+C_{\alpha \beta }{e}^{i\mathbf{QR}},$$

$$(m_{st}{)}^2=3C.$$

$$\chi =\frac{-4C_1}{J+2J{\alpha }_1{C}_1+8J{\alpha }_2{C}_2+8J_z({\alpha }_2{C}_{1z}-{\alpha }_1{C}_1)}.$$

$$\frac{-4C_1}{1+2{\alpha }_1{C}_1+8{\alpha }_2{C}_2+8(J_z/J)({\alpha }_2{C}_{1z}-{\alpha }_1{C}_1)}$$

$$=\frac{-2C_z}{12({\alpha }_2{C}_{1z}-{\alpha }_z{C}_z)+(J_z/J)(1+4{\alpha }_2{C}_{zz})}.$$

Full text

I 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 Manousakis91 ). According to the Mermin-Wagner theorem Mermin66 , 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 Chakravarty89 . 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  Bishop12 ; Li12 , in the $J_{1}$ – $J_{2}$ – $J_{3}$ Heisenberg model on the 2D honeycomb lattice four competing magnetic phases (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 bond-depending Kitaev interaction $K_{\alpha}$ Kitaev06 . Besides the spin-liquid 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., Chaloupka10 ; Chaloupka13 ). In a model with anisotropic $K_{\alpha}$ interactions, e.g., $K_{z}>K_{x}=K_{y}$, the LRO survives at finite temperature as shown in Ref. Vladimirov16 .

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 $\chi(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 Reger89 ; Low09 ; Jiang12 ; Castro06 , spin-wave theory Castro06 ; Weihong91 , series expansions around the Ising limit Oitmaa92 ; Oitmaa12 , Schwinger boson method Mattsson94 , and variational RVB wave function approach Noorbakhsh09 . In Ref. Tsirlin10 it was suggested to consider the $\beta$-Cu2V2O2 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}\simeq 26$ K.

Besides the honeycomb AFHM, the honeycomb Hubbard model has been studied too. In Ref. Peres04 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 Ref. Aria15 , using the two-particle self-consistent approach for the Hubbard model on the honeycomb lattice, the semimetal to spin-liquid 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 Vladimirov05 ; Vladimirov09 using the generalized mean-field approximation (GMFA), as has been done in our study of the compass-Heisenberg model on the square lattice Vladimirov15 . 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 self-consistently, 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. KY72 ; ST91 ; BB94 ; WI97 ; Ihle99 ; Siurakshina00 ; Siurakshina01 ; ISW99 ; YF00 ; BCL02 ; Schmalfuss06 ; HRI08 ; JIB08 ; Junger09 ; MKB09 ; BMS11 ; HRG13 ; Schmalfuss04 ; Schmalfuss05 ). In particular, let us mention an application of the GMFA to study the spin-$1/2$ AFHM on the stacked kagomé lattice Schmalfuss04 , where no LRO was found even for a strong interlayer coupling due to the frustration character of the AF interaction on this lattice.

In Sect. II 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 Sect. III. The numerical results and discussion are given in Sect. IV. The conclusion can be found in Sect. V.

II Model and dynamic spin susceptibility

We consider congruently stacked honeycomb layers shown in Fig. 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 ${{\bf\delta}_{j}}$, and sites on $B$ are connected to $A$ by vectors $-{{\bf\delta}_{j}}$:

[TABLE]

The basis vectors in the layer are ${\bf a}_{1}={{\bf\delta}_{3}}-{{\bf\delta}_{2}}=({a_{0}}/{2})(\sqrt{3},3)$ and ${\bf a}_{2}={{\bf\delta}_{3}}-{{\bf\delta}_{1}}=({a_{0}}/{2})(-\sqrt{3},3)$, the lattice constant is $a=|{\bf a}_{1}|=|{\bf a}_{2}|=\sqrt{3}a_{0}$, 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 ${\bf k}_{1}=({2\pi}/{3})(\sqrt{3},1)$ and ${\bf k}_{2}=({2\pi}/{3})(-\sqrt{3},1)$.

The Heisenberg model on this layered honeycomb lattice can be written as

[TABLE]

where $i$, $j$ count the unit cells, $\alpha$ and $\beta$ are the sublattice indexes, $\langle i\alpha,j\beta\rangle_{xy}$ and $\langle i\alpha,j\beta\rangle_{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) Zubarev60 :

[TABLE]

where $[A,B]=AB-BA$, $\,A(t)=\exp(iHt)A\exp(-iHt)$, and $\,\theta(x)$ is the Heaviside function. In the SRI theory all GF components $\nu=x,y,z$ are equivalent; therefore, we consider only one of them, e.g., $\nu=x$. The Fourier representation in $({\bf q},\omega)$-space is defined by the relation:

[TABLE]

where $N$ and ${\bf R}_{i}$ is the number and position of unit cells, respectively, and the GF matrix reads

[TABLE]

In the relaxation-function theory developed on the basis of the equation of motion method in Refs. Vladimirov05 ; Vladimirov09 we obtain the following representation of the DSS $\chi({\bf q},\omega)=-{G}^{\nu}({\bf q},\omega)$:

[TABLE]

Here, $F({\bf q})$ is the frequency matrix of spin excitations in the GMFA, where the approximation $-\ddot{S}^{\nu}_{{\bf q}\alpha}=[[S^{\nu}_{\bf q\alpha},H],H]=\Sigma_{\beta}F_{\alpha\beta}({\bf q})S^{\nu}_{{\bf q}\beta}$ is made, $\tau_{0}$ is the unity matrix, and $m({\bf q})$ is the moment matrix with components $m_{\alpha\beta}({\bf q})=\langle[i\dot{S}^{\nu}_{{\bf q}\alpha},S^{\nu}_{{\bf q}\beta}]\rangle=\langle[[S^{\nu}_{{\bf q}\alpha},H],S^{\nu}_{{\bf q}\beta}]\rangle$. The self-energy $\Sigma({\bf q},\omega)$ can be expressed exactly by a multispin GF (see Refs. Vladimirov05 ; Vladimirov09 ).

III Generalized mean-field approximation

We consider the GMFA for the DSS neglecting the self-energy $\Sigma({\bf q},\omega)$ in Eq. (8). Then, for a lattice with basis the zero-order DSS is given by  Schmalfuss04 :

[TABLE]

For the static spin susceptibility we obtain

[TABLE]

The direct calculation of the matrix elements $m_{\alpha\beta}({\bf q})$ yields

[TABLE]

where

[TABLE]

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({\bf q})$ in Eq. (9), we start from the second derivative $-\ddot{S}^{\nu}_{i\alpha}$ that is proportional to products of three spin operators on different lattice sites along nn sequences, e.g., $\langle iA,jB,kA\rangle$. We perform the decoupling of them as follows

[TABLE]

where the index $i+i_{z}$ denotes the unit cell with position ${\bf R}_{i+i_{z}}={\bf R}_{i}+{\bf a}_{3}$. Here the vertex renormalization parameters $\alpha_{1}$ and $\alpha_{2}$ are attached to the nn and the next-nn correlation functions $C_{1}$ and $C_{2}$ respectively. In the 3D case we introduce $\alpha_{z}$ associated with the nn interlayer correlation function $C_{z}$. For the next-nn correlation functions between the layers $C_{zz}=\langle S^{y}_{iA}S^{y}_{i+2i_{z}A}\rangle$ and $C_{1z}=\langle S^{y}_{iA}S^{y}_{j+i_{z}B}\rangle$ we attach the same vertex parameter $\alpha_{2}$ as for the next-nn within the layer. Using these decouplings we obtain the frequency matrix $F({\bf q})$:

[TABLE]

where

[TABLE]

with

[TABLE]

Since the matrices $m({\bf q})$ and $F({\bf q})$ commute, it is convenient to solve Eq. (9) by introducing the eigenvalues $m_{\pm}({\bf q})$ and the normalized eigenvectors $|E_{\pm}({\bf q})\rangle$ of the matrix (11):

[TABLE]

which are given by

[TABLE]

For the same eigenvectors $|E_{\pm}({\bf q})\rangle$ the spin-excitation frequencies are obtained as

[TABLE]

In this notation the DSS reads:

[TABLE]

where

[TABLE]

and $\langle\alpha|E_{\pm}\rangle\langle E_{\pm}|\beta\rangle=1/2$ for $\alpha=\beta$, otherwise $\langle\alpha|E_{\pm}\rangle\langle E_{\pm}|\beta\rangle=\mp\gamma_{1}({\bf q})/(2|\gamma_{1}({\bf q})|)$.

Introducing the Fourier representation for the correlation function as in Eq. (4),

[TABLE]

and using the spectral representation for the GF, we calculate the spin correlation function:

[TABLE]

The correlation functions $C_{\bf R}$ are written as:

[TABLE]

where $C_{\alpha\alpha}=-C_{\alpha\neq\beta}=C$, and the wave-vector ${\bf Q}$ characterizes the LRO. The condensation part $C$ appears in the ordered phase when $\omega_{+}({\bf q})$ condensates at ${\bf Q}$ which determines the LRO, $\omega_{+}({\bf Q})=0$. In the case of AF order in the two-sublattice model we have ${\bf Q}=(0,0,\pi)$, and the staggered magnetization $m_{st}$ is determined by

[TABLE]

Let us consider the uniform static susceptibility $\chi=\chi_{AA}(0)+\chi_{AB}(0)=\chi_{-}(0)$ and the staggered susceptibility $\chi_{st}=\chi_{AA}({\bf Q})-\chi_{AB}({\bf Q})=\chi_{+}({\bf Q})$. Expanding $\chi_{-}({\bf q})$ around ${\bf q}=0$ we get

[TABLE]

Considering ${\bf q}\Rightarrow 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\Rightarrow 0$: $\lim_{q_{x,y}\rightarrow 0}\chi_{-}({\bf q})|_{q_{z}=0}=\lim_{q_{z}\rightarrow 0}\chi_{-}({\bf q})|_{q_{x,y}=0}$ (see Refs. Ihle99 ; Siurakshina00 ; Schmalfuss04 ; Schmalfuss06 ; Junger09 ). This condition gives us the relation between the intralayer and interlayer correlation functions:

[TABLE]

We expand $\chi_{+}({\bf q})$ in the neighborhood of the AF vector ${\bf Q}$ and obtain $\chi_{+}({\bf Q}+{\bf k})=\chi_{+}({\bf Q})[1+\xi_{xy}^{2}(k_{x}^{2}+k_{y}^{2})+\xi_{z}^{2}k_{z}^{2}]^{-1}$, where for the intralayer correlation length $\xi_{xy}$ we get:

[TABLE]

At zero temperature in the 2D case or at $T=T_{N}$ in the 3D case, the LRO occurs when both the correlation length (31) and $\chi_{+}({\bf Q})$ diverge.

IV Results

To evaluate the spin-excitation spectrum and the thermodynamic properties, the correlation functions $C_{{\bf R},\alpha\beta}$ and the vertex parameters $\alpha_{1}$, $\alpha_{2}$, and $\alpha_{z}$ appearing in the spectrum $\omega_{\pm}({\bf 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 Eq. (27) for calculating the correlation functions, we have the sum rule $C_{{\bf R}=0,\alpha\alpha}=\langle{S}^{\nu}_{i\alpha}{S}^{\nu}_{i\alpha}\rangle=1/4\;$, the isotropy condition (30), and the LRO condition $\omega_{+}({\bf 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$ Low09 , where we take the latter value in the region $0<J_{z}\leq 1$. For finite temperatures we follow Refs. ST91 ; WI97 ; Siurakshina00 ; Siurakshina01 ; Junger09 and use the ansatz $r_{\alpha}(T)\equiv[\alpha_{2}(T)-1]/[\alpha_{1}(T)-1]=r_{\alpha}(0)$.

The results of our self-consistent calculations are presented in Table 1 and in Figs. 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 $\chi(0)$ rather well agree with the calculations by various methods. Note that the ground-state energy per site taken from Ref. Low09 , $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 Ref. Low09 , by $E_{0}/N_{s}=\frac{3}{2}e_{0}$. In Fig. 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 $\chi$ is shown in Fig. 3 for the 2D and 3D ($J_{z}=J$) cases. The increase of $\chi$ 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 $\chi$ reveals a maximum at a temperature of the order of the exchange interaction $J$. For $J_{z}=0$ we obtain the maximum value $\chi^{max}=0.122/J$ at $T^{max}=0.96J$. This is close to the QMC result of Ref. Low09 : $\chi^{max}=0.117/J$ at $T^{max}=0.72J$. In the 3D case with $J_{z}=J$ our results for the temperature dependence of $\chi$ in Fig. 3 are in a very good agreement with the QMC calculations in Ref. Low09 . The interlayer coupling $J_{z}>0$ lowers the susceptibility and slightly shifts the maximum of $\chi$ 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 $\beta$-Cu2V2O2 compound, $T^{max}_{exp}=50$K He2008 . Taking $J=60$K Tsirlin10 we get $T^{max}=58$K which is near to the experimental value.

In Fig. 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 Ref. Tsirlin10 yield the lower value $T_{N}^{QMC}=0.41J$. The inequality $T_{N}^{GMFA}>T_{N}^{QMC}$ was also found for the AFHM on the stacked square lattice Junger09 . Considering $\beta$-Cu2V2O2 with $J_{z}=0.2J$ and $J=$(60–66)K Tsirlin10 we have $T_{N}\simeq$(30–35)K that is close to the experimental value $T_{N}^{exp}=26$K He2008 . 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 Ref. Low09 .

The Néel temperature as a function of the interlayer coupling $J_{z}$ is shown in Fig. 5, where $T_{N}$ reveals the strongest decrease with decreasing $J_{z}$ for $J_{z}/J\ll 1$ and tends to zero logarithmically for $J_{z}\rightarrow 0$. The $J_{z}$ dependence of $T_{N}$ in the region $0\leq J_{z}\leq 0.2J$ can be approximated with the accuracy of $1\%$ by the formula

[TABLE]

where $\,A=1.48,\,B=1.25\,$. This behavior resembles the empirical formula proposed in Ref. Yasuda05 and the result of Refs. Junger09 ; Schmalfuss05 . Note that qualitatively the same law (32) with $\lim_{J_{z}\rightarrow 0}\partial T_{N}/\partial J_{z}=\infty$ was also found in the random phase approximation (see Ref. Junger09 ).

In Fig. 6 the influence of the interlayer coupling on the temperature dependence of the intralayer correlation length $\xi_{xy}$ is illustrated. In the 2D case, $\xi_{xy}$ diverges exponentially at $T=0$. For $J_{z}>0$, $\xi_{xy}$ diverges at $T_{N}$, since the gap $\omega_{+}({\bf Q})$ closes as $T$ approaches $T_{N}$ from above. In the vicinity of $T_{N}$, $\xi_{xy}^{-1}$ behaves as $T-T_{N}$ (corresponding to the critical index $\nu=1$), as it was also found by previous mean-field approaches (see Refs. Siurakshina00 ; Junger09 ).

The spectrum of spin excitations $\omega_{\pm}({\bf q})$ is shown in Figs. 7 and 8 along the symmetry directions $X(-1,0)\rightarrow\Gamma(0,0)\rightarrow Y(0,1)\rightarrow\Gamma^{\prime}(1,1)\rightarrow M(1/2,1/2)\rightarrow\Gamma$ of the BZ using the same notation as in Ref. Chaloupka13 . 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 Fig. 7(a) and  8(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\pi\xi^{-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 Ref. Junger05 . For $q<q_{c}$ the spin-wave picture breaks down. Considering $T=0.6J$ in Figs. 7(b) and 8(b), for $J_{z}=0(0.2J)$ we get $\xi^{-1}=0.86(0.41)$ (see Fig. 6). Correspondingly, in the whole $q$ range of Fig. 7(b) we have $q<q_{c}$, and the spin excitations may be named “paramagnon” excitations with the energies $\omega_{\pm}({\bf q})$. Note that such excitations have been measured, for example, by resonant inelastic x-ray scattering on high-Tc superconductors Tacon11 ; Tacon13 . In Fig. 8(b), due to the lower value of $q_{c}$ as compared with Fig. 7(b), 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 $\Gamma$ point, i.e., at the AF wave vector ${\bf Q}$ characterizing the LRO phase in the two-sublattice model (${\bf Q}=(0,0)$ at $J_{z}=0$ and ${\bf Q}=(0,0,\pi)$ at $J_{z}>0$). As can be seen in Figs.  7 and 8, the spin-excitation 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 Refs. Vladimirov09 . 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 $\omega_{+}({\bf q})$ by scattering experiments, where the intensity of the scattering is determined by the staggered susceptibility $\chi_{+}({\bf Q})$ given by Eq.(24), so that the gap $\omega_{+}({\bf Q})$ may be measured.

V 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).

Acknowledgements.

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.