# Gauge-Free Duality in Pure Square Spin Ice: Topological Currents and
Monopoles

Cristiano Nisoli
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Theoretical Division, Los Alamos National
Laboratory
Los Alamos, NM, 87545, USA

###### Abstract

We consider a pure square spin ice, that is a square ice where only nearest
neighbors are coupled. A gauge-free duality between the perpendicular and
collinear structure leads to a natural description in terms of topological
currents and charges as the relevant degrees of freedom. That, in turn, can be
expressed via a continuous field theory where the discrete spins are subsumed
into entropic interactions among charges and currents. This approach produces
structure factors, correlations, and susceptibilities for spins, monopoles,
and currents. It also generalizes to non-zero temperature the height formalism
of the disordered ground state. The framework can be applied to a zoology of
recent experimental results, especially realizations on quantum annealers and
can be expanded to include longer range interactions.

###### Contents

  1. I Introduction
  2. II Square Ice and Its Gauge-Free Duality
    1. II.1 Gauge-Free duality in 2D
    2. II.2 Square Spin Ice
    3. II.3 Heuristic Entropy and Ice-Like Correlations
  3. III Energy and States
  4. IV Field Theory: Exact Results
  5. V Field Theory: Approximations
    1. V.1 High Temperature
      1. V.1.1 Effective Energies and Entropic Interactions
      2. V.1.2 Correlations, Susceptibilities, and Height Functions
    2. V.2 Low Temperature
  6. VI Cases
    1. VI.1 Pure Degenerate Spin Ice ($\kappa=0$)
      1. VI.1.1 Charge Correlations
      2. VI.1.2 Spin Correlations
    2. VI.2 Weak Ferromagnetic Case ($\kappa\simeq 0^{+}$)
    3. VI.3 Weak Antiferromagnetic Case ($\kappa\simeq 0^{-}$)
  7. VII Naïve Kinetics
    1. VII.1 General Considerations
    2. VII.2 Relaxation Dynamics
      1. VII.2.1 Dynamic Susceptibilities, Conductivities, Dispersion
      2. VII.2.2 Real Space Picture
  8. VIII Conclusions
  9. A Conceptualizing the Dressing at Low $T$

##  I Introduction

Degenerate artificial square ice is perhaps the simplest two-dimensional
system in which to study disorder constrained by the ice ruleLieb (1967), and
its violations as monopole excitationsRyzhkin (2005); Castelnovo, Moessner,
and Sondhi (2008). While spin ice pyrochloresBramwell and Harris (2020);
Ramirez  _et al._ (1999); den Hertog and Gingras (2000); Bramwell and Gingras
(2001) had opened new vistas in the study of geometric frustration and
constrained disorder, more recently artificial spin ices—systems of
interacting, magnetic nanoislandsTanaka  _et al._ (2006); Wang  _et al._
(2006); Nisoli, Moessner, and Schiffer (2013); Heyderman and Stamps (2013);
Skjærvø  _et al._ (2019); Ortiz-Ambriz  _et al._ (2019)—have provided
controllable platforms that can be characterized at the constituent level.

Although the field has developed to include new forms of frustration and
geometries, allowing for the realization of magnets often not found in
natureMorrison, Nelson, and Nisoli (2013); Stamps (2014); Nisoli, Kapaklis,
and Schiffer (2017); Li  _et al._ (2018); Saccone  _et al._ (2019); Makarova
_et al._ (2021); Gliga  _et al._ (2017); Stopfel  _et al._ (2018), and
revealing novel phenomena absent in its crystal analogueSkjærvø  _et al._
(2019); Schiffer and Nisoli (2021), recent fabrication and characterization
advancesArnalds  _et al._ (2012); Farhan  _et al._ (2013); Kapaklis  _et al._
(2014) have brought back to the fore the austere simplicity of celebrated
early models, such as kagome and square ices.

The square geometry of spin ice was among Tanaka  _et al._ (2006) the first to
be realized artificially Wang  _et al._ (2006), when it was shown that non-ice
rule Bernal and Fowler (1933) vertices are suppressed after AC demagnetization
Wang  _et al._ (2006); Nisoli  _et al._ (2007, 2010). However, because in
nanoislands moments impinging perpendicularly in the vertex interact more
strongly than moments impinging collinearly, the degeneracy of the ice
manifold is lifted and antiferromagnetic (AFM) vertices are favored, leading
to a phase transition in the Ising classWu (1969); Anghinolfi  _et al._ (2015)
toward an ordered antiferromagnetic ground stateMorgan  _et al._ (2010); Porro
_et al._ (2013); Zhang  _et al._ (2013); Sendetskyi  _et al._ (2019).

Soon after, Möller and Moessner Möller and Moessner (2006) proposed to offset
the height of half of the nanoislands to regain the ice rule degeneracy
Giauque and Ashley (1933); Pauling (1935). The idea was recently
realizedPerrin, Canals, and Rougemaille (2016); Farhan  _et al._ (2019).
Meanwhile, ice rule degeneracy in square ice has been also demonstrated in
rectangular lattices Nascimento  _et al._ (2012); Ribeiro  _et al._ (2017), or
via cleverly placed interaction modifiers placed in the verticesÖstman  _et
al._ (2018), or nano-holes in a connected spin ice of nanowires Schanilec  _et
al._ (2019), or by rotation of the moments Macêdo  _et al._ (2018). Recently,
a pure square ice, i.e. a square ice with only nearest neighbor interaction,
was realized in a quantum annealer King  _et al._ (2021), and used to
demonstrate purely entropic monopole interactions. Then, various different
regimes can be achieved by tweaking specs and couplings, leading to
antiferromagnetic states, line states, or ice manifolds, which lead to
different spectral characterizations rou .

The scope of this work is to provide an unifying and expandable framework that
can cover many different square ice systems around the ice rule degeneracy
point Möller and Moessner (2006); King  _et al._ (2021), by considering
topological charges and currents as the relevant degrees of freedom and
subsuming the spin structure into effective, entropic interactions among them.
This approach flows naturally from a gauge-free duality of the system which is
absent in three dimensions (3D).

We limit ourselves to pure square ice, that is a sixteen vertex model where
interactions are limited to spins within the same vertex Lieb and Wu (1972);
Baxter (1982); Wu (1969). We consider no long-range interactions, and thus no
3D-Coulomb (i.e. $1/r$) interaction among monopoles or currents (we have
considered such case elsewhere Nisoli (2020a)). However, because of the
emergent nature of these objects, we show that they interact via 2D-Coulomb
(i.e. $\sim\ln r$) entropic interactions.

Within this model, we compute free energies, entropic interactions, structure
factors, susceptibilities, correlations, screening, relaxation dynamics. We
discuss strengths and limits of this approach. We show new results, but also
re-derive in a coherent framework results that were previously appreciated in
similar systems through a variety of methods. These had included
phenomenological approaches via coarse grained field, height models, or
analogies with chemical physics approaches Henley (2010); Isakov  _et al._
(2004); Garanin and Canals (1999); Henley (2005); Youngblood and Axe (1981);
Huse  _et al._ (2003); Henley (2011, 1997); Bramwell (2012); Twengström,
Henelius, and Bramwell (2020); Jaccard (1964); Ryzhkin (2005). We also
particularized some results that what we had already found on generic graphs
Nisoli (2020b).

##  II Square Ice and Its Gauge-Free Duality

There is in 2D a gauge-free duality absent in 3D pyrochlore ice. It is related
to the rather gravid mathematical fact that in 2D a Helmholtz decomposition
has no gauge freedom. That in turns follows from the fact that orthogonal
directions are uniquely defined in 2D. It is amply used in 2D continuum
theories, from fluid dynamics to the $XY$ model Kosterlitz and Thouless
(1973), and is behind the entire edifice of complex analysis.

###  II.1 Gauge-Free duality in 2D

In 2D we can always write a continuum vector field $\vec{S}$ in terms of
longitudinal and transverse potentials $h_{||}$, $h_{\\!\perp}$,

$\displaystyle\vec{S}=\vec{S}_{||}+\vec{S}_{\\!\perp}=\vec{\nabla}h_{||}-\hat{e}_{3}\wedge\vec{\nabla}h_{\\!\perp},$
(1)

where $\hat{e}_{3}=\hat{e}_{1}\wedge\hat{e}_{2}$, $\hat{e}_{1},\hat{e}_{2}$ is
an orthonormal basis of the plane, and we call
$\vec{S}_{||},\vec{S}_{\\!\perp}$ the longitudinal and perpendicular
components of the field. Unlike in the 3D case, where the perpendicular part
of $\vec{S}$ is the curl of a vector potential, there is no gauge freedom in
Eq (1).

If we define the charge and current distributions of the field as

$\displaystyle q[\vec{S}]\vcentcolon$
$\displaystyle=-\vec{\nabla}\cdot\vec{S}$ $\displaystyle
i[\vec{S}]\vcentcolon$
$\displaystyle=\hat{e}_{3}\cdot\vec{\nabla}\wedge\vec{S},$ (2)

then

$\displaystyle q$ $\displaystyle=-\Delta h_{||}$ $\displaystyle i$
$\displaystyle=-\Delta h_{\\!\perp}.$ (3)

If for a vector $\vec{w}$ we call

${}^{\perp}\\!\vec{w}\vcentcolon=\hat{e}_{3}\wedge\vec{w}$ (4)

the perpendicular of $\vec{w}$, we have then

$\displaystyle q[\vec{S}]$ $\displaystyle=i[^{\perp}\\!\vec{S}]$
$\displaystyle i[\vec{S}]$ $\displaystyle=-q[^{\perp}\\!\vec{S}],$ (5)

which expresses the duality between charges and currents, or longitudinal and
perpendicular components of the field, under a $\pi/2$ rotation.

This duality has a discretized analogue in square spin ice.

###  II.2 Square Spin Ice

Square spin ice (Fig. 1) is a set of classical, binary spins $\vec{S}_{e}$
aligned on the $N_{e}$ edges $e$ of a square lattice of $N_{v}=N_{e}/2$
vertices labeled by $v$. Spins form four vertex topologies often classified
Wang  _et al._ (2006) as t-I, …, t-IV, where t-I and t-II obey the ice rule
Bernal and Fowler (1933); Pauling (1935); Baxter (1982); Möller and Moessner
(2006) (i.e. have two spin pointing in, two pointing out).

In a degenerate square ice, ice rule vertices are degenerate and energetically
favored, and lead to a ground state Lieb (1967) of constrained disorder,
called the Ice Manifold. The latter is a Coulomb phase Henley (2010); Isakov
_et al._ (2004); Garanin and Canals (1999); Henley (2005), i.e a topological
state labeled, in lieu of an order parameter, by a height field Henley (2011);
Lamberty, Papanikolaou, and Henley (2013).

Figure 1: Top: the sixteen vertices of square ice can be divided into four
topologies, listed with degeneracy in parenthesis and topological charge.
Below, an ice rule obeying configuration of $\vec{S}$ (black) and its height
function $h_{\\!\perp}$ built from ${}^{\perp}\\!\vec{S}$ (gray). Also, the
coupling constants among spins.

Consider $Q_{v}[S]$, the topological charge of the vertex $v$, defined as the
number of spins pointing in the vertex minus those pointing out. Then, an ice
rule vertex $v$ has $Q_{v}=0$.

We can similarly define the topological current $I_{p}[S]$ of a minimal square
plaquette $p$, as the number of spins pointing clockwise around the edge of
the plaquette minus those pointing counterclockwise 111We call it current
because a magnetization $\vec{M}$ generates an electrical current density
$\vec{j}=\vec{\nabla}\wedge\vec{M}$.

For a spin configuration $\vec{S}$, consider its perpendicular configuration
${}^{\perp\\!}\vec{S}$ (Fig. 1), for which old plaquettes are now vertices and
old vertices are now plaquettes. We have then

$\displaystyle Q_{v}[\vec{S}]$ $\displaystyle=I_{v}[^{\perp}\\!\vec{S}],$
$\displaystyle I_{p}[\vec{S}]$ $\displaystyle=-Q_{p}[^{\perp}\\!\vec{S}]$ (6)

as in Eqs. (5).

For each configuration of spins that obeys the ice rule configuration, a
unique (up to a constant) height function $h_{\\!\perp}$ can be defined on the
plaquettes such that

${{}^{\perp}\\!\vec{S}_{e}}\cdot\hat{pp^{\prime}}=h_{\\!\perp
p^{\prime}}-h_{\\!\perp p},$ (7)

is true ($\hat{pp^{\prime}}$ is the unit vector pointing from plaquette $p$ to
$p^{\prime}$ separated by the edge $e$ as in Fig. 1). This follows from the
fact that ${}^{\perp\\!}S$ is “irrotational”: the line-sum of spins
${}^{\perp\\!}\vec{S}$ (grey in Fig 1) along a closed loop is zero.

Similarly, if a spin configuration has zero topological current on each
plaquette ($I_{p}[S]=0$ $\forall p$), a height function $h_{||v}$ can be
defined on the vertices by

$\vec{S}_{e}\cdot\hat{vv^{\prime}}=h_{||v^{\prime}}-h_{||v}$ (8)

where $e$ is the edge connecting the vertices $vv^{\prime}$.

Thus, in the ice manifold there are no charges, currents are disordered, and
${}^{\perp}\\!S$ is the discrete gradient of $h_{\\!\perp}$, which “labels”
the disorder of currents in absence of charges. Conversely, in the current-
free manifold, $S$ is the discrete of $h_{||}$, which labels the disorder of
charges in absence of currents.

The Eqs. (7, 8) are thus the discrete analogous to Eqs. (1) in the continuum.
A significant difference is that they are well defined only in the charge free
or current free state. We will show in the section V.1.2 how to generalize
them to any spin ensemble, inclusive of monopole and current excitations, and
thus at non-zero temperature.

###  II.3 Heuristic Entropy and Ice-Like Correlations

The height formalism can usefully, if heuristically, describe the pure ice
manifold. Height models are said to be in a “rough” (degenerate) of “flat”
(ordered) phase, a jargon derived from the theory of the roughening transition
historically associated with these models by various exact mappings van
Beijeren (1977); Chui and Weeks (1976).

The ice manifold of the square ice, i.e. the six vertex model, is known to be
equivalent to a dimer cover model Zinn-Justin (2009); Baxter (1982) and thus
in a rough phase Henley (2010, 2011, 1997). A widespreadvan Beijeren (1977);
Chui and Weeks (1976); Henley (1997, 2011) though by no means rigorously
justified ansatz ascribes to a configuration in the ice manifold an entropy
Raghavan, Henley, and Arouh (1997) that is quadratic in the height function.
In our language we can write

${\cal
S}[h_{\\!\perp}]=-\frac{1}{2\chi_{0}}\int\left(\vec{\nabla}h_{\\!\perp}\right)^{2}d^{2}x,$
(9)

where $h_{\\!\perp}$ is homogenized into a continuum field and $\chi_{0}$ is a
positive uniform susceptibility (see later). (Clearly, the same is true for
the current-free manifold by replacing $h_{\\!\perp}$ with $h_{\\!||}$.)

Eq. (9) can be understood in terms of the zero temperature partition function

$Z[\vec{H}]=\int\left[dh_{\\!\perp}\right]\exp\left({{\cal
S}[h_{\\!\perp}]+\hat{e}_{3}\cdot\int\vec{H}\wedge\nabla
h_{\\!\perp}d^{2}x}\right)$ (10)

where $\vec{H}$ is an external field.

Note that Eqs. (9,10) are not obviously unproblematic in a 2D (and thus gauge-
free) theory. In 3D we would be safe, as gauge invariance of the transverse
part of the field forbids the proliferation of relevant operators at the fixed
point (which, incidentally, is why an Higgs boson is needed in the standard
model). Eq. (9) merely happens to work in reproducing correlations that can
also in part be computed exactlySutherland (1968); Baxter (1982) (see also
refHenley (1997) and references therein for a discussion).

A series of interesting deductions come from Eqs. (9,10). For heigh function
correlations one immediately finds in reciprocal space

$\langle|\tilde{h}_{\\!\perp}(k)|^{2}\rangle=\frac{\chi_{0}}{k^{2}}.$ (11)

From that and the continuum limit of Eq (7) (or
$\vec{S}=-\hat{e}_{3}\wedge\vec{\nabla}h_{\\!\perp}$) one obtains the spin
correlator as 222We use dyadics: if $\vec{v}$ is a vector, the components of
$\vec{v}\vec{v}$ are $[\vec{v}\vec{v}]_{ij}\vcentcolon=v_{i}v_{j}$;
$\mathbbm{1}$ is the unitary matrix.

$\chi_{0}^{-1}\langle\vec{\tilde{S}}^{*}(k)\vec{\tilde{S}}(k)\rangle=\frac{{}^{\perp}\\!\vec{k}^{\perp}\\!\vec{k}}{k^{2}}=\mathbbm{1}-\frac{\vec{k}\vec{k}}{k^{2}}$
(12)

which is purely transversal.

Note that from Eq. (11), in real space, correlations of the height function
are logarithmic, or

$\langle
h_{\\!\perp}(x)h_{\\!\perp}(y)\rangle=-\frac{\chi_{0}}{2\pi}\ln(|x-y|).$ (13)

Spin correlations in real space can be obtained from Eq. (12) or more easily
as partial transversal derivatives (or
${}^{\perp}\vec{\nabla}=\hat{e}_{3}\wedge\vec{\nabla}$) of Eq. (11), obtaining

$\chi_{0}^{-1}\langle\vec{S}(x)\vec{S}(0)\rangle=\delta(x)\mathbbm{1}+\mathbbm{P}(x),$
(14)

where $\mathbbm{P}(x)$ is the kernel of the dipole-dipole interaction in 2D,
or

$\mathbbm{P}(x)=\frac{1}{2\pi}\left(\frac{\mathbbm{1}}{x^{2}}-2\frac{\vec{x}\vec{x}}{x^{4}}\right).$
(15)

Analogously, in the case of pyrochlore spin ice the spins correlations are the
kernel of the 3D dipolar interactionHenley (2010).

The spin correlations are therefore algebraic, making the ice-manifold a
critical phase of infinite correlation length. On the other hand the
correlation length for currents is zero: from Eqs. (11, 3) we have

$\langle|\tilde{i}(k)|^{2}\rangle={\chi_{0}}{k^{2}},$ (16)

which implies the infinitely localized screening of any pinned current.

In terms of currents, from Eq. (9) the entropy for the ice manifold can be
rewritten as

${\cal S}[i]=\frac{1}{2\chi_{0}}\int
d^{2}x~{}\\!i(x)\ln(x-y)i(y)~{}\\!d^{2}y,$ (17)

i.e. as a pairwise 2D-Coulomb interaction among the currents.

Our phenomenological picture is thus the following: in the ice manifold
charges are absent, disorder can be labeled by currents and their 2D-Coulomb
mutual interaction determines the entropy.

From Eqs. (9, 1), the entropy of a configuration in the ice manifold, can be
written in terms of $\vec{S}$ as

${\cal S}[\vec{S}]=-\frac{1}{2\chi_{0}}\int\vec{S}^{2}d^{2}x,$ (18)

with the constraint $\vec{\nabla}\cdot\vec{S}=0$.

Equation (18), unlike our previous formulas, is also valid in 3D. There, it
has been appreciated as Jaccard entropy in water ice Jaccard (1964); Nagle
(1966, 1978, 1979); Ryzhkin and Whitworth (1997), and later in pyrochlore spin
ice Huse  _et al._ (2003); Isakov  _et al._ (2004); Henley (2005, 2010) as
necessary to produce purely transverse correlations in the ice manifold. We
see therefore that in square ice the Jaccard entropy describes in fact a
2D-Coulomb interaction among disordered currents. The uniform susceptibility
$\chi_{0}$ is thus related to the so-called $\Phi$ constant Ryzhkin and
Whitworth (1997).

We note that in the previous deductions we have taken some cavalier liberties
with the boundary terms. With collaborators, we have already shown how fixing
the boundaries can induce a net charge in the bulk via a geometrical
expression of the Gauss’s Law King  _et al._ (2021). Taking the orthogonal of
the spin distribution, that implies a net current in the bulk if a current is
present on the boundaries. In a future work we will consider these interesting
topological effects for currents at the boundaries.

In the following we will show how to deduce a field theory for charges and
currents in which the heuristic Eqs. (1, 3) make sense, the intuitive height
function formalism [Eqs. (9-18)] finds a solid ground, and it is generalized
for $T>0$.

##  III Energy and States

The following Hamiltonian

$\displaystyle{\cal
H}[Q,I]=\frac{\epsilon}{2}\sum_{v}Q_{v}^{2}+\frac{\kappa}{2}\sum_{p}I_{p}^{2},$
(19)

reflects the current-charge duality by placing a cost or advantage on
topological currents and monopoles ($\epsilon$ and $\kappa$ are energies). In
terms of an Ising model, it is equivalent to a $J_{1},J_{2},J_{3}$ model where
(Fig. 1): $J_{1}=\epsilon-\kappa$ , $J_{2}=-\kappa$, $J_{3}=\epsilon$. By the
duality, the symmetry by orthogonalization corresponds to
$\epsilon\leftrightarrow\kappa$. The Hamiltonian describes various cases,
often close to the experimental reality.

If $\kappa=0$ and $\epsilon>0$, the ground state is the ice manifold of Fig. 1
(black arrows). Equivalently, by gauge-free duality, if $\epsilon=0$ and
$\kappa>0$, the ground state is an extensively degenerate ice manifold for
${}^{\perp}\\!\vec{S}$, or the grey arrows in Fig. 1.

If $\kappa=0$ and $\epsilon<0$, the ground state is the charge full state,
i.e. the ordered, antiferromagnetic tessellation of t-IV vertices. If
$\epsilon=0$ and $\kappa<0$ the ground state is the current full state, i.e.
the ordered, antiferromagnetic tessellation of t-I vertices, which is the
orthogonal of the charge full state.

For $\kappa>0$, $\epsilon>0$, both charges and currents are suppressed.
Because $J_{1}<J_{3}$, ferromagnetic t-II vertices are promoted over t-I.
Because $J_{2}<0$, t-II vertices want to align and the ground state is the
four-fold ferromagnetic state, made of t-II vertices ferromagnetically
aligned. (More loosely: as the ground state is current-free and charge-free,
we have $\Delta h_{||}=\Delta h_{\\!\perp}=0$ which implies a uniform
$\vec{S}$.) This case has not been investigated experimentally, though it can
certainly be realized in a quantum annealer King  _et al._ (2021). It might
approximate, however, experimental situations where t-II vertices can be
favoredPerrin, Canals, and Rougemaille (2016); Östman  _et al._ (2018);
Perrin, Canals, and Rougemaille (2019), leading to a line state that is
disordered but of sub-extensive entropy.

For $\kappa<0,\epsilon>0$, the ice rule is enforced at low $T$, currents are
promoted, $J_{1}>J_{3}$, and the ground state is an ordered antiferromagnetic
(AFM) tessellation of t-I vertices. Large $|\kappa|$ describes early square
ice realizationsWang  _et al._ (2006). Small $|\kappa|$ might approximate spin
ices that are designed so that ice rule vertices are degenerate Möller and
Moessner (2006); Perrin, Canals, and Rougemaille (2016) but where the dipolar
interaction still favors closed magnetic fluxes and thus promotes topological
currents and an ordered ground state.

When $\kappa=\epsilon$, we have $J_{1}=0$ and the set of vertical and
horizontal arrows become two decoupled systems. When $\kappa=\epsilon>0$ each
subsystem is ferromagnetic and the ground state is the four-fold fully
polarized state. When $\kappa=\epsilon<0$, each subsystem is antiferromagnetic
and the ground state is a four-fold antiferromagnetic state, which includes
the two orientations of the charge-full state and the two orientations of the
current-full state (one is the orthogonal of the other, as expected on the
symmetry line $\kappa=\epsilon$).

We will not study here the full monopole-currents model of Eq (19), which
leads to a rich phase diagram in the $\beta\epsilon\times\beta\kappa$ plane
(where $\beta=1/T$ and $T$ is temperature measured in units of energy) to be
compared with other models Levis  _et al._ (2013); Wu (1969). We will consider
only $\epsilon>0$ and $|\kappa|/\epsilon$ small and investigate how ice
manifold features are retained by small perturbations around the spin ice
point $\kappa=0$.

##  IV Field Theory: Exact Results

Because charges and currents represent an emergent description of pure square
ice, we deduce a field theory for which they are the relevant degrees of
freedom.

We generalize our previous approach for general graphs Nisoli (2020b, a) to
include currents. The partition function from Eq. (19) reads

$\displaystyle Z$ $\displaystyle=\sum_{S}\exp\left(-\beta{\cal
H}\right)\times$
$\displaystyle\exp\\!\left(\beta\sum_{e}\vec{S}_{e}\cdot\vec{H}_{e}+\beta\sum_{v}V_{q,v}Q_{v}+\beta\sum_{p}V_{i,p}I_{p}\right),$
(20)

and it is the generator of correlations

$\displaystyle\langle\vec{S}_{e_{1}}\dots\vec{S}_{e_{n}}\rangle=\partial_{\beta\vec{H}_{e_{1}}\dots\beta\vec{H}_{e_{n}}}\\!\\!\ln
Z$ $\displaystyle\langle Q_{v_{1}}\dots Q_{v_{n}}\rangle=\partial_{\beta
V_{q,v_{1}}}\dots\partial_{\beta V_{q,v_{n}}}\\!\\!\ln Z$
$\displaystyle\langle I_{p_{1}}\dots I_{p_{n}}\rangle=\partial_{\beta
V_{i,p_{1}}}\dots\partial_{\beta V_{i,p_{n}}}\\!\\!\ln Z.$ (21)

The fields $\vec{H},V_{q},V_{i}$ are measured in units of energy.

To obtain a continuum field theory we insert in the sum of (20) the tautology

$\displaystyle 1$ $\displaystyle=(2\pi)^{-2N_{v}}\prod_{v}\int
dq_{v}d\phi_{v}\exp\left[\mathbbm{i}\phi_{v}\\!\\!\left(q_{v}-Q_{v}\right)\right]$
$\displaystyle\times\prod_{p}\int
di_{p}d\psi_{p}\exp\left[\mathbbm{i}\psi_{p}\left(i_{p}-I_{p}\right)\right]$
(22)

and then sum over the spins, obtaining

$Z=\int\left[dqdi\right]\tilde{\Omega}[q,i]\mathbbm{e}^{-\beta{\cal
H}[q,i]+\sum_{v}q_{v}V_{q,v}+\sum_{p}i_{p}V_{i},p}$ (23)

where $\left[dqdi\right]=(2\pi)^{-N_{v}}\prod_{v}dq_{v}\prod_{p}di_{p}$.
$\tilde{\Omega}[q,i]$ is a generalized density of states for $q_{v},i_{p}$,
given by

$\displaystyle\tilde{\Omega}[q,i]=\int[d\phi
d\psi]\Omega[\phi,\psi]\mathbbm{e}^{\sum_{v}q_{v}\mathbbm{i}\phi_{v}+\sum_{p}i_{p}\mathbbm{i}\psi_{p}},$
(24)

and is therefore the functional Fourier transform of

$\Omega[\phi,\psi]=2^{N_{e}}{\prod_{{\hat{vv^{\prime}}}}}\cosh\left(-\mathbbm{i}\nabla_{vv^{\prime}}\phi-\mathbbm{i}\nabla_{pp^{\prime}}\psi+\beta
H_{vv^{\prime}}\right).$ (25)

(The product runs on all the edges $e={\hat{vv^{\prime}}}$ once, and
$\nabla_{vv^{\prime}}\phi\vcentcolon=\phi_{v^{\prime}}-\phi_{v}$,
$\nabla_{pp^{\prime}}\psi\vcentcolon=\psi_{p^{\prime}}-\psi_{p}$,
$H_{vv^{\prime}}\vcentcolon=\vec{H}_{e}\cdot\hat{vv^{\prime}}$, while
$\hat{pp^{\prime}}=-\hat{e}_{z}\wedge\hat{vv^{\prime}}$.)

Note that by construction $\langle Q_{v_{1}}\dots Q_{v_{n}}\rangle=\langle
q_{v_{1}}\dots q_{v_{n}}\rangle$, $\langle I_{p_{1}}\dots
I_{p_{n}}\rangle=\langle i_{p_{1}}\dots i_{p_{n}}\rangle$. Note also that the
use of $\vec{H}$, $V_{q},V_{i}$ is superabundant: $V_{q}\to
V_{q}+V^{\prime}_{q},V_{i}\to V_{i}+V^{\prime}_{i}$ is equivalent to
$\vec{H}\to\vec{H}+\vec{\nabla}V^{\prime}_{q}-\hat{e}_{3}\wedge\nabla
V^{\prime}_{i}$, but we keep it because it is useful.

We have gone from binary variables to a theory of continuous emergent
topological charges and currents constrained by an entropy

$S[q,i]=-T\ln\tilde{\Omega}[q,i]$ (26)

which conveys the effect of the underlying spin ensemble. Equivalently, in the
language of field theory, charges $q_{v}$ and currents $i_{p}$ interact
entropically via the fields

$\displaystyle V_{q}^{e}=\mathbbm{i}T\phi$ $\displaystyle
V_{i}^{e}=\mathbbm{i}T\psi,$ (27)

of generalized free energy

${\cal F}[\phi,\psi]=-T\ln\Omega[\phi,\psi].$ (28)

Note that ${\cal F}[\phi,\psi]$ can have imaginary values. As a consequence,
though the variables $\phi,\psi$ are themselves real, their expectation values
$\langle\phi\rangle$, $\langle\psi\rangle$ are imaginary, and thus the
entropic fields $\langle V_{q}^{e}\rangle,\langle V_{i}^{e}\rangle$ are real.
Indeed, by integrating over $q$ in Eq (23) and applying the second and third
equation in (21), one finds

$\displaystyle\langle V_{q}^{e}\rangle+V_{q}=\epsilon\langle q\rangle,$
$\displaystyle\langle V_{i}^{e}\rangle+V_{i}=\kappa\langle i\rangle,$ (29)

when $\kappa>0$. Similar Gaussian gymnastics also prove that

$\langle S_{vv^{\prime}}\rangle=\langle\tanh\left(\beta
H_{vv^{\prime}}-\mathbbm{i}\nabla_{vv^{\prime}}\phi-\mathbbm{i}\nabla_{pp^{\prime}}\psi\right)\rangle.$
(30)

We thus see that while $\mathbbm{i}\phi,\mathbbm{i}\psi$ correlate charges and
currents, the vector fields $-\mathbbm{i}T\nabla_{vv^{\prime}}\phi$,
$-\mathbbm{i}T\nabla_{pp^{\prime}}\psi$ correlate spins that would otherwise
be trivially paramagnetic.

In the spirit of field theory, we could therefore say the following:
$V_{q}^{e}$ is an entropic potential acting on charges and $V^{e}_{i}$ on
currents, and correspondingly there is an entropic field $\vec{B}^{e}$ acting
on the spins which in the long wavelength approximation is given by

$\displaystyle\vec{B}^{e}=\vec{B}^{e}_{||}+\vec{B}^{e}_{\\!\perp}=-\vec{\nabla}V_{q}^{e}-\hat{e}_{3}\wedge\vec{\nabla}V_{i}^{e}.$
(31)

##  V Field Theory: Approximations

###  V.1 High Temperature

When $\epsilon>0$ and $\kappa>0$ (see Section VI.3 for $\kappa\leq 0$)
integrating over charges and currents in Eq. (23) returns Gaussians in the
fields, of variance $\langle\phi^{2}\rangle=\epsilon/T$,
$\langle\psi^{2}\rangle=\kappa/T$: as temperature increases, the entropic
fields mediating spin correlations become smaller, as one would expect. In
fact, in the infinite temperature limit (but with $\beta H$ constant), the
Gaussian distributions in the field become Dirac deltas, and Eq. (23) becomes
the standard “paramagnetic” partition function

$Z=2^{N_{s}}\prod_{\hat{vv^{\prime}}}\cosh\\!\\!\left(\beta
H_{\hat{vv^{\prime}}}\right).$ (32)

Note that for $H=0$, the equation above returns the correct entropy per spin
at infinite temperature, $s=\ln 2$.

In the high $T$ limit, we can linearize Eq. (30) and write, in the long
wavelength limit,

$\displaystyle\langle\vec{S}\rangle=\beta\left(\vec{H}+\langle\vec{B}^{e}\rangle\right).$
(33)

Then, from $q=-\vec{\nabla}\cdot\vec{S}$,
$j=\hat{e}_{3}\cdot\vec{\nabla}\wedge\vec{S}$, and Eqs. (29) we find screened
Poisson equations for $q$, $i$

$\displaystyle-\xi_{||}^{2}\Delta q+q$ $\displaystyle=\beta\nu
q_{\text{ext}},$ $\displaystyle-\xi_{\perp}^{2}\Delta i+i$
$\displaystyle=\beta\nu i_{\text{ext}}$ (34)

(here ${\nu}q_{\text{ext}}=\vec{\nabla}\cdot\vec{H}$, where $q_{\text{ext}}$
is the external charge, $\nu i_{\text{ext}}=\vec{\nabla}\wedge\vec{H}$, where
$i_{\text{ext}}$ is the external current, $\nu$ is an energy) and therefore
“correlation lengths” 333When monopoles interact via a $1/r$ law, $\xi_{||}$
is no longer a correlation length, as the monopole-monopole interaction
destroys the screening at least in principle, as we have shown elsewhere
Nisoli (2020a). as

$\displaystyle\xi^{2}_{||}$ $\displaystyle={\epsilon/T}$
$\displaystyle\xi^{2}_{\\!\perp}$ $\displaystyle=\kappa/T.$ (35)

This expression for $\xi_{||}$ was already appreciated in other works, for
different geometries, via different meansGaranin and Canals (1999); Henley
(2005) and we had generalized it on a graph Nisoli (2020b). Note that for
$\kappa<0$, $\xi_{\\!\perp}$ would be imaginary, pointing to a periodicity
typical of the AFM ensemble, as we will discuss in section VI.3. Instead when
$\kappa=0$ the integral over $i$ returns a functional delta function on
$\psi$, and $\psi$ disappears from the equations, which for all purpose is
equivalent to taking $\xi_{\perp}=0$.

####  V.1.1 Effective Energies and Entropic Interactions

We now explore these heuristic deductions more precisely. Because we are
interested in the monopole liquid below the ice manifold threshold $T\simeq
2\epsilon$ but above other possible low $T$ transitions, a high $T$
approximation is a good starting point. Since it corresponds to small entropic
fields, we expand $\ln\Omega[\phi,\psi]$ at quadratic order and Fourier
transform via
$g_{x}=\int_{BZ}\tilde{g}(\vec{k})\mathbbm{e}^{-\mathbbm{i}\vec{k}\cdot
x}{d^{2}k}/{(2\pi)^{2}},$ where BZ is the Brillouin Zone, $g$ is a generic
field, and $x=v,p,l$ represents vertices, edges or plaquettes.

We obtain the approximated partition function at zero loop

$\displaystyle Z_{2}=\int[dqdi][d\phi
d\psi]\exp\\!\left(-\int_{\text{BZ}}\beta{\cal
H}_{2}(\vec{k})\frac{d^{2}k}{(2\pi)^{2}}\right)$

where

$\displaystyle{\cal H}_{2}[\tilde{q},\tilde{i},\tilde{\phi},\tilde{\psi}]$
$\displaystyle=\frac{\epsilon}{2}\left|\tilde{q}\right|^{2}+\frac{\kappa}{2}|\tilde{i}|^{2}+\frac{T\chi_{0}}{2}\gamma^{2}\\!\\!\left(|\tilde{\phi}|^{2}+|\tilde{\psi}|^{2}\right)$
$\displaystyle-\mathbbm{i}T\\!\left(\tilde{q}^{*}\tilde{\phi}+\tilde{i}^{*}\tilde{\psi}\right)$
$\displaystyle+\mathbbm{i}\chi_{0}\\!\left({}^{\perp}\\!\vec{\gamma}\tilde{\psi}^{*}-\vec{\gamma}\tilde{\phi}^{*}\right)\cdot\vec{\tilde{H}}-\chi_{0}\frac{\beta}{2}\left|\vec{\tilde{H}}\right|^{2}$
(37)

is an effective Hamiltonian for the new variables and the vector
$\vec{\gamma}$ has components

$\gamma_{\alpha}\vcentcolon=2\sin(k_{\alpha}/2)$ (38)

for $\alpha=x,y$, and contains informations on the square symmetry of the
lattice [in the long wavelength limit: $\vec{\gamma}\simeq\vec{k}+O(k^{3})$].

The first line of Eq. (37) contains the free energies for the uncoupled
charge, currents and entropic fields.

The second line contains the coupling between charges, currents and their
entropic fields: importantly, no $\phi\psi$ cross term survives at quadratic
order in the high $T$ approximation and charge and currents are independent at
second order, leading to so-called magnetic fragmentation Brooks-Bartlett  _et
al._ (2014); Petit  _et al._ (2016); Canals  _et al._ (2016), that is the
decoupling between charge-full and current-full ensembles.

The third line contains the coupling with the external field $\vec{H}$ (we
neglect here $V$). Because $-i\vec{\gamma}\cdot\vec{\tilde{w}}$,
$-i\vec{\gamma}\wedge\vec{\tilde{w}}$ are the generalized divergence and curl
respectively on the lattice in moment space for a generic field $\vec{w}$, the
entropic fields for currents (resp. charges) couple to the curl (resp.
divergence) of the external field.

The quantity $\chi_{0}$, already encountered in the previous section on the
opposite limit ($T=0$) is the uniform susceptibility. In Eq. (37) it is
$\chi_{0}=1$, but we include it nonetheless in the equation because it can be
$\chi_{0}\neq 1$ at low temperatures (see below).

Integrating $Z_{2}$ over $\tilde{\phi},\tilde{\psi}$ when $\tilde{H}=0$
returns the effective free energies for $q$ and $i$ in absence of external
field, which are decoupled at second order:

$\displaystyle{\cal H}_{2}[q,i]={\cal H}_{2}[q]+{\cal H}_{2}[i],$ (39)

where the two terms can be written as

$\displaystyle\beta{\cal H}_{2}[q]$
$\displaystyle=\frac{1}{2\chi_{0}}\\!\\!\left({\xi^{2}_{||}}+\gamma^{-2}\right)|\tilde{q}|^{2}$
$\displaystyle\beta{\cal H}_{2}[i]$
$\displaystyle=\frac{1}{2\chi_{0}}\\!\\!\left({\xi^{2}_{\perp}}+\gamma^{-2}\right)|{\tilde{i}}|^{2}.$
(40)

In field theory language, the first terms in Equations (40) are the “masses”
of the monopole or of the current. The second impliy that the underlying spin
ensemble mediates a pairwise entropic interaction among charges (or currents)
which is the 2D-Coulomb potential. In real space at large distances (where
$\gamma^{2}\simeq k^{2}$) the entropic interactions for charges and currents
are

$\displaystyle
V_{q}^{e}(v-v^{\prime})\simeq-\frac{T}{\chi_{0}}\frac{q_{v}q_{v^{\prime}}}{2\pi}\ln|v-v^{\prime}|$
$\displaystyle
V_{i}^{e}(p-p^{\prime})\simeq-\frac{T}{\chi_{0}}\frac{i_{p}i_{p^{\prime}}}{2\pi}\ln|p-p^{\prime}|.$
(41)

The origin of entropic interactions is in the emergent nature of charges and
currents, which exist in a spin vacuum. An assignation of charges and/or
currents to the system changes the number of ways in which the underlying spin
ensemble can be compatibly arranged. Thus the entropy associated to a fixed
distribution of charges and currents depends on their mutual position and
distance. While this is obvious, it is not obvious that changes in entropy can
be written via pairwise terms, as shown above, leading to an effective
pairwise interaction. We see that subsuming the effect of the underlying spins
leads to a “2D electrodynamics” formalism for charges and currents.

Integrating instead over $\tilde{q},\tilde{i}$ in Eq. (13)

####  V.1.2 Correlations, Susceptibilities, and Height Functions

From Eq. (39) and equipartition we obtain the correlations for charges and
currents

$\displaystyle\langle|\tilde{q}(\vec{k})|^{2}\rangle=\gamma(\vec{k})^{2}\tilde{\chi}_{||}(k)$
$\displaystyle\langle|\tilde{i}(\vec{k})|^{2}\rangle=\gamma(\vec{k})^{2}\tilde{\chi}_{\\!\perp}(k),$
(42)

where $\tilde{\chi}_{||}(k),\tilde{\chi}_{\\!\perp}(k)$ are defined as

$\displaystyle\tilde{\chi}_{||,\perp}(k)=\frac{\chi_{0}}{1+\xi_{||,\perp}^{2}\gamma(\vec{k})^{2}},$

Correctly, from Eq. (42) we have$\langle|\tilde{q}(\vec{0})|^{2}\rangle=0$,
since
$\langle|\tilde{q}(\vec{0})|^{2}\rangle=\langle\left(\sum_{v}q_{v}\right)^{2}\rangle=0$
(and same for currents). Figure 2 shows a heat map for
$\langle|\tilde{q}(\vec{k})|^{2}\rangle$.

Furthermore, if $\langle q^{2}\rangle$ is the average charge per vertex, or
$\langle q^{2}\rangle=\frac{1}{N_{v}}\sum_{v}\langle q_{v}^{2}\rangle$, then

$\langle
q^{2}\rangle=\int_{\text{BZ}}\langle|\tilde{q}(\vec{k})|^{2}\rangle\frac{d^{2}k}{(2\pi)^{2}}.$
(44)

From it, and Eq. (42), we obtain $\langle q^{2}\rangle\to 4$ for
$T\uparrow\infty$, which is correct. Indeed $\langle q^{2}\rangle=4$ is the
value deducible from a multiplicity argument ($2^{2}/2+4^{2}/8=4$). Because of
the duality, the same is true for currents.

Figure 2: Plots of the charge form factor
$\langle|\tilde{q}(\vec{k})|^{2}\rangle\xi_{||}^{2}/\chi_{0}$ for $\xi_{||}=1$
(top) and $\xi_{||}=5$ (bottom). White lines denote the Brillouin zone.
Identical plots hold for
$\langle|\tilde{i}(\vec{k})|^{2}\rangle\xi_{\\!\perp}^{2}/\chi_{0}$

.

Note that $\tilde{\chi}_{||}(k),\tilde{\chi}_{\\!\perp}(k)$ are the the
longitudinal and perpendicular static susceptibilities (multiplied by $T$,
thus accounting already for the Curie-Weiss law) for
${\tilde{S}}_{||}=\vec{\tilde{S}}\cdot\hat{\gamma}$ (the longitudinal, charge-
full, current-fee part of the magnetization) and
${\tilde{S}}_{\\!\perp}=\vec{\tilde{S}}\cdot^{\perp\\!\\!}\hat{\gamma}$ (the
perpendicular, charge-free, current-full magnetization), respectively
444$\hat{\gamma}$ is the unit vector of $\vec{\gamma}$. Indeed, by integrating
the Gaussian integral in Eq. (LABEL:Z2-n2), and remembering that
$\hat{\gamma}_{\alpha}\hat{\gamma}_{\alpha^{\prime}}+{{}^{\perp}\\!\hat{\gamma}_{\alpha}}\\!\\!{{}^{\perp}\\!\hat{\gamma}_{\alpha^{\prime}}}=\delta_{\alpha\alpha^{\prime}}$,
one obtains the total free energy as a function of the external field at
lowest order as

$\displaystyle\beta{\cal
F}_{2}[H]=-\frac{1}{2}(\beta\vec{\tilde{H}}^{*})\cdot\left({\hat{\gamma}\hat{\gamma}}{\tilde{\chi}_{||}}+{{{}^{\perp}\\!\hat{\gamma}}{{}^{\perp}\\!\hat{\gamma}}}{\tilde{\chi}_{\\!\perp}}\right)\cdot(\beta\vec{\tilde{H}}).$
(45)

From it, the spin correlations are

$\langle\tilde{S}^{*}_{\alpha}(\vec{k})\tilde{S}_{\alpha^{\prime}}(\vec{k})\rangle=\hat{\gamma}_{\alpha}\hat{\gamma}_{\alpha^{\prime}}\tilde{\chi}_{||}+{{}^{\perp}\\!\hat{\gamma}_{\alpha}}\\!\\!{{}^{\perp}\\!\hat{\gamma}_{\alpha^{\prime}}}\tilde{\chi}_{\\!\perp}.$
(46)

In the limit $T\to\infty$, Eq. (46) correctly returns
$\langle\tilde{S}_{\alpha}(\vec{k})\tilde{S}_{\alpha^{\prime}}(\vec{k})\rangle\to\delta_{\alpha\alpha^{\prime}}$,
or uncorrelated spins.

Equation (46) also implies that the magnetic susceptibilities in reciprocal
space can be obtained from experimentally measured spin spin correlations as

$\displaystyle\tilde{\chi}_{\\!\perp}(\vec{k})$
$\displaystyle=^{\perp}\\!\\!\hat{\gamma}\cdot\langle\vec{\tilde{S^{*}}}(\vec{k})\vec{\tilde{S}}(\vec{k})\rangle_{\text{exp}}\cdot^{\perp}\\!\\!\hat{\gamma}$
$\displaystyle\tilde{\chi}_{||}(\vec{k})$
$\displaystyle=\hat{\gamma}\cdot\langle\vec{\tilde{S^{*}}}(\vec{k})\vec{\tilde{S}}(\vec{k})\rangle_{\text{exp}}\cdot\hat{\gamma}$
(47)

and similarly, from (42)

$\displaystyle\langle|\tilde{q}(\vec{k})|^{2}\rangle=\vec{\gamma}\cdot\langle\vec{\tilde{S^{*}}}(\vec{k})\vec{\tilde{S}}(\vec{k})\rangle\cdot\vec{\gamma},$
$\displaystyle\langle|\tilde{i}(\vec{k})|^{2}\rangle=^{\perp}\\!\\!\vec{\gamma}\cdot\langle\vec{\tilde{S^{*}}}(\vec{k})\vec{\tilde{S}}(\vec{k})\rangle\cdot^{\perp}\\!\\!\vec{\gamma}.$
(48)

Finally, from the spin correlations we can obtain the magnetic structure
factor

$\Sigma_{m}(\vec{k})=^{\perp}\\!\\!\hat{k}\cdot\langle\vec{\tilde{S^{*}}}(\vec{k})\vec{\tilde{S}}(\vec{k})\rangle\cdot^{\perp}\\!\\!\hat{k},$
(49)

which at small $k$ corresponds to $\tilde{\chi}\\!_{\perp}$.

Taking the long wavelength approximation, Eq. (46) implies the following
effective free energy for the coarse grained spins

$\beta{\cal
F}_{2}[S]=\int_{\text{BZ}}\\!\left(\frac{1}{2\tilde{\chi}_{||}}{\tilde{S}}_{||}^{*}\cdot{\tilde{S}}_{||}+\frac{1}{2\tilde{\chi}_{\\!\perp}}{\tilde{S}}_{\\!\perp}^{*}\cdot{\tilde{S}}_{\\!\perp}\right)\frac{d^{2}k}{(2\pi)^{2}}.$
(50)

which in the long wavelength limit becomes

$\beta{\cal
F}_{2}[S]=\frac{1}{2\chi_{0}}\int\\!\left[{\xi^{2}_{||}}\left(\vec{\nabla}\cdot\vec{S}\right)^{2}+{\xi^{2}_{\\!\perp}}\left(\vec{\nabla}\wedge\vec{S}\right)^{2}+S^{2}\right]d^{2}x.$
(51)

The first two terms in the previous equation are energetic. The third is the
Jaccard entropy Jaccard (1964); Ryzhkin (2005) mentioned above. In the ice
manifold it is the only surviving term, since $\xi_{\\!\perp}=0$ and
$\nabla\cdot\vec{S}=0$, thus returning Eq. (18), previously found
heuristically.

Equation (51), when $\xi_{\\!\perp}=0$ reduces to the functional found
elegantly via methods of chemical physics for pyrochlore ice by Bramwell
Bramwell (2012); Twengström, Henelius, and Bramwell (2020). Indeed, when
expressed in terms of spins $\vec{S}$ rather than currents and charges $i,q$,
the two formalisms are expected to coincide.

Expressing the coarse grained $\vec{S}$ via height functions as in Eq. (1) we
obtain the free energy for the heigh functions in the long wavelength limit,
at quadratic oder:

$\displaystyle\beta{\cal F}_{2}[h_{||},h_{\\!\perp}]$
$\displaystyle=\frac{1}{2}\int\left[\left(\vec{\nabla}h_{||}\right)^{2}+\xi_{||}^{2}\left(\Delta
h_{||}\right)^{2}\right]d^{2}x$
$\displaystyle+\frac{1}{2}\int\left[\left(\vec{\nabla}h_{\\!\perp}\right)^{2}+\xi_{\perp}^{2}\left(\Delta
h_{\\!\perp}\right)^{2}\right]d^{2}x.$ (52)

When $\xi_{\perp}=0$, Eq. (52) represents the generalization for $T>0$ of the
heuristic Eq. (9). When $T=0$ $h_{||}$ is constant and Eq. (52) reduces to Eq.
(9). Moreover, in fact, Eq. (52) tends to Eq. (9) when $T\to 0$. For $T>0$
monopoles appear, and thus the longitudinal height function is also present.

###  V.2 Low Temperature

When $T\downarrow 0$, fluctuations of the entropic should fields diverge. One
could then perform a proper study of the higher orders expansion of
$\Omega(\phi,\psi)$ in terms of Feynman diagrams. Because of the complex shape
of $\Omega(\phi,\psi)$ such program is challenging. It might also be
uninteresting in real systems.

That is because at very small $T$ our model becomes insufficient for real
systems. Firstly, in magnetic systems long range interactions become relevant
at low temperature, and it can induce ordering at equilibrium Melko, den
Hertog, and Gingras (2001). Secondly, regardless of the range of interaction,
creation of monopoles is suppressed when $T\ll\epsilon$. Therefore, spin flips
that do not require large (compared to $T$) activation energy either pertain
to spins impinging on a monopole and that move the monopole—but monopoles are
exponentially few at low $T$—or to collaborative flips of entire loops, which
are exponentially unlikely in the size of the loop. Thus, one expects a
freeze-in as $T\downarrow 0$, in analogy with what is seen in pyrochlores
Snyder  _et al._ (2001), but our approach is at equilibrium.

We will therefore intend $T$ as small but not too small and proceed by
assuming that the functional form of an effective theory is quadratic, has the
same functional form as the theory above, but that interactions among
fluctuations lead to a “dressing” in the constants.

Because correlation lengths diverge at low $T$, from Eq. (44) we have e.g. for
charges

$\langle
q^{2}\rangle\sim\chi_{0}\xi_{||}^{-2}~{}~{}~{}~{}~{}~{}~{}\text{for}~{}T\downarrow
0,$ (53)

which implies a redefinition of constants in the theory at low $T$:
$\epsilon,\kappa$ are dressed as

$\displaystyle\epsilon\to\epsilon_{r}(T)$ $\displaystyle\sim T\chi_{0}/\langle
q^{2}\rangle$ $\displaystyle\kappa\to\kappa_{r}(T)$ $\displaystyle\sim
T\chi_{0}/\langle i^{2}\rangle~{}~{}~{}~{}~{}~{}~{}\text{for}~{}T\downarrow
0.$ (54)

Note that this dressing correspond to a quadratic theory that has the correct
mean square charge and current already from equipartition. Note also that Eq.
(53) has the form of a Debye screening length for a 2D-Coulomb potential whose
coupling constant is proportional to $T$, which is the case for our entropic
potential. We have shown elsewhere Nisoli (2020a) that Eq. (53) can follow
from a Debye-Hückel approach to the entropic potentials (see also Appendix I).

Note finally that $\chi_{0}$ must be finite at low $T$ (see also below), and
that $\langle q^{2}\rangle$ is decently approximated by the naive
$\overline{q^{2}}$ computed by assuming uncorrelated vertices, each with
proper multiplicity and Boltzmann weight. At low $T$,
$\overline{q^{2}}\sim(16/3)\exp(-2\epsilon/T)$. We have thus
$\xi_{||}\sim\exp\left({\frac{\epsilon}{T}}\right)$ for $T\downarrow 0$. A
similar exponential behavior for the correlation length was indeed suggested
experimentally by analyzing the pinch points in the structure factor of
pryrochlore iceFennell  _et al._ (2009). This exponential singularity points
to the topological nature of the $T=0$ ice-manifold.

In conclusion, at low $T$ we can use all the equations of the previous
sections, if expressed in terms of $\xi_{||}$, $\xi_{\\!\perp}$, where

$\xi_{||}^{2}\sim\begin{dcases}\epsilon/T~{}~{}~{}~{}~{}~{}\text{for}~{}~{}T/\mu\uparrow\infty\\\
\chi_{0}/\langle q^{2}\rangle~{}~{}\text{for}~{}~{}T/\mu\downarrow
0\end{dcases}$ (55)

and similarly

$\xi_{\\!\perp}^{2}\sim\begin{dcases}\kappa/T~{}~{}~{}~{}~{}\text{for}~{}~{}T/\kappa\uparrow\infty\\\
\chi_{0}/\langle i^{2}\rangle~{}~{}\text{for}~{}~{}T/\kappa\downarrow
0\end{dcases},$ (56)

and for $\chi_{0}$, $\chi_{0}\sim 1$ for $T\uparrow\infty$ while it remains
finite at $T=0$ (We show below that e.g for $T=0,\kappa=0$ we have
$\chi_{0}=2$.)

In Appendix A, we further motivate how this choice is mathematically
reasonable. Note also that no Kosterlitz-Thouless transition of monopole or
current unbinding is present Kosterlitz and Thouless (1973) because the
interaction is entropic. Transitions are driven by the interplay between
temperature and energy, but here interaction is itself entropic and thus
thermal.

Figure 3: Experimental results of entropic screening of a pinned monopole in a
square ice realized in a quantum annealer, from ref King  _et al._ (2021).
Here $J\propto\epsilon/T$ (see reference for details).

.

##  VI Cases

###  VI.1 Pure Degenerate Spin Ice ($\kappa=0$)

In this case, the ground state, aka ice manifold, is the degenerate six-vertex
modelLieb (1967), described in section II. In the purest form—i.e. without the
interference of long-range interaction—it was recently realized in a quantum
annealer King  _et al._ (2021), where the entropic effects were cleanly
studied. Nanomagnetic realizations Möller and Moessner (2006); Perrin, Canals,
and Rougemaille (2016); Farhan  _et al._ (2019) imply long range interaction,
whose ulterior effects we have described elsewhere Nisoli (2020a).

When $\kappa=0$, the functional integration over $\prod_{p}di_{p}$ in Eq. (23)
produces the delta functions $\prod_{p}\delta(\psi_{p})$. Further integration
over $\prod_{p}d\psi_{p}$ returns for the partition function

$Z=\int\left[dq\right]\tilde{\Omega}[q]\mathbbm{e}^{-\beta{\cal
H}[q]+\sum_{v}q_{v}V_{q,v}},$ (57)

where $\tilde{\Omega}[q]$ is the density of states for the charges, and it is
given by

$\displaystyle\tilde{\Omega}[q,i]=\int[d\phi]\Omega[\phi]\mathbbm{e}^{\sum_{v}q_{v}\mathbbm{i}\phi_{v}},$
(58)

and is therefore the functional Fourier transform of

$\Omega[\phi]=2^{N_{e}}{\prod_{{\hat{vv^{\prime}}}}}\cosh\left(-\mathbbm{i}\nabla_{vv^{\prime}}\phi-\beta\nabla_{pp^{\prime}}V_{i}+\beta
H_{vv^{\prime}}\right).$ (59)

In other words, the currents disappear from the picture, and the entropic
potential $-\mathbbm{i}T\psi$ is replaced in the equations by the external
potential $V_{i}$ acting on currents. Then, proceeding by the quadratic
approximation as before, and using the third of Eq. (21), one obtains for the
currents correlations

$\langle|\tilde{i}(\vec{k})|^{2}\rangle=\chi_{0}\gamma(\vec{k})^{2},$ (60)

which for small wave vectors reduces to Eq. (16).

####  VI.1.1 Charge Correlations

More generally, the reader will find that all the equations of the theory
developed above apply to the pure ice case by taking $\xi_{\perp}=0$. E.g.,
from Eq. (LABEL:taus), $\tilde{\chi}_{\\!\perp}(\vec{k})=\chi_{0}$, and the
perpendicular susceptibility in real space is a delta function, consistent
with zero correlation length.

Figure 4: Structure factors (divided by $\chi_{0}$) at high (top) and low
(middle) temperature plotted from Eq. (49), for pure degenerate square ice.
Figure 5: Structure factor for pure degenerate square ice across the
$k_{x}=2\pi$ line (top) shows sharpening as the correlation length increases,
leading to pinch points, and across the $k_{x}=\pi$ line (bottom).

Considering charge correlations and screening, Eq. (42) can be rewritten as

$\langle|\tilde{q}(k)|^{2}\rangle=\chi_{0}\xi_{||}^{-2}\left(1-\frac{1}{1+\xi_{||}^{2}\gamma^{2}}\right)$
(61)

where the first term Fourier-transforms to a Kronecker delta, whereas the
second term implies at large distance the charge correlation

$\displaystyle\langle
q_{v_{1}}q_{v_{2}}\rangle=-\frac{\chi_{0}}{2\pi\xi_{||}^{4}}K_{0}\left({|v_{1}-v_{2}|}/{\xi_{||}}\right).$
(62)

Note that the modified Bessel function $K_{0}$ is the screened 2D-Coulomb
potential. It is exponentially screened, or
$K_{0}(x)\sim\sqrt{\pi/2x}\exp(-x)$ making $\xi_{||}$ the
correlation/screening length. In 3D, we would have a screened 3D-Coulomb or
$\exp(-|v_{1}-v_{2}|/\xi_{||})/|v_{1}-v_{2}|$. This result is general: for
spin ice on a generic graph, for which Laplacian operators can also be
defined, entropic interactions lead to screened Coulomb correlationsNisoli
(2020b). Note that, from Eq (LABEL:taus), the kernel of the longitudinal
susceptibility functional is also a screened 2D-Coulomb, given by the the
modified Bessel function $K_{0}$ with correlation length $\xi_{||}$.

In artificial realizations it is possible to pin a charge $Q_{\mathrm{pin}}$
in $v_{0}$. Then, it is easy to show that the pinned charge generates a charge
distribution

$\displaystyle\langle q_{v}\rangle=\frac{\langle
q_{v}q_{v_{0}}\rangle}{\langle q^{2}\rangle}Q_{\mathrm{pin}}.$ (63)

Figure 6: Structure factors at different temperatures plotted from Eq. (49),
for the weak FM square ice. Note the line state emerging, the persistence of
pinch points, which are maxima. Figure 7: Structure factor for the weak FM
square ice across the $k_{x}=\pi$ line (top) shows line state features, and
across the $k_{x}=2\pi$ line (bottom) shows sharpening as the correlation
length increases, leading to pinch points.

Remarkably, the screening comes entirely from the entropic interaction. This
has been verified in a quantum annealer, where charges can be pinned King  _et
al._ (2021). Figure 3 reports experimental results of screening, which verify
Eq. (62). At high $\epsilon/T$ the curve becomes flat, consistent with Eqs.
(61, 62) as the correlation length exceeds the finite size of the sample 555
The situation becomes considerably more complex in an impure square ice. If
monopoles interact also via a real 3D-Coulomb law, various screening regimes
are predicted Nisoli (2020a), due to the interplay of the screening length and
the Bjerrum length..

We can gain some knowledge of screening at small distances by approximating
around the $K=(\pm\pi,\pm\pi)$ points of the BZ. There, $\gamma(k)^{2}$ is
maximum and $\gamma(\vec{k}+K)^{2}=8-k^{2}$. This leads to a screening
function

$\langle|\tilde{q}(K+k)|^{2}\rangle\simeq\xi_{||}^{-2}\left(1-\frac{\xi_{||}^{-2}}{\xi_{p}^{-2}-k^{2}}\right)$
(64)

where $\xi_{p}^{2}=\xi_{||}^{2}/(1+8\xi_{||}^{2})$. Then for small $\Delta
v=v-v^{\prime}$ the charge correlation (or equivalently screening) has a sign
alternation with the Manhattan distance on the graph and an envelope function
$E(|v-v^{\prime}|/\xi_{p})$ of periodicity $\xi_{p}/2\pi$, of the form

$\displaystyle\langle q_{v}q_{v^{\prime}}\rangle=(-1)^{\Delta v_{x}+\Delta
v_{y}}E(|\Delta v|/\xi_{p})+\text{constant}.$ (65)

Such sign alternation with the Manhattan distance for screening at short
distance was also verified in experiments King  _et al._ (2021) on a quantum
annealer666Note, however, that monopoles are characterized not only by a
charge, but also by a net magnetic moment, which was not taken into account in
our formalism. Therefore, at short distances the screening would be
anisotropic, just as monopoles are. In a future work we will incorporate that
degree of freedom to study at short length or equivalently around the $K$
points of the BZ..

####  VI.1.2 Spin Correlations

Considering now the spin correlations, Eq. (46) particularizes now to

$\displaystyle\chi_{0}^{-1}\langle\tilde{S}^{*}_{\alpha}(\vec{k})\tilde{S}_{\alpha^{\prime}}(\vec{k})\rangle$
$\displaystyle=\delta_{\alpha\alpha^{\prime}}-\hat{\gamma}_{\alpha}\hat{\gamma}_{\alpha^{\prime}}\left(1-\tilde{\chi}_{||}/\chi_{0}\right)$
$\displaystyle=\delta_{\alpha\alpha^{\prime}}-\hat{\gamma}_{\alpha}\hat{\gamma}_{\alpha^{\prime}}\xi_{||}^{2}|\tilde{q}(k)|^{2}$
$\displaystyle=\delta_{\alpha\alpha^{\prime}}-\frac{\gamma_{\alpha}\gamma_{\alpha^{\prime}}}{\xi_{||}^{-2}+\gamma^{2}}.$
(66)

In real space at distances larger than lattice discretization we have

$\chi_{0}^{-1}\langle
S_{\alpha}(\vec{x})S_{\alpha^{\prime}}(\vec{0})\rangle\simeq-2\pi\partial^{2}_{\alpha\alpha^{\prime}}K_{0}(|\vec{x}|/\xi_{||}),$
(67)

which is not algebraic at non-zero temperature. At long wavelength, the spin
correlations in real space in the pure ice manifold are thus the kernel of a
screened, 2D dipolar interaction, i.e. obtained as partial derivatives of the
screened 2D-Coulomb interaction ($K_{0}$), thus generalizing the heuristic Eq.
(14).

Spin correlations become algebraic in the (unrealistic) $T=0$ manifold, where
$\xi_{||}^{-1}=0$. Then, for the pure ice manifold we have the familiar
transverse form

$\displaystyle\chi_{0}^{-1}\langle\tilde{S}^{*}_{\alpha}(\vec{k})\tilde{S}_{\alpha^{\prime}}(\vec{k})\rangle_{\text{IM}}={{}^{\\!\perp}\\!\hat{\gamma}_{\alpha}}\\!\\!{{}^{\\!\perp}\\!\hat{\gamma}_{\alpha^{\prime}}}=\delta_{\alpha\alpha^{\prime}}-\frac{\vec{\gamma}_{\alpha}\vec{\gamma}_{\alpha^{\prime}}}{\gamma^{2}}.$
(68)

That is because, the longitudinal susceptibility $\tilde{\chi}_{||}$ goes to
zero in the ice manifold, as it is associated to monopoles, which disappear.
Thus only the transverse part of Eq. (46) remains. Indeed, temperature adds a
longitudinal part to the spin correlations, such that the difference between
correlations with and without temperature is

$\langle\tilde{S}^{*}_{\alpha}(\vec{k})\tilde{S}_{\alpha^{\prime}}(\vec{k})\rangle-\langle\tilde{S}^{*}_{\alpha}(\vec{k})\tilde{S}_{\alpha^{\prime}}(\vec{k})\rangle_{\text{IM}}=\chi_{0}\hat{\gamma}_{\alpha}\hat{\gamma}_{\alpha^{\prime}}\tilde{\chi}_{||}.$
(69)

In Fig. 4 we plot the structure factor $\Sigma_{m}(\vec{k})$ in units of
$\chi_{0}$ from Eq. (66), demonstrating the formation of pinch points as $T$
is reduced, shown in Fig. 5, top. Note that they compare very favorably with
structure factors obtained experimentally in a quantum annealer (Fig. 2 of ref
King  _et al._ (2021)) in the case of a line state, as well as from
simulations rou ; Brunn  _et al._ (2021). Interestingly, Fig. 5 bottom show
light bumps in the $K$ points of the BZ, that are also see, though more
pronounced, in experimental and numerical results, and correspond to small AFM
domains of t-I in the Ice Manifold, whose net AFM staggered order parameter is
nonetheless zero (no symmetry breaking).

Finally, spin correlations allow us to say something about $\chi_{0}$. Because
spins have values $S_{\alpha}=\pm 1$, then $1=\langle
S^{2}_{\alpha}\rangle=N_{v}^{-1}\langle\sum_{e}S^{2}_{\alpha,e}\rangle=\int_{\text{BZ}}\langle|\tilde{S}_{\alpha}(\vec{k})|^{2}\rangle_{\text{IM}}\frac{d^{2}k}{(2\pi)^{2}}$,
and summing over $\alpha=x,y$, from Eq. (68) we obtain

$\chi_{0}=2~{}~{}~{}~{}\text{for}~{}~{}~{}~{}T=0.$ (70)

Proceeding in the same way, but using the second line of Eq. (66) we obtain
for generic temperature

$\chi_{0}=1+\frac{1}{2}\xi_{||}^{2}\langle q^{2}\rangle,$ (71)

which relates the uniform susceptibility to the correlation length and the
mean square charge, at all temperatures. Then from Eqs. (55, 71) we find again
$\chi_{0}\to 2$ for $T\downarrow 0$ and $\chi_{0}\to 1$ fo $T\uparrow+\infty$.

Figure 8: Plot of $\langle|\tilde{i}(\vec{k})|^{2}\rangle\chi_{0}^{-1}$ along
the line $k_{x}=k_{y}$ for different negative values of $\kappa/T$
($\kappa/T=0$ in red, $\kappa/T=0.08$ in blue, $\kappa/T=0.1$ in black)
showing the growing divergence at the $K$ points $k_{x}=k_{y}=\pi$.

Figure 9: Structure factors at different temperatures plotted from Eq. (49),
for the weakly AFM square ice when $\xi_{\perp}=0.1\xi_{||}\mathbbm{i}$. Note
the formation of peaks at the $K$ points of the Brillouin zone indicating FM
ordering.

###  VI.2 Weak Ferromagnetic Case ($\kappa\simeq 0^{+}$)

This case has an ordered, fourfold, ferromagnetic ground state, corresponding
to the four possible full polarizations of the system. We are interested here
in small $\kappa/\epsilon$ (hence “weak”), at temperatures above the ordering,
whose scale is set by $\kappa$, but below the crossover into the ice manifold,
whose scale is set by $\epsilon$. While this case does not describe at low
temperature square ice with heigh offset exceeding the critical offset Möller
and Moessner (2006); Perrin, Canals, and Rougemaille (2016), which instead
leads to a disordered, albeit subextensively so, line liquid, it is likely a
good approximation for it at intermediate temperatures.

Because of the gauge-free duality, everything said above for charge
correlations applies here both to charges and currents. Thus, on top of the
previous equations, there are also screened correlations for currents as

$\displaystyle\langle
i_{p}i_{0}\rangle=-\frac{1}{2\pi\xi_{\perp}^{2}}K_{0}\left({|p|}/{\xi_{\\!\perp}}\right).$
(72)

In particular, a pinned current generates a screening from other currents,
similar to the screening of a charge.

While the considerations in the previous subsections on charge correlations
still apply, the spin correlations are now different and are given by the full
Eq. (46). The structure factor (Figz. 6, 7) remains reminiscent of spin ice.
As temperature is reduced it also shows the features of the line state (see
Fig. 2e in ref rou  and Fig 2 in ref King  _et al._ (2021)). However, it also
show a maximum at the center of the BZ that signals incipient ferromagnetic
ordering, absent in the standard line state, where instead it is replaced by
maxima on the $k_{x},k_{y}$ axes.

These features compare very favorably with those in structure factors obtained
experimentally in a quantum annealer (Fig. 2 of ref King  _et al._ (2021)) in
the case of a line state, as well as from simulations in refs rou ; Brunn  _et
al._ (2021).

###  VI.3 Weak Antiferromagnetic Case ($\kappa\simeq 0^{-}$)

When $k<0$, ${\cal H}_{2}$ in Eq. (39) is not bounded from below when $T\leq
T_{c}^{\text{afm}}=8|\kappa|$ (because $\gamma^{2}$ reaches its maximum
$\gamma^{2}=8$ on the $K$ corners of the BZ). This merely testifies to the
expected ordering criticality in the AFM case. Currents are promoted by a
negative $\kappa$ and the AFM state is an ice rule state that maximizes
currents, made of a tessellation of t-I vertices. There is a second order
phase transition to such ordered state. Because we have employed a mean field
theory, $T_{c}^{\text{afm}}$ is not the actual critical temperature, but
merely a useful parameter in the context of our framework.

For $T>T_{c}^{\text{afm}}$, $\tilde{\chi}_{\\!\perp}(k)$ has a maximum on the
$K$ corners of the BZ. Thus, we expand around $K$,
$\gamma(K+\vec{k})^{2}\simeq 8-k^{2}$ in $\tilde{\chi}_{\\!\perp}(k)$ and from
Eq (42) we obtain for large $|p|$

$\displaystyle\langle
I_{p}I_{0}\rangle\simeq\frac{(T/\kappa)^{2}}{2\pi}(-1)^{p_{x}+p_{y}}K_{0}\left({|p|}/{\xi_{\text{afm}}}\right),$
(73)

which expectedly alternates sign on adiacent plaquettes.

The AFM correlation length is given by

$\xi_{\text{afm}}^{2}=|\kappa|/(T-T_{c}^{\text{afm}})$ (74)

and is a measure of the size of the antiferromagnetic domains. It diverges at
the ordering criticality, though with the wrong critical exponent. Since one
expects a 2D Ising transitionWu (1969); Anghinolfi  _et al._ (2015), the
exponent should be $1$. Here it is $1/2$, or the mean-field exponent, as we
have employed a quadratic mean-field approximation.

Because $\kappa/\epsilon$ is small, $T^{\text{afm}}$ is much smaller than the
crossover temperature for the ice regime (typically $\sim 2\epsilon$). Thus,
above $T_{c}^{\text{afm}}$ the divergence-free and divergence-full fields
behave independently and features of the essentially transversal IM structure
factor, such as pinch points and charge correlations, are still present. We
plot the structure factor in Fig. 9. Note there the. growing maxima at the $K$
points of the BZ, corresponding to AFM ordering.

Again, these features compare very favorably with those in structure factors
obtained experimentally in a quantum annealer (Fig. 2 of ref King  _et al._
(2021)) in the case of a line state, as well as from simulations in refs rou ;
Brunn  _et al._ (2021).

##  VII Naïve Kinetics

The formalism above can be used for general considerations about kinetics on
time scales comparable to the relaxation times or shorter Chaikin, Lubensky,
and Witten (1995). We anticipate here the main deductions, leaving a more in-
depth study to future work on “slow” electrodynamics of spin ice.

###  VII.1 General Considerations

Unlike equilibrium thermodynamics, which merely concerns itself with samplings
of the phase space, dynamics is generally system-specific. Therefore, many
different models of kinetics could be introduced, in or out of equilibrium,
closer to the constitutive properties of a specific material or more general.
In this subsection, we consider notions of “kinematics” that hold true in any
such model.

As noted already Ryzhkin (2005), simple considerations should readily convince
that the flow density vector 777We use here “flow density vector” to denote
what is more generally called “current density vector”, to avoid confusion
with the notion of currents $i$ previously introduced. for the charges, in or
out of equilibrium, is always given by

$\vec{J}^{q}=\dot{\vec{S}},$ (75)

in perfect analogy with the relation between electrical current and the
dielectric polarization vector in standard electrodynamics.

Reasoning in the long wavelength limit, and taking the divergence we obtain
the conservation equation for the charge

$\dot{q}=-\vec{\nabla}\cdot\vec{J}^{q},$ (76)

which is merely the derivative of Eq. (2). Similarly, taking instead the curl,
we have

$\dot{i}=\hat{e}_{3}\cdot\vec{\nabla}\wedge\vec{J}^{q}=-\vec{\nabla}\cdot\vec{J}^{i}$
(77)

with

$\vec{J}^{i}=\hat{e}_{3}\wedge\vec{J}^{q}.$ (78)

We found that the flow density vectors for charges $q$ and currents $i$ are
orthogonal, as expected by gauge-free duality. Thus a flow of monopoles always
implies a flow of currents, and vice versa.

Even a distracted look at the physical system should convince that there
cannot be steady states of charge flow. After the application of a uniform
field, the spins reorient. This can be interpreted as flow of monopoles and
currents. At equilibrium all the spins have reached the new configuration, the
system is static and there is no more flow. Direct current is therefore only
possible during that relaxation, and it is not steady 888Depending on the
system, the aging might be slow rather than exponential Castelnovo, Moessner,
and Sondhi (2010); Levis and Cugliandolo (2013), in which case tiny currents
might be detected after a long time..

###  VII.2 Relaxation Dynamics

As mentioned, a simple relaxation dynamics at equilibrium would not apply to
real materials, as it is clear both experimentally Bovo  _et al._ (2013);
Dusad  _et al._ (2019); Goryca  _et al._ (2021) and conceptually Castelnovo,
Moessner, and Sondhi (2010); Levis and Cugliandolo (2013); Nisoli (2020c).
Depending on the system, however, deviations can be limited to kinetics much
faster than relaxation Nisoli (2020c).

A fundamental aspect is the choice of the degrees of freedom that relax. In
our system it seems reasonable to choose the spins as “real variables”, and
not charges or currents, nor, as done by e.g. C. Henley Henley (1997) in the
context of the similar dimer model, the height functions. Thus we proceed with
the relaxation equation for the spins $S$ which we write as

$\tau_{0}\dot{\vec{S}}=-\beta\frac{\delta{\cal F}[\vec{S}]}{\delta\vec{S}}$
(79)

where $\tau_{0}$ is a characteristic time of the kinetics, and changes with
temperature.

####  VII.2.1 Dynamic Susceptibilities, Conductivities, Dispersion

From the Eqs. (79, 50) we obtain for the longitudinal and transverse spin
components the equations of motion in reciprocal space.

$\displaystyle\tau_{0}\dot{\tilde{S}}_{||,\perp}(\vec{k})=-\tilde{\chi}_{||,\perp}(\vec{k})^{-1}\tilde{S}_{||,\perp}(\vec{k})+\beta\tilde{H}_{||,\perp}.$
(80)

Thus, the longitudinal and transverse relaxation times are

$\tau_{||,\perp}(\vec{k})=\tau_{0}\tilde{\chi}_{||,\perp}(\vec{k}).$ (81)

and thus the relaxation frequencies (from $\nu=1/\tau$) obey the dispersion
relations

$\nu_{||,\perp}(\vec{k})=\frac{1}{\chi_{0}\tau_{0}}+D_{||,\perp}\gamma(\vec{k})^{2},$
(82)

where $D_{||},D_{\\!\perp}$ are diffusion constants (see below) given by

$D_{||,\perp}=\frac{\xi^{2}_{||,\perp}}{\tau_{0}\chi_{0}}.$ (83)

[Also, from Eq. (83, 81) one finds Eq. (7.5) of ref Bramwell (2012).]

From Eq. (LABEL:taus) we see that when $\epsilon>0,\kappa>0$, the relaxation
time is minimal at the $K$ points of the Brillouin zone, which corresponds to
processes that involve few, neighboring spins, e.g. the creation/annihilation
of a monopole pair. When $\kappa=0$ the transverse relaxation time
$\tau_{\\!\perp}(\vec{k})$ is flat and equal to $\tau_{0}$ because processes
concerning currents do not change the energy.

From Eq. (80) we find the dynamic susceptibilities (defined for dimensional
convenience as $\tilde{S}=\tilde{\chi}\beta\tilde{H}$) are

$\displaystyle\tilde{\chi}_{||,\perp}(\omega,\vec{k})=\frac{\chi_{||,\perp}(\vec{k})}{1-\mathbbm{i}\omega\tau_{0}\chi_{||,\perp}(\vec{k})}$
(84)

as already appreciated via other methods in pyrochlores spin ice Bramwell
(2012).

The dynamic susceptibilities in Eq. (84) are typical of an exponential
relaxation to equilibrium Topping and Blundell (2018). They are wrong at high
frequency. Indeed, Eq. (84) and the fluctuation-dissipation theorem imply that
the power spectrum is a Lorentzian, and correspond to a brown noise, i.e.
scaling as $1/\omega^{2}$ at large $\omega$. Instead, recent experiments
report a “color” of the noise that depends on temperature and varies between
brown and pink Dusad  _et al._ (2019); Goryca  _et al._ (2021), suggesting
instead a generalized Debye function Cole and Cole (1941), further
underscoring the expected limits at frequencies of a mean field relaxation
kineticsNisoli (2020c).

From the first of Eq. (80), and remembering $\tilde{q}=\gamma\tilde{S}_{||}$,
$\tilde{i}=\gamma\tilde{S}_{\\!\perp}$ we obtain the relaxation equation for
the charge as

$\displaystyle\tau_{0}\dot{\tilde{q}}(\vec{k})=-\left[1+\xi_{||}^{2}\gamma(\vec{k})^{2}\right]\tilde{q}(\vec{k})-\mathbbm{i}\beta\vec{\gamma}\cdot\tilde{H}$
$\displaystyle\tau_{0}\dot{\tilde{i}}(\vec{k})=-\left[1+\xi_{\\!\perp}^{2}\gamma(\vec{k})^{2}\right]\tilde{q}(\vec{k})-\mathbbm{i}\beta^{\perp\\!}\vec{\gamma}\cdot\tilde{H},$
(85)

Finally, from Eqs. (75-84) we have, for the longitudinal and transverse
conductivity, defined as $J^{q}=\sigma H$, the expression

$\displaystyle\tilde{\sigma}_{||,\perp}(\omega,\vec{k})=\frac{1}{T\tau_{0}}\frac{-\mathbbm{i}\omega\tau_{||,\perp}(\vec{k})}{1-\mathbbm{i}\omega\tau_{||,\perp}(\vec{k})},$
(86)

which are zero at $\omega=0$, as it should be: there are no direct currents.

####  VII.2.2 Real Space Picture

Equations (85) become, in real space at length scales larger than the lattice
constant, diffusion-relaxation equations

$\displaystyle\dot{q}=D_{||}\Delta q+\left(\beta\nu
q_{\text{ext}}-q\chi_{0}^{-1}\right)\tau_{0}^{-1}$
$\displaystyle\dot{i}=D_{\\!\perp}\Delta i+\left(\beta\nu
i_{\text{ext}}-i\chi_{0}^{-1}\right)\tau_{0}^{-1}.$ (87)

The first term in Eqs. (87) accounts for the diffusion of charges or currents.
The diffusion flow density vectors are

$\displaystyle\vec{J}^{q,\text{diff}}=-D_{||}\vec{\nabla}q$
$\displaystyle\vec{J}^{i,\text{diff}}=-D_{\\!\perp}\vec{\nabla}i,$ (88)

for charges and currents respectively.

The second term in Eqs. (87), in absence of external charges, accounts for
pair annihilation of charges and relaxes exponentially the charge to the
equilibrium value $\langle q\rangle=0$. In presence of external charges the
system relaxes to the screening Equations (34).

A similar diffusion-relaxation equation can be written for $\vec{S}$. From Eq.
(80) we have

$\displaystyle\dot{S_{\alpha}}$
$\displaystyle=D_{||}\nabla_{\alpha}\vec{\nabla}\cdot\vec{S}+D_{\\!\perp}\\!\left[\vec{\nabla}\wedge\left(\vec{\nabla}\wedge\vec{S}\right)\right]_{\alpha}\\!\\!$
$\displaystyle+\left(\beta
H_{\alpha}-S_{\alpha}\chi_{0}^{-1}\right)\tau_{0}^{-1}$ (89)

which can also be written as

$\displaystyle\dot{S_{\alpha}}$
$\displaystyle=-\nabla^{\alpha^{\prime}}J^{S}_{{\alpha^{\prime}}{\alpha}}+\left(\beta
H_{\alpha}-S_{\alpha}\chi_{0}^{-1}\right)\tau_{0}^{-1}$ (90)

in terms of the flow density tensor for magnetization

$\displaystyle
J^{S}_{{\alpha^{\prime}}{\alpha}}=-D_{||}\nabla_{\alpha}S_{\alpha^{\prime}}-D_{\\!\perp}\left(\nabla_{\alpha}S_{\alpha^{\prime}}-\nabla_{\alpha^{\prime}}S_{\alpha}\right).$
(91)

However magnetization is not conserved, but because of the second term
$-S_{\alpha}$, it decays exponentially to the equilibrium value $\langle
S_{\alpha}\rangle=0$ imposed by the $Z_{2}$ symmetry of the problem (if
$\vec{H}=\vec{0}$).

Equation (89) can also be written in terms of charges and currents

$\displaystyle\dot{\vec{S}}$
$\displaystyle=-D_{||}\vec{\nabla}q-D_{\\!\perp}\hat{e}_{3}\wedge\vec{\nabla}i+\left(\beta
H_{\alpha}-S_{\alpha}\chi_{0}^{-1}\right)\tau_{0}^{-1},$ (92)

from which Eqs. (87) can be deduced directly by taking the divergence and the
curl. Then, in Eq. (92) we recognize the diffusion flow density vectors of Eq.
(88). We can therefore write for the total flow density vector of the charge

$\vec{J}^{q}=\vec{J}^{q,\text{diff}}+\hat{e}_{3}\wedge\vec{J}^{i,\text{diff}}+\left(\beta
H_{\alpha}-S_{\alpha}\chi_{0}^{-1}\right)\tau_{0}^{-1}.$ (93)

The second term is clearly divergence-free. From Eq. (78) the flow density
vector of the currents is obtained as its rotation.

The third term in Eqs. (92,93) is a drift term to which magnetization is
subtracted. Something similar was previously found in pyrochlore spin icevia
considerations of chemical physics of electrolytes Ryzhkin (2005); Bramwell
(2012). At equilibrium, for a constant uniform field $\vec{H}$, that term is
zero, currents and charges are uniform, and thus there is no longer monopole
flow. Again, direct current of monopoles is only possible during relaxation.

When $\kappa=0$, $D_{\perp}=0$ and the picture simplifies. Currents no longer
diffuse, but relax uniformly. However from Eq. (78) there is a divergence free
flow of currents, if monopoles are flowing.

##  VIII Conclusions

We have illustrated the duality between charges and currents in pure spin ice,
that is square spin ice coupled only at the vertex level. We have built on it
a field theory where elementary currents and monopoles are the degrees of
freedom, while spins are subsumed into entropic interactions. In pure spin
ice, this leads to a 2D electrodynamics formalism where 2D-Coulomb
interactions are entropic.

Within this framework, we have deduced free energies, static and dynamic
susceptibilities, relaxation times, form factors, structure factors, for the
three cases of: degenerate spin ice, line state, antiferromagnetic square ice.
Our purely analytical results compare well with a wealth of experimental and
numerical data. They generalize to thermal states previous heuristic
approaches at zero $T$ variously developed in the context of the dimer model,
and based on the height function formalism. They also accord well with
chemical physics approaches to 3D pyrochlores.

##  Appendix A Conceptualizing the Dressing at Low $T$

To understand the mathematical origin of the dressing, consider that
$\Omega[\phi,\psi]$ is periodic in the gradient of the entropic fields. Such
periodicity comes from the sum over the spin ensemble of the Fourier transform
of Dirac deltas. The latter enforce the discrete nature of charges and
currents. Thus, $\Omega[\phi,\psi]$ has two roles: one is to convey the
entropic interaction, and the other is to preserve the information that
charges are discrete. In our high temperature limit, however, we had taken
$\nabla\phi,\nabla\psi$ small, thus losing the periodicity needed to constrain
the magnitude of charges. In this scenario, the low temperature dressing [Eqs.
(54)] within an effective quadratic theory takes care of that constraint
already at the level of equipartition, while maintaining a formalism of
continuum charge distribution.

To separate the two effects in $\Omega$ consider the following, simple
calculation. We performe it in all generality in coordination $z$ on a
bipartite lattice of alternating vertices $A-B$. We take for simplicity
$\kappa=0$ (and thus integration over currents enforces $\psi=0$) and
$\vec{H}=0$.

$\Omega[\phi]$ in Eq. (25) can be written as products on $A$ vertices $v_{a}$,
or

$\Omega[\phi]=2^{N_{l}}\prod_{v_{a}}\omega[\phi_{v_{a}}]$ (94)

with

$\omega[\phi_{v_{a}}]=\prod_{v_{b}\in\partial
v_{a}}\cos(\phi_{v_{a}}-\phi_{v_{b}}),$ (95)

where $\partial v_{a}$ is the set of vertices $B$, connected to $v_{a}$. The
next step is to perform the mean field approximation
$\cos(\phi_{v_{a}}-\phi_{v_{b}})\simeq\cos(\phi_{v_{a}}-\overline{\phi}_{v_{a}})$,
where $\overline{\phi}_{v_{a}}$ is the mean of the fields $\phi_{v_{b}}$ on
vertices $v_{b}$ neighboring $v_{a}$, so that
$\omega[\phi_{v}]\simeq\cos(\phi_{v}-\overline{\phi}_{v})^{z}$. That leads to

$\int
d\phi_{v}e^{-\mathbbm{i}\phi_{v}q_{v}}\omega[\phi_{v}]=e^{\mathbbm{i}\overline{\phi}_{v}q_{v}}\tilde{\omega}(q_{v})$
(96)

where

$\tilde{\omega}(q)=2^{-z}\sum_{n=0}^{z}{{z}\choose{n}}\delta(q-q_{n})$ (97)

restricts $q$ to the only possible charges $q_{n}=(2n-z)$ for
$n=-z,-z+1,\dots,z-1,z$, each with proper multiplicity ${{z}\choose{n}}$.

We have thus obtained a mean field approach that, unlike the high $T$
approach, separates in $\Omega[\phi]$ the effect of the entropic interaction,
now expressed via the mean field $\bar{\phi}$, from the enforcement of
discreteness of charges, which is necessary at very low $T$ where excitations
are sparse.

Because at low $T$ the correlation length is much larger than the lattice
constant, it is physically intuitive to take $i\overline{\phi}_{v}=\langle
i\phi_{v}\rangle=V^{e}_{v}$, the entropic field, which is real. Then, the
$A\leftrightarrow B$ symmetry and the previous equations imply that the
probability of a charge distribution $q_{v}$ can be factored as

$\displaystyle\rho[q]=\prod_{v}\rho_{v}(q_{v}),$ (98)

in terms of the probability of having a charge $q$ on a vertex $v$, given by

$\rho_{v}(q)=\tilde{\omega}(q)\exp\left(-\frac{\epsilon}{T}q^{2}+qV^{e}_{v}\right)/Z_{v}.$
(99)

In Eq. (99) correlations are transmitted among vertices by the entropic field
$V^{e}_{v}$. $Z_{v}$ normalizes $\rho_{v}(q)$ and depends on $v$ via
$V^{e}_{v}$, which in turns depends on the collective charge distribution and
is therefore non-local, or $V^{e}_{v}=V^{e}_{v}[q]$. Finally, while Eq. (99)
was here deduced from a mean field approximation of the exact Eq. (23), it can
however also be obtained directly à la Landau, so to speak, as we have shown
previously Nisoli (2014).

What we are missing now is an information on how the entropic field depends on
the charge distribution. Assuming that at low $T$ charges are sparse, that
correlation lengths are large, and that the coarse graining of the geometry is
isotropic (as it is the case for the square lattice), we can write the
equation at lowest order in charges, fields and derivatives

$\overline{q}=\chi_{0}\nabla^{2}V^{e}$ (100)

where $\overline{q}_{v}=\int dq\rho_{v}(q)q$. From it and Eq. (99) we leave it
as a simple exercise for the reader to show that linearization à la Debye-
Hückel leads to the same correlations and screening as in the previous
subsection, but with $\xi_{||}$ replaced by Eq. (53) (see also ref Nisoli
(2020a)).

###### Acknowledgements.

We thank Andrew King (D-Wave Systems) for useful discussions and Beatrice
Nisoli for proofreading. This work was carried out under the auspices of the
U.S. DoE through the Los Alamos National Laboratory, operated by Triad
National Security, LLC (Contract No. 892333218NCA000001).
Data Availability Statement Data sharing is not applicable to this article as
no new data were created or analyzed in this study.

## References

  * Lieb (1967) E. H. Lieb, “Residual entropy of square ice,” Physical Review 162, 162 (1967).
  * Ryzhkin (2005) I. Ryzhkin, “Magnetic relaxation in rare-earth oxide pyrochlores,” Journal of Experimental and Theoretical Physics 101, 481–486 (2005).
  * Castelnovo, Moessner, and Sondhi (2008) C. Castelnovo, R. Moessner, and S. L. Sondhi, “Magnetic monopoles in spin ice,” Nature 451, 42–5 (2008).
  * Bramwell and Harris (2020) S. T. Bramwell and M. J. Harris, “The history of spin ice,” Journal of Physics: Condensed Matter 32, 374010 (2020).
  * Ramirez  _et al._ (1999) A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. S. Shastry, “Zero-point entropy in ‘spin ice’,” Nature 399, 333–335 (1999).
  * den Hertog and Gingras (2000) B. C. den Hertog and M. J. Gingras, “Dipolar interactions and origin of spin ice in ising pyrochlore magnets,” Physical review letters 84, 3430 (2000).
  * Bramwell and Gingras (2001) S. T. Bramwell and M. J. Gingras, “Spin ice state in frustrated magnetic pyrochlore materials,” Science 294, 1495–501 (2001).
  * Tanaka  _et al._ (2006) M. Tanaka, E. Saitoh, H. Miyajima, T. Yamaoka, and Y. Iye, “Magnetic interactions in a ferromagnetic honeycomb nanoscale network,” Physical Review B 73, 052411 (2006).
  * Wang  _et al._ (2006) R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi, and P. Schiffer, “Artificial ’spin ice’ in a geometrically frustrated lattice of nanoscale ferromagnetic islands,” Nature 439, 303–6 (2006).
  * Nisoli, Moessner, and Schiffer (2013) C. Nisoli, R. Moessner, and P. Schiffer, “Colloquium: Artificial spin ice: Designing and imaging magnetic frustration,” Reviews of Modern Physics 85, 1473 (2013).
  * Heyderman and Stamps (2013) L. Heyderman and R. Stamps, “Artificial ferroic systems: novel functionality from structure, interactions and dynamics,” Journal of Physics: Condensed Matter 25, 363201 (2013).
  * Skjærvø  _et al._ (2019) S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, “Advances in artificial spin ice,” Nature Reviews Physics , 1–16 (2019).
  * Ortiz-Ambriz  _et al._ (2019) A. Ortiz-Ambriz, C. Nisoli, C. Reichhardt, C. J. Reichhardt, and P. Tierno, “Colloquium: Ice rule and emergent frustration in particle ice and beyond,” Reviews of Modern Physics 91, 041003 (2019).
  * Morrison, Nelson, and Nisoli (2013) M. J. Morrison, T. R. Nelson, and C. Nisoli, “Unhappy vertices in artificial spin ice: new degeneracies from vertex frustration,” New Journal of Physics 15, 045009 (2013).
  * Stamps (2014) R. L. Stamps, “Artificial spin ice: The unhappy wanderer,” Nature Physics 10, 623–624 (2014).
  * Nisoli, Kapaklis, and Schiffer (2017) C. Nisoli, V. Kapaklis, and P. Schiffer, “Deliberate exotic magnetism via frustration and topology,” Nature Physics 13, 200–203 (2017).
  * Li  _et al._ (2018) Y. Li, G. W. Paterson, G. M. Macauley, F. S. Nascimento, C. Ferguson, S. A. Morley, M. C. Rosamond, E. H. Linfield, D. A. MacLaren, R. Macêdo, _et al._ , “Superferromagnetism and domain-wall topologies in artificial “pinwheel” spin ice,” ACS nano 13, 2213–2222 (2018).
  * Saccone  _et al._ (2019) M. Saccone, K. Hofhuis, Y.-L. Huang, S. Dhuey, Z. Chen, A. Scholl, R. V. Chopdekar, S. van Dijken, and A. Farhan, “Dipolar cairo lattice: Geometrical frustration and short-range correlations,” Physical Review Materials 3, 104402 (2019).
  * Makarova  _et al._ (2021) K. Makarova, V. Strongin, I. Titovets, A. Syrov, I. Zinchenko, V. Samoylov, K. Hofhuis, M. Saccone, A. Makarov, A. Farhan, _et al._ , “Low-energy states, ground states, and variable frustrations of the finite-size dipolar cairo lattices,” Physical Review E 103, 042129 (2021).
  * Gliga  _et al._ (2017) S. Gliga, G. Hrkac, C. Donnelly, J. Büchi, A. Kleibert, J. Cui, A. Farhan, E. Kirk, R. V. Chopdekar, Y. Masaki, _et al._ , “Emergent dynamic chirality in a thermally driven artificial spin ratchet,” Nature materials 16, 1106–1111 (2017).
  * Stopfel  _et al._ (2018) H. Stopfel, E. Östman, I.-A. Chioar, D. Greving, U. B. Arnalds, T. P. Hase, A. Stein, B. Hjörvarsson, and V. Kapaklis, “Magnetic order and energy-scale hierarchy in artificial spin-ice structures,” Physical Review B 98, 014435 (2018).
  * Schiffer and Nisoli (2021) P. Schiffer and C. Nisoli, “Artificial spin ice: Paths forward,” Applied Physics Letters 118, 110501 (2021).
  * Arnalds  _et al._ (2012) U. B. Arnalds, A. Farhan, R. V. Chopdekar, V. Kapaklis, A. Balan, E. T. Papaioannou, M. Ahlberg, F. Nolting, L. J. Heyderman, and B. Hjörvarsson, “Thermalized ground state of artificial kagome spin ice building blocks,” Applied Physics Letters 101, 112404 (2012).
  * Farhan  _et al._ (2013) A. Farhan, P. M. Derlet, A. Kleibert, A. Balan, R. V. Chopdekar, M. Wyss, J. Perron, A. Scholl, F. Nolting, and L. J. Heyderman, “Direct observation of thermal relaxation in artificial spin ice,” Physical review letters 111, 057204 (2013).
  * Kapaklis  _et al._ (2014) V. Kapaklis, U. B. Arnalds, A. Farhan, R. V. Chopdekar, A. Balan, A. Scholl, L. J. Heyderman, and B. Hjörvarsson, “Thermal fluctuations in artificial spin ice,” Nature nanotechnology 9, 514–519 (2014).
  * Bernal and Fowler (1933) J. Bernal and R. Fowler, “A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions,” The Journal of Chemical Physics 1, 515–548 (1933).
  * Nisoli  _et al._ (2007) C. Nisoli, R. Wang, J. Li, W. McConville, P. Lammert, P. Schiffer, and V. Crespi, “Ground State Lost but Degeneracy Found: The Effective Thermodynamics of Artificial Spin Ice,” Phys. Rev. Lett. 98, 217203 (2007).
  * Nisoli  _et al._ (2010) C. Nisoli, J. Li, X. Ke, D. Garand, P. Schiffer, and V. H. Crespi, “Effective Temperature in an Interacting Vertex System: Theory and Experiment on Artificial Spin Ice,” Phys. Rev. Lett. 105, 047205 (2010).
  * Wu (1969) F. Y. Wu, “Critical behavior of two-dimensional hydrogen-bonded antiferroelectrics,” Phys. Rev. Lett. 22, 1174–1176 (1969).
  * Anghinolfi  _et al._ (2015) L. Anghinolfi, H. Luetkens, J. Perron, M. Flokstra, O. Sendetskyi, A. Suter, T. Prokscha, P. Derlet, S. Lee, and L. Heyderman, “Thermodynamic phase transitions in a frustrated magnetic metamaterial,” Nature communications 6 (2015).
  * Morgan  _et al._ (2010) J. P. Morgan, A. Stein, S. Langridge, and C. H. Marrows, “Thermal ground-state ordering and elementary excitations in artificial magnetic square ice,” Nat. Phys. 7, 75–79 (2010).
  * Porro  _et al._ (2013) J. Porro, A. Bedoya-Pinto, A. Berger, and P. Vavassori, “Exploring thermally induced states in square artificial spin-ice arrays,” New Journal of Physics 15, 055012 (2013).
  * Zhang  _et al._ (2013) S. Zhang, I. Gilbert, C. Nisoli, G.-W. Chern, M. J. Erickson, L. O?Brien, C. Leighton, P. E. Lammert, V. H. Crespi, and P. Schiffer, “Crystallites of magnetic charges in artificial spin ice,” Nature 500, 553–557 (2013).
  * Sendetskyi  _et al._ (2019) O. Sendetskyi, V. Scagnoli, N. Leo, L. Anghinolfi, A. Alberca, J. Lüning, U. Staub, P. M. Derlet, and L. J. Heyderman, “Continuous magnetic phase transition in artificial square ice,” arXiv preprint arXiv:1905.07246 (2019).
  * Möller and Moessner (2006) G. Möller and R. Moessner, “Artificial Square Ice and Related Dipolar Nanoarrays,” Phys. Rev. Lett. 96, 237202 (2006).
  * Giauque and Ashley (1933) W. Giauque and M. F. Ashley, “Molecular rotation in ice at 10 k. free energy of formation and entropy of water,” Physical review 43, 81 (1933).
  * Pauling (1935) 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 (1935).
  * Perrin, Canals, and Rougemaille (2016) Y. Perrin, B. Canals, and N. Rougemaille, “Extensive degeneracy, coulomb phase and magnetic monopoles in artificial square ice,” Nature 540, 410–413 (2016).
  * Farhan  _et al._ (2019) A. Farhan, M. Saccone, C. F. Petersen, S. Dhuey, R. V. Chopdekar, Y.-L. Huang, N. Kent, Z. Chen, M. J. Alava, T. Lippert, _et al._ , “Emergent magnetic monopole dynamics in macroscopically degenerate artificial spin ice,” Science advances 5, eaav6380 (2019).
  * Nascimento  _et al._ (2012) F. S. Nascimento, L. A. S. Mol, W. A. Moura-Melo, and A. R. Pereira, “From confinement to deconfinement of magnetic monopoles in artificial rectangular spin ices,” New Journal of Physics 14, 115019 (2012).
  * Ribeiro  _et al._ (2017) I. Ribeiro, F. Nascimento, S. Ferreira, W. Moura-Melo, C. Costa, J. Borme, P. Freitas, G. Wysin, C. de Araujo, and A. Pereira, “Realization of rectangular artificial spin ice and direct observation of high energy topology,” arXiv preprint arXiv:1704.07373 (2017).
  * Östman  _et al._ (2018) E. Östman, H. Stopfel, I.-A. Chioar, U. B. Arnalds, A. Stein, V. Kapaklis, and B. Hjörvarsson, “Interaction modifiers in artificial spin ices,” Nature Physics 14, 375 (2018).
  * Schanilec  _et al._ (2019) V. Schanilec, Y. Perrin, S. L. Denmat, B. Canals, and N. Rougemaille, “Artificial vertex systems by design,” arXiv preprint arXiv:1902.00452 (2019).
  * Macêdo  _et al._ (2018) R. Macêdo, G. Macauley, F. Nascimento, and R. Stamps, “Apparent ferromagnetism in the pinwheel artificial spin ice,” Physical Review B 98, 014437 (2018).
  * King  _et al._ (2021) A. D. King, C. Nisoli, E. D. Dahl, G. Poulin-Lamarre, and A. Lopez-Bezanilla, “Qubit spin ice,” Science 373, 576–580 (2021).
  * (46) .
  * Lieb and Wu (1972) E. Lieb and F. Wu, “Two-dimensional ferroelectric models,” (1972).
  * Baxter (1982) R. Baxter,  _Exactly solved models in statistical mechanics_ (Academic, New York, 1982).
  * Nisoli (2020a) C. Nisoli, “Equilibrium field theory of magnetic monopoles in degenerate square spin ice: Correlations, entropic interactions, and charge screening regimes,” Physical Review B 102, 220401 (2020a).
  * Henley (2010) C. L. Henley, “The coulomb phase in frustrated systems,” Annu. Rev. Condens. Matter Phys. 1, 179–210 (2010).
  * Isakov  _et al._ (2004) S. Isakov, K. Gregor, R. Moessner, and S. Sondhi, “Dipolar spin correlations in classical pyrochlore magnets,” Physical Review Letters 93, 167204 (2004).
  * Garanin and Canals (1999) D. Garanin and B. Canals, “Classical spin liquid: Exact solution for the infinite-component antiferromagnetic model on the kagomé lattice,” Physical Review B 59, 443 (1999).
  * Henley (2005) C. Henley, “Power-law spin correlations in pyrochlore antiferromagnets,” Physical Review B 71, 014424 (2005).
  * Youngblood and Axe (1981) R. Youngblood and J. Axe, “Polarization fluctuations in ferroelectric models,” Physical Review B 23, 232 (1981).
  * Huse  _et al._ (2003) D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, “Coulomb and liquid dimer models in three dimensions,” Physical review letters 91, 167004 (2003).
  * Henley (2011) C. L. Henley, “Classical height models with topological order,” Journal of Physics: Condensed Matter 23, 164212 (2011).
  * Henley (1997) C. L. Henley, “Relaxation time for a dimer covering with height representation,” Journal of statistical physics 89, 483–507 (1997).
  * Bramwell (2012) S. T. Bramwell, “Generalized longitudinal susceptibility for magnetic monopoles in spin ice,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 5738–5766 (2012).
  * Twengström, Henelius, and Bramwell (2020) M. Twengström, P. Henelius, and S. T. Bramwell, “Screening and the pinch point paradox in spin ice,” Physical Review Research 2, 013305 (2020).
  * Jaccard (1964) C. Jaccard, “Thermodynamics of irreversible processes applied to ice,” Physik der kondensierten Materie 3, 99–118 (1964).
  * Nisoli (2020b) C. Nisoli, “The concept of spin ice graphs and a field theory for their charges,” AIP Advances 10, 115102 (2020b).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, “Ordering, metastability and phase transitions in two-dimensional systems,” Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Lamberty, Papanikolaou, and Henley (2013) R. Z. Lamberty, S. Papanikolaou, and C. L. Henley, “Classical topological order in abelian and non-abelian generalized height models,” Physical review letters 111, 245701 (2013).
  * Note (1) We call it current because a magnetization $\mathaccentV{vec}17EM$ generates an electrical current density $\mathaccentV{vec}17Ej=\mathaccentV{vec}17E\nabla\wedge\mathaccentV{vec}17EM$.
  * van Beijeren (1977) H. van Beijeren, “Exactly solvable model for the roughening transition of a crystal surface,” Physical Review Letters 38, 993 (1977).
  * Chui and Weeks (1976) S. Chui and J. Weeks, “Phase transition in the two-dimensional coulomb gas, and the interfacial roughening transition,” Physical Review B 14, 4978 (1976).
  * Zinn-Justin (2009) P. Zinn-Justin, “Six-vertex, loop and tiling models: integrability and combinatorics,” arXiv preprint arXiv:0901.0665 (2009).
  * Raghavan, Henley, and Arouh (1997) R. Raghavan, C. L. Henley, and S. L. Arouh, “New two-color dimer models with critical ground states,” Journal of statistical physics 86, 517–550 (1997).
  * Sutherland (1968) B. Sutherland, “Correlation functions for two-dimensional ferroelectrics,” Physics Letters A 26, 532–533 (1968).
  * Note (2) We use dyadics: if $\mathaccentV{vec}17Ev$ is a vector, the components of $\mathaccentV{vec}17Ev\mathaccentV{vec}17Ev$ are $[\mathaccentV{vec}17Ev\mathaccentV{vec}17Ev]_{ij}\mathrel{\mathop{:}}=v_{i}v_{j}$; $\mathbbm{1}$ is the unitary matrix.
  * Nagle (1966) J. F. Nagle, “Lattice statistics of hydrogen bonded crystals. i. the residual entropy of ice,” Journal of Mathematical Physics 7, 1484–1491 (1966).
  * Nagle (1978) J. Nagle, “Configurational statistics,” Journal of Glaciology 21, 73–83 (1978).
  * Nagle (1979) J. Nagle, “Theory of the dielectric constant of ice,” Chemical Physics 43, 317–328 (1979).
  * Ryzhkin and Whitworth (1997) I. Ryzhkin and R. Whitworth, “The configurational entropy in the jaccard theory of the electrical properties of ice,” Journal of Physics: Condensed Matter 9, 395 (1997).
  * Perrin, Canals, and Rougemaille (2019) Y. Perrin, B. Canals, and N. Rougemaille, “Quasidegenerate ice manifold in a purely two-dimensional square array of nanomagnets,” Physical Review B 99, 224434 (2019).
  * Levis  _et al._ (2013) D. Levis, L. F. Cugliandolo, L. Foini, and M. Tarzia, “Thermal phase transitions in artificial spin ice,” Physical review letters 110, 207206 (2013).
  * Note (3) When monopoles interact via a $1/r$ law, $\xi_{||}$ is no longer a correlation length, as the monopole-monopole interaction destroys the screening at least in principle, as we have shown elsewhere Nisoli (2020a).
  * Brooks-Bartlett  _et al._ (2014) M. Brooks-Bartlett, S. T. Banks, L. D. Jaubert, A. Harman-Clarke, and P. C. Holdsworth, “Magnetic-moment fragmentation and monopole crystallization,” Physical Review X 4, 011007 (2014).
  * Petit  _et al._ (2016) S. Petit, E. Lhotel, B. Canals, M. C. Hatnean, J. Ollivier, H. Mutka, E. Ressouche, A. Wildes, M. Lees, and G. Balakrishnan, “Observation of magnetic fragmentation in spin ice,” Nature Physics 12, 746 (2016).
  * Canals  _et al._ (2016) B. Canals, I.-A. Chioar, V.-D. Nguyen, M. Hehn, D. Lacour, F. Montaigne, A. Locatelli, T. O. Menteş, B. S. Burgos, and N. Rougemaille, “Fragmentation of magnetism in artificial kagome dipolar spin ice,” Nature communications 7 (2016).
  * Note (4) $\mathaccentV{hat}05E\gamma$ is the unit vector of $\mathaccentV{vec}17E\gamma$.
  * Melko, den Hertog, and Gingras (2001) R. G. Melko, B. C. den Hertog, and M. J. Gingras, “Long-range order at low temperatures in dipolar spin ice,” Physical review letters 87, 067203 (2001).
  * Snyder  _et al._ (2001) J. Snyder, J. Slusky, R. Cava, and P. Schiffer, “How ‘spin ice’freezes,” Nature 413, 48–51 (2001).
  * Fennell  _et al._ (2009) T. Fennell, P. Deen, A. Wildes, K. Schmalzl, D. Prabhakaran, A. Boothroyd, R. Aldus, D. McMorrow, and S. Bramwell, “Magnetic coulomb phase in the spin ice ho2ti2o7,” Science 326, 415–417 (2009).
  * Note (5) The situation becomes considerably more complex in an impure square ice. If monopoles interact also via a real 3D-Coulomb law, various screening regimes are predicted Nisoli (2020a), due to the interplay of the screening length and the Bjerrum length.
  * Note (6) Note, however, that monopoles are characterized not only by a charge, but also by a net magnetic moment, which was not taken into account in our formalism. Therefore, at short distances the screening would be anisotropic, just as monopoles are. In a future work we will incorporate that degree of freedom to study at short length or equivalently around the $K$ points of the BZ.
  * Brunn  _et al._ (2021) O. Brunn, Y. Perrin, B. Canals, and N. Rougemaille, “Signatures of farther neighbor couplings in artificial square ice,” Physical Review B 103, 094405 (2021).
  * Chaikin, Lubensky, and Witten (1995) P. M. Chaikin, T. C. Lubensky, and T. A. Witten,  _Principles of condensed matter physics_ , Vol. 10 (Cambridge university press Cambridge, 1995).
  * Note (7) We use here “flow density vector” to denote what is more generally called “current density vector”, to avoid confusion with the notion of currents $i$ previously introduced.
  * Note (8) Depending on the system, the aging might be slow rather than exponential Castelnovo, Moessner, and Sondhi (2010); Levis and Cugliandolo (2013), in which case tiny currents might be detected after a long time.
  * Bovo  _et al._ (2013) L. Bovo, J. Bloxsom, D. Prabhakaran, G. Aeppli, and S. Bramwell, “Brownian motion and quantum dynamics of magnetic monopoles in spin ice,” Nature communications 4, 1–8 (2013).
  * Dusad  _et al._ (2019) R. Dusad, F. K. Kirschner, J. C. Hoke, B. R. Roberts, A. Eyal, F. Flicker, G. M. Luke, S. J. Blundell, and J. S. Davis, “Magnetic monopole noise,” Nature 571, 234–239 (2019).
  * Goryca  _et al._ (2021) M. Goryca, X. Zhang, J. Li, A. Balk, J. Watts, C. Leighton, C. Nisoli, P. Schiffer, and S. Crooker, “Field-induced magnetic monopole plasma in artificial spin ice,” Physical Review X 11, 011042 (2021).
  * Castelnovo, Moessner, and Sondhi (2010) C. Castelnovo, R. Moessner, and S. Sondhi, “Thermal quenches in spin ice,” Physical review letters 104, 107201 (2010).
  * Levis and Cugliandolo (2013) D. Levis and L. F. Cugliandolo, “Defects dynamics following thermal quenches in square spin ice,” Physical Review B 87, 214302 (2013).
  * Nisoli (2020c) C. Nisoli, “The colored noise of spin ice: Subdiffusion of magnetic monopoles in a coevolving vacuum,” arXiv preprint arXiv:2004.11965 (2020c).
  * Topping and Blundell (2018) C. Topping and S. Blundell, “Ac susceptibility as a probe of low-frequency magnetic dynamics,” Journal of Physics: Condensed Matter 31, 013001 (2018).
  * Cole and Cole (1941) K. S. Cole and R. H. Cole, “Dispersion and absorption in dielectrics i. alternating current characteristics,” The Journal of chemical physics 9, 341–351 (1941).
  * Nisoli (2014) C. Nisoli, “Dumping topological charges on neighbors: ice manifolds for colloids and vortices,” New Journal of Physics 16, 113049 (2014).