format: text

text:
# An investigation of the quantum 𝐽₁-𝐽₂-𝐽₃ model on the honeycomb lattice

## 1 Introduction

Frustrated quantum antiferromagnetic (AF) spin systems on low dimension (D) lattices have attracted a great deal of interest in recent years. Quantum fluctuations, largest for small values of the spin $`S`$ of the magnetic ions, low D and small coordination of the lattice, are expected to lead to novel magnetic behaviors. Their effects, have been preeminently seen in 1D. They have been investigated on a few 2D systems. The most studied models are the AF Heisenberg model on the triangular blp92 ; capriotti99  or kagomé lattice smlbpwe00  which are geometrically frustrated systems, the $`J_1J_2`$ model on the square lattice capriotti00 ; singh99 ; kos99  where frustration is introduced by 2nd neighbor interaction, the $`J_1J_2`$ model cj92 ; k93 ; lblp95a  and the multi-spin exchange model (MSE) mblw98 , on the triangular lattice.

Less studied oitmaa78 ; reger89 ; oitmaa92 ; weihong91 , spin models on the honeycomb lattice deserve attention due to the special properties of the lattice and since there are experimental realizations. A first feature of the lattice is that, like the square lattice, it is not geometrically frustrated for AF nearest neighbor interactions but has lower coordination. Thus quantum fluctuations are expected to be larger than for the square lattice. For this reason the spin-1/2 Heisenberg antiferromagnet on the honeycomb lattice, has been studied theoretically by various methods oitmaa78 ; reger89 ; oitmaa92 ; weihong91  which all predicted that Néel long range order (LRO) subsists but with an order parameter smaller than for the square lattice case. This also motivated a Schwinger-boson study of the effects of frustrating second neighbor interactions in the $`J_1J_2`$ model mattsson94 .

A major incentive to study frustrated magnets on the honeycomb lattice is the availability of experimental data in the family of compounds $`\mathrm{BaM}_2(\mathrm{XO}_4)_2`$ (M= Co, Ni; X= P, As) obtained, some years ago, by Regnault and Rossat-Mignod regnault-1 . The magnetic ions M have small spins (it is supposed to be $`S=1/2`$ for the Co oxide and $`S=1`$ for Ni), disposed in weakly coupled layers where they sit on a honeycomb lattice. The simplest model relevant to these quasi 2D compounds is a $`J_1J_2J_3`$ model on a honeycomb lattice with first, second and third neighbor interactions and either on site if ($`S=1`$) or XXZ (if $`S=1/2`$) anisotropy.

So far the $`J_1J_2J_3`$ model was only investigated within first order linear spin-wave theory (LSW) regnault-1 ; rastel79 , and to our knowledge the renormalization of the order parameter by quantum fluctuations has not be calculated even in this simplest approach. The experimental results motivated us to do this calculation in the large S limit and then attack the $`S=1/2`$ problem with exact diagonalizations (ED).

In this paper, as a first step, we consider the case of purely isotropic interactions. The Hamiltonian of the model reads:

$$=J_1\underset{<i,j>_1}{}𝐒_i.𝐒_j+J_2\underset{<i,k>_2}{}𝐒_i.𝐒_k+J_3\underset{<i,k>_3}{}𝐒_i.𝐒_k$$
(1)
where the first, second and third sums run on the first, second and third neighbor pairs of spins, respectively (see Fig.1). The coupling constants $`J_i`$ can be either AF ($`J_i>0`$) or ferromagnetic ($`J_i<0`$). Depending of the values of the parameters $`J_i`$, this model displays various classical ground-states: a collinear AF ground-state, two degenerate manifolds of planar helimagnetic ground-states with four or eight sublattices and a ferromagnetic ground-state. The classical phase diagrams of the isotropic and anisotropic XXZ models display only minor differences. In particular the classical ground-states are the same for the parameters believed to be relevant to $`\mathrm{BaCo}_2(\mathrm{AsO}_4)_2`$. In this paper we concentrate on the quantum effects in the isotropic model. The study of the quantum XXZ $`J_1J_2J_3`$ model with parameters appropriate to $`\mathrm{BaCo}_2(\mathrm{AsO}_4)_2`$ will be presented in a separate paper fouet00b .

We investigated quantum effects for selected values of the $`J_i`$ chosen to give a broad picture of the quantum effects encountered in the model. The restriction to isotropic interactions limits the number of parameters and enable to separate the effects of anisotropy. In addition to the $`J_1J_2`$ models on the square and triangular lattices, the present model may be compared with the $`J_1J_2J_3`$ model on the square lattice which has a similar variety of classical ground-states and to which a few studies have been devoted fkksrr90 ; moreo90 ; read-sachdev91 ; ceccatto93 ; ferrer93 ; ll96 .

This paper is divided into five parts. In section II, we recall the classical phase diagram of the model, obtained by Rastelli et al. rastel79 , and identify the degeneracies of the classical ground-states not considered by these authors, we also discuss the stability of the first order spin-wave approximation for these different phases. The ED results for the case of antiferromagnetic and ferromagnetic nearest neighbor coupling are presented in section III and IV respectively. In section V we draw conclusions, and enumerate some conjectures relative to the appearance of the various generic two-dimensional spin-liquids. We described in an appendix the various technical features specific to our present ED calculations on different samples.

## 2 Classical phase diagram and semi-classical deviations

### 2.1  Planar ground-state configurations

The classical model was studied by Rastelli et al. rastel79 . They searched for planar or uniformly canted configurations minimizing the classical energy $`E_{cl}`$. The former were found energetically favored over the latter. They represent spiral configurations, characterized by a wave-vector $`𝐐`$. The classical spin (of length $`S`$) sitting at cell $`𝐑`$ of the triangular Bravais lattice on sublattice $`\alpha `$ is given by:

$$𝐒_{𝐑,\alpha }=S(\mathrm{cos}(𝐐.𝐑+\varphi _\alpha )𝐮+\mathrm{sin}(𝐐.𝐑+\varphi _\alpha )𝐯)$$
(2)
where $`𝐮`$ and v are two orthogonal unit vectors defining the plane of the spins, $`\varphi _\alpha `$ can be chosen to be zero on one sublattice and will be noted $`\varphi `$ on the other. The set of spiral wave-vectors $`𝐐`$ minimizing the classical energy will be noted $`\{𝐐\}`$.

The phase diagram of planar solutions of this type is reproduced in Fig.2 ($`J_1>0`$) and Fig.3 ($`J_1<0`$). There is a mapping between the two phase diagrams: the transformation $`J_1J_1`$, $`J_3J_3`$ and $`𝐒_i𝐒_i`$ for $`i`$ on the black triangular sublattice of Fig.1 leaves the Hamiltonian unchanged, and maps the ground-state for $`J_1>0`$ on that for $`J_1<0`$.

There are six regions in each phase diagram: four collinear phases (I,II,IV,VI) and two spiral ones (III,V). In the collinear regions I and VI, the wave-vector of the magnetic order is $`𝐐=0`$, whereas in regions II and IV, $`\{𝐐\}=\{𝐊_𝐢\}`$, where $`𝐊_i`$ are the three inequivalent middles of edges of the Brillouin zone (see Fig.6). The phases are: In I, $`\varphi =\pi `$ ($`0`$) if $`J_1>0`$ ($`J_1<0`$), in VI, $`\varphi =0`$ ($`\pi `$) if $`J_1>0`$ ($`J_1<0`$), in II $`\varphi =\pi `$ ($`0`$) if $`J_1>0`$ ($`J_1<0`$), in IV $`\varphi =0`$ ($`\pi `$) if $`J_1>0`$ ($`J_1<0`$).

The separation lines I-III, I-V, II-III, IV-V represent continuous transitions; the others are first order phase transitions.

### 2.2 Non planar ground-states manifolds

Ansatz (Eq.2) is usually generic to find all allowed classical ground-states villain77 ; bertaut63 . It assumes that, up to the trivial degeneracy associated to a global spin rotation, the ground-state is unique or exhibits at most a discrete degeneracy. It is valid if a linear combination of the different $`𝐐`$ modes of the same set $`\{𝐐\}`$, can be excluded as violating the constraint $`|𝐒_{𝐑,\alpha }|=S`$ on every site. Exceptions occur for special sets $`\{𝐐\}`$, in particular if $`𝐐`$ is half or one fourth of a reciprocal lattice vector $`𝐆`$ villain77 . This is the case in region II and IV where $`\{𝐐\}=\{𝐆/\mathrm{𝟐}\}=\{𝐊_𝐢\}`$ (see Fig.6). Here, the linear combination of three $`𝐊_i`$ solutions

$$𝐒_{𝐑,\alpha }=\underset{i=1}{\overset{3}{}}S\mathrm{cos}(𝐊_i.𝐑+\varphi _\alpha )𝐮_i$$
(3)
with unnormalized $`𝐮_i`$, is submitted to the constraints $`𝐮_i.𝐮_j=\delta _{i,j}`$ and $`𝐮_i^2=1`$. The ground-state is a two dimensional manifold continuously connecting the three $`𝐊_i`$ solutions (there are nine degrees of freedom for choosing the three $`𝐮_i`$, minus three for global rotations of the spins and four constraints). This gives birth to the non planar ground-states manifolds described below.

In regions IV for $`J_1>0`$ or II for $`J_1<0`$ ( $`\varphi =0`$), the ground-state manifold is the set of four-sublattice ordered solutions such as $`𝐒_A+𝐒_B+𝐒_C+𝐒_D=0`$. This could be seen directly from the expression of the classical energy of these configurations:

$`E_{cl}`$ $`=`$ $`{\displaystyle \frac{2}{N}}\left(J_1+2J_2\right)\left(𝐒_A+𝐒_B+𝐒_C+𝐒_D\right)^2`$
$``$ $`{\displaystyle \frac{2}{N}}\left(J_1+2J_23J_3\right)\left(𝐒_A^2+𝐒_B^2+𝐒_C^2+𝐒_D^2\right).`$
In this equation, N represents the total number of spins of the sample. The generic four-sublattice configurations are shown in Fig.4, as well as the three collinear configurations, which appear as special cases of it, with $`𝐒_A=𝐒_B=𝐒_C=𝐒_D`$ and the two other combinations of parallel spins.

The situation is reminiscent of the $`J_1J_2`$ model on the triangular lattice cj92 ; k93 ; lblp95a  with $`1/8<J_2/J_1<1`$.

In regions IV for $`J_1<0`$ or in II for $`J_1>0`$ ($`\varphi =\pi `$), the ground-state manifold is the set of eight-sublattice solutions shown in Fig.5 where sublattices labelled by the same letter are paired (partner sublattice is over-headed by a bar) and

$$𝐒_{\overline{\alpha }}=𝐒_\alpha $$
(5)
for $`\alpha \{A,B,C,D\}`$ and

$$𝐒_A+𝐒_B+𝐒_C+𝐒_D=0.$$
(6)
This minimizes the classical energy:

$`E_{cl}`$ $`=`$ $`{\displaystyle \frac{8}{N}}J_2(𝐒_A+𝐒_B+𝐒_C+𝐒_D`$ (7)
$`+𝐒_{\overline{A}}+𝐒_{\overline{B}}+𝐒_{\overline{C}}+𝐒_{\overline{D}})^2`$
$`+`$ $`{\displaystyle \frac{8}{N}}\left(J_12J_2\right)\left(𝐒_A+𝐒_B+𝐒_C+𝐒_D\right)`$
$`\left(𝐒_{\overline{A}}+𝐒_{\overline{B}}+𝐒_{\overline{C}}+𝐒_{\overline{D}}\right)`$
$`+`$ $`{\displaystyle \frac{4}{N}}(3J_3J_1)[(𝐒_A+𝐒_{\overline{A}})^2+(𝐒_B+𝐒_{\overline{B}})^2`$
$`+(𝐒_C+𝐒_{\overline{C}})^2+(𝐒_D+𝐒_{\overline{D}})^2]`$
$``$ $`{\displaystyle \frac{4}{N}}(J_1+2J_2+3J_3)(𝐒_A^2+𝐒_B^2+𝐒_C^2+𝐒_D^2`$
$`+𝐒_{\overline{A}}^2+𝐒_{\overline{B}}^2+𝐒_{\overline{C}}^2+𝐒_{\overline{D}}^2).`$
It is highly probable that thermal fluctuations will stabilize the collinear solutions as it does in similar situations on the square and triangular lattice. We will show below that quantum fluctuations indeed do it.

Continuous degeneracy of the ground-state also occurs when $`𝐐=𝐆/4`$ in III and in V but since this happens only on lines (for instance if $`J_2=0.5`$ in III) and not in full regions, we shall skip their description which is not essential to our present goal. The transition line III-V between the two spiral regions is very special. It has an infinite degeneracy of spiral ground-states corresponding to:

$$\mathrm{cos}(𝐐.𝐭_1)+\mathrm{cos}(𝐐.𝐭_2)+\mathrm{cos}(𝐐.(𝐭_1𝐭_2))=1/8J_2^23/2$$
(8)
The lines of $`𝐐`$ solutions of Eq.8 are shown in Fig.6 for $`J_2=0.2,0.4,0.5`$.

### 2.3 Stability of the quasi-classical phase diagram in the large S limit: LSW results

The renormalization of the order parameter $`m^{}`$ in the first order spin-wave approximation is already large in the pure model ($`J_1=1,J_2=0,J_3=0`$): $`m^{}0.48`$ i.e. a value reduced to $`48\%`$ of its classical value (Fig.7).

For $`J_1>0`$, the interplay of quantum fluctuations and frustration quickly destroys Néel LRO: $`m^{}`$ goes to zero for $`J_20.1`$ or $`J_30.1`$. The helimagnetic LRO of zone V disappears: near the $`J_3=0`$ axis, this is mainly due to the large classical degeneracy of the ground-state and near the $`J_2=0`$ axis, the main cause is the vanishing of the spin-wave velocity at the point $`J_1=1,J_2=0,J_3=0.25`$. The zone V being very small, we conclude that the Néel helimagnetic ground-state does not survive in region V for antiferromagnetic $`J_1`$.

For $`J_1<0`$ the ground-state is ferromagnetic in zone I of figure 3. The ferromagnetic state is an exact eigenstate of the hamiltonian, quantum fluctuations don’t destroy it. However, the classical degeneracy on the boundary between zones III and V implies a whole branch of soft modes and a disappearance of the helimagnetic LRO, on and near this line (as in the AF $`J_1`$ case). The main difference with the AF first neighbor case is the persistence of Néel order in zone V near the $`J_1=1,J_2=0,J_3=0.25`$ point and in the vicinity of the $`J_2=0`$ axis.

The nature of the quantum spin-$`1/2`$ phase for small $`J_2`$ and $`J_3`$ appears an open problem that we will now attack with the help of exact diagonalizations (ED).

Exact diagonalizations were performed on samples of $`N=18,24,26,28,30,32`$ sites with appropriate boundary conditions (see Appendix). The technical problems encountered in such approaches have already been studied in details in previous references blp92 ; mblw98  and will not be described in details in this paper. We briefly discuss in Appendix the different characteristics of the studied samples with respect to the present model. We will now proceed to the analysis of the ED results.

## 3 $`J_1>0`$: AF nearest neighbor interactions

### 3.1 The purely Heisenberg model, and phase I of the quantum model

The AF Heisenberg point ($`J_1=1`$, $`J_2=J_3=0`$) was previously investigated by exact diagonalization calculations on small samples oitmaa78 , Monte Carlo simulations reger89 , series expansions around the Ising limit oitmaa92 , spin-wave theory up to second order weihong91 , and Schwinger-boson mean field theory mattsson94 . All concluded that the quantum system exhibits Néel LRO but with a large reduction of the order parameter due to quantum fluctuations (the Monte Carlo result, rather close to the spin-wave value, is $`m^{}=0.44\pm 0.06`$). Our ED results for sizes up to $`N=32`$ are consistent with this conclusion. The approach is however different. In antiferromagnets with linear Goldstone modes, the scaling law for $`m^{}`$ is an expansion in $`1/N^{1/2}`$ fisher89 ; hasen93 . For the sizes encountered in ED, the asymptotic $`1/N^{1/2}`$ law is never reached blp92  and the extrapolation of the order parameter to the thermodynamic limit remains uncertain. A qualitative idea of the finite size scaling of $`m^{}`$ can be obtained from the LSW results: they show that the asymptotic $`1/N^{1/2}`$ regime cannot be expected for sizes smaller than $`N400`$ (see Fig. 7).

Nevertheless confirmation of Néel LRO can be obtained thanks to characteristic features of the spectrum itself which have more favorable scaling behaviors:

* For a given sample size, the lowest eigenlevels in each sector of the total spin $`S`$ evolve as $`E_0(S,N)S(S+1)`$ up to $`S\sqrt{N}`$, as shown in Fig.8. They are the eigenlevels associated with the collective dynamics of the order parameter, the so-called Quasi Degenerate Joint States (QDJS)blp92 , which can be described by the effective Hamiltonian:

$$H_{coll.dyn.}=\frac{1}{2}\mathrm{\Delta }(N)𝐒^{\mathrm{\hspace{0.17em}2}}$$
(9)
where $`\mathrm{\Delta }(N)`$ is the finite-size difference in total energy between the absolute ground-state and the first triplet excitation.
* The number of QDJS and their symmetries are those expected for the projections of the classical Néel order on the irreducible representations (IR) of $`SU(2)𝒢`$ ($`𝒢`$:lattice symmetry group): There is just one state for each $`S`$ value since there is just one way to couple two spins of magnitude $`N/4`$ in a total spin $`S`$. These states are invariant under lattice translations (their wave-vector is $`𝐤=0`$), under $`2\pi /3`$ rotations around the center of an hexagon, they are even (odd) under inversion with respect to the center of an hexagon for even $`S`$ and $`N=4p`$ (respectively odd $`S`$ and $`N=4p+2`$).
* In the thermodynamic limit, the QDJS collapse on the singlet ground-state as $`1/N`$ (Fig. 10b). The QDJS remain distinct from the softest magnon excitation which collapse to the ground-state as $`1/N^{1/2}`$. The QDJS and the softest magnon are shown for the $`N=32`$ sample in Fig.9.
* The asymptotic $`1/N`$ behavior of $`\mathrm{\Delta }(N)`$ is not yet reached for our largest sizes as seen in Fig.10b. The next order term in the $`1/N^{1/2}`$ expansion reads fisher89 ; hasen93 :

$$\mathrm{\Delta }(N)=\frac{1}{4\chi N}(1\beta \frac{c}{\rho \sqrt{N}})+𝒪(\frac{1}{N^2})$$
(10)
where $`\chi `$ is the spin susceptibility, $`c`$ is the spin-wave velocity, $`\rho `$ the spin stiffness and $`\beta `$ is a number of order one. A fit of the spin-gaps to this scaling law is shown in Fig.10b. The importance of the term in $`1/N^{3/2}`$ is not unexpected since this term is $`c/\rho `$ and quantum fluctuations which reduce $`\rho `$ with respect to its classical value are strong (as already shown by the reduction of the order parameter).
* The energy per site $`e_0(N)=E_0(0,N)/N`$ of the ground-state scales as:

$$e_0(N)=e_{\mathrm{}}\frac{\alpha ^{}}{N^{3/2}}(1\beta ^{}\frac{c}{\rho \sqrt{N}})+𝒪(\frac{1}{N^{5/2}}).$$
(11)
for a Néel order. Fig.10a shows that the leading term of order $`𝒪(1/N^{3/2})`$ is enough to describe the size effects in this range.

A rapid analysis of the quantum phase diagram in region I does not reveal new phases, but both LSW calculations and ED confirm that a weak antiferromagnetic second or ferromagnetic third neighbor coupling are sufficient to kill LRO: the boundary between phase I and phases III and V is shifted upwards by quantum effects.

### 3.2 Region IV

In region IV the classical model presents a degenerate manifold of four-sublattice ordered ground-states. The finite-size spectra clearly show that this degeneracy is lifted by quantum fluctuations which favor a collinear two-sublattice order (see Fig.11): the low lying levels of these spectra, below the magnons excitations, exhibit a large family $`\{^4\stackrel{~}{E}\}`$ of QDJS associated to four-sublattice solutions. At the bottom of this family there appears a line of eigenlevels with definite symmetries: these levels constitute the family $`\{^2\stackrel{~}{E}\}`$ of QDJS states associated to a collinear symmetry breaking.

The situation, seen here, is very similar to the one previously studied for the $`J_1J_2`$ antiferromagnet on the triangular lattice lblp95a . The expected number of states $`{}_{}{}^{4}N_{S}^{}`$, $`{}_{}{}^{2}N_{S}^{}`$ in $`\{^4\stackrel{~}{E}\}`$ and $`\{^2\stackrel{~}{E}\}`$ and their space symmetries are easily determined for each value of the total spin $`S`$. The eigenstates of $`\{^4\stackrel{~}{E}\}`$ can be labelled by the five irreducible representations $`\mathrm{\Gamma }_i`$ ($`i=1,5`$) of $`S_4`$ (permutation group of four elements). The mapping between the space group operations on the four-sublattice solutions and permutations of $`S_4`$ is described in Table 1, together with the character Table of $`S_4`$. The four-sublattice order is invariant in two-fold rotations ($`_\pi `$): thus the eigenstates of $`\{^4\stackrel{~}{E}\}`$ belong to the trivial representation of $`C_2`$. Since it is also invariant under a two-step translation of the Bravais lattice they have either a wave-vector $`𝐤=\mathrm{𝟎}`$ or a wave-vector $`𝐊_i`$. $`\mathrm{\Gamma }_1,\mathrm{\Gamma }_2`$ and $`\mathrm{\Gamma }_3`$ belong to the $`𝐤=\mathrm{𝟎}`$ subspace, whereas $`\mathrm{\Gamma }_4`$ and $`\mathrm{\Gamma }_5`$ belong to the subspace $`\{𝐊_i\}`$. $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ are invariant under the three-fold rotations $`_{2\pi /3}`$ of $`C_3`$, whereas $`\mathrm{\Gamma }_3`$ is associated with the two-dimensional representation of $`C_{3v}`$. The number of replicas of $`\mathrm{\Gamma }_i`$ that should appear for each $`S`$ value can be computed as in ref lblp95a . The result are shown in Table 2 for the N=32 sample. Analysis of the two-sublattice order can be done similarly: the collinear solution has a three-fold degeneracy, the set of eigenstates $`\{^2\stackrel{~}{E}\}`$ maps on $`Z_3`$. It is characterized by the IR $`\mathrm{\Gamma }_1`$, $`\mathrm{\Gamma }_3`$ and $`\mathrm{\Gamma }_4`$. The number of repliquas are shown in Table 3 for the N=32 sample.

The “order out of disorder” phenomenon villain77  is clearly seen for $`J_1=1,J_2=0,J_3=2`$. In Fig.11 we show the lower part of the $`N=32`$ spectrum at this point. The lowest eigenstates in each $`S`$ sector are the states of $`\{^2\stackrel{~}{E}\}`$, describing collinear order. Further support to this assumption is given by the finite size effects of the energy gaps. As shown in Fig.12, a plot of the spin-gap of the $`N=24,32`$ samples versus $`1/N`$ is consistent with a vanishing value for $`N\mathrm{}`$. On the other hand the gap between the two states $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_3`$ of $`\{^2\stackrel{~}{E}\}`$ of the $`S=0`$ sector tends to close when the size goes to infinity, whereas the gaps between the levels of $`\{^2\stackrel{~}{E}\}`$ and the other levels of $`\{^4\stackrel{~}{E}\}`$ increase with $`N`$.

We have investigated the scaling behavior of the spectra at some other points of region IV not too close from the classical boundaries and found essentially the same behavior and a selection of collinear LRO by quantum fluctuations.

Closer to the boundary between region IV and V, the separation between the $`\{^4\stackrel{~}{E}\}`$ states and the magnons states decreases. This is an indication of a softening of the magnons and the neighborhood of a 2nd order phase transition towards another phase. The behavior of the spin-gaps at $`J_1=1,J_2=0.5,J_3=0.5`$ and $`J_1=1,J_2=0,J_3=1`$, similar to Fig.12, nevertheless indicates that these points of the quantum phase diagram are still in the collinear phase IV.

Various studies of the spectra of the $`N=18,26,28`$ samples under suitable boundary conditions confirm these results for the quantum phase IV, and indicate that the quantum boundary between phase IV and V is probably slightly shifted down relatively to the classical boundary shown in Fig.2.

Results of ED calculations (not shown) in region II for $`J_1<0`$ (which has the same classical manifold of degenerate ground-states as in IV for $`J_1>0`$) suggest a similar selection of the collinear solution there too. In conclusion up to a slight motion of the boundaries , the semi classical behavior in regions I, II, IV and VI, is not qualitatively affected by the strong quantum fluctuations of the spins 1/2.

### 3.3 Quantum phases between I and IV

The intermediate phases between the two collinear Néel phases cover region V and part of region III. In this part of the quantum phase diagram, $`SU(2)`$ symmetry is unbroken, and there is a gap to triplet excitations: these phases only support short range order in the spin-spin correlations (see Fig.13).

Our LSW and ED calculations indicate that this quantum region likely extends in regions I and IV <sup>1</sup><sup>1</sup>1LSW calculations predict non vanishing order parameters for the classical spiral solutions inside III away from the boundaries but a vanishing order parameter in the whole region V.. In this work we study region V, region III close to I, and the transition line III-V with ED calculations using TBC on $`N=18,24`$ samples (see Appendix) and PBC on the $`N=24`$ and $`N=32`$ samples. A thorough search of the ED spectra, sweeping the twist angles at the sample boundaries, did not yield evidence of incommensurate helical LRO, neither with the wave-vectors of the classical solutions nor at other wave-vectors. In all cases no tower of QDJS was found. The ED results corroborate the conclusion of LSW calculations that the classical spiral solutions are destabilized by quantum fluctuations. This seems a rather general statement in systems where the quantum fluctuations are strong enoughmoreo90 ; read-sachdev91 .

Is this quantum phase a quantum disordered one? To answer this question we performed extended ED calculations on $`N=24,32`$ samples on different points of the transition line III-V where the classical model has an infinite set of spiral ground-states, and the LSW calculation diverges. Along this line, we found evidence of two different phases both with a gap.

Let us begin by the phase around $`J_2=0.4`$: this point is very close to the point where the energy versus $`J_2`$ is the largest and may be considered as a point of maximum frustration. The spectrum of the $`N=32`$ sample is shown in Fig.14. This spectrum differs from the spectra of the collinear ordered system in IV:

* the lowest states are not IR of $`\{^2\stackrel{~}{E}\}`$

and features associated with Néel LRO are missing:

* The lowest eigen-energies for each $`S`$ value do not increase as $`S(S+1)`$ with $`S`$
* The lowest states in each $`S`$ sector are not separated from the others as the QDJS are separated from the magnons, instead there is a dense continuum of states in each $`S`$ sector except the $`S=0`$ one.
* Furthermore a plot of the spin-gap of the $`N=24,32`$ samples versus $`1/N`$, displayed in Fig.15, shows that the scaling law characteristic of a Néel ordered system is not obeyed and indicates a large spin-gap $`0.2`$ for $`N\mathrm{}`$.

Most likely however the system is not fully disordered but exhibits dimer LRO (see Fig.16). The dimer operator on a pair of sites $`(i,j)`$ is $`d_{i,j}=\left(1P_{i,j}\right)/2`$ where $`P_{i,j}=2(𝐒_i.𝐒_j+1/4)`$ is the spin permutation operator. This projector is greater (less) than $`0.25`$ when the spin-spin correlation is negative (positive), equal to 1 on a singlet and to 0 on a triplet. For $`J_1=1,J_2=0.4`$, on the $`N=32`$ sample, the first neighbor correlation is $`<d_{k,l}>=0.4899`$. On the symmetry breaking Spin-Peierls state (pure product of ordered dimers), the average value of the dimer operator is 1 on the bonds where there is a dimer, and $`1/4`$ on the other bonds. As the exact eigenstate does not break $`C_3`$ symmetry the number $`0.4899`$ should be compared with the result obtained on the symmetric superposition of the three Spin-Peierls states aligned along the three main directions of the lattice: in this symmetrized Spin-Peierls state this correlation is $`d_{k,l}^{\mathrm{\Psi }_{sym}}=0.5`$. The average value of the dimer operator in the exact ground-state is thus very close to the Spin-Peierls value.

The dimer-dimer correlation between a reference bond $`(i,j)`$ and the bond $`(k,l)`$ is $`D_{(i,j),(k,l)}=<d_{i,j}d_{k,l}><d_{i,j}><d_{k,l}>`$. As in ref mblw98 , we normalized $`D_{(i,j),(k,l)}`$ by its maximum value which is achieved when the two bonds are completely correlated. We thus measured dimer correlations by

$`p_{(i,j),(k,l)}`$ $`=`$ $`{\displaystyle \frac{D_{(i,j),(k,l)}}{<d_{k,l}><d_{i,j}><d_{k,l}>}}`$
$`=`$ $`{\displaystyle \frac{<d_{i,j}d_{k,l}><d_{i,j}><d_{k,l}>}{<d_{i,j}>\left(1<d_{k,l}>\right)}}`$
If $`p_{(i,j),(k,l)}=1`$ the presence of a dimer on bond $`(i,j)`$ implies the existence of a dimer on bond $`(k,l)`$; if $`p_{(i,j),(k,l)}=0`$ there is an absence of correlations between singlets on bonds $`(i,j))`$ and $`(k,l)`$. If $`p_{(i,j),(k,l)}`$ is negative, a singlet on bond $`(i,j)`$ induces a tendency towards ferromagnetic correlation on bond $`(k,l)`$.

The correlation pattern for dimer on first neighbor bonds is displayed in Fig.16.

The calcul of dimer-dimer correlations on the state $`\mathrm{\Psi }_{sym}^{S.P}`$ gives $`p_{(i,j),(k,l)}=+0.5`$ if (i,j) and (k,l) are parallel and $`p_{(i,j),(k,l)}=0.25`$ if (i,j) and (k,l) are non parallel bonds. The exact dimer-dimer correlations are not too far from these values and decay very slowly with distance. This is in favor of a columnar LRO of dimers with a $`C_3`$ symmetry breaking, previously proposed by Einarsson and coll. ej91 .

Moreover the degeneracy of the ground-state for $`N\mathrm{}`$ points to the same conclusion: Fig.15 indicates that the gap between the ($`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_3`$, $`S=0`$) lowest states, which are both $`𝐤=\mathrm{𝟎}`$ states, closes for $`N\mathrm{}`$, while the gap between these states and the upper levels increases with the size. $`\mathrm{\Gamma }_1`$ is non degenerate and $`\mathrm{\Gamma }_3`$ twice degenerate: this allows the building of the three columnar dimer patterns with a $`C_3`$ symmetry breaking and no translation breaking. In this picture the finite size ground-state is the symmetric combination of these three states. This degeneracy corresponds to a true symmetry breaking with a local non zero order parameter (dimer LRO): this is a Valence-Bond Crystal(VBC).

The honeycomb lattice could a priori accommodate a different kind of VBC with alternation of hexagons with three dimers and hexagons without dimers: this pattern breaks both $`C_3`$ and translational symmetry. In their large N approach, Read and Sachdev rs90  found that this structure might be the ground-state. In the range $`J_20.30.35`$, we find a short range structure roughly reminiscent of this arrangement. In fact the short range dimer-dimer correlations are even more symmetric than in this VBC crystal and would be more compatible with a crystal of hexagon-plaquettes in a symmetric $`S=0`$ statemo01 : for example the correlation (1-2)(7-6) should be negative and equal to $`0.25`$ in the Read and Sachdev VBC state whereas it is equal to $`0.1`$ in the pure hexagon-plaquette VBC. In fact in our $`SU(2)`$ model all correlations decrease noticeably with distance (Fig. 13, 17) and the pattern does not seem to propagate at large distances. The ground state in this range of parameter is probably a RVB spin liquid. This conclusion is qualitatively substantiated by the study of the energy gaps to the ground-state: plausibly none of them goes to zero at the thermodynamic limit, which would be consistent with the RVB hypothesis. Unfortunately the finite-size effects are rather chaotic: the ground-state energy of the $`N=18`$ and $`N=30`$ samples (samples on which the Read and Sachdev dimer pattern is not frustrated) are larger than the ground-state energy of the $`N=26`$ and $`N=32`$ samples, which frustrate it. This is probably related to the fact that the $`N=18`$ and $`N=30`$ samples do not allow the system to take full advantage of the second neighbor antiferromagnetic coupling, whereas the $`N=26`$ and $`N=32`$ samples do. But the building of singlets on second neighbor bonds tends to destroy VBC patterns and favor a RVB ground-state. All these arguments point in favor of an RVB phase at this coupling parameter: unfortunately the sizes that can be studied do not allow a quantitative determination of the gaps.

The quantum AF $`J_1J_2`$ model on the square lattice close to the point of maximum frustration exhibits dimer singh99  or plaquette capriotti00  LRO; the same kind of conclusion has been drawn for the $`J_1J_3`$ model ll96 , and for the MSE model on the square lattice cgb92 . These phases share qualitative properties with the phase identified for $`J_2=0.4`$: in each cases a collinear LRO is destabilized by frustration giving birth through a 2nd-order phase transition to a massive phase with dimer LRO. Such VBC phases appear in many models on bipartite lattices: Rokhsar and Kivelson rk88  in the Quantum Dimer approach (QD), Dombre and Kotliar dk89  for the Hubbard model, Read and Sachdev read-sachdev91  in the $`SU(N)`$ approach of the Heisenberg model found VBC phases. These phases, as the first one described here (for $`J_2=0.4`$), have a gap for all excitations, a discrete degeneracy of the ground-state, exponential decrease of the two points spin-spin correlations but LRO in higher correlation functions; they have only confined spinons.

Up to now we only know few spin-1/2 models exhibiting true RVB phases with a clear-cut gap: the MSE model on the triangular latticemblw98  and the Quantum Dimer model on the triangular latticems00 . More work is needed to know if the excitations of these different RVB states are similar and in particular if they sustain deconfined spinons excitations.

## 4 $`J_1<0`$: ferromagnetic nearest neighbor interactions

As already underlined above in Sect. 3.2 the classical collinear phases (F or AF) observed for large $`J_2`$ and $`J_3`$ are likely to survive to quantum fluctuations. We thus focus our study on the region of maximum frustration, $`0<|J_2|,|J_3|<0.5`$, corresponding to region V and part of region III of Fig.3. In this situation LSW calculations predict a non vanishing order parameter of the spiral solutions for values of $`J_2`$ and $`J_3`$ not too close to the transition lines. However extensive ED calculations performed on the $`N=18,24`$ samples with twisted boundary conditions do not yield any evidence of spiral LRO.

We thus studied samples up to $`N=32`$ spins to investigate the nature of the ground-state for a few sets of parameters. The most extensive calculations were done at the point $`J_2=.25`$, $`J_3=0`$ on the transition line III-V which may be considered as highly frustrated (for this $`J_2`$ value the ground-state energy is close to its maximum, and in the LSW approach quantum fluctuations destroy LRO). Strong indications that the model has a RVB spin-liquid ground-sate, were found at this point:

* The spectra do not exhibit a tower of QDJS as shown in Fig.18 for $`N=32`$, and $`E_0(S)`$ clearly does not evolve as $`S(S+1)`$ with $`S`$.
* A plot of the spin-gap versus $`1/N`$, shown in Fig.19, indicates that the spin-gap is small but finite when $`N\mathrm{}`$ <sup>2</sup><sup>2</sup>2 In view of Fig.19, one may object to our extrapolation to $`N\mathrm{}`$ on two numerical samples. In fact our conclusion is supported by examination of both the gap and the energy per bond of the samples with sizes 18, 24, 28, 32 with various boundary conditions (available on request at: fouet@lptl.jussieu.fr). This study shows that the variations of these quantities with the size is very small and mainly due to the frustration of the short range antiferromagnetic order between third neighbors due the boundary conditions (see below) and not to the cut-off in the low-energy long wave-length quantum fluctuations. Notice that the energy per bond on the two non frustrating sizes 24 and 32 does not increase with the system size but decreases by a very small amount ($`10^3`$). We thus conclude that we are in the cross-over regime for both sizes 24 and 32 and the extrapolation of the spin gap in Fig.19 is reasonable. .

But contrary to the case with positive $`J_1`$ and $`J_2=0.4`$:

* The correlations display a strong short range order but plausibly no LRO. The short range pattern is original: the first neighbor spin-spin correlation is ferromagnetic ($`<𝐒_i.𝐒_j>=0.10`$), the second (third) neighbor spin-spin correlations are antiferromagnetic ( $`<𝐒_i.𝐒_j>_{n.n.}`$$`=0.13`$, $`<𝐒_i.𝐒_j>_{n.n.n.}`$$`=.25`$ ), but no long range pattern does emerge from this picture. The dimer-dimer correlations equally show a strong short range pattern and apparently no LRO. Fig. 20 represents first-neighbor dimer-dimer correlations: they are much weaker than in the AF case (remark that triplet-triplet correlations are equal to singlet-singlet ones and compare with Fig. 16). Fig. 21 displays second-neighbor dimer-dimer correlations which also decrease with distance. The third-neighbor dimer-dimer correlations decrease even quicker. The strength of the short range correlations explains the finite size results on small frustrating sizes (see footnote <sup>4</sup>).
* The ground-state is probably unique in the thermodynamic limit. The first two singlet excitations are shown together with the first triplet excitation in Fig. 19. The spin-gap and the ground-state energy per spin display a very small sensitivity to the size for $`N=24,32`$. The gap to the third excitation seems more sensitive to the size but this might due to the fact that the different sizes do not accommodate the same wave vectors. In view of the results it seems probable that the singlet excitations will not collapse to the absolute ground-state in the thermodynamic limit.

We thus conjecture that the quantum ground-state of this system does not break any symmetries of the Hamiltonian or of the lattice: it is a quantum spin-liquid  where all excitations are gapped. Results obtained for the $`N=24`$ and $`32`$ samples for $`J_2=.5`$ also indicate a finite spin-gap when $`N\mathrm{}`$. This suggest that there is a spin-liquid phase in a finite range of parameters.

Such a quantum massive phase, without LRO, is highly reminiscent of the spin-liquid phase found in the MSE model on the triangular lattice mblw98 . Curiously enough it appears in the two cases in the vicinity of a ferromagnetic phase destabilized by antiferromagnetic couplings. There is a difference in the degeneracy of the ground-state in the two cases: whereas the ground-state on the honeycomb lattice is unique, it has a 4-fold degeneracy on the triangular lattice. This is easily understood as the honeycomb lattice is not a Bravais lattice and has two spin-1/2 per unit cell. Thus the uniqueness of the ground-state in this latter case does not contradict the Lieb-Schultz-Matthis-Affleck conjecture affleck88 , or the topological approach of Read and Chakrabortyrc89 .

## 5 Conclusions and conjectures

This study of the spin-1/2 $`J_1J_2J_3`$ model on the honeycomb lattice has brought the following new results:

* For small frustrations $`J_2/J_1`$ or $`J_3/J_1`$ less than $`0.15`$ or larger than $`1`$, the system remains essentially classical: when various kinds of LRO are possible, quantum fluctuations, as well as thermal fluctuations in the classical case, select the LRO with the most symmetric order parameter amongst the various possibilities.
* The classical symmetry between the phase diagram for ferromagnetic $`J_1`$ and the phase diagram for antiferromagnetic $`J_1`$ discussed in Sect. 2.1 is destroyed by quantum fluctuations.
* For the largest frustrations these models exhibit gapped phases.
* For an antiferromagnetic first neighbor coupling, a Valence Bond Crystal phase has been clearly evidenced around $`J_2/J_1=0.4`$.
* For an intermediate frustration $`J_2/J_1=0.3`$, an RVB spin-liquid appears between the Néel ordered phase and the VBC phase.
* For a ferromagnetic first neighbor coupling, the present results favor the hypothesis of a RVB spin-liquid phase in a large range of parameters. No VBC has been found in that case.

This study of the spin-1/2 $`J_1J_2J_3`$ model on the honeycomb lattice, when compared to similar approaches of $`SU(2)`$ Hamiltonians leads us to formulate some conjectures on the generic behavior of such models on different lattices.

* In 2D the pure S=1/2 Heisenberg model is Néel ordered on any bipartite lattices with coordination number $`3`$. It is disordered on the triangular-based kagomé lattice which has a coordination number equal to 4.
* Non-coplanar classical ground-states are not robust against quantum fluctuations in the isotropic models.
* Néel order or ferromagnetism disappears around the points of maximum classical frustration giving birth to phases with spin-gap and short range spin-spin correlations.
* Disappearance of a ferromagnetic phase due to antiferromagnetic frustrations leads generically to a spin-liquid phase, with short range correlations in all observables and a gap to all excitations.
* Disappearance of a collinear antiferromagnetic phase might lead to a VBC phase either directly ($`J_1J_2`$ model on the square lattice), or through an intermediate RVB phase (this study). The spin-liquid phase observed by Santoro et al in the spin-orbital modelssgpt99  might be rather similar to the RVB phase described here.

For completion we might add:
* Disappearance of a non-collinear phase (3-sublattice Néel phase) takes place through a phase with a spin gap but a continuum in the singlet sectorlblps97 ; lmsl00 .

Acknowledgements: Computations were performed at The Centre de Calcul pour la Recherche de l’Université Pierre et Marie Curie and at the Institut de Développement des Recherches en Informatique Scientifique of C.N.R.S. under contract 990076 .

## Appendix: Special properties of the studied samples, and boundary conditions

The ED calculations were performed as in refs blp92 ; mblw98  on systems of $`N=18,24,26,28,30,32`$ sites shown in Fig.22. The $`N=18,24,32`$ samples have the full point group symmetry of the lattice, whereas the $`N=26`$ sample misses axial but still has rotational $`C_3`$ symmetry, the $`N=28,30`$ have neither.

With periodic boundary conditions (PBC), all the samples are of course compatible with the $`𝐐=0`$ order in region I. In region IV, however only the $`N=24`$ and $`32`$ samples have the full symmetry of the classical order. The $`N=28`$ sample is compatible with one collinear solution but frustrates the two others as well as the non coplanar solutions. The other samples are frustrating but can allow a collinear order if twisted boundary conditions (TBC) are used. This is the case for the $`N=18`$ sample if a twist of $`\pi `$ is applied along $`𝐭_1`$ or $`𝐭_2`$.

To search for spiral order, we used TBC and sweep the whole range $`[0,2\pi ]`$ of twist angles $`\varphi _{1,2}`$ in the $`𝐭_1`$ and $`𝐭_2`$ directions. These specific boundary conditions are defined as:

$$𝐒(𝐑_i+𝐭_j)=e^{i\varphi _jS_z(𝐑_i)}𝐒(𝐑_i)e^{i\varphi _jS_z(𝐑_i)}.$$
(13)
This allows to look for boundary conditions which would not frustrate helical ground-states. This approach was found effective for the Heisenberg model on the triangular lattice to deal with samples frustrating the three-sublattice Néel order blp92 . For such samples the ground-state energy was found to reach its minimum for the twists which release the frustration: at that point the spectra recover the characteristic features of Néel order.

A VBC with the pattern of Read and Sachdev rs90  (considered in Sect. 3.3) fits in the $`N=30`$ sample but not on a $`N=32`$ sample.