L D C Jaubert

Owen Benton

Jeffrey G Rau

J Oitmaa

R R P Singh

Nic Shannon

Michel J P Gingras

Are Multiphase Competition and Order by Disorder the Keys to Understanding Yb 2 Ti 2 O 7 ?

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The wide interest in magnetic frustration largely stems from the diversity of unconventional phenomena it begets, ranging from anomalous Hall effect [1] to multiferroicity [2], to name only a few. The engine for this diversity is the inherent indecisiveness of frustrated magnets towards ordering, opening an avenue for new mechanisms to control their low-temperature properties. The understanding of such mechanisms not only enlarges our broad understanding of the principles via which matter organizes itself; it is also the necessary first step in order to ultimately control the exotic properties of frustrated magnets.

In the context of emergent degeneracies in pyrochlores [3,4], we present here a thorough analysis of the ordering mechanisms of a realistic microscopic model of multiphase competition. We believe our theory is the missing key to unlock many of the puzzling features of Yb 2 Ti 2 O 7 , whose ordering -or absence of -has been a matter of debate for nearly 15 years [5][6][7][8][9], discussed in the context of quantum spin liquid [10][11][12], Higgs mechanism [13] and magnetic monopoles [14,15]. Specifically, our theory accounts for the multi-step ordering process and field dependence of Yb 2 Ti 2 O 7 [16] while providing a natural setting to understand its sample dependence issue [8,17].

Starting from the most general model for anisotropic nearest-neighbor exchange on the pyrochlore lattice [19,20], we consider a range of parameters known to describe Yb 2 Ti 2 O 7 [10]. We show that a splayed ferromagnetic (SFM) phase such as the one observed in Multiphase competition of the anisotropic nearest-neighbor pyrochlore model of Eq. ( 1), as determined by classical Monte Carlo simulations for {Ji=2,3,4} = {-0.22, -0.29, 0} meV. The ψ2 and ψ3 phases are selected by order-by-disorder within the antiferromagnetic U(1) manifold and separated from the splayed ferromagnet (SFM) by a first-order transition, whose slope agrees with classical lowtemperature expansion (green line, see Appendices). Rotating all spins of a given ψ2 state by the same angle around their local easy-axes produces the entire U(1) manifold, including the ψ3 states [18]. The "splaying" angle of the spins in the SFM ground state is uniquely determined by the {Ji} parameters.

the sister compound Yb 2 Sn 2 O 7 [21,22] competes against a degenerate U(1) manifold for ordering (see Fig. 1). As in Er 2 Ti 2 O 7 , a strong contender for experimental realisation of order-by-disorder (ObD) [23][24][25][26][27], soft modes of excitations lift the U(1) degeneracy. This degeneracy lifting creates an additional competition between the so-called ψ 2 and ψ 3 states within the U(1) manifold. Using a wide range of techniques, we find that both thermal and quantum fluctuations enhance the stability of the U(1) manifold to the detriment of the SFM phase. As a consequence, quantum fluctuations bring Yb 2 Ti 2 O 7 close to the boundary between all three phases where multiple phase transitions take place. In addition to the theoretical interest of understanding how multiple ObD selections take place within a classically energetically unstable manifold [3,28,29], our work offers a natural explanation of the properties of Yb 2 Ti 2 O 7 , as recently observed in bulk measurements [16].

Model -We consider the anisotropic Hamiltonian

H = ij S i Jij S j with J =   J 2 J 4 J 4 -J 4 J 1 J 3 -J 4 J 3 J 1   (1)
where i, j are pyrochlore nearest-neighbors. All coupling matrices Jij can be deduced from J by appropriate symmetry transformations [10], with only four independent parameters {J i=1..4 } allowed by the symmetry of the pyrochlore lattice [19,20]. Our work focuses on the parameter line J 1 ∈ [-0.09 : 0] meV and {J i=2,3,4 } = {-0.22, -0.29, 0} meV, relevant to Yb 2 Ti 2 O 7 for J 1 = -0.09(3) [START_REF]The number between parentheses represents the error bar[END_REF] [10] and including the T = 0 boundary between SFM and U(1) phases at J 1 = -0.029 (see Fig. 1). We consider both quantum spins S = 1/2 and classical Heisenberg spins (| S| = 1/2) whose phase diagrams have been studied in Refs. [11, 12, 26, 31, and 32] and [3], respectively.

Classical thermal fluctuations -Starting at T = 0, the 6-fold degenerate SFM phase persists for J 1 < -0.029 before giving way to the one-dimensional degenerate U(1) manifold for J 1 > -0.029. At the boundary, new continuously degenerate ground states emerge which confer additional zero-mode fluctuations to ψ 3 configurations [3].

This zero-temperature framework sets the scene for the phase diagram of Fig. 1 computed by Monte Carlo (MC) simulations using parallel tempering [33,[START_REF] Geyer | Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface[END_REF] and over-relaxation [START_REF] Creutz | [END_REF]. The U(1) degeneracy does not survive thermal fluctuations and collapses predominantly in the ψ 2 configurations, with a small ψ 3 island around the boundary due to the above-mentioned soft modes of excitations. An interesting consequence of this ObD competition is that the ψ 2 /ψ 3 phases are precisely selected to optimize the entropy of the U(1) manifold for a given set of parameters and temperature. At finite temperature, this optimization puts the energetically selected SFM phase at a disadvantage and gives rise to multiple phase transitions for J 1 ∈ [-0.034 : -0.029].

Since the Hamiltonian of Eq. ( 1) supports a large variety of emergent degeneracies and potential ObD transitions at the boundaries between ordered phases [3,[START_REF] Benton | Classical and Quantum Spin Liquids on the Pyrochlore lattice[END_REF], we expect such a phenomenology to be a common feature of pyrochlores [29,[START_REF] Chern | [END_REF] and frustrated magnetism in general [28]. But since temperature is not the only source of fluctuations, how do quantum fluctuations fit in this picture ?

Quantum fluctuations at zero temperature -Since we know the competing classical phases (SFM, ψ 2 and ψ 3 ), it is of interest to analyze their stability in the semiclassical limit using linear spin-wave theory (LSW) [START_REF] Fazekas | Lecture Notes on Electron Correlation and Magnetism[END_REF]. However, when applied to classically unstable states, LSW usually becomes rather tedious as it requires higher order terms in the spin wave expansion. To circumvent this problem, we used the method developed in Ref. [START_REF] Coletta | [END_REF], which provides an upper bound of the semiclassical ψ 2 /ψ 3 energies for J 1 < -0.029 meV (see Appendix). Keeping in mind that this approach underestimates the stability of the U(1) manifold, LSW shows that the semiclassical T = 0 frontier is shifted by quantum zero-point fluctuations from -0.029 meV down to -0.062 meV (see Table I).

We now consider the full quantum nature of the spins. Since frustration is already at play in the constituting bricks of the pyrochlore lattice, namely the tetrahedra, exact diagonalization of a finite number of tetrahedra provides a good indication of the local influence of quantum effects. To preserve the symmetry of the pyrochlore lattice, we consider clusters of 4 and 16 spins, forming respectively 1 and 5 tetrahedra and allowing for standard ED calculations. Defining the order parameter M and correlator C = M 2 -M 2 of a given phase, then the quantity ∆C = C U(1) -C SFM is a direct measure of the SFM/U(1) competition, bringing the T = 0 frontier to J 1 ≈ -0.052(2) and -0.064(2) for N = 4 and N = 16 respectively, in agreement with the semi-classical results (see Table I).

Quantum fluctuations at finite temperature -Even if the joint analysis of thermal and quantum fluctuations is particularly challenging for such a frustrated problem, the build up of correlations when approaching the phase transition from high temperature remains accessible thanks to numerical linked-cluster computation [40][41][42] and high-temperature expansion [START_REF] Oitmaa | Series Expansion Methods for Strongly Interacting Lattice Models[END_REF] (see Appendix for technical details).

HTE confirms the shifting of the boundary down to J 1 = -0.06(3) meV and with transition temperatures lower than 500 mK (see inset of Fig. 2), in agreement with our classical simulations. As for NLC, at high temperature where quantum effects ultimately disappear, ∆C changes sign at the classical limit J 1 ≈ -0.03 meV, which can be understood from β 2 terms in hightemperature expansion. Then, as temperature decreases for J 1 < -0.03, instead of diverging towards SFM ordering, ∆C shows a clear upturn towards enhanced U(1) correlations (see Fig. 2). This upturn is adiabatically evolving to lower temperature as J 1 is decreased. NLC thus puts the U(1)/SFM frontier at J 1 = -0.070(3) meV and, together with HTE, confirms the quantum nature of the boundary shift observed in ED and LSW, but now observed at non-zero temperature (see Table I).

Multiphase competition in Yb 2 Ti 2 O 7 -The Hamiltonian (1) was shown to fit well inelastic neutron scattering data of Yb 2 Ti 2 O 7 under large field [10], allowing for the estimation of the exchange parameters {J i=1,2,3,4 } = {-0.09(3), -0.22(3), -0.29(2), 0.01(2)} meV. This set of parameters has been useful to understand the paramagnetic and high field regimes [3,10,12,13,41,[START_REF] Hayre | [END_REF]. But many questions remain open at low temperature and zero field. For instance, while single crystals show a broad peak at ≈ 185 mK in specific heat, powder samples typically display a sharp peak at ≈ 265 mK [8,9,17], suggesting the influence of structural disorder as in Tb 2 Ti 2 O 7 [45] or in Dy 2 Ti 2 O 7 [46]. But interestingly, recent bulk measurements have brought to light new generic features between powder and single crystals [16].

Based on these latter experiments [16], and thanks to the present analysis, we believe it is finally possible to flesh out a common framework for the powder and single crystals which magnetically order. In Ref. [16], both samples undergo a double phase transition with a high-temperature non-ferromagnetic transition followed by a low-temperature ferromagnetic first-order transition. Furthermore, the application of a magnetic field h does not destroy the lower-temperature ferromagnetic transition, but rather increases its temperature for h > 5 mT; the transition remaining first order up to h c ≈ 20 mT, before becoming continuous or vanishing, which is experimentally difficult to distinguish. These experimental results fit perfectly within the framework of our theory. The double transition is a direct consequence of the SFM/U(1) competition, as shown in Fig. 1 and for a similar range of temperature as in experiments. Furthermore, since the U(1) manifold is antiferromagnetic, it does not couple with h. The magnetic field thus only favors the SFM phase, and the first-order ferromagnetic transition is expected to persist until the U(1) phase is destroyed at h c . Monte Carlo simulations confirm this scenario with h c ≈ 15 mT, followed by a crossover for h > h c (see Fig. 3.b). We also mention as a further similarity the experimental presence of a reversible bump in the single crystal magnetisation at ≈ 180 mK, between the two transitions. According to our MC simulations, such feature could correspond to a ψ 2 /ψ 3 ObD transition, but with the caveat that the ψ 3 phase does not persist above 50 mK in our classical phase diagram. In that case, structural disorder may play an important role [17,45,46], since it is known to i) favor ψ 3 over ψ 2 order [48] and ii) to be stronger in single crystals than powder samples where no bump was observed [16].

It should be noted that even if the parametrization that we used, taken from Ref. [10], was done on samples different from the ones in Ref. [16], the quantum shifting of the SFM/U(1) boundary would bring this classical scenario within the experimental uncertainty of the J 1 = -0.09 ± 0.03 meV parametrization range. This quantum shifting is best illustrated in Fig. 3, where the structure factor calculated from classical simulations at J 1 = -0.0335 meV is almost identical to quantum (NLC) results at J 1 = -0.06 meV. As for neutron scattering , but rather an emphasis on the importance of quantum fluctuations. The temperature of 0.65 K in the right panels has been chosen such that the ratio between measurement and transition temperatures is the same as in Ref. [13]. The color scale is fixed from 0 to the maximum SF intensity, except for the panel at 0.3 K where the color scale was chosen for visual comparison with experiments [13]. b) When a field is applied along the [001] direction, the antiferromagnetic ψ2 phase gradually disappears in MC simulations. Our theory thus explains why the first-order transition only persists at low field in Yb2Ti2O7 [16]. We used J1 = -0.033 meV and the anisotropic g-tensor of Yb2Ti2O7 [47]: g ⊥ = 4.18 and g = 1.77.

measurements of Yb 2 Ti 2 O 7 (see Fig. 2.a of Ref. [13]), the comparison is noticeably better with classical simulations for J 1 = -0.0335 meV where a double transition takes place, than for J 1 = -0.09 meV, confirming one more time the relevance of our theory to Yb 2 Ti 2 O 7 .

Last but not least, our work brings Yb 2 Ti 2 O 7 as the missing link between Yb 2 Sn 2 O 7 and Yb 2 Ge 2 O 7 , whose ground states are respectively splayed ferromagnetic [21,22] and a yet not characterized antiferromagnet [49], which we tentatively associate with U(1). On the basis of our work, we anticipate that a natural path will take this series of compounds through a transition from SFM to U(1), achieved either under chemical pressure (Sn→Ti→Ge) or hydrostatic pressure, in analogy with spin ice [50] and Tb 2 Ti 2 O 7 [51] experiments, respectively.

Conclusion -Using a palette of complementary numerical methods, our work sets a theoretical benchmark common to both powder and single crystals of Yb 2 Ti 2 O 7 . We have shown how the multi-step ordering recently observed in bulk measurements [16] naturally arises from the competition between a SFM phase and a U(1) manifold mediated by order-by-disorder selection. On a more fundamental level, the underlying ordering mechanism should be understood as the joint consequence of both quantum and thermal fluctuations. The T = 0 frontier between the U(1) manifold and the splayed ferromagnet is shifted by quantum fluctuations, raising the interesting possibility to modulate frustration by tuning quantum fluctuations. In the general context of multiphase competition, order-by-disorder can be viewed as a freeenergy optimization process which reinforces the stability of the degenerate phase it acts upon, by naturally selecting the subset of configurations with higher entropy and/or quantum zero-point fluctuations.

In light of the numerous models and phases supported by pyrochlores, ranging from spin liquids and spin ices to (partially) ordered phases [11,[52][53][54][55][56], and subsequent boundaries between them [3,29,[START_REF] Chern | [END_REF], our present work is a paradigmatic example of why the properties of frustrated magnets should generically be understood as the sum of competing phases, rather than coming from a single controlling state. Some of these properties would indeed seem to be "coming from nowhere" in absence of a global phase diagram. We expect such competition between neighboring phases to be particularly relevant to some of the most difficult materials to characterize, such as Tb 2 Ti 2 O 7 [45] and Er 2 Sn 2 O 7 [3,57,58], and to exacerbate the sample dependence issues [17,45,46] by the proximity of phase boundaries. In that respect it is interesting to note that Er 2 Ti 2 O 7 , whose coupling parameters lie far away from any phase boundary [3,24,26], is one of the most robust rare-earth pyrochlore for reproducibility of experiments, while Yb 2 Ti 2 O 7 is essentially the antithesis. Hence, we expect the interplay between multiphase competition and disorder to become a very topical question, necessary to account for experiments in pyrochlores and frustrated magnetism, and promisingly rich in exotic physics.

APPENDIX Order parameters and correlators

As illustrated in Fig. 1 of the main text, three main phases are studied in this paper: the U(1) manifold, di-vided into ψ 2 and ψ 3 phases, and the splayed ferromagnet (SFM). Their order parameters can be found in the literature [3,24,25,29,[START_REF] Chern | [END_REF] and are repeated below for convenience.

In what follows, we adopt the convention of Ref. [10] regarding the labeling of the spins, defined by their positions relative to the centre of the tetrahedron

r 0 = 1 2 , 1 2 , 1 2 , r 1 = 1 2 , - 1 2 , - 1 2 , r 2 = - 1 2 , 1 2 , - 1 2 , r 3 = - 1 2 , - 1 2 , 1 2 . ( 2
)
The spins at the four sites are denoted as S i=0,1,2,3 = (S x i , S y i , S z i ) in the cubic coordinates, with

| S i | = 1/2.
For the U(1) manifold, the order parameter, per tetrahedron, is a two-dimensional vector.

m U(1) =     m α U(1) m β U(1)     =     1 2 √ 6 (-2S x 0 + S y 0 + S z 0 -2S x 1 -S y 1 -S z 1 + 2S x 2 + S y 2 -S z 2 + 2S x 3 -S y 3 + S z 3 ) 1 2 √ 2 (-S y 0 + S z 0 + S y 1 -S z 1 -S y 2 -S z 2 + S y 3 + S z 3 )     (3)
To differentiate between the ψ 2 and ψ 3 states within the U(1) manifold, we define the angle θ U(1) = arctan(m β U(1) /m α U(1) ); the function cos(6θ U(1) ) respectively equals to +1 and -1 for ψ 2 and ψ 3 states [3,29,[START_REF] Chern | [END_REF].

As for the splayed ferromagnetic phase, it is fully described by two three-dimensional order parameters [3]

m SFM1 =         m α SFM1 m β SFM1 m γ SFM1         =         1 2 (S x 0 + S x 1 + S x 2 + S x 3 ) 1 2 (S y 0 + S y 1 + S y 2 + S y 3 ) 1 2 (S z 0 + S z 1 + S z 2 + S z 3 )         (4) m SFM2 =         m α SFM2 m β SFM2 m γ SFM2         (5) =         -1 2 √ 2 (S y 0 + S z 0 -S y 1 -S z 1 -S y 2 + S z 2 + S y 3 -S z 3 ) -1 2 √ 2 (S x 0 + S z 0 -S x 1 + S z 1 -S x 2 -S z 2 + S x 3 -S z 3 ) -1 2 √ 2 (S x 0 + S y 0 -S x 1 + S y 1 + S x 2 -S y 2 -S x 3 -S y 3 )        
where m SFM1 is simply the global magnetization and m SFM2 accounts for the splayed nature of the ferromagnetism, i.e. the fact that the spins are not colinear. For a given SFM ground state, both m SFM1 and m SFM2 are finite.

In order to compare the growth of correlations between the two phases, we define the order-parameter correlators

C IJ = m I m J -m I m J (6)
where m I is a given order parameter. For the U(1) manifold, the use of the correlators is rather straightforward. However, for the SFM phase, since both m SFM1 and m SFM2 are finite, one needs to compute the matrix of correlators

    C SFM1,SFM1 C SFM1,SFM2 C SFM2,SFM1 C SFM2,SFM2     . (7)
By definition, this matrix is symmetric but a priori nondiagonal. Upon diagonalization, the maximum eigenvalue is kept as correlator of the SFM phase. Please note that by symmetry of the lattice, calculations can be simplified by considering only one component of the order parameters, say e.g. m α SFM1 and m α SFM2 .

Monte Carlo simulations

The Fig. 1 of the main text has been obtained by classical Monte Carlo simulations of Heisenberg spins on the pyrochlore lattice. The conventional cubic unit cell consists of 16 spins and the system size has N = 3456 spins. We used the standard Metropolis algorithm where a Monte Carlo step (MCs) is defined as N random singlespin-flip attempts. To improve the quality of the simulations, parallel tempering [33,[START_REF] Geyer | Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface[END_REF] and over-relaxation [START_REF] Creutz | [END_REF] were included in the simulations.

The error bars of Fig. 1 (main text) were obtained by using two different cooling procedures during the equilibration of the simulations. For the first "annealed" procedure, the initial configuration is chosen randomly; the system is then gradually cooled down from high temperature (fixed at 10 K) to the temperature of measurement T during 10 6 MCs; it is then equilibrated at the temperature T during 10 6 additional Monte Carlo steps. Since it is starting from high temperature, the annealed procedure tends to favor the phase with higher entropy, i.e. the U(1) manifold, and provides a lower boundary to the transition temperature. For the second "quenched" procedure, the initial configuration is fixed in the SFM phase; the system is then equilibrated at temperature T during 10 6 additional Monte Carlo steps. Since it starts in the ordered SFM phase, the quenched procedure favors the SFM phase and provides an upper boundary to the transition temperature. Following these equilibration procedures, measurements are done every 10 MCs during 10 7 MCs.

For the phase diagram of Fig. 3 b) (main text), the system size was 3456 sites with measurements during 10 6 MCs. For the structure factors of Fig. 3 a) (main text), no parallel tempering was necessary for simulations above the transition temperature. The system size was 128000 sites with measurements during 10 6 MCs.

Classical low-temperature expansion

Classical low-temperature expansion is an expansion in small fluctuations around an ordered state of classical spins. It enables calculation of the free energy of a given ordered phase up to leading term in temperature. We have used it to determine the low-temperature dependence of the phase boundary between the SFM and ψ 3 phases, for comparison with MC simulation as indicated by the green line in Fig. 1 of the main text. The method is a standard one, outlined in (e.g.) Ref. [59].

Linear spin wave theory

Linear spin wave theory (LSW) is a semi-classical method to study the stability of a classical phase in presence of the quantum zero-point energy. When the phase is a classical ground state, such as the SFM phase in the double-transition region of our paper, the method is rather straightforward. But if the phase is not a classical ground state, such as the U(1) manifold in the doubletransition region, the inclusion of zero-point energy requires a variation of the LSW theory, as proposed in Ref. [START_REF] Coletta | [END_REF] and outlined below.

At first, the approach is similar to standard spin wave expansion. We rewrite the spin operators S i in terms of Holstein-Primakoff bosons and perform a LSW , it can be diagonalized by Fourier transform followed by a Bogoliubov transformation which will return real, positive, frequencies and therefore a meaningful excitation spectrum. The semi-classical energy can then be written

E 0 semi-cl = E 0 cl 1 + 1 S + 1 2 kλ ω kλ , (8)
where E 0 cl is the classical ground state energy, ω kλ is the spin wave dispersion of band λ at wave-vector k.

However, if the chosen ordered state is not a classical ground state, then H

LSW may still vanish if the spin configuration we are expanding around is a saddle point of the classical energy. As shown in Ref. [24], this is the case for the ψ 2 and ψ 3 configurations. The last remaining point is the diagonalization of H

(2) LSW . Following the approach developped by Coletta et al. [START_REF] Coletta | [END_REF], we add to H

(2) LSW the positive definite term V = δ i a † i a i , where the a † i and a i operators are Holstein-Primakoff bosons. This additional term V does not change the classical energy and the parameter δ may be adjusted to the minimum value for which we can obtain a real, positive excitation spectrum. The energy calculated with the inclusion of V is

E δ semi-cl = E cl 1 + 1 S -N δ 2 + 1 2 kλ ω δ kλ (9)
where E cl is the classical energy of the saddle-point configuration and ω δ kλ is the spin wave dispersion calculated including the potential V . Since V is a positive definite operator, E δ semi-cl is an upper bound on the semiclassical energy of the saddle-point configuration.

Exact diagonalization

We consider the properties of the exchange Hamiltonian shown in Eq. ( 1) of the main text for a single tetrahedron (N = 4) and a single cubic unit cell of the pyrochlore lattice (N = 16) with periodic boundary conditions. Thanks to the small number of sites the full spectrum is analytically accessible for N = 4 and the lowest lying eigenstates can be easily found via standard Lanczos diagonalization for N = 16.

Numerical linked cluster expansion

At finite temperature, we have also performed a numerical linked cluster expansion (NLC) [40,42], defined as

P (L)/N = C⊂L L(C)W (C) ( 10
)
where P is some extensive quantity and N is the number of lattice sites. The sum runs over clusters C of the lattice, L(C) counts the number of clusters of type C per site and W (C) is the weight evaluated on the cluster. This weight is computed using inclusion-exclusion rule

W (C) = P (C) - C ′ ⊂C W (C ′ ) (11)
where P (C) is the quantity computed on the cluster C and the sum runs over proper subclusters of C. There is some freedom in choosing the classes of clusters to sum over in this expansion. We follow the approach of Ref. [41] and use tetrahedra as our building block. A linked cluster with n T tetrahedra will have at most 3n T + 1 sites so we are limited to n T ≤ 4 in the expansion. The properties P we will compute are defined on the tetrahedron and so conform well to this expansion.

For convenience, we reproduce in Table II the geometrical clusters used in Ref. [41], with the appropriate embedding constant. These include the 0 th order point up to 4 th order which includes clusters composed of four tetrahedra. Each cluster has an associated Hamiltonian H C obtained from the exchange Hamiltonian (Eq. (1) of the main text) which we diagonalize numerically. The largest n T = 4 clusters we consider have 13 sites and thus Hilbert spaces of dimension 2 13 = 8192. These remain amenable to full diagonalization and thus we can compute arbitrary thermodynamic quantities at finite temperature.

The series is organized into terms P n including up to n tetrahedra. Explicitly carrying out the expansion using the embedding constants one has

P 0 = +P (C 0 ) (12)
P 1 = -P (C 0 ) + 1 2 P (C 1 ) (13)
Note that P (C 0 ) does not appear directly past first order. Following Ref. [41] the Euler resummation method is used on the final two terms to accelerate convergence. If we define the differences S n = P n+1 -P n then the Euler approximants are given by

E 2 = P 2 (17)
E 3 = E 2 + 1 2 S 3 (18)
E 4 = E 3 + 1 4 (S 3 + S 4 ) (19)
The difference between 3 rd and 4 th order of expansion is the uncertainty of our NLC computations.

High Temperature Expansion

The High Temperature Expansions (HTE) are done up to order β 8 , as explained in the book of Ref. [START_REF] Oitmaa | Series Expansion Methods for Strongly Interacting Lattice Models[END_REF]. The series expansion are then analyzed using Padé approximants, i.e. based on rational functions. We constructed all near-diagonal Padé approximants with 8 or 7 terms in the series, i.e.

FIG. 1 .

 1 FIG. 1.Multiphase competition of the anisotropic nearest-neighbor pyrochlore model of Eq. (1), as determined by classical Monte Carlo simulations for {Ji=2,3,4} = {-0.22, -0.29, 0} meV. The ψ2 and ψ3 phases are selected by order-by-disorder within the antiferromagnetic U(1) manifold and separated from the splayed ferromagnet (SFM) by a first-order transition, whose slope agrees with classical lowtemperature expansion (green line, see Appendices). Rotating all spins of a given ψ2 state by the same angle around their local easy-axes produces the entire U(1) manifold, including the ψ3 states[18]. The "splaying" angle of the spins in the SFM ground state is uniquely determined by the {Ji} parameters.

FIG. 2 .

 2 FIG. 2. The difference of correlators ∆C = C U(1) -CSFM computed with NLC confirms the quantum shifting of the boundary towards more negative values of J1 than for the classical system, estimated at J1 = -0.070(3) meV. J1 is given in the caption while {Ji=2,3,4} = {-0.22, -0.29, 0}. Inset: ∆C as computed from HTE for J1 = -0.09 (red), -0.06 (cyan) and -0.03 (black). The breadth of each curve represents the uncertainty for NLC and HTE at low temperature.

FIG. 3 .

 3 FIG. 3. Application to Yb2Ti2O7: a) structure factor as measured by neutron scattering and b) phase diagram in a field. a)The spin-flip (SF, top panels) and non-spin-flip (NSF, bottom panels) calculated by quantum NLC and classical Monte Carlo are almost identical when J1 is shifted from -0.0335 meV to -0.06 meV (please note that the temperature has been renormalized by 5/4 for a better agreement). This agreement confirms the quantum shifting of the boundary at finite temperature. When approaching the phase transition, the MC structure factor for J1 = -0.0335 meV reproduces the characteristic features of Yb2Ti2O7 neutron scattering data measured at 300 mK (to be compared with Fig.2aof Ref.[13]). On the other hand, the comparison between experiments[13] and classical simulations using the parametrization of Ref.[10] (J1 = -0.09 meV) are noticeably less successful, especially around (220)[3]. This disagreement is not a criticism of the parametrization by Ross et al., but rather an emphasis on the importance of quantum fluctuations. The temperature of 0.65 K in the right panels has been chosen such that the ratio between measurement and transition temperatures is the same as in Ref.[13]. The color scale is fixed from 0 to the maximum SF intensity, except for the panel at 0.3 K where the color scale was chosen for visual comparison with experiments[13]. b) When a field is applied along the [001] direction, the antiferromagnetic ψ2 phase gradually disappears in MC simulations. Our theory thus explains why the first-order transition only persists at low field in Yb2Ti2O7[16]. We used J1 = -0.033 meV and the anisotropic g-tensor of Yb2Ti2O7[47]: g ⊥ = 4.18 and g = 1.77.

  The authors would like to thank Han Yan for collaboration on related projects. LDCJ, OB and NS are supported by the Okinawa Institute of Science and Technology Graduate University. RRPS is supported in part by NSF-DMR 1306048. The work at the U. of Waterloo was supported by the NSERC of Canada, the Canada Research Chair program (M.G., Tier 1) and by the Perimeter Institute (PI) for Theoretical Physics. Research at PI is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation.

P 2

 2 1 ) + 10P (C 2 ) -21P (C 3 ) +3P (C 4a ) + 6P (C 4b ) + 2P (C 4c )

2 TABLE

 2 II. Clusters used for the NLC expansion. A graphical representation is shown along with the embedding constant L(C).

  [4/4], [5/3], [3/5], [6/2], [2/6], [4/3], [3/4], [5/2], [2/5] where [m/n] stand for the powers of the polynomial in the numerator and the denominator. The resulting spread in Padé approximant values represents the error bars of our HTE computations.

  1/S expansion around the local spin configuration of the chosen ordered state. At the harmonic level, our Hamiltonian takes the form H ≈ H LSW = H

	(0) LSW +H LSW +H (1) LSW . H (2) LSW is simply (0)
	the classical energy of the ordered state around which we
	are expanding, while H LSW and H (1) LSW contain only linear (2)
	and quadratic terms in boson operators, respectively.
	If the chosen ordered state is a classical ground state,
	then H LSW must vanish. As for H (1)	(2)


L D C Jaubert

R Moessner

Multiferroicity in spin ice: towards a magnetic crystallography of Tb 2 Ti 2 O 7 in a field

We combine two aspects of magnetic frustration, multiferroicity and emergent quasi-particles in spin liquids, by studying magneto-electric monopoles. Spin ice offers to couple these emergent topological defects to external fields, and to each other, in unusual ways, making possible to lift the degeneracy underpinning the spin liquid and to potentially stabilize novel forms of charge crystals, opening the path to a "magnetic crystallography". In developing the general phase diagram including nearest-neighbour coupling, Zeeman energy, electric and magnetic dipolar interactions, we uncover the emergence of a bi-layered crystal of singly-charged monopoles, whose stability, remarkably, is strengthened by an external [110] magnetic field. Our theory is able to account for the ordering process of Tb2Ti2O7 in large field for reasonably small electric energy scales.

By providing mechanisms for strong magneto-electric coupling, frustration has become a key ingredient in multiferroics [1][2][3][4][5][6][7][8]. While the search for high-temperature multiferroics is appealing for technological application such as memory devices [8], frustration opens a window on novel fundamental properties of magnetic matter at low temperature where even weak perturbations can play an important role. This holds especially for the collective behaviour of spin liquids in a wide range of compounds from rare-earth [9] and copper [10][11][12] oxides to organic Mott insulators [13] or iridates [14,15].

In spin ice materials, the constraints imposed by frustration support an extensively degenerate ground state where magnetic fluxes are locally conserved [16]. Such flux conservation can be described as a divergence-free condition, categorizing the spin ice ground state as a Coulomb spin liquid by analogy with Maxwell's electromagnetism [17][18][19], where excitations take the form of classical magnetic monopoles (Fig. 2.c-d) [START_REF] Castelnovo | [END_REF].

In addition to their magnetic properties, it has been recently theorized that magnetic monopoles could also carry an electric dipole moment [21] (Fig. 2.c).

Here we shall investigate the multiple facets of such magneto-electric coupling, as an unexplored generic ordering process in rare-earth pyrochlores, able to lift the degeneracy of spin liquids and to manipulate topological excitations in frustrated magnets. Our results are double. First of all, we give a precise description of the mosaic of competing phases in our multiferroic spin ice model (Eq.( 1)). We show how a ferromagnetic double-layer structure of monopoles (DL) is stabilized by electric dipolar interactions and even enhanced by a magnetic field in the [110] direction. We then use this double-layer structure as a signature of multiferroicity in rare-earth oxides able to account for recent experiments on the spin liquid candidate Tb 2 Ti 2 O 7 in a field.

The model -we consider classical Ising spins S aligned with their local easy-axes on the pyrochlore lattice supporting electric moments P induced by magneto-electric coupling [21] (Fig. 2), interacting via nearest neighbour spin coupling, Zeeman energy, magnetic and electric dipolar interactions where i, j and α, β are respectively indices for magnetic spins on the pyrochlore lattice and electric dipoles on the diamond lattice. r m and r e = 3/2 r m are the respective nearest neighbour distances. The nearest neighbour vector e ij , magnetic S i and electric P α moments have unit length. The size of the moments µ and p is included in the energy-scale prefactors

H = J ij S i • S j -h • i S i (1)
+ D m r 3 m i>j S i • S j -3 S i • e ij S j • e ij
D m = µ 0 µ 2 4π r 3 m , D e = p 2 4πε 0 r 3 e , h = µ 0 µ H (2)
where µ 0 and ε 0 are respectively the vacuum magnetic permeability and electric permittivity and H is the external magnetic field.

The Hamiltonian is studied via classical Monte Carlo simulations, using parallel tempering, worm and singlespin-flip Metropolis algorithm. The dipolar energies have been computed with the Ewald summation [22,23], in absence of demagnetization factor in order to develop a sample-independent theory [24]. All spin configurations are given in the Appendix.

Monopole Double Layer -first of all, what happens in absence of magnetic interactions, i.e. D m = h = J = 0 ? We find that electric dipoles induce a bi-layer structure of single charges with zero polarization and saturated magnetization along the [110] axis (Fig. 2.a). Because the electric field is even under time reversal, the apparition of such magnetization has to be spontaneous. This configuration is unfrustrated at the nearest-neighbour level, whose contribution constitutes 96% of the total energy. To our knowledge, magneto-electric coupling is the first intrinsic mechanism favouring single magnetic charges down to zero temperature.

Local chemical potential Jin terms of monopoles, J plays the role of a chemical potential favoring the Coulomb spin liquid for J < 0 (vacuum of charge) and the AIAO double-charge crystal for J > 0 [25,26], making single charges gapped topological excitations in both cases. However, the previously observed J = 0 doublelayer structure turns out to be robust over a large range of values for J/D e ∈ [-2.08 : 0.69] (Fig. 1), raising the question on the nature of the mechanism able to stabilize such "excitations".

For instance when J < 0, creating a pair of single charges out of the Coulomb spin liquid costs |4J/3| while the energy gain is at most -2D e /3, making such monopole-pair creation unfavorable for J/D e < -1/2. An energetically stable cluster of bi-layered monopoles thus needs to get bigger and bigger as J decreases in order to minimize its surface-over-volume ratio. The need for this kind of nucleation process to seed and grow a cluster makes the transition first order and prevents full thermalization of the simulations in the vicinity of the extensively degenerate Coulomb spin liquid (see yellow/green hatched region in Fig. 1). As a consequence, for J ≈ -2D e , an experimental cooling down protocol would probably fall out-of-equilibrium; once the magnet enters the Coulomb spin liquid with a low density of monopoles, it will be difficult to nucleate a big enough cluster of magnetic charges to crystallize the double-layer structure. Such phenomena also exist for J > 0, but to a lesser extent because of the low entropy of the AIAO ordered phase (see red/green hatched region in Fig. 1).

Hence, the magneto-electric opportunity to stabilize monopole excitations comes at the cost of large (free) energy barriers and multiple metastable states, which can naturally account for strong out-of-equilibrium effects in pyrochlores. In order to build a comprehensive and experimentally relevant picture of the problem, let us now include magnetic dipolar interactions, before adding a magnetic field able to tune these energy barriers, and finally applying our theory to experiments.

Long range magnetic dipolar interactions D m -since the electric polarization is coming from magneto-electric coupling, it would be improper to neglect the D m energy scale, especially if we keep in mind rare-earth materials with potentially large magnetic moments. In spin ice, magnetic dipolar interactions are responsible for the effective Coulomb interactions between monopoles [START_REF] Castelnovo | [END_REF]. The property of "projective equivalence" [17] ensures the quasi-degeneracy of the Coulomb spin liquid in presence of magnetic dipolar interactions [23], which is only weakly lifted in favour of the 2 in -2 out long range ordered dipolar spin ice (ODSI) state for T ≪ D m [27,28]. This is why, by favoring 2 in -2 out configurations, opposing the proximity of same-charge monopoles and hindering long-range ferromagnetic order, the D m interaction seems very unfavorable to the DL phase. Indeed for negative J, the double-layer phase makes way for the ODSI at finite D m (Fig. 3). Remarkably, the projective equivalence, valid for magnetic dipolar interactions in the Coulomb spin liquid but not electric ones in the DL phase, makes the transition temperature one order of magnitude smaller from the DL to the ODSI phase.

However the nearest neighbour contribution of the D m term decreases the monopole chemical potential which is four times bigger for double charges than for single ones [START_REF] Castelnovo | [END_REF]. The counter-intuitive consequence is that magnetic dipolar interactions favour DL order for positive values of J where it was absent at D m = 0. Because magnetic Coulomb interactions in spin ice are four orders of magnitude smaller than between bare electric charges at the same distance, a relatively low-energy coupling such as D e is sufficient to counter-balance the repulsion between neighbouring magnetic charges, opening the path for novel crystal structures made of magnetic monopoles. A [110] magnetic field hin the absence of magneto-electric couplings, such a field polarizes half of the spins along the α-chains (Fig. 2.a) while magnetic dipolar interactions align the remaining spins along antiferromagnetically-ordered β-chains (see Appendix) [29][30][31]. This 2 in -2 out state is noted ODSI [110] .

By ordering the α-chains the field hinders thermal fluctuations which improved the thermalization of simulations, even if the presence of modulated phases could not be completely ruled out for all boundaries of the fourdimensional parameter space. The T = 0 phase diagram of the simulated phases was then computed by Ewald summation (see Appendix and Fig. 4). Because of the intrinsic quasi-degeneracy of the 2 in -2 out configurations, the boundary with the DL phase only barely shifts as the ODSI order quickly gives way to ODSI [110] upon increasing h. The AIAO phase, on the other hand, is suppressed by the [110] But is Tb 2 Ti 2 O 7 a good candidate for our theory ? We believe so for three reasons. Firstly, the Ising anisotropy of Tb 3+ ions [34,35] and the pinch points observed in polarized neutron scattering [36][37][38] are strong indications for underlying spin-ice physics. Also, Tb 2 Ti 2 O 7 possesses a giant magnetostriction [39,40], especially along the [110] direction [41], and the stability of the low temperature spin liquid phase has been recently ascribed to spin-phonon hybridization of the excitations [42], via dynamical Jahn-Teller coupling [43]. Last but not least, while sample dependence seems to be an issue in this compound, the double-layer structure in a large [110] field has been confirmed by independent experiments [32,33] and is thus a robust feature of Tb 2 Ti 2 O 7 . It should indeed be noted that the nature of the zero-and low-field phases of this material remain under debate -possibly because of light stuffing/dilution [44] -alternatively described as spin liquid [36,37,[45][46][47] or glassy [48][49][50] with antiferromagnetic correlations [50,51]. Such behaviour is reminiscent of another rare-earth pyrochlore, Yb 2 Ti 2 O 7 : while being remarkably well parametrized under a high magnetic field [52], the zero field properties of Yb 2 Ti 2 O 7 noticeably vary between samples [53][54][55][56], also possibly due to light stuffing [54]. Given the present low-field uncertainty, our goal here is to propose an alternative scenario for the high field region and to put a new benchmark on the 15-year-old puzzle that is Tb 2 Ti 2 O 7 .

Because of a complex single-ion crystal field [47,[57][58][59], the effective size of the Tb 3+ magnetic moments is not fixed, but µ = 6 µ B is a good estimate at low temperature [34] and high field [60]. With r m = 3.59 Å, Eq. 2 gives D m = 0.48 K. From our theory, the double-layer structure can then be stabilized for D 0 e > D m /1.66 ≈ 0.30 K, corresponding to an electric moment of ∼ 2.10 -31 C.m and a displacement of oxygen ions of ∼ 0.6 pm, which are reasonable estimates for multiferroics [8,21]. Using the parametrization of [61] with J = 2.7 K, our multiferroic spin ice model can also explain why the DL phase only appears at finite field [33,50] (Fig. 5). In light of the sample-dependence issue, it is difficult to push the comparison further to low field, where the phase might be antiferromagnetic but probably not AIAO [50], and should be separated from the DL structure by a collective paramagnet up to 2 Tesla. This is where magnetostriction comes into play as a potential complementary facet of our model.

The [110] direction for a magnetic field has been shown to maximize magnetostriction in Tb 2 Ti 2 O 7 , resulting in field-dependent oxygen displacement [41]. Hence, an external magnetic field stabilizes the DL structure not only by suppressing the AIAO order (Fig. 4), but possibly also by increasing the electric energy scale D e . In that case, an even wider range of parameters (e.g. J = 0.96 K [62]) and perturbations will lead to a double-layer phase at high field. The possibility to include disorder, anisotropic [43,63,64] or next-nearest-neighbour interactions [65] and quantum fluctuations [65][66][67] to our multiferroic spin ice model opens a rich diversity of potential phases to account for the yet uncertain phase of Tb 2 Ti 2 O 7 in low field. In particular it is tempting to speculate whether the glassy behaviour observed in some samples [48][49][50] might be a consequence of the dynamically difficult nucleation process discussed in this paper, especially if Tb 2 Ti 2 O 7 lies close to a spin liquid phase [36,37,[45][46][47]68]. spin liquid by creating interactions between topological excitations. In spin ice these excitations condense into a bi-layered monopole crystal strengthened by a [110] magnetic field, a non-trivial example of "magnetic crystallography". Our theory offers a simple and robust explanation for the ordering of Tb 2 Ti 2 O 7 in a large [110] field for a reasonably small electric energy scale D e .

Conclusion -in summary, we have shown how magneto-electric coupling can lift the degeneracy of a

If we look at the diverse physics emerging from itinerant electrons coupled to spin ice (anomalous and spontaneous Hall effects in Nd 2 Mo 2 O 7 [69] and Pr 2 Ir 2 O 7 [70], non-Kondo resistivity minimum [71][72][73][74] and a new kind of quantum criticality [75] in iridates), we should expect the coupling to an additional, ferroic, degree of freedom to bring a new flavor to spin ice and spin liquids, both at equilibrium and dynamically [76]. Experiments on Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 already suggest the presence of magneto-electric effects [77][78][79][80], which could be enhanced or even qualitatively modified by doping and chemical pressure [81][START_REF] Lin | 12th Joint MMM-Intermag Conference[END_REF][START_REF] Xu | [END_REF] or with the inclusion of an electric field. More generally, multiferroicity offers a promising mechanism to control topological defects in magnets. We hope our work will motivate further theoretical and experimental investigations of multiferroic effects in pyrochlores and spin liquids.

The authors are thankful to Pascal Quémerais for collaborations at an early stage of this project, and to Owen Benton, Bruce Gaulin, Isabelle Mirebeau, Yukitoshi Motome, Karlo Penc and Oleg Tchernyshyov for useful discussions. This work was supported by funding from the Theory of Quantum Matter unit of the Okinawa Institute of Science and Technology Graduate University.

APPENDIX Ground states energies

The T = 0 phase diagram of Figs. 1 and4 in the paper was calculated based on the following energies per number of spins N . The spin configurations are given in figures 8, 9 and 10.

AIAO:

E = +4.09 D m -J (3)
DL:

E = 0.0455 D m -h/ √ 6 -0.692 D e (4) ODSI: E = -1.95 D m + J/3 (5) ODSI [110] : E = -1.90 D m -h √ 6 + J/3 (6)
ML:

E = -0.370 D m -h √ 6 ( 7
)
One sees immediately that the DL phase wins over the ML one for D e > 0.60D m , and that it is stable for J/D e ∈ [-2.07; 0.692] when D m = h = 0. While the prefactors for the coupling J and the Zeeman term are straightforward to calculate, the remaining terms require Ewald summation for a precise estimate. However, many of them can actually be calculated analytically to a good approximation [START_REF] Castelnovo | [END_REF]84]. Since it provides a useful insight into the analogy between dipoles and monopoles, we briefly explain the method in the next section. Please note these results are obtained in absence of demagnetization factor. We recover the same values as Yoshida et al. [30] for the magnetic energies of ODSI and ODSI [110] when including the demagnetization factor of a sphere in vacuum.

Analytical calculation of the energies

Let Q = n q m denote the magnetic charge on a given diamond site where q m = 2µ/r e and n ∈ {-2, -1, 0, 1, 2}. In presence of magnetic dipolar interactions only, the energy cost p n to create a monopole of charge n is (see Supplementary Informations of [START_REF] Castelnovo | [END_REF])

p n = - 8 3 1 + 2 3 n 2 D m (8)
The magnetic interaction between charges is difficult to calculate for any random configuration, but if the system is charge ordered then it is possible to use the Madelung constant of the corresponding crystal structure in order to calculate its Coulomb energy [84]. Both AIAO and ML configurations are ordered in the zincblende structure with Madelung constant α zb = 1.638. With N/2 diamond sites, the Coulomb energy is

U c n = - 1 2 N 2 α zb n 2 µ 0 q 2 m 4πr e = -N α zb 2 3 2 3 n 2 D m . (9)
giving a total energy

E n = U c n -p n N 2 E n = 2N 3 2 + (2 -α zb ) 2 3 n 2 D m (10) = 1.53 n 2 D m (11)
where n = 1 for the ML and n = 2 for AIAO. This value should be compared with a vacuum of charges, i.e. the Coulomb spin liquid. Since the degeneracy of the Coulomb spin liquid is weakly lifted, the choice for a reference energy is somewhat arbitrary. With a lowest energy state at -1.95D m (ODSI) and the highest energy one at -1.85D m (fully saturated in the [001] direction), we choose E ref = -1.90D m as reference energy and obtain

∆E 1 = E ML -E ref = 1.53D m vsE 1 = 1.53D m ∆E 2 = E AIAO -E ref = 5.99D m vsE 2 = 6.12D m
which are in remarkably good agreement. As for the double layer structure, we could not find the Madelung constant for such charge ordering. So we calculated it using Ewald summation for Coulomb interactions and obtained α DL = 0.976. This gives E DL,Madelung = 1.89D m , to be compared with ∆E DL = E DL -E ref = 1.95D m , within 3% of error.

Finite temperature simulations

Now that the zero temperature boundary can be determined exactly from equations (3) to (7), let us turn our attention to the finite temperature phase diagram, and in particular to the double-layer phase. All simulations were done with parallel tempering, usually with 1 mK difference between parallel temperatures. To further help thermalization, we have developed a variant of the worm algorithm for dipolar spin ice [24], adapted for both electric dipolar interactions and the presence of singly charged monopoles. Measurements were made for a given number of Monte Carlo steps (MCs) noted t max .

We used two different equilibration processes to compute the error bars:

• (i) the system is slowly cooled down from high temperature to the temperature T of measurement during t max /10; then it is thermalized at temperature T during t max /10; for our model, this method provides a lower bound for the transition temperature.

• (ii) for any given set of parameters, the ensemble of ground states is known exactly (cf. the previous sections); the system is then quenched into one of the ground state configurations and thermalized at temperature T during t max /10; this method is biased towards ordering and offers an upper bound for the transition temperature. These two values provide the error bars plotted on Figs.

1 and 5 of the main text for a system of size N = 1024 spins. When not visible, the error bars are smaller than the symbols. For the 3-dimensional plot of Fig. 3 of the main text (no error bars), we only used the equilibration process (i) for a system of size N = 432, which provides a lower estimate of the stability of the DL phase. For Figs. 1, 3 and 5 of the main text, we used t max = 10 6 MCS, except in the double layer phase of Fig. 1 where t max = 10 7 MCS.

In Fig. 6, we show how these two processes converge to the same value for four different sets of parameters ordering in the double-layer phase (the most difficult ordering process in our simulations). Because of the double longrange interactions and the very long time of thermalization, big system sizes are difficult to simulate: Fig. 7 displays how the transition temperature converges to a finite value with increasing linear system size L.

Husimi tree

However close to the boundaries, especially with the Coulomb spin liquid, if the equilibration process (ii) always orders in its ground state configuration, the process (i) might not be able to find the true ground state and will be dominated by the neighbouring phase. This is what happens in the hatched regions of Fig. 1 in the main text. In that case, we cannot rely solely on simulations to determine the finite temperature phase diagram and an analytical approach becomes necessary.

The out-of-equilibrium region between the DL and AIAO is rather narrow, which is why we shall focus on the broader one between the Coulomb spin liquid (CSL) and the double layer, by estimating the free energy of the two phases for D m = h = 0 (cf. yellow/green hatched region in Fig. 1

of the main text).

Let us first consider the double layer phase. According to Eq. 4, its internal energy is U DL = -0.692D e . As for the entropy S DL , since the transition is strongly first order, fluctuations can be neglected in a first approximation when compared with the Coulomb spin liquid of extensive degeneracy. Thus the free energy is F DL = U DL -T S DL = -0.692 if we fix D e = 1.

As for the Coulomb spin liquid, since D m = h = 0, it corresponds to the canonical nearest neighbour spin ice model with electric interactions only between singlycharged monopoles. To obtain an upper and lower estimate of the boundary between the Coulomb spin liquid and the DL phase, we consider the two following cases

• 1) electric interactions between the dilute monopoles are modeled by an effective chemi-cal potential: the biggest interaction energy gained by a pair of monopoles being -2D e /3, we estimate an effective chemical potential of -D e /3 per monopole.

• 2) electric interactions are neglected: since they tend to lower the interacting energies of the dilute monopoles, this should give an upper estimate of the transition temperature.

The free energies of both cases can be calculated in the Husimi tree approximation All tetrahedra carry a single magnetic charge. The green arrows are the electric polarization carried by each monopole. The thick lines are the α-chains of spins along the [110] direction. There is no global polarization but there is a saturated magnetization along the α-chains in both configurations. For a 3 in -1 out configuration (resp. 3 out -1 in), the minority spin is the outward (resp. inward) spin. In the ML phase, each minority spin (circled in orange) is a minority spin for both adjacent monopoles.

F CSL 1 = - T 2

r 3 ij+ D e r 3 3 αβFIG. 1 .

 3331 FIG.1. When Dm = h = 0, the electric dipoles stabilize a monopole double layer (DL -green, see Fig.2.a), in competition with all in / all out order (AIAO -red), and the Coulomb spin liquid (yellow). The circles (crosses) are the transition (crossover) temperatures obtained from Monte Carlo simulations. In the hatched regions, even if T = 0 calculations confirm the energetic stability of the double-layer structure, the first order nature of the transition prevents full thermalization of the simulations. The solid lines are upper and lower mean field estimates of the boundary. See Appendix for details on simulations and calculations.

FIG. 3 .

 3 FIG. 3. Phase diagram of J/De and Dm/De obtained from Monte Carlo simulations. The vertical axis is the normalized transition temperature into the AIAO (red), double-layer (green) or ODSI (orange) phase. See Appendix for details.

FIG. 4 .

 4 FIG. 4. Zero temperature phase diagram in a [110] field h obtained from Ewald summation. The arrows show the evolution of the boundaries as h increases. The DL phase is further stabilized by h over the AIAO order, but restricted to small values of Dm by the apparition of a monolayer phase of monopoles (ML, in cyan).

  field and gives way to the DL phase. The electrically induced double-layer structure is thus strengthened by the [110] magnetic field to larger J and up to D m ≈ 1.66 D e . For D m > 1.66 D e , the effective Coulomb repulsion between same-sign monopoles breaks the DL phase in favour of a zincblende or monolayer (ML) structure of monopoles with saturated magnetisation along the [111] direction. Multiferroicity in rare-earth oxides -our results provide a clear signature of what to look for in experiments and remarkably, this double-layer structure has indeed been previously observed in Tb 2 Ti 2 O 7 under an external [110] magnetic field [32, 33] !

FIG. 5 .

 5 FIG. 5. Phase diagram parametrized for Tb2Ti2O7 for J = 2.7 K, Dm = 0.48 K and De = 0.32 K. The [110] magnetic field aligns the α chains along the [110] direction, destroying the AIAO order in favour for the double-layer structure. All error bars are smaller than the dots except for the red/green hatched region where simulations were difficult to equilibrate.

24 FIG. 6 . 24 FIG. 7 .

 246247 FIG.6. Convergence of the transition temperature as a function of measurement time tmax for both equilibration processes (i -green) and (ii -blue). For all figures De = 1, h = 0 and N = 432. We confirm that even if the process is very slow, simulations converge to the same transition temperature indicated by the dashed line.

log 4 e

 4 βDe/3 + e 2βJ + 3 e -2βJ/reproduce the asymptotic limits of the Coulomb spin liquid entropy, namely the Pauling (T → 0 + ) and paramagnetic (T → +∞) entropies. When compared with F DL , equations (12) and (13) provide respectively the lower and upper solid line of Fig.1of the main text for J ∈ [-2.08 : -1.62]. The estimated extent of the boundary shift due to electric interactions is consistent, within error bars, with our simulations.

FIG. 8 .

 8 FIG. 8.Monopole Double Layer (DL) (top) and Monopole Mono Layer in a [110] field (ML) (bottom). All tetrahedra carry a single magnetic charge. The green arrows are the electric polarization carried by each monopole. The thick lines are the α-chains of spins along the [110] direction. There is no global polarization but there is a saturated magnetization along the α-chains in both configurations. For a 3 in -1 out configuration (resp. 3 out -1 in), the minority spin is the outward (resp. inward) spin. In the ML phase, each minority spin (circled in orange) is a minority spin for both adjacent monopoles.


Karim Essafi
email: karim.essafi@oist.jp

Owen Benton
email: owen.benton@oist.jp

L D C Jaubert

A Kagome Map of Spin Liquids: from XXZ to Dzyaloshinskii-Moriya Ferromagnet

Motivated by a growing number of materials, the kagome lattice sits at the crossroad of modern research efforts in quantum spin liquids, chiral phases, emergent skyrmion excitations and anomalous Hall effect to name but a few. In light of this diversity, our goal in this paper is to build a unifying picture of the underlying magnetic degrees-of-freedom. We present a three-fold mapping on the kagome lattice which transforms the celebrated Heisenberg antiferromagnet and XXZ model onto two lines of time-invariant Hamiltonians. These two lines have anisotropic nearest-neighbour couplings and opposed Dzyaloshinskii-Moriya interactions. The mapping is exact for classical and quantum spins alike, i.e. it preserves the energy spectrum of the original Heisenberg and XXZ models. As a consequence, at the classical level, all phases have an extensive ground-state degeneracy. These ground states support a variety of phenomena such as ferromagnetically induced pinch points in the structure factor and the possibility for spontaneous scalar chirality. For quantum spin-1/2, the XXZ model has been recently shown to be a quantum spin liquid. Applying our three-fold mapping to the XXZ model gives rise to a connected network of quantum spin liquids, centered around a paragon of quantum disorder, namely the Ising antiferromagnet. We show that this quantum disorder spreads over an extended region of the phase diagram at linear order in spin wave theory, which overlaps with the parameter region of Herbertsmithite ZnCu3(OH)6Cl2. We conclude this work by discussing the consequences of our results in light of the chiral spin liquids found on kagome with further nearestneighbour interactions, and of the recently synthesized ternary intermetallic materials.

Competing interactions have proven able to stabilize extended phases where chirality could be encoded in the spin texture, i.e. coming from the collective behaviour of spins. This spin-chirality is responsible for phenomena as varied as the anomalous Hall effect 1-4 , multiferroicity 5 and possibly high-T c superconductivity 6 . In this context, kagome systems are fertile soil for exotic spin textures. Motivated by a growing number of materials [7][8][9][10][11][12] , the kagome lattice, whose name comes from a traditional Japanese woven bamboo pattern 13 , has attracted the attention of chemists, experimentalists and theorists alike. The classical kagome antiferromagnet is a canonical example of order-by-disorder 14 , a counter-intuitive mechanism where order is induced by fluctuations 15 . As for its quantum counterpart, it is one of the few models that has been confirmed to be a quantum spin liquid by a gamut of complementary approaches [16][17][18][19][20] . Recently, the kagome lattice has also been shown to support examples of the long-sought Kalmeyer-Laughlin chiral spin liquid [21][22][23][24] , a bosonic analogue of the fractional quantum Hall effect with anyonic excitations 25,26 .

Our present work sits at the frontier of these ideas of unconventional phenomena, spin liquids and chiral phases. We unveil a three-fold mapping between kagome spin liquids, summarized in Figs. 1 and4, which is exact both at the classical and quantum level. This mapping brings into a general framework the famous Heisenberg antiferromagnet and XXZ models, together with a continuously connected network of systems with Dzyaloshinskii-Moriya and anisotropic ferromagnetic couplings. All interactions are time-invariant and between nearest neighbours (see Hamiltonian (1)). For the end points of this connected network, time-reversal (T ) symmetry can be spontaneously broken in the classical ground state, giving rise to finite scalar chirality. The Heisenberg antiferromagnet maps onto a pair of systems characterized by ferromagnetic pinch points in their structure factor, signatures of algebraic correlations constrained by an effective local flux conservation. Interestingly for quantum spin-1/2, our work puts the Ising antiferromagnet at the centre of this connected network of quantum spin liquids, shedding a new light on the reluctance of this model to order 28 . On the experimental front, our phase diagram includes the Herbertsmithite compound ZnCu 3 (OH) 6 Cl 2 which sits at the tip of an extended region of quantum disorder at linear order in spin wave theory. Our work is also motivated by the experimental possibility to explore a broad range of anisotropic interactions in the recently synthesized rare-earth kagome materials Dy 3 Ru 4 Al 12 11 and Yb 3 Ru 4 Al 12 12 and in optical lattices 29,30 .

I. PRESENTATION

A. Model

We focus on the nearest-neighbour Hamiltonian with anisotropic XXZ and Dzyaloshinskii-Moriya interactions:

H = ij J ⊥ S ⊥ i • S ⊥ j + J z S z i S z j + D z • S i × S j . (1)
We shall first consider classical Heisenberg spins of unit length

| S i | = 1 with in-plane components S ⊥ i = (S x i , S y i ) HAF X -X + XXZ 0 FDM - FDM + XXZ X X Z - X X Z + B - i z x 0 y 0 x 1 y 1 x 2 y 2 global basis z x y B + i z x 0 y 0 x 1 y 1 x 2 y 2
FIG. 1. Three-fold mapping between kagome spin liquids -Top: We show the existence of an exact one-toone mapping, made of local proper spin rotations (bottom panels), between the celebrated Heisenberg antiferromagnet (HAF) and two novel spin liquids X ± with opposite vector chirality. By tuning the anisotropy coupling δ of Eqs. ( 5) and ( 6), our mapping directly extends onto the anisotropic XXZ model. Its chiral counterparts (named XXZ ± ) share the same extensive ground-state degeneracy as the XXZ model, until the end point δ = -1/2 (FDM ± ), which belongs to the ferromagnetic model with Dzyaloshinskii-Moriya interactions, and where chirality becomes scalar. Bottom: The local bases are rotated by ∓ 2π 3 around the z-axis when moving from

B ± 0 → B ± 1 → B ± 2 .
The z-axis are the same for all bases, which are right-oriented.

for a system of N spins. The sublattice indices and Cartesian bases are given in Fig. 1. For perfect kagome symmetry, the Dzyaloshinskii-Moriya vector is restricted along the unit vector z, orthogonal to the kagome plane 31,32 , using the clockwise convention for choosing the pairs of spins around the triangles.

B. Heisenberg antiferromagnet (HAF)

Parametrized by J ⊥ = J z = J > 0 and D = 0, the extensively degenerate ground-state manifold of the HAF is locally constrained by a magnetization flux conservation. This constraint appears clearly if the Hamiltonian is rewritten as

H HAF = J ij S i • S j = J 2 ∆ 2 i=0 S i 2 -N J, (2)
where the flux conservation takes the form of a null magnetization on all triangles ∆:

2 i=0 S i = 0.

II. THREE-FOLD MAPPING

The peculiarity of the HAF lies in the form of its Hamiltonian (2). The idea of this paper is to find a oneto-one mapping (automorphism) of the spin degrees-offreedom which gives a Hamiltonian that can be re-written in the same form, while conserving the kagome symmetry and the spin unit-length, without imposing any spurious constraints.

To ensure the spin unit length | S i | = 1, we consider local transformations Γ acting on each spin independently, i.e. transformations from the global basis to a local one, B i :

S Bi i = Γ Bi S i with | S Bi i | = 1.
Then for the transformation to be non-triviali.e. for B i to be non-uniform -and to respect translation invariance, we attach one basis B i to each kagome sublattice. As a result, there are only two transformations respecting the space group symmetry of the kagome lattice. They are made of local proper rotations as illustrated in Fig. 1. They transform the HAF into the following models which we name X ±

H X ± = J 2 ∆ 2 i=0 S B ± i i 2 -N J (3) = - J 2 i,j S ⊥ i • S ⊥ j -2S z i S z j ± √ 3 z • ( S i × S j ) (4)
where

J = J z = -2J ⊥ = ∓2D/ √ 3 > 0.
Since H HAF and H X ± have the same form (see Eqs. ( 2) and ( 3)), spin configurations connected by the one-to-one mappings Γ B ± i necessarily have the same energy in their respective Hamiltonians. Hence, the HAF, X -and X + models have the same energy spectrum and thus the same extensive ground-state degeneracy. But the spin rotation confer them very peculiar signatures when probed magnetically.

The ground state manifold of the Heisenberg antiferromagnet supports algebraic spin correlations 33 . In neutron scattering measurements, these correlations take the form of anisotropic diffuse scattering known as "pinch points" (see Ref. [34] for a review). As depicted in Fig. 2, pinch-point singularities are clearly visible in the structure factors of the X ± models. The striking similarity of the HAF and X ± structure factors is actually a quantitative illustration of their underlying equivalence. But because the planar spin components are respectively antiferromagnetically and ferromagnetically coupled in H HAF and H X ± , their collective fluctuations induce reversed spin correlations. This provides a noticeable example of pinch points induced by continuous ferromagnetic fluctuations.

(a) S ⊥ (q) • S ⊥ (-q) HAF (b) S ⊥ (q) • S ⊥ (-q) X ±
As T → 0 + , the X ± models are expected to undergo the same thermal order-by-disorder selection as the Heisenberg antiferromagnet 14 , with the additional flavour that the octupolar order 35 now bears a finite vector chirality.

III. A CONNECTED FAMILY OF SPIN LIQUIDS

Spin chirality takes multiple forms. The intrinsic noncollinearity of the spins is directly measured by the vector chirality χ ij = S i × S j . For triangular units, one can further define a scalar chirality χ ijk = S i • S j × S k which is a measure of the solid angle formed by the three spins. Vector chirality comes from the spin current involved in the strong magneto-electric coupling of some multiferroics 5 and the emergence of skyrmion excitations 37 . As for scalar chirality, it is often associated with the possibility for anomalous Hall effect when coupled to itinerant electrons 1-4 .

While vector-chirality is intrinsically induced by the Dzyaloshinskii-Moriya term, we do not expect any longrange scalar-chiral order in the X ± models since the HAF spin liquid does not break T symmetry. It is thus tantalizing to see if, by taking advantage of the present threefold mapping, it were possible to tune the Hamiltonians and induce scalar chirality spontaneously. Since our three-fold mapping does not affect the z-axis (see Fig. 1), decreasing J z has the same influence on the HAF and X ± Hamiltonians: it tunes the HAF into the XXZ model

H XXZ = J i,j S ⊥ i • S ⊥ j + δ S z i S z j (5)
which is mapped onto what we name the XXZ ± models

H XXZ ± = - J 2 i,j S ⊥ i • S ⊥ j -2δ S z i S z j ± √ 3 z • ( S i × S j ) (6)
with J > 0. This mapping is valid for all δ but for -1/2 < δ < 1, the XXZ ground state remains a sub-ensemble of the HAF one (δ = 1) where all spins lie in plane (χ ijk = 0). This model is equivalent to the three-coloring problem up to a global O(2) symmetry 38 , whose degeneracy is countable and extensive 39 . As illustrated in Fig. 3, the system is entirely paved with only two kinds of triangular configurations, A and Ā, with opposite vector-chirality 38 . The noticeable consequence of our mapping is that in the XXZ -and XXZ + ground states, the A and Ā configurations are respectively replaced by a collinear state F with zero chirality; the resulting imbalance ensures finite and opposite vector-chirality between the two XXZ ± ground states, while preserving their extensive degeneracy.

From this point of view, it is interesting to think of the XXZ ground state as coming from the cancelation of positive and negative Dzyaloshinskii-Moriya terms, once ferromagnetism has been taken out.

IV. CHIRAL SPIN LIQUIDS

On the other hand, for δ = -1/2 Dzyaloshinskii-Moriya interactions become perfectly balanced by isotropic ferromagnetic coupling

H FDM ± = - J 2 i,j S i • S j ± √ 3 z • ( S i × S j ) . (7)
As a consequence, for each triangle, both the Dzyaloshinskii-Moriya (DM) induced 32 and ferromagnetic (FM) ground-state configurations minimize the classical energy

S ={0,1,2} = (sin θ cos φ ± , sin θ sin φ ± , cos θ), (8)
DM : φ ± = φ ± 2π 3 ⇒ χ 012 = ± 3 √ 3 2 cos θ sin θ 2 F M : φ ± = φ ⇒ χ 012 = 0
where the ± index distinguishes the two FDM ± models. With respect to the XXZ ± models where θ was imposed to be zero, the global degeneracy of the FDM ± ground states is enhanced to O(3). Thus, while the S ⊥ degrees-of-freedom conserves the character of a classical spin liquid, with the extensive degeneracy and algebraic correlations of the three-coloring problem, S z = cos θ can now take a finite uniform value, conferring a finite scalar chirality to any ground-state configuration with θ = {0, π/2, π} (see Eq. ( 8)).

The three-fold mapping transforms the FDM ± models back into the XXZ 0 Hamiltonian of Eq. ( 5) with δ = -1/2, where the scalar chirality has vanished but the enhanced global O(3) degeneracy remains. It is noteworthy that the end point value δ = -1/2 takes an elegant meaning along the XXZ ± lines, namely that the ferromagnetic coupling becomes isotropic, which is hidden if only considering the XXZ model.

The emergence of scalar chirality in what is essentially a "simple" ferromagnet with Dzyaloshinskii-Moriya interactions is quite remarkable, with a rich potential for unconventional phenomena. For example the interplay between a chiral spin liquid and itinerant electrons is an up-and-coming topic 22,40 . Indeed, the FDM ± ground -c) In this context, the apparition of vector-chirality in the XXZ ± ground states is understood as the suppression of either the A or Ā configurations in favor of a collinear state (F). The same scenario holds for the XXZ0 and FDM ± ground states where the finite out-of-plane magnetization makes the chirality scalar. An example of spin configuration with finite scalar chirality is given in panel (d): the planar projection of the spins corresponds to configuration A.

state is neither fully ordered like a solid, or paramagnetic like a gas. In a pictorial way it is a magnetic liquid where strong correlations and fluctuations co-exist, which can then couple via double-exchange to another "fluid" made of itinerant electrons. While hopping on the scalar-chiral spin texture, the itinerant electrons pick up a Berry phase that might not only induce anomalous Hall conductivity 1,22,41,42 , but at the same time feed back into the strongly correlated spin texture to induce exotic magnetic order 22,[43][44][45] . This feedback actually does not require scalar chirality and would also be pertinent to the XXZ ± lines.

It should be noted that given the large value of D = √ 3J, an experimental realization of the FDM ± models per se would arguably be difficult in solid state physics, but on the other hand, particularly motivating in op-tical lattices. Indeed, the kagome geometry 29 and spin anisotropy 30 have been experimentally realized with ultracold atoms. There is also good hope that the active research on synthetic gauge fields might be able to produce synthetic Dzyaloshinskii-Moriya interactions [46][47][48] , with the caveat that the Dzyaloshinskii-Moriya vector should be out-of-plane here.

Last but not least, ferromagnetic insulators with Dzyaloshinskii-Moriya coupling have been studied in the context of magnon Hall effect, i.e. where a transverse heat current is induced by a temperature gradient. It is intriguing to notice that the FDM ± sits at the frontier between two different topological phases, indicating the closing of a gap between two magnon bands 49 . In light of our present work, and since the topological phase for D < √ 3J is the same down to D = 0 49 , it would be of great interest for future work to study the finite temperature physics of the Dzyaloshinskii-Moriya ferromagnet. This is especially true since chiral magnonic edge states and topological skyrmion excitations have been observed in simulations for D/J ∼ 0.4 50 .

V. QUANTUM FLUCTUATIONS

Our analysis has been so far focused on classical spins in order to precisely determine the nature of their classical ground states. However, it is important to keep in mind that our present three-fold mapping is also exact with quantum spins, since the local transformations are proper rotations, i.e. unitary matrices, and therefore preserve the commutation relations of the spin components.

We have investigated the consequences of quantum fluctuations for these spin liquids in the framework of linear spin wave theory within the parameters of Hamiltonian (1). These results are illustrated in Fig. 4. Approaching any of the HAF or X ± models (marked by dots), the linear spin wave Hamiltonian takes the same form, which simply confirms the equivalence of these three spin liquids in the presence of quantum fluctuations. The linear spin wave calculation also indicates the likelihood of quantum disorder around the center of Fig. 4. Approaching the white-triangle region, a flat band of excitations collapses to zero energy leading to a divergence in the quantum correction to the order parameter.

Within the phase diagram of Fig. 4, the J ⊥ = J z = J > 0 line has drawn substantial interest for its relevance to Herbertsmithite ZnCu 3 (OH) 6 Cl 2 where Dzyaloshinskii-Moriya interactions are not negligible (D/J ∼ 0.044 -0.08) 51,[53][54][55] . We reproduce the results of Refs. [32,56] done on this line of parameters, namely that classically 32 and up to linear order in spin wave theory 56 , magnetic order is stabilized for any finite value of D. However, higher order terms in quantum fluctuations studied by Exact 1) obtained from linear spin wave theory for Jz > 0 -The HAF and X ± spin liquids are marked by dots and the yellow, orange and purple regions are long-range ordered phases. The white triangle delimits a regime where quantum corrections to the order parameters diverge, indicating a possible extended region of quantum disorder. In particular, DMRG has shown that the entire XXZ model (solid and dashed red lines) is a quantum spin liquid. Our three-fold transformation maps this quantum spin liquid onto the XXZ ± models of Hamiltonian (6) (green and blue lines) for 0 < δ < 1 (solid) and δ > 1 (dashed), which are thus also quantum spin liquids. Experimentally, independent parametrizations of the Herbertsmithite compound 51,52 put it at the tip of the white-triangle region (black rectangle).

Diagonalization 57,58 , Schwinger-boson 59,60 and perturbative methods 61 have shown that quantum disorder actually persists over a finite region up to D/J ∼ 0.1, which includes Herbertsmithite. Our goal here is not to claim explaining the spin liquid nature of Herbertsmithite which has been extensively studied, but rather to set our theory on an experimental footing. In particular, it should be noted that at linear order in quantum fluctuations, the small XXZ anisotropy observed in Herbertsmithite 52,54 (J ⊥ /J z ≈ 0.9) brings this compound at the tip of the white-triangle region with quantum disorder.

Over the past year, the XXZ line (D = 0 and 0 < J z < J ⊥ ) has also received significant attention for the spin-liquid nature of its ground state for spins S = 1/2 [62][63][64][65] , and the complex quantum order-bydisorder mechanism that takes place for S > 1/2 62,66 . Noticeably for spin-1/2, the density matrix renormalization group (DMRG) approach indicates that the quantum spin liquid persists for 0 < J ⊥ /J z < 1 and D = 0 63 . The consequences on our work are multiple. First of all, our three-fold mapping makes the XXZ ± lines, as well as their continuations along J ⊥ = ±D/ √ 3 < 0, quantum spin liquids for spin-1/2 (see green and blue dashed lines in Fig. 4). Furthermore, the central point of our phase diagram D = J ⊥ = 0 has Ising anisotropy and is known for remaining a quantum paramagnet even for arbitrarily small transverse fields 28,[67][68][69] . As such, it has been described as a rare example of "disorder-by-disorder", a mechanism proposed by Fazekas and Anderson 70 where quantum fluctuations select a disordered sub-manifold of the classically degenerate ground state. Within the framework of our three-fold mapping, this remarkable resistance to order can be understood as the consequence of being at the intersection of three (dashed) lines of spin liquids, in a way reminiscent of what has been observed in pyrochlore systems 71 . However, please note that the spin-1/2 phase diagram is expected to be highly anisotropic around this central point, since it has been shown to order into a superfluid phase for D = 0 and J ⊥ < 0 63,72 , and thus also along the symmetric lines J ⊥ = ±D/ √ 3 > 0 around the origin according to our three-fold mapping.

VI. CONCLUSION

We have discovered a connected network of quantum spin liquids on the kagome lattice, which can be mapped onto each other via a three-fold transformation. One of the branches of this network is the anisotropic XXZ model, known to be a quantum spin liquid for spin-1/2 [62][63][64][65] (see the red lines in Figs. 1 and4), which includes the famous Heisenberg antiferromagnet. While every triad of Hamiltonians connected by this mapping have exactly the same energy spectrum at the classical and quantum level, their spin configurations are transformed. As a consequence, the three-fold mapping of the XXZ line gives rise to new spin liquids with intrinsic vector chirality because of Dzyaloshinskii-Moriya interactions. The Ising antiferromagnet sits at the centre of this map (see Fig. 4), which sheds a new light on the unique propensity of this model to remain disordered 28,[67][68][69] .

Beyond these three branches of quantum spin liquids, we have studied the stability of Hamiltonian (1) for J z > 0, up to linear order in spin wave theory. We have found an extended region of the phase diagram in Fig. 4 where quantum disorder prevails. The small XXZ anisotropy observed in Herbertsmithite 52,54 (J ⊥ /J z ≈ 0.9) brings this compound within the tip of this extended region.

At the classical level, the Heisenberg antiferromagnet maps onto two models where algebraic correlations take the form of ferromagnetic pinch points visible in the structure factor of Fig. 2. Keeping the Dzyaloshinskii-Moriya term constant, if one tunes the J z coupling of these models until they become isotropically ferromag-netic, the chirality can spontaneously become scalar.

Our work opens a wide range of exciting directions to follow, both by theorists and experimentalists. In light of the intense research on the Heisenberg antiferromagnet and XXZ models, here we propose two lines of systems with the same energy spectra, but different (chiral) magnetic signatures. With this new probe at hand and Figs. 1 and4 in mind, it would be of great interest to look for new insights as one approaches these models and their chiral counterparts from different angles in parameter space (J ⊥ , J z , D). In particular, the spreading of quantum disorder within the white triangle of Fig. 4 and in its vicinity shall conserve the three-fold symmetry, and be mediated by quantum order-by-disorder mechanisms as we vary the spin length S 62,66 .

The inclusion of 2 nd and 3 rd nearest-neighbour interactions J 2 = J 3 = J N N N is known to stabilize a chiral spin liquid at finite value J c N N N [22][23][24]27,63,65 . This value J c N N N has been shown by DMRG to decrease as the antiferromagnetic J z coupling vanishes. This means that the chiral spin liquid is getting closer to the nearest-neighbour XXZ model as J z goes from 1 to 0. It would thus be very tempting to extend this work to ferromagnetic coupling (J z < 0) towards the XXZ 0 and equivalent FDM ± points. Since the classical ground states of the FDM ± models can support spontaneous scalar chirality, the possible connection with the chiral spin liquids at finite J N N N is a particularly attractive open question.

Beyond kagome physics, our results suggests that systems supporting the "disorder-by-disorder" mechanism 70 , such as the Ising antiferromagnet here 28 , are good places to look for hidden spin liquids in the neighbouring parameter space.

On the experimental front, our work fits within the on-going effort for the experimental realization of frustrated systems in optical lattices 29,30 , and especially to produce tunable synthetic Dzyaloshinskii-Moriya interactions [46][47][48] .

We hope that our results will motivate experimental efforts for the synthesis and characterization of kagome materials with anisotropic nearest-neighbour interactions. The recently synthesized ternary intermetallic compounds Dy 3 Ru 4 Al 12 11 and Yb 3 Ru 4 Al 12 12 are very promising materials to start with, since the 4f orbitals of rare-earth ions are known to induce very anisotropic and short-range interactions. Furthermore the presence of itinerant electrons make them natural materials to probe the chirality of the underlying spin texture. Their crystal structure, however, corresponds to a distorted kagome lattice. To impose kagome symmetry is a chemistry challenge, but so was the case for Volborthite Cu 3 V 2 O 7 (OH) 2 • 2H 2 O, 14 years ago 7 , which nucleated the synthesis of a growing number of materials with essentially perfect kagome symmetry [8][9][10] . According to our three-fold mapping, the places to look for would be large antifer-romagnetic J z , as well as around the XXZ 0 point where only small Dzyaloshinskii-Moriya terms are required. In light of Refs. [50,49], the region neighbouring the FDM ± models is also very promising, even for smaller values of D and anisotropic J z . At the proximity of these highsymmetry points, especially the one at the centre of the white triangle, chemical, hydrostatic and uni-axial pressure might help the exploration of the phase diagram, as observed in rare-earth pyrochlore oxides 73 .

FIG. 2 .

 2 FIG. 2. Structure factor of the HAF (a) and X ± (b) spin liquids -The Fourier transforms of the spin correlations have been computed using the method developed by Henley 36 for Coulomb spin liquids in which the local constraints are enforced by a projection operator in reciprocal space. We have considered the planar spin components S ⊥ where "pinch point" structures are formed in the centre of the Brillouin zones, characteristic of the local flux conservation. The structure factors clearly illustrate the underlying equivalence of the models, and the difference of their in-plane fluctuations; antiferromagnetic in the HAF and ferromagnetic in the X ± models. Only one figure is shown for the X ± phases because they cannot be distinguished by the structure factor of the planar spin components. The colour scales are fixed by the maximum of intensity on each figure.

FIG. 3 .

 3 FIG. 3.Three-coloring ground states -The classical ground-state ensembles of the XXZ model and its chiral counterparts XXZ ± are equivalent to the three-coloring problem, up to a global O(2) symmetry. (a) The equivalence is transparent for the XXZ ground state38 , where each triangle possesses the three possible spin orientations rotated by 2π/3 from each other, and S z = 0. The color code of the spin orientations is given in the top-left legend. The two antiferromagnetic permutations A= {•, •, •} and Ā = {•, •, •} are possible, giving a zero vector-chirality on average. (b -c) In this context, the apparition of vector-chirality in the XXZ ± ground states is understood as the suppression of either the A or Ā configurations in favor of a collinear state (F). The same scenario holds for the XXZ0 and FDM ± ground states where the finite out-of-plane magnetization makes the chirality scalar. An example of spin configuration with finite scalar chirality is given in panel (d): the planar projection of the spins corresponds to configuration A.

FIG. 4 .

 4 FIG. 4. Phase diagram of Hamiltonian(1) obtained from linear spin wave theory for Jz > 0 -The HAF and X ± spin liquids are marked by dots and the yellow, orange and purple regions are long-range ordered phases. The white triangle delimits a regime where quantum corrections to the order parameters diverge, indicating a possible extended region of quantum disorder. In particular, DMRG has shown that the entire XXZ model (solid and dashed red lines) is a quantum spin liquid. Our three-fold transformation maps this quantum spin liquid onto the XXZ ± models of Hamiltonian (6) (green and blue lines) for 0 < δ < 1 (solid) and δ > 1 (dashed), which are thus also quantum spin liquids. Experimentally, independent parametrizations of the Herbertsmithite compound51,52 put it at the tip of the white-triangle region (black rectangle).

ACKNOWLEDGMENTS

It is a pleasure to thank Nic Shannon, Mathieu Taillefumier, Andreas Laeuchli and Andriy Nevidomskyy for useful discussions. This work was supported by the Okinawa Institute of Science and Technology Graduate University.