Field-control of symmetry-broken and quantum disordered phases in frustrated moir´e bilayers with population imbalance

,

The quantum simulation of strongly correlated fermions is established in ultracold atoms, and it was demonstrated that alkaline-earth atoms in optical lattices realise strongly correlated systems with a tunable number of flavors N and SU(N)-symmetric interactions [36][37][38][39].Moiré transition metal dichalcogenides (TMD) offer recent solid-state alternatives for the controlled study of strongly correlated electron systems including triangular-lattice Hubbard models [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54].In particular, it is possible to form an SU(4) pseudo-spin out of layer and real spin degrees of freedom in twisted AB-stacked bilayers or three-layer hetero-structures with insulating middle layer [55].In an experimental realisation with WSe 2 competing electronic states with correlated insulators at integer fillings were reported [56].An important tuning parameter in these experiments is given by a perpendicular electric field which controls the layer polarisation.For integer layer populations, Mott insulators are formed at strong coupling, while at imbalanced layer population inter-layer excitonic insulators (EI) can emerge [56,57].The Zeeman effect of an in-plane magnetic field and a population imbalance in cold-atom experiments acts analogously to such a polarising field.All these fields detune the population of pairs of flavours against each other.Hence, they can be used to interpolate between effective SU (4) and SU(2) symmetric models from balanced to full polarisation.This is particularly interesting for filling factors n = 1 or n = 3, where the SU(2) limit corresponds to the half-filled Hubbard model (as opposed to a band insulator for n = 2).Theoretically, however, the effect of layer/population imbalance is not well studied.
In this work, we investigate population-imbalanced ABstacked TMD bilayers and ultracold fermionic alkalineearth atoms via the SU(4) symmetric triangular-lattice Heisenberg model in a field.We map out the phase diagram as a function of the imbalance P z and nextto-nearest-neighbour coupling J ′ employing flavour-wave theory and exact diagonalisation.In previous studies of the SU(2) symmetric case, a quantum spin liquid (of debated nature) was found between a 120 • and a stripe magnetic phase when J ′ is increased [26,27,29,33,34].In the SU(4) limit, there is evidence for a transition from a quantum liquid to four-sublattice magnetic order upon increasing J ′ [18,21,55,58,59].We show that one can tune between these two limits via an external field and determine the different ground states and their excitations in between.We find that the SU(2) 120 • antiferromagnet (AFM) develops ferromagnetic "dopants" in the minority layer and simultaneous tripartite interlayer excitonic order when the field depopulates the halffilled majority layer.For larger J ′ , we obtain an evolution from the magnetic stripe order of the SU(2) limit into a four-sublattice state of the SU(4) limit with intermediate AFM and excitonic stripes.Furthermore, we demonstrate that a large part of the phase diagram is occupied by a strongly fluctuating phase (SFP) in which quantum fluctuations prevent any long-range order, and we argue that the SFP continuously connects the candidate spin liquids of the SU(2) and SU(4) limits.

III. THE CLASSICAL GROUND STATE
We first determine the mean-field phase diagram as function of J ′ and δ and consider the role of quantum fluctuations in the next section.In the SU(4) limit δ = 0 and for J ′ = 0 the classical ground state (GS) is extensively degenerate: any state with different flavours on neighbouring sites minimises the energy.A finite J ′ selects a four-sublattice ground state out of this manifold [58].Similarly, at J ′ = 0, we expect an infinitesimal field δ to select a three-sublattice state out of the manifold because it possesses the maximal polarisation |P z | = 1/3 (see Figure 1).The three-sublattice state with |P z | = 1/3 is still extensively degenerate because any site of the third sublattice can be spin up or down.When |P z | < 1/3 the ground state keeps being highly degenerate and can be obtained starting from the three-sublattice state by substituting e.g.flavors 3 or 4 with flavors 2 or 1 with the constraints of always having different flavours on neighboring sites, as shown in Figure 1.In the SU(2) limit for large polarisation δ/J ≫ 1, the effective half-filled triangular lattice possesses 120 • antiferromagnetic order for small J ′ and transitions to a stripe phase for J ′ > 1/8 [23,24].
To obtain the mean-field phase diagram for general δ between these limits, we perform a product state ansatz |Ψ⟩ = i |ψ i ⟩ i and minimize the classical energy E cl = ⟨Ψ| H |Ψ⟩ [60].As an ansatz we choose a state with in-plane spin order and homogeneous layer polarization P z given by Phase diagram in the plane δ vs J ′ within mean-field (dashed lines δc 1 and δc 2 ) and including quantum fluctuations (solid lines).We find three phases displaying long-range order which we label with the wave-vector triplet Q t s , Qp, Q b s + Qp (see text).The blue Γ-K-K' regime describes 120 • spin and exciton order with minority ferromagnetism.In the turquoise M1-M2-M3 region, spin and excitons form stripes leading to a four-sublattice order, and in the red region, they order with incommensurate wave vectors.The green area in the phase diagram that we named strongly fluctuating phase (SFP) indicates the region where quantum fluctuations suppress the order parameter ∆12 to zero.It connects candidate spin liquids of the SU(2) (δ/J ≫ 1) and SU(4) (δ = 0) limit.(Center) Mean-field polarisation curves as a function of the external field for different values of J ′ /J.The inset shows the Maxwell construction for J ′ /J = 0.075.Dashed line is a guide to the eye for the first-order jump.(Right) Ground state degenerate manifold for |Pz| = 0.183 < 1/3 given by three different curves made of Q-vector triplets Q t s (green), Q b s (blue), and Qp (red) that minimize the classical energy.We explicitly mark three examples of triplets s , Qp around K (circles) and K ′ (triangles) are degenerate.Arrows denote how triplets evolve with θ accounting for double valued (FM) order of the top spin in the scarcely populated layer, and tripartite 120 • order of the bottom spin in the densely populated layer, which reduces to the conventional 120 • AFM order in the SU(2) limit.Simultaneously, we have a four-component inter-layer excitonic order parameter, where the non-vanishing components are given by Ô1,0 , Ô2,3 , Ô1,1 , Ô1,2 .The pairs ( Ô1,0 , Ô2,3 ) and ( Ô1,1 , Ô1,2 ) also form a 120 • configuration.At P z = 1/3, the tripartite ΓKK state with homogeneous polarisation is degenerate with the manifold of three-sublattice state with flavourpolarised sites.But for increasing |P z |, the energy of the ΓKK state remains minimal O((|P z − 1/3|) 2 J), while the energy of the flavour-polarised state rises strongly O(J) (see right panel of Figure 1).
Our numerical minimisation confirms the analytical considerations and allows us to obtain the full mean-field phase diagram between limiting cases.In total, we find three broken-symmetry states divided by two critical lines δ c1 and δ c2 (see Fig. 2).For values of the field δ < δ c1 or J ′ /J > 1/8 the energy is minimised by the wave-vector triplet and those related to this by symmetry).This gives rise to a four-sublattice order where top and bottom spin Ŝt,b form the same stripe arrangement in one lattice direction given by M 1 , and the inter-layer exciton has two components Ô1,0 and Ô1,1 forming stripes in the other lattice directions defined respectively by M 2 and M 3 .This four-sublattice order is reduced to stripes in the SU(2) limit for J ′ > 1/8 [23,24] and recovers the SU(4)-symmetric case at δ = 0 [58].When δ > δ c2 the GS is determined by the tripartite ΓKK state.For δ c1 < δ < δ c2 we find a third GS with triplets of incommensurate wave vectors.The incommensurate order is classically stable for J ′ /J ≲ 0.11.
For further characterisation, we calculate the layer polarisation as function of the external field δ for different values of the NNN exchange intensity J ′ (see Fig. 2).In the four-sublattice state the polarisation increases linearly as a function of the field , where it jumps abruptly to a higher value and the order becomes tripartite (if 0.11 ≲ J ′ /J < 1/8) or incommensurate (if J ′ /J ≲ 0.11).In the incommensurate phase the polarisation has a nonlinear behavior for δ c1 < δ < δ c2 and at δ = δ c2 the system continuously transitions into the tripartite phase.In the tripartite phase the polarisation increases linearly as |P z | = (8/27)δ +1/3.We evaluate the size of the first-order jump using Maxwell's construction (see Figure 2), which displays a non-monotonic behaviour as function of J ′ .At J ′ = 0 the polarisation jumps at δ = 0 + from zero to |P z | = 1/3.In approaching this point, the slope of the polarisation in the incommensurate region becomes steeper with decreasing J ′ so that the layer polarisability κ = dPz dδ diverges at J ′ = 0.Such an instability of the incommensurate order is accompanied by the onset of a continuous manifold of degenerate ground states for the classically forbidden values of the polarisation |P z | < 1/3 [62].It is characterised by Q-vector triplets defined on three curves Q t s (θ), Q b s (θ), Q p (θ) around Γ and K (K ′ ) in the Brillouin zone parameterised by the angle θ as shown in Fig. 2 for the representative case of |P z | = 0.183.

IV. THE ROLE OF QUANTUM FLUCTUATIONS
We study the excitation spectrum and the stability of the different phases within flavour-wave theory.Before adding quantum fluctuations to the mean-field solution, it is useful to introduce a unitary transformation U i that brings |ψ i ⟩ → |1⟩, in the same spirit as it has been done for evaluating the spin waves of the 120 degrees AFM [63] and the non-homogenous quantum Ising model [64].In this new basis, the ground state assumes the form of a homogeneous fully polarized state where every site is in the state |1⟩.The Hamiltonian in the new basis reads where We rewrite the Hamiltonian in Eq.( 2) as E cl +δH where δH contains the quantum fluctuations in terms of a generalised Holstein-Primakoff trans- n , n ∈ {1, 2, 3}, and expansion in 1/M .In the harmonic approximation, we can calculate the characteristic flavour-wave spectrum for the different phases by diagonalising δH for various values of δ and J ′ .We also explicitly checked that the dynamical structure factors contain the same excitations [61].
We first discuss the results at J ′ = 0.In Figure 3, we show four different spectra as a function of the crystalline momentum for different values of the polarization greater than 1/3.For large enough values of δ, when the system is fully polarized |P z | = 1, the spin waves of the half-filled bottom layer coincide with the ones of the mono-layer 120 • AFM [23,24,63,65].They are gap-less at k = Γ, K and their energy increases linearly in the vicinity of those points.On top of these excitations, we find two degenerate gapped "ferromagnetic"-like spin waves, which encode the hopping of one electron from the bottom to the top layer.Their gap at the Γ point is controlled by 2|δ − δ s | for δ > δ s , where δ s = 9/4J is the minimum value of the field necessary for full polarization.For intermediate values of |P z | lying in the interval (1/3, 1), we have three distinct Goldstone modes at the Γ point, two of which have a linear dispersion with different velocities while the third one displays a quadratic behavior that is associated to the FM Q t s = 0 wave-vector of the top-layer spin order.At the K-point we observe only one linear Goldstone mode.Furthermore, the FM mode is strongly suppressed by decreasing the polarization toward the critical value of 1/3.At the same time, the two other branches approach each other.At the critical point |P z | = 1/3, we obtain two degenerate excitations with linear Goldstone modes at the Γ and K points, and the FM mode completely flattens to zero.As mentioned before, these zero-energy excitations are already present in the classical picture, where any spin in the minority layer of the state with flavor-polarised sites can be flipped without energy cost, and we observe that they survive upon inclusion of Gaussian quantum fluctuations.They strongly affect the system's properties and mark the end of long-range order that for It is also instructive to analyse the role of quantum fluctuations in the classically forbidden region of the parameter space |P z | < 1/3.In this region, the classical ground state is highly degenerate.However, quantum fluctuations can remove such a high degeneracy and select particular states which become energetically favored once quantum corrections are taken into account, a well known phenomenon known as quantum order-by-disorder [66][67][68][69][70].In order to determine the ground state, we calculate quantum corrections to the energy for every degenerate classical state by evaluating the zero-point energy of quantum fluctuations [61].We find that Q t s and Q b s selected by quantum fluctuations lie respectively on the Γ-M and K-M directions.Figure 4 displays the flavour-wave spectrum as a function of the crystalline momentum for P z = −0.183.We observe that the fluctuation energy vanishes along the Γ-M direction and that we have Goldstone modes at the M-point and at two incommensurate wave vectors lying respectively on the M -K and K-Γ directions.The nodal line in Γ-M direction is again a consequence of the high degeneracy of the classical ground state and disorders the system.As we show below the presence and proximity of the zero modes strongly renormalises the order parameter and suppresses it to zero in large parts of the phase diagram.
When J ′ is finite the degeneracy of the classical ground state is removed everywhere in the phase diagram even in the case of incommensurate order.Figure 4 shows the flavour-wave spectra for J ′ /J = 0.01 and |P z | = 0.183.We find that even small values of the NNN super-exchange are enough to remove the nodal lines appearing in the classically forbidden region at J ′ = 0. Furthermore, we observe that the Goldstone mode appearing at the Mpoint for J ′ = 0 now shifts to an incommensurate vector lying in the Γ-M direction, and that the extra Goldstone mode appearing along Γ − K for J ′ = 0 acquires a gap.The regularisation of the spectra introduced by a nonvanishing J ′ allows us to quantify the impact of quantum fluctuations on the order parameter and map out the phase diagram beyond the mean-field approximation.The classical ground state in the new basis defined in Eq.( 2) is given by the completely polarised state i |1⟩ i , with n 1 = 1 and n 2 = n 3 = n 4 = 0, where n α = 1 V i ⟨S α α (i)⟩ is the α-th flavor population.With the inclusion of quantum fluctuations, the density matrix n αβ = 1 V i S α β (i) acquires off-diagonal terms and it has a block-diagonal form given by a 1×1 (α, β = 1) and a 3×3 (α, β ∈ {2, 3, 4}) block [61].Then, it is natural to choose the basis that diagonalises the the density matrix.The eigenvalues are the occupation numbers of the minority flavors that, for simplicity, we still refer to as n 2 , n 3 and n 4 , and we sort in descending order, i.e. n 2 > n 3 > n 4 .The population of the majority flavor can be computed from the knowledge of the renormalised minority ones via n 1 = 1 − 4 α=2 n α .In the generic case, the four occupation numbers are nonvanishing and we can express the order parameter in terms of three components ∆ 1α = n 1 − n α with α = 2, 3, 4. We identify the region in the phase diagram enclosed by the contour ∆ 12 = 0 (see Figure 2) as a strongly fluctuating phase (SFP), where quantum fluctuations are so strong that n 1 does not represent the occupation of the majority flavours anymore [71].
The lower panels of Figure 4 show the order parameter ∆ 1α as a function of NNN super-exchange and polari-sation.As we explained before, at |P z | = 1 we have a fully polarised layer with one fermion per site, which is equivalent to the SU(2) Heisenberg model.In this case, the GS transitions from a 120 • AFM to a striped phase at J ′ /J = 1/8 on the mean-field level [23,24].We note that the occupation numbers n 3 = n 4 = 0 in the SU(2) limit, and the order parameter has only one component, ∆ 12 .It displays a cusp at the transition point J ′ /J = 1/8 and crosses zero at two different close-by points, namely J ′ a /J ≈ 0.1 and J ′ b /J ≈ 0.138.Between these points J a < J ′ < J b the order parameter vanishes and becomes negative.This means that the harmonic approximation cannot be trusted anymore because quantum fluctuations are so strong to the point of destroying the order parameter.Interestingly, this interval is quantitatively comparable to the range of values where a spin-liquid phase was predicted by Monte-Carlo [29,33] and DMRG [26,27] calculations.
As an example, we also show the occupation numbers as a function of the external field at fixed J ′ /J = 0.075 in Figure 4.For values of P z close to complete polarisation, the system is in the Γ-K-K phase and we find n 4 = 0 so that n 1 = ∆ 14 only mildly decreases as a function of decreasing δ.At δ ∼ J, the GS becomes incommensurate, n 4 > 0 and all ∆ 1α are lowered quite rapidly by decreasing δ crossing zero (∆ 1α = 0) thereafter within the incommensurate phase.For smaller values of δ, the harmonic approximation breaks down and the system enters the strongly fluctuating regime.
We summarize the impact of quantum fluctuations in the phase diagram in Fig. 2. We observe that the strongly fluctuating regime connects the putative spin liquid phases of the SU(4) and SU(2) limits when the layer population is varied with δ.Similarly, for larger J ′ the four-sublattice and stripe phases of the limiting cases are continuously connected.In contrast, the ΓKK phase, which reduces to the 120 • AFM for full polarisation, is destroyed by strong fluctuations before the SU(4) δ = 0 limit can be reached.For small J ′ /J the strongly fluctuating regime is preempted by an incommensurate phase.

V. EXACT SPECTRA
In order to gain additional insights about the underlying physics of the exact solution we perform Lanczos exact diagonalisation for a 12-site cluster with periodic boundary conditions (PBC) defined by the lattice vectors T 1 = 4a 2 − 2a 1 and T 2 = 2a 1 + 2a 2 , where a 1 = (1, 0) and ).In Fig. 5 we show the polarisation as a function of the external field.The finite size allows for seven different polarisation values, namely P z = 1, 5/6, 2/3, 1/2, 1/3, 1/6, 0. Therefore, the polarisation displays seven plateaus.The width of a given plateau matches the energy difference between the two ground states belonging to the two different polarisation sectors, i.e. ∆E Pz = E 0,Pz − E 0,Pz−1/6 .In the thermodynamic limit, the layer-polarisability is given by κ = (lim ∆Pz→0 ∆E Pz ∆Pz ) −1 .Hence, we can use the finite difference ∆E Pz as an approximate estimate of the inverse layer polarisability (see inset in Fig. 5).We observe that ∆E −1  Pz is strongly peaked at P z = 1/3, where its value is one order of magnitude larger than for the other polarisations.This is in line with the sharpening of the slope of the polarisation upon approaching |P z | = 1/3 in the mean-field calculation (Fig. 2).
Further insights about the tendency of the system to stabilise a long range order can be gained by studying the low-energy sector of the energy eigenspectrum.Mean-field theory plus quantum fluctuations suggest that the bottom (top) layer orders antiferromagnetically (ferromagnetically) for P z ≥ 1/3.We can check the tendency of the system to form an AFM order in the bottom layer by plotting the energy eigenvalues as a function of the total spin of the bottom layer, i.e.
, which is a conserved quantity.In fact, the low energy eigenstates of antiferromagnets form a structure known as the Anderson tower of states: well separated states that are proportional to S 2 /N [58,[72][73][74].On the other hand, ferromagnets do not display the same behavior.Therefore, for P z ≥ 1/3, we expect the spectrum to only form the Anderson tower when the energy eigenvalues are plotted as function of S 2 b , and not when it is plotted as a function of the total spin.In Figure 5, we show the energy spectra for different values of the polarisation sectors.For P z = 1, the spectrum of the bottom layer reproduces the one of the SU(2)-Heisenberg model [73].At P z = 2/3, the ground state is found in the sector with S b = 0 and S t = 1: we observe that the lowenergy excited states with S t = 1 show a linear behavior as a function of S b (S b + 1).However, when the spectrum is plotted as a function of the total spin eigenvalues the linear behavior is lost as expected.Furthermore, the GS is found in the maximal spin sector of the top layer, signalling the tendency towards a FM order in the top-layer.At P z = 1/3, we start to observe a smalldeviation from linear behavior of the spectrum as a function of S b (S b +1).Furthermore, we now find the GS in the minimal top-layer spin sector, i.e. S t = 0, probably suggesting that FM order is already lost for this value of the polarisation.At P z = 1/6 we observe a sizable deviation from linear behavior, suggesting that at this value of the polarisation the bottom layer does not order antiferromagnetically.Thus, the results of the exact diagonalisation provide strong indication that the ground state is quantum disordered for small enough polarisations.

VI. CONCLUSIONS
Motivated by realisations in cold atoms and moiré TMDs, we carried out a theoretical study of the paradigmatic SU(4)-Heisenberg model on the triangular lattice in presence of a polarising field δ, which controls a population imbalance of of flavor pairs.On the classical level, the model is strongly frustrated with an extensive ground-state degeneracy, which we argued can persist for finite fields up to polarisations |P z | ≤ 1/3.Through a combination of variational mean-field calculations, flavour-wave theory, and exact diagonalisation we determined the ground states and excitation spectra for different values of the field δ and nearest-neighbor coupling J ′ .We mapped out a rich phase diagram with commensurate and incommensurate long-range orders, as well as a strongly fluctuating phase that shows evidence for a quantum disordered ground state.
For small J ′ /J and large enough δ, we found a tripartite phase where the bottom-layer spin and the inter-layer exciton order in a 120 • fashion while the top layer is ferromagnetically ordered.Accordingly, its flavour-wave spectrum shows AFM and FM excitations.For large J ′ /J, we found a four-sublattice phase for all values of δ which displays a striped configuration for top spin, bottom spin, and excitonic order parameter.In between these two phases, we found a small sliver of incommensurate order and a large regime of a strongly fluctuating phase where quantum fluctuations suppress long-range order.This is mirrored by extended zero modes in the flavour-wave spectrum for J ′ → 0. Our ED calculations on a 12-site cluster provided supporting evidence for a transition from an ordered to a disordered phase between large and small values of the polarisation.Furthermore, we observed a strong increase in the polarisability upon approaching the transition.
We argued that in the case of full polarisation, where the system effectively models the SU(2)-symmetric triangular lattice, the SFP can be identified as the precursor of the spin liquid state found in DMRG and Monte Carlo calculations [26,27,29,33,34].Similarly, the SFP coincides with a putative quantum liquid in the SU(4) limit δ = 0 [18,21,55].Interestingly, the two limits are continuously connected via the SFP making future studies on possible liquid phases in this regime highly desirable.This could also shed new light on the debated nature of the SU(2) quantum spin liquid.
Given the rich phase structure, spin and charge configurations can be effectively manipulated via the external field.The tunability of cold-atom and moiré systems provides an ideal opportunity to investigate these quantum many-body phases experimentally.The polarising field can be readily controlled and the NNN exchange J ′ is not expected to be very large so that it seems possible to reach the SFP regime, in particular, because it occupies an extended region of the phase space.For example, the softening of the FM mode, that could be probed by measuring the dynamic structure factor [61], would yield direct evidence of approaching the SFP.Intriguingly, a recent experiment on twisted AB-stacked WSe 2 reported evidence for paramagnetic insulators at hole density n = 1 for zero and full polarisation, and a potential excitonic insulator for intermediate polarisation.While this is in accordance with our findings, further experimental and theoretical investigations are needed to elucidate the nature of the insulating states.Most notably, this includes the possible emergence of Mott insulators in the associated Hubbard model.
|3⟩ + e iQ b s •Ri |4⟩ is given by where P z = ⟨ψ| Pz |ψ⟩, and τ (ρ) are (next-)nearestneighbor vectors of the triangular lattice.Generally, the NN sums are minimised by Q i = ±K with 3 τ =1 cos Kτ = −3/2, while NNN sums prefer any of the three M-vectors When δ is large so that P z → −1, the leading order terms ∝ (1 − P z ) 2 = O(1) of NN and NNN sums compete.They become equal ( when J ′ /J = 1/8 reproducing the mean-field transition from 120 • AFM to stripe phase in the SU(2) case.Thus, when we consider the case of J ′ = 0 in the P z → −1 limit, we obtain Q b s = K.Then, the nextto-leading order terms ∝ (1 − P 2 z ) = O((1 + P z )) can all be simultaneously minimized by the configuration Q P = K and Q t s = 0, which is the tripartite state we describe in the main text.The total energy of this state is E J /(N J) = 3(1 − 3P z ) 2 /8, where N is the number of lattice sites.This becomes degenerate E J = 0 with the manifold of three-sublattice states that have flavourpolarised sites at P z = 1/3.However, to slightly increase P z ≳ 1/3 in the latter state, an energy of order J is needed to flip one site from top to bottom layer.In contrast, the energy cost of the homogeneous configuration is much smaller (1 − 3P z ) 2 ≪ 1.
Considering J ′ /J > 1/8 next, the leading order term fixes Q b s = M .It is again possible to minimise all secondorder terms O((1+P z )) via Q t s = M and Q P = M ′ , where M and M ′ are two inequivalent M-vectors.The energy of this homogeneous four-sublattice state is (E J + E ′ J )/N = 2(J + J ′ )P 2 z , which becomes degenerate with the foursublattice states with flavour-polarised sites in the SU(4) limit P z → 0.

Appendix B: Inclusion of quantum fluctuations
We shall assume that, after the unitary transformation, the ground state of Eq.( 2) is very close to the fully polarised state i |1⟩ i .Hence, we can use the following approximation [58,60] for the operators: , where b † α (i) and b α (i) are creation and annihilation operators following bosonic statistics, i.e. [b α (i), b † β (j)] = δ ij δ αβ and M is the classical expectation value of S 1  1 .Substituting Eq.(B1) in Eq.( 2) we have: where H n contains bosonic operators to the power of n.Within the harmonic approximation we truncate the series in Eq.(B2) to the second order.The first term in the expansion gives back the classical energy which reads: It is clear from the last equation that δ must be proportional to M in order to be consistent with Eq.(B1).
The second term in the expansion contains linear terms in the bosonic fields and reads: This term must vanish for stability reasons and this condition fixes the value of δ that, after we set M = 1, reads: It is worth to note that the quantity on the right hand side of Eq.(B5) does not depend on the index α for symmetry reasons.
After setting M = 1, the quadratic term reads: where we used the short hand notation for b α (i + τ ) which indicates the destruction operator of a boson with flavor α at site R i + τ .After expanding the bosonic fields in their Fourier components b α (i) = 1 √ N k e ikRi b kα , we can rewrite the quadratic Hamiltonian in momentum space as following: where After substituting these expressions in Eq.(B7) we can finally write the Hamiltonian in the following matrix form: where p k = (p k1 , p k2 , p k3 ), x k = (x k1 , x k2 , x k3 ), with x kα and p kα being conjugate variables obeying the commutation relation The Hamiltonian is a 6 × 6 matrix that has been represented in a block form in Eq.(B9), where each block represents a 3 × 3 matrix.
In particular, we have that: The additive constant in Eq.(B9) is given by C = − N 2 α χ αα .The spectrum of the quantum excitations is given by the symplectic spectrum of the Hamiltonian in Eq.(B9), that coincides with the eigenvalues ϵ kα of iJ H(k), with J = −i σ y ⊗ I 3×3 [75].
Within the harmonic approximation, the expression of the density matrix is given by: where α, β = 2, 3, 4 and the subscript BZ means that the quantum expectation values for the different crystalline momenta must be averaged over the Brillouin zone.The expectation value of the majority flavors, according to Eq.(B1), is then given by n 1 = 1 −

Order by disorder
Quantum fluctuations can remove the high degeneracy of the classical ground state in the regime |P z | < 1/3 and select particular states.To investigate this quantum order-by-disorder phenomenon [66][67][68][69][70] we calculate quantum corrections to the energy for every degenerate classical state.This is done by evaluating the zero-point energy of the quantum fluctuations that in our case reads E ZP (θ, ℓ) = kα ϵ kα (θ, ℓ), where θ and ℓ = K, K ′ specify the incommensurate order, and then finding its minimum as a function of the angle.In Fig. 6 (a), we show the zero-point energy as a function of the polar angle θ for P z = −0.183: the energy has six minima at ±π/6, ±5π/6 and ±π/2.Since E ZP (θ, ℓ) does not depend on the second index for symmetry reasons, every minimum is doubly degenerate.In Fig. 6 (b), a color plot shows how the zeropoint energy is distributed on the same contours as the ones shown in Fig. 2. The Q-vector triplet minimising the energy, given by θ = π/6 and ℓ = K, is plotted with white filled dots.Interestingly, Q t s and Q b s selected by quantum fluctuations lie respectively on the Γ-M and K-M directions.For convenience, we show again the flavor-wave spectrum as a function of the crystalline momentum for P z = −0.183 in Fig. 6.

Dynamical Structure Factor
Let us define the following retarded Green's function: where we introduced the six-dimensional vector ξ k = (p k , x k ).We can express G ab (k, t) in terms of the Green's function of the quasi-particles by means of a canonical transformation that diagonalises the Hamiltonian in Eq.(B9) and that preserves the canonical variables commutation relations.This is done via a transformation S k so that S † k J S k = J and S † k H k S k = diag(ϵ k , ϵ k ), where H k is the 6×6 matrix appearing in Eq.(B9), and ϵ k = (ϵ k1 , ϵ k2 , ϵ k3 ).Such a transformation defines a set of new canonical coordinates, namely ξ k = ξk S † k and ξ −k = S k ξ−k .In particular, let us define the following Hermitian matrix Jk = i H , that can be rotated into the following diagonal form U † k Jk U k = diag(ϵ k , −ϵ k ).It is useful to introduce the following unitary transformation R k = U k T , where T = exp(iσ x π/4) ⊗ I 3×3 .Finally, we can write the canonical transformation as following: Hence, Eq.(B12) can be rewritten in the following way: The Fourier transform of G, i.e.
∞ 0 Gab (k, t)e iωt dt is a 6×6 matrix with the following structure: where α = {1, 2, 3} and the block diagonal terms are given by: The Green's function evaluated along the imaginary Matsubara frequencies can be calculated starting from Let us note that, since the canonical transformation S k is not unitary, the trace appearing in the last equation is affected non trivially by the weights [S k ] ab in Eq.(B13).For this reason it is useful to evaluate the dynamical structure factor which measures the absorption intensity.Fig. 7 shows the dynamic structure factor for different values of the layer polarisation and for J ′ = 0. We clearly recover the flavour-wave spectrum.In particular, we observe that the ferromagnetic mode that flattens to zero at |P z | = 1/3 is clearly visible for all values of the polarisation.Therefore, the softening of this mode could be used as direct experimental evidence of the highdegeneracy of the excitation spectrum which we identified as the main mechanism leading to the suppression of LRO and the onset of possible spin-liquid phases.

3 FIG. 1 .
FIG.1.Sketch of different flavour configurations for J ′ = 0 and three values of the layer polarisation.We identify |1⟩ = |↑ t⟩, |2⟩ = |↓ t⟩, |3⟩ = |↑ b⟩, |4⟩ = |↓ b⟩ for ↑/↓ spin and t(op)/b(ottom) layer.At zero field, the energy is minimised by any configuration with unequal nearest neighbours.(Left) For an infinitesimal field, such degenerate, optimal configurations can be maintained while maximising the energy gain from the field when |Pz| = 1/3.A possible optimal configuration is given by a tripartite order where flavour |2⟩ is excluded.Flipping 1 → 2 randomly at any site leaves the energy unchanged.(Center) Possible optimal configuration at |Pz| < 1/3 obtained from the previous tripartite order where in two sites (highlighted in red) flavours |3⟩ or |4⟩ have been flipped into |2⟩.(Right) High-energy state (3J) with non-homogeneous polarisation obtained from the tripartite order by flipping at one site (highlighted in blue) flavour 1 into 3.

14 FIG. 4 .
FIG. 4. (Top Panels) Flavour-wave spectra for |Pz| = 0.183 and two different values of J ′ /J = 0, 0.01.(Lower-left) Renormalised order parameter as a function of J ′ /J at full polarisation.(Lower-right) Order parameter components as a function of the external field for fixed J ′ /J = 0.075.When ∆1α < 0 (gray shaded area in the plots), quantum fluctuations destroy the mean-field order.

S 3 Pz = 2 / 3 ,Pz = 1 / 6 ,FIG. 5 .
FIG. 5. Numerical results from Lanczos Exact Diagonalisation.(Left top) Polarisation as a function of the external field.The inset displays the inverse of the plateaux width as a function of polarisation, which is an estimate of the system's polarisability (see main text).The other panels show the low-energy eigenspectra of the Hamiltonian belonging to different polarisation sectors.Antiferromagnets form the Anderson tower of states proportional to S 2 /N .Dashed lines are guide to the eyes.

FIG. 6 .
FIG. 6.Quantum fluctuations for Pz = −0.183which lies in the classically forbidden region.(a) Zero-point energy plotted as a function of the polar angle θ.(b) Energy distribution on the classical contours shown in Figure 2. (c) Spin waves spectrum plotted as a function of the crystalline momentum.