Modified spin-wave theory with ordering vector optimization I: frustrated bosons on the spatially anisotropic triangular lattice

We investigate a system of frustrated hardcore bosons, modeled by an XY antiferromagnet on the spatially anisotropic triangular lattice, using Takahashi's modified spin-wave (MSW) theory. In particular we implement ordering vector optimization on the ordered reference state of MSW theory, which leads to significant improvement of the theory and accounts for quantum corrections to the classically ordered state. The MSW results at zero temperature compare favorably to exact diagonalization (ED) and projected entangled-pair state (PEPS) calculations. The resulting zero-temperature phase diagram includes a 1D quasi-ordered phase, a 2D Neel ordered phase, and a 2D spiraling ordered phase. We have strong indications that the various ordered or quasi-ordered phases are separated by spin-liquid phases with short-range correlations, in analogy to what has been predicted for the Heisenberg model on the same lattice. Within MSW theory we also explore the finite-temperature phase diagram. We find that the zero-temperature long-range-ordered phases turn into quasi-ordered phases (up to a Berezinskii-Kosterlitz-Thouless temperature), while zero-temperature quasi-ordered phases become short-range correlated at finite temperature. These results show that modified spin-wave theory is very well suited for describing ordered and quasi-ordered phases of frustrated XY spins (or, equivalently, of frustrated lattice bosons) both at zero and finite temperatures. While MSW theory, just as other theoretical methods, cannot describe spin-liquid phases, its breakdown provides a fast method for singling out Hamiltonians which may feature these intriguing quantum phases. We thus suggest a tool for guiding our search for interesting systems whose properties are necessarily studied with a physical quantum simulator.


I. INTRODUCTION
Lattice models of strongly interacting bosons have recently been implemented experimentally thanks to impressive developments with trapped ultra-cold atoms in optical lattice potentials [1,2]. A particularly appealing perspective in this arena is the study of strongly correlated lattice boson models with frustration, arising for instance from the coupling of bosons to a (artificial) magnetic field [3,4,5], or from a periodical shaking of the optical lattice [6]. Frustration in the intersite hopping amplitudes is formally equivalent to a description of the system in a rotating reference frame, which implies that the system is subject to the spontaneous appearance of vortices.
Such vortices can form ordered arrays (vortex crystals) coexisting with Bose condensation, which consequently takes place in a macroscopic wavefunction sustaining persisting circulating currents (see Ref. [7] and references therein); or they can even disrupt condensation completely, and lead to a disordered insulating state [8]. Such disordered states are notoriously difficult to study theoretically.
In the particular limit of a very strong interparticle repulsion, frustrated bosonic models can be exactly mapped onto S = 1/2 frustrated XY antiferromagnets [9]. In two dimensions these models are known to exhibit ground states with spiral order, representing the magnetic counterpart to the aforementioned Bose-condensed states with vortex arrays. In special circumstances the interplay between quantum fluctuations and frustration may lead to disordered spin-liquid states, which are in one-to-one correspondence with bosonic insulating phases. XY antiferromagnets can also be regarded as the limiting case of antiferromagnetic Hamiltonians with planar anisotropy in the couplings, relevant to the description of frustrated antiferromagnetic materials, and they can describe the physics of Cooper pairs in arrays of ultra-small Josephson junctions with magnetic frustration [10]. More recently we have proposed that frustrated XY antiferromagnets can be experimentally implemented by loading planarly trapped ions into an optical lattice [11]. useful quantum simulation, possibly outperforming any classical computation (see e.g. Ref. [12] and references therein). Indeed, as mentioned above, solving frustrated bosonic models amounts to solving a large class of frustrated antiferromagnets. Fundamental steps have very recently been taken experimentally towards the implementation of artificial magnetic fields in cold atom experiments via Raman schemes [13]. Hence exciting progress in this field is expected in the near future.
However, the difficulty of finding accurate theoretical descriptions of disordered quantum lattice models makes it hard to tell a priori which systems will present such interesting phases in an experiment. The most attractive aspect of quantum simulators is their potential ability to emulate model Hamiltonians whose phase diagram cannot be accurately predicted with current theoretical approaches. Therefore it would be highly desirable to dispose of a fast tool that can outline quantum-mechanical phase diagrams, point out disordered phases, and thus classify model Hamiltonians according to their potential interest for experimental quantum simulation. We propose that the methods presented here will serve this very purpose.
Planar systems of bosons in optical lattices, coupled to an artificial magnetic field, can be described by the Bose-Hubbard Hamiltonian where b i , b † i are bosonic operators, n i = b † i b i , and i, j represents pairs of nearest neighbor sites. The hopping amplitude t ij = −t ij exp (iA ij ), wheret ij ≥ 0, is spatially modulated by the line integral of the vector potential along the i, j bond, A ij = r i r j A(r) · dl, as well as by possible spatial anisotropies in the optical lattice (contained in the bare hopping amplitudest ij > 0). In the following we will focus on the limit of infinite repulsion U → ∞ and half filling n i = 1/2, under which the Bose-Hubbard model maps onto the S = 1/2 XY Hamiltonian [11]: where S α i is the α th component of the S = 1/2 spin operator acting on site i. In this work we focus on a triangular lattice with antiferromagnetic nearest-neighbor interactions, which can be seen as a ferromagnetic lattice (t ij > 0) with half a magnetic flux quantum threaded through each lattice plaquette. This magnetic flux can, for instance, be interpreted as flipping the signs of all hopping amplitudes of the bonds along the horizontal direction of the lattice ( Fig. 1). In this gauge the hopping is transformed to t ij =t ij for r i − r j = ±aτ 1 , where a is the lattice spacing and τ 1 = (1, 0), while t ij = −t ij for the other (diagonal) nearest-neighbor bond directions. A canonical transformation rotating the spins in odd rows by an angle π around the z axis flips the signs of all diagonal bonds, leading to a model of a spatially anisotropic triangular XY antiferromagnet with all t ij > 0. Such a model exhibits strong frustration on the triangular lattice. In the following we specialize to the case in which the hopping amplitudes t ij take two values, t 1 and t 2 depending on the orientation of the bond [parallel to τ 1 or along the diagonals ±aτ 2 and ±a (τ 2 − τ 1 ), with τ 2 = 1/2, √ 3/2 -see Fig. 1]. This particular case is relevant for the physics of a number of systems: it can describe neutral atoms trapped in triangular lattices formed with three lasers intersecting at 120 • angles [14], and one of which has a different intensity than the other two; another implementation of this model are neutral atoms in an isotropic triangular lattice in the absence of artificial magnetic fields but with elliptical time-dependent forcing of the lattice [6]; furthermore, it is relevant for triangular Wigner crystals of ions trapped in the minima of an optical lattice [11].
The aim of this work is the determination of the ground-state and finite-temperature phase diagrams of the S = 1/2 XY antiferromagnet (AF) -or, alternatively, of frustrated half-filled hardcore bosons -on a spatially anisotropic triangular lattice (SATL) by means of spin-wave theory. In particular, we show that Takahashi's modified spin-wave (MSW) theory [15] gives an adequate description of the main features of the ground-state and low-temperature phase diagram, while keeping the computational effort at a minimum. Spin-wave methods generally account for weak quantum fluctuations around the ordered state corresponding to the ground state in the classical limit S → ∞. We show how MSW theory can be extended to arbitrary reference states, whose ordering vector may be shifted with respect to the ground state in the classical limit [16].
Furthermore, we demonstrate how this procedure allows for a convenient calculation of the spin stiffness. In the specific case of the S = 1/2 XYAF on a SATL, the spin stiffness proves to be an effective tool supporting our search for spin-liquid phases. Supplementary results derived by projected entangled-pair states (PEPS) [17] and exact diagonalizations (ED) using the Lanczos method allow us to validate the zero-temperature results of the MSW method with ordering vector optimization. We moreover provide the finite-temperature phase diagram of the S = 1/2 XYAS on a SATL, and find a region where the breakdown of MSW theory indicates a disordered state.
Such a disordered phase is an ideal candidate for performing a useful quantum simulation because theoretical tools for studying it adequately are currently lacking. In an upcoming publication we will extend this method to the ground state phase diagrams of the Heisenberg SATL and the J 1 J 2 J 3 model [18].
For coherence with the theoretical technique used -spin-wave theory -and for a better comparison with existing results in the literature, the results of this paper will be generally expressed in the language of spin physics, but guidance will be provided on how to translate the magnetic observables into bosonic observables.

II. MODIFIED SPIN-WAVE FORMALISM
In the past MSW theory has been found to give a satisfactory qualitative account of the lowtemperature properties of spin systems, even frustrated or disordered ones. In this section we review its formalism, mainly following Xu and Ting [16] in the first steps, but considering XY interactions rather than Heisenberg interactions (for the latter see also Ref. [18]). This requires only minor modifications of the formulas, and it is expected that the validity of the spin-wave approach is even better justified for XY interactions since in this case the influence of quantum fluctuations is reduced by the anisotropy in spin coupling.
Our aim is to determine the phase diagram of the Hamiltonian of Eq. (2). A fundamental assumption of spin-wave theory as applied to the XY model is that the ground state shows longrange order (LRO) with the spins classically lying in the xy-plane; for a translationally invariant system, like the one under investigation, the ordered ground state is characterized by a well defined ordering vector Q. Under this assumption it is convenient to rotate the local reference system of each spin from (x, y, z) to (η, ζ, ξ) so that the ground state in the local reference frame has all spins aligned in the same direction. This amounts to the following transformation: Then S ζ i , which will be the quantization axis, lies parallel to the classical spin S i = (cos (Q · r i ) , sin (Q · r i ) , 0). The component S η i lies perpendicular to it in the xy plane, and S ξ i is perpendicular to the xy plane. Unlike in ordinary spin-wave theories we do not make any assumption on the ordering vector Q. In particular, it may well differ from the one exhibited by the system in the classical limit (Q cl ). Spin waves can be described by applying the Dyson-Maleev (DM) transformation [19,20], which maps the physical spins to interacting bosons, where S ± i ≡ S ξ i ± iS η i . The DM transformation is an exact mapping between spins and bosons as long as projectors are retained which keep the system in the physical subspace, i.e. the subspace where at each site only 2S DM bosons are present at most. It can be shown that these projectors have the form P = 1 + O[n/(2S)] 3 where n is the DM boson density [21]. Hence, under the assumption of diluteness of the DM boson gas, n/(2S) < 1 (in fact n = S, see below), we can safely neglect the P projectors.
If the spin Hamiltonian under investigation is obtained as the hardcore limit of the Bose-Hubbard Hamiltonian Eq. (1), it is important to distinguish the DM bosons from the physical b-bosons from which the effective spin Hamiltonian originated. Indeed the DM bosons at a site i quantify the deviation of the i th spin from the local direction in the xy plane set by the ordered structure with ordering vector Q. On the other hand, the physical bosons correspond to the spin deviations with respect to full alignement of the spin along the z axis. The particle-hole symmetry of the bosonic Hamiltonian, Eq. (1), in the hardcore limit leads to half filling of the physical bosons, which accidentally coincides with the average filling imposed by the Takahashi's constraint on the DM bosons for S = 1/2 (see next subsection). Yet all other properties are in general quite different.
A. Derivation of a mean-field Hamiltonian and Takahashi's constraint Applying Eqs. (3) and (4) to Eq. (2) one arrives at the bosonic Hamiltonian Here we have defined the correlators These correlators can be rewritten in terms of independent particles by first Fourier transforming where N is the number of sites, and then applying a Bogoliubov transformation Requiring that the Bogoliubov particles be non-interacting imposes α k α k = α † k α † k = 0. This condition also removes the anti-Hermitian part of the Hamiltonian.
The correlators are now with n k = α † k α k = 1/ (exp (ω k /T ) − 1) being the occupation number of Bogoliubov mode k at temperature T (with the Boltzmann constant k B set to unity). The dispersion relation ω k is determined self-consistently in the next section.
So far we have essentially formulated a standard Hartree-Fock theory for the gas of interacting DM bosons. A very important modification to this theory, due to Takahashi [15], is the introduction of the constraint of zero magnetization at each site, The implementation of this constraint amounts to effectively reducing the Hilbert space dimension The correct spin wave description is found by minimizing the free energy is the entropy of a set of harmonic oscillators. Minimizing with respect to θ k and ω k under the constraint of Eq. (10) yields a set of self-consistent equations, with where µ is the Lagrange multiplier for Eq. (10) corresponding to the chemical potential for changing the total magnetization.
In Eqs. (13) we have abbreviated F ij = F (r ij ), and G ij = G (r ij ). Note that in the classical limit S → ∞ one gets G ij , F ij ≈ S and Eqs. (13) become analogous to their LSW counterparts.
The spin-wave spectrum reads Inserting Eq. (13) into Eq. (14) shows that a finite µ entails a gap at k = 0. This is to be seen in contrast to LSW theory where the spectrum always has a gapless Goldstone mode at k = 0. The correlators at the minimum take the form At T = 0 where n k = 0 ∀k = 0, one finds that µ vanishes. This implies also the disappearance of the gap at k = 0 that may exist for finite temperature. A vanishing gap is a necessary requirement for the appearance of the Goldstone mode associated with magnetic LRO. It also enables Bose condensation of the DM bosons in the k = 0 mode. Separating out the contribution of the zero mode, a † k=0 a k=0 /N = a k=0 a k=0 /N ≡ M 0 (corresponding to the order parameter measuring the total spiraling magnetization in the quantization axis directions given by the ordering vector Q), one arrives at the zero-temperature equations and the constraint Eq. (10) becomes As mentioned above, the occupation of the zero mode M 0 corresponds to a Bose-Einstein condensate of the DM bosons in the minimum of the dispersion relation. This condensate is depleted by interactions of the DM bosons. The larger this depletion, the more DM bosons reside at momenta different from zero, thus decreasing magnetic LRO.

C. Optimization of the ordering vector
It is not a priori clear that the classical ordering vector Q cl correctly describes the LRO in the quantum system. To account for the competition between states with LRO at different ordering vectors Q we extend the MSW procedure by optimizing the free energy with respect to the ordering vector Q. This procedure, first introduced in Ref. [16], significantly improves the predictions of MSW theory. It amounts to finding the best ordered reference state with in-plane ordering vector Q (spiral state) whose free energy is minimized not at the classical level, but including the effect of quantum fluctuations self-consistently within MSW theory.
The minimization of F with respect to Q x and Q y yields two additional equations which must be added to the set of self-consistent equations, The values of F ij and G ij can now be calculated by solving self-consistently Eq. (18) together with Eqs. (10,(12)(13)(14)(15). At zero temperature, Eq. (10) and Eq. (15) have to be replaced by Eq. (17) and Eq. (16), respectively. Through Wick's theorem the knowledge of the quantities F ij and G ij allows for the computation of the expectation value of any observable.

D. Spin stiffness
The optimization of the ordering vector allows for a straightforward calculation of the spin stiffness. This additional information, complementary to the order parameter, helps us in identifying candidate regions for spin-liquid behavior.
The MSW Ansatz always returns only a single one of all the possible ordering vectors Q 0 as the optimal ordering vector. However, if the true ground state is only short-range ordered, we might expect the Q-minimum to be relatively shallow, and that a slight change of the ordering vector barely affects the free energy F. This means that the order is not very stable against twists of the spin configuration. We quantify the curvature of the minimum by the spin stiffness tensor evaluated at the optimized ordering vector Q 0 . In particular we will determine the parallel spin stiffness ρ ≡ 1 2 Tr ρ = 1 2 (ρ xx + ρ yy ) and the Gaussian spin stiffness The spin stiffness gives a measure of how stiff magnetic LRO order is with respect to distortions of the ordering vector and it provides a fundamental self-consistency check of our approach. In fact, finding a small spin stiffness casts doubt on the reliability of the spin-wave approach in describing such a strongly fluctuating state, and hence suggests that the true ground state might be quantum disordered.
Since a change in Q affects the correlators F ij and G ij , we must compute Υ self-consistently.
After finding the optimal Q 0 by the self-consistent procedure described in the previous sections, we calculate 1 N F (Q x , Q y ) self-consistently for several ordering vectors Q = Q 0 + ∆Q and fit a quadratic form to the results. Since the minimum in the free energy can be very shallow, this procedure can be somewhat affected by numerical noise. As an approximation to the true spin stiffness, the partial spin stiffness ρ partial αβ can be computed via the partial derivatives, i.e. without recalculating the self-consistent equations. It reads We define Υ partial analogously to Υ [Eq. (20)] as the determinant of the partial spin-stiffness tensor.
The system can lower its energy by adjusting F ij and G ij to the new ordering vector, and therefore Υ is always smaller than Υ partial . We will see later that in some cases the partial Gaussian spin stiffness Υ partial gives a good estimate of the real Gaussian spin stiffness Υ, but there are cases where it is considerably larger.
In the following we present the zero-temperature phase diagram of the anisotropic triangular lattice. This paradigmatic example will show that the spin stiffness is an important quantity that provides a deeper insight into the order properties of the system.

E. From spins to bosons
As we mentioned in the introduction, in the S = 1/2 case the spin Hamiltonian is equivalent to that of infinitely repulsive bosons at half filling with frustrated hoppings. Hence it is important to match spin observables with their bosonic counterparts. Following the Dyson-Maleev or the Holstein-Primakoff transformation, spin operator and physical (hardcore) b bosons obey the rela- [43]. Here the b operators obey anticommutation rules on site, The ordering vector Q corresponds to the finite momentum at which condensation occurs. The condensed state in the spiral phase is characterized by a pattern of persistent currents forming a crystal of vortices, whose geometrically correlated structure is captured by the spin chirality (see below). Finally the parallel spin stiffness corresponds to the superfluid density of the bosons, ρ s = ρ /S.

III. GROUND STATE PHASE DIAGRAM OF THE ANISOTROPIC TRIANGULAR LATTICE
In this section we compute the ground-state phase diagram of the spatially anisotropic triangular lattice (SATL) with nearest-neighbor (NN) XY interactions. We consider a wide range of α ≡ t 2 /t 1 , where t 1 denotes the bond strengths along the chains and t 2 the bond strengths along the diagonals (black and red bonds, respectively, in Fig. 1). The parameter α interpolates between an ensemble of decoupled one-dimensional chains at α = 0, the isotropic triangular lattice at α = 1, and the square lattice for α → ∞.
If we assume the spins to behave classically, the 2D-Néel order, present for α ≥ 2, starts to continuously deform into spiral order at α ≤ 2 [compare Fig. 2 (a)]. The spiral phase extends down to α = 0 where the chains decouple. In a previous publication we presented the quantum mechanical phase diagram as predicted by projected entangled-pair states (PEPS) calculations [11]; it is reproduced for convenience in Fig. 2 (b). According to this, both the square lattice limit (α → ∞) and the most frustrated case, the isotropic triangular lattice (α = 1), display magnetic LRO. In the limit of uncoupled chains (α = 0) the system displays quasi-LRO with algebraically decaying correlations. However, similarly to what has been found in the Heisenberg model [24], the system seems to feature spin-liquid phases with exponentially decaying correlations between different types of order or quasi-order. In Appendix A we provide a further spectral feature, coming from the exact diagonalization of a small cluster, which is consistent with the observation coming from PEPS calculations. Further distinct features of the quantum model are that the transition between 2D-Néel and spiral order is shifted by quantum fluctuations to considerably smaller values of α, and that the quasi-ordered 1D-like state extends over a whole region of finite α in the phase diagram. phase diagram from Ref. [11]. SL is short for spin liquid.
In the following we compare predictions of MSW theory with these PEPS results and exact diagonalizations (ED). This will allow us to validate the reliability of the MSW method. The following system geometries are considered for the three different methods: • PEPS: a rhombic lattice of 20 × 20 = 400 spins with open boundary conditions and bond dimension D = 2. PEPS is a powerful numerical tool which goes beyond mean-field theory, but which, for small bond dimension D, only partially accounts for the entanglement in the ground state. This limitation becomes particularly serious close to quantum phase transitions. However, in Ref. [11] it was demonstrated that D = 2 is already accurate enough to effectively capture the physics of the system.
• ED: Lanczos diagonalization of clusters of 24 and 30 spins (the latter is shown in Fig. 3), again with open boundary conditions (necessary in order to allow for the accommodation of arbitrary ordering vectors).
• MSW: rhombic lattices of 32 × 32 = 1024 and 64 × 64 = 4096 spins and in the infinite-lattice (thermodynamic) limit, under periodic boundary conditions. We find that at these lattice sizes all quantities have essentially reached the infinite lattice limit in most of parameter space. As it can be expected, the deviations from the thermodynamic limit are sizable only in the one-dimensional limit and at critical points. The thermodynamic limit is computed by replacing finite sums over the first Brillouin zone with integrals.
In the triangular lattice with nearest-neighbor interactions one finds in the MSW formalism that Eq. (18b) gives Q y = 0 and from Eq. (18a) we obtain the formula Here τ 1 = (1, 0) and τ 2 = 1/2, √ 3/2 are the primitive lattice vectors (here and in the rest of the work we take the lattice spacing a equal to unity). For F ij = G ij = S, attained when S → ∞, this reduces to the ordering vector of the LSW theory Q cl x , Q cl y = (2 arccos (−α/2) , 0) which coincides with the classical ordering vector.

A. Breakdown regions for MSW theory
As a first step in our analysis we investigate the parameter regions where LRO is to be expected, and the locations where MSW theory ceases to be applicable. To this end we first investigate if there appear imaginary modes in the dispersion relation, which would indicate instabilities. Afterwards we study the order parameter M 0 and the spin stiffness.

Imaginary frequencies and breakdown of convergence
Convergence in the self-consistent equations of MSW theory with ordering vector optimization, Eqs. (12-14, 16, 17, 22), cannot be achieved in selected regions of the ground-state phase diagram, namely for α 0.18 and for 1.35 α 1.66, as summarized in Fig. 4. (Interestingly, convergence is restored in the pure 1D limit, α = 0, for which the theory formulates surprisingly good predictions.) This breakdown of convergence corresponds to the appearance of an imaginary part in the spinwave frequencies, Eq. (14), signaling an instability of the ordered ground state. The breakdown of a self-consistent description of the system in terms of an ordered ground state is strongly suggestive of the presence of a quantum-disordered ground state in the exact behavior of the system. Hence, one can interpret these parameter regions as candidates for the spin-liquid phases predicted by PEPS calculations [11] [compare Fig. 2 (b)].

Order parameter and spin stiffness
A fundamental indication for the validity of spin-wave theories is generally given by the order parameter M 0 (Fig. 5) and the spin stiffness (Fig. 6). The influence of quantum fluctuations is strong where they are small, and the primary assumption that the system can be described by a semi-classical spin-wave state begins to falter in such a case. Since the MSW formalism only takes quantum fluctuations partially into account, a small order parameter and/ or spin stiffness also suggests that the true quantum ground state could be completely disordered.
Interestingly in the regions of largest spatial isotropy of the interactions, i.e. around α = 1 (isotropic triangular lattice) and at large α (isotropic square lattice), the order parameter M 0 coincides with that of LSW theory. As we will see later, at the points α = 0, 1, ∞ the ordering vector found with the present MSW approach also exactly matches the classical one, due to symmetry.
In the square lattice limit α → ∞, the order parameter attains the value M 0 = 0.435 in the thermodynamic limit, which is very close to M = 0.437 as extrapolated from quantum Monte Carlo calculations [25]. For the spin stiffness Ref. [25] obtained ρ /α = 0.270; the MSW method returns the only slightly larger value ρ /α = 0.272. It appears that in this special case the main quantum corrections are correctly captured by MSW theory. The large values of the order parameter (around 87% of the theoretical maximum) and of the spin stiffness support the assumption that the classical picture remains essentially valid in the large-α limit.
The loss of LRO in the 1D-limit (α → 0) is reflected in the breakdown of the order parameter M 0 , which occurs at a finite value of the inter-chain coupling, α ≈ 0.18 (note that within LSW theory the order parameter vanishes only for α → 0). This coincides with the appearance of imaginary spin-wave energies as discussed in the previous section. At small but finite α the spin stiffness ρ yy essentially vanishes, which is characteristic of a 1D-like state that consists of effectively decoupled A single XY chain can be solved exactly by Bethe-Ansatz equations, and by use of twisted boundary conditions one can obtain the exact solution for the spin stiffness ρ xx = 1/π ≈ 0.318 [26]. Our MSW result of ρ xx ≈ 0.308 lies surprisingly close. For one-dimensional models it is known that a non-zero spin stiffness is accompanied by quasi long-range correlations with powerlaw decay. The critical nature of the state in the 1D-like phase reflects itself also in the fact that finite-size effects play an important role.
The MSW order parameter can be compared with results from exact diagonalization (ED) and PEPS calculations. In both cases the Fourier transform of the spin-spin correlations, the static structure factor The comparison of ED and MSW results shows that MSW is quantitatively reliable in the Néel phase for α 1.66. For smaller α values the comparison is more problematic: in particular, while ED and PEPS confirm the existence of an ordered spiral region for α around 1, the magnitude of the order parameter appears to be largely overestimated by MSW theory, which is not surprising considering the partial account of quantum fluctuations by the MSW approach. In particular, MSW theory produces the counterintuitive prediction that the frustrated spiral phase for 0.18 < α < 1.35 has an order parameter which can be larger (around the isotropic α = 1 point) than that of the unfrustrated case of the square lattice (recovered for α → ∞). We observe that MSW theory is only moderately improving upon linear spin wave theory around the α = 1 point for what concerns the quantum fluctuations of the order parameter -in particular, its prediction for M 0 essentially coincides with that of LSW for α = 1.
In summary, from the MSW order parameter M 0 we can derive a loss of LRO at α 0.18, and the spin stiffness suggests a strong weakening of inter-chain correlations already at α 0.35.
The spin-stiffness also decreases strongly upon approaching the parameter region 1.35 α 1.66.
Together with the appearance of imaginary spin-wave frequencies for α 0.18 and 1.35 α 1.66 this strongly indicates the appearance of disordered phases in these regions. Due to its semiclassical character the MSW Ansatz is not adapted to properly describe these phases, and we must resort to methods such as PEPS which take quantum fluctuations into account more completely. However, in the rest of parameter space magnetic LRO order seems to survive quantum fluctuations.
In the next section we investigate in detail the nature of the ordered phases.

B. Ordering vector, spin-spin and chiral correlations
In this section we introduce several observables which reveal the type of order appearing in the system, and which will be used in the discussions of sections III C and III D.
The ordering vector is a direct outcome of MSW theory, and can be extracted from the ED data by determining the position of the peak of the static structure factor, Eq. (23). Figure 7 displays the x-component of the ordering vector Q (with Q y = 0). Three limiting values are known. For α = 0 intra-chain antiferromagnetic (Néel) order is described by Q = πx. For α → ∞ square-lattice Néel order is described by Q = 2πx. In the isotropic lattice (α = 1) the threefold symmetry forces the ordering vector to Q = 4π 3x . (The ED and PEPS results deviate at α = 1 because the required threefold symmetry is broken by the shape of the simulation clusters, Fig. 3.) The importance of optimizing the ordering vector is apparent in Fig. 7 when comparing the MSW results to the classical (and LSW) curve.
The spin-spin correlations of nearest neighbors shed further light on the order properties. We analyze them through the two-site total spin, In Fig. 8 we plot it for nearest neighbors. This quantity vanishes if the spins are in a singlet, which is equivalent to perfect anticorrelation, takes the value 3 4 if they are uncorrelated, and the value 1 if the spins form a triplet, which means perfect correlation. For PEPS and ED we report the values of K ij averaged over the central spins, where boundary effects are minimal. Spiral phases carry not only a magnetic order parameter, but also a chiral order parameter. In particular, a vector chirality [27] can be defined on an upwards pointing triangle with counter-clockwise labeled corners (i, j, k) as κ ∆ = 2 and on a downwards pointing triangle with counter-clockwise labeled corners (i, l, j) as κ ∇ = Long-range chirality correlations are defined as [28] where the triangle pairs (∆, ∇) and (∆ , ∇ ) share a τ 1 edge. In Fig. 9 we   between Néel and spiral correlations at the slightly larger value of α ≈ 1.44. In the case of MSW theory, the jump in the ordering vector is realized when going across the breakdown region, namely when passing from α ≈ 1.66 (which is the lower bound to the Néel phase within MSW theory) to α ≈ 1.35 (which represents the upper bound of the spiral phase). In particular, all three different approaches point to the fact that Néel order persists to much lower α than the classical value α = 2.
The persistence of Néel order over a larger parameter region than in the classical case is reminiscent of what is observed in other models. Indeed, quantum fluctuations generally stabilize states where spins are ordered collinearly (see e.g. Refs. [29,30]), a property that is reproduced by the MSW Ansatz with ordering vector optimization. The mechanism behind it is strongly connected to order-by-disorder phenomena [30].
The abrupt transition from a phase with dominant Néel correlations to a phase with dominant spiral correlations is confirmed by the spin-spin correlations as displayed in Fig. 8. Anti-correlation along the strong τ 2 -bonds and correlation along the weak τ 1 -bonds are characteristic of a 2D-Néel ordered state; these correlations decrease rapidly for α < α crit outside of the Néel-ordered phase.
The change of the type of order is further supported by the overlap | ψ α |ψ ∞ | of the new ground state with the 2D-Néel ordered state of α = ∞, which we plot in Fig. 10 for the ED: above the phase transition it still attains a finite and quite large value, while it vanishes identically below the phase transition. Finally, the onset of strong chiral correlations shows that the new phase is indeed a spirally correlated one (Fig. 9).
The breakdown region of MSW theory, 1.35 α 1.66, is strongly suggestive of the loss of magnetic LRO, corresponding to a spin-liquid state. This region of disordered behavior is only roughly consistent with that indicated by PEPS calculations [11] for the appearance of a shortranged spin-liquid phase, namely 1.2 α 1.4. Nonetheless it is tempting to associate the breakdown of MSW theory to this quantum-disordered phase.

D. Persistence of 1D quasi-LRO up to finite inter-chain couplings
In MSW theory, at α ≈ 0.18 the order parameter M 0 breaks down. This is an indication of the transition to a phase without magnetic LRO such as the one reproduced in the one-dimensional limit which suggests that, if we took quantum fluctuations into account in a more accurate way than in MSW theory, spiral (2D) order would be probably lost below α 0.35.
In the limit α → 0 we find that the intra-chain spin-spin correlation K i,i+τ 1 ≡ K τ 1 0.355 lies close to the exact result that can be obtained by use of the Jordan-Wigner transform, K τ 1 = − 1 π − 1 π 2 + 3 4 0.330. The inter-chain spin-spin correlations on the other hand vanish, which is equivalent to K i,i+τ 2 ≡ K τ 2 = 3 4 . We find that the ordering vector of the MSW theory compares very well to the one computed by ED in the entire range of α (Fig. 7). Especially the very weak dependence of Q on α near the 1D-limit is found in both approaches, contrary to classical and LSW theories which exhibit linear dependences on α. This means again that quantum fluctuations stabilize collinear order within the chains. This is confirmed by ED: the overlap 6 i=1 ψ α |ψ i 0 2 of the ground state with the subspace spanned by the six-fold degenerate [44] ground-states of α = 0 remains very large (almost 80 percent) up to α ≈ 0.5 (Fig. 10). Moreover, the small chiral correlations in both PEPS and ED show that for small α there is no spiral long-range order.

E. Momentum distribution of the hardcore bosons
After characterizing the zero-temperature phase diagram of the spin model in the previous sections, we now wish to make contact with the cold-atom implementation of such a model via hardcore bosons. The most common observable in cold-atom experiments is the momentum distribution [1], which exactly corresponds to the static structure factor for S = 1/2 spins (23) via the spin-to-hardcore-boson mapping described in Sec. II E. Figure 11 shows  The possibility of observing the BKT transition is a particular advantage of MSW theory.
The BKT phase with algebraic order is generally predicted by linear spin-wave (LSW) theory to remain stable at arbitrary temperatures. The non-linearities contained in MSW theory allow the disruption of quasi-LRO and the transition to the short-range-ordered (SRO) phase. However, vortex-antivortex excitations are not explicitly present in the theory, and therefore in principle T BKT cannot be accurately estimated.

A. Spin-spin correlations
An important observable for the analysis of a temperature-dependent phase diagram is the two-point correlation function  In reality, the thermal onset of a gap we observed within MSW is typically very gradual, and a clear identification of the transition point via the gap is generally problematic. This observation can be understood on the basis of a well-known fact: the chemical potential of the half-filled DM boson gas, which determines the existence of a gap, cannot vanish at finite temperature because of the absence of Bose-Einstein condensation in two dimensions. As a consequence, we find a finite gap at any finite temperature, which means that the correlations decay exponentially at long distances. This suggests that, strictly speaking, MSW theory is not able to describe the BKT transition. However, for temperatures much lower than the BKT transition (estimated as explained above) the gap is very small, being below our numerical precision. For all practical purposes such a small gap entails a decay of correlations which is not distinguishible from an algebraical decay.
Moreover in a selected region of the phase diagram (corresponding to the spiral phase) the gap is seen to increase drastically around the estimated BKT transition temperature, and correspondingly the correlation function is seen to decay much more rapidly above that temperature. Hence we conclude that MSW theory still accounts for one of the most salient features of the BKT transition, namely a discontinuous behavior of correlations as the temperature is increased.
Finally we can extract from the correlations the temperature at which the MSW formalism breaks down. It is characterized by the complete loss of all correlations, even to the nearest neighbor. This behavior, occurring at temperatures of the order of the coupling strength, is clearly an artifact of the method, since in real systems the complete loss of correlations occurs only at extremely large temperatures where spin-spin interactions become negligible.

B. The phase diagram
In this section we present the finite-temperature phase diagram of the anisotropic triangular lattice model obtained via the MSW method with ordering vector optimization [45]. We derive it from the observables introduced in the previous section. For reference we first present a summarizing sketch of the phase diagram in Fig. 13, which introduces the phase labels we will refer to in the following discussion. Table I lists the main properties of these phases.
First, at small α, there is a phase with properties similar to the algebraic 1D-Néel-like state found at T = 0 but with exponential decay of intra-chain correlations (phase A). Further, we generally find that the phase diagram contains two quasi-LRO regions: a region at intermediate α corresponding to spiral quasi-LRO (phase B), and another region at large α which is characterized by Néel quasi-LRO (phase D). These phases undergo BKT transitions to similar phases with short-range order (SRO), phases C and E, respectively. Moreover, between them lies a region where imaginary frequencies occur in the spin-wave dispersion relation, which can be interpreted as an indication for an extremely short-range-ordered phase (phase F). This general structure of the phase diagram is supported by all the observables we investigate. At large T the MSW method breaks down (see sec. IV A) and therefore does not allow for any interpretation in that domain (region G).
It is also important to note that our calculations cease to converge properly for too low temperatures, when the chemical potential becomes smaller than the accuracy of our numerical integrations.
Depending on the region of the phase diagram the lowest temperatures for which appropriate results could be derived vary from less than one-tenth of a percent to several percent of the coupling strengths. This pathology is not observed at T = 0 (as calculated in section III) because an exact vanishing of the chemical potential allows the special treatment of the zero-mode as captured in Eqs. (16) and (17). Except for some points, we typically calculated down to T /(t 1 + 2t 2 ) = 0.025.
Since the bond strengths are the only energy scales in the problem, it seems a reasonable assumption that our finite temperature results can be analytically continued down to T = 0 + without encountering discontinuities (except possibly at exactly T = 0, where -in contrast to any finite temperature -Bose-Einstein condensation of the DM bosons becomes possible). Nevertheless, this issue should be kept in mind in the following analysis.
The breakdown of the calculations for too low temperature can be clearly seen in Fig. 14, which displays the phase diagrams obtained from several observables.
A natural starting point for a thorough analysis of the temperature dependent phase diagram is given by the respective ground state phases. Proceeding from small to large α, we divide the analysis in sections corresponding to different ground state behavior.

1D-like phase (phase A)
The one-dimensional quasi-ordered ground state phase for which we found strong indications below α ≈ 0.18 becomes a short-range-ordered phase at finite temperature (phase A). It is characterized by vanishing correlations between neighboring chains already at zero temperature. The finite gap leads to exponentially decaying intra-chain correlations for all T . This is consistent with the expected finite-temperature behavior above a ground state with quasi-LRO.
The assumption that this low-α phase really describes decoupled chains is reinforced by the component ρ partial yy of the spin stiffness, which vanishes in this region, and by the ordering vector that takes on the 1D value (π, 0), similar to the equivalent of the ground state phase diagram.
Moreover, neighboring spins on different chains are uncorrelated whereas nearest-neighbors on the same chain are anti-correlated. It is remarkable that this phase is preferred over the quasi-ordered spiral phase with rising temperatures. In section III D we have seen that quantum fluctuations stabilize 1D-Néel quasi-order. The same mechanism is at work here: collinear spin correlations are stabilized by fluctuations, in this case thermal ones.
Note that in the 1D-like phase A the inter-and intra-chain correlations behave completely differently. In the rest of the phase diagram they follow one and the same pattern, since in a truly two-dimensional structure the correlations in one direction typically cannot disappear without affecting the correlations in the other one.

Spiral phases (phases B and C)
At intermediate inter-chain couplings 0.18 α 1.35 and low temperature we find a spiral phase with magnetic quasi-LRO (phase B). It can be seen as the finite-temperature continuation of the spirally ordered ground state phase. At larger temperatures a BKT transition to a phase with a spiral ordering vector but with an exponential decay of correlations occurs (phase C). Within our resolution of the phase diagram it seems that phase C disappears on the large-α side of the B-phase dome, and that phase B is delimited at large α by the F region where imaginary spin-wave frequencies appear. Furthermore, on the low-α side phase C becomes extremely narrow and is almost immediately followed by a transition to phase A described in the previous section. The broadest extent in temperature of C is around α = 1.
At the isotropic point α = 1 the BKT transition from B to C is approximately located at T BKT /(t 1 + 2t 2 ) = 0.0836. Quantum effects lower the transition temperature considerably from the classical value T cl BKT /(t 1 + 2t 2 ) = 0.165 found by classical Monte Carlo simulations [36]. A purequantum self-consistent harmonic approximation, developed in Ref. [37], gives T BKT /(t 1 + 2t 2 ) = 0.0625. The fact that MSW theory gives a significantly higher estimate is not surprising given that Ref. [37] takes vortex-antivortex excitations explicitly into account while MSW theory does not.
We also find that in the same domain of large frustration (α close to 1) where spiral quasi-order is most stable, the breakdown of the theory occurs at a lower critical temperature.
We note also a strong drop of the transition temperatures around α ≈ 0.4. In Fig. 13 this is marked by a dashed line which separates phase B from a phase B' with similar properties. We believe that this behavior is a numerical artifact and that in fact B and B' are one and the same phase.

Spin-liquid candidate region (phase F)
At the high-α side of the spiral phases we find an extended region where the spin-wave dispersion acquires imaginary modes. This means that MSW theory predicts an instability here. The width of this region in α stays approximately constant but it moves to smaller α with increasing temperature, leaving space to the collinear short-range-ordered phase (E). The region F extrapolates well down to the suspected spin-liquid phase between 1.35 α 1.66 at T = 0. Given that MSW is seen to break down at a putative spin-liquid phase at T = 0 due to its lack of order, a fortiori one can expect MSW to break down in the same parameter range at finite temperatures, because at finite T the theory would be required to describe not only the ground state but also the excitations on top of it.
Note also that the spin-stiffness decreases upon approaching this region, which could be interpreted as a precursor of a short-range-ordered phase.

2D-Néel states (phases D and E)
As expected from BKT theory, when going to finite temperatures the 2D-Néel ground state first changes into a low-T quasi-long-range ordered phase (phase D), which at a temperature T BKT undergoes a transition into a high-T short-range-ordered phase (phase E). Both are characterized by an ordering vector at the 2D-Néel value Q = (2π, 0). Furthermore neighboring spins which share a diagonal bond are strongly anticorrelated whereas neighboring spins which lie on the same chain are positively correlated.
The square XY lattice, which is reached as α ≡ t 2 /t 1 → ∞, has been extensively studied in the past. The classical BKT-temperature T cl BKT /t 2 = 0.695, which has been calculated by use of classical Monte Carlo simulations [38], is significantly lowered in the quantum limit to around T BKT /t 2 ≈ 0.35 (Quantum Monte Carlo calculations [39,40]). Our MSW results yield a BKT temperature of T BKT /t 2 ≈ 0.27 at α = 100, where the system has practically reached the square lattice limit [46]. Once again the disagreement with the MSW result is not surprising, given that this theory does not account properly for vortex-antivortex excitations. In particular the BKT line for α 1.6 is not very distinct, due to the location problems mentioned in Sec. IV A. Therefore its quantitative value should be interpreted with caution, especially in this parameter range. However, the qualitative behavior of the phase diagram seems to be described correctly.
In the following section we turn to observables which show more clearly how order in the different phases persists at finite temperatures.

C. Observables distinguishing between LRO and SRO
Here we focus on observables which help to distinguish between different types of regimes (i.e. quasi-long-range order or short-range order), namely the partial Gaussian spin stiffness Υ partial , the gap ∆, the entropy, and the occupation of the zero mode n k=0 .

Entropy, spin stiffness, and gap
The entropy [Eq. (11)] shows behavior consistent with the quasi-ordered character of the lowtemperature 2D-Néel phase D and the spiral phase B in the sense that it is smaller in phases with stronger order as can be seen in Fig. 14 (g). Correspondingly, in phases B and D Υ partial is large At the BKT transition no sharp change occurs in these observables as may be expected. However, the contour lines of the spin stiffness, the entropy, and the gap, all seem to be consistent with the shape of the T BKT curve.
We report the gap ∆ for two representative values of α in Figs. 15 (a) and (b). It evolves smoothly through the BKT transition in the Neel phase, whereas in the spiral phase it is seen to display a sharp increase right above the BKT transition up to the transition to phase A. On the contrary, for all phases a sharp increase of ∆ (accompanied by a sharp drop of Υ partial ) can be discerned at the breakdown temperature where correlations are completely lost.
In phase A (better visible for smaller values α which are not shown) we find that the gap is almost a linear function of temperature up to very close to the breakdown temperature. This is typical of critical systems where the temperature introduces the only energy scale, which is reflected in the gap. Similar to what was seen for the gap in the previous section, at the transition between phases C and A n k=0 drops from very large values to something of the order of 1. Afterwards it changes only slightly with T . It is useful to remember that the average mode occupation, k n k /N , equals S by virtue of constraint (10). This means that in the 1D-like short-range-ordered phase A n k=0 is relatively small but still larger than the occupation of the other modes.
The behavior of n k=0 is different for the 2D-Néel phase. At the BKT line that we extracted from the analysis of the two-point correlation functions C m τ 1 and C m τ 2 , n k=0 decreases strongly but smoothly. Up to the breakdown point its values are still several times larger than in the 1Dlike phase, however. This supports our identification of the BKT transition; but the smoothness of n k=0 also shows the reason why the observables of section IV C could not point out a sharp transition.

E. Discussion
The finite-temperature phase diagram is observed to be a natural extension of the groundstate phase diagram. We find that zero-temperature long-range ordered phases are reflected in finite-temperature phases with quasi-LRO, while phases with quasi-LRO at T = 0 turn into shortrange-ordered phases at any finite temperature. When temperatures are at or below a few percent of the coupling strengths, the main characteristics of the ground state phase diagram are retrieved, with a short-range 1D-like phase (A), and two quasi-ordered phases [one with spiral properties near the isotropic triangular limit (B) and one with 2D-Néel-like characteristics at large values of α (D)] which are separated by a potential spin liquid (F). This last phase was identified by (i) the breakdown of MSW theory, which indicates that the assumption of an underlying ordered state is invalid, and (ii) the lowering of the spin stiffness as this phase is approached. In A we show further evidence for spin-liquid behavior in this phase at T=0 gleaned from very limited exact-diagonalization results.
We have given a rough estimate for T BKT over the entire range of anisotropies. We find agreement in the rough magnitude of T BKT at points where estimations of the BKT temperatures computed by other methods exist. In our results the BKT transition is more clearly visible in the correlations for the spiral phase than for the 2D-Néel phase.

V. CONCLUSIONS
We have extended Takahashi's modified spin-wave theory by an optimization of the ordering vector, which allows to account for order that deviates from the classical one.
We have used this method to calculate the ground state phase diagram of the spatially anisotropic triangular lattice with S = 1/2 spins and XY interactions. We found the expansion of a quasi-ordered 1D-like phase to finite inter-chain couplings, a spiral phase, and a 2D-Néel phase. At the transition between the latter two the breakdown of MSW theory indicates the loss of LRO.
We have extended this phase diagram to finite temperatures and computed Berezinskii-Kosterlitz-Thouless transitions, although the results are to be interpreted only semi-quantitatively because MSW theory does not explicitly account for vortex-antivortex excitations. We find that the ground state phases clearly imprint their properties on the finite temperature phase diagram at low temperatures, with long-range ordered phases being replaced by quasi-ordered ones and quasi-ordered phases by short-range ordered ones.
Qualitative and even quantitative agreement with PEPS and ED calculations was found in the regions where magnetic LRO is to be expected. In particular it has been shown that MSW theory with ordering vector optimization is able to account satisfactorily for the main quantum corrections to the ordering vector. Furthermore, our calculations show that the spin stiffness is a useful observable for finding candidate regions for spin-liquid behavior in the ground state. Indeed the breakdown of MSW theory, or the very weak stiffness of the magnetic order it predicts, can be used as indications of the absence of long-range order in the exact ground state of the model.
We find two main recurrent features for strongly frustrated quantum-magnets in two dimensions: collinear order is considerably stabilized by quantum and/ or thermal fluctuations against spiral order, and ordered or quasi-ordered phases characterized by different forms of order (collinear vs. spiral) do not continuously connect to each other, but they rather seem to be separated by quantum with the geometry depicted in Fig. 3. The system Hamiltonian, Eq. (2), commutes with the total magnetization along the z axis, S tot z , so that excited states can be classified on the basis of this quantum number [47]. Figure 16 shows the excitation energies of the first excited states in each S tot z sector (up to S tot z = 11) with respect to the minimum energy in each sector, E 0 (S tot z ) [E 0 (S tot z = 0) corresponds to the ground state energy]. Upon varying α we observe a significant evolution of the low-energy spectrum of the system, which points at the widely different regimes explored by the system.
In particular, in the spirally and Néel ordered phases -exemplified in Fig. 16 by α = 1 and α = 2, respectively -we observe that in each S tot z sector there are a few states lying close to the minimum energy one, and separated from the other excited states by a large gap. According to a standard 'tower-of-states' argument [41], these low-lying states are expected to collapse to the ground state in the thermodynamic limit, giving rise to degenerate superpositions of all S tot z sectors, each breaking the U(1) rotational symmetry of the Hamiltonian and displaying spiral or Néel order. The higher-energy states will instead reproduce the true excitation spectrum in the thermodynamic limit.
This tower-of-states feature is on the contrary absent in other regions of the phase diagram, in which the energy levels in each S tot z sector are more homogeneously spaced. The absence of a low-lying multiplet of states separated from the higher energy states by a large gap is observed in models whose ground state is generally considered to be a spin liquid [42]. We therefore introduce an observable aimed at quantifying the extent to which the spectrum exhibits the expected features in presence of spontaneous symmetry breaking in the thermodynamic limit. We consider the average maximal level spacing∆ max , defined as namely∆ max is the maximal level spacing in each S tot z sector, averaged over the N S + 1 = 12 sectors considered. The maximal level spacing is extracted by considering the lowest 10 levels E i (S tot z ), which captures the behavior of the low-energy part of the spectrum. The above quantity is chosen so that it will be maximal in presence of a large separation between the low-lying tower of states and the higher-energy spectrum, while it will be minimal in presence of homogeneously spaced levels in each sector.
When plotting∆ max as a function of α, as shown in Fig. 16, we observe two very pronounced relative minima, at α ≈ 0.6 and α ≈ 1.4. Remarkably these two minima correspond to the two regions in parameter space where PEPS calculations predict the occurrence of a spin-liquid phase [11] (compare Fig. 2). Hence the lack of the tower-of-states feature in the spectra of a small cluster is consistent with the PEPS prediction.