# Theory of multiple-phase competition in pyrochlore magnets with anisotropic exchange with application to Yb2Ti2O7, Er2Ti2O7, and Er2Sn2O7

Han Yan, Owen Benton, L D C Jaubert, Nic Shannon

► To cite this version:

Han Yan, Owen Benton, L D C Jaubert, Nic Shannon. Theory of multiple-phase competition in pyrochlore magnets with anisotropic exchange with application to Yb2Ti2O7, Er2Ti2O7, and Er2Sn2O7. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2017, 95, pp.094422. 10.1103/PhysRevB.95.094422. hal-01542118

HAL Id: hal-01542118

https://hal.science/hal-01542118

Submitted on 19 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Theory of multiple-phase competition in pyrochlore magnets with anisotropic exchange, with application to Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> and Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>

Han Yan,<sup>1,2</sup> Owen Benton,<sup>1,3</sup> Ludovic Jaubert,<sup>1,4</sup> and Nic Shannon<sup>1,2,3</sup>

<sup>1</sup>Okinawa Institute of Science and Technology Graduate University, Onda-son, Okinawa 904-0395, Japan

<sup>2</sup>Clarendon Laboratory, University of Oxford, Parks Rd., Oxford OX1 3PU, UK

<sup>3</sup>H. H. Wills Physics Laboratory, University of Bristol, Tyndall Av, Bristol BS8-1TL, UK

<sup>4</sup>Rudolf Peierls Centre for Theoretical Physics, University of Oxford, 1-6 Keeble Rd, Oxford OX1 3NP, UK

(Dated: March 2, 2017)

The family of magnetic rare-earth pyrochlore oxides R<sub>2</sub>M<sub>2</sub>O<sub>7</sub> plays host to a diverse array of exotic phenomena, driven by the interplay between geometrical frustration and spin-orbit interaction, which leads to anisotropy in both magnetic moments and their interactions. In this article we establish a general, symmetry-based theory of pyrochlore magnets with anisotropic exchange interactions. Starting from a very general model of nearest-neighbour exchange between Kramers ions, we find four distinct classical ordered states, all with **q** = 0, competing with a variety of spin-liquids and unconventional forms of magnetic order. The finite-temperature phase diagram of this model is determined by Monte Carlo simulation, supported by classical spin-wave calculations. We pay particular attention to the region of parameter space relevant to the widely studied materials Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, and Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>. We find that many of the most interesting properties of these materials can be traced back to the “accidental” degeneracies where phases with different symmetries meet. These include the ordered ground state selection by fluctuations in Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, the dimensional-reduction observed in Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, and the lack of reported magnetic order in Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>. We also discuss the application of this theory to other pyrochlore oxides

PACS numbers: 74.20.Mn, 11.15.Ha, 75.10.Jm

## I. INTRODUCTION

Like high-energy physics, condensed matter is dominated by the idea of symmetry. Any physical property which *cannot* be traced back to a broken symmetry is therefore of enormous fundamental interest. In this context, the spin liquid phases found in frustrated magnets are a rich source of inspiration¹. Perhaps the most widely studied examples are the “spin ice” states in Ho₂Ti₂O₇ and Dy₂Ti₂O₇, classical spin-liquids famous for their magnetic monopole excitations². And there is now good reason to believe that a *quantum* spin-liquid phase, in which the magnetic monopoles are elevated to the role of “elementary” particles, could exist in spin-ice like materials where quantum effects play a larger role ³⁻¹⁴.

The extraordinary physics of spin ice stems from the combination of the geometrical frustration inherent to the pyrochlore lattice on which the magnetic rare earth ions R³⁺ reside, and the strongly anisotropic nature of the interactions between rare-earth ions¹⁵,¹⁶. This mixture of geometrical frustration and strong spin anisotropy is common to many pyrochlore materials, and gives rise to a wide array of interesting physical behaviors¹⁷.

The spin ices belong to a wider family of rare-earth pyrochlore oxides R<sub>2</sub>M<sub>2</sub>O<sub>7</sub> in which the magnetic ions have a doublet ground state, and highly-anisotropic interactions. The physical properties of these materials depend on the choice of rare-earth R³⁺ and transition metal M⁴⁺, and are fabulously diverse. In addition to spin ices, this family includes a wide range of systems that order magnetically, spin glasses and systems where local

<figure><img src="image_2.png" /><figcaption>FIG. 1: Classical ground-state phase diagram for a pyrochlore magnet with anisotropic exchange interactions. The model considered is the most general nearest-neighbour exchange Hamiltonian on the pyrochlore lattice $H_{ex}$ [Eq. (1)], with symmetric off-diagonal exchange $J_3 < 0$, and vanishing Dzyaloshinskii-Moriya interactions ($J_4 = 0$). There are four distinct ordered phases, illustrated in the insets of Fig. 3. Points correspond to published estimates of parameters for Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> [43], Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> [26], and Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub> [34], setting $J_4 = 0$. The white circle corresponds to the path through parameter space shown in Fig. 3.</figcaption></figure>

moments couple to itinerant electrons¹⁷⁻¹⁹. Materials of
<figure><img src="image_3.png" /><figcaption>FIG. 2: Suppression of classical order by quantum fluctuations in pyrochlore magnets with anisotropic exchange interactions, as described by $\mathcal{H}_{ex}$ [Eq. (1)], with $J_3 < 0$, $J_4 = 0$. Coloured regions show the four ordered phases illustrated in Fig. 3. White regions indicate where quantum fluctuations eliminate conventional magnetic order, within a linear spin-wave theory. Parameters $J_1/|J_3| = \rho \cos\theta$, $J_2/|J_3| = \rho \sin\theta$ are shown on a log-polar scale with $0 < \rho \lesssim 10^6$. The white circle corresponds to the path through parameter space shown in Fig. 3.</figcaption></figure>

current interest include Yb$_2$Ti$_2$O$_7$, which exhibits striking “rod-like” features in neutron scattering$^{20–23}$, and has been argued to undergo a Higgs transition into a ferromagnetically ordered state$^{24}$; Er$_2$Ti$_2$O$_7$, which appears to offer an elegant worked example of (quantum) order by disorder$^{25–29}$; and Er$_2$Sn$_2$O$_7$, which has yet to be seen to order at any temperature$^{30–34}$. Alongside continuing investigations into these materials, the last few years has seen the synthesis of a steady stream of new rare earth pyrochlore oxides, exhibiting both ordered$^{35–39}$ and disordered$^{53}$ low temperature states.

Given this “embarrassment of riches”, it seems reasonable to ask whether there is any common framework which can connect the properties of different rare-earth pyrochlore oxides, place new materials in context, and help guide the search for novel magnetic states. In this Article, we enlarge on the results in an earlier preprint$^{40}$ to develop a broad scenario for these materials, based on the concept of multiple-phase competition. We go on to show how this approach can be used to explain many of the interesting properties of Yb$_2$Ti$_2$O$_7$, Er$_2$Ti$_2$O$_7$ and Er$_2$Sn$_2$O$_7$.

Our starting point is the most general model of nearest-neighbour interactions compatible with the symmetries of the pyrochlore lattice$^{41–43}$

$$ \mathcal{H}_{\text{ex}} = \sum_{\langle ij \rangle} J_{ij}^{\mu\nu} S_i^{\mu} S_j^{\nu}, \quad (1) $$

where the sum on $\langle ij \rangle$ runs over the nearest-neighbour

bonds of the pyrochlore lattice, $S_i = (S_i^x, S_i^y, S_i^z)$ describes the magnetic moment of the rare-earth ion, and the matrix $J_{ij}^{\mu\nu}$ is a function of four independent parameters. Following the notation of Ross et al. [43], we identify these as “X-Y” ($J_1$), “Ising” ($J_2$), “symmetric off-diagonal” ($J_3$) and “Dzyaloshinskii-Moriya” ($J_4$) interactions. This model encompasses an extremely rich variety of different magnetic physics, including an exchange-based “spin-ice” ($J_1=-J_2=J_3=J_4 < 0$), and the Heisenberg antiferromagnet on a pyrochlore lattice ($J_1=J_2 > 0$, $J_3=J_4=0$), both of which are believed to have spin-liquid ground states$^{15,44}$. Nonetheless, materials such as Er$_2$Ti$_2$O$_7$, which is extremely well-described by a nearest-neighbour exchange model$^{25–27}$, do order magnetically$^{25}$.

The phase diagram of $\mathcal{H}_{ex}$ [Eq. (1)] for a quantum spin-1/2 has previously been studied using mean-field and spin-wave approximations, with many papers emphasising connections with spin ice$^{5,8,9,11,45}$. In this Article we take a different approach, starting from an analysis of the way in which different spin configurations break the point-group symmetries of the pyrochlore lattice. We show that, for classical spins, the problem of finding the ground state of $\mathcal{H}_{ex}$ [Eq. (1)] can be neatly separated into two steps: i) finding the ground state of a single tetrahedron and ii) understanding how the spin-configuration on that tetrahedron can be used to tile the pyrochlore lattice. The first step, in turn, reduces to understanding how the different interactions in the model transform under the symmetries $T_d$ of a single tetrahedron. The second step, summarized in a simple set of “Lego-brick” rules, enables us to encompass both ordered ground states, which break lattice symmetries, and spin-liquids, which do not.

This approach, augmented by spin-wave calculations and extensive classical Monte Carlo simulations, makes it possible both to establish a complete phase diagram for $\mathcal{H}_{ex}$ [Eq. (1)] as a function ($J_1, J_2, J_3, J_4$), and to link ground state properties to predictions for neutron-scattering experiments. In this article, taking our motivation from estimated parameters for Yb$_2$Ti$_2$O$_7$ [43], Er$_2$Ti$_2$O$_7$ [26] and Er$_2$Sn$_2$O$_7$ [34], we concentrate on ordered phases in the limit $J_3 < 0$, $J_4 = 0$. Here there is a competition between four different types of order: a Palmer-Chalker$^{46}$ phase ($\Psi_4$), a non-collinear ferromagnet (FM), a coplanar antiferromagnet ($\Psi_3$), and a non-coplanar antiferromagnet ($\Psi_2$). The way in which these phases relate to one another is illustrated in Fig. 1, Fig. 2 and Fig. 3.

Crucially, the same symmetry-based approach used to find ordered ground states also permits us to explore the way in which these physically distinct states are connected by the “accidental” degeneracies arising at boundaries between phases with different symmetry. The enlarged ground-state manifolds at these phase boundaries have far-reaching consequences, once quantum and thermal fluctuations are taken into account. The common theme which emerges is of systems “living on the edge”
<figure><img src="image_4.png" /><figcaption>FIG. 3: Finite-temperature phase diagram for a pyrochlore magnet with anisotropic exchange interactions. The model considered is $\mathcal{H}_{\text{ex}}$ [Eq. (1)], with $J_1 = 3|J_3|\cos\theta$, $J_2 = 3|J_3|\sin\theta$, $J_3 < 0$, and $J_4 \equiv 0$, corresponding to the white circle in Fig. 1. Points show finite temperature phase transitions found from classical Monte Carlo simulations. The four ordered phases, Palmer-Chalker ($\Psi_4$), non-collinear ferromagnetic (FM), coplanar antiferromagnetic ($\Psi_3$) and non-coplanar antiferromagnetic ($\Psi_2$), are illustrated at the top of the figure. Each of these phases is six-fold degenerate, with zero crystal momentum, and is completely specified by the spin configuration in a single tetrahedron.</figcaption></figure>

— the physical properties of materials showing one type of magnetic order being dictated by the proximity of another, competing, ordered phase.

Thus, in $Yb_2Ti_2O_7$, we find ferromagnetic order proximate to competing, “$\Psi_3$” and “$\Psi_2$” phases, which manifest themselves in the “rods” seen in neutron scattering. Meanwhile, in $Er_2Ti_2O_7$, we discover that the reason fluctuations select the well-established “$\Psi_2$” ground state$^{25-27}$, is proximity to a neighboring Palmer-Chalker phase, as illustrated in Fig. 4. And in the case of $Er_2Sn_2O_7$, we find that fluctuations of Palmer-Chalker order predominate, but that all forms of magnetic order are strongly suppressed by the proximity of a degenerate ground-state manifold connected to a neighbouring “$\Psi_2$” phase.

We note that the same approach of combining symmetry analysis and the “Lego-brick” rules can also be used to systematically search for unconventional ordered states and new (classical) spin-liquid phases on the pyrochlore lattice. This is a theme which will be developed elsewhere$^{47,48}$. The remainder of the present article is structured as follows:

In Section II we introduce a general model of nearest-neighbour exchange interactions on a pyrochlore lattice and, restricting to classical spins, establish the conditions under which the model has a magnetically ordered ground state. We also provide a complete classification of possible ordered states in terms of the irreducible representations of the tetrahedral symmetry group $T_d$.

In Section III we show that this symmetry analysis can be used to determine the classical ground state of $\mathcal{H}_{\text{ex}}$ [Eq. (1)] for arbitrary parameters ($J_1, J_2, J_3, J_4$). The nature of the ground states in the limit ($J_3 < 0, J_4 = 0$), which is of particular relevance to real materials, is explored in some detail, including analysis of the degenerate manifolds arising at the phase boundaries of the

model.

In Section IV we explore the spin-wave excitations associated with these ordered phases. This enables us to make predictions for neutron scattering, and to develop a ground state phase diagrams for classical and semiclassical spins, focusing again on the limit $J_3 < 0$ and $J_4 = 0$, as illustrated in Fig. 1. It also enables us to identify regions of the phase diagram where strong quantum fluctuations are liable to eliminate classical order entirely, as illustrated in Fig. 2.

In Section V we use classical Monte Carlo simulation to explore the finite-temperature phase transitions which separate each of the ordered phases from the high-temperature paramagnet. The results of this analysis are summarized in Fig. 3.

In Section VI we study the finite-temperature consequences of the enlarged ground-state manifolds arising at the boundary between different ordered phases. This is illustrated in Fig. 4.

In Sections VII, VIII and IX, we discuss the implications of these results for the rare-earth pyrochlore oxides $Er_2Ti_2O_7$, $Yb_2Ti_2O_7$ and $Er_2Sn_2O_7$, respectively. Other rare-earth pyrochlore magnets to which the theory might apply are discussed briefly in Section X.

We conclude in Section XI with a summary of our results, and an overview of some of the interesting open issues.

Technical details of calculations are reproduced in a short series of Appendices at the end of the Article: Appendix A provides details of the local coordinate frame throughout the Article, and the associated form of the $g$-tensor. Appendix B provides technical details of the linear spin-wave calculations described in Section IV. Appendix C provides technical details associated with the classical Monte Carlo simulations described in Section V.

As far as possible, we have endeavoured to make Sec-
<figure><img src="image_5.png" /><figcaption>FIG. 4: Selection of an ordered ground state by thermal fluctuations in Er$_{2}$Ti$_{2}$O$_{7}$. For $T < T_c$, fluctuations select six states with non-coplanar antiferromagnetic order ($\Psi_2$) from a one-dimensional manifold of degenerate ground states. The entropic selection of these six states can be traced to an enlarged ground-state manifold found on the boundary with the Palmer-Chalker phase ($\Psi_4$). Plots show the probability distribution of the order parameter $m_E = (m_{E_x}, m_{E_y})$, as described in Section VI.</figcaption></figure>

tions VII–X, describing the application of the theory to experiments on rare-earth pyrochlores, self-contained. Readers chiefly interested in these materials may safely omit the theoretical development in Sections III to VI of the Article.

## II. MICROSCOPIC MODEL OF ANISOTROPIC EXCHANGE

### A. Magnetism at the level of a single ion

Pyrochlore oxides, A$_{2}$B$_{2}$O$_{7}$, are a ubiquitous feature of igneous rocks throughout the world. This broad family of materials takes its name from the mineral “pyrochlore” [(Ca, Na)$_{2}$Nb$_{2}$O$_{6}$(OH,F)], which burns with a green ($\chi\lambda\omega\rho\delta\zeta$) fire ($\pi\nu\rho$), and shares its crystal structure with a great many other oxides, halides and chalcogenides. Here we concentrate on those pyrochlore oxides in which the B-cation is a non-magnetic transition metal, such as Ti$^{4+}$ or Sn$^{4+}$, while the remaining cation A$^{3+}$ is a rare-earth ion with a magnetic doublet ground state. These magnetic ions form a pyrochlore lattice, built of corner-sharing tetrahedra, which shares the same cubic

symmetry $Fd\bar{3}m$ as the parent material.

Even within this restricted group of rare-earth oxides, the interplay between strong spin-orbit coupling, and the crystal electric field (CEF) at the A-cation site, leads to a huge variation in the magnetic properties of the rare-earth ion. For example, Dy$^{3+}$ provides the strong Ising moment in the spin-ice Dy$_{2}$Ti$_{2}$O$_{7}$, while Er$^{3+}$ forms a moment with XY-like character in Er$_{2}$Ti$_{2}$O$_{7}$ [49].

The goal of this article is not to explore the intricate CEF ground states of rare-earth ions (see e.g. Ref. [50,51] for a discussion on this topic), but rather to understand the way in which the anisotropic exchange interactions between them shape the magnetism of rare-earth pyrochlore oxides. We therefore concentrate on materials in which the ground state of the rare-earth ion is a Kramers doublet, with an odd number of electrons, like Yb$^{3+}$ ([Xe]4f$^{13}$) or Er$^{3+}$ ([Xe]4f$^{11}$).

In this case, as long as the temperature is small compared with the lowest-lying CEF excitation, the magnetic ion can be described by a pseudospin-1/2 degree of freedom

$$[S^\mu, S^\nu] = i\epsilon_{\mu\nu\xi}S^\xi. \qquad (2)$$

It is important to note that, even with the restriction to Kramers doublets, there is more than one possibility for how $S^\mu$ will transform under space-group operations$^{52}$. In this article we will focus on the case where $S^\mu$ transforms like a magnetic dipole, which is the case appropriate to Yb$^{3+}$ and Er$^{3+}$ based pyrochlores. We note that an alternative “dipolar-octupolar” case may be realized in Dy$^{3+}$, Nd$^{3+}$ and Ce$^{3+}$ based pyrochlores$^{52–54}$.

Where $S^\mu$ transforms like a magnetic dipole, it will be associated with an effective magnetic moment

$$m_i^\mu = \sum_{\nu=1}^{3} g_i^{\mu\nu} S_\nu \qquad (3)$$

where $\mu, \nu = \{x, y, z\}$. Since the magnetic anisotropy of rare-earth ion is determined by the local CEF, the tensor $g_i^{\mu\nu}$ is site-dependent, as described in Appendix A. This has important consequences for the magnetic correlations measured in neutron scattering experiments, discussed below.

### B. Anisotropy in exchange interactions

The interplay between spin-orbit coupling and CEF leads to anisotropy in the interactions between rare-earth ions, just as it leads to anisotropy in the magnetic ground state of an individual ion$^{50}$. It is possible to make estimates of exchange interactions in a pyrochlore oxides from knowledge of the CEF ground state and low-lying excitations$^{51,55}$. However for the purposes of this article it is sufficient to consider the constraints on these interactions imposed by the symmetry of the lattice.
In the case of Kramers ions on a pyrochlore lattice, the most general form of nearest-neighbour exchange can be broken down into a sum over tetrahedra $t$

$$ \mathcal{H}_{\text{ex}} = \sum_{\langle ij \rangle} J_{ij}^{\mu\nu} S_i^{\mu} S_j^{\nu} = \sum_t \mathcal{H}_{\text{ex}}^{\text{tet}}[t], \quad (4) $$

where

$$ \mathcal{H}_{\text{ex}}^{\text{tet}}[t] = \sum_{i,j \in t} \mathbf{S}_i \mathbf{J}_{ij}^{[t]} \mathbf{S}_j. \quad (5) $$

Here, $\mathbf{S}_i = (S_i^x, S_i^y, S_i^z)$, and $\mathbf{J}_{ij}^{[t]}$ is a $3 \times 3$ matrix specific to the bond $ij$, within tetrahedron $t$. The exchange interactions $\mathbf{J}_{ij}$ do not, in general, possess any continuous spin-rotation invariance. Nonetheless, the form of exchange $\mathbf{J}_{ij}$ is strongly constrained by the symmetry of the bond $ij$, and the interactions on different bonds must also be related by lattice symmetries.

Once these constraints are taken into account$^{41}$, $\mathbf{J}_{ij}$ is a function of just four independent parameters and, for the six bonds which make up the tetrahedron shown in Fig. 5, can be written

$$ \begin{aligned} \mathbf{J}_{01} &= \begin{pmatrix} J_2 & J_4 & J_4 \\ -J_4 & J_1 & J_3 \\ -J_4 & J_3 & J_1 \end{pmatrix} & \mathbf{J}_{02} &= \begin{pmatrix} J_1 & -J_4 & J_3 \\ J_4 & J_2 & J_4 \\ J_3 & -J_4 & J_1 \end{pmatrix} \\ \mathbf{J}_{03} &= \begin{pmatrix} J_1 & J_3 & -J_4 \\ J_3 & J_1 & -J_4 \\ J_4 & J_4 & J_2 \end{pmatrix} & \mathbf{J}_{12} &= \begin{pmatrix} J_1 & -J_3 & J_4 \\ -J_3 & J_1 & -J_4 \\ -J_4 & J_4 & J_2 \end{pmatrix} \\ \mathbf{J}_{13} &= \begin{pmatrix} J_1 & J_4 & -J_3 \\ -J_4 & J_2 & J_4 \\ -J_3 & -J_4 & J_1 \end{pmatrix} & \mathbf{J}_{23} &= \begin{pmatrix} J_2 & -J_4 & J_4 \\ J_4 & J_1 & -J_3 \\ -J_4 & -J_3 & J_1 \end{pmatrix} \end{aligned} \quad (6) $$

where we label lattice sites and interactions following the conventions of Ross et al. [43]. The structure of these matrices imply that the different contributions to the interaction $\mathbf{J}_{ij}$ [Eq. (6)] can be approximately identified as

* $J_1 \to \text{"XY"}$ with respect to the local bond

* $J_2 \to \text{"Ising"}$ with respect to the local bond

* $J_3 \to \text{symmetric off diagonal exchange}$

* $J_4 \to \text{Dzyaloshinskii-Moriya}$

The anisotropic nearest-neighbour exchange model, $\mathcal{H}_{\text{ex}}$ [Eq. (1)], has been applied with considerable success to a number of pyrochlore oxides. In the case of Yb$_2$Ti$_2$O$_7$, $\mathcal{H}_{\text{ex}}$ [Eq. (1)] has been shown to give an excellent description of spin-wave spectra measured in magnetic field$^{43}$. Thermodynamic quantities, calculated from $\mathcal{H}_{\text{ex}}$ [Eq. (1)] using the parameters from [43], also gave a very good description of experiments$^{56,57}$. Parameters for Er$_2$Ti$_2$O$_7$ have been extracted from equivalent inelastic neutron scattering experiments$^{26}$, and from measurements of the field-dependence of magnetisation at low

<figure><img src="image_6.png" /><figcaption>FIG. 5: A single tetrahedron within the pyrochlore lattice, showing the convention used in labelling sites. The positions of the magnetic sites relative to the centre of the tetrahedron are defined in Appendix A.</figcaption></figure>

<table><caption>TABLE I: Estimates of the parameters for anisotropic near-neighbour exchange, taken from experiments on Yb$_2$Ti$_2$O$_7$ [43], Er$_2$Ti$_2$O$_7$ [26], and Er$_2$Sn$_2$O$_7$ [34]. Values are quoted for exchange interactions in both the crystal coordinate frame $\mathcal{H}_{\text{ex}}$ [Eq. (1)], and the local coordinate frame, $\mathcal{H}_{\text{ex}}^{\text{local}}$ [Eq. (8)], following the notation of Ross et al. [43]. An alternative set of parameters for Yb$_2$Ti$_2$O$_7$ has recently been proposed by Robert et al. [71].</caption><thead><tr><th></th><th>Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> [43]</th><th>Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> [26]</th><th>Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub> [34]</th></tr></thead><tbody><tr><td>J<sub>1</sub></td><td>-0.09 meV</td><td>0.11 meV</td><td>0.07 meV</td></tr><tr><td>J<sub>2</sub></td><td>-0.22 meV</td><td>-0.06 meV</td><td>0.08 meV</td></tr><tr><td>J<sub>3</sub></td><td>-0.29 meV</td><td>-0.10 meV</td><td>-0.11 meV</td></tr><tr><td>J<sub>4</sub></td><td>0.01 meV</td><td>-0.003 meV</td><td>0.04 meV</td></tr><tr><td>J<sub>zz</sub></td><td>0.17 meV</td><td>-0.025 meV</td><td>0</td></tr><tr><td>J<sub>&pm;</sub></td><td>0.05 meV</td><td>0.065 meV</td><td>0.014 meV</td></tr><tr><td>J<sub>&pm;&pm;</sub></td><td>0.05 meV</td><td>0.042 meV</td><td>0.074 meV</td></tr><tr><td>J<sub>z&pm;</sub></td><td>-0.14 meV</td><td>-0.009 meV</td><td>0</td></tr></tbody></table>

temperature$^{58}$. The model $\mathcal{H}_{\text{ex}}$ [Eq. (1)], using parameters taken from neutron scattering$^{26}$, has been shown to give good agreement with the observed spin wave spectrum in Er$_2$Ti$_2$O$_7$, consistent with quantum order-by-disorder [26]. Anisotropic nearest-neighbour exchange parameters for Er$_2$Sn$_2$O$_7$ have also been estimated from measurements of the magnetization curve$^{34}$.

Representative estimates of the exchange parameters ($J_1, J_2, J_3, J_4$) taken from experiment on Yb$_2$Ti$_2$O$_7$ [43], Er$_2$Ti$_2$O$_7$ [26], and Er$_2$Sn$_2$O$_7$ [34], are shown in Table I. The typical scale of interactions is $|J| \sim 0.1$ meV (i.e. $|J| \sim 1$ K), with typical uncertainty in estimates of order $\delta J \sim 0.02$ meV (i.e. $\delta J \sim 0.2$ K) [26,34,43]. In all of these cases, the symmetric off-diagonal exchange interaction $J_3$ is negative, while the Dzyaloshinskii-Moriya interaction $J_4$ is relatively small$^{154}$.
### C. Anisotropic exchange in a local frame

Since both the anisotropy in magnetic ground state of the rare-earth ion, and the anisotropy in its interactions, are dictated by the local CEF field, it is often convenient to describe them in a local coordinate frame

$$ \{\mathbf{x}_i^{\text{local}}, \mathbf{y}_i^{\text{local}}, \mathbf{z}_i^{\text{local}}\} $$

such that the axis $\mathbf{z}_i^{\text{local}}$ aligns with the $C_3$ symmetry axis of the local CEF on site $i$, as described in Appendix A. We introduce a SU(2) (pseudo) spin-1/2 in this local frame

$$ [S_i^\alpha, S_i^\beta] = i\epsilon_{\alpha\beta\gamma}S_i^\gamma, \quad (7) $$

where $\alpha, \beta, \gamma = \{\mathbf{x}_i^{\text{local}}, \mathbf{y}_i^{\text{local}}, \mathbf{z}_i^{\text{local}}\}$. Note that throughout the manuscript $S_i^\alpha$ will refer to the spin components in this local frame, while $S_i^\alpha$ refers to the spin components in the global, crystal, coordinate system.

In the local coordinate frame, the most general form of exchange interactions between Kramers ions on the pyrochlore lattice can be written$^{43}$

$$ \mathcal{H}_{\text{ex}}^{\text{local}} = \sum_{\langle ij \rangle} \left\{ \begin{aligned} & J_{zz} S_i^z S_j^z - J_{\pm} (S_i^+ S_j^- + S_i^- S_j^+) \\ & + J_{\pm\pm} [\gamma_{ij} S_i^+ S_j^+ + \gamma_{ij}^* S_i^- S_j^-] \\ & + J_{z\pm} [S_i^z (\zeta_{ij} S_j^+ + \zeta_{ij}^* S_j^-) + i \leftrightarrow j] \end{aligned} \right\} \quad (8) $$

where the matrix

$$ \zeta = \begin{pmatrix} 0 & -1 & e^{i \frac{\pi}{3}} & e^{-i \frac{\pi}{3}} \\ -1 & 0 & e^{-i \frac{\pi}{3}} & e^{i \frac{\pi}{3}} \\ e^{i \frac{\pi}{3}} & e^{-i \frac{\pi}{3}} & 0 & -1 \\ e^{-i \frac{\pi}{3}} & e^{i \frac{\pi}{3}} & -1 & 0 \end{pmatrix} \quad \gamma = -\zeta^* \quad (9) $$

encode the change in coordinate frame between different sublattices.

The relationship between the parameters in this local frame, ($J_{zz}$, $J_\pm$, $J_{\pm\pm}$, $J_{z\pm}$), and exchange parameters in the global frame of the crystal axes ($J_1$, $J_2$, $J_3$, $J_4$), is given in Table II. Corresponding estimated parameters from experiment on Yb$_2$Ti$_2$O$_7$ [43], Er$_2$Ti$_2$O$_7$ [26], and Er$_2$Sn$_2$O$_7$ [34], are shown in Table I.

### D. Proof of the existence of a classical ground state with $q=0$, 4-sublattice order

Finding the ground state of $\mathcal{H}_{\text{ex}}$ [Eq. (1)], for a quantum (pseudo)spin-1/2, and arbitrary exchange interactions ($J_1$, $J_2$, $J_3$, $J_4$), is a very difficult problem, in general only tractable as a mean-field theory$^{5,8,9,11,45}$. However, many rare-earth pyrochlores are known to have relatively simple ground states, with vanishing crystal momentum $q=0$, implying a 4-sublattice magnetic order$^{17}$. Here we show that, under the restriction that $\mathbf{S}_i$ is a classical variable, $\mathcal{H}_{\text{ex}}$ [Eq. (4)] always possesses a ground state of this

<table><caption>TABLE II: Relationship between the parameters of the anisotropic nearest-neighbour exchange model in the local coordinate frame, $\mathcal{H}_{\text{ex}}^{\text{local}}$ [Eq. (8)], and the exchange parameters in the crystal coordinate frame $\mathcal{H}_{\text{ex}}$ [Eq. (1)]. The notation used for the different components of the interaction follows Ross et al. [43].</caption><thead><tr><th>interaction in local<br>coordinate frame</th><th>exchange parameters<br>in global frame</th></tr></thead><tbody><tr><td>$J_{zz}$</td><td>$-\frac{1}{3}(2J_1 - J_2 + 2J_3 + 4J_4)$</td></tr><tr><td>$J_\pm$</td><td>$\frac{1}{6}(2J_1 - J_2 - J_3 - 2J_4)$</td></tr><tr><td>$J_{\pm\pm}$</td><td>$\frac{1}{6}(J_1 + J_2 - 2J_3 + 2J_4)$</td></tr><tr><td>$J_{z\pm}$</td><td>$\frac{1}{3\sqrt{2}}(J_1 + J_2 + J_3 - J_4)$</td></tr></tbody></table>

type. In Section II E, below, we explore the conditions under which this classical ground state is unique.

We begin with the simple observation that, since $\mathcal{H}_{\text{ex}}$ [Eq. (1)] is expressed as a sum over individual tetrahedra, any state which minimizes the energy of each individual tetrahedron must be a ground state. It is convenient to split this sum into two pieces

$$ \mathcal{H}_{\text{ex}} = \sum_{t \in A} \mathcal{H}_{\text{ex}}^{\text{A}}[t] + \sum_{t' \in B} \mathcal{H}_{\text{ex}}^{\text{B}}[t']; \quad (10) $$

where A and B refer to the two distinct sublattices of tetrahedra, with

$$ \mathcal{H}_{\text{ex}}^{\text{A}}[t] = \sum_{i,j \in t} \mathbf{S}_i \mathbf{J}_{ij}^{[\text{A}]} \mathbf{S}_j \quad (11) $$

$$ \mathcal{H}_{\text{ex}}^{\text{B}}[t'] = \sum_{i,j \in t'} \mathbf{S}_i \mathbf{J}_{ij}^{[\text{B}]} \mathbf{S}_j . \quad (12) $$

The interactions $\mathbf{J}_{ij}^{[\text{A}]}$ and $\mathbf{J}_{ij}^{[\text{B}]}$ are related by inversion about a single site $\mathcal{I}$

$$ \mathbf{J}_{ij}^{[\text{B}]} = \mathcal{I} \cdot \mathbf{J}_{ij}^{[\text{A}]} \cdot \mathcal{I}. \quad (13) $$

Since $\mathcal{I}^2 = 1$, we can write

$$ \begin{aligned} \mathbf{S}_i \cdot \mathbf{J}_{ij}^{[\text{A}]} \cdot \mathbf{S}_j &= \mathbf{S}_i \cdot \mathcal{I}^2 \cdot \mathbf{J}_{ij}^{[\text{A}]} \cdot \mathcal{I}^2 \cdot \mathbf{S}_j \\ &= \mathbf{S}_i \cdot \mathcal{I} \cdot \mathbf{J}_{ij}^{[\text{A}]} \cdot \mathcal{I} \cdot \mathbf{S}_j \end{aligned} \quad (14) $$

where we have used the fact that the spin $\mathbf{S}_i$ is invariant under lattice inversion. This implies

$$ \mathbf{J}_{ij}^{[\text{A}]} = \mathcal{I} \cdot \mathbf{J}_{ij}^{[\text{A}]} \cdot \mathcal{I} = \mathbf{J}_{ij}^{[\text{B}]} . \quad (15) $$

It follows that interactions for any tetrahedron $t$ is the same, regardless of which tetrahedral sublattice it belongs to, and we can safely write

$$ \mathcal{H}_{\text{ex}}^{\text{tet}} = \sum_{i,j \in t} \mathbf{S}_i \mathbf{J}_{ij} \mathbf{S}_j \quad (16) $$

where $\mathbf{J}_{ij}$ are given by Eq. (6) and the sum runs over all pairs of sites $i, j$ in a given tetrahedron $t$, which may now be of either sublattice.
The proof we are seeking follows directly from this result [Eq. (16)] : for classical spins, [$\mathcal{H}_{\text{ex}}^{\text{A}}, \mathcal{H}_{\text{ex}}^{\text{B}}$] = 0, and we can construct a ground state of $\mathcal{H}_{\text{ex}}$ by choosing any state which minimizes the energy of a single tetrahedron, and repeating it across all A-sublattice (or B-sublattice) tetrahedra. Since every spin is shared between one A- and one B-sublattice tetrahedron, and the Hamiltonians for A- or B-sublattices are equivalent [Eq. (16)], any such classical spin-configuration which minimizes the energy on one tetrahedral sublattice, simultaneously minimizes the energy on the other tetrahedral sublattice, and is a ground state of $\mathcal{H}_{\text{ex}}$ [Eq. (1)].

It follows that there *always* exists a classical, $\mathbf{q} = 0$ ground state of $\mathcal{H}_{\text{ex}}$ [Eq. (1)], with 4-sublattice long-range magnetic order, for arbitrary exchange interactions ($J_1, J_2, J_3, J_4$). This is true *even* in the presence of finite Dzyaloshinskii-Moriya interaction $J_4$.

Such a $\mathbf{q} = 0$ ground state has a finite, discrete degeneracy associated with the breaking of point-group and time-reversal symmetries (in the case of classical spins, time-reversal corresponds to the inversion of all spins $\mathbf{S}_i \rightarrow -\mathbf{S}_i$). This degeneracy must be at least 2 (time-reversal), and is typically 6 ($C_3$ rotations $\otimes$ time-reversal), for the ordered phases considered in this article.

### E. Conditions for the uniqueness of 4-sublattice order — the “Lego-brick” rules

The existence of a classical ground state of $\mathcal{H}_{\text{ex}}$ [Eq. (1)] with 4-sublattice magnetic order, for arbitrary exchange interactions ($J_1, J_2, J_3, J_4$) constitutes a enormous simplification, since it is much easier to determine the spin-configuration which minimizes the energy of a single tetrahedron (as described in Section III, below) than to find the ground state of the entire lattice. However, as we shall see, many of interesting properties of rare-earth pyrochlores follow from the fact that while such a classical ground state must exist, it need not be unique.

Establishing the uniqueness of a 4-sublattice ground state, up to the discrete degeneracy of the state itself, amounts to determining the number of ways in which the spin configurations which minimize the energy of a single tetrahedron can be used to tile the entire lattice. For many purposes, it is convenient to think of these as a set of “Lego-brick” rules for fitting together spin-configurations on a lattice. These rules allow us to determine the degeneracy, and nature, of the ground states of the whole lattice, using the ground states of a single tetrahedron.

The rules can be stated as follows:

1. If the spin on every site of the tetrahedron points in a different direction, in each of its classical ground-states, then the 4-sublattice ground state of the lattice is unique (up to global symmetry operations). In this case, the degeneracy of the ground state of the lattice is the same as that of a single tetrahedron.

2. If, within the set of ground states for a single tetrahedron, there are two states in which a single site has the same spin orientation, the 4-sublattice ground state of the lattice is not unique. In this case, the system undergoes a dimensional reduction into independent kagome planes, and the number of classical ground states is at least $O(2^L)$, where $L$ is the linear size of the system.

3. If, within the set of ground states for a single tetrahedron, there is a pair of states which have the same spin orientation on two sites, the 4-sublattice ground state is also not unique. However in this case, the number of classical ground states must grow as *at least* $O(2^{L^2})$, corresponding to dimensional reduction into independent chains of spins. In the special case of spin-ice, the corresponding classical ground-state degeneracy is extensive.

The first rule guarantees the uniqueness of a 4-sublattice ground state, where the spin on every site of the tetrahedron points in a different direction in each of the ground states of a single tetrahedron. Away from phase boundaries, this is true for all of the 4-sublattice $\mathbf{q} = 0$ ordered phases discussed in this article. However, it is clear from Rules 2 and 3 that if two of the ground states of a single tetrahedron share a common spin — i.e. the spin on a given site points in the same direction in more than one ground state — then it is *always* possible to construct other ground states, with finite $\mathbf{q}$.

To give a concrete example of how these “Lego-brick” rules work, let us assume that two different ground states for a single tetrahedron have identical orientation of the spin on site 0, but different orientation of the spins on sites 1, 2 and 3. In this case it is possible to divide the pyrochlore lattice into a set of parallel kagome planes, containing spins associated with sites 1,2 and 3 of a tetrahedron, separated by triangular-lattice planes associated with site 0. Since each successive kagome plane can take on one of two different spin configurations, the number of such ground states grows as $2^{N_K}$, where $N_K$ is the number of kagome planes, and encompasses all possible $\mathbf{q} \parallel [111]$. Dimensional reduction of this type occurs on the classical phase boundary between ordered FM and Palmer-Chalker phases discussed in Section III of this article.

An example where Rule 3 applies, and a set of independent chains emerges in the ground state manifold, is the phase boundary between the Palmer-Chalker phase and the non-coplanar antiferromagnet discussed in Section III G. However the “Lego-brick” rules permit even larger ground-state degeneracies, as is known from the “two-in, two-out” states, made famous by the spin-ice problem. In this case there are a total of 6 possible ground states for a single tetrahedron, but each possible spin orientation, on each site, belongs to 3 different ground states. According to Rule 3, the 4-sublattice classical ground-state — a ferromagnet — should not be unique, and the total number of classical ground states must grow as *at least* $O(2^{L^2})$. In fact, the
<figure><img src="image_7.png" /><figcaption>FIG. 6: The “Lego-brick” rules describing how the ground states of a single tetrahedron can be connected to tile the pyrochlore lattice. The two tetrahedra in the left and right panels represent a pair of tetrahedra in distinct ground-state configurations. Distinct spin orientations on the sites of each tetrahedron are denoted by letters **a**–**h**. Three cases are shown. In case (1) all of the ground states for a single tetrahedron have different spin orientations for any given site. This means that two tetrahedra in distinct ground states cannot be joined together because they do not share a common spin orientation on any site. In this case the 4-sublattice ground state of the lattice is unique with $q = 0$ order, up to global symmetry operations. In case (2) there are two ground states configurations for a tetrahedron which share a common spin orientation on a single site, here denoted **a**. These tetrahedra can be joined together by sharing the spin in orientation **a**. In this case the ground state of the lattice has a degeneracy of at least $O(2^L)$. Indeed, successive kagome layers of spins can be independently in BDC or FEG configurations. In case (3) there is a pair of ground states which share common spin orientations on two sites, here denoted **a** and **b**. These tetrahedra can be joined together by sharing the spin in orientation **a** or the spin in orientation **b**. In this case the ground state of the lattice has a degeneracy of at least $O(2^{L^2})$.</figcaption></figure>

ground-state degeneracy of spin ice is extensive, scaling as $\Omega_{\text{ice}} \sim (3/2)^{N/2}$, where N is the total number of sites in the lattice$^{61}$. This manifold of spin-ice states includes ground states with all possible **q**.

## F. Representation theory

Except in very specific limits, such as the Heisenberg model ($J_1 = J_2 = J, J_3 = J_4 = 0$), Hamiltonian $\mathcal{H}_{\text{ex}}$ [Eq. (1)] does not possess any continuous spin-

rotation symmetry. The key to unlocking its properties, therefore, is to understand how different ordered states, and indeed different spin-fluctuations, break the space-group symmetries of the pyrochlore lattice. This task is made easier by the realization that a classical ground state with $\mathbf{q} = 0$ exists for all possible ($J_1, J_2, J_3, J_4$), as discussed in Section IID. It is therefore possible to restrict discussion to the point-group symmetries of the lattice.

In what follows, we explore the consequences of applying representation theory for these point-group operations to a general model of anisotropic nearest-neighbor exchange on the pyrochlore lattice, $\mathcal{H}_{\text{ex}}$ [Eq. (1)]. This analysis serves a two-fold purpose: it reduces the Hamiltonian for a single tetrahedron $\mathcal{H}_{\text{ex}}^{\text{tet}}$ [Eq. (16)] to a diagonal form, and it provides a set of order parameters with which to characterise the $\mathbf{q} = 0$, 4-sublattice ordered phases found in real materials.

The point-group symmetry of the pyrochlore lattice is the cubic symmetry group $O_h = T_d \times I$. Here $T_d$ is symmetry group of a single tetrahedron, and $I = \{\epsilon, \mathcal{I}\}$, where $\epsilon$ is the identity and $\mathcal{I}$ corresponds to the lattice inversion introduced in Eq. (13). For classical spins, lattice inversion plays a benign role (cf. Section IID), and it is sufficient to consider $T_d$ alone. The group $T_d$ has 24 elements$^{62}$, corresponding to the symmetries of the tetrahedron: $8 \times C_3 - \frac{2\pi}{3}$ rotation around a [111] axis; $3 \times C_2 - \pi$ rotation around [100] axis; $6 \times S_4 - \frac{\pi}{2}$ rotation around a [100] axis followed by reflection in the same [100] plane; $6 \times \sigma_d$ — reflection in [011] plane; and $\epsilon$ — the identity$^{63}$.

The different ways in which classical ground states with $\mathbf{q} = 0$ break the symmetries of a tetrahedron can be fully characterised by introducing order parameters $\mathbf{m}_\lambda$ which transform with the non-trivial irreducible representations $\lambda = \{A_2, E, T_1, T_2\}$ of $T_d$. These order parameters are formed by linear combinations of spin components, and can be expressed in either global coordinate frame of $\mathcal{H}_{\text{ex}}$ [Eq. (1)] — cf. Table III — or in the local coordinate frame of $\mathcal{H}_{\text{ex}}^{\text{local}}$ [Eq. (8)] — cf. Table IV.

The anisotropic exchange Hamiltonian $\mathcal{H}_{\text{ex}}^{\text{tet}}$ [Eq. (16)] can be transcribed exactly in terms of same set of irreps.

$$ \mathcal{H}_{\text{ex}}^{\text{tet}} \equiv \frac{1}{2} \left[ a_{A_2} m_{A_2}^2 + a_E \mathbf{m}_E^2 + a_{T_2} \mathbf{m}_{T_2}^2 + a_{T_{1,A}} \mathbf{m}_{T_{1,A}}^2 + a_{T_{1,B}} \mathbf{m}_{T_{1,B}}^2 + a_{T_{1,AB}} \mathbf{m}_{T_{1,A}} \cdot \mathbf{m}_{T_{1,B}} \right], \quad (17) $$

where the coefficients

$$ \begin{aligned} a_{A_2} &= -2J_1 + J_2 - 2(J_3 + 2J_4) \\ a_E &= -2J_1 + J_2 + J_3 + 2J_4 \\ a_{T_2} &= -J_2 + J_3 - 2J_4 \\ a_{T_{1,A}} &= 2J_1 + J_2 \\ a_{T_{1,B}} &= -J_2 - J_3 + 2J_4 \\ a_{T_{1,AB}} &= -\sqrt{8}J_3 \end{aligned} \quad (18) $$

are completely determined by the parameters of $\mathcal{H}_{\text{ex}}$ [Eq. (1)]. Equivalent expressions for $a_\lambda$ can be found
<table><caption>TABLE III: Order parameters $\mathbf{m}_\lambda$, describing how the point-group symmetry of a single tetrahedron within the pyrochlore lattice is broken by magnetic order. Order parameters transform according to irreducible representations of the point-group $T_d$, and are expressed in terms of linear combinations of spin-components $\mathbf{S}_i = (S_i^x, S_i^y, S_i^z)$, in the global frame of the crystal axes — cf. $\mathcal{H}_{\text{ex}}$ [Eq. (1)]. Labelling of spins within the tetrahedron follows the convention of Ross et al. [43] — cf. Fig. 5. The notation $\Psi_i$ for ordered phases is taken from [64].</caption><thead><tr><th>order parameter</th><th>definition in terms of spin components</th><th>associated ordered phases</th></tr></thead><tbody><tr><td>$m_{A_2}$</td><td>$\frac{1}{2\sqrt{3}}(S_0^x + S_0^y + S_0^z + S_1^x - S_1^y - S_1^z - S_2^x + S_2^y - S_2^z - S_3^x - S_3^y + S_3^z)$</td><td>"all in-all out"</td></tr><tr><td>$\mathbf{m}_E$</td><td>$\left(\frac{1}{2\sqrt{6}}(-2S_0^x + S_0^y + S_0^z - 2S_1^x - S_1^y - S_1^z + 2S_2^x + S_2^y - S_2^z + 2S_3^x - S_3^y + S_3^z)\right)$</td><td>$\Psi_2$ and $\Psi_3$</td></tr><tr><td>$\mathbf{m}_{T_{1,A}}$</td><td>$\left(\frac{1}{2}(S_0^x + S_1^x + S_2^x + S_3^x)\right)$</td><td>collinear FM</td></tr><tr><td>$\mathbf{m}_{T_{1,B}}$</td><td>$\left(\frac{1}{2}(S_0^y + S_1^y + S_2^y + S_3^y)\right)$</td><td>non-collinear FM</td></tr><tr><td>$\mathbf{m}_{T_2}$</td><td>$\left(\frac{-1}{2\sqrt{2}}(S_0^x + S_0^z - S_1^y - S_1^z - S_2^y + S_2^z + S_3^y - S_3^z)\right)$</td><td>Palmer-Chalker ($\Psi_4$)</td></tr></tbody></table>

in terms of the parameters of $\mathcal{H}_{\text{ex}}^{\text{local}}$ [Eq. (8)]. We note that a similar analysis, applied to Er$_2$Ti$_2$O$_7$, appears in Ref. [42], with a different choice of basis for the two appearances of the $T_1$ irrep. Note that this is entirely different from the analysis appearing in Ref. [41] which classifies wavefunctions for spin-1/2 on a tetrahedron according to the irreps of $T_d$, not linear combinations of spin operators.

Symmetry permits a finite coupling $a_{T_{1,AB}} \neq 0$ between the two distinct $T_1$ irreps $\mathbf{m}_{T_{1,A}}$ and $\mathbf{m}_{T_{1,B}}$. This can be eliminated from $\mathcal{H}_{\text{ex}}^{\text{tet}}$ by a coordinate transformation

$$ \begin{aligned} \mathbf{m}_{T_{1,A'}} &= \cos \theta_{T_1} \mathbf{m}_{T_{1,A}} - \sin \theta_{T_1} \mathbf{m}_{T_{1,B}} \\ \mathbf{m}_{T_{1,B'}} &= \sin \theta_{T_1} \mathbf{m}_{T_{1,A}} + \cos \theta_{T_1} \mathbf{m}_{T_{1,B}} \end{aligned} \quad (19) $$

where

$$ \theta_{T_1} = \frac{1}{2} \arctan \left( \frac{\sqrt{8J_3}}{2J_1 + 2J_2 + J_3 - 2J_4} \right). \quad (20) $$

is the canting angle between spins and the relevant [100] axis in the ferromagnetic ground state. The Hamiltonian $\mathcal{H}_{\text{ex}}^{\text{tet}}$ then becomes

$$ \mathcal{H}_{\text{ex}}^{[T_d]} = \frac{1}{2} [a_{A_2} m_{A_2}^2 + a_E \mathbf{m}_E^2 + a_{T_2} \mathbf{m}_{T_2}^2 + a_{T_{1,A'}} \mathbf{m}_{T_{1,A'}}^2 + a_{T_{1,B'}} \mathbf{m}_{T_{1,B'}}^2]. \quad (21) $$

with coefficients given in Table V. We wish to emphasize that $\mathcal{H}_{\text{ex}}^{[T_d]}$ [Eq. (21)] is an *exact transcription* of $\mathcal{H}_{\text{ex}}^{\text{tet}}$ [Eq. (16)] and *not* a phenomenological Landau theory. As such, $\mathcal{H}_{\text{ex}}^{[T_d]}$ [Eq. (21)] is subject to the constraint that every spin has fixed length.

For the majority of the discussion in this article, we shall be concerned with classical vectors $\mathbf{S}_i$ representing a (pseudo)spin-1/2, with

$$ S = 1/2, \quad (22) $$

in which case

$$ |\mathbf{S}_i|^2 = 1/4. \quad (23) $$

For spins belonging to a single tetrahedron, we can use symmetry$^{155}$ to express this constraint as

$$ \begin{aligned} \mathbf{S}_0^2 + \mathbf{S}_1^2 + \mathbf{S}_2^2 + \mathbf{S}_3^2 &= 1 \\ \mathbf{S}_0^2 + \mathbf{S}_1^2 - \mathbf{S}_2^2 - \mathbf{S}_3^2 &= 0 \\ \mathbf{S}_0^2 - \mathbf{S}_1^2 + \mathbf{S}_2^2 - \mathbf{S}_3^2 &= 0 \\ \mathbf{S}_0^2 - \mathbf{S}_1^2 - \mathbf{S}_2^2 + \mathbf{S}_3^2 &= 0. \end{aligned} \quad (24) $$

The constraint of fixed spin-length, Eq. (24), plays a central role in determining the allowed classical ground states, below.

We note in passing that the addition of a single-ion anisotropy term $-D(\mathbf{S} \cdot \mathbf{z}_i^{\text{local}})^2$ can easily be included in the analysis by a simple modification of the coefficients $a_\lambda$ in $\mathcal{H}_{\text{ex}}^{[T_d]}$ [Eq. (21)], and so do not affect any of the conclusions reached about ground states, below. However, since interactions of this form contribute only to a trivial energy constant for a Kramers doublet, we will not pursue this point further here.
<table><caption>TABLE IV: Order parameters $\mathbf{m}_\lambda$, describing how the point-group symmetry of a single tetrahedron within the pyrochlore lattice is broken by magnetic order. Order parameters are irreducible representations of the point-group $T_d$, and are expressed in terms of linear combinations of spin-components $S_i = (S_i^x, S_i^y, S_i^z)$, in the local frame of the magnetic ions — cf. $H_{\text{ex}}^{\text{local}}$ [Eq. (8)]. For convenience, in this table, the local axes $(x^{\text{local}}, y^{\text{local}}, z^{\text{local}})$ are simply written $(x, y, z)$. Labelling of spins within the tetrahedron follows the convention of Ross et al. [43] — cf. Fig. 5.</caption><thead><tr><th>order parameter</th><th>definition in terms of local spin components</th></tr></thead><tbody><tr><td><b><i>m</i><sub>A<sub>2</sub></sub></b></td><td><sup>1</sup>&frasl;<sub>2</sub> (<i>S</i><sub>0</sub><sup><i>z</i></sup> + <i>S</i><sub>1</sub><sup><i>z</i></sup> + <i>S</i><sub>2</sub><sup><i>z</i></sup> + <i>S</i><sub>3</sub><sup><i>z</i></sup>)</td></tr><tr><td><b><i>m</i><sub>E</sub></b></td><td><sup>1</sup>&frasl;<sub>2</sub> (<i>S</i><sub>0</sub><sup><i>x</i></sup> + <i>S</i><sub>1</sub><sup><i>x</i></sup> + <i>S</i><sub>2</sub><sup><i>x</i></sup> + <i>S</i><sub>3</sub><sup><i>x</i></sup>)</td></tr><tr><td><b><i>m</i><sub>T<sub>1,A</sub></sub></b></td><td>(<sup>1</sup>&frasl;<sub>2&radic;3</sub>(-&radic;<span style="text-decoration:overline">2<i>S</i><sub>0</sub><sup><i>x</i></sup> + <i>S</i><sub>0</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">2<i>S</i><sub>1</sub><sup><i>x</i></sup> + <i>S</i><sub>1</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">2<i>S</i><sub>2</sub><sup><i>x</i></sup> - <i>S</i><sub>2</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">2<i>S</i><sub>3</sub><sup><i>x</i></sup> - <i>S</i><sub>3</sub><sup><i>z</i></sup></span></span></span>)<br>(<sup>1</sup>&frasl;<sub>12</sub>(&radic;<span style="text-decoration:overline">6<i>S</i><sub>0</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>0</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>0</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>1</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>1</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>1</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">6<i>S</i><sub>2</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>2</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>2</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>3</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>3</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>3</sub><sup><i>z</i></sup></span></span></span></span></span></span></span></span></span></span></sup>)<br>(<sup>1</sup>&frasl;<sub>12</sub>(&radic;<span style="text-decoration:overline">6<i>S</i><sub>0</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>0</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>0</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>1</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>1</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>1</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>2</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>2</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>2</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">6<i>S</i><sub>3</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>3</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>3</sub><sup><i>z</i></sup></span></span></span></span></span></span></span></span></span></sup>)</span></td></tr><tr><td><b><i>m</i><sub>T<sub>1,B</sub></sub></b></td><td>(<sup>1</sup>&frasl;<sub>2&radic;3</sub>(-&radic;<span style="text-decoration:overline">2<i>S</i><sub>0</sub><sup><i>x</i></sup> - <i>S</i><sub>0</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">2<i>S</i><sub>1</sub><sup><i>x</i></sup> - <i>S</i><sub>1</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">2<i>S</i><sub>2</sub><sup><i>x</i></sup> + <i>S</i><sub>2</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">2<i>S</i><sub>3</sub><sup><i>x</i></sup> + <i>S</i><sub>3</sub><sup><i>z</i></sup></span></span></span>)<br>(<sup>1</sup>&frasl;<sub>12</sub>(&radic;<span style="text-decoration:overline">6<i>S</i><sub>0</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>0</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>0</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>1</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>1</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>1</sub><sup><i>x</i></sup> + &radic;<span style="text-decoration:overline">6<i>S</i><sub>2</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>2</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>2</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>3</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>3</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>3</sub><sup><i>z</i></sup></span></span></span></span></span></span></span></span></span></sup>)<br>(<sup>1</sup>&frasl;<sub>12</sub>(&radic;<span style="text-decoration:overline">6<i>S</i><sub>0</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>0</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>0</sub><sup><i>z</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>1</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>1</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>1</sub><sup><i>x</i></sup> - &radic;<span style="text-decoration:overline">6<i>S</i><sub>2</sub><sup><i>x</i></sup> - 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>2</sub><sup><i>y</i></sup> + 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>2</sub><sup><i>z</i></sup> + &radic;<span style="text-decoration:overline">6<i>S</i><sub>3</sub><sup><i>x</i></sup> + 3&radic;<span style="text-decoration:overline">2<i>S</i><sub>3</sub><sup><i>y</i></sup> - 2&radic;<span style="text-decoration:overline">3<i>S</i><sub>3</sub><sup><i>z</i></sup></span></span></span></span></span></span></span></span></span></sup>)</span></td></tr><tr><td><b><i>m</i><sub>T<sub>2</sub></sub></b></td><td>(<sup>1</sup>&frasl;<sub>2</sub>(-<i>S</i><sub>0</sub><sup><i>y</i></sup> - <i>S</i><sub>1</sub><sup><i>y</i></sup> + <i>S</i><sub>2</sub><sup><i>y</i></sup> + <i>S</i><sub>3</sub><sup><i>y</i></sup>)<br>(<sup>1</sup>&frasl;<sub>4</sub>(<span style="text-decoration:overline">&radic;<span style="text-decoration:overline">3<i>S</i><sub>0</sub><sup><i>x</i></sup> + <i>S</i><sub>0</sub><sup><i>y</i></sup> - &radic;<span style="text-decoration:overline">3<i>S</i><sub>1</sub><sup><i>x</i></sup> - <i>S</i><sub>1</sub><sup><i>y</i></sup> + &radic;<span style="text-decoration:overline">3<i>S</i><sub>2</sub><sup><i>x</i></sup> + <i>S</i><sub>2</sub><sup><i>y</i></sup> - &radic;<span style="text-decoration:overline">3<i>S</i><sub>3</sub><sup><i>x</i></sup> - <i>S</i><sub>3</sub><sup><i>y</i></sup></span></span></span>)<br>(<sup>1</sup>&frasl;<sub>4</sub>(-<span style="text-decoration:overline">&radic;<span style="text-decoration:overline">3<i>S</i><sub>0</sub><sup><i>x</i></sup> + <i>S</i><sub>0</sub><sup><i>y</i></sup> + &radic;<span style="text-decoration:overline">3<i>S</i><sub>1</sub><sup><i>x</i></sup> - <i>S</i><sub>1</sub><sup><i>y</i></sup> + &radic;<span style="text-decoration:overline">3<i>S</i><sub>2</sub><sup><i>x</i></sup> - <i>S</i><sub>2</sub><sup><i>y</i></sup> - &radic;<span style="text-decoration:overline">3<i>S</i><sub>3</sub><sup><i>x</i></sup> + <i>S</i><sub>3</sub><sup><i>y</i></sup></span></span></span>)</span></td></tr></tbody></table>

### III. ANALYSIS OF CLASSICAL PHASE DIAGRAM FOR $T=0$

#### A. General considerations

Given the existence of a classical ground state with $\mathbf{q} = 0$, 4-sublattice order, it is easy to determine a ground-state phase diagram directly from the Hamiltonian $H_{\text{ex}}^{[T_d]}$ [Eq. (21)]. The method developed below is quite general and can be applied for arbitrary $(J_1, J_2, J_3, J_4)$. However, for concreteness, we concentrate on the limit

$$J_3 < 0, \quad J_4 \equiv 0, \qquad (25)$$

which is of particular relevance to known pyrochlore materials, leading to the phases shown in Fig. 1—Fig. 3.

The classical ground state of $H_{\text{ex}}^{[T_d]}$ [Eq. (21)] can be found by first identifying the irrep $\lambda^*$ for which $a_\lambda^*$ takes on the minimum value, and then imposing the constraint on the total length of the spin [Eq. (24)] on $\mathbf{m}_\lambda^*$, which implies that

$$m_{A_2}^2 + \mathbf{m}_E^2 + \mathbf{m}_{T_2}^2 + \mathbf{m}_{T_{1,A}'}^2 + \mathbf{m}_{T_{1,B}'}^2 \equiv \sum_{\lambda} \mathbf{m}_{\lambda}^2 = 1. \quad (26)$$

Such an approach is possible because each individual order parameter $\mathbf{m}_\lambda$ can reach a maximum value of unity within physical spin configurations

$$\max \mathbf{m}_{\lambda}^{2} = 1. \qquad (27)$$

This method of determining the classical ground state is completely general and, once generalized to the lattice, is not restricted to conventionally ordered states$^{47,48}$.

In the limit $J_3 < 0$, $J_4 \equiv 0$, the coefficients $a_\lambda$ with the lowest values are $a_E$, $a_{T_{1,A}'}$, or $a_{T_2}$, and the corresponding $\mathbf{q} = 0$ ordered ground states found have $E, T_1$

and $T_2$ symmetry. The boundaries between these phases occur where

$$a_{T_2} = a_E < a_{T_2}, a_{T_{1B}'}, a_{A_2} \\ \Rightarrow J_2 = J_1 > 0 \qquad (28)$$

$$a_{T_2} = a_{T_{1,A'}} < a_E, a_{T_{1B}'}, a_{A_2} \\ \Rightarrow J_2 = -J_1 > 0 \qquad (29)$$

$$a_E = a_{T_{1,A'}} < a_{T_2}, a_{T_{1B}'}, a_{A_2} \\ \Rightarrow J_2 = \frac{J_1(4J_1 - 5J_3)}{4J_1 - J_3} < 0 \qquad (30)$$

as illustrated in Fig. 7.

In what follows, we explore the classical ground states with $E, T_1$ and $T_2$ symmetry in some detail, paying particular attention to what happens on phase boundaries where more than one order parameter can take on a finite value. We will not consider the “all-in, all out” ground state, a simple Ising-type order with two degenerate ground states. “All-in, all out” order has a finite value of the order parameter $m_{A_2}$ [cf Table IV], which requires $a_{A_2}$ to be the lowest coefficient. This only occurs for $J_4 > 0$ and/or $J_3 > 0$, and so falls outside the scope of this Article.

#### B. Non-collinear FM with $T_1$ symmetry

We begin by considering what happens where interactions are predominantly ferromagnetic (i.e. $J_1, J_2 < 0$), as in Yb$_2$Ti$_2$O$_7$ [43] — cf. Table I. For most of this region, as might be expected, the classical configuration with the lowest energy is a state with a finite magnetisation. This is the ground state throughout the region bounded by $a_{T_{1,A'}} = a_{T_2}$ [Eq. (29)], and $a_{T_{1,A'}} = a_E$ [Eq. (30)] — cf. Fig. 7. Here the energy is minimised by setting

$$\mathbf{m}_{T_{1,A'}}^2 = 1 \qquad (31)$$
<table><caption>TABLE V: Coefficients $a_λ$ of the scalar invariants $|m_λ|^2$ appearing in $H_{ex}^{[T_d]}$ [Eq. (21)]. Coefficients are expressed as a function of $(J_1, J_2, J_3, J_4)$, the parameters of $H_{ex}$ [Eq. (1)]; and $(J_{zz}, J_±, J_{±±}, J_{z±})$, the parameters of $H_{ex}^{local}$ [Eq. (8)], with the canting angle $θ_{T_1}$ defined in Eq. (20). The classical ground states of $H_{ex}$ [Eq. (4)] can be found by identifying the coefficient(s) $a_λ$ with the lowest value, and imposing the constraint of fixed spin-length, Eq. (24), on the associated $m_λ$.</caption><thead><tr><th>coefficient<br>of |m<sub>λ</sub>|<sup>2</sup></th><th>definition in terms of<br>parameters of H<sub>ex</sub> [Eq. (1)]</th><th>definition in terms of<br>parameters of H<sub>ex</sub><sup>local</sup> [Eq. (8)]</th></tr></thead><tbody><tr><td>a<sub>A<sub>2</sub></sub></td><td>-2J<sub>1</sub> + J<sub>2</sub> - 2(J<sub>3</sub> + 2J<sub>4</sub>)</td><td>3J<sub>zz</sub></td></tr><tr><td>a<sub>E</sub></td><td>-2J<sub>1</sub> + J<sub>2</sub> + J<sub>3</sub> + 2J<sub>4</sub></td><td>-6J<sub>±</sub></td></tr><tr><td>a<sub>T<sub>2</sub></sub></td><td>-J<sub>2</sub> + J<sub>3</sub> - 2J<sub>4</sub></td><td>2J<sub>±</sub> - 4J<sub>±±</sub></td></tr><tr><td>a<sub>T<sub>1,A'</sub></sub></td><td>(2J<sub>1</sub> + J<sub>2</sub>) cos<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>-(J<sub>2</sub> + J<sub>3</sub> - 2J<sub>4</sub>) sin<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>+√2J<sub>3</sub> sin(2θ<sub>T<sub>1</sub></sub>)</td><td>(1/3)(4J<sub>±</sub> + 8J<sub>±±</sub> + 8√2J<sub>z±</sub> - J<sub>zz</sub>) cos<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>+ (1/3)(1J<sub>±</sub> + 2J<sub>±±</sub> - 4√2J<sub>z±</sub> - J<sub>zz</sub>) sin<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>+ (√2/3)(-2J<sub>±</sub> - 4J<sub>±±</sub> + 2√2J<sub>z±</sub> - J<sub>zz</sub>) sin(2θ<sub>T<sub>1</sub></sub>)</td></tr><tr><td>a<sub>T<sub>1,B'</sub></sub></td><td>(2J<sub>1</sub> + J<sub>2</sub>) sin<sup>2</sup>(θ<sub>T<sub>1</sub></sub>) -<br>(J<sub>2</sub> + J<sub>3</sub> - 2J<sub>4</sub>) cos<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>-√2J<sub>3</sub> sin(2θ<sub>T<sub>1</sub></sub>)</td><td>(1/3)(4J<sub>±</sub> + 8J<sub>±±</sub> + 8√2J<sub>z±</sub> - J<sub>zz</sub>) sin<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>- (1/3)(1J<sub>±</sub> + 2J<sub>±±</sub> - 4√2J<sub>z±</sub> - J<sub>zz</sub>) cos<sup>2</sup>(θ<sub>T<sub>1</sub></sub>)<br>- (√2/3)(-2J<sub>±</sub> - 4J<sub>±±</sub> + 2√2J<sub>z±</sub> - J<sub>zz</sub>) sin(2θ<sub>T<sub>1</sub></sub>)</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>FIG. 7: Classical ground state phase diagram of $H_{ex}$ [Eq. (4)] for $J_3 < 0$, $J_4 = 0$, as a function of $(J_1, J_2)/|J_3|$. In the absence of fluctuations, the ground states are a non-collinear FM transforming with the $T_1$ irrep of $T_d$; a one-dimensional manifold of states transforming with the E irrep of $T_d$; and the Palmer-Chalker phase, a coplanar antiferromagnet transforming with the $T_2$ irrep of $T_d$. All three phases have long-range 4-sublattice order. Analytical expressions for the boundaries between phases are given in Eq. (28-30), with coefficients $a_λ$ defined in Table V.</figcaption></figure>

and

$$m_{A_2} = m_E = m_{T_2} = m_{T_{1B'}} = 0 \quad (32)$$

The constraints on the spin lengths [Eq. (24)] further imply that

$$ \begin{aligned} m_{T_{1A'}}^y, m_{T_{1A'}}^z &= 0 \\ m_{T_{1A'}}^x, m_{T_{1A'}}^z &= 0 \\ m_{T_{1A'}}^x, m_{T_{1A'}}^y &= 0. \end{aligned} \quad (33) $$

It follows that there are 6 possible ground states

$$ \mathbf{m}_{T_{1A'}} = \begin{pmatrix} \pm 1 \\ 0 \\ 0 \end{pmatrix}, \quad \begin{pmatrix} 0 \\ \pm 1 \\ 0 \end{pmatrix}, \quad \begin{pmatrix} 0 \\ 0 \\ \pm 1 \end{pmatrix}. \quad (34) $$

Written in terms of spins, these are 6, non-collinear ferromagnetic (FM) ground states, with typical spin configuration

$$ \begin{aligned} S_0 &= S \left( \sin \theta_{T_1} / \sqrt{2}, \sin \theta_{T_1} / \sqrt{2}, \cos \theta_{T_1} \right) \\ S_1 &= S \left( -\sin \theta_{T_1} / \sqrt{2}, \sin \theta_{T_1} / \sqrt{2}, \cos \theta_{T_1} \right) \\ S_2 &= S \left( \sin \theta_{T_1} / \sqrt{2}, -\sin \theta_{T_1} / \sqrt{2}, \cos \theta_{T_1} \right) \\ S_3 &= S \left( -\sin \theta_{T_1} / \sqrt{2}, -\sin \theta_{T_1} / \sqrt{2}, \cos \theta_{T_1} \right) \end{aligned} \quad (35) $$

where $θ_{T_1}$ is given by Eq. (20) and, following Eq. (22), $S = 1/2$.

The magnetisation of this FM ground state, illustrated in Fig. 8, is parallel to a [001] axis, with spins canted away from this axis, in an “ice-like” manner. This state has been identified as the ground state in Yb<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>, where it was referred to as a “splayed FM” [36], and in most samples of Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> which order at low temperature<sup>24,65–71</sup> although a different form of canting has recently been claimed in Ref. 72. It is also the observed ordered state of the Tb based pyrochlore Tb<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub> [73,74]

### C. One-dimensional manifold of states with E symmetry

For a wide range of parameters, predominantly with antiferromagnetic “XY” interactions $J_1 > 0$, the classical ground state of $H_{ex}^{[T_d]}$ [Eq. (21)] is a one-dimensional manifold of states which transforms with the E irrep of $T_d$. These ground states occur in a region bounded by $a_E = a_{T_{1,A'}}$ [Eq. (30)] and $a_E = a_{T_2}$ [Eq. (28)] — cf Fig. 7 — and is characterised by spins lying in the
<figure><img src="image_9.png" /><figcaption>FIG. 8: Spin-configuration in the 4-sublattice non-collinear FM phase, transforming with the $T_1$ irrep of $T_d$: (a) viewed along the [001] axis; (b) viewed slightly off the [110] axis. The magnetisation is aligned with the [001] axis. Spins are canted into the plane perpendicular to this, with canting angle $\theta_{T_1}$. The canting is of an “ice-like” form such that the projection of the spin configuration onto a [001] plane has two spins oriented into the tetrahedron and two oriented out, as shown in (a).</figcaption></figure>

“XY” plane normal to the local [111] axis on each site [cf Eqs. (A2,A3,A4)].

For this range of parameters, the classical ground state energy can be minimised by setting

$$ \mathbf{m}_{\mathrm{E}}^{2}=1 \qquad (36) $$

and

$$ m_{A_2} = \mathbf{m}_{T_2} = \mathbf{m}_{T_{1A'}} = \mathbf{m}_{T_{1B'}} = 0. \qquad (37) $$

These solutions *automatically* satisfy the constraint on the total length of the spin Eq. (24) and are conveniently characterised by writing

$$ \mathbf{m}_{\mathrm{E}} = (\cos \theta_{\mathrm{E}}, \sin \theta_{\mathrm{E}}). \qquad (38) $$

It follows that the ground state is a continuous, one-dimensional manifold of states parameterised by the single angle $0 \le \theta_{\mathrm{E}} < 2\pi$. The spin configuration in this

<figure><img src="image_10.png" /><figcaption>FIG. 9: Example of a spin configuration within the one-dimensional manifold of states transforming with the E irrep of $T_d$: (a) viewed along [001] axis; (b) viewed slightly off the [110] axis. The manifold possesses 4-sublattice long-range order, with spins lying in the “XY” plane perpendicular to the local [111] axis at each site. The manifold is continuous, and can be parameterised with a single angle $\theta_E$. The manifold can be generated by a clockwise rotation of all spins around their respective local axes.</figcaption></figure>

manifold is given by

$$ \begin{align} \mathbf{S}_0 &= S \left( \sqrt{\frac{2}{3}} \cos(\theta_E), \sqrt{\frac{2}{3}} \cos(\theta_E + \frac{2\pi}{3}), \sqrt{\frac{2}{3}} \cos(\theta_E - \frac{2\pi}{3}) \right) \nonumber \\ \mathbf{S}_1 &= S \left( \sqrt{\frac{2}{3}} \cos(\theta_E), -\sqrt{\frac{2}{3}} \cos(\theta_E + \frac{2\pi}{3}), -\sqrt{\frac{2}{3}} \cos(\theta_E - \frac{2\pi}{3}) \right) \nonumber \\ \mathbf{S}_2 &= S \left( -\sqrt{\frac{2}{3}} \cos(\theta_E), \sqrt{\frac{2}{3}} \cos(\theta_E + \frac{2\pi}{3}), -\sqrt{\frac{2}{3}} \cos(\theta_E - \frac{2\pi}{3}) \right) \nonumber \\ \mathbf{S}_3 &= S \left( -\sqrt{\frac{2}{3}} \cos(\theta_E), -\sqrt{\frac{2}{3}} \cos(\theta_E + \frac{2\pi}{3}), \sqrt{\frac{2}{3}} \cos(\theta_E - \frac{2\pi}{3}) \right). \tag{39} \end{align} $$

with spins lying in the local “XY” plane [cf. Eqs. (A3) and (A4)].

### D. Non-coplanar antiferromagnet, $\Psi_2$, with E symmetry

It is now well-established that in Er$_2$Ti$_2$O$_7$ quantum fluctuations$^{26,27,45}$, classical thermal fluctuations at low temperature$^{25}$ and thermal fluctuations near the ordering temperature$^{28,78}$ fluctuations all act within the one-dimensional manifold of classical ground states described in Section III C, to select a non-coplanar antiferromagnet, $\Psi_2$, illustrated in Fig. (10). Structural disorder, meanwhile, favours the coplanar antiferromagnet, $\Psi_3$.
<figure><img src="image_11.png" /><figcaption>FIG. 10: Spin configuration in the 4-sublattice non-coplanar antiferromagnet, $\Psi_2$, selected by fluctuations from the one-dimensional manifold of states transforming with $E$: (a) viewed along [001] axis; (b) viewed slightly off the [110] axis. At the phase boundary with the Palmer-Chalker phase, each of the six $\Psi_2$ ground states can be transformed continuously into a Palmer-Chalker state without leaving the ground-state manifold.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11: Spin configuration in the 4-sublattice coplanar antiferromagnet, $\Psi_3$, selected by fluctuations from the one-dimensional manifold of states transforming with $E$: (a) viewed along [001] axis; (b) viewed slightly off the [110] axis. At the phase boundary with the non-collinear FM phase, each of the six $\Psi_3$ ground states can be transformed continuously into a non-collinear FM state, without leaving the ground-state manifold.</figcaption></figure>

[75,76], illustrated in Fig. (11). Together, this pair of states form a basis for the $E$ irrep of $T_d$ [64].

The $\Psi_2$ ground state is six-fold degenerate, with spins canted symmetrically out of the [100] plane. The six spin configurations for $\Psi_2$ states are given by Eq. (39) with $\theta_E = \frac{n\pi}{3}$, $n = 0, 1, 2, \dots, 5$. The $\Psi_2$ state is characterised by the primary order parameter $m_E$ [cf. Table III], and by $c_E > 0$, where (cf. Refs. [77,78])

$$c_E = \langle \cos 6\theta_E \rangle \quad (40)$$

Symmetry allows for fluctuations to induce a finite value of $m_{A_2}$ in the $\Psi_2$ state$^{79}$, but classically this must vanish as $T \to 0$ since the energy is minimised by $m_{A_2} = 0$ within the region of phase diagram which favours the $\Psi_2$ state.

## E. Coplanar antiferromagnet, $\Psi_3$, with $E$ symmetry

For parameters bordering on the non-collinear FM phase, fluctuations select a coplanar antiferromagnet, $\Psi_3$, from the one-dimensional manifold of states transforming with $E$. The $\Psi_3$ ground state is six-fold degenerate, with spins lying in a common [100] plane.

The six spin configurations for $\Psi_3$ states are given by Eq. (39) with $\theta_E = \frac{(2n+1)\pi}{6}$, $n = 0, 1, 2, \dots, 5$. These states are characterised by a finite value of the order parameter $m_E$ [cf Table III], and by $c_E < 0$ [cf. Eq. (40)]. An example of a typical spin configuration is shown in Fig. (11).

Taken together $\Psi_2$ and $\Psi_3$ form a complete basis for the $E$ irrep of $T_d$.

## F. Palmer-Chalker phase, $\Psi_4$, with $T_2$ symmetry

In a region bounded by $a_{T_2} = a_{T_{1,A'}}$ [Eq. (29)] and $a_{T_2} = a_E$ [Eq. (28)] — cf. Fig. 7 — the energy is min-

imised by setting

$$\mathbf{m}_{T_2}^2 = 1 \quad (41)$$

and

$$m_{A_2} = m_E = m_{T_{1A'}} = m_{T_{1B'}} = 0 \quad (42)$$

The constraints on the total length of the spin, Eq. (24), further imply that

$$\mathbf{m}_{T_2}^2 = 1 \quad (43)$$

$$m_{T_2}^y m_{T_2}^z = 0 \quad (44)$$

$$m_{T_2}^x m_{T_2}^z = 0 \quad (45)$$

$$m_{T_2}^x m_{T_2}^y = 0 \quad (46)$$

giving us a set of 6 ground states

$$\mathbf{m}_{T_2} = \begin{pmatrix} \pm 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \pm 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \pm 1 \end{pmatrix}. \quad (47)$$

Within these ground states, spins are arranged in helical manner in a common [100] plane, with a typical spin configuration given by (see Fig. 12).

$$\begin{align}
\mathbf{S}_0 &= S \left( \frac{1}{\sqrt{2}}, -\frac{1}{\sqrt{2}}, 0 \right) \nonumber \\
\mathbf{S}_1 &= S \left( -\frac{1}{\sqrt{2}}, -\frac{1}{\sqrt{2}}, 0 \right) \nonumber \\
\mathbf{S}_2 &= S \left( \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right) \nonumber \\
\mathbf{S}_3 &= S \left( -\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right) \tag{48}
\end{align}$$

This phase is the “Palmer-Chalker” phase, first identified as the ground state of a model with antiferromagnetic nearest neighbour Heisenberg interactions and long-range
<figure><img src="image_13.png" /><figcaption>FIG. 12: Spin configuration in the 4-sublattice Palmer-Chalker phase, $\Psi_4$, transforming with the $T_2$ irrep of $T_d$: (a) viewed along [001] axis; (b) viewed slightly off the [110] axis. At the phase boundary with the $\Psi_2$ phase, each of the six Palmer-Chalker ground states can be transformed continuously into a $\Psi_2$ state.</figcaption></figure>

dipolar interactions on the pyrochlore lattice$^{46}$. The 6 degenerate ground states in this phase are described by the basis vectors $\Psi_{4,5,6}$ [64], which are all equivalent by symmetry. For brevity we refer to this phase as $\Psi_4$, but it should be remembered that all three basis vectors $\Psi_{4,5,6}$ are equivalent ground states.

The Palmer-Chalker states are superficially similar to the $\Psi_3$ states [Section III E], being coplanar, antiferromagnet configurations with all spins lying in a [100] plane. However, their symmetry properties are quite different (as is expressed in the fact that they transform according to different irreps of $T_d$.) A simple example of these different symmetry properties is their behaviour under the 3 $C_2$ rotations around $\langle 100 \rangle$ axes. The $\Psi_3$ configurations are invariant under all such rotations. The PC states are invariant under one such rotation (the one around the axis normal to all the spins), but reverse all spin orientations under the other two.

### G. Boundary between Palmer-Chalker phase and the one-dimensional manifold of states with E symmetry

The boundary between the Palmer-Chalker phase and the one-dimensional manifold of states with E symmetry occurs when $a_E = a_{T_2}$ [cf. Eq. (28)]. In this case, $H_{ex}^{[T_d]}$ [Eq. (21)] is minimized by setting

$$ \mathbf{m}_E^2 + \mathbf{m}_{T_2}^2 = 1 \qquad (49) $$

and

$$ m_{A_2} = \mathbf{m}_{T_{1A'}} = \mathbf{m}_{T_{1B'}} = 0. \qquad (50) $$

Substituting

$$ \mathbf{m}_E = m_E(\cos \theta_E, \sin \theta_E), \qquad (51) $$

<figure><img src="image_14.png" /><figcaption>FIG. 13: Structure of the ground state manifold at the boundary between the Palmer-Chalker (PC) phase and the one-dimensional manifold of states with E symmetry. The black circle denotes the manifold of E-symmetry ground states, including the six $\Psi_2$ ground states (black dots). At the boundary with the PC phase, this manifold branches at the $\Psi_2$ states, to connect with three, additional, one-dimensional manifolds. These manifolds in turn interpolate to the six Palmer-Chalker ground states with $T_2$ symmetry (red dots). An exactly equivalent picture holds on the boundary between the non-collinear ferromagnet (FM), and the one-dimensional manifold of states with E symmetry. However in this case the different manifolds intersect at the $\Psi_3$ states.</figcaption></figure>

and imposing the constraint Eq. (24), we find

$$ 2m_E m_{T_2}^x \sin(\theta_E) - m_{T_2}^y m_{T_2}^z = 0 $$

$$ 2m_E m_{T_2}^y \sin \left( \theta_E - \frac{2\pi}{3} \right) - m_{T_2}^x m_{T_2}^z = 0 $$

$$ 2m_E m_{T_2}^z \sin \left( \theta_E + \frac{2\pi}{3} \right) - m_{T_2}^x m_{T_2}^y = 0. \qquad (52) $$

It is easy to show that there are no solutions to Eqs. (52) where more than one component of $\mathbf{m}_{T_2}$ is finite. There are, however, three distinct one-dimensional manifolds which connect pairs of Palmer-Chalker states
to the one-dimensional manifold of E-symmetry states:

$$ \mathbf{m}_E = \cos(\alpha) \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \quad \mathbf{m}_{\mathcal{T}_2} = \sin(\alpha) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \quad (53) $$

$$ \mathbf{m}_E = \cos(\beta) \begin{pmatrix} -\frac{1}{2} \\ \frac{\sqrt{3}}{2} \end{pmatrix}, \quad \mathbf{m}_{\mathcal{T}_2} = \sin(\beta) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \quad (54) $$

$$ \mathbf{m}_E = \cos(\gamma) \begin{pmatrix} -\frac{1}{2} \\ -\frac{\sqrt{3}}{2} \end{pmatrix}, \quad \mathbf{m}_{\mathcal{T}_2} = \sin(\gamma) \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \quad (55) $$

where the angles $\alpha, \beta$ and $\gamma$ run from 0 to $2\pi$.

A typical spin configuration for one of the three connecting manifolds is

$$ \mathbf{S}_0 = S \sqrt{\frac{2}{3}} \left( -\cos(\alpha), \cos\left(\alpha + \frac{\pi}{3}\right), \cos\left(\alpha - \frac{\pi}{3}\right) \right) $$

$$ \mathbf{S}_1 = S \sqrt{\frac{2}{3}} \left( -\cos(\alpha), -\cos\left(\alpha + \frac{\pi}{3}\right), -\cos\left(\alpha - \frac{\pi}{3}\right) \right) $$

$$ \mathbf{S}_2 = S \sqrt{\frac{2}{3}} \left( \cos(\alpha), \cos\left(\alpha - \frac{\pi}{3}\right), -\cos\left(\alpha + \frac{\pi}{3}\right) \right), $$

$$ \mathbf{S}_3 = S \sqrt{\frac{2}{3}} \left( \cos(\alpha), -\cos\left(\alpha - \frac{\pi}{3}\right), \cos\left(\alpha + \frac{\pi}{3}\right) \right) \quad (56) $$

where $\alpha = 0, \pi$ correspond to the $\Psi_2$ ground states with $\theta_E = 0, \pi$, and $\alpha = \pi/2, 3\pi/2$ correspond to two of the six Palmer-Chalker ground states. These manifolds are illustrated in Fig. 13.

By application of the “Lego-brick” rules given in Section II E, the ground state degeneracy on this phase boundary must be at least $O(2^{L^2})$. This follows from the fact that, on the boundary, the Palmer-Chalker ground states $\Psi_4$ share two spin orientations with the neighbouring $\Psi_3$ configuration. This, in turn, is connected with the $O(L^2)$ number of zero modes appearing in spin wave expansions around the $\Psi_2$ configurations at this boundary [cf. Section IV B].

We note that a special case of this phase boundary is studied for pure XY spins (i.e. with infinite easy-plane anisotropy) in Ref. [80], finding the same structure of ground state manifolds.

### H. Boundary between the non-collinear ferromagnet and the one-dimensional manifold of states with E symmetry

The boundary between the non-collinear ferromagnet and the one-dimensional manifold of states with E symmetry occurs when $a_E = a_{\mathcal{T}_{1,A'}}$ [cf. Eq. (30)]. In this case, $H_{\text{ex}}^{[T_d]}$ [Eq. (21)] is minimised by setting

$$ \mathbf{m}_{\mathcal{T}_{1,A'}}^2 + \mathbf{m}_E^2 = 1 \quad (57) $$

and

$$ m_{A_2} = \mathbf{m}_{\mathcal{T}_{1,B'}} = \mathbf{m}_{\mathcal{T}_2} = 0. \quad (58) $$

Defining, for the sake of brevity, the quantities

$$ \begin{aligned} \mu(\theta_{\mathcal{T}_1}) &= (\sqrt{2} \cos(\theta_{\mathcal{T}_1}) - \sin(\theta_{\mathcal{T}_1})) \\ \nu(\theta_{\mathcal{T}_1}) &= (\sin(\theta_{\mathcal{T}_1})^2 + \sqrt{2} \sin(2\theta_{\mathcal{T}_1})) \end{aligned} \quad (59) $$

and imposing the constraint Eq. (24) we obtain

$$ \begin{aligned} 2m_E m_{\mathcal{T}_{1,A'}}^x \cos(\theta_E) &= -\frac{\mu(\theta_{\mathcal{T}_{1,A'}})}{\nu(\theta_{\mathcal{T}_{1,A'}})} m_{\mathcal{T}_{1,A'}}^y m_{\mathcal{T}_{1,A'}}^z \\ 2m_E m_{\mathcal{T}_{1,A'}}^y \cos\left(\theta_E - \frac{2\pi}{3}\right) &= -\frac{\mu(\theta_{\mathcal{T}_{1,A'}})}{\nu(\theta_{\mathcal{T}_{1,A'}})} m_{\mathcal{T}_{1,A'}}^x m_{\mathcal{T}_{1,A'}}^z \\ 2m_E m_{\mathcal{T}_{1,A'}}^z \cos\left(\theta_E + \frac{2\pi}{3}\right) &= -\frac{\mu(\theta_{\mathcal{T}_{1,A'}})}{\nu(\theta_{\mathcal{T}_{1,A'}})} m_{\mathcal{T}_{1,A'}}^x m_{\mathcal{T}_{1,A'}}^y \end{aligned} \quad (60) $$

where $\theta_{\mathcal{T}_1}$ is the (fixed) canting angle [Eq. (20)], $\theta_E$ is the (variable) angle within the U(1) manifold [Eq. (38)]. For the parameters considered here, the quantities $\mu(\theta_{\mathcal{T}_1})$ and $\nu(\theta_{\mathcal{T}_1})$ are always finite.

Arguments identical to those developed for the boundary with the Palmer-Chalker phase, give us three further 1D manifolds in addition to that associated with the E phase. However the intersections of the manifolds are now located at $\theta_E = \frac{(2n+1)\pi}{6}$, corresponding to the $\Psi_3$ states. This explains the model’s general entropic preference for $\Psi_3$ states in the region proximate to the ferromagnetic phase.

A typical spin configuration for one of the three connecting manifolds, parameterised by an angle $\eta$ is

$$ \begin{aligned} \mathbf{S}_0 &= S \left( \cos(\theta_{\mathcal{T}_1}) \sin(\eta), \frac{1}{\sqrt{2}}(-\cos(\eta) + \sin(\eta) \sin(\theta_{\mathcal{T}_1})), \right. \\ & \qquad \left. \frac{1}{\sqrt{2}}(\cos(\eta) + \sin(\eta) \sin(\theta_{\mathcal{T}_1})) \right) \\ \mathbf{S}_1 &= S \left( \cos(\theta_{\mathcal{T}_1}) \sin(\eta), \frac{1}{\sqrt{2}}(\cos(\eta) - \sin(\eta) \sin(\theta_{\mathcal{T}_1})), \right. \\ & \qquad \left. \frac{1}{\sqrt{2}}(-\cos(\eta) - \sin(\eta) \sin(\theta_{\mathcal{T}_1})) \right) \\ \mathbf{S}_2 &= S \left( \cos(\theta_{\mathcal{T}_1}) \sin(\eta), \frac{1}{\sqrt{2}}(-\cos(\eta) - \sin(\eta) \sin(\theta_{\mathcal{T}_1})), \right. \\ & \qquad \left. \frac{1}{\sqrt{2}}(-\cos(\eta) + \sin(\eta) \sin(\theta_{\mathcal{T}_1})) \right) \\ \mathbf{S}_3 &= S \left( \cos(\theta_{\mathcal{T}_1}) \sin(\eta), \frac{1}{\sqrt{2}}(\cos(\eta) + \sin(\eta) \sin(\theta_{\mathcal{T}_1})), \right. \\ & \qquad \left. \frac{1}{\sqrt{2}}(\cos(\eta) - \sin(\eta) \sin(\theta_{\mathcal{T}_1})) \right). \end{aligned} \quad (61) $$

Here $\eta = 0$ corresponds to the $\Psi_3$ ground state with $\theta_E = \pi/2$, and $\eta = \pi/2$ to one of the six FM ground states.

We note that an equivalent ground-state manifold was discussed by Canals and coauthors$^{81}$, and later by Chern$^{77}$, in the context the Heisenberg antiferromagnet with Dzyaloshinskii-Moriya interactions on the pyrochlore lattice. This case corresponds to a single point on the phase boundary considered here.
# I. Boundary between the Palmer-Chalker phase and the non-collinear ferromagnet

The boundary between the Palmer-Chalker phase and the non-collinear ferromagnet occurs when $a_{T_2} = a_{T_{1,A'}}$ [cf. Eq. (29)]. In this case, $H_{ex}^{[T_d]}$ [Eq. (21)] is minimised by setting

$$ \mathbf{m}_{T_2}^2 + \mathbf{m}_{T_{1,A'}}^2 = 1 \qquad (62) $$

and

$$ m_{A_2} = \mathbf{m}_E = \mathbf{m}_{T_{1,B'}} = 0. \qquad (63) $$

Imposing the constraint Eq. (24) we obtain

$$ \begin{aligned} & -m_{T_2}^y m_{T_2}^{\tilde{z}} + (\sin(\theta_{T_1})^2 + \sqrt{2}\sin(2\theta_{T_1}))m_{T_{1A'}}^y m_{T_{1A'}}^{\tilde{z}} \\ & \qquad + (\sqrt{2}\cos(\theta_{T_1}) - \sin(\theta_{T_1})) (\mathbf{m}_{T_{1A'}} \times \mathbf{m}_{T_2})_x = 0 \\ & -m_{T_2}^x m_{T_2}^{\tilde{z}} + (\sin(\theta_{T_1})^2 + \sqrt{2}\sin(2\theta_{T_1}))m_{T_{1A'}}^x m_{T_{1A'}}^{\tilde{z}} \\ & \qquad + (\sqrt{2}\cos(\theta_{T_1}) - \sin(\theta_{T_1})) (\mathbf{m}_{T_{1A'}} \times \mathbf{m}_{T_2})_y = 0 \\ & -m_{T_2}^x m_{T_2}^y + (\sin(\theta_{T_1})^2 + \sqrt{2}\sin(2\theta_{T_1}))m_{T_{1A'}}^x m_{T_{1A'}}^y \\ & \qquad + (\sqrt{2}\cos(\theta_{T_1}) - \sin(\theta_{T_1})) (\mathbf{m}_{T_{1A'}} \times \mathbf{m}_{T_2})_z = 0 \end{aligned} \qquad (64) $$

where $\theta_{T_1}$ is defined in Eq. (20).

In general, the ground state manifold on the boundary of the Palmer-Chalker phase is locally two-dimensional. To establish this, we consider small deviations from a given solution

$$ \begin{aligned} \mathbf{m}_{T_2} &= \mathbf{m}_{T_2}^0 + \delta\mathbf{m}_{T_2} \\ \mathbf{m}_{T_{1A'}} &= \mathbf{m}_{T_{1A'}}^0 + \delta\mathbf{m}_{T_{1A'}} \end{aligned} \qquad (65) $$

and expand the constraint Eq. (64) to linear order in $\delta\mathbf{m}$. Generally, we find two linearly-independent solutions for $(\delta\mathbf{m}_{T_2}, \delta\mathbf{m}_{T_{1A'}})$, and the manifold in the vicinity of $(\mathbf{m}_{T_2}^0, \mathbf{m}_{T_{1A'}}^0)$ is two-dimensional.

However if we expand around a state $(\tilde{\mathbf{m}}_{T_2}^0, \tilde{\mathbf{m}}_{T_{1A'}}^0)$ where both order parameters are aligned with the same cubic axis, e.g.

$$ \tilde{\mathbf{m}}_{T_2}^{0y} = \tilde{\mathbf{m}}_{T_2}^{0z} = \tilde{\mathbf{m}}_{T_{1A'}}^{0y} = \tilde{\mathbf{m}}_{T_{1A'}}^{0z} = 0 \qquad (66) $$

one of the Eqs. (64) is satisfied trivially, leaving only three constraints on six variables. It follows that the manifold is locally three-dimensional in the vicinity of $(\tilde{\mathbf{m}}_{T_2}^0, \tilde{\mathbf{m}}_{T_{1A'}}^0)$.

This set of ground states on the tetrahedron includes multiple states where one of the spins has the same direction. Applying the “Lego-brick” rules, described in Section II E, this means that neighbouring kagome planes can be effectively decoupled in the ground state and there is a ground state degeneracy on the lattice of at least $O(2^L)$.

# J. All-in, all-out order with $A_2$ symmetry

In the Section III B—III F we have explicitly discussed the different types of ordered, classical ground state which occur for $J_3 < 0, J_4 = 0$. For a more general choice of parameters, with $J_4 > 0$ or $J_3 > 0$, it is possible to find situations where the lowest parameter in $H_{ex}^{\text{tet}}$ [Eq. (21)] is $a_{A_2}$. In this case, the ground state will have 4-sublattice order with a finite value of the order parameter $m_{A_2}$ [Table III]. As can be seen from the definition of the order parameter, this type of order is particularly simple, with all spins aligned along the local [111] axes, and all spins pointing either into, or out of, tetrahedra on the A-sublattice. This type of order is commonly referred to as “all-in, all-out” and is observed in some pyrochlore magnets, including Nd$_2$Zr$_2$O$_7$ [115].

Since this type of order does not occur for $J_3 < 0, J_4 = 0$, we will not discuss it further in here. However we note that the order parameter $m_{A_2}$ is a scalar, and that finite-temperature phase transitions are therefore expected to fall into the Ising universality class, in the absence of a first order phase transition.

# IV. THEORY OF CLASSICAL AND QUANTUM SPIN WAVE EXCITATIONS

In order to complete the classical phase diagram described in Section III, it is necessary to understand how quantum and/or classical fluctuations select between the one-dimensional manifold of states described by $\mathbf{m}_E$. At low temperatures, this can be accomplished by exploring the way in which spin-wave excitations contribute to the free energy. Knowledge of the spin-wave excitations also makes it possible to make predictions for inelastic neutron scattering, discussed below, and to benchmark the results of the classical Monte Carlo simulations described in Section V.

In what follows, we describe a general theory of classical, and quantum spin-wave excitations about the different ordered states described in Section III. In Section IV A we establish a classical, low-temperature spin-wave expansion, which makes it possible to determine the boundary between the $\Psi_2$ and $\Psi_3$ ground states for classical spins in the limit $T \to 0$ [cf. Fig. 1]. In Section IV C, we develop an equivalent quantum theory, within the linear spin-wave approximation, which allows us to estimate the boundary between the $\Psi_2$ and $\Psi_3$ ground states for quantum spins in the limit $T \to 0$ [cf. Fig. 2]. We find that the high degeneracies at classical phase boundaries, described in Section III, strongly enhance quantum fluctuations, and in some cases eliminate the ordered moment entirely. In Section IV D we show how both classical and quantum spin-wave theories can be used to make predictions for inelastic neutron scattering. Readers interested in the relationship between classical and quantum spin-wave theories are referred to the discussion in Ref. [82].
A. Classical spin-wave expansion

To obtain the low-energy excitations around the or-
dered ground states of $\mathcal{H}_{\text{ex}}$ [Eq. 4] we use a description in
terms of classical spin waves, analogous to that described
in Ref. [83]. We begin by defining a local coordinate
system by introducing a set of orthogonal unit vectors
$\{\mathbf{u}_i, \mathbf{v}_i, \mathbf{w}_i\}$ for each of the four sublattices $i = 0, 1, 2, 3$
[cf. Fig. 5]. The local “z-axis”, $\mathbf{w}_i$, is chosen to be aligned
with the spins in a given four-sublattice ground state

$$
\mathbf{S}_i = S \mathbf{w}_i \quad \forall i \qquad (67)
$$

The remaining unit vectors, **u**<sub>i</sub> and **v**<sub>i</sub>, are only deter-
mined up to a rotation about **w**<sub>i</sub>, and any convenient
choice can be made.

Using this basis, the fluctuations of the spin $\mathbf{S}_{ik}$ on
sublattice $i$ of tetrahedron $k$ can be parameterized as

$$
\begin{align}
\mathbf{S}_{ik} &= \begin{pmatrix} \frac{\sqrt{S} \delta u_{ik}}{\sqrt{S} \delta v_{ik}} \\ \frac{\sqrt{S} \delta v_{ik}}{\sqrt{S^2 - S \delta u_{ik}^2 - S \delta v_{ik}^2}} \end{pmatrix} \nonumber \\
&\approx \begin{pmatrix} \frac{\sqrt{S} \delta u_{ik}}{\sqrt{S} \delta v_{ik}} \\ S - \frac{1}{2} \delta u_{ik}^2 - \frac{1}{2} \delta v_{ik}^2 \end{pmatrix}. \tag{68}
\end{align}
$$

Substituting Eq. (68) into $\mathcal{H}_{\text{ex}}$ [Eq. (4)] we obtain

$$
\begin{align}
\mathcal{H}_{\text{ex}} &= \sum_{\text{tet}} \sum_{k \ i<j} \mathbf{S}_{ik} \cdot \mathbf{J}_{ij} \cdot \mathbf{S}_{jk} \nonumber \\
&= \mathcal{E}_0 + \mathcal{H}_{\text{ex}}^{\text{CSW}} + \dots \tag{69}
\end{align}
$$

where

$$
\mathcal{E}_0 = \frac{NS^2}{4} \sum_{i,j=0}^{3} \mathbf{w}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{w}_j \quad (70)
$$

is the classical ground-state energy of the chosen
4-sublattice state, and

$$
\mathcal{H}_{\text{ex}}^{\text{CSW}} = \frac{S}{2} \sum_k \sum_{i,j=0}^{3} \left[ \begin{aligned}
& -\frac{1}{2}(\delta u_{ik}^2 + \delta u_{jk}^2 + \delta v_{ik}^2 + \delta v_{jk}^2) (\mathbf{w}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{w}_j) \\
& + \delta u_{ik} \delta u_{jk} (\mathbf{u}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{u}_j) + \delta v_{ik} \delta v_{jk} (\mathbf{v}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{v}_j) \\
& + \delta u_{ik} \delta v_{jk} (\mathbf{u}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{v}_j) + \delta v_{ik} \delta u_{jk} (\mathbf{v}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{u}_j)
\end{aligned} \right] \tag{71}
$$

describes the leading effect of (classical) fluctuations
about this state. Performing a Fourier transformation,
we find

$$
\mathcal{H}_{\text{ex}} = \frac{NS^2}{4} \sum_{i,j=0}^{3} \mathbf{w}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{w}_j \\
+ \frac{1}{2} \sum_{\mathbf{q}} \tilde{u}(-\mathbf{q})^T \cdot \mathbf{M}(\mathbf{q}) \cdot \tilde{u}(\mathbf{q}) \quad (72)
$$

Here $\tilde{u}(\mathbf{q})$ is the vector

$$
\tilde{u}(\mathbf{q}) = \left( \delta u_0(\mathbf{q}), \delta u_1(\mathbf{q}), \delta u_2(\mathbf{q}), \delta u_3(\mathbf{q}), \right. \\
\left. \delta v_0(\mathbf{q}), \delta v_1(\mathbf{q}), \delta v_2(\mathbf{q}), \delta v_3(\mathbf{q}) \right)^T, \quad (73)
$$

and $\mathbf{M}(\mathbf{q})$ the 8 × 8 matrix

$$
\mathbf{M}(\mathbf{q}) = 2S \begin{pmatrix} \mathbf{M}^{11}(\mathbf{q}) & \mathbf{M}^{12}(\mathbf{q}) \\ \mathbf{M}^{21}(\mathbf{q}) & \mathbf{M}^{22}(\mathbf{q}) \end{pmatrix} \quad (74)
$$

built from 4 × 4 blocks

$$
\mathbf{M}_{ij}^{11}(\mathbf{q}) = \cos(\mathbf{q} \cdot \mathbf{r}_{ij}) \left( \mathbf{u}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{u}_j - \delta_{ij} \sum_l (\mathbf{w}_l \cdot \mathbf{J}_{lj} \cdot \mathbf{w}_j) \right) \quad (75)
$$

$$
\mathbf{M}_{ij}^{12}(\mathbf{q}) = \mathbf{M}_{ji}^{21}(\mathbf{q}) = \cos(\mathbf{q} \cdot \mathbf{r}_{ij}) (\mathbf{v}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{u}_j) \quad (76)
$$

$$
\mathbf{M}_{ij}^{22}(\mathbf{q}) = \cos(\mathbf{q} \cdot \mathbf{r}_{ij}) \left( \mathbf{v}_i \cdot \mathbf{J}_{ij} \cdot \mathbf{v}_j - \delta_{ij} \sum_l (\mathbf{w}_l \cdot \mathbf{J}_{lj} \cdot \mathbf{w}_j) \right) \quad (77)
$$

where $i, j \in \{0, 1, 2, 3\}$ and $\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i$ [cf. Eq. (A1)].

The matrix $\mathbf{M}(\mathbf{q})$ [Eq. (74)] can be diagonalized by a
suitable orthogonal transformation, $\mathbf{U} = (\mathbf{U}^T)^{-1}$ to give

$$
\mathcal{H}_{\text{ex}}^{\text{CSW}} = \frac{1}{2} \sum_{\mathbf{q}} \sum_{\nu=1}^{8} \kappa_{\nu\mathbf{q}} v_{\nu\mathbf{q}} v_{\nu-\mathbf{q}} \quad (78)
$$

where the eight normal modes of the system are given by

$$
v(\mathbf{q}) = \mathbf{U} \cdot \tilde{u}(\mathbf{q}) \qquad (79)
$$

with associated eigenvalues $\kappa_{\nu}(\mathbf{q})$. Since $\mathcal{H}_{\text{ex}}^{\text{CSW}}$ [Eq. (78)]
is quadratic in $v_{\nu\mathbf{q}}$, the associated partition function can
be calculated exactly

$$
\begin{align}
\mathcal{Z}_{\text{ex}}^{\text{CSW}} &= \left(\frac{1}{\sqrt{2\pi}}\right)^{2N} \exp\left(-\frac{\mathcal{E}_0}{T}\right) \int \left[\prod_{\nu=1}^{8} \prod_{\mathbf{q}} dv_{\nu\mathbf{q}}\right] \nonumber \\
&\quad \exp\left(-\frac{1}{2} \frac{\sum_{\nu=1}^{8} \sum_{\mathbf{q}} \kappa_{\nu\mathbf{q}} v_{\nu\mathbf{q}} v_{\nu-\mathbf{q}}}{T}\right) \nonumber \\
&= \exp\left(-\frac{\mathcal{E}_0}{T}\right) \prod_{\nu=1}^{8} \prod_{\mathbf{q}} \left(\sqrt{\frac{T}{\kappa_{\nu\mathbf{q}}}}\right). \tag{80}
\end{align}
$$

It follows that, for $T \to 0$, the free energy of the system
is given by

$$
\mathcal{F}_{\text{ex}}^{\text{low-T}} = \mathcal{E}_0 + \frac{T}{2} \sum_{\nu \mathbf{q}} \ln \kappa_{\nu \mathbf{q}} - NT \ln T + \mathcal{O}(T^2). \quad (81)
$$
Where the $O(T^2)$ corrections arise from the higher order, spin wave interaction, terms neglected in Eq. (69).

Within this classical, low-T expansion, the eigenvalues $\kappa_{\nu}(\mathbf{q})$ correspond to independent, low-energy modes, which determine the physical properties of the states, and have the interpretation of a classical spin-wave spectrum. However the classical spectrum $\kappa_{\nu}(\mathbf{q})$ should not be confused with the quantum spin-wave dispersion $\omega_{\nu}(\mathbf{q})$, measured in inelastic neutron scattering experiments and discussed in Section IV C.

## B. Ground-state selection within the one-dimensional manifold of states with E symmetry

Energy alone does not select between the one-dimensional manifold of states with E symmetry$^{25,84}$. However quantum fluctuations$^{26,27,45}$, thermal fluctuations at low temperature$^{25,84}$, thermal fluctuations near the ordering temperature$^{28,78}$, structural disorder$^{75,76}$ and structural distortion$^{75}$ are effective in selecting an ordered ground state. In what follows, we use knowledge of the free energy within a classical spin-wave theory, $\mathcal{F}_{\text{ex}}^{\text{low-T}}$ [Eq. (81)], to determine which of the possible E-symmetry ground states is selected by thermal fluctuations in the limit $T \to 0$. A parallel treatment of the quantum problem is given in Ref. [45] and a classical analysis applied to the limiting case of pure XY spins (i.e. with infinite easy plane anisotropy) is given in Ref. [80].

As a first step, it is helpful to write down a minimal, symmetry allowed, form for the free energy in terms of the components of

$$ \mathbf{m}_E = \begin{pmatrix} m_E \cos(\theta_E) \\ m_E \sin(\theta_E) \end{pmatrix} \qquad (82) $$

(cf. [79]). Keeping only those terms which respect the lattice symmetries, and going up to 6th order in $m_E$

$$ \mathcal{F}_E = \mathcal{F}_0 + \frac{1}{2} a \, m_E^2 + \frac{1}{4} b \, m_E^4 + \frac{1}{6} c \, m_E^6 + \frac{1}{6} d \, m_E^6 \cos(6 \theta_E) + \mathcal{O}(m_E^8) \qquad (83) $$

where $\mathcal{F}_0$ is an unimportant constant.

Note that Eq. (83) does not contain the symmetry allowed coupling to $m_{A_2}$:

$$ m_{A_2} m_E^3 \cos(3\theta_E) $$

which appears in [79]. This is for two reasons: (i) we are considering the $T \to 0^+$ limit where a finite value of $m_{A_2}$ is energetically unfavourable and will be very small; (ii) $m_{A_2}$ can, in any case, be integrated out to arrive at Eq. (83) with renormalized coefficients for the 6th order terms.

It follows from Eq. (83) that i) a suitable order parameter for symmetry breaking within this manifold is $c_E = \cos 6\theta_E$ [77,78], cf. Eq. (40), and that ii) the two states spanning $\mathbf{m}_E$, $\Psi_2$ and $\Psi_3$, are distinguished only

at sixth-order (and higher) in $m_E$ [42]. The quantity $c_E$ is a secondary order parameter, in the sense that a finite expectation value of $c_E$ is induced by coupling to the primary order parameter $\mathbf{m}_E$. These facts have important consequences for the finite temperature phase transition into the paramagnet, as discussed below.

For $T \to 0$, we can parameterise $\mathcal{F}_E$ [Eq. (83)] from $\mathcal{F}_{\text{ex}}^{\text{low-T}}$ [Eq. (81)]. Since $\mathcal{H}_{\text{ex}}^{[\text{Td}]}$ [Eq. (21)] is quadratic in $\mathbf{m}_E$, all other terms in the free energy must be of purely entropic origin. Moreover, symmetry requires that the entropy associated with the E-symmetry states will vary as

$$ S_E(\theta_E) = N \sum_{n=0,1,2,\dots} s_n \cos(6n\theta_E) \qquad (84) $$

The sign of the coefficients $s_n$ then determines the ground state selected by fluctuations. Taking the derivative of Eq. (81) with respect to $T$ allows us to explicitly calculate $S_E(\theta_E)$

$$ \frac{S_E(\theta_E)}{N} = -\frac{1}{N} \frac{\partial \mathcal{F}_{\text{ex}}^{\text{low-T}}}{\partial T} = \ln T + 1 - \frac{1}{2N} \sum_{\mathbf{q}} \ln(\det(\mathbf{M}_{\theta_E}(\mathbf{q}))) \qquad (85) $$

[cf. Ref. [83]], where $\mathbf{M}_{\theta_E}(\mathbf{q})$ is the 8 × 8 matrix defined in Eq. (74), calculated by expanding around a state with a particular value of $\theta_E$. These results are illustrated in Fig. 14. Equivalent calculations, carried out numerically for all parameters associated with E-symmetry ground states, lead to the phase boundary between $\Psi_2$ and $\Psi_3$ shown in Fig. 1. For parameters appropriate to Er$_2$Ti$_2$O$_7$ [26], we find that fluctuations select a $\Psi_2$ ground state, in keeping with published work$^{26,27,42,84}$.

We can now learn more about how ground state selection works by realising that, for some choices of parameters, the operation connecting different E-symmetry ground states becomes an exact symmetry of the Hamiltonian. This is most easily seen in a coordinate frame tied to the local [111] axis, as described in Section II C. Considering $\mathcal{H}_{\text{ex}}^{\text{local}}$ [Eq. (8)], for the simple choice of parameters

$$ (J_{zz}, J_{\pm}, J_{\pm\pm}, J_{z\pm}) = (0, J, 0, 0) \qquad J > 0 $$

the ground state belongs to E and the Hamiltonian reduces to that of an “XY” ferromagnet. In this case the entire one-dimensional manifold of E-symmetry states are connected by an explicit symmetry of the Hamiltonian (rotation around the local $\langle 111 \rangle$ axes). It follows that order-by-disorder is ineffective, and the ground state retains its U(1) symmetry — for a related discussion, see [45].

To gain insight into the phase diagram for $J_3 < 0$, $J_4 \equiv 0$ [cf. Fig. 1], we expand about a point in parameter space

$$ \begin{array}{l} (J_{zz}, J_{\pm}, J_{\pm\pm}, J_{z\pm}) = (-2J, J, 0, 0) \qquad J > 0 \\ \implies (J_1, J_2, J_3, J_4) = (2J, -2J, 0, 0). \end{array} $$
<figure><img src="image_15.png" /><figcaption>FIG. 14: Variation of entropy per spin within the one-dimensional manifold of states with symmetry E. Entropy $S(\theta_E)$ has been estimated using the low-temperature expansion [Eq. (85)], for a range of values of $J_2$, with the entropy of the $\Psi_3$ state subtracted as a reference, i.e. $\Delta s_{\theta_E} = [S(\theta_E) - S(\pi/6)]/N$. The parameters $J_1 = 0.115$meV and $J_3 = -0.099$meV were fixed at values appropriate to Er$_2$Ti$_2$O$_7$ [26], setting $J_4 = 0$. In all cases, $\Delta s_{\theta_E}$ repeats with period $\pi/3$. For a choice of $J_2$ appropriate to Er$_2$Ti$_2$O$_7$ [$J_2/J_1 = -0.49$ — solid purple line], entropy takes on its maximum value for $\theta_E = \frac{n\pi}{3}$, with $n = 0, 1, 2, 3, 4, 5$, corresponding to the six $\Psi_2$ ground states. The extreme variation in entropy at the boundary of the Palmer-Chalker phase [$J_1 = J_2$ — dashed black line], reflects the presence of an $O(L^2)$ set of zero modes in the spectrum of $\Psi_2$ ground state. None the less, the entropy difference between $\Psi_2$ and $\Psi_3$, $\Delta s_{\pi/3} \approx 0.18$ remains finite. For sufficiently negative $J_2$ (dashed blue line, dotted green line) $\Delta s_{\pi/3} < 0$, and fluctuations select the $\Psi_3$ state. All results have been calculated from Eq. (85), with the sum evaluated numerically by a Monte Carlo method. Statistical errors are smaller than the point size.</figcaption></figure>

At this point the ground state manifold is formed from linear combinations of $E$ and $A_2$ symmetry states and the entire ground state manifold is connected by an exact symmetry of the Hamiltonian, so once again there is no order by disorder. For $J_3 < 0$, states with a finite value of $m_{A_2}$ are removed from the ground state manifold and fluctuations select a ground state from amongst the $E$ states. It follows that, for $J_3 \to 0^-$, $J_4 \equiv 0$, the phase boundary between the $\Psi_2$ and $\Psi_3$ states should tend to the line $J_2/|J_3| = -J_1/|J_3|$ [cf. Fig. 1].

To see which phase is preferred for finite $J_3$, we expand the difference in entropy $S_E(\theta_E)$ between the $\Psi_2$ and $\Psi_3$ ground states

$$ \Delta s_{\pi/3} = \frac{S_E(\pi/3) - S_E(\pi/6)}{N} \qquad (86) $$

in powers of $J_{\pm\pm}$ and $J_{z\pm}$. We do this by writing the matrix $M(q)$ [Eq. 74] as

$$ M(q) = M_0(q) + \epsilon X(q) \qquad (87) $$

where $M_0(q)$ is the matrix associated with the high-symmetry point, and $X(q)$ that associated with the perturbation, and noting that

$$ \ln(\det(\mathbf{M}_0 + \epsilon \mathbf{X})) = \ln(\det(\mathbf{M}_0)) + \sum_{n=1}^{\infty} (-1)^{(n+1)} \frac{\epsilon^n}{n} \mathrm{Tr} \left[ (\mathbf{X} \cdot \mathbf{M}_0^{-1})^n \right]. \qquad (88) $$

We then expand in powers of $J_{\pm\pm}$ and $J_{z\pm}$.

We find that the leading correction to $\Delta s_{\pi/3}$ is

$$ \Delta s_{\pi/3} \approx a \left( \frac{J_{\pm\pm}}{J_{\pm}} \right)^3 \qquad (89) $$

where $a = 0.0045$. It follows that, for sufficiently small $J_3$, the phase boundary between $\Psi_2$ and $\Psi_3$ should tend to the line $J_{\pm\pm} = 0$, with the $\Psi_2$ phase favoured for $J_{\pm\pm} > 0$ and $\Psi_3$ favoured for $J_{\pm\pm} < 0$. Numerical evaluation of Eq. (85), in the limit $J_3 \to 0$, yields results in agreement with these arguments [cf. Fig. 1].

On the line $J_{\pm\pm} = 0$ itself, we find that the leading correction to the difference in entropy is

$$ \Delta s_{\pi/3} \approx b \left( \frac{J_{z\pm}}{J_{\pm}} \right)^6 \qquad (90) $$

with $b = -5.3 \times 10^{-5}$. Hence the $\Psi_3$ state is weakly preferred, and the phase boundary will bend towards positive $J_2/|J_3|$, as observed in Fig. 1 of the main text. Since $J_{z\pm}$ is a term which drives ferromagnetic out of plane fluctuations, a negative sign for $b$ is consistent with the argument that $\Psi_3$ is better connected to the ferromagnetic phase, and hence has a softer spectrum for ferromagnetic out-of-plane fluctuations. On the other hand, symmetries of the $\Psi_2$ states have been shown to allow for small antiferromagnetic out-of-plane ordering, of the “all in, all out” type$^{79}$.

In the limit $|J_3| \gtrsim (|J_1|, |J_2|)$, numerical evaluation of Eq. (85) yields the more complex, reentrant behaviour, as
seen in Fig. 1. This behaviour occurs over a very narrow region of parameter space, and is discussed in detail in Ref. [45] for the case of quantum, as opposed to thermal, order by disorder.

### C. Quantum spin-wave theory

Quantum spin-wave theories for $\mathbf{q} = 0$, 4-sublattice classical ground states of $\mathcal{H}_{\text{ex}}$ [Eq. (1)] have been discussed by a number of authors, with attention focused on comparison with inelastic neutron scattering in applied magnetic field$^{26,43}$, and the way in which quantum fluctuations select between the one-dimensional manifold of states with E symmetry$^{26,27,45}$. To date, all calculations have been carried out in the linear spin-wave approximation

$$ \mathcal{H}_{\text{ex}} \approx \mathcal{E}_0 \left(1 + \frac{1}{S}\right) + \sum_{\mathbf{q}} \sum_{\nu=0}^{3} \omega_{\nu}(\mathbf{q}) \left(b_{\nu}^{\dagger}(\mathbf{q})b_{\nu}(\mathbf{q}) + \frac{1}{2}\right) + \dots \quad (91) $$

where $\mathcal{E}_0$ is the classical ground state energy defined in Eq. (70), $\omega_{\nu}(\mathbf{q})$ the spin-wave dispersion, and $b_{\nu}$ a set of non-interacting bosons describing spin-wave modes with band index $\nu = 0, 1, 2, 3$.

Since the anisotropic exchange model $\mathcal{H}_{\text{ex}}$ [Eq. (1)] does not, in general, possess any continuous symmetry, spin-wave excitations about a 4-sublattice classical ground states will generally be gapped, and quantum effects are small. However the enlarged ground-state manifolds occurring where different symmetry ground states meet can lead to “accidental” degeneracies in the spin-wave spectrum, and large quantum fluctuations about the ordered state. The effect of these fluctuations on classical order may be estimated by calculating the correction to the ordered moment on sublattice $i$, $\langle S_i^w \rangle$ defined in Eq. B1. Details of calculations are given in Appendix B. In all of the 4-sublattice phases described in this text, $\langle S_i^w \rangle$ is the same for all sublattices $i = 0, 1, 2, 3$.

In Fig. 15 we show the effect of quantum fluctuations on the classical, zero-temperature ground states of $\mathcal{H}_{\text{ex}}$ [Eq. (1)]. For parameters which are deep within the ordered phases, the ordered moment approaches its full classical value. However the enlarged ground-state degeneracies on classical phase boundaries lead to additional zero-modes in the spin-wave spectrum, and correspondingly larger corrections. Corrections to the ordered moment diverge logarithmically approaching the boundary with the Palmer-Chalker phase from $\Psi_2$, where there are entire planes of zero modes, and approaching the boundary with the ferromagnetic phases from the Palmer-Chalker phase.

In Fig. 2 [Section I], we show the quantum phase diagram of $\mathcal{H}_{\text{ex}}$ [Eq. (1)], within linear spin-wave theory.

Regions where quantum fluctuations eliminate the ordered moment entirely, are shaded white. The effect is strongest where the degeneracy of the classical ground state is highest, i.e. approaching the Heisenberg line

$$ J_1/|J_3| = J_2/|J_3| \to \infty, $$

and in the vicinity of the special point$^{47}$

$$ J_1/|J_3|, J_2/|J_3| \to 0. $$

The absence of an ordered moment within linear spin-wave theory can indicate a region where conventional magnetic order breaks down entirely, and typically underestimates the extent of any unconventional order$^{85-88}$. It therefore seems reasonable to suggest that the vicinity of these phase boundaries will be favorable places to find novel quantum ground states and quantum spin liquids. We will return to this point in our discussion of Er$_2$Sn$_2$O$_7$, in Section IX.

### D. Cross section in neutron scattering

Inelastic neutron scattering experiments measure the dynamical structure factor

$$ S(\mathbf{q}, \omega) = \sum_{\alpha, \beta=1}^{3} \sum_{i,j=0}^{3} \left( \delta_{\alpha\beta} - \frac{q_{\alpha}q_{\beta}}{q^2} \right) \langle m_{\alpha}^{i}(-\mathbf{q}, -\omega) m_{\beta}^{j}(\mathbf{q}, \omega) \rangle \quad (92) $$

where the projection operator

$$ \left( \delta_{\alpha\beta} - \frac{q_{\alpha}q_{\beta}}{q^2} \right) $$

reflects the fact the neutron interacts with the components of the spin transverse to the momentum transfer $\mathbf{q}$, and

$$ m_{\alpha}^{i}(\mathbf{q}) = \frac{1}{\sqrt{2\pi}} \sqrt{\frac{4}{N}} \sum_{\beta=1}^{3} g_{i}^{\alpha\beta} \int dt \left( \sum_{\mathbf{R}_{i}} e^{-i\omega t} e^{i\mathbf{q} \cdot \mathbf{R}_{i}} S_{i}^{\beta}(\mathbf{R}_{i}, t) \right) \quad (93) $$

is the Fourier transform of the magnetic moment associated with the rare-earth ions, for a given sublattice $i = 0, 1, 2, 3$. The associated g-tensor $g_i^{\alpha\beta}$ is defined in Appendix A.

The equal-time structure factor measured in energy-integrated, quasi-elastic neutron scattering, is given by

$$ S(\mathbf{q}, t=0) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} d\omega \, S(\mathbf{q}, \omega) \\ = \sum_{\alpha, \beta=1}^{3} \sum_{i,j=0}^{3} \left( \delta_{\alpha\beta} - \frac{q_{\alpha}q_{\beta}}{q^2} \right) \\ \times \langle m_{\alpha}^{i}(-\mathbf{q}, t=0) m_{\beta}^{j}(\mathbf{q}, t=0) \rangle. \quad (94) $$
<figure><img src="image_16.png" /><figcaption>FIG. 15: Fraction of the classical moment achieved in ordered phases of a pyrochlore magnet with anisotropic exchange interactions, within linear spin-wave theory. Away from the phase boundaries, the ordered moment is close to its full classical value. All results are obtained within a linear spin-wave analysis of $H_{ex}$ [Eq. (1)] — cf. Appendix B — with $J_1/|J_3| = 3 \cos \theta$, $J_2/|J_3| = 3 \cos \theta$, $J_3 < 0$, $J_4 \equiv 0$, corresponding to the white circle in Fig. 1 and Fig. 2.</figcaption></figure>

For many purposes it is also convenient to resolve the equal-time structure factor $S(\mathbf{q}) = S(\mathbf{q}, t = 0)$ into spin-flip (SF) and non-spin flip (NSF) components, or comparison with experiments carried out using polarised neutrons. For neutrons with polarisation along a direction $\hat{\mathbf{n}} \perp \mathbf{q}$, these are given by

$$S^{\text{NSF}}(\mathbf{q}) = \sum_{i,j=0}^{3} \langle (\mathbf{m}^i(-\mathbf{q}) \cdot \hat{\mathbf{n}})(\mathbf{m}^j(\mathbf{q}) \cdot \hat{\mathbf{n}}) \rangle \quad (95)$$

$$S^{\text{SF}}(\mathbf{q}) = \sum_{i,j=0}^{3} \frac{1}{q^2} \langle (\mathbf{m}^i(-\mathbf{q}) \cdot (\hat{\mathbf{n}} \times \mathbf{q})) \times (\mathbf{m}^j(\mathbf{q}) \cdot (\hat{\mathbf{n}} \times \mathbf{q})) \rangle \quad (96)$$

In this article, where we quote results for SF and NSF components of $S(\mathbf{q})$, we follow the conventions of Fennell et al. [89] and consider $\hat{\mathbf{n}} = (1, -1, 0)/\sqrt{2}$.

Connection to theory is made by using simulation, or spin-wave theory to evaluate the correlations of the magnetic moments $m_{\alpha}^i(\mathbf{q}, t)$. Equal-time correlations $S(\mathbf{q})$ [Eqs. (94-96)] can be calculated directly using classical Monte Carlo simulation described in Section V, or using the classical spin wave theory described in Section IV A. The quantum spin-wave theory described in Section IV C and Appendix B, gives access to the full dynamical structure factor $S(\mathbf{q}, \omega)$ [Eq. (92)], as measured by inelastic neutron scattering.

In the case of the classical spin-wave theory discussed in Section IV A, equal time correlations $\langle m_{\alpha}^i(-\mathbf{q})m_{\beta}^j(\mathbf{q}) \rangle$ are expressed in terms of the correlations of the spin wave modes $v_{\nu}(\mathbf{q})$ using Eqs. (68) and (79). The required

correlation functions $\langle v_{\nu\mathbf{q}}v_{\lambda-\mathbf{q}} \rangle$ can then be obtained directly from Eq. (80)

$$\langle v_{\nu\mathbf{q}} v_{\lambda-\mathbf{q}} \rangle = \delta_{\nu\lambda} \frac{T}{\kappa_{\nu\mathbf{q}}}. \quad (97)$$

In Figs. 20, 21, 22, explicit comparison is made between the equal-time structure factor $S(\mathbf{q})$ calculated within classical Monte Carlo simulation, and from a low-temperature classical spin wave theory. We find excellent, quantitative agreement between the two approaches. This confirms classical spin-wave theory as a useful link between the exact analytical zero-temperature theory developed in Sections II and III and the finite-temperature simulations presented in Section V.

In the case of the quantum spin-wave theory described in Section IV C and Appendix B, it is necessary to re-express the spin correlation functions $\langle S_i^{\alpha}S_j^{\beta} \rangle$ in terms of the spin-wave operators $b_{\nu}^{\dagger}$, $b_{\nu}$ [cf. Eq. (91)] This can be accomplished using the Bogoliubov transformations described in Appendix B. Within linear spin-wave theory, the correlations of $b_{\nu}$ take on a simple form

$$\langle b_{\nu}^{\dagger}(\mathbf{q}, \omega) b_{\nu'}^{\dagger}(-\mathbf{q}, -\omega) \rangle = \langle b_{\nu}(\mathbf{q}, \omega) b_{\nu'}(-\mathbf{q}, -\omega) \rangle = 0 \quad (98)$$

$$\begin{aligned} \langle b_{\nu}(\mathbf{q}, \omega) b_{\nu'}^{\dagger}(\mathbf{q}, \omega) \rangle & = \delta_{\nu\nu'} \delta(\omega - \omega_{\nu}(\mathbf{q})) + \langle b_{\nu}^{\dagger}(\mathbf{q}, \omega) b_{\nu'}(\mathbf{q}, \omega) \rangle \\ & = \delta_{\nu\nu'} \delta(\omega - \omega_{\nu}(\mathbf{q}))(1 + n_B(\omega_{\nu}(\mathbf{q}))) \end{aligned} \quad (99)$$

where $\delta(\omega - \omega_{\nu}(\mathbf{q}))$ is the Dirac delta function enforcing conservation of energy and $n_B(\omega)$ is the Bose-Einstein
distribution

$$n_B(\omega) = \frac{1}{\exp\left(\frac{\omega}{T}\right) - 1} \quad (100)$$

## V. FINITE-TEMPERATURE PROPERTIES

The symmetry analysis and “Lego-brick” rules described in Section II, analysis of ground-state energy described in Section III, and spin-wave theory described in Section IV, together make it possible to determine the ordered ground-states of the anisotropic exchange model $H_{ex}$ [Eq. (1)], in the limit $T \to 0$. The resulting classical ground-state phase diagram is shown in phase diagram Fig. 1. However the interesting and unusual properties of rare-earth pyrochlore oxides all come from experiments carried out at finite temperature, and often relate to paramagnetic, rather than ordered phases. We have therefore used classical Monte Carlo simulations to explore the physics of $H_{ex}$ at finite temperature. The main conclusions of these simulations are summarized in the finite-temperature phase diagram Fig. 3.

We note that, although the systems we wish to describe are fundamentally quantum in nature, the classical simulations can be expected to give a qualitatively correct description of the physics in the high temperature, paramagnetic phase, where thermal fluctuations destroy quantum coherence. The classical description should also work, at least in some respects, at low temperature as long as the ground state of the quantum system is in a classically ordered phase. Ultimately, the best justification for our approach will come a posteriori from agreement with neutron scattering experiments [cf. Fig. 24].

In what follows we document the nature of phase transitions from the paramagnet into each of the ordered phases [Sections V A—V D], and explore how the enlarged ground-state degeneracies at classical phase boundaries manifest themselves at finite temperature [Section VI]. Technical details of simulations are given in Appendix C.

### A. Finite-temperature transition from the paramagnet into the Palmer-Chalker phase

The most revealing feature of any broken-symmetry state is usually its finite-temperature phase transition. In Fig. 16 we show simulation results for the finite-temperature phase transition from the paramagnet into Palmer-Chalker phase with $T_2$ symmetry. Simulations were carried out for parameters

$$(J_1, J_2, J_3, J_4) = (0, 0.3, -0.1, 0) \text{ meV}$$

deep within the Palmer-Chalker phase. Clear evidence for a phase transition can be found in the anomalies in both the specific heat $c_h(T)$ [Fig. 16(a)], and

<figure><img src="image_17.png" /><figcaption>FIG. 16: Finite-temperature phase transition from the paramagnet into the Palmer-Chalker phase [$\Psi_4$], as determined by classical Monte Carlo simulation of $H_{ex}$ [Eq. (1)], for parameters $J_1 = 0$ meV, $J_2 = 0.3$ meV, $J_3 = -0.1$ meV, $J_4 = 0$ meV. a) Temperature dependence of the specific heat $c_h(T)$. b) Temperature dependence of the order-parameter susceptibility, $\chi_{T_2}(T)$. c) Temperature dependence of the order parameter, $|m_{T_2}(T)|$. d) Probability distribution of the energy $E$ evaluated at the transition temperature for a cluster of size $L = 12$. The black dashed line in (a)-(c) indicates a first-order phase transition at $T_{T_2} = 305 \pm 5$ mK. Simulations were performed for clusters of $N = 16L^3$ spins, with $L = 4, 6, 8, 12$.</figcaption></figure>

the order-parameter susceptibility $\chi_{T_2}(T)$ [Fig. 16(b)] at $T_{T_2} = 305 \pm 5$ mK.

Symmetry permits a continuous phase transition between the paramagnet and the Palmer-Chalker phase. However for this parameter set, fluctuations drive the transition first order, as is evident from the discontinuity in the value of the order parameter $m_{T_2}$ for $T = T_{T_2}$ [Fig. 16(c)], and the double peak in the probability dis-
tribution for the energy [Fig. 16(d)].

## B. Transition from the paramagnet into the non-colinear ferromagnetic phase

In Fig. 17 we show simulation results for the finite-temperature phase transition from the paramagnet into the non-colinear ferromagnet, for parameters appropriate to Yb₂Ti₂O₇ [43], setting $J_4 = 0$

$$ (J_1, J_2, J_3, J_4) = (-0.09, -0.22, -0.29, 0) \text{ meV} $$

Anomalies in both the specific heat $c_h(T)$ [Fig. 17(a)] and order-parameter susceptibility $\chi_{T_1}(T)$ [Fig. 17(b)] at $T_{T_1} = 455 \pm 5$ mK, provide clear evidence of a phase transition.

At low temperatures, the temperature-dependence of the order parameters $m_{T_{1,A}}$ and $m_{T_{1,B}}$ [Fig. 17(c)-(d)] converges on the values expected from the zero-temperature analysis of Section III B, and with a slope predicted by a low-temperature expansion about the FM ground state (not shown).

The single peak in the probability distribution for the energy [Fig. 17(e)] suggests that for these parameters, the thermal phase transition from paramagnet to non-collinear FM in a classical model is probably continuous and at most very weakly first order.

## C. Transition from the paramagnet into the $\Psi_2$ phase

In Fig. 18 we show simulation results for the finite-temperature phase transition from the paramagnet into the $\Psi_2$ phase, for parameters appropriate to Er₂Ti₂O₇ [26], setting $J_4 = 0$

$$ (J_1, J_2, J_3, J_4) = (0.11, -0.06, -0.1, 0) \text{ meV} $$

This shows a number of interesting features.

Anomalies in both the specific heat $c_h(T)$ [Fig. 18(a)] and order-parameter susceptibility $\chi_E(T)$ [Fig. 18(b)] at $T_E = 505 \pm 5$ mK offer clear evidence of a phase transition.

Both the smooth evolution of the primary order parameter, $m_E$ [Fig. 18(c)], and the single peak in the probability distribution for the energy [Fig. 18(e)] suggests that the phase transition seen in simulation is at most weakly first-order. For the clusters simulated, we confirm that it is possible to obtain a fairly good collapse of data for $\chi_E(T)$ [Fig. 18(b)] using 3D XY exponents$^{78}$.

However there are only a discrete number of $\Psi_2$ ground states, and a finite value of $|m_E|$ alone does not imply $\Psi_2$ order. Evidence for the $\Psi_2$ ground state comes from the secondary order parameter $c_E = \cos 6\theta_E > 0$ [Fig. 18(d)]. Here simulation results are strongly size-dependent, but suggest a slow crossover into the $\Psi_2$ state, occurring at

<figure><img src="image_18.png" /><figcaption>FIG. 17: Finite-temperature phase transition from the paramagnet into the non-collinear ferromagnet (FM), as determined by classical Monte Carlo simulation of $H_{ex}$ [Eq. (1)], for parameters appropriate to Yb₂Ti₂O₇$^{43}$, i.e. $J_1 = -0.09$ meV, $J_2 = -0.22$ meV, $J_3 = -0.29$ meV setting $J_4 = 0$ meV. a) Temperature dependence of the specific heat $c_h(T)$. b) Temperature dependence of the order-parameter susceptibility, $\chi_{T_1}(T)$. c) Temperature dependence of the order parameter, $|m_{T_{1,B}}(T)|$. d) Temperature dependence of the order parameter, $|m_{T_{1,A}}(T)|$. e) Probability distribution of the energy $E$ evaluated at the transition for a system of size $L = 12$. The black dashed line in (a)-(d) indicates a phase transition at $T_{T_1} = 455 \pm 5$ mK. Simulations were performed for clusters of $N = 16L^3$ spins, with $L = 4, 6, 8, 12$.</figcaption></figure>
a $T^* \ll T_E$, without any accompanying feature in $c_h(T)$ [Fig. 18(a)].

On the basis of Landau theory we anticipate that any finite value of $m_E = |\mathbf{m}_E|$ will induce degeneracy breaking in $\theta_E$, and that the entropic selection within the E manifold should therefore occur concurrently with the onset of magnetic order. Depending on the sign of the relevant coupling,

$$ \delta\mathcal{F}_E = \frac{1}{6} d m_E^6 \cos 6\theta_E \quad (101) $$

the system will then enter either a $\Psi_2$ or a $\Psi_3$ ground state.

However, the free-energy barrier separating the $\Psi_2$ and $\Psi_3$ ground states is very small, and this in turn sets a very large length-scale for the selection of the $\Psi_2$ ground state. Based on the low-temperature expansion $\mathcal{F}_{\text{ex}}^{\text{low-T}}$ [Eq. (81)], we estimate that clusters of $N \sim 10^9$ sites should be able to clearly resolve which state is selected at the magnetic ordering temperature.

### D. Transition from the paramagnet into the $\Psi_3$ phase

In Fig. 19 we show simulation results for the finite-temperature phase transition from the paramagnet into the $\Psi_3$ phase, for parameters

$$ (J_1, J_2, J_3, J_4) = (0, -0.3, -0.1, 0) \text{ meV} $$

close to the border with the non-collinear ferromagnet. Anomalies in both the specific heat $c_h(T)$ [Fig. 19(a)] and order-parameter susceptibility $\chi_E(T)$ [Fig. 19(b)] at $T_E = 395 \pm 5$ meV offer clear evidence of a phase transition. Both the smooth evolution of the primary order parameter, $\mathbf{m}_E$ [Fig. 19(c)], and the single peak in the probability distribution for the energy [Fig. 19(e)] suggest that this phase transition is continuous.

Evidence for the $\Psi_3$ ground state comes from the finite value of the secondary order parameter $c_E = \cos 6\theta_E < 0$ [Fig. 19(d)]. This secondary order parameter shows only a slow onset, consistent with a crossover into the $\Psi_3$ state, and is very strongly size-dependent. As with the $\Psi_2$ state considered above, we infer that, with increasing system size, the temperature associated with this crossover scales towards $T = T_N$, and that in the thermodynamic limit, a single phase transition takes place from the paramagnet into the $\Psi_3$ state.

### E. Comparison between Monte Carlo simulation and classical spin wave theory

Here, to demonstrate the validity of our theory, we compare the structure factors, as calculated from the classical spin wave theory $\mathcal{H}_{\text{ex}}^{\text{CSW}}$ [Eq. (78)] and Monte Carlo simulation, for three different parameter sets: the

<figure><img src="image_19.png" /><figcaption>FIG. 18: Finite-temperature phase transition from the paramagnet into the non-coplanar antiferromagnet $\Psi_2$, as determined by classical Monte Carlo simulation of $\mathcal{H}_{\text{ex}}$ [Eq. (1)], for parameters appropriate to Er$_2$Ti$_2$O$_7^{26}$, i.e. $J_1 = 0.11$ meV, $J_2 = -0.06$ meV, $J_3 = -0.1$ meV setting $J_4 = 0$ meV. a) Temperature dependence of the specific heat $c_h(T)$. b) Temperature dependence of the order-parameter susceptibility, $\chi_E(T)$. c) Temperature dependence of the order parameter, $|\mathbf{m}_E(T)|$. d) Temperature dependence of the secondary order parameter, $\cos 6\theta_E$. e) Probability distribution of the energy $E$ evaluated at the transition temperature for a system of size $L = 12$. The black dashed line indicates a continuous phase transition at $T_N = 505 \pm 5$ meV. Simulations were performed for clusters of $N = 16L^3$ spins, with $L = 4, 6, 8, 12$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 19: Finite-temperature phase transition from the paramagnet into the coplanar antiferromagnet $\Psi_3$, as determined by classical Monte Carlo simulation of $H_{ex}$ [Eq. (1)], for parameters $J_1 = 0$ meV, $J_2 = -0.3$ meV, $J_3 = -0.1$ meV, $J_4 = 0$ meV. a) Temperature dependence of the specific heat $c_h(T)$. b) Temperature dependence of the order-parameter susceptibility, $\chi_E(T)$. c) Temperature dependence of the order parameter, $|\mathbf{m}_E(T)|$. d) Temperature dependence of the secondary order parameter, $\cos 6\theta_E$. e) Probability distribution of the energy $E$ evaluated at the transition temperature for a system of size $L = 12$. The black dashed line indicates a continuous phase transition at $T_E = 395 \pm 5$ mK. Simulations were performed for clusters of $N = 16L^3$ spins, with $L = 2, 4, 6, 12$.</figcaption></figure>

parameters of $Yb_2Ti_2O_7$ as found in Ref. [43] where the classical ground state is ferromagnetic, the parameters of $Er_2Ti_2O_7$ as found in Ref. [26] where we expect the order by disorder mechanism to favour the $\Psi_2$ states and one set of parameters where the order by disorder mechanism favours the $\Psi_3$ states. We find excellent, quantitative agreement between the two methods.

In Fig. 20 we show the structure factor $S(q)$ calculated both from classical spin wave theory and from Monte Carlo simulation at $T = 0.05$K, in the NSF, SF and total scattering channels (see Eqs. (94), (95) and (96)). We have used the experimentally determined parameters for the g-tensor$^{90}$ $g_z = 1.77$, $g_{xy} = 4.18$ and exchange integrals$^{43}$ $J_1 = -0.09$meV, $J_2 = -0.22$meV and $J_3 = -0.29$meV, setting $J_4 = 0$. The excellent, quantitative agreement between spin wave theory and simulation demonstrates the excellent equilibration of the simulations down to 0.05K for the parameters of $Yb_2Ti_2O_7$.

$S(q)$ is also useful for studying the entropic ground state selection within the one-dimensional manifold of state with $E$ symmetry. For a given set of parameters we may compare the diffuse scattering calculated in spin wave theory in expansions around the $\Psi_3$ and $\Psi_2$ phases with the diffuse scattering calculated in simulations.

Such a comparison is shown in Fig. 21, for exchange parameters appropriate to $Er_2Ti_2O_7$ ($J_1 = 0.11$ meV, $J_2 = -0.06$ meV and $J_3 = -0.10$ meV, setting $J_4 = 0$) and temperature $T = 0.36$K. From the entropy calculations shown in Fig. 14 we expect the $\Psi_2$ state to be preferred for these values of the exchange parameters. Comparison of the distribution of weight in the vicinity of the (1, 1, 1), (3, 3, 3) and (1, 1, 3) reciprocal lattice vectors between the Monte Carlo data and the spin wave expansions around the $\Psi_2$ and $\Psi_3$ phases supports this conclusion.

Similarly, in Fig. 22 we show a comparison of the diffuse scattering between Monte Carlo simulations and spin wave expansions around the $\Psi_2$ and $\Psi_3$ phases for exchange parameters approaching the non-collinear ferromagnetic phase ($J_1 = 0$, $J_2 = -1.0$ meV and $J_3 = -0.1$ meV, $J_4 = 0$), at $T = 0.4$K. Calculations of the entropy within spin wave theory show that the $\Psi_3$ state should be preferred by fluctuations for these parameters, and this is confirmed by the comparison of the structure factors, in particular by the presence of bright rods along the $\langle 111 \rangle$ directions.
<figure><img src="image_21.png" /><figcaption>FIG. 20: Comparison between results for equal-time structure factor $S(q)$ obtained in classical Monte Carlo (MC) simulation and classical low-temperature expansion (spin wave theory) for parameters appropriate to $Yb_2Ti_2O_7^{43}$. a) Total scattering in the $(h, h, l)$ plane within MC simulation. b) Associated scattering in the non spin-flip (NSF) channel. c) Associated scattering in the spin-flip (SF) channel. d) Total scattering in the $(h, h, l)$ plane within a spin-wave expansion about the ferromagnetic ground state. e) Associated scattering in the NSF channel. f) Associated scattering in the SF channel. Rods of scattering in the $(h, h, h)$ direction, associated with a low-energy spin-wave excitation [cf. Fig. 25], are visible in both SF and NSF channels. All results were obtained at $T = 0.05$ K, for exchange parameters $J_1 = -0.09$meV, $J_2 = -0.22$meV, $J_3 = -0.29$meV, setting $J_4 = 0$. SF and NSF channels are defined with respect to a neutron with polarisation in the $(1, -1, 0)$ direction, as in Ref. 89. $S(q)$ has been calculated using the experimentally measured g-tensor for $Yb_2Ti_2O_7^{43,90}$, with $g_z = 1.77$, $g_{xy} = 4.18$. In order to avoid saturating the colour scale, the intensity associated with Bragg peaks at reciprocal lattice vectors has been subtracted.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIG. 21: Comparison between results for equal-time structure factor $S(q)$ obtained in classical Monte Carlo (MC) simulation and low-temperature expansion (classical spin wave theory) for parameters appropriate to $Er_2Ti_2O_7$. a) Total scattering in the $(h, h, l)$ plane within MC simulation. b) Associated scattering in the non spin-flip (NSF) channel. c) Associated scattering in the spin-flip (SF) channel. d) Total scattering in the $(h, h, l)$ plane within a spin-wave expansion about a $\Psi_2$ ground state. e) Associated scattering in the NSF channel. f) Associated scattering in the SF channel. g) Total scattering in the $(h, h, l)$ plane within a spin-wave expansion about a $\Psi_3$ ground state. h) Associated scattering in the NSF channel. i) Associated scattering in the SF channel. Careful comparison of the distribution of scattering in the vicinity of the (1, 1, 1), (3, 3, 3) and (1, 1, 3) reciprocal lattice vectors supports the conclusion that the $\Psi_2$ state is preferred for these exchange parameters, in agreement with experiment and the calculations described in the text. All results were obtained at $T = 0.36$ K, for exchange parameters $J_1 = 0.11$ meV, $J_2 = -0.06$ meV, $J_3 = -0.10$ meV, setting $J_4 \equiv 0$ and g-tensor parameters $g_z = 2.45$ and $g_{xy} = 5.97$ [26]. For clarity, intensity associated with Bragg peaks at reciprocal lattice vectors has been subtracted.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIG. 22: Comparison between results for equal-time structure factor $S(q)$ obtained in classical Monte Carlo (MC) simulation and low-temperature expansion (classical spin wave theory) in the ordered phase for parameters $J_1 = 0$, $J_2 = -1.0$ meV, $J_3 = -0.10$ meV, $J_4 \equiv 0$, approaching the non-collinear ferromagnet from within the E-symmetry phase. a) Total scattering in the $(h, h, l)$ plane within MC simulation showing strong rod-like features in [111] directions. b) Associated scattering in the non spin-flip (NSF) channel. c) Associated scattering in the spin-flip (SF) channel. d) Total scattering in the $(h, h, l)$ plane within a classical spin-wave expansion about a $\Psi_3$ ground state, showing strong rod-like features in [111] directions. e) Associated scattering in the NSF channel. f) Associated scattering in the SF channel. g) Total scattering in the $(h, h, l)$ plane within a classical spin-wave expansion about a $\Psi_2$ ground state. h) Associated scattering in the NSF channel. i) Associated scattering in the SF channel. Comparison of the scattering supports the conclusion that the $\Psi_3$ ground state is found in simulation, in agreement with the results of the low-T expansion [cf. Fig. 14]. An isotropic g-tensor $g_z = 1$, $g_{xy} = 1$ has been assumed. For clarity, intensity associated with Bragg peaks at reciprocal lattice vectors has been subtracted.</figcaption></figure>
## VI. LIVING ON THE EDGE: THE INFLUENCE OF GROUND STATE MANIFOLDS ON FINITE-TEMPERATURE PHASE TRANSITIONS

The major assertion of this article is that many of the interesting properties of pyrochlore magnets — for example the rods of scattering observed in Yb$_2$Ti$_2$O$_7$, and the order-by-disorder selection of a Ψ$_2$ ground state in Er$_2$Ti$_2$O$_7$, see Sections VII and VIII — are the direct consequence of competition between different ordered phases, and in particular, of the high ground-state degeneracy where phases with different symmetry meet. While the arguments for enlarged ground state manifolds at T = 0 are easy to understand, it is far less obvious that this degeneracy should make itself felt at finite temperature, especially where it is not protected by symmetry.

We can test the internal consistency of these ideas by using the probability distribution of the order parameter

$$ \mathbf{m}_E = m_E (\cos \theta_E, \sin \theta_E) $$

[cf. Eq. (38)] to deconstruct the order-by-disorder selection of Ψ$_2$ and Ψ$_3$ ground states in finite-temperatures simulations of $H_{ex}$ [Eq. (4)]. The probability density function $P(\mathbf{m}_E)$ is sensitive both to the formation of a one-dimensional manifold of states with E symmetry — which manifests itself as a ring in $P(\mathbf{m}_E)$ — and to the selection of an ordered ground state within this manifold — which will appear as six degenerate maxima within the ring.

$P(\mathbf{m}_E)$ also enables us to study the evolution of the ground state manifolds at the boundaries between phases with competing symmetry — in this case either with $T_2$ or with $T_{1,A'}$. At these phase boundaries, $\mathbf{m}_E$ takes on a new, constrained set of values, characteristic of the way in which different manifolds connect. For example Eqs. (53–55) predicts that, on the boundary with the Palmer-Chalker phase, the one-dimensional manifold of states with $|\mathbf{m}_E| = 1$ acquires “spokes” in the directions

$$ \theta_E = \left\{ 0, \frac{\pi}{3}, \frac{2\pi}{3}, \pi, \frac{4\pi}{3}, \frac{5\pi}{3} \right\} $$

connecting $\mathbf{m}_E = 0$ with the six Ψ$_2$ ground states. Observation of such a “spoked wheel” pattern in $P(\mathbf{m}_E)$ at finite temperature would therefore confirm that the zero-temperature degeneracies were still operative.

In an exactly parallel manner, we find that on the boundary with the ferromagnetic phase it is possible to deform the ground state continuously from Ψ$_3$ to a corresponding ferromagnetic state. For this reason Ψ$_3$ should be favoured approaching the boundary with the ferromagnet [cf. Fig. 1].

In Fig. 23 we present results for $P(\mathbf{m}_E)$ and $S(\mathbf{q})$ taken from simulations of $H_{ex}$ [Eq. (4)] for three sets of parameters

$$ (A) \quad (J_1, J_2, J_3, J_4) = (0, -0.3, -0.1, 0) \text{ meV} $$

where we expect a Ψ$_3$ ground state, while approaching the border with the non-collinear FM [Fig. 23 (a, d, g, j)];

$$ (B) \quad (J_1, J_2, J_3, J_4) = (0.11, 0.06, -0.1, 0) \text{ meV} $$

where we expect a Ψ$_2$ ground state, approaching the border with the Palmer-Chalker phase [Fig. 23 (b, e, h, k)]; and

$$ (C) \quad (J_1, J_2, J_3, J_4) = (0.11, 0.11, -0.1, 0) \text{ meV} $$

exactly on the T = 0 border of the Palmer-Chalker phase [Fig. 23 (c, f, i, l)].

The results for $S(\mathbf{q})$ shown in Fig. 23(a-c), demonstrate the diffuse structure expected in the paramagnet in each case: (A) Fig. 23(a) — rods of scattering, reminiscent of those observed in Yb$_2$Ti$_2$O$_7$ [20–24]; (B) Fig. 23(b) — a diffuse web of rings, reminiscent to that observed in experiments on Er$_2$Ti$_2$O$_7$ [see Fig. 24(a-d) for more details], also ordering in Ψ$_2$; (C) Fig. 23(c) — “bow-tie” patterns reminiscent of the pinch points observed in the Heisenberg antiferromagnet on a pyrochlore lattice [44]. Indeed, the Heisenberg antiferromagnet corresponds to the parameters of figure 23(c) with $J_3 = 0$.

The results for $P(\mathbf{m}_E)$ strongly validate our understanding of the problem in terms of degenerate ground state manifolds, even at finite temperatures. On the border of the Palmer-Chalker phase (C) $P(\mathbf{m}_E)$ shows a diffuse spoked wheel both above $T_c$ [Fig. 23(f)] and at $T_c$ [Fig. 23(i)], confirming that the connection implied by the T = 0 ground state manifold survives even for temperatures above the phase transition. The remnant of this spoked wheel pattern is even seen at $T_c$ when we tune away from the phase boundary (parameter set (B)) [Fig. 23(h)]. This spoked wheel pattern in $P(\mathbf{m}_E)$ is strong evidence that the manifolds continue to be operative at $T_c$ and play a role in ordered state selection. This perspective is consistent with the Ginzburg-Landau theory presented in Ref. 79. And it gives new insight into how the terms in the Ginzburg-Landau theory which lift the U(1) degeneracy are generated — namely by fluctuations into these low-energy manifolds.

The reason why we see the spoked wheel pattern at $T_c = 0.26$ K for parameter set (B) and not at $T_c = 0.39$ K for parameter set (A) is probably a consequence of the strong finite size dependence of the entropic selection between Ψ$_2$ and Ψ$_3$.

To conclude, we should add that where these two phase boundaries approach one another, the soft modes associated with the two different sets of manifolds compete. This leads to the complicated, reentrant behaviour seen in Fig. 1, and studied for quantum spins in [45].
<figure><img src="image_24.png" /><figcaption>FIG. 23: Influence of ground-state degeneracy on finite-temperature phase transitions, as revealed by the probability distribution of the order parameter $\mathbf{m}_E = m_E (\cos \theta_E, \sin \theta_E)$ [Eq. (38)]. Results are taken from simulation of $H_{ex}$ [Eq. (4)], with three different sets of exchange parameters, which we label as parameter sets (A), (B) and (C). Parameter set (A), used to calculate panels (a, d, g, j), corresponds to a $\Psi_3$ ground state, approaching the non-collinear FM with $T_c = 0.39$ K; Parameter set (B), used to calculate panels (b, e, h, k), corresponds to a $\Psi_2$ ground state, with $T_c = 0.26$ K; Parameter set (C), used to calculate panels (c, f, i, l), corresponds to a $\Psi_2$ ground state, on the border of the Palmer-Chalker phase with $T_c = 0.065$ K. (a)-(c) quasi-elastic scattering $S(\mathbf{q})$ in the paramagnetic phase $T > T_c$. (d)-(f) corresponding results for the probability density function, $P(\mathbf{m}_E)$. (g)-(i) $P(\mathbf{m}_E)$ at the transition temperature $T = T_c$. (j)-(l) $P(\mathbf{m}_E)$ in the ordered phase $T < T_c$. Parameter set (A): For a finite-size system, the onset of $\Psi_3$ occurs progressively, through (g) the emergence of a one-dimensional manifold of states with finite $|\mathbf{m}_E|$, and then (j) the entropic selection of $\theta_E$ corresponding to one of six distinct $\Psi_3$ ground states. (a) The connection with the non-collinear FM is evident in $S(\mathbf{q})$, with rods of scattering strongly reminiscent of those seen in Yb$_2$Ti$_2$O$_7$. Parameter set (B): The same process occurs, but in this case $P(\mathbf{m}_E)$ shows that $\Psi_2$ ground states are favoured at low temperatures (k) and even at the transition (h). Parameter set (C): On the boundary of the Palmer-Chalker phase, the ground state manifold includes additional manifolds of states which mix $\mathbf{m}_E$ and $\mathbf{m}_{T_2}$. These are evident (f,i) in the “spoked wheel” seen in $P(\mathbf{m}_E)$ at $T \ge T_c$, and drive the entropic selection of the $\Psi_2$ ground state. (c) The high degeneracy at this phase boundary is also evident in the “bow-tie” structure in $S(\mathbf{q})$. Further details of simulations and the parameters corresponding to (A), (B) and (C) are given in the text.</figcaption></figure>
<figure><img src="image_25.png" /></figure>

FIG. 24: Correlations in the high-temperature paramagnetic phase, as revealed by the quasi-elastic structure factor $S(q)$. (a-c) results for parameters interpolating from (a) Er$_2$Ti$_2$O$_7$ [26] to (c) the boundary of the Palmer-Chalker phase ($\Psi_4$). The diffuse scattering characteristic of the $\Psi_2$ phase evolve into sharp features reminiscent of pinch points when bordering the $\Psi_4$ phase. Results are taken from classical Monte Carlo simulations carried out for (a) $J_2 = -0.06$ meV, $T = 750$ mK; (b) $J_2 = 0.06$ meV, $T = 390$ mK; (c) $J_2 = 0.11$ meV, $T = 100$ mK. In all cases, $J_1 = -0.11$ meV, $J_3 = -0.1$ meV, $J_4 \equiv 0$, and $S(q)$ has been calculated using g-tensor parameters appropriate to Er$_2$Ti$_2$O$_7$ [26]. (d) detail of $S(q)$ for parameters appropriate to Er$_2$Ti$_2$O$_7$ at $T = 616$ mK, plotted with a colour scale chosen to match Fig. 14 of [91] and with the same temperature ratio $T/T_c = 1.22$ ($T_c = 505$ mK in simulations). (e-g) results for parameters interpolating from Yb$_2$Ti$_2$O$_7$ [cf. Ref. (43)], to the border of the $\Psi_3$ phase. The rods of scattering along [111] directions, interpreted as evidence of dimensional reduction in Yb$_2$Ti$_2$O$_7$ [22], evolve into weakly-dispersing, low-energy excitations in the neighbouring $\Psi_3$ phase. Results are taken from classical Monte Carlo simulations of $H_{ex}$ [Eq. (1)] for (e) $J_1 = -0.09$ meV, $T = 750$ mK; (f) $J_1 = -0.04$ meV, $T = 400$ mK; (g) $J_1 = -0.0288$ meV, $T = 450$ mK. In all cases, $J_2 = -0.22$ meV, $J_3 = -0.29$ meV, $J_4 \equiv 0$, and $S(q)$ has been calculated using g-tensor parameters appropriate to Yb$_2$Ti$_2$O$_7$ [90].

## VII. APPLICATION TO Er$_2$Ti$_2$O$_7$

Early heat capacity measurement of Er$_2$Ti$_2$O$_7$ revealed a phase transition at $T_c = 1.25$ K, releasing an entropy $\Delta s \approx 0.97k_B \ln 2$ per spin, consistent with the ordering of the ground state doublet of Er [18]. Later, neutron scattering studies revealed the nature of the low temperature order, finding it to correspond to the $\Psi_2$ configurations illustrated in Fig. 10 [25,26,64].

The selection of the $\Psi_2$ states in Er$_2$Ti$_2$O$_7$ out of the 1D manifold of states transforming with $E$ symmetry has been identified as a textbook example of order-by-disorder$^{25}$, with quantum zero-point fluctuations$^{26,27,45}$, low-temperature classical thermal fluctuations$^{80}$ and thermal fluctuations near the ordering temperature$^{28,78}$ all favouring $\Psi_2$ order. A corollary of this conclusion is that there should be a small, fluctuation induced, gap at $q = 0$ in the spin wave spectrum. And, such a gap has now been observed in inelastic neutron scattering$^{92}$. We note that, an alternative scenario has been proposed in which the selection of $\Psi_2$ comes instead from virtual

fluctuations into higher crystal field levels$^{29,42,93}$. Thus it may in fact be that the selection of $\Psi_2$ has multiple contributions- both from harmonic fluctuations of the ground state order and from virtual crystal field fluctuations.

Here, using the exchange parameters for Er$_2$Ti$_2$O$_7$ taken from [26], we confirm that thermal fluctuations of classical Heisenberg spins select a $\Psi_2$ phase at finite temperature. Estimates from our Monte Carlo simulations give $T_c \approx 500$ mK [Fig. 18], somewhat lower than both experiments and those obtained in a high temperature series expansion of the quantum spin model$^{28}$. But within the paramagnetic phase, our simulations of the spin structure factor $S(q)$ are in excellent agreement with neutron scattering measurements$^{91}$ showing the build-up of long-range order, as can be seen from the comparison between Fig. 24(d) of this work and Fig. 14 of Ref. [91].

A question which remains is *why* fluctuations should favour $\Psi_2$ in the case of Er$_2$Ti$_2$O$_7$ and more generally why they should favour either $\Psi_2$ or $\Psi_3$ for a given set of exchange parameters $\{J_i\}$. Our work provides the answer
<figure><img src="image_26.png" /><figcaption>FIG. 25: Spin-wave dispersion calculated within a classical, low-temperature expansion, showing dimensional reduction of a subset of excitations. (a) Excitations of the FM ground state, for exchange parameters appropriate to Yb$_2$Ti$_2$O$_7$, i.e. $J_1 = -0.09$meV, $J_2 = -0.22$meV, $J_3 = -0.29$meV, setting $J_4 = 0$. The ferromagnet possesses a flat band in the $(h, h, h)$ ($\Gamma \to L$) direction at energy $\Delta \approx 0.22$ meV, which gives rise to rods in the equal time structure factor (cf. Fig. 20). (b) Excitations of the $\Psi_3$ ground state, for exchange parameters on the boundary between the $\Psi_3$ and FM phases, i.e. $J_1 = -0.029$meV. $J_2 = -0.22$meV, $J_3 = -0.29$meV with $J_4 = 0$. The $\Psi_3$ phase, on the phase boundary also possesses a quasi-flat band along $(h, h, h)$, which in this case is gapless at the $\Gamma$ point of the Brillouin zone. This leads us to suggest that the low-energy rod-like features observed in the paramagnetic phase of Yb$_2$Ti$_2$O$_7$ arise from its proximity in parameter space to the $\Psi_3$ phase and the low-energy modes which are present on the phase boundary (see discussion in Section VIII).</figcaption></figure>

to this question and, in so doing, underlines how the properties of a frustrated magnet are strongly influenced by competing phases.

The mechanism by which the $\Psi_2$ states are selected is inherited from the phase boundary with the neighbouring Palmer-Chalker phase. At this boundary three additional continuous sets of ground states appear connecting the 6 Palmer-Chalker ground states, to the 1D manifold of E symmetry states [Fig. 13]. The points in configuration space at which these sets of ground states meet are none other than the $\Psi_2$ configurations. Due to their favoured position at the junctions of the ground state manifold the $\Psi_2$ states gain additional soft modes and are selected by fluctuations in the region approaching the boundary with the Palmer-Chalker phase. The consequences of these connected manifolds are visible even in finite temperature simulations as shown in Fig. 23.

It is worth noting that an exactly parallel mechanism selects the $\Psi_3$ states for parameters proximate to the

ferromagnetic phase. In the region proximate to both Palmer-Chalker and Ferromagnetic phases, a complicated re-entrant behaviour is observed$^{45}$ [Fig. 1].

Our work shows that the preference of fluctuations for $\Psi_2$ ordering in Er$_2$Ti$_2$O$_7$ is a property inherited from a nearby phase boundary where the $\Psi_2$ states sit at the junctions of a connected ground state manifold. In the context of this result it is interesting to ask how the spin correlations evolve as the exchange parameters are tuned from those appropriate to Er$_2$Ti$_2$O$_7$ to the boundary of the Palmer-Chalker phase. This is illustrated using Monte Carlo simulations of the spin structure factor $S(q)$ in Fig. 24(a-c). For the parameters appropriate to Er$_2$Ti$_2$O$_7$, our simulations reproduce the smooth features observed in experiment. As the phase boundary is approached these smooth features evolve into sharp, pinch-point like features, associated with the large ground state manifold on the phase boundary. In the limit $J_3 \to 0$ these pinch-point like features become the pinch points
<figure><img src="image_27.png" /><figcaption>FIG. 26: Prediction for the spin-wave excitations of the ferromagnetically-ordered ground state of Yb$_2$Ti$_2$O$_7$, calculated within linear spin-wave theory for the parameters given by Ross et al. [43]. (a) Spin-wave dispersion for a single ferromagnetic domain with magnetisation parallel to [100], showing a minimum gap to excitations $\Delta \sim 0.3$ meV, occurring in the zone centre. (b) Spin-wave dispersion averaged over the 6 possible ferromagnetic domains. (c) Path in reciprocal space used in making plots. Calculations were carried out for the anisotropic exchange model $\mathcal{H}_{\text{ex}}$ [Eq. (1)], with results convoluted with a Gaussian of full width at half maximum 0.014 meV to mimic finite experimental resolution. The relative intensity of scattering is shown in false color. Details of calculations are given in Section IV C and Appendix B.</figcaption></figure>

associated with the Coulomb phase of the O(3) Heisenberg model on the pyrochlore lattice$^{44}$.

Our study of the ground state selection in Er$_2$Ti$_2$O$_7$ emphasises that the properties of frustrated magnets can be strongly influenced by the soft modes appearing on nearby boundaries. This same essential insight also manifests itself – in rather different ways – in the study of Yb$_2$Ti$_2$O$_7$ and Er$_2$Sn$_2$O$_7$, to which we now turn.

VIII. APPLICATION TO Yb$_2$Ti$_2$O$_7$

Like its sister compound Er$_2$Ti$_2$O$_7$, Yb$_2$Ti$_2$O$_7$ was first identified as undergoing a finite temperature ordering transition in the heat capacity study by Blöte et al., nearly 50 years ago$^{18}$. That study revealed a sharp anomaly in the heat capacity at $T_c = 0.214$K, with a corresponding release of entropy of $\Delta s = 0.97k_B \ln 2$ per spin. Since then, the presence of this phase transition in stoichiometric Yb$_2$Ti$_2$O$_7$ has been debated in the literature, with different groups, with different samples, reporting differing results for the presence or absence of magnetic order. Nevertheless, it is now widely accepted that stoichiometric Yb$_2$Ti$_2$O$_7$ undergoes a thermodynamic phase transition into a state with finite magnetisation at a temperature $T_c \sim 0.2$K [24,65–67,69–72,94,95].

For the parameters given by Ross et al. [43], the the-

ory developed in Section III predicts that Yb$_2$Ti$_2$O$_7$ has a **q** = 0 ground state, with non-collinear order ferromagnetic order. This “splayed ferromagnet” is consistent with the interpretation of neutron scattering experiments given in [66,69–71]. An ordered ground state of this type would normally be expected to support coherent, dispersing spin-wave excitations, with a finite gap coming from the anisotropy of exchange interactions, as illustrated in Fig. 26.

Curiously, however, gapped, coherent spin-waves have yet to be observed in Yb$_2$Ti$_2$O$_7$, with a succession of experiments reporting a broad, gapless continuum at low temperatures$^{21,69–71}$. The origin of this gapless continuum remains a puzzle, although the presence of competing classical ordered phases must ultimately impact on quantum excitations$^{68}$. It is also important to recall that the ordered ground state breaks only the point-group symmetries of the anisotropic exchange model $\mathcal{H}_{\text{ex}}$ [Eq. (1)], and so spin-waves cannot be interpreted as Goldstone modes. It follows that interaction effects may play an important role, even at low orders in $1/S$. And it is interesting to note that the broad continuum observed in experiment has more in common with semi-classical simulations of the spin-excitation of the paramagnetic phase of Yb$_2$Ti$_2$O$_7$ [71,96], than with the linear spin-wave excitations of the ground state, as shown in Fig. 26. This point will be discussed further elsewhere$^{97}$.
While the nature of the ground state of Yb$_2$Ti$_2$O$_7$ has proved controversial, and the associated excitations remain to be understood, all neutron-scattering experiments agree about the signature feature of its paramagnetic phase — striking “rod”-like structures along the 〈111〉 directions of reciprocal space. First observed more than ten years ago$^{20}$, these rods of scattering have since been interpreted as evidence of dimensional reduction$^{21,22}$ and, in the context of $\mathcal{H}_{\text{ex}}$ [Eq. (1)], as evidence of significant anisotropic exchange interactions$^{23,43,49}$. They are a robust feature of $S(\mathbf{q})$, as calculated from $\mathcal{H}_{\text{ex}}$ [Eq. (1)] within both the (semi-classical) random phase approximation$^{23,24}$, and classical Monte Carlo simulations$^{68}$ [Fig. 24]. However, despite their ubiquity, the origin of these rods of scattering remains mysterious.

To understand the origin of the rods of scattering we must once again look to the influence of the phase boundaries. The classical ground states of $\mathcal{H}_{\text{ex}}$ reduce to a set of independent kagome planes on the boundary between FM and Palmer-Chalker phases, and to a set of independent chains on the boundary between the $\Psi_2$ and Palmer-Chalker phases. However the rods of scattering seen in Yb$_2$Ti$_2$O$_7$ occur for parameters where the ground state of $\mathcal{H}_{\text{ex}}$ is expected to be ordered and fully three-dimensional$^{43}$. Indeed, our classical Monte Carlo simulations predict that Yb$_2$Ti$_2$O$_7$ orders at 450 mK [cf. Fig. 17], a little higher than the $T_c \approx 200$ mK found in experiment.

Within the scenario of multiple-phase competition, rods of scattering can be traced back to dimensionally-reduced excitations, due to quasi-degenerate lines of low-lying spin wave excitations, which evolve into low-lying excitations of $\Psi_3$ on the boundary between the FM and the $\Psi_3$ phases [cf. Fig. 25]. This progression is also clear in the evolution of $S(\mathbf{q})$ from parameters appropriate to Yb$_2$Ti$_2$O$_7$ [Fig. 24(e)] to the border of the $\Psi_3$ phase [Fig. 24(g)], also shown in Ref. [68].

Seen in this light, the observation of rods of scattering in Yb$_2$Ti$_2$O$_7$ is a consequence of the proximity of competing ordered states- in this case the $\Psi_2$ and $\Psi_3$ states. The importance of these competing ground states in driving the unusual physics of Yb$_2$Ti$_2$O$_7$ has been underlined in recent work$^{68,71}$. In Ref. [68] it was shown that both quantum and thermal fluctuations bring the phase boundary between ferromagnetic and E symmetry phases closer to the parameter regime appropriate to Yb$_2$Ti$_2$O$_7$. In this sense, Yb$_2$Ti$_2$O$_7$ really is a material “living on the edge” between differing magnetic orders, and the rods of scattering are a manifestation of this.

## IX. APPLICATION TO Er$_2$Sn$_2$O$_7$

Like Er$_2$Ti$_2$O$_7$ and Yb$_2$Ti$_2$O$_7$, the magnetic ions in Er$_2$Sn$_2$O$_7$ have a Kramers doublet ground state$^{30,98}$, and their interactions are believed to be well-described by $\mathcal{H}_{\text{ex}}$ [Eq. (1)] [34]. Correlations reminiscent of the Palmer-Chalker phase have been observed in neu-

<figure><img src="image_28.png" /><figcaption>FIG. 27: Ordered moment in the region of parameter space relevant to Er$_2$Sn$_2$O$_7$, as calculated in linear spin wave theory. On the Palmer-Chalker ($\Psi_4$) side of the phase boundary, the quantum correction to the ordered moment is small. However, the correction diverges on the approach the phase boundary from the non-coplanar antiferromagnet ($\Psi_2$), indicating the possibility of a region of quantum disorder between these two phases. Experimental estimates of exchange parameters in Er$_2$Sn$_2$O$_7$ [34] place it close to this phase boundary, making it a good candidate for the observation of quantum spin-liquid physics. Calculations were carried out for the anisotropic exchange model $\mathcal{H}_{\text{ex}}$ [Eq. (1)], as described in Section IV C and Appendix B, with parameters $J_2 = 0.08$ meV, $J_3 = -0.11$ meV taken from [34], setting the Dzyaloshinskii-Moriya interaction $J_4 = 0$. The error bars on the estimated value of $J_1$ are taken from Ref. [34].</figcaption></figure>

tron scattering$^{34}$, and magnetization measurements show some evidence of spin-freezing at low temperatures$^{34}$. Nonetheless, Er$_2$Sn$_2$O$_7$ shows no evidence of magnetic order, in thermodynamic measurements$^{33,34}$, $\mu$SR$^{31}$, or neutron scattering$^{33,34}$, down to a temperature of 20 mK [31].

The exchange parameters determined for Er$_2$Sn$_2$O$_7$ in Ref. [34] would place it extremely close to the phase boundary between the Palmer-Chalker and $\Psi_2$ states [Fig. 1]. Classical Monte Carlo simulations with this parameter set predict a phase transition into the Palmer-Chalker state at $T_c \approx 200$ mK. However, we can once more gain further insight by looking at the behaviour of the model $\mathcal{H}_{\text{ex}}$ approaching the phase boundary.

As the phase boundary is approached, the ground-state value of the ordered moment, as calculated in linear spin wave theory, is reduced by quantum fluctuations. This is illustrated in Fig. 27. Approaching the boundary from the Palmer-Chalker side, this quantum correction is small. However, approaching the boundary from the $\Psi_2$ side, the correction is logarithmically divergent. Since spin wave theory typically underestimates quantum effects, this divergence is a likely indicator of a region of quantum disorder between the Palmer-Chalker and $\Psi_2$
region of the phase diagram [cf. Fig. 2]. The placement of Er₂Sn₂O₇ immediately adjacent to this classical phase boundary thus makes it a prime candidate for the observation of quantum spin liquid physics.

We can see therefore that the competition between Palmer-Chalker and Ψ₂ ordering in Er₂Sn₂O₇ enhances quantum fluctuations in that material and may even stabilize a quantum disordered state. Such a scenario would be consistent with the lack of observed magnetic order in Er₂Sn₂O₇ and would make Er₂Sn₂O₇ the first example of a pyrochlore spin liquid with dominantly XY-like interactions. While the recent years have seen considerable theoretical advances in the understanding of quantum spin liquid states occurring close to the Ising (spin ice) limits of $H_{ex}$ [5–11], the limit of dominant XY interactions has been much less explored for quantum spins. A deeper understanding of Er₂Sn₂O₇ calls for further work in this direction.

## X. OTHER PYROCHLORE MAGNETS

The family of rare-earth pyrochlore oxides R₂M₂O₇ is very diverse¹⁷, and it is becoming even more so, with high-pressure synthesis techniques allowing for many new combinations of rare earth, R, and transition metal, M to be realized⁹⁹. In this Section we briefly explore the properties of further examples of pyrochlore magnets. We will mostly restrict our comments to materials based on the Kramers ions such as Er³⁺, Yb³⁺, where the model $H_{ex}$ [Eq. (1)], with couplings of the form of $J_{01}$ [Eq. (6)], offers a completely a completely general description of nearest-neighbor interactions on the pyrochlore lattice. That said, many of the same physical phenomena arise in pyrochlore magnets based on non-Kramers ions and, where these have a doublet ground state, interactions may also be described using $H_{ex}$ [8,51]. With this in mind, we also make a few brief comments about materials based on the non-Kramers ions Tb³⁺ and Pr³⁺.

Yb₂Sn₂O₇ [35,36,100] and Yb₂Pt₂O₇ [39] have both been identified as having ferromagnetic ground states, and may therefore be placed in the non-collinear ferromagnet ($T_1$) region of our phase diagram. Meanwhile, the ground state of another Yb based system, Yb₂Ge₂O₇ has been shown in neutron scattering experiments to belong to the manifold of E symmetry states, although any ground state selection between Ψ₂ and Ψ₃ configurations has yet to be determined¹⁰¹,¹⁰². The progression, as a function of decreasing transition metal ionic radius $Sn→Pt→Ti→Ge$ for the Yb₂M₂O₇ compounds thus tunes across the phase boundary between the ferromagnetic and Ψ₂/Ψ₃ regions in Fig. 1 [68,101]. The spin excitations above these ordered states remain a puzzle, however, with a recent systematic study of the Sn, Ti and Ge compounds showing an absence of coherent spin waves in all three materials, and a continuity of the inelastic neutron scattering spectrum across the finite temperature ordering transition in each case⁷⁰.

Amongst Er based pyrochlores, Er₂Ge₂O₇ has been observed to order antiferromagnetically at $T_N = 1.41$K [38]. Neutron scattering experiments reveal this to belong to the Ψ₂/Ψ₃ region of the phase diagram in Fig. 1. The behaviour of the intensity of the magnetic Bragg peaks under external magnetic field suggests that fluctuations select a Ψ₃ ground state out of the E symmetry manifold for Er₂Ge₂O₇ [38]. Er₂Pt₂O₇, also orders antiferromagnetically at $T_N = 0.3$K. If the variation of exchange constants with the size of the transition metal ion M is monotonous, then it would be expected for this material to lie near the boundary between the Ψ₂ and Palmer-Chalker regions of the phase diagram in Fig. 1. However, this assumption could fail in the case of Er₂Pt₂O₇ since Pt⁴⁺ ion differs from {Ti⁴⁺, Sn⁴⁺, Ge⁴⁺} in that it possesses a partially filled d-shell³⁹,¹⁰³.

The Gd based pyrochlores Gd₂M₂O₇ have attracted significant research interest over a period of nearly two decades¹⁷,¹⁰⁴. The physics of these pyrochlores is somewhat different from (e.g.) Yb and Er based systems because the Gd³⁺ ions have vanishing orbital angular momentum $L=0$. The interactions of the $S = 7/2$ Gd spins are thus quite isotropic, and the nearest neighbour anisotropic exchange which is the focus of this article is a less important consideration than further neighbour interactions- including dipole-dipole interactions. Nevertheless, combining the effect of nearest neighbour antiferromagnetic exchange and the nearest neighbour part of the dipole-dipole interaction our theory does predict a Palmer-Chalker ground state which is consistent with observations on Gd₂Sn₂O₇ [105,106].

In Gd₂Ti₂O₇ further neighbour interactions drive a complex phenomenology involving multiple phase transitions and “partially ordered” states¹⁰⁷⁻¹¹⁰. The precise nature of the magnetic ground state remains a matter of discussion in the literature¹¹¹, and is beyond the scope of the present study. Recently, another Gd pyrochlore, Gd₂Pb₂O₇ has been synthesised presenting an antiferromagnetic ordering transition at $T=0.81$K into an as yet unidentified ground state¹¹².

Nd based pyrochlores have attracted considerable recent attention and provide examples of “all-in, all out” ordering on the pyrochlore lattice¹¹³⁻¹¹⁶. In particular, it has recently come to light that Nd₂Zr₂O₇ exhibits the novel phenomenon of “moment fragmentation” in which the spin correlations simultaneously show the pinch points of a Coulomb phase and the Bragg peaks of an ordered state¹¹⁷⁻¹¹⁹. However, we note that the ground doublet of the Nd ions is of the dipolar-octupolar type¹¹⁵,¹¹⁶. In this case, the anisotropic exchange interactions would have a different form to Eq. (6) [52]. Recently, a Ce based pyrochlore Ce₂Sn₂O₇, also believed to belong to the group of dipolar-octupolar pyrochlores, has been synthesised which appears to show an absence of magnetic order down to $T = 0.02$K [53], suggesting it as a promising candidate spin-liquid system.

There are also many pyrochlore systems where the rare-earth ion R is a non-Kramers ion, such as Tb³⁺ or
Pr³⁺ [17]. Many of these systems also exhibit two-fold degenerate crystal field ground states, and a pseudospin-1/2 description of the magnetic degrees of freedom may be appropriate⁵¹,¹²⁰,¹²¹ (although in the Tb³⁺ systems the picture is complicated by the relatively small gap to crystal field excitations¹²²⁻¹²⁴). In this case the nearest neighbour bilinear exchange Hamiltonian for these pseudospins takes the form of Eq. (8), with the constraint that J<sub>z</sub> = 0 [8,51] — cf. Table II. The pseudospin **S** must also be interpreted differently in this case with the part of the pseudospin perpendicular to the local [111] directions corresponding to a quadrupolar degree of freedom. Thus, an “easy plane” order of the pseudospins (as occurs in the **E** and T<sub>2</sub> regions of our phase diagram) actually corresponds to quadrupolar order for non-Kramers ions.

Amongst the Tb pyrochlores, Tb<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub> is known to exhibit ferromagnetic order<sup>73,74,125</sup> and therefore belongs to the T<sub>1</sub> region of the phase diagram. Tb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, meanwhile, has been studied for a long time as a candidate spin-liquid<sup>126,127</sup> exhibiting power-law spin correlations<sup>128,129</sup>. The apparent spin-liquid behaviour has been linked with quantum spin ice physics<sup>55</sup>, but it has recently been proposed that an alternative form of spin-liquid physics may be at work<sup>47</sup>. Recent studies have also revealed the presence of competing ordered states, with quadrupolar order (corresponding to a Palmer-Chalker like configuration of the pseudospins<sup>130,131</sup>) and antiferromagnetic **q** = (1/2, 1/2, 1/2) order<sup>132,133</sup> observed depending on the sample details and the experimental cooling protocol. Tb<sub>2</sub>Ge<sub>2</sub>O<sub>7</sub> has recently come to light as another system of interest for spin liquid physics, with no long range order observed down to 20mK and short range ferromagnetic correlations<sup>134</sup>. This may suggest that Tb<sub>2</sub>Ge<sub>2</sub>O<sub>7</sub> lives in a region of disorder, proximate to the ferromagnetic (T<sub>1</sub>) phase, such as that proposed in Ref. [47].

At the same time, the Pr based pyrochlores, Pr<sub>2</sub>M<sub>2</sub>O<sub>7</sub>
(M=Sn, Zr, Hf, Pb) have emerged as promising can-
didates for quantum spin ice physics, exhibiting an
absence of magnetic order and dynamic, spin-ice like
correlations<sup>112,135-137</sup>.

XI. CONCLUSION

Rare-earth pyrochlore oxides offer a veritable treasure
trove of novel physical phenomena, ranging from classi-
cal and quantum spin liquids, to dimensional reduction,
and phases governed by order-by-disorder effects. In this
Article we have established a general theory of multiple-
phase competition in materials with anisotropic exchange
interactions on the pyrochlore lattice, and shown how
it can be can be applied to three specific materials :
Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> and Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>. The recurring
theme throughout this analysis is of materials “living on
the edge”, in the sense of having properties which are
dictated by the competition between neighbouring forms

of magnetic order.

Starting from a very general model of interactions between nearest-neighbour spins on the pyrochlore lattice, H<sub>ex</sub> [Eq. (1)], we have used an analysis based on point-group symmetry to establish the exact, classical, ground-state phase diagram [Sections II and III]. As a by-product, we provide a complete classification of possible four-sublattice ordered states, according to the way they lift the symmetry of a single tetrahedron. Moreover, using the “Lego-brick” rules developed in Section II E, we are able to identify the conditions under which the classical ground state manifold undergoes a dimensional reduction into independent planes or chains of spins, opening the door to new physical phenomena.

We have given particularly careful consideration to the
ground-state manifolds in the limit where the symmetric
off diagonal exchange J₃ < 0 and the Dzyaloshinskii-
Moriya interaction J₄ = 0 [Section III]. Based on the
experimental parameterisations of exchange interactions
for Kramers pyrochlores²⁶,³⁴,⁴³, this limit is of particular
experimental relevance. We have elucidated the nature
of the expanded ground state manifolds which occur at
the phase boundaries of this model, and it is these which
drive much of the physics of the surrounding regions of
parameter space.

Stepping out of the ground state manifold, we have
given, in Section IV, calculations of the spin wave exci-
tations in the ordered phases. Among other things, this
allows us to determine the ground state selection by both
quantum and thermal fluctuation and to identify regions
of the phase diagram where classical order will be melted
by quantum fluctuations.

We have also studied the finite-temperature properties
of the anisotropic exchange model [Eq. (1)] using clas-
sical Monte Carlo simulations, presented in Sections V
and VI. These simulations make it possible to determine
the finite temperature phase diagram [Fig. (3)] and to
show how the expanded ground state manifolds on the
phase boundaries manifest themselves at finite tempera-
ture [Fig. 23].

The implications of our theory for three specific
pyrochlore materials — Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> and
Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub> — are expounded in Sections VII, VIII and
IX. We find that the influence of nearby phase bound-
aries accounts for the ground state selection by fluctu-
ations in Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, the apparent dimensional reduction
in the paramagnetic phase of Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, and the sup-
pression of magnetic order in Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>. The unusual
properties of these three materials can be understood as
“living on the edge” — having properties controlled by
the competition between different ground states.

As discussed in Section X, the family of rare earth py-
rochlore magnets is a large one, extending well beyond
the three materials covered in Sections VII to IX. In par-
ticular, recent work has seen the synthesis of rare earth
pyrochlores R<sub>2</sub>M<sub>2</sub>O<sub>7</sub> with M =Ge, Pt, Pb, Os, Zr, Hf
[38,39,112,137–139]. This work suggests the possibility
to move around the phase diagram shown in Fig. 1 and
Fig. 2, using chemical or physical pressure. This would be particularly interesting in cases where systems related by a change of transition metal ion M live on opposite sides of a classical phase boundary. One might then hope to tune through a region of strong quantum fluctuations using substitution of the transition metal ion. Such an opportunity would seem to present itself for Er₂M₂O₇ with M =Sn, Ti and Yb₂M₂O₇ with M =Ge, Ti [68].

Our analysis may also be useful in the study of related systems, such as the rare-earth spinel CdEr₂Se₄ [140], where the Er ions also form a pyrochlore lattice. Looking further afield, a modification of our theory could be used in the understanding of “breathing” pyrochlore compounds, where the tetrahedra of the pyrochlore lattice alternate in size¹⁴¹⁻¹⁴⁴. It is also interesting to note recent neutron scattering experiments on NaCaCo₂F₇, a pyrochlore material with quenched exchange disorder¹⁴⁵. The observed diffuse scattering in that material is rather similar to that predicted by our Monte Carlo simulations in Fig. 24(b). This may spring from a connection between the low-energy configurations found in the clean limit of the anisotropic exchange model, studied in this Article, and the low-energy configurations of the disordered system.

From a theoretical perspective, our work also highlights the importance of large, classical, ground-state degeneracies which are not related to the well-studied examples of spin ice, or of the Heisenberg antiferromagnet on a pyrochlore lattice. These degeneracies, which emerge in a number of different limits of $H_{ex}$ [Eq. (1)], could lead to novel forms of classical or quantum spin liquid, as well as entirely new forms of classical and quantum order⁴⁸. One such case, where fluctuations lead to a spin-liquid described by a rank-2 tensor field with a continuous gauge symmetry, has been developed in Ref. [47]. However there are many other regions of parameter space where strong fluctuations persist to low temperature¹⁴⁶, and the majority of these have yet to be fully explored. It seems that the study of rare earth pyrochlore magnets with anisotropic exchange interactions may have many more surprises yet in store.

## Acknowledgments

The authors are pleased to acknowledge helpful conversations with Bruce Gaulin, Michel Gingras, Alannah Hallas, Edwin Kermarrec, Isabelle Mirebeau, Sylvain Pectit, Karlo Penc, and Kate Ross, and a critical reading of the manuscript by Mathieu Taillefumier. This work was supported by the Theory of Quantum Matter Unit of the Okinawa Institute of Science and Technology Graduate University.

## Appendix A: g-tensor in local and global coordinate frames

The local crystal-electric field (CEF), acting on a given magnetic ion, affects both the character of its ground state, and the nature of its exchange interactions with other magnetic ions. For this reason, it is often convenient chose a coordinate frame

$$ \{\mathbf{x}_i^{\text{local}}, \mathbf{y}_i^{\text{local}}, \mathbf{z}_i^{\text{local}}\} $$

which is tied to the local CEF on site i. We can accomplish this by choosing $\mathbf{z}_i^{\text{local}}$ to be parallel with the [111] axis on site i, i.e. the local axis with $C_3$-symmetry

For the tetrahedron shown in Figure 5, the magnetic ions labelled $S_0$, $S_1$, $S_2$ and $S_3$ occupy positions

$$ \begin{aligned} \mathbf{r}_0 &= \frac{a}{8} (1, 1, 1) & \mathbf{r}_1 &= \frac{a}{8} (1, -1, -1) \\ \mathbf{r}_2 &= \frac{a}{8} (-1, 1, -1) & \mathbf{r}_3 &= \frac{a}{8} (-1, -1, 1), \quad (\text{A1}) \end{aligned} $$

relative to the centre of the tetrahedron, in units such that the cubic, 16-site unit cell of the pyrochlore lattice occupies a volume $V = a^3$. The local [111] axes on these sites are given by

$$ \begin{aligned} \mathbf{z}_0^{\text{local}} &= \frac{1}{\sqrt{3}}(1, 1, 1) & \mathbf{z}_1^{\text{local}} &= \frac{1}{\sqrt{3}}(1, -1, -1) \\ \mathbf{z}_2^{\text{local}} &= \frac{1}{\sqrt{3}}(-1, 1, -1) & \mathbf{z}_3^{\text{local}} &= \frac{1}{\sqrt{3}}(-1, -1, 1). \end{aligned} \quad (\text{A2}) $$

In defining $(\mathbf{x}_i^{\text{local}}, \mathbf{y}_i^{\text{local}})$ we follow the conventions of Ross et al. [43], and make the convenient choice

$$ \begin{aligned} \mathbf{x}_0^{\text{local}} &= \frac{1}{\sqrt{6}}(-2, 1, 1) & \mathbf{x}_1^{\text{local}} &= \frac{1}{\sqrt{6}}(-2, -1, -1) \\ \mathbf{x}_2^{\text{local}} &= \frac{1}{\sqrt{6}}(2, 1, -1) & \mathbf{x}_3^{\text{local}} &= \frac{1}{\sqrt{6}}(2, -1, 1), \end{aligned} \quad (\text{A3}) $$

such that all $\mathbf{y}_i^{\text{local}}$ lie in a common plane

$$ \begin{aligned} \mathbf{y}_0^{\text{local}} &= \frac{1}{\sqrt{2}}(0, -1, 1) & \mathbf{y}_1^{\text{local}} &= \frac{1}{\sqrt{2}}(0, 1, -1) \\ \mathbf{y}_2^{\text{local}} &= \frac{1}{\sqrt{2}}(0, -1, -1) & \mathbf{y}_3^{\text{local}} &= \frac{1}{\sqrt{2}}(0, 1, 1). \end{aligned} \quad (\text{A4}) $$

In this local coordinate frame, the magnetic moment

$$ m_i^{\alpha} = \sum_{\beta=1}^{3} g_{\text{local}}^{\alpha\beta} S_i^{\beta} \quad (\text{A5}) $$

is connected to the (pseudo) spin-1/2 operator $S_i^\alpha$ [Eq. (7)], through a g-tensor with a diagonal simple form

$$ \mathbf{g}_{\text{local}} = \begin{pmatrix} g_{xy} & 0 & 0 \\ 0 & g_{xy} & 0 \\ 0 & 0 & g_z \end{pmatrix} \quad (\text{A6}) $$
<table><caption>TABLE VI: Estimates of the components of the g-tensor in the local frame $\mathbf{g}_{\text{local}}$ [Eq. (A6)], taken from experiment on Yb$_{2}$Ti$_{2}$O$_{7}$ [90], Er$_{2}$Ti$_{2}$O$_{7}$ [26], and Er$_{2}$Sn$_{2}$O$_{7}$ [34].</caption><thead><tr><th></th><th>Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub></th><th>Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub></th><th>Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub></th></tr></thead><tbody><tr><th>g<sub>xy</sub></th><td>4.18</td><td>5.97</td><td>7.52</td></tr><tr><th>g<sub>z</sub></th><td>1.77</td><td>2.45</td><td>0.05</td></tr></tbody></table>

where $\alpha, \beta = \{\mathbf{x}_{\text{local}}, \mathbf{y}_{\text{local}}, \mathbf{z}_{\text{local}}\}$, and $\mathbf{g}_{\text{local}}$ is independent of the site considered. Estimates of $g_{xy}$ and $g_z$, taken from experiment on Yb$_{2}$Ti$_{2}$O$_{7}$ [90], Er$_{2}$Ti$_{2}$O$_{7}$ [26], and Er$_{2}$Sn$_{2}$O$_{7}$ [34] are shown in Table VI. For rare-earth ions with Ising character, such as Dy$^{3+}$ in Dy$_{2}$Ti$_{2}$O$_{7}$, $g_z > g_{xy}$, while for the rare-earth ions considered in this paper with easy-plane character, $g_z < g_{xy}$.

The g-tensor in the coordinate frame of the crystal axes, $\mathbf{g}_i$ [Eq. (3)], can be found by rotating $\mathbf{g}_{\text{local}}$ [Eq. (A6)] back into the global coordinate frame $\mu, \nu = \{\mathbf{x}, \mathbf{y}, \mathbf{z}\}$. Since the required rotation depends on the lattice site, the resulting g-tensor is sublattice-dependent

$$ \mathbf{g}_0 = \begin{pmatrix} g_1 & g_2 & g_2 \\ g_2 & g_1 & g_2 \\ g_2 & g_2 & g_1 \end{pmatrix} \qquad \mathbf{g}_1 = \begin{pmatrix} g_1 & -g_2 & -g_2 \\ -g_2 & g_1 & g_2 \\ -g_2 & g_2 & g_1 \end{pmatrix} \\ \mathbf{g}_2 = \begin{pmatrix} g_1 & -g_2 & g_2 \\ -g_2 & g_1 & -g_2 \\ g_2 & -g_2 & g_1 \end{pmatrix} \qquad \mathbf{g}_3 = \begin{pmatrix} g_1 & g_2 & -g_2 \\ g_2 & g_1 & -g_2 \\ -g_2 & -g_2 & g_1 \end{pmatrix} \qquad (\text{A7}) $$

where

$$ g_1 = \frac{2}{3}g_{xy} + \frac{1}{3}g_z \qquad g_2 = -\frac{1}{3}g_{xy} + \frac{1}{3}g_z. \qquad (\text{A8}) $$

## Appendix B: Linear spin-wave theory for a general 4-sublattice ground state

A general framework for linear spin-wave theory on the pyrochlore lattice is set out in [43], following the pattern that can be found in [147]. For completeness here we reproduce the technical steps needed to apply such a theory to the 4-sublattice, $\mathbf{q} = 0$ classical ground states discussed in Section III.

As with the classical spin-wave theory developed in Section IV A, it is convenient to work in a local basis, in which spins are quantised such that their local z-axis is aligned with the classical ground state. Following Eq. (68), we label these local axes

$$ \{\mathbf{u}_i, \mathbf{v}_i, \mathbf{w}_i\} $$

and quantize fluctuations about the classical ground state by introducing Holstein-Primakoff bosons

$$ S_i^w = S - a_i^\dagger a_i \qquad (\text{B1}) $$

$$ S_i^+ = S_i^u + iS_i^v = (2S - a_i^\dagger a_i)^{1/2} a_i \approx \sqrt{2S} a_i \qquad (\text{B2}) $$

$$ S_i^- = S_i^u - iS_i^v = a_i^\dagger (2S - a_i^\dagger a_i)^{1/2} \approx \sqrt{2S} a_i^\dagger \qquad (\text{B3}) $$

where $[a_i, a_j^\dagger] = \delta_{ij}$.

Substituting these expressions in $\mathcal{H}_{\text{ex}}$ [Eq. (4)] and Fourier transforming them, we obtain

$$ \mathcal{H}_{\text{ex}} \approx \mathcal{E}_0 + \mathcal{H}_{\text{ex}}^{\text{LSW}} + \dots \qquad (\text{B4}) $$

where $\mathcal{E}_0$ is the classical ground state energy defined in Eq. (70), and

$$ \mathcal{H}_{\text{ex}}^{\text{LSW}} = \frac{1}{2} \sum_{\mathbf{q}} \tilde{A}^{\dagger}(\mathbf{q}) \cdot \mathbf{X}(\mathbf{q}) \cdot \tilde{A}(\mathbf{q}) \qquad (\text{B5}) $$

describes quantum fluctuations at the level of linear spin wave theory. Here $\tilde{A}^\dagger(\mathbf{q}), \tilde{A}(\mathbf{q})$ are eight-component vectors of operators

$$ \begin{aligned} \tilde{A}^{\dagger}(\mathbf{q}) = & (a_0^{\dagger}(\mathbf{q}), a_1^{\dagger}(\mathbf{q}), a_2^{\dagger}(\mathbf{q}), a_3^{\dagger}(\mathbf{q}), \\ & a_0(-\mathbf{q}), a_1(-\mathbf{q}), a_2(-\mathbf{q}), a_3(-\mathbf{q})) \end{aligned} \qquad (\text{B6}) $$

and $X(\mathbf{q})$ is an $8 \times 8$ matrix written in block form as

$$ \mathbf{X}(\mathbf{q}) = 2S \begin{pmatrix} \mathbf{X}^{11}(\mathbf{q}) & \mathbf{X}^{12}(\mathbf{q}) \\ \mathbf{X}^{21}(\mathbf{q}) & \mathbf{X}^{22}(\mathbf{q}) \end{pmatrix} \qquad (\text{B7}) $$

$$ \mathbf{X}_{ij}^{11}(\mathbf{q}) = \cos(\mathbf{q} \cdot \mathbf{r}_{ij}) \begin{pmatrix} \mathbf{c}_i \cdot \mathbf{J}^{ij} \cdot \mathbf{c}_j^* - \delta_{ij} \sum_l \mathbf{w}_l \cdot \mathbf{J}^{lj} \cdot \mathbf{w}_j \end{pmatrix} \qquad (\text{B8}) $$

$$ \mathbf{X}_{ij}^{12}(\mathbf{q}) = \mathbf{X}_{ji}^{21*} = \cos(\mathbf{q} \cdot \mathbf{r}_{ij}) \begin{pmatrix} \mathbf{c}_i \cdot \mathbf{J}^{ij} \cdot \mathbf{c}_j \end{pmatrix} \qquad (\text{B9}) $$

$$ \mathbf{X}_{ij}^{22}(\mathbf{q}) = \cos(\mathbf{q} \cdot \mathbf{r}_{ij}) \begin{pmatrix} \mathbf{c}_i^* \cdot \mathbf{J}^{ij} \cdot \mathbf{c}_j - \delta_{ij} \sum_l \mathbf{w}_l \cdot \mathbf{J}^{lj} \cdot \mathbf{w}_j \end{pmatrix} \qquad (\text{B10}) $$

where

$$ \mathbf{c}_i = \frac{1}{\sqrt{2}} (\mathbf{u}_i + i\mathbf{v}_i). \qquad (\text{B11}) $$

The spin-wave Hamiltonian $\mathcal{H}_{\text{ex}}^{\text{LSW}}$ [Eq. (B5)] can be diagonalized by a suitable Bogoliubov transformation. We accomplish this following the method outlined in Ref. [148] by introducing new Bose operators $[b_i, b_j^\dagger] = \delta_{ij}$, such that

$$ \begin{aligned} B^{\dagger}(\mathbf{q}) = & (b_0^{\dagger}(\mathbf{q}), b_1^{\dagger}(\mathbf{q}), b_2^{\dagger}(\mathbf{q}), b_3^{\dagger}(\mathbf{q}), \\ & b_0(-\mathbf{q}), b_1(-\mathbf{q}), b_2(-\mathbf{q}), b_3(-\mathbf{q})) \\ = & \tilde{A}^{\dagger}(\mathbf{q}) \cdot \mathbf{U}^{\dagger}(\mathbf{q}) \end{aligned} \qquad (\text{B12}) $$

The condition that these operators are Bosonic may be written as

$$ [B_i(\mathbf{q}), B_j^\dagger(\mathbf{q}')] = \sigma_{ij} \delta_{\mathbf{q}\mathbf{q}'} \qquad (\text{B13}) $$

where

$$ \hat{\sigma} = \begin{pmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} \end{pmatrix}. \qquad (\text{B14}) $$
is an 8 × 8 matrix (written in block form) and leads to a pseudo-unitary condition on $\mathbf{U}(\mathbf{q})$

$$ \mathbf{U}^{-1}(\mathbf{q}) = \hat{\sigma} \cdot \mathbf{U}^{\dagger}(\mathbf{q}) \cdot \hat{\sigma}. \qquad (\text{B15}) $$

Substituting in Eq. (B5), we obtain

$$ \begin{aligned} \mathcal{H}_{\text{ex}}^{\text{LSW}} &= \frac{1}{2} \sum_{\mathbf{q}} B^{\dagger}(\mathbf{q}) \cdot \mathbf{U}^{-1\dagger}(\mathbf{q}) \cdot \mathbf{X}(\mathbf{q}) \cdot \mathbf{U}^{-1}(\mathbf{q}) \cdot B(\mathbf{q}) \\ &= \frac{1}{2} \sum_{\mathbf{q}} B^{\dagger}(\mathbf{q}) \cdot \hat{\sigma} \cdot \mathbf{U}(\mathbf{q}) \cdot \hat{\sigma} \cdot \mathbf{X}(\mathbf{q}) \cdot \mathbf{U}^{-1}(\mathbf{q}) \cdot B(\mathbf{q}). \end{aligned} \qquad (\text{B16}) $$

The object $\mathbf{U}(\mathbf{q}) \cdot \hat{\sigma} \cdot \mathbf{X}(\mathbf{q}) \cdot \mathbf{U}^{-1}(\mathbf{q})$ is a similarity transformation on the matrix $\hat{\sigma} \cdot \mathbf{X}(\mathbf{q})$, and for correctly chosen $\mathbf{U}(\mathbf{q})$, will be a diagonal matrix containing the eigenvalues of $\hat{\sigma} \cdot \mathbf{X}(\mathbf{q})$. We then arrive at

$$ \mathcal{H}_{\text{ex}}^{\text{LSW}} = \frac{1}{2} \sum_{\mathbf{q}} B^{\dagger}(\mathbf{q}) \cdot \hat{\sigma} \cdot \begin{pmatrix} \omega_{\nu}(\mathbf{q}) & 0 \\ 0 & -\omega_{\nu}(\mathbf{q}) \end{pmatrix} \cdot B(\mathbf{q}). \qquad (\text{B17}) $$

Collecting all terms, reordering operators and inserting into Eq. (B4) we obtain the result quoted in Section IV C

$$ \mathcal{H}_{\text{ex}} \approx \mathcal{E}_0 \left(1 + \frac{1}{S}\right) + \sum_{\mathbf{q}} \sum_{\nu=0}^{3} \omega_{\nu}(\mathbf{q}) \left(b_{\nu}^{\dagger}(\mathbf{q})b_{\nu}(\mathbf{q}) + \frac{1}{2}\right) + \dots \qquad (\text{B18}) $$

The dispersion $\omega_{\nu}(\mathbf{q})$ of the four branches of spin waves can be found by numerical diagonalization of $\hat{\sigma} \cdot \vec{X}(\mathbf{q})$.

## Appendix C: Classical Monte Carlo simulation

The Monte Carlo simulations described in this paper are based on the Metropolis algorithm with parallel tempering$^{149,150}$ and over-relaxation$^{151}$. The spins are modelled as classical vectors of length $|S_i| = 1/2$ and locally updated using the standard Marsaglia method$^{152}$. We consider cubic clusters of linear dimension $L$, based on the 16-site cubic unit cell of the pyrochlore lattice, and containing $N = 16L^3$ sites. A Monte Carlo step (MCs) is defined as $N$ attempts to locally update a randomly chosen spin, and $t_{max}$ (measured in MCs) is the total Monte Carlo time over which data are collected.

Equilibration is performed for each temperature in two successive steps. First the system is slowly cooled down from high temperature (random initial spin configuration) to the temperature of measurement $T$ during $t_{max}/10$ MCs. Then, the system is equilibrated at

temperature $T$ during additional $t_{max}/10$ MCs. After equilibration, Monte Carlo time is set to zero and measurements start and go on for $t_{max} \sim 10^5 - 10^7$ MCs.

All thermodynamical observables have been averaged over Monte Carlo time every 10 MCs, except for calculations of the equal-time structure factor $S(\mathbf{q})$, where data points were taken every 100 MCs for efficiency. The parallel tempering method implies simultaneously simulating a large number of replicas of the system in parallel, with each replica held at a different temperature. The program then regularly attempts to swap the spin configurations of replicas with neighbouring temperatures, in such a way as to maintain detailed balance$^{149,150}$. Simulating ~120 replicas, with swaps attempted every 100 MCs appears to offer a good compromise between efficiency and decorrelation for $L = 6$.

In the case of the over-relaxation method, after each Monte Carlo step, two further sweeps are made of the entire lattice. Each spin feels an effective field due to the interaction with its six nearest neighbours; any rotation around this axis conserves the energy and is thus an acceptable move respecting detailed balance. To avoid rotating successive neighbouring spins, we first update all spins of sublattice 0, then sublattice 1, 2 and finally 3. The first iteration of all $N$ spins is deterministic, i.e. we rotate them by the maximum allowed angle; while for the second iteration, a random angle of rotation is chosen for each spin. The generation of so many random numbers is of course time consuming but is recommended for better equilibration$^{153}$. We note that convergence of the specific heat $c_h \to 1$ for $T \to 0$ is a good indication of the equilibration of ordered phases at low temperatures.

The main results of Monte Carlo simulations are summarised in the finite-temperature phase diagram Fig. 3, which spans all four of the ordered phases discussed in the article. This phase diagram was determined from simulations for 64 different parameter sets, equally spaced on the circle defined by $\sqrt{J_1^2 + J_2^2} = 3|J_3|$ illustrated by the white circle in Fig. 1, with $J_3 = -0.1$ meV and $J_4 = 0$. Transition temperatures for each phase were extracted from the relevant order-parameter susceptibilities, as described in Section V.

Simulations were performed for a cluster of $N = 3456$ spins ($L = 6$), and data averaged over 10 independent runs during $t_{max} = 10^6$ MCs. Parallel tempering was used, typically with 121 replicas, at temperatures equally-spaced from 0 to 1.2 K. However, close to the boundaries between phases with different symmetries, the large number of competing ground states makes simulations difficult to equilibrate. Here, additional data points with better statistics were sometimes necessary, typically with 201 temperatures on a smaller temperature window, with $t_{max} = 10^7$ MCs and $N = 8192$ (i.e. $L=8$).
1 L. Balents, *Spin liquids in frustrated magnets*, Nature **464**, 199 (2010).

2 C. Castelnovo, R. Moessner, and S. L. Sondhi, *Spin Ice, Fractionalization, and Topological Order*, Annu. Rev. Cond. Mat. Phys. **3**, 35 (2012).

3 M. Hermele, M. P. A. Fisher, and L. Balents, *Pyrochlore photons: The U(1) spin liquid in a S=1/2 three-dimensional frustrated magnet*, Phys. Rev. B **69**, 064404 (2004).

4 A. Banerjee, S. V. Isakov, K. Damle, and Y.-B. Kim, *Unusual Liquid State of Hard-Core Bosons on the Pyrochlore Lattice*, Phys. Rev. Lett. **100**, 047208 (2008).

5 L. Savary and L. Balents, *Coulombic Quantum Liquids in Spin-1/2 Pyrochlores*, Phys. Rev. Lett. **108**, 037202 (2012).

6 N. Shannon, O. Sikora, F. Pollmann, K. Penc, and P. Fulde, *Quantum Ice: A Quantum Monte Carlo Study*, Phys. Rev. Lett. **108**, 067204 (2012).

7 O. Benton, O. Sikora, and N. Shannon, *Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice*, Phys. Rev. B **86**, 075154 (2012),

8 S.-B. Lee, S. Onoda and L. Balents, *Generic quantum spin ice*, Phys. Rev. B **86**, 104412 (2012).

9 L. Savary and L. Balents, *Spin liquid regimes at nonzero temperature in quantum spin ice*, Phys. Rev. B **87**, 205130 (2013).

10 M. J. P. Gingras and P. A. McClarty, *Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets*, Rep. Prog. Phys. **77**, 056501 (2014).

11 Z. Hao, A. G. R. Day, and M. J. P. Gingras, *Bosonic many-body theory of quantum spin ice*, Phys. Rev. B **90**, 214430 (2014).

12 L. D. Pan, S. K. Kim, A. Ghosh, C. M. Morris, K. A. Ross, E. Kermarrec, B. D. Gaulin, S. M. Koohpayeh, O. Tchernyshyov and N. P. Armitage, *Low-energy electrodynamics of novel spin excitations in the quantum spin ice Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Nat. Commun.* **5**, 4970 (2014).

13 L. D. Pan, N. J. Laurita, K. A. Ross, B. D. Gaulin and N. P. Armitage, *A measure of monopole inertia in the quantum spin ice Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Nature Phys.* **12**, 361-366 (2016).

14 Y. Tokiwa, T. Yamashita, M. Udagawa, S. Kittaka, T. Sakakibara, D. Terazawa, Y. Shimoyama, T. Terashima, Y. Yasui, T. Shibauchi and Y. Matsuda, *Possible observation of highly itinerant quantum magnetic monopoles in the frustrated pyrochlore Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Nat. Commun.* **7**, 10807 (2016).

15 M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, *Geometrical Frustration in the Ferromagnetic Pyrochlore Ho<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Phys. Rev. Lett.* **79**, 2554 (1997).

16 R. Moessner, *Relief and generation of frustration in pyrochlore magnets by single-ion anisotropy*, Phys. Rev. B **57**, R5587(R), (1998).

17 J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, *Magnetic pyrochlore oxides*, Rev. Mod. Phys. **82**, 53 (2010).

18 H. W. J. Blöte, R. F. Wielinga and W. J. Huiskamp, *Heat capacity measurements on rare-earth double oxides R<sub>2</sub>M<sub>2</sub>O<sub>7</sub>, Physica* **43**, 549 (1969).

19 G. W. Chern, *Novel Magnetic Orders and Ice Phases in Frustrated Kondo–Lattice Models*, Spin **5**, 1540006 (2015).

20 P. Bonville, J. A. Hodges, E. Bertin, J.-P. Bouchaud, P. Dalmas de Reotier, L.-P. Regnault, H. M. Ronnow, J.-P. Sanchez, S. Sosin, and A. Yaouanc, *Transitions and Spin Dynamics at Very Low Temperature in the Pyrochlores Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> and Gd<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>, Hyperfine Interact.* **156/157**, 103 (2004).

21 K. A. Ross, J. P. C. Ruff, C. P. Adams, J. S. Gardner, H. A. Dabkowska, Y. Qiu, J. R. D. Copley, and B. D. Gaulin, *Two-Dimensional Kagome Correlations and Field Induced Order in the Ferromagnetic XY Pyrochlore Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Phys. Rev. Lett.* **103**, 227202 (2009).

22 K. A. Ross, L. R. Yaraskavitch, M. Laver, J. S. Gardner, J. A. Quilliam, S. Meng, J. B. Kycia, D. K. Singh, T. Proffen, H. A. Dabkowska, and B. D. Gaulin, *Dimensional evolution of spin correlations in the magnetic pyrochlore Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Phys. Rev. B **84**, 174442 (2011).

23 J. D. Thompson, P. A. McClarty, H. M. Ronnow, L. P. Regnault, A. Sorge, and M. J. P. Gingras, *Rods of Neutron Scattering Intensity in Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>: Compelling Evidence for Significant Anisotropic Exchange in a Magnetic Pyrochlore Oxide*, Phys. Rev. Lett. **106**, 187202 (2011).

24 L.-J. Chang, S. Onoda, Y. Su, Y.-J. Kao, K.-D. Tsuei, Y. Yasui, K. Kakurai, and M. R. Lees, *Higgs transition from a magnetic Coulomb liquid to a ferromagnet in Yb<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Nat. Commun.* **3**, 992 (2012).

25 J. D. M. Champion, M. J. Harris, P. C. W. Holdsworth, A. S. Wills, G. Balakrishnan, S. T. Bramwell, E. Cizmar, T. Fennell, J. S. Gardner, J. Lago, D. F. McMorrow, M. Orendac, A. Orendacova, D. McK. Paul, R. I. Smith, M. T. F. Telling, and A. Wildes, *Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>: Evidence of quantum order by disorder in a frustrated antiferromagnet*, Phys. Rev. B **68**, 020401(R), (2003).

26 L. Savary, K. A. Ross, B. D. Gaulin, J. P. C. Ruff, and L. Balents, *Order by Quantum Disorder in Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Phys. Rev. Lett.* **109**, 167201 (2012).

27 M. E. Zhitomirsky, M. V. Gvozdikova, P. C. W. Holdsworth, and R. Moessner, *Quantum Order by Disorder and Accidental Soft Mode in Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>, Phys. Rev. Lett.* **109**, 077204 (2012).

28 J. Oitmaa, R. R. P. Singh, B. Javanparast, A. G. R. Day, B. V. Bagheri and M. J. P. Gingras, *Phase transition and thermal order-by-disorder in the pyrochlore antiferromagnet Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub>: A high-temperature series expansion study*, Phys. Rev. B **88**, 220404(R), (2013).

29 J. G. Rau, S. Petit and M. J. P. Gingras, *Order by virtual crystal field fluctuations in pyrochlore XY antiferromagnets*, Phys. Rev. **93**, 184408 (2016).

30 K. Matsuhira, Y. Hinatsu, K. Tenya, H. Amitsuka and T. Sakakibara, *Low-Temperature Magnetic Properties of Pyrochlore Stannates*, J. Phys. Soc. Jpn. **71**, 1576, (2002).

31 J. Lago, T. Lancaster, S. J. Blundell, S. T. Bramwell, F. L. Pratt, M. Shirai and C. Baines, *Magnetic ordering and dynamics in the XY pyrochlore antiferromagnet: a muon-spin relaxation study of Er<sub>2</sub>Ti<sub>2</sub>O<sub>7</sub> and Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>, J. Phys.: Condens. Matter **17**, 979, (2005).

32 Masae Shirai, *Experimental Investigations of Frustrated Magnets*, Ph.D. Thesis, University of London, (2007).

33 P. M. Sarte, H. J. Silverstein, B. T. K. van Wyk, J. S. Gardner, Y. Qiu, H. D. Zhou and C. R. Wiebe, *Absence of long-range magnetic ordering in the pyrochlore compound Er<sub>2</sub>Sn<sub>2</sub>O<sub>7</sub>, J. Phys.: Condens. Matter **23**, 382201, (2011).