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

HAL
open science

# Magnetic monopole dynamics in spin ice

L D C Jaubert, P C W Holdsworth

► To cite this version:

L D C Jaubert, P C W Holdsworth. Magnetic monopole dynamics in spin ice. Journal of Physics: Condensed Matter, 2011, 23, pp.164222. 10.1088/0953-8984/23/16/164222. [hal-01541666]

HAL Id: hal-01541666

https://hal.science/hal-01541666

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
# Magnetic Monopole Dynamics in Spin Ice

L. D. C. Jaubert* and P. C. W. Holdsworth†

October 6, 2010

*Max-Plack-Institut für Physik komplexer Systeme, 01187 Dresden, Germany.

†Université de Lyon, Laboratoire de Physique, École Normale Supérieure de Lyon, 46 Allée d'Italie, 69364 Lyon cedex 07, France.

October 6, 2010

## Abstract

One of the most remarkable examples of emergent quasi-particles, is that of the "fractionalization" of magnetic dipoles in the low energy configurations of materials known as "spin ice", into free and unconfined magnetic monopoles interacting via Coulomb's $1/r$ law [Castelnovo et. al., Nature, **451**, 42-45 (2008)]. Recent experiments have shown that a Coulomb gas of magnetic charges really does exist at low temperature in these materials and this discovery provides a new perspective on otherwise largely inaccessible phenomenology. In this paper, after a review of the different spin ice models, we present detailed results describing the diffusive dynamics of monopole particles starting both from the dipolar spin ice model and directly from a Coulomb gas within the grand canonical ensemble. The diffusive quasi-particle dynamics of real spin ice materials within "quantum tunneling" regime is modeled with Metropolis dynamics, with the particles constrained to move along an underlying network of oriented paths, which are classical analogues of the Dirac strings connecting pairs of Dirac monopoles.

## 1 Introduction

Spin ice materials[1, 2] form part of a series of rare-earth oxide insulator $R_2M_2O_7$ with space group $Fd\bar{3}m$ where $R^{3+}$ is a magnetic ($Dy^{3+}$ and $Ho^{3+}$) and $M^{4+}$ a non-magnetic ($Ti^{4+}$ or $Sn^{4+}$) ion. The cations sit on the vertices of two interpenetrating pyrochlore lattices formed by corner-sharing tetrahedra (see figure 1). The well-established members of this family are $Dy_2Ti_2O_7$, $Ho_2Ti_2O_7$ and $Ho_2Sn_2O_7$, the latter being less studied as it is only available in polycrystalline form. Other potential candidates are $Dy_2Sn_2O_7$ or $Pr_2Sn_2O_7$, while related members of the series with rather different properties include $Tb_2Ti_2O_7$ and $Tb_2Sn_2O_7$. The total angular momentum $J = L + S$ is a good quantum number for rare earth elements and the electronic ground states are determined by Hund's rules.

The free ion $Dy^{3+}$ (resp. $Ho^{3+}$) then has a 16-fold (resp. 17-fold) degeneracy that is lifted by the surrounding crystal field; the corresponding energy levels have been estimated by neutron time-of-flight spectroscopy [3], revealing an almost pure ground state doublet.

As the first excited states are separated from the doublet by an energy scale of $\Delta \approx 2 - 300K$ (see Table 2.), the low temperature behaviour of these materials can be approximated by classical Ising spins with a very large magnetic moment ($\approx 10\mu_B$), oriented along the line connecting the centers of two
<table><caption>Table 1: **Ground state electronic configuration** of the atom [R={Dy,Ho}], the free ion [R<sup>3+</sup>], and the corresponding quantum numbers *S*, *L* and *J*</caption><thead><tr><th></th><th>[R]</th><th>[R<sup>3+</sup>]</th><th>S</th><th>L</th><th>J</th></tr></thead><tbody><tr><td>Dy</td><td>[Xe] 4f<sup>10</sup> 6s<sup>2</sup></td><td>[Xe] 4f<sup>9</sup></td><td>5/2</td><td>5</td><td>15/2</td></tr><tr><td>Ho</td><td>[Xe] 4f<sup>11</sup> 6s<sup>2</sup></td><td>[Xe] 4f<sup>10</sup></td><td>2</td><td>6</td><td>8</td></tr></tbody></table>

<table><caption>Table 2: **Crystal field level:** The ground states |J, Jz⟩ (GS), the Landé factor g<sub>J</sub>, the estimated magnetic moment μ = g<sub>J</sub> J μ<sub>B</sub> and the energy level of the first excited state Δ. (g<sub>J</sub> = 1 + (J(J+1) - L(L+1) + S(S+1)) / (2J(J+1)))</caption><thead><tr><th></th><th>GS</th><th>g<sub>J</sub></th><th>μ</th><th>Δ (meV)</th><th>Δ (K)</th></tr></thead><tbody><tr><td>Dy</td><td>|15/2, ±15/2⟩</td><td>4/3</td><td>≈ 10 μ<sub>B</sub></td><td>~ 0.033</td><td>~ 380</td></tr><tr><td>Ho</td><td>|8, ±8⟩</td><td>5/4</td><td>≈ 10 μ<sub>B</sub></td><td>~ 0.020</td><td>~ 240</td></tr></tbody></table>

neighbouring tetrahedra and forming a local set of body centered cubic axes:

$$ \mathbf{S}_1 = \frac{\sigma_1}{\sqrt{3}} \begin{pmatrix} -1 \\ -1 \\ +1 \end{pmatrix}, \quad \mathbf{S}_2 = \frac{\sigma_2}{\sqrt{3}} \begin{pmatrix} +1 \\ +1 \\ +1 \end{pmatrix}, \quad \mathbf{S}_3 = \frac{\sigma_3}{\sqrt{3}} \begin{pmatrix} +1 \\ -1 \\ -1 \end{pmatrix}, \quad \mathbf{S}_4 = \frac{\sigma_4}{\sqrt{3}} \begin{pmatrix} -1 \\ +1 \\ -1 \end{pmatrix}. \quad (1) $$

$\sigma_i = \pm 1$ is an Ising pseudo-spin with chosen convention that $\sigma_i = 1$ corresponds to the moment pointing out of a *down* tetrahedron (see figure 1).

The physics of spin ice is captured to a good approximation by an effective model with nearest neighbour ferromagnetic interactions between the moments which, together with the strong crystal field, gives rise to a frustrated ferromagnetic system:

$$ \mathcal{H} = -3J_{\text{eff}} \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j = J_{\text{eff}} \sum_{\langle i,j \rangle} \sigma_i \sigma_j. \quad (2) $$

Harris & al. [1] showed that the low energy phase space of this system is identical to the Pauling model for ice [4], which shows extensive ground state entropy in which states satisfy the Bernal-Fowler ice rules and which has been measured, both in ice ([5]) and spin ice ([6]). The nearest neighbour spin ice (NNSI) Hamiltonian also maps onto an Ising antiferromagnet for the pseudo spins [7, 8] (see equation (2)); a model first derived by Anderson to describe spinel materials [9]. The equivalence between these three nearest neighbour models is illustrated in figure 2.

The Bernal-Fowler ice rules requiring two magnetic moments to enter and two to leave each tetrahedron constitute a topological constraint, with the result that the Pauling states make up an ensemble of gauge invariant topological sectors with U(1) symmetry. A consequence of such constraints in this and in
<figure><img src="image_3.png" /><figcaption>Figure 1: **Pyrochlore lattice:** The spins are located on the corner of every tetrahedra and are fixed along their local [111] axis represented by the dashed lines. There are two types of tetrahedra that we shall call the *down* tetrahedra (bottom left with the four spins) and the *up* tetrahedra. Each *down* tetrahedra is connected to four *up* ones, and vice-versa. The cube represents a unit cell with 8 tetrahedra (four of each kind) and 16 spins, and defines the [100] (x), [010] (y) and [001] (z) axes. We introduce the length of the unit cell $a \approx 10$ Å as well as the distance between nearest neighbour $r_{nn} = \frac{\sqrt{3}}{4} a \approx 3.5$ Å and between the centres of two connected tetrahedra $r_d = \frac{\sqrt{3}}{4} a \approx 4.3$ Å. The smallest close loop encompasses 6 spins (see green dotted loop).</figcaption></figure>

related systems is the development of algebraic correlations [10, 11, 12, 13, 14]: consider the magnetisation of a tetrahedron $\alpha$: $M(\mathbf{r}_\alpha) = \sum_{i \in \alpha} S_i$ and define a coarse-grained field $M(\mathbf{r}) = \frac{1}{V} \sum_{\alpha \in V} M(\mathbf{r}_\alpha)$ averaged over a volume $V$ large enough to make $M(\mathbf{r})$ a smoothly varying function. The proof of algebraic correlations requires a two-stage argument based on the frustrated 2 in - 2 out ground states. On this manifold of states, the Gibbs free energy $G$ is purely entropic and can be expressed as a functional of $\{M(\mathbf{r})\}$. For a tetrahedron respecting the ice-rules, $M(\mathbf{r}_\alpha)$ can only take 6 values proportional to $(\pm 1, 0, 0), (0, \pm 1, 0)$ or $(0, 0, \pm 1)$. Hence computing the entropy density for a volume $V$ is similar to the problem of a random polymer made of monomers fixed on a cubic lattice and whose extremities are separated by $M(\mathbf{r})$. A large value of $|M(\mathbf{r})|$ requires most of the spins/monomers to be aligned along this direction, allowing only a relatively small number of configurations. On the other hand, zero magnetisation fixes a minimum number of degrees of freedom. Hence the Gibbs free energy and its Fourier transform are expressed to lowest order and up to a constant as

$$ G[\{M(\mathbf{r})\}] = \frac{T}{V} \int d^3\mathbf{r} \frac{1}{2} K |M(\mathbf{r})|^2 \quad (3) $$

$$ \Leftrightarrow \tilde{G}[\{\tilde{M}(\mathbf{k})\}] = T \sum_{\mathbf{k}} \frac{1}{2} K |\tilde{M}(\mathbf{k})|^2 \quad (4) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Models equivalent to nearest neighbour spin ice: The different mappings between Anderson’s model for antiferromagnetic spinels, spin ice and water ice: the white spheres are the hydrogen atoms and the red one is the oxygen. A spin pointing inside (resp. outside) the tetrahedron corresponds to a down (resp. up) spin in the Anderson model and a short covalent bond (resp. long H-bond) for water ice.</figcaption></figure>

which gives a temperature independent Gaussian probability $e^{-G/T}$ and $\langle \tilde{M}_{\mu}(-\mathbf{k})\tilde{M}_{\nu}(\mathbf{k}) \rangle = \delta_{\mu\nu}/K$. However this argument neglects the non-local influence of the ice-rules. If we regard the spins as local fluxes of magnetisation then these constraints impose a flux conservation for each tetrahedron (2 spins in and 2 spins out); from a coarse-grained point of view, it is equivalent to a divergence free condition.

$$ \nabla \cdot \mathbf{M}(\mathbf{r}) = 0 \quad \Leftrightarrow \quad \mathbf{k} \cdot \tilde{M}(\mathbf{k}) = 0 \qquad (5) $$

which imposes $\tilde{M}(\mathbf{k})$ to be orthogonal to $\mathbf{k}$ and thus

$$ \langle \tilde{M}_{\mu}(-\mathbf{k}) \tilde{M}_{\nu}(\mathbf{k}) \rangle = \frac{1}{K} \left( \delta_{\mu\nu} - \frac{k_{\mu} k_{\nu}}{|\mathbf{k}|^2} \right) \qquad (6) $$

where the Greek indices $\mu, \nu$ are Cartesian coordinates labels. After a final Fourier transform back to real space, we obtain the desired correlations

$$ \langle M_{\mu}(0) M_{\nu}(\mathbf{r}) \rangle \propto \frac{1}{K} \frac{\delta_{\mu\nu} - 3 \hat{r}_{\mu} \hat{r}_{\nu}}{r^3} \quad \text{where} \quad \hat{r}_{\mu} = \frac{r_{\mu}}{|\mathbf{r}|} \qquad (7) $$

This result deserves a few comments. First of all, we have lost the details of the pyrochlore lattice as soon as we used a coarse-grained field and this result can be generalised to models others than spin ice [10, 11, 12, 13, 14], the only hypothesis being

* a *Coulomb phase* gauge theory imposed by equation (4);

* a divergence free condition, consequence of the ice-rules constraints (5).

The magnetisation field $\mathbf{M}(\mathbf{r})$ is then analogous to a magnetic field without monopoles which explains the specific form of the correlations similar to those generated by a *dipole-dipole interaction*. One might think that the $1/r^3$ dependence, in a 3d crystal should lead to a logarithmic divergence in the susceptibility and logarithmic peaks in scattering function. This is not, however the case as the dipolar angular dependence of the correlation function results in a cancellation of the divergent part when integrated over all space. The correlations can be observed experimentally however through the development of narrowing *pinch points* [15] in the diffuse scattering profile [16, 17, 18].

It is clear that the NNSI model already contains extremely rich phenomenology related to the constrained Coulomb phase with $U(1)$ gauge symmetry, which is the basis for monopole physics. Most notably, single spin flip excitations create a pair of topological defects; pairs of tetrahedra with 3 in -
1 out and 3 out - 1 in . The defects are expected to have an effective entropic interaction [14], but once created can move at zero energy cost. They therefore posses some but not all the characteristics of neutral magnetic monopole pairs. In order to see the emergence of true magnetic quasi-particles one has to go beyond the NNSI through the inclusion of dipolar interactions, as is considered in the next section.

## 2 The dipolar spin ice model

Dipolar interactions are of particular importance [19, 20] for the rare earth ions Ho$^{3+}$ and Dy$^{3+}$ at the heart of spin ice physics: as the exchange coupling between the ions is due to 4f electrons, buried behind 5s and 5p layers$^{1}$, it is very weak (~1K) in comparison with other ferromagnets (for example $J \sim 10^2$ K for metallic iron). The magnetic moments in spin ice being very large, the dipolar interactions, that are negligible on short length scale for iron, are on the same energy scale as the exchange here. A quantitative description has therefore to take both interactions into account while to an excellent approximation moments can be described by Ising degrees of freedom orientated along the body centers of the tetrahedra. This separation of scales leads to the dipolar spin ice (DSI) Hamiltonian [20]

$$ \mathcal{H} = -J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j + D r_{nn}^3 \sum_{i>j} \left[ \frac{\mathbf{S}_i \cdot \mathbf{S}_j}{|\mathbf{r}_{ij}|^3} - \frac{3 (\mathbf{S}_i \cdot \mathbf{r}_{ij}) (\mathbf{S}_j \cdot \mathbf{r}_{ij})}{|\mathbf{r}_{ij}|^5} \right], \quad (8) $$

where $J, D$ and $r_{nn} \approx 3.5$ Å are respectively an antiferromagnetic exchange coupling, the dipole-dipole coupling and the nearest neighbour distance between rare earth ions, which has been shown to provide a comprehensive quantitative description of spin ice materials. Here, $J$ is a free parameter, while $D$ is defined through the dipolar interaction

$$ D = \frac{\mu_0 \mu^2}{4\pi r_{nn}^3} \approx 1.4K, \quad (9) $$

for Dy$_{2}$Ti$_{2}$O$_{7}$. A numerical estimate of $J$ has been obtained by comparison between Monte Carlo simulations and experimental data for the specific heat in Dy$_{2}$Ti$_{2}$O$_{7}$ [20] or of neutron scattering on Ho$_{2}$Ti$_{2}$O$_{7}$ [21] and Ho$_{2}$Sn$_{2}$O$_{7}$ [22]. Uncertainties in the estimated values for the Stanates remain high; around 50% for Ho$_{2}$Sn$_{2}$O$_{7}$ because of the unavailability of single crystals. The dipolar contribution turns out to be indispensable for understanding certain magnetic field induced transitions [23, 24, 25, 26] as well as the quantitative feature of magnetic monopole excitations [27]. However, a remarkable property of this Hamiltonian is that the extensive degeneracy of an effective nearest neighbour model is almost maintained, as the nearest neighbour contribution to the dipolar interaction is the same for all Pauling states and the long ranged part of the interaction is almost perfectly screened [28, 29]. The dipolar interaction between neighbouring spins 1 and 2 defined in figure 2 is given by:

$$ \begin{aligned} & D r_{nn}^3 \left[ \frac{\mathbf{S}_1 \cdot \mathbf{S}_2}{|\mathbf{r}_{12}|^3} - \frac{3 (\mathbf{S}_1 \cdot \mathbf{r}_{12}) (\mathbf{S}_2 \cdot \mathbf{r}_{12})}{|\mathbf{r}_{12}|^5} \right] \\ & = \frac{5D}{3} \sigma_1 \sigma_2. \end{aligned} \quad (10) $$

By symmetry, the result is the same for all pairs of spins in a tetrahedron. Hence the NN term of the dipolar interaction is equivalent to an exchange term and the nearest neighbour spin ice (NNSI) Hamiltonian is recovered as follows

$$ \mathcal{H} = -3J_{\text{eff}} \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j = J_{\text{eff}} \sum_{\langle i,j \rangle} \sigma_i \sigma_j \quad \text{where} \quad J_{\text{eff}} = D_{nn} + J_{nn} = \frac{5D}{3} + \frac{J}{3} \quad (11) $$
<table><caption>Table 3: **Nearest neighbour interactions:** Exchange, dipolar and effective anti-ferromagnetic couplings in Kelvin.</caption><thead><tr><th></th><th>D</th><th>D<sub>nn</sub></th><th>J</th><th>J<sub>nn</sub></th><th>J<sub>eff</sub></th></tr></thead><tbody><tr><td>Dy [20]</td><td>1.41</td><td>2.35</td><td>-3.72</td><td>-1.24</td><td>≈ 1.1</td></tr><tr><td>HoTi [21]</td><td>1.41</td><td>2.35</td><td>-1.65</td><td>-0.52</td><td>≈ 1.8</td></tr><tr><td>HoSn [22]</td><td>1.41</td><td>2.35</td><td>~ 1.0</td><td>~ 0.33</td><td>~ 2.7</td></tr></tbody></table>

Values for the spin ice materials is summarised in table 3. below.

The origin of this approximate screening, the so called projective equivalence, [29] can be seen physically by considering a further abstraction of the DSI to the dumbbell model [30], which highlights the remarkable symmetry properties of dipoles on the pyrochlore lattice. In this model the point dipoles of the DSI are extended into magnetically charged dumbbells, as shown in figure (3), whose ends sit on the centers of the tetrahedra, forming a diamond lattice of magnetically neutral sites for all Pauling states. Projective equivalence is clearly exact for this model and the error incurred in DSI, which turn out to be quadrupolar corrections only[29] is due to the difference between point and extended dipoles.

<figure><img src="image_5.png" /><figcaption>Figure 3: **Dumbbell model:** Each magnetic dipole is seen as two charges sitting on the vertices of the diamond lattice (dashed lines). 2 *in* - 2 *out* is a vacuum whereas 3 *in* - 1 *out* is a positive charge.</figcaption></figure>

## 3 Magnetic monopoles and classical "Dirac strings"

Within the dumbbell model, the ensemble of Pauling states provide a quasi-particle vacuum for neutral pairs of topological defects. Once created, the 3 *in* - 1 *out* and 3 *out* - 1 *in* defects break magnetic charge neutrality on each dual lattice, so that, when separated, the defects should interact via Coulombs law. Hence, one can see by construction that fractionalization of the dipolar dumbbells into deconfined magnetic Coulomb charges occurs when the ice rules are broken: in this sense these quasi-particles are magnetic monopoles. To show that, thanks to projective equivalence, the same is true to an excellent

<footnote n="1">¹The 6s layer is empty because Ho³⁺ and Dy³⁺ lost three electrons.</footnote>
approximation in the DSI, we have followed [30] and simulated a DSI system. The long range interactions have been treated using the Ewald method [31, 32, 23], where the slowly converging sum of all dipolar interactions between spins is replaced by two strongly convergent series. Periodic boundary conditions ensures correct convergence. Starting with a random 2 *in* - 2 *out* configuration, we can create and force the diffusion of a single pair of monopoles by flipping spins and then compute the energy of the system for each configuration. After averaging over a large number of initial microstates and paths of diffusion, we obtain the potential of interaction between two monopoles of opposite charges shown in figure 5. For a distance *r* = *x r*<sub>d</sub> where *r*<sub>nn</sub> is the distance between two vertices on the diamond lattice (see figure 1), the Coulomb potential can be written [30]

$$V(r) = -\frac{\mu_0 Q^2}{4\pi r} = \frac{V_{\min}}{x}, \quad \text{where} \quad Q = \frac{2\mu}{r_d}, \qquad (12)$$

and where

$$V_{\min} = -\frac{\mu_0 4 \mu^2}{4\pi \frac{3}{2} r_{nn}^2} \frac{1}{\sqrt{\frac{3}{2} r_{nn}}} = -\frac{8}{3} \sqrt{\frac{2}{3}} D \approx -3.07 \text{ K}, \qquad (13)$$

with, from equation (9), *D* ≈ 1.41 K. *V*<sub>min</sub> is the energy gained, with respect to the vacuum, by the creation of a nearest neighbour pair of monopoles. This potential *V* for the dumbbell model is compared with our numerical results for DSI in figure 5. As we have periodic boundaries, the pair of monopoles under consideration will also feel the presence of its *image charges* separated by at least one system size. While this is convenient for simulating a real system (as is done later in this paper), it leads to a finite size bias to the Coulomb energy here, where we consider a single pair of monopoles. This bias remains negligible for small *x* but can be observed at larger *x* and scales away with system size. This is illustrated in 5 for two system sizes *L* = 4 (green crosses) and *L* = 8 (red dots). For the larger system size the 1/*x* behaviour is clearly observed and the data respects the energy scale fixed by *V*<sub>min</sub> to a good approximation. The data plotted here are averaged over at least 100 simulations, so that the small differences between our results and Coulomb’s law reflect the limits of accuracy of the dumbbell model, which are of order 1/*r*<sup>5</sup>. The only fitting parameter in the figure is the reference energy scale for the DSI, which translates the numerical data along the *y*-axis. Similar results are obtained in reference [30]. However while they computed the energy of separation between monopoles for a single configuration, we have performed a statistical average, which produces some statistical noise, but confirms that this is a general property of the ensemble of constrained Pauling states.

<table><caption>Table 4. summarises the main parallels between the dumbbell and DSI models.</caption><thead><tr><th>Dipolar Spin Ice</th><th>Dumbbell / Monopoles</th></tr></thead><tbody><tr><td>2 in - 2 out ground state</td><td>vacuum</td></tr><tr><td>3 in - 1 out defects (local excitations)</td><td>quasi-particles (monopoles)</td></tr><tr><td>dipolar interactions between magnetic dipoles</td><td>Coulomb interactions between magnetic monopoles</td></tr><tr><td>pyrochlore lattice</td><td>diamond lattice</td></tr><tr><td>canonical ensemble</td><td>grand canonical ensemble</td></tr></tbody></table>

<caption>Table 4: Mapping from the dipolar spin ice to the dumbbell model</caption>
<figure><img src="image_6.png" /><figcaption>Figure 4: **Emergence of monopoles** a) The magnetic ions ($Ho^{3+}$ or $Dy^{3+}$) lie on the sites of the pyrochlore lattice and are constrained to the bonds of the diamond lattice (dashed lines). Local topological excitations 3 *in* - 1 *out* or 3 *out* - 1 *in* correspond to magnetic monopoles with positive (blue sphere) or negative (red sphere) charges respectively. b) The diamond lattice provides the skeleton for the network of Dirac strings with the position of the monopoles restricted to the vertices. The orientation of the Dirac strings shows the direction of the local field lines in **H**.</figcaption></figure>

The theoretical basis for magnetic monopoles in the DSI model is thus clearly established. We particularly stress that the scale of the Coulomb interaction is set by *magnetic* constants: the permeability of free space, $\mu_0$, and the nearest neighbour distance on the diamond lattice. Hence the quasi-particles behave quantitatively as *magnetically* charged particles. However, two fundamental differences still exists between these particles and Dirac monopoles: firstly they correspond to divergences in the magnetic intensity **H**, or magnetic moment **M**, rather than in the magnetic induction, so that $\nabla \cdot \mathbf{B} = \nabla \cdot (\mathbf{H}+\mathbf{M}) = 0$ and they do not require a modification of Maxwell's equations. Secondly, their magnetic charge is not quantified. The magnetic charge on a Dirac monopole is quantified through a quantum conjugation relation with the fundamental electronic charge [33, 34]. This is clearly not the case in spin ice and in fact the charge can be continuously varied by modifying the nearest neighbour distance between spins by applying external pressure [30]. However, on all length scales above the atomic scale, a 3 *in* - 1 *out* defect appears to be a local sink in **M** and therefore a source of field lines in **H**. Furthermore we can assign to them a positive or negative charge when immersed in a magnetic field [35] that is conserved through the creation and annihilation of neutral monopole pairs. As the Coulomb potential is not confining<fnref n="2" /> in 3d, these quasi-particles provide an example of fractionalisation in high dimensions: the point dipoles of the DSI appear to be separated into free poles of opposite sign that are free to wander independently through the system. They are therefore classical analogues of Dirac monopoles and spin ice provides the first 3d experimental realisation of deconfined magnetic monopoles.

In the theory of Dirac, monopoles of a neutral pair of charges are connected by a tensionless Dirac string of overturned dipoles, taken in the continuum limit. The string, which corresponds to nodes in the wave function for the pair of particles, should be unobservable. Classical analogues of Dirac strings exist in DSI and in the dumbbell model [30], which correspond to the magnetic moments in the

<footnote n="2">²having a Coulomb gas on a *lattice* prevents the divergence for $r \to 0^+$. The energy of deconfinement is finite and equals to $|V_{\min}|$ on average.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 5: **Coulomb potential:** Effective interaction (in Kelvin) between a pair of monopoles in the dipolar spin ice model as a function of the distance between them in units of nearest neighbour distance between monopoles. The system sizes are 4 × 4 × 4 unit cells, i.e. 1024 spins (green crosses) and 8 × 8 × 8 unit cells, i.e. 8192 spins (red dots). The solid line is the Coulomb potential expected from the dumbbell model $V(x) = V_{min}/x$. As we are on the diamond lattice, the smallest possible distance between monopoles is $x = 1$. The data are the average of over more than 100 samples.</figcaption></figure>

constrained 2 *in* - 2 *out* states. That is, by creating and separating a pair of monopoles one flips a series of moments, defining one possible Dirac string connecting them. Starting from a disordered 2 *in* - 2 *out* state, the path taken to create the monopoles and place them at positions A and B is not unique and the arbitrariness provides a classical analogue of the quantum uncertainty of Dirac monopole trajectories in space. Further to this analogy, within the dumbbell model, as long as the dumbbells are taken to have point magnetic charges at each end (the dumbbells can also be thought of as infinitesimally thin magnetic needles, touching at the centres of the tetrahedra), the magnetic flux associate with each dipole is confined to infinitely narrow elements joining the charges on each dumbbell, so that the network of classical Dirac strings is also invisible. This of course is a mathematical abstraction; in the DSI, we have point dipoles on the vertices of the tetrahedra and because of this difference (because of corrections to projective equivalence in fact), the string network *is* observable. Indeed, recent experiments have reported the observation of "Dirac strings"[17, 36]. As far as dynamics are concerned, the Dirac string network plays a crucial role: moving a monopole involves a spin flip, so that the quasi-particles leave a wake of re-arranged links in the network of Dirac strings as they move, with 3 *in* - 1 *out* and 3 *out* - 1 *in* defects moving in opposite directions over the network (see figure 4). This directionality makes the dynamics highly constrained compared with the diffusing electric charge of an equivalent electronic Coulomb gas. At equilibrium and in zero field this constraint simply renormalizes the diffusion constant but in a rapid quench the constraints lead to glassy behaviour [37] reminiscent of other kinetically constrained models[38]. We will return to the dynamics in finite field in the conclusion.

Given the accessibility of these magnetic quasi-particles, the development of an experimental signature is of vital importance and interest and has recently given rise to considerable experimental activity[17, 36,
18, 39]. In particular, recent experiments have measured the charge carried on free deconfined particles via the magnetic Wien effect [40, 41]. Here, we will concentrate on magnetic relaxation, available from measurements of magnetic response in weak ac magnetic fields [42, 43, 44, 45]. In the following section we will show how these data can be used as a direct signature for diffusive monopole dynamics.

# 4 Dynamics in spin ice

In figure (6) we show the evolution of the magnetic relaxation time scale for Dy$_2$Ti$_2$O$_7$ measured from bulk a.c. susceptibility measurement [45]. The data show three distinct features: a high temperature regime above 12 K, which, has been associated with thermal excitations of higher crystal field levels [46, 44, 47], an intermediate temperature plateau regime between 3 and 10K and a regime with rapidly increasing time scale below 3 K. Qualitatively similar behaviour has been observed in Holmium Titanate[48]. We will use this data to illustrate that the low and intermediate temperature dynamics of spin ice can be interpreted in terms of the stochastic dynamics of monopole quasi-particles. Given the separation of energy scales in spin ice materials we follow the assumption that the dynamics in the low temperature and the plateau regions are controlled by quantum tunneling events [44] taking the magnetic moments directly from one Ising position to another through the existence of non-zero off-diagonal components of the dipolar interaction from neighbouring spins, giving a small mixing of the magnetic states. Within this picture an excitation requires thermal assistance to make an Ising spin flip, but the time scale for the flip should be given by rate constant $\tau_0$, which is essentially temperature independent. The high temperature limit of this model is then an Ising paramagnetic system evolving with temperature independent rate constant, giving a plateau in the evolution of the relaxation time scale. The deviation from the plateaux above 10 K is due to thermal mixing of the crystal field levels [46, 44, 49, 50, 51, 52] and will not be discussed further here.

From figure (6) the tunneling rate, $\tau_0$ for Dy$_2$Ti$_2$O$_7$, appears to be around $10^{-3}$ s, while for Ho$_2$Ti$_2$O$_7$, it is on the scale of $\tau_0 \approx 10^{-8}$ s [48]. These extremely slow, yet exponentially varying microscopic time scales from one compound to another, seem consistent with a quantum tunneling process exponentially dependent on the energy levels of the rare-earth ions CEF. Unfortunately, the complexity of these CEF prevents quantitative estimates for the experimental inversion rate and the time scale $\tau_0$ remains a fitting parameter in our theory. The quantum tunneling on these microscopic time scales should lead to stochastic dynamics on all length scales and if a single micropscopic process dominates the dynamics, one should be able to model the relaxation using an Arrhenius law. For a more complex system, for example one with long range interactions, one should be able to use the Metropolis dynamics of a Monte Carlo simulation to model the evolution of the system. These techniques and the validity of the model, with temperature independent $\tau_0$ are discussed in detail in the next section.

## 4.1 Arrhenius law and deconfined quasi-particles

The first quantitative analysis of spin ice dynamics [42] relates to Ho$_2$Sn$_2$O$_7$ and Ho$_2$Ti$_2$O$_7$. Here the relaxation time, $\tau(T)$ was fitted with an Arrhenius law in the freezing regime ($T \approx 1$ K) with characteristic energy barrier, $E_f = 19.6$ K and 27.5 K respectively. These numbers already suggest that a simple Arrhenius law cannot give a complete description of the dynamics, as they do not correspond to any of the energy scales appearing in the DSI:

* $\sim 200$ K between the ground state doublet and the first excited crystal field levels;

* $\sim 4J_{eff} = 7.2$ K for Ho$_2$Ti$_2$O$_7$ and 4.4 K for Dy$_2$Ti$_2$O$_7$, due to single spin flips within the NNSI model;

* $\sim 3$ K for the limit of infinite separation between monopoles (see figure 5).
<figure><img src="image_8.png" /><figcaption>Figure 6: **Experimental relaxation time** τ for Dy₂Ti₂O₇ vs temperature T from Snyder & al. [45], courtesy of Peter Schiffer. The starting point of this work has been shed new light on this result. Above 12 K the spins are not exactly Ising anymore. The quasi-plateau between 2 and 12 K is divided in two parts: the “high temperature” of the spin ice region above 4 K where one needs to consider double defects 4 *in* or 4 *out* and a large monopole density, and the “low temperature” regime below 6 K where our monopole model should hold. Below ~2 K, the system becomes frozen.</figcaption></figure>

Quantitative studies of Dy₂Ti₂O₇ rule out such a simple thermally activated process [43, 45], as shown in figure 7. Here, the data of Snyder & al. are fitted in the low temperature regime with an energy barrier of E_f = 6.7 J-eff ~ 6.6 K, with J-eff = 1.11 K, the value estimated for Dy₂Ti₂O₇ [20]. The Arrhenius scaling is quantitatively good below 2 K, but fails completely to reproduce the quasi-plateau region. As we have argued above, if the tunneling picture is correct, the plateau region, through its simplicity, offers an anchor point for any fit and any procedure that does not fit here cannot contain the physics of the problem.

Rather than fitting the low temperature freezing region one should therefore start by fitting the data through the plateau region. By doing this we see the first evidence for deconfined quasi-particle excitations. Following the idea that the quasi-plateau is in fact the *tail* of an exponential, thermally activated process, we have plotted in figure 7 different Arrhenius laws τ₀ exp(Eₚ/T) where the energy barrier is varied. The time scale τ₀ is fixed by scaling to the experimental time at 4 K; a temperature at which we expect the density of double defects (4 *in* and 4 *out*) to be negligible, leading to a direct relationship between spin flips and monopole pair creation. From the figure one can see the main result of this section: the experimental data are poorly reproduced by an energy barrier Eₚ = 4 J-eff, the cost of a single spin flip within the NNSI model, but are quantitatively well fitted with Eₚ = 2 J-eff over the region between 2.5 K and 5 K (see the lower panel of figure 7). This means that the lowest excitation responsible for the dynamics of spin ice materials is *not* the energy cost of a single spin flip, but only half of that *i.e.* the energy cost of a single topological defect. Within this effective temperature range, the NNSI model seems sufficient to capture the physics, as suggested by polarized neutron scattering [53]. Once a pair of defects is created, they can propagate freely and separately in the system, leading to magnetic relaxation from independent (deconfined) objects.
<figure><img src="image_9.png" /><figcaption>Figure 7: **Deconfined monopoles:** The experimental data of Snyder & al. [45] (× and dashed line) are compared with different Arrhenius functions. **Top:** The quasi-plateau region is in quantitative agreement with a thermally activated process with energy barrier $E_p = 2 J_{\text{eff}}$ (red line) whereas the spin freezing is well reproduced with $E_p \approx 6 J_{\text{eff}}$ (green line), but no unique function can fit the whole temperature window. **Bottom:** The characteristic excitation is the creation of a unique defect ($E_p = 2 J_{\text{eff}}$, red line) rather than a single spin flip ($E_p = 4 J_{\text{eff}}$, grey line). The fit is improved at higher temperature if we include the energy scales due to double defects in a modified Arrhenius law (blue line).</figcaption></figure>
As shown in figure 7, some corrections to Arrhenious behaviour develop above 5 K, which we show in the next section, is due to the excitation of 4 *out* and 4 *in* double defects, by constructing a variant of the Arrhenius law, including all relevant energy scales for single spin flips within the context of the NNSI model.

## 4.2 Multi-energy Arrhenius law

We use an independent tetrahedron picture, assuming that the probability of finding a tetrahedron in a certain state is simply given by the Boltzmann weights of the two tetrahedra involved in any spin flip. We label the different configurations 2 *in* - 2 *out*, 3 *in* - 1 *out* (or 3 *out* - 1 *in*) and 4 *in* (or 4 *out*) 2, 3 and 4 respectively. In table IV.2 we show a list of all energy changes resulting from a single spin flip shared between two tetrahedra, from which we can calculate the mean energy barrier $\Delta E$.

<table><caption>Table 5: Energy scales resulting from a single spin flip: Let us explain this table using the second row as an example. It means that if we flip a spin that initially belongs to a tetrahedron 2 *in* - 2 *out* and a tetrahedron 3 *in* - 1 *out*, then there is a probability of 3/4 simply to inverse the position of the defect (3|2) at no energy cost, and a probability of 1/4 to create an additional pair of defects (3|4) that cost an energy $\delta E = +8 J_{\text{eff}}$. The probabilities are a statistical average, whether the spin under consideration is one of the three *in* spins of the 3 *in* - 1 *out* tetrahedron ($\delta E = 0$), or the fourth *out* spin ($\delta E = 8 J_{\text{eff}}$).</caption><thead><tr><th>Initial</th><th>Final</th><th>δE</th><th>Probability</th></tr></thead><tbody><tr><td>2|2</td><td>3|3</td><td>+ 4 J<sub>eff</sub></td><td>1</td></tr><tr><td rowspan="2">2|3</td><td>3|2</td><td>0</td><td>3/4</td></tr><tr><td>3|4</td><td>+ 8 J<sub>eff</sub></td><td>1/4</td></tr><tr><td>2|4</td><td>3|3</td><td>- 4 J<sub>eff</sub></td><td>1</td></tr><tr><td rowspan="3">3|3</td><td>2|2</td><td>- 4 J<sub>eff</sub></td><td>9/16</td></tr><tr><td>2|4</td><td>+ 4 J<sub>eff</sub></td><td>6/16</td></tr><tr><td>4|4</td><td>+12 J<sub>eff</sub></td><td>1/16</td></tr><tr><td rowspan="2">3|4</td><td>4|3</td><td>0</td><td>1/4</td></tr><tr><td>2|3</td><td>- 8 J<sub>eff</sub></td><td>3/4</td></tr><tr><td>4|4</td><td>3|3</td><td>-12 J<sub>eff</sub></td><td>1</td></tr></tbody></table>

There are two points worth noting: the configurations (2|3) and (3|2) are different, which gives a factor of 2 for their Boltzmann weight. More importantly, a single spin flip that actually creates a pair of defects will effectively only cost one half of $\delta E$, because of the fractionalisation of the excitations, but the process (3|3)→(2|4) really costs $\delta E = +4 J_{\text{eff}}$ as it simply moves a charge. We find $\Delta E$:

$$ \begin{aligned} \Delta E(T) &= \frac{2(6e^{2\beta J_{\text{eff}}})^2 + 4(2\frac{1}{4}6e^{2\beta J_{\text{eff}}})8 + 6(\frac{1}{16}8^2) + 4(\frac{6}{16}8^2)}{(6e^{2\beta J_{\text{eff}}})^2 + (2\frac{1}{4}6e^{2\beta J_{\text{eff}}})8 + (\frac{1}{16}8^2) + (\frac{6}{16}8^2)} J_{\text{eff}} \\ &= 3 \frac{6e^{4\beta J_{\text{eff}}} + 8e^{2\beta J_{\text{eff}}} + 10}{9e^{4\beta J_{\text{eff}}} + 6e^{2\beta J_{\text{eff}}} + 7} J_{\text{eff}} \xrightarrow{T \to 0} 2 J_{\text{eff}}, \end{aligned} \quad (14) $$

from which we recover the creation energy of a single defect in the low temperature limit. The modified Arrhenius law $\tau_0 \exp(\Delta E(T)/T)$ is plotted in figure 7 (blue curve). It is almost identical to the thermally activated process with $E_p = 2 J_{\text{eff}}$ (red line) below 4 K, but differs at higher temperature and follows
the experimental data more closely up to ≈ 8 K. Quantitative differences develop above this limit as one expects: above this temperature the thermally activated processes of the high temperature regime should make their presence felt. Using the function $τ_0 \exp(ΔE/T)$ clearly extends the range over which this phenomenological model can describe the data and we believe that it convincingly confirms the validity of the quantum tunneling model to describe the dynamics up to 8 – 10K. We believe that this phenomenological approach provides a very strong argument in favour of magnetic relaxation assisted by the thermally activated creation of deconfined quasi-particles in spin ice and hence in favour of **fractionalisation** of dipoles into free effective monopoles!

However as shown in 7, any Arrhenius function ultimately fails below 2-3 K, underestimating the time scale in the low temperature region. Along the quasi-plateau the density of monopoles is high enough for Debye screening to occur with the result that effective models based on the NNSI provide a convincing description. As the temperature is lowered, their density is reduced and a quantitative description can only be achieved by taking explicitly into account the long range Coulomb interactions. In this case, simple Arrhenius behaviour is excluded and we must revert to numerical simulation.

# 5 Spin freezing

We have simulated the dynamics of both the DSI model and a Coulomb gas of magnetic monopoles in the grand canonical ensemble occupying the sites of the diamond lattice, using a Metropolis Monte Carlo algorithm. The Coulomb particles are excitations out of the quasi-particle vacuum provided by the dumbbell model and their dynamics is constrained by the underlying network of Dirac strings. In both cases we treat the long range interactions, either dipolar or Coulombic using the Ewald method [31, 32].

## 5.1 Dipolar spin ice

Simulation of the DSI model is relatively straightforward. We implement both nearest neighbour and dipolar interactions, using the numerical values of $J$ and $D$ for $Dy_2Ti_2O_7$ given in table I.3. We extract the characteristic relaxation time $τ$, from the auto-correlation function

$$C(t) = \frac{1}{N} \sum_{i} S_{i}(0) \cdot S_{i}(t), \qquad (15)$$

where $N$ is the total number of spins ($N = 16 \times 8^3 = 8192$ here) and $S_i(t)$ is the value of the unit vector representing the spin at the Monte Carlo time $t$. One Monte Carlo time step, for a system of size $N$ is defined as a single sweep through $N$ moments chosen at random. For the initial conditions we take a random configuration, which we let evolve at high temperature $T = 10$ K until equilibrium, defined here as the time beyond which $C(t)$ decays below 0.01. This defines $t = 0$ and $C(t)$ is then computed and stored for a given temperature. The Monte Carlo time is re-set to zero when $C(t)$ decays beyond 0.001 and the temperature $T$ is lowered by $δT$. The process is repeated until we reach $T = 0.7$ K when numerical equilibration becomes difficult. The resulting auto-correlation is averaged over 200 samples in order to give a smoothly decaying function down to $C(t) \lesssim 0.03$. The outcome is plotted in figure 8 for $T = 1.5$ and 3 K.

$C(t)$ relaxes almost exponentially, but with a small difference between short and long time scales. Different time scales, $τ_{i=1,2}$, can be extracted by fitting the correlation function $C(t)$ to an exponential function, either over the interval $C = 1 - 0.8$ ($τ_1$) or over the interval $C = 0.3 - 0.03$ ($τ_2$)(see dashed lines in figure 8). The characteristic relaxation time $τ$, is taken to be the average of $τ_1$ and $τ_2$. The Monte Carlo time is again scaled to that measured by Snyder & al. by equating the scales at 4 K, $τ(T = 4$ K) = 2.99 ms, for which the density of double defects is negligible ($\lesssim 1\%$), giving a quantitative
<figure><img src="image_10.png" /><figcaption>Figure 8: **Time correlation function C(t) for DSI** for T = 1.5 (green) and 3 K (red) on a semi-log scale. The x-axis is the Monte Carlo time t. The dashed lines are exponential fits for short and long relaxation times.</figcaption></figure>

equivalence between our Metropolis dynamics and the relaxation dynamics of the experiment, with one Monte Carlo step = 2.5 ± 0.5 ms.

Preliminary results (not shown) were for the NNSI model using the $J_{eff}$ from Table 3. The data agree quantitatively with the analytic calculation, $\tau_0 \exp(\Delta E(T)/T)$, (see figure (7)) over the temperature range 0 – 10 K showing that our rather simple approach accurately describes the stochastic dynamics of the NNSI model. Any deviation from this behaviour is therefore explicitly due to the long range interactions of the DSI. Our results for $\tau$ vs $T$ for the DSI are compared with the data of Snyder & al. [45] in figure (9). We find remarkably good agreement between the numerical results and experiment over the temperature range between 1 K and 10 K. In particular there is a substantial improvement over the Arrhenius law of figure (7) below 3 K as we move into the low temperature regime with rapidly increasing time scale. Although when one looks more closely there are still some systematic differences between the data from the real and the numerical experiment, we believe that the general agreement shown in figure (9) represents a considerable success for both the DSI and stochastic quantum tunneling hypothesis, in which the tunneling rate is temperature independent. Further, the tunneling time scale is approximately equal to the Monte carlo time scale fixed at 4 K, giving $\tau_0 \sim 2.5$ milliseconds.

We start the investigation of exactly how the Coulomb interactions modify the spin ice dynamics by calculating the evolution of topological defects as a function of temperature in both the NNSI and DSI models. The evolution is shown in figure 10. In both cases the density of 3 *in* - 1 *out* and 3 *out* - 1 *in* defects shoots rapidly to zero below 1 – 2 K, which corresponds to the low temperature regime, while the double defect concentration goes to zero below the 5 – 6 K range. This confirms the picture that the plateau region is due to the proliferation of a high density gas, while the low temperature region with increasing time scale is driven by an exponential fall in the monopole concentration. It is therefore clear from this picture that monopole physics should explain the agreement between real and numerical experiments shown in figure (9). A second point to notice is that the monopole density falls to zero at a noticeable higher temperature for the DSI than for the NNSI model which must surely be related to the brutal slowing down of the dynamics compared with our analytical Arrhenius law calculation. To investigate this slowing down further we move to the dumbbell model and a description of the problem
<figure><img src="image_11.png" /><figcaption>Figure 9: Relaxation time scale from DSI The experimental data (×) are from Snyder & al. [45] (●). The system size is $N = 16000$ spins. The details of the simulations are given in the text. The bars do not represent the error bars on the value of $\tau_1$; rather they show the extreme values $\tau_1$ and $\tau_2$ at each temperature. In order to see the improvement due to dipolar interactions, the modified Arrhenius formula 14 is also plotted in blue.</figcaption></figure>

uniquely in terms of monopoles.

## 5.2 Monopoles and strings

Within the monopole picture spin flips out of and into the constrained manifold correspond to creation/annihilation of pairs, while all other flips correspond to monopole hopping between nearest neighbour sites of the diamond lattice. Hence in a Monte Carlo step one considers two nearest neighbour sites at time $t$. The possible outcomes at time $t+1$ are:

* If there are no quasi-particles on these two sites, then we consider the creation of a pair of opposite charges.

* If there is only one monopole, it can either move to the other site or stay where it is;

* If there are two opposite charges, they can annihilate;

* If there are two charges of the same sign, no movement is possible; as we have disallowed monopole defects carrying double charges.

### 5.2.1 Chemical potential $\mu$

Moving into the monopole picture means we change the Gibbs ensemble from which we select microstates. In the magnetic picture the independent variables are $T$, magnetic field $\vec{B}$ (chosen to be zero throughout
<figure><img src="image_12.png" /><figcaption>Figure 10: **Density of monopoles** for the dipolar spin ice model (•) and their expected Boltzmann weight for the nearest neighbour spin ice model (solid line) with the numerical values of Dy₂Ti₂O₇. In red (resp. blue) are plotted the concentration of 3 *in* - 1 *out* and 3 *out* - 1 *in* defects (resp. 4 *in* and 4 *out*), carrying a magnetic charge ±1 (resp. ±2). The system size is *N* = 16000 spins.</figcaption></figure>

this paper) and the number of spins *N*, which is also fixed to the volume of the sample. In the monopole picture, the number of quasi particles varies, so the independent variables are *T*, $\vec{B} = 0$, the volume *V* and $\mu$ the chemical potential for monopole pair creation. The relevant thermodynamic potential is therefore the Grand Potential $\Phi = U - ST - N\mu$, for which we have to define $\mu$. We are, of course free to choose any value of $\mu$ and arbitrary values will take us away from spin ice physics, but as the present work is motivated by the experimental results of Snyder et. al. we need to choose the value of $\mu$ that best corresponds to Dy₂Ti₂O₇. In a first series of simulations we have estimated $\mu = \mu_1$ numerically by calculating the difference between the Coulomb energy gained by creating a single pair of neighbouring magnetic monopoles of opposite sign $\Delta U_{\text{mono}} < 0$ in a vacuum and that required to produce a unique pair of topological defects out of the 2 *in* - 2 *out* manifold in the DSI, $\Delta U_{\text{def}} > 0$, giving a configurationally averaged estimate $\mu_1 \equiv \Delta U_{\text{mono}} - \Delta U_{\text{def}} = -8.92$ K. Note that while $\Delta U_{\text{mono}}$ is simply $V_{\text{min}}$ from equation (13), $\Delta U_{\text{def}}$ has to be computed by simulation, as projective equivalence is not exact for the DSI. In a second series of simulations, $\mu_2(T)$ was taken as the value required to reproduce the same defect concentration as in a simulation of dipolar spin ice at temperature *T*, as shown in figure 10. We found that $\mu_2 = \mu_1$ within numerical error at *T* = 0.7K and varied by ≈3% only, over the temperature range 1 - 4 K, showing that our procedures are consistent and underlining the fact that the chemical potential used is *not* a free parameter. Rather, it is taken from detailed comparison with the DSI that best describes Dy₂Ti₂O₇, with parameters from Table 3. Note that in references [27, 54] we have used an unconventional definition of $\mu$, with a sign change compared to the above value. Here, we make a definition compatible with conventional statistical mechanics such that, in a Monte Carlo move we compute the change in $\delta\tilde{U} = \delta U_{\text{mono}} - \mu\delta N$. As the value given, $\mu \approx -8.9$ K, is for a pair of particles, $\delta N = \pm 1$ or zero.
It is interesting to dwell further on $\mu$ and its composite parts at finite temperature, as their evolution with $T$ turns out to be at the origin of the rapid increase in the time scale at low temperature. We can make a third estimate for the chemical potential; $\mu(T) = \Delta U_{\text{mono}}(T) - \Delta U_{\text{def}}(T)$ using the same procedure as for $\mu_1$, but at finite temperature, where a pair of monopoles is created in a microstate containing a finite concentration of monopoles corresponding to finite temperature $T$. To do this, we simulated the DSI model, with $N = 8192$ spins in parallel with a 'slave' Coulomb gas, i.e. for each single spin flip the corresponding monopole configuration on the diamond lattice is updated accordingly, making it possible to compare the energy of the spin configuration with its mirror in the dumbbell model at each Monte Carlo step. As we are interested in the creation of monopoles, we have artificially suppressed double defects in the DSI for all temperatures. Hence after equilibration at a given temperature $T$ and for each accepted creation of a pair of defects, we compute the energy difference between the Monte Carlo move in the DSI model $\Delta U_{\text{def}}$ and the evolution in the slave Coulomb gas $\Delta U_{\text{mono}}$, which gives a value of $\mu$ for this specific move. By averaging these values over the Monte Carlo time and for different initial configurations, one finds the temperature evolution of these energies as well as their standard deviation, plotted in figure 11.

<figure><img src="image_13.png" /><figcaption>Figure 11: **Creation of a pair of quasi-particles:** Energy required in the dipolar spin ice model ($\Delta U_{\text{def}}$, red ●) or gained for a Coulomb gas ($\Delta U_{\text{mono}}$, blue ■) and the resulting chemical potential $\mu(T)$ (▲) as a function of temperature. The vertical bars are *not* error bars but the standard deviation of these quantities. $\mu(T)$ tends to the limit $\mu_1 = -8.92$ K (lower dashed line) at very low temperature. The other dashed line is the analytical prediction of the energy gained by creating a pair of monopoles in a vacuum $V_{\text{min}} = -3.07$ K (see equation (12)). The temperature scale below 3 K is enlarged for a better display.</figcaption></figure>

Much of the physics of spin ice at low temperature is included in this figure: the standard deviation
of both $\Delta U_{\text{def}}(T)$ (red) and $\Delta U_{\text{mono}}(T)$ (blue) are large above 1 K, up to ±40% of the mean values. However, these quantities are highly correlated and their difference, giving rise to $\mu(T)$ (black) is very precisely defined with a standard deviation of ±4% at most. Furthermore, the extracted $\mu(T)$ follows closely the values for $\mu_2$, which was estimated more phenomenologically. This agreement, despite the large fluctuations, illustrates the success of the monopole picture and of the dumbbell model in describing the physics of spin ice: the energy of separation of two monopoles is given by the magnetic Coulomb potential, including many body effects for a finite particle concentration and the energy of creation (and annihilation, not shown here) is consistent with a Coulomb gas in the grand canonical ensemble. On average, large energy cost $\Delta U_{\text{def}}$ (top of the red bars) correspond to small energy gain $\Delta U_{\text{mono}}$ (top of the blue bars) and vice-versa.

The data in figure (11) show two regimes of behaviour: from 10 K down to 2 K, corresponding to the plateau region for relaxation times and below 2 K where the relaxation time rapidly increases. In the higher temperature regime the mean values, $\langle\Delta U_{\text{def}}\rangle$ and $\langle\Delta U_{\text{mono}}\rangle$ are roughly constant and $\langle\Delta U_{\text{def}}\rangle = \langle\Delta U_{\text{mono}}\rangle - \langle\mu(T)\rangle$ is close to the energy cost of a single spin flip in the NNSI model $E = 4 J_{\text{eff}} = 4.44$ K. Thus, apart from the fact that there is a quasi-continuous range of energy values to create monopoles instead of a unique one and that the fluctuations in the energy scales are surprisingly large, the results are consistent with the NNSI model, which explains why the NNSI provides a reasonable description of spin ice down to 2 K. This equivalence between models with long and short range interactions is a direct consequence of Debye screening [36, 14]. That is, throughout this regime, $\langle\Delta U_{\text{mono}}\rangle$ (blue curve) is in excess of the energy of deconfinement of a pair of monopoles in a vacuum, $V_{\text{min}} = -3.07$ K (lower dashed line, see equation (12)) and from the fluctuations one can see that pair creations providing less energy are rare. Hence, a newly created pair is more stable than in a vacuum with the additional energy coming from the build up of correlations within the Coulomb gas, favouring a neighbouring charge cloud of opposite charge around each particle. As the temperature is reduced, screening becomes more marked so that low energy events become increasingly rare, leading a maximum in the mean Coulomb energy gain on pair creation, at around 2.5 K. However, below this temperature the monopole concentration drops rapidly (see figure 10) so that the Debye screening is starved and there is a marked evolution of the energies in figure (11). This appears to occur as the smallest values of $\Delta U_{\text{mono}}$ disappear brutally and $\langle\Delta U_{\text{mono}}\rangle$ increases rapidly towards $V_{\text{min}}$. This evolution is accompanied by a surprising increase in $\langle\Delta U_{\text{def}}\rangle$ (red line), ensuring that $\mu(T)$ remains temperature independent. The increase in the mean energy required to create a pair of topological defects as temperature decreases is clearly due to the absence of many body screening effects from the mobile quasi-particles. The brutal non-Arrhenius spin freezing as one leaves the plateau region at low temperature is therefore an avalanche effect: fewer monopoles hinders the Debye screening, which reduces the number of creations of quasi-particles, which leads to less monopoles, ... If this slowing down is sharp but continuous, it is because the system remains thermally activated. That is, as long as the material remains “hot” enough to allow the creation of a pair of monopoles out of the vacuum, i.e. a single spin flip of energy $\Delta U_{\text{def}} = \mu_1 + V_{\text{min}} \approx 5.8$ K, the system will remain dynamic through the creation and diffusion of monopoles. Our simulations hit this limit for $T \approx 0.6$ K and one might expect equilibration in experiments to become dramatically difficult below this temperature. This is exactly what Snyder & al. observed at 0.65 K for Dy$_2$Ti$_2$O$_7$ in FC-ZFC measurements [45], although magnetocaloric measurements [55] and recent muon spin resonance experiments[41] appear to measure the dynamical properties of spin ice materials at considerably lower temperature where the intrinsic monopole concentration, as predicted by the DSI or the dumbbell model would be considerably lower.

### 5.2.2 Dirac strings

The monopoles hop between nearest neighbour sites via the Metropolis Monte Carlo algorithm, giving diffusive dynamics, but with the local constraint discussed in the introduction: in the spin model a 3 in - 1 out topological defect can move at low energy cost by flipping one of the three in spins, the direction of the out spin being barred by an energy barrier of $8 J_{\text{eff}}$ that would result in a double defect and which
is set to infinity in the Coulomb gas. An isolated monopole can therefore hop to 3 out of 4 of its nearest neighbour sites only, dictated by an oriented network of constrained trajectories corresponding to the ensemble of classical Dirac strings [30] of connected dipoles [33]. The positively charged monopoles move in one sense along the network while the negative charges move in the opposite direction (see figure 4.b). The network is dynamically re-arranged through the evolution of the monopole configuration. In fact the characteristic time scale that we compare with experiment comes from the evolution of the network of Dirac strings rather than from the monopoles themselves. That is, the magnetic relaxation experiment of Snyder et. al. [45] ultimately measures spin dynamics, due monopole diffusion, rather than giving a direct measure of the monopoles themselves. We define the string network by an integer $\sigma = \pm 1$ giving the orientation of the Dirac string along each bond of the diamond lattice. The Dirac string network is an image of the underlying spin structure for which we can define the auto-correlation function

$$C(t) = \frac{1}{N} \sum_{i} \sigma_{i}(0)\sigma_{i}(t). \qquad (16)$$

We start the simulation with the Dirac string network in an ordered configuration and let it evolve at $T = 4$ K until equilibrium is reached. A high density of monopoles is generated, which move following the criteria defined above and time step is defined as a sweep through $N$ links of the monopole vacuum. We compute $C(t)$ averaged over 200 simulations following the same procedure as for the calculation of spin-spin autocorrelation function in the DSI. The Monte Carlo time was again scaled to the experimental time at $T = 4$ K and no other fitting parameters were used.

## 5.3 Comparison with experiment

In accordance with our analysis of figure 11, the characteristic time scale for relaxation of the Dirac string network does indeed reproduce the sharp spin freezing observed experimentally [27]. Results are shown in figure (12) in the interval below 4 K. Data from the Coulomb gas simulation is in very close agreement with the Monte Carlo data for the DSI, confirming that the low temperature dynamics of the DSI is driven by the dynamics of Coulombic quasi particles, so that a Coulomb gas in the grand canonical ensemble really is the correct description of the system at low temperature. Also shown in the figure is the modified Arrhenius formula 14 describing freezing in the NNSI. On this temperature scale the difference between data for the NNSI and that for the DSI or the Coulomb gas is considerable, confirming that the sharp, non-Arrhenius freezing at low temperature is due to the presence of long range interactions between quasi-particles. The most significant effect here is the gradual suppression of Debye screening below $\sim 2$ K, provoking the avalanche effect which accelerates the disappearance of the monopoles below this temperature. In addition we also expect the diffusive dynamics to be slowed by the formation of neutral pairs, bound by the Coulomb interaction. However, as the Coulomb interaction is not confining in three dimensions the monopole mobility should remain finite. Below 1.5 K the estimate of the time scale for the Coulomb gas lies slightly below that for the DSI. This small difference provides a measure of the approximation involved in passing from the spin system to the charged particles.

We show simulation results using both our definitions of the chemical potential. The data sets are very close to each other, although a small systematic difference appears around 1 K, with data for the fixed chemical potential, $\mu_1$ giving a slightly lower estimate for $\tau$. This small deviation seems consistent with the small evolution in $\mu(T)$ vs $T$ observed in figure (11).

Finally, when looked at with this level of precision, all numerical data falls below the experimental data at low temperature. This small but systematic difference may indicate the limits of applicability of both the DSI and our model of stochastic dynamics although more experiments need to be done to quantify this remark. More recent analysis of the static properties of Dysprosium Titanate suggest that the DSI can be improved by including further neighbour exchange[56], which will change fine details of the dynamics. The approximation of having a single, temperature independent microscopic time scale,
$\tau_0$, will also have its limitations which will also be manifested in the fine detail of the difference between simulation and experiment.

<figure><img src="image_14.png" /><figcaption>Figure 12: Spin freezing in Dy$_2$Ti$_2$O$_7$: The experimental data (×) are from Snyder & al. [45]. The Arrhenius law (blue line) represents the free diffusion of topological defects in the nearest neighbour model. The relaxation time scale of the Dirac string network for the Coulomb gas has been obtained for fixed chemical potential (▲) and with $\mu$ varying slowly to match the defect concentration in dipolar spin ice (■), whereas the dipolar spin ice relaxation time is given by •. As in figure 9 and as defined in the text, the bars do not represent the error bars on the value of $\tau$; rather they show the extreme values $\tau_1$ and $\tau_2$ at each temperature. The system size is 8 × 12$^3$ = 13824 sites on the diamond lattice for the Coulomb gas and 8192 spins on the pyrochlore lattice for DSI simulations.</figcaption></figure>

# 6 Conclusion

The recent discovery of monopole quasi-particles in spin ice models and materials provides a rare occasion to think outside the magnetic box, while studying a set of magnetic materials. The theoretical[30, 36], numerical[27, 54, 37] and experimental[17, 36, 41, 18] work so far presented in this field make this point very clear. Here, the collective quasi-particle excitations are localized in real space, rather than reciprocal space, despite the absence of quenched disorder. The quasi-particle vacuum has hidden internal degrees of freedom with $U(1)$ gauge symmetry and thanks to the unique role played by the dipole interactions on the pyrochlore lattice structure, the excitations have both the topological and energetic properties of magnetic charges. The correct description at low temperature is therefore that of a Coulomb (lattice) gas in the grand canonical ensemble, whose properties are almost identical to those of its electronic equivalent, an electrolyte of positive and negatively charged ions whose dynamics can be considered as purely stochastic.
There are some exceptions to this analogy, notably spin ice materials cannot support a dc monopole current in the steady state[27]. This is because of the dynamical constraint that positive and negative charges move in opposite directions over the underlying Dirac string network, changing the direction of the network as they pass. The system only posses a finite number of monopole trajectories in each direction before complete ordering of the network occurs. Or, in magnetic language, the sample can provide a finite number of spin flips in a given direction before becoming magnetically ordered. This fact seems difficult to reconcile with recent experiments[41], where the magnetic charge has been measured through application of Onsager’s Wien effect [40] to Dysprosium Titanate. In an electrolyte, in the low charge concentration limit, the electric current density, $\sigma(E)$ increases linearly with electric field $E$, with the result that the ratio $\sigma(E)/\sigma(0)$ is a universal number containing $E$, $T$, charge $q$ and fundamental constants only. By treating muon relaxation rate in an equivalent way to the electric current, Bramwell et. al. have been able to measure the magnetic charge $Q \approx 4.6\mu_B A^{-1}$ in good agreement with the theoretical prediction of Castelnovo et. al. [30]. How can this be, given that no permanent current can exist? One possible explanation is that, in weak fields and on short time scales, the constraints do not qualitatively affect the dynamics. One can define a drift length scale for field $X$, $\ell_X = k_B T/QX$, where $X$ is the magnitude of a magnetic or electric field. Above this scale one should observe the deterministic drift of the ions in the electrolyte and the effects of the constraints in the magnetic problem. In the field and temperature range of the Wien effect experiment ($T = 100 - 200$ mK, $B \approx 2$ mT), $\ell_B \sim 70$ Å. Given this large scale it seems reasonable that the magnetic and electric systems should show the same behaviour within the experimental set up. Much work is required here to understand these points in detail.

Spin ice dynamics at low temperature should be compared with the dynamics of a spin glass. In the latter, the glassy behaviour is the result of a hierarchy of time scales, as manifested, for example in the phenomenological random energy model [57]. Here the evolution steps occur with an exponentially varying set of delay times. For spin ice, at least in zero field things are very different. We have modeled the dynamics of the NNSI using a single time scale associated with the quantum tunneling of an Ising spin $\tau_0$. This scale dictates two processes; a hopping move of an existing deconfined monopole by one lattice spacing, or the creation (annihilation) of a nearest neighbour pair of particles, whose characteristic time is dressed by the Boltzmann weight, $\tau = \tau_o \exp(-\beta\Delta E)$. At low temperature this becomes the rate determining step, so that once created, the diffusing particles are actually extremely efficient at decorrelating the system. For example, from figure (7), the characteristic Arrhenius time scale at 500 mK is $\tau \sim 1$ sec, while the microscopic time scale, $\tau_0$ is estimated to be about 2.5 milliseconds. That is, the system decorrelates in 400 microscopic time steps thanks to the free diffusion of the deconfined particles. In the presence of Coulomb interactions, the decorrelation times are longer as the microscopic processes are dressed by a continuous range of energy scales coming from the Coulomb interaction, but the basic picture is provided by the simpler nearest neighbour model where free particles move on an essentially flat energy surface, making them efficient in decorrelating the system.

In this paper we have interpreted a magnetic relaxation experiment in terms of monopole dynamics. A more fundamental approach is to design experiments that explicitly and directly test for monopole dynamics. Up to press the only experiments which attempt this are those using the Wien effect [41, 39], in which monopole mobility can be probed by exploiting the chemical equilibrium existing between bound and free particles in a field. An even more direct approach would be an equivalent of the “Stanford” superconducting coil experiment [58, 30], whose goal was to measure the passage of free Dirac monopoles through a superconducting coil. Such an experiment could in principle detect the passage of a single magnetic quasi-particle but, given that the charges have no mass and therefore have diffusive, rather than Newtonian dynamics, at least in the absence of an external field and that the monopole charge in Dysprosium Titanate is approximately 8000 times smaller than a Dirac Monopole[30], this experiment presents serious challenges. A starting point could be a many body equivalent. In fact, such an experiment has been set up in another context by Hérrison and Ocio[59] to look at the violation of the fluctuation dissipation relation in spin glass systems. In this experiment the magnetic fluctuations are directly
measured and compared with magnetic response. The analysis presented in this paper suggests that the
resulting magnetic noise in spin ice would be due to the diffusive motion of monopoles. Unfortunately,
for the moment this set up is limited to the 4 K temperature range which remains in the high monopole
density regime. The goal for a monopole experiment would be first to descend below one Kelvin and then
to reduce the scale of the sample sufficiently so that the resulting noise was due to a mesoscopic, or even
microscopic number of quasi-particles.

In conclusion, the discovery of magnetic monopoles in spin ice has exposed a rich vein of physics,
in which a Coulomb gas, equivalent to a weak electrolyte appears at low temperature. This emergent
'magnetolyte'[41] shows many of the static and dynamic properties of the electric system within a very
idealized environment. Much experimental, theoretical and numerical work remains to be done in order
to fully exploit it.

Acknowledgements

It is a pleasure to thank S. T. Bramwell and R. Moessner for stimulating discussions and for related collaborations. We have, in addition had useful conversations with L. Balents, L. Berthier, B. Canals, C. Castelnovo, T. Fennell, S. Giblin, M. Gingras, V. Kaiser, K. Matsuhira and P. Schiffer. PCWH thanks the Max Planck Institute for Complex Systems in Dresden for financial support through the study program, "Unconventional Magnetism in High Fields" and apologises to the organizers for not exceeding 10 mT.

References

[1] M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey. Geometrical frustration in the ferromagnetic pyrochlore $h_{2}ti_{2}O_{7}$. *Physical Review Letters*, 79(13):2554–2557, SEP 29 1997.

[2] M. J. P. Gingras. *Introduction to Frustrated Magnetism*, volume 164 of *Solid State Science*, chapter Spin Ice. Springer (Ed. Lacroix, Mendels & Mila), 2011.

[3] S. Rosenkranz, A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. S. Shastry. Crystal-field interaction in the pyrochlore magnet $h_{2}ti_{2}O_{7}$. *Journal of Applied Physics*, 87(9):5914–5916, MAY 1 2000.

[4] L. Pauling. The structure and entropy of ice and of other crystals with some randomness of atomic arrangement. *Journal of the American Chemical Society*, 57:2680–2684, JUL-DEC 1935.

[5] W. F. Giauque and J. W. Stout. The entropy of water and the third law of thermodynamics - the heat capacity of ice from 15 to 273 degrees k. *Journal of the American Chemical Society*, 58:1144–1150, JUL-DEC 1936.

[6] A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. S. Shastry. Zero-point entropy in 'spin ice'. *Nature*, 399:333–335, MAY 27 1999.

[7] R. Moessner. Relief and generation of frustration in pyrochlore magnets by single-ion anisotropy. *Physical Review B*, 57(10):R5587–R5589, MAR 1 1998.

[8] S. T. Bramwell and M. J. Harris. Frustration in ising-type spin models on the pyrochlore lattice. *Journal of Physics-Condensed Matter*, 10(14):L215-L220, APR 13 1998.

[9] P. W. Anderson. Ordering and antiferromagnetism in ferrites. *Physical Review*, 102(4):1008–1013, 1956.
[10] D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi. Coulomb and liquid dimer models in three dimensions. *Physical Review Letters*, 91(16):167004, OCT 17 2003.

[11] M. Hermele, M. P. A. Fisher, and L. Balents. Pyrochlore photons: The u(1) spin liquid in a s=1/2 three-dimensional frustrated magnet. *Physical Review B*, 69(6):064404, FEB 2004.

[12] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi. Dipolar spin correlations in classical pyrochlore magnets. *Physical Review Letters*, 93(16):167204, OCT 15 2004.

[13] C. L. Henley. Power-law spin correlations in pyrochlore antiferromagnets. *Physical Review B*, 71(1):014424, JAN 2005.

[14] C. L. Henley. The "coulomb phase" in frustrated systems. *Annual Review of Condensed Matter Physics*, 1(1):179-210, 2010.

[15] R. W. Youngblood and J. D. Axe. Polarization fluctuations in ferroelectric models. *Physical Review B*, 23(1):232-238, 1981.

[16] T. Fennell, S. T. Bramwell, D. F. McMorrow, P. Manuel, and A. R. Wildes. Pinch points and kasteleyn transitions in kagome ice. *Nature Physics*, 3(8):566-572, AUG 2007.

[17] 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. Magnetic coulomb phase in the spin ice $ho_2ti_{2o_7}$. *Science*, 326(5951):415-417, October 2009.

[18] Hiroaki Kadowaki, Naohiro Doi, Yuji Aoki, Yoshikazu Tabata, Taku J. Sato, Jeffrey W. Lynn, Kazuyuki Matsuhira, and Zenji Hiroi. Observation of magnetic monopoles in spin ice. *Journal of the Physical Society of Japan*, 78(10), October 2009.

[19] R. Siddharthan, B. S. Shastry, A. P. Ramirez, A. Hayashi, R. J. Cava, and S. Rosenkranz. Ising pyrochlore magnets: Low-temperature properties, "ice rules", and beyond. *Physical Review Letters*, 83(9):1854-1857, AUG 30 1999.

[20] B. C. den Hertog and M. J. P. Gingras. Dipolar interactions and origin of spin ice in ising pyrochlore magnets. *Physical Review Letters*, 84(15):3430-3433, APR 10 2000.

[21] S. T. Bramwell, M. J. Harris, B. C. den Hertog, M. J. P. Gingras, J. S. Gardner, D. F. McMorrow, A. R. Wildes, A. L. Cornelius, J. D. M. Champion, R. G. Melko, and T. Fennell. Spin correlations in $ho_2ti_{2o_7}$: A dipolar spin ice system. *Physical Review Letters*, 87(4):047205, JUL 23 2001.

[22] H. Kadowaki, Y. Ishii, K. Matsuhira, and Y. Hinatsu. Neutron scattering study of dipolar spin ice $ho_2sn_{2o_7}$: Frustrated pyrochlore magnet. *Physical Review B*, 65(14):144421, APR 1 2002.

[23] R. G. Melko and M. J. P. Gingras. Monte carlo studies of the dipolar spin ice model. *Journal of Physics-Condensed Matter*, 16(43):R1277-R1319, NOV 3 2004.

[24] S. Yoshida, K. Nemoto, and K. Wada. Ordered phase of dipolar spin ice under [110] magnetic field. *Journal of the Physical Society of Japan*, 73(7):1619-1622, JUL 2004.

[25] R. Higashinaka and Y. Maeno. Field-induced transition on a triangular plane in the spin-ice compound $dy_2ti_{2o_7}$. *Physical Review Letters*, 95(23):237208, DEC 2 2005.

[26] J. P. C. Ruff, R. G. Melko, and M. J. P. Gingras. Finite-temperature transitions in dipolar spin ice in a large magnetic field. *Physical Review Letters*, 95(9):097202, AUG 26 2005.
[27] L. D. C. Jaubert and P. C. W. Holdsworth. Signature of magnetic monopole and dirac string dynamics in spin ice. *Nature Physics*, 5(4):258-261, APR 2009.

[28] M. J. P. Gingras and B. C. den Hertog. Origin of spin-ice behavior in ising pyrochlore magnets with long-range dipole interactions: an insight from mean-field theory. *Canadian Journal of Physics*, 79(11-12):1339-1351, NOV-DEC 2001.

[29] S. V. Isakov, R. Moessner, and S. L. Sondhi. Why spin ice obeys the ice rules. *Physical Review Letters*, 95(21):217201, NOV 18 2005.

[30] C. Castelnovo, R. Moessner, and S. L. Sondhi. Magnetic monopoles in spin ice. *Nature*, 451:42-45, JAN 3 2008.

[31] S. W. de Leeuw, J. W. Perram, and E. R. Smith. Simulation of electrostatic systems in periodic boundary-conditions .1. lattice sums and dielectric-constants. *Proceedings of the Royal Society of London A*, 373(1752):27-56, 1980.

[32] D. Frenkel and B. Smit. *Understanding Molecular Simulations*. Academic Press Elsevier Science, 2002.

[33] J. D. Jackson. *Classical Electrodynamics*. Wiley, New-York, 1999.

[34] P. A. M. Dirac. Quantised singularities in the electromagnetic field. *Proceedings of the Royal Society of London A*, 133:60-72, 1931.

[35] I. A. Ryzhkin. Magnetic relaxation in rare-earth oxide pyrochlores. *Journal of Experimental and Theoretical Physics*, 101(3):481-486, 2005.

[36] D. J. P. Morris, D. A. Tennant, S. A. Grigera, B. Klemke, C. Castelnovo, R. Moessner, C. Czternasty, M. Meissner, K. C. Rule, J. U. Hoffmann, K. Kiefer, S. Gerischer, D. Slobinsky, and R. S. Perry. Dirac strings and magnetic monopoles in the spin ice $dy_2ti_2o_7$. *Science*, 326(5951):411-414, October 2009.

[37] C. Castelnovo, R. Moessner, and S. L. Sondhi. Thermal quenches in spin ice. *Physical Review Letters*, 104(10), March 2010.

[38] P Mayer, S Leonard, L Berthier, JP Garrahan, and P Sollich. Activated aging dynamics and negative fluctuation-dissipation ratios. *Physical Review Letters*, 96(3), January 2006.

[39] S. R. Giblin, S. T. Bramwell, P. C. W. Holdsworth, D. Prabhakaran, and I. Terry. Creation and measurement of long-lived magnetic monopole currents in spin ice. *unpublished*, 2010.

[40] L. Onsager. Deviations from ohm's law in weak electrolytes. *Journal of Chemical Physics*, 2:599, 1934.

[41] S. T. Bramwell, S. R. Giblin, S. Calder, R. Aldus, D. Prabhakaran, and T. Fennell. Measurement of the charge and current of magnetic monopoles in spin ice. *Nature*, 461(7266):956-959, October 2009.

[42] K. Matsuhira, Y. Hinatsu, K. Tenya, and T. Sakakibara. Low temperature magnetic properties of frustrated pyrochlore ferromagnets $ho_2sn_2o_7$ and $ho_2ti_2o_7$. *Journal of Physics-Condensed Matter*, 12(40):L649-L656, OCT 9 2000.

[43] K. Matsuhira, Y. Hinatsu, and T. Sakakibara. Novel dynamical magnetic properties in the spin ice compound $dy_2ti_2o_7$. *Journal of Physics-Condensed Matter*, 13(31):L737-L746, AUG 6 2001.
[44] G. Ehlers, A. L. Cornelius, M. Orendac, M. Kajnakova, T. Fennell, S. T. Bramwell, and J. S. Gardner. Dynamical crossover in 'hot' spin ice. *Journal of Physics-Condensed Matter*, 15(2):L9-L15, JAN 22 2003.

[45] J. Snyder, B. G. Ueland, J. S. Slusky, H. Karunadasa, R. J. Cava, and P. Schiffer. Low-temperature spin freezing in the $dy_2ti_{2o}_7$ spin ice. *Physical Review B*, 69(6):064414, FEB 2004.

[46] J. Snyder, J. S. Slusky, R. J. Cava, and P. Schiffer. How 'spin ice' freezes. *Nature*, 413:48-51, SEP 6 2001.

[47] G. Ehlers, A. L. Cornelius, T. Fennell, M. Koza, S. T. Bramwell, and J. S. Gardner. Evidence for two distinct spin relaxation mechanisms in 'hot' spin ice $ho_2ti_{2o}_7$. *Journal of Physics-Condensed Matter*, 16(11):S635-S642, MAR 24 2004.

[48] J. P. Clancy, J. P. C. Ruff, S. R. Dunsiger, Y. Zhao, H. A. Dabkowska, J. S. Gardner, Y. Qiu, J. R. D. Copley, T. Jenkins, and B. D. Gaulin. Revisiting static and dynamic spin-ice correlations in $ho_2ti_{2o}_7$ with neutron scattering. *Physical Review B*, 79(1):014408, JAN 2009.

[49] John P. Sutter, Satoshi Tsutsui, Ryuji Higashinaka, Yoshiteru Maeno, Olaf Leupold, and Alfred Q. R. Baron. Relaxation in the spin ice $dy_2ti_{2o}_7$ studied using nuclear forward scattering. *Physical Review B*, 75(14):140402, APR 2007.

[50] J. Shi, Z. Tang, B. P. Zhu, P. Huang, D. Yin, C. Z. Li, Y. Wang, and H. Wen. Dynamical magnetic properties of the spin ice crystal $dy_2ti_{2o}_7$. *Journal of Magnetism and Magnetic Materials*, 310(2):1322-1324, MAR 2007.

[51] J. Lago, S. J. Blundell, and C. Baines. µsr investigation of spin dynamics in the spin-ice material $dy_2ti_{2o}_7$. *Journal of Physics-Condensed Matter*, 19(32):326210, AUG 15 2007.

[52] K. Kitagawa, R. Higashinaka, K. Ishida, Y. Maeno, and M. Takigawa. High-temperature spin relaxation process in $dy_2ti_{2o}_7$ probed by ti-47-nqr. *Physical Review B*, 77(21), June 2008.

[53] L. J. Chang, Y. Su, H. Schneider, R. Mittal, M. R. Lees, G. Balakrishan, and Th Brückel. Magnetic correlations in pyrochlore spin ice as probed by polarized neutron scattering. *Journal of Physics: Conference Series*, 211:012013, 2010.

[54] L. D. C. Jaubert. *Topological Constraints and Defects in Spin Ice*. PhD thesis, Ecole Normale Supérieure de Lyon, SEP 2009.

[55] M. Orendac, J. Hanko, E. Cizmar, A. Orendacova, M. Shirai, and S. T. Bramwell. Magnetocaloric study of spin relaxation in dipolar spin ice $dy_2ti_{2o}_7$. *Physical Review B*, 75(10):104425, MAR 2007.

[56] T. Yavors'kii, T. Fennell, M. J. P. Gingras, and S. T. Bramwell. $Dy_2ti_{2o}_7$ spin ice: A test case for emergent clusters in a frustrated magnet. *Physical Review Letters*, 101(3):037204, JUL 18 2008.

[57] B. Derrida. Random energy model - an exactly solvable model of disordered systems. *Physical Review B*, 24(5):2613-2626, 1981.

[58] B. Cabrera. First results from a superconductive detector for moving magnetic monopoles. *Physical Review Letters*, 48:1378-1381, 1982.

[59] D. Herisson and M. Ocio. Fluctuation-dissipation ratio of a spin glass in the aging regime. *Physical Review Letters*, 88(25), June 2002.