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

HAL
open science

# Spin ice Thin Film: Surface Ordering, Emergent Square ice, and Strain Effects

L.D.C. Jaubert, T. Lin, T.S. Opel, P.C.W. Holdsworth, M.J.P. Gingras

► To cite this version:

L.D.C. Jaubert, T. Lin, T.S. Opel, P.C.W. Holdsworth, M.J.P. Gingras. Spin ice Thin Film: Surface Ordering, Emergent Square ice, and Strain Effects. Physical Review Letters, 2017, 118, pp.207206. 10.1103/PhysRevLett.118.207206. [hal-01528212](https://hal.science/hal-01528212)

HAL Id: hal-01528212

https://hal.science/hal-01528212

Submitted on 28 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Spin Ice Thin Film: Surface Ordering, Emergent Square Ice and Strain Effects

L. D. C. Jaubert,¹,² T. Lin,³ T. S. Opel,³ P. C. W. Holdsworth,⁴ and M. J. P. Gingras³,⁵,⁶

¹Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0495, Japan

²CNRS, Université de Bordeaux, LOMA, UMR 5798, 33400 Talence, France

³Department of Physics and Astronomy, University of Waterloo,
200 University Avenue West, Waterloo, Ontario, N2L 3G1, Canada

⁴Université de Lyon, Laboratoire de Physique, École normale supérieure de Lyon,
CNRS, UMR5672, 46 Allée d'Italie, 69364 Lyon, France

⁵Perimeter Institute for Theoretical Physics, 31 Caroline North, Waterloo, Ontario, N2L 2Y5, Canada

⁶Canadian Institute for Advanced Research, Toronto, Ontario, M5G 1Z8, Canada

(Dated: March 13, 2017)

Motivated by recent realizations of Dy₂Ti₂O₇ and Ho₂Ti₂O₇ spin ice thin films, and more generally by the physics of confined gauge fields, we study a model spin ice thin film with surfaces perpendicular to the [001] cubic axis. The resulting open boundaries make half of the bonds on the interfaces inequivalent. By tuning the strength of these inequivalent “orphan” bonds, dipolar interactions induce a surface ordering equivalent to a two-dimensional crystallization of magnetic surface charges. This surface ordering may also be expected on the surfaces of bulk crystals. For ultrathin films made of one cubic unit cell, once the surfaces have ordered, a square ice phase is stabilized over a finite temperature window. The square ice degeneracy is lifted at lower temperature and the system orders in analogy with the well-known F-transition of the 6-vertex model. To conclude, we consider the addition of strain effects, a possible consequence of interface mismatches at the film-substrate interface. Our simulations qualitatively confirm that strain can lead to a smooth loss of Pauling entropy upon cooling, as observed in recent experiments on Dy₂Ti₂O₇ films.

Highly frustrated magnets host an astonishing array of exotic many-body phenomena, taking us far from the conventional paradigms of collective magnetic behavior [1]. The formulation of the local frustrated constraints in terms of gauge fields has revolutionized our perspective on these systems. Depending on the system, such gauge structure can take the form of electromagnetism [2–5] with photon and magnetic-monopole excitations [5–13], or be similar to quantum chromodynamics [14] and linearized general relativity [15], as well as support phase transitions lying outside the Ginzburg-Landau-Wilson framework [16–19].

As vividly exposed in classic texts on electromagnetism [20], boundaries dramatically influence the behavior of gauge fields. It is therefore natural to ask what may be the role of boundary conditions in frustrated magnets described by gauge theories. Classical spin ice [21] presents itself as a nascent paragon to lay the foundations of such concepts, motivated by recent experimental studies of spin ice thin films of Dy₂Ti₂O₇ [22, 23] and Ho₂Ti₂O₇ [24, 25], along with the promising possibilities offered by spin ice heterostructures [26, 27] and pyrochlore-iridates thin films [28–34]. This work was further motivated by the repeated reports of exotic phenomena arising at the interfaces of oxide heterostructures [35].

In this paper, we study the dipolar spin ice model [36] in thin film (slab) geometry (slab-DSI) defined by Eq. (1) below. As illustrated in Fig. 1, the slab geometry renders the nearest-neighbor bonds on the surface inequivalent: either a bond belongs to a “bulk tetrahedron”, or it is an “orphan bond”, belonging to a “virtual tetrahedron” which has been “cleaved away” at the interface. Each

virtual tetrahedron may therefore carry a magnetic surface charge that can propagate to and from the surface into the bulk. The orphan bonds act as effective chemical potentials for surface charges, allowing the

<figure><img src="image_3.png" /><figcaption>FIG. 1. Spin ice thin film, thickness $L_h = 1$ in the [001] direction. Open boundaries lead to orphan bonds (green) on the surfaces. Each spin amounts to a pair of $\oplus$ and $\ominus$ effective magnetic charges $(\pm q)$ [8]. Negative $\delta_o$ favors magnetic surface charges $Q = \pm 2q$. Below the surface ordering temperature, at $T < T_{so}$, the central tetrahedra (orange) develop into a square ice phase, as seen by the projection of the spins in plane. Shown, one of the ground states below the F-transition ($T < T_F$). $r_{nn}$ is the nearest-neighbor distance.</figcaption></figure>
slab-DSI to undergo a surface charge-ordering transition. This is monopole crystallization [37] in two-dimensions, driven by the magnetic Coulomb potential between surface charges, which originates from the dipolar interactions. This crystallization is limited to the microscopic surface layer without penetration into the bulk. The thinnest slab where this phenomenology can be explored is one cubic unit-cell thick, containing three layers of tetrahedra. For such thickness, below the surface ordering temperature, $T_{so}$, the central layer emerges in the form of a constrained square ice system with extensive degeneracy and correlations characteristic of a two dimensional Coulomb phase [38]. Dipolar interactions ultimately lift the degeneracy of the Coulomb phase similarly to the venerable antiferroelectric F-model [39]. To conclude, we explore the possibility for strain effects to lift the residual Pauling entropy, as suggested by experiments on Dy$_{2}$Ti$_{2}$O$_{7}$ films [22].

**Slab Dipolar Spin Ice** – We consider thin films of the pyrochlore lattice of corner-sharing tetrahedra and whose surfaces are normal to the [001] cubic axis [Fig. 1]. The magnetic moments are described by classical Ising pseudospins $S_i \equiv \sigma_i \hat{z}_i$ [40], where $\hat{z}_i$ is the local easy-axis and $\sigma_i = \pm 1$ [41]. The spins interact via nearest-neighbor couplings $J$, dipolar interactions $D$, surface effects $\delta_o$ on orphan bonds, and strain effects $\delta_d$ on all bonds in the $xy$ planes parallel to the surfaces over the whole system:

$$ \mathcal{H} = J \sum_{\langle ij \rangle} \sigma_i \sigma_j + \delta_o \sum_{\langle ij \rangle_{\text{orphan}}} \sigma_i \sigma_j + \delta_d \sum_{\langle ij \rangle_{xy}} \sigma_i \sigma_j \quad (1) \\ + Dr_{nn}^3 \sum_{i>j} \left( \frac{\hat{z}_i \cdot \hat{z}_j}{r_{ij}^3} - \frac{3(\hat{z}_i \cdot \vec{r}_{ij})(\hat{z}_j \cdot \vec{r}_{ij})}{r_{ij}^5} \right) \sigma_i \sigma_j. $$

$\delta_o$ and $\delta_d$ are perturbations in addition to the coupling $J$.

The bulk DSI model ($\delta_o = \delta_d = 0$) is characterized by a low energy band of states with extensive entropy [42, 43] in which each tetrahedron satisfies the ice rules with two spins pointing inwards and two spins pointing outwards (2-in/2-out) [21]. This is a Coulomb phase with Pauling entropy [38]. Topological excitations out of the Coulomb phase carry either a single (3-in/1-out or 3-out/1-in) or a double (4-in or 4-out) gauge charge. Dipolar interactions endow these topological excitations with an effective (monopole) magnetic charge $\pm Q$. This is the dumbbell model where each moment is recast as a pair of positive/negative magnetic charges $\pm q$ sitting in the centers of the adjoining tetrahedra such that $Q = 2q$ [Fig. 1] [8]. Within the dumbbell model, the low energy ice states are exactly degenerate. The ordering in bulk-DSI [44, 45] is due to corrections to this dumbbell description [8, 43].

Slab-DSI differs from bulk-DSI [46] because of the open boundaries in the [001] direction, the presence of orphan bonds ($\delta_o \neq 0$) and strain effects ($\delta_d \neq 0$). Given the rapid variation of exchange constants with distance and

<figure><img src="image_4.png" /><figcaption>FIG. 2. Specific heat $C_h$ for films ($\delta_o = +4$ K in orange and $-4$ K in violet) of different thickness ($L_h = 1$ (▲), 2 (●), 3 (□)) and for bulk DSI (bulk-DSI, black). Inset: Low temperature $C_h$ for $\delta_o = -4$ K. Temperature axes: log (main), linear (inset). Throughout, $J = -1.24$ K and $D = 1.41$ K. The transition $T_F$ approaches the bulk $T_c \approx 0.18$ K and becomes first order [44, 45] in the limit $L_h \to \infty$; the vertical blue line in inset.</figcaption></figure>

geometry [40], we consider here a full range of parameter values, of the order of the bulk couplings themselves. Nevertheless, the dumbbell model remains a useful description if one also considers surface charges [Fig. 1]. To set an experimental context, we use $J = -1.24$ K and $D = 1.41$ K, as in a minimal model of Dy$_{2}$Ti$_{2}$O$_{7}$ [36].

**Method** – We use Monte Carlo simulations with parallel tempering [47, 48] and a loop algorithm adapted to account for dipolar interactions [44, 45] and monopoles [37, 49]. The system size is $L \times L \times L_h$ in units of the 16-site cubic unit cell, with thickness $L_h < L = 8$. We take open boundary conditions along the [001] cubic axis and periodic ones along the other two cubic axes. As in other simulations of dipolar spin ice [45], we use the Ewald summation method to describe the long-range dipolar interactions. The slab geometry is implemented here by inserting an empty space, as large as needed, between replicas in the [001] direction [50, 51]. We chose an empty space of $L_e = 1000$ unit cells. No difference was observed over the temperature range considered upon varying $L_e$ for $L_e \gg L > L_h$.

**Surface ordering** – We first set $\delta_d = 0$. By varying $\delta_o$ with respect to a threshold value $\delta_o^c = -J - 1.945(2)D$, one can favor either antiferromagnetic or ferromagnetic configurations on orphan bonds. For $\delta_o > \delta_o^c$, this extra surface energy scale manifests itself in the specific heat as a Schottky anomaly at $T \sim \delta_o$ [Fig. 2]. Bulk physics is then rapidly recovered as $L_h$ increases. For $\delta_o < \delta_o^c$, a sharp feature appears in the specific heat at an in-
<figure><img src="image_5.png" /><figcaption>FIG. 3. Two-step ordering in thin films of thickness $L_h = 1$ and $\delta_o = -4$ K: (a) specific heat and entropy, (b) order parameters and density $\rho$ of single charges (blue) inside the thin film, i.e. not on the surfaces. The latter is noticeably higher than for the bulk-DSI model (black). The surface ordering at $T_{so} \approx 900$ mK stabilizes a square ice model in the middle of the slab, as confirmed by the entropy in the intermediate region delimited by the dashed lines. Dipolar interactions ultimately lift the square-ice entropy at $T_F \approx 300$ mK (orange) in favor of the same antiferromagnetic ground state as the F-model [39]. Temperature axes: log.</figcaption></figure>

intermediate temperature $T_{so}$ whose height decreases with increasing $L_h$, suggesting a surface rather than a bulk transition, as discussed below.

One can think of $\delta_o$ as a chemical potential shift compared with the bulk, inciting charge on the surface. Indeed, the ferromagnetic (resp. antiferro.) alignment of the pseudospins on an orphan bond amounts to a net magnetic charge, $\pm Q$ (resp. vacuum) placed above the orphan bond at its mid point [Fig. 1]. The charge sites form a square array with a lattice constant of $2r_{nn}$. The dipolar interaction between spins generates an effective Coulomb interaction between the surface charges [8], so that each surface of the thin film can be thought of as a (square) Coulomb lattice gas. The surface charges can propagate into the bulk, which increases the density of monopoles above that of the bulk-DSI at intermediate temperature [Fig. 3(b)]. At low temperature, when varying $\delta_o$ below $\delta_o^c$, the surface state transforms from charge vacuum to charge crystal with a checkerboard pattern in order to minimize the Coulomb potential [Figs. 1 and 3]. A suitable order parameter is $M_{so} = |\bar{M}^b| + |\bar{M}^t|$ and $\bar{M}^{b,t} = \sum_\alpha Q_\alpha^{b,t} \epsilon_\alpha / (2L^2)$, where $Q_\alpha^{b,t}$ is a bottom (b) or top (t) surface charge and $\epsilon_\alpha = \pm 1$ accounts for the bipartite nature of the square array. Finite size effects are discussed in the Supplemental Materials [52].

Such surface order is not limited to thin films and

could also occur on a sufficiently pristine [100] surface of bulk crystals. The enhancement of monopole density at intermediate temperature should also persist within a thin layer below the surface of spin ice compounds, which could be manipulated by a magnetic field [10]. Near-surface dynamical properties could possibly be probed by β-NMR [53, 54] or µSR with slow muons [55].

**Emergent square ice for $L_h = 1$ –** Below $T_{so}$, the system enters a temperature regime in which most tetrahedra respect the ice rules. Considering a tetrahedron on the top layer, the orientation of the two surface spins becomes fixed by the surface ordering [Fig. 1]. The two lower spins are not fixed, but are coupled by the ice rules, forming a composite $\mathbb{Z}_2$ degree of freedom with projection along $\pm[110]$. The same holds for the tetrahedra on the bottom layer, with composite spins projected along $\pm[\bar{1}10]$. These composite projections now form the famous 6-vertex model [56] whose vertices correspond to the middle-layer tetrahedra respecting the ice rules [Fig. 1]. This is confirmed by an entropy plateau for intermediate temperature at $S_{sq} = \frac{1}{4}(\frac{3}{4}\ln\frac{4}{3})$ [Fig. 3(a)], which is the exact square ice entropy [57] (the composite spins account for 1/4 of the original degrees of freedom), and is visible as pinch points in the structure factor [52].

**Ground states** – Just as in bulk-DSI [44, 45], the degeneracy of the Coulomb phase is eventually lifted at sufficiently low temperature, here $T_F$, to give a long-range ordered ground state [Figs. 1 and 2(inset)]. However, the 12-fold degenerate ground states of bulk-DSI are inequivalent in slab geometry, as only 4 of them support surface charge order. When $\delta_o < \delta_o^c$, these 4 states are quasi-degenerate ground states; the energy remains quasi-degenerate if all spins on either one or both surfaces of the film are flipped. When $\delta_o > \delta_o^c$, orphan bonds favor 4 of the bulk-DSI ground states with no surface charges.

When $L_h = 1$, the transition at $T_F$ is equivalent to the $\mathbb{Z}_2$ symmetry breaking of the antiferroelectric F-model [39] [Figs. 1 and 3(b)]. The order parameter is $M_F = \sum_i \sigma_i \eta_i / N'$ where $\eta_i = \pm 1$ transcribes the $\mathbb{Z}_2$ ordering of the $N' = 8L^2$ spins of the middle layer of tetrahedra [Fig. 1]. The transition is characterized by a very smooth specific heat, as opposed to the sharp first-order singularity for thicker films where the square-ice mapping does not hold [inset of Fig. 2].

**Substrate and strain effects** – Since films are grown on a substrate, the surfaces, and their respective orphan-bond couplings $\delta_o^{\text{bot}}$, $\delta_o^{\text{top}}$, will be asymmetric unless the top surface is capped with the substrate material. Surfaces with different couplings will order at different temperatures and one could have frustrated ordering for $\delta_o^{\text{top}} < \delta_o^c < \delta_o^{\text{bot}}$ (or vice-versa). In particular, one might expect a regime of parameters where the charge order on one surface hinders the ordering of the entire film at
<figure><img src="image_6.png" /><figcaption>FIG. 4. Strain effects smoothly lifts the Pauling entropy for -0.2 K $\lesssim$ $\delta_d$ < 0 down to 0.4 K, while promoting bulk-DSI order for $\delta_d$ > 0 (e.g. $\delta_d$ = +0.2 K). Here $\delta_o$ = 0. Thickness: $L_h$ = 5 (5 nm) as in the thinnest sample of Ref. [22]</figcaption></figure>

lower temperature.

If the lattice parameters of the substrate are incom-
mensurate with those of the spin-ice film, lattice dis-
tortions will invariably be induced on the surface which
could penetrate into the film, inducing a “bond distortion
field”. In experiments on Dy₂Ti₂O₇ films [22] with sur-
faces normal to the [110] axis, such a substrate-induced
strain has been suggested as a possible mechanism for the
loss of residual Pauling entropy [22] compared to bulk
samples. We also note that for films of ionic crystals,
termination with a polar surface often leads to a “po-
lar catastrophe” with structural repercussions over the
entire film [58].

We introduce bond distortions through $\delta_d$ in Eq. (1)
which is applied uniformly throughout the sample.
For $\delta_d > 0$, the distortion enhances the degeneracy
lifting of the Coulomb phase driven by the dipolar
interactions [44, 45], such that the ordering occurs
at a higher temperature than in the bulk (see results
for $\delta_d$ = +0.2 K in Fig. 4). As for negative $\delta_d$, it
induces a topological “KDP” transition whose ground
state is ferromagnetic [59–61]. The description of the
diverse ground states resulting from the competition
with dipolar interactions, as a function of $\delta_d$ and $L_h$,
is beyond the scope of this letter. Nevertheless, in a
first theoretical step towards understanding the loss
of Pauling entropy observed in experiments [22], we
simulated a thin film with $L_h$ = 5 corresponding to the
thinnest sample of Ref. [22] [Fig. 4]. Indeed, for weak
distortion -0.2 K $\lesssim$ $\delta_d$ < 0, we find that the entropy
dips smoothly below the Coulomb phase plateau as the
temperature is reduced, without any phase transition
down to 0.4 K. Our model is thus in qualitative agree-
ment with the results of Ref. [22], even if it fails to
reproduce the precise shape of the specific heat, whose

Schottky peak is broader in experiments. The difference
could be due to the difference of surface orientation,
a strain gradient in the film, or the possible presence
of disorder, common in rare-earth pyrochlores [62, 63].
Because the KDP transition is topological in nature,
the addition of a dedicated *directed* worm algorithm
might be necessary for thermalization down to the
lowest temperatures [61, 64]. The absence of such global
coherent updates, which is likely in experiments, makes
the study of out-of-equilibrium effects in strained thin
films a promising open problem.

**Discussion** – In this work, we have provided a the-
oretical benchmark for spin ice thin films. A surface ter-
mination perpendicular to the [001] direction generates
strains and surface effects captured by orphan bonds.
The latter plays the role of a chemical potential for emer-
gent magnetic monopoles on the surfaces. We observe
surface ordering, driven by long-ranged dipolar interac-
tions, resulting in the crystalization of the magnetic sur-
face charges. Over a finite temperature range, this sur-
face ordering does not penetrate into the bulk, provid-
ing a magnetic example of partial wetting (see [Fig. 2]).
For ultra-thin films ($L_h$ = 1), the surface charge order
stabilizes an effective square-ice phase in the intervening
layers of spins. This offers a solid-state route of this ex-
otic phase, recently reported in experiments on confined
water ice [65] and artificial lattices designed by nano-
lithography [66–68]. Last but not least, our simulations
also confirm that strain coming from interface mismatch
can lead to a smooth loss of Pauling entropy in Dy₂Ti₂O₇
films, as proposed in Ref. [22]. Our model thus offers a
useful first step and guide for future *ab initio* calcula-
tionsaimed at developing a quantitative understanding
of thin film experiments.

It is interesting to compare the partial wetting ob-
served here to that occurring in water ice. Ice crystals
show a thin surface layer stabilized down to -40°C with
remarkably high charge mobility [69-72]. Although often
described as a liquid layer, its conductivity is orders of
magnitude larger than that of bulk liquid water. This
is admittedly a vastly more complex problem than our
model spin ice [73-75]. However, on cleaving a surface
from an ice crystal, one might expect surface charges to
be induced in an analogous manner. While the putative
charge crystal could then melt under the high Coulombic
pressure, it has also been proposed that the high surface
conductivity could come through a super-ionic mecha-
nism in an *ordered* ionic array [70], reminiscent of the
surface ordering we found here.

In conclusion, we believe the theory presented here
is just the tip of the iceberg of possibilities for spin ice
films and surfaces of single crystals. We hope our work
motivates further efforts in the investigation of surface
and confinement phenomena in strongly correlated
systems described by an emergent gauge theory.
We thank Bruce Gaulin, Gabriele Sala, Etienne Lantagne-Hurtubise, Rob Kiefl, Andrew MacFarlane, Jeff Rau, Nic Shannon and Tommaso Roscilde for useful discussions. We acknowledge the hospitality of the Ecole Normale Supérieure de Lyon, Université Lyon 1 and CNRS (LDCJ and MJPG) and from the Okinawa Institute of Science and Technology Graduate University (PCWH and MJPG). LDCJ is supported by the Okinawa Institute of Science and Technology Graduate University. The work at the U. of Waterloo was supported by the NSERC of Canada, the Canada Research Chair program (M.G., Tier 1) and by the Perimeter Institute (PI) for Theoretical Physics. Research at the Perimeter Institute is supported by the Government of Canada through Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science.

[1] C. Lacroix, P. Mendels, and F. Mila, eds., *Introduction to Frustrated Magnetism* (Springer, Heidelberg, 2011).

[2] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Physical Review Letters **93**, 167204 (2004).

[3] C. L. Henley, Phys. Rev. B **71**, 014424 (2005).

[4] M. J. Lawler, New Journal of Physics **15**, 043043 (2013).

[5] M. Hermele, M. P. A. Fisher, and L. Balents, Physical Review B **69**, 64404 (2004).

[6] A. Banerjee, S. V. Isakov, K. Damle, and Y. B. Kim, Phys. Rev. Lett. **100**, 047208 (2008).

[7] O. Benton, O. Sikora, and N. Shannon, Phys. Rev. B **86**, 075154 (2012).

[8] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature **451**, 42 (2008).

[9] I. A. Ryzhkin, Journal of Experimental and Theoretical Physics **101**, 481 (2005).

[10] L. D. C. Jaubert and P. C. W. Holdsworth, Nature Physics **5**, 258 (2009).

[11] D. J. P. Morris, D. A. Tennant, S. A. Grigera, B. Klemke, C. Castelnovo, R. Moessner, C. Czternasty, M. Meissner, K. C. Rule, J. U. Hoffmann, K. Kiefer, S. Gerischer, D. Slobinsky, and R. S. Perry, Science **326**, 411 (2009).

[12] T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. Morrow, and S. T. Bramwell, Science **326**, 415 (2009).

[13] H. Kadowaki, N. Doi, Y. Aoki, Y. Tabata, T. J. Sato, J. W. Lynn, K. Matsuhira, and Z. Hiroi, Journal of the Physical Society of Japan **78**, 103706 (2009).

[14] O. Cépas and A. Ralko, Phys. Rev. B **84**, 020413 (2011).

[15] O. Benton, L. D. C. Jaubert, H. Yan, and N. Shannon, Nature Communications **7** (2016).

[16] L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. **100**, 067207 (2008).

[17] S. Powell and J. T. Chalker, Physical Review B **78**, 24422 (2008).

[18] D. Charrier, F. Alet, and P. Pujol, Phys. Rev. Lett. **101**, 167205 (2008).

[19] S. Powell, Physical Review B **84**, 94437 (2011).

[20] J. D. Jackson, *Classical Electrodynamics* (Wiley, New York, 1999).

[21] M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Phys. Rev. Lett. **79**, 2554 (1997).

[22] L. Bovo, X. Moya, D. Prabhakaran, Y.-A. Soh, A. Boothroyd, N. Mathur, G. Aeppli, and S. Bramwell, Nature Communications **5**, 3439 (2014).

[23] O. Petrenko, Nature Materials **13**, 430 (2014).

[24] D. P. Leusink, F. Coneri, M. Hoek, S. Turner, H. Idrissi, G. Van Tendeloo, and H. Hilgenkamp, Applied Physics Letter Materials **2**, 032101 (2014).

[25] K. Kukli, M. Kemell, M. C. Dimri, E. Puukilainen, A. Tamm, R. Stern, M. Ritala, and M. Leskel, Thin Solid Films **565**, 261 (2014).

[26] T. Sasaki, E. Imai, and I. Kanazawa, Journal of Physics: Conference Series **568**, 052029 (2014).

[27] J.-H. She, C. H. Kim, C. J. Fennie, M. J. Lawler, and E.-A. Kim, ArXiv e-prints (2016), arXiv:1603.02692 [cond-mat.str-el].

[28] X. Hu, A. Rüegg, and G. A. Fiete, Phys. Rev. B **86**, 235141 (2012).

[29] B.-J. Yang and N. Nagaosa, Phys. Rev. Lett. **112**, 246402 (2014).

[30] E. J. Bergholtz, Z. Liu, M. Trescher, R. Moessner, and M. Udagawa, Phys. Rev. Lett. **114**, 016806 (2015).

[31] X. Hu, Z. Zhong, and G. A. Fiete, Scientific Reports **5**, 11072 (2015).

[32] T. C. Fujita, Y. Kozuka, M. Uchida, A. Tsukazaki, T. Arima, and M. Kawasaki, Scientific Reports **5**, 9711 (2015).

[33] T. C. Fujita, M. Uchida, Y. Kozuka, S. Ogawa, A. Tsukazaki, T. Arima, and M. Kawasaki, arXiv:1601.02710 (2016).

[34] K. Hwang and Y. B. Kim, ArXiv e-prints (2016), arXiv:1602.01836 [cond-mat.str-el].

[35] H. Y. Hwang, M. Iwasa, Y. amd Kawasaki, B. Keimer, N. Nagaosa, and Y. Yokura, Nature Materials **11**, 103 (2010).

[36] B. C. den Hertog and M. J. P. Gingras, Phys. Rev. Lett. **84**, 3430 (2000).

[37] M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, Phys. Rev. X **4**, 011007 (2014).

[38] C. L. Henley, Annual Review of Condensed Matter Physics **1**, 179 (2010).

[39] E. H. Lieb, Phys. Rev. Lett. **18**, 1046 (1967).

[40] J. G. Rau and M. J. P. Gingras, Phys. Rev. B **92**, 144417 (2015).

[41] The single-ion crystal field doublets of Dy$^{3+}$ and Ho$^{3+}$ are defined through a spectral weight in which $m_J = \pm J$ dominate and whose angular momenta are $J = 15/2$ and $J = 8$, respectively. As $m_J$ is maximal, we may assume that the Ising nature of the moment at the surface remains largely protected through the absence of high rank crystal field perturbations that would renormalize the doublet [40].

[42] M. J. P. Gingras and B. C. den Hertog, Canadian Journal of Physics **79**, 1339 (2001).

[43] S. V. Isakov, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. **95**, 217201 (2005).

[44] R. G. Melko, B. C. den Hertog, and M. J. P. Gingras, Phys. Rev. Lett. **87**, 067203 (2001).

[45] R. G. Melko and M. J. P. Gingras, Journal of Physics-Condensed Matter **16**, R1277 (2004).
[46] M. J. P. Gingras, "Introduction to frustrated magnetism," (Springer (Ed. Lacroix, Mendels & Mila), 2011) Chap. Spin Ice.

[47] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986).

[48] C. J. Geyer, Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface, p. 156 (1991).

[49] T. Lin, X. Ke, M. Thesberg, P. Schiffer, R. G. Melko, and M. J. P. Gingras, Phys. Rev. B 90, 214433 (2014).

[50] I.-C. Yeh and M. L. Berkowitz, The Journal of Chemical Physics 111, 3155 (1999).

[51] A. Bródka, Chemical Physics Letters 400, 62 (2004).

[52] See Supplemental Material for an illustration of the structure factor at low temperature and the finite-size effects in Monte-Carlo simulations.

[53] G. D. Morris, Hyperfine Interactions 225, 173 (2014).

[54] W. A. MacFarlane, Solid State Nuclear Magnetic Resonance 68-69, 1 (2015).

[55] P. Bakule and E. Morenzoni, Contemporary Physics 45, 203 (2004).

[56] R. J. Baxter, *Exactly solved models in statistical mechanics* (Dover Publications, 2007).

[57] E. H. Lieb, Phys. Rev. Lett. 18, 692 (1967).

[58] N. Nakagawa, H. Y. Hwang, and D. A. Muller, Nature Materials 5, 204 (2006).

[59] J. C. Slater, Journal of Chemical Physics 9, 16 (1941).

[60] E. H. Lieb, Phys. Rev. Lett. 19, 108 (1967).

[61] L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. 105, 087201 (2010).

[62] K. A. Ross, T. Proffen, H. A. Dabkowska, J. A. Quil- liam, L. R. Yaraskavitch, J. B. Kycia, and B. D. Gaulin,

Phys. Rev. B 86, 174424 (2012).

[63] G. Sala, M. J. Gutmann, D. Prabhakaran, D. Pomaran- ski, C. Mitchelitis, J. B. Kycia, D. G. Porter, C. Castel- novo, and J. P. Goff, Nature Materials 13, 488 (2014).

[64] S.-C. Lin and Y.-J. Kao, Phys. Rev. B 88, 220402 (2013).

[65] G. Algara-Siller, O. Lehtinen, F. C. Wang, R. R. Nair, U. Kaiser, H. A. Wu, A. K. Geim, and I. V. Grigorieva, Nature (London) 519, 443 (2015).

[66] G. Moller and R. Moessner, Phys. Rev. Lett. 96, 237202 (2006).

[67] I. Gilbert, G.-W. Chern, S. Zhang, L. OBrien, B. Fore, C. Nisoli, and P. Schiffer, Nature Physics 10, 670675 (2014).

[68] Y. Perrin, B. Canals, and N. Rougemaille, Nature 540, 410413 (2016).

[69] N. H. Fletcher, Philosophical Magazine 7, 255 (1962).

[70] I. A. Ryzhkin and V. F. Petrenko, Phys. Rev. B 65, 012205 (2001).

[71] Y. Li and G. A. Somorjai, The Journal of Physical Chemistry C 111, 9631 (2007).

[72] I. A. Ryzhkin and V. F. Petrenko, Journal of Experimental and Theoretical Physics 108, 68 (2009).

[73] N. H. Fletcher, Philosophical Magazine Part B 66, 109 (1992).

[74] V. Buch, H. Groenzin, I. Li, M. J. Shultz, and E. Tosatti, Proceedings of the National Academy of Sciences 105, 5969 (2008).

[75] D. Pan, L.-M. Liu, G. A. Tribello, B. Slater, A. Michaelides, and E. Wang, Phys. Rev. Lett. 101, 155703 (2008).
Supplemental Material:
“Spin Ice Thin Film: Surface Ordering, Emergent Square Ice and Strain Effects”

L. D. C. Jaubert,¹,² T. Lin,³ T. S. Opel,³ P. C. W. Holdsworth,⁴ and M. J. P. Gingras³,⁵,⁶

¹Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0495, Japan

²CNRS, Univ. Bordeaux, LOMA, UMR 5798, F-33400 Talence, France

³Department of Physics and Astronomy, University of Waterloo,
200 University Avenue West, Waterloo, Ontario, N2L 3G1, Canada

⁴Université de Lyon, Laboratoire de Physique, École normale supérieure de Lyon,
CNRS, UMR5672, 46 Allée d'Italie, 69364 Lyon, France

⁵Perimeter Institute for Theoretical Physics, 31 Caroline North, Waterloo, Ontario, N2L 2Y5, Canada

⁶Canadian Institute for Advanced Research, Toronto, Ontario, M5G 1Z8, Canada

(Dated: February 18, 2017)

STRUCTURE FACTOR

The structure factor corresponds to the Fourier trans-
form of the spin-spin correlations,

$$
S(\mathbf{q}) = \frac{1}{N} \left| \sum_{j=1}^{N} \sigma_j e^{i\mathbf{q} \cdot \mathbf{r}_j} \right|^2 \quad (1)
$$

$S(\mathbf{q})$ is an elegant observable to illustrate the two-step ordering discussed in the main text. Square-ice physics gives rise to pinch-point singularities at [200] and equivalent points [see white arrows in Fig. 1]. The pinch points are present at $T = 600$ mK, but absent at 250 mK, i.e. below the F-transition. Also, the scattering at $\vec{q} = [\bar{2}\bar{1}0]$ corresponds to diverging Bragg peaks at 250 mK, but only to short-range correlations at 600 mK (see insets).

<figure><img src="image_7.png" /><figcaption>FIG. 1. Structure factor in the [hk0] Fourier plane for $L_h = 1$, $\delta_o = -4$ K and $\delta_d = 0$. The insets (bottom left of each panel) are the same data than the main panel, but plotted on a broader color scale (×6).</figcaption></figure>

FINITE-SIZE SCALING

As visible in the order parameters of Fig. 2, the transi-
tions become sharper as the size of the system increases.
The maximum of the low-temperature peak in specific

heat is independent of the system size, as expected for
the F-transition observed when L<sub>h</sub> = 1. As for the high-
temperature peak, its maximum is not constant as ex-
pected for a transition in the 2D Ising universality class,
but slowly increases with system size. We interpret this
increase as the possible consequence of the strong, alge-
braic, correlations in between the two surfaces (see the
pinch points in Fig. 1), which influence the critical prop-
erties of the surface ordering. Further work is required
to precisely address this question.

<figure><img src="image_8.png" /><figcaption>FIG. 2. Finite-size effect in the two-step ordering in thin films of thickness $L_h = 1$, $\delta_o = -4$ K and $\delta_d = 0$: (a) specific heat divided by temperature and order parameters for (b) the surface ordering and (c) the F-transition. The temperature axis is in log scale. The monopole density and entropy are essentially independent of the system size.</figcaption></figure>