Semiclassical evidence of columnar order in the fully frustrated transverse field Ising model on the square lattice

# Semiclassical evidence of columnar order in the fully frustrated transverse field Ising model on the square lattice

Tommaso Coletta Institute of Theoretical Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland    Sergey E. Korshunov L. D. Landau Institute for Theoretical Physics RAS, 142432 Chernogolovka, Russia    Frédéric Mila Institute of Theoretical Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
July 15, 2019
###### Abstract

We investigate the zero-temperature phase diagram of the fully frustrated transverse field Ising model on the square lattice both in the classical limit and in the presence of quantum fluctuations. At the classical level (the limit of infinite spin ), we find that upon decreasing the transverse field this model exhibits a phase transition from the fully polarized state into an eight-fold degenerate translational symmetry breaking state. This phase can be identified to correspond to plaquette order in the dimer language and remains the lowest-energy state in the entire range of fields below the critical one, . The eight-fold degenerate solution which corresponds to columnar order in the dimer language is a saddle point of the classical energy. It is degenerate with the plaquette solution at and is only slightly higher in energy in the whole interval . The effect of quantum fluctuations is investigated in the context of a large S expansion both for the plaquette and columnar structures. For this purpose we employ an approximate method allowing to estimate from above the fluctuation-induced correction to the energy of a configuration which at the classical level is a saddle point of the energy, not a local minimum, and find that the harmonic quantum fluctuations show a clear tendency to overcome the energy difference between the two states. For relatively high fields the transition from the plaquette to the columnar state takes place at values of large enough to be in the domain of validity of the harmonic approximation.

###### pacs:
05.50.+q, 71.10.-w, 75.10.Jm

## I Introduction

Fully frustrated transverse field Ising models Moessner et al. (2000); Moessner and Sondhi (2001a) (FFTFIMs) are characterized by the interplay of two ingredients: a frustrated Ising-like interaction between neighboring spins 1/2 and a field term coupled to the transverse component of the spin which is the source of quantum dynamics in the system. In magnetic systems, the term fully frustrated implies that there is an odd number of antiferromagnetic (AF) bonds on each elementary plaquette, which makes it impossible to satisfy simultaneously all bonds on any of them. Such a model was first introduced by Villain under the name of ‘odd model’.Villain (1977) However, the investigation of fully frustrated Ising models has started long before that: the exact solution of the classical AF Ising model on the triangular lattice (evidently belonging to this class) has been constructed already in 1950.Wannier (1950); Houtappel (1950)

Starting from a classical Ising model, one can obtain its quantum version by including a transverse field term. Thus, the FFTFIM is described by the Hamiltonian

 H=∑⟨i,j⟩Jijσziσzj−Γ∑iσxi, (1)

where are spin operators defined on sites of some regular lattice and the first sum runs over the pairs of nearest neighbors on this lattice. The magnitude of the coupling is the same on all bonds of the lattice while its sign is positive for an odd number of bonds of each plaquette and negative for the remaining bonds. All Ising models on the same lattice satisfying this rule can be transformed into each other by a gauge transformation and accordingly are equivalent. Villain (1977) This leaves one with the freedom to choose the most convenient gauge, that is, which bonds (satisfying the “odd rule“) to consider as antiferromagnetic. References Moessner et al., 2000 and Moessner and Sondhi, 2001a review the properties of FFTFIM on a number of periodic two-dimensional lattices.

One of the main interests in FFTFIMs comes from the fact that they are closely related to the quantum dimer model (QDM) of Rokhsar and KivelsonRokhsar and Kivelson (1988) first introduced in the context of high cuprate superconductors. Moessner, Sondhi and Chandra Moessner et al. (2000) showed that in the limit FFTFIMs can be mapped on purely kinetic QDMs (i.e. QDMs whose potential energy is equal to zero) on the dual lattice. A detailed review of the mapping between the FFTFIM on the honeycomb lattice and the purely kinetic QDM on the triangular lattice can be found in Ref. Misguich and Mila, 2008. In contrast to the triangular-lattice QDM, for which the ground state symmetry at is well known, Moessner and Sondhi (2001b); Ralko et al. (2005); *Ralko06; *Ralko07 the nature of the ground state in the purely kinetic square-lattice QDM has long remained controversial, with different numerical studies predicting columnar,Sachdev (1989); Syljuåsen (2006) plaquetteLeung et al. (1996) or mixed typeRalko et al. (2008) order.

In the present work, the competition between the plaquette and columnar orders in the FFTFIM on the square lattice is systematically investigated in the context of a large expansion. The main results of our work were briefly announced earlier in Ref. Wenzel et al., 2012 (by Sandro Wenzel and the authors of the present article) focused mainly on large scale quantum Monte Carlo simulations of the FFTFIM on the square lattice.

Reference Wenzel et al., 2012 argues in favor of the existence of only two phases in the zero-temperature phase diagram of this model as the ratio is varied. For low fields, a state which in terms of the closely related quantum dimer model has columnar character is stabilized. 111A similar conclusion was reached earlier by Jalabert and Sachdev Jalabert and Sachdev (1991) for the three-dimensional classical model (constructed by stacking the fully frustrated Ising model) which, at finite temperatures, can be expected to behave as the zero-temperature quantum FFTFIM. However, as pointed out in Ref. Wenzel et al., 2012, the numerical analysis of Ref. Jalabert and Sachdev, 1991 did not resolve the phase of the order parameter used in this work, without which it is impossible to distinguish the columnar, plaquette and mixed phases from each other. As the field is increased beyond a critical value, the ordering disappears and the system becomes uniformly polarized.

The present work confirms the conclusion on the columnar nature of the ordering in the FFTFIM on the square lattice by showing that, to first order in , quantum fluctuations have the general tendency to stabilize columnar order. This conclusion has to be taken as only indicative for weak transverse field because of strong quantum fluctuations and anharmonic contributions, but for not too weak field it has more solid grounds and therefore can be considered as an analytical confirmation of the results of Ref. Wenzel et al., 2012 .

The paper is organized as follows. Sec. II discusses the classical version of the model corresponding to the limit of infinite spin . The polarized state with all classical spins aligned along the field extends down to . Our analytical analysis both in the critical region just below and in the low field limit, , demonstrates that the classical ground state corresponds to a plaquette structure in the dimer language. The results of the numerical energy minimization for a single plaquette are then used to prove that the same conclusion is valid in the whole interval . On the other hand, the spin configurations corresponding to columnar order in the dimer language are identified to be saddle points of the classical energy. In the vicinity of , the degeneracy between the columnar and plaquette states is lifted only when the expansion of the classical energy is extended up to the eighth order, the difference in energies scaling like .

Sec. III studies the effect of quantum fluctuations around the plaquette and columnar structures. Harmonic fluctuations in the plaquette state are treated via the standard linear spin-wave approximation. The same approach is inapplicable for analyzing the fluctuations in the columnar state as it is a saddle point of the classical energy, not a local minimum. This forces us to restrict ourselves with finding an upper estimate for the zero-point energy of the state with columnar structure. Within the framework of this approximation, the harmonic fluctuations turn out to stabilize the columnar state over the plaquette one in the entire field range for . The amplitude of harmonic fluctuations as well as the relevance for the QDM are also discussed in this section. A short conclusion is given in Sec. IV.

## Ii Classical limit

In this section we discuss the classical ground state of the model for different values of the ratio . Throughout this paper by classical version of the model we mean not the classical fully frustrated Ising model one obtains when putting , but the model whose Hamiltonian can be obtained by replacing in Eq. (1) the Pauli matrices by the corresponding components of vectors of unit lengths. Note that at the conceptual level this simple transformation can be interpreted as consisting of two separate steps. Firstly, the quantum FFTFIM discussed in the introduction is generalized to the case of spins of arbitrary magnitude by rewriting Eq. (1) as

 H=1S2∑⟨i,j⟩JijSziSzj−ΓS∑iSxi, (2)

secondly one takes the classical limit, , in which the fluctuations are suppressed. In this limit, the normalized spin operators are replaced by three-dimensional vectors of unit length. Below we refer to these variables as classical spins. In terms of classical spins the original Hamiltonian reduces to

 H=∑⟨i,j⟩Jijnzinzj−Γ∑inxi (3)

which has exactly the same structure as (1), but with the spin operators replaced by classical spins with .

In contrast to the analogous model on the triangular lattice, the fully frustrated Ising model on the square lattice does not have any “natural” gauge, hence it is convenient to discuss the structure of the ordered phases in this model in terms of some gauge-invariant variables. For this purpose, following Refs. Misguich and Mila, 2008,Wenzel et al., 2012 and Coletta et al., 2011, we introduce the quantity

 dij=12(1+JijJnzinzj)∈[0,1], (4)

which is equal to 0 if the energy of the bond is negative (that is, is equal to ) and to 1 if this bond is frustrated, that is, has positive energy . In terms of the quantum dimer model Rokhsar and Kivelson (1988) to which the spin-1/2 quantum Ising model (1) is equivalent Moessner et al. (2000) in the limit , the variables defined in Eq. (4) correspond to the occupation number of dimers on the bond of the dual lattice (on which the dimer model is defined) crossing the bond on the Ising model lattice. For this reason, the variables can be called dimer densities.

In the remaining part of this section we analyze the structure of the ground state of the classical model (3) on the square lattice for different values of the ratio . When it is necessary to specify a particular gauge, we always use the simplest periodic gauge in which the antiferromagnetic bonds occupy every second vertical line of the lattice whereas all other bonds are ferromagnetic. This pattern of couplings is depicted in Fig. 1 where dashed lines represent and solid lines . In the chosen gauge the lattice is split into two nonequivalent sublattices and accordingly its unit cell has a rectangular shape and includes two sites, see Fig. 1(a).

### ii.1 Large Γ/J case and critical modes

In this section we discuss the structure of the classical ground state of the model (3) for large values of the ratio . Since for any , the energy of any spin with can be decreased by rotating it around the axis until vanishes, in the lowest energy states all spins must be oriented in the plane. Accordingly, in all configurations that we have to compare, the classical spins can be parametrized by their projections on the axis,

 ni=(√1−(nzi)2,0nzi) (5)

which allows one to rewrite the expression for the classical energy given by Eq. (3) as

 E=∑⟨i,j⟩Jijnzinzj−Γ∑i√1−(nzi)2. (6)

It is evident that in the infinite transverse field limit the minimum of the energy is achieved when . However, upon decreasing , the spins can acquire a small component in the direction. Expanding the square roots in Eq. (6) in powers of yields

 E=Epol+∑⟨i,j⟩Jijnzinzj+Γ∑i[(nzi)22+(nzi)48+(nzi)616+5(nzi)8128+…], (7)

where is the energy of the fully polarized state, being the total number of lattice sites. Introducing the Fourier transforms of via [with index labeling the sublattice to which site belongs, see Fig. 1(a)], the energy per site of the system with periodic boundary conditions can be rewritten as

 E=−Γ−J4∑q,m,m′nz−q,m[M(q)−ΓJ]m,m′nzq,m′+Γ16T4+Γ32T6+5Γ256T8+… (8)

where denotes the sum

 Tμ=2∑m=1μ∏ν=1(∑qνnzqν,m)⋅∑Gδ∑μτ=1qτ,G , (9)

where the last sum runs over the reciprocal lattice of the rectangular lattice defined by the two-site unit cell of figure 1(a). in Eq. (8) is the matrix

 M(q)=(−2cosqz1+e−i2qx1+ei2qx2cosqz). (10)

The eigenvalues of are

 λ±(qx,qz)=±√2√2+cos2qx+cos2qz. (11)

The largest eigenvalues are given by . Therefore, the quadratic form obtained by truncating Eq. (8) to the second order in is positive definite as long as . This defines the critical value of the transverse field at which the fully polarized state becomes unstable. The wavevectors of the critical modes, and , indicate that an instability occurs towards a phase with a unit cell which is twice as large and contains four sites, see Fig. 1(b).

These critical modes are real and can be parametrized as

 (nzqA,1,nzqA,2)=nzAu(+)A,(nzqB,1,nzqB,2)=nzBu(+)B, (12)

where the real coefficients and are their amplitudes and

 u(+)A=(sinπ8,cosπ8),u(+)B=(cosπ8,sinπ8) (13)

are the normalized eigenvectors of and associated to the eigenvalue . Since the degeneracy of the two critical modes is not related to symmetry, one can expect it to be removed in the higher orders of the expansion of the classical energy (6).

The minimization of

 E(4)0=Γ−Γc4[(nzA)2+(nzB)2]+3Γ64[(nzA)2+(nzB)2]2, (14)

the sum of the second- and fourth-order contributions of the critical modes to Eq. (8), fixes only the combination

 (nzA)2+(nzB)2=83Γ(Γc−Γ), (15)

whereas the ratio between and remains undefined. It follows from Eq. (15) that, in the leading order, and .

Taking into account the sixth-order contributions to the energy one obtains

 E(6)0=E(4)0+5Γ256[(nzA)2+(nzB)2]3. (16)

As , depends only on . The dependence on the relative strength of the two amplitudes appears only in the eighth-order contribution to the energy

 E(8)0=E(6)0+5Γ8192[17(nzA)8+84(nzA)6(nzB)2+70(nzA)4(nzB)4+84(nzA)2(nzB)6+17(nzB)8]. (17)

The minimization of Eq. (17) fixes the ratio . This is readily seen when parametrizing the amplitudes of the critical modes as

 (18)

With this choice, the sum of the second- to eighth-order contributions to the energy becomes

 E(8)0=Γ−Γc4n2z+3Γ64n4z+5Γ256n6z+5Γ8192n8z(35−cos8ϕ)2 (19)

which, for any value of , is minimal when . This selects eight degenerate solutions with . To the lowest order in , the corresponding real space spin configurations are

 nz1=nzsin[π8(2p+1)],nz2=nzsin[π8(2p+3)],nz3=nzsin[π8(2p+5)],nz4=nzsin[π8(2p+7)], (20)

where , and where the indices that run from to keep track of the four sublattices [see Fig. 1 (b)]. In terms of the dimer density distribution, the eight solutions correspond to the four states with so-called columnar structure. The example of such a state presented in Fig. 2(b) corresponds to the solution with in Eq. (20). These four columnar coverings are equivalent, that is, they are related to each other by symmetry. Each of them corresponds to two different values of which differ by . In terms of classical spins the shift of by corresponds to reflecting all spins with respect to the plane.

This approach is based on the assumption that all other modes contribute to the energy expansion only in the higher orders. However, it is known from the analysis of the FFTFIM on the honeycomb lattice, Coletta et al. (2011) that when it is necessary to go beyond the fourth order of the Ginzburg-Landau expansion, by restricting the calculation just to the subspace of the critical modes of , one may miss some contributions which are equally essential for the structure selection. We shall now show that in the considered model, some second-, fourth- and sixth-order terms involving noncritical (gapped) modes also make contributions that are essential for determining the optimal value of .

The momenta of the critical modes enforce a four-sublattice structure of the ground state which allows for the excitation of other modes compatible with the same periodicity. Since the only wave vectors consistent with this periodicity are those of the critical modes, the momenta of the relevant noncritical modes will also be and . The most important terms coupling the critical modes with extra modes are expected to be linear in the amplitudes of these extra modes and of a higher order in the amplitudes of critical modes.

Let us denote the Fourier coefficients associated to these extra modes by

 (¯nzqA,1,¯nzqA,2)=¯nzAu(−)A,(¯nzqB,1,¯nzqB,2)=¯nzBu(−)B, (21)

where the real coefficients and are their amplitudes and

 u(−)A=(−cosπ8,sinπ8),u(−)B=(−sinπ8,cosπ8) (22)

are the normalized eigenvectors of and associated to the eigenvalue . The terms in the energy functional to sixth-order that are linear and harmonic in the extra modes are

 E(6)1=Γc+Γ4[(¯nzA)2+(¯nzB)2]+Γ16[(nzA)3¯nzA+3(nzA)2nzB¯nzB−3(nzB)2nzA¯nzA−(nzB)3¯nzB]+3Γ32[(nzA)2(¯nzA)2+3(nzA)2(¯nzB)2+3(nzB)2(¯nzA)2+(nzB)2(¯nzB)2+4nzAnzB¯nzA¯nzB]+3Γ64[(nzA)5¯nzA+5(nzA)4nzB¯nzB−5(nzB)4nzA¯nzA−(nzB)5¯nzB]. (23)

The dominant contribution to the amplitudes of the subcritical modes is captured by the variation of the first line of Eq. (23) with respect to and . This yields

 ¯nzA≈116[−(nzA)3+3(nzB)2nzA],¯nzB≈116[(nzB)3−3(nzA)2nzB], (24)

indicating that the lowest order contribution to these quantities is . Note that in Eqs. (24) has been replaced by . This has been done because we are focused only on the leading order contributions in . Injecting solution (24) into (23) one obtains

 E(6)1=−Γcn6z512+3Γcn8z(−8+cos8ϕ)8192. (25)

Hence is of same order as . Combining the terms which are of order coming from and and which depend on yields

 18192Γcn8z2cos8ϕ. (26)

Since the prefactor of in this expression is positive, the minimum of the energy is achieved when . This defines eight degenerate solutions with . The corresponding real space spin configurations can be obtained replacing by in Eq. (20). However, the energetic stability of the solution in comparison with other solutions is ensured by the form of the higher-order corrections to (20) which are proportional to . All eight degenerate solutions corresponding to have one of the spins in the four-sublattice structure aligned with the transverse field. All nearest neighbors to the polarized spin have the same component which is greater than that of next nearest neighbors [Fig 2(a) presents an example of a state with such a structure corresponding to the solution with ]. In the dimer representation, these solutions correspond to four equivalent plaquette structures which can be obtained from one another by translation. The rotational symmetry of the plaquette state is whereas that of the columnar state is .

In this approach only the leading contributions to the amplitudes and in powers of are considered. Hence, one may wonder if the higher order corrections to these amplitudes may lead to other dependent terms relevant for the structure selection. We verified that this is not the case and that to order all dependent terms are captured by the expression in Eq. (26). This was done by observing that the energy difference between the columnar and plaquette structures obtained by numerical calculation of their energies in the vicinity of is equal to the prefactor in expression Eq. (26) multiplied by two since for the columnar state and for the plaquette state .

From this calculation we deduce that just below the critical field the classical solution with the lowest energy has the plaquette structure and that the energy difference between the plaquette and columnar states is of the order of per site. In the next section we analyze the properties of the ground states in the opposite limit of weak fields.

### ii.2 Weak field case

In the weak field limit, the structure of the ground states can be also studied analytically. This allows one to understand in simple terms why part of the spins are completely polarized along the field. At zero field, the energy of a single square plaquette is given by

 E=p∑m=1Jm,m+1nzmnzm+1 (27)

where subscript numbers the spins along the perimeter of the plaquette, and .

It has been proven in Ref. Coletta et al., 2011 that for the minimum of this expression is achieved in the states in which one of the spins (for example, the one with ) has an arbitrary orientation, whereas the remaining spins are all directed along the axis with

 nzm+1=−sign(Jm,m+1)nzm=±1, (28)

which corresponds to minimizing the energy on each of the bonds connecting them (as in an open chain of length ). This proof does not rely on the particular value of and is applicable for any , in particular for . The ground state energy of a single plaquette is therefore equal to , which for gives .

On the honeycomb lattice (), it is possible to minimize simultaneously the energy of each of the three plaquettes sharing a given site only if the spin on this site is directed along (). For this reason, the ground state manifold of the model (3) with coincides with that of the discrete Ising model with the same Hamiltonian but Coletta et al. (2011)

In contrast to that, on the square lattice, the ground state manifold of the model (3) with is essentially wider than that of its discrete version. Different grounds states can be obtained by taking a ground state of the discrete model and rotating in an arbitrary way some of the spins for which the sum describing their interaction with the neighboring spins is equal to zero. However, two rotated spins cannot be the nearest neighbors of each other.

The application of a weak transverse field strongly suppresses the degeneracy of the ground states. It is evident that in order to decrease the energy, all spins which in the absence of the field are free to rotate now have to be directed along the field. This gives a negative contribution to the energy equal to times the number of such spins. Therefore, in the states minimizing the energy this number has to be as large as possible. Since each polarized spin sits in the middle of a cell along whose perimeter the spins have to be parallel (or antiparallel) to the axis, the maximal fraction of spins polarized by the field is equal to one quarter. In such a case, on each plaquette one of the four spins is polarized by the field. In particular, this can be realized in the state with the four-sublattice structure. In the next section we show that for any field the global minimum of energy can also be achieved in the framework of the four-sublattice ansatz.

### ii.3 Reduction to the single-plaquette problem

The aim of this section is to show that the minimum of the energy for any field can be achieved in a state with the four-sublattice structure. To demonstrate this, it is convenient to rewrite the classical energy (3) as

 E=∑αEα, (29)

where

 Eα=−J2[nzj1(α)nzj2(α)+nzj2(α)nzj3(α)+nzj3(α)nzj4(α)−nzj4(α)nzj1(α)]−Γ44∑m=1√1−[nzjm(α)]2, (30)

the index runs over all four-site plaquettes of the lattice, and denotes the site belonging to plaquette and to sublattice number [, see Fig. 1(b)]. Each term depends only on four variables which are associated to the four sites belonging to plaquette . Below is often called the energy of plaquette .

It is evident that if it is possible to minimize simultaneously all terms in the sum (29), this will give the absolute minimum of energy. Since Eq. (30) has exactly the same structure for all plaquettes, this aim is easily achieved by minimizing for a single plaquette and then assuming that the state has the four-sublattice structure in which all variables defined on the sites belonging to sublattice have the same value.

The reasoning above does not prove that all ground states have to have the four-sublattice structure. To check if this is really so it is necessary first to find what spin configurations minimize , which is the topic of the next section.

### ii.4 Single plaquette energy minimization

In Sec. II.2 we have shown that in the limit the energy is minimized when on each plaquette one of the spins is fully polarized along the field. On the other hand, the results of Sec. II.1 suggest that the same property holds also when approaches . It seems plausible that the state with such a structure minimizes the energy for any . To check this, let us first find the explicit form of the spin configuration minimizing the energy of a single plaquette defined by Eq. (30) under the assumption that one of the four spins is fully polarized. Let us denote by the site where the spin is polarized along the field [that is, ], by and its two neighbors and by the site diagonally opposite to . Under this assumption the energy of a single plaquette can be written as

 Eα=12nzl(Jljnzj+Jlknzk)−Γ4(1+nxj+nxk+nxl). (31)

The spin configuration minimizing (31) at any field is given by:

 (32)

with the signs of and determined by the sign of ,

 sign(nzj)=−sign(Jjlnzl),       sign(nzk)=−sign(Jklnzl). (33)

The solution described by Eqs. (32) and (33) is valid regardless of the position of the fully polarized spin with respect to the antiferromagnetic bond as long as and denote the sites neighboring the polarized spin and the site diagonally opposite to it.

It follows from Eqs. (32) and (33) that in such a state the force acting on the spin at site is equal to zero,

 Jijnzj+Jiknzk=0. (34)

This implies that the same energy can be obtained whatever the value of , and not only for . We have verified with the help of a numerical minimization of as a function of four variables that the solutions found above are not local but global minima of . In total, for a given plaquette there are eight spin configurations minimizing its energy, which are related to each other by symmetries. The factor four comes from the possibility to choose a site at which the spin is fully polarized and the additional factor two comes from the possibility to choose the sign of .

### ii.5 Classical ground states and their dimer representation

After finding eight spin configurations minimizing the energy of a single plaquette we can immediately construct eight ground states of the model on the infinite square lattice. As it has been already mentioned in Sec. II.3, this can be achieved by choosing one of these configuration for one particular plaquette and after that assuming that the ground state has a four-sublattice structure, that is, on each of the four sublattices [see Fig. 1(b)] all spins have the same orientation. It turns out to be impossible to construct any other state with the same energy, because on each of the bonds the two spins belonging to it always have different orientations. Therefore, the choice of one of the eight configurations on one plaquette uniquely determines which configurations have to be chosen on neighboring plaquettes, and so on, which reproduces nothing else but a state with the four-sublattice structure. Accordingly, the ground state manifold is restricted to the eight four-sublattice states.

In terms of the gauge-invariant variables defined by Eq. (4) (which in the quantum case can be interpreted as dimer densities on the bonds of the dual lattice crossing the corresponding bonds of the original lattice) these ground states correspond to the plaquette dimer structure exhibiting a rotational symmetry. In this structure the dimer densities on the bonds of the dual lattice surrounding the sites of the original lattice on which the spins are fully polarized by the field are all equal to , while the dimer densities of all remaining bonds of the dual lattice are again all equal to each other but are smaller than . Hence the bonds with the highest dimer densities form a regular pattern of square plaquettes, as shown in Fig. 2(a), which explains the origin of the widely used term “plaquette phase”. The plaquette pattern is characterized by a four-fold degeneracy related to translations, however each pattern corresponds to two different ground states in terms of classical spins. These two states are transformed into each other by a reflection of all the spins with respect to plane (in other terms, by changing the signs of all ).

The columnar state corresponds to the dimer density pattern of the type depicted in Fig. 2 (b). In the framework of the classical model, this state appears if the analytical analysis in the vicinity of takes into account only the critical modes. The columnar pattern has a four-fold degeneracy with one factor two related to translations and another one to rotations by 90 degrees. Naturally, in terms of spins, the degeneracy is doubled, because gauge-invariant variables are not sensitive to a simultaneous change of signs of all . All eight spin realizations of the columnar state have a periodicity compatible with the four-sublattice structure [see Fig. 1 (b)]. However, for each of them, it is possible to choose the gauge in such a way that the spin configuration becomes compatible with a larger translation group [that is, has a two-sublattice structure]. In particular, the choice of gauge depicted in Fig. (1) is the one realizing the doubling of the translation group of the columnar solution represented in Fig. 2 (b). For other columnar structures this gauge is obtained from that of Fig. (1) by a translation of one lattice parameter in the direction, or by a rotation by .

Within the 4-sublattice ansatz the columnar solution corresponds to a saddle point of the classical energy lying between two plaquette states. More precisely, in the phase space of four-sublattice structures there exists an almost degenerate circle of low energy states which has eight equivalent minima (the plaquette states) and between them eight degenerate maxima corresponding to the columnar structures. In the vicinity of , this family of low energy states can be parametrized by Eqs. (20) treating as a continuous variable. Insofar the emergence of the eight-fold degeneracy in the FFTFIM was discussed Moessner and Sondhi (2001a) only in relation with its appearance in the framework of the Ginzburg-Landau expansion for the free energy of the classical three-dimensional version of this model.Blankschtein et al. (1984) Our analysis has revealed that such a degeneracy has an even more evident origin.

The energy barrier separating “neighboring” plaquette states (that is, the difference in energy between the plaquette and columnar states) is always relatively small. In the vicinity of its value (per site) is of the order of and is much smaller than the characteristic energies of the plaquette and columnar state (counted off from that of the paramagnetic one) which both are of the order of . In the low field limit the height of the barrier is of the order of and is small in comparison with , the energy scale characterizing the manifold of the four-sublattices states. Between the analytically tractable limits, we verified numerically that the energy difference between the columnar and plaquette states never exceeds per site. This value is achieved at , where the characteristic energy scale is still of the order of . In the following section we show that fluctuations around the considered periodic structures change the classical picture and stabilize the columnar solution over the plaquette one.

## Iii Semiclassical approach

### iii.1 Linear spin-wave approximation

Having discussed the classical phase diagram of the model we now examine the effect of quantum fluctuations on the competition between the plaquette and columnar states. Quantum fluctuations are investigated in the context of the large expansion (spin-wave approximation). Without loss of generality we consider, out of the eightfold degenerate plaquette solutions, the one having the same orientation of the spins on the sublattices and and fully polarized spins on sublattice [see Fig. 2(a)], with the signs of and chosen to be positive. The columnar solution used in the calculations is the one having the simplest structure in the gauge used in this work and consisting of just two sublattices, see Fig. 2(b). These structures are the starting classical solutions used in our spin-wave calculation.

The first step in the construction of the spin-wave expansion consists in rotating the spin operators on each site around the axis,

 Sxj=nzm(j)Sx′j+nxm(j)Sz′j,Syj=Sy′j,Szj=−nxm(j)Sx′j+nzm(j)Sz′j, (35)

Here denotes the number of the sublattice to which site belongs, whereas and are the two components of the classical spin lying in the plane and belonging to the sublattice. The purpose of this rotation is to achieve a situation where on each site the axis is aligned with the classical spin. After that the spin operators in the rotated frame are mapped to bosons via the standard Holstein-Primakoff (HP) transformation Holstein and Primakoff (1940) according to

 Sx′j+iSy′j=√2S ⎷1−a†jaj2Saj,Sx′j−iSy′j=√2Sa†j ⎷1−a†jaj2S,Sz′j=S−a†jaj. (36)

The spin-wave approximation is based on replacing the square root in Eqs. (36) by its expansion in powers of , the number of HP particles divided by the value of the spin being the small expansion parameter. When truncating the expansion to harmonic level, the original Hamiltonian (2) expressed in terms of Holstein-Primakoff bosons can be split intro three contributions

 H=H(0)+H(1)+H(2), (37)

where is the classical energy of the system, and respectively contain only terms which are linear and quadratic in bosonic operators. Furthermore, the coefficients in are proportional to the derivatives of the classical energy with respect to the spin orientations. For this reason, both for the plaquette solution, which is the ground state of the classical energy, and for the columnar solution, which is a saddle point of the classical energy.

After Fourier transformation the quadratic bosonic Hamiltonian of the plaquette structure can be expressed as

 HP=NEP+∑q[→a†q¯HPq→aq+Δ], (38)

where is the classical energy per site of this configuration. Its transverse field dependence is given by

 EP=−Γ4−√(4J2+Γ2)(16J2+Γ2)8J. (39)

In Eq. (38) and is the matrix

 ¯HPq=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝AE⋆q0D⋆q0E⋆q0D⋆qEqCG⋆q0Eq0G⋆q00GqAFq0Gq0FqDq0F⋆qBDq0F⋆q00E⋆q0D⋆qAE⋆q0D⋆qEq0G⋆q0EqCG⋆q00Gq0Fq0GqAFqDq0F⋆q0Dq0F⋆qB⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (40)

with coefficients

 A=12S(2Jnz1nz2+Γnx1),B=Γ2S,C=12S(4Jnz1nz2+Γnx2),Dq=J4Snx1(1+ei2qz),Eq=−J4Snx1nx2(1+ei2qx),Fq=−J4Snx1(1+ei2qx),Gq=−J4Snx1nx2(1+ei2qz). (41)

In the above expressions, and with are the components of the classical spins in the plaquette state considered in this spin-wave calculation, see Fig. 2(a). Their values are given by Eqs. (32) with the identification and . The term in Eq. (38) is defined by .

Similarly, the quadratic bosonic Hamiltonian of the columnar structure is of the form

 HC=NEC+∑q[→a†q¯HCq→aq+~Δq], (42)

where is the classical energy per site of the columnar structure. The energy per site of the structure having the symmetries of the columnar state [Fig. 2 (b)] takes the form

 EC(n1,n2)=−Jnz1nz2−J2(nz2)2+J2(nz1)2−Γ2(nx1+nx2), (43)

and, accordingly, in Eq. (42) denotes the minimum of with respect to and . In Eq. (42) and is the matrix

 (44)

with coefficients

 ~Aq=12S[2Jnz1(nz2−nz1)+Γnx1+J(nx1)2cosqz],~Bq=12S[2Jnz2(nz2+nz1)+Γnx2−J(nx2)2cosqz],~Cq=J2S(nx1)2cosqz,~Dq=−J2S(nx2)2cosqz,~Eq=−J4Snx1nx2(1+ei2qx). (45)

Naturally, the classical spin components in Eqs. (45) are those minimizing . For all spins are almost parallel to the axis, the deviations from it depending on as

 nx1≈(Γ/J)1/3 ,nx2≈Γ/4J, (46)

while just below all spins are almost parallel to the axis, the components behaving as

 (47)

The term in Eq. (42) is defined by .

In the two classical structures that are compared, the fluctuation Hamiltonians (38) and (42) do not contain terms which are linear in the Holstein-Primakoff bosons and therefore are purely quadratic. In the case of the plaquette state, the correction to the classical energy is obtained by diagonalizing the fluctuation Hamiltonian. This is done via a standard Bogoliubov transformation of (38) which yields the dispersion relations and the correction to the classical energy.

The computation in the case of the columnar state is more involved. This time the starting classical configuration is not a minimum but a saddle point of the classical energy. This results in a quadratic fluctuation Hamiltonian which is not positive definite. An attempt to diagonalize the quadratic Hamiltonian via a Bogoliubov transformation in this case would yield a spectrum which is not well defined for some values of momenta. For the regions in the Brillouin zone near , the spectrum would not be real but would take non-physical complex values. In such a situation, the correction to the classical energy cannot be computed. The fact that the spectrum is not well defined results from the truncation of the spin-wave approximation to the harmonic level. If a state with columnar structure is to become the ground state when quantum fluctuations are fully taken into account, the higher-order terms of the spin-wave expansion must enforce the spectrum to have real frequencies.

In order to avoid the explicit inclusion of the higher-order terms into the analysis, we proceed in the following way. We add to the Hamiltonian , Eq. (42), the positive term,

 V=δS∑i[S−Sz′(i)i], (48)

describing the presence on each site of an auxiliary field oriented along the direction of the classical spin at this site, that is, along . In Eq. (48), the summation is taken over all sites and parametrizes the auxiliary field strength. In terms of Holstein-Primakoff bosons, takes the form . The addition of shifts the coefficients and in the columnar-state spin-wave Hamiltonian, Eq. (45), up by the same value,

 ~Aq→~Aq+δ2S,~Bq→~Bq+δ2S, (49)

but does not change any other coefficients. Accordingly the matrix and the term become

 ¯HCq→¯HCq+δ2SI,~Δq→~Δq−δS. (50)

Hence the effect of the auxiliary field is to shift all eigenvalues of up by .

The value of this auxiliary field is adjusted to obtain a fluctuation Hamiltonian which is non-negatively defined, allowing it to be diagonalized by a Bogoliubov transformation. This is done by choosing for any given ratio such that the lowest eigenvalue of is equal to zero. The resulting spectrum has real and positive frequencies with soft modes only at the wavevector . The advantages of this approach are the following: the addition of to the Hamiltonian does not change the classical energy of the state considered [ is left unchanged] and allows to obtain dispersion relations which are physically meaningful. Furthermore, is a strictly positive contribution to the Hamiltonian, hence the corrections to the energy of the columnar state computed with this approach provide an upper bound for the energy of this state at order .

Now we can compare the energy of the plaquette state corrected by the inclusion of the harmonic fluctuations with an estimate from above for the energy of the columnar state calculated to the same order in . We find that as a function of and , the upper bound for the energy of the columnar state is lower than the energy of the plaquette state in a significant parameter range, see Fig. 3.

This approach shows that quantum fluctuations to order easily overcome the energy difference between the classical energies of the plaquette and columnar structures. This is not surprising since, as already discussed in Sec. II.5, the difference of classical energies of the two states is very small compared to the energy of the plaquette state. The auxiliary field required to obtain real and positive frequencies for the columnar state ranges from to approximately and is maximal for .

It has to be noted that at both columnar and plaquette states are the ground states of the Hamiltonian and remain degenerate for all values of . The approach developed in this section suggests that for any finite value of as soon as is turned on the harmonic fluctuations favor the columnar state over the plaquette state. In fact, at low transverse fields, we find that the critical value of at which the transition between the two phases occurs scales like (see inset of Fig. 3). Naturally, for small (including ) the higher-order corrections in may be important. However, when one disregards them, our analysis predicts that at the columnar state is stabilized for all fields below .

### iii.2 Fluctuations and the applicability of the harmonic approximation

In this section we present the analysis of the amplitude of fluctuations and discuss the transverse field domain for which the harmonic approximation is justified. In fact it is reasonable to wonder whether fluctuations completely melt the classical order of the structures considered. In spin problems, one can distinguish longitudinal and transverse fluctuations (definitely related to each other). By longitudinal fluctuations we refer to the suppression of the spin projection along its average direction. Since the spin projection along the classical direction is given by