Competing spin liquids and hidden nematic order in spin ice with frustrated transverse exchange

Frustration in magnetic interactions can give rise to disordered ground states with subtle and beautiful properties. The spin ices Ho2Ti2O7 and Dy2Ti2O7 exemplify this phenomenon, displaying a classical spin liquid state, with fractionalized magnetic--monopole excitations. Recently there has been great interest in closely-related"quantum spin ice"materials, following the realization that anisotropic exchange interactions could convert spin ice into a massively-entangled, quantum, spin liquid, where magnetic monopoles become the charges of an emergent quantum electrodynamics. Here we show that even the simplest model of a quantum spin ice, the XXZ model on the pyrochlore lattice, can realise a still-richer scenario. Using a combination of classical Monte Carlo simulation, semi--classical molecular--dynamics simulation, and analytic field theory, we explore the properties of this model for frustrated transverse exchange. We find not one, but three, competing forms of spin liquid, as well as a phase with hidden, spin-nematic, order. We explore the experimental signatures of each of these different states, making explicit predictions for inelastic neutron scattering. These results show an intriguing similarity to experiments on a range of pyrochlore oxides.


I. INTRODUCTION
The search for spin liquids -disordered phases of magnets which support entirely new forms of magnetic excitation -has become one of the defining themes of modern condensed-matter physics [1][2][3]. In this context, the pyrochlore lattice, a corner-sharing network of tetrahedra found in a wide range of naturally-occurring minerals, has proved an amazing gift to science. Pyrochlore magnets play host to a variety of unconventional forms of magnetic order, and include systems which have not been observed to order at any temperature [4]. Perhaps the most celebrated of these is the "spin ice" found in the Ising magnets Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 [5]; a classical spin liquid, described by an emergent U(1) lattice gauge theory with magnetic monopole excitations [6].
The main message of this Article is that even the simplest model of a quantum spin ice -the XXZ model on a pyrochlore lattice -has far more to offer than spin ice alone. Working in the classical limit, accessible to large-scale simulation, we find that frustrated transverse exchange gives rise to not one, but three, distinct, spin-liquid regimes [ Fig. 1]. We explore the way in which these spin liquids compete, identify the different gauge groups associated with each spin liquid, and make explicit predictions for their experimental signatures [Fig. 2]. We find that one of these spin liquids posses a highly-unusual U(1)×U(1) gauge structure and, as an added bonus, undergoes a phase transition into a state with hidden, spinnematic, order. We also use molecular dynamics simulations to characterise the excitations of this spin-nematic phase [Fig. 5]. The portrait which emerges has striking similarities with the behavior of a number of pyrochlore materials.
The simplest model able to capture quantum effects in a spin ice [7][8][9][10][11][12] is the XXZ model on the pyrochlore lattice where S i = (S x i , S y i , S z i ) is a (pseudo)spin-half operator describing the two states of the lowest energy doublet of a magnetic ion subject to a strong crystal electric field (CEF). The symmetry of the lattice requires that the quantization axis of each spin (here, S z i ) lies on a local [111] axis, as defined in Appendix A.
Ising interactions, J zz > 0, favor states obeying the "ice rules" in which two spins point into, and two spins point out of, each tetrahedron on the lattice. The transverse term, J ± , introduces dynamics about these spin-ice configurations and, for larger, positive, values of J ± /J zz , drives the system into a state with easy-plane order [17,18,[47][48][49][50]. The physical meaning of this easy-plane order depends on the nature of the magnetic ion. For Kramers ions like Yb 3+ and Er 3+ all components of S relate to a magnetic dipole moment [11], and the ordered phase is an easy-plane antiferromagnet. However for non-Kramers ions such as Pr 3+ and Tb 3+ [9,44], or "dipolar-octupolar" Kramers ions like Nd 3+ or Ce 3+ [51], the easy-plane order may have quadrupolar (octupolar) character. In what follows, we consider explicitly the case of Kramers ions. However, suitably reinterpreted, these results also have important implications for non-Kramers ions.
For J ± > 0, H XXZ [Eq. (1)] is unfrustrated, in the sense that it is free of sign problems in Quantum Monte Carlo (QMC) simulation. In this case, the phase diagram is already well-established [14,21,23]. For J ± /J zz 0.05, QMC simulations find a crossover from a conventional paramagnet into a classical spin-liquid (spin ice) at a temperature T * /J zz ∼ 0.2, and a second crossover into a quantum spin liquid (QSL) at a much lower temperature T * QSL /J zz ∼ (J ± /J zz ) 3 . In the low temperature quantum spin liquid regime, the magnetic monopoles of classical spin ice become dynamic, fractional, spin excitations (spinons), while the spectrum of the model also includes gapless photons [13,16]. For J ± /J zz 0.05, the U(1) QSL gives way to easy-plane antiferromagnetic order (AF ⊥ ), in which spins lie in the plane perpendicular to the local S z -axis [14,21,23].
Very little is known about the properties of H XXZ for J ± < 0 [10,18,44]. On perturbative grounds, it is expected that the ground state for |J ± |/J zz 1 will also be a U(1) QSL [13], albeit one with a modified spinon dispersion [18,52]. Gauge mean-field calculations suggest that this QSL persists over a broad range of parameters, −4.13 J ± /J zz < 0 [18]. But the nature of competing ordered -or disordered -phases for J ± < 0 remains an open question.
There are many reasons to believe that the proper-ties of the quantum spin ice model, H XXZ [Eq. (1)] for frustrated coupling J ± < 0, could be even richer than for J ± > 0. In particular, for J ± /J zz = − 1 2 , H XXZ [Eq. (1)] is equivalent (up to a site-dependent spin-rotation), to the Heisenberg antiferromagnet (HAF) on a pyrochlore lattice. Like spin ice, the HAF is known to support a classical spin liquid [53][54][55][56], and it has also been argued to support a QSL ground state [57][58][59][60]. And, crucially, both the classical and quantum spin liquids in the HAF have a qualitatively different character from those found in spin ice. This sets up a competition between two different kinds of spin liquid, namely spin ice for J ± /J zz ≈ 0, and a state homologous to the HAF for J ± /J zz ≈ − 1 2 . It also opens the door for yet more novel magnetic phases for J ± /J zz < − 1 2 .

II. PHASE DIAGRAM DETERMINED BY CLASSICAL MONTE CARLO SIMULATIONS
Since the quantum spin ice model, H XXZ [Eq. (1)], is inaccessible to QMC for J ± < 0, we instead study its finite-temperature properties using classical Monte Carlo (MC) simulation -the results are summarised in the phase diagram Fig. 1. For J ± > 0, this phase diagram is very similar to that previously found in QMC simulations [14,21,23] -at a qualitative level, the only significant difference is the absence of a QSL below T * QSL /J zz ∼ (J ± /J zz ) 3 0.005. At a quantitative level, we find changes in numerical values of the crossover temperature associated with the spin ice regime, T * 1 , and the position of the zero-temperature boundary between SI and AF ⊥ . These changes can be ascribed to the fact that the magnetic monopoles (spinons) are not quantized in classical simulations and do not develop phase coherence [62]. Further details of classical MC simulations for J ± > 0 will be presented elsewhere [63].
We now turn to the frustrated case, J ± < 0. At low temperatures, spin-ice correlations persist up to J ± /J zz = − 1 2 [10,44], as illustrated in Fig. 2a. Upon reaching J ± /J zz = − 1 2 the system becomes thermodynamically equivalent to a HAF. This high-symmetry point gives rise to a new form of spin liquid at finite temperature, labelled a pseudo-Heisenberg antiferromagnet (pHAF) in Fig. 1. Once again, this spin liquid has algebraic correlations, as shown in Fig. 2c, but with qualitatively different character from spin ice [ Fig. 2a]. These correlations persist up to a crossover temperature T * 3 associated with the Curie-law crossover (CLC) in the magnetic susceptibility [64].
While the correlations measured in the equal-time structure factor S(q) are also different from those found in the HAF [54,55,65], the two models are equivalent up to a local coordinate transformation. And, by analogy with earlier work on the HAF [55,56,66], the spin liquid pHAF can be described by a U(1)×U(1)×U(1) gauge theory.
The situation for J ± /J zz < − 1 2 is even more interest-ing. Below a second crossover scale, T * 2 < T * 3 , identifiable by a reduction in the fluctuations of the z-components of the spins [see Appendix B], the pHAF gives way to an easy-plane spin liquid, labelled SL ⊥ in Fig. 1. Spin correlations in this regime have algebraic character, with pinch-points in S(q) [Fig. 2b]. However these correlations are qualitatively different from those in either spin ice [ Fig. 2a], or the pHAF [ Fig. 2c]. At a still lower temperature, T SN < T * 2 , the system undergoes a thermodynamic phase transition, marked by a clear anomaly in the specific heat. None the less, this phase transition does not give rise to any magnetic Bragg peaks in S(q) and, at least as far as dipolar spin correlations are concerned, the system remains disordered.
While the new phase for T < T SN -labelled SN ⊥ in Fig. 1 -does not exhibit any conventional magnetic order, it does posses a hidden, spin-nematic, order. The ordered state does not break translational symmetry, but breaks the U(1) symmetry of H XXZ [Eq. (1)] by selecting an axis within the local xy-plane. Such an order can be described by the bond-based order parameter [67][68][69] where the sum on ij runs over the nearest-neighbour bonds of the lattice, and This type of easy-plane order is formally identical to the spin-nematic phases found in a range of frustrated magnets in applied magnetic field [69][70][71]. And in common with these systems, the associated Landau theory lacks a cubic term, and therefore permits a continuous phase transition. Simulations suggest that the phase transition at T = T SN is indeed continuous for J ± /J zz − 1 2 , becoming first-order approaching the high-symmetry point J ± /J zz → − 1 2 . Further details of the thermodynamics of this transition are given in Appendix B.

III. THEORY OF THE EASY-PLANE SPIN LIQUID
Spin correlations in spin ice (SI) can be described using a U(1) lattice gauge theory [6,56,72], which gives rise to characteristic "pinch-points" in the spin structure factor S(q) [Fig. 2a]. Meanwhile, for classical spins, spin correlations in the Heisenberg AF on the pyrochlore latticeand by extension in the pHAF -can be described using a U(1)×U(1)×U(1) gauge theory [53][54][55][56]. The pHAF has qualitatively different pinch-points from spin ice, as illustrated in Fig. 2c. It is clear that the correlations of the easy-plane spin liquid, SL ⊥ [ Fig. 2b] are very different from either spin ice [ Fig. 2a] or the pHAF [ Fig. 2c]. None the less, the presence of pinch points suggests that  SL ⊥ , too, may be described by some form of gauge theory.
We can develop a field-theory for the spin-liquid SL ⊥ by applying the methods developed in [50,73]. The starting point of this approach is to recast the spins S i in H XXZ [Eq. (1)] in terms of five order-parameter fields defined on each tetrahedron r. These objects m λ (r) describe the different kinds of four-sublattice magnetic order consistent with the point group symmetry of the pyrochlore lattice. Definitions of each field m λ in terms of the spins S i are given in Appendix D.
The most general exchange Hamiltonian on the pyrochlore lattice can be transcribed exactly in terms of m λ [49]. This greatly simplifies the determination of classical ground states and, where classical ground states form an extensive manifold, one can use this approach to determine the local constraints which control the resulting spin-liquid [50,73]. In the case of SL ⊥ , for T → 0, we have The spin fluctuations at low temperature are thus dominated by the fluctuations of the remaining fields m T2 (r) and m T1planar (r). These fields have significance as the order-parameters of the competing four-sublattice magnetic orders which would be induced by the symmetryallowed perturbation where γ ij are complex phase factors arising from the change in coordinate frame between different lattice sites [7, 9-11, 47, 74]. For this reason, the spin-liquid SL ⊥ falls very naturally into the "multiple-phase competition" scenario for pyrochlore magnets [49,50,75,76]. In Fig. 3, we show the classical ground-state phase diagram of anisotropic exchange model This contains three distinct regions of 4-sublattice order : the easy-plane ordered phases described by the fields m E (denoted AF ⊥ in Fig. 1), m T1planar , and m T2 (Palmer-Chalker state [77]). These border a region of spin ice (denoted SI in Fig. 1), dominated by fluctuations of m T1ice . We note that a closely-related phase diagram has been derived for non-Kramers ions [10,44]; in this case easy-plane order must be interpreted in terms of the quadrupole moment of the magnetic ion. The non-trivial correlations in the spin-liquid SL ⊥ arise from the fact that neighbouring tetrahedra share a spin, so that the fields m λ (r) on neighbouring tetrahedra are not independent of one another. This point, combined with Eq. (5), imposes spatial constraints on the fluctuations of m T2 (r) and m T1planar (r). After coarse graining to extract the long wavelength physics these constraints may be written in terms of two, independent, vector fluxes  4)], as described in [49]. The minimal model of a quantum spin ice H XXZ [Eq. (1)] exists on the line J±± = 0 -for J± < − 1 2 (white line), two phases with 4-sublattice easy-plane order meet, and the resulting enlarged ground-state manifold gives rise to the easy plane spin liquid SL ⊥ , and spin-nematic phase SN ⊥ . A closely-related mean-field phase diagram for non-Kramers ions is given in [10,44].
which each separately obey their own Gauss' law We can therefore write and the theory has two, independent, U(1) gauge degrees of freedom. The free energy associated with the fluctuations of these fields is of entropic origin [66]. The only choice of Gaussian free-energy consistent with both the point group symmetry and the U (1) symmetry of H XXZ is where the coefficient λ can be determined through fits to simulation, or a large-N expansion [55,73]. It follows from the existence of the conserved fluxes B 1 and B 2 and the free-energy Eq. (11) that SL ⊥ is a Coulomb phase with algebraic correlations [66]. The validity of this description is demonstrated in Fig. 4 where we compare analytic calculations of the flux structure factor   At finite temperature, we anticipate that the spin liquid SL ⊥ will be perturbatively stable against terms such as δH ±± [Eq. (6)], which retain the point-group symmetry of the lattice but lift the U (1) symmetry of the spins. In this case the free energy will be modified : This form of free energy will still lead to pinch points in S αβ Bµ (q) and S(q), but these will take on a more anisotropic character.

IV. DYNAMICS IN THE SPIN-NEMATIC PHASE
For temperatures, T < T SN the easy-plane spin-liquid (SL ⊥ ) gives way to a phase with hidden spin-nematic order, labelled SN ⊥ in Fig. 1. As far as the dipole moments (a) Dynamical structure factor for spins, S(q, ω) To better understand the dynamics of the spinnematic phase, we have calculated the dynamical structure factor S(q, ω), within a semi-classical moleculardynamics (MD) simulation, using the methods described in [78]. Relevant definitions are given in Appendix C.
For ω/J zz 0.2, S(q, ω) presents a featureless, nondispersing continuum [Fig. 5a]. Relics of dispersing excitations are visible in S(q, ω) at higher energies, but these are explicitly not Goldstone modes, and have nothing to do with the hidden spin-nematic order. Examining the evolution of S(q, ω) as a function of temperature, we find that results for S(q, ω) in the spin-nematic phase for T < T SN , are very similar to those found in the spin liquid SL ⊥ for T > T SN .
Incoherent, non-dispersing structure of the type shown in Fig. 5a is reminiscent of theoretical predictions [79][80][81] and experimental measurements [44,82,83], for a wide range of different spin liquids. In a quantum spin liquid the presence of a non-dispersing continuum reflects the fact that, unlike conventional spin waves (magnons), single elementary excitations of a spin liquid cannot be created by local processes. It follows that, when a neutron scatters from a spin liquid, the energy, momentum and angular momentum (spin) transferred are not absorbed by a single excitation with a well-defined energy and momentum, but rather shared between multiple excitations [84]. In the semi-classical limit studied here, it is probably unsafe to attribute such a continuum to fractionalized excitations [63]. However, the fact that S(q, ω) only records dipolar spin correlations obscures a more important fact -the spin-nematic order which is present for T < T SN which breaks a continuous, U(1), symmetry of the Hamiltonian. And, by Goldstone's theorem, it must, therefore, also support gapless Goldstone modes.
In order to resolve this conundrum, it is necessary to examine the dynamical correlations of the quadrupole moments of spin. In Fig. 5b we present MD simulation results for the dynamical susceptibility χ Q site ⊥ (q, ω), which measures fluctuations of the on-site quadrupolar moments, which are well defined for classical spins [cf.
Appendix B]. A sharp excitation, with dispersion can now be resolved near to the zone centers with q rl = (0, 0, 0), (1, 1, 1), (2,2,2). These are the same zone centers for which the Bragg peaks associated with the hidden spin-nematic order SN ⊥ would occur in a quadrupolar structure factor, which might, in principle, be measured in resonant X-ray experiments [85]. At present, relatively little is known about the dynamical properties of spin-nematic states. Field-theoretic analysis [71,[86][87][88], based on the symmetry of the order parameter, predicts that spin-nematics support gapless Goldstone modes, visible in χ Q ⊥ (q, ω). This Goldstone mode has dispersion ω ∝ |q| [cf. Eq. (15)], and at zero temperature the associated intensity diverges as ∼ 1/ω for ω → 0 [71]. The same behaviour is seen in "flavour-wave" calculations and QMC simulations of spin-1 models constructed to support quadrupolar order [89][90][91]. The dynamics of the spin-1/2 frustrated ferromagnetic spin chain have also been studied using DMRG, and reveal a broad continuum of excitations at high energies [92]. However, because of the absence of long-range order, no Goldstone modes can be resolved. A continuum of excitations at high energies is also found in calculations for two-dimensional frustrated ferromagnets, within a slave-particle mean-field picture [93].
Our MD simulations of SN ⊥ clearly reveal a linearlydispersing Goldstone mode, with intensity which diverges for ω → 0 [cf. inset to Fig. 5b]. The form of this divergence is ∼ 1/ω 2 , rather than ∼ 1/ω. This follows from the fact that simulations are carried out at finite temperature, and probe thermal rather than quantum fluctuations. Most striking, however, is the broad continuum of excitations visible in both spin-and quadrupole structure factors. And it is also interesting to note that the way in which the Goldstone mode "dissolves" into this continuum bears some resemblence to what is seen in QMC simulations of a spin-1 model, at higher temperature [91]. Overall, the picture which emerges from MD simulation is consistent with all known facts about spin-nematics, and should provide a reliable guide for comparison with experiment.
Further details of the spin dynamics in the spinnematic phase, and specifically the characterization of the Goldstone mode are given in Appendix F. We note that the U(1) symmetry of H XXZ [Eq. (1)] is not a necessary condition for spin-nematic order to exist. However, if this symmetry were broken, the low-energy (pseudo-)Goldstone mode associated with SN ⊥ would acquire a small gap.

V. DISCUSSION
Spin liquids [1][2][3] and spin nematics [67][68][69] are prime examples of unconventional states of matter, and have many unusual and interesting properties. The experimental search for these exotic states has a long history, with many twists and turns, and not a few dead ends. Given this, finding both in one simple, canonical, and experimentally-motivated model is remarkable. It is therefore worth considering the possibilities for observing the unconventional spin liquid SL ⊥ , and the spinnematic SN ⊥ , in the type of rare-earth pyrochlore magnet described by Eq. (1).
In the case of SN ⊥ , it is important to make a distinction between the type of spin-nematic order considered in this manuscript, which is driven by fluctuations, and the quadrupolar or octupolar order which can arise directly from the ground states of rare-earth ions. Here we particularly have in mind the non-Kramers ions Pr 3+ [9,18,44] and Tb 3+ [9,39,45,94], and Kramers doublets of dipolar-octupolar character, such as Nd 3+ [51], and Ce 3+ [95]. For these ions, quadrupolar or octupolar order may occur as a"classical" order of the transverse part of the pseudospins S i . The multipolar character of the order follows from the symmetry of the crystalfield ground-state (doublet) of the magnetic ion, which is described by S i . Where multipolar order of this kind occurs, experiments which probe the dynamics of dipoles will see a gapped response and a sharp excitation spectrum. In contrast, in the easy-plane spin-nematic SN ⊥ , dipole moments remain in a spin-liquid like state, with strong fluctuations at low temperature and a broad, gapless response coexisting with the hidden nematic order [ Fig. 5a].
Where, then, might we observe these unusual magnetic states? Further experimental work will be necessary to definitively answer this question, but there are already a few trails to follow. In particular, the Prbased pyrochlores have the recommended single-ion and interaction anisotropies [9,18,44]. Coupling parameters of Pr 2 Zr 2 O 7 for example have been suggested to sit in the AF ⊥ phase of Fig. 1 [44], although it seems that the coupling of structural disorder to the non-Kramers doublets plays a significant role [24,34,96,97]. Since chemical pressure has already proven to be a useful tool to move a family of compounds across a phase diagram [42,49,76,98,99], Pr 2 X 2 O 7 (X=Sn,Hf,Pb) are promising candidates to investigate, with ferromagnetic correlations consistent with positive J zz and no dipole order yet observed [27,32,33,41,100].
The notion of hidden order also resonates with the elusive physics of Yb-based pyrochlores. As far as we know, Yb pyrochlores lie in a different regime of magnetic interactions than the H XXZ model of Eq. (1). Specifically, experiments on Yb 2 Ti 2 O 7 point to an unfrustrated value of J ± > 0 [11,101,102], and to an important role for other competing exchange interactions. In the light of this, the properties of that particular material seem to be connected with a different phase boundary from that associated with SL ⊥ [46,49,76]. That being said, some of the similarities between our results and the Ybpyrochlores are striking: a gapless continuum of spin ex-citations, oblivious to thermodynamic phase transitions [42,[103][104][105] [Fig. 5(a)], and robust in temperature up to a broad feature in specific heat [42] (in the present Article, between pHAF and SL ⊥ ). And while the magnetic order in Yb-pyrochlores is, at least partially, an order of dipolar moments [28,37,38,[106][107][108], recent experiments have indicated that the primary order parameter may be "hidden", and distinct from a standard dipole order [42]. Thus, while the specific case developed in this Article probably does not apply to the Yb-pyrochlores, related physics may be at play.
Furthermore, since the spin-nematic phase SN ⊥ is found within the spin liquid SL ⊥ [ Fig. 1], this work provides a prototype for the peaceful co-existence of emergent gauge fields and long-range order. In this sense, SN ⊥ is an interesting new addition to the other phases where gauge fluctuations and broken symmetries co-exist, such as the Coulombic ferromagnet [17,109], and states with magnetic moment fragmentation [110], as recently observed in Nd 2 Zr 2 O 7 [43,111] and Ho 2 Ir 2 O 7 [112].
We also note that many other magnetic systems outside the rare earth oxides R 2 X 2 O 7 feature moments located on a pyrochlore lattice. Of particular interest are materials such as NaCaCo 2 F 7 and NaSrCo 2 F 7 [113,114] which boast XY like interactions with much higher energy scales than observed in the rare-earth oxides. If such a case could be found with frustrated transverse coupling J ± < 0 then it would render the physics discussed here accessible at a much more amenable temperature range.
In almost all spin-liquid candidates, the role of quenched structural and chemical disorder is an important issue [34,39,96,97,[115][116][117][118][119][120]. Depending on the type and strength of disorder, its consequences can vary. It is worth noting however, that disorder is not necessarily deleterious to spin liquid physics. It is known, for example, that weak disorder in non-Kramers pyrochlores, which leads to splittings in the low energy non-Kramers doublet can actually play a role in promoting a U (1) spin liquid ground state [24]. The spin-liquid states discussed in this manuscript do not depend on the translational symmetry of the Hamiltonian, but rather on the emergent gauge symmetries which arise from the local constraints in the ground state [Eq. (5)]. Thus, as long as the disorder is not so strong that these constraints are strongly violated, the essence of the spin liquids should be maintained in the presence of disorder, at least at finite temperature. For sufficiently strong disorder, or sufficiently low temperature, disorder may lead to orderby-disorder or glassiness [121,122]. A quantitative study of the effects of disorder on the phase diagram in Fig. 1 is a large undertaking and is beyond the scope of the present work but may be an interesting direction for future consideration.

VI. SUMMARY AND CONCLUSIONS
"Quantum spin ice", in which magnetic ions on a pyrochlore lattice interact through highly-anisotropic exchange interactions, have become an important paradigm in the search for quantum spin liquids. In this Article we have used large-scale classical Monte-Carlo simulation to explore the physics of the canonical model of a quantum spin ice, the XXZ model on a pyrochlore lattice H XXZ [Eq. (1)]. We find that this model has far more to offer than spin ice alone, supporting three distinct types of spin liquid, each with a different emergent gauge symmetry. Each of these spin liquids has a different signature in neutron scattering [ Fig. 2]. And the states found include a completely new form of spin liquid, described by a U(1)×U(1) gauge theory. At low temperatures this novel spin-liquid undergoes a phase transition to a state with hidden spin-nematic order [ Fig. 1], but retaining algebraic correlations of the spin dipoles. We have studied the excitations of this phase using state-of-the-art dynamical simulations, revealing a sharply defined Goldstone mode which would be hidden from conventional neutron scattering techniques.
So far as experiment is concerned, the main lesson of these results is that "quantum spin-ice" materials, can play host to a great many different spin-liquid and (hidden-)order phases, even where they are described by a Hamiltonian as simple as H XXZ [Eq. (1)]. This reinforces the point that pinch-points in pyrochlore magnets need not imply spin ice [66,73]. The existence of a sharp Goldstone mode in the nematic phase SN ⊥ also serves as a salutary reminder that broad, non-dispersing continua of excitations can hide a multitude of secrets [ Fig. 5].
From the theoretical point of view, this work identifies a new spin liquid, a novel spin nematic phase, and opens an interesting new perspective on the way in which different spin liquids can compete. The effect of quantum fluctuations on the phase diagram shown in Fig. 1 for J ± < 0 remains a subject for future study. However, experience with QMC simulation of H XXZ [Eq. (1)] for J ± > 0 suggests that quantitative values of the crossover temperature T * 2 and T * 3 may be substantially renormalized, but that the qualitative structure of the phase diagram should remain the same down to very low temperatures [14,21,23,26]. The high-symmetry point, J ± /J zz = −1/2 is also a high-symmetry point for quantum spins, and so remains the anchor for the spin liquid pHAF. None the less, the fate of this U(1)×U(1)×U(1) spin liquid for quantum spins at T = 0 remains an open question [57,58,60,123]. And, to the best of our knowledge, quantum analogues of the new spin liquid, SL ⊥ , which has a U(1)×U(1) gauge structure, remain unexplored [124]. However it seems reasonable to speculate that quantum effects will enhance, rather than suppress, the fluctuations which drive SL ⊥ and pHAF, and that the phase SN ⊥ will survive as hidden quantum spin-nematic order, within a quantum spin liquid. And preliminary numerical results for the spin-1/2 model at high temper-atures are entirely consistent with the topology of the phase diagram shown in Fig. 1 [125]. All of these questions open exciting avenues for future research.

ACKNOWLEDGEMENTS
We are pleased to acknowledge many helpful discussions with Karim Essafi, Jaan Oitmaa and Rajiv Singh. This work was supported by the Theory of Quantum Matter Unit of the Okinawa Institute of Science and Technology Graduate University (OIST). Numerical calculations were carried out using HPC facilities provided by OIST.

Appendix A: Definition of local coordinate frame
We describe the local-coordinate frame which is defined for four spins on a pyrochlore tetrahedron S 0 , S 1 , S 2 , S 3 occupying positions where a is the length of a cubic, 16-site unit cell of the pyrochlore lattice. The pseudospins in the global, crystal, coordinate frame S i relate to the pseudospins in the local frame S i [Eq. (1)] as and We have used this relationship between the local coordinate frame of Eq. (1) and the crystal coordinate frame to plot a representation of the bond quadrupolar order in real space in Fig. 1(b). The ellipsoid on each bond ij in Fig. 1(b) has principal axes aligned along the cubic axes of the pyrochlore lattice, with the length of each principal axis given by where c = 0.08 is chosen to make the figure readable and S α i is the component of spin i in the global frame, along crystal axis α = x, y, z.  Thermodynamics of the QSI in the region of spin-nematic order. (a) Specific heat cV (T ), showing an upturn followed by a shallow maximum in the region of the crossover into the spin liquid SL ⊥ at T * 2 /Jzz ∼ 10 −1 , and small peak associated with the onset of spin-nematic order at TSN/Jzz ≈ 10 −2 . (b) Correlation function T χ T 1 ice (T ) used to determine the crossover temperature into the spin-liquid SL ⊥ . χ T 1 ice is the susceptibility of the field m T 1 ice defined in Appendix D. The crossover temperature T * 2 is defined as the point at which T χ T 1 ice drops below its infinite temperature limit with temperature set in linear scale for J ± /J zz > − 1 2 and 256 temperatures in logarithmic scale for J ± /J zz ≤ − 1 2 . The phase boundary of the antiferromagnetically ordered (AF ⊥ ) phase, T N , was extracted from the susceptibility of the relevant order parameter, m E , as defined in Appendix D.
The crossover scale for the spin-ice regime (SI), T * 1 , was estimated from the Schottky-like peak in the heat capacity.
The crossover scale T * 3 for the spin-liquid pHAF was estimated from the Curie-Law crossover shown in Fig. 6.
For J ± < − 1 2 , the crossover scale T * 2 is associated with a weakening of the correlations of the local z-components of the spins. This can be observed by measuring the susceptibility, χ T1Ice (T ), of the field m T1Ice , defined in Appendix D]. Decreasing the temperature for −1 < J± Jzz < −0.5 the quantity T χ T1Ice (T ) first increases during the crossover from the paramagnet to pHAF and then drops as the system enters SL ⊥ . We define the crossover temperature T * 2 as the point at which the quantity T χ T1ice (T ) drops below its infinite temperature value This is illustrated in Fig. 7b. For quantum, spin-1/2 moments, the onset of spinnematic order is heralded by the bond-based order parameter Eq. (2). However for the purpose of classical simulation, it is sufficient to consider the simpler, sitebased order parameter Please note the prefactor of 4 in the definition of the spin nematic site order parameter [Eq. (B2)] to account for the spin length |S i | = 1 2 . The onset of spin-nematic order in simulations can be observed in either the site-based [cf. Fig. 7d] or bond-based order parameters [cf. Fig. 7c]. However, for simplicity, the values of the spin-nematic ordering temperature T SN shown in Fig. 1, were extracted from the peak in the order-parameter susceptibility associated with the site-based order parameter, Eq. (B2) [cf. Fig. 7e] Fig. 7 is obtained using 300 temperatures in logarithmic scale coverring 3 orders of magnitude, parallel tempering every 100 Monte Carlo steps, simulated annealing and thermalization at temperature T for 100000 Monte Carlo steps each. Measurements consist of 100000 different samples with 10 Monte Carlo steps between each sample. Error bars were estimated by comparing the results of three independent runs of the simulation. In Fig. 2 we show predictions for neutron scattering experiments, based on the equal-time (i.e. energyintegrated) structure factor where the dynamical structure factor S(q, ω) is defined through and the g βγ i is the g-tensor written in the local coordinate frame [49]. For simplicity, we have here taken g βγ i = 2δ βγ for all of the calculations in this paper. R αβ i is a rotation matrix which rotates from the local coordinate frame on site i, to the global, crystal coordinate frame. The definition of the local coordinate frame is given in Appendix A. Results for S(q) are shown in the left half-panels of Fig. 2. These results were taken from classical MC simulations of H XXZ at a given temperature, with further averaging provided by numerically integrating the semi-classical equations of motion for the spins. This secondary molecular-dynamics (MD) simulation was carried out using methods described in Ref. [78]. It is also useful to decompose the structure factor into the spin-flip (SF) and non spin-flip (NSF) channels measured in polarised neutron-scattering experiments.
wheren is the direction of polarization of the neutron magnetic moment. Following Fennell et al. [61], we takê n = (1, −1, 0)/ √ 2. Simulation results for S SF (q) and S NSF (q) are shown in the right half-panels of Fig. 2.
We have also used MD simulation to calculate the dynamical structure factor S(q, ω). Results for S(q, ω) within the spin-nematic phase of the quantum spin ice model are shown in Fig. 5a. Further details of the calculation of dynamical properties can be found in Appendix F.

Appendix D: Definitions of local order parameter fields
The definitions of the local order parameter fields m λ which appear in the theory of the spin liquid SL ⊥ [Section III] are given in Table I.
Definition in terms of spins within tetrahedron mA 2 Table I. Order-parameter fields m λ , derived from irreducible representations (irreps) of the tetrahedral point-group T d . Spin components Si = (S x i , S y i , S z i ) are written in the local frame of the magnetic ions, see Appendix A for a definition of this coordinate frame. The convention for the labelling of the spins with an tetrahedron is given in Appendix A.
Here we give the definitions in terms of the spins written in the local coordinate frame S i (defined in Appendix A), cf. Ref. 49 where definitions are given in the global, crystal basis.  Table I.
The tetrahedra of the pyrochlore lattice may be divided into two sets A and B. The centres of each set of tetrahedra each form an FCC lattice.
To calculate S αβ Bµ we use Eq. (12) where B µ (q) is defined as the lattice Fourier transform of B µ (r) over only the A sublattice of tetrahedra.
where N u.c. is the number of unit cells in the system. Simulations were carried out using local spin updates, augmented by over-relaxation, within a parallel tempering scheme with 300 temperatures distributed on a log scale between T = 0.003 J zz and T = 0.1 J zz . Thermalisation was accomplished through a process of simulated annealing, with 10 4 Monte Carlo steps (MCs) of annealing from high temperature to temperature T , followed by 10 4 MCs of thermalization at temperature T , and 10 5 MCs of measurements at temperature T . Spin configurations were sampled every 100 MCs during the measurements, giving an ensemble of 1000 samples. showing the expected behaviour (q) = q at small q. (b) Intensity I(q) of the peak as a function of momentum q. The dashed line shows the expected behaviour at finite temperature, I(q) ∝ 1/q 2 . Results are taken from moleculardynamics simulations of a cluster of N = 65536 spins, for J±/Jzz = −1.0, T = 0.002Jzz. Momentum q is measured relative to q = (0, 0, 0).

Appendix F: Dynamics of excitations in the spin-nematic phase
To study the Goldstone mode associated to the development of spin-nematic order, we calculate the dynami-cal correlation function.
where fluctuations of spin-nematic order are given by and the order parameter Q site ⊥ (r i , t) is defined through Eq. (B2).
χ Q site ⊥ (q, ω) is calculated numerically from 200 sample configurations extracted from Monte Carlo simulations on a system of linear size L = 16. We used 20000 steps for the simulated annealing spaced by 10 Monte Carlo steps between each simulated annealing step. The other parameters for the thermalization and parallel tempering are identical to the parameters used to calculate the phase diagram [Appendix B].
The ensemble of configurations obtained from Monte Carlo is then evolved in time according to the equation of motion, is the effective exchange field acting on site i, J ij is the anisotropic exchange interaction tensor and the sum in Eq. (F4) runs over the neighbors of i. The numerical integration of this nonlinear equation of motion proceeds as described in Ref. [78]. In Fig 8 we plot the dispersion of the Goldstone mode found in molecular dynamics (MD) simulations of H XXZ within the spin-nematic phase for J ± /J zz = −1. The dispersion was extracted from the position of low-energy dispersing peak in χ Q ⊥ (q, ω), as shown in the inset to Fig. 5b.