Spin ice thin films: Large-N theory and Monte Carlo simulations

We explore the physics of highly frustrated magnets in confined geometries, focusing on the Coulomb phase of pyrochlore spin ices. As a specific example, we investigate thin films of nearest-neighbor spin ice, using a combination of analytic large-N techniques and Monte Carlo simulations. In the simplest film geometry, with surfaces perpendicular to the [001] crystallographic direction, we observe pinch points in the spin-spin correlations characteristic of a two-dimensional Coulomb phase. We then consider the consequences of crystal symmetry breaking on the surfaces of the film through the inclusion of orphan bonds. We find that when these bonds are ferromagnetic, the Coulomb phase is destroyed by the presence of fluctuating surface magnetic charges, leading to a classical Z_2 spin liquid. Building on this understanding, we discuss other film geometries with surfaces perpendicular to the [110] or the [111] direction. We generically predict the appearance of surface magnetic charges and discuss their implications for the physics of such films, including the possibility of an unusual Z_3 classical spin liquid. Finally, we comment on open questions and promising avenues for future research.


I. INTRODUCTION
The emergence of gauge structures in strongly correlated systems has proven to be an essential thread in the fabric of modern condensed matter physics [1][2][3][4]. In the prototypical example of a gauge theory -electromagnetism -boundary conditions can play a key role in the physics [5]. Indeed, realizations of gauge theories in systems with confined geometries can lead to rich and varied phenomena as, for example, in the Casimir effect [6,7]. In the same spirit, questions pertaining to surface effects in emergent gauge theories of strongly correlated systems, have only recently begun to be addressed [8]. A paradigmatic example where such an emergent U(1) gauge theory arises is in spin ice materials [9,10], a class of highly frustrated three-dimensional magnets. Given the level of maturity of research on spin ice [11,12], with many theoretical successes and several well-understood experimental examples, it is a natural system to explore the effects of confined geometries in emergent gauge theories.
In the prototypical spin ice materials Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 , the magnetic moments reside on the sites of a pyrochlore lattice, which is formed of corner-sharing tetrahedra, as shown in Fig. 1. At low temperatures, these magnetic moments are forced by the crystalline electric field [13] to point either in or out of any given tetrahedron. The strongly frustrated interactions between the magnetic moments then give rise to a local "2-in/2-out" constraint on every tetrahedron in the ground state, a close of analogue of the arrangement of protons in common water ice [11,12]. This constraint can be rewritten in a form similar to Gauss' law in electromagnetism [14], giving rise to an emergent Coulomb phase [15,16] in these materials. This phase is characterized by algebraic spin correlations with fractionalized excitations taking the form of emergent magnetic monopoles [16,17]. Furthermore, quantum models of spin ice materials have been suggested as promising platforms to realize a related quantum spin liquid phase [18]. Research on spin ice materials has so far mainly focused on bulk properties [12]. In water ice, some interesting physics has been found by looking at the effects of confined geometries. For example, by squeezing water between two sheets of graphene [19], one finds that the water molecules form a square lattice, reminiscent of the six-vertex model. The investigation of spin ice films in such confined geometries -thin films, in particular -is now developing [20][21][22]. Very recently, the first films of Dy 2 Ti 2 O 7 [20] and Ho 2 Ti 2 O 7 [22] were grown, with Ref. [20] reporting a vanishing residual entropy at low temperature, in strong contrast with bulk physics [23]. The theoretical work to date has tackled a variety of issues; for example, Refs. [24,25] considered heterostructures involving spin ice materials, while Jaubert et al. [26] investigated the dipolar spin ice model in a thin film geometry using Monte Carlo methods. However, to the best of our knowledge, there currently is no theoretical understanding of even the simplest minimal model, that of nearest-neighbor spin ice films.
In this paper, we explore the physics of nearest-neighbor spin ice films, considering the fate of the three-dimensional Coulomb phase as well as the effects of different surface terminations. We take a two-pronged approach: we use the analytical large-N method, which has been successful in applications to bulk spin ice [27] and films of ferromagnets [28], and validate its predictions for nearest-neighbor spin ice films using Monte Carlo simulations. We focus our investigation on the simplest highly symmetric film geometry, with surfaces perpendicular to the crystallographic [001] direction. We find that: (i) The characteristic pinch points found in the spin-spin correlations [14,15] of bulk spin ice remain intact for momenta parallel to the surfaces, a signature of a two-dimensional Coulomb phase (a classical U(1) spin liquid). (ii) The direct space spin-spin correlations oscillate as a function of depth in the sample, with an amplitude that increases with decreasing temperature. (iii) By including orphan bonds to capture some of the crystal symmetry breaking of the film surfaces, we find that the Coulomb phase and its associated pinch points disappear when the exchange on the orphan bonds is ferromagnetic, yielding a classical Z 2 spin liquid [29]. These results are summarized in the phase diagram shown in Fig. 2. Finally, building on this understanding, we extend these results to discuss the surface states of films with cleaved surfaces perpendicular to the [110] or [111] direction. From general considerations, we predict the appearance of surface magnetic charges in the ground state, akin to the monopoles realized as excitations in bulk spin ice. We discuss some implications of these surface charges, offering guidance for future studies on spin ice films with such geometries.
The rest of the paper is organized as follows: in Sec. II, we detail our model and then, in Sec. III, develop the large-N formalism used to investigate spin ice films. In Sec. IV, we apply the large-N method to films with surfaces perpendicular to the [001] direction, and compare these results to those obtained from Monte Carlo simulations. In Sec. V, we discuss the topological order of the classical U(1) and Z 2 spin liquids found in these films. Sec. VI briefly addresses other cleaving geometries, while Sec. VII offers concluding remarks and comments on possible avenues for future research. In Apps. A and B, we provide details of the large-N theory for bulk spin ice and [001] films. In App. C, we discuss the numerical solution of the large-N saddle point equations. Finally, in App. D, we provide details of the Monte Carlo algorithm used to simulate Ising (N = 1) spin ice films.

II. MODEL
A. Nearest-neighbor spin ice model To set the stage, we first review the essential features of the nearest-neighbor (nearest neighbor (NN)) spin ice model and then, in Sec. II B, we minimally extend it to the context of

films.
Recall that in classical spin ice [13], the magnetic moments are represented by pseudospins S i = σ iẑi living on pyrochlore lattice sites labeled by i, where σ i = ±1 is a classical Ising variable andẑ i are unit vectors along the local quantization axes (see App. A). In this work, we consider the simplest spin ice model, which only takes into account NN Ising exchange interactions. In the bulk this is the celebrated pyrochlore Ising antiferromagnet model [30] where J > 0, and the sum runs over NN bonds of the pyrochlore lattice. This Hamiltonian has a degenerate ground state manifold where every tetrahedron respects the local ice rules, i.e. the sum of Ising spins on every tetrahedron is zero. This realizes a classical U(1) spin liquid, with an extensive ground-state degeneracy that remains down to zero temperature, thus giving a nonzero residual entropy. The structure of the ground-state manifold can be formulated in terms of a coarse-grained effective "magnetic" field B(r) defined on each tetrahedron as B(r) ≡ (−1) r i∈r σ iẑi , where the sign, (−1) r , depends on the sublattice of the tetrahedron in the dual diamond lattice (see Ref. [15] for a review). In terms of B, the ice-rule constraint amounts to a divergence-free condition, ∇ · B = 0. Excitations above the spin ice ground-state manifold appear as pointlike sources or sinks of the field B, behaving effectively as deconfined magnetic monopoles [17]. At low temperatures, an analogue of classical magnetostatics thus emerges and induces a cooperative paramagnetic state dubbed a "Coulomb phase" [15]. The divergence-free constraint also implies dipolar spin-spin correlations which manifest themselves as sharp "pinch points" in reciprocal space [14,15].
Even though it appears greatly simplified compared to dipolar spin-ice (dipolar spin ice (DSI)) materials [31] such as Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 , the NN model captures much of the essential physics of the Coulomb phase shared with more realistic DSI models [17,[31][32][33]. Although they are the best examples, dipolar interactions are not the only route to realizing spin ice. Rare-earth magnets where super-exchange is dominant could potentially host more faithful realizations of NN spin-ice [Eq. (1)] due to the short-range nature of the exchange physics. For example, the Pr 2 M 2 O 7 family [34][35][36], recently discussed as quantum spin-ice [18] candidates, are expected to have NN Ising exchange that is significantly larger than the magnetostatic dipolar interactions [37,38].

B. Film geometries and boundary conditions
In order to model NN spin ice films, one must first define the boundary conditions, such as choosing a cleaving plane along which to cut the pyrochlore lattice, exposing free surfaces to a putative vacuum. For simplicity, we consider a free standing film and ignore complications arising from the presence of a substrate [20,22] Fig. 3. For simplicity, we consider thicknesses corresponding to an integer number L of conventional cubic unit cells, comprising 4L spin layers (which we label by l) where the top chains are along [110] and the bottom chains run along [110]. The primitive unit cell of the film thus comprises 4L layers and 8L spins (see App. B). The associated conventional unit cells for film thicknesses L = 1, 3, 5 are shown in Fig. 3.
As noted in Ref. [26], this cleaving renders some of the surface bonds locally inequivalent to those in the bulk. Generically, one expects the bonds that join the remaining spins of a cut-off tetrahedra, which we call orphan bonds following Ref. [26], to have an Ising coupling, J O , different than the other bonds in the film. We thus consider the following minimal model [110] [110] . . . [110] [110] [110] . . .  where the first sum, i j , runs over all NN bonds while the second sum, i j ∈ O, runs only over the orphan bonds [39].

III. METHODS
With our model of spin ice thin films defined, we now outline the methods used to tackle these systems. We first discuss the large-N method and review its application to bulk spin ice. Next, we discuss the modifications needed for an application to spin ice thin films. Finally, we discuss the Monte Carlo methods used to simulate Ising (N = 1) spin ice films directly.

A. Large-N method in bulk spin ice
The Hamiltonian for NN spin ice, Eq. (1), can be investigated using an analytically tractable approximation scheme, the so-called large-N expansion. This method allows one to obtain semi-quantitative spin-spin correlation functions at low temperature [27]. Consider the classical partition function, Z, for the Hamiltonian in Eq. (1), written as where the sum over i and j runs over all pyrochlore lattice sites, and we have defined V i j = 1 when i and j are nearestneighbors, and V i j = 0 otherwise. Replacing the classical Ising spins σ i by continuous variables s i and enforcing the unit spin length constraint leads to the partition function In this form, the length constraints render the partition function as intractable as the original model. The large-N approach circumvents this problem by extending these real variables, s i , to N-component vectors s i subject to the constraint The interaction between the spins is also extended to be O(N) symmetric, with the resulting Z N partition function now taking the form Clearly, if we set N = 1, we recover the original Ising model with Z 1 ≡ Z. The constraints at each site i can be enforced using constraint fields µ i , so Z N becomes where Ds ≡ j ds j and Dµ ≡ j dµ j . Integrating out the s fields yields where we have defined the diagonal matrix µ with elements µ i j ≡ δ i j µ i . In this form, it is clear that as N → ∞, the partition function is dominated by the saddle points of the exponential.
In the saddle-point solutions, µ is purely imaginary, so we consider the real variable λ ≡ iµ. The saddle-point equations are then given by The correlation functions between the s i can be readily obtained from Eqs. (6,7) by taking a derivative with respect to V i j . One obtains where a, b = 1, · · · N label the spin components. Note that by invoking this correlation function, the saddle-point condition can be interpreted as an average length constraint on the spins s i with For bulk spin ice, the translation and rotation symmetries of the lattice enforce that the λ i are site independent with λ i ≡ λ 0 . We can then block diagonalize the interaction matrix V using a Fourier transform which leads to where α, β = 1, . . . , 4 index the four sublattices of the pyrochlore lattice, and the explicit form of the matrix V(q) is given in App. A. The average length constraint [Eq. (10)] can then be expressed as where n is the total number of spins in the system. The previous derivation, leading to Eqs. (11) and (12), is equivalent (in outcome) to the so-called self-consistent Gaussian approximation or spherical approximation [40]. Perhaps surprisingly, this large-N treatment has been found to provide semi-quantitative correlation functions when compared to Monte Carlo simulations of pyrochlore Ising (N = 1) and Heisenberg (N = 3) antiferromagnets [27]. However, this method does not provide a good description of the physics for N = 2 due to the manifestation of order-by-disorder [41][42][43]. In principle, a 1/N expansion [27] around this exactly solvable point would allow one to obtain more precise correlation functions for spins with finite N, but given the success of the N → ∞ results, this seems unnecessary. The method has also been extended to include features found in more realistic models of spin ice; this includes further-neighbor exchange interactions [44,45] and dipolar interactions [46].

B. Large-N method for spin-ice films
The application of the large-N method to spin ice models in film geometries introduces additional complications to the methodology. In particular, the breaking of the translational symmetry in the finite direction of the film forbids a uniform constraint field, λ 0 (as is the case in bulk spin ice). However, one can still take advantage of the translational symmetry in the plane parallel to the surfaces, defining a constraint field on each layer l. Such layer-dependent constraint fields were used previously in Ref. [28] to study Casimir effects in films of ferromagnets.
Proceeding as in the previous section, and integrating out the spin variables s i , we obtain the large-N partition function where the matrix λ now has layer-resolved elements, λ i j = λ l δ i j , and the layer index l is an implicit function of the lattice site i. Using the translational symmetry in the plane, the saddlepoint solution [Eq. (10)] can be defined on each layer, where n l is the number of spins on layer l. Spin-spin correlations in reciprocal space are obtained from Eq. (9) after performing a Fourier transform in the plane, where r runs over the n c primitive unit cells of the film, q ⊥ are in-plane wave vectors, and r α are basis vectors locating each sublattice α within the unit cell. Note that we identified the pyrochlore lattice sites as i ≡ (r, α). One ultimately finds where The numerical values of the constraint fields λ l [47] are obtained by enforcing the saddle-point conditions, given in Eq. (14), which can be expressed as for each layer l. We note that the framework provided by Eqs. (16) and (17) is completely general and does not suppose a particular choice of surface geometry, which appears in the definition of the unit cell and through the structure of the matrix M(q ⊥ ). We also note that this analysis would carry through for spin ice films that include further-neighbor or dipolar interactions [26], in their paramagnetic phases. One simply needs to compute the Fourier transform V(q ⊥ ) of the corresponding interaction matrix in the chosen film geometry (see Refs. [44][45][46] for details in the bulk case). The inclusion of orphan bonds, as in Eq. (2), is also straightforward, with the corresponding interaction matrix V(q) given in App. B.

C. Monte Carlo simulations
To confirm that the large-N method correctly captures the physical behavior of the Ising (N = 1) films at low temperatures, as it does in the bulk case [27], we perform a classical Monte Carlo simulation of the model of Eq. (2) for the appropriate film geometries. To avoid issues with equilibriation, we use a non-local Monte Carlo update. Specifically, we adapt the cluster algorithm of Ref. [48] to the film geometry. To implement the surfaces in the [001] direction we consider a periodic system of cubic cells with dimensions L ⊥ × L ⊥ × L. This can be modified into the appropriate film geometry by cutting the bonds that pass through a plane with normalẑ and changing the remaining two bonds to carry the orphan coupling, J O . This modifies the weights used for the surface tetrahedra in the cluster algorithm of Ref. [48], but otherwise leaves the algorithm unaffected for any choice of J O /J (see App. D for further details). This cluster algorithm is closely related to the standard loop or worm algorithm used in spin ice simulations [49,50], similar to the relationship between the Swensden-Wang [51] and Wolff [52] cluster algorithms used in unfrustrated Ising models. For our purposes, one advantage of this formulation is the availability of an improved estimator [48] for the spin-spin correlation functions that allows to access larger system sizes at lower computational cost. Typically, accurate spin-spin correlations can be obtained with samples generated using only 10 3 steps of the cluster algorithm when employing this improved estimator.
For the single-layer films (L = 1), we considered systems up to L ⊥ = 64, while for the thicker films (L = 3 and L = 5) we considered sizes up to L ⊥ = 32. For a given system size n = 16L 2 ⊥ L, we expect finite size effects to become important when the monopole density ∼ e −2J/T [53] drops below ∼ 1/n. Below the crossover temperature T * /J ∼ 1/ log n one expects the system to be confined to the ice manifold itself and recover the T = 0 behavior. For the lattice sizes of interest, the simulations become finite-size limited for temperatures less than T * /J ∼ 0.1 − 0.2. When considering the direct-space spin-spin correlators, we used smaller sizes of L ⊥ = 16, but with a larger number of samples, typically of order 10 5 , to ensure small statistical errors.

IV. [001] SPIN ICE FILMS
To begin our exploration of spin ice films, we consider what is perhaps the simplest geometry: films with surfaces perpendicular to the [001] crystallographic direction. We start with equivalent orphan and bulk bonds, J O = J, and then consider the more general and richer case, J O J. In both cases, we apply the large-N method described in Sec. III B, and confirm its results via Monte Carlo simulations as described in Sec. III C.

A. Equivalent orphan and bulk bonds
We first consider films with the orphan bond coupling equal to that of the bulk (J O = J). We start with the thinnest films (L = 1), where we expect the most pronounced effects when compared to the bulk case. To proceed, Eq. (17) must be solved numerically in order to obtain the set of constraint fields λ l as a function of T/J, which is needed to access all other observables. This is accomplished by applying the Newton-Raphson descent algorithm (see App. C). Although there are four spin layers, the symmetry of the slab leads to only two distinct constraint fields, as shown in Fig. 4(d). As T/J → 0, the constraint fields for the middle layers approach the expected value for bulk spin ice, λ 0 = 1/2 [27], whereas the constraint fields for the surface layers converge to a significantly lower value. As a result, the in-plane spin-spin correlations acquire a layer-resolved character, with stronger correlations at the free surfaces than in the middle of the slab. As an example of this behavior, we show in Fig. 4(a) the direct space correlation  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 4 5 6 7 8 9 10 11 12 13 14 15 16 17  To explore the fate of the key signature of the Coulomb phase, the presence of pinch points [14][15][16], we consider the spin-spin correlation function, S (q), of the films in reciprocal space where q ≡ (q ⊥ , q z ) = 2π(hx + kŷ + lẑ) is a three-dimensional wave vector, expressed using Miller indices [hkl]. In terms of the spin-spin correlation matrix M(q ⊥ ), we obtain In Fig. 5, we show S (q) for various temperatures T/J in two high-symmetry planes: [hk0] and [hhl]. Strikingly, the characteristic pinch-points remain intact in the [hk0] scattering plane, parallel to the surfaces. However, as expected due to the finite extent in theẑ direction, they are washed out in scattering planes with a non-zero q z component, normal to the film. We also observe "scattering rods" in the q z direction, and "secondary" pinch-points near [110] and equivalent wavevectors.
All these features are reproduced in Monte Carlo simulations for the same geometry. The only substantive difference between the Monte Carlo and large-N results lies in the temperature dependence of the build up of spin-ice correlations, as found in the bulk case [27]. For the large-N case, due to the continuous nature of the spins, the correlation functions only approach the asymptotic T = 0 result algebraically [14], while for discrete Ising spins this occurs exponentially [54]. This can be seen explicitly in the pinch points; as ∼ √ T/J → 0, their width decays as ∼ √ T/J for the large-N case, as com- We now examine how the properties uncovered above change as the thickness is increased towards the bulk limit. Surprisingly, the numerical solution for the constraint fields λ l  shows oscillations as a function of depth in the sample, as is illustrated in Fig. 4(e,f) for films with L = 3 and L = 5. These oscillations have increasing amplitude with decreasing temperature and a characteristic damping length scale which appears independent of thickness, indicative of a surface effect. They are also seen in the layer-resolved direct-space correlations C l and are well reproduced in the Monte Carlo simulations, as shown in Fig. 4(b,c), keeping in mind the difference in temperature dependence expected between the large-N theory and the Ising model. We have verified that a large-N treatment of thin films of a pyrochlore Ising ferromagnet [i.e. Eq. (1) with J → −J] with the same geometry does not show oscillations in the constraint fields λ l , but rather a monotonic behavior from the surface to the middle of the sample, as long as the system remains in the high-temperature paramagnetic phase. These results thus suggest that the oscillations are directly related to the geometrical frustration.
We also compute the spin-spin correlation functions S (q) [given by Eq. (19)] at T/J = 0.1, deep in the Coulomb phase, for thicknesses of L = 1, 3 and 5 (see Fig. 6). In the [hk0] plane, pinch-points are always present, with an increased contrast for thinner films. In the [hhl] plane, the washed-out pinch-points are progressively restored with increasing thickness, as the system crosses over from a two-dimensional to three-dimensional Coulomb phase. The restoration of the "three-dimensional" pinch-points in the [hhl] plane is set by the thickness of the film L which cuts off the Coulomb correlations in theẑ direction. Roughly, we expect quasi-two-dimensional behaviour when the correlation length, ξ ∼ e J/T , becomes of the order of the film thickness, L, giving a crossover temperature T ∼ J/ log L. As with the L = 1 films and bulk results, the Monte Carlo and large-N results agree well, aside from the aforementioned temperature dependence (algebraic compared to exponential) of the build up of Coulomb correlations.

B. Inequivalent orphan and bulk bonds
We now move to the more general case and consider the influence of differing orphan and bulk bonds (J O J), as defined in Eq. (2), on the physics of NN spin ice films. At low temperature (T |J O |, J), we find that the system remains paramagnetic with the physics depending only on the sign of J O (not on its magnitude). This is expected because the "bulk" ice rules are always compatible with minimizing the energy of the orphan bonds. We show the spin-spin correlations function, S (q), for three representative values J O /J = +1, 0, −1 in Fig. 7. The case J O /J = 1 corresponds to the results of the previous section, with sharp pinch points characteristic of a two-dimensional Coulomb phase. However, these pinch points completely disappear for J O /J = 0 and J O /J = −1, revealing only broad features in S (q) in the low-temperature regime.
The preservation or destruction of the Coulomb phase (depending on the sign of J O /J) can be understood in terms of a simple picture of the ground state manifold for arbitrary film thicknesses. Note that for J O /J > 0, the flux lines of the field B (see Sec. II A) run along the surface and are then redirected back into the bulk of the film [see 8(a)]. This choice of orphan bond coupling is thus compatible with the (bulk) constraint ∇ · B = 0 and the Coulomb phase, now two-dimensional, is preserved. However, as discussed in Ref. [26], when J O /J < 0, the orphan bonds host pairs of aligned spins which correspond to surface magnetic charges. These surface charges serve as endpoints to the "flux-lines" of the effective magnetic field B [55] [see Fig. 8(c)]. Just as a finite density of thermally populated monopoles endows the pinch points in bulk spin ice with a finite width [15,53], the finite density of these fluctuating surface charges broaden the pinch points in spin ice films. However, unlike the bulk case, these charges are present even at T = 0 and thus destroy the Coulomb phase. We note that the destruction of the Coulomb phase is ultimately averted in Ref. [26] by the ordering of the surface charges due to the long-range dipolar interactions. This static ordering inhibits thermal fluctuations of the surface charges and restores the twodimensional Coulomb phase below the ordering temperature -a behavior reminiscent of magnetic fragmentation [56]. In other words, the state we find for J O /J < 0 can be viewed as the "parent" state out of which the ordering of Ref. [26] arises.
Finally, we consider the special case J O /J = 0. Since the orphan bonds provide no energetic constraint, their spin configuration is essentially random. This leads to both the presence and absence of surface charges, as illustrated in Fig. 8(b). The finite density of these surface charges, while not maximal (as for J O /J < 0), is sufficient to destroy the Coulomb phase. This results in the absence of pinch points in S (q), as observed in Fig. 7(b).
We have verified that when L is increased, the spin-spin correlation functions S (q) approach the bulk result for any J O /J (not shown). In particular, the gradual recovery of the pinchpoints indicates that algebraic correlations are present up to a length scale set by the film thickness, L, even in the presence of surface charges. As mentioned above, this is analogous to the case of thermally activated charged defects (monopoles) in bulk spin ice which cut off the algebraic correlations at a length scale set by their average separation [53]. We have also investigated the depth dependence of the constraint fields λ l and the layer-resolved correlators C l ; we find that these oscillate as a function of layer index l (for L > 1) for all values of J O /J considered (not shown).

V. CLASSICAL TOPOLOGICAL ORDER IN [001] FILMS
The previous section has exposed, via large-N and Monte Carlo results, how boundary conditions affect the Coulomb phase present in the parent bulk system. In this section, we relate these results to the topological order that characterizes different classical spin liquids [29,57]. In particular, we argue that the low-temperature state of films with J O /J < 0, while not a Coulomb phase, is nonetheless a non-trivial collective paramagnet -a classical Z 2 spin liquid [29].
First recall that in bulk spin ice, the spins can be mapped to dimers on the dual diamond lattice [58]. Specifically, we identify σ = +1 with the presence of a dimer on the corresponding bond of the diamond lattice (similarly, σ = −1 is identified with the absence of such a dimer). The ice rules then   correspond to the requirement that two dimers touch at each diamond lattice site. To move within the ice manifold, one uses "loop moves" that reverse all the spins on an alternating spin loop [49,50]. In the dimer picture, this corresponds to swapping the occupied and unoccupied bonds on this loop. Any "short" loop (not spanning the system) thus preserves three winding numbers, corresponding to the number of dimers crossing three planes oriented in thex,ŷ orẑ directions. These winding numbers are topological invariants that characterize the classical U(1) spin liquid (Coulomb phase), and can only be changed by "large" loops that wind around the periodic directions of the system [59].
How does this physics change for films? First, since our system is two-dimensional, we can define at most two distinct winding numbers. Second, in addition to the usual "bulk" loop moves, we can also construct zero-cost moves that involve the surface spins. When J O /J > 0, the constraint of anti-aligned orphan bond spins only allows the construction of loops that run along the orphan bonds, as shown in Fig. 8(a). The two in-plane winding numbers thus remain topological invariants and we have a two-dimensional Coulomb phase, a classical U(1) spin liquid.
The J O /J < 0 case is more interesting. Here, to construct a zero-cost move, we must consider open strings of alternating spins that end in pairs on the orphan bonds, since preserving the surface constraint requires flipping both orphan bond spins. One can view such a pair of strings as a usual loop, where the alternation pattern is reversed when an orphan bond is encountered [see Fig. 8(c)]. Therefore, contributions from the two strings add up, and these moves can change the winding numbers only by even amounts, leading to the destruction of the aforementioned classical U(1) topological order. However, not all is lost; one can still define two Z 2 topological invariants [29], corresponding to the two "winding parities", that can only be changed by moves that wrap around the system. This leads us to identify the paramagnetic phase found for J O /J < 0 films as a classical Z 2 spin liquid, consistent with the absence of pinch-points exposed in Sec. IV [60]. A direct consequence of such (classical) Z 2 topological order is the predicted presence of fractionalized magnetic moments bound to vacancies in the film [29].
Finally, we turn to the J O /J = 0 case. At this point the system is neither a Z 2 nor a U(1) spin liquid. We can see this noting that strings of alternating spins terminating on the orphan bonds can be flipped at zero energy cost [see Fig. 8(b)]. Flipping these strings can change the winding numbers arbitrarily, even for short strings. We thus conclude that the J O /J = 0 case does not support topological sectors. Given that this case sits at the critical point between the U(1) and Z 2 spin liquids, its properties have more general implications for the regime where |J O | T J. In this limit, the orphan bonds are effectively at high temperature, so the system will behave more like the J O /J = 0 point, rather than the classical U(1) or Z 2 spin liquids (see the phase diagram in Fig. 2).

VI. MAGNETICALLY CHARGED SURFACES IN [110] AND [111] FILMS
In the previous section, we showed how specific boundary conditions (the orphan bond exchange J O ) have, through the formation of fluctuating surface charges, dramatic effects on the properties of the film. With an understanding of the "simple" case of [001] films, we now proceed to briefly discuss more complicated geometries, specifically films with surfaces perpendicular to the [110] and [111] directions. These geometries are obtained by cutting one (or three) spins per surface tetrahedron, as shown in Fig. 9. The resulting slab is comprised of alternating kagome and triangular layers for [111] films, but has a somewhat more complicated geometry for [110] films.
These two geometries differ drastically from the [001] films which, having two spins remaining per surface tetrahedron, can still (in principle) respect the divergence-free condition ∇ · B = 0 defining the Coulomb phase, by having one spin pointing in and one pointing out. Whether this is energetically favorable depends on the value of J O /J, as explained in Sec. IV. This is, however, impossible for [110] or [111] spin ice films, where the surface tetrahedra have either one or three spins remaining -in other words, the boundary conditions imposed on the B field at the surfaces are fundamentally incompatible with the divergence-free condition of the bulk. As a result, the surface tetrahedra must host a charge of the B field, irrespective of the value of J O /J, as illustrated in Fig. 9. In [110] films, these effective magnetic charges sit on parallel "zig-zag" chains running on the surfaces [see Fig. 9(a)], whereas in [111] films, they live on a triangular lattice [see Fig. 9(b)]. As discussed in Sec. IV, when these charges can fluctuate, they serve as free, zero-cost end-points for the (effective) magnetic flux lines. The presence of such zero-cost end points will generically destroy the two-dimensional Coulomb correlations. We thus expect any two-dimensional pinch-points to be destroyed for [111] or [110] films.
However, we do not expect a classical Z 2 spin liquid in these geometries. For J O /J > 0, the surface monopoles are at the end of only one flux line. Thus, there is no constraint on changes of the winding numbers incurred by local updates. In contrast, for J O /J < 0 the surface triangles consist of three aligned spins, and thus three strings must terminate at each triangle. Borrowing the arguments of Sec. V, this would imply the presence of a classical Z 3 spin liquid [61,62] as the winding numbers can only change in multiples of three. Whether the physics discussed above extends to zero temperature, or is preempted by some ordering phenomena requires detailed numerical study, which we leave to future work.
While this picture of surface charges applies to the Ising case (N = 1), the N ≥ 3 cases should be qualitatively different. Consider as an example a [111] film terminating on a kagome layer; while three Ising spins cannot add to zero on orphan triangles, thus requiring surface charges, three N = 3 (or higher) vectors can sum to zero. This implies that such surface charge defects (i.e. a non-zero sum of spins) are not required for the N = 3 (or higher) case. We thus expect that for [110] and [111] films with orphan triangles, the large-N result will be qualitatively different than the Ising case, with the twodimensional Coulomb phase not (necessarily) destroyed [63].

VII. DISCUSSION
We now discuss some potential extensions and implications of our work. In particular, we outline applications to continuous spin systems, possible extensions to dipolar spin ice films and the surfaces of bulk single crystal dipolar spin ices. In addition, we speculate on more theoretical aspects, such as the effects of a magnetic field and the physics of quantum spin ice films [18].
In view of making concrete contact with experimental realizations, we discuss some aspects of films of spin ice materials. There are several complications in connecting the ideas discussed in this work to the physics of these compounds. First, crystal symmetry breaking at the interface between spin ice and vacuum could weaken or even destroy the Ising nature of the spins and their interactions [13] near the surfaces. For canonical spin ices such as Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 , this may not be a serious concern. Since the crystal field ground doublets of Dy-and Ho-based pyrochlores are predominantly maximalrank (mostly J z = ±15/2 and J z = ±8 respectively), strong perturbations to the crystal field would be necessary to generate significant effects on the single-ion or two-ion properties [13]. One might then expect that the induced transverse (quantum) exchange at the surfaces would be small for both compounds and that the splitting of the non-Kramers doublet expected in Ho 2 Ti 2 O 7 might be negligible. For the case of the Pr 2 M 2 O 7 family mentioned in Sec. II, these effects are likely to be more drastic. Indeed, the presence of large random transverse fields [64] due to weak structural disorder appears to be a feature even in bulk samples. Given that the crystal field doublet in these non-Kramers compounds lacks the "axial protection" present in Ho 2 Ti 2 O 7 [13], we expect these transverse fields to be further enhanced at any surfaces. These complications could be minimized through more clever engineering of the films. For example, one might consider heterostructures composed of a thin layer of a spin ice material sandwiched between layers of non-magnetic pyrochlore materials having the same crystal structure and similar lattice constant, such as La 2 M 2 O 7 , Lu 2 M 2 O 7 and Y 2 Ti 2 O 7 .
However, one serious difficulty with all of the proposals for spin ice thin films and heterostructures is the effect of substrate-induced strain. This could be due to a lattice constant mismatch, or simply to slight chemical bonding differences at the interface. The presence of such strain will generically strengthen some of the bonds and remove the ground state degeneracy and residual entropy [26,65]. In non-Kramers compounds, this could also induce a transverse field at each site (depending on the film geometry and strain direction); as in the case of surfaces, this could be negligible for Ho 2 Ti 2 O 7 , but significant in the Pr 2 M 2 O 7 family. How this strain can be minimized, so that the intrinsic physics of spin ice films can be exposed, is an important but exciting challenge in the fabrication of these systems.
In dipolar spin ice such as Dy 2 Ti 2 O 7 or Ho 2 Ti 2 O 7 , the longrange tail of dipolar interactions should have significant effects on the physics discussed here. For bulk spin ice, these differences are suppressed due to the structure of the spin-ice manifold -the dipolar interactions are effectively short-range when acting on ice states [32,33]; this is the so-called self-screening or projective equivalence. Consequently, the splitting of the ice manifold is small, and the associated ordering due to the dipolar interactions only occurs at low temperature compared to the bare scale of the dipolar interactions [49]. However, when monopoles or surface charges are present, the dipolar interaction is significant, promoting the entropic Coulomb interaction between the defects into an energetic one. Indeed, it was found in Ref. [26] that for [001] surfaces with ferromagnetic orphan bonds, this attraction induces a phase transition to long-range order of the surface charges into a checkerboard pattern. For films with [110] or [111] surfaces, the analogous physics is likely to be even richer. For example, in the [111] geometry, the charges live on a triangular lattice, with either a single monopole or anti-monopole at each site. While each surface is not required to be neutral (due the presence of the other surface), this will be favored energetically. One thus expects a one-component Coulomb gas at half-filling (the other component being treated as background) on a triangular lattice. At the nearest neighbor level, this is equivalent to an anti-ferromagnetic triangular lattice Ising model and is highly frustrated, with a macroscopically degenerate set of ground states [66]. The effective long-range Coulomb interactions will presumably lift this degeneracy but, as in dipolar spin ice, only weakly, due to the approximate local charge neutrality.
Some of the physics discussed here is expected to carry over from the thin film context to that of exposed surfaces of bulk crystals of spin ice materials. For example, the presence of fluctuating surface charges in certain geometries may have screening effects on the fields from monopoles in the bulk.
As discussed in Sec. VI, the physics of different surface terminations depends strongly on the nature of the spins in question. While we have mainly discussed Ising spins here, it is known that the large-N method also works well for O(3) spins [27,44]. Examples of pyrochlore anti-ferromagnets with such continuous, classical spins include certain spinels [67], as well as the recently discovered chemically disordered fluorine pyrochlores [68,69]. Thin films of such compounds may be an interesting playground to explore analogues of the physics discussed in this work. Other interesting avenues in this direction include the case of O(2) spins, where the large-N method is known to be unreliable due to the appearance of order-by-disorder [41][42][43]. Experimentally, there are several promising candidates; for example the bulk XY pyrochlores are known to exhibit rather exotic behaviors, from the order-by-disorder physics of Er 2 Ti 2 O 7 [70][71][72] to the unusual physics of the Yb 2 M 2 O 7 family, Yb 2 Ge 2 O 7 in particular [73][74][75] (see Ref. [76] for a review). Further, the limit of an anti-ferromagnetic XY model has been found to harbor several exotic phases, including a spin liquid at intermediate temperature and a "hidden" quadrupolar order at low temperature [77]; the effect of a film geometry would likely lead to rich physics.
On the more theoretical front, there are many fundamental open questions about spin ice films. For example, bulk spin ice shows a complex phase diagram in an external magnetic field, with the physics strongly dependent on the field direction [78][79][80]; the effect of magnetic fields on films is likely to be similarly rich. The effects of transverse exchange also raises a host of interesting questions: in bulk spin ice this induces tunnelling between the ice states and stabilizes a U(1) quantum spin liquid, quantum spin ice [18]. However, the quantum analogue of the two dimensional Coulomb phase is fundamentally unstable [81], meaning that a direct two-dimensional analogue of quantum spin ice does not exist [82]. The fate of quantum spin ice films is thus unclear; possibilities include magnetically ordered states, valence bond solids [82] or potentially a quantum Z 2 spin liquid, as found in related models on the kagome lattice [83]. How the resulting state in this quantum case depends on the film geometry and the choice of orphan bond exchange presents many directions to pursue in future studies. Given the ability to readily grow high quality rareearth pyrochlore oxide films [20,22], one might expect such theoretical investigations to motivate a range of experimental studies which will, in return, undoubtedly fuel new sets of theoretical questions.
The study of films of pyrochlore magnets is a nascent field with many open questions, and rapid development could be expected in the near future. We our hope that this work will help shed light onto these systems and provide guidance and motivation for upcoming experimental and theoretical work.
where the indexing matches that of the sublattice vectors.

Interaction matrix
The explicit form of the interaction matrix V(q) for NN bulk spin ice is where the term proportional to the identity matrix makes the Lagrange multiplier λ in the large-N theory λ consistent with the stiffness parameter of the Coulomb phase. The so-called adjacency matrix A(q) is given by where c αβ ≡ cos[q · (r α − r β )].
and the sublattice vectors r α are given by for the first eight spins, and r α+8k = r α + kẑ, k = 1, .., L − 1 for the remaining spins. The corresponding conventional unit cell, showing the structure of stacked layers in theẑ direction (each layer made of parallel chains in the [110] or [110] direction, alternatively) is shown in Fig. 3.

Interaction matrix
Similar to the bulk theory, we define The adjacency matrix A(q ⊥ ) is a 8L × 8L matrix with a tridiagonal 8 × 8 block structure -that is, only the diagonal blocks and next-to-diagonal blocks are non-trivial. For spins in the same cubic unit cell, the diagonal 8 × 8 block reads where e αβ ≡ exp[iq ⊥ · r αβ ], and r αβ is the nearest-neighbor vector connecting sublattices α and β. The next-to-diagonal 8 × 8 blocks are given by: 0 · · · 0 e 06 e 07 0 · · · 0 e 16 e 17 . . . . . . 0 0 0 One can check that the complete adjacency matrix A(q) is indeed Hermitian. All other couplings are zero.

Including orphan bonds
When including orphan bonds, with an Hamiltonian given by Eq. (2), the interaction matrix becomes where the adjacency matrix corresponding to the orphan bond couplings, A O (q ⊥ ), has only the following non-zero matrix elements: This method decomposes the pyrochlore lattice into tetrahedral clusters; we label the sixteen states of each tetrahedron as S Q µ where Q = 0, ±1, ±2 is the charge and the index µ runs over the number of distinct states with the given charge. For Q = 0 one has six states, for Q = +1 or −1 one has four states each and for Q = +2 or −2 one has a single state each (see Fig. 10 for an illustration). Generally, we can write the partition function of a nearest neighbor model on the pyrochlore lattice as where I is a tetrahedron and ω(S Q I µ I ) is the statistical weight of the configuration on that tetrahedron. For bulk nearest neighbour spin ice one can define the weights where z ≡ e −2βJ , since the energy depends only on the charge Q I of a given tetrahedron, Films in the [001] direction can be implemented simply using this formalism. First, consider (bulk) nearest neighbor spin ice composed of L ⊥ × L ⊥ × L conventional cubic unit cells with periodic boundary conditions. The desired film geometry can then be realized by cutting the bonds that pass thorough a fixed plane with normalẑ, changing the remaining bonds on those cut tetrahedra to carry the orphan coupling J O . For J O /J > 0 we define the weights, ω + O (S Q µ ), on such "orphan" tetrahedron to be where z O ≡ e −2β|J O | . Similarly, for J O /J < 0 the weights, ω − O (S Q µ ), can be defined as after a constant shift of the energy. The remaining non-orphan tetrahedra simply have the bulk weights given in Eq. (D2) independent of J O /J. Next we define the probabilistic graph assignments that define the clusters, following the framework of Refs. [84,85]. To this end, we decompose a weight as where G is a graph defined on a tetrahedron and the ∆(S Q µ , G) = 0 or 1 are compatibility factors, with zero being incompatible and one being compatible. For our purposes, these graphs consist of isolated spins or possible pairings of two spins, as illustrated in Fig. 11(a). For the orphan tetrahedra, we do not include graphs that connect spins across the cut; the allowed graphs are shown in Fig. 11(b). The bulk tetrahedra and the J O /J > 0 orphan tetrahedra graphs are defined to be compatible with a state if the pairs of spins joined take on opposite values while, for the J O /J < 0 graphs, the two joined spins must be equal.
With these definitions, one can solve Eq. (D5) to obtain the graph weights W(G). As in the case of the partition function weights, we denoted the orphan tetrahedra to have graph weights as W ± O (G). A solution for the bulk tetrahedra is given in Ref. [48] as At low temperature, T J, one has z → 0 and thus only the three "icelike" graphs (G 2,a ) have non-zero assignment probability. In this limit, the algorithm (for the bulk case) reduces to a variant of the usual loop algorithm [49]. For the orphan tetrahedra one finds a solution These weights are positive and satisfy the required sum rules for any choice of J O [48]. At low temperature, T J O , one has z O → 0, with only the two pair graphs (G 2,0 O ) having nonzero probability. For J O /J > 0 this corresponds to continuing the usual loops along the surface, while for the J O /J < 0 case it corresponds to a loop where the alternation pattern is reversed when an orphan bond is encountered, as discussed in Sec. V [see Fig. 8(c)]. For J O /J = 0, one has z O = 1 and thus only the free graph (G 0,0 O ) has non-zero assignment probability. This corresponds to allowing the strings of alternating spins to end at the orphan bonds. Note that these probabilities factorize; we could also define the weights on the orphan bonds separately at each surface, rather than using a combined orphan tetrahedron.
The method then proceeds as usual, as discussed in Refs. [84,85,48]. A Monte Carlo step consists of first assigning graphs to each tetrahedron following the probabilities, W(G), given in Eq. (D6) and (D7). When assigned to the whole lattice, these graphs form clusters, in this case strings and loops, which are identified and then flipped or not flipped with equal probability [51].