Modified spin-wave theory with ordering vector optimization II: Spatially anisotropic triangular lattice and $J_1J_2J_3$ model with Heisenberg interactions

We study the ground state phases of the $S=1/2$ Heisenberg quantum antiferromagnet on the spatially anisotropic triangular lattice and on the square lattice with up to next-next-nearest neighbor coupling (the $J_1J_2J_3$ model), making use of Takahashi's modified spin-wave (MSW) theory supplemented by ordering vector optimization. We compare the MSW results with exact diagonalization and projected-entangled-pair-states calculations, demonstrating their qualitative and quantitative reliability. We find that MSW theory correctly accounts for strong quantum effects on the ordering vector of the magnetic phases of the models under investigation: in particular collinear magnetic order is promoted at the expenses of non-collinear (spiral) order, and several spiral states which are stable at the classical level, disappear from the quantum phase diagram. Moreover, collinear states and non-collinear ones are never connected continuously, but they are separated by parameter regions in which MSW breaks down, signaling the possible appearance of a non-magnetic ground state. In the case of the spatially anisotropic triangular lattice, a large breakdown region appears also for weak couplings between the chains composing the lattice, suggesting the possible occurrence of a large non-magnetic region continuously connected with the spin-liquid state of the uncoupled chains.


I. INTRODUCTION
Low-dimensional frustrated quantum spin systems can display an intriguing interplay between order and disorder: classical order has been shown to be quite resilient in two or three dimensions [1][2][3][4]; frustration, however, can lead to the melting of magnetic long-range order (LRO) and the emergence of quantum disordered states like valence-bond solids or resonating valence bond states [5,6]. Understanding such magnetically disordered quantum phases is important for the search for fractionalized excitations in two dimensions [5], as well as for the understanding of the behavior of layered magnetic insulators/metals in which magnetism is disrupted by charge doping, leading to dramatic phenomena such as superconductivity at high critical temperature [7][8][9].
In this work, we investigate the S = 1/2 Heisenberg antiferromagnetic Hamiltonian on two-dimensional frustrated lattices making use of Takahashi's modified spin-wave (MSW) theory [29], supplemented with the optimization of the ordering vector [30]. In a previous paper [31], we have shown that (for the SATL with XY interactions) this approach provides a significant improvement over conventional spin-wave theory (as well as over conventional MSW theory), as it allows to correctly account for the dramatic quantum effects occurring to the form of order which appears in frustrated quantum antiferromagnets, and for the quantum corrections to the stiffness of the ordered phase. In particular, a very low stiffness, or the complete breakdown of the theory, provide strong signals that the true ground state might be quantum disordered; hence, this method serves as a viable approach to finding candidate models potentially displaying spin-liquid behavior. For a more detailed description of the formalism we refer the reader to Ref. [31].
Here, we apply this MSW theory with ordering vector optimization to the Heisenberg SATL, as well as to the square lattice with nearest, next-to-nearest and next-to-next-tonearest neighbor couplings (the J 1 J 2 J 3 model [27,[32][33][34]). Both models feature a very complex T = 0 phase diagram, with spirally and collinearly ordered regions, whose ordering vector is subject to strong quantum corrections with respect to the classical (S → ∞) limit.
They also feature extended breakdown regions for MSW theory, pointing at the possible spinliquid nature of the true ground state of the system. Comparison with numerical results coming from exact diagonalization and projected-entangled-pair-state (PEPS) calculations show that MSW theory correctly accounts for some of the most salient features of the quantum phase diagram of these systems, and that it hence represents a very versatile tool to probe the robustness (or the breakdown) of a semi-classical description of the ground state of frustrated quantum magnets.
The remainder of this paper is organized as follows: Section II presents the ground state phase diagram of the SATL with nearest-neighbor Heisenberg interactions; in Section III, we calculate the ground state phase diagram of the J 1 J 2 J 3 model; finally, in Section IV we present our conclusions. The technical aspects of MSW theory applied to Heisenberg antiferromagnets are presented in the Appendix.

II. MSW THEORY ON THE SPATIALLY ANISOTROPIC TRIANGULAR LATTICE WITH NEAREST-NEIGHBOR HEISENBERG-BONDS
The triangular lattice with Heisenberg interactions has been considered as one of the first candidate systems for quantum-disordered behavior in the ground state [5]. Recently, the phase diagram of the spatially anisotropic triangular lattice (SATL) up to values of α ≡ t 2 /t 1 = 1 has been studied by Yunoki and Sorella using variational quantum Monte Carlo methods [16]. They find that the gapless spin-liquid phase of the isolated chains (t 2 = 0) persists also at finite coupling up to a critical value α ≈ 0.65, followed by a gapped spin liquid; for α ≈ 0.8 the gap closes and the system undergoes an ordering transition to spiral order, continuously connected with the 3-sublattice order of the isotropic Heisenberg antiferromagnet (α = 1). This scenario is still controversial, however: studies based on low-energy effective field theory for the description of the coupled chains in the case α < 1 indicate that the system might still exhibit long-range antiferromagnetic order even for very weak coupling among the chains. This form of order results from high-order perturbation theory in the inter-chain coupling, and it is necessarily very weak, given that numerical methods cannot detect it. Its observation is clearly beyond the capabilities of our MSW approach. Coming from the large-α limit, series expansions by Weihong et al. indicate that 2D-Néel order -appearing on the square lattice defined by the dominant t 2 -couplings -persists down to α 1.43, followed by a phase without magnetic order in the interval 1.1 α 1.43 [15]. Below this region the authors find incommensurate spiral order connecting continuously to the isotropic point α = 1. In Ref. [35], qualitative similar results have been obtained using the Schwinger-boson approach. The resulting phase diagram differs strongly from the classical one, which is characterized by spiral order for 0 < α < 2, and by Néel order for α ≥ 2. The classical phase diagram is contrasted with the quantum mechanical one (composed from Refs. [15] and [16]) in Fig. 1  Ref. [19] showed that, for a finite inter-chain coupling, spinons tunnel between chains in bound pairs with S = 1 (so-called triplons), so that the fractionalization in two dimensions is strictly speaking not present. Ref. [19] argues that the spinons in Cs 2 CuCl 4 are descendants of the excitations of the individual 1D chains and not characteristic of any exotic 2D state.
This further reinforces the idea of a quasi one-dimensional behavior up to relatively high inter-chain interactions mentioned in the previous paragraph.
A. MSW predictions for the ground-state phase diagram In this section, we discuss the ground-state phase diagram resulting from the predictions of MSW theory for the S = 1/2 SATL with nearest-neighbor (NN) Heisenberg interactions.
In order to assess the validity of MSW results, we compare them with exact diagonalizations (ED). Using the Lanczos method, we compute the ground state of small clusters of 14, 24, and 30 spins. The considered geometry for the 30-spin system can be found in Fig. 2. Convergence in the self-consistent equations of MSW theory with ordering vector optimization, Eqs. (A3-A7, A9), cannot be achieved in selected regions of the ground state phase diagram, namely for α 0.65 and for 1.14 α 1.3. (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 spin-wave frequencies, Eq. (A5), 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-liquids predicted from Refs. [15,16] [compare Fig. 1 (b)]. Both for α < 1 and α > 1, we find that the breakdown region of MSW appears to be fully contained within the region of SL behavior (either gapped or gapless) estimated in Refs. [15,16]. Hence MSW theory is seen to possibly underestimate the width of the quantum-disordered regions in the phase diagram, which is to be expected due to the partial account of quantum fluctuations given by MSW theory.   [16], also plotted in Fig. 3. For comparison, we also show the curve that they obtain with a projected-BCS (p-BCS) wave-function. In the isotropic triangular lattice, the MSW energy compares also favorably to the data from the Green's function Monte Carlo method with stochastic reconfiguration (GFMCSR) from Ref. [40], but both energy and order parameter (see section II A 3) lie closest to the variational quantum Monte Carlo calculation from Weber et al. [38], who used a mixture of a BCS wave-function and a wave function with spiral order as their starting point (BCS+spiral).
At α = 0 the MSW value E 0 (α = 0) = −0.4647 is relatively close to the exact result  given in brackets [16,38]. FN is short for lattice fixed node and FNE for the improved FN effective Hamiltonian method [16]. Furthermore included are the exact result of the Heisenberg chain in the thermodynamic limit, the ED results for the 30-spin cluster, and the ED results from finite size extrapolations of calculations on clusters of up to 40 spins [37]. Also given are the estimates of LSW theory from Ref. [16] and the Green's function Monte Carlo method with stochastic reconfiguration (GFMCSR) [40]. The last column gives the staggered magnetization or, in the case of MSW theory, the population of the zero mode M 0 .
of the one-dimensional case, −(ln 2 − 1/4) = −0.44315. However, it is located below the exact value. This apparent puzzle is easily resolved by noticing that the MSW method is not variational due to the incomplete inclusion of the kinematic constraint (see Appendix).
We also notice that the ground state energies derived from ED of the 30-site system lie very close to the values from the other methods except in the 1D phase. This could be attributed to the small system size: if the interpretation is correct that for small α the Heisenberg SATL is in a 1D-like phase with algebraic correlations, it is natural that finite size effects play a very important role in the critical 1D phase. This would explain the strong deviation of the ED energy in that parameter region.
On the square lattice (α → ∞), Takahashi showed already twenty years ago the extremely good performance of MSW theory [29]: its ground state energy per spin is −0.6699, which We further display the value obtained in the isotropic limit by Ref. [38] by use of a VMC method with a mixture of a BCS and a spiral ordered wave-function (BCS+spiral), and the exact result of the 1D limit. The numbers in the labels of the curves are the respective system sizes considered.

Order parameter and spin stiffness from MSW theory.
Our next step is to determine the regions where the presence of a finite order parameter M 0 and spin stiffness Υ reveal magnetic long-range order (LRO). Even when M 0 and Υ are finite, a caveat is still in order: a finite order parameter with a very small stiffness might suggest that taking quantum fluctuations more completely into account than in MSW theory could lead to a completely disordered state.
The order parameter M 0 , drawn in Fig. 4, shows that magnetic LRO is present in the intervals 0.65 < α < 1.14 and α > 1. 3. This is to be contrasted with linear SW (LSW) The numbers in the labels of the curves are the respective system sizes considered in the calculations.
theory, which predicts the breakdown of magnetic order only for α 0.3 [42]. However, in the isotropic case, α = 1, MSW theory predicts a stronger order parameter than what is predicted by LSW, as well as by most of the other numerical estimates, which are presented in Table I. In the square lattice limit, α → ∞, on the other hand, both MSW and LSW theory attain a staggered magnetization of 0.303, which compares favorably with the most recent estimates M 0 = 0.311, based upon diagonalizations of small clusters of up to 40 spins [37].
The MSW order parameter drops drastically upon approaching the region 1.14 α 1.3 and when reaching the region α 0.65 from above, the regions where the self-consistent description breaks down, further corroborating the assumption that in these regions magnetic LRO disappears in the true quantum ground state. This assumption is strongly reinforced by considering the Gaussian spin stiffness (Fig. 5): It vanishes at α = 0.65 and it drops significantly when approaching α = 1.14 from below.

Spin and chirality correlations from MSW theory
Now we describe the ordered phases found by the MSW Ansatz for the Heisenberg SATL in more detail. To this end, we analyse the following quantities 1. The ordering vector Q (Fig. 6). Three limiting values for the ordering vector 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 ; 2. The spin-spin correlations (Fig. 7). We analyze the spin-spin correlations of nearest neighbors through the two-site total spin, 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; 3. The mean chiral correlations (Fig. 8). Spiral phases carry not only a magnetic order parameter, but also a chiral order parameter. In particular, a vector chirality can be defined on an upwards pointing triangle with counter-clockwise labeled corners (i, j, k) as [45] and on a downwards pointing triangle with counter-clockwise labeled corners (i, l, j) as Chirality correlations are defined as [46] where the triangle pairs (∆, ∇) and (∆ , ∇ ) share a τ 1 ≡ (1, 0) edge. In Fig. 8  The numbers in the labels give the system sizes.
A comparison of these quantities shows a spiral phase at around 0.65 α 1.14 and a 2D-Néel ordered phase for α 1.3. Moreover, when approaching α ≈ 0.65 from above, the ordering vector, the spin-spin correlations and the ground state energy approach their respective 1D values. This is an indication that below α ≈ 0.65 the true ground state of the system may enter a 1D-like spin-liquid phase. Nonetheless, the vanishing of the spin stiffness for α → 0.65 + is not consistent with the onset of a gapless 1D spin-liquid phase, for which the spin stiffness should remain finite. Hence, the MSW results rather suggest that the phase appearing below α = 0.65 is a gapped spin liquid, and that the gapless 1D spin-liquid phase, connected continuously with the limit α = 0, is only attained for even smaller α. This seems consistent with the prediction of Ref. [16] that a gapped spin-liquid phase separates the spirally ordered phase from the 1D-like gapless disordered one.

Order parameter and correlations in comparison with exact diagonalization.
In the case of ED, the static structure factor  Finally, we notice that at α = 1 the ED results deviate from the isotropic value Q x = 120 • because the required threefold symmetry is broken by the shape of the simulation cluster, Fig. 2.
The nearest-neighbor spin-spin correlations T ij , Eq. (2), [71] are in qualitative agreement with the MSW results as well (Fig. 7). In particular, they show 1D-like behavior at small α, a spiral phase in an intermediate parameter range around the isotropic limit α = 1, and a 2D-Néel structure at large α. In summary, despite the significant deviations in the magnitude of the order parameter, both ED and MSW theory give a coherent picture, both qualitatively and quantitatively, of the evolution of the nature of spin-spin correlations upon increasing the α parameter, going from quasi-1D to spiral to Néel.

B. Discussion
Despite its limitations, the MSW approach with ordering vector optimization reproduces faithfully the main characteristics of the phase diagram as sketched in Fig. 1 (b), and thus remarkably improves on the results that were previously obtained for this model with The classical phase diagram of the J 1 J 2 J 3 model [34,[47][48][49] is sketched in Fig. 9 (b).
One identifies: IV) A second spiral phase, this time with ordering vector Q = (q, q); q → π/2 for J 3 → ∞, attaining the limit of two decoupled and Néel-ordered J 3 −sublattices.
This phase diagram is believed to change considerably in the quantum limit [27,[32][33][34]: In phase II quantum fluctuations select the columnar ordered states with Q = (π, 0) or Q = (0, π) from all the possible classical states. Furthermore, the Néel phase I increases in size considerably and Néel order persists up to the vicinity of the line (J 2 + J 3 ) /J 1 = 1/2.
In the vicinity of this line, the classical order is believed to be destabilized and to be replaced by a non-magnetic state. The controversy about the exact nature of the ground state in this highly frustrated region, however, is still not settled. In particular, it has been suggested that it could have the nature of a columnar valence bond crystal [50] with both translational and rotational broken symmetries, of a plaquette state with no broken rotational symmetry [27], or of a spin liquid with all symmetries restored [51][52][53][54][55].
In the following, we investigate the quantum model using the modified spin-wave (MSW) formalism, and compare it to recent results from projected entangled-pairs states (PEPS) calculations. The MSW lattice size is again N = 32 × 32. In most of parameter space, a lattice of N = 32 × 32 spins is essentially already converged to the infinite lattice, except close to a quantum critical point.
In Ref. [56], some of us reported numerical calculations of the J 1 J 2 J 3 model based on the projected entangled-pair state (PEPS) variational Ansatz for varying lattice sizes. In the following, we will focus on the extrapolations to the thermodynamic limit, except if stated otherwise.
We first discuss in more detail the special cases of the J 1 J 2 model (i.e., J 3 = 0) and the J 1 J 3 model (i.e., J 2 = 0). Both models have been studied before within the MSW formalism [30,[57][58][59][60][61]. On the one hand, we confirm existing results on the J 1 J 2 case, for which the optimization of the ordering wave-vector returns only two possible values (corresponding to Néel order [Q = (π, π)] or columnar order [Q = (π, 0) or Q = (0, π)]), and we give further insight into the spin stiffness and the dimer-dimer correlation functions. On the other hand, we analyze the J 1 J 3 model with optimization of the ordering wavevector, which proves crucial to correctly capture the quantum effects on the classical spiraling phases appearing in this case [30]. Finally, we give an overview of the entire quantum ground state phase diagram of the J 1 J 2 J 3 model.  magnetization that where obtained in Ref. [37] from diagonalization of small clusters. In agreement with other methods, e.g., exact diagonalization (ED) [37,62,63] or Schwinger bosons [64], MSW theory finds Néel order with Q = (π, π) at small J 2 /J 1 and columnar order with Q = (π, 0) or Q = (0, π) at large J 2 /J 1 (see Fig. 10). As it is well known from previous studies, there is a region between 0.56 In agreement with the MSW prediction, Q PEPS is located at the Néel value (π, π) up to J 2 /J 1 = 0.6, while above this it lies at the value of columnar order (π, 0).
We find a remarkable correspondence of the ground state energy per spin between the MSW prediction and ED results extrapolated to the infinite lattice from Ref. [37] [ Fig. 11 (a)]. Moreover, the noticeable kink associated with the Néel-to-columnar transi- tion of MSW theory at J 2 /J 1 = 0.6 is exhibited as well by the 40-sites system from Ref. [37].
Therefore, ED confirms that J 2 /J 1 = 0.6 marks a transition point, although in the true ground state such a transition might connect the columnar state to a quantum-disordered state. A similarly good agreement is found with the PEPS results extrapolated to the infinite size limit.
As shown in Fig. 11 (b), at small J 2 /J 1 , i.e., deep in the Néel phase, the finite size extrapolation of the ED staggered magnetization from Ref. [37] lies very close to the MSW results.
As it is well known [29], in the unfrustrated square lattice limit (J 2 = 0) the MSW value M 0 = 0.303 is only slightly smaller than M = 0.311 from ED. For the PEPS calculations an analogous quantity can -similar to section II A 5 -be derived from the peak height of the static structure factor, Eq. (5). We show its finite size extrapolation in Fig. 11 (b). In the Néel phase PEPS agrees very well with MSW theory, considerably better than ED, which decreases faster towards the strongly frustrated region. In the entire columnar phase, however, PEPS and ED data lie closer together, while MSW overestimates the order parameter.
Around the transition, however, agreement between PEPS and MSW theory is very good.
The PEPS data suggest that the magnetically disordered region, predicted by ED to occur in the range 0.35 J 2 /J 1 0.66, is either much smaller or does not occur at all.
The MSW spin stiffness ρ ≡ (ρ xx + ρ yy ) /2, however, while being finite for any considered value of the ratio J 2 /J 1 , is strongly suppressed in the region 0.3 J 2 /J 1 0. 6 [Figs. 11 (c) and (d)], suggesting as usual that accounting for quantum fluctuations beyond the MSW approximation could lead to the disappearance of magnetic order. A suppression of spin stiffness is also observed in previous results coming from ED of finite clusters [62] or from the Schwinger boson approach [64,65]. As a consequence, even though MSW admits a stable solution with magnetic order for any J 2 /J 1 value, for J 2 /J 1 = 0.6 it exhibits a clear transition from soft Néel order to a stiff columnar order, suggesting that this transition could actually separate the columnar state from a quantum disordered phase.
1. Dimer correlations in the J 1 J 2 model.
The nature of the state in the transition region between Néel and columnar order, where magnetic order is strongly reduced, can be further investigated through the study of the dimer-dimer correlations where k and l, and i and j are pairs of neighboring spins. Figure 12 sketches the expectation for the dimer-dimer correlations in (a) a columnar valence bond crystal and (b) a columnar magnetic state.
In Fig. 13, we show the spatially resolved dimer-dimer correlations from MSW theory.
Below J 2 /J 1 = 0.6 the dimer-dimer correlations have a structure compatible with a Néel state (namely they are positive and nearly equal for all bond pairs), while above J 2 /J 1 = 0.6 the dimer-dimer correlations acquire the expected structure in a columnar state, with opposite signs for the correlations between dimers of the same spatial orientation (both horizontal and both vertical) and between dimers of opposite orientations. Nonetheless, for J 2 /J 1 0.7, remarkably MSW theory shows a short-range modulation in the strength of the dimer correlations whose structure is compatible with that of a valence bond crystal. Although MSW theory is not appropriate to characterize non-magnetic states such as a valence bond crystal, it is remarkable to observe that it identifies a columnar valence-bond structure as the dominant form of dimer correlations at short range; this indication is consistent with, e.g., the results of PEPS [66], which also point towards columnar valence-bond order in the non-magnetic region of the J1J2 model.
C. Ground state properties of the J 1 J 3 model We now turn to the J 1 J 3 model. Classically, this model has a transition from Néel to spiral order at J 3 = 0.25J 1 . Recent PEPS calculations show that for S = 1/2 Néel order persists up to approximately J 3 /J 1 = 0.3 [56]. Above this point the peak of the structure factor is still at the Néel ordering vector (π, π) but its height vanishes in the thermodynamic limit, which suggests a complete loss of magnetic LRO. A different type of LRO arises anew at approximately J 3 /J 1 = 0.6 with an ordering vector Q = (q, q) that tends to (π/2, π/2) in the limit of large J 3 (see Fig. 14). For large enough J 3 the nature of the ordered phase becomes similar to that of the classical limit.
The optimization of the ordering wave-vector within MSW calculations shows that, for small J 3 /J 1 , Néel order persists up to J 3 /J 1 = 0.39 (see Fig. 14), confirming the assumption that quantum fluctuations stabilize Néel order against spiral order with respect to the classical limit. Coming from the opposite limit of J 3 ∼ J 1 , we observe a spiral phase with continuously varying pitch vector Q = (q, q), where q approaches π/2 for J 3 /J 1 → ∞, and increases up to q ≈ 0.7π for J 3 /J 1 → 0.52 + . In the region 0.39 < J 3 /J 1 < 0.52, convergence of the MSW calculations breaks down, which points at a possible spin-liquid phase, in agreement with the predictions from PEPS calculations. Fig. 15 (a) shows the PEPS energy extrapolated to the thermodynamic limit. Agreement to the MSW results is again found to be extremely good.  The indication of a disordered phase drawn from the break down of MSW theory is further corroborated by the order parameter M 0 [ Fig. 15 (b)], which decreases strongly for J 3 /J 1 → 0.39 − and for J 3 /J 1 → 0.52 + , and by the spin stiffness [ Fig. 15 (c) and (d)], which is drastically reduced when approaching the above two limits. In particular, the Gaussian spin stiffness Υ is already strongly reduced for J 3 /J 1 0.3. These results are consistent with the vanishing of the spin stiffness at J 3 /J 1 = 0.35 that was found by ED of a system of 20 sites in Ref. [67].
A destabilization of magnetic order at around J 3 /J 1 0.3 seems to be confirmed by the PEPS order parameter, Fig. 15 (b), which vanishes in the range 0.3 J 3 /J 1 0.5. Note that, again, we find that the PEPS order parameter deep in the Néel phase is similar to the MSW data, but that in the spiral phase MSW data for the order parameter lie well above the PEPS ones.
In our calculations, despite using the same equations as in Ref. [30], we find a considerably larger breakdown region. However, the region where our calculations do not yield a result is very stable, i.e., it does not depend much on system size nor on the exact algorithm for solving the self-consistent MSW equations.
The precise nature of the state in the candidate region for quantum-disordered behavior cannot be determined reliably by the use of MSW theory. From an analysis of the dimerdimer correlations in the convergence regions, we can find no indications of any exotic disordered quantum state; on the contrary, PEPS results indicate a plaquette state in the region of maximal frustration J 3 ≈ J 1 /2 [56].
D. Ground state phase diagram of the J 1 J 2 J 3 model

MSW results
After having investigated the two limiting cases of the J 1 J 2 and the J 1 J 3 models, we consider more generally the J 1 J 2 J 3 model over the relevant parameter range 0 ≤ J 2 /J 1 , J 3 /J 1 ≤ 1. As already seen in the case of the J 1 J 3 model, we observe a sizable parameter range over which the convergence of MSW theory breaks down, and which is then pointed out as a candidate region for non-magnetic behavior. We notice that, while convergence is achieved for any J 2 /J 1 ratio at J 3 = 0, a region of convergence breakdown opens up by adding a A low value marks a destabilization of magnetic LRO.
Q cl = (q, π) with continuously varying q, in favor of the columnar phase and of a potentially quantum-disordered phase.

Comparison to PEPS calculations
In Fig. 18, we display the peak height of the static structure factor, Eq. (5), from a PEPS calculation on a 8 × 8 lattice with auxiliary dimension D = 3. We observe a broad asymmetric v-shaped region in which the magnetic order, quantified by the height of the peak in the structure factor, is strongly suppressed. We notice that this region is strongly reminiscent of (albeit broader than) the breakdown region of MSW theory. In particular, the asymmetry is due to the fact that the bottom of the "v" lies at J 2 /J 1 > 0. 5 Similarly to what happens in the above spin-wave calculations, a pronounced peak at the Néel ordering vector (π, π) appears if both J 2 /J 1 and J 3 /J 1 are small, while at large J 2 /J 1 but small J 3 /J 1 the structure factor is peaked at the columnar ordering vector (π, 0). For large J 3 /J 1 , finally, the peak is located at (q, q), where q tends to π/2.

IV. CONCLUSION
In this work, we made use of Takahashi to a magnetically ordered ground state even though the true ground state is disorderedalthough in this case it will probably feature a small value for the order parameter, or a small stiffness, suggesting that the magnetic order is not robust when dealing with quantum fluctuations more accurately; on the other hand, the breakdown of MSW theory seems to be a strong indication that the true ground state is disordered.
In particular, in the case of the SATL, MSW theory completely breaks down for sufficiently weak couplings between the chains composing the lattice, suggesting that the system remains in a disordered 1D-like state even when the chains are coupled, as already predicted by recent variational approaches. A further disordered phase appears when the inter-chain couplings exceed the intra-chain ones: this phase is sandwiched in between the spiral phase of the nearly isotropic triangular lattice and the Néel phase appearing at large interchain couplings. In the case of the J 1 J 2 J 3 lattice, a large breakdown region separates the Néelordered region for small J 2 and J 3 , from the columnar-ordered region for J 2 > J 1 /2 and small J 3 , and from the spiral phase at large J 3 . Hence, a general conclusion that we can draw from the study of these two models is that collinearly ordered phases (Néel and columnar) and spiral phases cannot be connected adiabatically -at least at the MSW level -but they are always separated by a breakdown region; this is a signal that in the true ground state collinear and spiral phases might always be divided by an intermediate quantum-disordered phase.
Quantitative comparisons with more accurate methods (exact diagonalization, and variational Ansatzes based on projected BCS states and projected entangled-pair states) reveal that MSW theory with ordering wave-vector optimization goes well beyond linear spin-wave theory in dealing with quantum effects, and it correctly accounts for the quantum correction to the ordering wave-vector of the ordered phases, and for the strong suppression (or total cancellation) of magnetic order in correspondence with the candidate regions for quantumdisordered behavior. Given its flexibility and its modest numerical cost, MSW theory serves therefore as a unique tool for the identification of novel quantum phases in strongly frustrated quantum Heisenberg antiferromagnets.
where a i (a † i ) destroys (creates) a Dyson-Maleev boson at site i, S is the length of the spin, and Q the ordering vector. Here, we neglected the kinematic constraint which restricts the Dyson-Maleev-boson density n to the physical subspace n < 2S, given by the length of the spins S. Moreover, we dropped terms with six boson operators, which are of order O[n/(2S) 3 ] and are negligible for n/(2S) < 1. Using Wick's theorem [70], and defining the correlators a † i a j = F (r ij ) − 1 2 δ ij and a i a j = a † i a † j = G (r ij ), the expectation value E ≡ H can be written as (1 + cos (Q · r ij )) (A2) (1 − cos (Q · r ij )) .
After Fourier transforming, a k = 1 √ N i a i e −ik·r i , where N is the number of sites, and a subsequent Bogoliubov transformation, α k = cosh θ k a k − sinh θ k a † −k , and α † −k = − sinh θ k a k + cosh θ k a † −k , we minimize the free energy under the constraint of vanishing magnetization at each site, a † i a i = S [29]. (This guarantees that the kinematic constraint is satisfied in the mean.) This yields a set of self-consistent equations, where µ is the Lagrange multiplier for the constraint. The spin-wave spectrum reads 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 condition for magnetic LRO. It also enables Bose condensation 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 magnetic order parameter), one arrives at the zero-temperature equa- and the constraint of vanishing magnetization at each site becomes 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 F with respect to the ordering vector Q. This yields two additional equations which must be added to the set of self-consistent equations, In the SATL with NN interactions these simplify to Q y = 0 and Q x = 2 arccos − α 2 where τ 1 = (1, 0) and τ 2 = 1/2, √ 3/2 are the lattice vectors.
The values of F ij and G ij can now be calculated by solving self-consistently Eqs. (A4-A8).
Through Wick's theorem the knowledge of the quantities F ij and G ij allows the computation of the expectation value of any observable.
a. Spin stiffness The optimization of the ordering vector allows a straightforward calculation of the spin stiffness, which gives a measure of how stiff magnetic LRO order is with respect to distortions of the ordering vector, and thus provides a fundamental selfconsistency 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.
The spin stiffness tensor is defined as ρ αβ = 1 N d 2 F dQαdQ β Q=Q 0 , evaluated at the optimized ordering vector Q 0 . From this we can extract the parallel spin stiffness ρ ≡ 1 2 (ρ xx + ρ yy ) and the Gaussian spin stiffness Υ = det ρ.
Since a change in Q affects the correlators F ij and G ij , we must compute Υ selfconsistently. After finding the optimal Q 0 by the self-consistent procedure described above, we calculate 1 N F (Q x , Q y ) self-consistently for several fixed 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 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 Υ as the determinant of the partial spin-stiffness tensor.