Sovann En

Frédéric Jurie. UNSUPERVISED DEEP HASHING WITH STACKED Bruno Crémilleux

Frédéric Jurie

UNSUPERVISED DEEP HASHING WITH STACKED CONVOLUTIONAL AUTOENCODERS

Keywords: learning based hashing, unsupervised learning, convolutional autoencoder

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

This paper addresses the problem of unsupervised feature learning, with the motivation of producing compact binary hash codes that can be used for indexing images. Such a representation is essential for many computer vision problems. For instance, in image retrieval, such representations enable efficient computations and low storage when dealing with large datasets [START_REF] Xia | Unsupervised deep hashing for large-scale visual search[END_REF][START_REF] Junbo | Stacked what-where auto-encoders[END_REF][START_REF] Wang | Semisupervised hashing for scalable image retrieval[END_REF][START_REF] Lin | Deep learning of binary hash codes for fast image retrieval[END_REF]. For image classification, it offers a higher level of semantic information compared to engineered feature extraction techniques [START_REF] Masci | Stacked convolutional auto-encoders for hierarchical feature extraction[END_REF][START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF][START_REF] Diederik | Auto-encoding variational bayes[END_REF].

Broadly speaking, there are two main ways to design hashing functions [START_REF] Wang | A survey on learning to hash[END_REF]: locality sensitive hashing (LSH) and learning to hash. LSH consists in combining dataindependent hashing functions, while learning to hash relies on a supervised learning stage for producing hashing functions adapted to the input data, for a specific task. The former approach generally does not scale well, making it not suitable for large scale nearest neighbor search [START_REF] Wang | Learning to hash for indexing big data -a survey[END_REF]. In addition, because of the locality preserving structure, LSH may not be able to preserve the semantic similarity.

Learning to hash techniques [START_REF] Wang | Semisupervised hashing for scalable image retrieval[END_REF][START_REF] Liu | Supervised hashing with kernels[END_REF][START_REF] Xia | Supervised hashing for image retrieval via image representation learning[END_REF] have been proposed to overcome such drawbacks, relying on data annotations to train hashing functions which can map semantically similar examples to similar binary codes. Although they can help to address the aforementioned shortcomings, such approaches heavily rely on the availability of training data. This is the reason why a new trend in hashing is to consider semi-supervised or even unsupervised learning strategies [START_REF] Xia | Unsupervised deep hashing for large-scale visual search[END_REF][START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF][START_REF] Zhang | Ssdh: semi-supervised deep hashing for large scale image retrieval[END_REF][START_REF] Jegou | Product quantization for nearest neighbor search[END_REF]. We build on this line of works.

Autoencoders are good candidates for unsupervised learning, due to their ability to learn data manifolds from unlabeled data (see e.g., [START_REF] Zhang | Augmenting supervised neural networks with unsupervised objectives for large-scale image classification[END_REF]). These networks are mostly based on regular autoencoders, made of fully connected layers, while recent state-of-the-art results are now also considering the use of convolutional autoencoders [START_REF] Lecun | Deep learning[END_REF]. Inspired by the recent works of [START_REF] Junbo | Stacked what-where auto-encoders[END_REF][START_REF] Makhzani | Winner-takeall autoencoders[END_REF], this paper takes advantage of convolutional deep autoencoders to design a novel and efficient way to design deep-hashing networks. While the objective of these aforementioned works is to learn hierarchical representations, our work is to focus on the production of compact binary codes.

In this paper, we learn to represent images by compact and discriminant binary codes, through the use of stacked convolutional autoencoders, relying on their ability to learn meaningful structure without the need of labeled data [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF]. The layers are stacked to form a network which maps an input image hierarchically into a latent space, through a highly nonlinear function (encoder part), and reconstruct hierarchically the original data (decoder part). Regular stacked convolutional autoencoders does not produce binary codes. We force the central layer, the one containing the code, to contain only binary values by introducing a regularization term in the loss function. This is one of the key contributions of this paper.

The capability of our network to produce semantic preserving codes is validated on the MNIST dataset on an image retrieval and image classification tasks. Experimental results illustrate the competitiveness of our network compared to state-of-the-art approaches.

The rest of the paper is as follows: Section 2 reviews the related works, Section 3 presents the method and the architecture, Section 4 gives experimental validation on image retrieval and classification while Section 5 concludes the paper.

RELATED WORK

Pioneer works in hashing [START_REF] Gionis | Similarity search in high dimensions via hashing[END_REF] are based on the idea of Locality Sensitive Hashing (LSH), one of the most popular choices [START_REF] Wang | Learning to hash for indexing big data -a survey[END_REF]. Variants of LSH such as the multi-probe LSH [START_REF] Dong | Modeling lsh for performance tuning[END_REF] or the non-linear hashing projection of [START_REF] Dasgupta | Fast locality-sensitive hashing[END_REF] have also been investigated. Such approaches usually use hand-crafted features which does not provide an optimal solution for hashing.

Different from LSH, which relies on a random projection, learning based hashing takes advantages of modern learning algorithms to find an adapted hashing function on a given dataset. The hash functions can be developed in three different ways: supervised, semi-supervised and unsupervised learning. They differ from one to another in terms of the loss function and the constraints used. In supervised learning, the most important constraint is the preservation of semantic information. This is usually achieved by forcing the distance of the hash codes of similar (resp. dissimilar) data to be smaller (resp. bigger) than a margin [START_REF] Junbo | Stacked what-where auto-encoders[END_REF][START_REF] Lin | Deep learning of binary hash codes for fast image retrieval[END_REF][START_REF] Liu | Supervised hashing with kernels[END_REF][START_REF] Xia | Supervised hashing for image retrieval via image representation learning[END_REF][START_REF] Liu | Deep supervised hashing for fast image retrieval[END_REF][START_REF] Weiss | Spectral hashing[END_REF]. Semisupervised hashing uses both labeled and unlabeled date. For instance, [START_REF] Wang | Semisupervised hashing for scalable image retrieval[END_REF][START_REF] Wang | Semisupervised hashing for large-scale search[END_REF] designed a loss function to minimize the distance between examples of the same class and maximize the variance between unlabeled data. In [START_REF] Zhang | Ssdh: semi-supervised deep hashing for large scale image retrieval[END_REF], in addition to pairwise similarity losses, the authors model the neighborhood structure of the labeled and unlabeled data. This structure acts as a regularizer so that the neighborhood of the data will be preserved in the hashing space. Unsupervised learning does not exploit any high level semantic information. For instance, [START_REF] Gong | Iterative quantization: A procrustean approach to learning binary codes for largescale image retrieval[END_REF] tries to minimize the quantization errors produced by the binarization of the features. [START_REF] Jegou | Product quantization for nearest neighbor search[END_REF][START_REF] Gong | Iterative quantization: A procrustean approach to learning binary codes for largescale image retrieval[END_REF] design their hashing function based on the idea of preserving the neighborhood of the data in the Euclidean space.

Recent breakthrough in deep learning gradually attracted lots of researchers to propose new hashing techniques. [START_REF] Xia | Unsupervised deep hashing for large-scale visual search[END_REF] employs classical autoencoder to exploit the underlying structure of the data in an unsupervised manner followed by a RBM layer with constraints to reduce the dimension of the hamming space. [START_REF] Dilokthanakul | Deep unsupervised clustering with gaussian mixture variational autoencoders[END_REF] designed a new hashing technique based on variational autoencoders. Variational autoencoders assume the input data can be generated via a random probability distribution, which will be learned by the network.

Closely related to hashing techniques, [START_REF] Junbo | Stacked what-where auto-encoders[END_REF] proposes a stacked what-where autoencoder based on convolutional autoencoders in which the necessity of switches (what-where) in the pooling/unpooling layers is highlighted. Another closely related work is the one of [START_REF] Makhzani | Winner-takeall autoencoders[END_REF]. The authors utilize convolutional autoencoders but with an aggressive sparsity constraints. On each feature map, only the nb maximal activation is kept while the rest is put into lifetime sparsity (zero value). These closely related works, however, are designed for automatic features learning while our objective is to produce compact binary codes for hashing purpose.

MODEL ARCHITECTURE

Our model is based on a stacked convolutional autoencoder mapping input images into a compact latent space, through an encoder network, and reconstructing the original image through a decoder network. The network comprises of 3 convolutional autoencoders and one classical autoencoder. Each convolutional autoencoder can be summarized as: conv -relu-batch normalization-pooling (encoder) and unpooling -conv -relu. The classical autoencoder is based on fully connected layer and can be compactly described as: fcxx-fcnb-fcxx. Figure 1 illustrates the whole architecture.

Convolutional Encoder-Decoder

More formally, let us assume that our network has L convolutional layers. The hierarchy of outputs given by the encoder is denoted as x i ∈ R Wi×Hi×Ci , for the i th -level, where W × H is the size of the image and C the number of channels. On the other hand, the reconstructions given by the decoder are denoted as xi ∈ R Wi×Hi×Ci for the i th -level. x 0 is the input image and x0 is the output of the network. The activation of the k th layer of the encoder is defined by x k = max pool(σ(x k-1 * W + b)), while the activation of the k th layer of the decoder is given by : xk-1 = unpool(σ(x k * W + b)), where W / W are the weights of the convolution/deconvolution layer, and b/ b the biases. * denotes the convolutional operation. σ is the activation functions which are in our case rectified linear units except for the hashing layer in the middle of the network, and max pool is the max pooling operator (with a stride of 2 pixels in our case).

The parameters of the networks are learned with SGD, following a standard autoencoder optimization procedure where the loss function is the mean square error between input and reconstruction: minimize

W,b 1 = 1 N N j=1 xj 0 -x j 0 2 2 (1)
with N the amount of training samples. For simplicity, we drop all the subscripts of the weights W and biases b.

Fully connected central layer

Between the encoder and the decoder, we introduce a regular autoencoder made of 3 fully connected layers to further reduce the dimension of the features. The hidden layer of this autoencoder constitutes the code used to represent the images.

It should be noticed that we also experimented with fully convolutional encoder-decoder network, without noticing much difference in terms of performance networks with a fully connected autoencoders. In practice, the fully convolutional network is faster to train and easier to optimize. The first and last FC layer are also equipped with relu. The decoder network is responsible to output the feature maps and gradually transform it into its original dimensions followed by 1x1 convolutional layer with sigmoid activation.

As for the fully convolutional encoder-decoder network, the optimization of the fully connected central autoencoder is based on the mean square error as given in Equation (1).

Producing sparse binary codes

While our goal is to learn a mapping function transforming input images into compact binary codes, the codes produced by the network given above are vectors of real values. In this section, we show how to introduce regularization terms making these codes binary, sparse and well distributed.

The binarization of the codes cannot be simply obtained by thresholding the real values of the codes. Such a thresholding would introduce a non-differentiable operation in the network, and standard back-propagation techniques would not be usable. This is why we opted for introducing a binary relaxation term [START_REF] Liu | Deep supervised hashing for fast image retrieval[END_REF]. This binary relaxation term aims at forcing the activation of the middle layer to values as close as possible to binary values

[-1, 1]. Let b j = (b j 1 , b j 2 , .
. . , b j d ) T be the activation vector of the hashing layer, for image j, d being the dimension of the code. The binary relaxation is done using: minimize

W,b 2 = 1 N N j=1 d i=1 |b j i | -1 (2)
In order to minimize Eq. ( 2), the activation b i is pushed toward -1 or 1. The total loss function is then obtained by computing the weighted sum of 1 end 2 : minimize

W,b = 1 N N j=1 x0 -x 0 2 2 + α N N j=1 d i=1 |b j i | -1 (3)
Interestingly, autoencoders allow to represent the data in low-dimensional spaces. If the data can be reconstructed well from such low dimensional spaces, it means that these latent spaces contain enough high-level semantic information. In addition, making the codes sparse allows to index them more efficiently [START_REF] Xia | Unsupervised deep hashing for large-scale visual search[END_REF][START_REF] Wang | Semisupervised hashing for scalable image retrieval[END_REF][START_REF] Makhzani | Winner-takeall autoencoders[END_REF]. Assuming that we want only nb bits to be 1 and all the rest equal to -1, Eq. ( 3) can be rewritten as: minimize

W,b = 1 N N n=1 xn -x n 2 2 + α N N j=1 i∈imax j b j i -1 + i∈imin j b j i + 1
where imax j (resp. imin j ) are the index of the first nb (resp. last d -nb) maximum values of b j .

Theoretically, the nb bits should be distributed uniformly in the hashing layer. However, in practice, we notice that the network tends to focus on the same bits to encode all images. This makes our network having similar hash codes when adding the binarization term in the loss. Following [START_REF] Gong | Iterative quantization: A procrustean approach to learning binary codes for largescale image retrieval[END_REF], and to produce more balanced bits, we maximize the variance of each bit and force the bits to be pairwise uncorrelated. This can be done by setting the correlation matrix to be as close as possible to the identity matrix:

1 N b T .b = I (4)
It is worth mentioning that all the regularizations mentioned earlier are only applied to the hashing layer. However, as noted in [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF][START_REF] Makhzani | Winner-takeall autoencoders[END_REF], without any regularization or denoising, autoencoder can barely learn any useful information. This is the reason why the relu activation and the pooling layer are Table 1. Retrieval mAP@1000 as a function of the number of bits of the hash code, on the MNIST dataset. XXX-VGGF denotes the hashing methods built on the top of features produced by the pre-trained VGG-F network. Our method produces its best overall performance with 32 filters in the convolutional layers. used in the network. To some extent, these two layers can act as a good regularization to perturb the network not to learn to memorize all the information. In addition, the relu activation also serves as a powerful non-linear model for the network to uncover the underlying structure in the data manifold.

EXPERIMENTATIONS

Hereafter, we describe the dataset used in our experimentations, give some implementation details and finally present our results and comparisons with state-of-the-art.

Datasets:

The MNIST dataset consists of 70K 28×28 greyscale images of handwritten digits from 0 to 9. It consists of 60k training images and 10k testing images. Following [START_REF] Xia | Unsupervised deep hashing for large-scale visual search[END_REF][START_REF] Xia | Supervised hashing for image retrieval via image representation learning[END_REF][START_REF] Wang | Semisupervised hashing for large-scale search[END_REF], we randomly select 1k images (100 per class from the testing set) as our queries and keep the rest to learn the autoencoder. For image classification, we randomly selected a number of labeled data to train our classifier and test on the whole testing set. All the images are scaled so that their raw pixel intensities are in the range [0, 1].

Implementation: Each convolution layer is set to have {16, 32, 64} filters with kernel size of 5x5. Each pooling and unpooling are set to reduce the dimension of the feature maps by 2 and vice versa (i.e. their stride is of [START_REF] Junbo | Stacked what-where auto-encoders[END_REF][START_REF] Junbo | Stacked what-where auto-encoders[END_REF]). Our model can be summarized as conv64-conv64-conv64-fc128-fcX-fc128-conv64-conv64-conv64 (X ∈ {12, 24, 36, 48}). The network is first trained in a layer-wise fashion before being fined tuned as a whole. Next, we apply the binary relaxation constraints by gradually increasing α value from 1e -4 until the hash layer contains binary values (with a precision of 1e -3 ). Then, we fine tune again the network by gently increasing the β term (balancing bits regularizer) to make the hash codes well distributed. [START_REF] Junbo | Stacked what-where auto-encoders[END_REF] 0.118 0.058 0.042 0.025 CONV [START_REF] Masci | Stacked convolutional auto-encoders for hierarchical feature extraction[END_REF] ---0.076 Our method 0.12 0.06 0.05 0.05

Image Retrieval Task

Table 1 shows the retrieval mAP@1000 obtained with our method, on the MNIST dataset, and compares it with state of the art methods. The mAP@1000 refers to the calculation of mAP on the 1000 first returned results for each query [START_REF] Xia | Unsupervised deep hashing for large-scale visual search[END_REF][START_REF] Xia | Supervised hashing for image retrieval via image representation learning[END_REF]. The results are taken from the most recent works available on the MNIST dataset [START_REF] Zhang | Ssdh: semi-supervised deep hashing for large scale image retrieval[END_REF]. The proposed approach outperforms all the state-of-the-art of unsupervised method for hashing. The significant gain in performance fully justifies the competitiveness of our method compared to other recent unsupervised learning methods.

Image Classification Task

This section evaluates how our compact representation performs on classification tasks. Our model acts as a feature extractor, and we employ a simple k-NN classifier to show the robustness and discriminative power of our features. Table 2 shows the classification errors obtained with our method, on the MNIST dataset, and compares it with stateof-the-art approaches. It should be noticed that the compared methods are trained end-to-end using labeled data, while our method act as an unsupervised feature learning method followed by a k-NN. Despite the fact that we do not use any annotations when learning the encoders, the performance is better or comparable to the state-of-the-art supervised techniques. This shows that our model successfully extract meaningful information from unlabeled data.

CONCLUSIONS

We proposed a novel architecture based on deep convolutional autoencoders to learn compact binary hash codes. The convolutional layers help us to hierarchically learn an embedding containing high-level semantic information. To turn the central layer into binary codes, a constraint based on binary relaxation is integrated together with a bit balancing regularizer. Experimental results demonstrate the competitiveness of our approach over state-of-the-art methods in image retrieval. Image classification experiments show that our model can also act as a good unsupervised feature learning, comparable to end-to-end networks without requiring any supervision.

Fig. 1 .

 1 Fig.1. Illustration of our network architecture. The left, middle and right blocks corresponds to the encoder, the fully connected layers and the decoder part. Each block in the encoder part consists of a convolutional, relu, batch normalization and a pooling layer. The fully connected layers, in the middle, learn to map the output of the encoder into a low-dimensional and binary spaces. The first and last FC layer are also equipped with relu. The decoder network is responsible to output the feature maps and gradually transform it into its original dimensions followed by 1x1 convolutional layer with sigmoid activation.

Table 2 .

 2 Classification error rate on the MNIST dataset using our network as a feature extractor, as a function of the numbers of labeled examples used to train the classifier. The dimension of the features in our case is of 12.

	Number of training examples
	100	300	600	1000
	SWWA


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

Motivated by recent realizations of Dy2Ti2O7 and Ho2Ti2O7 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 Dy2Ti2O7 films.

Highly frustrated magnets host an astonishing array of exotic many-body phenomena, taking us far from the conventional paradigms of collective magnetic behavior [START_REF]Introduction to Frustrated Magnetism[END_REF]. 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 [START_REF] Isakov | [END_REF][3][4][5] with photon and magnetic-monopole excitations [5][6][7][8][9][10][11][12][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][17][18][19].

As vividly exposed in classic texts on electromagnetism [START_REF] Jackson | Classical Electrodynamics[END_REF], 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 [START_REF] Harris | [END_REF] 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 2 Ti 2 O 7 [22,23] and Ho 2 Ti 2 O 7 [24,25], along with the promising possibilities offered by spin ice heterostructures [26,27] and pyrochlore-iridates thin films [28][29][30][31][32][33][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 Open boundaries lead to orphan bonds (green) on the surfaces. Each spin amounts to a pair of ⊕ and ⊖ effective magnetic charges (±q) [8]. Negative δo favors magnetic surface charges Q = ±2q. Below the surface ordering temperature, at T < Tso, 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 < TF ). rnn is the nearest-neighbor distance.

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 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 ≡ σ i ẑi [40], where ẑi is the local easy-axis and σ i = ±1 [START_REF]crystal field doublets of Dy 3+ and Ho 3+ are defined through a spectral weight in which mJ = ±J dominate and whose angular momenta are J = 15/2 and J = 8, respectively[END_REF]. The spins interact via nearest-neighbor couplings J, dipolar interactions D, surface effects δ o on orphan bonds, and strain effects δ d on all bonds in the xy planes parallel to the surfaces over the whole system:

H = J ⟨ij⟩ σ i σ j + δ o ⟨ij⟩ orphan σ i σ j + δ d ⟨ij⟩xy σ i σ j (1)
+ Dr 3 nn i>j ẑi • ẑj r 3 ij - 3 (ẑ i • ⃗ r ij ) (ẑ j • ⃗ r ij ) r 5 ij σ i σ j .
δ o and δ d are perturbations in addition to the coupling J.

The bulk DSI model (δ o = δ d = 0) is characterized by a low energy band of states with extensive entropy [START_REF] Gingras | [END_REF]43] in which each tetrahedron satisfies the ice rules with two spins pointing inwards and two spins pointing outwards (2-in/2-out) [START_REF] Harris | [END_REF]. 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 ±Q. This is the dumbbell model where each moment is recast as a pair of positive/negative magnetic charges ±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 [START_REF] Gingras | Introduction to frustrated magnetism[END_REF] because of the open boundaries in the [001] direction, the presence of orphan bonds (δ o ̸ = 0) and strain effects (δ d ̸ = 0). Given the rapid variation of exchange constants with distance and

FIG. 2.
Specific heat C h for films (δ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 δo = -4 K. Temperature axes: log (main), linear (inset). Throughout, J = -1.24 K and D = 1.41 K. The transition TF approaches the bulk Tc ≈ 0.18 K and becomes first order [44,45] in the limit L h → ∞; the vertical blue line in inset.

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 [START_REF] Swendsen | [END_REF][START_REF] Geyer | Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface[END_REF] and a loop algorithm adapted to account for dipolar interactions [44,45] and monopoles [37,[START_REF] Lin | [END_REF]. The system size is L×L×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 ≫ L > L h . termediate 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 δ 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, ±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 [ o , the surface state transforms from charge vacuum to charge crystal with a checkerboard pattern in order to minimize the Coulomb potential [Figs. 1 and3]. A suitable order parameter is

M so = |M b | + |M t | and M b,t = α Q b,t α ϵ α /(2L 2 )
, where Q b,t α is a bottom (b) or top (t) surface charge and ϵ α = ±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 Z 2 degree of freedom with projection along ±[110]. The same holds for the tetrahedra on the bottom layer, with composite spins projected along ±[ 110]. These composite projections now form the famous 6-vertex model [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF] 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 = 1 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 longrange 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 δ o < δ c o , 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 δ o > δ c o , 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 Z 2 symmetry breaking of the antiferroelectric F -model [39] [Figs. 1 and 3(b)]. The order parameter is M F = i σ i η i /N ′ where η i = ±1 transcribes the 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 orphanbond couplings δ bot o , δ top o , 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

δ top o < δ c o < δ bot o
(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 Here δo = 0. Thickness: L h = 5 (5 nm) as in the thinnest sample of Ref. [22] lower temperature.

If the lattice parameters of the substrate are incommensurate with those of the spin-ice film, lattice distortions will invariably be induced on the surface which could penetrate into the film, inducing a "bond distortion field". In experiments on Dy 2 Ti 2 O 7 films [22] with surfaces 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 "polar catastrophe" with structural repercussions over the entire film [58].

We introduce bond distortions through δ d in Eq. ( 1) which is applied uniformly throughout the sample. For δ 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 δ d = +0.2 K in Fig. 4). As for negative δ d , it induces a topological "KDP" transition whose ground state is ferromagnetic [59][60][61]. The description of the diverse ground states resulting from the competition with dipolar interactions, as a function of δ 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 δ 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 agreement 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 theoretical benchmark for spin ice thin films. A surface termination perpendicular to the [001] direction generates strains and surface effects captured by orphan bonds. The latter plays the role of a chemical potential for emergent magnetic monopoles on the surfaces. We observe surface ordering, driven by long-ranged dipolar interactions, resulting in the crystalization of the magnetic surface charges. Over a finite temperature range, this surface ordering does not penetrate into the bulk, providing 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 exotic phase, recently reported in experiments on confined water ice [65] and artificial lattices designed by nanolithography [66][67][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 2 Ti 2 O 7 films, as proposed in Ref. [22]. Our model thus offers a useful first step and guide for future ab initio calculationsaimed at developing a quantitative understanding of thin film experiments.

It is interesting to compare the partial wetting observed 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][70][71][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][74][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 mechanism 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.

STRUCTURE FACTOR

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

S(q) ≡ 1 N N j=1 σ j e iq•r j 2 (1)
S(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 ⃗ q =[ 21 0] corresponds to diverging Bragg peaks at 250 mK, but only to short-range correlations at 600 mK (see insets).

FINITE-SIZE SCALING

As visible in the order parameters of Fig. 2, the transitions 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 h = 1. As for the hightemperature peak, its maximum is not constant as expected 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, algebraic, correlations in between the two surfaces (see the pinch points in Fig. 1), which influence the critical properties of the surface ordering. Further work is required to precisely address this question.

FIG. 1 .

 1 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 ⊕ and ⊖ effective magnetic charges (±q)[8]. Negative δo favors magnetic surface charges Q = ±2q. Below the surface ordering temperature, at T < Tso, 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 < TF ). rnn is the nearest-neighbor distance.

  FIG.3. Two-step ordering in thin films of thickness L h = 1 and δo = -4 K: (a) specific heat and entropy, (b) order parameters and density ρ 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 Tso ≈ 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 TF ≈ 300 mK (orange) in favor of the same antiferromagnetic ground state as the F -model[39]. Temperature axes: log.

Fig. 3 .

 3 (b)]. At low temperature, when varying δ o below δ c

4 3 4

 44 ln 4 3 [Fig.3(a)], which is the exact square ice entropy[START_REF] Lieb | [END_REF] (the composite spins account for 1/4 of the original degrees of freedom), and is visible as pinch points in the structure factor[52].

FIG. 4 .

 4 FIG. 4. Strain effects smoothly lifts the Pauling entropy for -0.2 K δ d < 0 down to 0.4 K, while promoting bulk-DSI order for δ d > 0 (e.g. δ d = +0.2 K).Here δo = 0. Thickness: L h = 5 (5 nm) as in the thinnest sample of Ref.[22]

FIG. 1 .

 1 FIG. 1. Structure factor in the [hk0] Fourier plane for L h = 1, δo = -4 K and δ 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).

FIG. 2 .

 2 FIG. 2. Finite-size effect in the two-step ordering in thin films of thickness L h = 1, δo = -4 K and δ 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.

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.

Supplemental Material: "Spin Ice Thin Film: Surface Ordering, Emergent Square Ice and Strain Effects"