# Phase diagram of the spin-1 Heisenberg model with three-site interactions on the square lattice

###### Abstract

We study the spin antiferromagnetic Heisenberg model on the square lattice with, in addition to the nearest-neighbor interaction, a three-site interaction of the form . This interaction appears naturally in a strong coupling exansion of the two-orbital, half-filled Hubbard model. For spin 1/2, this model reduces to a Heisenberg model with bilinear interactions up to third neighbors, with a second-neighbor interaction twice as large the third-neighbor one, a very frustrated model with an infinite family of helical classical ground states in a large parameter range. Using a variety of analytical and numerical methods, we show that the spin-1 case is also very frustrated, and that its phase diagram is even richer, with possibly the succession of seven different phases as a function of the ratio of the three-site interaction to the bilinear one. The phases are either purely magnetic phases with collinear order, or of mixed magnetic and quadrupolar character with helical order.

## I Introduction

The quest for exotic phases in magnetic quantum systems has driven a lot of research over the past years. The appearance of new phases is often related to the inclusion of terms beyond the first neighbor Heisenberg interaction. For spin , it can be for example longer range interactions leading to frustrationbook_FM (), Dzyaloshinskii-Moriya interactionsDzyaloshinsky (); Moriya () or plaquette interactionsColdea (). These terms also appear in spin systems. However, the fact that the local Hilbert space is of higher dimension for allows for new types of interaction. Starting from a two-orbital Hubbard model at half-filling with Hund’s coupling, a strong coupling expansion leads to an effective spin-1 model. At second order, only a Heisenberg interaction between first neighbor appears. At fourth order, on top of the usual terms that also appear in the case starting from the single band Hubard modelTakahashi (); macdonald (), namely second-neighbor and plaquette interactions, two new interactions appear, the biquadratic interaction, and a three-spin interaction PhysRevB.76.132412 (); MVMM () of the form . The biquadratic interaction has been intensively studied over the past few years, both in one dimensionlai (); Sutherland (); Takhtajan2 (); Babujian2 (); AKLT (); Chubukov (); Fath1 (); Fath2 (); Jolicoeur (); Laeuchli (); Salvatore () and two dimensionsKawashima (); BLBQtsunetsugu (); BLBQmila (); tamas (). However to the best of our knowledges, the effect of the three-spin interaction has not been investigated in 2D, which has led us to study the following model on a square lattice:

(1) | |||||

where are spin operators, sums over first neighbors and sums over all possible configurations where and are first neighbors of and are different from each other.

It is worthwhile noticing that in the case of spin-1/2, the three-spin interaction reduces to a second-neighbor interaction. On a chain, this leads to the model which has a transition to a dimerized ground state at J1J2transition (). In Ref. MVMM, , it has been shown that the three-spin interaction generalizes this result to higher spins and that, for any , the ground state shows dimerization in some region of the phase diagram. In two dimensions and for spin-1/2, the Hamiltonian of eq. (1) reduces to a Hamiltonian with Heisenberg interactions , and to first, second and third neighbor with . A discussion of this model can be found in Ref. [J1J2J3rastelli, ; J1J2J3chandra, ; J1J2J3, ; J1J2J3ferrer, ; J1J2J3gochev, ; J1J2J3yang, ; J1J2J3m, ; J1J2J3Arlego, ]. The line is of particular interest because it lies at the classical transition between a helical phase with pitch vector and another helical phase with pitch vector . Starting from an antiferromagnetic phase with Néel order, the system undergoes at a phase transition to a phase with an infinite number of helical ground states whose pitch vectors are defined by the condition:

(2) |

At the level of linear spin-wave theory, quantum fluctuations have been shown to completely disorder this phase. This is the first hint that frustration might be very high in the case of the Hamiltonian of Eq. (1) as well. As we shall see, the situation is even more complicated, and the competition between several phases leads to a very rich phase diagram.

The paper is organized as follows. In section II, we determine the different phases in the classical limit and in particular we discuss the degeneracy of the phase appearing in the limit . In section III we discuss the mean-field phase diagram based on a product of local wave-functions and the difference between this phase diagram and the classical phase diagram. We add quantum fluctuations to the system in the context of a semiclassical expension around the classical solutions in section IV. To confirm the previous results, we give some exact diagonalization results in section V. We finish the paper by a conclusion in section VI

## Ii Classical phase diagram

In this section, we consider the spins as three dimensional arrows of length . The zero temperature phase diagram is obtained by minimizing the energy of Hamiltonian (1) under this assumption. The classical phase diagram can be formally established by decomposing the Hamiltonian into a sum of local Hamiltonians in the following way:

where is the number of sites, and and sum over the nearest neighbors. Minimizing independently all these local Hamiltonians is a sufficient, though not necessary, condition to have a global minimum. We will see that we were able to independently minimize the local Hamiltonian for all values of the ratio .

For symmetry reasons, the central spin of the local Hamiltonian can be chosen to point in the direction. Since the energy only depends on the relative angle between the central spin and its four neighbor (and not on the angle between the neighbors themselves, as would be the case with for example a second neighbor interaction), all the spins can be assumed to point in the same plane, say the plane. Under these assumptions, the classical energy can be written:

where are the angles of the four neighbors of the central spin with respect to . This energy can be easily minimized, and we find three different types of minima, depending on the value of . These three different local minima can be extended to the entire lattice and constitute the three phases that appear in the classical phase diagram. We will now describe in detail the three different phases.

### ii.1 Néel phase

In the limit where goes to zero, the classical ground state is of Néel type. The energy per site is given by .

### ii.2 Up-up-down-down phase

The first transition that appears when increasing leads to the phase depicted in Fig. 1. This phase realizes a compromise between the Heisenberg interaction and the three-body interaction. The three-body interaction favors phases where we have a local up-up-down or down-down-up structure. In the intermediate phase, this is realized by having up-up-down-down chains which are coupled antiferromagnetically so that the Heisenberg interaction is still partially satisfied. The energy per site is then given by . The transition point can be established by comparing the classical energy of the two phases. We see that the transition to this phase appears already at . In the following, we will refer to this phase as the up-up-down-down (uudd) phase.

### ii.3 phases

If , the minimization of the local Hamiltonian is trivial. Because of the rotational invariance of the Hamiltonian, the central spin can point in any direction. If we choose it to point up, then the minimum is realized by having two neighbors pointing up, and two neighbors pointing down. This still holds if the central spin points down. A global minimum can be found if we manage to build a configuration where every spin has two neighbors pointing up and two neighbors pointing down. An easy way to build phases which respect this constraint is to draw lines of spins up and lines of spins down so that two lines of spins up (or spins down) never touch each other. Several examples are given in Fig. 2. In the first one (1) we draw straight lines. This can be deformed by adding domain walls where all lines form an angle, as in (2). In (3), we put all possible domain walls. A second possibility shown in (4) is to draw squares. This phase can also be transformed by adding bigger squares around the square (5). Finally, we can also mix phases with closed loops and phases with lines as seen in (6).

The local constraint seems weak. However, the residual entropy scales as and not as (where L is the linear size of the system). To prove it, let us construct a classical ground state as described in Fig. 3. We start from an arbitrary cell of four spins (in purple in Fig. 3). Then, by turning around this cell following a spiral, we put new spins who respect the constraint imposed by the spins in the inner side of the spiral. The spin along the sides are completely constrained by the inner spins (in black on Fig. 3). At most two spins when the path form an angle might^{1}^{1}1We say "might" and not "can" because two neighbors are enough to constrain one spin if they are pointing in the same direction be choosen freely (in white in Fig 3). The number of free spins will therefore scales at best as and the number of possible configuration as . This is only an upper bound for the degeneracy. Now, one can
build at least ground state in a system of linear size by adding from to domain walls to the phase with straight lines of Fig. 2 (1) which give us a lower bound. Therefore, the residual entropy in the thermodynamic limit lies between and and scales as the linear size of the system .

For classical spins, the Heisenberg interaction does not lift the degeneracy between these different phases, since it is zero for all of them (by definition of the constraint). In the classical case, we therefore expect all these phases to coexist, even for finite . The energy per site of this phase is . The transition to the intermediate phase takes place at .

Overall, the classical phase diagram contains three regions, the last one showing a big degeneracy. It is sketched in Fig. 4.

## Iii Mean-field phase diagram

While classical phases give in general some first intuition about what the quantum phase diagram could be, they still show some down sides. One of them is the fact that if one considers quantum spin systems, the local order parameter can be of quadrupolar type instead of magnetic type. This kind of local order is known to be of crucial importance when biquadratic interactions enter the gameintro_spin_nematic (). Since the three-body interaction involves the square of some local operator, it is legitimate to wonder if, in this case as well, the quadrupolar component plays a role. To tackle this difficulty, it is useful to do a mean-field phase diagram instead of a classical phase diagramintro_spin_nematic ().

To be more precise, the mean-field approach consists in reducing the Hilbert space to states which are products of local wave functions, i.e. to consider only states of the form:

The local wave function for a spin is in general described by three complex numbers. However, the condition that the norm is one, and the freedom to fix the phase, reduce the freedom to only four real parameters. The local wave function can be chosen as:

where , , and are real numbers. It is easy to check that the norm of the local wave function is one and therefore . The mean-field approach consists then in minimizing the mean-field energy , which is a function of variables.

One of the main reasons to use this approach is to detect quadrupolar phases, as shown in the case of biquadratic interactionsBLBQmila (); tamas (). However, even for purely magnetic states, the phase diagram can be different from the classical one if one considers terms beyond the Heisenberg interaction. This effect can be traced back to the fact that the mean-value of a product of operators is the product of the mean-values of these operators only if these operators commute with each other. Since the spin operators on different sites commute, for a Hamiltonian which is linear in spin operator on all sites, the classical and the mean field values are the same. For example, for a Heisenberg Hamiltonian:

By contrast, if one considers a Hamiltonian with higher order on-site terms such as , which appear in the biquadratic and in the three body interactions, the decomposition as a product of mean field terms cannot be performed and the classical and mean-field energies are different, even for coherent states. As an example, let us compare the mean field and the classical energy of the biquadratic interaction in the state:

The two energies are indeed different. In the model presented here, this difference between the classical and the mean-field energy for a given configuration plays a more important role than the possibility to have quadrupolar order, and it lies at the root of the appearance of helical phases which are not present in the classical case.

The strategy to find the phase diagram is the following. First, we minimize the mean-field energy of finite clusters, up to , leaving all variables free. Based on the results, we do an educated guess to reduce the number of free variables to a number which does not scale with the system size (for example the angle between nearest-neighbor sites). This allows us to find the energy in the thermodynamic limit with only a few angles left free. We then have to check that the trial wave function gives an energy per site which is smaller or equal to the energy of finite systems. However, we can never be sure that a bigger cluster would not lead to a better energy, and therefore that we are not missing some phases.

Obviously, the classical phases are included in the Hilbert space of the mean-field phase. They can therefore appear in the phase diagram. Indeed, the three classical phases are realized in some region of the phase diagram. However, new phases appear between the classical phases. Based on the results obtained from the numerical simulations, we assume that the system has either a unit cell of or is an helical phase obtained by the rotation of a cell. We also assume that the spins lie all in the plane. Finally, we assume that the length of the quadrupole is the same on every site. This provides us with an energy with only four independent variables, one for the relative angle inside the unit cell, two for the rotation of the unit cell in the and directions, and a last one for the length of the spin. Using this assumption, we can find an expression for the energy in the thermodynamic limit. However, the exact form of the energy is still too complicated to be minimized by an analytical treatment. Nevertheless, it can be minimized numerically with a very high precision. In the case of a second order phase transition, we can also predict the critical value where the phase changes to a helical phase by taking the Taylor expansion of the energy around the different classical phases. The transition point is given by the condition that the quadratic term vanishes.

### iii.1 First helical phase

Starting from the Néel state, we have a first instability at . This state is formed by exactly antiferromagnetic chains along one direction, while along the other direction, the angle between the spins is given by . Close to the transition, . The quadrupolar order is not driving the transition. However, once the helical phase appears, the spins develop a small quadrupolar component. The phase ends with a first order transition to the up-up-down-down phase, around 0.087. This phase is depicted in Fig. 5 (a).

### iii.2 Second helical phase

The mechanism leading to the second helical phase is similar. We start from an up-up-down-down phase, and we start to tilt the spins in the direction where they are coupled antiferromagnetically. The transition point is at . This is also a second order phase transition, and the deviation with respect to the classical phase will also behave as . As in the previous case, a small quadrupolar component appears. The phase ends with a first order transition to a third helical phase at . This phase is depicted in Fig. 5 (b).

### iii.3 Third helical phase

### iii.4 Phase diagram

Overall, the three classical phases show up at some point of the mean-field phase diagram, and on top of that, three helical phases appear between the classical phases. In Fig. 6 we show the difference between the mean-field energy and the classical energy. The mean-field phase diagram can be seen in Fig. 7.

## Iv Quantum fluctuations

While the first two phases of the classical phase diagram are expected to show up in the quantum case, with of course some renormalization of the boundary, the third phase is more subtle. We expect some lifting of the degeneracy by quantum fluctuations, i.e. by a process of order by disorder. To see which phase is selected by the quantum fluctuations, we use linear spin wave theory.

The Hamiltonian for any classical phase can be decomposed as a sum over the magnetic unit cell of local Hamiltonians:

where runs over the different magnetic cells, runs over the elements of one unit cell, and and run over the neighbors of site .

We can then perform a Holstein-Primakov transformationPhysRev.58.1098 () and expand it in power of . For a collinear phase, this leads to:

where is the length of the spin, is the pitch vector of the phase, is the position of the spin , and and represent Holstein-Primakov bosons.

Under this transformation, keeping only the term of order (classical energy) and the term of order (two-boson interactions), the Hamiltonian reads:

It is then well known that, because of the periodicity of the lattice, the Fourier transform of this expression will decouple the Hamiltonian. Taking the following convention:

where is the number of sites and is the position of site , the Hamiltonian reads:

The different phases have different pitch vectors and different magnetic cells, leading to different bosonic Hamiltonians.

### iv.1 Correction to the energy

Since for we have an infinite number of classical ground states, we cannot perform the spin wave calculation for all of them. In particular, we can only treat phases which are periodic. We restrain ourselves to four phases which look paradigmatic to us. The four phases with the associated unit cells are depicted in Fig 8. The phase formed by straight lines of up spins alternating with lines of down spins is the ground state of the model for large and we start our analysis by this phase. We can change this phase by adding domain walls on every other diagonal (which we refer to as the three-three phase), or on every diagonals (which we refer to as the step phase). Finally, the plaquette phase is formed by square of 4 spins pointing all in the same direction alternating with square of spins pointing all in the opposite direction.

The correction to the energy due to quantum fluctuations is depicted in Fig. 9. We see that quantum fluctuations indeed lift the degeneracy between most of the phases already at the level. After the transition from the up-up-down-down phase, quantum fluctuations select the plaquette phase. Around , the system undergoes another phase transition due to quantum fluctuations presumably to the phase. The phase selected after the transition is however not so clear from linear spin wave theory. Indeed, the phase is degenerate with the three-three phase. We can however speculate that, at higher order, the interaction between the domain walls lift the degeneracy, probably in favor of the phase without any domain wall, considering the fact that the phase with a high density of domain walls is not good energetically.

### iv.2 Reduction of local moment

It is also useful to consider the average number of bosons to have an estimate of the amplitude of quantum fluctuations. We have computed the fluctuations for all the classical ground states which are present in the spin-wave phase diagram. The helical phases which appear in the mean-field phase diagram are not minimum of the classical energy, and a conventionnal spin-wave theory cannot be performed around them. The calculation of the quantum fluctuations in these cases would require more sophisticated tools which are beyond the scope of this paper.

For a spin-1 system, the projection of the spin along the classical order axis is given by . If the approximation becomes unphysical, and it is well controlled in the limit . As seen in Fig. 10, the fluctuations remain quite small for all four phases, a good indication that the classical ground state is a good approximation of the true ground state in these regions. At for the phase, and at for the plaquette phase, the quantum fluctuations diverge (not shown), but within the region where the classical ground states are also minima of the mean-field phase diagram, they remain finite and never exceed .

## V Exact diagonalization

The exact diagonalization approach is very limited for spin system because of the size of the Hilbert space which grows as . The only cluster that is compatible with all the classical phases discussed in section IV and which is accessible with this method is a cluster. This cluster has some additional symmetries depicted in Fig. 12 (a) which make it difficult to treat. As shown in Fig.12 (b), if one picks a reference point, there are only four different types of sites which are not equivalent by symmetry arguments.

Despite this fact, we can still check that the obtained results are compatible with exact diagonalization. To do so, we compute the spin-spin correlation function with respect to a reference site. Since we do not expect the quantum ground state to break the symmetry of the Hamiltonian, we have to compare the quantum quantity to a superposition of classical configurations. More precisely, we have to average the quantities over all the classical configurations which are equivalent up to a symmetry of the cluster. In Fig. 13 we plot the correlation function for different values of . We clearly see at least three different phases. In particular, the second diagonal spin (2,2)^{2}^{2}2By this notation, we mean the site at a distance . starts with a ferromagnetic correlation, it then becomes antiferromagnetic and is again ferromagnetic for large .

The value for obviously corresponds to a Néel state. At , the fifth neighbor (2,2) is strongly antiferromagnetic. This is in agreement with the up-up-down-down phase. Moreover, the second neighbors (1,1) and (2,0) are very close to zero, which is also in agreement with the up-up-down-down phase where half the second neighbors are ferromagnetic and the other half are antiferromagnetic. Finally, the first neighbor is antiferromagnetic, reflecting the fact that in the classical phase, three of them are antiferromagnetic, and one is ferromagnetic. We are therefore confident that the up-up-down-down phase is stabilized in the thermodynamic limit.

The exact diagonalization for does not provide a lot of information. In a cluster, the phase and the plaquette phase are equivalent by the symmetry arguments presented in Fig. 12. It is therefore hopeless to try to distinguished between these two phases with this cluster. Moreover, the average over the first, second, third and fourth neighbor is the same for all the phases of Fig. 8^{3}^{3}3This is true only for this cluster. In the thermodynamic limit, for example the second neighbor is different from one phase to the other.. We therefore have to think about another quantity if we want to make the difference between the plaquette and the phase and the other phases. One quantity which is different is the product of spins around a loop: . This quantity should be positive in the and in the plaquette phase, and should be negative in the
step
phase. It should be negative but small in the case of the three-three phase. We find a value around for which clearly points toward the plaquette or the phase.

## Vi Conclusion

The phase diagram of the spin antiferromagnetic Heisenberg model on the square lattice with three-site interaction turns out to be extremely rich and complex, with, according to the present results, the possible succession of seven phases, as summarized in Fig. 14. Already for classical spins the situation is more complicated than for the equivalent spin-1/2 model, which has only one transition between a Néel phase and a highly degenerate phase, whereas the present model has two classical phases at small , a Néel phase and an up-up-down-down phase, followed by a highly degenerate ground state. Besides, a mean-field treatment based on a factorized wave-function suggests that the phase diagram is actually more complicated, with the appearance of three helical phases with a predominantly magnetic order parameter and a tiny quadrupolar component.

Let us now have a critical look at the results. The sequence and stability of the classical phases is in our opinion quite robust: all the phases are stable with respect to semiclassical fluctuations in their range of stability, and the order-by-disorder selection of the plaquette phase for not too large and of the phase for large is plausible in view of the linear-spin wave results, even if, as usual, the final answer would require to go beyond linear-spin wave theory.

What is less clear however is the fate of the helical phases. Indeed it is well known that, quite generally, quantum fluctuations tend to favor collinear structures because the harmonic spectrum is softer, and it is not excluded that the helical phases shrink to the point where they disappear altogether when quantum fluctuations are included, something we have not attempted to do in this paper. If that turned out to be the case however, the system is likely to develop quantum spin liquid phases at the transition between the Néel and the up-up-down-down phases as well as at the transition between the up-up-down-down and the plaquette phases since the correction to the local moment diverges in the up-up-down-down phase upon approaching 1/12 and in the plaquette phase upon approaching 1/4. This very interesting alternative, as well as the other open issues summarized above, are left for future investigation.

## Vii Acknowledgments

F. Michaud would like to thank T. Coletta for useful discussions about spin waves. The authors acknowledge the Swiss National Fund and MaNEP.

## References

- (1) C. Lacroix, P. Mendels, and F. Mila, Introduction to Frustrated Magnetism: Materials, Experiments, Theory, Springer Series in Solid-State Sciences (Springer, na, 2011).
- (2) I. Dzyaloshinsky, Journal of Physics and Chemistry of Solids 4, 241 (1958).
- (3) T. Moriya, Phys. Rev. 120, 91 (1960).
- (4) R. Coldea et al., Phys. Rev. Lett. 86, 5377 (2001).
- (5) M. Takahashi, Journal of Physics C: Solid State Physics 10, 1289 (1977).
- (6) A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys. Rev. B 37, 9753 (1988).
- (7) R. Bastardis, N. Guihéry, and C. de Graaf, Phys. Rev. B 76, 132412 (2007).
- (8) F. Michaud, F. Vernay, S. R. Manmana, and F. Mila, Phys. Rev. Lett. 108, 127202 (2012).
- (9) C. K. Lai, Journal of Mathematical Physics 15, 1675 (1974).
- (10) B. Sutherland, Phys. Rev. B 12, 3795 (1975).
- (11) L. A. Takhtajan, Physics Letters A 87, 479 (1982).
- (12) H. M. Babujian, Nuclear Physics B 215, 317 (1983).
- (13) I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799 (1987).
- (14) A. V. Chubukov, Journal of Physics: Condensed Matter 2, 1593 (1990).
- (15) G. Fáth and J. Sólyom, Phys. Rev. B 44, 11836 (1991).
- (16) K. Buchta, G. Fáth, O. Legeza, and J. Sólyom, Phys. Rev. B 72, 054433 (2005).
- (17) U. Schollwöck, T. Jolicoeur, and T. Garel, Phys. Rev. B 53, 3304 (1996).
- (18) A. Läuchli, G. Schmid, and S. Trebst, Phys. Rev. B 74, 144426 (2006).
- (19) S. R. Manmana, A. M. Läuchli, F. H. L. Essler, and F. Mila, Phys. Rev. B 83, 184433 (2011).
- (20) K. Harada and N. Kawashima, Phys. Rev. B 65, 052403 (2002).
- (21) H. Tsunetsugu and M. Arikawa, Journal of the Physical Society of Japan 75, 083701 (2006).
- (22) A. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 97, 087205 (2006).
- (23) T. A. Tóth, A. M. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 105, 265301 (2010).
- (24) K. Nomura and K. Okamoto, Journal of the Physical Society of Japan 62, 1123 (1993).
- (25) E. Rastelli, L. Reatto, and A. Tassi, Journal of Physics C: Solid State Physics 19, 6623 (1986).
- (26) P. Chandra, P. Coleman, and A. I. Larkin, Journal of Physics: Condensed Matter 2, 7933 (1990).
- (27) A. Moreo, E. Dagotto, T. Jolicoeur, and J. Riera, Phys. Rev. B 42, 6283 (1990).
- (28) J. Ferrer, Phys. Rev. B 47, 8769 (1993).
- (29) I. G. Gochev, Phys. Rev. B 51, 16421 (1995).
- (30) J. Yang, D.-K. Yu, and J.-L. Shen, Physics Letters A 236, 589 (1997).
- (31) M. Mambrini, A. Läuchli, D. Poilblanc, and F. Mila, Phys. Rev. B 74, 144422 (2006).
- (32) M. Arlego and W. Brenig, Phys. Rev. B 78, 224415 (2008).
- (33) K. Penc and A. M. Laeuchli, in Introduction to Frustrated Magnetism, Vol. 164 of Springer Series in Solid-State Sciences, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, na, 2011), pp. 331–362.
- (34) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).