Fabien Alet

Michel Gingras

Peter Holdsworth

Claire Lhuillier

Grégoire Misguich

Les remerciements sont généralement la première (sinon la seule!) partie qui est lue dans une thèse. C'est malheureusement bien souvent également la dernière à être rédigée. Je ne ferai pas exception à cette dernière règle. La légèreté du ton n'ôte toutefois rien à la sincérité de ces lignes.

Je voudrais tout d'abord commencer Tommaso, sans qui cette thèse n'aurait même pas pu commencer. J'ai énormément appris grâce à lui, que ce soit en tant que directeur de thèse ou en tant que professeur durant le Master. Je me souviendrai encore longtemps de son cours sur les jonction Josephson qui a grandement facilité mon choix (judicieux) de directeur de thèse l'année suivante. Je tiens à le remercier pour son aide précieuse tout au long de ces quelques années, et plus particulièrement pour son enthousiasme, sa grande disponibilité et les longues discussions que nous avons pu avoir.

Je tiens à remercier les membres du jury

Introduction 1.1 Consequences of frustration in many-body systems

Frustration [3,[START_REF] Houtappel | Order-disorder in hexagonal lattices[END_REF] is a phenomenon that arises from competing interactions, whose energies cannot be simultaneously minimized. The origin of this competition can be a strong disorder (e.g. in spin glass systems [START_REF] Toulouse | Theory of frustration effect in spin-glasses .1[END_REF]), or the geometry of the lattice (e.g. in structurally ordered magnetic materials). As it can preclude the ordering of the system at low temperature, frustration is a source of degeneracy and disorder which, in turn, lead to new interesting physics.

In ordered materials, geometrical frustration has inspired numerous works over the past four decades, both theoretical and experimental. These studies have paved the way for the understanding of new exotic phases. At the classical level, a typical example is spin ice [START_REF] Gingras | Spin ice[END_REF]. Its ground state is macroscopically degenerate, leading to the absence of magnetic order; moreover its low-energy excitations have non-trivial properties, detectable in the low-temperature behavior of the system. In the case of quantum spins, the ground-state degeneracy is typically finite, but the ground state(s) might still be magnetically disordered. If it preserves all the symmetries of the frustrated Hamiltonian, it is then called a quantum spin liquid [2,6]. In classical spin ice this property is realized on average over the degenerate ground-state manifold, so that it can be called a classical spin liquid [START_REF] Henley | The "Coulomb Phase" in frustrated systems[END_REF]. Quantum spin liquids stem from systems with a classical ground state with large degeneracy. In these systems, quantum fluctuations select a unique ground state (or a finite number of degenerate ground states). These different spin liquid phases are often associated with emerging gauge theories at low energy.

Unfrustrated systems

In spin systems, frustration is associated with antiferromagnetic interactions [START_REF] Wannier | Antiferromagnetism. the triangular ising net[END_REF][START_REF] Chalker | Geometrically frustrated antiferromagnets: Statistical mechanics and dynamics[END_REF]. Before describing how the behavior of frustrated antiferromagnets is exotic, we first briefly recall the properties of unfrustrated spin systems. To this end, let us consider a simple Heisenberg antiferromagnet with coupling J > 0, on a cubic lattice. Because the lattice is bipartite, the ground state of the system is a Néel state in which all spins of a sublattice point in the same direction and spins of the other sublattice point in the opposite direction. This ground state is unique up to a global rotation of the spins in the system. The low-energy excitations above this state are long wave-length spin-waves. These Goldstone modes [START_REF] Goldstone | Field theories with superconductor solutions[END_REF] are a consequence of the symmetry breaking in the ground state. Their energy ω(k) varies linearly with k around the Néel points k N = (±π, ±π) and, in particular, vanishes at k N . These excitations (along with higher-energy non linear excitations) are responsible for the transition to a paramagnetic state as the temperature is increased. In dimension d > 2, the finite bandwidth of the excitations leads to a transition at a finite critical temperature T c to an ordered Néel phase [4].

Frustrated systems

In contrast, geometrical frustration can introduce a large degeneracy in the ground state. The common example to illustrate this effect is the case of three Heisenberg spins on a triangle. It is impossible to find a configuration of this system in which each spin is antiparallel to the other two. The solution to minimize the energy is then to arrange the spins in a way that they make an angle of 120 • with each other. This condition allows for the construction of two inequivalent states, with opposite "handedness" (up to a global rotation, see figure 1.1). The degeneracy observed at the level of a few spins might lead in a many-body system to an exponential growth of the degeneracy of the ground state with the system size (see section 1.2.1). ? =⇒ Figure 1.1: Heisenberg spins cannot be antiferromagnetically aligned on a triangle. The minimization of the energy of the configuration leads to the two configurations on the right. Furthermore, this large degeneracy of the ground state is often accompanied by a large density of lowfrequency modes above the different classical ground states. These modes will be thermally activated and may be able to destroy at a very low temperature any form of order which appears in the ground state.

Experimentally, a first evidence of frustration can be obtained from the measurements of the magnetic susceptibility χ. Indeed, at high temperature χ -1 varies linearly with the temperature

χ -1 ∝ T -Θ CW (1.1)
where Θ CW is the Curie-Weiss temperature. This temperature contains information about the interactions.

In particular, for antiferromagnetic interaction J > 0 (and hence for frustrated systems), the Curie-Weiss temperature is negative. In absence of frustration, it is close to the critical temperature T c (Θ CW = -3J/2 for the S = 1/2 cubic lattice Heisenberg model for a transition temperature T c /J = 0.946 ± 0.001. [START_REF] Sandvik | Critical temperature and the transition from quantum to classical order parameter fluctuations in the three-dimensional heisenberg antiferromagnet[END_REF]). In frustrated systems, the critical temperature is much lower than |Θ CW | (or even equal to 0). It is then clear that the importance of the frustration can be estimated experimentally from the frustration ratio f ≡| Θ CW |/T c [START_REF] Ramirez | Strongly geometrically frustrated magnets[END_REF]. Fundamental insight can be obtained from the susceptibility and, in general, from the magnetic structure factor which contains information on the spin-spin correlations. The structure factor is experimentally reconstructed via neutron scattering, and it provides the characteristic signature of the low-temperature phase of spin ice systems. These systems are the topic of the next section.

Spin ice 1.2.1 Degeneracy of the ground state of ice and spin ice

The possibility of experimentally observing degenerate classical ground states has been identified in the early 1930's. From their remarkable heat capacity measurements, Giauque and co-workers found that the common solid phase of ice (hexagonal I h ) possess a residual entropy S 0 in the limit of low temperatures. They determined that cooling a piece of ice from 0 • C down to ∼ 10 K leaves it with an entropy S 0 = 0.82 ± 0.05 Cal/deg • mol. This result, in apparent contradiction with the third law of thermodynamics, was soon after explained by Pauling, who estimated that the proton disorder in ice leads to a residual entropy S 0 = 0.81 Cal/deg • mol [START_REF] Pauling | The structure and entropy of ice and of other crystals with some randomness of atomic arrangement[END_REF], very close to the experimental value.

• • • • • • • • • • O • H ↑ ↑ ↓ ↓ ↓ ↑ ↓ ↑ ↑ ↓ ↑ ↑ ↓ ↑ ↓ ↓

Ice rules and residual entropy

The residual entropy stems from simple principles. In the I h phase of ice, the oxygen atoms crystallize in a diamond lattice (see figure 1.2). Each of these atoms has four nearest neighbors connected by hydrogen bonds. The hydrogen bond is asymmetric (indeed, the O -H -O bond length is 2.76 Å whereas the covalent O -H bond length is 0.96 Å, less than half of this). On such a bond, the proton sits closer to one of the oxygens, to which it is covalently bound. The low-energy proton configuration of ice must satisfy the Bernal-Fowler ice rules [START_REF] Bernal | A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions[END_REF]. The first ice rule states that the water molecule are connected by hydrogen bonds, so that there is exactly one proton on each O -O bond. The second ice rule asserts that exactly two of the four protons linked to a given oxygen atom sit close to this atom and two sit away. Indeed, breaking this rule would amount to creating a HO -or a H 3 O + molecule, and the energy of the covalent O -H bond (∼ 10 eV) precludes this from happening.

Pauling used these rules to determine the residual entropy of water. His argument is the following. Consider a piece of ice containing N water molecules, and hence 2N hydrogen atoms. According to the first ice rule, each of these protons has two possible positions. There are then 2 2N such configurations for the whole system. Of the 2 4 = 16 possible configurations around a single oxygen, only 6 satisfy the second ice rule. Neglecting the correlations between tetrahedra, we get an upper bound for the number Ω of ice-rule configurations :

Ω ≤ 6 16 N 2 2N = 3 2 N . (1.2)
This gives an entropy per molecule s 0 = k B log Ω/N = k B log(3/2), and hence an entropy per mole S 0 = R log(3/2) = 0.81 Cal/deg • mol. The entropy has since been calculated in a more controllable and more precise way [START_REF] Nagle | Lattice statistics of hydrogen bonded crystals. i. the residual entropy of ice[END_REF], leading to S 0 = 0.8145±0.0002 Cal/deg•mol. Both estimates are in very good agreement with Giauque measurements. Many other occurrences of systems with residual entropy have been found since then [START_REF] Kozliak | Residual entropy, the third law and latent heat[END_REF], and their link with discrete spin systems has been characterized theoretically [START_REF] Griffiths | Microcanonical ensemble in quantum statistical mechanics[END_REF][START_REF] Leff | Proof of the third law of thermodynamics for ising ferromagnets[END_REF]1,[START_REF] Chow | Residual entropy and validity of the third law of thermodynamics in discrete spin systems[END_REF].

In particular, residual entropy is found at low temperature in many spin systems that have a crystal structure similar to that of the I h phase of ice, as we will discuss in the following.

Equivalent spin systems Anderson model

In his study of cation ordering in the inverse spinels materials, Anderson introduced the antiferromagnetic Ising model on the pyrochlore lattice (see figure 1.2). This lattice is composed of corner-sharing tetrahedra arranged in a diamond lattice, with an Ising spin attached to each corner. This model is frustrated, as all 6 links of a tetrahedron cannot be simultaneously in an antiferromagnetic configuration. The lowest energy configurations of a tetrahedron are such that only two links are frustrated, leading to the Ising version of the second ice rule (the first one is here enforced by definition of the spins) : each tetrahedron contains two up spins and two down spins. This model is formally equivalent to that of the proton configuration in ice, and should therefore exhibit the same properties. This was the first three dimensional frustrated antiferromagnet ever studied. The corresponding Heisenberg model was introduced by Villain [START_REF] Villain | Insulating spin glasses[END_REF] who predicted that it would remain disordered down to zero temperature, calling this phase a collective paramagnet [START_REF] Villain | Insulating spin glasses[END_REF], a fact that has since been verified [START_REF] Reimers | Absence of long-range order in a three-dimensional geometrically frustrated antiferromagnet[END_REF][START_REF] Moessner | Low-temperature properties of classical geometrically frustrated antiferromagnets[END_REF].

Three-dimensional spin ice

The experimental study of spin ice really started in 1997, when Harris and co-workers measured the structure factor of Ho 2 Ti 2 O 7 using neutron scattering and did not find any sign of long-range order down to temperatures much lower than the Curie-Weiss temperature Θ CW [START_REF] Harris | Geometrical frustration in the ferromagnetic pyrochlore ho 2 ti 2 O 7[END_REF]. This material has a pyrochlore structure, with magnetic Ho 3+ located on the corners of the tetrahedra, like in the Anderson model. However, the spins interact here with ferromagnetic Heisenberg interaction. The frustration comes from the strong crystal field that forces the spins of the Ho ions onto the m J = |±8 doublet, leading to a magnetic moment of order 10µ B . This large moment is a generic property of rare-earth magnetic materials. As a consequence, the spins align along the easy-axis connecting the centers of neighboring tetrahedra (see figure 1.2). The correspondence to the I h phase of ice is even more striking than in the case of the Anderson model. Indeed, the direction in which a spin points is directly related to the position of the corresponding proton in ice relatively to the center of the bond, whence the name spin ice.

The smoking-gun evidence of the equivalence between rare-earth pyrochlore compounds and ice was given in Ref. [START_REF] Ramirez | Zero-point entropy in 'spin ice[END_REF], in which precise measurements of the specific heat of Dy 2 Ti 2 O 7 were used to determine its residual entropy per tetrahedron S 0 = 0.93 Cal/deg • mol, which is consistent with Pauling's estimate of the entropy of ice (figure 1.3). More recent measurements give an even closer value [START_REF] Ke | Nonmonotonic zero-point entropy in diluted spin ice[END_REF]. Even though, even more recent measurements show that careful annealing of the magnetic state of Dy 2 Ti 2 O 7 is able to access a lower entropy state, possibly displaying long-range order due to residual dipolar interactions [START_REF] Pomaranski | Absence of pauling/'s residual entropy in thermally equilibrated dy2ti2o7[END_REF].

Two-dimensional ice and 6-vertex model

In this section, we concentrate on a two-dimensional version of spin ice, namely square ice. This model consists of a set of Ising spins sitting on the bonds of a square lattice (see figure 1.4). Similarly to the pyrochlore ice, the energy is determined by the configuration of the spins converging at a vertex. Particular variants of this model are relevant to describe the properties of copper formate tetrahydrate (CFT) and copper formate tetradeuterate (CFTD) [START_REF] Youngblood | Correlations in ice-rule ferroelectrics[END_REF]. Neutron diffraction [START_REF] Dorner | Inelastic neutron scattering from ice and other proton-containing substances[END_REF] and heat capacity [START_REF] Matsuo | Heat capacities of copper(ii) formate tetrahydrate and tetradeuterate: A comparative study of phase transitions in layer hydrate crystals[END_REF] showed that the latter is indeed a frustrated proton ice system. When restricting to configuration respecting the ice-rule (2in/2out) on each vertex, one gets the celebrated 6-vertex model [START_REF] Baxter | Exactly solvable models in statistical physics[END_REF]. This model was exactly solved by Lieb [START_REF] Lieb | Residual entropy of square ice[END_REF], and its spin correlations were computed in [START_REF] Sutherland | Correlation functions for two-dimensional ferroelectrics[END_REF][START_REF] Youngblood | Correlations in ice-rule ferroelectrics[END_REF][START_REF] Youngblood | Polarization fluctuations in ferroelectric models[END_REF].

This model shares many properties of three-dimensional spin ices. It exhibits a residual entropy per spin S 0 /k B = 3 2 log 4 3 ≈ 0.857 Cal/deg•mol, very close to the one of pyrochlore spin ice. This comes from the fact that the 6-vertex model can be viewed as a two-dimensional arrangement of "crushed" corner-sharing tetrahedra, with the same coordination as in the pyrochlore lattice (for this reason the two-dimensional 6-vertex model and other related models are often referred to as planar pyrochlore lattice models). Surprisingly, Pauling's approximation for the calculation of the entropy is still quantitatively correct in two dimensions. Figure taken from [START_REF] Ramirez | Zero-point entropy in 'spin ice[END_REF].

Topological properties and collective paramagnet Loops in ice systems

A key property of spin ice systems stemming from the ice rules concerns the way different ground state configurations are connected to each other. Let us consider two ice-rule configurations {σ i } and σ i on the square ice, and compute σi = σ iσ i on each bond (represented in figure 1.4). Because of the 2in-2out constraint, the non-zero elements of { σi } compose a set of closed strings (loops) of head-to-tail spins (if the boundary conditions are periodic, otherwise the strings could start and end at the boundary). The system goes from {σ i } to σ i (and vice-versa) by flipping all the loops in { σi }. This property has been widely used to sample the ice-rule states in Monte Carlo simulations (see section 3.4).

{σ i } {σ i } { σi }≡{σ i }-{σ i } Figure 1
.4: Two ice-rule configurations {σ i } and {σ i } and their difference { σi }≡{σ i }-{σ i }.

Topological sector

The ice-rule states can be classified into topological sectors, depending on their winding number w. This number has two components in two dimensions, one for each spatial direction (w x , w y ). It is defined as the total magnetization m x/y along a vertical cut (for w x ) or a horizontal cut (for w y ) line, as described in figure 1.5. This number is the same for all horizontal (vertical) lines, making it a topological invariant of the configuration. w can only be changed by the flip of a winding loop, i.e. a loop that goes around the system. Such a loop is topologically non trivial as it is not contractible (i.e. cannot be reduced to a point by successive local moves, such as the flip of the four spins of a given plaquette). The configuration space is then composed of topological sectors connected to each other by the flip of winding loops. Flipping only topologically trivial loops will force the system to be confined to a given topological sector.

w y = 1 w x = -3 Each component of the winding number is obtained as the total magnetization of the spin cutting a vertical (w x ) or horizontal (w y ) line. In this example we have w x = -3 and w y = 1. It is easily verified that it is indeed the same number for each line/column.

Collective paramagnetism

Because the winding number of an ice-rule state is the same along each cut, it is directly related to the total magnetization M = (M x , M y ). Indeed, M = (L × w x , L × w x ) = Lw x . The fluctuations in the winding number are then linked to the structure factor at q = (0, 0)

w 2 -w 2 = 1 L 2 M 2 -M 2 ≡ S (0, 0) . (1.3)
Most of the key features of the Coulomb phase can still be computed correctly without leaving the (0, 0) topological sector. Indeed, we observe that restricting the system to the (0, 0) topological sector only changes the structure factor of the 6-vertex model at the origin (see figure 1.6). This is precisely because, in this case, the winding number fluctuations vanish by construction. At finite temperature, the ice rules can be locally violated. Therefore the winding number ceases to be well defined. Nevertheless, the fluctuations due to winding loops flip can be tracked through the calculation of the structure factor at the origin. This height is independent of the system size, up to vanishing corrections [START_REF] Sutherland | Correlation functions for two-dimensional ferroelectrics[END_REF][START_REF] Youngblood | Correlations in ice-rule ferroelectrics[END_REF]: S (0, 0) = σ 2 0 + O(1/L) . (1.4) This property holds for both two-dimensional and three-dimensional ice systems. For the 6-vertex model, we find with, by a Monte Carlo simulations, σ 2 0 ≈ 1.89. Hence the low-temperature susceptibility per spin χ should reads

χ ≡ β L 2 M 2 -M 2 ≈ 1.89β. (1.5)
For an Ising paramagnet, the Curie law gives χ x = β. The susceptibility is here 1.89 times larger. The system is then qualified as being a collective [START_REF] Villain | Insulating spin glasses[END_REF] or a cooperative [START_REF] Moessner | Low-temperature properties of classical geometrically frustrated antiferromagnets[END_REF] paramagnet. This is a generic property of the Coulomb phase [START_REF] Moessner | Magnets with strong geometric frustration[END_REF], and of classical spin liquids in general [2,6]. Ryzhkin predicted σ 2 0 = 2 Figure 1.6: S q ≡ M x (-q)M x (q) + M x (-q)M y (q) (which is proportional to the structure factor of the projection of the magnetization onto the x + ŷ direction) averaged on the complete ice states manifold (left) and on the (0, 0) topological sector(right). There seems to be no visible difference between the two figures. However, there is a dramatic change at the pinch point. The difference between the two structure factors (middle) is zero every where but at the pinch points. The black square is the first Brillouin zone.

for pyrochlore ice [START_REF] Ryzhkin | Magnetic relaxation in rare-earth oxide pyrochlores[END_REF], as found three dimensional spin ice compounds [START_REF] Jaubert | Topological-sector fluctuations and curie-law crossover in spin ice[END_REF]. This value is remarkably close to the two-dimensional values.

Coulomb phase

Spin correlations in spin ice

In the language of spin ice, the ice rules impose that on each tetrahedron two spins are pointing outwards and two are pointing inwards (2in-2out rule). This condition has a very specific consequence on the correlations in the system. In particular, this simple local constraint implies that the spin-spin correlations decay algebraically, and that it takes, in three dimensions, the functional form of a dipole-dipole interaction. This property has a very salient manifestation in the structure factor. Its expression can be derived via a field theory obtained by coarse-graining the magnetization in the system. To do so, let us note M i the magnetization of site i, defined as a unit vector joining the two adjacent tetrahedra of i and oriented according to the pseudospin on i (this unit vector can be thought of as the spin on site i). We now consider the total magnetization M(r) = 1

V D i∈D M i of a domain D centered at position r. D is chosen large enough so that M(r) is a smooth function.

For an ideal spin ice model, with a perfect degeneracy of all the ice states, the Gibbs free energy G is purely entropic. For a large value of M, most of the links within D are polarized in the same direction, whereas for M = 0 they can take any orientation so long as their sum vanishes. It is then clear that the number of ice configuration of D with fixed magnetization M is maximal for M = 0. Consequently, if the magnetization fluctuations are sufficiently large (i.e. if the correlations between domains decay fast enough), the central limit theorem holds, so that the distribution of magnetizations approaches a Gaussian. The free energy of the total system then reads

G [{M(r)}] ≈ T d d r 1 2 KM(r) 2 (1.6)
which, in Fourier space, corresponds to G {M(q)} ≈ T d d q 1 2 KM(q) 2 .

(1.7)

K is here an empirical prefactor that characterizes the width of the magnetization distribution. This energy is easily minimized by a field M (G) (r) whose correlations read M (G) µ (-q)M (G) ν (q) = 1 K δ µν , which corresponds to a paramagnetic state.

If all vertices within D obey the ice rules, then as many links crossing the boundary of D are oriented outward as inward. This amounts to imposing a divergence-free condition on M ∇•M= 0 (1.8) In Fourier space, this translates into q •M(q) = 0 .

(1.9)

Using standard projection method [START_REF] Temam | Une méthode d'approximation de la solution des équations de navier-stokes[END_REF][START_REF] Chorin | Numerical solution of the navier-stokes equations[END_REF], we find the field M(r) that minimizes (1.7) under the icerule constraint. It is obtained as the orthogonal part of the solution M (G) (r) to the same Gaussian problem without the constraint M(q) = M (G) (q) -M (G) (q) • q q .

(1.10)

The correlations between magnetization components are then given by [START_REF] Youngblood | Correlations in ice-rule ferroelectrics[END_REF] M µ (-q)M ν (q) = 1 K δ µν -q µ q ν q 2 (1.11)

and they now describe a state which differs substantially from a trivial paramagnet. This expression has a singular behavior around the q = 0 point. More precisely, along the q x = 0 line M x (-q)M x (q) ≡ 0, whereas along the q y = 0 line M x (-q)M x (q) ≡ 1. This property appears in all correlation functions. The structure factor then looks like it has been "pinched" at the origin (see figure 1.7, where the structure factor is represented for a two-dimensional system). Hence this particular feature has been named a "pinch M x (-q).M x (q) M y (-q).M y (q) M x (-q).M y (q) point" in the ice and spin ice literature. Pinch points were first measured in CFDT (a two-dimensional ice) by Youngblood and Axe [START_REF] Youngblood | Polarization fluctuations in ferroelectric models[END_REF][START_REF] Youngblood | Polarization fluctuations in ferroelectric models[END_REF]. In three dimensions they were first measured in Ho 2 Ti 2 O 7 [START_REF] Bramwell | Spin ice state in frustrated magnetic pyrochlore materials[END_REF], and they compare with the result of Monte Carlo calculations. Recent measurements on the same compound show even clearer pinch points [START_REF] Fennell | Magnetic coulomb phase in the spin ice ho2ti2o7[END_REF] (figure 1.8). Note that in 2d the height of this pinch point is given by the winding-number fluctuations defined in (1.4).

Transforming (1.11) back to real space, one obtains an algebraic decay of the correlations for sufficiently large distance r

C µν (r) = M µ (0).M ν (r) ∝ 1
Kr d δ µνdr µ rν .

(1.12) Pinch points structures appear at the [START_REF] Polyakov | Quark confinement and topology of gauge theories[END_REF] points. Figures taken from [START_REF] Fennell | Magnetic coulomb phase in the spin ice ho2ti2o7[END_REF].

Because of the coarse-graining approach, all details about the underlying lattice have been lost in the derivation of the spatial correlations. The result is then fairly general and can be applied to a whole family of models that possess the same key ingredients : a disordered coarse-grained vector field that satisfies a divergence-free condition [START_REF] Henley | The "Coulomb Phase" in frustrated systems[END_REF]. The quadratic form of the free energy and the divergence-free conditions resemble what is found for a magnetic field vacuum or for an electric field in vacuum. The possibility of associating a vector potential A to M such that M = ∇×A allows to identify in this system an emergent gauge structure. In particular, the theory with a free energy of the form 1.6 supplemented with the Gauss law condition (∇•M = 0) is a paradigmatic example of a classical gauge theory, and in particular of a lattice gauge theory, when M is ascribed back to the original lattice degrees of freedom. As we will see in the following, the emergent gauge field is not able to confine the gauge charges and therefore the lowtemperature phase of the system corresponds to a deconfined phase of a lattice gauge theory [START_REF] Polyakov | Quark confinement and topology of gauge theories[END_REF], the so-called Coulomb phase. A fundamental difference with ordinary lattice gauge theories is that the field retains a finite value on each site, even in the ground state of our system. This is due to the discreteness of the lattice degrees of freedom : it is therefore what could be called a "frustrated" lattice gauge theory, with exponential degeneracy of the ground state.

Monopole excitations

At finite temperature, the ice-rule are not strictly enforced, and defects can occur, corresponding to thermally activated single-spin flips. Replacing the spins by elongated dipoles with opposite charges sitting at the center of the tetrahedra/vertices leads to the so-called dumbbell model [START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF]. In this model, the ice rules impose the neutrality of each vertex. In this language, the fundamental spin excitations are then pairs of monopoles of opposite charge, i.e. vertices with a net charge, obtained from the flip of a dipole. Once created, the two monopoles can move away from each other and thus fractionalize [START_REF] Laughlin | Nobel lecture: Fractional quantization[END_REF] the spin-flip excitation [START_REF] Castelnovo | Spin ice, fractionalization, and topological order[END_REF].

If one considers the dipolar spin-spin interactions in spin ice, the monopoles inherit an effective interaction which has a Coulombic form. Because of the effective expression of the free energy (1.6), there is an additional Coulombic coupling between the monopoles that is temperature dependent [START_REF] Castelnovo | Debye-Hückel theory for spin ice at low temperature[END_REF].

Artificial spin ice

In the previous section, we presented various remarkable properties of pyrochlore spin ice. In particular, it is a rather standard procedure to determine experimentally the residual entropy and/or the magnetic structure factor of pyrochlore spin-ice compounds. However, the actual detection of monopoles and the measurement of their properties are still very difficult [START_REF] Bramwell | Measurement of the charge and current of magnetic monopoles in spin ice[END_REF][START_REF] Dunsiger | Spin ice: Magnetic excitations without monopole signatures using muon spin rotation[END_REF][START_REF] Sala | Magnetic coulomb fields of monopoles in spin ice and their signatures in the internal field distribution[END_REF][START_REF] Bovo | Brownian motion and quantum dynamics of magnetic monopoles in spin ice[END_REF]. A way to observe the monopoles unambiguously would be to probe the magnetization of the system locally in real time. It is a very difficult task in a three-dimensional system, but it is conceivable in two-dimensional systems.

Dipolar nanoarrays

The hunt for a two-dimensional analogue to spin ice has recently gained momentum after the work of Wang and collaborators [START_REF] Wang | Artificial 'spin ice' in a geometrically frustrated lattice of nanoscale ferromagnetic islands[END_REF] (figure1.9), which realized artificial spin ice [START_REF] Nisoli | Artificial spin ice: Controlling geometry[END_REF]. This system consists of lithographically fabricated single-domain ferromagnetic nanostructures that behave like giant Ising spins, interacting via dipolar couplings. The lithographic techniques developed for these systems enable the design of all kinds of planar lattices (honeycomb, kagome, square...). Furthermore, the configuration of the system can be directly measured, e.g. by magnetic force microscopy (MFM), which allows for direct measurements of the dynamics. The realization of these systems was rendered possible by the tremendous progress in the field of nanomagnetism (for a review on this point, see [START_REF] Bader | Colloquium : Opportunities in nanomagnetism[END_REF]).

However, most of these artificial spin ice have an intrinsic limitation : there is a high energy barrier that prevents the thermal flip of the magnetization of a single domain. Given the dimensions of the nanoarrays (∼ 100 nm), this energy is of the order of 10 4 -10 5 K. This precludes the thermal equilibration of the sample at room temperature. Indeed, it is the presence of mobile monopoles (created by the flip of domains) that allows the system to scan the ice-rule manifold. Heating the system in order to recover some dynamics would not help, because the material would melt before reaching the necessary temperature to induce fluctuations in the domain orientations.

To circumvent this problem, complex protocols involving varying external magnetic fields have been designed [START_REF] Wang | Artificial 'spin ice' in a geometrically frustrated lattice of nanoscale ferromagnetic islands[END_REF]. These protocols have efficiently allowed to approach the ground state, but it remains difficult to observe the finite temperature properties of the systems. Nonetheless, it has been recently proven experimentally [START_REF] Kapaklis | Melting artificial spin ice[END_REF], that the thermal-fluctuation regime is accessible in "soft" artificial spin ice. In a series of recent works soft magnetic alloys were used to design magnetic domains with tunable Curie temperature T C . In this system, it is possible to set T C below the room temperature, hence enabling thermal fluctuations [START_REF] Morgan | Thermal ground-state ordering and elementary excitations in artificial magnetic square ice[END_REF][START_REF] Porro | Exploring thermally induced states in square artificial spin-ice arrays[END_REF].

Very recently, magnetic fluctuations have been observed in real time in small hexagonal lattice clusters [START_REF] Farhan | Exploring hyper-cubic energy landscapes in thermally active finite artificial spin-ice systems[END_REF], and thermalization processes have allowed to reach the ground state of artificial square ice [START_REF] Zhang | Crystallites of magnetic charges in artificial spin ice[END_REF]. These works suggest that thermodynamic properties of square ice might soon be measured in artificial systems.

Other implementations

The quest for 2D artificial ice has reached other domains of physics as well. In particular, colloids have been used to simulate the antiferromagnetic Ising model on the triangular lattice [START_REF] Han | Geometric frustration in buckled colloidal monolayers[END_REF]. In this work, the Ising degrees of freedom were encoded in the position of a colloidal particle, similarly to the position of the proton on the O-O bond in water ice. A possible implementation of ice physics has been proposed in the context of optically trapped colloids [START_REF] Libál | Realizing colloidal artificial ice on arrays of optical traps[END_REF], based on recent experimental progress in the design of the traps and the tuning of the interactions [8]. This proposal consists of an array of double-well optical traps containing each a singly charged colloid. The colloids interact with each other through (screened) Coulombic interaction, their positions in each trap encoding the Ising degrees of freedom of the ice model. Here the height of the barrier inside a trap sets the temperature below which the system freezes. It is typically of order 1 -10 K, allowing for thermal fluctuations at room temperature.

The idea of realizing ice physics via charged particles in an array of double-well traps can be scaled down to the atomic size by using trapped ions. The Ising degree of freedom can be again encoded either in the position of the ion in a double-well trap [START_REF] Retzker | Double well potentials and quantum phase transitions in ion traps[END_REF][START_REF] Schmied | Quantum phases of trapped ions in an optical lattice[END_REF]. Alternatively, the internal state of the ion can serve to encode this Ising spin, with spin-spin couplings induced by lasers [START_REF] Johanning | Quantum simulations with cold trapped ions[END_REF][START_REF] Clark | A two-dimensional lattice ion trap for quantum simulation[END_REF][START_REF] Kim | Quantum simulation of frustrated ising spins with trapped ions[END_REF][START_REF] Bermudez | Frustrated quantum spin models with cold coulomb crystals[END_REF][START_REF] Korenblit | Quantum simulation of spin models on an arbitrary lattice with trapped ions[END_REF][START_REF] Schneider | Experimental quantum simulations of many-body physics with trapped ions[END_REF]. In particular, the interactions between the ions can be tuned in such systems, allowing for the simulation of a large variety of models [START_REF] Golovach | Dynamic response of one-dimensional bosons in a trap[END_REF]. More importantly, quantum effects can be tuned in these systems, therefore they are designated candidates to study the physics of quantum (spin) ice [START_REF] Molavian | Dynamically induced frustration as a route to a quantum spin ice state in tb 2 ti 2 o 7 via virtual crystal field excitations and quantum many-body effects[END_REF], i.e. the behavior of spin ice systems in presence of quantum fluctuations.

Spin liquids

As discussed in the previous chapter, Classical frustrated magnets that have a large ground state degeneracy and an extensive residual entropy seem to violate the third law of thermodynamics. Nonetheless the extensive degeneracy requires fine tuning of the Hamiltonian and, in experimental systems, the ground-state degeneracy will be lifted by the inevitable presence of disorder and/or interactions not accounted for in the model Hamiltonian. However, from a purely theoretical point of view the third law of thermodynamics can be restored by quantum mechanics. The resulting finitely degenerate ground state may break some symmetries of the Hamiltonian via a quantum order-by-disorder mechanism, or it may remain as symmetric as the classical parent ground-state manifold and hence realize a quantum spin liquid [START_REF] Balents | Spin liquids in frustrated magnets[END_REF].

Valence-bond states

In standard unfrustrated magnetic systems in dimension d > 1, the ground state exhibits long-range magnetic order, and hence breaks the spin rotational or inversion symmetry and in some cases even the translational symmetry. The situation is richer in frustrated antiferromagnetic systems. In this case, the ground state can be disordered and respect all the symmetries of the Hamiltonian. This kind of state was first conjectured by Anderson [2], when looking for the ground state of the antiferromagnetic S = 1/2 Heisenberg model on the triangular lattice. A convenient way to represent a state that does not break the global SU (2) symmetry is in terms of singlets (or valence bonds). A singlet ground state can always be exactly written as a superposition of possible singlets coverings of the lattice [START_REF] Liang | Some new variational resonating-valence-bond-type wave functions for the spin-1 2 antiferromagnetic heisenberg model on a square lattice[END_REF]. The idea of reducing the Hilbert space of the spin system to that of singlet coverings of the lattice was later taken up by Rokhsar and Kivelson [START_REF] Rokhsar | Superconductivity and the quantum hard-core dimer gas[END_REF], launching a whole line of investigation on the so-called dimer models [START_REF] Fisher | On dimer solution of planar Ising models[END_REF][START_REF] Nagle | Exact-solutions for a variety of dimer models[END_REF][START_REF] Raman | Quantum dimer models and exotic orders[END_REF].

Depending on the system considered, the valence-bond ground state can break or preserve the translational symmetry. In the former case, the translationally ordered ground state is called a valence-bond solid (VBS). The VBS breaks a discrete symmetry, and therefore does not have gapless Goldstone modes, but rather features a gap in its excitation spectrum. In this case, the excitations have integer spin and do not fractionalize.

If the ground state preserves the translational symmetry, then it is known as a resonating valence-bond (RVB) state. It was first introduced by Anderson in the triangular antiferromagnetic Heisenberg model [2], an it initiated the concept of quantum spin liquid. In this case the ground state is a superposition of a large number of short-range singlet coverings of the lattice (see figure 1.10). This sate has very exotic properties, common to all spin liquids, some of which we now briefly present.

+ + + ...

Features of a spin liquid

In the context of valence-bond states, the spin liquid is easily defined as the alternative to a solid state. Intuitively, spin liquids are thought as states with no long-range order. Indeed, symmetry-wise the spin-liquid phase is continuously connected to the paramagnetic phase. However, fundamental topological properties distinguish the two states. The general definition of a quantum spin liquid involves subtler concepts, describing the quantum nature of these phases. In particular, the integer-spin excitations of a spin liquid can generally decay rapidly into spin-1/2 quasi-particles, known as spinons, that may be deconfined, in the sense that the charges of the underlying lattice gauge theory are deconfined. This is the case for the RVB state, and it can be used as a characterization of a spin liquid. Anyhow, as far ad the scope of our discussion is concerned, spin liquids can be thought as phases with an underlying emergent gauge theory.

Lattice gauge theories 1.4.1 Comparison with spin ice systems

The link between spin liquids and spin ice is rather natural. Indeed, as discussed in section 1.2.4, the Coulomb phase does not break any symmetry and it supports fractionalized monopole excitations. These excitations are a trace of an underlying classical gauge theory, emerging from the ice-rule constraint, as discussed in section 1.2.4. When considering quantum spin ice, one might then conjecture the existence of a quantum spin liquid phase supporting fractionalized excitations if the quantum version of the underlying lattice gauge theory admits a deconfined phase. In the following we will discuss the relevant lattice gauge theory for classical and quantum spin ice, which are of the U(1) type.

Compact lattice quantum electrodynamics

Spin ice models are defined by the local constraint that the magnetization must be divergenceless. This stems from the ice rules and has fundamental consequences on the low-energy properties. In particular, quantum spin ice can be rewritten as a lattice gauge theory as we will explore in details in chapter 6. The topic of lattice gauge theory was introduced by Wegner [START_REF] Wegner | Duality in generalized ising models and phase transitions without local order parameters[END_REF] who constructed an Ising model with a gauge (=local) symmetry, and found a phase transition "without an order parameter". We will here briefly describe the construction and the properties of the lattice gauge theory which is relevant for a given limit of quantum spin ice in two dimensions, namely the (2 + 1)-dimensional U(1) compact lattice gauge theory or equivalently the two-dimensional compact lattice quantum electrodynamics (cQED) and show that they lead to confinement and to massive photons, and constrast this to the three-dimensional case where a deconfined phase is instead found. The interested reader is referred to Kogut's review [START_REF] Kogut | An introduction to lattice gauge theory and spin systems[END_REF], or Polyakov's book [START_REF] Polyakov | Gauge fields and strings[END_REF] for more details.

Compact U(1) lattice gauge theory

x âx (x) Êy (x) × p • • • • • • • • • • • • Figure 1.11: Compact QED
fields. The operators âµ (x) and ʵ (x) act on the link (x,µ) of a square lattice.

In the case of d-dimensional classical electrodynamics, the action in continuous space is

S = 1 4µ 0 d d x (∂ µ A ν -∂ ν A µ ) 2 (1.13)
where A is the vector potential. There are at least two lattice gauge theories that lead to this action in the continuum limit. The most natural one is, a priori,

S = 1 4µ 0 x,µ,ν F 2 x,µν F x,µν = A x,ν + A x+µ,ν -A x+ν,µ -A x,µ -∞≤A x,µ ≤ +∞ (1.14)
where the vector potential A x,µ lives on the links of a d-dimensional hypercubic lattice while F x,µν lives on the center of the plaquettes. But one can also choose the compact prescription

S = 1 2µ 0 x,µ,ν (1 -cos F x,µν ) -π ≤ A x,µ ≤ +π. (1.15)
It describes a compact gauge field compact because the values of A µ (r) are defined modulo 2π. Both models yield the same limit, but the last one (compact U(1) lattice gauge theory) naturally produces the charge quantization [START_REF] Polyakov | Gauge fields and strings[END_REF], whereas the latter has to be imposed in the first version (non-compact U(1) lattice gauge theory). Furthermore, compact U(1) lattice gauge theory turns out to be the relevant lattice gauge theory for quantum spin ice. Let us now describe some of its properties. Consider the partition function Z defined by the compact action (1.15). It turns out that it can be factorized into two parts

Z≡ π -π x,µ dA x,µ exp        - 1 2µ 0 x,µ,ν (1 -cos F x,µν )        = Z SW Z c (1.16)
where

Z SW = π -π x,µ dA x,µ exp        - 1 4µ 0 x,µ,ν F 2 x,µν        (1.17)
is the partition function of a non-compact U(1) lattice gauge theory, and

Z c = π -π x,µ dA x,µ exp        - π 2µ 0 r,r q(r)q(r ) |r -r | -k r q 2 (r)        (1.18)
where the q(r) are a set of gauge charges interacting via a Coulombic interaction, sitting at the center of the cubes. (they are the instanton solutions of this theory). As the Coulombic interaction is not confined in three dimension, these particles are unbound at all temperature and the theory describes a Coulomb gas of deconfined charged particles.

In two dimensions, the Coulombic interaction between the gauge charges becomes a logarithmic interaction (indeed, this interaction comes from the Green function of the Laplacian), which is confining. Therefore, in the ground states the gauge charges are confined in this case, and can only be deconfined at finite temperature.

Compact lattice quantum electrodynamics

We now turn to the (2 + 1)-dimensional compact lattice quantum electrodynamics (cQED). Its construction is very similar to that of the U(1) lattice gauge theory. Indeed, the vector potential A x,µ of (1.15) has to be replaced by the eigenvalue of an operator Âx,µ living on the links of the hypercubic lattice. We introduce then its canonical conjugate Êx,µ (the electric field) satisfying

µ (x), Êν (x ) = i δ µν δ xx (1.19)
We define then the lattice curl Φp of A as its flux around plaquette p

Φp = Ây (x + ŷ) -Ây (x) - Âx (x + x) -Âx (x) = Ây (x + ŷ) -Ây (r) -Âx (x + x) + Âx (x) (1.20)
whose eigenvalues corresponds to the F x,µν of (1.15). Its action on an eigenstate e µ (x) of ʵ (x) is

e ±i Φp e µ (x) = e µ (x) ± 1 . (1.21)
This quantization of the electric field comes from the compactness of the theory. The action of the cQED is very similar to that of (1.15), but with a spatially anisotropic coupling constant [START_REF] Kogut | An introduction to lattice gauge theory and spin systems[END_REF]:

S = 1 2µ 0 x,ν (1 -cos F x,0ν ) -t x,i,ν cos F x,iν . (1.22)
where the i and ν are respectively summed over the spatial directions and all space-time directions. This leads to the Hamiltonian

H = 1 2µ 0 x Ê(x) 2 -t p cos Φp . (1.23)
As a consequence, the ground state of this theory is a vacuum of field (e µ (x ≡ 0). In the low energy limit, the cosine can be expanded, and the Hamiltonian becomes quadratic in Φp , which is the lattice equivalent of the magnetic field, so that we recover the standard energy of an electromagnetic field.

Frustrated compact quantum electrodynamics from square ice

Starting from the square ice model, one can construct a frustrated compact QED (fcQED) [START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF]. In the compact lattice QED, the electric field is quantized. This quantization was obtained from the compactness of the theory. Here, the construction is the opposite. The electric field of the fcQED will correspond to the discrete magnetization M i = ±1/2 of the square ice model, introduced in section 1.2.4, which is quantized by construction. The electric field is defined so that its eigenvalues satisfy e µ (x) = M µ (x) + 1/2 = ±0, 1, and e -µ (x + μ) = -e µ (x). In this language, the ice rules reads

∆ µ e µ (x) ≡ e x (x + x) -e x (x) + e y (x + ŷ) -e y (x) (1.24)
This is Gauss's law, for the discrete divergence ∆. The electric field ʵ (x) = Ŝ µ (x) + 1/2 = ±1 and the vector potential Âx,µ are then defined exactly like in the cQED. We can then enlarge the electric field Hilbert space, allowing for any integer value of e µ (x) instead of just 0 or 1, but the other values will be energetically penalized by a term of the form

H 0 = 1 2g x,µ ʵ (x) -1/2 2 .
(1.25)

The path-integral formulation of this model yield a Hamiltonian similar to that of the unfrustrated cQED

H = 1 2g x,µ ʵ (x) -1/2 2 -t p cos Φp . (1.26)
The frustration comes from the 1/2 term in the electric field part of the Hamiltonian, that forces the electric field to take values 0 or 1. Without this term, this is exactly the compact QED. Ref. [START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF] concludes that the same Polyakov mechanism responsible for the confinement in ordinary cQED is also at play in its frustrated version. The gauge charges can only become deconfined at finite temperatures. Similarly to the cQED, three dimensional fcQED is deconfined for small quantum fluctuations [START_REF] Hermele | Pyrochlore photons: The u(1) spin liquid in a s = 1 2 three-dimensional frustrated magnet[END_REF].

Spin liquids in spin ice

As already mentioned, the ground-state manifold of spin ice systems can be identified with a "classical" spin liquid. Observing quantum effects in spin ice materials is a current challenge. Indeed, in Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 that the gap between the crystal field ground state doublet (i.e. the effective Ising states) and the first excited states is of order 300 K [START_REF] Gingras | Quantum spin ice: A search for gapless quantum spin liquids in pyrochlore magnets[END_REF]), and therefore they admit a classical treatment. However, in some spin ice materials, the this gap is smaller and quantum effects are more important. For instance, it has been suggested that Er 2 Ti 2 O 7 should order at low temperature by a quantum order-by-disorder mechanism [START_REF] Zhitomirsky | Quantum order by disorder and accidental soft mode in er 2 ti 2 o 7[END_REF][START_REF] Savary | Order by quantum disorder in er 2 ti 2 o 7[END_REF]. The first material for which the name quantum spin ice was used [START_REF] Molavian | Dynamically induced frustration as a route to a quantum spin ice state in tb 2 ti 2 o 7 via virtual crystal field excitations and quantum many-body effects[END_REF] was Tb 2 Ti 2 O 7 , in which spin-ice like correlations as well as transverse fluctuations might have been observed [START_REF] Gingras | Quantum spin ice: A search for gapless quantum spin liquids in pyrochlore magnets[END_REF]. Besides, recent experimental and theoretical work on Yb 2 Ti 2 O 7 has proven that quantum effects are important, even though the ground state is supposed to be magnetically ordered [7].

Different theories have been recently developed to describe these new quantum spin ice materials. In [START_REF] Benton | Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice[END_REF], Benton and collaborators developed a lattice QED theory to describe the physics of quantum pyrochlore spin ice. In particular, they found a "ghostly" gapless photon mode, with a vanishing density of states as ω → 0 (see figure 1.12). Furthermore, they calculated the structure factor within their theory. They found that quantum fluctuations have a tendency to reduce the height of the pinch points at T = 0 and hence reduce the topological sector fluctuations. However, at finite temperature the pinch points are partially recovered (figure 1.13).

An alternative approach to quantum spin ice, the gauge mean field theory (gMFT), has been recently introduced [START_REF] Savary | Coulombic quantum liquids in spin-1/2 pyrochlores[END_REF][START_REF] Savary | Spin liquid regimes at nonzero temperature in quantum spin ice[END_REF] in the study of S = 1/2 pyrochlore spin ices, and it will be discussed in chapter [START_REF] Benton | Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice[END_REF] 6 in the context of two-dimensional spin ice. In particular gMFT allows to take into account the spinon excitations (namely the gauge charges). For pyrochlore spin ice, this theory constructs a rich phase diagram with various exotic phases, including a quantum spin liquid phase.

Plan of the thesis

In this thesis, we will explore different versions of square ice. On the one hand we study the effects of dipolar interactions on the classical Coulomb phase. In particular we are interested in determining the optimal experimental parameters to observe the Coulomb phase in artificial spin ice, be it realized in magnetic nanoarrays or trapped particles. We characterize the different phases and the ordering transition induced by the dipolar interactions. On the other hand we will examine the effects of quantum fluctuations introduced by a transverse field on an Ising representation of square ice. In particular, making use of a novel Monte Carlo scheme, we determine the non-perturbative phase diagram of quantum square ice, and link it to the results of spin-wave theory, of degenerate perturbation theory and of gauge mean-field theory.

The structure of this manuscript is as follows. The first part is concerned with the classical dipolar square ice. In the first chapter we introduce the model and discuss the spectrum of its normal modes. The second chapter deals with the principles of Monte Carlo simulation, and in particular it describes the update algorithms we used. In the third chapter the results are detailed and the phase diagram of dipolar spin ice is given, including the analysis of the different phase transitions.

The second part is a study of quantum square ice, described as a transverse-field Ising model (TFIM) on the checkerboard lattice. The first chapter is a classical and semi-classical analysis of the model, using linear-spin wave theory. The second chapter applies perturbation theory and gauge mean-field theory to the study of the low-field properties of the model. The third chapter describes a new quantum Monte Carlo (QMC) algorithm that was designed to study both the quantum square ice and a related model, the quantum link model. Finally, the last chapter is the discussion of the quantum Monte Carlo results. In particular, the phase diagram of the TFIM obtained from the QMC simulations is compared to the conclusion of the effective gauge theories derived in the second chapter.

Part I

Spin ice physics with trapped charged particles

Chapter 2

Dipolar square ice

Experimental motivations

Inspired by experiments and proposals of artificial ices [START_REF] Wang | Artificial 'spin ice' in a geometrically frustrated lattice of nanoscale ferromagnetic islands[END_REF][START_REF] Libál | Realizing colloidal artificial ice on arrays of optical traps[END_REF], we study the phase diagram of square ice with dipolar interactions.

The dipolar nanoarrays realization of artificial ice are approaching the physics of ice. However, in all experiments realized up to now, the ice-rule manifold degeneracy was strongly lifted. We consider the possibility of modifying the geometry of these systems in order to reduce this energy splitting and recover, at least partially, the physics of the degenerate square ice.

ε + δ ε

Figure 2.1: Energy splitting of the ice vertices in the planar geometry. The vertices of the "current" kind (left) are more energetic than the vertices of the "vortex" kind (right). The energy difference δ between the two kinds of ice vertices sets the ordering temperature of the system.

We consider a realization of the square ice model consisting of elongated dipoles with point-like opposite charges on its ends, siting along the links of a square lattice. These dipoles interact via the Coulombic interactions of the charges, which becomes a dipolar interaction at long distance. This model is particularly relevant for a realization of square ice using trapped ions or colloids. However, the essential physics of dipolar nanoarrays artificial ice should captured by this model. In any strictly two-dimensional realization of artificial spin ice, arranging the optical traps or the nanoarrays within a plane leads to a different in-vertex interaction between parallel dipoles and between perpendicular dipoles. This lifts the degeneracy of the ice-rule vertices, and hence also the exponential degeneracy of the square ice ground state. The question is then to find a work around to recover ice physics in an artificial ice set up, and to determine what are the optimal parameters of such a system and how sensitive they are to imperfect tuning.

The model

In order to recover the degeneracy and the interesting physics of ice, it has been suggested [START_REF] Moller | Artificial square ice and related dipolar nanoarrays[END_REF] to introduce a height displacement h between the two sublattices of the system. The vertical spins of the vertex model are then shifted to a height h above the horizontal spins (see figure 2.2). Given the lattice spacing 2a and the trap half-width λ, tuning the value of h then allows to recover the energy degeneracy of the ice-rule vertices.

h 2a λ J 1 J 2 Figure 2
.2: Array of double-well traps. The system is a set of double well traps of half width λ. There is a height shift h between the "horizontal" and the "vertical" traps. The system forms a square lattice of lattice spacing 2a with the traps on its links. Here we represent a possible vertex configuration that obeys the ice rules. On the left the "actual" system, on the right its representation in the vertex model. The position of the colloid in the trap encodes the pseudospin degree of freedom. We note J 1 (blue dashed lines) the interaction between orthogonal pseudospins and J 2 (red solid lines) the interaction between pseudospins.

Mapping onto an Ising model

The interaction between a given pair of particles only depends on the position of each particle in its trap. Having two possibilities for each particle, this gives 4 possible energies. The pairwise interaction can then be cast into a 2 × 2 Hamiltonian coupling the occupancies of the wells of a trap with those of the other. Let us note n i,p = 0, 1 the occupation of the well p of site i, and V pp ij the interaction energy between two particles on sites (i, p) and (i, p ). The Hamiltonian of the system can then be written as

H = 1 2 ij pp V pp ij n i,p n j,p . (2.1)
This in turn can be rewritten as an Ising Hamiltonian (i.e. a quadratic Hamiltonian coupling ±1 degrees of freedom). We introduce the pseudospin on site i, given by σ i = n i,2n i,1 . The constraint of having exactly one particle per trap (i.e. n i,2 + n i,1 = 1) gives

n i,1 = 1 -σ i 2 n i,2 = 1 + σ i 2 (2.2)
Casting this into the Hamitonian (2.1) we get a quadratic spin Hamiltonian

H = 1 2 ij J ij σ i σ j J ij = 1 4 pp (-1) p+p V pp ij . (2.3)
Indeed, the linear term disappears, up to a boundary term. It is given by

1 8 i, j p,p V pp ij (σ i -σ j ) = 1 8 i, j         p,p V pp ij - p,p V pp ji         σ i ≡ 0 . (2.4)
The last equation comes form the symmetry

V pp ij = V p p
ji . This result was expected, as a linear term would act as a magnetic field and break the symmetry of the original model. There is also a constant term, that we drop because it only shifts the origin of the energies.

Neutral system

In the case of trapped colloids, the Coulombic interactions are in fact screened, leading to an effective Yukawa interaction. However, inspired by the nanoarray systems and the possibility to extend the colloid proposal to trapped ions, we choose unscreened interactions.

Neutral system

We consider an ensemble of trapped charged particles (from now on referred to as "ions"). Each ion sits in one of the minima of a double-well trap. The physics of the system is unchanged when introducing a neutralizing background. As the interaction between the ions is Coulombic, in presence of this background, the Hamiltonian (2.1) becomes

H = 1 2 ij pp V pp ij (n i,p -1/2)(n j,p -1/2) ≡ 1 8 ij pp V pp ij ρ i,p ρ j,p (2.5)
where we have introduced the charge ρ i,p = 2n i,p -1 = ±1 at {p, i}. By construction, each trap is charge neutral (i.e. ρ i,2 + ρ i,1 = 0), so we get, with a pseudospin defined as

σ i = ρ i,2 -ρ i,1 /2: ρ i,1 = -σ i ρ i,2 = σ i . (2.6)
The Hamitonian is then again a quadratic spin Hamiltonian

H = 1 2 ij J ij σ i σ j J ij = 1 16 pp V pp ij (-1) p+p . (2.7)
Here, the interaction between the charges is simply a Coulombic potential

V pp ij = V 0 r i,p -r jp . (2.8)
This leads to a dipolar interaction between the pseudospins. Indeed, let us consider two double-traps centered at r i and r j . The charges are located at r i ± λd i /2 and r j ± λd j /2, with d i , d j ∈{ x, ŷ} depending on the orientation of the traps. We note r ij = r ir j . The interaction between the corresponding pseudospins writes then

J ij = V 0 16 α=±1 α =±1 (-1) p+p r ij + λ(αd i -α d j )
.

(2.9)

At long distance r ij = r ij λ, this becomes a dipolar interaction

J ij = V 0 16r ij α=±1 α =±1 (-1) p+p r 2 ij + 2λr ij • (αd i -α d j ) + λ 2 (αd i -α d j ) 2 = - V 0 λ 2 16r 3 ij         d i • d j -3 (r ij • d i )(r ij • d j ) r 2 ij         + O         λ 3 r 3 ij         , (2.10)

Nearest neighbors interaction : the restoration of the ice states degeneracy

Let us now restrict ourselves to the interactions between pseudospins of the same vertex, neglecting any further coupling. We do not need any detail about the interactions at this point. We just note J 1 and J 2 the interactions between orthogonal and collinear pseudospins respectively (c.f. figure 2.2). The value of J 1 depends on the height offset h. Hence tuning h amounts to tuning the ratio J 2 /J 1 .

There is then a particular value of the height h c (λ) for which J 1 = J 2 , recovering the full degeneracy of the ice-rule state. For any other value of h, the degeneracy of the ground state is reduced (see figure 2.3). This reduction is different whether h is larger or smaller than h c (λ).

h < h c (λ).
In this case, the "vortex" vertices are favored. There are only two ways to arrange a square lattice of such lattices, reducing the degeneracy to a finite value. The two ground states are the two staggered arrangements, later on referred to as "vortex states".

h > h c (λ).
In the case of a large offset, the vertices of "current" type are favored (see figure 2.1), i.e. vertices with a net magnetization. The ground state degeneracy remains divergent with the system size. Indeed, any state with ferromagnetically ordered vertical and horizontal lines minimizes the energy. For a system of L × L vertices, there are 2 2L such "current states" (2 for each of the lines). The entropy is then subextensive and we do not expect any ordering at low temperature.

This phase can be easily understood in the limit h →∞ . In this case, the coupling J 1 between the two sublattices vanishes. The system then splits into two sets of independent spin chains with ferromagnetic interactions. The ratio between the two prominent interactions (J 2 ,J 1 ) varies smoothly from 0.48 to 1.94 as h is increased from 0 to 1.2 a. For J 2 > 2J 1 , the excited 3 in-1 out vertices become less energetic than the vortex vertices. The system is then far from the ideal square ice description.

h J 2 /J 1 0 h c 1.2 1.5 1 0.5 0 h < h c (λ) Vortex → 2 states h > h c (λ) Current → 2 2L

Objectives

Considering this kind of interaction, we are interested in finding the range of temperature in which the ideal square ice features can be seen. As mentioned in section 2.2.3, the energy splitting of the ice states sets a temperature scale T * above which the degeneracy is effectively recovered. If this temperature is lower than the energy cost of a defect ∆, placing the system at an intermediate temperature T * < T < ∆ should allow to observe ice physics in a non ideal system. Indeed, it has been recently pointed out that the survival of the Coulomb phase features in dipolar 3D spin ice was due to thermal effects [START_REF] Sen | Coulomb phase diagnostics as a function of temperature, interaction range, and disorder[END_REF]. The value of T * comes from two kinds of deviations from ideal ice.

First, if h is not perfectly tuned, the degeneracy of the single-vertex energies is lifted. The induced splitting energy δ is included in T * .

Secondly, even if h is tuned so that the six ice-rule vertices are degenerate, the contribution tails of the interaction is not the same for every ice-rule sates. Therefore the longer-range interactions introduce another energy scale to include in T * .

The question is then to determine the range of temperature ∆T ≈ [T * , ∆] in which our system behavior mimics that of an ideal square ice. In particular, we want to check that there is a significative range of values for h for which ∆T T * , allowing for imperfect tuning and hence possible experimental realizations.

Expected ground states

Following the analysis of the short-range model, depending on the height offset we expect two kinds of ground states. The results are not changed by the long range interactions in the case h < h c . The ground state is still the two-fold degenerate vortex state. However, if h > h c , the degeneracy among the 2 2L current states is lifted. The ground states are selected by the third-neighbour interaction, between collinear spins belonging to adjacent vertices (see figure 2.4). Taking a dipolar approximation for J ij , we get J 2 = -2J 3 . This interaction is antiferromagnetic. It will then favor the current states with alternating orientations of the horizontal lines and of the vertical lines, like the traffic flow in Manhattan. Whence their name : the "Manhattan" states [48]. There are four of them. As described in chapter 4, this four-fold degeneracy has interesting implications on the ordering transition. J 3 Figure 2.4: Including the third-neighbor interaction J 3 is crucial to predict the low temperature phase in the h > h c (λ) case. This interaction is antiferromagnetic. The ground state it selects is the Manhattan state, which is four-fold degenerate. The magnetic unit cell (dotted red square) contains four geometric unit cells (dashed blue square), and hence 8 pseudospins.

Normal modes and bandwidth

A more quantitative way to search for an ordered phase is to look at the Fourier transform of the interactions. It transforms the quadratic Ising Hamiltonian in a set of harmonic oscillators known as the normal modes of the model. These modes are labelled by their wave vector q and have a frequency (an energy) ω q . The presence of a minimum in the spectrum ω q indicates an ordered state of lowest energy with the ordering wave vector q. Setting the system in the corresponding state will minimize its energy. Typical examples are the antiferromagnetic Ising and Heisenberg models on the square lattice. They share the same normal modes with a minimum at the (±π, ±π) points (figure 2.5), which corresponds to their ground state : the Néel state. There are minima at the equivalent (±π, ±π) points with energy -2J, corresponding to the Néel states which is the ground state of the model. The maximum at (0, 0) corresponds to the ferromagnetic state which is the highest energy state.

It should be noted that this approach only accurately describes systems with continuous variables. Indeed, let us consider the Fourier states σ q with energy ω q . The normal modes wave vectors lie in the first Brillouin zone of the reciprocal lattice. For the sake of simplicity, we will suppose that the interaction has the same symmetries as the lattice and that this lattice is a Bravais lattice. Under these assumptions, two spins σ i and σ j separated by a distance R j -R i =∆ R are related by σ j = e iq•∆R σ i . Hence for an Ising system q • ∆R/π has to be an integer. On a square lattice, this means that only Fourier states satisfying q = π(m x , m y ) with m x and m y integers preserve the symmetries of the interactions. If the lattice is decorated, the same argument holds with q half-integer multiple of the corresponding reciprocal lattice primitive vectors.

The lattice of the square ice model is not a Bravais lattice. Both for the model with short-range (J 1 -J 2 ) interactions only, and for the model with long-range interactions, the unit cell contains two spins (one vertical and one horizontal, see figure 2.4). Let us rewrite the Hamiltonian in a form that makes the Bravais lattice explicit. We note R, R the sites of the Bravais lattice, and denote by a, b the two sublattices. We get

H = 1 2 ij J ij σ i σ j = 1 2 RR a,b J a,b (R, R )σ a (R)σ b (R ) (2.11)
where J is the 2 × 2 matrix containing the interaction between the two spins of each site R and R . The periodic boundary conditions impose that J only depends on the distance ∆R = R -R between the sites. For each wave vector q there are then 2 normal modes α 1,2 q of energy ω 1,2 q . These are the eigenvectors and eigenvalues of the Fourier transform J q of J(∆R). The Hamiltonian then reads

H = q α,β J a,b q σ a q σ b -q = q,s ω s q α s q 2 .
(2.12)

In the case of the ideal square ice, the lower band is perfectly flat (figure 2.7). This flatness is directly related to the extensive degeneracy of the ice states manifold, and reflects the Fourier space representation of the zero-energy excitation modes in real space, which are flips of closed loops. First identified by Anderson [5] in the pyrochlore lattice, these excitations allow the system to connect the ice states with each other.

In presence of longer range interactions, the lower band of excitations becomes corrugated (see figure 2.7). Two important informations can be deduced from this corrugations.

Bandwidth First, as inferred in the previous section, we expect any lift of the degeneracy of the ice states to induce ordering at a transition temperature which is roughly given by the magnitude of the energy splitting. This energy scale is expected to be of the same order of magnitude as the width δω of the lower band of the normal modes. Indeed, if the depth of a minimum q 0 of the spectrum is larger than the temperature, thermal fluctuations will not be sufficient to allow the system to escape the vicinity of q 0 , and explore the entire manifold of ice-rule states.

h = 0 . 2 a h = 0 . 6 a h = h c h = h c L R h = a h = 1 . 2 a R/a h =0.2 a h =0.6 a h =1.2 a h = a h = h c h = h LR c Figure 2
.6: Evolution of the width δω of the lower band of excitations with respect of the cutoff range R of the interactions, for various values of h. This width is always lowered in the long range limit. Its saturation value depends on h. Minimizing δω sets the optimal value of h. The smaller band width is obtained for h LR c ≈ 0.85a (see figure 2.8) which is slightly larger than h c .

As expected, introducing the long-range tails induces corrugations in the normal modes bands. In particular, the lower band, which is flat in the ideal square ice, acquires a finite width δω. Interestingly, if one introduces a cutoff of the interactions at distance R it turns out that, as R is increased, the value of δω first increases and then decreases, saturating to a smaller value (see figure 2.6). The lower band is then flatter with long range interactions than with truncated interactions. It seems then that the tails of the interactions partially screen the rest of the couplings. This means that one has to be very careful when truncating the interactions in such models. If the chosen range is too small, one may find artificial ordering due to a larger bandwidth than in reality. This phenomenon has already been observed in Ising pyrochlore magnets [START_REF] Gingras | Origin of spin-ice behavior in ising pyrochlore magnets with long-range dipole interactions: an insight from mean-field theory[END_REF], where it was used to explain the observation of spin ice physics in presence of long range interactions. This effect is present for all values of the offset h. Tuning h amounts to reducing the saturation value of this width. The smallest band width is obtained for a value h LR c ≈ 0.85a of the offset that is slightly larger than the short range value h c ≈ 0.805a.

Ground states As already mentioned, from the position of the minima in the band dispersion we can read out the ordering wave vector of the ground state. The so-obtained low temperature phases are the same as in a short-range model, and are consistent with the analysis of the first three interactions (see e.g. figure 2.4).

In figure 2.6 we represent this lower band width δω in function of the offset h for both the nearest neighbor (interactions limited to spins of the same vertex) and the long range model. Converting δω into a temperature gives us a sketch of the expected phase diagram of the model (figure 2.8). ). The range of the interaction is increased from top to bottom (nearest neighbor, 1, 2 and 64 lattice spacings). For h = h c , the lower band is perfectly flat when retaining nearest-neighbour interactions only. For all values of h, the width of the band reduces in the long range limit. In particular, away from h c the global form of the spectrum is recovered in the long range limit. The long range and nearest neighbor spectra are particularly close around the minima. For h = h LR c , the minima at the vortex state -(±π, ±π) -and the Manhattan state -(0, ±π) and (±π, 0) -have the same depth, signaling the degeneracy between them.

For h < h c , we find as expected a minimum at the (π, π) point in both short and long range limits. This is consistent with the Néel ordering of the unit cells in the "vortex" phase for which the spins of the two sublattices are given by σ I,1 ,σ I,2 = σ 0 (-1) I (1, 1). There are no significative difference between the nearest neighbor and the long range spectra.

For h > h c , the situation depends on the range of the interactions. In the short range limit, there are two lines of minima for q = (q x , 0) and q = (0, q y ). This corresponds to the 2 2L states with ferromagnetic horizontal and vertical lines (see figure 2.3). In the long range limit, the spectrum looks rather similar. But there is a significant modification of the two lines of minima. They are slightly bent to give minima at the (0, ±π) and (±π, 0) points and a local maximum at the origin (see figure 2.7). This maximum is present for all ranges, but it is significantly reduced in the long range limit. This changes promote some of the nearest neighbor ground states, namely the Manhattan state (see figure 2.4).

The bandwidth δω of the lower band of excitations is represented in figure 2.8. Because of the tails of the interaction, it never vanishes. However, it is greatly suppressed around h = h c . If the temperature is below δω, we expect the system to be ordered, whereas at higher temperature it should be disordered.

Based on this analysis, we now know that we can expect to find ice physics in our model at finite temperature, and we know which ordered phases we can expect at low temperature. Any further quantitative analysis will imply the use of numerical tools. In particular, the precise temperature transition and the type of transitions that occur cannot be obtained from this simple approach. That is why we turned to a Monte Carlo analysis of the system, as this method is particularly indicated to study the thermodynamics of lattice spin models. This is the topic of the next chapter.

Chapter 3

Monte Carlo

Monte Carlo is one of the most used numerical methods in statistical physics. The basic idea is to take advantage of randomness and importance sampling in guided random walks to evaluate rapidly thermodynamical averages. Indeed, even before the name Monte Carlo was used, people have had considered random sampling as a way to evaluate areas or, more generally, integrals. The conceptual step of transferring these idea to the study of physical systems was made more than half a century ago by Metropolis and Ulam [START_REF] Metropolis | The monte carlo method[END_REF]. Because of the relative simplicity of the concepts involved, it has been successfully adapted to a wide range of problems. This method is particularly powerful when dealing with lattice and/or spin systems. As we are interested in Ising models, this is the designated method. In this chapter, we will briefly describe the principles of Monte Carlo along with some of the most used algorithms, namely the Metropolis spin flip and the Wolff cluster updates. There are numerous books and reviews on the topics. So we will be very brief and invite the interested reader to consult for example the books by Newman and Barkema [START_REF] Newman | Monte Carlo Methods in Statistical Physics[END_REF] or by Krauth [START_REF] Krauth | Cluster monte carlo algorithms[END_REF][START_REF] Krauth | Introduction to monte carlo algorithms[END_REF] (his use of the concept of "a priori probability" is particularly clarifying). We will nonetheless give more details about the key algorithm for ice models which is the loop algorithms. The last section of the chapter slightly deviates from the Monte Carlo discussion, as it is focused on how one can extrapolate the finite size Monte Carlo results to the thermodynamic limit. To this end, we introduce the basic ideas of scaling in critical systems and describe its consequences on finite size scaling.

Principles

The Monte Carlo technique is dedicated to simulations of stochastic processes (i.e. whose evolution is non-deterministic). During a Monte Carlo simulation, the system explore different microstates i k , k = 1,...,M. The goal is to relate the average of a quantity on these explored states Ā to its thermodynamical average A.

Estimators and averages

Let us consider a system with N s accessible microstates i of respective weights w i (typically their Boltzmann weight). The thermodynamical average of a quantity A is then given by

A ≡ 1 Z N s i A(i) w i (3.1)
where Z = N s i w i is the partition function of the model. Note that every microstate has to be taken into account in this sum. It is clear that it will not be possible to count them all one by one because there are simply too many of them (consider for example the 2 N states of a system of N Ising spins). The idea of Monte Carlo calculations is to evaluate such a sum using only a small subset of the microstates of the system. The cost of this sampling is the introduction of a statistical error on the obtained results. In order to have full control on the output result of the algorithm, it is crucial to evaluate correctly this error. This is briefly addressed in the following section.

Estimators

If the M microstates are generated using the simple sampling method (i.e. with a uniform probability to obtain each state), the estimator of the desired average is simply

A M ≡ M k A(i k ) w i k M k w i k . (3.2)
This quantity clearly tends to A in the limit M→N s . However, sampling randomly some of these states will not help in general. Indeed, for MN s (which is almost always the case) it is most likely that the sampled states will not be representative and the calculated averages will be spurious. Fortunately, most of the states give a vanishingly small contribution. And one can then try to sample the configuration space in a clever way that makes sure that the system visit mainly states with a higher weight. Let then note p(i) the probability for the microstate i to be visited by the chosen sampling procedure. The estimator for A then has to be modified as

A M = M k A(i k ) w i k /p i k M k w i k /p i k . (3.3)
Choosing p i k ∝ w i k will greatly simplify the estimator so it becomes a simple average :

A M = 1 M M k A(i k ) (3.4)
This prescription for p i is the one chosen in Monte Carlo calculations. However, in order to adjust p i to the proper value w i one has to turn the random sampling into a Markov process (namely to implement the importance sampling). But how can one generates these states and, above all, how does one ensure that the attributed probability of visiting them is correct? This is controlled by two key ingredients of any Monte Carlo Algorithms that directly come from (3.4): the detailed balance condition and the ergodicity principle. They are described in subsection 3.2.2. Moreover, in such a Markov process, each configuration is generated from the previous one. Consecutive states might then be correlated, a fact which complicates the estimate of the statistical error. This is the purpose of the following discussion.

Sampling the correct distribution Master equation and detailed balance

We consider a system with a finite number of microstates i with weights w i that can a priori vary with time. Starting from a random initial state, we let the system span the configuration space, following a Markov chain defined by the transition probabilities P(i → j). These transition probabilities have to be chosen carefully because they drive the time evolution of the weights of the state through the so-called Master equation

dw i (t) dt = j P( j → i) w j (t) -P(i → j) w i (t) . (3.5)
With constant transition probabilities, this equation corresponds to a stationary Markovian process. In this limit, the transition rates only depend on the present configuration, the system has no memory of its previous steps. This implies that the system will eventually reach an equilibrium state for which the weights are constant over time (i.e. dw i (t)/ dt = 0, ∀i). We are interested in this thermal equilibrium regime. A way to obtain this from equation (3.5) is to set the transition probabilities such as

P( j → i) w j = P(i → j) w i ⇐⇒ P( j → i) P(i → j) = w j w i . (3.6)
This condition is at the heart of most Monte Carlo algorithms for statistical mechanics system at equilibrium. It is known as the detailed balance condition. We could solve (3.6) with some P( j → i) = P(i → j) = 0, but if w i 0 and w j 0 this would prevent in general to reconstruct the statistical sum in its entirety because some states with finite weight might become unreachable. One needs then to satisfy the ergodicity condition :

In order for the correct statistical sum to be reconstructed, any state with finite weight has to be reachable from any other state with finite weight in a finite time.

A priori probability

When studying the thermodynamics of a system by a Monte Carlo algorithm, one often has to design clever ways to change the configuration of the system. The choice of these "update moves" is what makes an algorithm efficient. Devising them demands some insight into the system and they are often very specific to a given problem. A way to make sure that an update move does not violate the detailed balance condition is to split the update probability into two parts

P(i → j) = p(i → j)A(i → j) (3.7)
where p(i → j) is the probability to propose state j as the arrival state of the transition starting from state i (the a priori probability) and A(i → j) is the probability to accept the update. One could say that the choice of A ensures the ergodicity while the choice of p guarantees the detailed balance. The detailed balance condition then becomes

p( j → i) p(i → j) = w j w i A(i → j) A( j → i) . (3.8)
In all the algorithms that were involved in this work, we set the return move probability to be the same as the onward move A(i → j) = A( j → i). For the acceptance probability, we chose Metropolis prescription [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] :

p( j → i) = min 1, w j w i .
(3.9)

Error estimate and autocorrelation time

If the microstates over which the sum (3.1) is sampled are fully uncorrelated, then the A(i k ) are identical independant random variables and the central limit theorem applies.The error δA on A M is estimated by using the time series as a sample of which random variable A i k , and hence the error on the sum of random variables is (δA

) 2 = 1 M A 2 M -A 2 M . (3.10)
However, in a Monte Carlo simulation, the successive states are often correlated. Let us now evaluate the error in such a case. In order to do so we consider n independant series

A µ (µ = 1,...,n) of M measurements A µ (t)(t = 1,...

,M). Each series of measurements yields an average value

A µ M
with a statistical error

δA µ = A µ M
-A (A being the thermodynamical average that would be obtained from an infinite series of measures). The reasoning is valid in the general case, but in order to remain as simple as possible, we will suppose that A = 0. We are interested in the variance ∆A 2 of the partial averages. We get

∆A 2 ≡ (δA µ ) 2 n =         1 M M t=1 A µ (t)         2 n = 1 M         A µ (t) 2 Mn + 2 M t=1 1 - t M A µ (0)A µ (t) n         . (3.11)
In the limit of n →∞ , the averages becomes thermodynamical averages and we will drop the µ index. We then introduce the autocorrelation function of A

φ A (t) = A(0)A(t) A 2 . (3.12)
This function decays to zero on a time scale τ A (one can define τ A more precisely as

τ A = t φ A (t)dt).
Neglecting the linear term t/M of (3.11), we get in the limit M→∞,

∆A 2 = 1 M A 2 (1 + 2τ A ) . (3.13)
Comparing with (3.10), it is clear that the error is larger for large correlation times. The reduction of this correlation time is then a crucial point in the optimization of a Monte Carlo algorithm. In practice, one produces a single time series. Knowing the value of τ A (or at least a rough estimate, i.e. using (3.12) from a preliminary series), the series of A(i k ) is divided into bins of length M b τ A . One computes then the series of averages A b over each bin b which therefore approximate statistically independent time series. The statistical error on the average over the whole series of measures is then

∆A ≈ 1 N b A 2 b N b -A b 2 N b (3.14)
where N b = M/M b is the number of bins. Note that this value is independent of the choice of the size of the bins, as long as it is much larger than the correlation time τ A . However, as we are dealing with finite series of values, we want the number of bins to be as large as possible, so that the standard error is most accurately estimated.

Standard algorithms

The update moves used in a Monte Carlo calculation can be either local or non-local. We used both in our calculations. Before turning to the ones specifically implemented for ice models, let us introduce two of the most paradigmatic algorithms : the Metropolis single spin flip algorithm 1 and the Wolff cluster algorithms. We will present them in the context of a ferromagnetic Ising spins system but they can be generalized to any spin model.

Single spin flip

The single spin flip update generates the local magnetization fluctuations due to the thermal bath. It proceeds as follows 1. Choose one of the N spins at random.

2.

Compute the energy cost ∆E to flip this spin.

Flip the spin with probability

p = exp(-β∆E) (i.e. generate a random number r ∈ [0, 1], if r < p,
flip the spin). β is here the inverse temperature.

The a priori probability is here A(i → j) ≡ 1/N. It is easy to check that this algorithm respects the detailed balance condition.

Wolff cluster algorithm

When the energy landscape in configuration space becomes too corrugated, and particularly in critical systems, it often becomes extremely unlikely to link two statistically relevant states using only successive local moves. In particular, the number of single moves needed to significantly rearrange the configuration varies at a critical point as power L z of the system size, and exponentially with the system size in an ordered phase. Because of this freezing, the correlation time also diverges with the system size. According to (3.13), this implies that the statistical error becomes so large that calculated averages are unreliable in this case.

To prevent the system from freezing in a particular configuration, one can use non local updates that directly link states differing by an extended update of spins. The cluster updates belong to this class of non local updates. There are many of them, but we will present the Wolff algorithm [START_REF] Wolff | Collective monte-carlo updating for spin systems[END_REF] which is the one we used. It consists in constructing sequentially a cluster of parallel spins and flipping it. Let us consider a system of N Ising spins with a ferromagnetic coupling J > 0. The algorithm is the following :

1. Let C = ∅ be the growing cluster, and S = ∅ the spins that may be added to C.

Randomly choose a spin i and add it to C.

3. Add to S all the nearest neighbors of i that are not already in C.

4. If S ∅, take the first element j of S, else go to (8).

5. Add j to C with probability max 0, 1 -e -2βJS i S j 6. Remove j from S.

7. If j was added to C, go to (3) with i = j, otherwise go to (4).

8. Flip the spins in the cluster C.

1 it must be the most famous update there is, as the paper presenting it [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] has more than 15000 citations as of today long loop short loop Figure 3.1: Loop moves. On the left, a long loop (in red) is constructed starting from the circled spin. This loop is winding in the y direction (it changes the vertical magnetization) but not in the x direction. On the right, the same loop is constructed in the short loop algorithm. Once the red loop is closed, the construction is stopped and all the previous spins (in green) are erased from the loop.

A few remarks are in order. First, note that adding a spin to the cluster does not change the energy of the system. It is not adding the spin that costs energy. It is a crucial point that one must keep in mind when generalizing this algorithm. This algorithm always updates the system. It might flip only the initial spin or the entire system, but it is important to note that the rejection probability is zero. However, for frustrated system this algorithm is not very efficient, as it has be shown to percolate, even for systems larger than the correlation length [START_REF] Leung | Percolation properties of the wolff clusters in planar triangular spin models[END_REF]. For such systems, other algorithms have been designed.

Loop algorithms

Resorting to cluster updates partially solve this problem by maintaining some fluctuations in the system when local updates are no longer sufficient. However, such algorithms as Wolff's cluster algorithm are blind to any local constraint that may exist in the system. As often in Monte Carlo, the most efficient algorithms are specifically designed for a particular model. Spin-ice systems are locally constrained to satisfy the ice rules at T = 0. A specific algorithm has been devised to sample efficiently the ice-rule manifold. This key idea of the algorithm is to construct closed loops of head-to-tail spins [START_REF] Evertz | Cluster algorithm for vertex models[END_REF][START_REF] Barkema | Monte carlo simulation of ice models[END_REF][START_REF] Yanagawa | Calculations of correlation functions for 2-dimensional square ice[END_REF][START_REF] Rahman | Proton distribution in ice and kirkwood correlation factor[END_REF], reversing one in and one out spin on each vertex it crosses, and hence preserving the ice rules. We used different variants of this algorithm, generalizing it to finite temperature and long-range interactions.

Long loop algorithm

The most straightforward version of the loop algorithm is what we will call the "long-loop" algorithm. Let us start with the case of the ideal square ice model with periodic boundaries. In this case, the six ice vertices have the same energy and the loop move is microcanonical.

Principles

We place ourselves at T = 0. This means that there are no defect in the system, i.e. all vertices obey the ice rules. The construction of the loop is just a "directed" random walk that eventually comes back to the starting point (see figure 3.1) :

1. Randomly choose a spin i and add it to the loop. Let v 0 be the vertex in which i is pointing out.

2. Consider the vertex v in which i is pointing in. If one of the two out spins of v has already been visited, add the other to the loop, otherwise pick one randomly (with even probability) and add it to the loop. Note j the chosen spin.

3. If j belongs to v 0 , add it to the loop and go to (4). Otherwise, go to (2) with i= j.

4. Flip all the spins in the loop.

Each spin is then visited at most once and each vertex at most twice. The energy of the configuration is not changed by this move. Hence the Metropolis prescription for the transition probability ensure that the move is always accepted. Note that for a loop to be able to change the net magnetization of the system, it has to be winding around the system. Otherwise, it makes as many steps to the left as to the right (and as many upward as downward), hence leaving the magnetization unchanged.

Scaling and topological sector fluctuations

This algorithm allows to sample the ice states efficiently. However, the length of the constructed loop diverges with the system linear dimension L of the system. More precisely, the average length of the generated loops scales as L 5/3 [START_REF] Barkema | Monte carlo simulation of ice models[END_REF][START_REF] Saleur | The antiferromagnetic potts model in two dimensions: Berker-kadanoff phase, antiferromagnetic transition, and the role of beraha numbers[END_REF]. In particular, a non negligible fraction of these loops will be winding around the system. The topological sectors, labeled by their winding number w (defined in section 1.2.3), can only be changed by the flip of such a winding loop. However, these loops take a long time to generate and slow down the simulation for sizable system sizes. An efficient alternative is available to sample the ice states within a given topological sector.

Short loop algorithm

The first ever implementations of a loop algorithm [START_REF] Yanagawa | Calculations of correlation functions for 2-dimensional square ice[END_REF][START_REF] Rahman | Proton distribution in ice and kirkwood correlation factor[END_REF] already used an alternative version of the algorithm. The idea is that, even though the length of the loops generated by the long-loop algorithm diverges with the system size, most loops intersect with themselves and can be subdivided in shorter loops (see figure 3.1). The algorithm is very similar to the long-loop algorithm. Instead of waiting for the loop to close at its starting point, we flip all the spins between two crossings of the same vertex, hence selecting a short-loop that can be flipped like any loop generated by the long-loop algorithm, as described in figure 3.1. This algorithm basically removes the tail before the first self-crossing of the growing loop. This algorithm is much faster than the previous one, as the average length of the generated loop does not scale with the system size (It is close to a constant value l sl ≈ 4.74). It allows to access any ice state of a given topological sector, as any non-winding loop can be decomposed into smaller non self-intersecting loops.

Generalization to finite temperature and longer-range interactions

At T = 0, a we chose our Monte Carlo steps for ice models to consist of a few long loop moves (so that the topological sector can fluctuate) and many more short loop moves. Basically, for a system of linear size L, we will take per Monte Carlo step L 2 /L 5/3 = L 1/3 long loops and L 2 /4 short loops so that statistically every spin should be visited at least once. For the generalized ice models, or at finite temperature, the algorithms have to be slightly modified.

Finite temperature

At finite temperature, defects are present in the system. With ∆ the energy cost of flipping a spin and thus creating a pair of monopoles, the density of monopoles is proportional to e -β∆ , with β the inverse temperature. We create such defects in the system with spin-flip updates (L 2 of them). The construction of loops as described previously only works with defect-free vertices (i.e. ice-rule vertices). Consequently, we add the following ingredients to the construction of the loops :

• The loop can only start on an ice-rule vertex. Otherwise the growth simply stops.

• If the loop reaches a 3in/1out vertex, the growth stops, but the loop is flipped nonetheless. This has the effect of "teleporting" the defect to the starting vertex. The detailed balance is ensured as this is a microcanonical move.

• If the loop reaches a 4in/0out or a 4out/0in vertex, the growth stops without any update of the configuration. Maintaining the loop updates at finite temperature allows to keep the ergodicity even when the temperature is small (compared to spin flip gap ∆=4J). At low temperature, the acceptance rate is 1 by construction. At finite temperature, the move can be rejected because of presence of defects. These defects proliferate as the temperature is increased, thus reducing the acceptance of these updates. The acceptance rate remains nonetheless reasonable in a wide range of temperatures (see figure 3.2). Note that this rate is independent of the size of the system. The "teleportation" move speeds up the -otherwise diffusivedynamics of the monopoles.

Similarly to the T = 0 case, the short loops length does not vary with system size (see figure 3.2). It is only slightly reduced by the presence of defects. The length of the long loops however becomes size independent down to temperature much smaller than the spin flip gap where the scaling with system size is only recovered (see its ratio by L 5/3 in the inset of lower left panel). Thus the flip of a winding loop becomes very rare, and the total magnetization of the system is not changed by this update anymore. In this regime, the fluctuations of magnetization are induced by single-spin dynamics.

Long range interactions

acceptance rate acceptance rate The loop algorithms can also be extended to the case of biased and/or longer range interactions [START_REF] Melko | Long-range order at low temperatures in dipolar spin ice[END_REF][START_REF] Melko | Monte carlo studies of the dipolar spin ice model[END_REF]. In the presence of interactions beyond nearest neighbors and/or a splitting in the energy of the ice-rule vertices, the loop move ceases to be microcanonical. Therefore the detailed balance is guaranteed by accepting a posteriori the flip of the loop with a probability given by the induced energy change ∆E : p = max 1, exp(-β∆E) . We might expect this to dramatically lower the acceptance rate of the algorithm.

l ll l sl T/∆ 0 T/∆ 0 2.8
Nonetheless, the acceptance rate of both loop algorithms remains significative in a sizable range of parameters (temperature and height offset, see figure 3.3) which correspond to a phase where ice physics is not suppressed either by the finite temperature or the long range interactions. The dramatical drop of this rate at low temperature corresponds to the transition to an ordered phase, where the energetic cost of the flip of a loop becomes too large. However, except for a small temperature range around T = 0.9J 2 (which corresponds to the ordering transition), the average length of the long loops is size independent. (Note that J 2 is independent of h and will therefore be used as the energy unit). Thus the fluctuations of the total magnetization are dominated by single-spin flips as the winding loops updates are suppressed for any sufficiently large systems. The suppression of the non-local fluctuations is the key difference between the long-range and/or energy-biased ice models and the short-range and energy-unbiased ice models. This difference is related to the appearance of long-range order at low temperature.

Finite-size scaling

Any system simulated by a computer has by definition a finite size L. However, when investigating a statistical physics problem, in particular a phase transition, one is interested in the thermodynamic limit (L →∞). Fortunately, close to phase transitions there exist systematic ways to extract the critical behavior in the thermodynamic limit. These methods, known as finite-size scaling, are based on the scaling hypothesis for critical phenomena. There are many references on both critical phenomena [START_REF] Goldenfeld | Lectures on Phase Transitions and the Renormalization Group[END_REF][START_REF] Stanley | Introduction to Phase transitions and Critical Phenomena[END_REF] and finite size scaling [START_REF] Cardy | Finite-Size Scaling[END_REF], so we will only sketch the principles of the method here.

Order parameter

A phase transition is characterized by an order parameter m, whose value can fluctuate between -1 and 1. Above the transition temperature, the average value of m is zero. Below the transition temperature, the fluctuations are suppressed and the order parameter takes a finite value ±m 0 , signaling the ordering of the system and the corresponding symmetry breaking.

In the thermodynamic limit, this spontaneous symmetry breaking occurs because the symmetry-broken phases differs by an extensive number of local changes. Therefore at low temperature fluctuations between these states are suppressed.

In a finite system however, the system has a finite probability to go from one ordered state to another. Indeed, non-local moves are precisely designed to allow the simultaneous update of a large portion of the system. As a consequence, the average value of the order parameter remains zero even in the ordered phase, as the system fluctuates between the different low temperature states.

Instead of calculating the average value m of the order parameter, it is then preferable to estimate the average of its magnitude |m|. This quantity goes from zero in the disordered phase to 1 in the ordered phase. In the thermodynamic limit, it converges towards the magnitude of the average of the order parameter |m| ----→ L→∞ |m| .

(3.15)

Scaling

When a system approaches a continuous phase transition, the fluctuations of its energy and of its order parameter diverge. In particular, fluctuations occur at all length scales. Indeed, the only relevant length scale in a critical system is the correlation length. At the transition, this length diverges, rendering all other scales irrelevant. One must then develop a scale-free description of the system. This property of scale invariance is at the heart of the study of critical phenomena. It can be encoded in one of the corner stones of the study of critical phenomena, the scaling hypothesis. It describes how the singular part of the free energy g s (associated with the observed divergences in the thermodynamical quantities) changes when the length scale is changed. For a thermal transition at a critical temperature T c this hypothesis states that g s transforms as

g s (τ, h) = l -d g s (τl y τ , hl y h ) (3.16)
through a change of scale r → l × r, with τ = |T/T c -1| the reduced temperature, h the source field of the order parameter, and d the dimension of the system [START_REF] Widom | Equation of state in neighborhood of critical point[END_REF][START_REF] Widom | Surface tension and molecular correlations near critical point[END_REF][START_REF] Domb | On critical behaviour of ferromagnets[END_REF]. The exponents y τ and y h are not specified in this hypothesis. However, equation (3.16) implies that all thermodynamical quantities will scale in a similar way. A phase transition will be characterized by a power law behavior of several quantities descending from the scaling hypothesis of the free energy. The power laws are characterized by critical exponents that are universal in the sense that they do not depend on the specifics of the system, but only on a few salient properties (its symmetries, the range of the involved interactions, the dimension). A given set of exponent defines the universality class of a phase transition. In fact, the different critical exponents are not independent. Indeed, in most cases, two of them are enough to characterize the universality class.

Using the standard notation for the correlation length ξ, the order parameter m, the specific heat C and the susceptibility χ, their critical behavior reads

ξ(τ, h = 0) ∝ τ -ν C(τ, h = 0) ∝ τ -α χ(τ, h = 0) ∝ τ -γ (3.17) and m(τ, h = 0) ∝ τ β m(τ = 0, h) ∝ h 1/δ . (3.18)
These exponents are related through the scaling relations (for more details, see e.g. [START_REF] Fisher | Renormalization group theory: Its basis and formulation in statistical physics[END_REF]) :

νd = 2 -α 2β + γ = 2 -αβ (δ -1) = γ. (3.19)
The scaling hypothesis holds in the thermodynamic limit. Indeed, for a finite system, all quantities remain analytical, and the singularities are rounded. This rounding has been first pointed out in the two dimensional Ising model [START_REF] Fisher | Scaling theory for finite-size effects in critical region[END_REF][START_REF] Ferdinand | Bounded and inhomogeneous ising models .i. specific-heat anomaly of a finite lattice[END_REF] for which the correction was known exactly. The scaling hypothesis has been later applied to a generic finite size system [START_REF] Fisher | Scaling theory for finite-size effects in critical region[END_REF] and it is now commonly used to determine the universality class of a model from Monte Carlo calculations.

Finite size scaling

In a finite system, the correlation length is naturally bounded by the linear system size L. In general, any thermodynamical quantity is in this case a function of ξ/a and L/a, where a is a characteristic microscopic length of the system. The finite-size scaling hypothesis assumes that, close to the critical point, any (microscopic) length other than the correlation length drops out. The only dimensionless quantity one can then construct from ξ and L is their ratio. The expressions obtained for an infinite system are then extended to finite systems by including x = ξ/L as a parameter.

Thermodynamical quantities

Let A be a quantity that scales as τ -a in the thermodynamic limit. In a finite system, the dependence in ξ/L has to be included in a universal function φ A , so that

A = τ -a φ A (ξ/L) = ξ a/ν φ A (ξ/L) = L a/ν x -a/ν φ A (x) , (3.20)
where φ A (0) 0 (so that in the thermodynamic limit A = φ A (0) ∝ ξ a/ν ). One can then write A as

A = L a/ν φA (L 1/ν τ) (3.21) with φA (x) ∝ x -a φ A (x -ν ).
In particular, at T = T c , the correlation length grows linearly with the linear system size L. In practice, one plots AL -a/ν versus L 1/ν τ for different system sizes. If the exponents are correctly chosen, the curves collapse onto a single universal curve φA (x).

Critical temperatures

The value of the critical temperature T c enters the definition of τ. Unfortunately, this value is generally not known in advance. One could try to determine T c as the value that gives the "best" collapse of the curves of different sizes. Yet a joint determination of T c and of the critical exponents might have too many free parameters to be adjusted correctly. It is better to minimize the number of parameters determined simultaneously by a single finite-size scaling collapse.

In fact, one can get one more piece of information out of the previous reasoning. In a finite size system, the quantities diverging at the critical point show rounded peaks at a fictitious critical temperature T c (L) which converges to the true critical temperature T c only when L →∞. To estimate the scaling of T c (L), one can use the condition that for T = T c (L), ξ(Tc(L)) ∝ L and hence

|T c (L) -T c | ∝ L -1/ν . (3.22)
Fitting the obtained values of T c (L) to this expression allows then to extract T c and ν.

A more precise way to determine the critical temperature of a system is to compute the Binder cumulant [START_REF] Binder | Finite size scaling analysis of ising-model block distribution-functions[END_REF] of the order parameter. It is defined as

U 4 = 1 - m 4 3 m 2 2 . (3.23)
Assuming a gaussian distribution of the order parameter around m 0 with standard deviation σ, one gets U 4 = 0 in the ordered phase (m 0 = 0) and U 4 = 2/3 + O(σ 2 ) in the ordered phase (±m 0 0). U 4 is a universal function of ξ/L in the critical region. In particular, curves for the Binder cumulant for different system sizes cross as a function of temperature at a fixed point value U * (in the sense of renormalization groups). This crossing occurs at T c (L), defined as (3.22). The value of U * depends on the universality class of the transition and on the boundary conditions.

Corrections to scaling

For most systems, the finite-size scaling works very well. However, the data sometimes fall off the expected curves. Indeed, there are sometimes logarithmic corrections beyond this power law dependency of thermodynamic observables Instead of being proportional to τ -a , these quantity are proportional to τ -a | log τ| â. Most of the times this correction is barely visible and can be simply neglected. In some cases, in particular if a = 0 (like for the specific heat of the two dimensional Ising model), this correction has to be added in order to explain numerical or experimental observations. A lot of effort has been put recently into determining the proper derivation of these correction and the value of the logarithmic exponent of some models [START_REF] Kenna | Scaling relations for logarithmic corrections[END_REF][START_REF] Palchykov | Critical phenomena on scale-free networks: Logarithmic corrections and scaling functions[END_REF]. For a detailed review of the recent developments on the subject, see [START_REF] Kenna | Universal scaling relations for logarithmic-correction exponents[END_REF]. We now briefly describe what consequences this corrections have on the finite size scaling.

Logarithmic corrections

In the following, we will use standard notations for the critical exponents (α, β, ν, . . .). The corresponding logarithmic critical exponent will be noted with the same greek letter with a "hat" ( α, β, ν, . . .). With logarithmic corrections, the scaling of the thermodynamic quantities becomes

ξ(τ, h = 0) ∝ τ -ν | log τ| ν C(τ, h = 0) ∝ τ -α | log τ| α χ(τ, h = 0) ∝ τ -γ | log τ| γ (3.24)
and

m(τ, h = 0) ∝ τ β | log τ| β m(τ = 0, h) ∝ h 1/δ | log h| δ . (3.25)
An important modification to the standard scaling approach is how the correlation length at criticality vary with system size

ξ(τ = 0, h = 0) ∝ L(log L) q (3.26)
If q > 0 then the correlation length can effectively be larger than the size of the system. These logarithmiccorrection exponents [START_REF] Kenna | Scaling relations for logarithmic corrections[END_REF] are related to each other by scaling relations that are very similar to the ones for the standard exponents :

d(q -nu) = α -δ α,0 2 β -γ = α -δ α,0 β(δ -1) = δ δ -γ (3.27)
Let us now derive the finite-size expression implied by these logarithmic factors. Consider a quantity A with exponents a and â (i.e. A ∝ τ -a | log τ| â). Like in the standard case, we need first to express τ as a function of ξ. To the leading order, we get

τ ∝ ξ -1/ν | log ξ| ν/ν .
(3.28)

In the thermodynamic limit, to the leading order in ξ, one gets then

A ∝ ξ a/ν | log ξ| â-aν/ν . (3.29)
For a finite size system, we correct again A by a factor that only depends on ξ/L. The dependence of A with L and τ is then slightly modified

A = L a/ν (log L) â-aν/ν φ A L 1/ν τ| log τ| -ν/ν . (3.30)
The proper collapsing plot is then AL -a/ν (log L) -â+aν/ν versus L 1/ν τ| log τ| -ν/ν . Similarly to conventional scaling, the critical temperature can be extracted from the fit of T c (L) to

T c (L) -T c ∝ L -1/ν (log L) ν/ν . (3.31)
An example : the two dimensional four-state Potts model An example of model for which the logarithmic corrections are important is the q = 4 two dimensional Potts model [START_REF] Salas | Logarithmic corrections and finite-size scaling in the two-dimensional 4-state potts model[END_REF]. The Potts models are a family of generalizations of the Ising model (see [START_REF] Wu | The potts model[END_REF] for a review on these). It consists of a set of N variables σ i on a lattice that can take their values in {1, q}, where q is an integer that defines the model. The Hamiltonian of the system is then

H = -J i, j δ σ i ,σ j (3.32)
where δ is the Kroenecker function.

On a two-dimensional lattice, this model presents a phase transition at finite temperature for any q ≥ 2. For q = 2, this model is precisely the Ising model. The only non zero logarithmic critical exponent is α = 1. For q > 4, the transition is first order. The case q = 4 is then peculiar. For this model, logarithmic corrections are present in all quantities. In table 3.1, both sets of critical exponents (for standard and logarithmic scalings) are listed for the q = 2 (Ising) and q = 4 Potts models.

q ν, ν α, α β, β γ, γ q 2 (Ising) 1,0 0,1 1/8,0 7/4,0 0 4 2/3, 1/2 2/3,1 1/12
, -1/8 7/6,3 /4 0 Table 3.1: Critical exponents for the Ising (q = 2 Potts) and q = 4 two-dimensional Potts models [START_REF] Kenna | Universal scaling relations for logarithmic-correction exponents[END_REF]. Only the specific heat presents logarithmic corrections for the Ising model.

Equivalence of the vertex and Ising models

We have described square ice systems as vertex models. In this picture, the notion of closed loops is very intuitive, and the ice rules correspond explicitly to a divergenceless field. However, in this picture there are two inequivalent types of pseudospins (horizontal and vertical), even without considering the symmetries of the interactions. The structure factor has then several possible definitions (correlations between parallel pseudospins, between orthogonal pseudospins,. . . ). Besides, square ice can also be written as an Ising model on the checkerboard lattice. In this model, the spin structure factor is uniquely defined. As we are particularly interested in the properties of the structure factor we choose to study the properties of square ice as an Ising model.

In this section we describe how the properties of the vertex model translate into the language of the Ising model. From now on we will refer to the pseudospins of the vertex model as arrows in order to avoid confusion with Ising spins. Let us first describe the mapping which leads from the vertex model to the Ising model, as illustrated in figure 3.4. Starting from a 6-vertex configuration (Fig. 3.4(a)), one maps the sign of the projections of the arrows along, e.g., the y-axis onto Ising spins (pointing up for a positive projection and down otherwise -Fig. 3.4(b)). Flipping the Ising spins of every other row (Fig. 3.4(c)), gives zero (Ising-spin) magnetization on each vertex if the corresponding vertex configuration is a 6-vertex one obeying the 2-in/2-out ice rule (a similar mapping is obtained by flipping every other column). In the following, we will indicate as vertices (denoted by ) the squares with additional diagonal couplings, and as plaquettes (denoted by ) the squares without diagonal couplings. The links of the underlying square lattice correspond to J 1 couplings of the square ice model, whereas the J 2 couplings become diagonal couplings on half of the squares. In particular vortex vertices map onto Néel vertices for the Ising spins (with ferromagnetic links on the diagonals), while current vertices are mapped onto collinear vertices (with antiferromagnetic links on the diagonals).

Through this mapping, the vortex and the Manhattan states transform respectively into a Néel and a collinear states (see figure 3.5). In the Néel state, the diagonal links are frustrated (i.e. the spins it connects are ferromagnetically aligned) whereas in the Collinear state they are unfrustrated. In the "current" states of the vertex model, the horizontal and vertical lines of parallel arrows were ferromagnetically aligned. After the mapping onto an Ising model, these lines become diagonals with a staggered spin arrangement.

Structure factor

When considering the structure factor of the vertex model, one has to distinguish the correlations between arrows along parallel directions from those between arrows along orthogonal directions (see section 1.2.4). In the case of the Ising model on the checkerboard lattice, the structure factor is uniquely defined, and it involves the correlation between the components of the vertex-model arrows along the x + ŷ direction

C(r) = C xx (r) + C xy (r) + C yx (r) + C yy (r) (3.33)
vortex Néel Manhattan (π/2, π/2) collinear Figure 3.5: Ising model equivalents to the Vortex and the Manhattan states. The latter becomes a collinear state with a 8-site unit cell with a (π/2, π/2) ordering vector and the former a Néel state. All the current states map onto Ising states with antiferromagnetically ordered diagonals along the J 2 direction (dashed lines).

Using (1.11), the expected expression for these correlations in the 6-vertex model around the origin reads [START_REF] Youngblood | Correlations in ice-rule ferroelectrics[END_REF]

S q ≈ q→0 1 K (q x -q y ) 2 q 2 (3.34)
This predicts the presence of a pinch point located at the origin of the first Brillouin zone. The structure factor obtained with a Monte Carlo calculation of the 6-vertex ice model at T = 0 is shown in figure 3.6. There are other pinch points at the (q x , q y ) ∈ πZ 2 , that are equivalent to the one at the origin. The mapping from the vertex model to the Ising model consists of a π/4 rotation, followed by the flip of half of the lines of spins. The last step corresponds to multiplying all spins by a factor e iπy where (x, y) are the coordinates of the spins. The structure factor of the Ising model will then be obtained from that of the square ice model with the same rotation, followed by a translation of πŷ (see figure 3.6). The pinch point at the origin transforms into a pinch point at the (0,π) point. Its expression is given by

S I q = q=(0,π) 1 K q 2 x q 2
x + (q yπ) 2 .

(3.35)

The other equivalent pinch points are present at the (0, -π) and (±π, 0) points. In all the following discussion, we will use the Ising model description of the lattice and structure factor.

Loop updates and winding number

In the Ising representation, the ice states are connected to each other by flip of closed loops of alternating spins running on the links of the checkerboard lattice. These loops connect sites of the square lattice of vertices (see figure 3.7). Figure 3.7: Loop updates on the checkerboard lattice. The ice states are connected by the flip antiferromagnetic loops connecting sites of the dual lattice (green dots). An example of such a loop is represented in green. Flipping this loop does not change the magnetization of any vertical or horizontal line. For instance, it crosses the purple line one a down spin and an up spin. Flipping these spins will note change the magnetization of this line.

If the loop is not winding around a lattice direction, it does not change the magnetization of a given column or row. To explain this property, let us consider a loop crossing twice the y = y 0 row (see figure 3.7). A step of the loop gives a ± displacement in the y direction. These displacement must add up to zero between the two crossings. It crosses then the line y 0 an even number of times. As the spins of the loop are in a staggered arrangement, the two crossings must have an opposite orientation. Hence flipping them will not change the magnetization of the line. Only a non local winding-loop update can then change the magnetization of a line. Furthermore, because all vertices are neutrally charged, two adjacent rows (or columns) have opposite total magnetization in an ice configuration. The topological sector is then characterized by the total magnetizations w x and w y of a reference raw and a reference column.

At T = 0, one can use the fluctuations in the magnetization of a fixed line to determine whether or not topologically non-trivial updates occur in the system. If these fluctuations vanish, then the system is locked in a topological sector. Note that the winding numbers correspond to the (0,π) and (π, 0) magnetizations, and hence the fluctuations in the winding number are given by the spin structure factor at the same points. These points are precisely the location of the pinch points.

Chapter 4

Phase diagram of the dipolar square ice

We now go back to the dipolar square ice. The aim of this chapter is the quantitative reconstruction of its phase diagram and, in particular, our goal is to characterize the disordered phase at low (but finite) temperature, and the transition from this phase to the ordered phase(s). The model we study is the square array of dipoles made of two point-like charges, as described in section 2.2.2. In particular, we characterize how these results change when truncating the dipolar interactions.

Monte Carlo

Updates

In order to accelerate dynamics in our Monte Carlo calculations, we used a combination of various update moves. To simulate the thermal bath and ensure the proper density of defects, we used Metropolis single spin-flip updates. As for any ice system, we also used the long-and short-loop algorithms to prevent the system from freezing at low temperature. Indeed, a finite density of defects is responsible for the fluctuations of the system among the various configurations satisfying locally the ice rules (modulo the presence of dilute defects). These fluctuations comes from the fact that the defects flip lines of spins on their path as they diffuse through the system. Our algorithm is able to "teleport" them over finite lengths and hence accelerates their dynamics. In the ideal square ice, the lower energy excitations are creations of a pair of defects on neighboring vertices annihilating each other on another vertex and flipping a closed loop in the process.

As we expect two different regimes at low temperature, we used an additional move adapted to the ground state. Wolff cluster moves accelerate the relaxation towards the proper Néel ordering at small h. For large h on the other hand, there is a whole family of low energy states : the current states (see figure 2.3). These states are connected with each other through the flip of a straight line of L spins, changing the winding number. Such a "straight loop" can be obtained from the loop algorithm, but only if all the spins along the loop form a staggered pattern. Furthermore, any ice vertex encountered by the loop offers a binary choice for the exit spin, and therefore the probability to create a straight loop in the ordinary loop algorithm is of order 2 -L . The loop algorithm will then be very ineffective in this regime (see figure 3.3). In order to improve the convergence, we include a "single line" flip move. It tries to flip all the spins of a random line with a probability given by the Metropolis prescription. This update accelerates the convergence to the ordered state.

Monte Carlo Steps

Let us now turn to the thermodynamics of the model. A typical Monte Carlo step is composed of N = L 2 spin-flips, N/4 short-loop updates, 2N/L 5/3 = 2L 1/3 long-loop updates and L cluster updates or 2L line flip updates. The results presented here were typically obtained from series of 10000 bins of 200 Monte Carlo steps each. In order to obtain the correct low temperature ground state, and to prevent the system from freezing into a metastable state, we perform a linear annealing of the system, gradually lowering the temperature of the system from 4J 2 , and allowing an initial thermalization of 10000 MC steps at each temperature.

Phase diagram

We choose to take a trap half-width λ = 0.2 a, and we vary the height offset h from 0 to 1.2 a. This corresponds to varying the ratio between the two greatest couplings J 2 /J 1 from 0.48 to 1.94. The phase diagram obtained from Monte Carlo simulations is represented in figure 4.1. From now on, all energies are expressed in units of J 2 , which does not change with h. At high temperature, the system is paramagnetic. At low temperature, the system orders, into a Néel state for small h, and into a (π/2, π/2) collinear state for large h. The ordering temperature is much lower than the spin-flip gap above the ground state. Therefore there is a clear temperature range between the paramagnetic phase in which the defects proliferate (signaled by the single-spin flip gap, represented by the blue line in figure 4.1), and the ordered phase, in which the density of defects vanishes exponentially with the inverse temperature. The phase boundary between the ordered states are very close to the normal mode bandwidth calculated in section 2.2.4.

The properties of this phase are very similar to those of the Coulomb phase, and correspond to a collective paramagnetic phase, as indicated in figure 4.1. In particular, the pinch points are clearly visible in the structure factor (see section 4.2.1).

The ordering transitions are second order. The details about these transitions, in particular their universality class, are given in the next section. The transition from the Néel to the collinear phase as h is increased at fixed low temperature should be first order, as these phases have different symmetries. Indeed, we observe discontinuities in the order parameter in this case.

Coulomb phase

Considerable theoretical work has been devoted to describing the regime where defects are absent, or sufficiently rare to be unable to mask the underlying constrained manifold of ice states and the Coulomb phase features. It has been widely studied from different perspectives (in presence of long range interactions [START_REF] Isakov | Why spin ice obeys the ice rules[END_REF][START_REF] Siddharthan | Ising pyrochlore magnets: Low-temperature properties, "ice rules," and beyond[END_REF][START_REF] Castelnovo | Magnetic monopoles in spin ice[END_REF], finite temperature [START_REF] Sen | Coulomb phase diagnostics as a function of temperature, interaction range, and disorder[END_REF] or quantum fluctuations [START_REF] Balents | Fractionalization in an easy-axis kagome antiferromagnet[END_REF][START_REF] Tewari | Emergence of artificial photons in an optical lattice[END_REF][START_REF] Hermele | Pyrochlore photons: The u(1) spin liquid in a s = 1 2 three-dimensional frustrated magnet[END_REF]).

The presence of pinch points in the magnetic structure has been used several times as a prognostic tool for its identification (see section 1.2.4 for more details). We find here that this feature is indeed present in the dipolar square ice, the structure factor becoming that of an ordered phase through the transition line, with the apparition of Bragg peaks (see below).

These pinch points are not suppressed if the density of defects is sufficiently small. Fortunately, albeit finite, the density of defects is exponentially small for temperatures below the spin-flip gap. The defects are then unable to destroy the correlations in the system. In fact, they are responsible for the fluctuations between the ice states by hopping from site to site along (closed) loops.

Another important property of the Coulomb phase is the deconfinement of the excitations. It is deeply connected to the fluctuations of the winding number of the system. Even though the winding number is a T = 0 property, we can give it a definition at finite temperature. As explained in section 3.6.3, the height of the pinch point S q (π, 0) is a measure of the winding-number fluctuations. This is well defined at any temperature, and serves as a criterion to identify the confinement of the deconfinement of the monopoles in the system. Indeed, if the winding number fluctuates thanks to the flip of winding loops, then the defects are deconfined, as they can propagate through the whole system following such loops. If the winding number is fixed however, these loops are suppressed, and the excitations are confined.

Structure factor

The magnetic structure factor is defined as

S z q = 1 L 2 i e iq.r i σ i 2 (4.1)
with σ i = ±1 the the spin at position r i of the lattice. Experimentally, spin ice materials have been identified by the existence of pinch points in their magnetic structure factor [START_REF] Bramwell | Spin ice state in frustrated magnetic pyrochlore materials[END_REF][START_REF] Fennell | Magnetic coulomb phase in the spin ice ho2ti2o7[END_REF]. Because they are by definition vanishingly thin, the pinch points demand high accuracy measurements to be unarguably identified. Similarly, the resolution on the structure factor calculated on a finite system is limited by its size. Indeed, the pinch point has a width of order ξ -1 ∼ n -1/2 , with ξ the correlation length and n the density of defects. Hence a finite-size system is effectively in an ideal Coulomb phase if ξ>L, so that the pinch point remains resolution limited. For the optimal value of h, the structure factor is very similar to that of the ideal model in a wide range of temperatures (figure 4.2). The pinch points are clearly the dominant feature of the structure factor. The system is then precisely in the regime where the temperature is large enough to overcome the energy barriers betweens ice states, but small enough to prevent the proliferation of defects and the destruction of the Coulomb phase. At lower temperatures, Néel Bragg peaks appear. The evolution of the structure factor

h > h LR c
If the height offset is too large, the range of temperatures over which the structure factor develops pinchpoint-like features but no Bragg peaks is reduced, but remains finite. As the temperature is lowered, the structure factor develops scattering ridges, corresponding to the different collinear states that are, in this case, less energetic than the other ice states (figure 4.4). As a comparison, the structure factor obtained when restricted to the collinear states only is represented in figure 4.5. The same ridges are found on the lines q x ± q y ∈ π(2Z + 1). When the temperature is further reduced, Bragg peaks eventually appear at (±π/2, ±π/2), corresponding to the expected q = (π/2, π/2) collinear ground state.

h < h LR c
For small height offset, the only ground states are the Néel states. This is the relevant case for experimental two dimensional spin ices. Starting from a Néel state, the lowest excitation is the flip of a plaquette. It turns four vortex vertices into current vertices. It costs an energy ∆ p = 8(J 1 -J 2 ). The spin-flip gap above aNéel state is ∆=2(2J 1 -J 2 ). If J 2 < 2J 1 /3 (which occurs at h 0.51 a), the lowest excitation becomes a single-spin flip. The fluctuations are then governed by spin flips and not loop flips, and the ice rules are no longer relevant. Hence monopoles proliferate in the magnetically disordered phase, and there is only a weak trace of pinch point in the structure factor (figure 4.6).

Pinch point width

We define the width δq of the pinch point at (π, 0) as the half width at half maximum of the structure factor along the q y = 0 line. As the temperature is increased, the structure factor becomes flatter. In particular, its minimum, located at q = 0 increases. Therefore we define the height of the pinch point as S (π, 0) -S (0, 0). In the ideal square ice at T = 0, this width tends to zero in the thermodynamic limit as the structure factor becomes a Dirac delta in this direction, as can be seen in figure 4.7. Here, however, the width is finite, because of the long range interactions and the finite temperature (figure 4.7). Nevertheless, its dependence in the temperature is very similar to that of the ideal square ice for a wide range of temperature above the ordering transition. This corresponds to the thermal Coulomb phase, in which very few defects are present.

The sudden increase of the pinch point width at low temperature signals the onset of the Bragg peaks and the transition to the ordered ground state. When the width reaches π, the pinch point disappears as its extent reaches the entire Brillouin zone.

Despite the confinement of the defects induced by the long-range interactions, the correlations of the system have a similar form as in the ideal case. This is possible thanks to the finite temperature, that effectively smooths out the small energy difference between the ice states induced by the long-range interactions. This is consistent with other studies of long-range interacting ice models at finite temperature [START_REF] Sen | Coulomb phase diagnostics as a function of temperature, interaction range, and disorder[END_REF].

Density of defects

The presence of defects at finite temperature is unavoidable. Their density can even be rather high, as can be seen in figure 4.8. However, the presence of defects does not necessarily indicate that the Coulomb phase is destroyed. Indeed, the higher density of defects is observed for the h ≥ h LR c cases, for which the pinch points are more clearly visible in the structure factor than for the h ≤ h LR c cases.

The presence of defects changes the algebraic decay of the correlations of the ideal Coulomb phase to an exponential decay. This will primarily affect the long distance behavior, and hence the long wave length structure factor. This will then have a similar effect as a finite resolution [START_REF] Youngblood | Correlations in ice-rule ferroelectrics[END_REF]. As the temperature is increased, the correlation length becomes smaller, and hence the pinch point becomes wider. This increase of the pinch point width with temperature is observed in our results, so long as it is in the thermal Coulomb phase. This density is independant of the system size. This indicates that the tail of the long range interaction has practically no effect on the creation of a defect. The energy cost of a spin flip then only comes from closeby spins.

T /J 2 , T /J T /J 2 , T /J T /J 2 , T /J h =0.2 a h =1.2 a h =0.

Transition to the Néel state

For small height offset h, the low temperature phase is the Néel state, with a doubly degenerate ground state. We expect then a second-order phase transition belonging to the Ising universality class.

The order parameter m of the Néel phase is the square root of the (π, π) structure factor

m = 1 N i (-1) x i +y i σ i . (4.
2)

The specific heat and the susceptibility were measured as the fluctuations of the energy and of the order parameter

C = β 2 E 2 -E 2 χ/N = β m 2 -|m| 2 . ( 4.3)
The boundary between this phase and the paramagnetic phase in the phase diagram 4.1 were determined with finite-size scaling analysis of the thermodynamic quantities using the Ising critical exponents. This is indeed the correct universality class, as can be seen in the collapse plots in figure 4.9.

Transition to the (π/2, π/2) collinear state

For large h, the dominant interaction is the interaction J 2 along the diagonals of the vertices. However, this coupling is not sufficient to select a finite number of ground states. In order to predict the correct ground state, one has to include the coupling J 3 on the diagonals of the plaquettes. The predicted ground state is then the (π/2, π/2) collinear state. It has a 8-site unit cell and is four-fold degenerate (see figure 2.4). The two universality classes are compatible with a four-fold degenerate ground state : the 4-state clock model and the 4-state Potts model.

Universality classes

Clock models The clock models are a family of models that interpolate between the Ising and the XY models. The variables are rigid rotors S i = (cos ϑ i , sin ϑ i ) sitting on the sites of a lattice. These rotors are only allowed to point in q discrete directions of the plane ϑ i = 2πi/q, i = 1,...,q. The parameter q defines The clock models. The spins rotors point in any of the q angles defined by the model. For q = 2, one recovers the Ising model, whereas the XY model is obtained in the q →∞limit. (right) Spectrum of the single-bond energy for q = 4.

(T -T c )/L (T -T c )/L (T -T c )/L (T -T c )/L
T c = 3.25J 2 . q =2 • • q =4 • • • • q =6 • • • • • • q →∞ • • • ε = 0 ε = -J ε = J q = 4
which clock model is considered. For q = 2, these rotors are Ising spins, whereas for q →∞they become continuous XY spins. The Hamiltonian of this model is

H (q) = -J i, j S i .S j = -J i, j cos(ϑ i -ϑ j ) . (4.4)
For q = 4, the interaction energy between two sites i and j can take three values : ±J if ϑ i = ±ϑ j , and 0 if S i and S j are orthogonal. Above a configuration of minimal energy -J, there are two excited states (modulo a global rotation of all the rotors). For the q = 4, the clock model belongs to the Ising universality class [START_REF] Taroni | Universal window for two-dimensional critical exponents[END_REF]. Its critical exponents are given in table 3.1.

Potts models

The Potts models are another generalization of the Ising model. The variables are here integers σ that can take any values 1, 2 ...,q. The Hamiltonian reads

H (q) = -J i, j δ σ i ,σ j (4.5)
where δ is the Kronecker function. The Potts model. The variables can take any value between 1 and q. These can be viewed as q different colorings. For q = 4, there is a single minimal energy configuration with three excitations for each pair of interacting sites.

q = 2 q = 3 q = 4 ••• ε = 0 ε = -J q = 4
As in the clock model, for q = 2 one recovers the Ising model. In two dimensions, for q > 4, the ordering transition is first order. For q = 4, as in the clock model, there is a single minimal energy configuration for a pair of sites (modulo a global relabeling of the 4 single site configurations). However, the lowest excitation is here three-fold degenerate (for a fixed value of one of the two sites) with no higher energy configuration. This is an important difference from the clock model. The universality class of this model has important logarithmic correction to scaling (see section 3.5.4)

In view of the 4-fold degeneracy of the ground state of the dipolar ice for h > h LR c , the universality classes of the 4-state clock model and the 4-state Potts model represent two strong candidates for the scaling analysis of the Monte Carlo data. Yet we should keep in mind that other relevant models with a 4-fold degenerate ground state, such as the José-Kadanoff-Kirkpatrick-Nelson model [START_REF] José | Renormalization, vortices, and symmetrybreaking perturbations in the two-dimensional planar model[END_REF] with q = 4 as well as the J 1 -J 2 Ising model on the square lattice [START_REF] Kalz | Location of the potts-critical end point in the frustrated ising model on the square lattice[END_REF], feature continuously varying exponents (with fixed β/ν = 1/8).. Here we take the least ambitious goal of testing the consistency of our data with two universality classes only (q = 4 Potts and Ising). While continuously varying exponents would in principle suggest to perform a fully unbiased scaling analysis of our data, the long-range interactions that we deal with limit in practice the range of sizes that we can study, preventing us from obtaining unambiguous estimates of the critical exponents from an unbiased analysis.

Discriminating between the two universality classes is a subtle exercise, as they both share the same β/ν = 1/8 and γ/ν = 7/4 ratios. Furthermore, the logarithmic corrections of the q = 4 Potts model render the finite size scaling analysis difficult for accessible system sizes. However, the results are more consistent with the Potts universality class.

The ground state state has a standard (π/2, π/2) structure, so that the most natural order parameter is

m = 1 L 2
x,y e iπ(x+y)/2 σ(x, y)

. ( 4

.6)

As for the Néel transition the specific heat and the susceptibility were measured as the fluctuations of the energy and of the order parameter

Critical temperature

We extract the transition temperature from the extrapolation of the position of the peaks of the specific heat and of the susceptibility, using both the Potts and clock ansätze. For the clock model, we expect the "critical" temperature to scale as

T c (L) = T c (∞) + aL -1/ν = T c (∞) + aL -1 (4.7)
whereas the Potts model predicts The two resulting critical temperatures are very close (figure 4.12). For the Potts scaling, we obtain T (P) c /∆ 0 = 0.901 ± 0.002 whereas for the Ising the estimated transition temperature is T (I) c /∆ 0 = 0.899 ± 0.004. These estimates are very close, hence further analysis is needed to discriminate between the two universality classes.

T c (L) = T c (∞) + aL -1/ν (log L) ν = T c (∞) + aL -3/2 (log L)
0 2 2¤ 2 2¥ 2 2¦ 2 2 § 2 2£ L -3/2 (log L) 3/4 L -1

Binder cumulant

The transition temperature obtained from the crossing of the Binder cumulant curves (described in section 3.5.3) is T (I) c /∆ 0 = 0.9107 ± 0, 01 (see figure 4.13). It is consistent with both of the previous temperatures, and hence is not sufficient to discriminate between the universality class. The difference is clearer when one tries to collapse the Binder cumulant curves onto each other using both sets of exponents. The Potts universality class gives a better collapse (figure 4.13).

Finite size scaling for the thermodynamic observables

Using the critical temperatures obtained with the corresponding scaling, we try to collapse the heat capacity, the order parameter and the susceptibility onto universal functions using the scaling exponents of the Potts and clock universality classes. The correct set of exponent should give a better collapse. We concentrate on the h = 1.2 a geometry and consider systems of linear size L ranging from 24 to 48. The results are shown in figure 4.14. The 4-state Potts ansatz leads to a much better collapse of the data.

In particular, the height of the heat capacity peak is correctly rescaled using the 4-state Potts critical exponents. As β/ν and γ/ν are the same in both Ising an 4-state Potts universality classes, the difference is less obvious when looking at the order parameter or the susceptibility.

C/log L χ L -7/4 (log L) 1/8 CL -1 (log L) 3/2 χ L -7/4 mL 1/8 mL 1/8 (T -T c ) L 3/2 (log L) -3/4 (T -T c ) L 0 0.5 1

J 1 -J 2 -J 3 model

As mentioned in section 2.2.4, the (π/2, π/2) collinear state is selected by the third neighbor interaction J 3 . One might then expect that a J 1 -J 2 -J 3 model would capture the universal feature of the transition, given that is has the same symmetries and the same ground state as the system with long-range interactions. The Hamiltonian in this case writes

H = J 1 i, j σ i σ j + J 2 i, j σ i σ j + J 3 i, j σ i σ j .
(4.9)

For the chosen set of geometrical parameters (λ = 0.2 a and h = 1.2 a), the couplings are J 2 ≈ 1.94 J 1 and

J 3 = -J 2 /2 ≈ 0.97 J 1 . J 1 J 2 J 3
Figure 4.15: The interactions are truncated at J 3 . The couplings are here J 1 on the links of the square lattice, J 2 on the diagonal of the odd plaquettes and J 3 on the diagonal of the odd plaquettes. This is the minimal set of couplings to promote the Manhattan state as a ground state.

We have seen in section 2.2.4 that the band width of the normal modes is larger in the J 1 -J 2 -J 3 model than if the full range of interactions is included. We expect then the transition temperature to be higher in the former case. We checked this prediction, and studied the transition in the J 1 -J 2 -J 3 model with Monte Carlo calculations. The update algorithm is the same as that used for the model with long-range interactions, as described in section 3.3. A finite-size scaling analysis is conducted on systems with sizes ranging from L = 24 to L = 128.

Critical temperature

We proceed as in the previous case, and extract the critical temperature from the extrapolation of the position of the peaks of the specific heat and the susceptibility (figure 4.16). Both Potts and Ising scaling are consistent with the data within error bars. For the Potts scaling, we obtain T (P) c /∆ 0 = 1.5487 ± 0.0005 whereas for the Ising the estimated transition temperature is T (I) c /∆ 0 = 1.543 ± 0.001. As anticipated, the transition temperature is slightly higher than in the dipolar case because the screening due to the tails of the interactions is absent here. This is indeed due to the fact that the band width of the normal modes is larger in this case, as discussed in section 2.2.4.

Binder cumulant

The Binder cumulants of the Manhattan order parameter have here a clear crossing around T c = 1.545 ± 0.002. This value is a bit closer to the estimates from the clock than for the Potts scaling. Rescaling the curves according to each universality class gives a qualitatively good collapse in both case (figure 4.17). However, the collapse is apparently better for the Ising exponents close to T c . This suggests that the universality class of the transition might change from 4-state Potts to Ising when the range of interactions is truncated.

T c (L) L -3/2 (log L) 3/4
L -1

U 4 U 4 U 4 T (T -T c )L 3/2 (log L) -3/4 (T -T c )L 0.

Finite-size scaling for the thermodynamic observables

We perform the same analysis of section 4.2.4 on the J 1 -J 2 -J 3 model. Much larger system sizes are numerically accessible in this case. The collapse of the thermodynamic quantities are represented in figure 4.18. The collapse of the heat capacity improves much faster with the system size if one uses the Ising critical exponents. Furthermore, the scaling of the susceptibility is significantly better in the Ising

From reduced to long range interactions

The above results seems to indicate that the universality class of the ordering transition to the (π/2, π/2) collinear state changes from Ising to 4-state Potts as the range of the interactions is increased. The main difference between the clock and the Potts models lies in the degeneracy of the first excited levels. In the 4-state clock model there are two different levels of excited pairs of interacting sites, whereas all the excited pairs of sites in the Potts model have the same energy. This difference can be related to the difference in the bandwidth of the dipolar and the J 1 -J 2 -J 3 models. A future work would be to examine the transition to the (π/2, π/2) collinear state for 1/r α interactions with α = 3, 4, 5,....There might be a change in the universality class of the transition for α larger than a particular value α c .

Conclusion

We find that tuning the geometry of the system allows to increase the temperature range for which dipolar square ice exhibits similar features in the structure factor as the ones of the nearest-neighbor square ice. In particular, we find that the pinch points survive the presence of the long-range interactions. This is possible because the ordering transition occurs at a critical temperature T c smaller than the spin-flip gap ∆ 0 that sets the typical temperature above which the ice rules cease to dominate the physics and the Coulomb phase is destroyed.

Furthermore, it turns out that changing the geometry of the system (in particular the height offset h) is an efficient way to select both the ground state and the universality class of the ordering transition.

However, in the geometry that is relevant in the nanoarrays experiments (h a), the bias towards the Néel state is so large that the Coulomb phase is never approached at any temperature. A way to circumvent the problem could be to reduce the range of the interactions. It could be possible if one used trapped colloids [START_REF] Libál | Realizing colloidal artificial ice on arrays of optical traps[END_REF] where the interactions would be screened by the environment, or trapped ions [START_REF] Kim | Quantum simulation of frustrated ising spins with trapped ions[END_REF] for which the interactions can be more easily tuned. In the latter case, quantum effects would be important. The next step is then to investigate the fate of the different ice features once the quantum effects are included in the analysis.

Part II

Exotic phases in quantum spin ice

In this part, we investigate the properties of the square ice model in presence of quantum fluctuations. In particular, we are interested in the low temperature phases of the transverse field Ising model on the checkerboard lattice.

Inspired by previous works on the model [START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF], we look for signs of a Coulomb phase as well as both thermal and quantum order-by-disorder mechanisms that could destroy this phase.

This part is organized as follows. First we analyze the model in its S →∞limit, and perform a linear spin-wave expansion. Then we develop effective lattice gauge theories in order to predict the different phases of the phase diagram. And finally we turn to numerical techniques. We describe a new efficient algorithm to simulate strongly constrained quantum systems and use them to determine the phase diagram of the transverse field Ising model on the checkerboard lattice.

Quantum square ice : the transverse field Ising model on the checkerboard lattice

Quantum fluctuations via a transverse field

There are various ways to introduce quantum fluctuations in an Ising model. In the case of three dimensional spin ice, couplings between transverse components of the spins have been considered in the literature [START_REF] Savary | Coulombic quantum liquids in spin-1/2 pyrochlores[END_REF][START_REF] Savary | Spin liquid regimes at nonzero temperature in quantum spin ice[END_REF][START_REF] Shannon | Quantum ice: A quantum monte carlo study[END_REF][START_REF] Benton | Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice[END_REF] and it has been proven that they can lead to different phases when the quantum effects dominates the physics.

Another possibility is to add a transverse magnetic field coupled to the spins. This also leads to an offdiagonal term in the Hamiltonian. If the original Ising degrees of freedoms were encoded in the position of an ion in a double-well trap, this transverse magnetic field corresponds to the tunneling amplitude between the two wells of a trap. The Hamiltonian for an ion in a double trap is then

H t = -t(b † l b r + b † r b l ) (4.
σz = 2b † r b r -1 = 1 -2b † l b l σx = b † l b r + b † r b l σy = i(b † l b r -b † r b l ) , (4.11)
and the tunneling Hamiltonian turns into a transverse field term

H t = -t σx . (4.
12)

The same Hamiltonian can be obtained by encoding the spin in the internal state of a trapped ion. In this case, the tunneling term is obtained by Raman coupling between the internal states [START_REF] Schneider | Experimental quantum simulations of many-body physics with trapped ions[END_REF].

From the square ice system to the transverse field Ising model on the checkerboard lattice Adding the transverse field term to the Hamiltonian of the square ice model, we get

H = 1 2 i, j J ij σz i σz j -Γ i σx i . (4.13)
This is the Hamiltonian of a transverse field Ising model (TFIM) with a transverse field Γ. When J ij is restricted to the links of a checkerboard lattice (see figure 4. [START_REF] Binder | Finite size scaling analysis of ising-model block distribution-functions[END_REF]), we will call this model the quantum square ice, whose properties are the topic of the second part of this thesis. The spins lie on the sites of the square lattice (black dots), and their z components are coupled with couplings J 1 (solid lines) and J 2 (dashed lines).

Chapter 5

Quantum square ice : large S limit and spin-wave analysis

In this chapter, we analyze the large S limit of quantum square ice. We test for possible order-by-disorder phenomena in the S →∞limit, and subsequently for quantum order-by-disorder in the semi-classical context of linear spin-waves.

S →∞limit

We first analyze the transverse field Ising model on the checkerboard lattice in its classical limit. The Hamiltonian is the same as in the quantum case (4.13), but instead of N Ising spins we allow for fluctuations in the transverse components of the spins. To this end we consider N = L × L classical continuous spins of length S that can point in any direction of a sphere. The Hamiltonian is

H = - i, j J ij S z i S z j -Γ i S x i S i = (S x i , S y i , S z i ), |S i | 2 = S 2 .
(5.1)

Ground states

At T = 0(i.e. in the ground state), the y component of the spins vanishes. Indeed, they only bring an entropic contribution to the free energy. If Γ=0, then the ground state configuration is S z i = σ i S , where σ i ± 1 is an Ising variable. A transverse field will perturb this configuration by canting each spin in the fielddirection by an angle ϑ i , so that S z i = S σ i cos ϑ i and S x i = S sin ϑ i , with ϑ i ∈ [0, π/2]. Furthermore, in any ice-rule state, each spin sees the same mean field along z, so that the canting angle is the same ϑ i = ϑ for all spins. This reduces the z component but does not change its sign, and adds a uniform x component to all spins (see figure 5.1). We note E 0 ({σ

i }) = 1 2 S 2 i, j
J ij σ i σ j < 0, the energy of an Ising configuration {σ i } in zero field. The energy of the corresponding canted configuration is then

E Γ ({σ i }) = E 0 ({σ i }) cos 2 ϑ -NΓS sin ϑ. (5.2)
Minimizing the energy with respect to ϑ, we find the field dependance of the canting angle :

sin ϑ = min NSΓ -2E 0 ({σ i }) , 1 ≡ min Γ Γ c , 1 Γ c ≡- 2E 0 NS (5.3)
The system becomes fully polarized in the x direction if the field is larger than a critical value that depends on the configuration. Its energy is then just E Γ ({σ i }) = -NSΓ. If the field is less than Γ c , the energy of the system writes

E Γ ({σ i }) = E 0 ({σ i })        1 + Γ Γ c 2        (5.4)
Hence, considering two Ising configurations {σ i } and

σ i , if E 0 ({σ i }) ≤ E 0 σ i , then E Γ ({σ i }) ≤ E Γ σ i for all Γ.
The energy hierarchy of the Ising configurations in zero field case is preserved by the canting.

From now on, we will restrict the discussion to in-vertex couplings. We note J 1 the coupling between adjacent sites of the square lattice, and J 2 the coupling on the interacting diagonals of the checkerboard lattice (see figure 4.19). They both are negative so that the interactions are antiferromagnetic. The Hamiltonian in this case reads

H = J 1 S 2 cos 2 ϑ i, j> σ i σ j + J 2 S 2 cos 2 ϑ i, j> σ i σ j -ΓS sin ϑ i σ i (5.5)
As discussed in section 2.2.3 (in the case of the vertex model), the ground state depends on the ratio J 2 /J 1 .

J 2 < J 1
If the diagonal coupling is smaller, the zero field ground state is the (2-fold degenerate) Néel state. In these states, all diagonal links are ferromagnetically arranged, whereas the horizontal and vertical links are in an antiferromagnetic configuration. All spins have the same energy and the total energy is E 0 (Néel) = N(J 2 -2J 1 )S 2 (5.6)

The critical field that separates the Néel phase from the fully polarized phase is then

Γ c (Néel) = 2S (2J 1 -J 2 ) (5.7)
Consequently, for 0 < J 2 < J 1 , the classical ground state of the system is a canted Néel state up to Γ=2S (2J 1 -J 2 ) where it becomes a fully polarized states (see figure 5.2).

J 2 > J 1
If the diagonal coupling is the largest, the ground states are the "collinear" states. In these states, all the interacting diagonals are antiferromagnetically aligned. There are 2 2L of these states as they are obtained from one another by the flip of diagonals. All the diagonal links are then in an antiferromagnetic configuration, whereas only half of the J 1 links are in the same configuration. The energy of the system then writes E 0 (coll.) = -NJ 2 S 2 (5.8)

The critical field that separates the collinear states from the fully polarized phase is then Γ c (coll.) = 2SJ 2 (5.9)

For J 1 < J 2 , the classical ground state of the system is the canted collinear state up to Γ=Γ c where it becomes a fully polarized states (see figure 5.2). The ground state degeneracy remains sub extensive (i.e. proportional to the system linear size L) throughout this phase as the energy hierarchy of the ice-rule states remains the same.

J 2 = J 1 ≡ J
The expression for the energy of the ground states and the critical field for both phases connect smoothly for J 2 = J 1 (E 0 = -NJS 2 and Γ c = 2SJ). At this point, all the ice-rule states have the same energy. The only difference from the zero field phase is that the spins are canted towards the field by a uniform angle.

Finite temperature

The transverse field preserves the ice-rule states degeneracy for J 1 = J 2 . At finite temperature however, a particular ordered state might be favored by the entropic effects, leading to an order-by-disorder effect. In order to determine whether this is the case, we performed Monte Carlo calculations of the model in the classical limit. The energy of a given configuration {(ϑ i ,ϕ i )} i=1,...,N is then given by (5.1). It reads

E = JS 2 i, j> i, j> cos ϑ i cos ϑ j -ΓS i sin ϑ i cos ϕ i (5.10)
At this point, we can take S = 1 without any loss of generality.

MC approach

We focus our study to the degenerate point J 2 = J 1 . We used the same type of algorithm as in the Ising case, with a few minor modifications.

Spin flip update

The Ising spin flip update is replaced by a continuous spin update. The orientation of a spin i is changed from (ϑ i ,ϕ i ) to (ϑ i ,ϕ i ) with a probability given by the Metropolis heat bath prescription.

Loops update

The loop algorithm for Ising spin ice [START_REF] Barkema | Monte carlo simulation of ice models[END_REF] is generalized to the case of continuous spins in the following manner: a loop is built in the same way as for Ising spins, using the sign of the z component as effective Ising spin variable; Only the z components of the spins are flipped by this update, but because its magnitude is not necessarily the same on each site, the update is not microcanonical, and the loop move has to be reweighted with the energy change E 0 σ i -E 0 ({σ i }) in order to satisfy the detailed balance condition.

In zero field, the energy cost of the flip of a spin in one of the ground states is ∆=4S 2 J. At finite Γ, this spin flip gap becomes

∆ Γ = 4JS 2 cos 2 ϑ =∆        1 - ΓS ∆ 2        (5.11)
Similarly to the square ice case, we expect the density of defects to be negligible for T ∆ Γ , namely

T/∆ 1 -(Γ/∆) 2
(5.12)

In this case, it must be possible to observe the ground state of the system in a Monte Carlo calculation. We place ourselves at T = 5 × 10 -3 J = 1.25 × 10 -3 ∆ which guarantees to observe the physics of the ground state up to Γ ∆=4 J.

At this temperature, we study the evolution of the z component structure factor and the transverse polarization, and look for possible signs of order.

Structure factor

In the Coulomb phase, the main features of the structure factor are the pinch points at (±π, 0) and (0, ±π), and peaks at (±π, ±π) that diverge logarithmically with the system size. All these features can be captured following the path in figure 5 The evolution of the structure factor S q with Γ is represented in figure 5.3. As the transverse field increases, the overall intensity decreases. The discontinuity at the pinch point is gradually smoothed and the logarithmic peak at the Néel point is rounded up similarly. However, for a large the transverse-field range, this evolution can be attributed to the reduction of the magnitude of the z component at T = 0. Indeed, if one rescales the value of S q by a factor cos 2 ϑ, with ϑ the canting angle defined at (5.3) the curves at different fields collapse onto each other, up to a field Γ 1.75 J. At this point, the spin-flip gap becomes sufficiently small to allow the thermal fluctuations to break the ice rules.

There is no sign of long-range order at this temperature for any value of the transverse field. In particular, the height of the peaks at the Néel points vanishes in the thermodynamic limit (figure 5.4).This can be inferred from the finite-size scaling of the order parameter m s estimated as

m 2 s = (1/L 4 ) ij (-1) i+ j S z i S z j ( 5

Linear spin wave analysis

The previous results can rule out the presence of thermally induced order-by-disorder. Here we analyze the possibility of ordering via quantum fluctuations at the harmonic level, as described by linear spin-wave theory. [START_REF] Coletta | Phase diagram of the fully frustrated transversefield ising model on the honeycomb lattice[END_REF][START_REF] Misguich | Quantum dimer model on the triangular lattice: Semiclassical and variational approaches to vison dispersion and condensation[END_REF]. We will treat separately the spectrum of excitations above the various classical reference states of the system : Néel, collinear, and fully polarized.

Spin-boson transformation

Let us consider a generic classical ground state with long-range magnetic order, and with a magnetic unit cell containing n spins. We denote S l,p the p-th spin (p = 1 ...n) of the l-th cell. As seen in section 5.1.1, in the classical limit an applied transverse field rotates the spins around the y-axis by an angle ϑ.W e introduce a local rotation of the spin configuration, Sl,p = σ p R y (σ p ϑ)S l,p (as illustrated in figure 5.5), where σ p = 1(-1) if the spin in zero field has positive (negative) projection along the z axis, and R y (±ϑ) is the rotation matrix of an angle ±ϑ around the y axis. The ground state is a simple ferromagnetic state for the Sl,p spins, namely S z l,p = S everywhere. We then consider small quantum fluctuations around this classical reference state, by transforming the quantum spins to bosons via a linearised Holstein-Primakoff transformation [START_REF] Holstein | Field dependence of the intrinsic domain magnetization of a ferromagnet[END_REF], valid in the limit of a small number of bosons n l,p 2S :

S z S x S ⇒ S z S x S z S x σ = 1 S z S x S ⇒ S z S x S z S x σ = -1
S z l,p = S -a † l,p a l,p S x l,p ≈ S 2 a † l,p + a l,p . (5.14)
Here a l,p and a † l,p are bosonic operators, satisfying [a l,p , a † l,p ] = 1 and [a ( †) l,p , a ( †) l,p ] = 0.

Harmonic Hamiltonian for ordered ice-rule states

The Hamiltonian is then expanded up to quadratic order in the bosonic operators (the linear terms vanish by definition of the angle ϑ, chosen to minimize the classical energy). In the following, we will concentrate the discussion on reference classical states which in zero field are ice-rule states with long-range order, whose ordered structure allows to build a spin-wave theory. We introduce the number ν 1 of frustrated J 1 -bonds and ν 2 of frustrated J 2 -bonds on a vertex. In all reference ice-rule states we consider these numbers are the same for each vertex. For an ice state ν 1 + ν 2 = 2, and ν 2 can only take two values : either ν 2 = 0 (for the collinear vertices) or ν 2 = 2 (for the Néel vertices). In the following, we will consider states containing only Néel vertices or only collinear vertices. In zero field, the energy of such a state writes

ε cl = -S 2 [J 2 -ν 2 (J 2 -J 1 )] 1 +
Hnnn = 2S (1 -ν 2 ) cos 2 ϑ l,p a † l,p a l,p + S 2 sin 2 ϑ lp,l p a † l,p a † l ,p + a † l,p a l ,p + h.c. HΓ = -Γ sin ϑ l,p a † l,p a l,p (5.16)
Remarkably, the spin-wave Hamiltonian depends only on the frustration parameters ν 1 and ν 2 , while it is completely independent of the other details of the spin configuration of the unit cell. The frustration parameters ν 1 ,ν 2 distinguish among Néel states and collinear states, but they are not able to distinguish among different collinear states. Even more remarkably, at the square-ice point J 1 = J 2 the dependence on the frustration parameter drops. Therefore, as we will discuss further, quantum corrections at the harmonic level are not able to lift the degeneracy between ordered ice-rule states, regardless of the size of their magnetic unit cell. This result can be extended even to disordered ice-rule states, which can be regarded as ordered ones with an infinite unit cell.

To diagonalize the spin-wave Hamiltonian of equation (5.16) we first introduce the Fourier transform of the bosonic operators and then perform a Bogoliubov transformation as described in Appendix B. The Hamiltonian then becomes

H = N ε cl - ε 0 2 + k,p ω k,p b † k,p b k,p + 1 2
(5.17

)
where

ε 0 = 2S [J 2 -ν 2 (J 2 -J 1 )
] is the zero field energy per spin of the state under consideration.

Néel state

Let us first consider the Néel state (ν 2 = 2), defined for Γ < 2S (2J 1 -J 2 ). Its unit cell contains n = 2 spins (as described in figure 3.5). The diagonalization of the spin-wave Hamiltonian shows that the spectrum of the magnon excitations is gapped whenever the classical Néel state is definedi.e. if Γ < 2S (2J 1 -J 2 ). Its lower band is plotted in figure 5.6. It has minima at (0, 0) and at the four corners of the Brillouin zone. This corresponds to the structure of the classical Néel state.

Collinear states

As already mentioned in section 5.2.2, all collinear states admit the same frustration parameter ν 2 = 0, and hence the same spin-wave Hamiltonian. This means that they possess the same spectrum of harmonic spin-wave excitations (but folded into a smaller Brillouin zone, the larger the unit cell), and that zero-point quantum fluctuations cannot lift the degeneracy among them. We will then specify the discussion to the particular collinear state, the "Manhattan state", represented earlier. Its unit cell contains n = 8 spins. While not being the simplest of all collinear states, this state is relevant because it can be energetically stabilized against other collinear states by e.g. dipolar interactions, as we discussed in the first part of this thesis.

The magnon dispersion relation, obtained by diagonalizing the spin-wave Hamiltonian, is shown in figure 5.7. It shows a finite gap, and two lines of degenerate modes of minimum-energy along the axes of the first Brillouin zone of the magnetic lattice (1/8 of the Brillouin zone of the geometric lattice). These degenerate modes traveling with momentum (k x , ±k x ) for all k x values can be associated with deconfined monopole pairs, obtained by flipping a finite string of spins along a J 2 -diagonal of the checkerboard lattice. Given the degeneracy of all collinear states, not perturbed by quantum fluctuations, these pairs are deconfined along the J 2 -diagonals, and their energy is independent of momentum as long as it satisfies the constraint of diagonal motion.

ω min /J 1 J 2 /J 1 (ω q -ω min )/J 1 J 2 /J 1 =1.01 J 2 /J 1 =1.1 J 2 /J 1 =1.5 J 2 /J 1 =2 - - 0 9

Quantum square ice

At the square ice point J 2 = J 1 = J, the spin-wave expansion around any canted ice-rule state produces the same excitation spectrum and zero-point fluctuations. Hence harmonic fluctuations are not able to lift the classical degeneracy of the ice states. In particular, the elementary excitations are the same as in the classical case, namely deconfined monopole pairs moving with arbitrary momentum. This shows up in the spin-wave dispersion, whose lower band is perfectly flat, and whose first excited band has a gap that equals that of a classical spin flip ∆=4JS 2 .

Polarized states

For large Γ, the classical reference state is fully polarized along the field. Quantum fluctuations however reduces the total polarization, as the magnetization does not commute with the Ising part of the Hamiltonian. We build the spin-wave expansion of the Hamiltonian around the classical polarized state in a similar way to what we did for the ice-rule states. Even though the classical polarized state is uniform, the spin-wave unit cell contains two sites (exactly as in the case of the Néel state), due to the fact that the checkerboard lattice is not a Bravais lattice. The bosonic excitations correspond to deviations of the spins from full polarization along the x axis.

The magnon spectrum is shown in figure 5.8. It displays softer modes at the four corners of the squarelattice Brillouin zone (for J 1 > J 2 ), and along the edges of the checkerboard-lattice Brillouin zone (for J 1 < J 2 ). These modes become gapless when approaching the critical field Γ c , signaling the instability of the fully polarized state to a Néel state (for J 1 > J 2 ) and to degenerate collinear states (for J 1 < J 2 ).

Results of the spin wave analysis

Linear spin-wave theory gives easily access to the the internal energy E = H LS W and to the order parameter m = S z /S . Their expressions are derived in Appendix C. Making use of these two quantities, the phase diagram of our system can be constructed in the harmonic approximation.

ΓS/J 1 J 2 /J 1 =0.7 -ΓS/J 1 =1.30 J 2 /J 1 =0.7 -ΓS/J 1 =1.43 J 2 /J 1 =1.3 -ΓS/J 1 =1.43 J 2 /J 1 =1.3 -ΓS/J 1 =1.30

Construction of the phase diagram

For each value of J 2 we determine the excitation spectrum, the internal energy and the order parameter of each of the three reference states (Néel, collinear and polarized). The reference state displaying the minimum energy as well as a positive order parameter, is then identified as the ground sate of the quadratic Hamiltonian.

Classically, the Néel and the collinear states can be distinguished by the correlation between the nextto-nearest neighbors (i.e. connected by a J 2 link). In the classical limit, C (2) = S z i S z j i, j = (ν 2 -1)S 2 : for the Néel state C (2) Néel = S 2 , whereas for all collinear states C (2) collinear = -S 2 . This property provides a further criterion to determine the ground state : a reference state can be considered as a stable ground state only if harmonic quantum fluctuations do not alter the sign of C 2 compared with the classical case.

A last criterion is the real value of the frequencies. Indeed, in principle the spin-wave expansion can produce imaginary or complex frequencies, but those signal the instability of the phase.

The phase diagram (figure 5.9) is then constructed as follows. The stability of each phase is verified by determining its order parameter, the sign of its C (2) correlations and its energy. The ground state is then chosen as the less energetic of the stable states. If all states fail to meet the stability criterion the fluctuations become too strong and the spin wave theory breaks down, suggesting that long-range order is absent in the true ground state.

Ground state phase diagram

A sizable region of the above-cited instability of the spin wave theory is found around the classical phase boundary between the polarized state and each of the ordered states. This corresponds to the black region in figure 5.9, where all three reference phases have a negative order parameter, as can be seen in figure 5.10, hence signaling the onset of a quantum disordered state.

However, a number of anomalies appear in the linear spin wave theory close to the transition to the field-induced paramagnetic phase. These suggest a significantly larger region of the phase diagram where long-range order may be lost because of quantum fluctuations.

First of all, we observe a qualitative deviation from the classical limit when approaching the square ice limit (J 1 = J 2 ) in a strong field (0.7 Γ J 1 ). In this range (indicated by the dense-hatched region in figure 5.9), the energy hierarchy between the Néel and collinear states is reversed with respect to the classical

Polarized

Néel Collinear case. This occurs despite the fact that quantum fluctuations are stronger for the energetically favored phase, so that its order parameter is smaller (namely m Néel > m collinear when E Néel > E collinear , and viceversa, see figure 5.11). This strong quantum effect of energy hierarchy inversion suggests that classical order might be unstable around the hierarchy inversion region in figure 5.9 when considering quantum fluctuations beyond linear spin-wave theory. The real ground state of the system may then be an intermediate phase which cannot be described within the linear spin-wave approximation.

Polarized Néel Collinear J 2 /J 1 ΓS/J 1 J 2 /J 1 ΓS/J 1 m i =0 • C (2) sign • E hierarchy • m x =1
ΓS/J 1 ΓS/J 1 ΓS/J 1 3 0 1 2 3 /J 1 J 2 /J 1 =2 1.99 2
Another strong quantum effect is also revealed close to the classical phase boundaries. While classically E Néel and E collinear are monotonic functions of J 2 /J 1 , they become non-monotonic around the above mentioned field range (figure 5.11). In particular E Néel grows with increasing J 2 /J 1 until it reaches a maximum, beyond which it starts to decrease; from a classical point of view this is quite surprising. According to the Hellmann-Feynman theorem, the next nearest neighbor correlations C (2) are given by the derivative of the energy with respect to J 2 , namely C (2) = ∂H/∂J 2 . Consequently, a change of sign in the derivative of E corresponds to a change of sign in C (2) , which means that the harmonic ground state is dramatically different from the reference state. The locus of the maxima in the energy as a function of J 2 /J 1 represents the lower bound of the hatched region in figure 5.9 . A further element of inconsistency of spin-wave theory is offered by looking at the transverse magnetization, given by the Hellman-Feynman theorem

m x = - ∂H ∂Γ . (5.18)
This quantity increases with the transverse field up to a transverse field Γ 1 where it saturates to 1 (and even becomes larger than one, see figure 5.12). This value of the transverse field is very close to the boundary of the region of inversion of the C (2) correlations, as represented in the phase diagram 5.9. This unphysical magnetization signals again that an ordered reference state does not lead to consistent results. On the other hand, the results obtained for the polarized phase (Γ J 1 ) do not present the same problem. However, they do not connect to the low-field results. In particular, the transverse magnetization is found to vanish for Γ ≥ Γ c (while it saturates to 1 for Γ ≤ Γ c ).

If one instead considers the projection of the order parameter onto the transverse field mx = m sin ϑ (neglecting 1/S corrections to the canting angle [START_REF] Zhitomirsky | Magnetization curve of a square-lattice heisenberg antiferromagnet[END_REF][START_REF] Coletta | Semiclassical approach to ground-state properties of hard-core bosons in two dimensions[END_REF]), the results on both sides of Γ c connect more smoothly. The values thus obtained are very close to the actual magnetization (obtained by derivation of the energy) in the (classical) low-field limit and in the polarized phase (where they are strictly identical). However, this estimate is incorrect as it violates Hellmann-Feynman theorem (5.18) and leads to unphysical results (negative susceptibility) in a range of fields realizing unphysical values (larger than the saturation value) for the transverse magnetization estimated as the derivative of the energy.

ΓS/J 1 ΓS/J 1 ΓS/J 1 0 0.5 1 J 2 /J 1 =0 J 2 /J 1 =0.98 J 2 /J 1 =1 0 0.5 1 0 1 2 3 J 2 /J 1 =1.02 0 1 2 3 J 2 /J 1 =1.1 0 1 2 3 J 2 /J 1 =2 m x m x Figure 5
.12: Field-induced magnetization as a function of the field for different values of J 2 /J 1 .

The highlighted regions are as in figure 5.10. The dotted lines correspond to the projection of the order parameter onto the field direction. These curves deviate from the average induced magnetization, obtained from the derivatives of the average energies with respect to the magnetic field Γ (see text).

Conclusion

In the transverse-field Ising model on the checkerboard lattice harmonic quantum fluctuations are shown to be unable to lift the classical degeneracy, be it exponential in system size at the square ice point J 1 = J 2 , or exponential in the system linear size when J 2 > J 1 . This implies that spin wave theory is inconsistent in determining which classical ground state is selected by quantum effects, and that non-linear effects have to play a central role in the lifting of the degeneracy and the selection of the exact ground state. The study of non-linear quantum effects in the case J 2 = J 2 will be developed in the next chapter.

Chapter 6

Lattice gauge theories in the transverse field Ising model on the checkerboard lattice

In the present chapter, we study the effect of a small transverse field on the Ising model on the checkerboard lattice, using perturbation theory. The effective Hamiltonian restricted to ice-rule states is endowed with a local (gauge) symmetry, stemming from the ice rules (which are equivalent to Gauss's for an emergent lattice gauge field). Hence perturbation theory produces a family of compact lattice gauge theories for S = 1/2 variables. The first non-trivial term appears at 4-th order, and it provides a realization of compact quantum electrodynamics (QED) with S = 1/2 variables [START_REF] Hermele | Pyrochlore photons: The u(1) spin liquid in a s = 1 2 three-dimensional frustrated magnet[END_REF][START_REF] Moessner | Planar pyrochlore, quantum ice and sliding ice[END_REF][START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF], also known as the U(1) quantum link model in the high-energy physics literature [START_REF] Chandrasekharan | Quantum link models: A discrete approach to gauge theories[END_REF][START_REF] Orland | Lattice gauge magnets -local isospin from spin[END_REF]. When pushing the perturbation expansion up to the eighth order, a diagonal term that favors a selected ice-rule state, namely the Néel states, appears. As we will see in chapter 8, the competition between the 4-th order term and higher order terms leads to a complex phase diagram, including at least two symmetry breaking phases.

We will now derive the perturbation Hamiltonian for our system and show how it compares to lattice gauge theories.

Degenerate perturbation theory

In zero field, the ground state manifold of the transverse-field Ising model is the ice-rule manifold. In presence of a transverse field Γ, we construct the degenerate perturbation theory on this manifold, building an effective Hamiltonian as a power series in Γ.

Principles

Let us first discuss the general formalism of degenerate perturbation theory (DPT). We start with a Hamiltonian H 0 with a degenerate ground state of energy ε 0 . We note P 0 the projector onto the groundstate manifold G, so that P 0 H 0 = H 0 P 0 = ε 0 P 0 (6.1)

A perturbation V is then added to the Hamiltonian, that does not commute with H 0 . The total Hamiltonian writes

H λ = H 0 + λV (6.2)
with λ a small dimensionless parameter. In the presence of the perturbation, the ground states |ψ (0) i of H 0 , of energies ε 0 , evolve into the eigenstates |ψ (λ) i of H λ . We note P λ the projector onto those perturbed states. The idea of DPT is to derive an effective Hamiltonian H eff that is restricted to G, whose eigenstates | ψ(λ)

i have energies ε (λ) i which reproduce the exact eigenenergies, and are as close to the actual eigenstates as possible |ψ (λ) i (i.e. minimizing |ψ (λ) i -| ψ(λ) i ). This effective Hamiltonian is usually expressed as a power series in λ of operators acting on G.

Resolvent and projector onto the perturbed states

We use Kato's resolvent method [START_REF] Kato | On the convergence of the perturbation method .1[END_REF][START_REF] Kato | On the convergence of the perturbation method[END_REF] to obtain the projector onto the perturbed eigenstates P λ .T o this end, we introduce the resolvent R λ defined as

R λ (z) = (H λ -z) -1 = (H 0 + λV -z) -1 (6.3)
so that

P λ = - 1 2iπ C dzR λ (z) (6.4)
with the contour C chosen to enclose only the eigenvalues of the perturbed states ψ (λ) i . This method allows to express P λ in terms of the unperturbed projector P 0 . The resolvent is first expressed as a power series in

λ R λ = -R 0 n≥0 λ n (-1) n [VR 0 ] n . (6.5) Writing (H 0 -z) = (H 0 -ε 0 ) 1 -(ε 0 -z)(ε 0 -H 0 ) -1
and taking the inverse, one then gets

R 0 = - n≥0 (ε 0 -z) n (ε 0 -H 0 ) -1-n . (6.6)
We introduce then

S 0 = -P 0 S n = (1 -P 0 )(ε 0 -H 0 ) -1 n (6.7) so that R 0 = (ε 0 -z) -1 P 0 + (1 -P 0 )R 0 = - n≥0 (ε 0 -z) n-1 S n (6.8)
Injecting this expression into the contour integral gives the projector onto the perturbed states

P λ = - n≥0 λ n n k 0 ,...,k n =1 k 0 +...+kn=n S k 0 V ...VS k n ≡- n≥0
λ n P n . (6.9)

Des Cloizeaux's expansion

We now need to construct the effective perturbed eigenstates | ψ(λ) i . There are many ways to construct theses states [START_REF] Klein | Degenerate perturbation-theory[END_REF], and they can lead to slightly different Hamiltonians. We choose des Cloizeaux prescription [START_REF] Cloizeaux | Extension d'une formule de Lagrange à des problèmes de valeurs propres[END_REF][START_REF] Takahashi | Half-filled hubbard-model at low-temperaturealf-filled hubbard-model at lowtemperature[END_REF], that symmetrically eliminates the overlap between these states. In this case, the effective perturbed states are obtained from the unperturbed states with the transformation Γ λ defined as

| ψ(λ) i =Π λ |ψ (0) i Π λ = P λ P 0 (P 0 P λ P 0 ) -1/2 . (6.10)
This transformation conserves the norm (Π † λ Π λ = 1 so that ψλ

i | ψλ i = ψ λ i |ψ λ i ).
It connects an unperturbed state (hence the right P 0 factor) to a perturbed state (with the left P λ factor).

With this prescription, the effective Hamiltonian is

H eff =Π † λ H λ Π λ . (6.11)
The remaining step is the expansion of H eff in powers of λ. To do this, we use the series expansion of (1x -1/2 ) (P 0 P λ P 0 ) -1/2 = P 0 + ∞ n=1 (2n -1)!! (2n)!! [P 0 (P 0 -P λ ) P 0 ] n . (6.12)

Expansion of the effective Hamiltonian

The generic expansion of the effective Hamitonian contains terms of the form P 0 VS k 1 V ...VS k n VP 0 , with k 1 + ...+ k n = n, with a numerical prefactor a k λ n+1 . Following Klein notation [START_REF] Klein | Degenerate perturbation-theory[END_REF], we define

(k 1 ,...,k n ) ≡P 0 VS k 1 V ...VS k n VP 0 .
(6.13)

Effective Hamiltonian for the TFIM

We now apply these results to the Hamiltonian of the transverse field Ising model on the checkerboard lattice (4.13). We consider now identical couplings J 1 = J 2 ≡ J, so that the Hamiltonian reads

H = 4J i, j Ŝ z i Ŝ z j -2Γ i Ŝ x i = J i, j σz i σz j -Γ i σx i . (6.14)
We will take the spin operators as Pauli matrices.

The ground-state manifold G is the ice-rule manifold, the perturbation is the transverse-field term V = i S x i , and the perturbation parameter is λ ≡ Γ. The terms of the form (6.13) become

(k 1 ,...,k n ) = (-1) n i 1 ,...,i n P 0 S x i 1 S k 1 S x i 2 ...S x i n S k n S x i n+1 P 0 . (6.15)
Such term flips the spins i 1 ,...,i n . If n is odd, it vanishes; indeed, ice-rule states are connected by flips of closed loops of alternating spin orientation and it is impossible to construct such a loop with an odd number of spins. Furthermore, if there exist p such that k p = 0, then (k 1 ,...,k p-1 , 0, k p+1 ,...,k n ) = (k 1 ,...,k p-1 )(k p+1 ,...,k p-1 ) . (

Therefore, (k 1 ,...,k p-1 ) and (k p+1 ,...,k p-1 ) have to satisfy the same role of containing an even number of spin flips. The first few of these terms are order 2 : (1) order 4 : ( 111), ( 102), (201) order 6 : (11111), ( 20111), ( 11201), ( 11102), (

Most of these terms have no contribution either. Indeed, any term containing a zero will factorize as in (6.16).

Because of the translational invariance, any term (6.15) is at least extensive. Therefore any factorized term will be non-extensive and would lead to a non-extensive contribution to the energy, and the contribution of all the terms containing a k p = 0 term must necessarily add up to zero, as we have explicitly verified up to n = 6

The only remaining terms are those of the form k p ≡ 1, ∀p. For those terms, the prefactor is simply a 1...1 ≡ 1. The second order term, coming from (1) = P 0 V ...VP 0 VP 0 is simply the classical correction to the energy :

H (2) = - Γ 2 ∆ N = -NJ sin 2 ϑ. (6.18)
Where ∆=4J is the spin-flip gap, and ϑ is the classical canting angle of the spins towards the transverse field, as discussed in 5.2. We will drop this term alongside with all other similar trivial terms that correspond to virtual flip of an open string of spins (namely terms that flip and unflip the same string of spins) as they only lead to a global shift in the energy of all unperturbed ground states. The next two sections are devoted to the non-trivial higher-order terms.

Order 4 perturbation theory : quantum link model

The first non-trivial term appears at fourth order. It corresponds to the minimal move connecting two different ice-rule state, which is the flip of a Néel plaquette (represented in figure 6.1). The effective Hamiltonian to this order is then

H (4) = -Γ 4 i 1 ,i 2 ,i 3 ,i 4 P 0 S x i 1 S 1 S x i 2 S 1 S x i 3 S 1 S x i 4 P 0 ≡-a 4,4 Γ 4 ∆ 3 F (6.19)
where

F

Resulting Hamiltonian

The fourth order perturbation expansion leads to a model of compact lattice QED (fcQED) for a discrete (S = 1/2) gauge field [START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF][START_REF] Hermele | Pyrochlore photons: The u(1) spin liquid in a s = 1 2 three-dimensional frustrated magnet[END_REF] -also known as U(1) quantum link model or U(1) gauge magnet in the high-energy physics literature [START_REF] Chandrasekharan | Quantum link models: A discrete approach to gauge theories[END_REF][START_REF] Orland | Lattice gauge magnets -local isospin from spin[END_REF]):

H (4) = -K 4 F + const. ( 6.20)
The coupling constant here has value

K 4 = 20Γ 4 /∆ 3 = (5/2)Γ 4 /J 3 .
fcQED can be suspected to undergo confinement due to the Polyakov mechanism (as described in section 1.3), valid for compact QED in d = 2 [START_REF] Polyakov | Gauge fields and strings[END_REF][START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF]. This expectation is indeed verified by numerics, using exact diagonalization [START_REF] Shannon | Cyclic exchange, isolated states, and spinon deconfinement in an xxz heisenberg model on the checkerboard lattice[END_REF][START_REF] Banerjee | The (2 + 1) -dU(1) Quantum Link Model Masquerading as Deconfined Criticality[END_REF] or quantum Monte Carlo [START_REF] Syljuåsen | Resonating plaquette phase of a quantum six-vertex model[END_REF], finding that the Coulomb phase is removed from the ground state of fcQED in favor of a gapped plaquette valence-bond solid (pVBS -see figure 6.2). This phase breaks the translational symmetry but not the spin inversion symmetry. A caricature of this state has half of the plaquettes in a resonating state between both Néel configurations (|ψ

= |N + | N / √ 2) as repsresented in figure 6.2.
As we will see in the next section, considering a field of arbitrary strength leads to a richer physics, related to a more complex emerging lattice gauge theory. Figure 6.2: Resonating plaquette valence-bond-solid (pVBS). Half of the plaquette are in a resonating state between both Néel configurations. This phase breaks translational symmetry but not the spin inversion symmetry.

6-th and 8-th perturbation theory

The quantum link model is expected to correctly predict the ground-state and low-temperature properties of the TFIM for small transverse fields. In order to get insight on the physics of the system at a larger field, we push the perturbation up to 8th order. At this order the first non-trivial diagonal term appears in the effective Hailtonian. Such a term turns out to be essential to understand the results from quantum Monte Carlo simulations of the TFIM, as discussed in chapter 8.

6-th order terms

Similarly to the 4-th order case, the 6-th order perturbation terms correspond to flips of closed loops. This order of perturbation gives a correction to the plaquette term and yield another loop move, namely the flip of a loop of 6 spins, as described in appendix D. The 6th order part of the Hamiltonian is then of the form H (6) /∆=-a 6,4 γ 6 Fa 6,6 γ 6 l∈L 6 F 6l + const. (6.21) where F 6l is an operator that flips the loop l if it belongs to L 6 (which is the set of all closed loop of 6 spins that connects to another ice-rule state) and vanishes otherwise. For the sake of compactness, we have introduced perturbation parameter γ =Γ/∆.

8th order terms

The 8-th order terms introduce corrections to the plaquette and 6-loop moves obtained at lower order, and a new term that corresponds to the flip of a closed loop of 8 spins. Additionally, a diagonal term appears. This term corresponds to the forward and backward flip of a plaquette. In particular, this term favors flippable plaquettes (i.e. aNéel plaquette, that can be flipped without leaving the ice-rule manifold). This move is diagonal and yields a finite correction for all plaquette configurations (see appendix D for more details). This favors Néel order, which is the state with the most flippable plaquettes [START_REF] Moessner | Ising models of quantum frustration[END_REF].

The corresponding Hamiltonian is then of the form

H (8) /∆=-a 8,4 γ 8 F -a 8,6 γ 8 l∈L 6 F 6l -a 8,8 γ 8 l∈L 8 F 8l -a 8,8 F 2 -a 8,8
F F + const. (6.22) Figure 6.3: Possible non trivial eighth-order moves. It can be the flip of a closed loop of length 8 (black), the flip of two plaquettes sharing as spin (red) or a vertex (green), or the flip of two distant plaquettes(yellow). It can also be the virtual flip of a loop (blue), which is a diagonal term that gives a non-zero contribution for any plaquette configuration.

Total Hamiltonian

Putting all the terms together, the total Hamiltonian has the form

H eff /∆=-K 4 F -K 6 l∈L 6 F 6l -K 8 l∈L 8 F 8l -K 8 F 2 -K 8 F F + const. + O(γ 10 ) . (6.23)
Let us consider the model obtained by retaining in the Hamiltonian only dominant off-diagonal and the diagonal terms :

H = -K 4 F -K 8 F 2 (6.24)
This has been studied in Refs. [START_REF] Shannon | Cyclic exchange, isolated states, and spinon deconfinement in an xxz heisenberg model on the checkerboard lattice[END_REF][START_REF] Banerjee | The (2 + 1) -dU(1) Quantum Link Model Masquerading as Deconfined Criticality[END_REF] using exact diagonalization, and a transition from Néel to pVBS is predicted to occur for a critical ratio K 8 /K 4 ≈ 0.37. The ratio between these two coefficients can be controlled in the TFIM via the transverse field, 0.7 ≈ (a 8,8 /a 4,4 )(Γ/∆) 4 . The corresponding value of the field is Γ/J ≈ 0.64. We will return to this point in chapter 8, comparing with numerics on the TFIM.

From the quantum link model to the frustrated compact quantum electrodynamics

The quantum link model can be written as a frustrated compact quantum electrodynamics (described in 1.4.2). Starting from the spin operators, we introduce an angular operator Âi with compact eigenvalues, so that

Ŝ + i = e iA Ŝ - i = e -iA † Ŝ z i , Â j = iδ ij (6.25)
Using these operators, (6.20) becomes

H (4) = -2K 4 cos Âi 1 -Âi 2 + Âi 3 -Âi 4 = -2K 4 cos Φ (6.26)
where the spins i 1 , i 2 , i 3 , i 4 are the four spins of a plaquette, labelled clockwise, and Φ is the lattice curl of Â, as defined in (1.20). This Hamiltonian has precisely the form of that of a compact QED, and within the ice-rule manifold, (lattice) Gauss's law is verified by Ŝ , so that the effective gauge theory is a frustrated compact QED, as introduced in section 1.4.2.

Gauge mean field theory

Degenerate perturbation theory constructs a "pure" gauge theory. Indeed, by construction the presence of defects is not allowed in this approach. The relationship between quantum spin ice and lattice gauge theories can be further fleshed out using gauge mean-field theory (gMFT) that, on the other hand, provides a matter field representing the monopoles. This method was recently developed to study the low-temperature phase diagram of pyrochlore spin ice and explain the properties of Yb 2 Ti 2 O 7 in terms of that of a quantum spin ice [START_REF] Savary | Spin liquid regimes at nonzero temperature in quantum spin ice[END_REF]. The idea is to artificially separate the quantum dynamics in the ice-rule manifold from that of the defects. The former one is treated as a gauge field and the latter as a matter field.

Principles

In the gMFT approach, the S = 1/2 degrees of freedom are formally split into a "matter" part, and a gauge part. The matter part is the spinon (monopole) field living on the centers r of the vertices, which is represented by a bosonic field of integer modulus Φ r = e iφ r . The gauge part is a S = 1/2 gauge field s α rr , with α = x, y, z, living on the lattice sites between two vertices r and r (see figure 6.4).

S z i → s z i S + i → Φ † r s + rr Φ r (6.27)
with i the site shared by the vertices r and r . We introduce the charge operator Q r = r |r,r s z rr of vertex r, where the sum runs over all adjacent sites of the vertex. The operators Q r and Φ r obey the commutation relations

Φ r = e iφ r φ r , Q r = i . (6.28)
In principle, Q r = 0, ±1, ±2. The Hilbert space can be enlarge by allowing it to take any integer value, provided that the fluctuation in Q r around 0 are sufficiently small. Using these notations, the Hamiltonian of the TFIM (4.13) becomes

H = 2J r Q 2 r -Γ r,r Φ † r s x rr Φ r . (6.29)
Φ r Φ r s rr α Figure 6.4: gMFT parameters. The matter field Φ r lives on the vertices lattice and the gauge field s α rr lives on the sites of the lattice, between the r and r vertices.

Mean-field decoupling

The dynamics of the gauge field and of the spinon field are then decoupled in a mean-field fashion :

S x rr =Φ † r s x rr Φ r ≈Φ † r Φ r s x rr + s x rr Φ † r Φ r -s x rr Φ † r Φ r (6.30)
so that the Hamiltonian reads

H≈H Φ + H s + const. (6.31)
with

H Φ = -Γ rr s x rr cos(φ r -φ r ) + 2JS 2 r Q 2 r (6.32) H s = -Γ rr cos(φ r -φ r ) s x rr . (6.33)
Such decoupling implies that the gauge field cannot confine the matter field. Therefore the gMFT description can only be valid at finite (albeit low) temperatures, namely in the deconfined phase.

Confinement/deconfinement transition

The spin Hamiltonian is readily minimized by s x = 1/2, so that the spinon Hamiltonian reduces to a quantum rotor Hamiltonian on the square lattice

H Φ = - Γ 2 rr cos(φ r -φ r ) - J 2 r ∂ 2 ∂φ r 2 (6.34)
This Hamiltonian admits a numerical treatment using path-integral Monte Carlo [START_REF] Wallin | Superconductor-insulator transition in two-dimensional dirty boson systems[END_REF], and presents a Mott insulator/superfluid transition. As long as we are in the Mott regime, the fluctuations in Q r are limited. However, in the superfluid phase, these fluctuations become too strong and this theory breaks down. In particular, the enlargement of the spinon Hilbert space ceases to be justified. On the other hand, the quantum paramagnet phase could be identified with the spinon condensed phase. But, as we will see, the model under consideration exhibits a confining phase at intermediate fields, separating the spinon condensate phase from the deconfined phase. Therefore the gMFT can locate the confined phase even though it is not valid in this regime. A readily accessible quantity from the quantum rotor model is the kinetic energy, corresponding to the the transverse magnetization. It is given by It compares rather well to the linear-spin-wave results at low field, for which both tend to the classical value of the magnetization (m x = sin ϑ =Γ/2J), also given by second order perturbation theory (6.18).

S x ≈s x cos(φ r -φ r ) = 1 2 cos(φ r -φ r )≡ 1 2 m x ( 6

Chapter 7

The membrane algorithm for quantum spin ice Frustrated spin models are known to be difficult to tackle numerically because of the so-called sign problem [START_REF] Loh | Sign problem in the numerical simulation of many-electron systems[END_REF][START_REF] Henelius | Sign problem in monte carlo simulations of frustrated quantum spin systems[END_REF]. However, frustrated Ising models only exhibit "classical" frustration, namely the frustration in one spin component only, and then they are not affected by the sign problem. We use world-line Monte Carlo. This approach introduces an additional dimension to the system accounting for quantum fluctuations, while the spins only interact through physical Ising interactions in the original spatial direction (see [START_REF] Scalettar | World-line quantum monte carlo[END_REF] for more details on the method). This allows to generalize the loop algorithms, that have proven to be so efficient in classical models [START_REF] Barkema | Monte carlo simulation of ice models[END_REF], to the quantum domain.

In this chapter, we describe this mapping alongside with the algorithm we developed.

Trotter-Suzuki mapping

When adding the transverse-field term to a d-dimensional Ising model, the Hamiltonian acquires an off-diagonal part in the S z i basis. All the difficulties of the problem lie then in the diagonalization of this Hamiltonian. The Trotter-Suzuki transformation offers a clever way to get an effective, but rather accurate, description of the model [START_REF] Trotter | Critical behavior of random transverse-field ising spin chains[END_REF][START_REF] Suzuki | Generalized trotter's formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems[END_REF]. In this section, we consider a generic d-dimensional transverse field Ising model with N spins, whose Hamiltonian reads

H = i, j J ij σz i σz j -Γ i σx i . (7.1)
We will note H 0 = i, j J ij S z i S z j the Ising part of the Hamiltonian (which is by definition diagonal in the S z i basis), and H Γ =i S x i the transverse-field part.

Principles

Quantum Monte Carlo simulations aim at sampling the (quantum) partition function of the model under investigation. Through the Trotter-Suzuki mapping, we construct an effective classical system that shares the same partition function as the original quantum model. The partition function of the TFIM is given by

Z = Tr e -β Ĥ = Tr e -β Ĥ0 +βΓ ĤΓ . ( 7.2)
The key point of the method is the Trotter formula [START_REF] Trotter | Critical behavior of random transverse-field ising spin chains[END_REF][START_REF] Suzuki | Generalized trotter's formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems[END_REF], that gives the exponential of the sum of two non-commuting operators A and B :

e x(A+B) = lim M→∞ e x M A e x M B M (7.3) In particular, if M is finite, e x M A e x M B M = e x M (A+B) +O x 2 M 2 M = e x(A+B) +O x 2 M . (7.4)
We will use this formula to express the partition function of our d-dimensional quantum model as a (d + 1)dimensional classical partition function. This approach is fairly general, but in the case of transverse-field Ising models, the corresponding classical model takes the form of a (d + 1) Ising model.

Application to transverse-field Ising models

Trotter dir. The intralayer coupling (solid lines) becomes J/M, and the interlayer coupling (dashed lines) is λ(β, ) as defined in the text.

Effective partition function

We work in the S z basis, and we will note |S 0 ≡ S z 1,0 ,...,S z N,0 the basis vectors. We truncate the Trotter formula (7.3) at a finite M to obtain an approximation of the partition function. We get

Z≈ |S 0 S 0 | e -β M Ĥ0 e -βΓ M ĤΓ M |S 0 = |S p S p M p=1 e -β M Ĥ0 e -βΓ M ĤΓ S p + 1 (7.5)
where the last equation was obtained by introducing (M -1) completeness relations

S p S p S p = , with
|S M ≡ |S 0 . These M copies of the system correspond to the M layer of the effective classical system in the additional direction. We will call these layers Trotter layers and the corresponding direction the Trotter dimension. According to (7.4), the error induced by the finiteness of M in the expansion (7.5) is controlled by the Trotter parameter

= βΓ M . (7.6)
Given that H Γ is a sum of single-spin operators, we rewrite Z as :

Z = |S p M p=1 e -βE p (S ) N i=1 S z i,p e -Ŝ x i S z i,p+1 (7.7)
where

E p (S ) = 1 M ij J ij S z i,p S z j,
E(S) = ijp J ij M S z ip S z jp -λ ip S z ip S z ip+1 (7.12)
its energy. The total energy is naturally decomposed into an in-plane energy E and a transverse energy E ⊥ , which are defined as

E (S) = ijp J ij M S z ip S z jp E ⊥ (S) = -λ ip S z ip S z ip+1 ≡-λC z (S) . (7.13)
The Trotter dimension has been introduced to include quantum fluctuations. Indeed, fluctuations in the Trotter direction correspond to quantum fluctuations. The coupling in this direction λ diverges logarithmically to -∞ when T →∞or Γ → 0, which corresponds to the classical limit. In this case, all layers will be identical, and hence quantum fluctuations will be suppressed in the physical model. The other limit (λ → 0) corresponds to the zero temperature or infinite field limit of the physical model, in which quantum fluctuations are very large. In the effective system, the Trotter layers will be completely decoupled, and therefore fluctuations in the Trotter direction will be very large as well.

In the end the partition function Z is expressed in terms of the effective partition function Z of a classical Ising model with energy E(S):

Z = e MNν Z . (7.14)
Because of the e MNν prefactor, and because the coupling λ in the Trotter direction depends on β and Γ, thermodynamic observables of the quantum model (energy, heat capacity,. . . ) have to be extracted with care. Their expressions are derived in the next section.

Observables

We calculate the thermodynamic observables as derivatives of ln Z using (7.14).

Transverse magnetization

We note M x = Nm x the transverse magnetization of the system. Its expression is given by

Nm x ≡ 1 β ∂log Z ∂Γ = MN β ∂ν ∂Γ + 1 β ∂log Z ∂Γ . (7.15)
Using (7.10), we get

m x = M β ∂ν ∂Γ + ∂λ ∂Γ C z Z = 1 sinh(2) cosh(2) - 1 NM C z Z (7.16)
where . Z is the average of the classical model with partition function Z.

Internal energy

The mean energy E of the system is given by

E = Nε≡- ∂log Z ∂β = -MN ∂ν ∂β - ∂log Z ∂β . (7.17)
Using the expression for λ and ν, we get

Nε = E Z - NΓ sinh(2) cosh(2) - 1 MN C z Z . (7.18)
It it consistent with (7.16), so that

ε = 1 N E Z -Γm x . (7.19)
The diagonal part of the energy is simply obtained by averaging it on each of the layers.

Structure factor

As we have seen in chapter 4 in the case of the classical case, the features of the structure factor are characteristic properties of spin ice in their collective paramagnet phase. Therefore we are interested in determining it evolution with the transverse field. It is given by

S q = Z -1 Tr          1 N i e iq•r i σz i 2 e -β Ĥ          . (7.20)
We introduce then the M -1 completeness relations to get

S q = Z -1 |S p M p=1 1 N i e iq•r i S z i,p 0 2 e -βE p (S ) N i=1 S z i,p e -Ŝ x i S z i,p+1 = 1 N i e iq•r i S z i,p 0 2 Z . (7.21)
Because the system is invariant by translation in the Trotter direction, the layer p 0 can be chosen arbitrarily. We take then the average over all layers as the estimator of the structure factor :

S q = 1 M p 1 N i e iq•r i S z i,p 2 Z . (7.22)

Membrane algorithm for the TFIM

The d + 1 classical model with effective partition function Z is then studied using Monte Carlo. This allows to evaluate the . Z averages that appear in the expression for the thermodynamic quantities of the original d-dimensional quantum model. From now on, we restrict the Ising coupling J ij to that of the TFIM.

The precision of the calculations is set by the Trotter parameter. We observe that for ∼ 10 -2 -10 -3 the quantities of interest have converged to their → 0 value within statistical error bar (namely the Trotter error, i.e. the error induced by the truncation (7.4), is comparable or smaller than the statistical error). . For a given temperature and a given transverse field, fixing sets the Trotter dimension M of the system. Therefore, reducing the temperature increases the effective system size, and hence also computing times. We are aiming at the ground state properties of the system. However, we cannot reach very low temperature with reasonable Trotter error. Nevertheless, if a zero temperature phase survives up to sufficiently high temperature, it will be reachable with our simulations.

At low temperature and in zero transverse field, the relevant update moves in a Monte Carlo simulation of ice models are closed-loop flips, as discussed in section 3.4. We generalize this algorithm in the presence of a transverse field with the membrane algorithm, that flips clusters of spins obtained from the growth of a loop in the Trotter direction. p 0 ? ? ? ? Figure 7.2: Membrane update. First a loop is grown in layer p 0 . This loop can be closed (red) or not (blue). It is then propagated in the Trotter dimension, the loops in each layer being treated as the spins in a one dimensional Wolff algorithm.

Principle

The membrane algorithm works in three steps.

• First, a layer p 0 is chosen randomly. In this layer, a loop l p 0 = {S z i 1 ,p 0 ,...,S z i m ,p 0 } is created in the exact same way as in the classical case (as detailed in section 3.4). We use both short-and long-loop updates, with the possibility of open loops (or teleportation moves), as the temperature is not zero. By construction, flipping l p 0 cost no intralayer energy.

• Then the membrane is grown in the Trotter dimension. Starting from layer p 0 , the loop l p 0 propagates à la Wolff in the Trotter dimension (p = p 0 ± 1, p 0 ± 2,...), trying to flip the same loop in the successive neighboring layers. The interaction energy between the loops {S z i 1 ,p ,...,S z i m ,p } and {S z i 1 ,p+1 ,...,S z i m ,p+1 } of two successive layers p and p + 1 is

E(l p , l p+1 ) = -λ m α=1 S z i α ,p S z i α ,p+1 (7.23)
Hence the loop l p+1 is added to the membrane with a probability max 0, 1 -e 2βE(l p ,l p+1 ) .

• The layers are not all in the same configuration, hence propagating the loop through a given layer may create monopoles. Therefore the update is not microcanonical and costs intralayer energy. We accumulate the total intralayer energy variation ∆E of the flip of all the loops added to the membrane. Once the growth of the membrane stops, the membrane is then flipped with a probability p = min 1, exp(-β∆E ) . This last step guarantees that the algorithm satisfies the detailed balance condition.

This algorithm samples efficiently the low-energy states of the TFIM and hence allows to access the low-temperature properties of the model, as we will see in the next section.

Properties of the algorithm Acceptance rate

When βΓ 1, the acceptance rate of the move is close to 1. Indeed, λ diverges logarithmically as Γ → 0, and the system becomes perfectly correlated in the Trotter dimension. Therefore, as the loop is flippable in the starting layer, it must be so in all layers and the intralayer energy cost ∆E vanishes. For larger transverse fields, the system decorrelates in the Trotter direction, and hence the acceptance rate drops.

More precisely, let M be a membrane containing N M spins, and grown upon a loop L of length l L . Once it has been grown, M is flipped with a probability

P flip = min          1, exp          - 2βJ
M (i,k)∈M j∈N i σ i,k σ j,k                   (7.24)
where N i represents the set of lattice sites neighboring the site i, and the primed sum indicates that one has to exclude the sites belonging to the membrane. The probability P flip has value 1 in the classical limit → 0, λ →∞ , in which all the layers display the same configuration, and hence a microcanonical loop/string on a layer is equally microcanonical on every other layer. For a finite transverse field, on the other hand, the flip probability will be typically reduced due to the presence of discontinuities in the imaginary-time propagation, associated with defect vertices (namely monopoles) appearing in isolated layers. A naive estimate of the scaling of the membrane flip probability gives

P flip ∼ exp[-β(J/M)N M n m ] = exp[-βJln m ],
where n m is the density of (free) monopoles in the system. Such a scaling would imply that the probability is inevitably suppressed exponentially as the temperature is decreased. Yet a more careful estimate leads to a much slower decay of the acceptance rate. Indeed, at very low temperatures, βJ 1, the thermal monopole density n m is exponentially suppressed, while the monopoles induced by quantum fluctuations are bound, as discussed in chapter 8. Hence their effect on the suppression of the flip probability is not as simple as their density n m appearing in the previous scaling formula. In particular a simple estimate (coming from perturbation theory) of the typical size of a bound monopole pair gives l pair ∼|log(Γ/(2J))| -1 . We can therefore imagine that the flip probability of the membrane M will be affected by bound monopole pairs only if such monopole pairs cross the loop/string, hence if they fall within a region of size l L ×l pair . Making use of (7.16), we estimate the density of discontinuites (i.e. of field-induced pairs of monopoles) at low field as n pairs ∼ m x . This means that the exponential suppression of P flip due to bound monopole pairs can be estimated as P flip ∼ exp(-β(J/M)l L l pair Mm x ). If l L ∼ O(1) (using short loops), we find that the suppression is still exponential. However, there is a factor in the exponential that strongly reduces the decay with β of the acceptance rate compare to the previous naïve estimate.

The acceptance rate decreases with the length of the starting loop. Therefore in our calculations we used both long and short loops as starting points for the growth of membranes. The latter one creates shorter loops, whose average length is independent of the system size. Using this update assures that the algorithm retains a reasonable acceptance rate, even for large systems.

The long-loop algorithm is needed to ensure the sampling of different topological sectors. Indeed, although the short-loop algorithm can create winding loops, this will happen with a probability that drops exponentially with the system size. Therefore using only this kind of update would typically sample a single topological sector.

Figure 7.3 represents the acceptance rate of both algorithms as a function of the transverse field at T/J = 5 × 10 -3 and as a function of the temperature at Γ/J = 0.4, for different system sizes. As expected, the acceptance rate of the short-loop version of the algorithm depends weakly on the system size.

As we will see in chapter 8, the region 0.2 Γ/J 0.4 corresponds to the appearance of an ordered phase, well distinct from the Coulomb phase. Hence the drop in the acceptance rate of the membrane is due to the appearance of a strong energy separation between the ice-rules, similarly to what we observed in the Néel phase of the classical case in section 3.4.

The acceptance rate increases again for Γ/J 0.6, where the system enters a quantum paramagnetic phase. In this regime, the dominant update from membrane algorithm is the flip of very short open loops that are limited to one or only a few layers. The membrane algorithm reduces then to a single spin flip algorithm, which is sufficient to simulate efficiently the quantum paramagnet phase.

Size of the clusters

A large acceptance rate ensures that the algorithm updates the system. However, like any cluster algorithm, its strength lies in the possibility of flipping a large number of spin in a single move. In order to ensure the efficiency of the membrane algorithm, we have to check the size of the clusters it flips, both in real space dimensions as well as in the Trotter dimension.

T /J T /J Γ/J Γ/J In-layer loops In zero field, the loops generated by the membrane algorithm should be similar to the ones generated by the classical loop algorithms at the same temperature. We expect then the length of the long loops to scale as L 5/3 [START_REF] Barkema | Monte carlo simulation of ice models[END_REF], with L the linear system size.This is indeed the case, as can be seen in figure 3.2. Furthermore, as the transverse field is increased, the correlations in the Trotter direction decrease and the dependence of the acceptance rate on the length of the loop significantly increases (figure 7.3). Therefore we expect the average length of the loops to be reduced as the transverse field increases. This is exactly what we observe in figure 7.4.

Trotter direction

The average dimension l z of the membrane in the Trotter direction signals to what extent the algorithm is able to decorrelate the system in this direction. At small fields, the system is strongly correlated in the Trotter dimension, as expected in the classical limit. Therefore the membranes that the algorithm generates should span the entire Trotter dimension (i.e. l z ≡ M). This is indeed seen in figure 7.5.

At fixed temperature, the average length in the Trotter dimension of the membranes is almost the same for the long-loop and the short-loop algorithms. It is close to one for small fields and rapidly drops to zero at Γ ≈ 0.4J, corresponding to the onset of the quantum paramagnet phase.

Effective Hamiltonian : Quantum link model

The membrane algorithm can be adapted to other ice models. In particular, it can be used to study the quantum link model, described in chapter 6. The Hamiltonian is given by

H = -K 4 F (7.25)
with F the operator introduced in equation (6.19), that flips a plaquette of spins if it is in a Néel configuration and vanishes otherwise. The path-integral formalism of this model is slightly different from that of the TFIM, and we will describe it below.

Effective partition function

We can decompose the Hamiltonian into two parts F o = odd F and F e = even F , acting respectively on the odd and even plaquette sublattices, so that each term is a sum of terms that commute with each other. We proceed then as in the TFIM case and derive an effective expression for the partition function.

We note K = βK 4 , so that the partition function of the model reads There is no interaction within a single layer, but all layers must be in an ice-rule state. In a given even (resp. odd) layer p 0 , the even (resp. odd) plaquettes interact with the p 0 + 1 layer whereas the odd (resp. even) plaquettes interact with the p 0 -1 layer (see figure 7.6). Note that a spin configuration will only have a positive weight if all pairs of interacting plaquettes are either identical or flippable and antiparralel.

Z = Tr e -K(F e +F o ) . ( 7
w(, ) =                  sinh() if | ∈ |N , | N and | = | cosh() if | ∈ |N , | N and | = | 1 if | is
We now need an efficient way to sample the states of the effective model, compatible with the constraint of working within the ice-rule manifold.

Algorithm

We use two kinds of updates to simulate this model.

Long loop membranes

The first update is a long-loop membrane algorithm, along the line of what we described in section 7.2. It is essential for the correct sampling of the topological sectors. In the case at hand, long-loop membranes cannot have open ends in the Trotter dimension, unless they have a single plaquette as a support (see next paragraph). Indeed, two successive replicas cannot differ by the flip of a loop which extends over more than a single plaquette, as such a flip would require the flip of even as well as odd plaquettes, which is impossible in a single propagation step according to the Trotter-Suzuki decomposition of (7.27). Therefore long-loop membranes are automatically propagated along the whole Trotter dimension and their flip is accepted with a probability p = max {1, W(S )/W(S)}, where S is the initial configuration and S the configuration after the flip.

Plaquette membranes

The second and most important update is then a plaquette membrane update. This move introduces fluctuations in the Trotter direction. It is very similar to the membrane update of the TFIM, with a few crucial modifications to adapt it to the form of the interactions in the present model.

Let us consider an even layer 2p (the same reasoning holds for an odd layer), and let (1),2p , (2),2p , (3),2p and (4), 2 p be the four neighboring plaquettes of 2p within layer 2p. Flipping 2p without flipping 2p + 1 will flip one spin in each of the (i),2p and none in each of the (i),2p+1 , leading to a configuration with zero weight. The fundamental moves of this algorithm are then simultaneous flips of two non-interacting adjacent plaquettes. For each i ∈{1, 2, 3, 4}, if (i),2p and (i),2p+1 are flippable and in an opposite configuration (i),2p = (i),2p+1 , this move also creates a configuration with zero weight. If these plaquettes were identical in both layers, flipping 2p and 2p+1 will switch them between flippable and non-flippable, connecting two configurations with non-zero weight. Fortunately, the weight of having two interacting opposite plaquette configurations is sinh which is much smaller then the weights of identical plaquettes (1 or cosh ). These configurations will then be very rare and hence they do not prevent plaquette flips from being accepted at a sizable rate.

The plaquette membrane algorithm works like a regular one dimensional Wolff algorithm, where at each step two plaquettes are added to the cluster.

• First, randomly pick a layer p 0 and a plaquette p 0 in this layer. We will consider here the case where p 0 and are even. If p 0 is not flippable, stop the algorithm. Otherwise add p 0 and p 0 -1 to the membrane, and store the ratio δw p 0 = w( p 0 , p 0 +1 )/w( p 0 , p 0 +1 ) of the weights of the four neighboring pairs of plaquettes p 0 , p 0 +1 after the move by the same weights before the move. • Then propagates the flip of double plaquettes to layers p 0 + 1 and p 0 + 2 with a probability p = max 0, 1w( p 0 , p 0 +1 )/w( p 0 , p 0 +1 ) . and go on

• Once the propagation has stopped, flip the membrane with probability min

       1, p δw p        .

Properties of the algorithm

In the following we document the temperature-and size-dependence of the efficiency of the algorithm. To study the quantum link model, we used both the long-loop membrane and the plaquette membrane algorithms. This allows to sample efficiently the configurations of the system. We describe here this efficiency, which is characterized by the acceptance rate of the updates, the size of the membranes and the fluctuations between topological sectors.

Acceptance rate

Let us first consider the acceptance rate of the algorithm (figure 7 acceptance rate barely depends on the system linear size L or the coupling K. It is affected by the proportion of flippable plaquette in each layer(∼ cosh ), and the proportion of antiparallel pairs of interacting plaquettes (∼ tanh ), which are both function of only.

The long-loop membrane acceptance rate decays with K, but remains sizable ( 0.2) for a wide range of couplings. However, this acceptance rate drops at a maximal coupling that decreases with L (see figure 7.7). Indeed, in the ice-rule manifold the average length of a long-loop is proportional to L 5/3 , with L the linear system size. The size of the membrane would then scale as L 5/3 M = L 5/3 K/. With a constant density of antiparallel pairs of interacting plaquettes, the probability to hit one of them is then ∼ L 5/3 , therefore we expect a drop of the acceptance rate for K ∼ L -5/3 , which is what we observe.

Size of the membranes

In-layer loops The fact that the average length of long loops scales as l ll ∝L 5/3 is linked to the algebraically decaying correlations of the Coulomb phase. If the monopoles are confined, we expect this scaling to change, and the average length of the loops to become system-size independent. At low coupling, the classical L 5/3 scaling is recovered (see figure 7.8), whereas at large coupling the average loop length saturates at a constant value l ll ≈5 for any system size, as expected. At low couplings, we recover the classical limit and the whole Trotter dimension is updated at each move. As the temperature is decreased, this update becomes more and more local as the fluctuations set in. Topological sector fluctuations A fundamental property of the classical Coulomb phase is the presence of long-range fluctuations that connect different topological sectors. As discussed in section 1.2.3, the height of the pinch points observed in the structure factor is precisely given by the topological sector fluctuations. The system to undergo a transition to a confined phase [START_REF] Polyakov | Gauge fields and strings[END_REF] in which these fluctuations are suppressed. In the high temperature limit we recover the classical value of the pinch point height (∼ 1.88). As the temperature decreases, the topological sector fluctuations strongly decrease. For large system sizes, the fluctuations are abruptly suppressed when K reaches K * (L) and the long-loop membrane updates are suppressed. This effect is not physical, but it is uniquely related to the drop in the acceptance rate of the long-loop membrane move. The correct behavior of the winding number fluctuations can be recovered in principle by increasing the number of membrane moves per Monte Carlo step so as to maintain the number of accepted moves per MC step fixed even though this becomes prohibitive for sufficiently large systems. Low temperature thermodynamics of quantum square ice

In this chapter, we discuss the phase diagram of quantum square ice, realized by the transverse-field Ising model on the checkerboard lattice. The phase diagram is obtained by a combination of quantum Monte Carlo simulations on the transverse-field Ising model as well as on the effective Hamiltonian obtained from perturbation theory at 4-th order in the transverse field, namely the quantum link model. Our method allows to construct the non-perturbative phase diagram of the model, at finite temperature. We compare the numerical results to the predictions of the analytic approaches described in chapters 5 and 6.

Phase diagram of the quantum link model

Let us first analyze the model obtained from the fourth order perturbation theory, the quantum link model defined in section 6.2.1, whose hamiltonian is given by (6.20). This model is expected to undergo a confinement/deconfinement transition, corresponding to the onset of a low-temperature symmetry-breaking phase. In the following, we note K = βK 4 the coupling constant, where β is the inverse temperature.

Numerical results

Exact diagonalization results for the 4×4 system (detailed in appendix E) illustrate that the pVBS ground state of the model is to be found in the w = 0 sector. To treat larger system sizes, we make use of quantum Monte Carlo based on the algorithm described in section 7.3. For the following results, we set the Trotter parameter to = 10 -2 . To ensure convergence towards the correct topological sector, each set of data was obtained via an annealing procedure starting from K = 0.01, and approaching the desired temperature in steps of ∆K ≈ 0.1, consisting of 50000 MC steps each. For a system of linear size L, each MC step consists of L 2 /4 plaquette membrane flips and L long-loop membranes flips.

Coulomb to pVBS phase transition

As mentioned earlier (see section 1.2.3), a signature of the confinement of the system is the suppression of the fluctuations of its winding number (figure 7.10). These fluctuations decrease rapidly when the temperature is lowered, signaling the confinement of the excitations. By construction of the Hamiltonian, no monopoles are present in the system. At sufficiently high temperature, the system is disordered and fluctuates between the ice states. Therefore, in this regime pinch points are clearly visible in the magnetic structure factor (figure 8.1). As the temperature is lowered, the height of the pinch point, corresponding to the winding number fluctuations, is reduced, and the features of the structure factor related to ice-rule states are lost. On the other hand, at low temperature the structure factor develops peaks at the Néel points. These are not Bragg peaks, as they have a finite width. They are a consequence of the fact that the fluctuations concentrate in the w = 0 topological sector, and in particular many states in this sector are obtained from the flip of only a few plaquettes starting from a Néel state. In fact, in order to get a state with zero Néel order parameter out of a Néel state, one has to flip exactly half of the spins or half of the plaquettes of one of the plaquette sublattices. Therefore, even if the system is locally very close to a Néel state, the correlation length remains finite. We define the flippability f of a plaquette state | as

K =0.3958 K =0.01 K =1.0710 K =1.2639 K =2.0355 K =0.6852 K =0.8781 K =3 Figure 8
.1: Magnetic structure factor of the S z component of the spins for the quantum link model. At high temperature, the system is in a quantum Coulomb phase, as it is disordered while no monopoles are present (by construction of the Hamiltonian). Pinch points are indeed clearly visible. As the coupling is increased, the height of the pinch points (the winding number fluctuations) is progressively reduced while rounded peaks appear at (π, π).

f = | 2F 2 -1 | (8.1)
where F is the operator defined in (6.19) that flips the spins on plaquette . The flippability is 1 if the plaquette is in one of the two Néel states, and -1 otherwise. We then calculate the flippability structure factor. The pVBS phase (described in section 6.2.1) is characterized by a finite mean staggered flippability f 0,π and f π,0 . When the coupling K is increased (or the temperature is lowered) the system undergoes a transition from a Coulomb phase to a pVBS phase. This can be detected by tracking the flippability correlations at wavevectors (0,π) and (π, 0) through this transition. More precisely, we define the order parameter of the pVBS phase as

m 2 pV BS = (L/2) -4 , (-1) + f f (8.2)
To get an accurate estimate of the critical temperature T c , we calculate the Binder cumulant of the staggered flippability

U 4 = 1 - m 4 pV BS m 2 pV BS 2 (8.3)
From the crossing point of the Binder cumulant curves (figure 8.2), we extract the critical coupling K c = 1.38 (2), which corresponds to a transition temperature T c = 0.72(1)K 4 . This result provides the first quantitative estimate of the thermal confinement/deconfinement transitionof the U(1) quantum link model, as described qualitatively e.g. in [START_REF] Castro Neto | Ice: A strongly correlated proton system[END_REF].

Photon gap

One can extract an estimate of the excitation gap from the properties of the algorithm. Indeed, a plaquette membrane creates (or annihilates) two discontinuities in the imaginary-time propagation (i.e. two pairs of flippable plaquettes of opposite magnetization). The average Trotter length Mz pl of the membranes updated by the algorithm can be then related to the imaginary-time correlations of the flip operators F (0)F (τ), where F (τ) is acting on plaquette at imaginary time t = iτ. It is given by

F (τ)F (0) = 1 Z m e -βE m m| e iH(-iτ) F e -iH(-iτ) F |m = 1 Z m,n e -βE m -(E n -E m )τ |m| F |n| 2 (8.4)
where the |m are the eigenstates of the Hamiltonian, with energy E m . Be ∆ = E 1 -E 0 the gap above the ground state in the spectrum. Assuming that the plaquette flip operator F connects the ground state to the first excited state (i.e. 0| F |1 0), at low temperature, one has then

F (0)F (τ) ≈ β→∞ n e -(E n -E 0 )τ |0| F |n| 2 ≈ τ→∞ e -∆ τ |0| F |1| 2 (8.5)
The correlation length in the Trotter dimension at T = 0, ξ z , is then given by the inverse of the gap :

ξ z =∆ -1 (8.6)
The correlation length ξ z can be estimated via the characteristic distance between imaginary time discontinuities, namely the average length of plaquette membranes in imaginary time, βz pl . Setting ξ z ∼ βz pl , we obtain

∆ = ξ -1 z ∼ β→∞ (βz pl ) -1 .
(8.7)

Figure 8.3 represents T/(z pl ) as a function of K. This quantity appears to saturate to a finite value ( 1) when K 1. This is a clear proof of the gapped nature of the pVBS ground state. The lowest excitation of this frustrated lattice QED model for S = 1/2 is therefore a gapped photon, as it could be anticipated by extending Polyakov's argument on confinement in two dimensions [START_REF] Polyakov | Quark confinement and topology of gauge theories[END_REF] to the case of a frustrated compact lattice gauge theory.

Phase diagram of the transverse-field Ising model on the checkerboard lattice

We now turn to the original transverse-field Ising model, whose Hamiltonian is the one given in (6.14). Its phase diagram is represented in figure 8.4. Note the logarithmic scale on the y axis, emphasizing the fact that the physics we discuss here occurs at very low temperature (much lower that the classical single-spin flip gap ∆=4J).

At low temperature and low transverse field, the system is expected to be in a pVBS phase. This corresponds to the ground state of the quantum link model, that approximates correctly the square ice model at small Γ. Increasing the field at low temperature, the system is found to enter a magnetically ordered phase, namely the Néel phase. Increasing instead the temperature at low transverse field melts the pVBS phase and the Néel phase, and the system enters a quantum Coulomb phase. This phase is very similar to the classical Coulomb phase, the main difference being that in the quantum Coulomb phase the monopoles propagate coherently whereas in the classical one they have a thermally induced diffusive motion. Both phases smoothly connect to each other for temperatures T ∼ Γ, marking a smooth crossover. At large transverse-field, the system enters a quantum paramagnetic phase, where the polarization of the spins along the field induces strong fluctuations in the S z i component of the spins.

Plaquette valence-bond solid phase

At low field and low temperature, the system is in a pVBS state (as described in section 6.2.1). This state breaks the translational symmetry but not the spin inversion symmetry. In particular, the excitations are gapped and the defects (monopoles) are confined in this phase, accordingly to Polyakov's mechanism for two-dimensional compact QED. The transition line T pV BS separating the pVBS from the quantum Coulomb phase was drawn using the estimate of the critical temperature T c = 0.72(1)K 4 of the quantum link model

Magnetic structure factor and flippability structure factor

In order to track the Coulomb phases, as well a the magnetically ordered Néel state, we compute the magnetic structure factor of the system, as defined in (7.20). The results are given in figures 8.5 and 8.6. The pVBS phase has a staggered flippability, whereas the Néel phase is characterized by a uniform flippability. In order to get insight into the difference between these states, and in particular in order to determine whether there can be a coexistence of the two phases, we compute the flippability structure factor, defined as

S q = 1 N 2 , e iq•(r -r ) f f (8.8)
where r is the position of the bottom left corner of plaquette . The results are given in figures 8.5 and 8.8.

Néel phase

At intermediate transverse fields (0.25 Γ 0.55), the magnetic structure factor develops Néel Bragg peaks when the temperature is decreased, as can be seen in figures 8.5 and 8.6. The Néel phase boundaries were determined by computing the Néel order parameter, defined as with L the linear system size. In order to ascertain the presence of Néel order in the thermodynamic limit, we fit the system size dependance of m s to a cubic poylomial in 1/L, for various temperatures and trans- verse fields. At the phase boundary, the extrapolated order parameter thus obtained changes from zero to a finite value (see figure 8.9). Nevertheless, the extrapolated order parameter is very small (m s ≈ 0.02 for both Γ=0.360J and Γ=0.547J). This indicates that this phase must be richer than a simple magnetically

m s = 1 L 2 i (-1) i σ z i (8.9)
.9: Scaling of the Néel order parameter at T/J = 10 -2 ; solid lines are fits to cubic polynomials.

ordered Néel phase. Indeed, this can also be viewed in the flippability structure factor. The latter quantity indeed displays a sharp Bragg peak for q = (0, 0) (not represented in the figures), as expected in a uniformly flippable Néel phase, but this peak coexists with seemingly sharp Bragg peaks at (π, 0) and (0,π), corresponding to a staggered modulation of the flippability. This suggests that this Néel phase still possess pVBS order. Hence, the state at hand is far from being a classical state with weak quantum fluctuations. This is coherent with the fact that spin-wave theory fails to predict this ordered phase; Néel order is induced by a highly non-linear quantum order-by-disorder mechanism as indicated by degenerate perturbation theory.

The understanding of the Néel phase can be obtained from the 8-th order perturbation theory, as described in 6.2.2. Restricting the perturbation Hamiltonian to (6.24) we obtain the simplest model describing the competition between the pVBS phase and a maximally flippable state [START_REF] Moessner | Ising models of quantum frustration[END_REF], namely the Néel state. This model has been studied in [START_REF] Shannon | Cyclic exchange, isolated states, and spinon deconfinement in an xxz heisenberg model on the checkerboard lattice[END_REF] and [START_REF] Banerjee | The (2 + 1) -dU(1) Quantum Link Model Masquerading as Deconfined Criticality[END_REF] using exact diagonalization, and a transition from Néel to pVBS is predicted to occur for a critical ratio K 8 /K 4 ≈ 0.37. The ratio between these two coefficients can be controlled in the TFIM via the transverse field, and the above ratio between K 8 and K 4 is achieved at Γ/J ≈ 0.64, a value which lies reasonably close to the field range in which the Néel order is seen to appear in figure 8.4, considering that the effective Hamiltonian of the TFIM contains many other terms at 6-th and 8-th order beside those retained in (6.24).

Nonetheless, the fact that the Néel and pVBS states coexist in this phase would imply that the quantum phase transition occuring between the pVBS and the Néel phase at T = 0 (not accessible to our simulations), might be a continuous, second order transition, unlike what happens in the model (6.24) for which a first order transition is predicted [START_REF] Banerjee | The (2 + 1) -dU(1) Quantum Link Model Masquerading as Deconfined Criticality[END_REF][START_REF] Shannon | Quantum ice: A quantum monte carlo study[END_REF][START_REF] Chandrasekharan | Quantum link models: A discrete approach to gauge theories[END_REF].

Quantum paramagnet

At larger transverse fields, the system approaches nearly full polarization in the x direction as the field is increased. This phase phase could be a priori thought of as completely disordered in the S z i : a trivial quantum paramagnet. However, as can be seen in figure 8.6, when increasing the transverse field beyond the Néel phase, the correlations of the Coulomb phase seem to be recovered. This property has yet to be explained, as it goes beyond the range of applicability of the theories at hand (gauge mean-field theory and degenerate perturbation theory), and this phase might be richer than just a trivial quantum paramagnet.

Coulomb phase

We end our discussion of the phase diagram with the Coulomb phase. At "high" temperature, i.e. when the temperature is comparable with the spin-flip gap ∆, this corresponds to the thermal Coulomb phase described in 4.2.1. The monopoles/spinons are created and diffuse in the system via thermal fluctuations.

For T pV BS ≤ T Γ ∆, the spin-flip gap prevents the apparition of thermal defects. Indeed, for βJ ∼ 10 -2 , the density of such defects is of order e -100 , which is such an astronomically small number that any system of realistic finite size will be free of thermally induced defects in practice. Nonetheless the transverse field induces ground-state defects (vacuum spinon-pair fluctuations) which are identifiable with bound spinon pairs as long as the field is highly off-resonant compared to the spinon gap, namely Γ 4J. On the other hand, the monopole/spinon excitations living at energies Γ ∼ 4J differ substantially from their classical spin-ice counterparts as they can propagate coherently through the system with hopping Γ. While they may form a bound state due to gauge-field confinement, their binding energy is expected to be of the order of T pV BS . Hence for T pV BS T Γ the spinons are effectively deconfined thanks to thermal fluctuations, but thermal fluctuations in turn are too weak to affect the ballistic dynamics of the spinons.

This corresponds therefore to a quantum version of the Coulomb phase, or quantum Coulomb phase, attained at finite but exceedingly low temperature.

At T = 5×10 -3 J, the Coulomb phase survives up to Γ 0.25. This can be seen in the magnetic structure factor that retains all of its features, including the pinch points, indicating that the nature of correlations is the same as in the classical Coulomb phase up to astronomically large length scales (the average separation between thermally induced defects). In the deconfined quantum Coulomb phase, the excitations can be interpreted in terms of spinons of the gMFT described, in section 6.4.

When reaching the Néel/pVBS phase, the pinch points acquire a finite width, and their height is strongly reduced (figure 8.10). This effect directly signal the confinement transition. Both the Néel state and the pVBS state have a very large overlap with the w = 0 of ice-rule state, and they are compatible with the strong suppression of the pinch points height in this phase, as observed in figure 8.10. This suppression corresponds to a confinement transition of the spinons. This effect has already been seen in quantum pyrochlore spin ice [START_REF] Benton | Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice[END_REF], where the pinch points height is reduced by the quantum fluctuations, but restored by the thermal fluctuations. The situation is similar here, with the fundamental difference that pinch points are recovered via a sharp deconfinement transition. It is exactly the case here, as thermal and quantum fluctuations are similarly competing.

In figure 8.11, the transverse magnetization obtained from the Monte Carlo calculations is compared to the predictions of the gauge mean field theory described in chapter 6. This shows that gMFT is quantitatively The gMFT holds at low field despite the fact that it completely ignores the dynamics of the gauge field (6.33), and hence completely misses its confining effect on the spinons. The confining effects are then clearly not relevant in this regime. Using the quantum rotor picture developed in 6.4.3, we can regard the quantum Coulomb phase as a bosonic Mott insulator phase of the spinon field, and its elementary excitations (pairs of spinons) as particle hole pairs in a bosonic Mott insulator, with a gap corresponding exactly to the spin-flip gap. The ground state of a Mott insulator experiences quantum nucleation of bound particle-hole pairs, which corresponds to the vacuum spinon-pair fluctuations in quantum square ice. The result of such fluctuations is a finite kinetic energy of the spinons, and a corresponding finite transverse magnetization. Furthermore, the particle-hole excitations of a bosonic Mott insulator are gapped, deconfined and form a continuum [START_REF] Golovach | Dynamic response of one-dimensional bosons in a trap[END_REF]. Therefore the Mott insulator analogy offered by gMFT further corroborates the picture in which the excitation spectrum for the matter sector of the quantum Coulomb phase consists of a continuum of deconfined spinons.

Conclusion

We have found that quantum square ice, in its TFIM realization, has a rich phase diagram. Using quantum Monte Carlo based on the membrane algorithm, we were able to efficiently simulate the model throughout the temperature-field plane. We confirm the existence of confined phases at low temperature, as expected for the compact QED in two dimensions. We find evidences of two different confined phases with different symmetries, namely the pVBS phase and the Néel phase.

Moreover, we find a sizable region where the spinon excitations are thermally deconfined, and the spinspin correlations are those of a Coulomb phase. This quantum phase realizes a (thermally) deconfined U(1) spin liquid.

Conclusions and perspectives

The recurring theme of this thesis is the Coulomb phase. Indeed, we investigated its robustness in different realistic implementations of square ice. The first one is the dipolar square ice, whose study was inspired by experiments in the recent and rapidly developing topic of artificial spin ice [START_REF] Nisoli | Artificial spin ice: Controlling geometry[END_REF]. These experiments were designed to observe the physics of ice in two-dimensional systems with controllable geometry. However, the geometries that were realized experimentally up to now have prevented any direct observation of the Coulomb phase, in particular because dipolar interactions lead to a strong lift of the degeneracy of the ground state with respect to ideal square ice. We have here characterized the thermodynamical properties of the dipolar square ice model, when the height between the two dipole sublattices is varied. We reconstructed the phase diagram and we identified an optimal parameter range for the observation of the Coulomb phase.

Our phase diagram may serve as a guide for the experimental realization of artificial spin in the context of magnetic nanoarrays as well as for other experimental implementations. In particular, our model is well suited for the description the properties of spin ice realized with trapped nanoparticles [START_REF] Libál | Realizing colloidal artificial ice on arrays of optical traps[END_REF] as well as with trapped ions [START_REF] Schmied | Quantum phases of trapped ions in an optical lattice[END_REF][START_REF] Schneider | Experimental quantum simulations of many-body physics with trapped ions[END_REF]. For the specific case of a trapped-ion realization of square ice, the form of the interactions can be tuned, and it would be important for experiments to determine wether there is an optimal choice of interactions that facilitates the experimental observation of the Coulomb phase.

Besides, we find that the dipolar square ice orders into a four-fold degenerate ground state when its geometry is set to favor the alignment of parallel dipoles on each vertex. Interestingly, it turns out that the universality class of this transition depends on the range of the interactions. This phenomenon might be linked to the difference in the bandwidth of the normal mode spectrum of the model when truncating the interactions at a finite distance. In order to shed lights on this peculiar phenomenon, one might want to tune the range of the interactions continuously from short to long by decreasing the decay exponent of algebraically decaying interactions. While this might appear as a purely theoretical construction, it is actually accessible to the experiments on trapped ions as shown in [START_REF] Britton | Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of spins[END_REF]. Moreover, it would be interesting to investigate the phase diagram of dipolar spin ice in the case of dipoles perpendicular to the plane [START_REF] Mengotti | Dipolar energy states in clusters of perpendicular magnetic nanoislands[END_REF].

Drawing on this experience with classical dipolar spin ice, we studied the effects of quantum fluctuations in the ideal square ice realized by a microscopic Hamiltonian, namely the transverse-field Ising model on the checkerboard lattice. Transverse-field Ising models are relevant for the description of trapped-ion experiments, and therefore an experimental realization of the microscopic model Hamiltonian can be envisaged in that context. We constructed the phase diagram of this model, looking, among other phases, for a quantum Coulomb phase. To this end, we developed a novel Monte Carlo scheme that promotes the loop algorithm, commonly used in the study of classical spin ice, to a membrane algorithm for the (2 + 1)-dimensional effective classical model reproducing the same partition function as the quantum model. This algorithm efficiently samples the partition function of both the transverse-field Ising model and of effective Hamiltonians obtained via degenerate perturbation theory for weak fields. Among these Hamiltonians stands the quantum link model, realizing a lattice model of quantum electrodynamics for S = 1/2 variables.

In the phase diagram, we find a (thermally) deconfined phase that possesses the same spin correlations as the Coulomb phase up to extremely large length scale. This phase is expected to support fractionalized spinon excitations. Our next step will hence be the quantitative characterization of the excitations of such a phase. This can be done by computing the spin correlations in imaginary time, from which the gap of the excitations (if finite) can be directly extracted. The calculation of the complete dynamical structure factor (i.e. the spectrum of the excitations) can be obtained from these very correlations via a Wick rotation. On the other hand, the analytical continuation implied by the Wick rotation is an ill-posed numerical problem, that demands a high precision in the calculations of the correlations [START_REF] Jarrell | Bayesian inference and the analytic continuation of imaginary-time quantum monte carlo data[END_REF]. The exploration of different approaches to analytical continuations -some of them proposed in the recent literature [START_REF] Prokof | Spectral analysis by the method of consistent constraints[END_REF] -is work in progress.

Our results show that the transverse-field Ising model, despite its simplicity, realizes fundamental confinement/deconfinement transitions associated with an underlying lattice gauge theory. In this sense it represents a strong candidate for the quantum simulation of the phase diagram of lattice gauge theories possibly more accessible than other proposals formulated in the context of cold atoms [START_REF] Banerjee | Atomic quantum simulation of dynamical gauge fields coupled to fermionic matter: From string breaking to evolution after a quench[END_REF][START_REF] Zohar | Simulating compact quantum electrodynamics with ultracold atoms: Probing confinement and nonperturbative effects[END_REF][START_REF] Tagliacozzo | Optical abelian lattice gauge theories[END_REF].

Furthermore, the phase diagram of the TFIM turns out to be very rich, and many salient features still need to be clarified. Among other things, the nature of the quantum phase transitions between the groundstate ordered phases has to be clarified, as well as the nature of the deconfinement transition induced by the temperature. Moreover, the paramagnetic phase obtained for the largest fields we explored is found to exhibit spin correlations surprisingly similar to those of the Coulomb phase. In particular, the possible appearance of a Coulomb phase at high fields does not allow a standard image in terms of a gauge theory. Therefore it would be interesting to test whether this putative phase is at all subject to the same constraints obeyed by U(1) lattice gauge theories in two dimensions (namely the presence of a low-temperature confinement).

Nonetheless, the membrane algorithm we developed could be adapted to a wide range of similar models. It is particularly suited to models that admit a loop gas representation. For example, it could be applied to the Anderson model of Ising spins on a pyrochlore lattice and check its low-energy properties. As a threedimensional model, it is not limited by Polyakov's mechanism, and it is known to to be deconfined [START_REF] Polyakov | Gauge fields and strings[END_REF]. Efficient Monte Carlo simulation would allow to compare the non-perturbative results to the rich phase diagram predictied by field theories. [START_REF] Hermele | Pyrochlore photons: The u(1) spin liquid in a s = 1 2 three-dimensional frustrated magnet[END_REF][START_REF] Shannon | Cyclic exchange, isolated states, and spinon deconfinement in an xxz heisenberg model on the checkerboard lattice[END_REF] The quadratic Hamitonian can then be written in the compact form (1) k ,...,l (n) k ) = L k . The transformation matrix T k is then obtained as [START_REF] Wessel | Quantum fluctuations and excitations in antiferromagnetic quasicrystals[END_REF][START_REF] Mucciolo | Excitations and quantum fluctuations in sitediluted two-dimensional antiferromagnets[END_REF][START_REF] Blaizot | Quantum Theory of Finite Systems[END_REF]]

H 2 = N n p       ε cl,p - hp 2       + 1 2 k A † k M k A k (B.
U k such that U † k Z † k ΣZ k U k = diag(l
T k = Z k U k |L k | -1 /2 .
(B.9)

In particular, the eigenmodes ω As mentioned above, in the particular case of the Néel and collinear states of the checkerboard lattice Ising model studied in this work, the classical energies ε cl,p are all equal to ε 0 . The mean energy of the system then reads

E = N(ε cl -ε 0 /2) + 1 2 k,p ω k,p = Nε cl + ε 0 4 k,p          1 + 4 λ k,p ε 0 -1          (C.1)
As ∆ k is proportional to Γ 2 , so are its eigenvalues. We can then expand the mean energy per spin in powers of Γ 2 . We will introduce rescaled eigenvalues λk,p defined as λ k,p = ε 0 Γ 2 λk,p .

ε = ε cl + ε 0 N k ∞ m=1 α m λm k,p Γ 2m = ε cl + ε 0 N k,m α m Tr ∆ k ε 0 m (C.2)
where the α m are defined by

√ 1 + 4x -1 = 4 ∞ m=1
α m x m . In the case of a Néel or a collinear case, the trace of ∆ k averages to zero in the Brillouin zone. Thus the first non-zero correction to the classical energy in Eq. (C.2) is of fourth order in Γ.

In both cases, if we introduce the ratio J = J 2 /J 1 , we get

ε = ε cl + p>1 c p (J) J 1 S (ν 2 -J) 2 p Γ 2p (C.3)
where the c p coefficients only depend on J, and not on the considered state. It is then obvious that the expansion becomes independent of the classical state if J = 1. Thus the classical degeneracy of the ice model is not lifted by harmonic fluctuations.

Similarly, if all the ω k,p are real (which is the case whenever spin-wave theory holds), we have where the β m are defined by

√ 1 + 4x + 1 √ 1 + 4x -2 = 4 ∞ m=1 β m x m .
Appendix D

Calculation of some of the coefficients in the perturbation expansion

We derive here some of the numerical coefficients of the perturbation theory. As the spin operators on different sites commute with each other, this amounts to counting the number of different sequences of the single spin flips that constitute the total move, weighted by the intermediate energies.

D.1 Simple loop moves

The simplest terms that arise from the perturbation theory at order n are the flips of closed loops of length n. The numerical prefactor of these terms is determined here. As the action of these terms vanishes if flipping the loop does not lead to a new ice-rules state, we will only consider configurations in which the considered loop connects two ice-rules states. The different spin flip sequences can be classified by the sequence of the energies of the intermediate states. A S x i operators creates or annihilates a pair of defects, or simply propagates one of them. Hence the change of energy after a single spin flip is either ±∆ or 0.

The contribution of a given energy sequence {ε i } = {p i }∆, i = 0,...,n for a loop of length n (with p i integers) is then α {p i } (P i ∆) -n with P i = i p i and α {p i } the number of different spin flip sequences that yield the same energy sequence {ε i }. All one has to do is then to determine α {p i } .

D.1.1 4th order

For a plaquette flip, the energy sequence is either (0, ∆, ∆, ∆, 0) or (0, ∆, 2∆, ∆, 0) (see figure D.1). (0, ∆, ∆, ∆, 0)

• • • • • • • • ∆ ∆ ∆ ∆ ∆ 2∆ 0 0 0 0 • • • • • • • • • • • • S x i 1 S x i 2 S x i 3 S x i 4 S x i 1 S x i 2 S x i 3 S x i 4
For this move, there are 4 choices for the starting point, and two choices for the second and the third step (the next spin flip has to be on one of the neighbors of the current flip). This gives a contribution 4 × 2 × 2/∆ 3 = 16 ∆ -3 (0, ∆, 2∆, ∆, 0) For this move, there are 4 choices for the starting point. The second spin flip has to be on the opposite site of the plaquette, and the third one can be on either of the two remaining sites. This gives a contribution 4 × 1 × 2/(2∆ 3 ) = 4 ∆ -3

Total prefactor

Summing the all the contributions, one gets the total prefactor a 4,4 ∆ -3 = 20 ∆ -3 (D.1)

D.1.2 6th order

The same reasonning holds for any loop size. For the 6-th order loops, the calculation remains humanly tractable. The moves corresponding to the energy sequences (0, ∆, ∆, ∆, ∆, ∆, 0) and (0, ∆, 2∆, 3∆, 2∆, ∆, 0) are represented in figure D.2. (0, ∆, ∆, ∆, ∆, ∆, 0)

• • • • • • • • • • • • • • • • • • • • • S x i 1 S x i 2 S x i 3 S x i 4 S x i 5 S x i 6 ∆ ∆ ∆ ∆ ∆ 0 0 • • • • • • • • • • • • • • • • • • • • • S x i 1 S x i 2 S x
There are 6 choices for the starting point of the loop, and then 2 choices at each of the 4 following steps, so that this term gives a contribution 6 × 2 4 /∆ 5 = 96 ∆ -5 . (0, ∆, 2∆, 3∆, 2∆, ∆, 0) There are 6 choices for the starting point of the loop, and then 2 choices at the following step between its two next-to-nearest neighbors. The third one is deterministic. The next step is the flip of one of the three in-between spins, and then one of the two remaining spins. The total contribution of this is then 6 × 2 × 3 × 2/(12∆ 5 ) = 6 ∆ -5 . (0, ∆, 2∆, ∆, ∆, ∆, 0), (0, ∆, ∆, 2∆, ∆, ∆, 0), (0, ∆, ∆, ∆, 2∆, ∆, 0) For all these moves, there are 6 choices for the starting point of the loop, and then 2 choices at each step but the one that lower the energy from 2∆ to ∆ and the last one which are deterministic. The contribution of these terms is then 3 × 6 × 2 3 /∆ 5 = 96 ∆ -5 .

Other sequences

We construct similarly the contribution of the remaining coefficients (0, ∆, 2∆, 2∆, ∆, ∆, 0), (0, ∆, ∆, 2∆, 2∆, ∆, 0) → 48 ∆ -5 (0, ∆, 2∆, ∆, 2∆, ∆, 0) → 6 ∆ -5 (0, ∆, 2∆, 2∆, 2∆, 2∆, 0) → 48 ∆ -5

Total factor

Summing up all these contributions, we get the total prefactor a 6,6 ∆ -5 = 252 ∆ -5 (D.2) D.1.3 8th order and higher

For the 8th order loops, the calculation involve several dozens of terms. We resort then to computer assisted calculations. We obtain where γ ≡ Γ/∆ and the operator F 6l (resp. F 8l ) is the flip operator of a loop of length 6 (resp. 8). that vanishes if it acts onto a states in which this loop does not connect to another ice-rules state.

D.2 Virtual flip of a loop

At order 8 appears the first non-trivial diagonal term. It flips twice each of the four spins of a plaquette. This gives a quantum correction to every classical states that depends on the configuration of its plaquette. Thanks to the symmetries of the lattice, there are only four different energy corrections, depending only on the total magnetization m p of the plaquette (±4S , ±2S or 0) and its flippability (see figure D.3). The method to derive these corrections is a same as in the case of the off-diagonal terms.

Using the notations of figure D.

D.3 Combined loop and virtual flip moves moves

The last family of terms are combinations of virtual flips and loop moves. A few examples of such moves are represented in figure D.4. This terms give a non extensive part (that cancels out with the one obtained from the trivial terms) and a relevant extensive part. The first of these terms appear at order 6.

• • • • • • • • • • • • • • • Figure D.4:
Possible non trivial sixth-order plaquette moves. The flip of the plaquette is accompanied by the flip of a spin that can be non neighboring the plaquette (blue), neighbor of the plaquette (green) or part of the plaquette (red). In the latter case, the spin is flipped three times by the move. All these moves give O(γ 6 ) corrections to the fourth-order plaquette term.

D.3.1 6th order corrections

The sixth order terms that arise are corrections to the fourth-order plaquette term. They can be of three different kind (see figure D.4). The fifth spin, that is flipped twice, can be either part of the plaquette, a neighbor of the plaquette or any other spin. The resulting correction to the Hamiltonian is H (6) 4 = -

Flippability

We define the flippability f of a plaquette state | as

f = | 2F 2 -1 | (E.4)
Where F is the operator that flips the spin on plaquette . The flippability is 1 if the plaquette is in one of the two Néel states, and -1 otherwise. We then calculate the flippability structure factor. This allows to distinguish between the Néel and the pVBS phase. Indeed, in the former the flippability is uniform ( f ≡ 1) whereas it is staggered in the latter ( f = (-1) ).

For the 4 × 4 system, the ground state has indeed staggered flippability and almost no mean flippability :

f π,π = 0| 2F 2 -1 (-1) |0 = 4/5 f 0,0 0| 2F 2 -1 |0 = 1/5 (E.5)
They are not exactly 1 and 0 because of the contributions of the two Néel-states components of the ground state.

Glace bidimensionnelle classique et quantique : phases de Coulomb et phases ordonnées

La frustration -c'est-à-dire la présence d'interactions de nature compétitive -donne lieu à des effets de grande complexité en physique. La glace -aussi bien la phase bien connue de l'eau, que ses analogues magnétiques, dites glaces de spin -en offre un exemple remarquable. Pour des interactions à courte portée et des degrés de liberté classiques, son état fondamental est infiniment dégénéré, et comporte en outre des corrélations à longue portée induites par une contrainte locale, caractéristiques de la phase dite de Coulomb. Ses excitations élémentaires correspondent au retournement d'un dipôle qui se fractionnalise en deux monopôles. Dans cette thèse nous nous intéressons à la stabilité de cette phase de Coulomb dans la glace bidimensionnelle -réalisée aussi bien comme glace de protons dans des composés organiques, que comme glace de spin dans des systèmes nanomagnétiques.

Dans le cas classique, les interactions dipolaires -présentes dans les systèmes expérimentaux -déstabilisent la phase de Coulomb dans son état fondamental. Cependant, une déformation de la simple géométrie planaire permet de récupérer cette phase dans un régime où différents états ordonnés entrent en compétition. Dans le cas quantique, les fluctuations dues à un champ magnétique transverse induisent une brisure de symétrie dans l'état fondamental qui, à basse température, cède la place à une phase de Coulomb quantique, réalisant un liquide de spin quantique avec excitations fractionnalisées. Nos résultats sont obtenus à l'aide de méthodes analytiques (analyse harmonique classique et quantique et théorie de perturbations) aussi bien que numériques (Monte Carlo) fondées sur des algorithmes originaux.

Mots clés : Magnétisme frustré, glaces de spin, phase de Coulomb, Monte Carlo, transitions de phase, liquide de spin

Classical and quantum two dimensional ice : Coulomb and ordered phases

Frustration -namely the presence of competing interactions -gives rise to highly complex effects in physics. Ice -be it the well known water ice or its magnetic equivalent, the so-called spin-ice -offers a remarkable example in this context. For short-range interactions and classical degrees of freedom, its ground state is infinitely degenerate, and exhibits long-range correlations induced by a local constraint, characterizing the so-called Coulomb phase. Its elementary excitations correspond to the flip of a dipole fractionalizing into two monopoles. In this thesis we are interested in the stability of this Coulomb phase in the two-dimensional ice -realized in the form of a proton ice in organic compounds as well as of spin ice in nanomagnetic systems.

In the classical case, dipolar interactions -present in the experimental systems -destabilize the Coulomb phase in the ground state. However, a slight deformation of the simple planar geometry allows to recover this phase in a regime where different ordered states compete with each other. In the quantum case, fluctuations due to a transverse magnetic field induce a symmetry breaking in the ground state, that melts at low temperature into a quantum Coulomb phase, realizing a quantum spin liquid with fractionalized excitations. Our results were obtained with analytical (classical and quantum normal-mode analysis and perturbation theory) as well as numerical techniques (Monte Carlo) based on original algorithms.

Key words : Frustrated magnetism, spin ice, Coulomb phase, Monte Carlo, phase transitions, spin liquid

10 CHAPTER

 10

Figure 1 . 2 :

 12 Figure 1.2: (left) Pyrochlore structure of ice. The white and black spheres are oxygen and hydrogen atoms. (middle) Anderson model. (right) Equivalent structure of spin ice. The colors (blue/red) are the corresponding Ising variables (±1).

Figure 1 . 3 :

 13 Figure 1.3: Specific heat (top) and entropy(bottom) of Dy 2 Ti 2 O 7 .Figure taken from[START_REF] Ramirez | Zero-point entropy in 'spin ice[END_REF].

Figure 1 . 5 :

 15 Figure 1.5: Determination of the winding number.Each component of the winding number is obtained as the total magnetization of the spin cutting a vertical (w x ) or horizontal (w y ) line. In this example we have w x = -3 and w y = 1. It is easily verified that it is indeed the same number for each line/column.

Figure 1 . 7 :

 17 Figure 1.7: Structure factor of the square ice model, from equation (1.11).

Figure 1 . 8 :

 18 Figure 1.8: Structure factor of Ho 2 Ti 2 O 7 obtained from Monte Carlo simulations (left) and neutron scattering (right). Pinch points structures appear at the (111) points. Figures taken from [43].

Figure 1 . 9 :

 19 Figure 1.9: (left) AFM image of a permalloy array. (right) MFM image taken from the same array. Figures taken from [151].

Figure 1 . 10 :

 110 Figure 1.10: Short range dimer coverings superposed in the resonating valence-bond state on the triangular lattice.

Figure 1 .

 1 Figure 1.12: Dynamic Structure factor of quantum spin ice as calculated by Benton et al., showing the existence of a gapless photon mode. Figure from [16].

Figure 1 . 13 :

 113 Figure 1.13: Static structure factor of classical spin ice at T = 0(left), of quantum spin ice at T = 0(middle), and of quantum spin ice at finite temperature (right). Figures from[START_REF] Benton | Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice[END_REF]

Figure 2 . 5 :

 25 Figure 2.5: Normal modes of the square lattice antiferromagnetic Ising model with coupling J > 0.There are minima at the equivalent (±π, ±π) points with energy -2J, corresponding to the Néel states which is the ground state of the model. The maximum at (0, 0) corresponds to the ferromagnetic state which is the highest energy state.

Figure 2 . 7 :

 27 Figure 2.7: Lower band of the normal modes spectrum of the model. The value of h is increased from left to right (0.2a, h c ≈ 0.805a, h LR c ≈ 0.85a and 1.20a). The range of the interaction is increased from top to bottom (nearest neighbor, 1, 2 and 64 lattice spacings). For h = h c , the lower band is perfectly flat when retaining nearest-neighbour interactions only. For all values of h, the width of the band reduces in the long range limit. In particular, away from h c the global form of the spectrum is recovered in the long range limit. The long range and nearest neighbor spectra are particularly close around the minima. For h = h LR c , the minima at the vortex state -(±π, ±π) -and the Manhattan state -(0, ±π) and (±π, 0) -have the same depth, signaling the degeneracy between them.

Figure 2 . 8 :

 28 Figure2.8: Evolution of the width δω of the lower band of excitations in the short range (dotted black line) and long range (solid red line) interaction limits with h. This curve delimits the expected phase diagram of the model. We expect here three phases : The disordered paramagnetic phase, separated from the ordered phases by δω(h), and the two ordered phases (the vortex and the Manhattan states) separated by a vertical line at the minimum of δω.

Figure 3 . 2 :

 32 Figure 3.2: In all panels, data for ideal ice systems of size L = 32 (green), L = 64 (blue) and L = 128 (purple) are shown. (top panels) Acceptance rate of the long loop (left) and short loop (right) updates as a function of the temperature. (bottom panels) Average length of the long (left) and short (right) loops. The inset in the bottom left panel is the ratio by L 5/3 of the long-loop average length, showing that the T = 0 scaling is indeed recovered at low temperature.

Figure 3 . 3 :

 33 Figure 3.3: Efficiency of the loop algorithms in the dipolar square ice, for various geometries introduced in section 2.2 (blue : h = 0.2 a, red : h = h LR c , green : h = 1.2 a) and system sizes (squares : L = 16, circles : L = 24, triangles : L = 32). (top panels) Acceptance rate of the short (left) and long (right) loop updates as a function of temperature. (lower panels) Average length of the flipped short (left) and long (right) loops.

Figure 3 . 4 :

 34 Figure 3.4: Mapping between a vertex configuration and an Ising-spin configuration. Starting from the initial lattice, rotated to a π/4 angle, one first takes the y projection of the oriented links and then flip one every two horizontal lines.

Figure 3 . 6 :

 36 Figure 3.6: (left) Structure factor of the pseudospins for the six vertex model in units of the inverse distance between neighboring vertices. The black square is the magnetic Brillouin zone. (right) Structure factor for the Ising model on the checkerboard lattice. The first Brillouin zone is the black square.

Figure 4 . 1 :

 41 Figure 4.1: Phase diagram of the dipolar square ice, obtained for a trap half-width λ = 0.2 a.the two panels only differ in the range of temperature represented. The red dots are the Monte Carlo data obtained from the finite-size scaling analysis of the order parameter. The black dashed line is the normal mode band width. The green dashed line is the energy difference between the Néel and the (π/2, π/2) collinear state. The dotted lines visible in the right panel are the spin-flip gap (blue) above the ground state. These two last lines set a characteristic temperature scale for the crossover from the Coulomb phase to the paramagnetic phase.

Figures 4. 2

 2

  , 4.4 and 4.6 present the evolution with temperature of the structure factor of the dipolar square ice, at different points in the parameter space (h = 0.2 a, h = h LR c and h = 1.2 a). In all cases, pinch point like features are visible at sufficiently high temperature. h ≈ h LR c

Figure 4 . 2 :

 42 Figure 4.2: Magnetic structure factor of a dipolar square ice of size L = 64 at λ = 0.2a and h = h LR c . for the ideal model is represented in figure 4.3. The pinch points survive at temperatures much larger than the spin-flip gap. By definition, there is no ordering at low temperature in this case.

Figure 4 . 3 :

 43 Figure 4.3: Magnetic structure factor of an ideal square ice (six-vertex model with J 1 = J 2 = J) of size L = 64.

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Magnetic structure factor of a dipolar square ice of size L = 32 at λ = 0.2a and h = 1.2 a > h LR c .

Figure 4 . 6 :

 46 Figure 4.6: Magnetic structure factor of a dipolar square ice of size L = 32 at λ = 0.2 a and h = 0.2 a < h LR c .

Figure 4 . 7 : 4 T /J 2 Figure 4 . 8 :

 474248 Figure 4.7: Half-width at half maximum of the pinch point as a function of temperature for different height offset (left : h = 0.2 a, middle : h = h LR c , right : h = 1.2 a), for a size L = 64 (middle panel) and L = 32 (other panels). The dashed lines represent the transition temperatures. The black dotted curve is the same quantity for the ideal square ice with coupling J as a function of T/J with a linear size L = 64.

Figure 4 . 9 :

 49 Figure 4.9: Collapse plot of the Binder cumulant (top left), the heat capacity (top right), the Néel order parameter (bottom left) and the Néel susceptibility (bottom right) for the dipolar square ice at λ = 0.2 a and h = 0.2 a for system sizes ranging from L = 12 to L = 52. The value of the critical temperature used is T c = 3.25J 2 .

Figure 4 .

 4 Figure 4.10:The clock models. The spins rotors point in any of the q angles defined by the model. For q = 2, one recovers the Ising model, whereas the XY model is obtained in the q →∞limit. (right) Spectrum of the single-bond energy for q = 4.

Figure 4 .

 4 Figure 4.11:The Potts model. The variables can take any value between 1 and q. These can be viewed as q different colorings. For q = 4, there is a single minimal energy configuration with three excitations for each pair of interacting sites.

Figure 4 . 12 :

 412 Figure 4.12: Scaling of the position of the peaks of the heat capacity (squares) and the susceptibility (red), as a function of L -1/ν (log L) ν, with the Potts exponents (left) and the clock (Ising) exponents (right).

Figure 4 . 13 :

 413 Figure 4.13: (left) Binder cumulant of the Manhattan order parameter for system sizes ranging from L = 24 to L = 48. The different curves intersect around T c ≈ 0.91J 2 . The onset panel is a closer view around T c . Rescaled Binder cumulant using Ising (middle) and Potts critical exponents (right).

Figure 4 . 14 :

 414 Figure 4.14: Rescaled specific heat capacity (top), order parameter (middle) and susceptibility (bottom) Potts (left) and Ising (right) critical exponents. The data are obtained for system sizes ranging from L = 24 to L = 48. When present, logarithmic corrections have been taken into account.

Figure 4 . 16 :

 416 Figure 4.16: Scaling of the position of the peaks of the heat capacity (squares) and the susceptibility (red) of the J 1 -J 2 -J 3 model, as a function of L -1/ν (log L) ν, with the Potts exponents (left) and the clock (Ising) exponents (right). The extrapolated critical temperatures are the ones used in the finite size scaling of the other quantities.

Figure 4 .

 4 Figure 4.17: (above) Binder cumulant of the Manhattan order parameter for system sizes ranging from L = 40 to L = 128 (left). The different curves clearly intersect around T c ≈ 1.545 ± 0.002. The inset is just a closer view around T c .( below) Rescaled Binder cumulant using Ising (middle) and Potts (right) critical exponents.

8 (Figure 4 . 18 :

 8418 Figure 4.18: Rescaled specific heat capacity (top), order parameter (middle) and susceptibility (bottom) Potts (left) and Ising (right) critical exponents, in arbitrary units. The data are obtained for system sizes ranging from L = 24 to L = 128.

10 )

 10 where b † l,r is the creation operator of a particle in the left (right) pocket of the trap. Imposing that each trap contains exactly one ion (b † l b l + b † r b r = 1), the creation and annihilation operators are converted into Pauli matrices

Figure 4 . 19 :

 419 Figure 4.19: Quantum square ice represented as a Transverse field Ising model on the checkerboard lattice. The spins lie on the sites of the square lattice (black dots), and their z components are coupled with couplings J 1 (solid lines) and J 2 (dashed lines).

Figure 5 . 1 :

 51 Figure 5.1: Classical canting angle ϑ induced by the transverse field Γ.

Figure 5 . 2 :

 52 Figure 5.2: Classical ground state phase diagram of the TFIM. The zero field phases are smoothly connected to the polarized states as the z components are reduced to 0 while the {σ i } configuration remains the same.

Figure 5 . 3 :

 53 Figure 5.3: (left) Structure factor along three straight lines joining high symmetry points of the Brillouin zone O, P and N for transverse field intensities varying from 0 to 2. These points are represented in the middle panel. The background colormap is the structure factor in the T = 0, Γ=0 case, as a guide to the eye, and the dotted line is the first Brillouin zone of the square lattice. (right) The same structure factor divided by cos 2 ϑ. This renormalization accounts for the reduction of the intensity up to values of Γ 1.75 J.

Figure 5 . 4 :

 54 Figure 5.4: Finite-size scaling of the Néel order parameter in the continuous spin (S →∞ ) limit for different values of Γ, at T/J = 5 × 10 -3 . The order parameter extrapolates to 0 for the entire range of transverse field values -solid lines are fits to cubic polynomials.

Figure 5 . 5 :

 55 Figure 5.5: Local basis used for the construction of the bosonic operators, depending on the orientation of the classical reference configuration.

Figure 5 . 6 :

 56 Figure 5.6: Lowest band of the magnon spectrum for Γ=J 1 S /2 and various values of J 2 /J 1 , around the Néel state. For the purpose of readability, the bands have been offset by the energy of their lower edge. The lower band edge (corresponding to the minimum excitation gap) is plotted in the lower panel.

Figure 5 . 7 :

 57 Figure 5.7: Lowest band of the magnon spectrum for Γ=J 1 S /2 and various values of J 2 /J 1 , around any collinear state. For the purpose of readability, the bands have been offset by the energy of their lower edge. The lower band edge (corresponding to the minimum excitation gap) is plotted in the lower panel.

Figure 5 . 8 :

 58 Figure 5.8: Upper panels: lowest band of the magnon spectrum around the polarized state for Γ > Γ c (left column) and for Γ=Γ c (right column). Lower panels: Minimum excitation gap: (left) as a function of J 2 /J 1 for ΓS /J 1 = 1.3 (red line) and for ΓS /J 1 = 1.43 (blue line); (right) as a function of ΓS /J 1 for J 2 /J 1 = 0.7 and 1.3 (the curves for both cases coincide).

Figure 5 . 9 :

 59 Figure 5.9: (left) Phase diagram of the S = 1/2 transverse-field Ising model on a checkerboard lattice from spin-wave theory. The hatched, and dense-hatched regions correspond to the region of inversion of the next-to-nearest neighbor correlations and to the energy-hierarchy inversion region, respectively. The thick dotted line corresponds to the saturation of the transverse magnetization (see text). Below the thin dotted lines one can find -in the classical limittwo states which are (local) energy minima and which have order in the z spin components: one state with Néel order, and the other with collinear order. Above this line one of the two states becomes the polarized state. Therefore collinear and Néel ordered states can only be compared energetically below the dotted line. (right) Zoom around the classical tri-critical point. In the black domain, the order parameters of all considered phases vanish.

Figure 5 . 10 :

 510 Figure 5.10: Order parameter vs field for different values of J 2 . It vanishes at the critical value of the field and a small gap opens for J 2 close to J 1 where none of the three states has a finite order parameter. The insets are closer views around the critical value of the field. The hatched and dense-hatched regions correspond those identified in the phase diagram (figure 5.9).

1 Figure 5 . 11 :

 1511 Figure 5.11: Average energy (upper panel) and order parameter (lower panel) associated with the Néel (red solid line) and collinear (blue solid line) reference states. The dashed lines correspond to the classical energies. Both panels are for Γ/J 1 = 0.9.

Figure 6 . 1 :

 61 Figure 6.1: The fourth-order moves are flips of a plaquette of spins. This flip can be executed in any sequence of single-spin flips, potentially creating a virtual state of energy 2∆, with two pairs of defects. This move swaps Néel states (|N or | N) and vanishes on any other plaquette state.

Figure 6 . 5 :

 65 Figure 6.5: Transverse magnetization predicted by the linear spin waves (dashed line) and the gauge mean field theory (solid line). Both approaches gives the same susceptibility at low field as in the classical limit.

Figure 7 . 1 :

 71 Figure 7.1: Trotter-Suzuki mapping for the TFIM on the checkerboard lattice. One start with a two dimensional classical Ising model with coupling J (left). The Trotter-Suzuki mapping leads to a three-dimensional Ising model with the same connectivity in each of the M layers as in the starting model (right). The intralayer coupling (solid lines) becomes J/M, and the interlayer coupling (dashed lines) is λ(β, ) as defined in the text.

  p in the intralayer energy of the p-th Trotter layer. The transverse field term can be transformed into an Ising interaction between layers. We introduce λ and ν so that S As λ<0, the intralayer coupling is ferromagnetic. The partition function for the model is then

Figure 7 . 3 :

 73 Figure 7.3: Acceptance rate of the membrane algorithm as a function of the transverse field at β/J = 5 10 -3 (left) and as a function of the temperature at Γ/J = 0.4, grown from short loops (top) or long loops (bottom). The different colors are different system sizes (green : L = 8, blue : L = 12, purple : L = 16 and red : L = 20).

3 Figure 7 . 4 :

 374 Figure 7.4: Average in-plane dimension (i.e. loop length) of the membranes as a function of the transverse field at β/J = 5 10 -3 (left) and as a function of the temperature at Γ/J = 0.4, grown from short loops (top) or long loops (bottom). The different colors are different system sizes (green : L = 8, blue : L = 12, purple : L = 16 and red : L = 20).

Figure 7 . 5 :

 75 Figure 7.5: Proportion of the Trotter dimension updated by the membrane algorithm as a function of the transverse field at β/J = 5 10 -3 (left) and as a function of the temperature at Γ/J = 0.4, grown from short loops (top) or long loops (bottom). The different colors are different system sizes (green : L = 8, blue : L = 12, purple : L = 16 and red : L = 20). Note that for the range of temperature considered at Γ/J = 0.4 the acceptance rate is very close to one.

  not flippable and | = | 0 otherwise (7.29) Therefore the partition function of a L × L quantum link model maps onto that of a classical model of L × L × 2M Ising spins with a partition function given by

1 Figure 7 . 6 :

 176 Figure 7.6: Trotter-Suzuki mapping for the quantum link model. Even plaquettes interact between layers 2p and 2p + 1 (blue), and odd plaquettes between layers 2p -1 and 2p (red). (right) The fundamental moves of the plaquette membrane algorithm are simultaneous flips of two adjacent non-interacting plaquettes.

Figure 7 . 7 :

 77 Figure 7.7: Acceptance rate of the long-loop membrane algorithm (left) and the plaquette membrane algorithm (right) as a function of K. The different colors are different system sizes (green : L = 8, blue : L = 16, purple : L = 24 and red : L = 32). The plaquette algorithm acceptance rate barely depends on K or on the system size L. The acceptance rate of the long-loop membrane algorithm drops at K ≈ K * (L) ∝ L -5/3 , which is plotted in the inset.

3 Figure 7 . 8 :

 378 Figure 7.8: Average in-plane length l ll of the longloop membranes as a function of K. The different colors are different system sizes (green : L = 8, blue : L = 16, purple : L = 24 and red : L =[START_REF] Chorin | Numerical solution of the navier-stokes equations[END_REF]. At low K, the classical scaling is recovered and l ll ∝L 5/3 , as can be seen in inset where l ll /L 5/3 is represented.

Figure 7 . 9 :

 79 Figure 7.9: Portion z pl of the Trotter direction updated by the plaquette algorithm as a function of K. The different colors are different system sizes (green : L = 8, blue : L = 16, purple : L = 24 and red : L = 32). As K is increased, the system decorrelates in the Trotter dimension, with corresponding reduction of the length of the plaquette membranes.

Figure 7 . 10 :

 710 Figure 7.10: Winding number fluctuations as a function of K. The different colors are different system sizes (green : L = 8, blue : L = 16, purple : L = 24 and red : L = 32). It depends only slightly on the system size. For large system size, the jumps from a small value to zero are due to the finiteness of the Monte Carlo simulations that become much less efficient out of the quantum Coulomb phase

Figure 8 . 2 :

 82 Figure 8.2: (left) Staggered flippability of the quantum link model as a function of the coupling for linear system sizes L = 16, 20, 24. There is clearly a transition to a state with finite staggered flippability. (left) Binder Cumulant of the staggered flippability for the same parameters. The inset is the plot of the crossing point K c (L) as a function of 1/L. The crossings extrapolate to K c = 1.38(2).

1 Figure 8 . 3 :

 183 T/z pl as a function of K = βK 4 , for system sizes ranging from L = 8 to L = 32.

Figure 8 . 4 :

 84 Figure 8.4: Phase diagram of quantum square ice. The boundaries of the pVBS phase (T pVBS ) and of the canted Néel phase (T N ) are transition lines. The blue dots are Monte Carlo data. The dashed line T =Γmarks a crossover from coherent to diffusive spinon/monopole dynamics at the energy scale set by the transverse field Γ.

Figure 8 . 5 :

 85 Figure 8.5: Magnetic structure factor of the S z component of the spins for the TFIM at Γ/J = 0.5 and L = 24.

Figure 8 . 6 :

 86 Figure 8.6: Magnetic structure factor of the S z component of the spins for the TFIM at T/J = 10 -3 and L = 24. At low field, the structure factor is very similar to that of the classical model. In particular, clear pinch points are visible. At intermediate field (0.25 Γ/J 0.55), Néel peaks appear.

Figure 8 . 7 :

 87 Figure 8.7: Flippability structure for the TFIM at Γ/J = 0.5 and L = 24. The dominant peak at (0,0) is not represented so that the other features remain clearly visible.

Figure 8 . 8 :

 88 Figure 8.8: Flippability structure factor of the TFIM at T/J = 10 -3 and L = 24. The dominant peak at (0,0) is not represented so that the other features remain clearly visible.

32 Figure 8 . 10 :

 32810 Figure 8.10: (left) Scans in the static structure factor at T/J = 5 × 10 -3 and L = 24, showing the evolution of the pinch-point width. (right) Pinch point height (i.e. S q (π, 0)) as a function of the transverse field at T/J = 5.10 -3 .

Figure 8 .

 8 Figure 8.11: Transverse magnetization of TFIM for T/J = 5 × 10 -3 and L = 16, compared with the gauge mean-field theory (gMFT) prediction. The vertical dashed line marks the transition from Mott insulator (MI) to superfluid (SF) in the corresponding quantum rotor model.

  126APPENDIX B. SPIN-WAVE THEORY FOR GENERAL ISING HAMILTONIANS IN A TRANSVERSE FIELDThe quadratic Hamiltonian then readsH 2 = E cl + lr l ) pp a † lp + a lp a † l p + a l pwhere hp = 2ε cl,p /S +Γsin ϑ p J(∆r) pp = J(∆r) pp sin ϑ p sin ϑ p (B.5)We then introduce the Fourier transform of the bosonic operators and of the interactiona k,p = 2 N l e ik.r l a l,p J(k) = l e -ik.∆r J(∆r) . (B.6)

7 ) where hp = 2 Scl 0 n 0 n ε cl ∆ k = 1 2 J

 722 ε cl,p +Γsin ϑ p A † k = (a † k,1 ,...,a † k,n , a -k,1 ,...,a -k,n ) M k = ∆ k ∆ k ∆ k ∆ k -ε (k) + J(k) † (B.8)This Hamiltonian can be diagonalized by a n-mode Bogolyubov transformation. This consists in finding the transformationA k = T k B k , with B k = (b † k,1 ,...,b † k,n , b -k,1 ,...,b -k,n ) T , such that A † k M k A k = p ω (p) k b † k,p b k,p and [b k,p , b † k,p ] = 1 and [b ( †) k,p , b ( †)k,p ] = 0. We introduce the matrix Σ, given byΣ= I n 0 n 0 n -I n ,the matrix Z k of the right eigenvectors of ΣM k , and the unitary matrix

k

  are the eigenvalues of ΣM k . If the matrices ∆ k and ε cl commute (which is the case for the Néel and collinear states of the checkerboard Ising model studied here, having ε cl = ε 0 I n ), the eigenmodes ω (p) k can be expressed in terms of the eigenvalues λ (p) k of ∆ k in the form ω (p) k = ε cl,p 2 1 + 4 λ (p) k ε cl,p . (B.10)

128APPENDIX 1 - 1 NS

 11 C. ENERGY AND MAGNETIZATION IN THE CASE OF N ÉEL AND COLLINEAR STATES m =

Figure D. 1 :

 1 Figure D.1: The two types of spin flips sequences that give rise to a plaquette flip. When acting on a Néel plaquette, the total move connects two ice-rules states. The energy of each intermediate state is a multiple of the spin flip gap. (left) Only one pair of defects is created during the move. At each step the loop growth on either of its sides. (right) Moves that create two pairs of defects.

Figure D. 2 :

 2 Figure D.2: Two types of spin flips sequences that give rise to a 6-loop flip. (top) Only one pair of defects is created during the move. At each step the loop growth on either of its sides. (bottom) Moves that create three pairs of defects.

D. 1 . 4 4 F -525γ 6 6lF 6l -3432γ 8 8lF

 14468 Partial HamiltonianThe Hamiltonian obtained from these terms writesH/∆=-20γ 8l + O(γ) 10 ) (D.5)

Figure D. 3 :

 3 Figure D.3: The four different types of plaquette and their quantum corrections due to the eight-order non-trivial diagonal term.

35512 27 γ 6 F (D. 10 )D. 3 .2 8 - 8 ) 4 , 6 = -a 8 ,4 γ 8 F -a 8 5 4 K 4 .

 6103884688854 th order corrections At the eighth order, the same kind of term is obtained for the 4 and 6 plaquettes. The corresponding correction is then of the form H (EXACT DIAGONALIZATION FOR THE L = 4 SYSTEM This is indeed the ground state obtained when diagonalizing the 90 × 90 Hamiltonian on the entire ice-states manifold. It has energy ε 0 = -√

  states

	Figure 2.3: Ground state for the nearest neighbor interactions. For any height offset different
	from h c (λ), the ground state degeneracy is lifted. (left) If the offset is too small, the ground state
	is only two-fold degenerate. It only contains vortex vertices. (right) For large offset, the ground
	state degeneracy grows exponentially with the system linear size, leading to a sub extensive
	entropy. The different ground states are connected through flips of vertical or horizontal lines.

  .26) We use again the Trotter formula (7.3) and work in the S z basis, introducing two sets of closure relations for each Trotter step. The partition function then reads |N and | N the two Néel states for a plaquette, and | the plaquette state obtained by flipping all the spins of . With these notations, we have

	Z≈	|S	S |	e -F e e -F o	M |S =	|S 2p ,|S 2p+1		S 2p	M p=1	e -F e	S 2p+1	S 2p+1	e -F e	S 2p+2	with = K/M . (7.27)
	The basis can be decomposed into plaquette vectors | = of , so that |S = even | = odd |. The partition function is then expressed as S z i 1 , S z i 2 , S z i 3 , S z i 4 , with {i 1 , i 2 , i 3 , i 4 } the four sites
							M
				Z =	|S p	p=1	even	2p	e -F			2p+1	odd	2p+1	e -F	2p+2	(7.28)
	We now just need to evaluate w( 2p , 2p+1 ) ≡	2p	e -F	2p+1	. We note

Remerciements

Appendix A

Long-range interactions with periodic boundary conditions

As discussed in section 2.2.2, the tails of the long-range interactions have some non-negligible effects. In a finite-size system however, the range of the interactions is by definition limited. This is in particular the case of any Monte Carlo calculation. The interaction range is de facto truncated at the maximal distance between two sites of the sytem. Furthermore, it is standard procedure in a numerical analysis to use periodic boundary conditions. Such a choice allows to separate finite-size effects from boundary effects. The system can be viewed as an infinite tiling of space with the same periodic pattern. To comply with this image in the presence of long-range interactions, one has to replace the interaction energy of the finite-size system by that of the infinite periodic system by replacing the interaction between two sites (1 and 2) at a given distance with the interaction between site 1 and all the replicas of site 2 in the infinite periodic tiling.

Fortunately, the dipolar interactions are integrable in two dimensions. It is then possible to truncate the interactions to a large but finite range. The effective interaction between two spins is obtained by adding n shells of copies of the system around the original one There are then (2n + 1) 2 replicas for each spin, and the effective interaction of a spin i with a spin j (red arrow on the left) is then, up to a factor (2n + 1) 2 the sum of the interactions between i and all the copies of j (red arrows on the right).

First, we add n shells containing periodic copies of the system. The pth shell contains 8p copies. We note σ i,p the spin i in the copy p of the system and J pq ij the interaction between σ i,p and σ j,q . The total effective energy is then

Using the periodicity of the system (σ i,p ≡ σ i ), we get

where the last equality comes from the periodicity of the large system. Thus the effective coupling we used in the Monte Carlo calculations reads, for two spins i and j

In practice, we take n so that nL ∼ 100. This ensures that the truncation does not miss any physically relevant part of the tail of the interactions.

Appendix B

Spin-wave theory for general Ising Hamiltonians in a transverse field

We present here several general formulas to study the quadratic quantum fluctuation in a generic transverse field Ising system. We consider a generic classical ground state with long-range magnetic order, and with a magnetic unit cell containing n spins. We denote S l,p the p-th spin (p = 1 ...n) of the l-th cell. The most general Hamiltonian supporting such a ground state has the form

Here r l is the position of a reference site in the l-th unit cell, and J(∆r) is a n × n matrix containing the couplings between spins in unit cells at a distance ∆r.

In the classical limit an applied transverse field rotates the p-th spins around the y-axis by an angle ϑ p . We introduce a local rotation of the spin configuration, Sl,p = σ p R y (σ p ϑ p )S l,p , where σ p = ±1 is the orientation of the spin in zero field. The rotation has the effect of reducing the S = ∞ ground state to a perfectly ferromagnetic one. The classical energy of the p-th spin of each cell has the expression

[J(∆r)] pp σ p cos ϑ p -S Γ sin ϑ p , so that the total classical energy can be written as

where we have introduced the matrix [ε cl ] p,p = ε cl,p δ p,p . We then consider small quantum fluctuations around this classical reference state, by transforming the quantum spins to bosons via a linearised Holstein-Primakoff transformation [START_REF] Holstein | Field dependence of the intrinsic domain magnetization of a ferromagnet[END_REF] valid in the limit of a small number of bosons n l,p 2S :

Here a l,p and a † l,p are bosonic operators, satisfying [a l,p , a † l,p ] = 1 and [a ( †) l,p , a ( †) l,p ] = 0. The angles ϑ p are chosen so that the classical reference state is stable. Thus the linear terms in the bosonic operators vanish.

Appendix C

Energy and magnetization in the case of Néel and collinear states Appendix E Exact diagonalization for the L = 4 system Let us consider the case of a finite system with L = 4. There are only 90 ice-rule states for this system. Therefore the Hamiltonian is easily diagonalized. As the winding number w is conserved, the Hamiltonian is block diagonal, each block corresponding to a topological sector of fixed w. In the w = 0 sector, there are 30 states that can be labeled in terms of the distance d N measuring the minimal number of plaquette flips connecting them to the Néel state (they are represented in figure E.1).

Ground state

The ground state of the model lies in the w = 0 sector. In this sector there are 2 Néel states (d N = 0), 16 states with d N = 1 and 12 states with d N = 2. We construct the symmetric modes of each class :

Each of these states preserves the translational invariance of the Hamiltonian. We then look for the ground state |0 as a linear combination of the a i

Using this ansatz to solve the Schrödinger, we get