# From Néel long-range order to spin-liquids in the multiple-spin exchange model

## I INTRODUCTION

The two-dimensional triangular lattice antiferromagnet (2D-TLA) was firstly proposed to be a candidate for the disordered (or spin-liquid) ground-state of the spin-$`\frac{1}{2}`$ Heisenberg model (Anderson et al. in 70’s ). Different approaches failed to support this conjecture, but favor a ground-state with Néel long range order (LRO) . Nevertheless, the lattice frustration on the 2D-TLA attracts great interest of theorists, providing a challenge for exotic antiferromagnets. Recently, the multiple-spin exchange model has been extensively studied as an alternative to the Heisenberg model, showing a rich structure of ground-states . In this model, the ground-state can be ferromagnetic (FM), anti-ferromagnetic (AFM) with Néel LRO, or a spin-liquid (SL). A prospective phase diagram has been given by Misguich et al., who considered two-, four-, and five-spin exchange interactions on the 2D-TLA . They found that a large enough four-spin exchange interaction drives th FM phase into a SL phase. They did not study how the AFM Néel LRO is destroyed by the four-spin exchange interaction, and how the transition between Néel LRO and the short range RVB phase takes place. This question is the main object of this paper.

Unhappily there is no exact method allowing the study of the zero temperature phases of such frustrated systems. For a finite system, however, one can always, in principle, represent the eigenstates in the complete basis of spin configurations. This allows one to have a real touch on the exact ground-states of small size systems through numerical computations. The huge number of spin configurations ($`2^N`$) becomes a great obstacle on the way of numerical simulations. On the most recent computers, the largest sample that may be handled in exact diagonalizations has $`6\times 6`$ sites. On the triangular lattice, the Quantum Monte Carlo method is plagued by the well-known sign problem, but a new technique called Stochastic Reconfiguration allows handling samples up to $`12\times 12`$ sites . All these calculations point to Néel LRO, with a sublattice magnetization of the order of $`40\%`$ of the saturated value . Series expansions give a reduced ($`20\%`$) but non-zero sublattice magnetization .

On the other hand, in the case of short range correlations, the situation is more straightforward, as soon as the available sizes are of the order of, or larger than the correlation length. This is fortunately the case in the SL phase found in the $`J_2J_4`$ model ($`J_20,J_4>0`$) by Misguich et al. .

In this work, we use exact diagonalizations to obtain the exact eigenenergies versus wave vectors and total spin for 12, 16, 19, 21, 24 and 27 site samples of the $`J_2J_4`$ model ($`J_2=1,J_4>0`$) on the 2D-TLA. In the classical limit, the AFM ground-state of the 2D-TLA can be described as a three-sublattice structure, with spins of different sublattices making angles of $`2\pi /3`$. Periodic boundary conditions are compatible with the three-sublattice structure for samples with 12, 21, 24 and 27 sites, but not for the 16 and 19 site samples. Therefore, we use twisted boundary conditions for the 16 and 19 site samples  and periodic boundary conditions for the 12, 21, 24 and 27 site ones.

## II THE MODEL: THE MULTIPLE-SPIN EXCHANGE HAMILTONIAN

The Hamiltonian of the multiple-spin exchange model is given by

$$H=\underset{n}{}(1)^nJ_n(P_n+P_n^1),J_n>0,n2$$
(1)   
where $`J_n`$ are the $`n`$-spin exchange tunneling probabilities (exchange coefficients), $`P_n`$ and $`P_n^1`$ are the $`n`$-spin exchange operators and their inverse operators, respectively. The alternative sign in the summation over $`n`$ in Eq. 1 comes from the permutation of fermions. In general, the exchange coefficients decrease with increasing $`n`$. The two-spin exchange term gives exactly the Heisenberg Hamiltonian up to a constant, since one has

$$P_2=2𝐒_i𝐒_j+\frac{1}{2}$$
(2)   
where $`𝐬_i`$ and $`𝐬_j`$ are spins localized at site $`i`$ and $`j`$, respectively. The three-spin exchange operator is exactly equivalent to a sum of two-spin exchange operators . Thus, the three-spin exchange term of equation (1) can be absorbed into the two-spin exchange term, as long as $`J_2`$ is replaced by the effective two-spin exchange coefficient $`J_2^{\mathrm{eff}}=J_22J_3`$. Therefore, except for the two-spin exchange, the next most important term is the four-spin exchange. A pure positive two-spin exchange (i.e., the Heisenberg Hamiltonian) on the 2D-TLA gives an AFM phase with Néel LRO, and, as shown in Ref. , a pure four-spin exchange gives a SL phase. In this paper we use the specific properties of the spectra of these different kinds of phases to study the transition from one phase to the other, when the relative weight of the four-Spin Exchange $`J_4`$ increases relatively to the antiferromagnetic two-spin coupling. In the following, all the energies are measured in units of $`J_2^{\mathrm{eff}}=1`$.

Little is known on this region of the phase diagram. Previous works are based on a classical approximation , semi-classical spin-wave calculations  or mean-field Schwinger-boson results . The classical result predicts a transition from the 3-sublattice Néel state to a 4-sublattice tetrahedral state at $`J_4=0.24`$. Both quantum approaches indicate that the four-spin exchange strongly enhance fluctuations in the 3-sublattice Néel phase. Kubo et. al  found that the sublattice magnetization vanishes for $`J_4>0.17`$. In the Schwinger-boson approach , the Néel state is destroyed when $`J_4>0.25`$. These two techniques have a general tendency to underestimate the effects of quantum fluctuations on ordered phases. The exact diagonalization analysis presented here indeed shows that the Néel long-ranged order disappears for a smaller value of $`J_4`$ (the critical value is estimated to be in the interval $`J_4^C0.07\mathrm{}0.1`$).

## III CRITERION TO DISCRIMINATE BETWEEN NÉEL LRO AND SL PHASE: THE SPIN GAP ?

### A Finite-size energy spectrum of the Néel LRO phase

In the classical limit, a $`N`$-site 2D-TLA sample with Néel LRO is characterized by a three-sublattice structure with spin $`N/6`$ on each sublattice. Coupling of these three $`N/6`$-spins, gives total spin $`S`$ with $`min\{2S+1,N/2S+1\}`$ degeneracy . In an isotropic antiferromagnet (as the collinear AFM which has equal spin susceptibilities and spin wave velocities) the finite-size total energy depends on the total spin $`S`$ (to first order in $`1/N`$) as:

$$E_S=E_0+\frac{1}{2N\chi }S(S+1),$$
(3)   
where $`E_0=Nϵ_0`$ is the energy of the ground-state in the thermodynamic limit and $`\chi `$ is the isotropic magnetic susceptibility of the sample. In the anisotropic case this equation should be rewritten:

$$E_S=E_0+\frac{1}{2N\chi _{}}S(S+1)+\frac{1}{2N}(\frac{1}{\chi _{}}\frac{1}{\chi _{}})S_3^2,$$
(4)   
where $`S_3`$ is the component of the total spin $`S`$ on the internal symmetry axis of the spin system, and $`\chi _{}`$ and $`\chi _{}`$ are the magnetic susceptibilities on the internal symmetry axis and on the perpendicular plane, respectively. In the broken symmetry picture the symmetry axis is perpendicular to the plane of the spins and $`\chi _{}`$ (respectively $`\chi _{}`$) measures the spin fluctuations orthogonal to (respectively in) the spin plane. Eq. 3 (respectively Eq. 4), is the dynamical equation of a rigid rotator (respectively of a quantum top).

Eqs. (3, 4) show that the slopes of the total energy versus $`S(S+1)`$ and $`S_3^2`$ approach zero as $`1/N`$ does when $`N\mathrm{}`$. $`S_3`$ is an internal quantum number dynamically generated, which is not under control in a finite-size study. But the total spin $`S`$ is a good quantum number and the $`N^1`$ scaling of the $`S(S+1)`$ dependence of the total energy versus sample size is interesting because it is more rapid than the scaling law of the order parameter (which goes as $`N^{1/2}`$.

### B Finite-size scaling in the SL phase

In a SL phase, contrarily to the Néel LRO phase, the spin gap (i.e. the difference in total energy between ground-states in the $`S=1`$ sector and in the $`S=0`$ sector) does not collapse to zero in the thermodynamic limit. The finite-size scaling law in this second situation is not known exactly, insofar as the ”massive” phase is not characterized precisely. Heuristically, we expect the finite-size spin gap to decrease exponentially to a finite value $`\mathrm{\Delta }(\mathrm{})`$ with the characteristic length $`\xi `$ of the spin-spin correlations. For samples of linear size $`L`$ smaller than the correlation length and in the cross-over regime, there are not enough quantum fluctuations to destroy the sublattice magnetization and the system probably behaves as it were classical (i.e. with a spin gap decreasing as $`N^1`$). The following heuristic law might be used to interpolate between the two behaviors:

$$\mathrm{\Delta }(L)=\mathrm{\Delta }(\mathrm{})+\frac{\beta }{L^2}\times exp(L/\xi )$$
(5)   
### C Quantum critical regime

The use of this heuristic law (Eq. 5) encounters a severe difficulty as soon as the disordered system approaches a quantum critical point: in such a situation the correlation length $`\xi `$ diverges, the gap closes to zero and on a finite-size sample it is impossible to discriminate between such a situation and isotropic Néel LRO (Eq. 3).

### D Numerical results

In view of this difficulty we have done a pedestrian finite-size scaling of the spin gap by using the simplest linear $`1/N`$ behavior which probably gives a lower bound of the gap in a SL outside critical points. The physical reason is that we expect the finite-size corrections of the gap value to be smaller in a system with a finite correlation length (no long-range order) than in a LRO Néel phase where it vanishes as $`1/N`$<sup>*</sup><sup>*</sup>*This assumption would be invalid at a critical point where the uniform susceptibility vanishes: In such a case the gap might close as $`1/\sqrt{N}`$. The analysis of section IV shows that the behavior of the system changes rather abruptly from a Néel like spectrum to a spectrum with a very large number of low lying singlets below the first $`S=1`$ state (and potentially a $`T=0`$ residual entropy in the singlet sector). This happens before any decrease to zero of the spin velocity or of the homogeneous spin susceptibility. So, if the true thermodynamical system has indeed a critical point between the two phases, the sizes we are looking at are too small to scrutinize the critical regime. The results extrapolated to $`N\mathrm{}`$ are shown in Fig. 1. Strictly speaking the spin gap never extrapolates to zero except for a pure $`J_2`$ where it is equal to zero within its error bar (At $`J_4=0`$, data of Bernu et al.  for $`N=36`$ are added to the present results, see Fig. 1.). Nevertheless these data already show three distinct ranges for the parameter $`J_4`$: for very small $`J_4`$ (below 0.075) Néel LRO is plausible but should be confirmed by another approach. For $`J_4`$ larger than $`0.1`$ a gap certainly opens rapidly with increasing $`J_4`$ and then decrease for $`J_4>0.175`$. The spin gap criterion cannot give more insight on the phase diagram. We will now move to the analysis of the symmetries of the low lying levels of the spectra to characterize more precisely these three phases.

## IV SYMMETRIES OF THE LOW LYING LEVELS IN A NÉEL ORDERED PHASE

### A Theoretical background

Firstly we show the low energy spectrum of the pure Heisenberg model on the 21 site sample (Fig. 2). In order to emphasize the low energy structure we have displayed the low energy spectrum minus a rigid rotator energy $`\alpha S(S+1)`$.

Let us first concentrate on the lowest part of the energy spectrum in each S sector (solid and open triangular symbols in the figure). This family of levels forms on a finite-size lattice the quantum counterpart of the semi-classical Néel state. These specific states are in the trivial representation of the invariance group of the three-sublattice Néel ordered solution . To be definite:

* The Néel ground-state breaks the 1-step translation but is invariant in a 3-step translation: as a consequence the only wave vectors appearing in this family of QDJS (for Quasi-degenerate-joint-states defined by Bernu et al. ) are respectively the center $`𝐤=(0,0)`$ and corners $`\pm 𝐤_\mathrm{𝟎}`$ of the Brillouin zone.
* These QDJS belong specifically to the trivial representation of $`C_{3v}`$ as the Néel state itself (i.e. they are invariant in a $`2\pi /3`$ rotation, and in a reflection symmetry).
* As the $`\pi `$ rotation symmetry of the lattice is broken in this particular ground-state, these QDJS appear either in the odd or even representation of the 2-fold rotation group (see Ref.  for more details).
* The numbers and characteristics (quantum numbers) of the QDJS in each $`S`$ sector are precisely fixed by theory : for the 21 sites spectrum displayed in Fig. 2 (as for all sizes that have been studied up to now) the numbers of low lying levels and their quantum numbers correspond exactly to the above-mentioned theoretical predictions.

### B Numerical results

The dynamical law given by Eq. 4 is still imperfectly obeyed for the 21 site sample: in particular the generation of the internal symmetry is still imperfect but nevertheless the spectrum of a quantum top could already be anticipated.

Above these levels with specific properties, there appear eigenstates with wave vectors belonging to the inside of the Brillouin zone (simple dashes in Fig. 2). A group of such eigenstates with different total spin represents a magnon excitation of the Néel ground-state. As the antiferromagnetic magnons have a linear dispersion law, the softest magnon energy scales as the smallest wave vector accommodated in the Brillouin zone of the finite-size sample. Thus, for sizes large enough, these levels collapse to the ground-state as $`1/\sqrt{N}`$, more slowly than the QDJS which collapse as $`1/N`$ to the thermodynamic Néel ground-state energy. That is the reason of the appearance in Fig. 2 of a quasi gap between the QDJS and the magnon excitations.

This hierarchy of low lying levels is a very strong constraint on the finite-size samples spectra. It is perfect for sizes up to 27 and for $`J_4`$ smaller or equal to 0.075 and totally absent An illustration of the striking qualitative differences between spectra in the Néel region and in the spin-liquid phase has been previously published: two $`N=27`$ spectra at $`J_4=0`$ and $`J_4=0.1`$ are displayed in Fig. 2 and Fig. 3 of Ref. . for $`J_4`$ larger or equal to 0.1.

* This result, associated to the spin gap behavior, consistently proves that for $`J_4`$ larger or equal to 0.1 the system is in a spin-liquid state with rather short range spin-spin correlations.
* For $`J_40.075`$ the structure of the low lying eigen-levels of the spectra are compatible with Néel LRO. BUT as discussed above, it is indeed impossible to precisely point a quantum critical transition within this approach. In view of the spectra, we might speculate that the transition is second order and that it takes place between 0.07 and 0.1.
* From $`J_4=0`$ to $`J_40.1`$ we see a softening of the spin-wave velocity consistent with the gradual decrease of the Néel LRO (Fig. 3). However, Fig. 4 shows that a large number of singlet states are already present at low energy when $`J_4=0.1`$. Therefore, the system is certainly no longer in the ordered phase at $`J_40.1`$. We can conclude from Fig. 3 that the spin wave velocity does not vanish at the critical point (even if the precise location of the transition cannot be determined). This has been previously suggested by various analytical approaches . The sizes studied are nevertheless too small to check Azaria’s prediction  of an $`O(4)`$ symmetry of the effective field theory at the critical point.

## V THE LOW ENERGY EXCITATIONS OF THE SL PHASE

Contrarily to our expectations, the SL phase which appears immediately after the disappearance of the Néel ordered phase is not the phase studied by Misguich et al. . It is indeed a SL phase, with a gap and short range spin-spin correlations. But, as might be seen in Fig. 4, this phase exhibits a very large number of singlets in the magnetic gap and seems in this respect similar to the spin-liquid phase of the Heisenberg model on the kagomé lattice . However, we are not aware of any exponential degeneracy (i.e. $`\mathrm{exp}(N)`$) in the classical MSE model Momoi et al.  found a ground-state degeneracy in the classical MSE model but this degeneracy only grows as the exponential of linear size of the system ($`\mathrm{exp}\sqrt{N}`$). Moreover, this classical degeneracy was only found in a region where $`J_2<0`$ is ferromagnetic ($`\frac{1}{4}J_4/|J_2|\frac{3}{4}`$)., as is the case for the classical kagomé antiferromagnet.

We suspect that a much larger four-spin exchange parameter will be needed to recover the SL phase studied by Misguich et al.. These data point to the existence of this new phase in a finite range of parameters $`0.075J_40.25`$. However, one cannot disregard the hypothesis that these properties are in fact those of a critical point, with a critical region enlarged by finite-size effects. More work with different methods is needed to clarify this point.

## VI Magnetization Plateaus

In an external magnetic field $`B`$ along the $`z`$ axis, the total energy of the state with component $`S_z`$ of the total spin is given by:

$$E_B=E_SS_zB$$
(6)   
The magnetization is determined by the minimum of $`E_B`$ respective to $`S_z`$, which requires $`E_B/S_z=0`$. Therefore, in the isotropic case with Néel LRO, one has

$$m=2\chi B1/N$$
(7)   
where $`m=2S/N`$, is the polarization relative to the saturated magnetization $`N/2`$ and $`\chi `$ is indeed the magnetic susceptibility.

We noticed that when the four-spin exchange interaction increases, deviation from Eq. 3 occurs about $`S_z=N/6`$ and $`N/4`$. This is in agreement with the earlier mean-field calculation of Kubo and Momoi  who predicted magnetization plateaus at $`m=1/3`$ and $`m=1/2`$ in the $`J_2J_4`$ model. We present the magnetization curve of the 24-site sample in Fig. 5. A small plateau exists at $`1/3`$ magnetization for $`J_4=0`$, and its width first increases slowly as $`J_4`$ increases and then decreases from around $`J_4=0.125`$. This $`1/3`$ plateau has also been found in previous studies (see  and references therein) in the pure three-sublattice Néel ordered system. The $`m=1/2`$ magnetization appears at about $`J_4=0.1`$ and its width increases with $`J_4`$. Finally, there still exists a plateau at about $`m=1/2`$ in samples with odd numbers of sites, but it distributes over the two closest positions to the $`1/2`$ magnetization. The $`m=1/3`$ and $`m=1/2`$ plateaus correspond to the classical $`uud`$ and $`uuud`$ ordering of spins (see Momoi et al.  for the four-sublattice $`uuud`$ state). These phases are more ”classical” than the zero-field phase, but yet show a decrease of the sublattice magnetization from the classical saturation values.

## VII Conclusion

In conclusion, we studied the transition between the Néel ordered state and a Spin Liquid state of the multiple-spin exchange model by means of the exact diagonalization method. The pure three-sublattice Néel ordered phase is gradually destroyed by quantum fluctuations when increasing the 4-spin exchange coupling. The spin wave velocity decreases but apparently remains finite at the transition. The quantum disordered phase tuned by the 4-spin exchange coupling is different from the pure $`J_4`$ phase studied by Misguich et al.. It exhibits low energy singlet excitations, reminding of kagomé SL. This result opens many interesting questions that cannot be answered in the present framework: is this phase a new generic SL phase or a finite-size manifestation of a quantum critical regime ? Are these singlet excitations the “resonon” modes invoked by Rokhsar and Kivelson  ? In agreement with previous studies , we find two magnetization plateaus at $`1/3`$ and $`1/2`$ of the full magnetization. These plateaus are associated to the semi-classical $`uud`$ and $`uuud`$ ordering structures. A finite-size scaling on much larger sizes is needed to draw a definite conclusion on this magnetic phase diagram.