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

L.D.C. Jaubert, Owen Benton, Jeffrey G. Rau, J. Oitmaa, R.R.P. Singh, Nic

Shannon, Michel J.P. Gingras

► To cite this version:

L.D.C. Jaubert, Owen Benton, Jeffrey G. Rau, J. Oitmaa, R.R.P. Singh, et al., Are Multiphase Competition and Order by Disorder the Keys to Understanding Yb 2 Ti 2 O 7?. Physical Review Letters, 2015, 115, pp.267208. 10.1103/PhysRevLett.115.267208. hal-01542075

HAL Id: hal-01542075

https://hal.science/hal-01542075v1

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 CC BY 4.0 - Attribution - International License
# Are multiphase competition & order-by-disorder the keys to understanding Yb$_{2}$Ti$_{2}$O$_{7}$?

L.D.C. Jaubert,¹ Owen Benton,¹ Jeffrey G. Rau,² J. Oitmaa,³

R.R.P. Singh,⁴ Nic Shannon,¹ and Michel J.P. Gingras²,⁵,⁶

¹Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0395, Japan

²Department of Physics and Astronomy, University of Waterloo,
200 University Avenue West, Waterloo, Ontario, N2L 3G1, Canada

³School of Physics, The University of New South Wales, Sydney 2052, Australia

⁴Department of Physics, University of California, Davis, California 95616, USA

⁵Perimeter Institute for Theoretical Physics, 31 Caroline North, Waterloo, Ontario, N2L 2Y5, Canada

⁶Canadian Institute for Advanced Research, Toronto, Ontario, M5G 1Z8, Canada

(Dated: May 21, 2015)

If magnetic frustration is most commonly known for undermining long-range order, as famously illustrated by spin liquids, the ability of matter to develop new collective mechanisms in order to fight frustration is no less fascinating, providing an avenue for the exploration and discovery of unconventional properties of matter. Here we study an ideal minimal model of such mechanisms which, incidentally, pertains to the perplexing quantum spin ice candidate Yb$_{2}$Ti$_{2}$O$_{7}$. Specifically, we explain how thermal and quantum fluctuations, optimized by order-by-disorder selection, conspire to expand the stability region of an accidentally degenerate continuous symmetry $U(1)$ manifold against the classical splayed ferromagnetic ground state that is displayed by the sister compound Yb$_{2}$Sn$_{2}$O$_{7}$. The resulting competition gives rise to multiple phase transitions, in striking similitude with recent experiments on Yb$_{2}$Ti$_{2}$O$_{7}$ [Lhotel *et al.*, Phys. Rev. B **89** 224419 (2014)]. Considering the effective Hamiltonian determined for Yb$_{2}$Ti$_{2}$O$_{7}$, we provide, by combining a gamut of numerical techniques, compelling evidence that such multiphase competition is the long-sought missing key to understanding the intrinsic properties of this material. As a corollary, our work offers a pertinent illustration of the influence of chemical pressure in rare-earth pyrochlores.

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

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

Starting from the most general model for anisotropic nearest-neighbor exchange on the pyrochlore lattice [19, 20], we consider a range of parameters known to describe Yb$_{2}$Ti$_{2}$O$_{7}$ [10]. We show that a splayed ferromagnetic (SFM) phase such as the one observed in

<figure><img src="image_2.png" /><figcaption>FIG. 1. Multiphase competition of the anisotropic nearest-neighbor pyrochlore model of Eq. (1), as determined by classical Monte Carlo simulations for {$J_{i=2,3,4}$} = {-0.22, -0.29, 0} meV. The $\psi_2$ and $\psi_3$ phases are selected by order-by-disorder within the antiferromagnetic $U(1)$ manifold and separated from the splayed ferromagnet (SFM) by a first-order transition, whose slope agrees with classical low-temperature expansion (green line, see Appendices). Rotating all spins of a given $\psi_2$ state by the same angle around their local easy-axes produces the entire $U(1)$ manifold, including the $\psi_3$ states [18]. The “splaying” angle of the spins in the SFM ground state is uniquely determined by the {$J_i$} parameters.</figcaption></figure>

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

### Model
We consider the anisotropic Hamiltonian

$$ \mathcal{H} = \sum_{\langle ij \rangle} \vec{S}_i \cdot \bar{\mathcal{J}}_{ij} \vec{S}_j \quad \text{with} \quad \bar{\mathcal{J}} = \begin{pmatrix} J_2 & J_4 & J_4 \\ -J_4 & J_1 & J_3 \\ -J_4 & J_3 & J_1 \end{pmatrix} \quad (1) $$

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

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

This zero-temperature framework sets the scene for the phase diagram of Fig. 1 computed by Monte Carlo (MC) simulations using parallel tempering [33, 34] and over-relaxation [35]. The U(1) degeneracy does not survive thermal fluctuations and collapses predominantly in the $\psi_2$ configurations, with a small $\psi_3$ island around the boundary due to the above-mentioned soft modes of excitations. An interesting consequence of this ObD competition is that the $\psi_2/\psi_3$ phases are precisely selected to optimize the entropy of the U(1) manifold for a given set of parameters and temperature. At finite temperature, this optimization puts the energetically selected SFM phase at a disadvantage and gives rise to

<table><caption>TABLE I. Critical value of the exchange parameter $J_1$ [meV], separating the splayed ferromagnet (SFM) from the U(1) manifold, as estimated from different methods (Monte Carlo, linear spin-wave, exact diagonalization, numerical linked-cluster and high-temperature expansion) at zero temperature and upon cooling from high temperature. Quantum and thermal fluctuations work together to stabilize the U(1) manifold over the SFM phase. In particular, quantum effects bring the SFM/U(1) boundary within error bars of Yb₂Ti₂O₇ parameters [10].</caption><thead><tr><th rowspan="2"></th><th colspan="1">Classical</th><th colspan="4">Quantum</th></tr><tr><th>MC</th><th>LSW</th><th>ED</th><th>NLC</th><th>HTE</th></tr></thead><tbody><tr><th>$T > 0$</th><td>-0.0340(5)</td><td>n/a</td><td>n/a</td><td>-0.070(3)</td><td>-0.06(3)</td></tr><tr><th>$T = 0$</th><td>-0.0289(1)</td><td>-0.062</td><td>-0.064(2)</td><td>n/a</td><td>n/a</td></tr></tbody></table>

multiple phase transitions for $J_1 \in [-0.034 : -0.029]$. Since the Hamiltonian of Eq. (1) supports a large variety of emergent degeneracies and potential ObD transitions at the boundaries between ordered phases [3, 36], we expect such a phenomenology to be a common feature of pyrochlores [29, 37] and frustrated magnetism in general [28]. But since temperature is not the only source of fluctuations, how do quantum fluctuations fit in this picture ?

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

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

*Quantum fluctuations at finite temperature* – Even if the joint analysis of thermal and quantum fluctuations is particularly challenging for such a frustrated problem, the build up of correlations when approaching the phase transition from high temperature remains accessible thanks to numerical linked-cluster computation [40–42] and high-temperature expansion [43] (see Appendix for technical details).

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

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

Based on these latter experiments [16], and thanks
to the present analysis, we believe it is finally possible
to flesh out a common framework for the powder and
single crystals which magnetically order. In Ref. [16],
both samples undergo a *double phase transition* with a
high-temperature non-ferromagnetic transition followed
by a low-temperature ferromagnetic first-order transi-
tion. Furthermore, the application of a magnetic field
*h* does not destroy the lower-temperature ferromagnetic
transition, but rather increases its temperature for *h* > 5

mT; the transition remaining first order up to $h_c \approx 20$ mT, before becoming continuous or vanishing, which is experimentally difficult to distinguish. These experimental results fit perfectly within the framework of our theory. The double transition is a direct consequence of the SFM/U(1) competition, as shown in Fig. 1 and for a similar range of temperature as in experiments. Furthermore, since the U(1) manifold is antiferromagnetic, it does not couple with *h*. The magnetic field thus only favors the SFM phase, and the first-order ferromagnetic transition is expected to persist until the U(1) phase is destroyed at $h_c$. Monte Carlo simulations confirm this scenario with $h_c \approx 15$ mT, followed by a crossover for $h > h_c$ (see Fig. 3.b). We also mention as a further similarity the experimental presence of a reversible bump in the single crystal magnetisation at $\approx 180$ mK, between the two transitions. According to our MC simulations, such feature could correspond to a $\psi_2/\psi_3$ ObD transition, but with the caveat that the $\psi_3$ phase does not persist above 50 mK in our classical phase diagram. In that case, structural disorder may play an important role [17, 45, 46], since it is known to i) favor $\psi_3$ over $\psi_2$ order [48] and ii) to be stronger in single crystals than powder samples where no bump was observed [16].

It should be noted that even if the parametrization
that we used, taken from Ref. [10], was done on samples
different from the ones in Ref. [16], the quantum shift-
ing of the SFM/U(1) boundary would bring this classi-
cal scenario within the experimental uncertainty of the
$J_1 = -0.09 \pm 0.03$ meV parametrization range. This
quantum shifting is best illustrated in Fig. 3, where the
structure factor calculated from classical simulations at
$J_1 = -0.0335$ meV is almost identical to quantum (NLC)
results at $J_1 = -0.06$ meV. As for neutron scattering

<figure><img src="image_3.png" /><figcaption>FIG. 2. The difference of correlators $\Delta C = C_{U(1)} - C_{SFM}$ computed with NLC confirms the quantum shifting of the boundary towards more negative values of $J_1$ than for the classical system, estimated at $J_1 = -0.070(3)$ meV. $J_1$ is given in the caption while $\{J_{i=2,3,4}\} = \{-0.22, -0.29, 0\}$. Inset: $\Delta C$ as computed from HTE for $J_1 = -0.09$ (red), $-0.06$ (cyan) and $-0.03$ (black). The breadth of each curve represents the uncertainty for NLC and HTE at low temperature.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. *Application to Yb₂Ti₂O₇*: a) structure factor as measured by neutron scattering and b) phase diagram in a field. a) The spin-flip (SF, top panels) and non-spin-flip (NSF, bottom panels) calculated by quantum NLC and classical Monte Carlo are almost identical when *J*₁ is shifted from −0.0335 meV to −0.06 meV (please note that the temperature has been renormalized by 5/4 for a better agreement). This agreement confirms the quantum shifting of the boundary at finite temperature. When approaching the phase transition, the MC structure factor for *J*₁ = −0.0335 meV reproduces the characteristic features of Yb₂Ti₂O₇ neutron scattering data measured at 300 mK (to be compared with Fig. 2a of Ref. [13]). On the other hand, the comparison between experiments [13] and classical simulations using the parametrization of Ref. [10] (*J*₁ = −0.09 meV) are noticeably less successful, especially around (220) [3]. This disagreement is *not* a criticism of the parametrization by Ross et al., but rather an emphasis on the importance of quantum fluctuations. The temperature of 0.65 K in the right panels has been chosen such that the ratio between measurement and transition temperatures is the same as in Ref. [13]. The color scale is fixed from 0 to the maximum SF intensity, except for the panel at 0.3 K where the color scale was chosen for visual comparison with experiments [13]. b) When a field is applied along the [001] direction, the antiferromagnetic ψ₂ phase gradually disappears in MC simulations. Our theory thus explains why the first-order transition only persists at low field in Yb₂Ti₂O₇ [16]. We used *J*₁ = −0.033 meV and the anisotropic g-tensor of Yb₂Ti₂O₇ [47]: *g*⊥ = 4.18 and *g*∥ = 1.77.</figcaption></figure>

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

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

*Conclusion* – Using a palette of complementary nu-
merical methods, our work sets a theoretical benchmark
common to both powder and single crystals of Yb₂Ti₂O₇.
We have shown how the multi-step ordering recently ob-
served in bulk measurements [16] naturally arises from
the competition between a SFM phase and a U(1) mani-
fold mediated by order-by-disorder selection. On a more
fundamental level, the underlying ordering mechanism

should be understood as the joint consequence of both
quantum and thermal fluctuations. The *T* = 0 frontier
between the U(1) manifold and the splayed ferromagnet
is shifted by quantum fluctuations, raising the interest-
ing possibility to modulate frustration by tuning quan-
tum fluctuations. In the general context of multiphase
competition, order-by-disorder can be viewed as a free-
energy optimization process which reinforces the stability
of the degenerate phase it acts upon, by naturally se-
lecting the subset of configurations with higher entropy
and/or quantum zero-point fluctuations.

In light of the numerous models and phases supported
by pyrochlores, ranging from spin liquids and spin ices
to (partially) ordered phases [11, 52–56], and subsequent
boundaries between them [3, 29, 37], our present work
is a paradigmatic example of why the properties of frus-
trated magnets should generically be understood as the
sum of competing phases, rather than coming from a sin-
gle controlling state. Some of these properties would in-
deed seem to be “coming from nowhere” in absence of a
global phase diagram. We expect such competition be-
tween neighboring phases to be particularly relevant to
some of the most difficult materials to characterize, such
as Tb₂Ti₂O₇ [45] and Er₂Sn₂O₇ [3, 57, 58], and to exac-
erbate the sample dependence issues [17, 45, 46] by the
proximity of phase boundaries. In that respect it is in-
teresting to note that Er$_{2}$Ti$_{2}$O$_{7}$, whose coupling param-
eters lie far away from any phase boundary [3, 24, 26], is
one of the most robust rare-earth pyrochlore for repro-
ducibility of experiments, while Yb$_{2}$Ti$_{2}$O$_{7}$ is essentially
the antithesis. Hence, we expect the interplay between
multiphase competition and disorder to become a very
topical question, necessary to account for experiments in
pyrochlores and frustrated magnetism, and promisingly
rich in exotic physics.

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

APPENDIX

Order parameters and correlators

As illustrated in Fig. 1 of the main text, three main
phases are studied in this paper: the U(1) manifold, di-

vided into $\psi_2$ and $\psi_3$ phases, and the splayed ferromag-
net (SFM). Their order parameters can be found in the
literature [3, 24, 25, 29, 37] and are repeated below for
convenience.

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

$$
\begin{align}
\vec{r}_0 &= \left(\frac{1}{2}, \frac{1}{2}, \frac{1}{2}\right), & \vec{r}_1 &= \left(\frac{1}{2}, -\frac{1}{2}, -\frac{1}{2}\right), \nonumber \\
\vec{r}_2 &= \left(-\frac{1}{2}, \frac{1}{2}, -\frac{1}{2}\right), & \vec{r}_3 &= \left(-\frac{1}{2}, -\frac{1}{2}, \frac{1}{2}\right). \tag{2}
\end{align}
$$

The spins at the four sites are denoted as $\vec{S}_{i=0,1,2,3} = (S_i^x, S_i^y, S_i^z)$ in the cubic coordinates, with $|\vec{S}_i| = 1/2$.

For the U(1) manifold, the order parameter, per tetra-
hedron, is a two-dimensional vector.

$$
\vec{m}_{\text{U(1)}} = \begin{pmatrix} m_{\text{U(1)}}^{\alpha} \\ m_{\text{U(1)}}^{\beta} \end{pmatrix} = \begin{pmatrix} \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) \\ \frac{1}{2\sqrt{2}}(-S_0^y + S_0^z + S_1^y - S_1^z - S_2^y - S_2^z + S_3^y + S_3^z) \end{pmatrix} \quad (3)
$$

To differentiate between the $\psi_2$ and $\psi_3$ states within the U(1) manifold, we define the angle $\theta_{\text{U(1)}} = \arctan(m_{\text{U(1)}}^{\beta}/m_{\text{U(1)}}^{\alpha})$; the function $\cos(6\theta_{\text{U(1)}})$ respectively equals to +1 and -1 for $\psi_2$ and $\psi_3$ states [3, 29, 37].

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

$$
\vec{m}_{\text{SFM1}} = \begin{pmatrix} m_{\text{SFM1}}^{\alpha} \\ m_{\text{SFM1}}^{\beta} \\ m_{\text{SFM1}}^{\gamma} \end{pmatrix} = \begin{pmatrix} \frac{1}{2}(S_0^x + S_1^x + S_2^x + S_3^x) \\ \frac{1}{2}(S_0^y + S_1^y + S_2^y + S_3^y) \\ \frac{1}{2}(S_0^z + S_1^z + S_2^z + S_3^z) \end{pmatrix} \quad (4)
$$

$$
\begin{align}
\vec{m}_{\text{SFM2}} &= \begin{pmatrix} m_{\text{SFM2}}^{\alpha} \\ m_{\text{SFM2}}^{\beta} \\ m_{\text{SFM2}}^{\gamma} \end{pmatrix} \nonumber \\
&= \begin{pmatrix} \frac{-1}{2\sqrt{2}}(S_0^y + S_0^z - S_1^y - S_1^z - S_2^y + S_2^z + S_3^y - S_3^z) \\ \frac{-1}{2\sqrt{2}}(S_0^x + S_0^z - S_1^x + S_1^z - S_2^x - S_2^z + S_3^x - S_3^z) \\ \frac{-1}{2\sqrt{2}}(S_0^x + S_0^y - S_1^x + S_1^y + S_2^x - S_2^y - S_3^x - S_3^y) \end{pmatrix} \tag{5}
\end{align}
$$

where $\vec{m}_{\text{SFM1}}$ is simply the global magnetization and
$\vec{m}_{\text{SFM2}}$ accounts for the splayed nature of the ferromag-
netism, i.e. the fact that the spins are not colinear. For
a given SFM ground state, both $\vec{m}_{SFM1}$ and $\vec{m}_{SFM2}$ are finite.

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

$$C_{IJ} = \langle m_I m_J \rangle - \langle m_I \rangle \langle m_J \rangle \quad (6)$$

where $m_I$ is a given order parameter. For the U(1) manifold, the use of the correlators is rather straightforward. However, for the SFM phase, since both $\vec{m}_{SFM1}$ and $\vec{m}_{SFM2}$ are finite, one needs to compute the matrix of correlators

$$\begin{pmatrix} C_{SFM1,SFM1} & C_{SFM1,SFM2} \\ C_{SFM2,SFM1} & C_{SFM2,SFM2} \end{pmatrix}. \quad (7)$$

By definition, this matrix is symmetric but a priori non-diagonal. Upon diagonalization, the maximum eigenvalue is kept as correlator of the SFM phase. Please note that by symmetry of the lattice, calculations can be simplified by considering only one component of the order parameters, say e.g. $\vec{m}_{SFM1}^{\alpha}$ and $\vec{m}_{SFM2}^{\alpha}$.

## Monte Carlo simulations

The Fig. 1 of the main text has been obtained by classical Monte Carlo simulations of Heisenberg spins on the pyrochlore lattice. The conventional cubic unit cell consists of 16 spins and the system size has $N = 3456$ spins. We used the standard Metropolis algorithm where a Monte Carlo step (MCs) is defined as $N$ random single-spin-flip attempts. To improve the quality of the simulations, parallel tempering [33, 34] and over-relaxation [35] were included in the simulations.

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

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

## Classical low-temperature expansion

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

## Linear spin wave theory

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

At first, the approach is similar to standard spin wave expansion. We rewrite the spin operators $\vec{S}_i$ in terms of Holstein-Primakoff bosons and perform a $1/S$ expansion around the local spin configuration of the chosen ordered state. At the harmonic level, our Hamiltonian takes the form $\mathcal{H} \approx \mathcal{H}_{LSW} = \mathcal{H}_{LSW}^{(0)} + \mathcal{H}_{LSW}^{(1)} + \mathcal{H}_{LSW}^{(2)}$. $\mathcal{H}_{LSW}^{(0)}$ is simply the classical energy of the ordered state around which we are expanding, while $\mathcal{H}_{LSW}^{(1)}$ and $\mathcal{H}_{LSW}^{(2)}$ contain only linear and quadratic terms in boson operators, respectively.

If the chosen ordered state is a classical ground state, then $\mathcal{H}_{LSW}^{(1)}$ must vanish. As for $\mathcal{H}_{LSW}^{(2)}$, it can be diagonalized by Fourier transform followed by a Bogoliubov transformation which will return real, positive, frequencies and therefore a meaningful excitation spectrum. The semi-classical energy can then be written

$$E_{\text{semi-cl}}^{0} = E_{\text{cl}}^{0} \left( 1 + \frac{1}{S} \right) + \frac{1}{2} \sum_{\vec{k}\lambda} \omega_{\vec{k}\lambda}, \quad (8)$$

where $E_{\text{cl}}^{0}$ is the classical ground state energy, $\omega_{\vec{k}\lambda}$ is the spin wave dispersion of band $\lambda$ at wave-vector $\vec{k}$.

However, if the chosen ordered state is not a classical ground state, then $\mathcal{H}_{LSW}^{(1)}$ may still vanish if the spin
configuration we are expanding around is a saddle point of the classical energy. As shown in Ref. [24], this is the case for the $\psi_2$ and $\psi_3$ configurations. The last remaining point is the diagonalization of $\mathcal{H}_{\text{LSW}}^{(2)}$. Following the approach developed by Coletta et al. [39], we add to $\mathcal{H}_{\text{LSW}}^{(2)}$ the positive definite term $V = \delta \sum_i a_i^\dagger a_i$, where the $a_i^\dagger$ and $a_i$ operators are Holstein-Primakoff bosons. This additional term $V$ does not change the classical energy and the parameter $\delta$ may be adjusted to the minimum value for which we can obtain a real, positive excitation spectrum. The energy calculated with the inclusion of $V$ is

$$E_{\text{semi-cl}}^{\delta} = E_{\text{cl}} \left(1 + \frac{1}{S}\right) - N \frac{\delta}{2} + \frac{1}{2} \sum_{\vec{k}\lambda} \omega_{\vec{k}\lambda}^{\delta} \quad (9)$$

where $E_{\text{cl}}$ is the classical energy of the saddle-point configuration and $\omega_{\vec{k}\lambda}^{\delta}$ is the spin wave dispersion calculated including the potential $V$. Since $V$ is a positive definite operator, $E_{\text{semi-cl}}^{\delta}$ is an upper bound on the semiclassical energy of the saddle-point configuration.

## Exact diagonalization

We consider the properties of the exchange Hamiltonian shown in Eq. (1) of the main text for a single tetrahedron ($N = 4$) and a single cubic unit cell of the pyrochlore lattice ($N = 16$) with periodic boundary conditions. Thanks to the small number of sites the full spectrum is analytically accessible for $N = 4$ and the

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

The series is organized into terms $P_n$ including up to $n$ tetrahedra. Explicitly carrying out the expansion using

lowest lying eigenstates can be easily found via standard Lanczos diagonalization for $N = 16$.

## Numerical linked cluster expansion

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

$$P(\mathcal{L})/N = \sum_{\mathcal{C} \subset \mathcal{L}} L(\mathcal{C})W(\mathcal{C}) \quad (10)$$

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

$$W(\mathcal{C}) = P(\mathcal{C}) - \sum_{\mathcal{C}' \subset \mathcal{C}} W(\mathcal{C}') \quad (11)$$

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

the embedding constants one has

$$P_0 = +P(C_0) \quad (12)$$

$$P_1 = -P(C_0) + \frac{1}{2}P(C_1) \quad (13)$$

$$P_2 = -\frac{3}{2}P(C_1) + P(C_2) \quad (14)$$

$$P_3 = +\frac{3}{2}P(C_1) - 5P(C_2) + 3P(C_3) \quad (15)$$

$$P_4 = -\frac{1}{2}P(C_1) + 10P(C_2) - 21P(C_3) \\ +3P(C_{4a}) + 6P(C_{4b}) + 2P(C_{4c}) \quad (16)$$

Note that $P(C_0)$ does not appear directly past first order. Following Ref. [41] the Euler resummation method is used on the final two terms to accelerate convergence. If we define the differences $S_n = P_{n+1} - P_n$ then the Euler
$C L(C)$

$C_0 \quad 1 \cdot$

$C_1 \quad \frac{1}{2}$

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

$C_3 \quad 3$

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

$C_{4b} \quad 6$

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

$C L(C)$

$C_2 \quad 1$

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

$C_{4a} \quad 3$

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

$C_{4c} \quad 2$

<figure><img src="image_10.png" /><figcaption>TABLE II. Clusters used for the NLC expansion. A graphical representation is shown along with the embedding constant $L(C)$.</figcaption></figure>

approximants are given by

$$E_2 = P_2 \qquad (17)$$

$$E_3 = E_2 + \frac{1}{2}S_3 \qquad (18)$$

$$E_4 = E_3 + \frac{1}{4}(S_3 + S_4) \qquad (19)$$

The difference between 3rd and 4th order of expansion is the uncertainty of our NLC computations.

## High Temperature Expansion

The High Temperature Expansions (HTE) are done up to order $\beta^8$, as explained in the book of Ref. [43]. The series expansion are then analyzed using Padé approximants, i.e. based on rational functions. We constructed all near-diagonal Padé approximants with 8 or 7 terms in the series, i.e. [4/4], [5/3], [3/5], [6/2], [2/6], [4/3], [3/4], [5/2], [2/5] where [m/n] stand for the powers of the polynomial in the numerator and the denominator. The resulting spread in Padé approximant values represents the error bars of our HTE computations.

[1] Y. Taguchi, Y. Oohara, H. Yoshizawa, N. Nagaosa, and Y. Tokura, Science **291**, 2573 (2001).

[2] M. Mostovoy, Phys. Rev. Lett. **96**, 067601 (2006).

[3] H. Yan, O. Benton, L. D. C. Jaubert, and N. Shannon, arXiv:1311.3501 (2013).

[4] B. Canals, M. Elhajal, and C. Lacroix, Physical Review B **78**, 214431 (2008).

[5] J. A. Hodges, P. Bonville, A. Forget, A. Yaouanc, P. Dalmas de Réotier, G. André, M. Rams, K. Królas, C. Ritter, P. C. M. Gubbens, C. T. Kaiser, P. J. C. King, and C. Baines, Phys. Rev. Lett. **88**, 077204 (2002).

[6] Y. Yasui, M. Soda, S. Iikubo, M. Ito, M. Sato, N. Hamaguchi, T. Matsushita, N. Wada, T. Takeuchi, N. Aso, and K. Kakurai, Journal of the Physical Society of Japan **72**, 3014 (2003).

[7] 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, Phys. Rev. Lett. **103**, 227202 (2009).

[8] A. Yaouanc, P. Dalmas de Réotier, C. Marin, and V. Glazkov, Phys. Rev. B **84**, 172408 (2011).

[9] R. M. D'Ortenzio, H. A. Dabkowska, S. R. Dunsiger, B. D. Gaulin, M. J. P. Gingras, T. Goko, J. B. Kycia, L. Liu, T. Medina, T. J. Munsie, D. Pomaranski, K. A. Ross, Y. J. Uemura, T. J. Williams, and G. M. Luke, Phys. Rev. B **88**, 134428 (2013).

[10] K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Physical Review X **1**, 021002 (2011).

[11] L. Savary and L. Balents, Physical Review Letters **108**, 037202 (2012).

[12] L. Savary and L. Balents, Phys. Rev. B **87**, 205130 (2013).

[13] L.-J. Chang, S. Onoda, Y. Su, Y.-J. Kao, K.-D. Tsuei, Y. Yasui, K. Kakurai, and M. R. Lees, Nature Communications **3**, 992 (2012).

[14] L. 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, Nature Communications **5** (2014), 10.1038/ncomms5970.

[15] L. Pan, N. J. Laurita, K. A. Ross, E. Kermarrec, B. D. Gaulin, and N. P. Armitage, arXiv (2015), 1501.05638.
[16] E. Lhotel, S. R. Giblin, M. R. Lees, G. Balakrishnan, L. J. Chang, and Y. Yasui, Phys. Rev. B **89**, 224419 (2014).

[17] K. A. Ross, T. Proffen, H. A. Dabkowska, J. A. Quil- liam, L. R. Yaraskavitch, J. B. Kycia, and B. D. Gaulin, Phys. Rev. B **86**, 174424 (2012).

[18] A. Poole, A. S. Wills, and E. Lelievre-Berna, Journal of Physics: Condensed Matter **19**, 452201 (2007).

[19] S. H. Curnoe, Physical Review B **75**, 212404 (2007).

[20] S. H. Curnoe, Physical Review B **78**, 094418 (2008).

[21] A. Yaouanc, P. Dalmas de Réotier, P. Bonville, J. A. Hodges, V. Glazkov, L. Keller, V. Sikolenko, M. Bartkowiak, A. Amato, C. Baines, P. J. C. King, P. C. M. Gubbens, and A. Forget, Phys. Rev. Lett. **110**, 127207 (2013).

[22] Z. L. Dun, E. S. Choi, H. D. Zhou, A. M. Hallas, H. J. Silverstein, Y. Qiu, J. R. D. Copley, J. S. Gardner, and C. R. Wiebe, Phys. Rev. B **87**, 134408 (2013).

[23] J. D. M. Champion and P. C. W. Holdsworth, Journal of Physics: Condensed Matter **16**, S665 (2004).

[24] L. Savary, K. A. Ross, B. D. Gaulin, J. P. C. Ruff, and L. Balents, Phys. Rev. Lett. **109**, 167201 (2012).

[25] M. E. Zhitomirsky, M. V. Gvozdikova, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. **109**, 077204 (2012).

[26] A. W. C. Wong, Z. Hao, and M. J. P. Gingras, Phys. Rev. B **88**, 144402 (2013).

[27] J. Oitmaa, R. R. P. Singh, B. Javanparast, A. G. R. Day, B. V. Bagheri, and M. J. P. Gingras, Phys. Rev. B **88**, 220404 (2013).

[28] C. Pinettes, B. Canals, and C. Lacroix, Phys. Rev. B **66**, 024422 (2002).

[29] P. A. McClarty, P. Stasiak, and M. J. P. Gingras, Phys. Rev. B **89**, 024425 (2014).

[30] The number between parentheses represents the error bar on the last digit, i.e. $J_1 = -0.09 \pm 0.03$ here.

[31] S. B. Lee, S. Onoda, and L. Balents, Phys. Rev. B **86**, 104412 (2012).

[32] Z. Hao, A. G. R. Day, and M. J. P. Gingras, Phys. Rev. B **90**, 214430 (2014).

[33] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. **57**, 2607 (1986).

[34] C. J. Geyer, Computing Science and Statistics: Proceed- ings of the 23rd Symposium on the Interface, p. 156 (1991).

[35] M. Creutz, Physical Review D **36**, 515 (1987).

[36] O. Benton, *Classical and Quantum Spin Liquids on the Pyrochlore lattice*, Ph.D. thesis, University of Bristol (2014).

[37] G.-W. Chern, ArXiv e-prints (2010), arXiv:1008.3038 [cond-mat.str-el].

[38] P. Fazekas, *Lecture Notes on Electron Correlation and Magnetism* (World Scientific, Singapore, 1999).

[39] T. Coletta, M. E. Zhitomirsky, and F. Mila, Phys. Rev. B **87**, 060407 (2013).

[40] M. Rigol, T. Bryant, and R. R. P. Singh, Physical Review E **75**, 061118 (2007).

[41] R. Applegate, N. R. Hayre, R. R. P. Singh, T. Lin, A. G. R. Day, and M. J. P. Gingras, Phys. Rev. Lett. **109**, 097205 (2012).

[42] B. Tang, E. Khatami, and M. Rigol, Computer Physics Communications **184**, 557 (2013).

[43] J. Oitmaa, C. Hamer, and W. Zheng, *Series Expansion Methods for Strongly Interacting Lattice Models* (Cambridge University Press, Cambridge, UK, 2006).

[44] N. R. Hayre, K. A. Ross, R. Applegate, T. Lin, R. R. P. Singh, B. D. Gaulin, and M. J. P. Gingras, Phys. Rev. B **87**, 184423 (2013).

[45] T. Taniguchi, H. Kadowaki, H. Takatsu, B. Fåk, J. Ollivier, T. Yamazaki, T. J. Sato, H. Yoshizawa, Y. Shimura, T. Sakakibara, T. Hong, K. Goto, L. R. Yaraskavitch, and J. B. Kycia, Phys. Rev. B **87**, 060408 (2013).

[46] G. Sala, M. J. Gutmann, D. Prabhakaran, D. Pomaran- ski, C. Mitchelitis, J. B. Kycia, D. G. Porter, C. Castel- novo, and J. P. Goff, Nature Materials **13**, 488 (2014).

[47] J. A. Hodges, P. Bonville, A. Forget, and G. Andre, Canadian Journal of Physics **79**, 1373 (2001).

[48] V. S. Maryasin and M. E. Zhitomirsky, Phys. Rev. B **90**, 094412 (2014).

[49] Z. L. Dun, M. Lee, E. S. Choi, A. M. Hallas, C. R. Wiebe, J. S. Gardner, E. Arrighi, R. S. Freitas, A. M. Arevalo- Lopez, J. P. Attfield, H. D. Zhou, and J. G. Cheng, Phys. Rev. B **89**, 064401 (2014).

[50] H. D. Zhou, J. G. Cheng, A. M. Hallas, C. R. Wiebe, G. Li, L. Balicas, J. S. Zhou, J. B. Goodenough, J. S. Gardner, and E. S. Choi, Phys. Rev. Lett. **108**, 207206 (2012).

[51] I. Mirebeau, I. N. Goncharenko, P. Cadavez-Pares, S. T. Bramwell, M. J. P. Gingras, and J. S. Gardner, Nature **420**, 54 (2002).

[52] J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Reviews of Modern Physics **82**, 53 (2010).

[53] M. J. P. Gingras and P. A. McClarty, Reports on Progress in Physics **77**, 056501 (2014).

[54] M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, Phys. Rev. X **4**, 011007 (2014).

[55] B. Javanparast, Z. Hao, M. Enjalran, and M. J. P. Gin- gras, Phys. Rev. Lett. **114**, 130601 (2015).

[56] J. G. Rau and M. J. P. Gingras, arXiv (2015), 1503.04808.

[57] P. M. Sarte, H. J. Silverstein, B. T. K. V. Wyk, J. S. Gardner, Y. Qiu, H. D. Zhou, and C. R. Wiebe, Journal of Physics: Condensed Matter **23**, 382201 (2011).

[58] S. Guitteny, S. Petit, E. Lhotel, J. Robert, P. Bonville, A. Forget, and I. Mirebeau, Phys. Rev. B **88**, 134408 (2013).

[59] N. Shannon, K. Penc, and Y. Motome, Physical Review B **81**, 184409 (2010).