Non-Coplanar Model States in Quantum Magnetism Applications of the High-Order Coupled Cluster Method

Coplanar model states for applications of the coupled cluster method (CCM) to problems in quantum magnetism are those in which all spins lie in a plane, whereas three-dimensional (3D) model states are, by contrast, non-coplanar ones in which all the spins do not lie in any single plane. Here we extend the CCM to non-coplanar / 3D model states and we present results for three cases: (a) the spin-half one-dimensional Ising ferromagnet in an applied transverse magnetic field (as an exactly solvable test model to use as a yardstick for the viability and accuracy of our new methodology); (b) the spin-half triangular-lattice Heisenberg antiferromagnet in the presence of an external magnetic field; and (c) the spin-$S$ triangular-lattice {\it XXZ} antiferromagnet in the presence of an external magnetic field, for the cases $\frac{1}{2} \leq S \leq5 $. For 3D model states the sets of algebraic CCM equations for the ket- and bra-state correlation coefficients become complex-valued, but ground-state expectation values of all physical observables are manifestly real numbers, as required. Excellent correspondence is seen with the results of other methods, where they exist, for these systems. CCM results demonstrate explicitly that coplanar ordering is favoured over non-coplanar ordering for the triangular-lattice spin-half Heisenberg antiferromagnet at all values of the applied external magnetic field, whereas for the anisotropic {\it XXZ} model non-coplanar ordering can be favoured in some regions of the parameter space. Specifically, we present a precise determination of the boundary (i.e., the critical value of the {\it XXZ} anisotropy parameter $\Delta$) between a 3D ground state and a coplanar ground state for the {\it XXZ} model for values for the external magnetic field near to saturation, for values of the spin quantum number $S \leq 5$.


I. INTRODUCTION
The coupled cluster method (CCM) [1][2][3][4][5][6][7][8][9][10] is a powerful method of quantum many-body theory that has long been used to study strongly interacting and highly frustrated quantum spin systems with great success . The introduction of the "high-order" CCM [15][16][17] for these systems has led to a step-change in its accuracy. The high-order CCM employs very high orders of approximation schemes, for which the equations that determine all multispin correlations retained at any given level are both derived and subsequently solved by using massively parallel computational tools [15][16][17]37]. The CCM is now fully competitive with the best of other approximate methods, especially for systems of N spins on the sites of a lattice in two (see, e.g., Refs. [27,31,32,36] and references cited therein) or three (see, e.g., Refs. [33,34] and references cited therein) spatial dimensions. Unlike several other approximate quantum many-body methods that are limited in their range of applicability by frustration (i.e., where bonds in the Hamiltonian compete against each other to achieve energy minimisation), the CCM has been applied previously even to highly frustrated and strongly correlated quantum spin systems with much success. Thus, recently it has been demonstrated, for example, that the high accuracy needed to investigate the quantum ground-state selection of competing states of the kagome antiferromagnet is provided by high-order CCM calculations [27,32,33]. Another advantage of the high-order CCM is that it is very flexible. For example, both in principle and in practice, the CCM technique can treat essentially all Hamiltonians containing either single spin operators and/or products of two spin operators, on any crystallographic lattice, and for any spin quantum number S.
We note too that, unlike most alternative techniques, the CCM can be applied from the outset in the thermodynamic limit, N → ∞, at every level of approximate implementation, thereby obviating the need for any finite-size scaling of the results.
In all practical implementations of the CCM the many-body correlations present in the exact (ground or excited) state of the system under investigation are expressed with respect to a suitable model (or reference) state, as we explain in more detail below in Sec. II.
Coplanar model spin states used in the CCM are those states in which all spins lie in a plane, whereas three-dimensional (3D) model spin states are non-coplanar states in which the spins do not lie in any one plane. We remark that, until now, only coplanar model spin states have been used in all prior CCM calculations in the field of quantum magnetism for reasons that we now explain.
Thus, an important ingredient used in all practical applications of the CCM to spin-lattice systems  is to rotate the local spin axes of all spins in the model state such that they appear (mathematically only) to be point in the "downwards" z-direction. One can always choose a set of rotations that leads to terms in the Hamiltonian that contain only real-valued coefficients, with respect to the new set of local spin axes, for the coplanar model states.
By contrast, three-dimensional (3D) (non-coplanar) model states inevitably lead to terms in the new Hamiltonian after rotation of the local spin axes that contain complex-valued coefficients. These cases are more difficult to treat both analytically and computationally.
Of course, all macroscopic physical parameters calculated within the CCM, such as the ground-state energy and magnetic order parameter, still have to be real numbers because the transformations of local spin axes are unitary and the resulting Hamiltonian is still Hermitian. Nevertheless, the intervening multispin correlation coefficients are necessarily complex-valued quantities.
In this article, we explain how we can carry out CCM calculations for such Hamiltonians that contain terms in the Hamiltonian after rotation of local spin axes with complex-valued coefficients. We show that the amendments to the existing CCM code for spin-lattice models [37] to be able to treat such cases is non-trivial. In order to illustrate the new technique, we present three separate applications to models of considerable interest in quantum magnetism. As a first test of the new methodology we present results in Sec. III for the exactly solvable one-dimensional Ising model in a transverse external magnetic field, and an explicit analytical calculation of the lowest-order implementation of the CCM is presented in detail for this model in Appendix A. Secondly, in Sec. IV we then describe results for the spin-half Heisenberg model on the triangular lattice at zero temperature in the presence of an external magnetic field, in which we make explicit use of the 3D "umbrella" state as our CCM model state, and an explicit derivation of the Hamiltonian after the rotations of the local spin axes for this state is presented in Appendix B. Lastly, in Sec. V, the phase diagram of the spin-half XXZ model on the triangular lattice at zero temperature, also in the presence of an external magnetic field (near saturation), is examined. Here we again employ the 3D "umbrella" state as a possible CCM model state, and we show how its use now leads to an improved quantitative description and understanding of this model. We conclude with a brief summary of our results in Sec. VI.

II. THE COUPLED CLUSTER METHOD (CCM)
A. Ground-State Formalism As the methodology of the CCM has been discussed extensively elsewhere , only a brief overview of the method is presented here. The ground-state Schrödinger equations are given byĤ in terms of the HamiltonianĤ, and where formally, for normalisation, we require The bra and ket states for our N-spin system (with each of the spins carrying the spin quantum number S) are parametrised independently in the forms within the normal coupled cluster method, in terms of the multi-configurational CCM (creation and destruction) correlation operators,Ŝ andŜ, respectively. The index I here is a set-index that denotes a set of lattice sites, I = {i 1 , i 2 , · · · , i n ; n = 1, 2, · · · 2SN}, in which each site may appear no more than 2S times, for reasons we describe below. We shall be interested specifically in the case of infinite systems, N → ∞. Note thatĈ + 0 ≡1 is defined to be the identity operator in the many-body Hilbert space, the operatorsĈ + I andĈ − I ≡ (Ĉ + I ) † are respectively, ∀I = 0, multispin creation and destruction operators, for clusters of up to N spins, which are defined more fully below, and S I andS I are the CCM ground-state ket-and bra-state (c-number) multispin correlation coefficients, respectively. We use model states (denoted |Φ for the ket state and Φ| for the bra state) as references states for the CCM. The ket state |Φ is required to be a fiducial vector (or cyclic vector) with respect to the complete set of mutually commuting, multispin creation operators {Ĉ + I }. Equivalently, the set of states {Ĉ + I |Φ } is a complete basis for the ket-state Hilbert space. Furthermore, |Φ is also defined to be a generalised vacuum state with respect to the set of operators We note that, with these conditions fulfilled, the exact ground-state ket-and bra-state wave functions, |Ψ and Ψ |, respectively, now satisfy the normalisation conditions We now define the ground-state energy functional,H ≡ Ψ |Ĥ|Ψ = Φ|Ŝe −ŜĤ eŜ|Φ , such that the CCM ket-and bra-state equations are given by extremisingH with respect to all of the CCM multispin correlation coefficients, With these equations satisfied, the CCM ground-state energy is now given by Equation (9) in these rotated local spin-space frames, where i k denotes an arbitrary lattice site, | ↓ i k is the "downward-pointing" state of a spin on site i k with spin quantum number S (i.e., defined The CCM formalism would be exact if all possible multispin cluster correlations could be included in the operatorsŜ andŜ. However, this is normally impossible to achieve practically. In most cases, systematic approximation schemes are used to truncate the respective summations in Eqs. (3) and (4) for these operators, by restricting the sets of multispin configurations {I} to some manageable subset within some hierarchical scheme that becomes exact in the limit that all configurations are retained. In the present paper we use two schemes that are denoted as the SUBn-n and LSUBn schemes, respectively. The more general SUBn-m scheme retains all correlations involving only n or fewer spin flips (with respect to the respective model state |Φ ) that span a range of no more than m contiguous lattice sites. By contrast, in the localised LSUBn scheme all multispin correlations over all distinct locales on the lattice defined by n or fewer contiguous sites are retained. Each spin flip is defined to require the action of a spin-raising operatorŝ + in acting just once, and a set of lattice sites is said to be contiguous if every site of the set is a nearest neighbour (in some specified lattice geometry) to at least one other member of the set. The LSUBn and SUBn-n schemes are thus identical only for the limiting case when S = 1/2. For higher spin quantum numbers S, the LSUBm scheme is equivalent to the SUBn-m scheme if and only if n = 2Sm. Spin-cluster configurations I that are equivalent under the space-and point-group symmetries of the crystallographic lattice (as well as of both the Hamiltonian and the model state under consideration) are counted only once by explicitly incorporating those symmetries into the calculation, and these clusters are referred to as "fundamental clusters". The number of such fundamental clusters used for the ground-state expansions for |Ψ and Ψ | at the respective nth-order level of (either LSUBn or SUBn-n) approximation to Eqs. (3) and (4) is denoted by N f (n).
Although, formally, the CCM correlation operatorsŜ andŜ of Eqs. (3) and (4) must obey the condition which is implied by Hermiticity, in practice this may not be exactly fulfilled at finite levels of (LSUBn or SUBn-n) approximate implementation, due to the independent parametrisations of the two operators. However, this minor drawback of the CCM is far outweighed in practice by the two huge advantages that the method exactly obeys both the Goldstone linked-cluster theorem and the very important Hellmann-Feynman theorem at all levels in the approximation hierarchies. The former implies that we can work from the outset in the required thermodynamic limit of an infinite number of spins, N → ∞, while the latter implies that the expectation values of all physical parameters are calculated within the CCM on the same footing as the energy and in a fully self-consistent manner.
Unlike in many other competing formulations of quantum many-body theory, the CCM thus never needs any finite-size scaling of the results obtained with it. Indeed, the sole approximation that is ever made within any application of the CCM is to extrapolate the results obtained for any physical parameter within the (LSUBn or SUBn-n) approximation hierarchy used to the limit n → ∞ where the method becomes exact. By now, a great deal of experience has been acquired on how to perform such CCM extrapolations, and we allude here to one such calculation in Sec. III, and invite the reader to consult the literature cited above for further details.

B. Computational Aspects for 3D Model States
The CCM equations (7) and (8) may be readily derived and solved analytically at low orders of approximation. A full explanation of how this is carried out for the LSUB1 approximation for the spin-half ferromagnetic Ising chain in a transverse magnetic field, which we study in Sec. III is given in Appendix A. Highly intensive computational methods [15][16][17] are essential at higher orders of LSUBn or SUBn-n approximation because the number N f (n) of fundamental clusters (and so therefore also the computational resources necessary to store and solve them) scales approximately exponentially with the order n of the approximation scheme being used. There are four distinct steps to perform in carrying out high-order CCM calculations for the ground state for "3D model states," each of which has a counterpart in the "standard" CCM code [37] that pertains only to coplanar states. (As a short-hand only, we shall refer to any case that results in the Hamiltonian containing terms with coefficients that are complex-valued after rotation of the local spin axes to be a "3D model state," although clearly these are some essentially artificial cases, e.g., the transverse Ising model presented below, where the model state might be coplanar.) The first step is to read in "CCM script files" that define the basic problem to be solved.
We remark that the derivation of Hamiltonians after rotation of local spin axes for the 3D model states is non-trivial because we must carry out at least two sets of rotations. An example of this process is given for the spin-half triangular-lattice Heisenberg model in the presence of an external magnetic field in Appendix B. As a consequence the resulting CCM script file is much longer than for coplanar model states because we now have terms in the Hamiltonian with both real and imaginary coefficients.
The second step involves the enumeration of all connected clusters (also called "lattice animals") and all disconnected clusters that are to be retained at a given level of LSUBn or SUBn-n approximation for a given lattice and spin quantum number S that are distinct under the lattice, model state and Hamiltonian symmetries (and perhaps that also satisfy some such conservation rule as s z T = 0, whereŝ z T ≡ N i=1ŝ z i , which would pertain, for example, to all models whose Hamiltonians contain only spins interacting pairwise via isotropic Heisenberg exchange interactions). This step is no more difficult for 3D model states than for coplanar model states, although clearly this step is itself highly non-trivial to perform computationally.
The third step involves deriving and storing the basic CCM ground-state equations. In order to find these equations, we first partition the multispin cluster configuration pertaining to the set index I for the operatorĈ − I in the ket-state equation Φ|Ĉ − I e −ŜĤ eŜ|Φ = 0 of Eq. (7) into the products of "high-order CCM operators" [15][16][17]. There are a huge number of partitions potentially and each term in a new potential contribution to the ket-state equations must be tested for suitability, i.e., all subclusters are checked against a list of fundamental clusters after any appropriate space-and point group symmetries (plus any applicable conservation laws) have been employed. This is arguably the most difficult step in carrying out any high-order CCM calculation, and effectively we must run this code twice for the 3D model states: once for the terms in the Hamiltonian with real coefficients and again for the terms with imaginary coefficients.
The fourth step is to solve the ground-state ket and bra equations and to obtain the ground-state expectation values. For 3D model states, complex-number algebra must be implemented for all subroutines that solve the ket-and bra-state equations (solved by "direct iteration" for 3D model states), and also in those subroutines that determine expectation values such as the ground-state energy of Eq. (9) or other expectation values (i.e.,Ā = Φ|Ŝe −ŜÂ eŜ|Φ ). This is achieved by using options in the C++ compiler.
Although the process of updating the existing high-order CCM code [37] for coplanar states so as to be able now also to utilise 3D model states (resulting in Hamiltonians with terms involving complex-valued coefficients) is therefore straightforward in principle, this process is actually considerably less so in practice because the CCM code [37] itself is extensive and complex. In order to validate the new code, we show in Sec. III that analytical low-order LSUB1 results derived in Appendix A for the transverse Ising model are replicated by the new 3D CCM code. Similarly, for the same model, for higher orders of LSUBn approximation with n ≤ 12, we also show that the new code exactly replicates the corresponding results obtained using the "standard" code [37]. Both of these results are excellent tests of the new code. Furthermore, all results for each of the three models considered in Secs. III, IV and V are in excellent agreement with the results of other methods (where they exist).
Finally, we remark again that the creation of the CCM script files is more complicated for 3D model states than for coplanar model states.

TERNAL MAGNETIC FIELD
We take as a first example to demonstrate the feasibility and the accuracy of the new CCM approach an exactly solvable model, namely the one-dimensional (1D) spin- 1 2 Ising ferromagnet in a transverse magnetic field [38]. The corresponding Hamiltonian is given bŷ where the index k runs over all lattice sites on the linear chain (with site N + 1 equivalent to The strength of the applied external transverse magnetic field is given by λ. Clearly, in the case λ = 0 with no field applied, the spins are ferromagnetically aligned along the z-direction. Similarly, for high enough values of λ it is clear that the spins will align along the transverse (x-)direction. The CCM model state that we choose for this system is one in which all spins point in the downwards z-direction. Thus, the model state is expected to be better for low values of λ, particularly those below the phase transition that separates the two regimes where the spins are respectively canted to align along some Classically, the spins are canted at an angle α from the (say, downwards) z-direction in the presence of the transverse magnetic field λ. It is trivial to see that the classical groundstate energy E cl g is minimised for α = sin −1 λ for λ ≤ 1. There is then a classical phase transition at λ = λ cl c ≡ 1, such that for λ ≥ λ cl c , the spins are all aligned in the direction of the transverse field, with α = 1 2 π. We thus have that the classical ground-state energy per spin is given by Similarly, the classical values of the magnetisations in the z-direction (i.e., the Ising direction), M z , and in the transverse x-direction (i.e., the field direction), M trans. , are trivially found to be given by and M trans.  The quantum spin-1 2 version of the model can be exactly solved [38]. Thus we also have available to us the corresponding exact expressions for the classical parameters given above in Eqs. (14)- (16), against which we can compare our CCM results. Particular interest attaches to the model due to the fact that the classical phase transition at λ = λ cl c ≡ 1 is now shifted to the point λ = λ c ≡ 1 2 . The exact ground-state energy per spin is given by [38] E which expression is nonanalytic at the quantum phase transition point λ c = 1 2 . It is simple to check that in the two extremes λ → 0 and λ → ∞, Eq. (17) reduces respectively to the two limiting values, E g (λ = 0)/N = − 1 4 and E g (λ → ∞)/N → − 1 2 λ, exactly as for the classical case given by Eq. (14). This is also just as expected, since in these two limits the fully aligned ferromagnetic states are also eigenstates of the quantum Hamiltonian. Precisely at the quantum phase transition point, Eq. (17) yields the value E g (λ = 1 2 )/N = − 1 π . The classical result of Eq. (14) is compared with its exact counterpart of Eq. (17) in Fig. 1.
The corresponding exact result for the magnetisation in the (Ising) z-direction, M z , is given by [38] which now exhibits the phase transition at λ = λ c ≡ 1 2 much more clearly than Eq. (17) for the ground-state energy. Once again, the classical and exact results for M z , from Eqs. (15) and (18) respectively, are compared in Fig. 2. Finally, the exact result for the transverse magnetisation (i.e., in the field direction) is given by [38] which is again nonanalytic at the quantum phase transition point, λ c = 1 2 , where it takes the value M trans. (λ = 1 2 ) = 1 π . It is easy to confirm that when λ varies from zero to ∞, M trans. from Eq. (19) varies smoothly from zero to 1 2 , as shown in Fig. 3 where it is also compared to its classical counterpart of Eq. (16).
For present purposes we now wish to illustrate how the CCM can be applied when we carry out a unitary transformation of the local spin axes that leads to terms in the Hamiltonian with complex-valued coefficients. We use the unitary rotation of the local spin axes (now for all sites k on the linear chain) given bŷ which simply is equivalent to rotating the transverse field from the x-to the y-direction, while leaving the spins aligned in the (negative) z-direction. That leads to an alternative representation of the model given by the Hamiltonian where the term with an imaginary coefficient now appears in the transverse external field part of the Hamiltonian. However, we remark again that the eigenvalue spectrum for this Hamiltonian should not change compared to that for the Hamiltonian of Eq. (13). Note too that this rotation of the spins in the xy-plane does not affect the model state for this system, namely, one in which all spins point in the downwards z-direction.
CCM LSUB1 calculations for both Hamiltonians of Eqs. (13) and (21)  We note that convergence of the LSUBn sequences of approximants becomes worse for larger values of λ ( > ∼ 0.5), exactly as expected, since this region is precisely where the model state becomes a poorer starting point for the CCM calculations, due to the quantum phase transition that occurs at λ c = 1 2 . Nevertheless, it is clear by inspection of Figs. 1 and 3 that results for the ground-state energy and also M trans. compare extremely well with the exact results of Ref. [38] for all values of λ, especially for the higher-order LSUBn approximations with n > ∼ 6. Although the results for M z in Fig. 2 also compare well, by inspection, with the exact results of Ref. [38] in the region where M z is known to be non-zero from these exact calculations, (i.e., λ < 1 2 .), the agreement is now much poorer outside this region (i.e., λ > 1 2 ) for any of the LSUBn approximants shown. However, even in this case, the agreement is found to become excellent when the LSUBn sequence of approximants is extrapolated to the exact limit, n → ∞, as alluded to in Sec. II A. Thus, a very well-tested extrapolation scheme for use in such cases where the system undergoes a quantum phase transition (see, e.g., Ref. [36] and references cited therein) is where M z (n) is the nth-order CCM approximant (i.e., at the LSUBn or SUBn-n level) to M z . Thus, in Fig. 2 we also show the extrapolation using Eq. (22) as the fitting formula, together with the LSUBn approximants (for both the even values of n shown and the un- shown odd values) with 6 ≤ n ≤ 12 as the input data, to determine the extrapolated value µ 0 , which is plotted. Clearly, the extrapolation now agrees extremely well with the exact result, even in the very sensitive region very close to the critical value λ c , the value of which itself is now also predicted rather accurately.

IV. SPIN-HALF TRIANGULAR-LATTICE HEISENBERG ANTIFERROMAG-NET IN AN EXTERNAL MAGNETIC FIELD
We now consider the spin-1 2 triangular-lattice Heisenberg antiferromagnet in a magnetic field. The Hamiltonian that we will use here is given bŷ where the index i runs over all N lattice sites on the triangular lattice and the sum over the index i, j indicates a sum over all nearest-neighbour pairs, with each pair being counted once and once only. The strength of the applied external magnetic field is again given by λ. The triangular lattice is itself tripartite, being composed of three triangular sublattices, denoted as A, B and C, the sites of which we denote respectively as A n , B n , and C n . If the original lattice has a distance a between nearest-neighbour sites, the corresponding distance on each of the sublattices A, B and C is √ 3a.
It is easy to see that the classical spin-S model corresponding to Eq. (23) (see, e.g., Ref. [39]) has an infinitely (and continuously) degenerate family of ground states, red with the associated order parameter space being isomorphic to the 3D rotation group SO(3).
Thus, one may readily rewrite the classical energy per spin for this model in the form where λ = λẑ andẑ is a unit vector in the z-direction, and is defined to be the sum of the three spins on the kth elementary triangular plaquette on the lattice with nearest-neighbour vertices A ∆ k , B ∆ k , and C ∆ k . Equation (24) shows clearly that the energy is minimised when each of the squared terms in the sum over elementary triangular plaquettes is either zero (which is possible for λ ≤ 9S) or minimised [viz., to take the value (3S − 1 3 λ) 2 for λ > 9S]. Thus, we find rather simply that the classical ground-state energy per spin is given by For λ ≤ 9S, the ground state is clearly infinitely (and continuously) degenerate, since any configuration of spins that satisfies S ∆ k = 1 3 λ on all 2N elementary triangular plaquettes will yield the same energy. Furthermore, this condition immediately yields that the comparable classical value for the lattice magnetisation M, where M ≡ N i=1 S i = Mẑ (i.e., in the direction of the field), in the ground state is given by from which we also see that the magnetisation saturates at the value λ = λ s = 9S of the magnetic field strength.
In the zero-field case (λ = 0) the energy-minimising condition (viz., that S ∆ k = 0 on all 2N elementary triangular plaquettes) simply becomes the condition for the usual 120 degeneracy as that of the latter, due to the condition that S ∆ k = 1 3 λ on all 2N elementary triangular plaquettes. Thus, on each plaquette, each of the three spins has two orientational degrees of freedom, and the above condition simply reduces the overall degrees of freedom from six to three. The trivial degeneracy of the λ = 0 case is now, however, quite non-trivial in the λ = 0 case, since the local 120 • triangular-plaquette structures can become quite deformed by the application of the external field, even into non-coplanar configurations, as we now discuss.
From our discussion above, in principle, depending on the magnetic field strength, any of the five ground-state spin configurations sketched in Fig. 4 may appear. While the states I, II, III and IV are coplanar states, the "umbrella" state V is a 3D non-coplanar state. Although on the classical level both coplanar and non-coplanar sates are energetically degenerate, as we have noted above for λ ≤ λ s , thermal fluctuations tend to favour the coplanar configurations [39][40][41][42].
By minimising the energy, it is easy to show that the classical spin-S model described by Eq. (23) has a (coplanar) ground state of type I in Fig. 4 for λ < 3S, with a canting angle α given by At zero field (λ = 0) state I simply becomes the usual 120 • three-sublattice Néel state, while precisely at the value λ = 3S the state I becomes the collinear state II shown in Fig. 4, and as λ is increased further the ground state now smoothly transforms into state III shown in Fig. 4. The canting angles α and β are found to be given by and which may readily be shown to satisfy the condition, 2 cos α = cos β, which ensures that state III does not acquire any lattice magnetisation transverse to the applied field. When the field strength takes the value λ = 3S the angles are α = 1 2 π and β = − 1 2 π, which is again just equivalent to state II. As λ is then increased, up to the saturation value λ = λ s = 9S, the angle α first decreases to its minimum value, α = 1 3 π, at λ = 3 √ 3S, after which it again increases smoothly back to the value α = 1 2 π at λ = 9S. At the same time, as λ is increased beyond the value 3S, the angle β increases from − 1 2 π to 1 2 π at λ = 9S, taking the value β = 0 in between, precisely at the point λ = 3 √ 3S where α becomes a minimum. For all values λ > λ s = 9S the ground state is the fully saturated ferromagnetic state (viz., state III with α = 1 2 π = β). One may readily show that the classical ground-state energy of Eq. (23), for both states I and III at the respective values of their minimising canting angles and for the fully saturated state, is just that given previously in Eq. (25). Furthermore, the corresponding classical value for the lattice magnetisation (i.e., in the direction of the field) in the ground states I and III is given by our previous result of Eq. (26).
Although the state IV shown in Fig. 4 is not utilised as a CCM model state in any further application in this Section to the spin-1 2 case of the Hamiltonian of Eq. (23), it will be considered later in Sec. V. Hence, for completeness, we note that state IV has a minimum energy for the classical spin-S case of the present Hamiltonian for a value of the canting angle α shown in Fig. 4 given by Hence, unlike the situation for state I, which undergoes a smooth transformation to state III at a value, λ = 1 3 λ s , of the external field strength, (which then itself smoothly varies as λ is further increased up to the value λ s , at which point it becomes the fully saturated ferromagnetic state), state IV simply varies smoothly from the 120 • three-sublattice Néel state at zero field, λ = 0, to the fully saturated ferromagnetic state at λ = λ s .
For the classical case, as we have already noted, a non-coplanar state of the form of state V of Fig. 4 is degenerate in energy with states I and III above in their respective regimes. Thus, one readily finds that state V has a minimum energy for the classical spin-S Hamiltonian of Eq. (23) for a value of the out-of-plane angle θ given by Thus, with that value of θ, state V also yields a value for the energy identical to that of Eq. (25). Clearly, the lattice magnetisation is then also given by Eq. (26).
According to the CCM scheme briefly outlined in Sec. II we have to perform a passive rotation of the local spin axes of the spins such that all spins appear to point downwards for all five model states I, II, III, IV and V in Fig. 4. For the coplanar model states I, II, and III that procedure has been explained and discussed in detail in Ref. [24]. A similar rotation is also necessary for the coplanar model state IV, that does not play a role in this section, but which will be used in Sec. V. For the non-coplanar model state V, that comprises spins that make an angle θ to the plane perpendicular to the external field, the derivation of the Hamiltonian after rotation of the local spin axes is given in Appendix B, from which we note that the final result is given bŷ where the sums over i B,C,A → j C,A,B represent a shorthand notation to include the three sorts of "directed" nearest-neighbour bonds on each basic triangular plaquette of side a on the triangular lattice, which join sites i B and j C going from the B-sublattice to the C-sublattice, sites i C and j A going from the C-sublattice to the A-sublattice, and sites i A and j B going from the A-sublattice to the B-sublattice (in those directions only and not reversed). We see that this Hamiltonian now contains terms with both real and imaginary coefficients.
Clearly, when using any classical configuration of spins as a CCM model state, such as those shown in Fig. 4, there is no reason to expect that the quantum spin-S version of the model, with a finite value of the spin quantum number S, will take the same values of the angle parameters that characterise it as the classical version (i.e., in the S → ∞ limit), even in the case that the quantum ground state (at least partially) preserves the classical ordering inherent in the model state. For this reason a first step in using any such model state in a CCM calculation is to optimise the angle parameters that characterise the spin configuration.
To do so we simply choose those parameters that minimise the ground-state energy at each (either LSUBn or SUBn-n) level of approximation that we undertake, performing a separate such optimisation at each level.
Typical such CCM results for the ground-state energy of the spin-1 2 Hamiltonian described by Eq. (23) from using the non-coplanar state V as the model state are shown in Fig. 5 as a function of the out-of-plane angle θ and for various values of the external field strength in the range 0 ≤ λ ≤ λ s , for the particular case of the LSUB5 level of approximation. The groundstate energy is found to be a real number for all values of λ and θ. We note in particular that all purely imaginary contributions to the energy sum identically to zero. Figure 5 demonstrates the general result that the ground-state energy has a well-defined minimum with respect to the angle for all values of λ, at each LSUBn level of approximation. The angle that minimises the ground-state energy is plotted as a function of λ in Fig. 6  External Field, λ is coplanar) when the external field is zero (λ = 0). Also as required, all spins point in the direction of the field (i.e., shown by θ/π = 1 2 ) at the saturation field, λ s = 9 2 . The angle that minimises the energy varies continuously as a function of λ for model state V, excepting the limiting point at "saturation", λ s . Furthermore, we see from Fig. 6 that the non-coplanar energy-minimising configuration of spins converges very rapidly as the LSUBn approximation index n is increased, for all values of λ.
Ground-state energies for model state V are shown in Fig. 7, in which results for the coplanar model states I to III from Ref. [24] are also shown for comparison. Again, LSUBn results for the energy converge rapidly with increasing levels of the truncation index n for all values of λ. We see very clearly that ground-state energies for the coplanar model states lie lower than those of the 3D "umbrella" state (model state V) for all values λ, which is in agreement with the results of other methods [39][40][41]. Note that the energy difference between the coplanar and the non-coplanar states is particularly large in the plateau region around λ = 1.5, as can be seen clearly from the inset in Fig. 7.
Naturally, the (physical) lattice magnetisation is defined in terms of the spin directions before all rotations of the local spin axes have been carried out. Thus, the lattice magnetisation is given in terms of the "unrotated" coordinates as: After the rotations of the local spin axes for model state V, which led to the expression of Eq. (32), have been completed, this expression is given by Previous initial results for model state V [25] used computational differentiation and the Hellmann-Feynman theorem to evaluate this lattice magnetisation of Eq. (33). Here we evaluate it directly by finding both the ket-and bra-state correlation coefficients, which are now complex-valued, and then evaluating the expectation value explicitly, although we note that this method provides identical results (within the precision allowed by numerical differentiation) to that of the former technique, as required.
Once again, the results for the lattice magnetisation shown in Fig. 8   have now shown explicitly lie lower in energy than the 3D "umbrella" state, are in excellent agreement with experimental results for the magnetic compound Ba 3 CoSb 2 O 9 (a spin-half triangular-lattice antiferromagnet) and exact diagonalisations [61]. Note also that previous CCM results for the coplanar model states also indicate that a similar plateau occurs over the range 2.82 < ∼ λ < ∼ 3.70 for the for the spin-one triangular-lattice antiferromagnet, and this theoretical result has subsequently been established experimentally for the compound Ba 3 NiSb 2 O 9 (a spin-one triangular-lattice antiferromagnet) [30].

V. SPIN-S TRIANGULAR-LATTICE XXZ ANTIFERROMAGNET IN AN EX-TERNAL MAGNETIC FIELD
In recent investigations [62][63][64][65][66] of the anisotropic triangular-lattice XXZ model where the indices have the same meaning as in Eq. (23), it has been shown that for an easy-plane anisotropy (i.e., ∆ < 1) the 3D "umbrella" state V discussed in Sec. IV can become energetically favoured over the coplanar states, so as to form the true ground state under certain conditions that we now elaborate.
The corresponding phase diagram in the ∆-λ plane is rich, (see, e.g., Ref. [62]). Moreover, the phase boundary between the coplanar and non-coplanar ground states strongly depends on the spin quantum number S. We note that in the classical limit S → ∞ for ∆ < 1 the non-coplanar "umbrella" state is always energetically favoured over the planar states to form the ground state, as we elaborate further below, whereas for the extreme quantum case S = 1 2 there is a wide region of values of the anisotropy parameter ∆ and the field strength λ where coplanar states are favoured. We note further that since the energy differences between competing ground states can be very small, accurate and self-consistent calculations are hence required to be able to distinguish between them reliably.
By comparison with the derivation of Eq. (24) for the Hamiltonian of Eq. (23) in Sec. IV, it is clear that we may write the classical energy per spin for the current anisotropic triangularlattice XXZ model in the form It is evident that the second sum in Eq. (36) can now potentially favour non-coplanar states in the case of easy-plane anisotropy (i.e., when ∆ < 1). One readily finds, by making use of Eq. (36), that state V of the form shown in Fig. 4 has a minimum energy for the classical spin-S Hamiltonian of Eq. (35) for a value of the out-of-plane angle θ given by where λ s ≡ 3(1 + 2∆)S is the value of the field strength that reaches saturation (i.e., the fully aligned ferromagnetic state) for this state V. With this value of θ one may readily show that state V yields a value for the classical ground-state energy per spin given by which also replicates Eq. (25) at isotropy (i.e., when ∆ = 1). Furthermore, one may readily show that state V, with the energy-minimising value of the out-of-plane angle θ of Eq. (37), yields a classical value for the lattice magnetisation given by We may compare the above results for the non-coplanar state V with those of the coplanar states. For example, one may readily show, again by making use of Eq. (36), that state IV of the form shown in Fig. 4 has a minimum energy for the classical spin-S Hamiltonian of Eq. (35) for a value of the canting angle α given by where λ s = 3(1 + 2∆)S as before. Once again, this result is in accord with our previous result of Eq. (30) at the isotropic point, ∆ = 1. With this value of α one may show that state IV yields a value for the classical ground-state energy per spin given by which also replicates Eq. (25) at isotropy (i.e., when ∆ = 1). One may also show that state IV, with the energy-minimising value of the canting angle α of Eq. (40) yields a classical value for the lattice magnetisation given by which may be compared with the corresponding result of Eq. (39) for state V. Finally, at the classical level, one may readily show from Eqs. (38) and (41) that the difference between the minimum energy per spin for the "umbrella" state V and that for the coplanar state IV is given explicitly by It is evident from Eq. (43) that the "umbrella' state V always lies lower in energy than the coplanar state IV for all values of the anisotropy parameter ∆ < 1, as we have already asserted, in the classical limit S → ∞. Equation (43) shows clearly, however, that the energy difference between the states decreases both as λ approaches the saturation value λ s and as ∆ approaches unity (i.e., near the Heisenberg isotropic limit). One expects on rather general grounds that quantum fluctuations will favour phases with coplanar over noncoplanar configurations of spins. One also expects that the effects of quantum fluctuations will increase monotonically as the spin quantum number S is decreased smoothly from the large-S classical limit. Hence, for small positive values of (λ s − λ) it is, a priori, likely that the value ∆ 1 (S) of the anisotropy parameter at which a possible quantum phase transition occurs between the "umbrella" state forming the ground state (for ∆ < ∆ 1 ) to a coplanar state forming the ground state (for ∆ > ∆ 1 ) would decrease smoothly as S decreases from the classical value ∆ 1 (∞) = 1 at S → ∞ towards a lower value ∆( 1 2 ) at the extreme quantum limit, S = 1 2 . Thus, to demonstrate the capability of the new 3D CCM code to detect such anisotropydriven quantum transitions between coplanar and non-coplanar states we now consider the Other investigations [62][63][64][65][66] have shown that for strong magnetic fields with a strength λ infinitesimally below the saturation field strength λ s the states III, IV or V shown in Fig. 4 appear as the ground state, depending on the value of the anisotropy parameter ∆. In the XY limit, ∆ = 0, the 3D "umbrella" state V is always the ground state, for all values of the -2x10 -6 -  Therefore, we focus particular attention here on the transition between the 3D "umbrella" state V and the coplanar state IV, which we have already discussed above in the classical limit, S → ∞.
In Fig. 9 we in the lower panel of Fig. 9), does not actually indicate a reentrance of the "umbrella" state V, since for decreasing λ the coplanar state III, which we have not considered here, actually becomes the ground state (see, e.g., Refs. [62,63,65,66]). Based on curves such as those shown in Fig. 9 we derive the critical points ∆ 1 (S) with an accuracy of ±0.015, as shown in Fig. 10. We compare our CCM data with the large-S approach of Starykh et al. [63] which yields ∆ 1 (S) = 1 − 0.53/S, as well as with numerical data obtained by the dilute Bose gas expansion [65]. As can be clearly seen, there is good agreement of our results with both those of Ref. [65] and those of Ref. [63], except where the latter results from the large-S expansion naturally fail for the smaller values of S. Finally, we note too that our results also agree extremely well with those from a recent study that used the numerical cluster mean-field plus scaling method [66] to investigate the cases with S ≤ 3 2 . To conclude, we have shown that the new CCM code for 3D non-coplanar ground states provides accurate data that, in particular, allow us to examine with confidence the quantum selection of competing ground states of the frustrated triangular-lattice XXZ antiferromagnet in an external magnetic field.

VI. SUMMARY
Coplanar model states for applications of the CCM to problems in quantum magnetism are those states in which all spins lie in a plane, whereas 3D model states are non-coplanar states in which the spins do not lie in any plane. The first step in applying the CCM to lattice quantum spin systems is always to rotate the local spin axes (i.e., on each lattice site) so that all spins in our model state appear mathematically to points in the downwards (i.e., negative) z-direction. We have shown explicitly here that this process leads inevitably to terms in the rotated Hamiltonian with complex-valued coefficients for 3D model states, by contrast with the case for coplanar model states, where all of the respective terms carry (or can be made to carry) real-valued coefficients. Since these rotations represent unitary transformations the rotated Hamiltonians in each case are still Hermitian. Nevertheless, for the case of 3D model states, even though the expectation values of all physical operators are hence guaranteed to be real-valued quantities, the intervening CCM bra-and ket-state multispin cluster coefficients from which they will be calculated will be complex-valued at all levels of approximation.
In this paper we have demonstrated that the existing high-order CCM code [37] can be extended appropriately to be able to be applied to such rotated Hamiltonians with complexvalued coefficients that arise from the utilisation of 3D model states. An explicit derivation of such a Hamiltonian after all rotations of spin axes for a 3D model state was given for the triangular-lattice Heisenberg antiferromagnet in an external magnetic field. Although in such cases the CCM ket-and bra-state equations for the multispin cluster correlation coefficients are now complex-valued quantities, we have demonstrated explicitly for several models of interest that all expectation values are real numbers. This was also shown explicitly in analytical LSUB1 calculations for the one-dimensional spin- 1 2 Ising ferromagnet in a transverse external magnetic field.
Due to the length and complexity of the CCM code, its extension from coplanar to 3D model states is a non-trivial task. Furthermore, the task of defining the problem to be solved by the code in CCM "script files" becomes more difficult. Hence, this is an important advance for practical applications of the CCM, and one that greatly extends its range of applicability.
Finally, excellent correspondence with the results of other methods was found for all of the cases considered here, namely, (a) the 1D spin- 1 2 Ising ferromagnetic chain in a transverse external magnetic field; (b) the 2D spin-1 2 Heisenberg antiferromagnet on a triangular lattice in the presence of an external magnetic field; and (c) the 2D spin-S triangular-lattice XXZ antiferromagnet in the presence of an external magnetic field, for the cases 1 2 ≤ S ≤ 5. For the first case of the transverse Ising model, which simply provides a testbed for which we can artificially introduce terms in a rotated Hamiltonian pertaining to it that contain an imaginary coefficient, we showed that CCM results agree well with exact results. With the extended CCM methodology thereby validated we turned to two frustrated 2D models defined on a triangular lattice, for both of which a real 3D non-coplanar configuration of spins is physically relevant.
Our CCM results for the spin-1 2 triangular-lattice Heisenberg antiferromagnet show that coplanar ordering is favoured over non-coplanar ordering for all values of the applied external magnetic field, which agrees with the results of other approximate methods [39][40][41]. The boundary between a 3D ground state (state V) and a coplanar ground state (state IV) was also obtained for the XXZ model on the triangular lattice for values of the external magnetic field near to saturation and for spin quantum number S ≤ 5. The CCM calculations in this case were computationally intensive for this frustrated model, especially for high spin quantum numbers. However, we note that only a very few other approximate methods can deal effectively and accurately with the combination of strong frustration and higher spin quantum number. The differences in ground-state energies between the two states is also very small in the limit of field saturation. Hence, an accurate delineation of the phase boundary is an extremely delicate task. Despite this inherent difficulty excellent correspondence was seen with the results of other approximate methods. These results thereby constitute a useful advance in the understanding of this model, as well as providing an excellent quantitative test of high-order CCM using 3D model states. We first carry out a CCM LSUB1 calculation for the spin- 1 2 1D transverse Ising model of Eq. (13). We recall that the CCM model state contains spins that point in the downward z-direction. We note also that this is the exact (albeit trivial) ground state when λ = 0 and that all CCM correlation coefficients should therefore tend to zero in this limit, λ → 0. The LSUB1 ket-state operatorŜ is thus given bŷ By making use of the usual SU(2) spin commutation relations, [ŝ + l ,ŝ − l ′ ] = 2ŝ z l δ l,l ′ and [ŝ z l ,ŝ ± l ′ ] = ±ŝ ± l δ l,l ′ , and by also using the nested commutator expansion of Eq. (10), it is readily proven that the CCM similarity transforms of the operatorsŝ + l ,ŝ z l , andŝ − l are given by e −Ŝŝ+ l eŜ =ŝ + l ; e −Ŝŝz l eŜ =ŝ z l + aŝ + l ; (A2) at the LSUB1 level of approximation. Thus we see that the corresponding CCM LSUB1 similarity transform of the Hamiltonian of Eq. (13) is given by By making use of the relation,ŝ z |Φ = − 1 2 |Φ , for the present model state |Φ , the groundstate energy is now readily evaluated from Eq. (9) as Furthermore the ground-state LSUB1 equation for the coefficient a is given from Eq. (7) as follows, This quadratic equation has the physical solution where we have discarded the (unphysical) solution with the negative sign of the square root since a must be zero when λ = 0, as noted above. Using Eqs. (A4) and (A6), the ground-state energy is thus given by at the CCM LSUB1 level of approximation. As tests of this equation, we note that E g /N = −1/4 when λ = 0 and that E g /N → −λ/2 as λ → ∞, which are the correct results in these two limiting cases. For comparison both with the exact result of Eq. (17) and with those from higher-order CCM LSUBn approximations with n > 1, the LSUB1 result of Eq. (A7) is also shown in Fig. 1.
At the same CCM LSUB1 level of approximation the bra-stateŜ operator is given bŷ such that we may now evaluate the LSUB1 ground-state energy expectation value functional, H, asH We see immediately that which simply gives us the LSUB1 ket-state equation of Eq. (A5) again, as required. Furthermore, we see that the LSUB1 bra-state coefficientã can similarly now also be obtained as follows, The magnetisation in the Ising direction, can now easily be evaluated at the LSUB1 level of approximation as The transverse magnetisation is given by which is readily evaluated at the LSUB1 level of approximation as Once again, for comparison both with the corresponding exact results of Eqs. (18) and (19) and with those from higher-order CCM LSUBn approximations with n > 1, the LSUB1 results of Eqs. (A13) and (A16) are also shown in Figs. 2 and 3, respectively.
We now carry out a CCM LSUB1 calculation for the 1D transverse Ising model of Eq. We wish to compare LSUB1 results for this new Hamiltonian to those for the "unrotated" case with Hamiltonian given by Eq. (13). Indeed, we show explicitly that macroscopic quantities for this Hamiltonian do not change (and remain real), as they must, for the LSUB1 approximation even though the CCM correlation coefficients may now be complexvalued (i.e., contain real and imaginary components).
The LSUB1 approximation is again given by Eq. (A1) and the similarity transformed spin operators are given as before in Eq. (A3). Thus we see that The ground-state energy equation is now given by and the ground-state energy LSUB1 equation is given by Φ|ŝ − l e −ŜĤ eŜ|Φ = 0 ⇒ a + iλ 2 + iλa 2 2 = 0 .
Again, this quadratic equation has the physical solution given by where we have now discarded the (unphysical) solution with the positive sign of the square root since a must be zero when λ = 0, as noted previously. Thus, we see that the new LSUB1 ket-state correlation coefficient a is now complex-valued (indeed, pure imaginary) for all values of λ(> 0). Substitution of this value for a from Eq. (A20) into Eq. (A18) immediately yields that the ground-state energy at the CCM LSUB1 level of approximation is given by which is indeed seen to be identical to that given by Eq. (A7). The bra-state S operator is again given by Eq. (A8), such that we may explicitly evaluate the ground-state energy functional, We see immediately that again as required. Furthermore, we see that The magnetisation in the Ising direction is again given by Eq. (A12), which at the LSUB1 level of approximation (now using Eqs. (A20) and (A24)) is given by which is now explicitly real and in agreement with Eq. (A13), both as required. The transverse magnetisation should now, of course, be evaluated with respect to the y-direction (in terms of spin coordinates after rotation of the local spin axes) and it is given by Furthermore, in this step, we do not rotate the axes for spins at sites i C on the C-sublattice, since they are already pointing in the downwards direction: We may now use Eqs. (B4)-(B6) to rewrite the nearest-neighbour part of the Hamiltonian of Eq. (B3) that connects sites i B on the B-sublattice and j C on the C-sublattice in the direction going from site i B to site j C in the (unitarily equivalent) rotated form