# Long range Néel order in the triangular Heisenberg model

\[

## Abstract

We have studied the Heisenberg model on the triangular lattice using several Quantum Monte Carlo (QMC) techniques (up to 144 sites), and exact diagonalization (ED) (up to 36 sites). By studying the spin gap as a function of the system size we have obtained a robust evidence for a gapless spectrum, confirming the existence of long range Néel order. Our best estimate is that in the thermodynamic limit the order parameter $`m^{}=0.41\pm 0.02`$ is reduced by about 59% from its classical value and the ground state energy per site is $`e_0=0.5458\pm 0.0001`$ in unit of the exchange coupling. We have identified the important ground state correlations at short distance.

\]

Historically the antiferromagnetic spin-1/2 Heisenberg model on the triangular lattice was the first proposed Hamiltonian for a microscopic realization of a spin liquid ground state (GS) :

$$\widehat{H}=J\underset{i,j}{}\widehat{𝐒}_i\widehat{𝐒}_j,$$
(1)   
where $`J`$ is the nearest-neighbors antiferromagnetic exchange and the sum runs over spin-$`1/2`$ operators. At the classical level the minimum energy configuration is the well known $`120^{}`$ Néel state. The question whether the combined effect of frustration and quantum fluctuations favors disordered gapped resonating valence bonds (RVB) or long range Néel type order is still under debate. In fact, there has been a considerable effort to elucidate the nature of the GS and the results of numerical , and analytical works are controversial. From the numerical point of view, ED, which is limited to small lattice sizes, provides a very important feature: the spectra of the lowest energy levels order with increasing total spin, a reminiscence of the Lieb-Mattis theorem for bipartite lattices, and are consistent with the symmetry of the classical order parameter . However, other attempts to perform a finite size scaling study of the order parameter indicate a scenario close to a critical one or no magnetic order at all.

The variational Quantum Monte Carlo (VMC) allows to extend the numerical calculations to fairly large system sizes, at the price to make some approximations, which are determined by the quality of the variational wavefunction (WF). Many WF have been proposed in the literature and the lowest GS energy estimation was obtained with the long range ordered type. In particular, starting from the classical Néel state, Huse and Elser  introduced important two and three spin correlation factors in the WF:

$$|\psi _\mathrm{V}=\underset{x}{}\mathrm{\Omega }(x)\mathrm{exp}\left(\frac{\gamma }{2}\underset{i,j}{}v(ij)S_i^zS_j^z\right)|x,$$
(2)   
where $`|x`$ is an Ising spin configuration specified by assigning the value of $`S_i^z`$ for each site and

$$\mathrm{\Omega }(x)=T(x)\mathrm{exp}\left[i\frac{2\pi }{3}\left(\underset{i\mathrm{B}}{}S_i^z\underset{i\mathrm{C}}{}S_i^z\right)\right]$$
(3)   
represents the three sublattices (say A, B and C) classical Néel state in the $`xy`$-plane multiplied by the three spin term

$$T(x)=\mathrm{exp}\left(i\beta \underset{i,j,k}{}\gamma _{ijk}S_i^zS_j^zS_k^z\right),$$
(4)   
defined by the coefficients $`\gamma _{ijk}=0,\pm 1`$, appropriately chosen to preserve the symmetries of the classical Néel state, and by an overall factor $`\beta `$ as discussed in Ref. . Since the Hamiltonian is real and commutes with the $`z`$-component of the total spin, $`\widehat{S}_{\mathrm{tot}}^z`$, a better variational WF on a finite size is obtained by taking the real part of Eq. (2) projected onto the $`S_{\mathrm{tot}}^z=0`$ subspace.

For the two body Jastrow potential $`v(r)`$ it is also possible to work out an explicit Fourier transform $`v_q`$, based on the consistency with linear spin wave (SW) results and a careful treatment of the singular modes coming from the $`SU(2)`$ symmetry breaking assumption. This analysis gives $`v_q=1\sqrt{1+2\gamma _q/1\gamma _q}`$ for $`q0`$ and $`0`$ otherwise, where $`\gamma _q=\left[cos(q_x)+2cos(q_x/2)cos(\sqrt{3}q_y/2)\right]/3`$ and the $`q`$-momenta are the ones allowed in a finite size with $`N`$-sites. For a better control of the finite size effects we have chosen to work with clusters having all the spatial symmetries of the infinite system .

In the square antiferromagnet (AF) the classical part by itself determines exactly the phases (signs) of the GS in the chosen basis, the so called Marshall sign. For the triangular case the exact phases are unknown and the classical part is not enough to fix them correctly. Therefore, one has to introduce the three-body correlations of Eq. (4). Although these do not provide the exact answer, they allow to adjust the signs of the WF in a non trivial way without changing the underlying classical Néel order. To this respect it is useful to define an average sign of the variational WF relative to the normalized exact GS $`|\psi _0`$ as

$$s=\underset{x}{}|\psi _0(x)|^2\mathrm{sgn}\left(\psi _\mathrm{V}(x)\psi _0(x)\right),$$
(5)   
with $`\psi (x)=x|\psi `$.

We have compared the variational calculation with the exact GS obtained by ED on the $`N=36`$ cluster. For completeness we have considered the more general Hamiltonian with exchange easy-plane anisotropy $`\alpha `$, ranging from the $`XY`$ case ($`\alpha =0`$) to the standard spin isotropic case ($`\alpha =1`$). As shown in Tab. I, in the variational approach the most important parameter, particularly for $`\alpha 1`$, is the one, $`\beta `$, controlling the triplet correlations. Though the overlap of our best variational WF with the exact GS is rather poor, the average sign $`s`$ is in general very much improved by the triplet term. Our interpretation is that short range many body correlations are very important to reproduce the relative phases of the GS on each Ising configuration. The optimal parameters for our initial guess $`\psi _\mathrm{V}`$ of the GS $`\psi _0`$ are expected to be very weakly size-dependent but they are very difficult to determine accurately for large sizes. For $`\alpha =1`$ and $`N=36`$, where ED is still possible, our best guess for the GS WF - with the maximum overlap and average sign - is slightly different from the one determined with the optimization of the energy. Since the forthcoming calculations, which significantly improve the VMC, are more sensitive to the accuracy of the WF rather than to the one of the GS energy, henceforth we have chosen to work with $`\beta =0.23`$ for all the system sizes.

One way to get more accurate GS properties is to use the Green Function MC technique (GFMC). As in the fermionic case, for frustrated spin systems this numerical method is plagued by the well-known sign problem. Recently, to alleviate the above mentioned instability, the Fixed-Node (FN) GFMC scheme  has been introduced as a variational technique, typically much better than the conventional VMC. As shown in Fig. 1, and also pointed in Ref. , for frustrated spin systems, this technique does not represent a significative advance compared to VMC, leading therefore to results biased by the variational ansatz.

In order to overcome this difficulty we have used a recently developed technique: GFMC with Stochastic Reconfiguration (SR) , which allows to release the FN approximation, in a controlled but approximate way, yielding, as shown in Fig. 1 much lower energies, even for the largest sizes where ED is not possible. In the appropriate limit  of large number of walkers and high frequency of SR, the residual bias introduced by the SR depends only on the number $`p`$ of operators used to constrain the GFMC Markov process. These constraints, analogously to the FN one, allow simulations without numerical instabilities. In principle the exact answer can be obtained, within statistical errors, provided $`p`$ equals the huge Hilbert space dimension. Practically it is necessary to work with small $`p`$ and an accurate selection of physically relevant operators is crucial. As can be easily expected, the short range correlation functions $`\widehat{S}_i^z\widehat{S}_j^z`$ and $`(\widehat{S}_i^+\widehat{S}_j^{}+\widehat{S}_i^{}\widehat{S}_j^+)`$ contained in the Hamiltonian give a sizable improvement of the FN GS energy when they are put in the SR procedure. In order to be systematic we have included in the SR the short range correlations generated by $`\widehat{H}^2`$ (see Fig. 2), averaged over all spatial symmetries commuting with the Hamiltonian. This local correlations are particularly important to obtain quite accurate and reliable estimates not only of the GS energy but also of the mixed average  of the total spin square $`\widehat{S}_{\mathrm{tot}}^2`$ and of the order parameter $`m^2`$ (defined as in Ref. ). These quantities are easily estimated within the GFMC technique and compared with the exact values computed by ED for $`N=36`$ in Tab. II. In particular it is interesting that, starting from a variational WF with no definite spin, the GS singlet is systematically recovered by means of the SR technique. Furthermore, as it is shown in Fig. 1, the quality of our results is similar to the variational one obtained by P. Sindzingre et al. , using a long range ordered RVB wavefunction. The latter approach is almost exact for small lattices, but the sign-problem is already present at the variational level, and the calculation has not been extended to high statistical accuracy or to $`N>48`$.

Having obtained an estimate for the GS energy, at least an order of magnitude more accurate than our best variational guess, it appears possible to obtain physical features, such as a gap in the spin spectrum, that are not present at the variational level. For instance in the frustrated $`J_1J_2`$ spin model, with the same technique and a similar accuracy, a gap in the spin spectrum was found in the thermodynamic limit, starting with a similar ordered and therefore gapless variational WF .

In the isotropic triangular AF, the gap to the first spin excitation is rather small. Furthermore, for the particular choice of the guiding WF (2), the translational symmetry of the Hamiltonian is preserved only if projected onto subspaces with total $`S_{\mathrm{tot}}^z`$ multiple of three. Such an $`S=3`$ excitation belongs to the low-lying states of energy $`E_S`$ and spin $`S`$ of the ordered quantum AF, behaving as $`E_SE_0S(S+1)/N`$. If instead $`E_SE_0`$ remains finite for $`S=3`$ and $`N\mathrm{}`$, this implies a disordered GS. For all the above reasons we have studied the gap to the spin $`S=3`$ excitation as a function of the system size. As it is shown in Fig. 3, for the lattice sizes for which a comparison with ED data is possible, the spin gap estimated with the SR technique is nearly exact. The importance to extend the numerical investigation to clusters large enough to allow a more reliable extrapolation is particularly evident in the same figure in which the $`N=12`$ and 36 exact data extrapolate linearly to a large finite value. This behavior, is certainly a finite size effect and it is corrected by the SR data for $`N48`$, suggesting, strongly, a gapless excitation spectrum.

As we have seen GFMC allows to obtain a very high statistical accuracy on the GS energy, but does not allow to compute directly GS expectation values $`\psi _0|\widehat{O}|\psi _0`$. A straightforward way is to perturb the Hamiltonian with a term $`\lambda \widehat{O}`$ , calculate the energy $`E(\lambda )`$ in presence of the perturbation and, by Hellmann-Feynman theorem, estimate $`\psi _0|\widehat{O}|\psi _0=dE(\lambda )/d\lambda |_{\lambda =0}`$ with few computations at different small  $`\lambda `$’s. A further complication for non exact calculations like the FN or SR, is that if the off-diagonal matrix elements $`O_{x^{},x}`$ of the operator $`\widehat{O}`$ (in the chosen basis) have the opposite sign of the product $`\psi _\mathrm{V}(x^{})\psi _\mathrm{V}(x)`$, they cannot be handled exactly within FN because these matrix elements change the nodes of $`\psi _\mathrm{V}`$. A way to circumvent this difficulty if to split the operator $`\widehat{O}`$ in three contributions: $`\widehat{O}=\widehat{D}+\widehat{O}^++\widehat{O}^{}`$, where $`\widehat{O}^+`$ ($`\widehat{O}^{}`$) is the operator with the same off-diagonal matrix elements of $`\widehat{O}`$ when they have the same (opposite) signs of $`\psi _\mathrm{V}(x^{})\psi _\mathrm{V}(x)`$, and zero otherwise, whereas $`\widehat{D}`$ is the diagonal part of $`\widehat{O}`$. Then we can add to the Hamiltonian a contribution that does not change the nodes: $`\widehat{H}(\lambda )=\widehat{H}\lambda (\widehat{D}+2\widehat{O}^+)`$ for $`\lambda >0`$ and $`\widehat{H}(\lambda )=\widehat{H}\lambda (\widehat{D}+2\widehat{O}^{})`$ for $`\lambda <0`$. It is easy to show that $`\underset{\lambda 0}{lim}(E(\lambda )E(\lambda ))/2\lambda =\psi _0|\widehat{O}|\psi _0`$.

We plot in Fig. 4 $`m^2`$ estimated with this method using the FN and SR techniques. For the order parameter the inclusion of many short range correlations in the SR is not very important (see Tab. II). Then, in order to minimize the numerical effort, we have chosen to put in the SR conditions the first four correlation functions shown in Fig. 2, the order parameter itself and $`\widehat{S}_{\mathrm{tot}}^2`$. While the FN data extrapolate to a value not much lower than the variational result, the SR calculation provides a much more reliable estimate of the order parameter with no apparent loss of accuracy with increasing sizes. In this way we obtain for $`\widehat{m}^{}`$ a value well below the linear and the second order (which has actually a positive correction ) SW predictions. This is partially in agreement with the conclusions of the finite temperature calculations  suggesting a GS with a small but nonzero long range AF order and with series expansions  indicating the triangular antiferromagnetic Heisenberg model to be likely ordered but close to a critical point. However in our simulation, which to our knowledge represents a first attempt to perform a systematic finite size scaling analysis of the order parameter, the value of $`\widehat{m}^{}`$ remains sizable and finite, consistent with a gapless spectrum. This features could be also verified experimentally on the K/Si(111):B interface  which has turned out recently to be the first realization of a really bidimensional triangular AF.

Though there is classical long range order, both the VMC and the SR approach show the crucial role of GS correlations defined on the smallest four spin clusters: in the variational calculation they are important to determine the correct relative phases of the GS WF whereas in the latter more accurate approach this correlations allow to obtain very accurate results for the energy and the spin gap and to restore the spin rotational invariance of the finite size GS.

Useful communications with M. Boninsegni and P. W. Leung are acknowledged. We also thank M. Calandra and A. Parola for help and fruitful discussions. This work was supported in part by INFM (PRA HTSC and PRA LOTUS), CINECA grant and CONICET (A.E.T.).