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

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.

Degenerate artificial square ice is perhaps the simplest two-dimensional system in which to study disorder constrained by the ice rule 1 , and its violations as monopole excitations 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 icessystems 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][15][16][17][18][19][20][21] , and revealing novel phenomena absent in its crystal analogue 12,22 , recent fabrication and characterization advances [23][24][25] 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 8 the first to be realized artificially 9 , when it was shown that nonice rule 26 vertices are suppressed after AC demagnetization 9,27,28 . 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 class 29,30 toward an ordered antiferromagnetic ground state [31][32][33][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 been also demonstrated in rectangular lattices 40,41 , or via cleverly placed interaction modifiers placed in the vertices 42 , or nano-holes in a connected spin ice of nanowires 43 , or by rotation of the moments 44 . Recently, a pure square ice, i.e. a 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, which lead to different spectral characterizations 46 .
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 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 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 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 case 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, 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 2,50-60 . We also particularized some results that what we had already found on generic graphs 61 .

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 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 in terms of longitudinal and transverse potentials h || , h ⊥ , whereê 3 =ê 1 ∧ê 2 ,ê 1 ,ê 2 is an orthonormal basis of the plane, and we call S || , S ⊥ the longitudinal and perpendicular components of the field. Unlike in the 3D case, where the perpendicular part of 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 If for a vector w we call the perpendicular of w, we have then which expresses the duality between charges and currents, or longitudinal and perpendicular components of the field, under a π/2 rotation. This duality has a discretized analogue in square spin ice.

B. Square Spin Ice
Square spin ice (Fig. 1) is a set of classical, binary spins 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 9 as t-I, . . . , t-IV, where t-I and t-II obey the ice rule 26,35,37,48 (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 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 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 t-I (2), Q=0 t-II (4), Q=0 t-III (6),Q=±2 t-IV (2), Q=±4 pointing clockwise around the edge of the plaquette minus those pointing counterclockwise 64 . For a spin configuration S, consider its perpendicular configuration ⊥ S (Fig. 1), for which old plaquettes are now vertices and old vertices are now plaquettes. We have then as in Eqs. (5).
For each configuration of spins that obeys the ice rule configuration, a unique (up to a constant) height function h ⊥ can be defined on the plaquettes such that is true (pp 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 is "irrotational": the line-sum of spins ⊥ 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 ∀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 is the discrete gradient of h ⊥ , 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 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) of "flat" (ordered) phase, a jargon derived from the theory of the 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 though 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 a positive uniform susceptibility (see later). (Clearly, the same is true for the current-free manifold by replacing h ⊥ with h || .) Eq. (9) can be understood in terms of the zero temperature partition function where 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 exactly 48,69 (see also ref 57 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 From that and the continuum limit of Eq (7) (or S = −ê 3 ∧ ∇h ⊥ ) 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 ⊥ ∇ =ê 3 ∧ ∇) 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 spins 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,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. 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 S as with the constraint ∇ · S = 0. Equation (18), unlike our previous formulas, is also valid in 3D. There, it has been appreciated as Jaccard entropy in water ice 60,[71][72][73][74] , and later in pyrochlore spin ice 50,51,53,55 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 χ 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 the Gauss's Law 45 . 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.

III. ENERGY AND STATES
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): 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 of Fig. 1 (black arrows). Equivalently, by gauge-free duality, if = 0 and κ > 0, the ground state is an extensively degenerate ice manifold for ⊥ S, or the grey 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, > 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 currentfree and charge-free, we have ∆h || = ∆h ⊥ = 0 which implies a uniform S.) This case has not been investigated experimentally, though 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, > 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 and thus promotes topological currents and an ordered ground state.
When κ = , we have J 1 = 0 and the set of vertical and horizontal arrows become two decoupled systems. When κ = > 0 each subsystem is ferromagnetic and the ground state is the four-fold fully polarized state. When κ = < 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 κ = ).
We will not study here the full monopole-currents model of Eq (19), which leads to a rich phase diagram in the β × βκ plane (where β = 1/T and T is temperature measured in units of energy) to be compared with other models 29,76 . We will consider only > 0 and |κ|/ 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, 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 and then sum over the spins, obtaining is a generalized density of states for q v , i p , given bỹ and is therefore the functional Fourier transform of (25) (The product runs on all the edges e =vv once, and 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 generalized free energy Note that F[φ, ψ] can have imaginary values. As a consequence, though the variables φ, ψ are themselves real, their expectation values φ , ψ are imaginary, and thus the entropic fields V e q , V e i 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φ, iψ correlate charges and currents, the vector fields −iT ∇ vv φ, −iT ∇ pp ψ correlate spins that would otherwise be trivially paramagnetic.
In the spirit of field theory, we could therefore say the following: V e q is an entropic potential acting on charges and V e i on currents, and correspondingly there is an entropic field B e 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 Section VI C for κ ≤ 0) integrating over charges and currents in Eq. (23) returns Gaussians in the fields, of variance φ 2 = /T , ψ 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, 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, Then, from q = − ∇ · S, j =ê 3 · ∇ ∧ S, and Eqs. (29) we find screened Poisson equations for q, i (here νq ext = ∇ · H, where q ext is the external charge, νi ext = ∇ ∧ H, where i ext is the external current, ν is an energy) and therefore "correlation lengths" 77 as This expression for ξ || was already appreciated 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 VI C. Instead when κ = 0 the integral over i returns a functional delta function on ψ, and ψ disappears from the equations, which for all purpose is equivalent to taking ξ ⊥ = 0.

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 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 where is an effective Hamiltonian for the new variables and the vector γ has components for α = x, y, and contains informations on the square symmetry of the lattice [in the long wavelength limit: 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 so-called magnetic fragmentation [78][79][80] , that is the decoupling between charge-full and current-full ensembles.
The third line contains the coupling with the external field H (we neglect here V ). Because −i γ · w, −i γ ∧ w are the generalized divergence and curl respectively on the lattice in moment space for a generic field w, the entropic fields for currents (resp. charges) couple to the curl (resp. divergence) of the external field.
The quantity χ 0 , already encountered in the previous section on the opposite limit (T = 0) is the uniform susceptibility. 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φ,ψ whenH = 0 returns the effective free energies for q and i in 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 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 γ 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 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.
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χ || (k),χ ⊥ (k) are the the longitudinal and perpendicular static susceptibilities (multiplied by T , thus accounting already for the Curie-Weiss law) for S || = S ·γ (the longitudinal, charge-full, current-fee part of the magnetization) andS ⊥ = S · ⊥γ (the perpendicular, charge-free, current-full magnetization), respectively 81 . Indeed, by integrating the Gaussian integral in Eq. (36), and remembering thatγ αγα + ⊥γ α ⊥γ α = δ αα , one obtains the total free energy as a function of the external field at lowest order as From it, the spin correlations are In the limit T → ∞, Eq. (46) correctly returns S α ( k)S α ( k) → δ αα , 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χ ⊥ . 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 (51) 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 ∇ · S = 0, thus returning Eq. (18), previously found heuristically.
Equation (51), when ξ ⊥ = 0 reduces to the functional found elegantly via methods of chemical physics for pyrochlore ice by Bramwell 58,59 . Indeed, when expressed in terms of spins S rather than currents and charges i, q, the two formalisms are expected to coincide.
Expressing the coarse grained 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: When ξ ⊥ = 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 → 0. For T > 0 monopoles appear, and thus the longitudinal height function is also present.

B. Low Temperature
When T ↓ 0, fluctuations of the entropic should fields diverge. One could then perform a proper study of the higher orders expansion of Ω(φ, ψ) in terms of Feynman diagrams. Because of the complex shape of Ω(φ, ψ) 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 82 . Secondly, 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-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 ↓ 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, 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 which implies a redefinition of constants in the theory at low T : , κ are dressed as 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 49 that Eq. (53) can follow from a Debye-Hückel approach to the entropic potentials (see also Appendix I).
Note finally that χ 0 must be finite at low T (see also below), and that q 2 is decently approximated by the naive q 2 computed by assuming uncorrelated vertices, each with proper multiplicity and Boltzmann weight. At low T , q 2 ∼ (16/3) exp(−2 /T ). We have thus ξ || ∼ exp T 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 pryrochlore 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 the previous sections, if expressed in terms of ξ || , ξ ⊥ , where and similarly and for χ 0 , χ 0 ∼ 1 for T ↑ ∞ while it remains finite at T = 0 (We show below that e.g for T = 0, κ = 0 we have χ 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 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 section II. In the purest form-i.e. without the interference of longrange 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 we have described elsewhere 49 .
When κ = 0, the functional integration over p di p in Eq. (23) produces the delta functions p δ(ψ p ). Further integration over p dψ p returns for the partition function whereΩ[q] is the density of states for the charges, and it is given byΩ and is therefore the functional Fourier transform of (59) In other words, the currents disappear from the picture, and the entropic potential −iT ψ 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 which for small wave vectors reduces to Eq. (16).

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. E.g., from Eq. (43),χ ⊥ ( k) = χ 0 , and the perpendicular susceptibility in real space is a delta function, consistent with zero correlation length.
Considering charge correlations and screening, Eq. (42) can be rewritten as where the first term Fourier-transforms to a Kronecker delta, whereas the second term implies at large distance the charge correlation Note that the modified Bessel function K 0 is the screened 2D-Coulomb potential. It is exponentially screened, or K 0 (x) ∼ π/2x exp(−x) making ξ || the correlation/screening length. In 3D, we would have a screened 3D-Coulomb or 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 a screened 2D-Coulomb, given by the 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,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 γ( k + K) 2 = 8 − k 2 . This leads to a screening function where ξ 2 p = ξ 2 || /(1 + 8ξ 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 Such sign alternation with the Manhattan distance for screening at short distance was also verified in experiments 45 on a quantum annealer 86 .

Spin Correlations
Considering now the spin correlations, Eq. (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, 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 ξ −1 || = 0. Then, for the pure ice manifold we have the familiar transverse form That is because, the longitudinal susceptibilityχ || 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 In Fig. 4 we plot the structure factor Σ m ( k) 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 ) in the case of a line state, as well as from simulations 46,87 . 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 χ 0 . Because spins have values S α = ±1, (2π) 2 , and summing over α = x, y, from Eq. (68) we obtain Proceeding in the same way, but using the second line of Eq. (66) we obtain for generic temperature 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 χ 0 → 2 for T ↓ 0 and χ 0 → 1 fo T ↑ +∞.

B. Weak Ferromagnetic Case (κ 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 κ/ (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 at low temperature square ice with heigh offset exceeding the critical offset 35,38 , 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 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 46 and Fig 2 in ref 45 ). 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 45 ) in the case of a line state, as well as from simulations in refs 46,87 .
C. Weak Antiferromagnetic Case (κ 0 − ) When k < 0, H 2 in Eq. (39) is not bounded from below when T ≤ T afm c = 8|κ| (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 ordered state. Because we have employed a mean field theory, T afm c is not the actual critical temperature, but merely a useful parameter in the context of our framework.
For T > T afm c ,χ ⊥ (k) has a maximum on the K corners of the BZ. Thus, we expand around K, γ(K + k) 2  inχ ⊥ (k) and from Eq (42) we obtain for large |p| which expectedly alternates sign on adiacent plaquettes. The AFM correlation length is given by 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 transition 29,30 , 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 κ/ is small, T afm is much smaller than the crossover temperature for the ice regime (typically ∼ 2 ). Thus, above T afm c the divergence-free and divergencefull 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 45 ) in the case of a line state, as well as from simulations in refs 46,87 .

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 88 . We anticipate here the main deductions, leaving a more in-depth study to future work on "slow" electrodynamics of spin ice.

A. 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 2 , simple considerations should readily convince that the flow density vector 89 for the charges, in or out of equilibrium, is always 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 chargeq which is merely the derivative of Eq. (2). Similarly, taking instead the curl, we havė with 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 90 .

B. Relaxation Dynamics
As mentioned, a simple relaxation dynamics at equilibrium would not apply to real materials, as it is clear both experimentally [91][92][93] and conceptually [94][95][96] . Depending on the system, however, deviations can be limited to kinetics much faster than relaxation 96 .
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 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.

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.
Thus, the longitudinal and transverse relaxation times are and thus the relaxation frequencies (from ν = 1/τ ) obey the dispersion relations where D || , D ⊥ are diffusion constants (see below) given by [Also, from Eq. (83, 81) one finds Eq. (7.5) of ref 58 .] From Eq. (43) we see that when > 0, κ > 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 κ = 0 the transverse relaxation time τ ⊥ ( k) is flat and equal to τ 0 because processes concerning currents do not change the energy.
From Eq. (80) we find the dynamic susceptibilities (defined for dimensional convenience asS =χβH) arẽ as already appreciated via other methods in pyrochlores spin ice 58 . The dynamic susceptibilities in Eq. (84) are typical of an exponential relaxation to equilibrium 97 . 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/ω 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 98 , further underscoring the expected limits at frequencies of a mean field relaxation kinetics 96 .

Real Space Picture
Equations (85) become, in real space at length scales larger than the lattice constant, diffusion-relaxation equationsq The first term in Eqs. (87) accounts for the diffusion of charges or currents. The diffusion flow density vectors are J q,diff = −D || ∇q 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 q = 0. In presence of external charges the system relaxes to the screening Equations (34).
A similar diffusion-relaxation equation can be written for S. From Eq. (80) we havė which can also be written aṡ in terms of the flow density tensor for magnetization However magnetization is not conserved, but because of the second term −S α , it decays exponentially to the equilibrium value S α = 0 imposed by the Z 2 symmetry of the problem (if H = 0). Equation (89) can also be written in terms of charges and currentṡ 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 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 2,58 . At equilibrium, for a constant uniform field 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 κ = 0, D ⊥ = 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.

A. 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, Ω[φ, ψ] 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 ∇φ, ∇ψ 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 Ω 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 κ = 0 (and thus integration over currents enforces ψ = 0) and H = 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(φ va − φ v b ) cos(φ va − φ va ), where φ va is the mean of the fields φ v b on vertices v b neighboring v a , so that whereω (q) = 2 −z z n=0 z n δ(q − q n ) (97) restricts q to the only possible charges q n = (2n − z) for n = −z, −z +1, . . . , z −1, z, each with proper multiplicity z n .
We have thus obtained a mean field approach that, unlike the high T approach, separates in Ω[φ] the effect of the entropic interaction, now expressed via the mean fieldφ, 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φ v = iφ v = V e v , 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. (99) correlations are transmitted among vertices by the entropic field V e v . Z v normalizes ρ 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 99 .
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 where q v = dqρ 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 ξ || replaced by Eq. (53) (see also ref 49 ).