Order by disorder and phase transitions in a highly frustrated spin model on the triangular lattice

A. Honecker, D. C. Cabra, H.-U. Everts, P. Pujol, 5 and F. Stauffer Institut für Theoretische Physik, Georg-August-Universität Göttingen, 37077 Göttingen, Germany Departamento de F́ısica, Universidad Nacional de La Plata F́ısica and Instituto de F́ısica La Plata, C.C. 67, (1900) La Plata, Argentina Institut für Theoretische Physik, Leibniz Universität Hannover, 30167 Hannover, Germany Laboratoire de Physique Théorique, Université de Toulouse, UPS, (IRSAMC), 31062 Toulouse, France CNRS, LPT (IRSAMC), 31062 Toulouse, France RBS Service Recherche & Développement, Strasbourg – Entzheim, 67836 Tanneries Cedex, France (Dated: August 26, 2011)


I. INTRODUCTION
The study of frustrated quantum magnets is a fascinating subject that has stimulated many studies within the condensed matter community in recent years. [1][2][3] Such systems are assumed to be the main candidates for a rich variety of unconventional phases and phase transitions such as spin liquids and critical points with de-confined fractional excitations. 4 Frustration can also play an important rôle in classical systems. The phenomenon of order-by-disorder 5,6 is the perfect example where the interplay of frustration and fluctuations produces the emergence of unexpected order. Order-by-disorder implies that a certain low-temperature configuration is favored by its high entropy, not by its low energy. Order-bydisorder can also occur in a quantum system, where a naïve argument suggests that quantum fluctuations play the same rôle as thermal fluctuations in the classical system, albeit there are counterexamples where their rôle is in fact quite different. 7 A particularly illustrative example is provided by the spin 1/2 antiferromagnet on the kagome lattice.
A spin gap appears to be present both at zero magnetization 2,[8][9][10][11][12][13][14] (see, however, Refs. [15][16][17] and at 1/3 of the saturation value where it gives rise to a plateau in the magnetization curve. 7,[18][19][20][21] One would be tempted to believe that the nature of the ground state is similar in both cases. However, whether the ground state at zero field is ordered or not is still under debate and also the existence of a plateau in the isotropic spin 1/2 Heisenberg model at magnetization 1/3 has been questioned recently. 22,23 Nevertheless, the existence of a plateau at magnetization 1/3 is quite clear for easy-axis exchange anisotropies 7,19 and, using a correspondence with a quantum dimer model on the honeycomb lattice, 24 the ground state is identified as an ordered array of resonating spins. 7,25 In this paper we study an effective model that arises in the strong trimerization limit of the spin 1/2 kagome antiferromagnet. 26 This model has played an important rôle in analyzing the zero-field properties of the kagome antiferromagnet, 27,28 but here we will focus on magnetization 1/3 of the Heisenberg model, corresponding to full polarization of the physical spin degrees in the effective model. Thus, we are left with the chirality degrees of freedom of the original antiferromagnet which we will treat as 'spin' variables. In this sense our spin system can be considered as a purely orbital model similar to compass models recently considered in the literature (see, e.g., Refs. [29][30][31][32][33][34][35][36][37]. As an experimental realization of this model a system of spin-polarized fermions trapped in a trimerized optical kagome lattice at 2/3 occupancy has been suggested. [38][39][40] Beyond the particular realizations of our model its very rich physics which results from the interplay between classical and quantum fluctuations and frustration makes it an interesting model in its own right. As will be shown in this paper, a Hamiltonian with an (unusual) discrete symmetry but with a continuous degeneracy of the classical ground state, as it would be expected for a Hamiltonian with a continuous symmetry, is just one aspect of the rich phenomena emerging from this model.
The present paper is organized as follows: in section II we present the Hamiltonian and the symmetries of the classical and spin 1/2 cases. The Hamiltonian can be defined for both signs of the coupling constant J. We deliberately discuss in parallel the two cases throughout the entire paper to point out their similarities and differences. The spin 1/2 case is then treated in section III by means of exact diagonalization techniques, and we argue that a finite-temperature phase transition takes place. Since exact diagonalization can access only small lattices, we move to the classical model in section IV. We study in detail the manifold of lowest-energy configurations and their corresponding spin-wave spectra. The effect of soft modes in the order-by-disorder selection mechanism is argued to be the origin for the phase transition of the J > 0 case, in contrast to the J < 0 case, where the transition is of a more conventional purely energetic origin. In section V we apply Monte-Carlo techniques to the classical model and determine the transition temperature for J > 0 and J < 0. We also analyze the universality class of the transition, however, only for J < 0 since the transition temperature for J > 0 turns out to be so low that it is difficult to access. Finally section VI is devoted to some concluding remarks and comments.

A. Hamiltonian
We will study the Hamiltonian where with the third root of unity ω = e 2πi/3 . The sums in (1) run over the bonds of a triangular lattice, each corresponding to one of the three distinct directions of the lattice, as sketched in Fig. 1(a). The Hamiltonian (1) arises as an effective Hamiltonian for the trimerized kagome lattice, sketched in Fig. 1(a) behind the triangular lattice. Our notation follows the derivation from the half-integer spin Heisenberg model for the case where the remaining magnetic degrees of freedom are polarized. 26 In this case, there are two pseudospin states of opposite chirality for each triangle, as shown in Fig. 1(b). Plain first-order perturbation theory of the S z -S z interactions between two triangles yields Eq. (1) where the exchange constant J is proportional to the inter-triangle exchange constant of the kagome lattice and would thus typically assumed to be antiferromagnetic (J > 0). Note that we have chosen a convenient normalization of J. A similar derivation starting from a Fermi gas with two atoms per trimer also leads to the Hamiltonian (1). 38 Due to the two possible chiralities on each triangle, the pseudo-spin operators S i should be considered as quantum spin-1/2 operators. The derivations 26,38 also suggest a positive J > 0 to be more natural. In this paper we will relax these constraints and, for reasons that will become clear later, consider also classical spins, i.e., unit vectors S i , and the case J < 0.
Note that the Hamiltonian (1) is not symmetric under reflections of the lattice. Our conventions agree with those of Ref. 26 where this Hamiltonian appeared first, while some more recent works 38-41 use a reflected convention for the chirality. Note furthermore that our conventions for J differ by a factor 4 from previous studies of the model (1). 39,40 It may also be useful to represent the Hamiltonian (1) in a more compact form 42 where the vectors e i are indicated in Fig. 1 Models which are very similar to (1) have recently been studied in the context of spin-orbital models (see, e.g., Refs. 29-37).

B. Symmetries
The Hamiltonian (1) has the following symmetries on an infinite lattice: 1. Translations T x , T y along the two fundamental directions of the lattice.
2. Simultaneous rotation R 2π/3 of the lattice and all spins around the z-axis by angles of 2π/3 (the latter rotation amounts to a cyclic exchange of T A i , T B i , and T C i ). 3. A rotation by π around the z-axis in spin space: P : S x i → −S x i , S y i → −S y i while keeping the lattice fixed.

A spatial reflection combined with rotation of all
spins around a suitable axis in the x-y-plane by an angle π. One particular choice is I : combined with a reflection of the lattice along the dashed line in Fig. 1(a). 5. For spin 1/2, there is another symmetry implemented by Conservation of Q means that the number of spins pointing up (or down) along the z-axis is a good quantum number modulo two. This conservation law is most easily verified by observing that the interaction terms in (1) always invert a pair of spins in an eigenbasis of S z .
The choice of factors in (4) ensures that Q 2 = 1. Furthermore, one has R 3 2π/3 = P 2 = I 2 = 1. R 2π/3 and P together generate the abelian group Z 6 ∼ = Z 3 × Z 2 , as described for instance in Ref. 39. The combination of R 2π/3 , P , and I generates the dihedral group D 6 , which is non-abelian (I R 2π/3 I = R −1 2π/3 ). Finally, R 2π/3 and I generate the symmetric group S 3 which can be traced to the point-group symmetry of the underlying kagome lattice. The operators R 2π/3 , P , and I leave the Hamiltonian (1) invariant irrespective of the value of the spin quantum number. Thus, the group D 6 is a symmetry also of the classical variant of the model.
The symmetries P and Q are not present in the underlying kagome lattice, hence they should be specific to the lowest-order approximation. 26,38 Indeed, at least in the derivation from the Heisenberg model one observes that already the next correction 41 breaks the symmetries P and Q. Now let us consider the consequences of the combination of I and Q for the spin-1/2 model on a finite lattice with N sites. Then the relation I S z j I = −S z j leads to Since the eigenvalues of Q are q = ±1, this implies that I is an isomorphism between the subspace with q = −1 and the subspace with q = 1 for odd N and spin 1/2.

III. QUANTUM SYSTEM: EXACT DIAGONALIZATION FOR SPIN 1/2
In this section we will present numerical results for the Hamiltonian (1) with spin 1/2. We impose periodic boundary conditions and use the translational symmetries T x , T y in order to classify the states by a momentum quantum number k. We only consider lattices which do not frustrate a potential three-sublattice order, i.e., only values of N that are multiples of three. For the system sizes N considered already in Refs. 39 and 40, we will use the same lattices. In particular, the N = 12, 21, and 24 lattices are shown in Fig. 9 of Ref. 40. Furthermore, we will consider the N = 27 lattice which can be found, e.g., in Fig. 1 of Ref. 43.
Let us briefly discuss the consequences of the other symmetries mentioned above. We did not make explicit use of P , although it is present for any lattice. However, the symmetry Q (which is also present for any lattice, see (4)) is easily implemented if we work in an S z -eigenbasis. Concerning the rotation R 2π/3 , it is not possible to find lattices for all N such that it is a symmetry. If R 2π/3 is a symmetry, we use it to select one representative k for all equivalent momenta. Finally, the presence of the symmetry I is more delicate. We have performed computer checks and found that most of the lattices under consideration have a suitable spatial reflection symmetry, ensuring that I is a symmetry. The only exception is the N = 21 lattice where there is no such reflection symmetry. Nevertheless, we find the same spectra in the subspaces with q = −1 and q = 1 also for N = 21. Therefore, for N odd we can choose representatives for all symmetry sectors in the subspace with q = 1.
For N ≤ 21 the translational symmetries and Q lead to matrices with dimension up to 49 940 and we can obtain all eigenvalues. Dimensions increase up to 2 485 592 for N = 27. In this case, we have used the Lanczos method to compute the n lowest eigenvalues in each sector. Mainly for reasons of CPU time, we restrict to n ≈ 70 (150) for N = 27 (24) and J > 0.
A. Low-lying spectra Let us first look at the spectra. In order to take degeneracies into account, in Fig. 2 we show the integrated density of states, i.e., the number of states with energy less or equal than ∆E above the ground state.
Panel (a) of Fig. 2 shows results for J < 0. These results extend previously presented results 40 for N = 12, 18, and 21 to higher energies and larger N . One observes that there are at most 8 states for energies ∆E < ∼ 3 |J| with a substantial density of states setting in at higher energies. This suggests a thermodynamic gap ≈ 3 |J|. Fig. 2(b) shows the density of states for J > 0, extending previously published results for N = 18, 21, and 24. 39,40 In this case, we observe a large density of states at substantially lower energies than for J < 0. This large density of states is reminiscent of the large density of nonmagnetic excitations observed in the spin-1/2 Heisenberg antiferromagnet on the kagome lattice, both at zero magnetic field 10,11 and on the one-third plateau. 7 In particular the N = 27 data presented in Fig. 2(b) shows a large density of states for ∆E > ∼ 0.02 J. On the other hand, one observes at most 8 states with ∆E < ∼ 0.012 J in Fig. 2(b) for a given system size N . From these observations we infer that a gap is at most on the order of 0.02 J if present at all.
Since an ordered ground state breaks the symmetry group D 6 , such a ground state should be six-fold degenerate. Indeed, classical and semiclassical considerations predict a six-fold degeneracy in an ordered state (see section IV below). However, there is no clear separation of 6 low-lying states from the remainder of the spectrum, neither for J < 0 (see Fig. 2(a)), and even less so for J > 0   Fig. 2(b)). The considered lattice sizes may be too small to observe the expected low-energy structure of the spectrum. However, correlation functions exhibit pronounced 120 • correlations already on these small lattices. 39,40 B. Specific heat The specific heat C can be expressed in terms of the eigenvalues of the Hamiltonian. Since we have all eigenvalues for N ≤ 21, it is straightforward to obtain the specific heat for all temperatures and both signs of J. Fig. 3 shows the results of the specific heat per site C/N . The case J < 0 is shown in panel (a). There is a finitesize maximum slightly above T ≈ |J|. The large finitesize effects which are still observed here are consistent with a phase transition around T ≈ |J| in which case C should become non-analytic for N → ∞. Because of a possible phase transition, we have tried to obtain a low-temperature approximation to the specific heat for N > 21, J < 0 by keeping low-energy states. However, for N = 24 even 12 462 low-lying states going up to energies as high as ∆E < ∼ 12.6 |J| turned out to yield a specific  heat which has sufficiently small truncation errors only for temperatures T < ∼ 0.9 |J|. This result for N = 24 (also included in Fig. 3(a)) clearly does not include the maximum of the specific heat C.
Now we turn to the case J > 0 which is shown in panel (b) of Fig. 3. In this case, finite-size convergence at high temperatures is much faster than for J < 0. This fast convergence indicates that there is no phase transition associated to the high-temperature maximum (the position of this maximum is at T ≈ 2.1225 J with a value C ≈ 0.105717 N for N = 21). The reduced finite-size effects and the smaller value of C at the maximum reflect that there is a substantially smaller energy scale for J > 0 as compared to J < 0.
For J > 0, a second peak emerges in the specific heat at low temperatures, see Fig. 3(b). In order to investigate this in more detail, we use again the low-temperature approximation for the specific heat obtained from the low-lying part of the spectrum. For N = 24 and 27, we have used a total of 7 029 and 3 906 eigenvalues, respectively (the N = 24 data is included in Fig. 3(b), but it is difficult to see there since it is valid only at very low temperatures). Fig. 4 shows the specific heat divided by temperature (panel (a)) and the entropy per site (panel (b)) in the low-temperature region for J > 0 and system sizes N = 18, 21, 24, and 27. Our result for the specific heat with N = 21 obtained from the full spectrum agrees with a previous result for N = 21 based on approximately 2 000 low-lying states. 40 The finite T = 0 limit of the entropy for N = 21 and 27 in Fig. 4(b) corresponds to the two-fold degeneracy of the ground state for these system sizes, see Fig. 2(b). Although the maximum value of C/T increases with increasing N , there are nonsystematic finite-size corrections to the position of this maximum. Thus, we can only conclude that if there is a finite-temperature ordering transition for J > 0, it should have a very low transition temperature T c < ∼ J/100. shows that there is a remarkably large entropy S/N ≈ 0.2 . . . 0.3 associated to the finite-size lowtemperature peak of the specific heat. This is comparable to the entropy associated to the degeneracy of the classical ground states, see section IV C below. Therefore, this observation lends further support to the interpretation 40 of the low-energy states for S = 1/2 in terms of the classical ground states for J > 0.

IV. CLASSICAL SYSTEM: LOWEST-ENERGY CONFIGURATIONS AND SPIN-WAVE ANALYSIS
We will now proceed with a discussion of the lowenergy, low-temperature properties of the classical variant of the model (1), i.e., we will assume that the S i are unit vectors. We will parametrize the spin at site i by angles α i and γ i : Because the z-components do not contribute to the energy, configurations with extremal energy should have spins lying in the x-y-plane (γ i = 0). The energy E({α i }) for a given set of angles {α i }, γ i = 0 is obtained from (1) by identifying with Ω = 2 π/3. We will further be interested in small fluctuations The energy can then be expanded as Here, E zz is a diagonal quadratic function of the out-ofplane fluctuations˜ i which, to quadratic order, decouples from the relevant degrees of freedom i . The eigenvalues f i of the symmetric matrix M i,j correspond to the spinwave modes. The fact that {α i } describes a ground state implies f i ≥ 0. We will call a mode with f i = 0 'pseudo-Goldstone mode'.
A. Ground states with small unit cell for J < 0 Let us first consider the case J < 0. Then a ground state is given by a certain three-sublattice configuration where the angles α i between different sublattices differ by multiples of 2 π/3. 39,40 Fig. 5 shows such a low-temperature configuration as a snapshot which was taken during a Monte-Carlo simulation (details to be given in section V below). The energy of such states E < class. = 6 J N is invariant under global rotations of the spin configuration in the x-y-plane, i.e., there is a oneparameter family of ground states (note that this invariance under a continuous group is not a symmetry of the Hamiltonian). We parametrize this global rotational degree of freedom by an angle α of the spins on one sublattice. Using a three-site unit cell, we can exploit invariance of this ground state under translations to represent the matrix (8) in Fourier space by the following 3 × 3 matrix where A = e i k1 sin 2 (α) + e i k2 sin 2 (α + Ω) (11) Let us analyze now the effect of the fluctuations by computing the free energy associated with (8). To this end we can compute the partition function where f α i ( k) are the eigenvalues of (9) and Z zz is the Gaussian integral over the N quadratic variables corresponding to the out-of-plane fluctuations. Performing the Gaussian integral we get which yields the free energy as where F zz = − ln Z zz /β. The low-temperature behavior is therefore determined by the following contribution of the fluctuations to the free energy where M is the matrix (9). The result of the integral (15) is shown in Fig. 6. F < (α) is a 2 π/3-periodic function since the spin angles on the different sublattices differ by 2 π/3. Hence, it is sufficient to consider α ∈ [0, 2 π/3]. One observes that F < and thus the lowtemperature limit (14) of the free energy has minima at α = (2 n + 1) π/6, n = 0, 1, . . . , 5. This implies that the 120 • classical ground-state configuration locks in at these angles for T → 0. This lock-in can indeed be verified in histograms of Monte-Carlo simulations, see Ref. 42 for a planar variant of this model and also Fig. 11 below. Lock-in of the classical ground-state configuration at α = (2 n + 1) π/6 follows also from the semiclassical approach. 40 In this approach the ground-state energy is given by The three eigenmodes fi for the 120 • ground state with J < 0, assuming a lock-in of the ground state at α = (2 n + 1) π/6. Note that the shaded surfaces extend only over the first Brillouin zone.
where ω < i ( k, α) are the three sheets of spin-wave (SW) frequencies obtained from the linear Holstein-Primakoff approximation and where the sum over k runs over the magnetic Brillouin zone.
B. Ground states with small unit cell for J > 0 Now we turn to the case J > 0. There is a first ground state 39,40 described by α i = α with an arbitrary angle α. This 'ferromagnetic' state has energy E ferro class. = −3 J N . However, for J > 0 there is another ground state with a small unit cell, 39,40 again with three sublattices where the angles α i between different sublattices differ by multiples of 2 π/3. Also the energy E > class. = −3 J N is invariant under global rotations. The latter state is illustrated by the global structure of Fig. 8 which shows a snapshot of a low-temperature configuration taken during a Monte-Carlo simulation (details again to be given in section V below). Note that the sense of orientation around a triangle, i.e., the chirality of the spins in Fig. 8 is exactly The ferromagnetic state is the simplest case for the computation of fluctuations since M i,j is diagonalized by a Fourier transformation. One finds the modes with the k i defined in Eq. (11). As for the case J < 0, the classical frequencies f ferro class. are proportional to the squares of the SW frequencies ω ferro obtained from a linear Holstein-Primakoff approximation: 40 f ferro class. = ω ferro 2 / 12 S 2 J . By computing the contribution of the modes f ferro class.
to the free energy we find minima for α = n π/3, n = 0, 1, 2, . . ., so that the spins in the ferromagnetic structure lock in to the lattice directions of the triangular lattice. For the lock-in values of α, f ferro class. depends only on one of the k i and has a line of zeros in the perpendicular direction in momentum space.
The three-sublattice state leads to the following 3 × 3 matrix in Fourier space: For α = n π/3, diagonalization of (18) leads to three completely flat branches In particular the lowest branch f > class.,1 = 0 corresponds to a branch of soft modes. In real space these soft modes correspond to the rigid rotation of one single triangle. 40 Note that there is no such flat branch of soft modes for a value of α which is not an integer multiple of π/3.
When computing the contribution of fluctuations around these configurations (α = n π/3) to the free energy, one finds that one third of the modes are quartic instead of quadratic. This yields a free energy of the form: where, again, F zz ∼ N ln β/(2 β) corresponds to the trivial contribution of out-of-plane fluctuations and (compare also Ref. 6) At low temperatures, this term dominates the free energy. The flat branch of soft modes reduces the coefficient of ln β/β from N/2 as in the case of only quadratic modes (compare (14)) to N/3 + N/12 = 5 N/12. This implies two things: firstly, the angles of the 120 • state should lock in at α = n π/3 for low temperatures. Secondly, a thermal order-by-disorder mechanism should favor the 120 • state over the ferromagnetic state for T → 0. As in the case of the ferromagnetic state one finds that the relation f > class.,i = (ω > i ) 2 / 12 S 2 J , where ω > i , i = 1, 2, 3, are the SW frequencies obtained from a linear Holstein-Primakoff approximation, 39,40 holds for arbitrary values of α. Using the results for ω ferro ( k, α) and ω > i ( k, α) to calculate semiclassical ground-state energies E ferro semclass. (α) and E > semclass. (α) in the same manner as in (16) one finds that both are minimal at α = n π/3 and that E > semclass. (n π/3) < E ferro semclass. (n π/3). Thus the semiclassical approach is fully consistent with the classical findings.  hand, CPU time for a similar enumeration of all 6 N such configurations becomes prohibitively big for N > ∼ 15. On the other hand, all known ground states turn out to have mutual angles which are multiples of 2 π/3. Furthermore, we eliminate a global rotational degree of freedom by fixing one arbitrary angle α 0 = π. Then there remain just 3 N −1 configurations with α i ∈ {π/3, π, −π/3} to be enumerated. Direct enumeration of these 3 N −1 configurations can be carried out with reasonable CPU time for N ≤ 27, but becomes quickly impossible for larger N . We have therefore performed such enumerations for N ≤ 27, using the same lattices as in section III.
The number of ground states D N determined in this manner is given in Table I for J > 0. Note that the ordered states which we have described in section IV B are just two of the D N states, but there are many further ground states which can be interpreted as defects in and domain walls between the ordered states. 39,40 Indeed, also closer inspection of the snapshot shown in Fig. 8 reveals the presence of defects in the three-sublattice structure. The 240 = 6 D 12 states described previously 39,40 for N = 12 are recovered by a global rotation of the angles such that α 0 takes on the six values α 0 = n π/3. The last column of Table I lists ln (D N ) /N . The fact that these numbers stay almost constant indicates a finite ground-state entropy per site slightly above 0.3 in the thermodynamic limit.
It is straightforward to derive the N × N matrix M i,j defined in (8) for any ground state and diagonalize it.
Among the eigenmodes f i , one can then identify the g pseudo-Goldstone modes f i = 0 and in turn count the number n g of ground states with g pseudo-Goldstone modes. These numbers are also given in Table I in the form D N = g n g . One observes that all ground states have at least one pseudo-Goldstone mode. There are at most N/3 pseudo-Goldstone modes, and there is only one ground state with this maximal number of pseudo-Goldstone modes corresponding to the three-sublattice state described in section IV B (apart from N = 15; however, this lattice is special in that it has a period three translational symmetry T 3 y = 1). There are N/3 ground states which differ from the per-fect three-sublattice state by a rigid rotation of the spins on certain triangles by an angle 2 π/3, and another N/3 ground states where the spins on a different set of triangles are rotated by −2 π/3. These ground states have two pseudo-Goldstone modes less than the three-sublattice state. The data in Table I show that these 2 N/3 configurations with N/3 − 2 pseudo-Goldstone modes account for all states with the second largest number of pseudo-Goldstone modes for N ≥ 21. This indicates that ground states deviating from the homogeneous three-sublattice state are obtained at the expense of pseudo-Goldstone modes. At finite temperature such "inhomogeneous" ground states are then penalized by an entropic cost because of the reduced number of pseudo-Goldstone modes. This indicates that ground states with bigger deviations from the three-sublattice ground state have a higher free energy for small T since they have less soft modes. Thus, the global minimum of the free energy is the perfectly ordered state with many close-by configurations which deviate only locally from the perfect order. These arguments predict a thermal order-by-disorder selection of the 120 • state among the macroscopic number of ground states for T → 0 and J > 0.
The above enumeration procedure can also be performed for J < 0. In fact, in this case we have carried it out twice, first with α i ∈ {π/3, π, −π/3} and α 0 = π, and then again with all α i shifted by π/6 in order to match the lock-in predicted in section IV A. In sharp contrast to the large degeneracy found for J > 0, we confirm that the ground state is unique (up to global rotations of all angles) for J < 0 and thus identical to the three-sublattice state described in section IV A. Diagonalization of the corresponding N × N matrix M yields one pseudo-Goldstone mode, in complete agreement with the f i shown in Fig. 7 which have just one zero, namely f < class.,1 ( k = 0) = 0.

V. CLASSICAL SYSTEM: MONTE CARLO ANALYSIS
Section IV already contained some discussion of the low-temperature properties in the classical limit. The results of section III lead us to suspect a finite-temperature transition in the quantum system at least for J < 0. Indeed, a finite-temperature phase transition is allowed 44 since the model has only discrete and no continuous symmetries (compare section II). Since such a finitetemperature phase transition should be a classical phase transition, we may hope to gain insight into its universality class by studying the classical counterpart of the model. This motivates us to present results of a Monte-Carlo simulation of the classical model, treating the spins S i as classical O(3) vectors. Simulations were performed on square lattices with diagonal bonds (topologically equivalent to the triangular lattice) and periodic boundary conditions. First, we have used a standard single-spin flip  46 In order to determine error bars, we have used between 100 and 400 independent simulations for J < 0. For J > 0, the standard single-spin flip algorithm turns out to be no longer ergodic for temperatures T < ∼ 0.1 |J|. Such problems are in fact expected in view of the large ground-state degeneracy which we discussed in section IV C. In this region we have therefore used the parallel tempering Monte Carlo method (also known as exchange Monte Carlo -see Refs. 47-50 and references therein). In this framework, n simulations are performed in parallel, each at a different temperature using the standard Metropolis algorithm. Periodically, the exchange between the configurations of two simulations consecutive in temperature is proposed and accepted depending on the energy balance of such a move. A careful choice of the temperature points allows each configuration to shuffle through the entire temperature range during the simulation, greatly reducing the probability of getting stuck in a local minimum of the free energy. In principle, this allows an efficient exploration of the phase space, while not having to wait for rethermalization of the systems after each configuration switch. Strategies to optimize the choice of the temperature grid have been proposed (see, e.g., Ref. 50), but we simply opted for constructing a fine-grained temperature set using the rule of thumb that the probability for a configuration switch to be accepted should always be at least around 70% to 80%. The resulting grid consists of at least 96 points for T /J ∈ [0.01, 0.7], where of course most points lie in the low-temperature region. After an initial thermalization, observables are sampled until convergence of their error bars is observed (the typical duration of a simulation being at least 10 7 Monte-Carlo sweeps per system). We would like to insist that even if parallel tempering is adequate to the task of studying the low-temperature properties of such a highly degenerate frustrated magnet, it is still by no means easy to obtain physically relevant data at such a low temperature for continuous spherical spins as we shall see later.
A. Specific heat Fig. 9 shows results for the specific heat of the classical model for J < 0 (panel (a)) and J > 0 (panel (b)). Statistical errors should be at most on the order of the width of the lines. Although all lattice sizes are bigger than those used previously for the quantum model, there are remarkable similarities of the specific heat of the classical model shown in Fig. 9 with the specific heat of the quantum model, see Fig. 3. For J < 0, a singularity seems to develop in the specific heat for temperatures around T ≈ 1.5 |J|, signaling a phase transition. For J > 0, there is also a broad maximum at 'high' temperatures T ≈ 0.3 J. The finite-size effects for the latter maximum are small indicating that this does not correspond to a phase transition. In this case, the interesting features of the specific heat lie in the low-temperature region, as displayed in Fig. 10(b). As the system size increases, one can see that a small peak builds up in the specific heat for T ≈ 0.02 J. This seems to indicate that a phase transition might occur around that temperature, two orders of magnitudes smaller than for J < 0. We are unfortunately not on a par with the J < 0 data, as the CPU requirements are too steep to secure relevant data for systems larger than 27 × 27 sites even though the specific heat is a comparably robust quantity, and it is clear that other observables are needed to conclude on the existence of this phase transition.
An important difference between the s = 1/2 and the classical model arises at low temperatures: the specific heat of the quantum system has to vanish upon approaching the zero temperature lim T →0 C/N = 0, while due to the remaining continuous degrees of freedom, the specific heat of the classical system approaches a finite value for  T → 0. For J < 0, the equipartition theorem predicts an N/2 contribution to the specific heat for each transverse degree of freedom which yields lim T →0 C/N = 1, in excellent agreement with the results depicted in the panel (a) of Fig. 9. For J > 0, one must take into account the fact that the flat soft-mode branch of the threesublattice state is expected to contribute only N/12 to the specific heat. Thus for J > 0 one should expect lim T →0 C/N = 11/12 = 0.916666.... As can be seen in Fig. 10(b), we observe a specific heat lower than one in the low-temperature region along with a downward trend as T goes to 0 for all the system sizes studied. However, according to the data which we have at our disposal, it seems that one would have to go to very low temperatures T < 10 −3 J in order to verify the prediction for the zero-temperature limit.
Returning to the finite-temperature transition for J < 0, Fig. 10(a) shows a zoom into the relevant temperature range, including data for up to N = 90 × 90 spins. At these bigger system sizes, the position of the maximum continues to shift to lower temperatures and the 0 π/3 2 π/3 π 4 π/3 5 π/3 2 π φ i  Fig. 10(a) demonstrate that the maximum value of the specific heat starts to decrease as one goes to system sizes beyond N = 45 × 45. This implies that the exponent α which characterizes the divergence of the specific heat at the critical temperature is very small or maybe even negative.
B. Sublattice order parameter, Binder cumulant, and transition temperature for J < 0 According to subsection IV A, we expect that the phase transition observed for J < 0 is a transition into a threesublattice ordered state. This ordering is indeed exhibited at least at a qualitative level by snapshots of Monte-Carlo simulations at low temperatures (compare Fig. 5). In addition, one observes in Fig. 5 that the spins are lying essentially in the x − y-plane for low temperatures.
Furthermore, we expect a lock-in of the spins to one of 6 symmetrically distributed directions in the plane at low temperatures (compare Fig. 6). The latter prediction is indeed verified by the histogram of the angles of the in-plane component of the spins φ i at low temperatures shown in Fig. 11. Note that the histogram is rather flat for the smaller lattices (in particular the N = 6 × 6 lattice) and sharpens noticeably as the lattice size increase to N = 36 × 36 (the largest lattice which we have considered in this context). The fact that the lock-in occurs only on large lattices can be attributed to the replacement of the sum over k in (15) by an integral being a good approximation only for large lattices.
To test for the expected three-sublattice order, we in- troduce the sublattice order parameter where the sum runs over one of the three sublattices L of the triangular lattice. Fig. 12 shows the behavior of the square of this sublattice order parameter for J < 0. One observes that the sublattice order parameter indeed increases for T < 2 |J| and goes indeed to m 2 s = 1 for T → 0, as expected for a three-sublattice ordered state. Inclusion of larger lattices (up to N = 90×90) allows one to restrict the ordered phase to T < ∼ 1.7 |J|. However, more accurate estimates for the transition temperature can be obtained in a different manner.
A useful quantity to determine the transition into an ordered state accurately is the 'Binder' cumulant 45,51,52 associated to the order parameter (23) via We have chosen the prefactor in (24)  such that U s = 2/5 in this case. Hence, for an ordered state we expect U s ≈ 0.4 for T < T c . Fig. 13 shows results from classical Monte-Carlo simulations for the Binder cumulant U s of the system with J < 0. First, the broad temperature range shown in Fig. 13(a) confirms that indeed U s ≈ 0.4 in the ordered low-temperature phase and U s ≈ 0 for high temperatures. The transition temperature can now be accurately extracted from the crossings of the Binder cumulants at different sizes N . 45,51,52 For this purpose, Fig. 13(b) zooms in to the relevant temperature range, including bigger system sizes N . Although the crossings between any pair of system sizes N 1 and N 2 fall into a narrow temperature window, there still remains a small residual dependence on the sizes N 1 and N 2 considered. In order to perform an extrapolation N → ∞, we have analyzed the crossings between neighboring system sizes N 2 > N 1 ≥ 9 × 9. This leads to the estimate for the thermodynamic limit N → ∞.
C. Nature of the phase transition for J < 0 Having determined the transition temperature for J < 0, one would like to clarify the universality class of the phase transition. On the one hand, there is no evidence for any latent heat in the specific heat at T < c , see Fig. 10(a), i.e., the ordering transition appears to be continuous for J < 0. On the other hand, a negative dip in the Binder cumulant for T > T c , as observed in Fig. 13(a) is sometimes taken as evidence for a first-order transition (see, e.g., Ref. 53). In order to distinguish better between the two scenarios we use histograms of the energy E of the microstates realized in the Monte-Carlo procedure. [54][55][56] We have collected such histograms for several system sizes and temperatures. Fig. 14 shows two representative cases on the N = 90 × 90 lattice, namely T = 1.7025 |J| which corresponds to the maximum of the specific heat for the 90 × 90 lattice (compare Fig. 10(a)) and T = 1.566 |J|, the estimated critical temperature of the infinite system, see Eq. (25). We always find bell-shaped almost Gaussian distributions, which are characteristic for continuous transition. We never observed any signatures of a splitting of this single peak into two, as would be expected for a first-order transition. [54][55][56] Hence, the transition appears to be continuous and we will now try to characterize its universality class further in terms of critical exponents.
We start with the correlation length exponent ν which can be extracted from the finite-size behavior of the Binder cumulant: close to T c , the Binder cumulant should scale with the linear size of the system L as 45,51,52  The slope of the Binder cumulant yields the correlation length exponent ν (top panel), the sublattice magnetization ms yields the exponent β, and the specific heat C yields the exponent α. Lines show the fits which have been used to estimate the exponents. Note that the scale is doublelogarithmic in all three panels. Fig. 13(b) shows that there is very little curvature in the Binder cumulants U s as a function of temperature T close to the estimated critical temperature (25). Therefore d U s /dT can be extracted without much sensitivity to the error of (25). The result is shown by the symbols in the top panel of Fig. 15. Now we can determine ν by fitting this to (26). Since inclusion of the correction term renders the fit unstable, we use only the leading term (i.e., we set b = 0 in (26)). A fit (which is shown by the line in the top panel of Fig. 15) then yields We now turn to the order parameter exponent β which can be extracted from the finite-size behavior of the sublattice order parameter m s . The sublattice order parameter should have a scaling behavior (see for example Refs. 45 and 57) Specialization of (28) to T = T c yields The middle panel of Fig. 15 shows the Monte-Carlo results for m 2 s at three temperatures which cover the estimate (25) for T < c and its error bars. The fits of these results to (28) which are shown by the lines in the middle panel of Fig. 15 lead to Finally, we turn to the specific heat exponent α. We proceed in the same manner as for the order parameter exponent β and make again a scaling ansatz for the specific heat: 45,57 The lower panel of Fig. 15 shows the specific heat results at the estimate (25) for T < c . One observes that this does not follow a power law very well. Indeed, it is known that non-scaling contributions to the specific heat can be important. 57 However, including a constant in the ansatz (31) does not lead to a stable fit. We therefore fit only the data for L = 27, 36, 45, and 90 (lines in the lower panel of Fig. 15). This procedure leads to the estimate Note that the error bar includes just the error of the fit. In view of the deviations from a simple power law, this is probably too optimistic. Indeed, α could very well be (slightly) negative, as is suggested by the fact that the maximal value of C in Fig. 10(a) decreases when the system size increases from N = 45 × 45 to 90 × 90. Even if the error bars in (25), (27), and (32) should be too optimistic, it remains safe to conclude that we find a rather larger correlation length exponent ν > ∼ 3. It should be noted that in combination with a specific heat exponent α ≈ 0, we then find that the hyperscaling relation 45 (with the spatial dimension d = 2) is strongly violated. On the other hand, we could use the relation (33) to estimate α, in particular if we expect it to be negative (compare, e.g., Ref. 58 for a similar situation). Insertion of (27) into (33) yields a very negative exponent α/ν = −1.52 ± 0.04. Again, this reflects the large exponent ν. In fact, a large correlation length exponent ν ≈ 4 has been found in other two-dimensional disordered systems. 59,60 However, in those cases the large exponent corresponds to approaching the critical point via a fine-tuned direction in a two-dimensional parameter space and there is a second, substantially smaller correlation length exponent. 59,60 Thus, we are left with not completely unreasonable, but definitely highly unusual values of the critical exponents ν and α.

D. Critical temperature
As mentioned earlier, for J > 0, the specific heat alone is only mildly conclusive regarding the existence of a low-temperature phase transition to a three-sublattice ordered state. This statement requires to be supported by the analysis of other observables. The sublattice magnetization (23) will tell us whether significant order is developing in the low-temperature region or not. Our results shown in Fig. 16 show that three-sublattice order is indeed developing, although an appreciable order develops only at temperatures that are so low that they become increasingly difficult to access with increasing system size. To take a closer look at the low-temperature ordered state, we took some snapshots of the system during the simulation for N = 12×12 spins. A typical configuration is reproduced in Fig. 8. While the global structure corresponds indeed to a 120 • three-sublattice ordered state, we also observe the presence of defects. The presence of these defects is neither a surprise nor in contradiction with the existence of the phase transition, as they are in fact necessary ingredients of the order-by-disorder mechanism. Note also that the spins in Fig. 8 lie essentially in the x−y-plane and -up to small fluctuations-are aligned with the lattice.
As for J > 0, the Binder cumulant (24) allows us both to further support our conclusions concerning the lowtemperature ordered state and to obtain an estimate of T > c . Fig. 17 shows that all the curves for the different system sizes cross in a region around T ≈ 0.015J. First, this is a strong argument in favor of the existence of the ordering transition. We used the smallest three system sizes (N = 6 × 6, 9 × 9, and 12 × 12) for which we have the best statistics to obtain an estimate for the transition temperature: The error bars on the data are unfortunately too large to get precise values for the critical exponents and thus prevent us from investigating the nature and the universality class of the transition. However, the fact that the low-temperature ordered state breaks the same symmetries irrespectively of the sign of J suggests that the universality class for J > 0 is the same as for J < 0.

VI. CONCLUSIONS
Although the Hamiltonian (1) may seem unusual in the context of frustrated magnetism, it is instructive in many respects and illustrates the rich phenomenology often present in this subject. Either seen as the strong dimerization limit of the kagome lattice of spin 1/2 in a magnetic field, or as a possible illustration of an orbital model, 29-37 the underlying physics associated to this Hamiltonian is extremely interesting for both signs of the coupling constant J.
In the quantum case (for spin 1/2) we show, by studying the low-energy spectra using the Lanczos method, that a thermodynamic gap of the order of 3 |J| is present for J < 0, while for J > 0 the gap, if present, would be at most of the order of 0.02 J. The six-fold degeneracy of the would-be ordered ground state which is predicted by semiclassical considerations is not observed in our numerical results, probably due to the small lattices considered. These results illustrate very well how the deep quantum (S = 1/2) regime differs from the large S spinwave predictions. The specific heat curves point to a phase transition around T ≈ |J| for J < 0, while a lower temperature peak shows up in the positive J case. This last peak could be due to an ordering phase transition at a very low temperature T c ≤ J/100. In both cases one is tempted to envisage a finite-temperature phase transition whose nature could be understood by the analysis of the classical model.
The analysis of the classical model has revealed to be also quite interesting and instructive. For J < 0 the lowest-energy configuration consists in an in-plane antiferromagnetic arrangement of the spins with given chirality accompanied by a 'spurious' continuous rotational degeneracy which does not correspond to any symmetry of the Hamiltonian. This pseudo degeneracy is lifted by entropy at finite temperature giving rise to an ordering at low temperature as observed by Monte Carlo data which locate the transition temperature at T c /|J| ∼ 1.566 ± 0.005. Inspection of the histograms of the energy close to the transition temperature gave no evidence of a first-order transition. Hence, we analyzed it within the scenario of a continuous transition and estimated unusual values for the critical exponents α and ν, strongly violating the hyperscaling relation. It should nevertheless be mentioned that we cannot exclude the existence of a crossover scale which exceeds the lattice sizes accessible to us. The fact that lock-in of the spin components to the lattice requires a certain length scale may point in this direction. An unambiguous determination of the universality class of the transition would require improved methods. A first possibility is to restrict the degrees of freedom to the in-plane configurations 42 which are realized in the low-temperature limit. Even more efficiency could be gained by additionally restricting each spin variable to the 6 spin directions which are stabilized in the zero-temperature limit. However, it remains to be investigated whether the second modification changes the universality class of the transition.
For J > 0 the situation is even more interesting. Although a 'spurious' rotational degeneracy is also present for the antiferromagnetic 120 • configuration (with the opposite chirality than the one for J < 0), the manifold of lowest-energy configurations is more complex. There exist local discrete 'flips' of triangles which bring one from the homogeneous antiferromagnetic lowest-energy configuration to another configuration with the same energy.
The mechanism that gives rise to ordering is again understood by analyzing the entropic spectra over each of these configurations. The homogeneous antiferromagnetic configuration has a whole branch of soft modes in its classical spin-wave spectrum. Flipping one triangle to jump to another lowest-energy configuration also destroys one soft mode. One is then left with a scenario that can be understood with an Ising-type low-temperature expansion picture of the system, where each 'flippable' triangle plays the rôle of an Ising spin. The difference is in the fact that flipping one spin on an otherwise perfectly ordered background costs no energy but an entropy, or if one wishes a temperature-dependent pseudo energy. An ordering transition will also take place, as in a normal energetic system, but at a much smaller temperature. This transition temperature is observed in the Monte-Carlo analysis to be at around T c /|J| ≈ 0.0125, two orders of magnitude smaller than in the J < 0 case.
For J > 0, one may also wonder how the 120 • ordered state of the model (1) relates to the structure of the magnetization 1/3 state of the homogeneous kagome lattice. 7,25 In order to address this question, we need to associate a variational wave function to the classical 120 • ordered state which is indicated in Fig. 8. First, we associate quantum wave functions to the three classical spin directions as in Fig. 18(a). The phase factors are chosen in order to yield a convenient representation in terms of spin configurations of a triangle after insertion of Fig. 1(b) for the chirality pseudo spins. Insertion into the 120 • wavefunction for a triangle of the triangular lattice on the left side of Fig. 18(b) then yields the expression in terms of the 8 spin configuration of a ninesite unit of the underlying kagome lattice shown on the right side of Fig. 18(b). Note that the two terms on the first line of the right side of Fig. 18(b) amount exactly to the variational wave function for the magnetization 1/3 state of the homogeneous kagome lattice, 7,25 as it follows from a mapping to a quantum-dimer model on the honeycomb lattice. 24 Thus, the present results for the strongly trimerized kagome lattice may be smoothly connected to the 1/3 plateau state of the homogeneous kagome lattice.
To conclude, the model (1) has turned out to be a very interesting laboratory to understand the emergence of a hierarchy of energy scales originating from different levels of order by disorder. The emergence of such a hierarchy in a classical model is related to similar hierarchies in quantum systems as for example the huge difference between the magnetic and non-magnetic gaps in the kagome spin 1/2 system at magnetization 1/3. 7 Moreover, if the transitions observed in this work can be confirmed to be continuous, the exponents will probably correspond to exotic models, like parafermionic conformal field theories. 61

ACKNOWLEDGMENTS
We are grateful to P. C. W. Holdsworth, H. G. Katzgraber, F. Mila, and M. E. Zhitomirksy for useful discussions, comments and suggestions. This work has been supported in part by the European Science Foundation through the Highly Frustrated Magnetism network. D.C.C. is partially supported by CONICET (PIP 1691) and ANPCyT (PICT 1426). Furthermore, A.H. is supported by the Deutsche Forschungsgemeinschaft through a Heisenberg fellowship (Project HO 2325/4-2).