# The multiple symmetry sustaining phase transitions of spin ice

V. Raban  Université de Lyon, ENS de Lyon, Université Claude Bernard, CNRS,
Laboratoire de Physique, F-69342 Lyon, France  Laboratoire Charles Coulomb
(L2C), Université de Montpellier, CNRS, Montpellier, France  C. T. Suen
Université de Lyon, ENS de Lyon, Université Claude Bernard, CNRS, Laboratoire
de Physique, F-69342 Lyon, France  L. Berthier  Laboratoire Charles Coulomb
(L2C), Université de Montpellier, CNRS, Montpellier, France  P. C. W.
Holdsworth  Université de Lyon, ENS de Lyon, Université Claude Bernard, CNRS,
Laboratoire de Physique, F-69342 Lyon, France

###### Abstract

We present the full phase diagram of the dumbbell model of spin ice as a
function of temperature, chemical potential and staggered chemical potential
which breaks the translational lattice symmetry in favour of charge crystal
ordering. We observe a double winged structure with five possible phases,
monopole fluid (spin ice), fragmented single monopole crystal phases and
double monopole crystal, the zinc blend structure. Our model provides a
skeleton for liquid-liquid phase transitions and for the winged structures
observed for itinerant magnets under pressure and external field. We relate
our results to recent experiments on Ho2Ir2O7 and propose a wide ranging set
of new experiments that exploit the phase diagram, including high pressure
protocols, dynamical scaling of Kibble-Zurek form and universal violations of
the fluctuation-dissipation theorem.

##  I Introduction

Over the last decade, spin ice models and materials Harris  _et al._ (1997);
Bramwell and Gingras (2001) have emerged as model systems for the study of
generalized electrostatics on a lattice Isakov  _et al._ (2004a); Castelnovo
_et al._ (2008); Ryzhkin (2005); Jaubert and Holdsworth (2009); Castelnovo
_et al._ (2011); Brooks-Bartlett  _et al._ (2014); Kaiser  _et al._ (2015).
The emergence of the electrostatics can best be seen by replacing the point
dipole moments of spin ice by infinitesimally thin magnetic needles, lying
along the axes linking the centres of adjoining tetrahedra Möller and Moessner
(2006)(see Fig. 1). Within this dumbbell approximation Castelnovo  _et al._
(2008), the pyrochlore lattice of magnetic moments transforms den Hertog and
Gingras (2000); Isakov  _et al._ (2005) into a diamond lattice of vertices for
magnetic charge. The needles carry magnetic flux and dumbbells of effective
magnetic charge which touch at the vertices. By construction the ensemble of
low energy “Pauling states” Pauling (1935) with two spins into and two out of
each tetrahedron are degenerate in this approximation, with charge neutrality
imposed at each vertex. These ground states form a vacuum from which magnetic
monopole quasi-particles are excited by reversing the orientation of a needle,
breaking the ice rules on a pair of neighbouring sites Castelnovo  _et al._
(2008). Double monopoles can also be created by reversing a second needle, for
a vertex with all needles in or all out. The emerging Coulomb fluid of
magnetic origin is often referred to as a magnetolyte Jaubert and Udagawa
(2019) in analogy with its electrical counterpart.

Figure 1: From spins to dumbbells: The point dipoles are extended to needles
touching at the diamond lattice centres. The needles carry magnetic flux and
charge $q=\pm m/a$ at each end. In a 2in-2out configuration (left) the vertex
is charge neutral. A 3in-1out (3out-1in) configuration carries a monopole
charge $Q=2m/a$ ($-Q=-2m/a$) (left). A 4in (4out) configuration carries a
double monopole charge $2Q=4m/a$ ($-2Q=-4m/a$).

In this paper we study the full phase diagram of the dumbbell model, including
a staggered chemical potential, $\Delta$, which breaks a $Z_{2}$ translational
symmetry of the diamond lattice, favouring monopole and double monopole
crystallisation into bi-partite ionic cristals. The staggered chemical
potential lifts the degeneracy between single and double monopoles at the
crystallisation transition in a manner compatible with the staggered internal
magnetic field offered by iridium ions in the spin ice material Ho2Ir2O7
Lefrançois  _et al._ (2017).

As shown in Fig. 2, the dumbbell model offers a rich phase diagram in the
three dimensional space of parameters $\Delta$, energy scale $\nu$ fixing the
monopole and double monopole chemical potentials: $\mu=-\nu$, $\mu_{2}=-4\nu$,
and temperature $T$. The central plane with $\Delta=0$ corresponds to the
standard spin ice phase diagram within this approximation Melko and Gingras
(2004); Guruciaga  _et al._ (2014), with a transition from spin ice to “all-
in-all-out” (AIAO) order that changes from first to second order in a multi-
critical region. In the monopole language AIAO order corresponds to an ionic
crystal of double monopoles with the zinc blend structure. Emerging from this
region, there is a double winged structure of phase boundaries that terminate
in continuous lines of critical end points. The five phases separated by the
boundaries are the Coulomb fluid (spin ice) phase, a fragmented monopole
crystal Brooks-Bartlett  _et al._ (2014); Borzi  _et al._ (2013) in which the
magnetic moments appear to break up into independent divergence full and
divergence free parts and the double monopole crystal AIAO phase.

As $\Delta$ breaks the translational symmetry all transitions, away from the
central plane, are symmetry sustaining. In this sense the transition from
monopole fluid to single monopole crystal is thermodynamically equivalent to
the liquid-gas transition and that from single to double monopole crystal is
equivalent to liquid-liquid transitions observed experimentally in supercooled
liquids Katayama  _et al._ (2000); Sastry and Austen Angell (2003); Brovchenko
_et al._ (2005). Entirely analogous sets of phase transitions also occur in
itinerant magnetic compounds under pressure and in the presence of an external
field Kaluarachchi  _et al._ (2017); Kotegawa  _et al._ (2011). A consequence
of our work is that we are able to offer a generic framework and minimal model
to generate such seemingly exotic behaviour, occurring in diverse domains of
physics and chemistry.

Figure 2: Dumbbell model phase diagram: (0) monopole fluid ($\phi_{0}$), (1)
monopole crystal ($\phi_{1}$), (2) double monopole crystal ($\phi_{2}$).
Surfaces show $1^{st}$ order and solid lines $2^{nd}$ order transitions and
dotted lines show the multi-critical region. The long dashed lines show the
extension of the $2^{nd}$ order lines to infinity. Chemical potentials for
single and double monopoles are $\mu=-\nu$, $\mu_{2}=-4\nu$ (see eqn. (3)).
For the phase diagram of the $S=2$ Blume Capel model $\bar{\Delta}$ replaces
$\Delta$ and $\bar{\nu}$ replaces $\nu$.

Inspired by the $S=2$ Blume-Capel model Plascak  _et al._ (1993), in the next
section we will provide concrete and quantitative evidence for the existence
of the double winged phase diagram shown in Fig. 2, introducing general Blume-
Capel models, providing a detailed explanation of the multi-critical region
and investigating one of the continuous set of critical end points that takes
the model from the spin ice monopole fluid to fragmented monopole crystal. In
section III we present dynamical finite size scaling results in the region of
the critical point and show that it exihibits dynamical Kibble-Zurek scaling
in the three dimensional Ising universality class. We also present results
showing the universal violation of the fluctuation-dissipation relation
consistent with this universality class. In section IV we relate our results
to the observed monopole driven phase transition for spin ice materials in a
magnetic field $\vec{H}$ in the $[111]$ direction showing that monopole
crystallisation thermodynamics leads to a quantitative prediction of the
$\vec{H},T$ phase diagram. In section V we give some discussion, putting our
results in the wider context of liquid-liquid phase transitions and the
temperature-field-pressure phase diagram of itinerant magnets. We conclude
this section, returning to frustrated magnets, in particular Ho2Ir2O7 and the
possibility of observing such a rich phase diagram and its consequences in
future experiments.

The Kelvin energy scale is used throughout, fixing Boltzmann’s constant to
unity. We also set the permeability of free space $\mu_{0}=1$ so that the
field $\vec{H}$ is measured in Tesla. We follow standard notation for spin ice
simulations and refer to a dimensionless length $L$, measured in cubic units.
Each cubic cell contains 16 spins (dumbbells) so that the number of tetrahedra
(monopole sites), $N_{0}=8L^{3}$. In this paper quantitative measures refer to
the spin ice material Dy2Ti2O7 (DTO) for which diamond lattice constant
$a=4.33$ Å, the nearest neighbour spin distance
$r_{nn}=\frac{\sqrt{3}a}{2}=3.74$ Å, and cube length
$a_{c}=\frac{4a}{\sqrt{3}}\approx 10$ Å(see Fig. (1)).

##  II Monopole Crystal Phase Diagram

The dumbbell model is an excellent approximation to the dipolar spin ice model
(DSI) which is characterised by short range exchange interactions and dipole
interactions which provide long range forces for the monopole quasi-particles
den Hertog and Gingras (2000); Henelius  _et al._ (2016). The dumbbell model
captures all features of the DSI except for a low temperature ordering
transition which indicates the lifting of the degeneracy of the Pauling
states. Above this energy scale, the DSI shows a phase transition on varying
the ratio of the exchange terms to dipolar interaction, taking the model from
the spin ice phase to the AIAO phase Melko  _et al._ (2001). The transition
appears to change from first to second order via a multi-critical point
Guruciaga  _et al._ (2014).

###  II.1 Blume-Capel models

Figure 3: Blume Capel S=1: (0) Paramagnet, (1) antiferromagnet. Surfaces show
$1^{st}$ and lines $2^{nd}$ order transitions and the point shows a tri-
critical point.

Such physics is generically provided Cardy (1996) by the Blume-Capel (BC1)
model Blume (1966); Capel (1966), developed by Blume, Emery and Griffiths
Blume  _et al._ (1971) to study mixtures of 3He and 4He. In this model Ising-
like degrees of freedom, which could be spins or occupation numbers for a
neutral two component lattice fluid, take on values, $S_{i}=0,\pm 1$. Contact
with spin ice corresponds to the antiferromagnetic case with spins on a
bipartite lattice such as square, cubic or diamond with energy function

${\cal{H}}_{BC}=J\sum_{ij}S_{i}S_{j}+\bar{\nu}\sum_{i}S_{i}^{2}-\bar{\Delta}\sum_{i=1,N_{0}}(-1)^{i}S_{i}\;,$
(1)

where $J<0$ is a coupling constant, $\bar{\nu}$ is the energy scale for
exciting a site $i$, $\bar{\Delta}$ is a staggered field that breaks the
$Z_{2}$ symmetry of the bipartite lattice. Although the function
${\cal{H}}_{BC}$ may generically be referred to as the Hamiltonian, for future
reference we take the Hamiltonian to be the many body term only. The
parameters $\bar{\nu}$ and $\bar{\Delta}$ can be interpreted as Lagrange
multipliers which allow the evolution from the canonical to less constrained
ensembles, so that the single site terms contribute to the free energy but not
the internal energy Landau and Lifshitz (1959). A suitable order parameter can
be defined

$\phi=\bigg{\langle}\frac{1}{N_{0}}\sum_{i=1,N_{0}}\phi_{i}\bigg{\rangle},$
(2)

where $\langle\dots\rangle$ is a thermal average. The term
$\phi_{i}=\epsilon_{i}S_{i}$ distinguishes the two sublattices with
$\epsilon_{i}=1$ on an $A$ site and $-1$ on a $B$ site.

For $\bar{\Delta}=0$, on increasing $\bar{\nu}$, the transition changes from
$2^{nd}$ order, in the Ising universality class, to $1^{st}$ order via a tri-
critical point. The staggered term $\bar{\Delta}$ is conjugate to $\phi$ and
therefore guarantees a winged structure, as shown in Fig. (3). The first order
transitions terminate along a line of critical end points for finite
$\bar{\Delta}$ and temperature. The winged phase boundaries and finite
temperature critical end points stretch out to
$\bar{\Delta}\rightarrow\infty$, as even when the site occupation is perfectly
partitioned with $S_{i}=1$ on $A$ sites only and $S_{i}=-1$ on $B$ sites only,
the interaction between the sublattices remains, allowing for a singular jump
in site occupation at finite temperature. As $\bar{\Delta}$ breaks the lattice
symmetry, the transitions at the critical end points are symmetry sustaining.
They are characterised by an emergent Ising like order parameter at each point
and in this sense are liquid-gas like.

The Blume-Capel model can be extended Plascak  _et al._ (1993) to higher
values of $S$. Of particular interest is $S=2$ (BC2) which greatly resembles
the dumbbell model of spin ice. The order parameter is now defined on the
interval $-2\leq\phi\leq 2$ and according to mean field Plascak  _et al._
(1993) and pair approximation calculations Lara and Plascak (1998) the BC2
model allows for two ordered phases corresponding to $|\phi|\sim 1$ (referred
to as $\phi_{1}$) and $|\phi|\sim 2$ ($\phi_{2}$) as well as the disordered
phase with $\phi\approx 0$ ($\phi_{0}$). As a consequence, adding a finite
staggered field, $\bar{\Delta}$ to the BC2 energy function will open out a
double winged structure as shown qualitatively in Fig. (2) for the dumbbell
model and discussed in detail below.

###  II.2 The dumbbell model

Returning to the dumbbell model, the charge on vertex $i$ of the diamond
lattice takes values $Q_{i}=0,\pm Q,\pm 2Q$ with $Q=2m/a$, $m$ the magnetic
moment associated with a spin and $a$ the lattice constant (see Fig. (1)),
from which one can define a site occupation variable
$\hat{n}_{i}=\frac{Q_{i}}{Q}=0,\pm 1,\pm 2$ in analogy with the BC2 model
variables $S_{i}$. A magnetic north (south) monopole carries charge
$+Q\,(-Q)$. Within the dumbbell approximation, the dipolar spin ice
Hamiltonian for excitations above the lowest energy 2in-2out states can be
written:

${\cal{H}}={u(a)\over{2}}\sum_{i\neq
j}\left({a\over{r_{ij}}}\right)\hat{n}_{i}\hat{n}_{j}+\nu\sum_{i}\hat{n}_{i}^{2}-\Delta\sum_{i=1,N_{0}}(-1)^{i}\hat{n}_{i},$
(3)

where $u(a)={\mu_{0}Q^{2}\over{4\pi a}}$ is the nearest neighbour Coulomb
energy scale for a pair of monopoles. The mapping thus re-formulates the spin
ice problem as a lattice Coulomb fluid in the grand ensemble Castelnovo  _et
al._ (2008); Jaubert and Holdsworth (2009); Jaubert  _et al._ (2011, 2013);
Brooks-Bartlett  _et al._ (2014); Kaiser  _et al._ (2013, 2018) with chemical
potential for monopole and double monopole creation $\mu=-\nu$ and
$\mu_{2}=-4\nu$ respectively. The chemical potential $\mu$ can be calculated
for each material from the parameters of the corresponding (DSI) and that for
double monopoles is constrained to $\mu_{2}=4\mu$ by the spin Hamiltonian.
Here we add a staggered chemical potential term $\Delta$ which lifts the
degeneracy for quasi-particles with charge $\pm Q$ (and with charge $\pm 2Q$)
on the sublattices $A$ and $B$, $\mu^{A}_{\pm}=\mu\pm\Delta$,
$\mu^{B}_{\pm}=\mu\mp\Delta$ and the convention is such that $\Delta>0$
reduces the energy scale for creation of monopoles (double monopoles) with
positive charge on $A$ sites and with negative charge on $B$ sites.

The Hamiltonian in eqn. (3) is a BC2 type energy function with long range
Coulomb interactions, with order parameter $\phi$ given by eqn. (2) and with
$\hat{n}_{i}$ replacing $S_{i}$. However, the BC2 and dumbbell models are
different as they have different configurational phase spaces and so have
different entropies. In the dumbbell model one must take into account the
fragmented spin background Brooks-Bartlett  _et al._ (2014), the so-called
Dirac strings Castelnovo  _et al._ (2008); Castelnovo and Holdsworth (2019),
which emerge in the electrostatics as a divergence free electric field giving
Coulomb phase correlations Isakov  _et al._ (2004a); Henley (2010) at low
temperature in the $\phi_{0}$ phase. These strings possess their own
configurational entropy independently of the charges. As a consequence, for
zero or finite monopole density and even in the monopole crystal phases, the
entropy remains different from that of a lattice Coulomb fluid and hence of
the BC2 model. The zero temperature limits for these entropies are well known.
The entropy number density of the Coulomb fluid phase is the Pauling entropy,
$s_{0}\approx\ln(3/2)=0.405$ per tetrahedron. The entropy of the fragmented
monopole crystal is that of an ensemble of hard core dimers on a diamond
lattice Nagle (1966); Brooks-Bartlett  _et al._ (2014),
$s_{1}\approx\ln(1.3)=0.262$ while that of the double monopole crystal is
zero. One can develop an expression for the entropy of both monopoles and
strings at the Pauling level of approximation Pauling (1935); Ryzhkin (2005)
which works well in the monopole fluid phase Kaiser  _et al._ (2018);
Castelnovo and Holdsworth (2019), but breaks down in the crystal phases. More
detailed analysis requires a return to the field theoretic description of the
charges and its ensuring lattice Helmholtz decomposition Brooks-Bartlett  _et
al._ (2014); Maggs and Rossetto (2002).

###  II.3 The double winged phase diagram

Figure 4: $\Delta-T$ plane: (Left) shaded plane through the full phase diagram
at fixed $\nu$ \- dotted black lines show the intercepts of the phase
boundaries with the plane. (Right) The fixed $\nu$ plane rotated to give
$\Delta$ vs. $T$. Green lines show the phase boundaries, the dotted black line
shows an isothermal trajectory in the plane.

The entropy terms make some quantitative difference but similar phase diagrams
can be expected for the two models as can be seen from thermodynamic
arguments. The monopole free energy can be written

$\Omega=N_{0}\big{(}u_{c}+\nu n+4\nu n_{2}-\Delta\phi-sT\big{)},$ (4)

where $u_{c}$ and $s$ are the Coulomb energy and entropy number densities. As
we are dealing with ionic crystals, the energy of the three phases are known
exactly at zero temperature Brooks-Bartlett  _et al._ (2014):
$u_{c}^{\phi_{0}}=0$, $u_{c}^{\phi_{1}}=-u(a)\alpha/2$,
$u_{c}^{\phi_{2}}=-2u(a)\alpha$, where $\alpha=1.638$ is the Madelung constant
for a diamond lattice. Hence there are zero temperature phase boundaries
between the three phases with $\Delta\geq 0$

$\displaystyle\phi_{0}$
$\displaystyle:\phi_{1},\;\;\Delta=-\frac{u(a)\alpha}{2}+\nu\;K$ (5)
$\displaystyle\phi_{1}$
$\displaystyle:\phi_{2},\;\;\Delta=-\frac{3u(a)\alpha}{2}+3\nu\;K.$

Notice that, as both the Coulomb energies and the chemical potentials scale
with the square of the charge ($\mu_{2}=-4\nu$) the five phases intercept the
$\Delta=0$ axis at the same point, $\nu^{\ast}=u\alpha/2$. For smaller $\nu$
the Coulomb energy of the double monopole crystal wins out corresponding to
spin ice models passing directly into the AIAO phase. However, as $\Delta$
couples linearly to the charge the wings spread out from this point in the
$T=0$ plane.

The finite temperature phase boundaries can be estimated from the Clapeyron
equation for equilibrium between phases $i$ and $j$:

$\frac{d\Delta}{dT}=-\frac{s_{i}-s_{j}}{\phi_{i}-\phi_{j}},$ (6)

where $\phi_{i}$ and $s_{i}$ are the order parameter and entropy densities of
phase $i$. At small temperature we can assume that both order parameter and
entropy are constant: $s_{0}=0.405$, $s_{1}=0.262$ and $s_{2}=0$, $\phi=2,1$
and $0$, giving intercepts and slopes for the phase boundaries in a $T-\Delta$
plane for fixed $\nu>\nu^{\ast}$. At higher temperatures this “fixed entropy
approximation” will break down and the lines should terminate in critical end
points as illustrated in Fig. (4).

Figure 5: Multiple monopole crystallization: (Upper) The order parameter
$\phi$ vs. $\Delta$ simulated from the dumbbell model at fixed $\nu=4.35$ K.
Simulations for $N_{0}=4096$ ($L=8$) and periodic boundaries. All values are
in Kelvin. (Lower) Probability density $P(\phi)$ for $\Delta=2.05$ K, $T=0.3$K
and $\nu=4.35$ K.

We have tested this proposition numerically. In Fig. (5) we show $\phi$ vs
$\Delta$ for simulations with $L=8$ for different temperatures for fixed
$\nu=4.35$ K and $u(a)=2.88$ K, values estimated for DTO 111The value of
$u(a)=2.88$ used here corresponds to a magnetic moment for the spins,
$m=9.87\mu_{B}$, as deduced from crystal field calculations Yavors’kii  _et
al._ (2008) rather than the $10\mu_{B}$ often quoted in the literature, which
gives $u(a)=3.07$K and a corresponding difference in the energy scale for the
phase diagram . For the lowest temperature, sharp steps are indeed observed in
$\phi$ from $\phi\approx 0$ to $\phi\approx 1$ and from $\phi\approx 1$ to
$\phi\approx 2$ at a value slightly greater than $2$ K and $6$ K respectively.
The data is consistent with two $1^{st}$ order phase transitions from
$\phi_{0}$ to $\phi_{1}$ and from $\phi_{1}$ to $\phi_{2}$. As the temperature
is increased the steps in $\phi$ become rounded, consistent with the model
passing through a critical end point with the transitions evolving to
crossovers at high temperature. The singular nature of the transition between
$\phi_{0}$ and $\phi_{1}$ at $T=0.3$ is confirmed in the lower panel where we
show the probability density $P(\phi)$ estimated during the simulation. The
distribution is sharply peaked near $\phi=1$ but shows a lower peak in
probability near $\phi=0$, consistent with fluctuations between metastable
states separated by a finite jump in order parameter space. The inequality in
the peak heights shows that for these parameters, the system has passed into
the ordered phase. The lower peak in distribution occurs at a small but finite
value of $\phi$, consistently with $\Delta$ breaking the $Z_{2}$ symmetry of
the lattice even in the $\phi_{0}$ phase. The five phases confirming the
double winged structure are indeed the Coulomb fluid (spin ice) phase
($\phi_{0}$), the two fragmented monopole crystal phases Brooks-Bartlett  _et
al._ (2014); Borzi  _et al._ (2013) ($\phi_{1}$) and the double monopole
crystal AIAO phases ($\phi_{2}$).

The position of the $1^{st}$ order transitions in parameter space can be
estimated using eqns. (5) and (6). Taking the DTO values for $\nu$, and $u(a)$
and $\Delta>0$ the zero temperature intercept of the two phase boundaries are

$\displaystyle\phi_{0}$ $\displaystyle:\phi_{1},\;\;\Delta=1.99\;K$ (7)
$\displaystyle\phi_{1}$ $\displaystyle:\phi_{2},\;\;\Delta=5.97\;K.$

Assuming complete jumps in the order parameter at the transition, one finds
for $T=0.3$ K

$\displaystyle\phi_{0}$ $\displaystyle:\phi_{1},\;\;\Delta=2.03\;K$ (8)
$\displaystyle\phi_{1}$ $\displaystyle:\phi_{2},\;\;\Delta=6.05\;K,$

in close agreement with the results of Fig. (5).

###  II.4 A critical end point

We have made a quantitative estimate of the position of one critical end
point, that for the transition from $\phi_{1}$ to $\phi_{0}$ for $\nu=4.35$ K.
This can be extracted from the crossings of the Binder cumulant Binder (1981),
$B_{4}$ for the emergent Ising like order parameter at the critical end point,
$\varphi=\phi-\phi_{c}$, where $\phi_{c}=\langle\phi(T_{c})\rangle$:

$B_{4}(T,\Delta)=\frac{\langle(\varphi)^{4}\rangle}{\langle(\varphi)^{2}\rangle^{2}}.$
(9)

The parameters $\Delta_{c}$, $T_{c}$ and $\phi_{c}$ were estimated using an
iterative procedure. A first estimate of $T_{c}$ and $\Delta_{c}$ was made by
following the evolution of $P(\phi)$ from a double to single peak
distribution. From here a more accurate estimate of $\Delta$ was found from
the maximum of the susceptibility for $\phi$. This estimate was found to be
invariant under small temperature changes and the result can be established
with high precision Hamp  _et al._ (2015). We find $\Delta_{c}=2.03745\pm
0.00005$ K. The evolution of $B_{4}$ with temperature for this $\Delta_{c}$ is
shown in Fig. (6) for system sizes $L=8,10,12$ and for $\phi_{c}=0.42$. A
crossing point is found for $T=0.36752\pm 0.00001$ K with
$B_{4}(T,\Delta)=2.03\pm 0.01$. The crossing value should be compared with
other Ising like systems: $B_{4}(T_{c})=1.60$ for the 3D Ising model Fenz  _et
al._ (2007) and $B_{4}(T_{c},H_{c})=1.86$ for spin ice with field $\vec{H}$
along the $[111]$ cubic axis Hamp  _et al._ (2015). We found that the value
depends on $\phi_{c}$, reducing to $\approx 1.60$ for $\phi_{c}=0.5$, with
crossing at $T=0.3672$ K, but in this case the crossing was not so accurately
defined. From this analysis we estimate $T_{c}=0.3675\pm 0.0003$ K. In Fig.
(6) we show the probability density function, $P(\phi)$ calculated at
$T_{c},\Delta_{c}$ which resembles qualitatively the universal function $P(M)$
for the magnetisation $M$ of the three dimensional Ising model at the critical
point Rummukainen  _et al._ (1998); Cardozo and Holdsworth (2016) and is
centred on $\phi=0.5$. The universality class of the critical point is
discussed further below through a dynamical finite size scaling analysis and
the measurement of the fluctuation-dissipation ratio.

Figure 6: Locating the critical point. (Top) The Binder cumulant $B_{4}$ for
the emergent order parameter $\varphi=\phi-\phi_{c}$ for fixed $\nu=4.35$ K
and $\Delta=\Delta_{c}=2.03745$ K and $\phi_{c}=0.42$ at the
$\phi_{1}:\phi_{0}$ phase boundary (see text). (Bottom) Probability density
$P(\phi)$ at the critical end point, $\nu=4.35$ K, $\Delta_{c}=2.03745$,
$T_{c}=0.36752$.

###  II.5 The multicritical region

How the wings meet in the multicritical region is a rather subtle question.
The intersection of the five phases on the $\Delta=0$ plane at a single penta-
critical point is unlikely, as the plane is characterised by two variables $T$
and $\nu$ only. This allows the system to tune to a tri-critical point in
which both the quadratic and quartic terms in an expansion of the free energy
in $\varphi$ are zero Cardy (1996). However, a penta-critical point would
require the annulation of the sixth order term which, without a third
parameter would be accidental. In the model studied here, emergent from the
DSI, the monopole and double monopole costs are fixed: $\mu_{2}=4\nu$.
Floating $\mu_{2}$ away from this value could allow the tuning necessary to
establish penta-criticality but the evidence presented below suggests that in
our case the wings meet in two stages which indeed maintains the tri-
criticality of the BC1 model.

Figure 7: Low temperature ($T\ll\nu$): (Left) $\nu=\nu^{\ast}$, (Right)
$\nu<\nu^{\ast}$. (0) monopole fluid ($\phi_{0}$), (1) monopole crystal
($\phi_{1}$), (2) double monopole crystal ($\phi_{2}$) .

The five phases do however meet at $T=0,\Delta=0,\nu=\nu^{\ast}$. A phase
boundary between $\phi_{2}$ and $\phi_{0}$ then rises from the five phase
intercept along

$4(\nu^{\ast}-\nu)=s_{0}T,$ (10)

where we have again assumed a constant entropy approximation, valid for
$T\ll\nu^{\ast}$. Within this approximation the Pauling entropy of the spin
ice vacuum gives a finite slope away from $T=0$ which takes the system away
from $\phi_{1}$ and this phase is suppressed everywhere in the $\Delta=0$
plane except the special point at $T=0$. This can be seen in detail by
analysis of the three different free energies. As a consequence, $\Delta-T$
planes for $\nu=\nu^{\ast}$ and $\nu\leq\nu^{\ast}$ take the form shown in
Fig. (7) at low temperature. In the latter case there is a finite temperature
order disorder transition between $\phi_{2}$ and $\phi_{0}$ along the
$\Delta=0$ axis, ensuring that the $\phi_{2}$, $\phi_{1}$ and $\phi_{0}$
phases meet at a triple point for finite $\Delta$. The slopes of the phase
boundaries, $\eta_{ij}$ can be estimated from eqn. (6):
$\eta_{21}=(s_{1}-s_{2})=0.262$, $\eta_{10}=(s_{0}-s_{1})=0.143$ and
$\eta_{20}=\frac{1}{2}(s_{0}-s_{2})=0.203$ and the triple point, which is
allowed because of the linear dependence between the three boundary curves,
occurs at

$T_{t}=\frac{2(\nu^{\ast}-\nu)}{2s_{1}-s_{0}}.$ (11)

As $|\nu-\nu^{\ast}|$ increases, the $\Delta=0$ transition temperature
increases until at the tricritical point the transition changes from $1^{st}$
to $2^{nd}$ order, at which point the line structure in Fig. (7b) will have
evaporated through critical end points.

Heating up to the critical end points should therefore lead to the wings
meeting in two stages with three separate tri-critical points, with all lines
meeting tangentially Taufour  _et al._ (2016). Two of these being the critical
termination of the triple points for finite positive and negative $\Delta$ and
the third, a classic tri-critical point separating ordered and disordered
phases for $\Delta=0$.

This scenario can be compared with that of the BC2 model. In this case, the
same five phase intercept occurs at $T=0$ but the $\phi_{2}$ : $\phi_{0}$
phase boundary now rises vertically as the entropy of both phases approach
zero as $T$ goes to zero. However, at the level of mean field and pair
approximation calculations Plascak  _et al._ (1993); Lara and Plascak (1998) a
small sliver of $\phi_{1}$ appears at higher temperatures, stabilised by the
entropy of spin fluctuations. The $\phi_{1}$ : $\phi_{2}$ boundary ends at a
critical point in the $\bar{\Delta}=0$ plane as shown in Fig. 8. This suggests
that the tri-critical point of the BC1 model is again maintained with this
time, separate intercepts onto the central plane for the two wings for
positive and for negative $\bar{\Delta}$.

Figure 8: Symmetry breaking in the $S=2$ Blume Capel model: The
$\bar{\Delta}=0$ plane of the $S=2$ Blume Capel model Lara and Plascak (1998).

The undershoot and overshoot of the wing interceptions in the dumbbell and BC2
models illustrates the accidental nature of penta-criticality for this set of
parameters and strongly suggests that a generalised model with independent
$\mu$ and $\mu_{2}$ could be tuned to include a penta-critical point.

##  III Dynamic Scaling at a Critical End Point

###  III.1 Critical slowing down

Along the lines of critical end points there are divergent time scales
associated with the diverging correlation lengths and critical slowing down.
In Fig. (9) we show the evolution of the auto-correlation function

$C_{\phi}(t)=\frac{\left<\phi_{i}(t)\phi_{i}(0)\right>-\left<\phi_{i}(t)\right>\left<\phi_{i}(0)\right>}{\left<\phi_{i}(0)^{2}\right>-\left<\phi_{i}(0)\right>^{2}},$
(12)

with Metropolis Monte Carlo time as the critical end point for $\nu=4.35$ K,
$\Delta_{c}=2.03745$ is approached along the temperature axis. Note that in
eqn. (12) we study the critical dynamics using the local spin-spin
autocorrelation function, which is distinct from the autocorrelation function
of the global order parameter $\phi$. The spin autocorrelation function is
statistically easier to access, but it also captures the critical slowing
down. The data shows decay of correlations at equilibrium for a system of size
$L=12$.

As the transition is approached from above the correlation time increases and
$C_{\phi}(t)$ develops a powerlaw decay with exponent $\alpha\approx-1/2$, out
to a maximum of the order of $10^{3}$ Monte Carlo steps per dumbbell. The best
power law is observed for a temperature $T\approx 0.369$ K, higher than the
$T_{c}=0.3675$ estimated from analysis of the Binder cumulant.

Within the critical region time scales and length scales are bridged via the
dynamical critical exponent Hohenberg and Halperin (1977) $z$. The correlation
time $\tau$ diverges with the correlation length $\xi$ as

$\frac{\tau}{\tau_{0}}=\left(\frac{\xi}{a}\right)^{z}.$ (13)

Hence, as the spatial correlation function for the local order parameter in
dimension $d$ scales with distance in the critical region as
$r^{-(d-1-\eta)}$, with $d-\eta$ the anomalous dimension of the universality
class, one expects dynamical scaling of the form $C_{\phi}(t)\sim
t^{-(d-1-\eta)/z}$. Taking $\eta\simeq 0.0363\ll 1$ and $z\approx 2$, which
should be the case for local dynamics in the three dimensional Ising
universality class, one finds an exponent $\alpha\approx 1/2$ as observed. The
shift in effective transition temperature away from the Binder crossing point
is expected and is due to finite size effects.

Figure 9: Critical correlations at equilibrium: Autocorrelation function for
the monopole crystal $C_{\phi}(t)$ vs. $t$ for temperatures close to the
critical temperature. Solid black line shows power law decay with exponent
$-1/2$.

The cut off of the power law is compatible with the finite size cut off of
$\xi$: $\xi_{max}<L^{\prime}=\frac{4La}{\sqrt{3}}$. Taking $z\approx 2$,
$L=12$ and microscopic time $\tau_{0}$ equal to one Monte Carlo time step
indeed gives a cut off to the critical scaling of the order of $10^{3}$
Metropolis time steps. Below the critical temperature the time correlation
function develops a plateau which decays at longer time scales. This is
consistent with a change of regime in the dense crystalline phase where decay
of correlations is due to the creation and propagation of monopole holes
Jaubert (2015).

###  III.2 Kibble-Zurek scaling

A more quantitative picture of the emergent universality class of the critical
end point can be achieved by following the dynamical Kibble-Zurek Kibble
(1980); Zurek (1985) scaling protocol proposed in [Hamp  _et al._ , 2015]. In
this scenario the field like scaling variable is swept in time through a cycle
with characteristic time scale $\tau_{Q}$:

$\lambda(t)=\frac{\Delta(t)-\Delta_{c}}{\Delta_{c}}=\lambda_{0}\sin\left(\frac{\pi
t}{2\tau_{Q}}\right),$ (14)

with temperature fixed at $T_{c}$. Far from the critical point the equilibrium
time scale is small compared with $\tau_{Q}$ so that the evolution is
adiabatic but as the critical point is approached the equilibrium time scale
diverges. As a consequence, at a given point in each cycle the system falls
out of equilibrium creating hysteresis loops in the thermodynamic observables,
whose magnitude depends on sweep time.

Figure 10: Hysteresis loops: $\phi$ . vs. $\Delta$ for different sweep rates.
System size $L=8$. The $\phi(t)$ is a configurational average of $500$
samples, each starting at equilibrium at $t=0$.

In Fig. 10 we show the evolution of $\phi$ with $\Delta$ at $T_{c}$ for a
system of size $L=8$ and $\lambda_{0}=0.1$ for different $\tau_{Q}$.
Hysteresis loops centred on $\phi=\phi_{c}$ indeed appear and their amplitude
falls to zero as $\tau_{Q}$ increases.

Following eqn. (13), the correlation time diverges along the field axes as
$\tau\sim\tau_{0}\lambda^{-z\nu_{\lambda}}$ with $\nu_{\lambda}$ the field
driven correlation length exponent, so that $\tau$ and sweep time are related
through eqn. (14). The crossover from adiabatic to out of equilibrium response
occurs around the point $\left|\frac{d\tau}{dt}\right|=1$, which fixes a
characteristic Kibble-Zurek time scale,
$t_{KZ}=\tau_{0}\left(\frac{\tau_{Q}}{\tau_{0}}\right)^{z\nu_{\lambda}/(z\nu_{\lambda}+1)}$.

In the critical region the fall from equilibrium of the emergent order
parameter $\varphi$ is captured by the dynamical scaling hypothesis Hohenberg
and Halperin (1977)

$\varphi(t)=\left(\frac{1}{t_{KZ}}\right)^{\frac{D_{\lambda}}{z}}{\cal{G}}\left(\frac{t}{t_{KZ}}\right),$
(15)

where $D_{\lambda}$ is the scaling dimension of the field and ${\cal{G}}(x)$
is a scaling function allowing for data collapse for different data sets. For
short range systems, up to and including the upper critical dimension
$D_{\lambda}=\frac{d\nu_{\lambda}-1}{\nu_{\lambda}}$, while in the Gaussian
regime $D_{\lambda}=1$.

In Fig. (11) we show the Kibble-Zurek scaling collapse for
$\Delta\varphi=\varphi(\lambda)\uparrow-\varphi(\lambda)\downarrow$, the
difference in order parameter values on an up and down swing of the cycle. We
find a convincing collapse using known values for the three dimensional Ising
universality class Pelissetto and Vicari (2002) and local stochastic dynamics
Wansleben and Landau (1991), $\nu_{\lambda}=0.4$, $z=2$. We do not have access
to large enough system sizes or high enough resolution on our data to
distinguish between three dimensional XY and Ising universality classes but
the collapse shown is superior to that found using Gaussian exponents. Hence,
as in [Hamp  _et al._ , 2015] for the critical point observed for spin ice in
a $[111]$ field, we can exclude the possibility of the long range Coulomb
interactions influencing the the universal fluctuations.

Figure 11: Kibble-Zurek scaling: $\Delta\varphi$ weighted by
$t_{KZ}^{D_{\lambda}/z}$ as a function of $\left(\frac{t}{t_{KZ}}\right)$ for
$\nu_{\lambda}=0.4$ and $z=2$.

It is worth remarking that, however accurate the data, the field scaling
Kibble-Zurek protocol cannot unambiguously establish Ising universality, as
the procedure accesses only one of the two static scaling dimensions,
$D_{\lambda}$; $\nu_{\lambda}$ being independent of the second dimension
$D_{\epsilon}$. This yields
$D_{\lambda}=\frac{1}{\delta\nu_{\lambda}}=\frac{\beta}{\nu}$, where exponents
have their usual meaning Cardy (1996), establishing weak universality only
Suzuki (1974). This in principle allows for variation of $\beta$ and $\nu$
within the weak universality constraint Taroni  _et al._ (2008). A thermal
Kibble-Zurek protocol would fix the two static exponents through the presence
of both $D_{\lambda}$ and $\nu$ although one would then have a three parameter
fit ($D_{\lambda},\nu,z)$ for a single expression.

###  III.3 Aging and fluctuation-dissipation ratio

A further remarkable consequence of the diverging time scale at the critical
point is that if the system is suddenly quenched from a high temperature to
$T_{c}$, it will not reach equilibrium within the time window offered by
experiments or simulation. As a result, systems quenched to criticality
display universal aging properties, reported in an extensive literature
Godrèche and Luck (2000); Godreche and Luck (2000); Berthier  _et al._ (2001);
Henkel  _et al._ (2001); Henkel and Pleimling (2003); Calabrese and Gambassi
(2002a, b); Mayer  _et al._ (2003, 2004); Calabrese and Gambassi (2005);
Prudnikov  _et al._ (2015) showing explicitly that the tools developed in the
context of materials with slow glassy dynamics are highly relevant for aging
critical dynamics.

Two important properties emerge from the out-of-equilibrium dynamics. First,
the time correlation function in eqn. (12) is no longer time translationally
invariant, so that one needs to explicitly follow the dependence on the time
spent at criticality since the quench. As a result the system slowly ages
towards equilibrium in a manner reminiscent of disordered glassy systems
Bouchaud  _et al._ (1998). Second, the fluctuation-dissipation theorem (FDT)
which, in equilibrium connects linear response functions to time correlations
functions, is no longer valid. In glassy materials, violations of the FDT have
been found to take simple forms with appealing physical interpretations
Cugliandolo and Kurchan (1993); Cugliandolo  _et al._ (1997); Bouchaud  _et
al._ (1998); Crisanti and Ritort (2003). Studies of FDT violations in systems
quenched to criticality show that the deviations from the equilibrium relation
contains direct information about the universality class of the model Godrèche
and Luck (2000); Calabrese and Gambassi (2002a).

Inspired by these studies, we consider a numerical protocol in which the
temperature is instantaneously varied from $T=1000~{}K$ to
$T=T_{c}=0.3685~{}K$, and denote $t_{w}$ the “waiting time” spent at $T_{c}$
since the quench. We then define

$\tilde{C}_{\phi}(t,t_{w})=\frac{\left<\phi_{i}(t_{w})\phi_{i}(t)\right>-\left<\phi_{i}(t_{w})\right>\left<\phi_{i}(t)\right>}{\left<\phi_{i}(t)^{2}\right>-\left<\phi_{i}(t)\right>^{2}},$
(16)

where the waiting time dependance is now made explicit. Indeed, as with other
critical systems, we find that the time decay of the spin auto-correlation
function is not just a function of $t-t_{w}$ but now depends explicitly of
both times. We also define the linear response function associated with the
time correlation function in eqn. (16) as

$\chi_{\phi}(t,t_{w})=\frac{\partial\langle\phi_{i}(t)\rangle}{\partial\Delta_{i}(t_{w})},$
(17)

where $\Delta_{i}$ is the field conjugate to the local order parameter
$\phi_{i}$. We introduce the normalised response function
$\tilde{\chi}_{\phi}=\chi_{\phi}/(\left<\phi_{i}(t)^{2}\right>-\left<\phi_{i}(t)\right>^{2})$,
such that the equilibrium FDT reads $\tilde{\chi}=(1-\tilde{C}_{\phi})/T$.

In the aging regime following a quench, the FDT is not expected to be
satisfied, and it can generically be rewritten as

$\tilde{\chi}_{\phi}(t,tw)=\frac{X(t,t_{w})}{T}\left(1-\tilde{C}_{\phi}(t,t_{w})\right),$
(18)

which defines the fluctuation-dissipation ratio $X(t,t_{w})$ Cugliandolo and
Kurchan (1993). Physically, eqn. (18) is appealing as it has the same
mathematical form as in equilibrium, with the difference that the thermal bath
temperature is replaced by an effective temperature $T/X(t,t_{w})$ Cugliandolo
_et al._ (1997).

Figure 12:  FDT violations Parametric plot of the response function against
the correlation function, for various fixed times $t$ and using $t_{w}$ as a
running parameter after a quench at the critical point $T_{c},\Delta_{c}$. The
system size is $N_{0}=13824$ ($L=12$), and data are averaged over 500
independent quenches. The inset show limiting values for long times. The solid
lines show the expected universal value $X_{\infty}=0.38$ for the 3D Ising
universality class.

In Fig. (12) we display FDT violations by representing
$T\tilde{\chi}_{\phi}(t,t_{w})$ as a function of
$1-\tilde{C}_{\phi}(t,t_{w})$, for a fixed time $t$ and using $t_{w}$ as a
running parameter in the plot Mayer  _et al._ (2004). We repeat these
measurements for a series of $t$ values. In order to achieve statistical
accuracy, we adapt the most efficient Monte Carlo tools presented in refs.
Chatelain (2003); Ricci-Tersenghi (2003); Berthier (2007) to the dumbell
model.

The relevance of this representation is obvious as the slope of these curves
is a direct measure of the fluctuation-dissipation ratio, by virtue of eqn.
(12). Close to the origin, corresponding to short time differences $t-t_{w}$,
the equilibrium FDT is obeyed and the parametric response-correlation plot is
linear with slope given by the temperature $T$. In contrast, clear deviations
from the FDT are observed in the opposite limit of large time differences
$t-t_{w}$, with a fluctuation-dissipation ratio $X(t,t_{w})<1$. The physical
interpretation is that small-scale (and thus fast) fluctuations rapidly reach
thermal equilibrium and display equilibrium FDT, whereas large and slow
critical fluctuations retain their non-equilibrium nature and display FDT
violations, as seen in other critical systems Godrèche and Luck (2000);
Berthier  _et al._ (2001).

The limiting value of the fluctuation-dissipation ratio defined as

$X_{\infty}=\lim_{t_{w}\to\infty}\lim_{t\to\infty}X(t,t_{w})$ (19)

takes a finite value, specific to a particular universality class. In the
inset of Fig. 12, we compare the limiting value of the fluctuation-dissipation
ratio measured in our simulations to the known value, $X_{\infty}=0.38$
measured for the three dimensional Ising model Godreche and Luck (2000);
Prudnikov  _et al._ (2015). We find an excellent agreement with our data,
which again supports the idea that the critical end point is in this
universality class.

##  IV Comparison with spin ice in a $[111]$ field

At present, the only experimentally observable phase transition driven by
monopole ordering is that observed with magnetic field placed along the
$[111]$ crystal axis Sakakibara  _et al._ (2003); Aoki  _et al._ (2004);
Higashinaka  _et al._ (2004), $\vec{H}=\frac{H}{\sqrt{3}}(1,1,1)$. A field of
modest strength selects a subset of Pauling states with the moments of the
spins lying parallel to the field axis aligned in the field direction. The
system maintains a finite entropy related to configurations of the three spins
of each tetrahedron with components lying in the kagome planes perpendicular
to the field direction Isakov  _et al._ (2004b) (see Fig. (13)). On increasing
the field at low temperature, a first order transition is observed to a fully
ordered state of 3in-1out/3out-1in tetrahedra. As the temperature increases
the transition line terminates in a critical end point. In Fig. (14) we
reproduce data from Figure 4 of [Sakakibara  _et al._ , 2003], which reports
experiments on DTO. The figure shows the estimated phase diagram. The order of
magnitude of the applied field is $0.9$ Tesla and the critical temperature is
around $0.35$ K.

The transition has previously been successfully interpreted as a liquid-gas
like critical end point of a monopole crystalisation transtion Castelnovo  _et
al._ (2008) and in this sense is a close cousin of the transition separating
$\phi_{0}$ and $\phi_{1}$ discussed above. The main difference is that the
external field couples to both of the fragmented components of the magnetic
moments Brooks-Bartlett  _et al._ (2014), providing a staggered chemical
potential for the monopoles and introducing a preference for Dirac strings
oriented with the field. The field therefore breaks both the magnetic symmetry
and the monopole translational symmetry.

Figure 13: Monopole creation in a $[111]$ field:
(Upper panel) connected tetrahedra perpendicular to the $[111]$ axis form
kagomé planes of spins which are bases for alternating up and down
tetetrahedra. The arrow shows the direction of an applied field.
(Lower panel) Flipping the spins indicated (left) creates monopole paires with
broken translational symmetry. North pole (+) -red disc, South pole (-) blue
disc.

At the transition monopole pairs are created in abondance by flipping spins in
the kagome planes as illustrated in Fig. (13). The direct action of the field
on the charges is to provide a chemical potential gradient,
$\vec{\nabla}\mu_{\pm}=\pm Q\vec{H}$ so that, in addition to the energy scale
for monopole creation in zero field, there is also a contribution depending on
the direction of movement in the field. The chemical potential gradient alone
does not therefore provide a staggered energy profile. If one of the north
monopoles of Fig. (13) were to continue moving along the $\hat{z}$ axis it
would pick up energy at each step in the same manner. However, the constraints
of spin ice forbid this: movement between the kagome planes is blocked as, on
the magnetic plateau the spins joining the planes point in the wrong direction
to allow monopole movement between planes via a single spin flip. Preparing
the ground for this move requires flips of loops of spins at high energy cost
Castelnovo  _et al._ (2010) so that the monopoles are essentially confined to
two-dimensional strips perpendicular to the field axis.

The chemical potential gradient does provide a staggered energy landscape
within this confined space. The difference in potential energy for a (north)
monopole on an $A$ or a $B$ site of a kagome plane is
$2\Delta=Q\vec{a}.\vec{H}$ where $\vec{a}=\frac{a}{\sqrt{3}}(-1,-1,1)$ is a
lattice vector spanning the two sites. This yields
$\Delta=\left(\frac{mH}{3}\right)$, which is just the Zeeman energy of the
spin flip in the presence of the field.

Given the similarities, we can repeat the thermodynamic arguments of section
II C for modified phases $\phi_{0}^{\prime}$, the plateau phase with entropy
per tetrahedron $s_{o}^{\prime}=0.1616$ Udagawa  _et al._ (2002) and ordered
monopole crystal phase $\phi_{1}^{\prime}$ with entropy zero. From this, using
$m=9.87\mu_{B}$ Yavors’kii  _et al._ (2008) we predict a field for the
transition at zero temperature $H(T=0)=0.90$ Tesla and an initial slope for
the phase boundary $\frac{dH}{dT}=0.073$ Tesla K-1. We note that the observed
critical temperature for DTO Sakakibara  _et al._ (2003) is very close to our
calculated value $T_{c}=0.367$ K for the $\phi_{0}:\phi_{1}$ critical end
point. Taking this value and using the constant entropy approximation we find
$H_{c}=0.927$ T.

Our predicted phase diagram, shown in Fig. (14) is in quite remarkable
quantitative agreement with reference [Sakakibara  _et al._ , 2003]. However,
a word of caution is probably in order. As the entropy of the phases
$\phi_{0},\phi_{0}^{\prime}$ and $\phi_{1},\phi_{1}^{\prime}$ are different
there is no reason to expect such quantitative agreement between the two
critical temperatures. Indeed simulations of the transition using the dipolar
spin ice model den Hertog and Gingras (2000), while still in excellent
qualitative agreement with the experimental data show a significantly higher
critical temperature, $T_{c}=0.587$ K Castelnovo  _et al._ (2008); Hamp  _et
al._ (2015). However, quantitative modelling of experiments with the DSI at
such low temperatures requires corrections in the form of further neighbour
exchange terms Henelius  _et al._ (2016), which could also have significant
effects on the critical end point Castelnovo  _et al._ (2008). In general
these extra terms reduce the ordering temperature for symmetry breaking among
the Pauling states, compared with the original DSI model Melko  _et al._
(2001). As the dumbbell model has no such ordering transition these
corrections may play in its favour, but one could be forced to concede an
element of good fortune in this remarkable agreement. It would clearly be of
interest to pursue this subject in future research.

Figure 14: Monopole phase transition in a $[111]$ field: We reproduce data
from Figure 4 of [Sakakibara  _et al._ , 2003], which shows the thermal
variation of the transition field of Dy2Ti2O7 for $H$ parallel to $[111]$.
Solid circles (squares) denote the data points obtained for increasing
(decreasing) field sweeps. Open triangles are the averaged critical field
$H_{c}$, which show a nearly linear temperature variation with
$\frac{dH_{c}}{dT}=0.08$ Tesla K-1 (dashed line) at low temperature. The added
blue circles are our estimates of the field strength at $T=0$ and at $T_{c}$
with blue line the phase boundary given by the constant entropy approximation.

##  V Discussion

We have shown that the dumbbell model of spin ice has a rich phase diagram
with the double winged structure shown in Fig. (2). A key to its existence is
the presence of a first order line for the spin ice - AIAO transition in the
$\Delta=0$ plane Guruciaga  _et al._ (2014). The first order nature of the
transition ensures that the singularity survives application of a symmetry
breaking field giving symmetry sustaining transitions and the emergence of the
wings. The first order transition becomes second order via a tri-critical
point as discussed in detail in section II.5. Tri-critical behaviour with
first and second order sectors is common in frustrated magnetic systems
Champion  _et al._ (2002); Zhitomirsky  _et al._ (2012); Shahbazi and
Mortezapour (2008); Sadeghi  _et al._ (2015) and is related to the entropy of
fluctuations provided by the frustrated geometry. In the case of spin ice one
must go beyond the nearest neighbour spin ice model to generate a first order
transition as within this approximation the monopoles are non-interacting.
Ordering in this case is due uniquely to entropic considerations Guruciaga
_et al._ (2014) and can only be second order Castelnovo  _et al._ (2008).
Including the dipolar interactions in the spin model provides the emergent
monopoles with an energy versus entropy trade off which drives the transition
first order. However, truncating the Coulomb interaction beyond nearest
neighbour monopoles would not give a quantitative change to the phase diagram.

We have studied both the dumbbell model and the related Blume-Capel model, the
BC2. The five phases of the wings meet either in two stages for dumbbell, or
not at all for BC2, entering the $\Delta=0$ plane at two different values of
$\nu$ and $T$. This undershoot or overshoot is consistent with there being
only two independent variables on the plane. This could however be changed by
freeing the double monopole chemical potential from the fixed value,
$\mu_{2}=-4\nu$ of the present model. By tuning $\mu_{2}$ it should be
possible to find a parameter set $\nu$, $T$, $\mu_{2}$ for which the five
phases meet at a single penta-critical point. This corresponds, at the mean
field level to all terms up to and including order $\phi^{6}$ being zero in an
expansion of the free energy.

###  V.1 Liquid-gas, liquid-liquid and symmetry sustaining transitions

Liquid gas phase transitions have two defining characteristics.

Firstly they correspond to crossing lines of phase equilibria in a temperature
like - field like phase diagram, between two phases with the same symmetry. If
the line of transitions terminates at a critical end point, it is then
possible to move analytically from one phase to another by contouring this
special point. As a consequence, the only thing that defines the two separate
phases is the transition itself. There is broken symmetry at the transition,
but it is emergent, separating phase space into high and low density sectors
with the same microscopic symmetry.

The second is that the sustained symmetry is the highest allowed by the
Hamiltonian. The generic case is that of a fluid that changes from low to high
density through the control of temperature and pressure, or chemical potential
while maintaining continuous translational symmetry. In the quantum case,
temperature could be replaced by a coupling constant and thermal fluctuations
by quantum fluctuations, allowing for the transition from a long range
entangled quantum liquid, such as a quantum spin liquid Balents (2010) to a
classical paramagnet, or spin gas phase Savary and Balents (2013).

The first criterion is ubiquitous thermodynamics and can be generated for any
first order transition by the application of a field conjugate to the order
parameter characterising the transition. The second is a non-universal
property of strongly correlated systems, with the existence of the liquid-gas
singularity dependent on the microscopic properties of the model Barrat and
Hansen (2003).

The transitions discussed in this paper satisfy the first criterion but not
the second. In general Coulomb fluids on a bi-partite lattice do not offer a
liquid gas transition with the full discrete symmetry of the $T-\nu,\Delta=0$
plane Kobelev  _et al._ (2002). Rather, such a transition is usurped by
sublimation from low density fluid to crystal with broken Z2 translational
symmetry, as we have seen here in detail for the diamond lattice. The
transitions are therefore liquid-gas like in a weak sense: they are symmetry
sustaining but do not maintain the highest translational symmetry offered by
the diamond lattice. However, as the two criteria are equivalent from a
thermodynamic point of view, Blume-Capel type models and therefore spin ice
can be considered as generic systems for studying symmetry sustaining
phenomena, often occurring in liquids.

In particular, there has been much work on systems showing liquid-liquid phase
transitions. In molten phosphorous Katayama  _et al._ (2000), silicon Sastry
and Austen Angell (2003) or water Brovchenko  _et al._ (2005) for example
pressure takes the fluid from a low to high density liquid state via a first
order transition that terminates in a critical end point. The high density
transition often appears in a supercooled state as it is again usurped by
crystallisation in thermodynamic equilibrium. A characteristic of these
systems is the capacity to accommodate two kinds of local packing, open
(tetrahedral) and close packed. This can be modelled using two hard core
repulsion length scales Franzese  _et al._ (2001) but it is also proposed as
an emergent phenomenon due to frustration and inhomogeneities in simple fluids
Tanaka (2000). The BC2 model and hence spin ice clearly provides a generic
skeleton for this science. If passage from the $\phi_{0}$ to $\phi_{1}$ is
equivalent a liquid-gas phase transition then that from $\phi_{1}$ to
$\phi_{2}$ on one side of the double winged phase diagram of Fig. (2) is
thermodynamically equivalent to a liquid-liquid transition. Detailed
comparison with the models presented here could therefore provide new insight
into the necessary conditions for liquid-liquid transitions including the
possibility of liquid-liquid tri-criticality.

Slightly nearer to home, similar physics is observed in magnetic itinerant
electron systems under pressure. Both LaCrGe3 Kaluarachchi  _et al._ (2017)
and UGe2 Kotegawa  _et al._ (2011) show double winged phase diagrams as a
function of temperature, pressure and applied field with two ferromagnetic
phases extending out to finite field values. The ferromagnetic phase
transitions are symmetry sustaining in exact analogy with the transitions
presented in this paper, so that the BC2 type models again provide a skeleton
for this structure. Interestingly these materials provide experimental
examples of the two possible multi-critical regions discussed in section II.5,
confirming the accidental nature of the wing connections in the phase diagram.
In LaCrGe3 the wings meet in two stages for each field direction, as proposed
for spin ice, giving three distinct tri-critical points for positive and
negative characteristic fields and for $H=0$. In UGe2 on the other hand, the
wings meet the central plane separately, as is apparently the case for the BC2
model.

In LaCrGe3 the winged phase transitions are extrapolated to terminate at zero
temperature and finite field at a series of quantum critical points. This
prediction should be contrasted with antiferromagnetic BC2 type models for
which the lines of finite temperature critical end points extend out to
$\Delta=\pm\infty$. Here, as $\Delta$ becomes large the partitioning of north
and south poles on $A$ and $B$ sublattices becomes perfect, but the collective
interaction between charges of opposite sign still drives a liquid-gas like
discontinuity in the sublattice monopole density at finite temperature. It
would certainly be interesting to do more studies for the ferromagnetic case
including transverse spin fluctuations, the quantum case being accessible via
quantum Monte Carlo simulation.

###  V.2 Future experiments in frustrated magnetism

Motivation for this work has come in large part from experiments on the spin
ice material Ho2Ir2O7 (HIO)Lefrançois  _et al._ (2017). In this material both
the Ho3+ ions and the Ir4+ ions carry a magnetic moment and they sit on
interpenetrating pyrochlore structures. The moments of the Ir4+ ions order on
the scale of $100$ K into an AIAO structure which provides internal magnetic
fields which act in turn on the Ho3+ magnetic moments. In the monopole picture
the internal fields translate into the staggered chemical potential studied
here and proposed in [Brooks-Bartlett  _et al._ , 2014] as a mechanism for
separating the $\phi_{1}$ and $\phi_{2}$ monopole crystal phases and accessing
the $\phi_{1}$ phase. As temperature is lowered through the 1K range the Ho
sublattice continuously develops AIAO order with the ordered moment saturating
at 50 $\%$ of the total moment. The leftover moment gives correlated diffuse
scattering consistent with a Coulomb phase and the measured characteristics of
the powder sample are indeed consistent with the fragmented $\phi_{1}$ phase.

As the Ho sublattice shows no phase transition none of the winged structure
is, as yet observable directly in experiment. However, it is worth noting that
a different material in this series, Tb2Ir2O7 (TIO) settles into a ground
state with full AIAO order, that is into the $\phi_{2}$ phase as defined above
Lefrançois  _et al._ (2015). Its sister material Tb2Ti2O7 (TTO) is in some
sense spin ice like, falling close to the spin ice AIAO ($\Delta=0$) phase
boundary den Hertog and Gingras (2000). Hence, although TTO remains an enigma
Rau and Gingras (2018), the fact that TIO fully orders is completely
consistent with our logic. If one could chemically tune the values of $\nu$
and $\Delta$ from TIO to HIO one would pass through the $\phi_{2}:\phi_{1}$
phase boundary on the way. Once at the values corresponding to HIO, heating up
one could hit a further phase boundary but no further transition is required
from symmetry arguments as the transitions are symmetry sustaining. Although
the planes of first order transitions do not lie perpendicular to the
$\nu-\Delta$ planes they do fall with a very steep slope, the inverse of eqn.
(6). Hence for an accidental value of $\Delta$ it is quite likely that a
thermal trajectory would maintain the system well away from the phase
boundaries. We propose that this is the case for HIO.

The above conclusion immediately begs the question of if it is possible to
shift the value of $\Delta$ experimentally. One possibility would be to put
materials such as HIO or its dysprosium counterpart under pressure. High
pressure would presumably change both the strength of the internal fields and
the monopole chemical potentials $\mu$ and $\mu_{2}$ which are combinations of
exchange and dipole interactions Castelnovo  _et al._ (2008). One might expect
that increasing the pressure would have the effect of increasing the scale of
the antiferromagnetic exchange, therefore reducing the scale of $\mu$, while
at the same time increasing the scale of $\Delta$, moving the system towards
the $\phi_{1}:\phi_{2}$ boundary, but the evolution could equally well be
counterintuitive and go in the opposite direction. One could also consider the
effects of chemical pressure through the chemical substitution of Ir4+ ions
with non-magnetic species such as Ti4+, or Ge4+ which has a smaller ionic
radius than its counterparts Zhou  _et al._ (2012). In order to hit one of the
phase boundaries, starting from HIO one would need to shift $\Delta$ and/or
$\nu$ on the Kelvin scale, that is on the scale of the exchange constants
themselves. These are challenging experiments that open the door to rich
theoretical and numerical problems and the present results provide a
motivating framework in which to work.

Given the steepness of the slope of the phase boundaries in Fig. (2), if one
did cross a first order plane by altering $\Delta$, further tuning to find the
critical end point to the plane should be straightforward, at least in
comparison, giving access to Kibble-Zurek scaling experiments as outlined in
section III and proposed for the critical point in a $[111]$ field Hamp  _et
al._ (2015). The prospect of doing Kibble-Zurek scaling experiments is
particularly appealing as critical slowing down gives very weakly diverging
time scales and so is difficult to access experimentally. For example, if the
microscopic time scale is a nanosecond, getting the divergence into the
millisecond range requires a correlation length of 1000 times the microscopic
length and a reduced temperature or field of order $10^{-6}$. Such high
precision can be avoided by finding systems with either long microscopic
length or time scales. Long microscopic length scales occur naturally in cold
atom systems, which has recently led to successful Kibble-Zurek type
experiments Corman  _et al._ (2014); Labeyrie and Kaiser (2016). Spin ice, on
the other hand is ideally suited because of its naturally long microscopic
time scales, for example around a millisecond for DTO Jaubert and Holdsworth
(2009) so that the plethora of critical points presented here could open the
door to many such dynamical experiments. Once accessed, both field like and
temperature like protocols are envisageable.

Our work also suggests that it could be interesting to extend to spin ice
materials the type of noise measurements that were previously performed in
spin glasses Hérisson and Ocio (2002, 2004) to simultaneously detect linear
susceptibilities and time correlation functions, in order to experimentally
access the fluctuation-dissipation ratio introduced in section III.3.

##  VI Conclusion

Spin ice materials and models have proven to be the source of rich emergent
science Bramwell and Gingras (2001); den Hertog and Gingras (2000); Isakov
_et al._ (2004a); Castelnovo  _et al._ (2008); Ryzhkin (2005); Fennell  _et
al._ (2009); Brooks-Bartlett  _et al._ (2014), widening the scope and interest
of frustrated magnetism and offering multiple avenues for novel research. In
particular the monopole picture, which simplifies a complex and strongly
interacting frustrated system to a level in which it can be addressed in
incomparable detail, has provided an unexpected controlled environment in
which to study Coulomb fluids both from a field theoretic and charge
perspective. We have exploited the full phase diagram of the emergent, on
lattice magnetolyte in which both monopoles and double charged monopoles play
important roles. In doing so, we have exposed a model system for multiple
phase transitions with wide ranging applications. These include fluids showing
liquid-liquid phase transitions Katayama  _et al._ (2000); Sastry and Austen
Angell (2003); Brovchenko  _et al._ (2005) and itinerant magnetic systems
under pressure Kaluarachchi  _et al._ (2017); Kotegawa  _et al._ (2011) as
well as extensive new applications within the field of frustrated magnetism.

## Acknowledgments

It is a pleasure to thank Claudio Castelnovo, Laurent de Forges de Parny,
Ludovic Jaubert, Vojtech Kaiser, Elsa Lhotel, Sylvain Petit, Lucile Savary and
Valentin Taufour for useful discussions. We also thank Vojtech Kaiser for
sharing numerical codes with us, Camille Scalliet for help with the figures
and Zenji Hiroi for authorising the reproduction of our Figure 14. This
project was supported in part (VR and PCWH) by ANR grant Listen Monopoles. THS
thanks ENS de Lyon for financial support during a research project. This work
was supported by a grant from the Simons Foundation (# 454933, L. Berthier).

## References

  * Harris  _et al._ (1997) M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Physical Review Letters 79, 2554 (1997).
  * Bramwell and Gingras (2001) S. T. Bramwell and M. J. P. Gingras, Science 294, 1495 (2001).
  * Isakov  _et al._ (2004a) S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Physical Review Letters 93, 167204 (2004a).
  * Castelnovo  _et al._ (2008) C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).
  * Ryzhkin (2005) I. A. Ryzhkin, Journal of Experimental and Theoretical Physics 101, 481–486 (2005).
  * Jaubert and Holdsworth (2009) L. D. C. Jaubert and P. C. W. Holdsworth, Nature Physics 5, 258 (2009).
  * Castelnovo  _et al._ (2011) C. Castelnovo, R. Moessner, and S. L. Sondhi, Physical Review B 84, 144435 (2011).
  * Brooks-Bartlett  _et al._ (2014) 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).
  * Kaiser  _et al._ (2015) V. Kaiser, S. T. Bramwell, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. 115, 037201 (2015).
  * Möller and Moessner (2006) G. Möller and R. Moessner, Phys. Rev. Lett. 96, 237202 (2006).
  * den Hertog and Gingras (2000) B. C. den Hertog and M. J. P. Gingras, Phys. Rev. Lett. 84, 3430 (2000).
  * Isakov  _et al._ (2005) S. V. Isakov, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 95, 217201 (2005).
  * Pauling (1935) L. Pauling, Journal of the American Chemical Society 57 (1935).
  * Jaubert and Udagawa (2019) L. D. C. Jaubert and M. Udagawa, eds.,  _Spin Ice_ (Springer, 2019).
  * Lefrançois  _et al._ (2017) E. Lefrançois, V. Cathelin, E. Lhotel, J. Robert, P. Lejay, C. V. Colin, B. Canals, F. Damay, J. Ollivier, B. Fåk, L. C. Chapon, R. Ballou, and V. Simonet, Nature Communications 8, 209 (2017).
  * Melko and Gingras (2004) R. G. Melko and M. J. P. Gingras, Journal of Physics: Condensed Matter 16, R1277 (2004).
  * Guruciaga  _et al._ (2014) P. C. Guruciaga, S. A. Grigera, and R. A. Borzi, Phys. Rev. B 90, 184423 (2014).
  * Borzi  _et al._ (2013) R. A. Borzi, D. Slobinsky, and S. A. Grigera, Phys. Rev. Lett. 111, 147204 (2013).
  * Katayama  _et al._ (2000) Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura, M. Yamakata, and K.-i. Funakoshi, Nature 403, 170 EP (2000).
  * Sastry and Austen Angell (2003) S. Sastry and C. Austen Angell, Nature Materials 2, 739 EP (2003).
  * Brovchenko  _et al._ (2005) I. Brovchenko, A. Geiger, and A. Oleinikova,  _The Journal of Chemical Physics_ , The Journal of Chemical Physics 123, 044515 (2005).
  * Kaluarachchi  _et al._ (2017) U. S. Kaluarachchi, S. L. Bud’ko, P. C. Canfield, and V. Taufour, Nature Communications 8, 546 (2017).
  * Kotegawa  _et al._ (2011) H. Kotegawa, V. Taufour, D. Aoki, G. Knebel, and J. Flouquet,  _Journal of the Physical Society of Japan_ , Journal of the Physical Society of Japan 80, 083703 (2011).
  * Plascak  _et al._ (1993) J. Plascak, J. Moreira, and F. s Barreto, Physics Letters A 173, 360 (1993).
  * Henelius  _et al._ (2016) P. Henelius, T. Lin, M. Enjalran, Z. Hao, J. G. Rau, J. Altosaar, F. Flicker, T. Yavors’kii, and M. J. P. Gingras, Phys. Rev. B 93, 024402 (2016).
  * Melko  _et al._ (2001) R. G. Melko, B. C. den Hertog, and M. J. P. Gingras, Phys. Rev. Lett. 87, 067203 (2001).
  * Cardy (1996) J. Cardy,  _Scaling and Renormalization in Statistical Physics_ (Cambridge University Press, 1996).
  * Blume (1966) M. Blume, Phys. Rev. 141, 517 (1966).
  * Capel (1966) H. W. Capel, Physica 32, 966 (1966).
  * Blume  _et al._ (1971) M. Blume, V. J. Emery, and R. B. Griffiths, Phys. Rev. A 4, 1071 (1971).
  * Landau and Lifshitz (1959) L. D. Landau and E. M. Lifshitz,  _Course of Theoretical Physics Vol. 5 (Statistical Physics Part 1)_ (Permagon, 1959).
  * Lara and Plascak (1998) D. P. Lara and J. A. Plascak, International Journal of Modern Physics B 12, 2045 (1998), https://doi.org/10.1142/S0217979298001198 .
  * Jaubert  _et al._ (2011) L. D. C. Jaubert, M. Haque, and R. Moessner, Physical review letters 107, 177202 (2011).
  * Jaubert  _et al._ (2013) L. D. C. Jaubert, M. J. Harris, T. Fennell, R. G. Melko, S. T. Bramwell, and P. C. W. Holdsworth, Phys. Rev. X 3, 011014 (2013).
  * Kaiser  _et al._ (2013) V. Kaiser, S. T. Bramwell, P. C. W. Holdsworth, and R. Moessner, Nature Materials 12, 1033 (2013).
  * Kaiser  _et al._ (2018) V. Kaiser, J. Bloxsom, L. Bovo, S. T. Bramwell, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. B 98, 144413 (2018).
  * Castelnovo and Holdsworth (2019) C. Castelnovo and P. C. W. Holdsworth,  _Spin Ice, Edited by L. D. C. Jaubert and M. Udagawa_ (Springer, 2019).
  * Henley (2010) C. L. Henley, Annual Review of Condensed Matter Physics 1, 179 (2010).
  * Nagle (1966) J. F. Nagle, Phys. Rev. 152, 190 (1966).
  * Maggs and Rossetto (2002) A. C. Maggs and V. Rossetto, Phys. Rev. Lett. 88, 196402 (2002).
  * Note (1) The value of $u(a)=2.88$ used here corresponds to a magnetic moment for the spins, $m=9.87\mu_{B}$, as deduced from crystal field calculations Yavors’kii  _et al._ (2008) rather than the $10\mu_{B}$ often quoted in the literature, which gives $u(a)=3.07$K and a corresponding difference in the energy scale for the phase diagram.
  * Binder (1981) K. Binder, Zeitschrift für Physik B Condensed Matter 43, 119 (1981).
  * Hamp  _et al._ (2015) J. Hamp, A. Chandran, R. Moessner, and C. Castelnovo, Phys. Rev. B 92, 075142 (2015).
  * Fenz  _et al._ (2007) W. Fenz, R. Folk, I. M. Mryglod, and I. P. Omelyan, Phys. Rev. E 75, 061504 (2007).
  * Rummukainen  _et al._ (1998) K. Rummukainen, M. Tsypin, K. Kajantie, M. Laine, and M. Shaposhnikov, Nuclear Physics B 532, 283 (1998).
  * Cardozo and Holdsworth (2016) D. L. Cardozo and P. C. W. Holdsworth, Journal of Physics: Condensed Matter 28, 166007 (2016).
  * Taufour  _et al._ (2016) V. Taufour, U. S. Kaluarachchi, and V. G. Kogan, Phys. Rev. B 94, 060410 (2016).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Jaubert (2015) L. D. C. Jaubert, SPIN 05, 1540005 (2015), https://doi.org/10.1142/S2010324715400056 .
  * Kibble (1980) T. Kibble, Physics Reports 67, 183 (1980).
  * Zurek (1985) W. H. Zurek, Nature 317, 505 EP (1985).
  * Pelissetto and Vicari (2002) A. Pelissetto and E. Vicari, Physics Reports 368, 549 (2002).
  * Wansleben and Landau (1991) S. Wansleben and D. P. Landau, Phys. Rev. B 43, 6006 (1991).
  * Suzuki (1974) M. Suzuki, Progress of Theoretical Physics 51, 1992 (1974).
  * Taroni  _et al._ (2008) A. Taroni, S. T. Bramwell, and P. C. W. Holdsworth, Journal of Physics: Condensed Matter 20, 275233 (2008).
  * Godrèche and Luck (2000) C. Godrèche and J. Luck, Journal of Physics A: Mathematical and General 33, 1151 (2000).
  * Godreche and Luck (2000) C. Godreche and J. Luck, Journal of Physics A: Mathematical and General 33, 9141 (2000).
  * Berthier  _et al._ (2001) L. Berthier, P. C. Holdsworth, and M. Sellitto, Journal of Physics A: Mathematical and General 34, 1805 (2001).
  * Henkel  _et al._ (2001) M. Henkel, M. Pleimling, C. Godreche, and J.-M. Luck, Physical review letters 87, 265701 (2001).
  * Henkel and Pleimling (2003) M. Henkel and M. Pleimling, Phys. Rev. E 68, 065101 (2003).
  * Calabrese and Gambassi (2002a) P. Calabrese and A. Gambassi, Physical Review E 66, 066101 (2002a).
  * Calabrese and Gambassi (2002b) P. Calabrese and A. Gambassi, Physical Review E 65, 066120 (2002b).
  * Mayer  _et al._ (2003) P. Mayer, L. Berthier, J. P. Garrahan, and P. Sollich, Physical Review E 68, 016116 (2003).
  * Mayer  _et al._ (2004) P. Mayer, L. Berthier, J. P. Garrahan, and P. Sollich, Physical Review E 70, 018102 (2004).
  * Calabrese and Gambassi (2005) P. Calabrese and A. Gambassi, Journal of Physics A: Mathematical and General 38, R133 (2005).
  * Prudnikov  _et al._ (2015) V. V. Prudnikov, P. V. Prudnikov, E. A. Pospelov, and A. N. Vakilov, Physics Letters A 379, 774 (2015).
  * Bouchaud  _et al._ (1998) J.-P. Bouchaud, L. F. Cugliandolo, J. Kurchan, and M. Mezard, Spin glasses and random fields , 161 (1998).
  * Cugliandolo and Kurchan (1993) L. F. Cugliandolo and J. Kurchan, Physical Review Letters 71, 173 (1993).
  * Cugliandolo  _et al._ (1997) L. F. Cugliandolo, J. Kurchan, and L. Peliti, Physical Review E 55, 3898 (1997).
  * Crisanti and Ritort (2003) A. Crisanti and F. Ritort, Journal of Physics A: Mathematical and General 36, R181 (2003).
  * Chatelain (2003) C. Chatelain, Journal of Physics A: Mathematical and General 36, 10739 (2003).
  * Ricci-Tersenghi (2003) F. Ricci-Tersenghi, Physical Review E 68, 065104 (2003).
  * Berthier (2007) L. Berthier, Physical review letters 98, 220601 (2007).
  * Sakakibara  _et al._ (2003) T. Sakakibara, T. Tayama, Z. Hiroi, K. Matsuhira, and S. Takagi, Phys. Rev. Lett. 90, 207205 (2003).
  * Aoki  _et al._ (2004) H. Aoki, T. Sakakibara, K. Matsuhira, and Z. Hiroi, Journal of the Physical Society of Japan 73, 2851 (2004), http://dx.doi.org/10.1143/JPSJ.73.2851 .
  * Higashinaka  _et al._ (2004) R. Higashinaka, H. Fukazawa, K. Deguchi, and Y. Maeno,  _Journal of the Physical Society of Japan_ , Journal of the Physical Society of Japan 73, 2845 (2004).
  * Isakov  _et al._ (2004b) S. V. Isakov, K. S. Raman, R. Moessner, and S. L. Sondhi, Phys. Rev. B 70, 104418 (2004b).
  * Castelnovo  _et al._ (2010) C. Castelnovo, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 104, 107201 (2010).
  * Udagawa  _et al._ (2002) M. Udagawa, M. Ogata, and Z. Hiroi, Journal of the Physical Society of Japan 71, 2365 (2002), http://dx.doi.org/10.1143/JPSJ.71.2365 .
  * Yavors’kii  _et al._ (2008) T. Yavors’kii, T. Fennell, M. J. P. Gingras, and S. T. Bramwell, Phys. Rev. Lett. 101, 037204 (2008).
  * Champion  _et al._ (2002) J. D. M. Champion, S. T. Bramwell, P. C. W. Holdsworth, and M. J. Harris, Europhysics Letters 57, 93 (2002).
  * Zhitomirsky  _et al._ (2012) M. E. Zhitomirsky, M. V. Gvozdikova, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. 109, 077204 (2012).
  * Shahbazi and Mortezapour (2008) F. Shahbazi and S. Mortezapour, Phys. Rev. B 77, 214420 (2008).
  * Sadeghi  _et al._ (2015) A. Sadeghi, M. Alaei, F. Shahbazi, and M. J. P. Gingras, Phys. Rev. B 91, 140407 (2015).
  * Balents (2010) L. Balents, Nature 464, 199 (2010).
  * Savary and Balents (2013) L. Savary and L. Balents, Phys. Rev. B 87, 205130 (2013).
  * Barrat and Hansen (2003) J.-L. Barrat and J.-P. Hansen,  _Basic Concepts for Simple and Complex Liquids_ (Cambridge University Press, 2003).
  * Kobelev  _et al._ (2002) V. Kobelev, A. B. Kolomeisky, and M. E. Fisher, The Journal of Chemical Physics 116, 7589 (2002), http://dx.doi.org/10.1063/1.1464827 .
  * Franzese  _et al._ (2001) G. Franzese, G. Malescio, A. Skibinsky, S. V. Buldyrev, and H. E. Stanley, Nature 409, 692 EP (2001).
  * Tanaka (2000) H. Tanaka, Phys. Rev. E 62, 6968 (2000).
  * Lefrançois  _et al._ (2015) E. Lefrançois, V. Simonet, R. Ballou, E. Lhotel, A. Hadj-Azzem, S. Kodjikian, P. Lejay, P. Manuel, D. Khalyavin, and L. C. Chapon, Phys. Rev. Lett. 114, 247202 (2015).
  * Rau and Gingras (2018) J. G. Rau and M. J. P. Gingras, arXiv e-prints , arXiv:1806.09638 (2018), arXiv:1806.09638 [cond-mat.str-el] .
  * Zhou  _et al._ (2012) 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, Physical Review Letters 108, 207206 (2012).
  * Corman  _et al._ (2014) L. Corman, L. Chomaz, T. Bienaimé, R. Desbuquois, C. Weitenberg, S. Nascimbène, J. Dalibard, and J. Beugnon, Phys. Rev. Lett. 113, 135302 (2014).
  * Labeyrie and Kaiser (2016) G. Labeyrie and R. Kaiser, Phys. Rev. Lett. 117, 275701 (2016).
  * Hérisson and Ocio (2002) D. Hérisson and M. Ocio, Physical Review Letters 88, 257202 (2002).
  * Hérisson and Ocio (2004) D. Hérisson and M. Ocio, The European Physical Journal B-Condensed Matter and Complex Systems 40, 283 (2004).
  * Fennell  _et al._ (2009) T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell, Science 326, 415 (2009).