The multiple symmetry sustaining phase transitions of spin ice

TL;DR

This paper maps the full phase diagram of spin ice's dumbbell model, revealing complex phase transitions and structures, and connects these findings to experimental observations and potential future experiments.

Contribution

It provides a comprehensive phase diagram of the dumbbell model of spin ice, including new phases and structures, and links theoretical results to experimental systems and protocols.

Findings

Identified a double winged phase structure with five phases.

Connected phase transitions to experimental systems like Ho$_2$Ir$_2$O$_7$.

Proposed new experimental protocols based on the phase diagram.

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 Ho$_2$Ir$_2$O$_7$ 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.

Figures

Equations

$${\mathcal{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,$$

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

$$\mathcal{H}=\frac{u(a)}{2}\sum _{i\ne j}\left(\frac{a}{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,$$

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

$${\varphi }_0$$

$$\frac{d\Delta }{dT}=-\frac{s_i-s_j}{{\varphi }_i-{\varphi }_j},$$

$${\varphi }_0$$

$$B_4(T,\Delta )=\frac{⟨(\phi {)}^4⟩}{⟨(\phi {)}^2{⟩}^2}.$$

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

$$T_t=\frac{2({\nu }^{\ast }-\nu )}{2s_1-s_0}.$$

$$C_{\varphi }(t)=\frac{⟨{\varphi }_i(t){\varphi }_i(0)⟩-⟨{\varphi }_i(t)⟩⟨{\varphi }_i(0)⟩}{⟨{\varphi }_i(0{)}^2⟩-{⟨{\varphi }_i(0)⟩}^2},$$

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

$$\lambda (t)=\frac{\Delta (t)-{\Delta }_c}{{\Delta }_c}={\lambda }_0\sin \left(\frac{\pi t}{2{\tau }_Q}\right),$$

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

$${\tilde{C}}_{\varphi }(t,t_w)=\frac{⟨{\varphi }_i(t_w){\varphi }_i(t)⟩-⟨{\varphi }_i(t_w)⟩⟨{\varphi }_i(t)⟩}{⟨{\varphi }_i(t{)}^2⟩-{⟨{\varphi }_i(t)⟩}^2},$$

$${\chi }_{\varphi }(t,t_w)=\frac{\partial ⟨{\varphi }_i(t)⟩}{\partial {\Delta }_i(t_w)},$$

$${\tilde{\chi }}_{\varphi }(t,tw)=\frac{X(t,t_w)}{T}\left(1-{\tilde{C}}_{\varphi }(t,t_w)\right),$$

$$X_{\infty }=\underset{t_w\to \infty }{\lim }\underset{t\to \infty }{\lim }X(t,t_w)$$

Full text

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.

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.

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

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

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

[TABLE]

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$

[TABLE]

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$:

[TABLE]

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 [math], 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).

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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.

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.

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

[TABLE]

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

[TABLE]

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}$.

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

[TABLE]

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

[TABLE]

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.

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}$:

[TABLE]

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.

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)

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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).

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

[TABLE]

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.

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.

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).