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

TL;DR

This paper investigates a simple quantum spin ice model on the pyrochlore lattice, revealing three competing spin liquid phases and hidden nematic order, with predictions for experimental signatures like neutron scattering.

Contribution

It uncovers multiple competing spin liquids and hidden nematic order in a basic quantum spin ice model, expanding understanding of frustrated magnetic systems.

Findings

Identified three distinct spin liquid phases.

Discovered a phase with hidden spin-nematic order.

Provided predictions for experimental neutron scattering signatures.

Abstract

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.

Figures

Tables

Table 1. Table 1: Order–parameter fields 𝐦 λ subscript 𝐦 𝜆 {\bf m}_{\lambda} , derived from irreducible representations (irreps) of the tetrahedral point-group T d . Spin components 𝖲 i = ( 𝖲 i x , 𝖲 i y , 𝖲 i z ) subscript 𝖲 𝑖 subscript superscript 𝖲 𝑥 𝑖 subscript superscript 𝖲 𝑦 𝑖 subscript superscript 𝖲 𝑧 𝑖 \mathsf{S}_{i}=(\mathsf{S}^{x}_{i},\mathsf{S}^{y}_{i},\mathsf{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 .

	Definition in terms of spins within tetrahedron
$m_{{𝖠}_{\mathrm{𝟤}}}$	$\frac{1}{2}\left({𝖲}_0^z+{𝖲}_1^z+{𝖲}_2^z+{𝖲}_3^z\right)$
${𝐦}_{𝖤}$	$\frac{1}{2}\left(\begin{array}{c}{𝖲}_0^x+{𝖲}_1^x+{𝖲}_2^x+{𝖲}_3^x\\ {𝖲}_0^y+{𝖲}_1^y+{𝖲}_2^y+{𝖲}_3^y\end{array}\right)$
${𝐦}_{{𝖳}_{\mathrm{𝟣},\mathrm{𝗂𝖼𝖾}}}$	$\frac{1}{2}\left(\begin{array}{c}{𝖲}_0^z+{𝖲}_1^z-{𝖲}_2^z-{𝖲}_3^z\\ {𝖲}_0^z-{𝖲}_1^z+{𝖲}_2^z-{𝖲}_3^z\\ {𝖲}_0^z-{𝖲}_1^z-{𝖲}_2^z+{𝖲}_3^z\end{array}\right)$
${𝐦}_{{𝖳}_{\mathrm{𝟣},\mathrm{𝗉𝗅𝖺𝗇𝖺𝗋}}}$	$\left(\begin{array}{c}\frac{1}{2}\left({𝖲}_0^x+{𝖲}_1^x-{𝖲}_2^x-{𝖲}_3^x\right)\\ \frac{1}{4}\left(-{𝖲}_0^x+\sqrt{3}{𝖲}_0^y+{𝖲}_1^x-\sqrt{3}{𝖲}_1^y-{𝖲}_2^x+\sqrt{3}{𝖲}_2^y+{𝖲}_3^x-\sqrt{3}{𝖲}_3^y\right)\\ \frac{1}{4}\left(-{𝖲}_0^x-\sqrt{3}{𝖲}_0^y+{𝖲}_1^x+\sqrt{3}{𝖲}_1^y+{𝖲}_2^x+\sqrt{3}{𝖲}_2^y-{𝖲}_3^x-\sqrt{3}{𝖲}_3^y\right)\end{array}\right)$
${𝐦}_{{𝖳}_{\mathrm{𝟤}}}$	$\left(\begin{array}{c}\frac{1}{2}\left(-{𝖲}_0^y-{𝖲}_1^y+{𝖲}_2^y+{𝖲}_3^y\right)\\ \frac{1}{4}\left(\sqrt{3}{𝖲}_0^x+{𝖲}_0^y-\sqrt{3}{𝖲}_1^x-{𝖲}_1^y+\sqrt{3}{𝖲}_2^x+{𝖲}_2^y-\sqrt{3}{𝖲}_3^x-{𝖲}_3^y\right)\\ \frac{1}{4}\left(-\sqrt{3}{𝖲}_0^x+{𝖲}_0^y+\sqrt{3}{𝖲}_1^x-{𝖲}_1^y+\sqrt{3}{𝖲}_2^x-{𝖲}_2^y-\sqrt{3}{𝖲}_3^x+{𝖲}_3^y\right)\end{array}\right)$

Equations

$${\mathcal{H}}_{\mathsf{XXZ}}=\sum _{⟨ij⟩}{J}_{\mathsf{zz}}{\mathsf{S}}_i^z{\mathsf{S}}_j^z-J_{\pm }\left({\mathsf{S}}_i^{+}{\mathsf{S}}_j^{-}+{\mathsf{S}}_i^{-}{\mathsf{S}}_j^{+}\right)$$

$${\mathbf{Q}}_{\perp }=\sum _{⟨ij⟩}\frac{1}{3N}\left(\begin{array}{c}{\mathsf{S}}_i^x{\mathsf{S}}_j^x-{\mathsf{S}}_i^y{\mathsf{S}}_j^y\\ {\mathsf{S}}_i^x{\mathsf{S}}_j^y+{\mathsf{S}}_i^y{\mathsf{S}}_j^x\end{array}\right),$$

$$F_{\text{SN⊥}}=a_2(T){\mathbf{Q}}_{\perp }^2+a_4{\mathbf{Q}}_{\perp }^4+\dots ,$$

$$\{{\mathbf{m}}_{\lambda }\}=\{m_{{\mathsf{A}}_2},{\mathbf{m}}_{\mathsf{E}},{\mathbf{m}}_{{\mathsf{T}}_1\mathsf{ice}},{\mathbf{m}}_{{\mathsf{T}}_1\mathsf{planar}},{\mathbf{m}}_{{\mathsf{T}}_2}\}$$

$$m_{{\mathsf{A}}_2}(\mathbf{r})=0,{\mathbf{m}}_{\mathsf{E}}(\mathbf{r})=0,{\mathbf{m}}_{{\mathsf{T}}_1\mathsf{ice}}(\mathbf{r})=0\forall \mathbf{r}$$

$$\delta {\mathcal{H}}_{\pm \pm }=\sum _{⟨ij⟩}{J}_{\pm \pm }\left[{\gamma }_{ij}{\mathsf{S}}_i^{+}{\mathsf{S}}_j^{+}+{\gamma }_{ij}^{\ast }{\mathsf{S}}_i^{-}{\mathsf{S}}_j^{-}\right],$$

$${\mathcal{H}}_{\mathsf{ex}}={\mathcal{H}}_{\mathsf{XXZ}}+\delta {\mathcal{H}}_{\pm \pm }.$$

$${\mathbf{B}}_1=\frac{1}{2}(2m_{{\mathsf{T}}_1\mathsf{planar}}^x,-\sqrt{3}{m}_{{\mathsf{T}}_2}^y-m_{{\mathsf{T}}_1\mathsf{planar}}^y,$$

$$\sqrt{3}{m}_{{\mathsf{T}}_2}^z-m_{{\mathsf{T}}_1\mathsf{planar}}^z)$$

$${\mathbf{B}}_2=\frac{1}{2}(2m_{{\mathsf{T}}_1\mathsf{planar}}^x,-m_{{\mathsf{T}}_2}^y+\sqrt{3}{m}_{{\mathsf{T}}_1\mathsf{planar}}^y,$$

$$-m_{{\mathsf{T}}_2}^z-\sqrt{3}{m}_{{\mathsf{T}}_1\mathsf{planar}}^z),$$

$$\nabla \cdot {\mathbf{B}}_1=0,\nabla \cdot {\mathbf{B}}_2=0.$$

$${\mathbf{B}}_1=\nabla \times {\mathbf{A}}_1,{\mathbf{B}}_2=\nabla \times {\mathbf{A}}_2$$

$${\mathcal{F}}_{{\text{SL}}_{\perp }}=\frac{T}{V}\int d^{3}r\lambda ({\mathbf{B}}_1^2+{\mathbf{B}}_2^2)$$

$$=\frac{T}{V}\int d^{3}r\lambda \left[(\nabla \times {\mathbf{A}}_1{)}^2+(\nabla \times {\mathbf{A}}_2{)}^2\right]$$

$$S_{{\mathbf{B}}_{\mu }}^{\alpha \beta }(\mathbf{q})=⟨B_{\mu }^{\alpha }(-\mathbf{q})B_{\mu }^{\beta }(\mathbf{q})⟩$$

$${\mathcal{F}}_{{\text{SL}}_{\perp }}\to {\mathcal{F}}_{{\text{SL}}_{\perp }}+\delta {\mathcal{F}}_{{\text{SL}}_{\perp }}$$

$$\delta {\mathcal{F}}_{{\text{SL}}_{\perp }}=\frac{T}{V}\int d^{3}r{\lambda }^{\prime }\{{\left(B_1^x\right)}^2-\frac{1}{2}\left[{\left(B_1^y\right)}^2+{\left(B_1^z\right)}^2\right]$$

$$-{\left(B_2^x\right)}^2+\frac{1}{2}\left[{\left(B_2^y\right)}^2+{\left(B_2^z\right)}^2\right]-\sqrt{3}\left[B_1^z{B}_2^z-B_1^y{B}_2^y\right]\}$$

$$\omega \approx v_Q|\mathbf{q}-{\mathbf{q}}_{\mathsf{rl}}|$$

$${\mathbf{r}}_0=\frac{a}{8}\left(1,1,1\right)$$

$${\mathbf{S}}_i={\mathbf{x}}_i^{\mathsf{local}}{\mathsf{S}}_i^x+{\mathbf{y}}_i^{\mathsf{local}}{\mathsf{S}}_i^y+{\mathbf{z}}_i^{\mathsf{local}}{\mathsf{S}}_i^z$$

$${\mathbf{z}}_0^{\mathsf{local}}=\frac{1}{\sqrt{3}}(1,1,1)$$

$${\mathbf{x}}_0^{\mathsf{local}}=\frac{1}{\sqrt{6}}(-2,1,1)$$

$${\mathbf{y}}_0^{\mathsf{local}}=\frac{1}{\sqrt{2}}(0,-1,1)$$

$$l_{\alpha }=4(c+⟨S_i^{\alpha }{S}_j^{\alpha }⟩)$$

$$T_2^{\ast }{\chi }_{{\mathsf{T}}_1\mathsf{ice}}(T_2^{\ast })=\underset{T\to \infty }{\lim }T{\chi }_{{\mathsf{T}}_1\mathsf{ice}}(T)$$

$${\mathbf{Q}}_{\perp }^{\mathsf{site}}=\frac{4}{N}\sum _i\left(\begin{array}{c}{{\mathsf{S}}_i^x}^2-{{\mathsf{S}}_i^y}^2\\ 2{\mathsf{S}}_i^x{\mathsf{S}}_i^y\end{array}\right).$$

$${\chi }_{Q_{\perp }^{\mathsf{site}}}=\frac{N}{T}\left(⟨{{\mathbf{Q}}_{\perp }^{\mathsf{site}}}^2⟩-⟨|{\mathbf{Q}}_{\perp }^{\mathsf{site}}|{⟩}^2\right).$$

$$S(\mathbf{q})=\int d\omega S(\mathbf{q},\omega ),$$

$$S(\mathbf{q},\omega )=\sum _{\alpha \beta }\left({\delta }_{\alpha \beta }-\frac{q_{\alpha }{q}_{\beta }}{q^2}\right)⟨m^{\alpha }(-\mathbf{q},\omega )m^{\beta }(\mathbf{q},\omega )⟩$$

$$m^{\alpha }(\mathbf{q},\omega )=\sum _{i\beta \gamma }{R}_i^{\alpha \beta }{g}_i^{\beta \gamma }\left(\int {\mathsf{S}}_i^{\gamma }(t)e^{i\omega t}dt\right)e^{i\mathbf{q}\cdot {\mathbf{r}}_i},$$

$$S_{\text{SF}}(\mathbf{q})=\frac{1}{q^2}\int d\omega ⟨|\mathbf{m}(\mathbf{q},\omega )\cdot (\hat{\mathbf{n}}\times \mathbf{q}){|}^2⟩$$

$$S_{\text{NSF}}(\mathbf{q})=\int d\omega ⟨|\mathbf{m}(\mathbf{q},\omega )\cdot \hat{\mathbf{n}}{|}^2⟩,$$

Full text

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 Anderson (1973); Lee (2008); Balents (2010). 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 [Gardner et al., 2010]. Perhaps the most celebrated of these is the “spin ice” found in the Ising magnets Ho2Ti2O7 and Dy2Ti2O7 [Bramwell and Gingras, 2001]; a classical spin liquid, described by an emergent U(1) lattice gauge theory with magnetic monopole excitations [Castelnovo et al., 2012].

As the understanding of spin ice has grown, so more attention has been given to the role of quantum effects. These are of particular relevance where a spin–ice arises through anisotropic exchange interactions in a pyrochlore magnet Curnoe (2007); Molavian et al. (2007); Onoda and Tanaka (2010, 2011); Ross et al. (2011); Gingras and McClarty (2014), and have the potential to convert classical spin ice into a massively–entangled, quantum, spin liquid, described by an emergent U(1) quantum electrodynamics [Hermele et al., 2004; Banerjee et al., 2008; Shannon et al., 2012; Benton et al., 2012; Savary and Balents, 2012; Lee et al., 2012; Gingras and McClarty, 2014; Hao et al., 2014; McClarty et al., 2015; Kato and Onoda, 2015; Chen, 2016; Shannon, 2017; Savary and Balents, 2017a; Chen, 2017a; Huang et al., ]. At the same time, great progress has been made in synthesizing and characterizing magnetic pyrochlore oxides. As well as revealing a number of candidates for quantum spin–ice behaviour Zhou et al. (2008); Chang et al. (2012); Fennell et al. (2012); Kimura et al. (2013); Sibille et al. (2015, 2016); Anand et al. (2016); Wen et al. (2017); Sibille et al. , these experiments have turned up many unusual and unexpected magnetic states in systems with strongly anisotropic exchange Dalmas de Réotier et al. (2006); Dun et al. (2013); Yaouanc et al. (2013); Taniguchi et al. (2013); Hallas et al. (2014, 2015, 2016a); Petit et al. (2016a, b); Takatsu et al. (2016); Hallas et al. .

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. 1d]. 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)$\times$U(1) gauge structure and, as an added bonus, undergoes a phase transition into a state with hidden, spin–nematic, 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 Curnoe (2007); Molavian et al. (2007); Onoda and Tanaka (2010, 2011); Ross et al. (2011); Gingras and McClarty (2014) is the XXZ model on the pyrochlore lattice

[TABLE]

where $\mathsf{S}_{i}=(\mathsf{S}^{x}_{i},\mathsf{S}^{y}_{i},\mathsf{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, $\mathsf{S}^{z}_{i}$) lies on a local [111] axis, as defined in Appendix A.

Ising interactions, $J_{\sf 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_{\pm}$, introduces dynamics about these spin–ice configurations and, for larger, positive, values of $J_{\pm}/J_{\sf zz}$, drives the system into a state with easy–plane order Onoda (2011); Lee et al. (2012); Savary and Balents (2012); Wong et al. (2013); Yan et al. (2017); Benton (2015). The physical meaning of this easy–plane order depends on the nature of the magnetic ion. For Kramers ions like Yb3+ and Er3+ all components of $\sf S$ relate to a magnetic dipole moment Ross et al. (2011), and the ordered phase is an easy–plane antiferromagnet. However for non–Kramers ions such as Pr3+ and Tb3+ Onoda and Tanaka (2010); Petit et al. (2016b), or “dipolar-octupolar” Kramers ions like Nd3+ or Ce3+ Huang et al. (2014), 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_{\pm}>0$, $\mathcal{H}_{\sf 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 Banerjee et al. (2008); Kato and Onoda (2015); Shannon (2017). For $J_{\pm}/J_{\sf zz}\lesssim 0.05$, QMC simulations find a crossover from a conventional paramagnet into a classical spin–liquid (spin ice) at a temperature $T^{*}/J_{\sf zz}\sim 0.2$, and a second crossover into a quantum spin liquid (QSL) at a much lower temperature $T_{QSL}^{*}/J_{\sf zz}\sim(J_{\pm}/J_{\sf 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  Hermele et al. (2004); Benton et al. (2012). For $J_{\pm}/J_{\sf zz}\gtrsim 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 $\mathsf{S}^{z}$–axis Banerjee et al. (2008); Kato and Onoda (2015); Shannon (2017).

Very little is known about the properties of $\mathcal{H}_{\sf XXZ}$ for $J_{\pm}<0$ Onoda and Tanaka (2011); Lee et al. (2012); Petit et al. (2016b). On perturbative grounds, it is expected that the ground state for $|J_{\pm}|/J_{\sf zz}\ll 1$ will also be a U(1) QSL Hermele et al. (2004), albeit one with a modified spinon dispersion Lee et al. (2012); Chen (2017b). Gauge mean–field calculations suggest that this QSL persists over a broad range of parameters, $-4.13\lesssim J_{\pm}/J_{\sf zz}<0$ [Lee et al., 2012]. But the nature of competing ordered — or disordered — phases for $J_{\pm}<0$ remains an open question.

There are many reasons to believe that the properties of the quantum spin ice model, $\mathcal{H}_{\sf XXZ}$ [Eq. (1)] for frustrated coupling $J_{\pm}<0$, could be even richer than for $J_{\pm}>0$. In particular, for $J_{\pm}/J_{\sf zz}=-\frac{1}{2}$, $\mathcal{H}_{\sf 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 Moessner and Chalker (1998a, b); Isakov et al. (2004); Henley (2005), and it has also been argued to support a QSL ground state Canals and Lacroix (1998, 2000); Burnell et al. (2009); Huang et al. (2016). 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_{\pm}/J_{\sf zz}\approx 0$, and a state homologous to the HAF for $J_{\pm}/J_{\sf zz}\approx-\frac{1}{2}$. It also opens the door for yet more novel magnetic phases for $J_{\pm}/J_{\sf zz}<-\frac{1}{2}$.

II Phase diagram determined by classical Monte Carlo simulations

Since the quantum spin ice model, $\mathcal{H}_{\sf XXZ}$ [Eq. (1)], is inaccessible to QMC for $J_{\pm}<0$, we instead study its finite–temperature properties using classical Monte Carlo (MC) simulation — the results are summarised in the phase diagram Fig. 1d. For $J_{\pm}>0$, this phase diagram is very similar to that previously found in QMC simulations Banerjee et al. (2008); Kato and Onoda (2015); Shannon (2017) — at a qualitative level, the only significant difference is the absence of a QSL below $T_{QSL}^{*}/J_{\sf zz}\sim(J_{\pm}/J_{\sf zz})^{3}\lesssim 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 $\text{AF}_{\perp}$. 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 [Benton, ]. Further details of classical MC simulations for $J_{\pm}>0$ will be presented elsewhere [Taillefumier et al., ].

We now turn to the frustrated case, $J_{\pm}<0$. At low temperatures, spin–ice correlations persist up to $J_{\pm}/J_{\sf zz}=-\frac{1}{2}$ Onoda and Tanaka (2011); Petit et al. (2016b), as illustrated in Fig. 2a. Upon reaching $J_{\pm}/J_{\sf zz}=-\frac{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. 1d. 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 Jaubert et al. (2013).

While the correlations measured in the equal–time structure factor $S({\bf q})$ are also different from those found in the HAF Moessner and Chalker (1998b); Isakov et al. (2004); Conlon and Chalker (2009), the two models are equivalent up to a local coordinate transformation. And, by analogy with earlier work on the HAF Isakov et al. (2004); Henley (2005, 2010), the spin liquid pHAF can be described by a U(1)$\times$U(1)$\times$U(1) gauge theory.

The situation for $J_{\pm}/J_{\sf zz}<-\frac{1}{2}$ is even more interesting. 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 $\text{SL}_{\perp}$ in Fig. 1d. Spin correlations in this regime have algebraic character, with pinch–points in $S({\bf 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_{\text{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({\bf q})$ and, at least as far as dipolar spin correlations are concerned, the system remains disordered.

While the new phase for $T<T_{\text{SN}}$ — labelled $\text{SN}_{\perp}$ in Fig. 1d — 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 $\mathcal{H}_{\sf XXZ}$ [Eq. (1)] by selecting an axis within the local $xy$-plane. Such an order can be described by the bond–based order parameter Andreev and Grishchuk (1984); Chubukov (1991); Shannon et al. (2006)

[TABLE]

where the sum on $\langle ij\rangle$ runs over the nearest–neighbour bonds of the lattice, and $\mathsf{S}_{i}=(\mathsf{S}^{x}_{i},\mathsf{S}^{y}_{i},\mathsf{S}^{z}_{i})$ is expressed in the local frame of site $i$ (cf. Appendix A).

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 [Shannon et al., 2006, 2010; Smerald et al., 2015]. And in common with these systems, the associated Landau theory

[TABLE]

lacks a cubic term, and therefore permits a continuous phase transition. Simulations suggest that the phase transition at $T=T_{\text{SN}}$ is indeed continuous for $J_{\pm}/J_{\sf zz}\lesssim-\frac{1}{2}$, becoming first-order approaching the high–symmetry point $J_{\pm}/J_{\sf zz}\to-\frac{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 Henley (2005); Isakov et al. (2005); Castelnovo et al. (2012), which gives rise to characteristic “pinch–points” in the spin structure factor $S({\bf q})$ [Fig. 2a]. Meanwhile, for classical spins, spin correlations in the Heisenberg AF on the pyrochlore lattice — and by extension in the pHAF — can be described using a U(1)$\times$U(1)$\times$U(1) gauge theory Moessner and Chalker (1998a, b); Isakov et al. (2004); Henley (2005). 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, $\text{SL}_{\perp}$ [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 $\text{SL}_{\perp}$, too, may be described by some form of gauge theory.

We can develop a field–theory for the spin–liquid $\text{SL}_{\perp}$ by applying the methods developed in Benton et al. (2016); Benton (2015). The starting point of this approach is to recast the spins $\mathsf{S}_{i}$ in $\mathcal{H}_{\sf XXZ}$ [Eq. (1)] in terms of five order–parameter fields

[TABLE]

defined on each tetrahedron ${\bf r}$. These objects ${\bf m}_{\lambda}(\mathbf{r})$ describe the different kinds of four-sublattice magnetic order consistent with the point group symmetry of the pyrochlore lattice. Definitions of each field ${\bf m}_{\lambda}$ in terms of the spins ${\sf S}_{i}$ are given in Appendix D.

The most general exchange Hamiltonian on the pyrochlore lattice can be transcribed exactly in terms of ${\bf m}_{\lambda}$ Yan et al. (2017). 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 Benton et al. (2016); Benton (2015). In the case of $\text{SL}_{\perp}$, for $T\to 0$, we have

[TABLE]

The spin fluctuations at low temperature are thus dominated by the fluctuations of the remaining fields ${\bf m}_{\sf T_{2}}({\bf r})$ and ${\bf m}_{\sf T_{1}planar}({\bf r})$. These fields have significance as the order-parameters of the competing four-sublattice magnetic orders which would be induced by the symmetry–allowed perturbation

[TABLE]

where $\gamma_{ij}$ are complex phase factors arising from the change in coordinate frame between different lattice sites Curnoe (2007); McClarty et al. (2009); Onoda (2011); Onoda and Tanaka (2010, 2011); Ross et al. (2011). For this reason, the spin–liquid $\text{SL}_{\perp}$ falls very naturally into the “multiple–phase competition” scenario for pyrochlore magnets Yan et al. ; Jaubert et al. (2015); Benton (2015); Yan et al. (2017).

In Fig. 3, we show the classical ground–state phase diagram of anisotropic exchange model

[TABLE]

This contains three distinct regions of 4–sublattice order : the easy–plane ordered phases described by the fields ${\bf m}_{\sf E}$ (denoted $\text{AF}_{\perp}$ in Fig. 1d), ${\bf m}_{\sf T_{1}planar}$, and ${\bf m}_{\sf T_{2}}$ (Palmer–Chalker state Palmer and Chalker (2000)). These border a region of spin ice (denoted SI in Fig. 1d), dominated by fluctuations of ${\bf m}_{\sf T_{1}ice}$. We note that a closely–related phase diagram has been derived for non–Kramers ions Onoda and Tanaka (2011); Petit et al. (2016b); 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 $\text{SL}_{\perp}$ arise from the fact that neighbouring tetrahedra share a spin, so that the fields ${\bf m}_{\lambda}(\mathbf{r})$ on neighbouring tetrahedra are not independent of one another. This point, combined with Eq. (5), imposes spatial constraints on the fluctuations of ${\bf m}_{\sf T_{2}}({\bf r})$ and ${\bf m}_{\sf T_{1}planar}({\bf r})$. After coarse graining to extract the long wavelength physics these constraints may be written in terms of two, independent, vector fluxes

[TABLE]

which each separately obey their own Gauss’ law

[TABLE]

We can therefore write

[TABLE]

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 Henley (2010). The only choice of Gaussian free-energy consistent with both the point group symmetry and the $U(1)$ symmetry of $\mathcal{H}_{\sf XXZ}$ is

[TABLE]

where the coefficient $\lambda$ can be determined through fits to simulation, or a large–N expansion Isakov et al. (2004); Benton et al. (2016).

It follows from the existence of the conserved fluxes ${\bf B}_{1}$ and ${\bf B}_{2}$ and the free-energy Eq. (11) that $\text{SL}_{\perp}$ is a Coulomb phase with algebraic correlations Henley (2010). The validity of this description is demonstrated in Fig. 4 where we compare analytic calculations of the flux structure factor

[TABLE]

based on Eq. (11) with the results of Monte Carlo simulation. Pinch point singularities are clearly seen in both analytic and numerical calculations. It is the same fluctuations of ${\bf B}_{1}$ and ${\bf B}_{2}$ which are responsible for the characteristic pinch–point structures in the (spin) structure factor measured by neutron scattering, as shown in Fig. 2b.

At finite temperature, we anticipate that the spin liquid $\text{SL}_{\perp}$ will be perturbatively stable against terms such as $\delta\mathcal{H}_{\pm\pm}$ [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 :

[TABLE]

This form of free energy will still lead to pinch points in $S_{{\bf B}_{\mu}}^{\alpha\beta}(\mathbf{q})$ and $S({\bf q})$, but these will take on a more anisotropic character.

IV Dynamics in the spin–nematic phase

For temperatures, $T<T_{\text{SN}}$ the easy–plane spin-liquid ($\text{SL}_{\perp}$) gives way to a phase with hidden spin–nematic order, labelled $\text{SN}_{\perp}$ in Fig. 1d. As far as the dipole moments of spins are concerned, the spin–nematic phase is disordered, and neutron scattering experiments would reveal algebraic correlations, as in the spin liquid $\text{SL}_{\perp}$. However the pinch points in $S({\bf q})$ [cf Fig. 2b] hide a great wealth of interesting spin excitations.

To better understand the dynamics of the spin–nematic phase, we have calculated the dynamical structure factor $S({\bf q},\omega)$, within a semi–classical molecular–dynamics (MD) simulation, using the methods described in [Taillefumier et al., 2014]. Relevant definitions are given in Appendix C. For $\omega/J_{\sf zz}\lesssim 0.2$, $S({\bf q},\omega)$ presents a featureless, non–dispersing continuum [Fig. 5a]. Relics of dispersing excitations are visible in $S({\bf q},\omega)$ 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({\bf q},\omega)$ as a function of temperature, we find that results for $S({\bf q},\omega)$ in the spin–nematic phase for $T<T_{\text{SN}}$, are very similar to those found in the spin liquid $\text{SL}_{\perp}$ for $T>T_{\text{SN}}$.

Incoherent, non-dispersing structure of the type shown in Fig. 5a is reminiscent of theoretical predictions Knolle et al. (2014); Punk et al. (2014); Bieri et al. (2015) and experimental measurements Han et al. (2012); Petit et al. (2016b); Paddison et al. (2017), 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 Savary and Balents (2017b). In the semi–classical limit studied here, it is probably unsafe to attribute such a continuum to fractionalized excitations Taillefumier et al. . However, the fact that $S({\bf q},\omega)$ only records dipolar spin correlations obscures a more important fact — the spin–nematic order which is present for $T<T_{\text{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 $\chi_{\sf Q_{\perp}^{\sf site}}({\bf q},\omega)$, which measures fluctuations of the on-site quadrupolar moments, which are well defined for classical spins [cf. Appendix B]. A sharp excitation, with dispersion

[TABLE]

can now be resolved near to the zone centers with ${\bf q}_{\sf 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 $\text{SN}_{\perp}$ would occur in a quadrupolar structure factor, which might, in principle, be measured in resonant X–ray experiments Savary and Senthil .

At present, relatively little is known about the dynamical properties of spin–nematic states. Field–theoretic analysis [Smerald and Shannon, 2013; Starykh and Balents, 2014; Smerald et al., 2015; Smerald and Shannon, 2016], based on the symmetry of the order parameter, predicts that spin–nematics support gapless Goldstone modes, visible in $\chi_{Q_{\perp}}({\bf q},\omega)$. This Goldstone mode has dispersion $\omega\propto|{\bf q}|$ [cf. Eq. (15)], and at zero temperature the associated intensity diverges as $\sim 1/\omega$ for $\omega\to 0$ [Smerald et al., 2015]. The same behaviour is seen in “flavour–wave” calculations and QMC simulations of spin–1 models constructed to support quadrupolar order Tsunetsugu and Arikawa (2006); Läuchli et al. (2006); Völl and Wessel (2015). 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 Onishi (2015). 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 Shindou et al. (2013).

Our MD simulations of $\text{SN}_{\perp}$ clearly reveal a linearly–dispersing Goldstone mode, with intensity which diverges for $\omega\to 0$ [cf. inset to Fig. 5b]. The form of this divergence is $\sim 1/\omega^{2}$, rather than $\sim 1/\omega$. 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 Völl and Wessel (2015). 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 spin–nematic phase, and specifically the characterization of the Goldstone mode are given in Appendix F. We note that the U(1) symmetry of $\mathcal{H}_{\sf 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 $\text{SN}_{\perp}$ would acquire a small gap.

V Discussion

Spin liquids Anderson (1973); Lee (2008); Balents (2010) and spin nematics Andreev and Grishchuk (1984); Chubukov (1991); Shannon et al. (2006) 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 $\text{SL}_{\perp}$, and the spin–nematic $\text{SN}_{\perp}$, in the type of rare–earth pyrochlore magnet described by Eq. (1).

In the case of $\text{SN}_{\perp}$, 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 Pr3+ Onoda and Tanaka (2010); Lee et al. (2012); Petit et al. (2016b) and Tb3+ Onoda and Tanaka (2010); Guitteny et al. (2013); Taniguchi et al. (2013); Takatsu et al. (2016), and Kramers doublets of dipolar–octupolar character, such as Nd3+ Huang et al. (2014), and Ce3+ Li and Chen (2017). For these ions, quadrupolar or octupolar order may occur as a“classical” order of the transverse part of the pseudospins ${\sf S}_{i}$. The multipolar character of the order follows from the symmetry of the crystal–field ground–state (doublet) of the magnetic ion, which is described by ${\sf 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 $\text{SN}_{\perp}$, 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 Pr-based pyrochlores have the recommended single-ion and interaction anisotropies Onoda and Tanaka (2010); Lee et al. (2012); Petit et al. (2016b). Coupling parameters of Pr2Zr2O7 for example have been suggested to sit in the AF*⟂* phase of Fig. 1d Petit et al. (2016b), although it seems that the coupling of structural disorder to the non-Kramers doublets plays a significant role Savary and Balents (2017a); Wen et al. (2017); Benton (2017); Martin et al. . Since chemical pressure has already proven to be a useful tool to move a family of compounds across a phase diagram Dun et al. (2014); Wiebe and Hallas (2015); Jaubert et al. (2015); Hallas et al. (2016a); Yan et al. (2017), Pr2X2O7 (X=Sn,Hf,Pb) are promising candidates to investigate, with ferromagnetic correlations consistent with positive $J_{\sf zz}$ and no dipole order yet observed Matsuhira et al. (2004); Zhou et al. (2008); Hallas et al. (2015); Sibille et al. (2016); Anand et al. (2016).

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 $\mathcal{H}_{\sf XXZ}$ model of Eq. (1). Specifically, experiments on Yb2Ti2O7 point to an unfrustrated value of $J_{\pm}>0$ Ross et al. (2011); Robert et al. (2015); Thompson et al. (2017), 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 $\text{SL}_{\perp}$ Jaubert et al. (2015); Yan et al. (2017); Hallas et al. . That being said, some of the similarities between our results and the Yb–pyrochlores are striking: a gapless continuum of spin excitations, oblivious to thermodynamic phase transitions Ross et al. (2009); Maisuradze et al. (2015); Hallas et al. (2016a); Gaudet et al. (2016) [Fig. 5(a)], and robust in temperature up to a broad feature in specific heat Hallas et al. (2016a) (in the present Article, between pHAF and $\text{SL}_{\perp}$). And while the magnetic order in Yb-pyrochlores is, at least partially, an order of dipolar moments Yasui et al. (2003); Chang et al. (2012); Yaouanc et al. (2013); Dun et al. (2013); Lhotel et al. (2014); Hallas et al. (2016b), recent experiments have indicated that the primary order parameter may be “hidden”, and distinct from a standard dipole order Hallas et al. (2016a). 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 $\text{SN}_{\perp}$ is found within the spin liquid $\text{SL}_{\perp}$ [Fig. 1d], this work provides a prototype for the peaceful co-existence of emergent gauge fields and long–range order. In this sense, $\text{SN}_{\perp}$ is an interesting new addition to the other phases where gauge fluctuations and broken symmetries co-exist, such as the Coulombic ferromagnet Savary and Balents (2012); Powell (2015), and states with magnetic moment fragmentation Brooks-Bartlett et al. (2014), as recently observed in Nd2Zr2O7 Petit et al. (2016a); Benton (2016) and Ho2Ir2O7 Lefrançois et al. (2017).

We also note that many other magnetic systems outside the rare earth oxides R2X2O7 feature moments located on a pyrochlore lattice. Of particular interest are materials such as NaCaCo2F7 and NaSrCo2F7 Ross et al. (2016, 2017) 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_{\pm}<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 Ross et al. (2012); Wen et al. (2017); Benton (2017); Mostaed et al. (2017); Taniguchi et al. (2013); Kermarrec et al. (2014); Han et al. (2016); Norman (2016); Li et al. (2017); Martin et al. . 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 Savary and Balents (2017a). 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 order–by–disorder or glassiness Bellier-Castella et al. (2001); Saunders and Chalker (2007). A quantitative study of the effects of disorder on the phase diagram in Fig. 1d 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 $\mathcal{H}_{\sf 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)$\times$U(1) gauge theory. At low temperatures this novel spin–liquid undergoes a phase transition to a state with hidden spin–nematic order [Fig. 1d], 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 $\mathcal{H}_{\sf XXZ}$ [Eq. (1)]. This reinforces the point that pinch–points in pyrochlore magnets need not imply spin ice Henley (2010); Benton et al. (2016). The existence of a sharp Goldstone mode in the nematic phase $\text{SN}_{\perp}$ 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. 1d for $J_{\pm}<0$ remains a subject for future study. However, experience with QMC simulation of $\mathcal{H}_{\sf XXZ}$ [Eq. (1)] for $J_{\pm}>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 Banerjee et al. (2008); Kato and Onoda (2015); Shannon (2017); Huang et al. . The high–symmetry point, $J_{\pm}/J_{\sf 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)$\times$U(1)$\times$U(1) spin liquid for quantum spins at $T=0$ remains an open question Canals and Lacroix (1998, 2000); Tsunetsugu (2001); Huang et al. (2016). And, to the best of our knowledge, quantum analogues of the new spin liquid, $\text{SL}_{\perp}$, which has a U(1)$\times$U(1) gauge structure, remain unexplored Ueda et al. (2016). However it seems reasonable to speculate that quantum effects will enhance, rather than suppress, the fluctuations which drive $\text{SL}_{\perp}$ and pHAF, and that the phase $\text{SN}_{\perp}$ 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 temperatures are entirely consistent with the topology of the phase diagram shown in Fig. 1d Singh and Oitmaa . 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 ${\sf S_{0}}$, ${\sf S_{1}}$, ${\sf S_{2}}$, ${\sf S_{3}}$ occupying positions

[TABLE]

where $a$ is the length of a cubic, 16-site unit cell of the pyrochlore lattice.

The pseudospins in the global, crystal, coordinate frame ${\bf S}_{i}$ relate to the pseudospins in the local frame ${\sf S}_{i}$ [Eq. (1)] as

[TABLE]

Where

[TABLE]

[TABLE]

and

[TABLE]

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. 1d(b). The ellipsoid on each bond $ij$ in Fig. 1d(b) has principal axes aligned along the cubic axes of the pyrochlore lattice, with the length of each principal axis given by

[TABLE]

where $c=0.08$ is chosen to make the figure readable and $S^{\alpha}_{i}$ is the component of spin $i$ in the global frame, along crystal axis $\alpha=x,y,z$.

Appendix B Details of the numerical determination of the phase diagram

The phase diagram shown in Fig. 1d was extracted from classical Monte Carlo (MC) simulations of the quantum spin ice model, $\mathcal{H}_{\sf XXZ}$ [Eq. (1)]. Spins were treated as classical vectors with fixed length $|{\sf S}_{i}|=\frac{1}{2}$. These simulations were carried out for a cubic cluster of 8192 spins, using a single spin flip algorithm combined with simulated annealing, parallel tempering and over–relaxation. The phase diagram is obtained using 75000 simulated annealing steps using a stepwise decrease of temperature starting from $T=10J_{\sf zz}$ down to the target temperature. Each annealing step consists of 10 Monte Carlo steps (a Monte Carlo step consists of a full sweep of the lattice combined with over–relaxation). The simulated annealing is followed by 1000 parallel tempering steps with 500 Monte Carlo steps in between, and then by 200000 Monte Carlo steps for thermalization at fixed temperature. Measurements consist of 200000 samples separated by 10 Monte Carlo steps and combined to parallel tempering every 50 measures. We use 256 different replicas with temperature set in linear scale for $J_{\pm}/J_{\sf zz}>-\frac{1}{2}$ and 256 temperatures in logarithmic scale for $J_{\pm}/J_{\sf zz}\leq-\frac{1}{2}$.

The phase boundary of the antiferromagnetically ordered (AF*⟂*) phase, $T_{N}$, was extracted from the susceptibility of the relevant order parameter, ${\bf m}_{\sf 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}^{\ast}$ for the spin–liquid pHAF was estimated from the Curie-Law crossover shown in Fig. 6.

For $J_{\pm}<-\frac{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, $\chi_{\sf T_{1}Ice}(T)$, of the field ${\bf m}_{\sf T_{1}Ice}$, defined in Appendix D]. Decreasing the temperature for $-1<\frac{J_{\pm}}{J_{zz}}<-0.5$ the quantity $T\chi_{\sf T_{1}Ice}(T)$ first increases during the crossover from the paramagnet to pHAF and then drops as the system enters $\text{SL}_{\perp}$. We define the crossover temperature $T_{2}^{\ast}$ as the point at which the quantity $T\chi_{\sf T_{1}ice}(T)$ drops below its infinite temperature value

[TABLE]

This is illustrated in Fig. 7b.

For quantum, spin–1/2 moments, the onset of spin–nematic 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, site–based order parameter

[TABLE]

Please note the prefactor of 4 in the definition of the spin nematic site order parameter [Eq. (23)] to account for the spin length $|{\sf S}_{i}|=\frac{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_{\text{SN}}$ shown in Fig. 1d, were extracted from the peak in the order–parameter susceptibility

[TABLE]

associated with the site–based order parameter, Eq. (23) [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.

Appendix C Definitions of dynamical structure factors

In Fig. 2 we show predictions for neutron scattering experiments, based on the equal–time (i.e. energy–integrated) structure factor

[TABLE]

where the dynamical structure factor $S({\bf q},\omega)$ is defined through

[TABLE]

and the $g^{\beta\gamma}_{i}$ is the $g$-tensor written in the local coordinate frame Yan et al. (2017). For simplicity, we have here taken $g^{\beta\gamma}_{i}=2\delta_{\beta\gamma}$ for all of the calculations in this paper. $R^{\alpha\beta}_{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({\bf q})$ are shown in the left half–panels of Fig. 2. These results were taken from classical MC simulations of $\mathcal{H}_{\sf 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. [Taillefumier et al., 2014].

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.

[TABLE]

where $\hat{\bf n}$ is the direction of polarization of the neutron magnetic moment. Following Fennell et al. [Fennell et al., 2009], we take $\hat{\bf n}=(1,-1,0)/\sqrt{2}$. Simulation results for $S_{\text{SF}}({\bf q})$ and $S_{\text{NSF}}({\bf q})$ are shown in the right half–panels of Fig. 2.

We have also used MD simulation to calculate the dynamical structure factor $S({\bf q},\omega)$. Results for $S({\bf q},\omega)$ 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 ${\bf m}_{\lambda}$ which appear in the theory of the spin liquid $\text{SL}_{\perp}$ [Section III] are given in Table 1.

Here we give the definitions in terms of the spins written in the local coordinate frame ${\sf S}_{i}$ (defined in Appendix A), cf. Ref. Yan et al., 2017 where definitions are given in the global, crystal basis.

Appendix E Numerical simulation of the correlations of the flux

Values of the flux field ${\bf B}_{\mu}(\mathbf{r})$ are calculated for each tetrahedron ${\bf r}$ according to Eq. (8) and the definitions of ${\bf m}_{\lambda}$ given in Table 1.

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 ${\bf S}_{{\bf B}_{\mu}}^{\alpha\beta}$ we use Eq. (12) where ${\bf B}_{\mu}(\mathbf{q})$ is defined as the lattice Fourier transform of ${\bf B}_{\mu}(\mathbf{r})$ over only the A sublattice of tetrahedra.

[TABLE]

where $N_{\sf 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_{\sf zz}$ and $T=0.1\ J_{\sf zz}$. Thermalisation was accomplished through a process of simulated annealing, with 104 Monte Carlo steps (MCs) of annealing from high temperature to temperature $T$, followed by 104 MCs of thermalization at temperature $T$, and 105 MCs of measurements at temperature $T$. Spin configurations were sampled every 100 MCs during the measurements, giving an ensemble of 1000 samples.

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 dynamical correlation function.

[TABLE]

where fluctuations of spin–nematic order are given by

[TABLE]

and the order parameter ${\bf Q}^{\sf site}_{\perp}({\bf r}_{i},t)$ is defined through Eq. (23).

$\chi_{Q_{\perp}^{\sf site}}(q,\omega)$ 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,

[TABLE]

where

[TABLE]

is the effective exchange field acting on site $i$, ${\sf J}_{ij}$ is the anisotropic exchange interaction tensor and the sum in Eq. (32) runs over the neighbors of $i$. The numerical integration of this nonlinear equation of motion proceeds as described in Ref. [Taillefumier et al., 2014].

In Fig 8 we plot the dispersion of the Goldstone mode found in molecular dynamics (MD) simulations of $\mathcal{H}_{\sf XXZ}$ within the spin–nematic phase for $J_{\pm}/J_{\sf zz}=-1$. The dispersion was extracted from the position of low–energy dispersing peak in $\chi_{Q_{\perp}}({\bf q},\omega)$, as shown in the inset to Fig. 5b.