The cubic Kugel–Khomskii model for triply degenerate t2g electrons

The cubic Kugel–Khomskii Hamiltonian for titanates describes spin and orbital superexchange interactions between d1 ions in an ideal cubic perovskite structure in which the three t2g orbitals are degenerate in energy, and electron hopping in the presence of large Coulomb interactions is constrained by cubic site symmetry. We review results for the unusual symmetry of this model and give a simple physical argument that explains why this symmetry prevents long-range spin order at non-zero temperatures. We also review the Landau theory of the disordered phase of this model, which gives rise to susceptibilities that are dispersionless along one wavevector axis. We present new results for the mean-field equations, which describe possible long-range order (in the presence of suitable stabilizing perturbations). We also analyse the role of thermal and quantum fluctuations and for the first time give a renormalization group analysis of this model in d spatial dimensions. Finally, we briefly review extensions of this model which are needed to describe real systems, such as LaTiO3.


Introduction
High-temperature superconductivity [1] and colossal magnetoresistance [2] have sparked much recent interest in the magnetic properties of transition metal oxides, particularly those with orbital degeneracy [3,4]. In many transition metal oxides, the d electrons are localized due to the very large on-site Coulomb interaction, U. In cubic oxide perovskites, the crystal field of the surrounding oxygen octahedra splits the d-orbitals into a twofold degenerate e g and a threefold degenerate t 2g manifold. In most cases, these degeneracies are further lifted by a co operative Jahn-Teller (JT) distortion [3], and the low-energy physics is well described by an effective superexchange spin-only model [5]- [7]. However, some cubic perovskites, such as the titanates (RTiO 3 , where R = La, Y, etc) [8,9], do not undergo a significant JT distortion, in spite of the orbital degeneracy [10], but the degeneracy is lifted by a rotation of the oxygen octahedra [11]. Although it seems likely that it will not be possible to experimentally realize the ideal cubic model, it is interesting to study it, for two main reasons. Firstly, this model is very interesting theoretically, and the concepts needed to understand it may apply elsewhere. Secondly, some real systems may be close enough to being cubic, so that their behaviour may be obtained perturbatively from this limit.
In systems such as the ideal cubic model, the effective superexchange model must deal with not only the spin degrees of freedom but also the degenerate orbital degrees of freedom [3,4,12]. The large degeneracy of the resulting ground states may then yield rich phase diagrams, with exotic types of order, involving a strong interplay between the spin and orbital sectors [4,8,9]. In the idealized cubic model for the titanates, there is one d electron in the t 2g degenerate manifold. Following Kugel and Khomskii (KK) [12], one starts from a Hubbard model with onsite Coulomb energy U and nearest-neighbour (nn) hopping energy t. For large U, this model can be reduced to an effective superexchange model, which involves only nn spin and orbital coupling, with energies of the order of = t 2 /U. This low-energy model has been the basis for some theoretical studies of the titanates. In particular, it has been suggested [13] that the cubic KK Hamiltonian gives rise to an ordered isotropic spin phase, and that an energy gap in the spin excitations can be caused by spin-orbit interactions [14]. Recently [15,16] (this will be referred to as I) we have presented rigorous symmetry arguments, which show several unusual properties of the cubic KK Hamiltonian. Perhaps the most striking symmetry is that with respect to rotation of the total spin of α orbitals (where α denotes one of the three t 2g orbitals) summed over all sites in a plane perpendicular to the α-axis. As shown in I, this symmetry implies that the system does not support long-range spin order at any non-zero temperature. In addition, this symmetry implies that in the disordered phase the wavevector-dependent spin susceptibility for α orbitals, χ α (q), is dispersionless in the q α -direction [17] (this will be referred to as II). This peculiar rotational symmetry can be broken by almost any perturbation. In II, we considered the effect of symmetry-breaking perturbations due to (a) spin-orbit interactions, (b) next-nearestneighbour (nnn) hopping and (c) Hund's rule coupling. According to the general symmetry argument of I, although long-range order at non-zero temperature is possible when spin-orbit interactions are included, the system still possesses enough rotation symmetry that the excitation spectrum should be gapless. This argument would imply that mean-field theory will produce a state which has a continuous degeneracy associated with global rotation of the spins even in the presence of spin-orbit interactions. Here we will review the mean-field theory introduced in II and interpret the results obtained therefrom in light of the general symmetry arguments. We extend the analysis to study the self-consistent equations of mean-field theory. The symmetry of these solutions confirms the results of I and II, but has some additional symmetry (specific to mean-field theory) which we here show is removed by fluctuations. In addition we review results for the cubic KK Hamiltonian when the symmetry breaking perturbations mentioned are included. In the presence of spin-orbit interactions, the staggered moments of different orbital states are not collinear, so that the net spin moment is greatly reduced from its spin-only value. The effect of including only nnn hopping is also interesting. Within the mean-field theory, this perturbation was found to stabilize a state having long-range staggered spin order for each orbital state, but the staggered spins of the three orbital states add up to zero. When only Hund's rule coupling is included, mean-field theory predicts stabilization of long-range spin and orbital order. In this paper, we go beyond mean-field theory to treat quantum and thermal fluctuations. We show that these fluctuations eliminate the aforementioned additional symmetry present in mean-field theory, since they favour spin-only order. As a result, a state with long-range order of both spin and orbital degrees of freedom can only occur when the strength of the Hund's rule coupling exceeds some critical value which we cannot estimate within the present formalism. Finally, we give a renormalization group (RG) analysis of the model when it is generalized to d spatial dimensions. In analogy with the lower critical dimension d < (at and below which long-range order is not permitted) being raised from the usual d < = 2 to d < = 3, the RG analysis shows that the upper critical dimension d > (above which fluctuations away from mean-field theory can be ignored) is raised from the usual d > = 4 to d > = 5.
Briefly this paper is organized as follows. In section 2 we discuss the KK Hamiltonian and give a hand-waving discussion of its unusual invariance with respect to rotation of the spin associated with a given orbital. The consequences of this, such as the impossibility (shown in I) of long-range spin order at any non-zero temperature, are discussed. In section 3 we review results obtained in II for the high-temperature phase. The wavevector-dependent spin susceptibilities that diverge as the temperature is lowered through a critical value have dispersionless directions, so that, unusually, mean-field theory provides no 'wavevector selection' at the mean-field transition. We review there the effect of several symmetry breaking perturbations, namely spin-orbit interactions, nnn hopping and Hund's rule coupling. In each of these cases 'wavevector selection' leads to the usual two-sublattice structure, but the qualitative nature of ordering depends on which perturbation is considered. In section 4 we consider the mean-field description of the ordered phase (assumed to be stabilized by some symmetry-breaking perturbation). The self-consistent solutions for the order parameters reflect the unusual rotational symmetry found in I. Here we discuss the additional symmetry introduced by the approximations of mean-field theory and in section 5 we show that fluctuations remove this spurious symmetry. In section 6 we give an RG analysis of this model when it is generalized to d spatial dimensions. We give specific results in the limit of 5 − dimensions. In section 7 we discuss the significant modifications in the model which are necessary to treat a real system, namely, LaTiO 3 . Finally, our conclusions are summarized in section 8.

The Kugel-Khomskii Hamiltonian and its novel symmetries
Following the seminal work of Kugel and Khomskii [12] (KK), we describe our system by a Hubbard Hamiltonian H H of the form where c † iασ creates an electron in the orbital labelled α in spin state σ on site i and N α (i) = σ c † iασ c iασ ,˜ α ≡ α − (U αα /2), where α is the crystal-field energy of the α orbital, t αβ (i, j) (which we assume to be real) is the matrix element for hopping between orbital α of site i and orbital β of site j and U αβ is the direct Coulomb interaction between electrons in orbitals α and β on the same site. Note that this Hamiltonian does not include the somewhat smaller Coulomb exchange (Hund's rule) terms [18,19]. We assume that the crystal field splits the five orbital d states into three low-energy t 2g states, d yz ≡ |x , d xz ≡ |y and d xy ≡ |z and that the two other e g states have high enough energy that we can neglect their presence. We confine our attention to the wide class of materials in which hopping between magnetic ions is mediated by intervening oxygen ions. In that case, following KK, we note that the hopping matrix t αβ between nn's is diagonal in orbital indices and also t αα = 0 if the nn bond lies along the α-axis, which has been called [18] the 'inactive' axis for α hopping. Note that these symmetries only hold for tetragonal or cubic site symmetry. They are broken by a rotation of the oxygen octahedra, as in LaTiO 3 .
When t U, KK reduced equation (1) at lowest order in t/U to an effective Hamiltonian for the manifold of states which remain when U → ∞ and the crystal-field splitting is assumed to be very large. This is the KK Hamiltonian, where = 2t 2 /U and ij ∈ α means that the sum over pairs of nearest neighbours is restricted to those for which r ij is along an α-axis. This Hamiltonian can be regarded as a many-band generalization of the Heisenberg Hamiltonian, which is obtained from the single-band Hubbard model with one electron per site, for which the exchange constant is J = 4t 2 /U. To make the analogy more apparent, the KK Hamiltonian is often written in terms of spin variables as where S i ≡ 1 2 σ i is the vector spin operator for an electron in site i, ij indicates that the sum is over pairs of nearest neighbour ions and the exchange 'constant' is now an orbital operator written in terms of spinless fermion operators as where a † iα creates an electron in orbital |α on site i (of either spin) and α = ij means that we sum α over the values x, y and z, except for the value of α which corresponds to the co-ordinate direction of the bond ij , which is the inactive axis for orbital α. The mathematical relation between the spinless fermion operators a † iα and the true electron creation operators c † iασ is where the superscript on the vector operator σ labels its Cartesian components.
To understand the symmetry of the KK Hamiltonian with respect to spin rotation, it is useful to review a symmetry which has been recognized for some time in the literature. Consider the hopping of electrons which are in the orbital state |x ∼ yz. As noted, such electrons cannot hop along the inactive axis, which in this case is the x-direction. Therefore, it is easy to see that the total number of x electrons in any specified plane perpendicular to the x-axis (we will call such a plane an x-plane) is constant [19]. The inactive axis prevents these electrons from hopping from one x-plane to another x-plane. A more subtle version of this argument indicates that the total quantum spin vector S summed over all electrons of orbital α in an α plane is also a constant. This conclusion is implicit in the result that the band structure of a non-magnetic system in the absence of a magnetic field does not depend on the choice of quantization axes for the electron spin. In the present context, the Hamiltonian will be invariant under spin rotation provided we rotate the spin at all sites, which are connected either directly or indirectly by hopping. In particular, since hopping of α electrons can only take place within a single α-plane, we see that if we rotate the spin for all α-orbital electrons in some arbitrarily selected α plane, the Hamiltonian will be invariant.
A direct consequence of this unusual rotational invariance is that long-range order of the spin of α-orbital electrons, if it exists, must be confined to a single α-plane. But since an isotropic model in two dimensions cannot support long-range order at any non-zero temperature, we expect that this model does not exhibit long-range spin order at non-zero temperature. This hand-waving argument gives an intuitive meaning to the rigorous mathematical argument given in I, which is based on the Mermin-Wagner theorem [20], and which allows us to conclude that the cubic KK model does not support long-range spin order at any finite non-zero temperature.
A striking application of this symmetry is found when one determines the ground state of a system of eight sites placed on the vertices of a cube. Since each site has three orbital and two spin states, the Hilbert space for one electron per site for such a system has 6 8 ≈ 1.7 × 10 6 states. Using a sparse matrix package, we determined (in I) that the ground state has a energy of −6.2716 . In this ground state, the total number of electrons in orbital, say x, is a good quantum number. It was found that the ground state had four electrons of one orbital, four electrons of a second orbital and no electrons of the third orbital. The ground state was threefold degenerate corresponding to the three choices for the orbital which had no electrons. Given this information we may use the symmetries described above to simplify the determination of the ground state wavefunction and energy. The energy spectrum and the corresponding wavefunctions can be found by diagonalizing a much smaller matrix if one works within a manifold defined by the conserved numbers applied to each face of the cube (see I). We actually have 18 conserved quantum numbers and there is an astonishing numerical simplification when maximal use is made of these symmetries, in which case we deal with a matrix of dimensionality 16!

Landau expansion for T > T c
The Landau expansion of the free energy is a multivariable expansion of the free energy in powers of the full set of possible order parameters. Our derivation is based on the variational principle according to which the exact free energy is obtained by minimizing the free-energy functional F(ρ) as a function of the trial density matrix ρ, The first term is the trial energy and the second is −T times the trial entropy, where T is the temperature. Mean-field theory is obtained by the ansatz that ρ is the product of single-site density matrices, ρ(i): and F(ρ) is then minimized with respect to the variables used to parametrize the density matrix, ρ(i). Since ρ(i) acts in the space of t 2g states of one electron, it is a 6 × 6-dimensional Hermitian matrix with unit trace. Its most general form is where σ is the Pauli matrix vector and A αβ (i), B x αβ (i), B y αβ (i) and B z αβ (i) are 3 × 3 Hermitian matrices, of which the first is traceless. The diagonal terms of the matrix A are parametrized for later convenience as The thermal average of any operator O(i) is given by tracing over the six states |α, σ , Thus, the thermally averaged occupation number for α electrons at site i is given by and hence the parameters a 1 (i) and a 2 (i) describe unequal population of the orbitals at site i. They thus represent orbital-order parameters, which may become non-zero at low temperatures.
As for the spin-order parameters, the µ-components of the thermally averaged spin vector, s α (i), for α electrons at site i, is Using Fourier representation with the definition where N is the total number of lattice sites, we find that at second order where there are just two distinct types of susceptibility: a diagonal one, given by where c α ≡ cos(q α a) and an off diagonal one, given by At high temperatures all the eigenvalues of the inverse susceptibility matrix are finite and positive. As the temperature is reduced, one or more eigenvalues may become zero, corresponding to an infinite susceptibility. Usually this instability will occur at the star 5 of some wavevectors, and this set of wavevectors describes the periodicity of the ordered phase near the ordering transition. This phenomenon is referred to as 'wavevector selection'. The eigenvector associated with the divergent susceptibility contains information on the qualitative nature of the ordering. Here, a central question which the eigenvector addresses is whether the ordering is in the spin sector, the orbital sector or both sectors. If the unstable eigenvector is degenerate, one can usually determine the symmetries which give rise to Goldstone (gapless) excitations. In the present case, we see that the instabilities first appear at kT = 2 /3 for the diagonal susceptibilities. A detailed analysis of the susceptibilities for unequal occupations of the three orbital states (see II) shows that it occurs only when the wavevector q assumes its antiferromagnetic value Q = (π, π, π)/a, which leads to a two sublattice structure called the 'G' state. The twofold degeneracy is the symmetry associated with rotations in occupation number space N x , N y and N z with the constraint that the sum of these occupation numbers is unity. In contrast, the inverse spin susceptibility χ −1 αα of equation (15) has a flat branch so that it vanishes for kT = 2 /3 for any value of q α , when the two other components of q assume the antiferromagnetic value π/a. This wavevector dependence indicates that correlations in the spin susceptibility of α-orbital electrons become long ranged in an α-plane, but different α-planes are completely uncorrelated. Note that beyond the fact that there is no wavevector selection in the spin susceptibility, one has complete rotational invariance in B γ αα (q) for the components labelled by γ independently for each orbital labelled α. This result reflects the exact symmetry of the Hamiltonian with respect to rotation of the total spin in the α orbital summed over all spins in any single α-plane [15]. If we restrict attention to the G wavevector q = Q, we have complete rotational degeneracy in the 11-dimensional space consisting of the nine B γ αα (Q) spin order parameters and the two a n (Q) occupational order parameters, leading at this level of approximation to O(11) symmetry! Most of this symmetry only holds at quadratic order in mean-field theory, and we expect that higher-order terms in the Landau expansion will generate anisotropies in this 11-dimensional space to lower the symmetry to the actual cubic symmetry of the system. Note that this O(11) symmetry includes the possibility of rotations which mix spin and orbital degrees of freedom. However, such a mixing is not a true symmetry of the system. As we shall see, the anisotropy that inhibits this mixing of spin and orbital degrees of freedom is not generated by the quartic terms in the free energy. Perhaps unexpectedly, as we show in section 5, this anisotropy is only generated by fluctuations not accessible to mean-field theory.
Usually, a dispersionless susceptibility is an artefact of mean-field theory and does not represent a true symmetry of the full Hamiltonian. In such a case, the continuous degeneracy is lifted by fluctuations, which can either be thermal fluctuations [21] or quantum fluctuations [22]. Here we have a rather unusual case in that the spin susceptibility has a dispersionless direction, resulting from an exact true symmetry of the quantum Hamiltonian that persists even in the presence of thermal and quantum fluctuations.
In view of the Mermin-Wagner argument given in I, it is more sensible to apply mean-field theory to a KK model having small perturbations which lead to long-range ordering at non-zero temperatures, as has been done in II and is briefly reviewed below. In all cases we have studied so far, and probably for most of the realistic cases, the G wavevector is selected by perturbations quadratic in the order parameters. Thus, we expect that in the asymptotic critical regime, we only need to study the set of variables at the G wavevector. For that purpose the model we need to consider will have a quadratic term perturbed by the specific perturbation (spin-orbit, nnn coupling, etc) under consideration with a fourth-order term calculated for the cubic KK model without any such perturbation.
As shown in II, terms which are quartic in the critical variables [of wavevector Q ≡ (π, π, π)/a] come from two sources. Firstly, there are cubic terms which involve two critical variables and one non-critical variable. Their contribution is obtained upon minimizing with respect to the non-critical variables (at zero wavevector). The second contribution comes from the bare quartic terms. The resulting total fourth-order free energy is This quartic term reduces the O(11) symmetry. However, since this term (and higher order ones as well) depends only on s 2 α , we have symmetry with respect to rotation of s α (Q) for each value of α. When s α (Q) is zero, we still have rotational invariance in the orbital sector since the free energy then only depends on (a 2 1 + a 2 2 ). The effect of cubic symmetry in the purely orbital sector is only felt in sixth order in the variables a i . Unexpectedly, perhaps, the quartic term also allows a rotation between the spin and orbital sectors. We will see this more clearly when we consider the self-consistent mean-field equations for T < T c .

Effect of perturbations
The idealized KK Hamiltonian has sufficient symmetry, so that there is no wavevector selection 6 within the mean-field theory and the exact symmetry of this model does not support long-range order at non-zero temperatures [15,16]. However, various perturbations which are inevitably present, even when there is no distortion from perfect cubic symmetry, may lead to a wavevector selection. Here we examine this possibility by studying the effect of several perturbations.

Spin-orbit interactions.
These interactions destroy the peculiar invariance with respect to rotating planes of spins of different orbitals independently, and lead to a wavevector selection from the susceptibility. A plausible guess is that the system will select the wavevector Q to allow simultaneous condensation of spins of all three orbitals. This is indeed confirmed by the detailed calculations in [17]. However, although the spin-orbit interaction selects the directions for the spin vectors s α of orbital α relative to one another, there is rotational invariance when all the s α 's are rotated together. This indicates that relative to the mean-field state there are zero frequency excitations, which correspond to rotations of the staggered spin. More generally, one can establish this rotational invariance to all orders in the spin-orbit coupling λ and without assuming the validity of mean-field theory [6,15]. In addition, the spin state induced by spinorbit coupling does not have the spins of the individual orbitals, s α , parallel to one another and thus the net spin, S, is greatly reduced by this effect. We find that the total spin is 1/ √ 3 of what it would be if s α were parallel to one another.

Further neighbour hopping.
For a perfectly cubic system, the nnn hopping process comes from the next-to-shortest exchange path between magnetic ions. It is shown in II that this perturbation to the quadratic terms in the Landau free energy selects the wavevector Q, and yields the transition temperature kT c = 2( + )/3, where = (t ) 2 /U and t is the effective hopping matrix element connecting next-nearest neighbours. However, the fourth-order term does not select a particular direction for order. Rather, the three orbital spin vectors make a 120 • angle with each other and lie in a single plane. Hence, the net staggered moment vanishes. There is long-range spin order, but not of any simple type. Such a magnetic structure, for which the local moment (summed over all orbits) vanishes, will be rather difficult to detect experimentally.
It is instructive to argue for the above results without actually performing the detailed calculations of [17]. We expect the effect of indirect exchange between nnn's to induce an antiferromagnetic interaction between the spins of different orbitals of nnn's. Note that the wavevector Q describes a two sublattice structure in which nnn's are on the same sublattice. Accordingly, as far as mean-field theory is concerned, an nnn interaction between different orbitals is equivalent to an antiferromagnetic interaction between spins of different orbitals on the same site. So the spins of the three orbitals form the same structure as a triangular lattice antiferromagnet [23], namely the spins of the three different orbitals are equal in magnitude and all lie in a single plane with orientations 120 • apart. This state still has global rotational invariance, but also, as does the triangular lattice antiferromagnet, it has degeneracy with respect to rotation of the spins of two orbitals about the axis of the spin of the third one. Thus we expect (and in I it is proved) that there should be gapless excitations corresponding to this rotational invariance.

Hund's rule coupling.
As is shown in [17], the Hund's rule coupling favours antiferromagnetic orbital ordering, as described by the order parameters a 1 (Q) and a 2 (Q). Since the mean-field temperature for spin and orbital ordering were degenerate in the absence of the Hund's rule coupling (scaled by η), we conclude that within mean-field theory the addition of an infinitesimal η gives rise to an ordering transition in which the ordered state shows longrange antiferromagnetic orbital order. However, as we show in section 5 for the bare KK model, fluctuations stabilize the spin-only states relative to orbital states. We thus conclude that when fluctuations are taken into account, it will take a finite amount of Hund's rule coupling to bring about orbital ordering. For spin ordering, the mean-field state is degenerate with respect to an arbitrary rotation.
It is interesting to examine the anisotropy in the mean-field solution for orbital order. Wavevector conservation dictates that we can have only products involving an even number of the variables a 1 (Q) and a 2 (Q). If we write a 1 (Q) = a cos θ Q and a 2 (Q) = a sin θ Q , then one can show [17] that the contribution to the free energy of the order of a 4 is independent of θ Q , but the term of the order of a 6 is of the form δF = a 6 [C 0 + C 6 cos(6θ Q + φ)]. This form indicates an anisotropy, so that the mean-field solution is not subject to a rotational degeneracy in a 1 − a 2 space. If C 6 is positive and φ = 0, these minima come from the six angles that are equivalent to θ Q = π/2 + nπ/3. For θ Q = π/2, a 1 = 0 and we have ordering involving only a 2 , so that N z = 1/3, The six minima of cos(6θ Q ) correspond to the six permutations of coordinate labels, which give equivalent ordering under cubic symmetry. Somewhat different states occur for negative C 6 , but different solutions reproduce the cubic-symmetry operations.

Spin-orbit interactions and
Hund's rule coupling. As we have seen, with only spin-orbit interactions we get a spin state which has a rotational degeneracy, and with only Hund's rule interactions, the ordered phase has orbital rather than spin ordering. When both interactions are present, there is a competition between these two types of ordering. To study this competition, we need to compare the minimum eigenvalue of the spin susceptibility matrix and the orbital susceptibility one. This calculation results in the phase diagram shown in figure 1, which is not quite the same as that found in [19] for zero temperature.
The analysis of the fourth-order terms in the free energy shows that they give rise to an anisotropy that orients the staggered spin along a (1, 1, 1) direction. Since fluctuations favour In the 'spin-only' phase for η = 0, the staggered moment orients along a (1, 1, 1) direction, but the staggered spin moments of different orbital states are not collinear, thus reducing the net staggered spin. For η = 0, the mean-field state has rotational degeneracy, so no easy direction of staggered magnetization is selected and the excitation spectrum is gapless. In the orbital phase one has the sixfold anisotropy associated with the equivalent choices for differently populating orbital levels in cubic symmetry, as is discussed in the text.
the spin-only state, the phase boundary shown in figure 1 will be shifted by fluctuations to larger positive η. In the regime of orbital ordering, there is a sixfold anisotropy in the variables a 1 (Q) and a 2 (Q), as discussed above.

Mean-field Hamiltonian
Since we expect mean-field order only in the 11 components of the staggered orbital and spinorder parameters, we now derive the self-consistent mean-field equation just for these order parameters. This will allow us to study the temperature dependence of these parameters for all temperatures. To obtain these self-consistent equations, we find it convenient to write the KK Hamiltonian as in equation (3) H = 2 where ij α denotes a nearest-neighbour bond between sites i and j along the α-axis. The on-site diagonal orbital and/or spin-order parameters involve only the averages ( Q(i) is the same as 2 s(i), as defined in equation (12).) Anticipating possible antiferromagnetic order, we also define the above order parameters separately for each sublattice A and B.
Since there is one electron on each site, one has the constraint and therefore one must have the inequalities We next consider the mean-field approximation in which a product of two operators associated with different sites is approximated as follows: With this approximation, we find the mean-field Hamiltonian, where N 0 is the number of sites on the lattice. To calculate the free energy, we need to calculate the single-site partition function, e.g.
For the electron in the α-orbital, the spin can be quantized along the direction of Q α (B). Thus, where we denote y ≡ 2β , and with a similar expression for Z(B). Finally, the free energy per site becomes In minimizing F , it is important to realize that Z(A) depends only on the magnitudes Q α (B) ≡ | Q α (B)|. Therefore, the directions of the Q α s appear only in the scalar product Q α (A) · Q α (B), implying that these two vectors must be collinear with each other. Indeed, a direct minimization yields The second equation , leaving a complete rotational freedom for the direction of Q α (A) for each α separately! Although we started with a single spin per site, we ended up with three spin-like order parameters, associated with the three possible orbital states. The definitions of Q α and the constraint (20) imply that and therefore The above rotational symmetry implies that each Q α (A) can point in any direction. Therefore, the total average spin S on each sublattice can assume any direction and its z-component S z can assume any magnitude up to 1/2. All of the resulting spin states (with different total spin) are degenerate with each other. This is a very surprising property of the KK model. Since this is an intrinsic symmetry of the quantum Hamiltonian, it cannot be removed by any (thermal or quantum) fluctuations.

Possible low-temperature solutions
At T = 0, the energy is just Therefore, the ground state is highly degenerate: given equation (21), the lowest energy is E = 0.
Since the three terms in the sum are non-negative, they have to vanish individually. This can happen either when n α (A)n α (B) = 0 (and therefore also when Q α (A) · Q α (B) = 0), due to the occupation of different orbitals on the two sublattices, or when However, it turns out that as one approaches T = 0 from low finite temperatures, thermal fluctuations prefer specific ground states, with a much smaller degeneracy. This process is sometimes called 'order out of disorder' [21]. Therefore, we next concentrate on finding the lowest free energies at low T . From equation (26) Equation (25) can also be written as can have an integer value between 1 and 5. However, the self-consistency equations restrict the product Z(A)Z(B) to values between 1 and 9, with the largest value found when Z(A) = Z(B) = 3. The lowest free energy at low T thus turns out to be f min ≈ −T log 3 ≈ −1.0986T . We now review some of the self-consistency solutions.
As discussed above, it is sufficient to find the solutions for the magnitudes of the Q α 's. Specifically, the second equation in (27) and therefore also  (1, 0, 0) for the second case. Interestingly, these two limit solutions are equivalent, after an appropriate interchange of the sublattices and a permutation of the indices. It is also interesting to note that in this limit, the two sublattices are not equivalent: if one sublattice chooses to be in one orbital state, then the other sublattice has an equal amount of the two other orbital state, and vice versa. This reflects the antiferro-orbital interaction, due to the Coulomb repulsion of neighbouring electrons. The corresponding solutions represent orbitally ordered states. At low T , this solution has a free energy f orb ≈ −1.5T log 2 ≈ −1.0397T , vanishing at T = 0 and slightly higher than f min . As discussed after equation (30), all the degenerate ground states have zero energy. However, these many states may differ in their entropy, and thus approach the ground state with different free energies.
We next classify the low-temperature solutions with some 'Magnetic order', for which Q α (A) reaches a finite non-zero limit as T → 0, for some α. For such an α, at low T , equation (33) implies that Q α (B) ≈ n α (B)(1 − 2e −2yQ α (A) ) ≈ n α (B), and thus Z α (A) ≈ 1. In fact, this also follows from the weaker condition yQ α (A) 1, which could also arise when Q α (A) slowly approaches 0 at low T . If Q α (A) is finite for all α, then we find the limit solution n α (A) = n α (B) = Q α (A) = Q α (B) = 1/3 and Z(A) = Z(B) = 3. For a lower free energy we need to choose Q α (A) = − Q α (B), ending up with the free energy f 1 ≈ −T log 3 ≈ −1.0986T . As mentioned above, this solution (which we call 'M 1 ') has the lowest free energy. However, note the small difference in free energy from the 'orbital' solution, which may imply metastable states at low T . The solution M 1 has no orbital ordering. Note that in addition to the state in which all the Q α (A)'s are parallel to each other, with magnitude Q α (A) ≡ 1/3, which yields | S z | = 1/2 on the two sublattices, there exists a continuum of other solutions with the same free energy but with the Q α (A)'s having the same magnitude but pointing in arbitrary directions, i.e. having only partial magnetic ordering.
We next consider (21) implies that also n 1 (A) = 0, n 2 (A) = 0. As above, we have Q 1,2 (B) ≈ n 1,2 (B). The first equation of (27)  Finally, we consider The results now depend on how many of n 1 (A) and n 2 (A) vanish, and we end up with two basic limits, which we call M 3 and M 4 . The corresponding limiting order parameters and free energies are given in table 1, together with all the other limits discussed so far. Note that state M 4 is degenerate with state M 1 , having the lowest free energy.
In addition to the simple solutions, listed in table 1, there are solutions which combine sublattice values from different columns in that table. As explained above, the solutions in the table all assumed that when yQ α (A) 1 then also Q α (A) remains finite at T = 0. An example with a different scenario would arise, e.g., if for some α one would have log y yQ α (A) 1 (this could happen if, e.g., yQ α (A) ≈ c log y with c 1). In this case, it is easy to check that Q α (A) ≈ n α (A) → 0, but Z α (A) → 0. An example in which this occurs for α = 2 yields a solution for which n α (B) and Q α (B) approach (1/3, Depending on the particular ways in which some of the Q's decrease to zero at low T , one can now imagine other limits, combining pairs of the states listed for the two sublattices in different rows in table 1. However, all the combinations which we looked at had higher free energies than M 1 . Note that the ground state limits, M 1 and M 4 , can be represented as antiferromagnetic, since Q α (A) = − Q α (B) and (n α (A) − 1/3) = −(n α (B) − 1/3). In contrast, Orb, M 2 and M 3 all represent non-equivalent sublattices, which cannot be described by purely ferro-or antiferromagnetic order parameters. However, all the orbital states can be represented as  This is similar to equation (11), with a 1 (A) = √ 6a/2 and a 1 (B) = − √ 6b/2. The equality of two out of the three n α 's can be understood by the orbital symmetry breaking: one sublattice chooses an orbital state, and then the other two are equivalent, both on the same and on the other sublattice. The particular ratios of a and b in the ground state still need some intuitive explanation.
Finally, we note that for all the solutions discussed here one has f − f min T , which may imply difficulties in settling into the 'true' ground state at low T .

Extension to finite temperature
The self-consistency equations are easily solved numerically at finite temperatures. For the disordered phase, the results in table 1 apply at all T . For the magnetic phase M 1 , it is easy to see that equations (32) is satisfied if we set n α (A) = n α (B) ≡ 1/3 and The numerical solution to this equation is shown in figure 2 (here and below, we choose = 1/2, i.e. measure energies in units of 2 ). The corresponding free energy, given by turns out to be the lowest free energy at all T . This free energy is compared with that of the disordered state in figure 3. Interestingly, the linear dependence on T , given in table 1, provides an excellent approximation for f 1 over most of the range. The solution M 4 is degenerate with M 1 at all temperatures, and is obtained from the substitutions n α (A) = 1 3 + q, 1 3 − q, 1 3 , with q being the solution of equation ( An example of such a solution is shown in figure 4. This solution has Q α = Q β = 0, while the n α 's are as in equation (34). The upper part of the figure shows the 'order parameters' a and b, and the lower part shows the difference between the free energy of this phase and f 1 . Clearly, this difference is quite small for all T . This situation is typical for many of the 'metastable' solutions.
Before concluding this section, we note that the finite temperature extensions of the solutions M 1 and M 4 , equations (35) and (37), also minimize the Landau free energy discussed in the previous section, as can be checked explicitly.   (32) (upper frame), and the difference between its free energy and that of the lowest state (lower frame).

Perturbation theory for fluctuations
The Hamiltonian for a bond in the z-direction is the same as in equation (18), We now write down the mean-field Hamiltonian, H 0 , for site i. Unlike the discussion in the previous section, we now quantize all spins along the z-direction, and define Q α (i) as the quantum average of the z-component of the spin, Q α (i) = 2S iz a † i,α a i,α . Again, we assume a two-sublattice model in which n α (i) ≡ a † iα a iα and Q α (i) may assume different values on the two sublattices associated with the wavevector (π, π, π)/a. Then we get where site i is on one sublattice and site j is on the other.
The perturbation is written as the sum of bond terms. For a bond along α, we write its contribution V(i, j) to the perturbation as where site i is on one sublattice and site j is on the other. To evaluate the free energy when expanded about the mean-field solution we write where Z 0 is the mean-field partition function and V = ij V(i, j). So we are considering the leading effect of fluctuations, which may select from among the degenerate states of H 0 .
The self-consistent solutions obtained above indicate that we should consider the following three classes of mean-field states, which might develop different free energy in the presence of fluctuations: where n α (A) (n α (B)) is the value of n α (k) when the site k is an A (B) site and the three values in the parentheses are for α = x, y and z, respectively. Clearly, solutions |1 and |2 correspond to solution M 1 in table 1, while solution |3 corresponds to solution M 4 there. We recall the self-consistent equations for the order parameters: where Z is the single-site partition function and so that and also that 1 We first note that a property of mean-field theory is that so that We divide the perturbation into terms V d which commute with H 0 and terms V o which do not commute with H 0 . Then we may write where where ρ ≡ exp(−βH 0 ) is the mean-field density matrix. In appendix A we show that equation (49) gives the same free energy for states of classes #1 and #2, but a higher free energy for states of class #3. This result is not surprising because although the system is symmetric with respect to spin rotation within each orbital, there is no fundamental symmetry which involves a transformation between spin and orbital variables. In any case, this calculation implies that fluctuations prefer magnetic ordering over orbital ordering, thus leaving only the nine magnetic order parameters and eliminating the two orbital ones.

Renormalization group
For the pure KK model in three dimensions, our general symmetry arguments imply no longrange order. Choosing a specific axis, say the z-axis, the spins associated with the z-orbital, Q z (i), would have no interactions along the z-axis. Therefore, the Hamiltonian related to these spins would split into a sum over non-interacting planes. Within each such plane, one would have a three-component Heisenberg spin problem, with no long-range order at any non-zero temperature.
As discussed above, many small perturbations break this ideal symmetry, yielding some form of long-range order. In an extreme case, such perturbations would generate dispersion of the spin correlations along all directions, reducing the problem to that of an 11-component order-parameter model. The simplest interaction that would generate such dispersion would be an additional direct hopping between nearest-neighbour Ti ions. In the extreme isotropic case, all the 11 diagonal propagators have the isotropic form χ α (q) = 1/(r α + q 2 ) (where q, measured from the ordering vector Q, is assumed small, and initially all the r α 's are degenerate, T c )). The quartic spin terms would then be as given in equation (17). In Fourier space, these would have the form where it is understood that in each product of four order parameters, these are taken at the four wave vectors q 1 , q 2 , q 3 and q 4 , and the wave vectors are over the Brillouin zone.
For this isotropic case, as well as for any model with dispersion along all directions, the Ginzburg criterion identifies the upper critical dimension as d u = 4, and we can now follow the standard renormalization group analysis in d = 4 − ε dimensions [24]. The generalization to d dimensions is somewhat ambiguous, since in our model the number of orbitals is actually equal to d. However, one way to proceed is to keep this number equal to three, while performing the momentum integral in general d dimensions. At the end one extrapolates d back to three anyway. A similar ambiguity occurs for the dipolar interactions, which also couple the dimension of the order parameter with that of the space [24]. Other possibilities, e.g. replacing the 11-dimensional order parameter space by one of the dimension 4d − 1, give similar qualitative results.
It is now straightforward to follow the standard renormalization group procedure [24]: integrate out shells in momentum space, and then rescale spin and space to return to the original form of the free energy. The resulting recursion relations for the r α 's immediately split between the nine spin variables (which remain equal to each other, denoted by r S ) and the remaining two orbital ones (denoted r O ). To leading order in the r α 's, where n(= 3) is the number of spin components, u ∝ and A is related to the area of a d-dimensional sphere [24]. After iteration, we have r O > r S , indicating ordering of the spin degrees of freedom at a higher temperature. Thus, thermal fluctuations remove the degeneracy between the orbital-and spin-order parameters. Close to the transition, the orbital degrees of freedom become non-critical, and they can be integrated out of the free energy, leaving only the nine spin degrees of freedom. This conclusion is the same as the one obtained from the perturbative treatment in section 5. The remaining quartic free energy is now a special case of the general nm-component model, with n = m = 3 (see e.g. [24]). It can then be parametrized in the following way, where α, β = 1, 2, . . . , m, and µ, ν = 1, 2, . . . , n. For m = n = 3 it is well known [24,25] that the critical behaviour is dominated by the decoupled fixed point, where u * = 0 and each of the spins s α exhibits the usual Heisenberg critical behaviour for a three-component vector model.
The actual spin values in the ordered phase are then determined by the small irrelevant values of u. In our case, the initial Hamiltonian is probably far away from this isotropic limit. An alternative starting point would be to return to the model discussed above, with the propagators as given in equations (15) and (16). Repeating the Ginzburg criterion for this case, the upper critical dimension becomes equal to d u = 5, leading to an ε expansion in d = 5 − ε dimensions. Again, the recursion relations for the r α 's split into 9 + 2, leaving only the spin degrees of freedom. For second order in u and v, the recursion relations for u and v now require integrals of the form where the integral is over a shell in q and where is the propagator generated from the Gaussian average s αµ (q)s αµ (−q) . In our case, G αµ (q) has no dispersion at all for q along the α-axis. As a result, the integrals I αµβν diverge for d < 5 when α = β, but are less divergent for α = β. A similar behaviour was encountered for the displacive phase transitions in some perovskites [26]. However, in that case the dispersion along the α-axis was generalized to have a leading Lifshitz-like behaviour, q 2L α , resulting in a divergence below d = 5 − 1/L. In our case, L = ∞, and the upper critical dimension is shifted to d u = 5. Following similar steps to those used for Lifshitz points, and keeping only the most divergent integrals I αµαµ ≡ I, the recursion relations become du d = εu − 4I[(nm + 4)u 2 + 2(n + 2)uv] + · · · , where ε = 5 − d. For leading order in ε, these recursion relations do not have any stable fixed point. However, the same arguments used for the dispersion-less case [24] can be repeated here: since u represents an energy-energy coupling between the decoupled spin models, the stability exponent of the decoupled fixed point is still λ = α/ν, where α and ν now represent the specific heat and the correlation length exponents of the three-component Heisenberg model-except that now this is the model with the dispersionless behaviour. If α remains negative for this model, as one might anticipate (α usually becomes more negative as one moves further below d u ), then one ends up again with the decoupled behaviour. Interestingly, for n < 32/9 there exist no other real fixed points, except the Gaussian (u * = v * = 0) and the decoupled ones. However, for n > 32/9 there appear more non-trivial fixed points. Note that the shift of the upper critical dimension from 4 to 5 is accompanied by a similar shift of the lower critical dimension from 2 to 3, as noted in our discussion of the Mermin-Wagner proof for the lack of long-range order at d = 3. Therefore, all the above discussion probably applies only for 3 < d < 5, and is thus somewhat academic. However, as mentioned earlier, many symmetry breaking terms will remove these effects, and restore long-range order. The same anisotropies will also restore the convergence of the integral I above four dimensions, and map the problem back to the 'usual' (but anisotropic) nm-vector problem. The anisotropies will then break the degeneracy of the nine spin components, and end up with a slow crossover towards the critical behaviour appropriate for those reduced symmetries in three dimensions.
Since the exact cubic KK model has no long-range order at d = 3, the above renormalization group discussion for d > 3 is somewhat academic. However, it may serve as a starting point for adding small perturbations and addressing the crossover between this extreme critical behaviour and that of the models with dispersion. We leave this detailed analysis for some future works.
We next discuss the renormalization group equations for some of the lower symmetry cases discussed in section 3. When we have both spin-orbit and Hund's coupling, the problem reduces to one with only three spin variables, denoted by ξ x , ξ y and ξ z . The quadratic terms now have dispersion, with Since usually c 1 = c 2 , the upper critical dimension returns to four. The recursion relations for the r α 's again yield r S < r O , and we end up with the usual three-component spin problem, yielding the usual G-type antiferromagnetic order of the ξ's.
When we have nnn hopping, then again fluctuations prefer only spin ordering. Now the spins are expressed in terms of two vectors ξ and ρ. The renormalization group analysis of such a Hamiltonian, in 4 − ε, was carried out in [27]. The resulting recursion relations had no stable fixed point, which probably implies a flow towards an unstable quartic term indicating a first-order transition.

A nearly-cubic real system: LaTiO 3
It has been proposed [13] that perhaps the cubic KK model studied above could be taken as a starting point for a successful interpretation of LaTiO 3 (LTO). From what we have seen, the cubic model has some very unusual symmetries and consequently some very unusual properties. On the other hand, LTO seems not very different from a conventional Heisenberg antiferromagnetic insulator. Indeed, the temperature dependence of the order parameter and the wavevector dependence of the spin-wave energy given by Keimer et al [27], and shown in figure 5, is well described by a nn superexchange coupling having the value 15.5 meV, accompanied by a weak ferromagnetic moment consistent with a small Dzyaloshinskii-Moriya interaction of about 1.1 meV [27]. The ordered magnetic moment has been found to be rather small, 0.46−0.57µ B [11,27], leading to several attempts to explain the magnetic properties of this material. Thus, although this system may have some properties that require an explanation, it is nevertheless, to a first approximation, a conventional antiferromagnet. Accordingly, its Hamiltonian must differ in important respects from that of the cubic KK model studied here.
In particular, the threefold degeneracy of the t 2g levels, which prevails at strictly cubic symmetry, which is frequently lifted in real materials by the JT distortion, is lifted in the antiferromagnetic Mott insulator LaTiO 3 (T N ≈ 146 K). The JT effect in LaTiO 3 is caused by the twisting of the Ti-O bonds with respect to each other. The distortion leads to a crystal field that splits the t 2g levels [11], yielding a crystal-field gap of about 0.2 eV between the orbitally non-degenerate ground state and the first excited level, a value that has been confirmed by a study of photoelectron spectroscopy [28]. Thus, instead of the KK model, which is based on triply degenerate orbitals, LaTiO 3 is actually described by an orbital singlet. The non-degenerate ground-state orbital due to the crystal-field calculations given in [11] is consistent with the orbital order found in NMR measurements of the Ti-3d quadrupole moment [30]. The presence of orbital order at low temperatures has also been concluded from measurements of the dielectric properties and the dynamical conductivity [31]. A model for the magnetism of LaTiO 3 , which is based on the crystal-field calculation given in [11] is presented in [32]. This calculation starts from a point-charge summation of the static crystal field for the Ti ions, employing a full Madelung sum over the crystal. A recent calculation of the crystal field, which includes also the covalency contribution and the spin-orbit coupling [33], has led to roughly the same t 2g splitting scheme, but with a larger spacing between the t 2g and the e g levels.
Schmitz et al [32] have used a Slater-Koster parametrization of the Ti-O hopping to obtain an effective hopping matrix between the d orbitals of nearest-neighbour Ti ions. Treating this Ti-Ti hopping and the on-site spin-orbit coupling as perturbations, they calculated the superexchange coupling between the non-degenerate crystal-field ground states of the Ti 3+ ions. Because of the rather low symmetry of the system, in treating the Ti 2+ ions which appear as intermediate states of the exchange processes, it has been necessary to take into account the full on-site Coulomb interaction matrix. This, together with the spin-orbit coupling, gave rise to antisymmetric and symmetric anisotropies of the effective single-bond spin Hamiltonian.
The single-bond spin Hamiltonian is the basis for the magnetic Hamiltonian, from which the magnetic order of the classical ground state follows. The latter is constructed by decomposing the entire Ti-lattice into four sublattices (each magnetic unit cell includes four Ti ions, just as the crystallographic unit cell). Assigning a fixed-magnitude magnetization (per site) to each sublattice, M i , the macroscopic magnetic Hamiltonian takes the form [32] in which I ij is the macroscopic isotropic coupling, D ij 's are the Dzyaloshinskii vectors (to leading order in the spin-orbit coupling), which are the macroscopic antisymmetric anisotropies, and ij are the macroscopic symmetric anisotropy tensors (of second order in the spin-orbit coupling). Although all four sublattice magnetizations are of equal magnitude, their directions are all different. In general, according to the space group Pbnm symmetries of LaTiO 3 , there are four possibilities for the symmetry of sublattice magnetizations of the classical ground state. The minimization of the magnetic Hamiltonian yielded the magnetic structure shown in figure 6. This ground state is predominantly a G-type antiferromagnet along the x-direction, as reported in [11]. The z components of the magnetizations order ferromagnetically, again in accordance with the experiment. Finally, the y components of the magnetizations order antiferromagnetically, in an A-type structure, which has not been detected experimentally. However, the experiment of [11] is not sensitive to a small moment along the y-axis, and hence the small A-type antiferromagnetic order along this direction does not contradict the data. Symmetry allows for such ordering, given the G-type order along x and the ferromagnetic order along z. Indeed, in theYTiO 3 -system, which has the same space group as LaTiO 3 , such an order has been detected [34], but with different magnitudes of the canting angles, which cause the ferromagnetic order to dominate.
Equation (58) indicates that LaTiO 3 is a system of localized spins with S = 1/2, governed by an anisotropic spin Hamiltonian. However, even in the effective spin-1/2 ground manifold, the effects of the orbital degrees of freedom are not completely absent. For instance, because spin-orbit perturbations connect the ground and excited crystal-field levels, one has the so-called Van Vleck (VV) susceptibility tensor, χ VV . It is reasonable to calculate χ VV using the single-ion crystal-field states, in which case the α − β component of χ VV is given by where µ B is the Bohr magneton and |0 and |e are the orbital ground and excited crystal-field states, respectively. The leading perturbation of the magnetic moment due to the spin-orbit interaction can be calculated using the perturbed single-ion states (indicated by a subscript λ) which are where the Greek indices label spin wavefunctions and λ = 0.018 eV [34]. Using these perturbed states, we calculate the moment (exclusively orbital at this order) which is induced by the spinorbit perturbation, λl · s. Within the manifold of spin states associated with the crystal-field ground states, we find the operator relation The conclusion is simple: the zero-point orbital magnetic moment is directly proportional to the VV zero-temperature susceptibility. This is to be expected if one views the spin-orbit interaction as an effective field acting on the spin operator. The static crystal field wavefunctions of Schmitz et al [32] give

Discussion and summary
The cubic KK model has some very unusual and interesting symmetries, which cause meanfield theory to have some unusual features. In particular, for the simplest KK Hamiltonian, we found that mean-field theory leads to criticality for the wavevector-dependent spin susceptibility associated with orbital α, which is dispersionless along the q α direction of the wavevector. This result is consistent with the previous observation [15] that the Hamiltonian is invariant against an arbitrary rotation of the total spin in the orbital α summed over all spins in any single plane perpendicular to the α-axis. This 'soft mode' behaviour prevents the development of long-range spin order at any non-zero temperature [15], even though the system is a three-dimensional one.
Any perturbation which destroys this peculiar symmetry will enable the system to develop long-range spin order. In particular, we investigate the role of (a) spin-orbit interactions, (b) second-neighbour hopping and (c) Hund's rule coupling in stabilizing long-range spin order.
In the presence of spin-orbit interaction, we find wavevector selection (because now the spin of different orbitals cannot be freely rotated relative to one another) into a two-sublattice antiferromagnetic state with a greatly reduced spin magnitude. Since experiment shows such a reduction [8], this mechanism may be operative to some extent. However, as noted previously [15], the excitation spectrum does not have a gap until further perturbations are also included. The mean-field solution is consistent with this conclusion, because the mean-field state which minimizes the trial free energy is degenerate with respect to a global rotation of the staggered spin.
The ordered state which results when nnn hopping is added to the bare KK Hamiltonian is quite unusual. In this state, each orbital has a staggered spin moment, but these three staggered spin moments form a 120 • state such that the total staggered spin moment (summed over the three orbital states) is zero! It is not immediately obvious how such long-range order would be observed. Finally, we show that when the bare KK Hamiltonian is perturbed by the addition of only Hund's rule coupling, the resulting ordered state may exhibit long-range antiferromagnetic orbital order.
One caveat concerning our result should be mentioned. All our results are based on a stability analysis of the disordered phase. If the ordering transition is a discontinuous one, our results might not reveal such a transition. This could explain why our results for the first appearance of ordering as the temperature is reduced give a phase diagram which is not quite the same as that found in [19] at zero temperature. In section 6 we presented results for the temperaturedependence of the various mean-field solutions. Further analysis of the ordered phase is needed to fit our results onto those of [19].
One general conclusion from our work is that it is not safe to associate properties of real experimental systems with properties of a model Hamiltonian unless one is absolutely sure that the real system is a realization (at least in all important aspects) of the model Hamiltonian. Here the ideal cubic KK Hamiltonian has properties that are quite different from those observed for systems it supposedly describes. What this means is that it is necessary to take into account local distortions from cubic symmetry, as we reviewed in section 7 for LaTiO 3 . Alternatively, perhaps our work will inspire experimentalists to find systems that are as close as possible to that of the ideal cubic KK Hamiltonian treated here. Such systems would have quite striking and anomalous properties, such as spin-spin correlations which have different extreme anisotropies depending on the associated orbital occupancy.

A.1. Diagonal terms
We first look at the diagonal terms in V , i.e. those which commute with H. For a bond along α the diagonal terms are Note that the diagonal terms are of the form where the sum is over pairs (i, j) of nearest neighbouring sites. This enables us to write ×( 2S jz a † jβ a jβ a † jβ a jβ − a † jβ a jβ 2S jz a † jβ a jβ ) (A.4)

A.2. Off-diagonal terms
Now we consider terms in V which do not commute with H. We work in a representation in which S iz is diagonal. The off-diagonal terms for a bond along α are Note that cross terms involving more than one bond vanish because the states of both sites are changed by this off-diagonal perturbation. We want to construct where µ = 2 λ, we only kept terms which survive the trace, and (a † jγ a jγ n γ (i) + 2S jz a † jγ a jγ Q γ (i)) S j+ a † jβ a jβ (A.7n) In the Q's and R's we used the fact that a † iβ a iβ a † iβ a iβ vanishes unless β = β. Then (noting that γ a † iγ a iγ = 1) we have  We always take averages over a mean-field density matrix, which means that averages are taken separately for each site. So we have (A. 12) and interchanging the roles of β andβ we also havē . (A.13) Now use Q β (j) = 1 3 tanh(2β Q β (j)). (A.14) This relation indicates that kT c = 2 /3 and that Using these relations we writē .

(A.17)
This establishes thatT =Z. Similarly, one can show thatR =Q. These also follow from the fact that the contributions from two terms, one involving VV † and the other involving V † V , are equal. Thus

A.3. Comparison of states #1 and #2
In states #1 and #2, n aρ = 1/3, so the only terms that can differentiate between these states are the off-diagonal term with β =β. In table A.2 we list the relevant factors. So we now look at the difference between the free energies of states #1 and #2 (which are indicated by superscripts) (1) − (U i U j + W i W j + Z i Z j ) (2) ]. (A. 19) Here we used U i = Y i , etc. Also we noted that the contributions from Z i Z j and T i T j are equal. Out of the six terms in the sum over β andβ, we will only keep the four identical terms for which Q i,β = Q i,β . Then After the integration is performed and after applying equation (A.15), we see that F (1) = F (2) . So, as guaranteed by spin symmetry, both states are equally affected by quantum fluctuations.  N(β,ββ) is given by the appropriate linear combination of the n's. In the columns headed by the values of the indices β andβ we give the value of dn/dQ(±1), which is used to construct N(β,ββ).
(β,β) (β,β) Term Pre n x, y x, z y, x y, z z, x z, y Term Pre n x, y x, z y, x y, z z, x z, y When N j (β,β) = ±4, then where we used equation (46). Now we count the contributions for δ off F . For state #3, we see from table A.3 that T 1 gives four cases of τ and two cases of ρ, T 2 gives four cases of τ and two of ρ, 2T 3 gives eight cases of τ and four of ρ. So from these terms we get 16 cases of τ and eight of ρ. From 2T 4  For state #1 note that n ± β (i) and n ± β (j) do not depend on β. So from T 1 we get six cases of ρ and from T 2 we also get six cases of ρ. From T 3 we get 12 cases of τ and from 2T 4 we get 12 cases of τ. So δ off F (1) = − This shows that F < 0, i.e. the lowest-order perturbation due to fluctuations causes the free energy of state #1 to be less than that of state #3.