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

HAL
open science

# Multiferroicity in spin ice: Towards magnetic crystallography of Tb 2 Ti 2 O 7 in a field

L. D. C. Jaubert, R. Moessner

► To cite this version:

L. D. C. Jaubert, R. Moessner. Multiferroicity in spin ice: Towards magnetic crystallography of Tb 2 Ti 2 O 7 in a field. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2015, 91, pp.214422. 10.1103/PhysRevB.91.214422. hal-01542080

HAL Id: hal-01542080

https://hal.science/hal-01542080

Submitted on 19 Jun 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Multiferroicity in spin ice: towards a magnetic crystallography of Tb$_{2}$Ti$_{2}$O$_{7}$ in a field

L.D.C Jaubert¹ and R. Moessner²

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

²Max-Planck-Institut für Physik komplexer Systeme, 01187 Dresden, Germany

(Dated: February 16, 2015)

We combine two aspects of magnetic frustration, multiferroicity and emergent quasi-particles in spin liquids, by studying magneto-electric monopoles. Spin ice offers to couple these emergent topological defects to external fields, and to each other, in unusual ways, making possible to lift the degeneracy underpinning the spin liquid and to potentially stabilize novel forms of charge crystals, opening the path to a “magnetic crystallography”. In developing the general phase diagram including nearest-neighbour coupling, Zeeman energy, electric and magnetic dipolar interactions, we uncover the emergence of a bi-layered crystal of singly-charged monopoles, whose stability, remarkably, is strengthened by an external [110] magnetic field. Our theory is able to account for the ordering process of Tb$_{2}$Ti$_{2}$O$_{7}$ in large field for reasonably small electric energy scales.

By providing mechanisms for strong magneto-electric coupling, frustration has become a key ingredient in multiferroics [1–8]. While the search for high-temperature multiferroics is appealing for technological application such as memory devices [8], frustration opens a window on novel fundamental properties of magnetic matter at low temperature where even weak perturbations can play an important role. This holds especially for the collective behaviour of spin liquids in a wide range of compounds from rare-earth [9] and copper [10–12] oxides to organic Mott insulators [13] or iridates [14, 15].

In spin ice materials, the constraints imposed by frustration support an extensively degenerate ground state where magnetic fluxes are locally conserved [16]. Such flux conservation can be described as a divergence-free condition, categorizing the spin ice ground state as a Coulomb spin liquid by analogy with Maxwell’s electromagnetism [17–19], where excitations take the form of classical magnetic monopoles (Fig. 2.c-d) [20].

In addition to their magnetic properties, it has been recently theorized that magnetic monopoles could also carry an electric dipole moment [21] (Fig. 2.c). Here we shall investigate the multiple facets of such magneto-electric coupling, as an unexplored generic ordering process in rare-earth pyrochlores, able to lift the degeneracy of spin liquids and to manipulate topological excitations in frustrated magnets. Our results are double. First of all, we give a precise description of the mosaic of competing phases in our multiferroic spin ice model (Eq.(1)). We show how a ferromagnetic double-layer structure of monopoles (DL) is stabilized by electric dipolar interactions and even enhanced by a magnetic field in the [110] direction. We then use this double-layer structure as a signature of multiferroicity in rare-earth oxides able to account for recent experiments on the spin liquid candidate Tb$_{2}$Ti$_{2}$O$_{7}$ in a field.

The model – we consider classical Ising spins $\vec{S}$ aligned with their local easy-axes on the pyrochlore lattice supporting electric moments $\vec{P}$ induced by magnetoe

electric coupling [21] (Fig. 2), interacting via nearest neighbour spin coupling, Zeeman energy, magnetic and electric dipolar interactions

$$ \mathcal{H} = J \sum_{\langle ij \rangle} \vec{S}_i \cdot \vec{S}_j - \vec{h} \cdot \sum_i \vec{S}_i + D_m r_m^3 \sum_{i>j} \frac{\vec{S}_i \cdot \vec{S}_j - 3 (\vec{S}_i \cdot \vec{e}_{ij}) (\vec{S}_j \cdot \vec{e}_{ij})}{r_{ij}^3} + D_e r_e^3 \sum_{\alpha>\beta} \frac{\vec{P}_\alpha \cdot \vec{P}_\beta - 3 (\vec{P}_\alpha \cdot \vec{e}_{\alpha\beta}) (\vec{P}_\beta \cdot \vec{e}_{\alpha\beta})}{r_{\alpha\beta}^3} \quad (1) $$

<figure><img src="image_3.png" /><figcaption>FIG. 1. When $D_m = h = 0$, the electric dipoles stabilize a monopole double layer (DL - green, see Fig. 2.a), in competition with all in / all out order (AIAO - red), and the Coulomb spin liquid (yellow). The circles (crosses) are the transition (crossover) temperatures obtained from Monte Carlo simulations. In the hatched regions, even if $T = 0$ calculations confirm the energetic stability of the double-layer structure, the first order nature of the transition prevents full thermalization of the simulations. The solid lines are upper and lower mean field estimates of the boundary. See Appendix for details on simulations and calculations.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Electrically induced ground state of our multiferroic spin ice model, made of alternative bi-layers of positive (blue) and negative (red) magnetic charges stacked along a [001] axis. The α-chains (indicated by thick bonds) carry a saturated [110] magnetisation. (b-d) There are three different kinds of configurations for a given tetrahedron: 2 in - 2 out (vacuum of charge forming the Coulomb spin liquid), 3 in - 1 out (single magnetic charges carrying an electric moment $\vec{P}$ whose direction is dictated by the minority (here outward) spin and independent of the sign of the magnetic charge because of time-reversal symmetry, forming the double-layer structure), and 4 in (double charges forming the AIAO order).</figcaption></figure>

where i, j and α, β are respectively indices for magnetic spins on the pyrochlore lattice and electric dipoles on the diamond lattice. $r_m$ and $r_e = \sqrt{3}/2 r_m$ are the respective nearest neighbour distances. The nearest neighbour vector $\vec{e}_{ij}$, magnetic $\vec{S}_i$ and electric $\vec{P}_\alpha$ moments have unit length. The size of the moments $\mu$ and $p$ is included in the energy-scale prefactors

$$D_m = \frac{\mu_0 \mu^2}{4\pi r_m^3}, \quad D_e = \frac{p^2}{4\pi\epsilon_0 r_e^3}, \quad \vec{h} = \mu_0 \mu \vec{H} \quad (2)$$

where $\mu_0$ and $\epsilon_0$ are respectively the vacuum magnetic permeability and electric permittivity and $\vec{H}$ is the external magnetic field.

The Hamiltonian is studied via classical Monte Carlo simulations, using parallel tempering, worm and single-spin-flip Metropolis algorithm. The dipolar energies have been computed with the Ewald summation [22, 23], in absence of demagnetization factor in order to develop a sample-independent theory [24]. All spin configurations are given in the Appendix.

**Monopole Double Layer** – first of all, what happens in absence of magnetic interactions, i.e. $D_m = h = J = 0$? We find that electric dipoles induce a bi-layer structure of single charges with zero polarization and saturated magnetization along the [110] axis (Fig. 2.a). Because

the electric field is even under time reversal, the application of such magnetization has to be spontaneous. This configuration is unfrustrated at the nearest-neighbour level, whose contribution constitutes 96% of the total energy. To our knowledge, magneto-electric coupling is the first intrinsic mechanism favouring single magnetic charges down to zero temperature.

**Local chemical potential J** – in terms of monopoles, J plays the role of a chemical potential favoring the Coulomb spin liquid for $J < 0$ (vacuum of charge) and the AIAO double-charge crystal for $J > 0$ [25, 26], making single charges gapped topological excitations in both cases. However, the previously observed $J = 0$ double-layer structure turns out to be robust over a large range of values for $J/D_e \in [-2.08 : 0.69]$ (Fig. 1), raising the question on the nature of the mechanism able to stabilize such “excitations”.

For instance when $J < 0$, creating a pair of single charges out of the Coulomb spin liquid costs $|4J/3|$ while the energy gain is at most $-2D_e/3$, making such monopole-pair creation unfavorable for $J/D_e < -1/2$. An energetically stable cluster of bi-layered monopoles thus needs to get bigger and bigger as J decreases in order to minimize its surface-over-volume ratio. The need for this kind of nucleation process to seed and grow a cluster makes the transition first order and prevents full thermalization of the simulations in the vicinity of the extensively degenerate Coulomb spin liquid (see yellow/green hatched region in Fig. 1). As a consequence, for $J \approx -2D_e$, an experimental cooling down protocol would probably fall out-of-equilibrium; once the magnet enters the Coulomb spin liquid with a low density of monopoles, it will be difficult to nucleate a big enough cluster of magnetic charges to crystallize the double-layer structure. Such phenomena also exist for $J > 0$, but to a lesser extent because of the low entropy of the AIAO ordered phase (see red/green hatched region in Fig. 1).

Hence, the magneto-electric opportunity to stabilize monopole excitations comes at the cost of large (free) energy barriers and multiple metastable states, which can naturally account for strong out-of-equilibrium effects in pyrochlores. In order to build a comprehensive and experimentally relevant picture of the problem, let us now include magnetic dipolar interactions, before adding a magnetic field able to tune these energy barriers, and finally applying our theory to experiments.

**Long range magnetic dipolar interactions $D_m$** – since the electric polarization is coming from magneto-electric coupling, it would be improper to neglect the $D_m$ energy scale, especially if we keep in mind rare-earth materials with potentially large magnetic moments. In spin ice, magnetic dipolar interactions are responsible for the effective Coulomb interactions between monopoles [20]. The property of “projective equivalence” [17] ensures the
quasi-degeneracy of the Coulomb spin liquid in presence of magnetic dipolar interactions [23], which is only weakly lifted in favour of the 2 in - 2 out long range ordered dipolar spin ice (ODSI) state for $T \ll D_m$ [27, 28].

This is why, by favoring 2 in - 2 out configurations, opposing the proximity of same-charge monopoles and hindering long-range ferromagnetic order, the $D_m$ interaction seems very unfavorable to the DL phase. Indeed for negative $J$, the double-layer phase makes way for the ODSI at finite $D_m$ (Fig. 3). Remarkably, the projective equivalence, valid for magnetic dipolar interactions in the Coulomb spin liquid but not electric ones in the DL phase, makes the transition temperature one order of magnitude smaller from the DL to the ODSI phase.

However the nearest neighbour contribution of the $D_m$ term decreases the monopole chemical potential which is four times bigger for double charges than for single ones [20]. The counter-intuitive consequence is that magnetic dipolar interactions favour DL order for positive values of $J$ where it was absent at $D_m = 0$. Because magnetic Coulomb interactions in spin ice are four orders of magnitude smaller than between bare electric charges at the same distance, a relatively low-energy coupling such as $D_e$ is sufficient to counter-balance the repulsion between neighbouring magnetic charges, opening the path for novel crystal structures made of magnetic monopoles.

<figure><img src="image_5.png" /><figcaption>FIG. 3. Phase diagram of $J/D_e$ and $D_m/D_e$ obtained from Monte Carlo simulations. The vertical axis is the normalized transition temperature into the AIAO (red), double-layer (green) or ODSI (orange) phase. See Appendix for details.</figcaption></figure>

A [110] magnetic field $h$ – in the absence of magneto-electric couplings, such a field polarizes half of the spins along the $\alpha$-chains (Fig. 2.a) while magnetic dipolar interactions align the remaining spins along antiferromagnetically-ordered $\beta$-chains (see Appendix) [29–31]. This 2 in - 2 out state is noted ODSI$_{[110]}$.

By ordering the $\alpha$-chains the field hinders thermal fluctuations which improved the thermalization of simulations, even if the presence of modulated phases could not be completely ruled out for all boundaries of the four-dimensional parameter space. The $T = 0$ phase diagram of the simulated phases was then computed by Ewald summation (see Appendix and Fig. 4).

<figure><img src="image_6.png" /><figcaption>FIG. 4. Zero temperature phase diagram in a [110] field $h$ obtained from Ewald summation. The arrows show the evolution of the boundaries as $h$ increases. The DL phase is further stabilized by $h$ over the AIAO order, but restricted to small values of $D_m$ by the apparition of a monolayer phase of monopoles (ML, in cyan).</figcaption></figure>

Because of the intrinsic quasi-degeneracy of the 2 in - 2 out configurations, the boundary with the DL phase only barely shifts as the ODSI order quickly gives way to ODSI$_{[110]}$ upon increasing $h$. The AIAO phase, on the other hand, is suppressed by the [110] field and gives way to the DL phase. The electrically induced double-layer structure is thus strengthened by the [110] magnetic field to larger $J$ and up to $D_m \approx 1.66 D_e$. For $D_m > 1.66 D_e$, the effective Coulomb repulsion between same-sign monopoles breaks the DL phase in favour of a zincblende or monolayer (ML) structure of monopoles with saturated magnetisation along the [111] direction.

*Multiferroicity in rare-earth oxides* – our results provide a clear signature of what to look for in experiments and remarkably, this double-layer structure has indeed been previously observed in Tb$_2$Ti$_2$O$_7$ under an external [110] magnetic field [32, 33] !

But is Tb$_2$Ti$_2$O$_7$ a good candidate for our theory ? We believe so for three reasons. Firstly, the Ising anisotropy of Tb$^{3+}$ ions [34, 35] and the pinch points observed in polarized neutron scattering [36–38] are strong indications for underlying spin-ice physics. Also, Tb$_2$Ti$_2$O$_7$ possesses a giant magnetostriction [39, 40], especially along the [110] direction [41], and the stability of the low temperature spin liquid phase has been recently ascribed to spin-phonon hybridization of the excitations [42], via dynamical Jahn-Teller coupling [43]. Last but not least, while sample dependence seems to be an issue in this compound, the double-layer structure in a large [110] field
has been confirmed by independent experiments [32, 33] and is thus a robust feature of Tb$_2$Ti$_2$O$_7$. It should indeed be noted that the nature of the zero- and low-field phases of this material remain under debate – possibly because of light stuffing/dilution [44] – alternatively described as spin liquid [36, 37, 45–47] or glassy [48–50] with antiferromagnetic correlations [50, 51]. Such behaviour is reminiscent of another rare-earth pyrochlore, Yb$_2$Ti$_2$O$_7$: while being remarkably well parametrized under a high magnetic field [52], the zero field properties of Yb$_2$Ti$_2$O$_7$ noticeably vary between samples [53–56], also possibly due to light stuffing [54]. Given the present low-field uncertainty, our goal here is to propose an alternative scenario for the high field region and to put a new benchmark on the 15-year-old puzzle that is Tb$_2$Ti$_2$O$_7$.

Because of a complex single-ion crystal field [47, 57–59], the effective size of the Tb$^{3+}$ magnetic moments is not fixed, but $\mu = 6 \mu_B$ is a good estimate at low temperature [34] and high field [60]. With $r_m = 3.59$ Å, Eq. 2 gives $D_m = 0.48$ K. From our theory, the double-layer structure can then be stabilized for $D_e^0 > D_m/1.66 \approx 0.30$ K, corresponding to an electric moment of $\sim 2.10^{-31}$ C.m and a displacement of oxygen ions of $\sim 0.6$ pm, which are reasonable estimates for multiferroics [8, 21]. Using the parametrization of [61] with $J = 2.7$ K, our multiferroic spin ice model can also explain why the DL phase only appears at finite field [33, 50] (Fig. 5). In light of the sample-dependence issue, it is difficult to push the comparison further to low field, where the phase might be antiferromagnetic but probably not AIAO [50], and should be separated from the DL structure by a collective paramagnet up to 2 Tesla. This is where magnetostriction comes into play as a potential complementary facet of our model.

The [110] direction for a magnetic field has been shown to maximize magnetostriction in Tb$_2$Ti$_2$O$_7$, resulting in field-dependent oxygen displacement [41]. Hence, an external magnetic field stabilizes the DL structure not only by suppressing the AIAO order (Fig. 4), but possibly also by increasing the electric energy scale $D_e$. In that case, an even wider range of parameters (e.g. $J = 0.96$ K [62]) and perturbations will lead to a double-layer phase at high field. The possibility to include disorder, anisotropic [43, 63, 64] or next-nearest-neighbour interactions [65] and quantum fluctuations [65–67] to our multiferroic spin ice model opens a rich diversity of potential phases to account for the yet uncertain phase of Tb$_2$Ti$_2$O$_7$ in low field. In particular it is tempting to speculate whether the glassy behaviour observed in some samples [48–50] might be a consequence of the dynamically difficult nucleation process discussed in this paper, especially if Tb$_2$Ti$_2$O$_7$ lies close to a spin liquid phase [36, 37, 45–47, 68].

**Conclusion** – in summary, we have shown how magneto-electric coupling can lift the degeneracy of a

<figure><img src="image_7.png" /><figcaption>FIG. 5. Phase diagram parametrized for Tb$_2$Ti$_2$O$_7$ for $J = 2.7$ K, $D_m = 0.48$ K and $D_e = 0.32$ K. The [110] magnetic field aligns the $\alpha$ chains along the [110] direction, destroying the AIAO order in favour for the double-layer structure. All error bars are smaller than the dots except for the red/green hatched region where simulations were difficult to equilibrate.</figcaption></figure>

spin liquid by creating interactions between topological excitations. In spin ice these excitations condense into a bi-layered monopole crystal strengthened by a [110] magnetic field, a non-trivial example of “magnetic crystallography”. Our theory offers a simple and robust explanation for the ordering of Tb$_2$Ti$_2$O$_7$ in a large [110] field for a reasonably small electric energy scale $D_e$.

If we look at the diverse physics emerging from itinerant electrons coupled to spin ice (anomalous and spontaneous Hall effects in Nd$_2$Mo$_2$O$_7$ [69] and Pr$_2$Ir$_2$O$_7$ [70], non-Kondo resistivity minimum [71–74] and a new kind of quantum criticality [75] in iridates), we should expect the coupling to an additional, ferroic, degree of freedom to bring a new flavor to spin ice and spin liquids, both at equilibrium and dynamically [76]. Experiments on Dy$_2$Ti$_2$O$_7$ and Ho$_2$Ti$_2$O$_7$ already suggest the presence of magneto-electric effects [77–80], which could be enhanced or even qualitatively modified by doping and chemical pressure [81–83] or with the inclusion of an electric field. More generally, multiferroicity offers a promising mechanism to control topological defects in magnets. We hope our work will motivate further theoretical and experimental investigations of multiferroic effects in pyrochlores and spin liquids.

The authors are thankful to Pascal Quémerais for collaborations at an early stage of this project, and to Owen Benton, Bruce Gaulin, Isabelle Mirebeau, Yukitoshi Mottome, Karlo Penc and Oleg Tchernyshyov for useful discussions. This work was supported by funding from the Theory of Quantum Matter unit of the Okinawa Institute of Science and Technology Graduate University.
# APPENDIX

## Ground states energies

The T = 0 phase diagram of Figs. 1 and 4 in the paper was calculated based on the following energies per number of spins N. The spin configurations are given in figures 8, 9 and 10.

$$ \text{AIAO: } E = +4.09 D_m - J \quad (3) $$

$$ \text{DL: } E = 0.0455 D_m - h/\sqrt{6} - 0.692 D_e \quad (4) $$

$$ \text{ODSI: } E = -1.95 D_m + J/3 \quad (5) $$

$$ \text{ODSI}_{[110]}: E = -1.90 D_m - h\sqrt{6} + J/3 \quad (6) $$

$$ \text{ML: } E = -0.370 D_m - h\sqrt{6} \quad (7) $$

One sees immediately that the DL phase wins over the ML one for $D_e > 0.60D_m$, and that it is stable for $J/D_e \in [-2.07; 0.692]$ when $D_m = h = 0$. While the prefactors for the coupling J and the Zeeman term are straightforward to calculate, the remaining terms require Ewald summation for a precise estimate. However, many of them can actually be calculated analytically to a good approximation [20, 84]. Since it provides a useful insight into the analogy between dipoles and monopoles, we briefly explain the method in the next section. Please note these results are obtained in absence of demagnetization factor. We recover the same values as Yoshida et al. [30] for the magnetic energies of ODSI and ODSI$_{[110]}$ when including the demagnetization factor of a sphere in vacuum.

## Analytical calculation of the energies

Let $Q = n q_m$ denote the magnetic charge on a given diamond site where $q_m = 2\mu/r_e$ and $n \in \{-2, -1, 0, 1, 2\}$. In presence of magnetic dipolar interactions only, the energy cost $p_n$ to create a monopole of charge n is (see Supplementary Informations of [20])

$$ p_n = -\frac{8}{3} \left( 1 + \sqrt{\frac{2}{3}} \right) n^2 D_m \quad (8) $$

The magnetic interaction between charges is difficult to calculate for any random configuration, but if the system is charge ordered then it is possible to use the Madelung constant of the corresponding crystal structure in order to calculate its Coulomb energy [84]. Both AIAO and ML configurations are ordered in the zincblende structure with Madelung constant $\alpha_{zb} = 1.638$. With $N/2$ diamond sites, the Coulomb energy is

$$ U_n^c = -\frac{1}{2} \frac{N}{2} \alpha_{zb} n^2 \frac{\mu_0 q_m^2}{4\pi r_e} = -N \alpha_{zb} \frac{2}{3} \sqrt{\frac{2}{3}} n^2 D_m. \quad (9) $$

giving a total energy $E_n = U_n^c - p_n \frac{N}{2}$

$$ E_n = \frac{2N}{3} \left( 2 + (2 - \alpha_{zb}) \sqrt{\frac{2}{3}} \right) n^2 D_m \quad (10) $$

$$ = 1.53 n^2 D_m \quad (11) $$

where n = 1 for the ML and n = 2 for AIAO. This value should be compared with a vacuum of charges, i.e. the Coulomb spin liquid. Since the degeneracy of the Coulomb spin liquid is weakly lifted, the choice for a reference energy is somewhat arbitrary. With a lowest energy state at $-1.95D_m$ (ODSI) and the highest energy one at $-1.85D_m$ (fully saturated in the [001] direction), we choose $E_{\text{ref}} = -1.90D_m$ as reference energy and obtain

$$ \Delta E_1 = E_{\text{ML}} - E_{\text{ref}} = 1.53 D_m \text{ vs } E_1 = 1.53 D_m $$

$$ \Delta E_2 = E_{\text{AIAO}} - E_{\text{ref}} = 5.99 D_m \text{ vs } E_2 = 6.12 D_m $$

which are in remarkably good agreement. As for the double layer structure, we could not find the Madelung constant for such charge ordering. So we calculated it using Ewald summation for Coulomb interactions and obtained $\alpha_{DL} = 0.976$. This gives $E_{DL,\text{Madelung}} = 1.89D_m$, to be compared with $\Delta E_{DL} = E_{DL} - E_{\text{ref}} = 1.95D_m$, within 3% of error.

## Finite temperature simulations

Now that the zero temperature boundary can be determined exactly from equations (3) to (7), let us turn our attention to the finite temperature phase diagram, and in particular to the double-layer phase. All simulations were done with parallel tempering, usually with 1 mK difference between parallel temperatures. To further help thermalization, we have developed a variant of the worm algorithm for dipolar spin ice [24], adapted for both electric dipolar interactions and the presence of singly charged monopoles. Measurements were made for a given number of Monte Carlo steps (MCs) noted $t_{\text{max}}$.

We used two different equilibration processes to compute the error bars:

* (i) the system is slowly cooled down from high temperature to the temperature T of measurement during $t_{\text{max}}/10$; then it is thermalized at temperature T during $t_{\text{max}}/10$; for our model, this method provides a lower bound for the transition temperature.

* (ii) for any given set of parameters, the ensemble of ground states is known exactly (cf. the previous sections); the system is then quenched into one of the ground state configurations and thermalized at temperature T during $t_{\text{max}}/10$; this method is biased towards ordering and offers an upper bound for the transition temperature.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Convergence of the transition temperature as a function of measurement time $t_{max}$ for both equilibration processes (i – green) and (ii – blue). For all figures $D_e = 1$, $h = 0$ and $N = 432$. We confirm that even if the process is very slow, simulations converge to the same transition temperature indicated by the dashed line.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7. Finite size scaling of the transition temperature for equilibration process (i). For all figures $D_e = 1$, $h = 0$ and $L = \{3, 4, 5, 6\}$. The number of pyrochlore sites is $N = 16L^3$. The dotted line is the best fit of the form $a + b/L^c$, where $a, b, c$ are fitting parameters.</figcaption></figure>
These two values provide the error bars plotted on Figs. 1 and 5 of the main text for a system of size $N = 1024$ spins. When not visible, the error bars are smaller than the symbols. For the 3-dimensional plot of Fig. 3 of the main text (no error bars), we only used the equilibration process (i) for a system of size $N = 432$, which provides a lower estimate of the stability of the DL phase. For Figs. 1, 3 and 5 of the main text, we used $t_{max} = 10^6$ MCS, except in the double layer phase of Fig. 1 where $t_{max} = 10^7$ MCS.

In Fig. 6, we show how these two processes converge to the same value for four different sets of parameters ordering in the double-layer phase (the most difficult ordering process in our simulations). Because of the double long-range interactions and the very long time of thermalization, big system sizes are difficult to simulate: Fig. 7 displays how the transition temperature converges to a finite value with increasing linear system size $L$.

## Husimi tree

However close to the boundaries, especially with the Coulomb spin liquid, if the equilibration process (ii) always orders in its ground state configuration, the process (i) might not be able to find the true ground state and will be dominated by the neighbouring phase. This is what happens in the hatched regions of Fig. 1 in the main text. In that case, we cannot rely solely on simulations to determine the finite temperature phase diagram and an analytical approach becomes necessary.

The out-of-equilibrium region between the DL and AIAO is rather narrow, which is why we shall focus on the broader one between the Coulomb spin liquid (CSL) and the double layer, by estimating the free energy of the two phases for $D_m = h = 0$ (cf. yellow/green hatched region in Fig. 1 of the main text).

Let us first consider the double layer phase. According to Eq. 4, its internal energy is $U_{DL} = -0.692D_e$. As for the entropy $S_{DL}$, since the transition is strongly first order, fluctuations can be neglected in a first approximation when compared with the Coulomb spin liquid of extensive degeneracy. Thus the free energy is $F_{DL} = U_{DL} - T S_{DL} = -0.692$ if we fix $D_e = 1$.

As for the Coulomb spin liquid, since $D_m = h = 0$, it corresponds to the canonical nearest neighbour spin ice model with electric interactions only between singly-charged monopoles. To obtain an upper and lower estimate of the boundary between the Coulomb spin liquid and the DL phase, we consider the two following cases

• 1) electric interactions between the dilute monopoles are modeled by an effective chemical potential: the biggest interaction energy gained by a pair of monopoles being $-2D_e/3$, we estimate an effective chemical potential of $-D_e/3$ per monopole.

• 2) electric interactions are neglected: since they tend to lower the interacting energies of the dilute monopoles, this should give an upper estimate of the transition temperature.

The free energies of both cases can be calculated in the Husimi tree approximation

$$F_{CSL 1} = -\frac{T}{2} \log \left[ \frac{4 e^{\beta D_e / 3} + e^{2\beta J} + 3 e^{-2\beta J / 3}}{2} \right] \quad (12)$$

$$F_{CSL 2} = -\frac{T}{2} \log \left[ \frac{4 + e^{2\beta J} + 3 e^{-2\beta J / 3}}{2} \right] \quad (13)$$

Both free energies reproduce the asymptotic limits of the Coulomb spin liquid entropy, namely the Pauling ($T \to 0^+$) and paramagnetic ($T \to +\infty$) entropies. When compared with $F_{DL}$, equations (12) and (13) provide respectively the lower and upper solid line of Fig. 1 of the main text for $J \in [-2.08 : -1.62]$. The estimated extent of the boundary shift due to electric interactions is consistent, within error bars, with our simulations.

[1] T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura, Nature 426, 55 (2003).

[2] N. Hur, S. Park, P. Sharma, J. Ahn, S. Guha, and S. Cheong, Nature 429, 392 (2004).

[3] H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev. Lett. 95, 057205 (2005).

[4] I. A. Sergienko and E. Dagotto, Phys. Rev. B 73, 094434 (2006).

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

[6] D. I. Khomskii, Journal of Magnetism and Magnetic Materials 30, 1 (1980).

[7] S.-W. Cheong and M. Mostovoy, Nature Materials 6, 13 (2007).

[8] Y. Tokura, S. Seki, and N. Nagaosa, Reports on Progress in Physics 77, 076501 (2014).

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

[10] M. P. Shores, E. A. Nytko, B. M. Bartlett, and D. G. Nocera, Journal of the American Chemical Society 127, 13462 (2005), http://dx.doi.org/10.1021/ja053891p.

[11] Y. Okamoto, H. Yoshida, and Z. Hiroi, Journal of the Physical Society of Japan 78, 033701 (2009), http://dx.doi.org/10.1143/JPSJ.78.033701.

[12] J. A. Quilliam, F. Bert, R. H. Colman, D. Boldrin, A. S. Wills, and P. Mendels, Phys. Rev. B 84, 180401 (2011).

[13] K. Kanoda and R. Kato, Annual Review of Condensed Matter Physics 2, 167 (2011).

[14] K. A. Modic, T. E. Smidt, I. Kimchi, N. P. Breznay, A. Biffin, S. Choi, R. D. Johnson, R. Coldea, P. Watkins-Curry, G. T. McCandless, J. Y. Chan, F. Gandara, Z. Islam, A. Vishwanath,
<figure><img src="image_10.png" /><figcaption>FIG. 8. Monopole Double Layer (DL) (top) and Monopole Mono Layer in a [110] field (ML) (bottom). All tetrahedra carry a single magnetic charge. The green arrows are the electric polarization carried by each monopole. The thick lines are the α-chains of spins along the [110] direction. There is no global polarization but there is a saturated magnetization along the α-chains in both configurations. For a 3 in - 1 out configuration (resp. 3 out - 1 in), the minority spin is the outward (resp. inward) spin. In the ML phase, each minority spin (circled in orange) is a minority spin for both adjacent monopoles.</figcaption></figure>

<caption>FIG. 9. Ordered Dipolar Spin Ice (ODSI) (top) and Ordered Dipolar Spin Ice in a [110] field (ODSI$_{[110]}$) (bottom). All tetrahedra are in the 2 in - 2 out configuration, i.e. without electric polarization. For the ODSI phase, all spins within a (001) plane point in the same direction. The magnetisation from one (001) plane to the next rotates by $\pi/2$. There is no global magnetisation. The ODSI$_{[110]}$ configurations have the same structure as for ODSI except that all spins on the α-chains (thick violet bonds) point in the same direction, giving rise to a saturated magnetization.</caption>

A. Shekhter, R. D. McDonald, and J. G. Analytis, Nature Communications 5 (2014), 10.1038/ncomms5203.

[15] T. Takayama, A. Kato, R. Dinnebier, J. Nuss, and H. Takagi, ArXiv e-prints (2014), arXiv:1403.3296 [cond-mat.str-el].

[16] M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Physical Review Letters 79, 2554 (1997).

[17] S. V. Isakov, R. Moessner, and S. L. Sondhi, Physical Review Letters 95, 217201 (2005).

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

[19] C. Henley, Annual Review of Condensed Matter Physics 1, 179 (2010) [arXiv:0905.0910].

[20] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).

[21] D. Khomskii, Nature Communications 3, 904 (2012).

[22] S. W. de Leeuw, J. W. Perram, and E. R. Smith, Proceedings of the Royal Society of London A 373, 27 (1980).

[23] B. C. den Hertog and M. J. P. Gingras, Physical Review Letters 84, 3430 (2000).

[24] R. G. Melko and M. J. P. Gingras, Journal of Physics-Condensed Matter 16, R1277 (2004).

[25] R. Moessner, Physical Review B 57, R5587 (1998).

[26] S. T. Bramwell and M. J. Harris, Journal of Physics-Condensed Matter 10, L215 (1998).

[27] R. Siddharthan, B. S. Shastry, A. P. Ramirez, A. Hayashi, R. J. Cava, and S. Rosenkranz, Physical Review Letters 83, 1854 (1999).
<figure><img src="image_11.png" /><figcaption>FIG. 10. **All in / All out** (AIAO). There is no electric polarization and no global magnetization.</figcaption></figure>

[28] R. G. Melko, B. C. den Hertog, and M. J. P. Gingras, Physical Review Letters **87**, 067203 (2001).

[29] R. Higashinaka, H. Fukazawa, and Y. Maeno, Physical Review B **68**, 014415 (2003).

[30] S. Yoshida, K. Nemoto, and K. Wada, Journal of the Physical Society of Japan **73**, 1619 (2004).

[31] T. Fennell, O. A. Petrenko, B. Fak, J. S. Gardner, S. T. Bramwell, and B. Ouladdiaf, Physical Review B **72**, 224411 (2005).

[32] J. P. C. Ruff, B. D. Gaulin, K. C. Rule, and J. S. Gardner, Phys. Rev. B **82**, 100401 (2010).

[33] A. P. Sazonov, A. Gukasov, I. Mirebeau, and P. Bonville, Phys. Rev. B **85**, 214420 (2012).

[34] M. J. P. Gingras, B. C. den Hertog, M. Faucher, J. S. Gardner, S. R. Dunsiger, L. J. Chang, B. D. Gaulin, N. P. Raju, and J. E. Greedan, Physical Review B **62**, 6496 (2000).

[35] H. Cao, A. Gukasov, I. Mirebeau, P. Bonville, C. Decorse, and G. Dhalenne, Phys. Rev. Lett. **103**, 056402 (2009).

[36] T. Fennell, M. Kenzelmann, B. Roessli, M. K. Haas, and R. J. Cava, Phys. Rev. Lett. **109**, 017201 (2012).

[37] S. Petit, P. Bonville, J. Robert, C. Decorse, and I. Mirebeau, Phys. Rev. B **86**, 174403 (2012).

[38] S. Guitteny, J. Robert, P. Bonville, J. Olivier, C. Decorse, P. Steffens, M. Boehm, H. Mutka, I. Mirebeau, and S. Petit, Phys. Rev. Lett. **111**, 087201 (2013).

[39] I. Aleksandrov, B. Lidskii, L. Mamsurova, M. Neigauz, K. Pigal'skii, K. Pukhov, N. Trusevich, and L. Chcherbakova, Journal of Experimental and Theoretical Physics **89**, 2230 (1985).

[40] L. Mamsurova, K. Pigal'skii, and K. Pukhov, Journal of Experimental and Theoretical Physics **43**, 584 (1986).

[41] J. P. C. Ruff, Z. Islam, J. P. Clancy, K. A. Ross, H. Nojiri, Y. H. Matsuda, H. A. Dabkowska, A. D. Dabkowski, and B. D. Gaulin, Phys. Rev. Lett. **105**, 077203 (2010).

[42] T. Fennell, M. Kenzelmann, B. Roessli, H. Mutka, J. Olivier, M. Ruminy, U. Stuhr, O. Zaharko, L. Bovo, A. Cervellino, M. K. Haas, and R. J. Cava, Phys. Rev. Lett. **112**, 017203 (2014).

[43] P. Bonville, A. Gukasov, I. Mirebeau, and S. Petit,

Phys. Rev. B **89**, 085115 (2014).

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

[45] J. S. Gardner, S. R. Dunsiger, B. D. Gaulin, M. J. P. Gingras, J. E. Greedan, R. F. Kiefl, M. D. Lumsden, W. A. MacFarlane, N. P. Raju, J. E. Sonier, I. Swainson, and Z. Tun, Physical Review Letters **82**, 1012 (1999).

[46] J. S. Gardner, A. Keren, G. Ehlers, C. Stock, E. Segal, J. M. Roper, B. Fak, M. B. Stone, P. R. Hammar, D. H. Reich, and B. D. Gaulin, Physical Review B **68**, 180401 (2003).

[47] P. Bonville, I. Mirebeau, A. Gukasov, S. Petit, and J. Robert, Phys. Rev. B **84**, 184409 (2011).

[48] G. Luo, S. Hess, and L. Corruccini, Physics Letters A **291**, 306 (2001).

[49] Y. Yasui, M. Kanada, M. Ito, H. Harashina, M. Sato, H. Okumura, K. Kakurai, and H. Kadowaki, Journal of the Physical Society of Japan **71**, 599 (2002).

[50] K. Fritsch, E. Kermarrec, K. A. Ross, Y. Qiu, J. R. D. Copley, D. Pomaranski, J. B. Kycia, H. A. Dabkowska, and B. D. Gaulin, Phys. Rev. B **90**, 014429 (2014).

[51] Y. Yasui, M. Kanada, M. Ito, H. Harashina, M. Sato, H. Okumura, and K. Kakurai, Journal of Physics and Chemistry of Solids **62**, 343 (2001).

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

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

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

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

[56] L.-J. Chang, M. R. Lees, I. Watanabe, A. D. Hillier, Y. Yasui, and S. Onoda, Phys. Rev. B **89**, 184416 (2014).

[57] J. Zhang, K. Fritsch, Z. Hao, B. V. Bagheri, M. J. P. Gingras, G. E. Granroth, P. Jiramongkolchai, R. J. Cava, and B. D. Gaulin, Phys. Rev. B **89**, 134410 (2014).

[58] V. Klekovkina and B. Malkin, Optics and Spectroscopy **116**, 849 (2014).

[59] A. J. Princep, H. C. Walker, D. T. Adroja, D. Prab- hakaran, and A. T. Boothroyd, ArXiv e-prints (2015), arXiv:1501.04927 [cond-mat.str-el].

[60] A. P. Sazonov, A. Gukasov, I. Mirebeau, H. Cao, P. Bonville, B. Grenier, and G. Dhalenne, Phys. Rev. B **82**, 174406 (2010).

[61] M. Enjalran and M. J. P. Gingras, Physical Review B **70**, 174426 (2004).

[62] I. Mirebeau, P. Bonville, and M. Hennion, Physical Review B **76**, 184436 (2007).

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

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

[65] P. McClarty, O. Sikora, R. Moessner, K. Penc, F. Poll- mann, and N. Shannon, ArXiv e-prints (2014), arXiv:1410.0451 [cond-mat.str-el].

[66] H. R. Molavian, M. J. P. Gingras, and B. Canals, Physical Review Letters **98**, 157204 (2007).
[67] H. R. Molavian, P. A. McClarty, and M. J. P. Gingras, ArXiv e-prints (2009), arXiv:0912.2957 [cond-mat.stat-mech].

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

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

[70] Y. Machida, S. Nakatsuji, S. Onoda, T. Tayama, and T. Sakakibara, Nature **463**, 210 (2013).

[71] S. Nakatsuji, Y. Machida, Y. Maeno, T. Tayama, T. Sakakibara, J. van Duijn, L. Balicas, J. N. Millican, R. T. Macaluso, and J. Y. Chan, Phys. Rev. Lett. **96**, 087204 (2006).

[72] M. Sakata, T. Kagayama, K. Shimizu, K. Matsuhira, S. Takagi, M. Wakeshima, and Y. Hinatsu, Phys. Rev. B **83**, 041102 (2011).

[73] M. Udagawa, H. Ishizuka, and Y. Motome, Phys. Rev. Lett. **108**, 066406 (2012).

[74] G.-W. Chern, S. Maiti, R. M. Fernandes, and P. Wölfle, Phys. Rev. Lett. **110**, 146602 (2013).

[75] L. Savary, E.-G. Moon, and L. Balents, Phys. Rev. X **4**, 041027 (2014).

[76] A. Sarkar and S. Mukhopadhyay, Phys. Rev. B **90**, 165129 (2014).

[77] T. Katsufuji and H. Takagi, Phys. Rev. B **69**, 064422 (2004).

[78] M. Saito, R. Higashinaka, and Y. Maeno, Physical Review B **72**, 144422 (2005).

[79] D. Liu, L. Lin, M. F. Liu, Z. B. Yan, S. Dong, and J. M. Liu, Journal of Applied Physics **113** (2013), 10.1063/1.4793704.

[80] C. P. Grams, M. Valldor, M. Garst, and J. Hemberger, Nature Communications **5** (2014), 10.1038/ncomms5853.

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

[82] L. Lin, Z. Y. Zhao, D. Liu, Y. L. Xie, S. Dong, Z. B. Yan, and J. M. Liu, Journal of Applied Physics **113** (2013), 10.1063/1.4794129, 12th Joint MMM-Intermag Conference, Chicago, IL, JAN 14-18, 2013.

[83] Z.-C. Xu, M.-F. Liu, L. Lin, H. Liu, Z.-B. Yan, and J.-M. Liu, Frontiers of Physics **9**, 82 (2014).

[84] M. Brooks-Bartlett, S. Banks, L. Jaubert, A. Harman-Clarke, and P. Holdsworth, Phys. Rev. X **4**, 011007 (2014).