Majorana chiral spin liquid in Mott insulating cuprates

The large thermal Hall conductivity recently detected in Mott insulating cuprates has been attributed to chiral neutral spin excitations. A quantum spin liquid with Majorana excitations, Chern number +/-4 and large thermal Hall conductivity is found to be an excited state of a frustrated Heisenberg model on the square lattice. Using a Majorana mean-field theory and exact diagonalizations, we explore two possible routes to achieve this chiral quantum spin liquid, an orbital effect of an applied magnetic field and spin orbit couplings as present in cuprates. In particular, we show how only the orbital magnetic field allows this topological phase to be the ground state, while it remains an excited state of the Majorana mean field under the Dzyaloshinskii-Moriya terms. We interpret the large thermal Hall effect observed in Mott cuprates from their close proximity to a transition to a Majorana chiral quantum spin liquid which can be induced by an external magnetic field.

Introduction. The pseudogap phase of the cuprates continues to provide unexpected behavior. Thermal Hall experiments have recently found a surprisingly large thermal Hall conductivity 1 . As doping is reduced below the critical doping δ < δ * , a negative thermal Hall conductivity signal is observed, which becomes the largest when reaching the Mott insulator δ ≈ 0. The large absolute values and T -dependence of the thermal Hall conductivity in the undoped cuprate, La 2 CuO 4 , is very similar to observations in several spin liquid frustrated materials such as volborthite 2 or α-RuCl 3 3 . The half-quantization thermal Hall effect 4  Phonons have been identified as the main heat carriers in the thermal Hall effect observed in the pseudogap 5 and Mott insulating phases 6 of cuprates. The large values of the thermal Hall conductivity indicates that phonons have non-zero chirality whose origin is yet to be explained. Since magnetic impurity effects and magnons have been discarded, a possible scenario is that the paramagnetic phase of cuprates is a quantum spin liquid whose chirality is imprinted on the phonons through the spin-phonon coupling. An intriguing theoretical possibility would be a chiral quantum spin liquid on the square lattice with Majorana fermions as elementary excitations, as in the Kitaev quantum spin liquid in the honeycomb lattice, relevant to α-RuCl 3 . Other related but different chiral quantum spin liquids with either bosonic 7,8 or fermionic spinons 9 have been proposed in undoped cuprates. The d-density wave 10 has non-fractional electron constituents for δ = 0 though.
Here, we show how a chiral spin liquid state with Ma-jorana excitations is a good candidate for explaining the thermal Hall effect in Mott insulating cuprates. Using a Majorana representation of the J 1 -J 2 Heisenberg model on the square lattice, we find that a chiral spin liquid state breaking time-reversal symmetry with Majorana excitations emerges spontaneously. Such state which we denote as Majorana π-QSL has a large associated Chern number ν = ±4, leading to absolute values of the thermal Hall conductivity of the order ∼ k 2 B / as experimentally observed. However, this chiral spin liquid state is only an excited state of the system: the well known Néel, collinear and/or non-chiral disordered states being favored for different ranges of magnetic frustration. Extending the Heisenberg model to account for the presence of an external magnetic field through the orbital effect and/or for a Dzyaloshinskii-Moriya (DM) spin orbit coupling, we explore whether the π-QSL becomes the absolute ground state of the system or not. Considering several compatible choices of DM vectors, we show how the π-QSL is lowered in energy by the DM but always remains an excited state of the system. Only the orbital magnetic field is able to turn the π-QSL into the absolute ground state, consistent with observations where the magnetic field is essential to externally trigger the transition. Hence, the combined effect of the DM and the magnetic field orbital term can drive the relevant model for Mott insulating cuprates into the Majorana π-QSL phase which displays the large thermal Hall effect experimentally observed.
Model and methods. The simplest relevant model to describe the magnetic properties of undoped cuprates is the J 1 -J 2 S = 1/2 Heisenberg model on the square lattice: The bond patterns of the MMFT lowest energy QSL, the π-QSL. The arrows indicate that the bond average is positive (negative) in the same (opposite) direction of a bond between two sites. (d) Excitation spectra of the gapless π-QSL for J2 = 0 (red) and of the gapped π-QSL arising for J2 = 0 (J2 = J1/2 in this plot). The Brillouin zone with high symmetry points is displayed as an inset.
where the first sum runs over nearest-neighbor sites and the second over next nearest neighbors.
Electrons can couple to an external magnetic field B through their orbital motion. A strong coupling expansion of the Hubbard model to O(t 3 /U 2 ) leads to the additional chiral term 11 : where denotes a triangle lying in the square plaquettes of area A with vertices ijk taken in an anticlockwise direction (see Fig.1(b)). The three-spin-exchange coupling is J φ = − 24t2t 2 1 U 2 sin(2πφ/φ 0 ), with φ = BA/2 and φ 0 = hc/e the flux quantum.
A second important ingredient is the spin-orbit coupling effect through the Dzyaloshinskii-Moriya (DM) interaction 12-14 defined as: where D ij are DM vectors. The compatible vectors D ij are generically defined by 4 unit vectors d i with i = 1, 2, 3, 4 pointing in the different bond directions of a 4-site unit cell, as given in Fig.1   ). In this representation, the spin operators are S α i = i 2 b α i c i , where α = x, y, z. A three-channel mean-field decoupling of the interaction terms S α i S β j in the Hamiltonian yields to where the first three terms are associated with magnetic orders while the next three with spin liquid formation. With such decoupling, our Majorana mean-field theory 18 (MMFT) allows to analyze on equal footing quantum spin liquids (QSL)( S i = 0) as well as magnetically ordered states. To do so, as done in [18], we solve numerically the set of self consistent equations, together with a physical constraint connected to the number of particles per site that our theory has to fulfill.
Chiral QSL with Majorana excitations. The J 1 -J 2 Heisenberg model on the square lattice sustains a chiral quantum spin liquid state with Majorana fermion excitations (see Supplemental Material [26]). Among the possible QSL ansatze, the π-flux QSL ansatz or π-QSL from now on, shown Fig.1(c) has the lowest energy in a finite range of J 2 , consistent with Lieb's theorem 19 .
The eight corresponding Majorana bands -there are two sites per unit-cell 20 , labelled A and B -are given by ±3γ(k), ±γ(k) (triply degenerate) with For J 2 = 0, all these bands touch at the two Dirac points (0, ±π/2) in the Brillouin zone as shown in Fig.1(d), the π-QSL is gapless in this case. Moreover, it minimizes the total energy with E = −0.3442J 1 much lower than the 0-flux QSL ansatz with E = −0.2462J 1 . The ±π Berry phases at the Dirac cones can lead to non-trivial topology if the gap opens in a non-trivial manner. Indeed, this occurs when turning on a J 2 = 0, as shown in Fig.1(d). The gap opening lowers slightly the energy from E = −0.3442J 1 (|c 1 | = 0.4790) for J 2 = 0 to E = −0.3443J 1 (|c 1 | = 0.4777, |c 2 | = 0.0516) for J 2 = 0.5 and to E = −0.3726J 1 (|c 1 | = 0.4189, |c 2 | = 0.2700) for J 2 = J 1 . Moreover, a direct calculation of the topological invariant Chern number ν on the occupied Majorana bands 18 (with negative energy) shows that this state is a gapped topological π-QSL with ν = ±4.
In Fig. 2(a) we show the thermal Hall conductivity of the π-QSL obtained from the expression: 21 where σ xy ( ) = nk Ω z n (k)θ( − n (k)). Note that the r.h.s. of Eq. 5 is half of the standard expression used for fermionic spinons since in our MMFT approach, thermal energy is transported by the Majoranas. In the limit T → 0, we find κ xy /T → −( As temperature is raised and thermal excitations of the Majorana fermions occur, there is a cancellation between the Chern numbers of the bands above and below zero energy and κ xy /T → 0. The gapped π-QSL found in the J 1 -J 2 Heisenberg model for J 2 = 0 is doubly-degenerate since there are two possible senses of the bond amplitudes or chiralities, T ijk = 0, with the same energy. Hence, the Chern number of the π-QSL can have two values, ν = ±4, depending on the sign of the chirality as shown in Fig.2(b). It may seem that spontaneous symmetry breaking can occur since any infinitesimal J φ → 0 ± selects either the ν = 4 or the ν = −4 π-QSL solution. However, our MMFT stability analysis below reveals that the π-QSL is only an excitation of the pure J 1 −J 2 Heisenberg model. Under a finite J φ , the π-QSL with the favorable sign of the chirality for the orientation of the applied B, becomes the ground state of the system. This is consistent with experiments on the cuprates which find no thermal Hall effect when no magnetic field is applied implying no spontaneous time-reversal symmetry breaking ground state. Orbital magnetic field and spin-orbit coupling effects. We have performed fully self-consistent MMFT calcula- ). The resulting phase diagram shown in Fig.3 (a) displays Néel, stripe and spin disordered phases. For J φ = 0 we qualitatively recover the phase diagram of the J 1 -J 2 Heisenberg model on the square lattice. Subject of an intense research activity since the discovery of cuprate superconductivity, recent state-of-the-art numerical works 22-25 find a gapless quantum spin liquid and a valence bond solid around J 2 /J 1 ∼ 0.5. Our MMFT is partially consistent with this since it predicts a spin disor- Coll.
DM state π-QSL? dered valence bond crystal (VBC) phase sandwiched between the magnetically ordered Néel and collinear phases. Note that, at the mean field level, all VBCs like the columnar or the staggered phases, as well as possible resonating placquette phases 20 are degenerate all consisting of disconnected singlets. On the other hand, under sufficient strong J φ , the MMFT finds that the π-QSL becomes the ground state in a broad region of the J 2 -J φ phase diagram. Although nonphysical large J φ values are needed, it is well known that mean-field theory overestimates broken symmetry states. Indeed, we show below how quantum fluctuation effects strongly suppress the critical J φ at which the π-QSL emerges in the phase diagram.
We now analyze the effect of the DM interaction parameterised by the D vectors(see Supplemental Material [26]). The DM is also found to favor the gapped π-QSL (with Chern number ν = ±4) but a contrario to the orbital effect, it is never the absolute ground state of the J 1 -J 2 -D model. We have considered all possible MMFT self-consistent solutions analyzing in detail the competition between the π-QSL phase, the Néel and collinear magnetically ordered states and spin disordered states 20 for five different DM vector choices [16]. We have found that the disordered states are always the ground state and that the effect of DM is to continuously deform both the GS and the QSL to eventually reach, when D dominates, a magnetic order compatible with the considered D ij . The typical resulting phase diagram is shown in Fig.3(b) for YBCO. Interestingly, this phase diagram is qualitatively similar to the phase diagram of the J 1 -J 2 -J φ model discussed previously. Our MMFT results suggest that even for very large DM, D > 0.5J 1 , not present in the cuprates (D < 0.1J 1 ), the π-QSL phase with an associated large thermal Hall effect is not the absolute ground state but still a robust and close excited state. Nevertheless, it can still play an important role in effectively bringing the system closer to the π-QSL phase. Thus, from our MMFT analysis we expect that, in the presence of the DM, the large thermal Hall effect could occur at a smaller J φ than for D = 0.
Beyond Majorana mean-field theory. An important question is whether the π-QSL can occur beyond the MMFT approach. In Fig. 4 we show the quantum fidelity | ψ ref |ψ 0 | where |ψ 0 is the ground state for a given set of parameters and |ψ ref is a reference state at known limiting values of the parameters (see caption of Fig. 4 for more details). For J φ = D = 0, we identify the Néel and collinear phases a well as a possible disordered phase around J 2 ∼ 0.5−0.7. We also find that the phases shown in Fig. 4 arising between the Néel and collinear phase under J φ or under D are quite different. ED calculations of the spin structure factor(see Supplemental Material [26]) and spin chirality T ijk suggest that, while the intermediate phase is spin disordered (blue region) and chiral, T ijk = 0, in the J φ model, in the D-model it is nonchiral, T ijk = 0 and possibly also spin disordered as discussed below (green region). It however has to connect with the pure D limit at which magnetic order in-duced by D eventually sets in. Based on the consistent comparison of the ED and the MMFT phase diagrams of Fig. 3, we associate the (blue) spin disordered phase region with the MMFT π-QSL and the (green) magnetic DM with MMFT DM magnetically ordered phase.
Within the intermediate spin disordered phase of the J φ -model, close to the Néel phase, we find a topological QSL. This QSL is chiral characterized by a non-zero Chern number, ν = 1, associated with a two-fold degenerate ground state well separated from a continuum of excited states, as reported previously. 27 Interestingly this phase survives down to very low J φ for J 2 ∼ 0.6 as shown in Fig. 4. In the D-model, the evidence for a distinct and possibly spin disordered phase is corroborated by calculations combining the quantum fidelity, the gap and the spin structure factor. Indeed, in contrast to the orbital magnetic field case, the spin phase induced by D is clearly delimited by a gap closing associated with level crossings (see Supplemental Material [26]) suggesting, to the best of our knowledge, a novel intermediate phase induced by spin-orbit coupling. We let the analysis of this phase as a future work.
The MMFT could be used to unveil the nature of the chiral QSLs found in other numerical studies of the Heisenberg model on the triangular lattice with four-spin terms 28 and with three-spin orbital magnetic field terms 35 and on spin models of Kagomé Mott insulators. [29][30][31][32] Connection with cuprate materials. In Since DM is only of a few meV, D/J 0.1, we would predict that the system under no magnetic field is in the Néel phase as indeed is observed in undoped cuprates. If a magnetic field of B ∼ 10 T is applied to the system, the flux term, J φ /J 1 ∼ 10 −4 − 10 −3 would be tiny. Based on the phase diagrams obtained here, in such parameter range appropriate for cuprate materials, J 2 /J 1 ∼ 0.1 and J φ → 0, the system would be immersed in the Néel phase and so no thermal Hall effect would be expected based on our analysis. However, since the MMFT finds (see Supplementary Material [26]) that the π-QSL can coexist with the Néel AF, the π-QSL amplitude in the GS can still lead to a large thermal Hall effect. Such hybrid π-QSL + Néel AF state may also account for the anomaly observed around (π, 0) in the magnon dispersion of square lattice antiferromagnets 34 which could result from teh decay of S = 1 magnons into a continuum of Majorana excitations associated with the π-QSL.
Conclusions. We report a novel chiral QSL state with Majorana excitations and Chern number ν = ±4 which occurs as an excited state of the Heisenberg model on the square lattice. The thermal Hall conductivity of such Majorana π-QSL state leads to large absolute values |κ xy /T | ∼ (k 2 B / ) consistent with experiments in Mott insulating cuprates. MMFT predicts that this Majorana π-QSL becomes the ground state under sufficiently large J φ consistent with a topological chiral QSL found in ED. The DM present in cuprates plays a secondary role eventually suppressing critical J φ 's towards more physically realistic values. The large thermal Hall effect observed in Mott insulating cuprates can be interpreted from their proximity to a Majorana π-QSL transition which can be triggered by an external magnetic field via its orbital effect. cuprate Mott insulators Nd2CuO4 and Sr2CuO2Cl2, Nat.
Comm. 11 5325 (2020).  We first provide the details of the Majorana mean-field theory (MMFT) performed on the J 1 − J 2 Heisenberg model on a square lattice. We introduce a Majorana representation of the spin operators [1]: where α = x, y, z. Although the Hilbert space is enlarged by the Majorana representation, the actual physical space can be recovered by imposing the set of constraints over the Majorana fermions: corresponding to a single occupancy constraint on the original fermions. Any bilinear of spins S α i S β j can be mean-field Hartree-Fock decoupled as: where the first three terms are associated with magnetic ordering while the rest with spin liquid formation. This decoupling is performed on the Heisenberg terms of the Hamiltonian (α = β).The constraints in Eq. (S2) are enforced at the mean-field level through the set of Lagrange multipliers, {λ i }. Since there is no contribution from the last three terms we neglect them in our present analysis. In order to be able to describe collinear phases foursite unit cells are used in actual calculations. Here, we describe the implementation of the MMFT equations assuming a two-site unit cell allowing for Néel and/or πflux QSL states. The diagonalization of the hamiltonian is performed in momentum space by first performing a Fourier transform of the Majoranas: where we have used here, for convenience, the notation: On the other hand, the operators in momentum space, c a (k), satisfy c a (k) = c a † (−k) due to the property: c a i = c a † i which leads to standard fermionic anticommutation relations: We have: where A, B are the two sublattices in which the original lattice is divided. Since the spin operators are expressed through four Majorana fermions The hamiltonian block associated with sites in the same sublattice, a = A, B, reads: where: with: α = x, y, z.
The hamiltonian block associated with the A − B interaction reads: where: The self-consistent loop proceeds as follows. A random guess of variational parameters, ic α ia c α j b , ib α ia b α j b is arXiv:2112.14991v1 [cond-mat.str-el] 30 Dec 2021 injected in H M M F (k) and the hamiltonian diagonalized. The set of Lagrange multipliers {λ i } is fixed by the single occupancy constraint in Eq. (S2) using a least square minimization. Finally, we find a new set of variational parameters obtained from the mean-field ground state which is a Slater determinant constructed out from the filled orbitals of the mean-field hamiltonian. These three steps are repeated until convergence with the desired tolerance is reached. We now analyze the π-flux solution to the J 1 -J 2 Heisenberg model discussing the role played by Néel antiferromagnetism. We also explore other possible quantum spin disordered solutions such as the valence bond crystal (VBC). We focus on the 0 < J 2 0.5 regime which is sufficient for our purposes.

A. π-QSL solution
For the π-flux QSL bond pattern shown in Fig. 1 of the main text, the Majorana mean-field hamiltonian is: y, z for the flux pattern in Fig. 1 of the main text which has a Chern number of ν = −4. Also we have: H 00 BB (k) = −H 00 AA (k) and H αα BB (k) = −H αα AA (k).

B. π-QSL + Néel antiferromagnetism
An important aspect of the physics of undoped cuprates though, which we have ignored up to now, is Néel antiferromagnetism (AF). For J 2 = 0 we expect a Néel state in the square lattice. We have searched for hybrid π-flux QSL + Néel AF solutions. Indeed a nonzero AF moment is present for 0 < J 2 < 0.43 as shown in Fig. S1. For J 2 < 0.25 the Néel-AF is fully saturated, M = 0, 5, and so the energy follows the classical linear de-pendence with J 2 : E N eel = −0.5 + 0.5J 2 (taking J 1 = 1) shown in Fig. S1. In contrast, for 0.24 < J 2 < 0.43, a π-QSL+Néel AF i. e. a state having M < 0.5 and | c i c j | = 0 arises. A topological transition occurs at J 2 ∼ 0.43 where Néel antiferromagnetism fades away and the pure π-QSL with |ν| = 4 wins. Fig. S1 indeed shows how the energy of the π-QSL+Néel AF converges to the energy of the pure π-QSL a J 2 /J 1 at sufficiently large J 2 /J 1 .
Although our MMFT π-QSL + Néel AF is always gapped, the gap for J 2 < 0.43 associated with the Néel order is trivial while the gap for J 2 > 0.43 in the π-QSL phase is topological with Chern number, |ν| = 4. Hence, even though a hybrid π-flux QSL + Néel AF actually occurs in the MMFT equations, such state is nontopological, in fact only the pure π-QSL is topological. However, this could be an artifact of the MMFT in which broken symmetry states are treated classically.
The hybrid π-QSL + Néel AF state is not the lowest energy state in the whole 0 < J 2 /J 1 < 0.5 range. The spin disordered valence bond crystal (VBC) consisting of singlets between n.n. spins, has an energy per site of E V BC = −0.375 for any J 2 /J 1 since the singlets are effectively decoupled. As shown in Fig. S1 the VBC state wins for J 2 /J 1 > 0.25, actually the onset of the hybrid π-QSL + Néel AF state occurs. Although our total energy analysis suggests that a direct transition from a Néel AF to a VBC occurs around J 2 ∼ 0.25J 1 , it also indicates that a π-QSL + Néel AF is a possible selfconsistent solution of the MMFT. Theories beyond the Majorana mean-field treatment including quantum fluctuations may stabilize a coexistent π-QSL + Néel AF state with non-trivial topological properties.

II. THREE-SPIN MAGNETIC ORBITAL TERM
We now apply the MMFT to the magnetic orbital contribution considered in the main text: where the sum is taken over all triangular placquettes of the square lattice. This kind of three-spin terms has been discussed in the context of bosonic and fermionic spinon theories of the Hubbard model in the large-U limit. [2][3][4] Actually fermonic mean-field theories provide an exact solution of the large-N limit of the model as has been shown in the SU(N) Heisenberg model. [5] In order to make direct contact with these approaches we apply the MMFT decoupling directly onto the large-N expression for the chiral three-spin terms. Our starting point is the full expression of the three-spin contribution expressed in terms of Abrikosov fermions [2]: exchanges three fermions around the triangular placquette. Note that Einstein notation for the sums over spin indices is assumed here. The mean-field treatment of this contribution consists on contracting only the fermionic operators with the same spin which leads to the large-N expression: The fermions satisfy: α f † iα f iα = 1. It can be shown that even for the not so large value N = 2 the contractions neglected do not change qualitatively the physics.
Defining χ ij = f † jα f iα the mean-field hamiltonian, H φ reads: where the sum is over all triangular placquettes with vertices ijk. Performing a Majorana transformation: leads to the decomposition in terms of the Majoranas: We consider mean-field solutions with non-zero bond fluxes i .e. χ ij = ±iR ij with R ij ∈ R leading to the terms (S15) where the final overall sign in front depends on the final bond flux orientation. Hence, the final hamiltonian, H M F φ , will only contain diagonal contributions. Since H M F φ as well as the Heisenberg hamiltonian is diagonal in the Majoranas, we have that: The eight three-spin placquettes renormalising the Heisenberg exchange couplings between an A spin and its nearest-neighbor B and also its next-nearest-neighbor B spin are shown in Fig. S2. While the three-spin terms: are involved in the 1 − 4 and 1 − 7 n.n.n. couplings, the terms: renormalize the n.n.n. 1−6 and 1−9 exchange couplings. On the other hand the terms: renormalize the n.n. 1 − 3 and 1 − 8 couplings while: the n.n. 1 − 2 and 1 − 5 couplings.
The three-spin terms evaluated on zero-flux mean-field solutions such as a zero-flux QSL or a VBC solution leads to cancellations so that H M F −H 00 AA (k) and H αα BB (k) = −H αα AA (k). From the meanfield hamiltonian derived in Eq. (S20) it is evident that the three-spin terms just modify the bare Heisenberg amplitudes of the π-QSL solution. It is worth noting that even a pure H φ model with infinitesimally small non-zero J φ will induce a π-QSL on the square lattice.

III. MAGNETIC ORDER FROM EXACT DIAGONALIZATION
The various magnetic orders in the models considered are investigated by evaluating the ED static spin structure factor: where N is the number of sites of the cluster. In Fig. S3 we compare the S(Q) calculated on a 4 × 4 cluster for different parameters of the J 1 − J 2 − D and J 1 − J 2 − J φ models. The Néel and collinear phases display a large clear structure at their corresponding (π, π) and (π, 0) ordering wavevectors, respectively. This behavior is in contrast to the structureless S(Q) of the possible spin disordered phase (around J 2 /J 1 = 0.6) or the weak structure around (π, 0) ((π, π/2)) of the J φ (D) dominated phase. The fact that these three last phases display comparable value of S(Q) indicate that they are disordered or very weakly magnetically ordered.

IV. DZYALOSHINSKII-MORIYA CONTRIBUTION
In this section, we discuss the phase diagram of the J 1 − J 2 −D model in the presence of a spin-orbit coupling that forces the spin rotations to be coupled with real-space symmetry transformations. The Dzyaloshinskii-Moriya (DM) term is defined as: where D ij are the DM vectors. As said in the paper, we have considered fives compatible DM vector choices, YBCO, LSCO-LTO , LSCO-LTT for the ones compatible with cuprates, and the ones expected to realize rashba-like and dresselhaus-like SO as observed for electronic systems, considered by [6]. The corresponding d i vectors as displayed in Fig. 1  The effect of the DM term can be considered by employing the MMFT. Since the DM term only contains two spin operators, using Eq. S3 it is straightforward to mean field decouple the following terms: For each case, we have systematically compared the energies of the 2 most relevant ansatze in the MMFT, namely the VBC state and the π-QSL state. This is displayed in Fig. S4.
Three important features can be noted. First, thr π-QSL is always a stable anstaz of our MMFT. Second, the VBC is always lower in energy than the π-QSL, whatever the choice of the DM vectors. Finally, the π-QSL is also always lowered w.r.t. D indicating that the effect of SO coupling is positive, in the sense that it continues to favor it even though it is an excited state.
To further understand the role of the DM term, we have solved the J 1 − J 2 − D Hamiltonian by using exact diagonalizations on a 16 site cluster and the Majorana mean field theory for the case compatible with YBCO compounds. In Fig.S6, we show a ternary plot with the constraint J 1 + J 2 + D = 1, of the first gap obtained by ED. We immediately see five distinct domains separated by zero gap lines (white regions).
Along the D = 0 line, we recover the Néel region (close to J 1 = 1), the intermediate disordered phase (0.4 < J 2 < 0.6) and the collinear phase at higher J 2 . In order to characterize these phases we have calculated the quantum fidelity | ψ ref |ψ 0 | in the whole parameter space by starting from 4 reference states |ψ ref chosen deep in their corresponding expected parameter regions. They are indicated by black points in Fig.S5 for each of the 4 panels, and the value of the quantum fidelity, ranging from 0 to 1, scales from white (0) to full colored (1) hexagons respectively. Strikingly, already on such a limited cluster, the fidelity gives strong indications about the phase boundary locations. For instance, from the overlap with the pure Néel state obtained at J 1 = 1, we find that the Néel region extends in the lower right part of the ternary plot (yellow points). A very small, yet finite, fidelity can be found in the DM region as J 1 decreases, but considering that the Néel region is unambiguously characterized by a fidelity close to one then, along the J 2 = 0 line, the Néel region ends around J 1 0.4. What is remarkable here is that the Néel, the DM and the collinear regions are well separated even when the transition is second order (no level crossings). As previously mentioned, starting from a different reference state allows us to span entire regions of this model with close-to-one fidelities. This is in particular the case for the intermediate region (blue centered domain of the upper-left panel of Fig. S5) delimited by a frontier of zero energy gap as can be observed in Fig.S6. This region (the blue one in Fig.S5) is gapped and possibly magnetically disordered since the spin structure factor, S(Q), is rather featureless displaying only weak structures as shown in Fig. S3. We have denoted this intermediate phase the DM state in the main text. Finally, at large D, the system enters a fourth phase which is magnetically ordered. This can be interpreted as a quantum version of the magnetic DM state obtained with the MMFT approach. To the best of our knowledge, such intermediate possibly disordered DM state has never been reported so far.