# Thermodynamics of the quantum easy-plane antiferromagnet on the triangular lattice

## Abstract

The classical XXZ triangular-lattice antiferromagnet (TAF) shows both an Ising and a BKT transition, related to the chirality and the in-plane spin components, respectively. In this paper the quantum effects on the thermodynamic quantities are evaluated by means of the pure-quantum self-consistent harmonic approximation (PQSCHA), that allows one to deal with any spin value through classical MC simulations. We report the internal energy, the specific heat, and the in-plane correlation length of the quantum XX0 TAF, for $`S=1/2`$, 1, 5/2. The quantum transition temperatures turn out to be smaller the smaller the spin, and agree with the few available theoretical and numerical estimates.

A renewed interest has recently focused on triangular antiferromagnets (TAF) . Indeed, they turned out to describe the magnetic behavior of several real compounds as, for example, the stacked antiferromagnet $`\mathrm{NaTiO}_2`$ , the organic superconductors of the family $`\kappa (\mathrm{BEDT}\mathrm{TTF})_2\mathrm{X}`$  and the K/Si(111):B interface .

In this paper we investigate the thermodynamic properties of the quantum XXZ Heisenberg antiferromagnet on the triangular lattice, defined by the following Hamiltonian

$$\widehat{}=\frac{J}{2}\underset{𝐢,𝐝}{}\left(\widehat{S}_𝐢^x\widehat{S}_{𝐢+𝐝}^x+\widehat{S}_𝐢^y\widehat{S}_{𝐢+𝐝}^y+\lambda \widehat{S}_𝐢^z\widehat{S}_{𝐢+𝐝}^z\right),$$
(1)   
where $`J`$ is the positive (antiferromagnetic) exchange constant, and $`(\widehat{S}_𝐢^x,\widehat{S}_𝐢^y,\widehat{S}_𝐢^z)`$ are the spin operators sitting on the sites $`𝐢`$ of a triangular lattice. They satisfy $`𝒮𝒰(2)`$ commutation relations and belong to the spin-$`S`$ representation $`\left|\widehat{𝐒}_𝐢^2\right|=S(S+1)`$. The interaction is restricted to nearest-neighbors and $`𝐝`$ runs over their relative displacements. The planar character of the system is due to the presence of the anisotropy $`\lambda [0,1)`$, energetically favoring configurations with the spins lying in the $`xy`$ plane (easy-plane). For $`\lambda =0`$ the spin components on the $`z`$ axis do not appear in the Hamiltonian and the model is known as XX0 or quantum XY.

The minimum energy configuration of the classical counterpart of the Hamiltonian (1) for every value of $`\lambda [0,1]`$ consists of coplanar spins forming $`\pm 2\pi /3`$ angles between nearest-neighbors and this leads to a $`\sqrt{3}\times \sqrt{3}`$ periodic Néel state. In contrast to the isotropic case, where the plane in which the $`2\pi /3`$ structure lies can take any direction in the spin space, in the XXZ model such structure must take place in the easy-plane. As a result, in the planar TAF the frustration causes an additional discrete two-fold degeneracy of the classical ground state, which is due to the chirality (or helicity), defined as the sign of rotation of the spins along the sides of each elementary triangle. The resulting degeneracy corresponds to the group $`SO(2)\times Z_2`$. As the Mermin and Wagner theorem only states that the sublattice magnetization must vanish at any non zero temperature, long-range order can occur as far as the chirality is concerned, and an Ising-like phase transition is indeed observed , in addition to the usual Berezinskii-Kosterlitz-Thouless (BKT) critical behavior associated to the rotation symmetry in the $`xy`$ plane.

In the quantum case the situation is far less clear. In fact, unlike the antiferromagnet on the square lattice where there is a general consensus about the ordered nature of the ground state even for $`S=1/2`$, in the frustrated cases the lack of exact analytical results is accompanied by difficulties in applying stochastic numerical methods, as their reliability is strongly limited by the well-known sign problem. Indeed only very recently a systematic size-scaling of the order parameter and of the spin gap has been performed using a new Quantum Monte Carlo technique , confirming the existence of Néel long-range order in the ground state as also suggested by the symmetry properties of the first excited states, evidenced by Bernu et al. .

An even less clear situation is that concerning the finite temperature behavior. In fact an early numerical work , limited to lattice sizes up to 27 sites, indicated for the $`S=1/2`$ XX0 model a phase diagram similar to the classical one, in contrast with the high temperature expansion produced by Fujiki and Betts  where no evidence for a phase transition was found. For the XXZ Hamiltonian, Momoi and Suzuki , applying an effective field theory, conjectured that the chiral phase transition should persist for every value of $`\lambda [0,1)`$, as in the classical case, and were able to estimate the transition temperature for $`\lambda =0`$, obtaining a value very close to that found in Ref. . Recently Suzuki and Matsubara  using a quantum transfer Monte Carlo method to study clusters up to 24 sites, have claimed instead the absence of the chiral order at any finite temperature for $`\lambda 0.6`$.

In this context where, at least up to now, quantum numerical methods cannot give satisfactory answers, the pure-quantum self-consistent harmonic approximation (PQSCHA)  can provide an effective instrument to investigate the thermodynamics of quantum spin systems, as far as their ground state is ordered. The method is based on the path-integral formulation of quantum statistical mechanics, and has been successfully applied recently to a variety of unfrustrated spin models, both one-  and two-dimensional .

By the PQSCHA the evaluation of thermal averages in the quantum model can be reduced to the calculation of classical-like averages over a Boltzmann distribution defined by an effective Hamiltonian, which contains the contribution of the pure-quantum part of the fluctuations (approximated within a self consistent harmonic scheme) in its renormalized interaction parameters, which are temperature and spin dependent. As a result one can get accurate results on the quantum spin system using classical computational methods, like the transfer-matrix in the one-dimensional case and classical Monte Carlo simulations in the two-dimensional one.

The first step of the derivation of the effective Hamiltonian for the easy-plane TAF, is to apply the unitary transformation which defines a spatially varying coordinate system pointing along the local Néel direction, namely

$$𝒰=\mathrm{exp}\left(\frac{2\pi }{3}i\underset{𝐢\mathrm{B}}{}\widehat{S}_𝐢^z\frac{2\pi }{3}i\underset{𝐢\mathrm{C}}{}\widehat{S}_𝐢^z\right),$$
(2)   
where B and C labels two of the three sublattices. Unlike in the bipartite lattices where the corresponding transformation maps the antiferromagnet into a model with an in-plane ferromagnetic exchange and an antiferromagnetic coupling along the $`z`$ axis, thus allowing the demonstration of the Lieb and Mattis theorem  and computability with standard quantum Monte Carlo methods, in the triangular case the transformed Hamiltonian shows an extra current-like term  which contains the effects of the frustration, whose form is quite similar to the chiral order parameter , i.e., the physical quantity undergoing the order-disorder phase transition present in the classical model for every value of $`\lambda <1`$.

From now on the derivation follows the same lines already described in Refs. . A point worth being recalled is the use of the Villain transformation  in order to represent the spin operators in terms of canonically conjugated variables, which is a necessary step in the derivation of the effective Hamiltonian. As it is well known, this spin-boson transformation preserves the commutation rules but neglects the so called kinematic interaction due to the limited spectrum of $`S^z`$, thus giving a better description when the system has a good easy-plane character and the spin states with large fluctuations of $`S^z`$ are less relevant to the thermodynamics. In the square lattice case  such approximation scheme turns out to be reliable up to some value of $`\lambda _\mathrm{M}<1`$ ($`\lambda _\mathrm{M}=0.58`$ in the extreme quantum case $`S=1/2`$), when the mapping with the Villain transformation breaks down and a different spin-boson transformation is needed. However, it provides accurate results for the critical temperatures even for $`\lambda =0.5`$; a similar behavior is also found in the case of the quantum TAF. Finally we remind that Weyl ordering, which is inherent to the PQSCHA, naturally leads to define an effective classical spin length as $`\stackrel{~}{S}=S+1/2`$ and thus to set the natural energy scale $`ϵ=J\stackrel{~}{S}^2`$. Therefore in the following we use the reduced temperature, $`t=k_BT/J\stackrel{~}{S}^2`$.

In the case of the easy-plane TAF the effective Hamiltonian has the form $`_{\mathrm{eff}}=\overline{}+G(t)`$, where $`G(t)`$ is an additive uniform term, formally identical to that obtained for the square lattice, which is unessential in the calculation of the thermal averages, while:

$$\overline{}=\frac{ϵ}{2}j_{\mathrm{eff}}\underset{𝐢,𝐝}{}\left(s_𝐢^xs_{𝐢+𝐝}^x+s_𝐢^ys_{𝐢+𝐝}^y+\lambda _{\mathrm{eff}}s_𝐢^zs_{𝐢+𝐝}^z\right),$$
(3)   
where $`(s_𝐢^x,s_𝐢^y,s_𝐢^z)`$ are unit vectors, i.e., classical spins. Within the PQSCHA quantum effects are embodied in the temperature and spin dependence of the renormalized

$`j_{\mathrm{eff}}(t,S,\lambda )`$ $`=`$ $`(1{\displaystyle \frac{1}{2}}D_{})^2e^{\frac{1}{2}𝒟_{}},`$ (4)   
$`\lambda _{\mathrm{eff}}(t,S,\lambda )`$ $`=`$ $`\lambda (1{\displaystyle \frac{1}{2}}D_{})^1e^{\frac{1}{2}𝒟_{}},`$ (5)   
with

$`D_{}`$ $`=`$ $`(2\stackrel{~}{S}N)^1{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{b_𝐤}{a_𝐤}}(f_𝐤),`$ (6)   
$`𝒟_{}`$ $`=`$ $`(\stackrel{~}{S}N)^1{\displaystyle \underset{𝐤}{}}(1\gamma _𝐤){\displaystyle \frac{a_𝐤}{b_𝐤}}(f_𝐤),`$ (7)   
where

$`a_𝐤^2`$ $`=`$ $`{\displaystyle \frac{z}{2}}(1\frac{1}{2}D_{})\mathrm{e}^{{\scriptscriptstyle \frac{1}{2}}𝒟_{}}(1+2\lambda _{\mathrm{eff}}\gamma _𝐤),`$ (8)   
$`b_𝐤^2`$ $`=`$ $`{\displaystyle \frac{z}{2}}j_{\mathrm{eff}}(1\gamma _𝐤),`$ (9)   
$`f_𝐤=a_𝐤b_𝐤/(2\stackrel{~}{S}t)`$ , $`(x)=\mathrm{coth}xx^1`$ is the Langevin function, $`\gamma _𝐤=z^1_𝐝\mathrm{cos}(𝐤𝐝)`$ , and $`𝐤`$ is a wavevector varying in the first Brillouin zone. $`D_{}(S,\lambda ,t)`$ and $`𝒟_{}(S,\lambda ,t)`$ represent the pure-quantum square fluctuations of the out-of-plane and in-plane components of the spins respectively. They are decreasing functions of temperature and spin, vanishing for $`t\mathrm{}`$ and $`S\mathrm{}`$, i.e., when the quantum part of the fluctuations is negligible with respect to the classical one.

From the above equations, we can infer that the PQSCHA approach is valid under the condition that second order terms in $`(D_{})^2`$ can be neglected. One can take the criterion that the renormalization effects of quantum fluctuations must not reduce much more than, say, $`50\%`$ the effective exchange integral. Such a strong renormalization only occurs for $`S=1/2`$ and $`t0.2`$, while for higher spin values the PQSCHA is reliable at any temperature.

In the XX0 model, $`\lambda _{\mathrm{eff}}=\lambda =0`$ and all the information about the quantum system is hence contained in the renormalization of the energy scale. In this case the critical properties of the quantum system at a temperature $`t`$ are essentially those of its classical counterpart at the effective temperature $`t_{\mathrm{eff}}=t/j_{\mathrm{eff}}(t,S)`$, and we have used the results of classical Monte Carlo (MC) simulations recently obtained for lattice sizes up to $`N=120\times 120`$  to calculate the corresponding quantum observables.

Indicating the classical averages with the effective Hamiltonian as $`\mathrm{}_{\mathrm{eff}}`$, the internal energy per spin can be calculated as

$$e(t,S,\lambda )=\frac{\widehat{}}{Nϵ}=\overline{}_{\mathrm{eff}}+\frac{z}{2}\lambda 𝒟_{},$$
(10)   
where $`𝒟_{}`$ can be expressed as $`D_{}`$, Eq. (6), with an extra factor $`\gamma _𝐤`$ in the summand. For $`\lambda =0`$ the above equation reduces to

$$e(t,S)=j_{\mathrm{eff}}(t)e_{\mathrm{cl}}(t_{\mathrm{eff}}),$$
(11)   
being $`e_{\mathrm{cl}}(t)`$ the internal energy per spin of the corresponding classical system. In Fig. 1 $`e(t,S)`$ is plotted for various values of the spin in the range of temperatures where the PQSCHA is expected to give reliable results.

The energy curves flatten and increase with decreasing $`S`$ due to the increased quantum fluctuations. As said before the $`S=1/2`$ curve is reported only in the valid temperature range. As a matter of fact the extrapolation to lowest temperatures gives the self-consistent spin-wave ground-state energy. The difference from the most refined estimates  can be mainly attributed to $`1/S^2`$ constant contributions coming from the Villain transformation  and to the use of the low coupling approximation (LCA) . This term is not significant for $`S1`$.

Consistently with this picture the finite size ($`N=120\times 120`$) peak of the specific heat (Fig. 2), obtained by numerical derivation of the internal energy, moves towards lower temperatures and decreases in height as $`S`$ decreases. However, since the quantum renormalizations are essentially size-independent, classical scaling with size  is conserved and a logarithmic divergence of the specific heat, connected with the Ising-like chirality phase-transition, is therefore expected in the thermodynamic limit. By direct derivation of Eq. (11) it is easily seen that, in order for the quantum specific heat to vanish in the zero temperature limit, within our approximation we must have $`dj_{\mathrm{eff}}/dt|e_{\mathrm{cl}}(0)|^1`$ as $`t0`$ a condition which is fulfilled for every value of $`S`$, as can be verified analytically from the explicit expressions of the renormalization parameters; this is shown in the inset of Fig. 1.

Most papers in literature do mainly concern with the chiral order-disorder transition, while the XXZ TAF also supports another kind of phase transition. In fact, the classical system  displays as well a BKT critical behavior. For this reason we have calculated the magnetic correlation length which governs the decay of the in-plane correlation functions in the high-temperature phase, whose expression within the PQSCHA reads

$$\widehat{S}_𝐢^x\widehat{S}_{𝐢+𝐫}^x+\widehat{S}_𝐢^y\widehat{S}_{𝐢+𝐫}^y=G(𝐫)s_𝐢^xs_{𝐢+𝐫}^x+s_𝐢^ys_{𝐢+𝐫}^y_{\mathrm{eff}},$$
(12)   
where $`𝐢`$ and $`𝐢+𝐫`$ belong to the same sublattice and $`G(𝐫)`$ is bounded and essentially constant for large $`𝐫`$. As a consequence the asymptotic behavior of the correlation functions in the critical region is the same of the effective classical spin system. In particular, for $`\lambda =0`$, the correlation length can be simply found as $`\xi (t)=\xi _{\mathrm{cl}}(t_{\mathrm{eff}})`$. The result is reported in Fig. 3: as expected in a BKT transition, it displays a divergence at a temperature $`t_{\mathrm{BKT}}`$ which decreases with decreasing $`S`$, as a result of enhanced quantum fluctuations.

Within the PQSCHA, the quantum renormalizations cannot modify the critical behavior of the effective classical system , so that both the chirality and the BKT critical temperatures can be connected to their classical counterparts by the self-consistent relation

$$\frac{t_{\mathrm{crit}}(S,\lambda )}{t_{\mathrm{crit}}^{\mathrm{cl}}\left(\lambda _{\mathrm{eff}}(t_{\mathrm{crit}},S,\lambda )\right)}=j_{\mathrm{eff}}(t_{\mathrm{crit}},S,\lambda ),$$
(13)   
which can be solved numerically. The obtained critical temperatures for $`\lambda =0`$ and $`\lambda =0.5`$ are reported for various values of the spin in Table I. In the $`S=1/2`$ case, although our theory begins to become unreliable when $`t0.2`$, we notice that for $`\lambda =0`$ the extrapolated value for the chiral critical temperature, $`t_\mathrm{c}=0.193`$, agrees remarkably well with those obtained by the size scaling on the QMC  data ($`t_\mathrm{c}=0.195(1)`$) and the effective field theory of Ref.  ($`t_\mathrm{c}0.20`$).

L. C. acknowledges S. Sorella for continuous and fruitful discussions.