Properties of the two-dimensional spin-1 2 Heisenberg model on a honeycomb lattice with interlayer coupling

The magnetic properties of the two-dimensional S = 1 2 (quantum) antiferromagnetic Heisenberg model on a honeycomb lattice with and without interlayer coupling are studied by means of a continuous Euclidean time Quantum-Monte-Carlo algorithm. The internal energy, the magnetic susceptibility and the staggered magnetization are determined in the full temperature range. For the two-dimensional system the groundstate energy/bond is found to be E 0 = −0.36303(13), and the zero temperature staggered magnetization mst = 0.2681(8). For coupled planes of honeycomb systems a phase transition from an ordered phase to a disordered phase is found at T/J = 0.695(10).


Introduction
Besides spin models on the two-dimensional square lattice and on the cubic lattice other lattice topologies are interesting both due to their principal importance as quantum systems and because of their phenomenological relevance.
In particular, the honeycomb lattice which is a two-dimensional bipartite lattice with minimal coordination number has attracted continuous attention over the years.Early studies of its finite temperature properties by means of high temperature expansions date back to 1973 [1,2].A Quantum-Monte-Carlo study was performed by Reger, Riera and Young [3] in 1989 using the discrete world line algorithm.Bishop and Rosenfeld [4] used a coupled cluster method to evaluate the ground-state properties of the XXZ-model on the honeycomb lattice and compared with the series expansions by Weihong, Oitmaa, and Hamer [5,6].An elaborate study of the two-dimensional honeycomb lattice in the presence of site-percolation can be found in [7].
Also, only recently InCu 2/3 V 1/3 O 3 has been discussed as a possible candidate for a spin-1 2 substance [8,9] on a honeycomb lattice.The advent of this and possible other upcoming new substances as well as the availability of high precision Monte-Carlo methods [10,11] prompted us to make a renewed study of spin systems on the honeycomb lattice.We study both the twodimensional honeycomb system focusing on issues which have not been treated in the literature and on honeycomb planes coupled in the third dimension.To our knowledge this type of threedimensional system has not been studied thus far by quantum Monte-Carlo.For our investigation the continuous time loop algorithm [11] was employed, which allows us to study spin systems with unprecedented precision even at very low temperatures.The short autocorrelation times of the method are achieved by nonlocal loop updates and no extrapolation in Trotter number is needed, since the algorithm works with a continuous Trotter time variable.
Beard and Wiese [11] originally suggested this algorithm for the two-dimensional square lattices, but it can also be used for three-dimensional systems and to some extent also to other lattice topologies.Simulation of a honeycomb structure is achieved by "viewing" the honeycomb lattice as a double dimerized model with Hamiltonian 2 and J y = J.I. e. the Hamiltonian equation ( 1) has different dimerizations in x-direction for y even and y odd, where J is the exchange coupling of the honeycomb system.The Hamiltonian equation (1) for one layer is depicted in figure 1.The coupling in z-direction between the layers is J z = J.Note that the cubic lattice in terms of equation ( 1) is readily given by δ = 0 and J x = J y = J z = J.It is well known that including such a dimerization-structure does not pose a problem to the Monte Carlo-algorithm.Using dimerizations δ = ±1 to calculate ladder properties is a well established technique.Also the data shown here for the one-and two-dimensional systems indicate that breaking higher dimensional structures down to lower dimensions by putting bonds to zero works perfectly with this algorithm.
The paper is organized as follows.In section 1 the ground-state energy and the internal energy of the honeycomb lattices are studied as well as the presence of a phase transition in the threedimensional system is demonstrated.Section 2 is devoted to the susceptibility and section 3 to the staggered magnetization.

The internal energy
The ground-state energy of the one-dimensional Heisenberg model is exactly known since the seminal works of H. Bethe [12,13] and L.Hulthén [14].However, there is no counterpart of exact results for higher-dimensional versions of the model.Therefore, over the past years a substantial effort has been made both with approximate analytical and numerical techniques to obtain a value of the ground-state energy of the two-and three-dimensional Heisenberg models.
At low temperatures for the one-dimensional model, Bethe ansatz [18] and conformal theory [19] predict a quadratic dependence of the internal energy U (T ) on temperature.Whereas for two-(three)-dimensional models, U (T ) depends in the leading order spin wave theory on the third (fourth) power of temperature [20,21].For models with linear spin wave dispersion, this behavior does not depend on the details of the coupling, but only on the dimension of the system.
To compare with the above predictions and to determine the coefficients of the leading power laws, low temperature data for the internal energy were calculated both for a honeycomb lattice, and a honeycomb system with three-dimensional coupling J z = J.For reasons of reference the data for a chain, a square and a cubic lattice are also presented.
The coefficients for the one-dimensional system follow from the analysis of [18] and [19].For the square lattice they were determined by Takahashi [22] using modified spin-wave theory and by Avoras and Auerbach using Schwinger boson mean-field theory [23].Our fits to the one-dimensional case are in very good agreement with these predictions.Also for a 32 × 32 square lattice a fit to the data in the temperature range 0.005 T /J 1.1 reproduces the results of Takahashi [22] with good accuracy.For the honeycomb lattice we find (lower inset in 2) for a 32 × 64 lattice and a fit range (T /J 0.8) U hex /J = −0.363+ 0.799(T /J) 3  (2) which corroborates the universal results of Oguchi and Kubo, [21,20] of a T 3 behavior, however the validity is confined to smaller temperatures as compared to the square lattice.
For the honeycomb lattice with interlayer-coupling we find a best fit with a somewhat higher exponent of the temperature, U 3dim hc /J = −0.3124+ 0.325 (T /J) 4.8  (3) hinting at higher order contributions in spin wave theory.Forcing a T 4 fit gives For a graphical depiction of the results see inset of figure 2.

Figure 2.
Monte-Carlo results for the internal energy U/J per bond as a function of temperature, for a 400 site chain (circles), a 32 × 32 square lattice (squares), a honeycomb lattice of size 32 × 64 (diamonds), a honeycomb lattice with three-dimensional coupling (triangles), and a cubic system of size 16 3 (crosses).Upper inset: Fit (equation ( 3)) to low temperature behavior of the three-dimensional system(red), and fit (equation ( 4)) using a T 4 law (green).Lower inset: Low temperature behavior of the internal energy for the honeycomb lattice equation (2).
Next, the internal energy and the specific heat of the honeycomb lattice with three-dimensional coupling in an intermediate temperature range are discussed.Analogous to the Heisenberg model on a three-dimensional cubic lattice, one expects that the three-dimensionally coupled honeycomb lattices undergo a second-order phase transition at a finite temperature T c .Since rigorously speaking phase transitions can be seen only in infinite size systems, in numeric calculations involving systems of limited size one is in general restricted to the study of size dependent quantities, e.g. the increase of thermodynamic quantities with growing system size or scaling laws, to show the presence of the singularity connected to the phase transition.
Here we estimate the critical temperature of the coupled honeycomb planes by analyzing the specific heat.As can be seen in figure 3 the height of the peak of the specific heat grows with the system size in a temperature range, where one does not expect finite size effects.(Note that for not too low temperatures and T < 0.6J and also for T > J the curves of the internal energy U (T ) for different system sizes fall together.)Thus, the unconventional finite size behavior can be interpreted as a manifestation of a phase transition.Also, the location of the maximum of the specific heat is shifted to lower temperatures with growing system size.From this behavior we estimate (see inset in figure 3) that the transition is located at T c /J = 0.695 (10).The temperature T c corresponds approximately to the steep decay of the staggered magnetization analyzed in the section 3, indicating that the transition is between an ordered phase at low temperatures and a magnetically disordered high-temperature phase with vanishing staggered magnetization.As expected, the high temperature behavior of the three curves approaches the same value E ∞ = 0. Since we considered the energy per bond, the fact that the internal energy for the honeycomb lattice lies between the results for the square lattice and the chain is an effect of the quantum fluctuations independent of coordination number.
To obtain the value of the ground-state energy of the honeycomb lattice for the infinite size system the energies U (T = 0.005J) (see figure 4) are extrapolated for system sizes N = 16 . . .36.Taking into account the previous result of equation (2) one can expect that the error due to the finite temperature should be ≈ (0.005) 3 which is smaller than the statistical error.
As in the case of the square lattice the ground-state energy scales in leading order with 1 N 3 .For a N × N lattice with N = 16 . . .36 a linear fit gives E hc 0 = −0.363035(13) for the ground-state energy.This agrees with the value E hc 0 = −0.3630(3) of [3].From the nonlinear σ model description of the two-dimensional Heisenberg model [24][25][26] it is known that for the square lattice, the next order of the finite size behavior is given by a 1/N 4 correction.However, a significant deviation from the 1/N 3 law cannot be extracted from the data.
As a technical point we remark that for this analysis a rectangular lattice with N × 2N sites has been chosen, which amounts to an equal number of unit cells of the honeycomb lattice in x and y-direction.We found, however, that extrapolated quantities do not depend on the shape of the system.This is demonstrated in some detail in case of the staggered magnetization, a quantity which is particularly sensitive to the size of the system, because it is closely related to the correlations.

Susceptibility
In this section the uniform susceptibility is discussed.Figure 5 displays an overview of the results for the chain, the square, the honeycomb lattice and the three-dimensional systems.The susceptibility of the honeycomb lattice reaches its maximum at T max ≈ 0.72J.Both low and high temperature series expansions fail to reach this value and earlier Monte-Carlo simulations did not concentrate on this point.Note that data for the two-dimensional honeycomb-model were already used in [8] to confirm that the magnetically active ions of InCu 2/3 V 1/3 O 3 form a honeycombstructure.
Figure 5. Susceptibility of the Heisenberg model on a honeycomb lattice and a honeycomb lattice with interlayer coupling Jz = J compared to the susceptibility of a square lattice, a onedimensional chain with and a cubic system.The dot dashed line shows the high-temperatureseries expansion of [1,2] for the honeycomb lattice.
For a comparison with experimental data the most prominent feature of a susceptibility curve is its maximum, which most often serves as a first estimate of the exchange coupling.It is, therefore, certainly of phenomenological interest, that T max for the honeycomb lattice lies close to the maximum of the susceptibility of the chain at T max ≈ 0.64J [27], but it is well distinguishable from the square lattice whose maximum is at T max ≈ 0.94J.As for the honeycomb lattice with three-dimensional coupling the maximum is located at T max ≈ 1.03J.For the cubic system T max ≈ 1.22(2)J.Thus, inter-layer couplings shift the maximum of the susceptibility to higher temperatures.Also the slow ascent of the curves is typical of three-dimensional systems.
Next we study the low temperature behavior of the susceptibility of the two-dimensional honeycomb lattice, which has been discussed extensively by spin wave theory and series expansions [5,6].In contrast to the susceptibility of the one-dimensional Heisenberg model, which approaches its zero temperature value Jχ(T = 0) = 0.101322 [28] with infinite slope for two-dimensional lattices spin-wave theory predicts a linear behavior of χ(T ) at low temperatures.We clearly see this linear behavior in our system (see figure 6).Note that the exponential drop of χ(T ) is a finite size effect, originating from the energy gap of the finite systems.But in the temperature range T /J > 0.04, where the results fall both for the 32 × 32 and the 24 × 24 size system, the curves definitely represent the thermodynamic limit.A linear fit to the data for 0.04 T /J 0.15 gives Jχ(T = 0) = 0.05188 (8), which agrees well with the result of series expansion [6], whereas spinwave theory in the first order gives Jχ(T = 0) = 0.111 and thus lies out of the scale of figure 6. Second order spin-wave theory is ≈ 40% too small.Figure 6.The susceptibility of the honeycomb lattice at low temperatures for three system sizes.The solid line is a linear fit to the Monte-Carlo data of the 32 × 32 lattice in the range between T = 0.04J and T = 0.15J.The star marks the series expansion for χ(T = 0) from [6], the triangle the result of second order spin-wave [5].
The failure of spin-wave theory for the honeycomb lattice can be intuitively understood by the fact that quantum fluctuations are larger for the honeycomb than for the square lattice.

The staggered magnetization
It is rigorously known [29] that for the two-dimensional model the thermal expectation value of the staggered magnetization operator x,y vanishes at finite temperatures.But also the zero-temperature magnetization 0| M st |0 , where |0 is the ground-state of the finite system, does vanish, since M st is an operator with momentum q = (π, π) whose matrix elements for all non-degenerate states with fixed momentum are zero.It is thus not obvious how to conclude on the existence of a non-vanishing staggered magnetization from finite size systems, working at low, but nonzero temperatures.
Instead, it is common use to calculate at the lowest possible temperatures the expectation value of the square of M st which is proportional to the static structure factor at q = (π, π).Calculations of this type have a long history [15] for the two-dimensional Heisenberg model, yielding lim N →∞ m st = 0.3070(3) [16,30] for the square lattice.The finite size effects of the staggered magnetization are very large as compared to, e.g.finite size effects of the ground-state energy.However (m st ) 2 shows an almost perfect linear behavior in 1/N (see figure 7) with some small correction proportional to 1/N 2 , which allows to extrapolate with good accuracy.
We thus obtain for the two-dimensional system lim N →∞ m st = 0.2681(8) (7) as result for the infinite size system.Figure 7 also shows data points for rectangular systems, which show a slightly different slope and extrapolate to m st = 0.2671 (10).
Another way to proceed is to calculate the correlation at maximal distance on the periodic lattice [3] and extrapolate it to N → ∞.From this fit one concludes that lim which is within the error equal to m st .The results are compatible with the previous m st = 0.22(3) of [3] and in excellent agreement with [7] who find m st = 0.2677 (6) by means of the stochastic series expansion Monte-Carlo.Again the error due to finite temperature effects is smaller than the statistical error (see insets of figure 7).Thus, the honeycomb lattice displays a reduced magnetic order as compared to the square lattice.However, the Néel-type structure is not destroyed by quantum fluctuations.
The staggered magnetization m st for a honeycomb system with three-dimensional coupling is shown in figure 8.In this case one observes a sudden decrease of m st with temperature, which is a sign of the phase transition observed in section 1.For an increasing system size, the slope of m st increases.To corroborate that the peak found in the specific heat corresponds to a transition between an ordered phase at T < T c and a disordered phase with vanishing staggered magnetization at T > T c we extrapolated the data for m st for system size N × 2N × N , as a function of 1 N .For T > T c we find a very convincing N − 3 2 law extrapolating to m st = 0 for N → ∞.The lower inset in figure 8 shows a linear fit of m st for T /J = 0.9, 1.0, 1.1.The linear decrease to zero as a function of N − 3 2 is almost perfect.The finite size dependence for T < T c is shown in the upper inset, where the data for T /J = 0.2, 0.4, 0.6 and a fit to the data with fit function aN b + c is presented.For all three temperatures, m st approaches a finite value.Since, however, the steepest drop in m st depends on the system size, this type of extrapolations can only be performed not too close to the critical temperature.
The above analysis was performed to demonstrate the existence of the two phases.With vanishing and non-vanishing staggered magnetization.It cannot be used as a reliable method to determine a critical exponent.We wish to postpone an analysis of critical exponents to a future project beyond the scope of this work, which includes an examination of spin-spin-correlation functions.A study of this type will hopefully give further insight in the interesting physics of three-dimensional quantum spin systems.Also we hope, that analysing systems with variable inter-plane couplings and anisotropies will be of experimental interest for S = 1 2 systems with honeycomb lattice structure, like InCu 2/3 V 1/3 O 3 .

Figure 3 .
Figure 3. Monte-Carlo results for the internal energy U (T ) per bond of coupled honeycomb planes as a function of temperature, for systems of different size.Upper inset: Specific heat of the same systems.Lower inset: Location of the peak as a function of 1/N with N = 6, 8, 10, 12, 16.

Figure 7 .
Figure 7.The staggered magnetization (mst) 2 for square systems of size N × N (squares) and rectangular system of size N/2×N (diamonds) and the correlations function at maximal distance C( N 2 , N 2 ) (circles) of the honeycomb lattice at T /J = 0.01 as a function of 1/N .The insets show the temperature dependence of (mst) 2 (left figure) and C(N/2, N/2) (right figure) for N = 16 (triangles), N = 24 (diamonds) and N = 32 (squares).