We consider pure square spin ice, that is, square ice, where only nearest neighbors are coupled. The 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 the height formalism of the disordered ground state to non-zero temperature. The framework can be applied to the zoology of recent experimental results, especially realizations on quantum annealers, and can be expanded to include longer range interactions.

## I. INTRODUCTION

Degenerate artificial square ice is perhaps the simplest two-dimensional system in which the disorder constrained by the ice rule^{1} and its violations as monopole excitations are studied.^{2,3} While spin ice pyrochlores^{4–7} had opened new vistas in the study of geometric frustration and constrained disorder, more recently, *artificial spin ice—*systems of interacting, magnetic nanoislands^{8–13}—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 nature^{14–21} and revealing the novel phenomena absent in its crystal analog,^{12,22} recent fabrication and characterization advances^{23–25} have brought back the austere simplicity of celebrated early models to the fore, such as kagome and square ice.

The square geometry of spin ice was among^{8} the first to be realized artificially,^{9} when it was shown that non-ice rule^{26} vertices are suppressed after AC demagnetization.^{9,27,28} However, because 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 class^{29,30} toward an ordered antiferromagnetic ground state.^{31–34}

Soon after, Möller and Moessner^{35} proposed to offset the height of half of the nanoislands to regain the ice rule degeneracy.^{36,37} The idea was recently realized.^{38,39} Meanwhile, ice rule degeneracy in square ice has also been demonstrated in rectangular lattices,^{40,41} via cleverly placed interaction modifiers placed in the vertices,^{42} with nano-holes in a connected spin ice of nanowires,^{43} or by rotation of the moments.^{44} Recently, pure square ice, i.e., square ice with only nearest neighbor interaction, was realized in a quantum annealer^{45} 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 of different spectral characterizations.^{46}

The scope of this work is to provide a unifying and expandable framework that can cover many different square ice systems around the ice rule degeneracy point^{35,45} 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 the gauge-free duality of the system, which is absent in three dimension (3D).

We limit ourselves to *pure* square ice, that is, a 16 vertex model, where interactions are limited to spins within the same vertex.^{29,47,48} We consider no long-range interactions and thus no 3D-Coulomb (i.e., 1/*r*) interaction among monopoles or currents (we have considered such cases elsewhere^{49}). However, because of the emergent nature of these objects, we show that they interact via 2D-Coulomb (i.e., ∼ln *r*) entropic interactions.

Within this model, we compute free energies, entropic interactions, structure factors, susceptibilities, correlations, screening, and relaxation dynamics. We discuss strengths and limits of this approach. We show new results but also re-derive results that were previously achieved 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.^{2,50–60} We also particularized some results that we had already found on generic graphs.^{61}

## II. SQUARE ICE AND ITS GAUGE-FREE DUALITY

In 2D, there is *gauge-free* duality absent in 3D pyrochlore ice. It is related to the rather gravid mathematical fact that in 2D, the Helmholtz decomposition has no gauge freedom. That in turn 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,^{62} and is behind the entire edifice of complex analysis.

### A. Gauge-free duality in 2D

In 2D, we can always write a continuum vector field $S\u20d7$ in terms of longitudinal and transverse potentials *h*_{‖} and *h*_{⊥},

where $e\u03023=e\u03021\u2227e\u03022$, $e\u03021,e\u03022$ are orthonormal bases of the plane, and we call $S\u20d7\Vert andS\u20d7\u22a5$ the longitudinal and perpendicular components of the field, respectively. Unlike in the 3D case, where the perpendicular part of $S\u20d7$ is the curl of a vector potential, there is no gauge freedom, as shown in Eq. (1).

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

then

If for a vector $w\u20d7$ we call

the *perpendicular* of $w\u20d7$, we then have

which express the duality between charges and currents, or longitudinal and perpendicular components of the field, under a *π*/2 rotation.

This duality has a discretized analog in square spin ice.

### B. Square spin ice

Square spin ice (Fig. 1) is a set of classical, binary spins $S\u20d7e$ aligned on the *N*_{e} edges *e* of a square lattice of *N*_{v} = *N*_{e}/2 vertices labeled *v*. Spins form four vertex topologies often classified^{9} as t-I, …, t-IV, where t-I and t-II obey the ice rule^{26,35,37,48} (i.e., have two spins pointing in and two pointing out).

In degenerate square ice, ice rule vertices are degenerate and energetically favored, and they lead to a ground state^{1} of constrained disorder, called the *ice manifold*. The latter is *a Coulomb phase*,^{50–53} i.e., a topological state labeled, *in lieu* of an order parameter, by a *height field*.^{56,63}

Consider *Q*_{v}[*S*], the topological charge of the vertex *v*, defined as the number of spins pointing to 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 that pointing counterclockwise.^{64}

For a spin configuration $S\u20d7$, consider its perpendicular configuration $S\u20d7\u22a5$ (Fig. 1), for which old plaquettes are now vertices and old vertices are now plaquettes. We then have

as in Eq. (5).

For each configuration of spins that obeys the ice rule, a unique (up to a constant) *height function* *h*_{⊥} can be defined on the plaquettes such that

is true ($pp\u2032\u0302$ is the unit vector pointing from plaquette *p* to *p*′ separated by the edge *e* as in Fig. 1). This follows from the fact that $S\u22a5$ is “irrotational:” the line-sum of spins $S\u20d7\u22a5$ (gray 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 *∀p*), a height function *h*_{‖v} can be defined on the vertices by

where *e* is the edge connecting the vertices *vv*′.

Thus, in the ice manifold, there are no charges, currents are disordered, and $S\u22a5$ is the discrete gradient of *h*_{⊥}, which “labels” the disorder of currents in the absence of charges. Conversely, in the current-free manifold, *S* is the discrete of *h*_{‖}, which labels the disorder of charges in the absence of currents.

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

### C. 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) or “flat” (ordered) phase, which is jargon derived from the theory of *roughening transition* historically associated with these models by various exact mappings.^{65,66}

The ice manifold of the square ice, i.e., the six-vertex model, is known to be equivalent to a dimer cover model^{48,67} and thus in a rough phase.^{50,56,57} A widespread,^{56,57,65,66} although by no means, rigorously justified ansatz ascribes to a configuration in the ice manifold, an entropy^{68} that is quadratic in the height function. In our language, we can write

where *h*_{⊥} is homogenized into a continuum field and *χ*_{0} is the positive uniform susceptibility (seen later). Clearly, the same is true for the current-free manifold by replacing *h*_{⊥} with *h*_{‖}.

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

where $H\u20d7$ is an external field.

Note that Eqs. (9) and (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 a Higgs boson is needed in the standard model). Equation (9) merely happens to work in reproducing correlations that can also in part be computed exactly^{48,69} (see also Ref. 57 and references therein for a discussion).

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

From that and the continuum limit of Eq. (7) (or $S\u20d7=\u2212e\u03023\u2227\u2207\u20d7h\u22a5$), one obtains the spin correlator as^{70}

which is purely transversal.

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

Spin correlations in real space can be obtained from Eq. (12) or more easily as partial transversal derivatives (or $\u2207\u20d7\u22a5=e\u03023\u2227\u2207\u20d7$) of Eq. (11), obtaining

where $P(x)$ is the kernel of the dipole–dipole interaction in 2D, or

Analogously, in the case of pyrochlore spin ice, the spin correlations are the kernel of the 3D dipolar interaction.^{50}

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) and (3), we have

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

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

The 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) and (1), the entropy of a configuration in the ice manifold can be written in terms of $S\u20d7$ as

with the constraint $\u2207\u20d7\u22c5S\u20d7=0$.

Equation (18), unlike the previous formulas, is also valid in 3D. There, it has been labeled as *Jaccard entropy* in water ice^{60,71–74} and later in pyrochlore spin ice^{50,51,53,55} 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 *χ*_{0} is thus related to the so-called Φ constant.^{74}

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 Gauss’s Law.^{45}

## III. ENERGY AND STATES

Section II C pertains to degenerate states, charge-free or current-free, but does not account for excitations and thermal fluctuations.^{1} For that, we need an energetic description.

The following Hamiltonian,

reflects the current-charge duality by placing a cost or advantage on topological currents and monopoles (*ϵ* and *κ* 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} = *ϵ* − *κ*, *J*_{2} = −*κ*, and *J*_{3} = *ϵ*. By the duality, the symmetry by orthogonalization corresponds to *ϵ* ↔ *κ*. The Hamiltonian describes various cases, often close to the experimental reality.

If *κ* = 0 and *ϵ* > 0, the ground state is the ice manifold shown in Fig. 1 (black arrows). Equivalently, by gauge-free duality, if *ϵ* = 0 and *κ* > 0, the ground state is an extensively degenerate ice manifold for $S\u20d7\u22a5$, shown through gray arrows in Fig. 1.

If *κ* = 0 and *ϵ* < 0, the ground state is the charge-full state, i.e., the ordered, antiferromagnetic tessellation of t-IV vertices. If *ϵ* = 0 and *κ* < 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 *κ* > 0 and *ϵ* > 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 fourfold ferromagnetic state made of t-II vertices ferromagnetically aligned (more loosely; as the ground state is current-free and charge-free, we have Δ*h*_{‖} = Δ*h*_{⊥} = 0, which implies a uniform $S\u20d7$.). This case has not been investigated experimentally, although it can certainly be realized in a quantum annealer.^{45} It might approximate, however, experimental situations where t-II vertices can be favored,^{38,42,75} leading to a line state that is disordered but of sub-extensive entropy.

For *κ* < 0 and *ϵ* > 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 |*κ*| describes early square ice realizations.^{9} Small |*κ*| might approximate spin ices that are designed so that ice rule vertices are degenerate^{35,38} but where the dipolar interaction still favors closed magnetic fluxes, thus promoting topological currents and an ordered ground state.

When *κ* = *ϵ*, we have *J*_{1} = 0, and the set of vertical and horizontal arrows becomes two decoupled systems. When *κ* = *ϵ* > 0, each subsystem is ferromagnetic, and the ground state is a fourfold fully polarized state. When *κ* = *ϵ* < 0, each subsystem is antiferromagnetic, and the ground state is a fourfold 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 *κ* = *ϵ*).

Here, we will not study the full monopole-current model of Eq. (19), which leads to a rich phase diagram in the *βϵ* × *βκ* plane (where *β* = 1/*T* and *T* is the temperature measured in units of energy) to be compared with other models.^{29,76} We will consider only *ϵ* > 0 and |*κ*|/*ϵ* to be small and investigate how ice manifold features are retained by small perturbations around the *spin ice point* *κ* = 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^{49,61} to include currents. The partition function from Eq. (19) reads

and it is the generator of correlations,

The fields $H\u20d7,Vq,Vi$ are measured in units of energy.

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

and then sum over the spins, obtaining

where $dqdi=(2\pi )\u2212Nv\u220fvdqv\u220fpdip$. $\Omega \u0303[q,i]$ is a generalized density of states for *q*_{v}and*i*_{p}, given by

and it is therefore the functional Fourier transform of

(the product runs on all the edges $e=vv\u2032\u0302$ once, and ∇_{vv′}*ϕ* : = *ϕ*_{v′} − *ϕ*_{v}, ∇_{pp′}*ψ* : = *ψ*_{p′} − *ψ*_{p}, and $Hvv\u2032:=H\u20d7e\u22c5vv\u2032\u0302$, while $pp\u2032\u0302=\u2212e\u0302z\u2227vv\u2032\u0302$).

Note that by construction, $\u27e8Qv1\u2026Qvn\u27e9=\u27e8qv1\u2026qvn\u27e9$, and $\u27e8Ip1\u2026Ipn\u27e9=\u27e8ip1\u2026ipn\u27e9$. Note also that the use of $H\u20d7$, *V*_{q}, and *V*_{i} is superabundant: $Vq\u2192Vq+Vq\u2032,Vi\u2192Vi+Vi\u2032$ is equivalent to $H\u20d7\u2192H\u20d7+\u2207\u20d7Vq\u2032\u2212e\u03023\u2227\u2207Vi\u2032$, 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

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

of the generalized free energy

Note that $F[\varphi ,\psi ]$ can have imaginary values. As a consequence, although the variables *ϕ*and*ψ* are themselves real, their expectation values ⟨*ϕ*⟩ and ⟨*ψ*⟩ are imaginary, and thus, the entropic fields $\u27e8Vqe\u27e9,\u27e8Vie\u27e9$ are real. Indeed, by integrating over *q* in Eq. (23) and applying the second and third equation in (21), one finds

when *κ* > 0. Similar Gaussian gymnastics also prove that

We thus see that while $i$*ϕ* and $i$*ψ* correlate charges and currents, the vector fields −$i$*T*∇_{vv′}*ϕ* and −$i$*T*∇_{pp′}*ψ* correlate spins that would otherwise be trivially paramagnetic.

In the spirit of field theory, we could therefore say the following: $Vqe$ is an entropic potential acting on charges, and $Vie$, on currents; correspondingly, there is an entropic field $B\u20d7e$ acting on the spins, which in the long wavelength approximation is given by

## V. FIELD THEORY: APPROXIMATIONS

### A. High temperature

When *ϵ* > 0 and *κ* > 0 (see Sec. VI C for *κ* ≤ 0), integrating over charges and currents in Eq. (23) returns Gaussians in the fields of variance ⟨*ϕ*^{2}⟩ = *ϵ*/*T* and ⟨*ψ*^{2}⟩ = *κ*/*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 *βH* constant), the Gaussian distributions in the field become Dirac deltas, and Eq. (23) becomes the standard “paramagnetic” partition function,

Note that for *H* = 0, Eq. (32) 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,

Then, from $q=\u2212\u2207\u20d7\u22c5S\u20d7$, $j=e\u03023\u22c5\u2207\u20d7\u2227S\u20d7$, and Eq. (29), we find screened Poisson equations for *q* and *I*,

(here, $\nu qext=\u2207\u20d7\u22c5H\u20d7$, where *q*_{ext} is the external charge, $\nu iext=\u2207\u20d7\u2227H\u20d7$, where *i*_{ext} is the external current, and *ν* is an energy), and therefore “correlation lengths”^{77} as

This expression for *ξ*_{‖} was already achieved in other works, for different geometries, via different means,^{52,53} and we had generalized it on a graph.^{61} Note that for *κ* < 0, *ξ*_{⊥} would be imaginary, pointing to a periodicity typical of the AFM ensemble, as we will discuss in section VIC. Instead, when *κ* = 0, the integral over *i* returns a functional delta function on *ψ*, and *ψ* disappears from the equations, which for all purpose are equivalent to taking *ξ*_{⊥} = 0.

#### 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* ≃ 2*ϵ* 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 Ω[*ϕ*, *ψ*] at quadratic order and Fourier transform via $gx=\u222bBZg\u0303(k\u20d7)e\u2212ik\u20d7\u22c5xd2k/(2\pi )2$, where BZ is the Brillouin zone, *g* is a generic field, and *x* = *v*, *p*, *l* represents vertices, edges, or plaquettes, respectively.

We obtain the approximated partition function at zero loop,

where

is an effective Hamiltonian for the new variables and the vector $\gamma \u20d7$ has components

for *α* = *x*, *y*, which contains information on the square symmetry of the lattice [in the long wavelength limit: $\gamma \u20d7\u2243k\u20d7+O(k3)$].

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 *ϕψ* cross term survives at quadratic order in the high *T* approximation, and charge and currents are independent at second order, leading to the so-called magnetic fragmentation.^{78–80}

The third line contains the coupling with the external field $H\u20d7$ (here, we neglect V). Because $\u2212i\gamma \u20d7\u22c5w\u0303\u20d7$ and $\u2212i\gamma \u20d7\u2227w\u0303\u20d7$ are the generalized divergence and curl, respectively, on the lattice in moment space for a generic field $w\u20d7$, the entropic fields for currents (resp. charges) couple to the curl (resp. divergence) of the external field.

The quantity *χ*_{0}, already encountered in Eq. (9) on the opposite limit (*T* = 0), is the uniform susceptibility with respect to β *H*. In Eq. (37), it is *χ*_{0} = 1, but we include it, nonetheless, in the equation because it can be *χ*_{0} ≠ 1 at low temperatures (see below).

Integrating *Z*_{2} over $\varphi \u0303,\psi \u0303$ when $H\u0303=0$ returns the effective free energies for *q* and *i* in the absence of external field, which are decoupled at second order,

where the two terms can be written as

In field theory language, the first terms in Eq. (40) are the “masses” of the monopole or of the current. The second terms imply 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 *γ*^{2} ≃ *k*^{2}), the entropic interactions for charges and currents are

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

#### 2. Correlations, susceptibilities, entropy, and height functions

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

where $\chi \u0303\Vert (k),\chi \u0303\u22a5(k)$ are defined as

Correctly, from Eq. (42), we have$\u27e8|q\u0303(0\u20d7)|2\u27e9=0$ since $\u27e8|q\u0303(0\u20d7)|2\u27e9=\u27e8\u2211vqv2\u27e9=0$ (and the same for currents). Figure 2 shows a heat map for $\u27e8|q\u0303(k\u20d7)|2\u27e9$.

Furthermore, if ⟨*q*^{2}⟩ is the average charge per vertex, or $\u27e8q2\u27e9=1Nv\u2211v\u27e8qv2\u27e9$, then

From it, and Eq. (42), we obtain ⟨*q*^{2}⟩ → 4 for *T* ↑ ∞, which is correct. Indeed, ⟨*q*^{2}⟩ = 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.

Note that $\chi \u0303\Vert (k),\chi \u0303\u22a5(k)$ are the longitudinal and perpendicular static susceptibilities (multiplied by *T*, thus accounting already for the Curie–Weiss law) for $S\u0303\Vert =S\u0303\u20d7\u22c5\gamma \u0302$ (the longitudinal, charge-full, current-fee part of the magnetization) and $S\u0303\u22a5=S\u0303\u20d7\u22c5\u22a5\gamma \u0302$ (the perpendicular, charge-free, and current-full magnetization), respectively.^{81} Indeed, by integrating the Gaussian integral in Eq. (36) and remembering that $\gamma \u0302\alpha \gamma \u0302\alpha \u2032+\gamma \u0302\alpha \u22a5\gamma \u0302\alpha \u2032\u22a5=\delta \alpha \alpha \u2032$, one obtains the total free energy as a function of the external field at lowest order as

From it, the spin correlations could be obtained as

In the limit *T* → ∞, Eq. (46) correctly returns to $\u27e8S\u0303\alpha (k\u20d7)S\u0303\alpha \u2032(k\u20d7)\u27e9\u2192\delta \alpha \alpha \u2032$, or uncorrelated spins.

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

and similarly, from (42),

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

which at small *k* corresponds to $\chi \u0303\u22a5$.

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

which in the long wavelength limit becomes

The first two terms in the previous equation are energetic. The third is the Jaccard entropy^{2,60} mentioned above. In the ice manifold, it is the only surviving term since *ξ*_{⊥} = 0 and $\u2207\u22c5S\u20d7=0$, thus returning to Eq. (18), previously found heuristically.

In Eq. (51), *ξ*_{⊥} = 0 reduces to the functional that was found elegantly via methods of chemical physics for pyrochlore ice by Bramwell.^{58,59} Indeed, when expressed in terms of spins $S\u20d7$ rather than currents and charges *i*and*q*, the two formalisms are expected to coincide.

Expressing the coarse grained $S\u20d7$ via height functions as in Eq. (1), we obtain the free energy for the height functions in the long wavelength limit, at quadratic order,

### B. Low temperature

When *T* ↓ 0, fluctuations of the entropic fields should diverge. One could then perform a proper study of the higher order expansion of Ω(*ϕ*, *ψ*) in terms of Feynman diagrams. Because ∈ Ω has the form of a locking potential, and at low temperature $\u2207\u20d7\varphi ,\u2207\u20d7\psi $ are no longer limited, vertex operators might intervene.^{99} Because of the complex shape of Ω(*ϕ*, *ψ*), such a 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. First, in magnetic systems, long-range interactions become relevant at low temperature, and they can induce ordering at equilibrium.^{82} Second, regardless of the range of interaction, creation of monopoles is suppressed when *T* ≪ *ϵ*. 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 or to collaborative flips of entire loops, which are exponentially unlikely because of the size of the loop. Thus, one expects a freeze-in as *T* ↓ 0, in analogy with what is seen in pyrochlores,^{83} 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 and has the same functional form as the theory mentioned above, but those interactions among fluctuations lead to a “dressing” in the constants.

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

which implies a redefinition of constants in the theory at low *T*: *ϵ*and*κ* are dressed as

Note that this dressing corresponds 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^{49} that Eq. (53) can follow from a Debye–Hückel approach to the entropic potentials (see also the Appendix).

Note finally that *χ*_{0} must be finite at low *T* (see also below) and that ⟨*q*^{2}⟩ is decently approximated by the naive $q2\u0304$ computed by assuming uncorrelated vertices, each with proper multiplicity and Boltzmann weight. At low *T*, $q2\u0304\u223c(16/3)exp(\u22122\u03f5/T)$. We thus have $\xi \Vert \u223cexp\u03f5T$ for *T* ↓ 0. A similar exponential behavior for the correlation length was indeed suggested experimentally by analyzing the pinch points in the structure factor of pyrochlore ice.^{84} 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 Sec. V A, if expressed in terms of *ξ*_{‖} and *ξ*_{⊥}, where

and similarly,

For *χ*_{0} and *χ*_{0} ∼ 1, *T* ↑ ∞ while it remains finite at *T* = 0 (we show below that, e.g., for *T* = 0and*κ* = 0, we have *χ*_{0} = 2).

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

## VI. CASES

### A. Pure degenerate spin ice (*κ* = 0)

In this case, the ground state, aka ice manifold, is the degenerate six-vertex model,^{1} described in Sec. II. In the purest form—i.e., without the interference of long-range interaction—it was recently realized in a quantum annealer,^{45} where the entropic effects were cleanly studied. Nanomagnetic realizations^{35,38,39} imply long-range interaction, whose ulterior effects have been described elsewhere.^{49}

When *κ* = 0, the functional integration over ∏_{p}*di*_{p} in Eq. (23) produces the delta functions ∏_{p}*δ*(*ψ*_{p}). Further integration of ∏_{p}*dψ*_{p} returns it to the following partition function:

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

which is therefore the functional Fourier transform of

In other words, the currents disappear from the picture, and the entropic potential −$i$*Tψ* is replaced in the equations by the external potential *V*_{i} acting on currents. Then, proceeding based on the quadratic approximation derived before and using the third of Eq. (21), for the current correlations, one obtains

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

#### 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 *ξ*_{⊥} = 0. For example, from Eq. (43), $\chi \u0303\u22a5(k\u20d7)=\chi 0$, and the perpendicular susceptibility in real space is a delta function, consistent with zero correlation length.

Considering charge correlations and screening, the first of Eqs. (42) can be rewritten as

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

Note that the modified Bessel function *K*_{0} is the screened 2D-Coulomb potential. It is exponentially screened, or $K0(x)\u223c\pi /2xexp(\u2212x)$, making *ξ*_{‖} the correlation/screening length, as presaged above. In 3D, we would have a screened 3D-Coulomb form,^{100} exp(−|*v*_{1} − *v*_{2}|/*ξ*_{‖})/|*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 correlations.^{61} Note that from Eq. (43), the kernel of the longitudinal susceptibility functional is also screened 2D-Coulomb, given by the modified Bessel function *K*_{0} with correlation length *ξ*_{‖}.

In artificial realizations, it is possible to pin a charge *Q*_{pin} in *v*_{0}. Then, it is easy to show that the pinned charge generates a charge distribution

Remarkably, the screening comes entirely from the entropic interaction. This has been verified in a quantum annealer, where charges can be pinned.^{45} Figure 3 reports experimental results of screening, which verify Eq. (62). At high *ϵ*/*T*, the curve becomes flat, consistent with Eqs. (61) and (62) as the correlation length exceeds the finite size of the sample.^{85}

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

where $\xi p2=\xi \Vert 2/(1+8\xi \Vert 2)$. Then for small Δ*v* = *v* − *v*′, the charge correlation (or equivalently screening) has a sign alternation with the Manhattan distance on the graph and an envelope function *E*(|*v* − *v*′|/*ξ*_{p}) of periodicity *ξ*_{p}/2*π* of the form

#### 2. Spin correlations

The spin correlations are now considered. Equation (46) particularizes now to

In real space at distances larger than lattice discretization, we have

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, which is obtained as partial derivatives of the screened 2D-Coulomb interaction (*K*_{0}), thus generalizing the heuristic equation (14).

Spin correlations become algebraic in the (unrealistic) *T* = 0 manifold, where $\xi \Vert \u22121=0$. Then, for the pure ice manifold, we have the familiar transverse form,

That is because the longitudinal susceptibility $\chi \u0303\Vert $ goes to zero in the ice manifold as it is associated with 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

In Fig. 4, we plot the structure factor $\Sigma m(k\u20d7)$ in units of *χ*_{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. 45) as well as from simulations.^{46,87} Interestingly, Fig. 5 bottom shows light bumps in the *K* points of the BZ that are also seen, although 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 *χ*_{0}. Because spins have values *S*_{α} = ±1, then $1=\u27e8S\alpha 2\u27e9=Nv\u22121\u27e8\u2211eS\alpha ,e2\u27e9=\u222bBZ\u27e8|S\u0303\alpha (k\u20d7)|2\u27e9IMd2k(2\pi )2$; summing over *α* = *x* and *y*, from Eq. (68), we obtain

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

### B. Weak ferromagnetic case (*κ* ≃ 0^{+})

This case has an ordered, fourfold, ferromagnetic ground state, corresponding to the four possible full polarizations of the system. Here, we are interested in small *κ*/*ϵ* (hence “weak”), at temperatures above the ordering, whose scale is set by *κ*, but below the crossover into the ice manifold, whose scale is set by *ϵ*. While this case does not describe square ice with height offset exceeding the critical offset at low temperature,^{35,38} which instead leads to a disordered, albeit sub-extensively 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

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

While the considerations in Sec. VI A 2 on charge correlations still apply, the spin correlations are now different and are given by the full equation (46). The structure factor (Figs. 6 and 7) remains reminiscent of spin ice. As temperature is reduced, it also shows the features of the line state (see Fig. 2e in Ref. 46 and Fig. 2 in Ref. 45). However, it also shows 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} and *k*_{y} axes.

### C. Weak antiferromagnetic case (*κ* ≃ 0^{−})

When *k* < 0, $H2$ in Eq. (39) is not bounded from below when $T\u2264Tcafm=8|\kappa |$ (because *γ*^{2} reaches its maximum *γ*^{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 *κ*, 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 an ordered state. Clearly, due to our mean field theory approximation, $Tcafm$ is not the actual critical temperature but merely a useful parameter in the context of our framework.

For $T>Tcafm$, $\chi \u0303\u22a5(k)$ has maxima on the *K* corners of the BZ. These maxima are visible in the correlation among currents, shown in Fig. (8), and grow as *T* approaches $Tcafm$.

Thus, we expand around *K*, $\gamma (K+k\u20d7)2\u22438\u2212k2$, in $\chi \u0303\u22a5(k)$, and from Eq. (42), for large |*p*|, we obtain

which expectedly alternates sign on adjacent plaquettes.

The AFM correlation length is given by

and it is a measure of the size of the antiferromagnetic domains. It diverges at the ordering criticality, although with the wrong critical exponent. Since one expects a 2D Ising transition,^{29,30} the exponent should be 1. Here, it is 1/2, as expected from a quadratic mean field approximation.

Because *κ*/*ϵ* is small, *T*^{afm} is much smaller than the crossover temperature for the ice regime (typically ∼2*ϵ*). Thus, above $Tcafm$, 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.

## VII. NAÏVE KINETICS

From the formalism mentioned above, we can deduce general considerations about kinetics on time scales of the order of mean field relaxation times.^{88} Here, we anticipate the main deductions, leaving a more in-depth study to future work on “slow” electrodynamics of spin ice. Unlike equilibrium thermodynamics, which merely sample the phase space, kinetics are system-specific. The nature of the moments, of their inversion mechanism, depends on the specific system, and 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 Sec. VII A, we consider notions of “kinematics” that hold true in any case. In the following ones, we explore simple relaxation dynamics.

### A. General considerations

As noted already,^{2} simple considerations should readily convince that the flow density vector^{89} for the charges, in or out of equilibrium, is given by

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,

Taking the curl instead, we have

with

Therefore, the flow density vectors for charges *q* and currents *i* are orthogonal, as expected from the gauge-free duality: 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. The dynamics of the relaxation toward the field can be interpreted in terms of 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 relaxation, and it is not steady.^{90}

### B. Relaxation dynamics

As mentioned, various aspects of kinetics in real spin ice materials escape simple relaxation dynamics, as it is clear both experimentally^{91–93} and conceptually.^{94–96} Depending on the system, however, deviations can be limited to processes much faster than relaxation.^{96}

A fundamental aspect of relaxation dynamics is the choice of the degrees of freedom that relax. In our system it seems reasonable to choose the spins and not the charges or the currents, nor, as done by, e.g., Henley^{57} 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

where *τ*_{0} is a characteristic time of the kinetics, and changes with temperature, and *ω*_{0} its associated frequency.

#### 1. Dynamic susceptibilities, conductivities, and dispersion

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

Thus, the longitudinal and transverse relaxation frequencies obey the dispersion relations,

where *D*_{‖}, *D*_{⊥} are diffusion constants (see below), given by

[Equations (81) have an analogue in pyrochlore ice, see Eq. (7.5) of Ref. 58.] Correspondently, the relaxation times are

When *κ* = 0, the transverse relaxation time $\tau \u22a5(k\u20d7)$ is flat and equal to *τ*_{0} because processes concerning currents do not change the energy.

From Eq. (80), the dynamic susceptibilities (defined for dimensional convenience as $S\u0303=\chi \u0303\beta H\u0303$) are

as already determined via other methods in pyrochlores spin ice.^{58}

The dynamic susceptibilities in Eq. (84) are wrong at high frequency. Indeed, Eq. (84) and the fluctuation-dissipation theorem imply that the power spectrum is a Lorentzian, corresponding to a brown noise, i.e., scaling as 1/*ω*^{2} at large *ω*. Instead, recent experiments report a “color” of the noise that depends on temperature and varies between brown and pink,^{92,93} suggesting instead a generalized Debye function,^{97} further underscoring the expected limits at frequencies of mean field relaxation kinetics.^{96}

From the first of Eq. (80), and remembering $q\u0303=\gamma S\u0303\Vert $ and $i\u0303=\gamma S\u0303\u22a5$, we obtain the relaxation equation for the charge as

Finally, from Eqs. (75)–(84), the longitudinal and transverse conductivity for charges, defined as *J*^{q} = *σH*, is given by

which is zero at *ω* = 0 as there are no direct currents.

#### 2. Real space picture

In real space at length scales larger than the lattice constant, Eqs. (85) become diffusion-relaxation equations,

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

for charges and currents, respectively.

The second term in Eq. (87) accounts for pair annihilation of charges leading to an exponential relaxation to the equilibrium value given by the screening equation (34).

A similar diffusion–relaxation equation can be written for $S\u20d7$. From Eq. (80), we have

which can also be written as

In terms of the flow density tensor for magnetization,

which, because of the second term −*S*_{α}, coming from the Jaccard entropy, relaxes exponentially.

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

from which Eq. (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). Then, from Eq. (75), we can therefore write the *total* flow density vector of the charge as

The second term is clearly divergence-free. From Eq. (77), the flow density vector of the currents is obtained.

The third term in Eqs. (92) and (93) is a drift term to which magnetization is subtracted. Something similar was previously found in pyrochlore spin ice via consideration of chemical physics of electrolytes.^{2,58} At equilibrium, for a constant uniform field $H\u20d7$, that term is zero, currents and charges are uniform, and thus, there is no longer monopole flow. We see once more that direct current of monopoles is only possible during relaxation.

When *κ* = 0, *D*_{⊥} = 0, and the picture simplifies. Currents no longer diffuse but relax uniformly. However, from Eq. (78), it can be seen that 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, and structure factors for the three cases of degenerate spin ice, line state, and 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.

## ACKNOWLEDGMENTS

We thank Andrew King (D-Wave Systems) and Roderich Moessner (Max Planck) 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), and funded by a grant of the DOE-LDRD office at LANL, Grant No. 20190430ER.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### APPENDIX: CONCEPTUALIZING THE DRESSING AT LOW *T*

To understand the mathematical origin of the dressing, consider that Ω[*ϕ*, *ψ*] 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, Ω[*ϕ*, *ψ*] plays 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 ∇*ϕ* and ∇*ψ* to be small, thus losing the periodicity needed to constrain the magnitude of charges. In this scenario, the low temperature dressing [Eq. (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 Ω, consider the following simple calculation. We perform it in all generality on a bipartite lattice of coordination *z* and of alternating vertices *A* − *B*. For simplicity, we take *κ* = 0 (and thus integration over currents enforces *ψ* = 0) and $H\u20d7=0$.

Ω[*ϕ*] in Eq. (25) can be written as products on *A* vertices *v*_{a}, or

with

where ∂*v*_{a} is the set of vertices *B*, connected to *v*_{a}. The next step is to perform the mean field approximation $cos(\varphi va\u2212\varphi vb)\u2243cos(\varphi va\u2212\varphi \u0304va)$, where $\varphi \u0304va$ is the mean of the fields $\varphi vb$ on vertices *v*_{b} neighboring *v*_{a} so that $\omega [\varphi v]\u2243cos(\varphi v\u2212\varphi \u0304v)z$. That leads to

where

restricts *q* to the only possible charges *q*_{n} = (2*n* − *z*) for *n* = −*z*, −*z* + 1, …, *z* − 1, *z*, each with proper multiplicity $zn$.

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

Because the correlation length is much larger than the lattice constant at low *T*, it is physically intuitive to take $i\varphi \u0304v=\u27e8i\varphi v\u27e9=Vve$, the entropic field, which is real. Then, the *A* ↔ *B* symmetry and the previous equations imply that the probability of a charge distribution *q*_{v} can be factored as

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

In Eq. (A6) correlations are transmitted among vertices by the entropic field $Vve$. *Z*_{v} normalizes *ρ*_{v}(*q*) and depends on *v* via $Vve$, which in turn depends on the *collective* charge distribution and is therefore non-local, or $Vve=Vve[q]$. Finally, while Eq. (A6) was deduced here from a mean field approximation of the exact equation (23), it can, however, also be introduced directly in the manner of *Landau*, so to speak, as we have shown previously.^{98}

What we are missing now is 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 as

## REFERENCES

We call it current because a magnetization *⃗M* generates an electrical current density *⃗j* = *⃗*∇ ∧ *⃗M*.

We use dyadics: if $v\u20d7$ is a vector, the components of $v\u20d7$$v\u20d7$ are [$v\u20d7$$v\u20d7$]_{ij} ≔ $vivj$; $1$ is the unitary matrix.

When monopoles interact via a 1/*r* law, *ξ*_{‖} is no longer a correlation length, as the monopole–monopole interaction destroys the screening at least in principle, as we have shown elsewhere.^{49}

$\gamma \u02c6$ is the unit vector of $\gamma \u20d7$.

_{2}Ti

_{2}O

_{7}

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 predicte,^{49}due to the interplay of the screening length and the Bjerrum length.

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.

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.

Depending on the system, the aging might be slow rather than exponential,^{94,95}in which case tiny currents might be detected after a long time.