Static properties of 2D spin-ice as a sixteen-vertex model

# Static properties of 2D spin-ice as a sixteen-vertex model

Laura Foini, Demian Levis, Marco Tarzia and Leticia F. Cugliandolo
Laboratoire de Physique Théorique et Hautes Energies
Université Pierre et Marie Curie, Paris, France
Laboratoire de Physique de la Matière Condensée
Université Pierre et Marie Curie, Paris, France
June 30, 2019
###### Abstract

We present a thorough study of the static properties of models of spin-ice type on the square lattice or, in other words, the sixteen-vertex model. We use extensive Monte Carlo simulations to determine the phase diagram and critical properties of the finite dimensional system. We put forward a suitable mean-field approximation, by defining the model on carefully chosen trees. We employ the cavity (Bethe-Peierls) method to derive self-consistent equations, the fixed points of which yield the equilibrium properties of the model on the tree-like graph. We compare mean-field and finite dimensional results. We discuss our findings in the context of experiments in artificial two dimensional spin-ice.

## 1 Introduction

Many interesting classes of classical and quantum magnetic systems are highly constrained. In particular, geometric constraints lead to frustration and the impossibility of satisfying all competing interactions simultaneously. In most cases this phenomenon gives rise to the existence of highly degenerate ground states [1, 2], often associated with an excess entropy at zero temperature. Anti-ferromagnets on a triangular lattice and spin-ice samples are materials of this kind.

In conventional spin-ice [3] (see [4, 5] for reviews) magnetic ions form a tetrahedral structure in , i.e. a pyrochlore lattice. This is the case, for instance, of the Dy ions in the DyTiO compound. Since the -electron spins have a large magnetic moment, they can be taken as classical variables at low enough temperature, say, K, behaving as Ising doublets, forced to point along the axes joining the centers of the tetrahedra shared by the considered spin. The origin of geometric frustration in these systems is twofold: it arises from the non-colinearity of the crystal field and the effective ferromagnetic exchange, and from long-range couplings between the spins [5]. Since within this Ising formulation the long ranged dipolar interactions are almost perfectly screened at low temperature [6, 7], in a simplified description only short-range ferromagnetic exchanges are retained [3]. Thus, topological frustration arises from the fact that the Ising axes in the unit cell are fixed and different, forced to point towards the centers of neighboring tetrahedra. The configurations that minimize the energy of each tetrahedron are the six states with two-in and two-out pointing spins.

The ground state of the system is more easily visualized by realizing that each tetrahedron in can be considered as a vertex taking one out of six possible configurations on a lattice of coordination four. With this mapping the magnetic problem described above becomes equivalent to the model of water ice introduced in [8, 9]. In this context, the entropy of the ground state satisfying all the local (two-in - two-out) ice constraints, with all six vertices taken with the same statistical weight, was estimated by Pauling with a simple counting argument [9]. Pauling’s result is very close to the earlier measurements performed by Giauque and Stout [10] on water ice and to the ground-state entropy of the magnetic spin-ice sample measured in the late 90s [11]. Experimentally, the Boltzmann weights of the six vertices can be tuned by applying pressure or magnetic fields along different crystallographic axes. Indeed, the extensions of Pauling’s ice model to describe more general ferroelectric systems lead to ‘ice-type models’ [12]. We will come back to this issue below.

The local constraint leads to many peculiar features that have been studied experimentally, numerically, and analytically. In the ground state, the total spin surrounding a lattice point is conserved and constrained to vanish according to the two-in – two-out rule. This fact has been interpreted as a zero-divergence condition on an emergent vector field [7, 13]. Spins are regarded as fluxes and, quite naturally, an effective fluctuating electromagnetism emerges, where each equilibrium configuration is made of closed loops of flux. This analogy can be used to derive power-law decaying spatial correlations of the spins [14, 15], with a parameter dependent exponent, that were recently observed experimentally via neutron scattering measurements [16]. Such criticality of the disordered ground state, also called spin-liquid or Coulomb phase, had been observed numerically first [17] (see [18] for a more general discussion). A detailed description of the Coulomb phase can be found in [19].

Thermal (or other) fluctuations are expected to generate defects, in the form of ice-breaking rule vertices. In the electromagnetic analogy described above a defect corresponds to a charge, defined as the number of outgoing arrows minus the number of ingoing arrows. As such, a tetrahedra with three-in and one-out spins contains a negative charge of magnitude two and the reversed configuration a positive charge also of magnitude two. In turn, the four-in cells carry a charge minus four and the four-out ones a charge plus four. Such vertices should be present in the samples under adequate conditions. The possibility of observing magnetic monopoles and Dirac strings as being associated to defects has been proposed by Castelnovo et al. [20] and investigated experimentally by a number of groups [16, 21].

Spin-ice can be projected onto Kagome planes by applying specially chosen magnetic fields or pressure. Recently, the interest in spin-ice physics has been boosted by the advent of artificial samples [22] on simpler square lattices. These are stable at room temperature and have magnetic moments that are large enough to be easily observed in the lab, thus giving access to the micro-states that can be directly visualized with microscopy. Following the same line of reasoning exposed in the previous paragraphs, such ice-type systems should be modeled by the complete sixteen-vertex model on a square lattice, where all kind of vertices are allowed.

The exact solution of the ice model (i.e., the six-vertex model) [23], and some generalizations of it in which a different statistical weight is given to the six allowed vertices [24, 25, 26], was obtained by Lieb and Sutherland in the late 60s using the transfer matrix technique with the Bethe Ansatz. Soon after, Baxter developed a more powerful method to treat the generic six- and eight- vertex models [27] and founded in this way the theory of integrable systems (in the eight-vertex model vertices with four in-going or four out-going arrows are allowed). The equilibrium phase diagrams of the six- and eight- vertex models are very rich: depending on the Boltzmann weight of the vertices the systems can be found in a variety of different phases, such as a quasi long-range ordered spin liquid phase (SL), two ferromagnetic (FM) and one (or two) antiferromagnetic (AF) phases, separated by different types of transition surfaces in the phase diagram. In the six-vertex case the SL phase is critical similarly to what is observed in the Coulomb phase.

From a theoretical perspective integrable vertex models are of notable interest. Their static properties can be mapped into spin models with many-body interactions, loop models, three-coloring problems, random tilings, surface growth, alternating sign matrices and quantum spin chains. The Coulomb gas method [28, 29] and conformal field theory techniques [29, 30] have added significant insight into the phase transition and critical properties of these systems as well. A comprehensive discussion of these mappings goes beyond the scope of the present work. For a review the interested reader may consult Refs. [31, 32]. Some comments on these alternative representations, when useful, will be made in the text.

Much less is known about the equilibrium (and dynamic [33, 34, 35, 36]) properties of the sixteen-vertex model in two and three dimensions. At present, the experimental interest in classical frustrated magnets of spin-ice type is focused on understanding the nature and the role of defects (magnetic monopoles) and their effects on the samples’ macroscopic properties [20, 16, 37]. This motivates us to revisit the generic vertex model with defects, i.e. the sixteen-vertex model, and complete its analysis. The special experimental simplicity of two dimensional samples suggests to start from the case. Moreover, it is worth trying to extend at least part of the very powerful analytic toolbox developed for the integrable six- and eight- vertex models to the (non-integrable) complete sixteen-vertex models, where all defects are allowed.

In this paper we describe the equilibrium properties of spin-ice systems with defects. We proceed in two directions. On the one hand, we study the static properties of the sixteen-vertex model on a square lattice by using Monte Carlo (MC) simulations in equilibrium. For simplicity, and in order to keep the model close to experimental realizations, we assume that the statistical weight of the defects remains relatively small (as will be explained carefully in the text). We establish its phase diagram and its critical properties, which we compare to the ones of the integrable models. On the other hand, we introduce a suitable mean-field approximation, defined on a well-chosen tree of vertices, and we adapt the cavity (Bethe-Peierls) method to derive self-consistent equations on such trees, the fixed points of which yield the exact solution of the mean-field model. This method allows us to describe all expected phases and to unveil some of their properties, such as the presence of anisotropic equilibrium fluctuations in the symmetry broken phases. We discuss the range of validity of the mean-field approximation, and we compare the analytical solution to the numerical results for the finite dimensional system. In the conclusions we summarize our results and we briefly discuss some experimental features that we can describe with our methods [38].

## 2 Vertex models

In vertex models the degrees of freedom (Ising spins, valued Potts variables, etc.) sit on the edges of a lattice. The interactions take place on the vertices and involve the spins of the neighboring edges. In this Section we recall the definition and main equilibrium properties of Ising-like vertex models in two dimensions defined on a square lattice.

### 2.1 Definitions

We focus on an square lattice with periodic boundary conditions. We label the coordinates of the lattice sites by . This lattice is bipartite, namely, it can be partitioned in two sub-lattices and such that the sites having even belong to , those having odd belong to , and each edge connects a site in to one in . The degrees of freedom sit on the links between two sites or, in other words, on the “medial” lattice whose sites are placed on the midpoints of the links of the original lattice. The midpoints are hence labeled by and . Thus, in the models we consider the degrees of freedom are arrows aligned along the edges of a square lattice, which can be naturally mapped into Ising spins, say . Without loss of generality, we choose a convention such that corresponds to an arrow pointing in the right or up direction, depending on the orientation of the link, and conversely corresponds to arrows pointing down or left.

### 2.2 The six-vertex model

In the six-vertex model (i.e., spin-ice) arrows (or Ising spins) sit on the edges of a coordination four square lattice and they are constrained to satisfy the two-in two-out rule [8, 9]. Each node on the lattice has four spins attached to it with two possible directions, as shown in Fig. 1. A charge can be attributed to each single vertex configuration, simply defined as the number of out-going minus the number of in-coming arrows. Accordingly, the six-vertex model vertices have zero charge. (Note that the charge is not the sum of the spins attached to a vertex. Such a total spin will be defined and used in Sec. 4.)

Although in the initial modeling of spin-ice all such vertices were equivalent, the model was later generalized to describe ferroelectric systems by giving different statistical weights to different vertices: with the energy of each of the vertices. is the inverse temperature. Spin reversal symmetry naturally imposes that for the first pair of ferromagnetic (FM) vertices, for the second pair of FM vertices, and for the anti-ferromagnetic (AF) ones, see Fig. 1. We have here introduced the conventional names , , and of the three fugacities corresponding to the Boltzmann weight of the three kinds of vertices. In the literature it is customary to parametrize the phase diagram and equilibrium properties in terms of and . This is the choice we also make in this paper. Particular cases include: the F model of anti-ferroelectrics, in which the energy of the AF vertices is set to zero and all other ones are taken to be equal and positive, i.e.  [39]; the KDP model of ferroelectrics, in which the energy of a pair of FM vertices is set to zero and all other ones are equal and positive, e.g.  [12]; the Ice model, in which the energies of all vertices are equal, i.e.  [23]. It is important to note, however, that in the context of experiments in artificial spin-ice type samples, vertex energies are fixed and the control parameter is the temperature. In [38] we used this alternative parametrization and we compared the model predictions to experimental observations [40, 41].

The six-vertex model is integrable. The free-energy density of the model with and periodic boundary conditions was computed by Lieb in the late 60s with the transfer matrix technique and the Bethe Ansatz to solve the eigenvalue problem [23]. The method was then extended to calculate the free-energy density of the F model [24], the KDP model [25], and also models with generic values of , , and , and periodic boundary conditions [26], and the same general case with antisymmetric [42] and domain wall boundary conditions [43]. The effect of the boundary conditions is indeed very subtle in these systems as some thermodynamic properties depend upon them [44] contrary to what happens in usual short-range statistical physics models. Very powerful analytic methods such as the Yang-Baxter equations have been employed in the analysis of these integrable systems [27].

An order parameter that allows one to characterize the different phases is the total direct and staggered polarization per spin111A different order parameter can also be introduced in the eight-vertex model. The latter can be reformulated as an Ising model with nearest, next-to-nearest and plaquette interactions between spins sitting on the sites of the dual lattice  [45]. The usual spontaneous magnetization of the corresponding Ising model defines the order parameter  [27].

 ⟨M±⟩=12(⟨|mx±|⟩+⟨|my±|⟩) (1)

with the horizontal and vertical fluctuating components given by

 L2mx±=∑(m,n)∈A1sm+1/2,n ±∑(m,n)∈A2sm+1/2,n, (2) L2my±=∑(m,n)∈A1sm,n+1/2 ±∑(m,n)∈A2sm,n+1/2. (3)

The angular brackets denote here, and throughtout the rest of this article, the statistical average. A finite value of indicates the spontaneous breaking of symmetry, while a finite value of corresponds to the spontaneous breaking of and translational symmetry.

The four equilibrium phases are classified by the anisotropy parameter

 Δ6=a2+b2−c22ab, (4)

and they are the following.

-Ferromagnetic (-FM) phase: ; i.e. . Vertices and are favored. Spin reversal symmetry is broken. The lowest energy state in the full FM phase is doubly degenerate: either all arrows point up and right or down and left [i.e. , with the magnetization density defined in eq. (1)]. In this phase the system is frozen as the only possible excitations involve a number of degrees of freedom of the order of . In all this phase the exact free-energy per vertex is given by [27]

 fFM=ϵ1. (5)

At () the system experiences a frozen-to-critical phase transition from the frozen FM to a disordered (D) or spin liquid (SL) phase that we discuss below.

-Ferromagnetic (-FM) phase: ; i.e. . This phase is equivalent to the previous one by replacing - by -vertices. The free-energy is and the phase transition towards the SL phase is also of frozen-to-critical type.

Spin liquid (SL) phase: ; i.e. , and . In this phase the averaged magnetization is zero, , and one could expect the system to be a conventional paramagnet (PM). However, the ice constraints are strong enough to prevent the full decorrelation of the spins even at finite temperature. The system is in a quasi long-range ordered phase with an infinite correlation length. At there is a Kosterlitz-Thouless (KT) phase transition between this critical phase and an AF phase with staggered order that is discussed below. Some particular points in parameter space belong to the SL phase as the spin-ice point for which . At this special point the ground state is macroscopically degenerate giving rise to the residual entropy [23]

 S/N=3/2 ln4/3 (6)

with the number of vertices in the sample.

The exact solution found by Baxter yields the free-energy density as a function of the parameters through a number of integral equations [27]. In Sec. 4 we evaluate it numerically and we compare it to the outcome of the Bethe approximation. Close to the FM transitions the free-energy density can be approximated by

 fSL≃max(ϵ1,ϵ3)−12kBT(b+ca−1)=max(ϵ1,ϵ3)−12kBT t2−α, (7)

with being the reduced distance from criticality, , and an exponent that plays the role of the one of the heat-capacity and takes the value here. The first derivative of with respect to the distance from the transition shows a step discontinuity at the SL-FM transition as it would in a first-order phase transition, even though the FM phase is frozen. This corresponds to a critical-to-frozen phase transition.

Antiferromagnetic (AF) phase: ; i.e. . Vertices and are favored. The ground state is doubly degenerate, corresponding to the configurations . The staggered order is not frozen, due to thermal fluctuations. This is confirmed by the exact expression of the staggered magnetization found by Baxter [46, 47]. The free-energy has an essential singularity at the critical temperature (towards the SL phase)

 fAF≃e−cst/√t, (8)

with cst being a constant and the distance from criticality, as it is typical for an infinite order phase transition.

The transition lines are straight lines (given by for the SL-FM and for the SL-AF) and they are shown in Fig. 2 as solid (red) lines. The dashed line along the diagonal represents the range of variation of the parameters in the F model. The horizontal dashed line is the one of the KDP model. The intersection of this two lines corresponds to the ice-model. Although the transitions are not of second order, critical exponents have been defined and are given in the first column of Table 1 for the SL-FM transition. The exponent is taken from the expansion of the free-energy close to the transition, see eqs. (7) and (8). We denote by the critical exponent associated to the order parameter as defined in eq. 1, i.e. the polarization of the arrows on the edges of . The ratios , and are defined using the (divergent in the SL phase) correlation length instead of as the scaling variable [48].

### 2.3 The eight-vertex model

The eight-vertex model is a generalization of the six-vertex model, first introduced to get rid of some of its very unconventional properties due to the hard ice-rule constraint (frozen FM state, quasi long-range order at infinite temperature, etc.) [49, 50]. In this model the allowed local configurations are those for which each vertex is surrounded by an even number of arrows pointing in or out, resulting in the addittion of the two vertices with four ingoing and four outgoing arrows shown in Fig. 3 to the ones in Fig. 1. The eight-vertex model can still be mapped into a loop model on the lattice in such a way that the integrability property is preserved. It was first solved by Baxter in the zero-field case (i.e., with symmetry) [46, 47]. In order to do that he introduced the celebrated Star-Triangle relations (now called Yang-Baxter equations).

The phase diagram is characterized by the anisotropy parameter

 Δ8=a2+b2−c2−d22(ab+cd) (9)

which becomes the six-vertex one when . This model sets into the following five phases depending on the weight of the vertices:

-Ferromagnetic phase (-FM): . Spin reversal symmetry is broken. This ordered phase is no longer frozen and .

-Ferromagnetic phase (-FM): . This phase is equivalent to the previous one replacing - by - vertices.

Paramagnetic phase (PM): . As soon as this phase is truly disordered, with a finite correlation length. The averaged magnetization vanishes .

-Antiferromagnetic phase (-AF): . Translational symmetry is broken. The configurations are dominated by -vertices with an alternating pattern of vertices of type and with defects; .

-Antiferromagnetic phase (-AF): . Translational symmetry is broken. The configurations are dominated by -vertices, with an alternating pattern of vertices and with defects. is also different from zero in this phase. This order parameter does not allow one to distinguish the -AF from the -AF.

The transition lines are given by for the PM-FM ones and for the PM-AF ones. The projection of the critical surfaces on the plane yields straight lines translated by with respect to the ones of the six-vertex model, in the direction of enlarging the PM phase, as shown by the dashed lines in Fig. 2222Beyond the parameter which is crucial in the classification of the phases of the eight-vertex model, a second parameter plays an important role [27] and can be put in relation through the ratio with the values of the critical exponents (see Table 1)..

The effect of the -vertices on the order of the different phase transitions is very strong. As soon as , the KT line between the -AF and the SL phases becomes ‘stronger’, i.e. second order, except at the intersection with the or planes when the transition is first order. On the contrary, the frozen-to-critical lines between the FMs and SL phases get “softer”, i.e. second order, and become KT transitions on the and planes. Finally, the separation between the -AF and disordered phases is second order for and it is of KT type on the and planes. As we will show in Sec. 3, this is consistent with our numerical results 333 In this work we will generically refer to the eight-vertex model as the model where all fugacities, and are different from zero, in contrast with the six-vertex model where (or other particular cases where any of the four Boltzmann weights is zero)..

The critical exponents can be found from the analysis of the free-energy density close to the transition planes. In the -AF regime (referred to as ’principal regime’) they depend explicitly on the weights of the vertices via the parameter  [46]. The critical behaviour close to the -FM-PM transition is obtained from the principal regime by replacing the parameter by and by . The critical exponents for the PM-FM transitions in this model are given in the second column of Table 1. The values of ,   and for the six-vertex model given in the first column are consistent with the eight-vertex model results as the limit (when ).

### 2.4 The sixteen-vertex model

The most general model obtained by removing the ice-rule is the sixteen-vertex model, in which no restriction is imposed on the value of the binary variables attached on each edge of the lattice, and all the vertex configurations can occur. The (eight) three-in one-out and three-out one-in vertices that are added to the ones already discussed are shown in Fig. 4. In order to preserve the symmetry (in absence of an external magnetic field that would break the rotational symmetry), the same statistical weight is given to all these ‘defects’ with charge and . In the figure, vertices are ordered in pairs of spin-reversed couples ( is the spin reversed of and so on) and the difference with the following couples is a rotation by ( is equal to apart from a -rotation and so on).

The new vertices naturally entail the existence of new phases. One can envisage the existence of a critical SL phase for and as this new four-vertex model is equivalent to the dimer model solved by Kasteleyn [51]. It is quite easy to see that -AF stripe order is also possible. For instance, one can build an ordered configuration with alternating lines of and vertices, or another one with alternating columns of and vertices. Phases of this kind should appear if one favors one pair of spin-reversed related vertices by giving them a higher weight than the others, and the transition to this phase should be continuous, since local fluctuations made by elementary loops around a square plaquette are allowed.

The sixteen-vertex model loses the integrability properties [32, 52, 53]. Nonetheless, some exact results are available for a few special sets of parameters. The equivalent classical Ising model has only nearest and next-nearest neighbor two-body interactions when  [32, 52]. In the -AF sector this condition leads to the generalised F model defined by , , and , with the constraint . The model has been solved for the special cases: (i) and (i.e. and ), the associated spin model simplifies into an AF Ising model with only nearest neighbor interactions. This model is known to exhibit a second-order phase transition444The transition occurs at . Using our conjectured (defined below) we get a critical temperature . with a logarithmic divergence of the specific heat (). (ii) Using a different approach it has been shown that for and (i.e. and ) the system also exhibits a second-order phase transition in the same universality class as (i). Notice that the exactly solvable F model is recovered in the limit and . In the same way, in the -FM sector this leads to the generalised KDP model [54] by setting , , , and again . For and (i.e. and ) the system exhibits a second-order phase transition with the same properties of its -AF analog discussed above. For and (i.e. and ) the system also exhibits a second-order phase transition in the same universality class as the previous case.

Since the phase diagram of the generic sixteen-vertex model is rather complex, and our aim is to consider and vertices as defects having a relative small statistical weight, in this work we will focus on the effect of the presence of vertices on the phases and the phase transition lines described in Fig. 2.

### 2.5 Numerical simulations

The numerical analysis of the equilibrium properties of vertex models has been restricted, so far, to the study of the six and eight-vertex cases. Single spin-flip updates break the six- and eight- vertex model constraints and cannot be used to generate new allowed configurations. Instead, as each spin configuration can be viewed as a non-intersecting (six-vertex) or intersecting (eight-vertex) loop configuration, stochastic non-local updates of the loops have been used to sample phase space [55, 56]. By imposing the correct probabilities all along the construction of non-local moves, cluster algorithms can be designed [56, 57]. Non-trivial issues such as the effect of boundary conditions have been explored in this way [56, 58].

Loop-algorithms, as usually presented in the context of Quantum Monte Carlo methods, exploit the world-line representation of the partition function of a given quantum lattice model [59]. It is well known that the six- and eight-vertex models are equivalent to the Heisenberg XXZ and XYZ quantum spin- chains, respectively [49, 60]. It is then not surprising to find the same kind of loop-algorithms in the vertex models literature. A configuration in terms of bosonic world lines of the quantum spin chain in imaginary time can be one-to-one mapped into a vertex configuration on the square lattice, so that the same loop algorithm samples equivalently the configurations of both models.

The loop algorithms could be modified to include three-in – one-out and one-in – three-out defects for the study of spin-ice systems [61, 62]. However, the simultaneous inclusion of four-in and four-out defects makes this algorithm inefficient compared to a MC algorithm with local updates. For this reason, we will use local moves in our numerical studies, as we explain in the next Section.

## 3 Numerical study of the sixteen-vertex model

In this Section we summarize the results obtained for the sixteen-vertex model using MC simulations. We first explain the numerical algorithm. Then, we discuss the phase diagram and the critical properties of the model. All our results are for a square lattice with linear size and periodic boundary conditions.

### 3.1 Numerical methods

We used two numerical methods to explore the equilibrium properties of the generic model; the Continuous time Monte Carlo (CTMC) method that that we briefly explain in Sec. 3.1.1 and in App. A and the Non-equilibrium relaxation method (NERM) that we equally briefly explain in Sec. 3.1.2.

#### 3.1.1 Continuous time Monte Carlo

The usual Metropolis algorithm, from now on called Fixed Step Monte Carlo algorithm (FSMC), is very inefficient to study the equilibrium properties of frustrated magnets. The dynamics freeze when and the acceptance probability for most of the single-spin flip updates is extremely small. We implemented an algorithm which overcomes this difficulty, the CTMC algorithm [63, 64]. This method is also known with different names: Bortz-Kalos-Lebowitz [63], n-fold way or kinetic MC. The basic aim of this algorithm is to get rid of the time wasted due to a large number of rejections when the physics of the problem imposes a very small acceptance ratio. It is extremely useful for the study of the long time behavior of systems with complicated energy landscapes and a large number of metastable states. The main idea behind the method is to sample stochastically the time needed to update the system and then do it without rejections. In our case, we chose to use single spin updates such that the two vertices connected by the spin are updated in the move. Details on this algorithm are given in App. A. Equilibrium can be achieved for relatively large samples. The finite size scaling analysis of the equilibrium MC data yield the thermodynamic properties of the generic model.

#### 3.1.2 Non-equilibrium relaxation method

The fact that dynamic scaling applies during relaxation at a critical point [65] suggested to use short-time dynamic measurements to extract equilibrium critical exponents with numerical methods [64, 66, 67]. With this method it is not needed to equilibrate the systems and only short-time scales are evaluated. Therefore, it is not necessary to use the CTMC version but a plain MC is sufficient. We parametrize the consecutive steps of the MC simulation with a parameter measured in MC step units, each of these corresponding to a sweep of the single spin flip MC algorithm over spins taken at random on the sample. Magnetized, , and non-magnetized, , configurations can be used as starting conditions and the critical relaxation

 M±(t)≃t−βe/(νz) F(tx0/zM0±) (10)

where is the dynamic critical exponent, and for and for , can be used to extract either the critical parameters or the critical exponents. This expression is expected to hold for and with the equilibrium correlation length.

### 3.2 Phase diagram and critical singularities

In this subsection we present a selected set of results from our simulations, and we describe the kind of phases and critical properties found. The strategy to study the different phase transitions is the following. We first chose the relevant order parameter, or , to study FM or AF phases. From the finite size scaling analysis of the corresponding fourth-order reduced Binder’s cumulant

 KM±=1−⟨M4±⟩3⟨M2±⟩2≃ΦK(tL1/ν) (11)

where is the distance from the critical point, we extracted the critical exponent . From the maximum of the magnetic susceptibility

 χ±=L2[⟨M2±⟩−⟨M±⟩2]≃Lγe/νΦχ(tL1/ν) (12)

we extracted , then as was already known. From the maximum of the specific heat

 CE=L−2[⟨E2⟩−⟨E⟩2]≃Lα/νΦC(tL1/ν) (13)

we extracted , then . Here, with the number of vertices of type and their energy. The direct measurement of is difficult, we thus deduced it from the scaling relation . Finally, we checked hyper-scaling, i.e. whether is satisfied by the exponent values obtained, that we summarize in Table 2 for the SL/PM-FM transition and two choices of parameters.

#### 3.2.1 The PM-FM transition

In order to reduce the number of parameters in the problem we studied the PM-FM transition for the special choice .

As the direct magnetization density defined in eqs. (1)-(3) is the order parameter for the PM-FM transition in the six-vertex model, we study this quantity to investigate the fate of the FM phase in the sixteen-vertex model. In Fig. 5 we show as a function of for and three values of the fugacity (all normalized by ). The data for the model (shown with colored points) demonstrate that the presence of defects tends to disorder the system and, therefore, the extent of the PM phase is enlarged for increasing values of . Moreover, the variation of the curves gets smoother for increasing values of suggesting that the transitions are second order (instead of frozen-to-critical) in presence of defects. The data displayed with black points for the same parameters are the result of the mean-field analysis of the model, which will be discussed in Sec. 4. In [33] we suggested that the equilibrium phases of the sixteen-vertex model with could be characterized by a generalization of the anisotropy parameter of the eight-vertex model recalled in eq. (9):

 Δ16=a2+b2−c2−(d+3e)22[ab+c(d+3e)]. (14)

In the same way as in the integrable cases, the proposal is that the PM phase corresponds to the region of the parameters’ space where , the FM phases corresponds to , and the AF ones to . It follows that the projection of the FM-transition hyper-planes onto the plane should then be parallel to the ones of the six- and eight-vertex models and given by (or equivalently ). As shown in Fig. 5, the numerical results are well fitted by Eq. (14) which is, however, not exact. As we will explain in detail in Sec. 4, our mean-field treatment of the model defined on a tree of single vertices predicts a similar shift of the transition lines given by . The parameter capturing all transition lines is, within this analytic approach,

 Δsv16=a2+b2−c2−d2+2(a+b−c−d)e2(cd+ab+e(a+b+c+d+2e)) (15)

(see Sec. 4 for the technical details). For the more sophisticated tree made of ‘plaquette’ units the transition lines are not parallel to the ones of the six- and eight-vertex models and an analytic form of the anisotropy parameter, , has not been found. Nevertheless, the transition lines can be computed numerically and their evaluation leads to the phase diagram depicted in Fig. 19.

Further evidence for the transition becoming second order comes from the analysis of the fourth-order cumulant defined in eq. (11). Raw data on such Binder’s cumulant across the FM-PM transition were shown in [33] where we showed that they intersect at a single point, as expected in a second order phase transition. In Fig. 6 we display raw data for (a) and scaled data for (b) as a function of . In both cases and, as above, we normalize all fugacities by . From the analysis of the scaling properties we extract for and for . Sets of data for linear system sizes are scaled quite satisfactorily by using for the small and for the large value of .

In order to complete the analysis of this transition we studied the magnetic susceptibility (12) associated to the direct magnetization and its finite size scaling. Figure 7 displays for , (a) and (b), and five linear sizes, . The data collapse is very accurate and it allows us to extract the exponent in both cases. The study of the maximum of displayed in the insets confirms this estimate for . We repeated this analysis for other values of and we found that in all cases critical scaling is rather well obeyed and, interestingly enough, is, within numerical accuracy, independent of . This is similar to what happens in the eight-vertex model where, as shown in Table 1, this ratio is independent of the parameters.

The ratio is obtained from the finite size analysis of the specific heat (not shown) that is consistent with (instead of for a first order phase transition). We found for and and a logarithmic divergence of the heat capacity, i.e. for (cf. Table 2), although it is very difficult to distinguish numerically a logarithmic divergence from a power-law one with a very small exponent.

The critical exponents extracted numerically at the second order phase transition with are compared to the ones of the six-vertex model and the Ising model in Table 2. It is interesting to note that for very small value of the exponents are rather close to the ones of the six-vertex model while for large value of they approach the ones of the Ising model. In terms of a Renormalization Group (RG) approach, this behavior suggests the existence of two fixed points, one in the plane describing the critical behavior of the six-vertex model, and another for , belonging to the Ising universality class. The first one might be unstable as soon as an infinitesimal amount of defects is allowed. A RG treatment of the model is required to confirm this guess.

The critical exponents obtained from the numerical analysis depend on the fugacity , namely, they vary along the transition lines, as it happens for the eight-vertex model. However, as shown in Table 2, the ratios of critical exponents , and do not depend qualitatively on the choice of the parameters and are equal to the ones of the Ising model.

The NERM yields values of the critical that are in agreement (within numerical accuracy) with the ones found with the conventional analysis. We do not show this analysis here.

#### 3.2.2 The c-AF-PM transition

We now focus on the transition between the -AF and PM phases. For the six-vertex model this is a KT transition while for the eight-vertex model it is of second order as soon as . In this case we chose to work with and with . We present data obtained with the NERM.

Figure 8 shows the relaxation of the staggered averaged magnetization at and different values of given in the caption. The power-law relaxation, typical of the critical point, is clearly identifiable from the figure. We extract the critical value for the -AF-PM phase transition. Moreover, the data demonstrate that the PM phase is not of SL-type as soon as a finite density of defects is allowed. Indeed, the relaxation of does not follow a power law in the PM phase, the decay being exponential for . Although this strategy gives a rather precise determination of , it is very hard to determine the value of the exponent , and hence of , with good precision as is extremely sensitive to the choice of .

The standard analysis of the -AF-PM transition is not as clean as for the FM-PM one. Figure 9 (a) shows the scaling plot of the Binder cumulant of the staggered magnetization for and, in this case, . From it one extracts . The susceptibility fluctuates too much to draw any stringent conclusion about the exponent . The analysis of the specific heat (not shown) suggests a logarithmic divergence . The dependence of the critical line on the fugacities is reasonably well described by .

Both and vertices make the disordered phase be a conventional PM, and the transitions be second order. Although does not break integrability while does, their effect, in these respects, are similar. The fact that the critical exponents in the eight-vertex model depend on the fugacities is known from the exact solution. In the sixteen-vertex model, where integrability is lost, this is not the case. Unfortunately, our numerical analysis does not allow us to draw stringent enough conclusions for the -AF-PM transition, since it is very hard to get precise measurement of the observables close to the critical lines. Possibly, a RG approach would be useful in this respect.

## 4 Bethe-Peierls mean-field approximation: Vertex models on a tree

In this Section we study the properties of the six-, eight- and sixteen-vertex models defined on the Bethe lattice by using the cavity method. First, we will consider a standard Bethe lattice of uniform connectivity (Sec. 4.2.1). In order to get more accurate results, we also study the model on a tree of plaquettes of vertices (Sec. 4.2.2), which account for local excitations and fluctuations induced on short scales by the presence of small loops.

### 4.1 The cavity method

The cavity method is a technique that allows one to calculate the average properties of statistical models defined on tree-like graphs in the thermodynamic limit [68, 69]. The method, which is equivalent to the Bethe-Peierls (BP) approximation, is based on the assumption that due to the tree-like structure of the lattice, in absence of a given site (the cavity), the neighbors of that site are not correlated and their marginal joint probability factorizes.

Removing one site from the graph creates rooted trees, each one being a tree where all the sites of the bulk have the same connectivity , apart from the root which has only neighbors. The evaluation of physical observables is based on the determination of the properties of the site at the root of a rooted tree. Thanks to the above-mentioned factorization property one can write relatively simple recursion equations for the marginal probabilities of the rooted sites. Such equations have to be solved self-consistently, the fixed points of which yield the free energy of the system along with all the thermodynamic observables.

The BP method can be interpreted either as the exact solution of the model on the tree, or as an approximate solution of the original model on the Euclidean lattice. In order to mimic an Euclidean lattice in dimensions the tree should have connectivity . Such an approach takes into account short-scale () correlations and it is expected to give more accurate quantitative results than the standard fully-connected mean-field approach. It is well-known, for instance, that for the Ising model the BP method yields the mean-field critical exponents with a better estimate of than the one obtained on the fully-connected graph.

### 4.2 The trees

The definition of the vertices requires the selection of a particular orientation of the edges adjacent to a given site. We will define “horizontal” and “vertical” edges, each one with two possible orientations. With this procedure we associate a statistical weight to each vertex configuration, even in such non-Euclidean geometry. In the recursion equations this partition will translate into four different species of rooted trees, depending on the “position” (left, right, down or up) of the missing edge.

#### 4.2.1 A tree of vertices

In the models we are interested in, each site is a vertex and its coordination, which is fixed and equal to four, is the number of vertices connected to it. In order to distinguish one type of vertex from another, and to identify all possible phases, we define the analogue of the two orthogonal directions of the Euclidean square lattice: each vertex has four terminals that we call “up” (u), “down” (d), “left” (l) and “right” (r). So far, the vertices were labeled by their positions. Here, for the sake of clarity, we label them with a single latin index, say , , , . Vertices are connected through edges and that link respectively the down extremity of a vertex with the up terminal of a neighboring vertex , or the left end of with the right end of its neighbor . The symbols and denote undirected edges, so that and . In this way, one creates a bipartition of the edges into horizontal (left-right ) and vertical (up-down ) edges. This notation gives a notion of which vertex is above () and which one is below () and similarly for the horizontal direction. See the sketch in Fig. 10.

Each edge is occupied by an arrow shared by two vertices. An arrow defined on the edge is the left arrow for the vertex and the right arrow for the neighboring vertex. A similar distinction holds for the vertical edges. Each arrow, as any binary variable, can be identified with a spin degree of freedom, taking values in . In this construction there are two kinds of spins, those living on horizontal edges, , and those sitting on vertical edges, . Without loss of generality, we choose a convention such that if the arrow points up and otherwise. Similarly, if the arrow points right and otherwise, for the horizontal arrows (see Fig. 10). This is the analog of the convention used for the spin sign assignment in the model (of course, alternative choices of the signs of the spins can be used).

The local arrow configuration defines the state of the selected vertex. With the spin definition given above, we can assign a total spin to each vertex , and define it as the sum of the spins attached to it, , where indicates the neighborhood of . With this assignment, the spin associated to each type of vertices is if the vertex is of type , if it is a one, if the vertex is of kind , for a vertex, and for the ones.

Consider now a site (vertex) at the root of a rooted tree, in absence of an edge with one of its neighbors, say site . There are four distinct rooted trees depending on whether the missing edge with vertex is the one on its left, right, up or down direction. By analogy with the case, one could interpret these rooted trees as the result of the integration of a transfer matrix in four possible directions. This can be emphasized by taking into account the particular direction of the missing edge at the root that we will indicate as follows: , with . For instance, an “up rooted tree” is the one in which the root has no connection to the up terminal of the vertex , i.e. the link is absent. As shown in Fig. 11, such a rooted tree is obtained by merging a left, an up and a right rooted tree (with root vertices , and respectively) with the addition of a new vertex through the links , , (pictorially the transfer matrix is moving down). Similarly, a “left rooted tree” is obtained by merging a down, a left and an up rooted tree, and so on. The Bethe lattice is finally recovered by joining an up, a left, a down and a right rooted tree with the insertion of the new vertex. Equivalently, given a tree, one creates rooted (cavity) trees by removing an edge.

#### 4.2.2 A tree of plaquettes

As we will show in the following, the results obtained by using the tree structure described above compare extremely well to the ones in in many respects. In order to further improve the approximation, in particular relatively to the nature of some transitions, we introduce a Bethe lattice of “plaquettes” (see Fig. 12), where the basic unit cell is not a single vertex, but a square of vertices. This tree is constructed by connecting each plaquette to other four plaquettes (without forming loops of plaquettes), in the same way as we constructed the tree of single vertices. In this procedure one has to be careful with the orientation of each plaquette and of the vertices on it, and with the order between the two edges outgoing from each side of the unit.

In the following we will refer to the first simpler geometry as the “single vertex problem” and to the second one as the “plaquette model”.

#### 4.2.3 Discussion

A mean-field approximation for the pyrochlore spin-ice system in , based on the same tree-like structure of single vertices described in Sec. 4.2.1, has been already employed in [70] and [61]. In these papers no distinction between the orientation (up, down, left, right) of the edges was made. This approach was apt to deal with the SL and FM phases [61] only. Conversely, our approach keeps track of the four different directions and allows us to simultaneously study all possible phases, including the AF ones. Within our approach it is also easier to remove the degeneracy of the vertices by introducing external magnetic fields or other kinds of perturbations.

The single-vertex BP approximation proposed in [70, 61] provides a very good qualitative and quantitative description of the transition towards the frozen FM phase (KDP problem) in pyrochlore spin-ice in . This is due to the absence of fluctuations in the frozen FM phase.555This is also due to the fact that coincides with the upper critical dimension of the problem in the absence of defects (i.e., the six-vertex model): the mean-field BP approximation is almost exact in apart from logarithmic corrections. However, such a single-vertex approximation is not precise enough to describe the unfrozen staggered AF order, due to thermal fluctuations. By constructing the tree of plaquettes we will manage to capture some of these fluctuations, and to obtain a more accurate description of the unfrozen phases.

Let us finally mention that the ice-rule for the six-vertex model, or the “parity” rule for the eight-vertex might be viewed as a particular form of constraint that forces to flip all the variables on an entire loop in order to go from one allowed configuration to another. Such form of constraints has in many respects a strong algorithmic impact, in particular for algorithms that work with local updates. For this reason it has been investigated in the context of combinatorial optimization, for a wide class of problems mainly defined on tree-like graphs, with techniques similar to those adopted here for the tree of single vertices [71].

### 4.3 The six- and eight-vertex models on a tree of vertices

In this Subsection we study the six- and the eight-vertex models on a tree of single vertices. We obtain the self-consistent recursive equations for the marginal probabilities along with their fixed points. We perform the stability analysis of the solutions, compute the free-energy of the different phases, and present the resulting phase diagram.

#### 4.3.1 Recursion equations

We call , with , the probability that the root vertex – in a rooted tree where is the missing edge – be of type . Such probabilities must satisfy the normalization condition

 ∑ζ∈χ8vψiα→jβζ=1 ,  ∀ ⟨iαjβ⟩ . (16)

In the recurrence procedure we will only be concerned with the state of the arrow on the missing edge. As a consequence, on the root vertex of a rooted tree with a missing edge , we define as being the probability that the arrow that lies on the missing edge takes the value . Then,

 ψui=ψid→juv1+ψid→juv3+ψid→juv6+ψid→juv8 , ψdi=ψiu→jdv1+ψiu→jdv3+ψiu→jdv5+ψiu→jdv7 , ψli=ψir→jlv1+ψir→jlv4+ψir→jlv6+ψir→jlv7 , (17) ψri=ψil→jrv1+ψil→jrv4+ψil→jrv5+ψil→jrv8 ,

and

 ψiα→jβ(−1)=1−ψiα→jβ(+1)=1−ψβi . (18)

Clearly, other parameterizations are possible for these probabilities. For instance, one could use an effective field acting on the spin , i.e. , and the recursive equations would be equivalently written in terms of the fields , as shown in Appendix B.

The operation of merging rooted trees allows us to obtain in a self-consistent fashion the set of probabilities associated to the new root in terms of those of the previous generation. In the bulk of the tree one does not expect such quantities to depend on the precise site, since all expected phases are homogeneous. As a result, the explicit reference to the particular vertex index of the probabilities defined in eq. (17) can be dropped. By so doing, the following four coupled self-consistent equations for the probabilities are obtained:

 ψu=^Ψu[a,b,c,d,ψu,ψd,ψl,ψr]=gu(a,b,c,d,ψu,ψd,ψl,ψr)/zu=1zu[aψlψuψr+b(1−ψl)ψu(1−ψr)+c(1−ψu)(1−ψl)ψr+dψl(1−ψu)(1−ψr)]ψl=^Ψl[a,b,c,d,ψu,ψd,ψl,ψr]=gl(a,b,c,d,ψu,ψd,ψl,ψr)/zl=1zl[aψdψlψu+b(1−ψd)ψl(1−ψu)+cψd(1−ψl)(1−ψu)+d(1−ψd)(1−ψl)ψu]ψd=^Ψd[a,b,c,d,ψu,ψd,ψl,ψr]=gd(a,b,c,d,ψu,ψd,ψl,ψr)/zd=1zd[aψrψdψl+b(1−ψr)ψd(1−ψl)+c(1−ψr)(1−ψd)ψl+dψr(1−ψd)(1−ψl)]ψr=^Ψr[a,b,c,d,ψu,ψd,ψl,ψr]=gr(a,b,c,d,ψu,ψd,ψl,ψr)/zr=1zr[aψuψrψd+b(1−ψu)ψr(1−ψd)+cψu(1−ψr)(1−ψd)+d(1−ψu)(1−ψr)ψd] (19)

where , , , and are the fugacities of the vertices, as defined in Sec. 2 (see Figs. 1 and 3), and are normalization constants that ensure the normalization of :

 zα=gα(a,b,c,d,ψu,ψd,ψl,ψr)+gα(a,b,c,d,1−ψu,1−ψd,1−ψl,1−ψr) . (20)

The functions have been defined in Eq. (19). The first term in this sum corresponds to the un-normalized contribution of a spin while the second term is for a spin . For the sake of simplicity in eq. (19) we considered the argument in the functions to be the same for all the directions , although in each equation the field defined along the opposite direction does not appear explicitly.

#### 4.3.2 Fixed points and free-energy

In order to allow for a fixed point solution associated to -AF and -AF staggered order we study eqs. (19) on a bipartite graph. We partition the graph into two distinct subsets of vertices and , such that each vertex belonging to is only connected to vertices belonging to and vice versa. This amounts to doubling the fields degrees of freedom , one for the sub-lattice and the other one for the sub-lattice , and to solving the following set of coupled equations:

 ψα1=^Ψα[a,b,c,d,ψu2,ψd2,ψl2,ψr2],ψα2=^Ψα[a,b,c,d,ψu1,ψd1,ψl1,ψr1] , (21)

with . The FM and PM phases are characterized by , while the AF phases by with .

Considering the solution associated to only, the fixed points are

 ψPM=(ψu=12,ψl=12,ψr=12,ψd=12) ,ψa-FM=(ψu=1,ψl=1,ψr=1,ψd=1) ,ψb-FM=(ψu=1,ψl=0,ψr=0,ψd=1) ,ψc-AF=(ψu=1,ψl=1,ψr=0,ψd=0) ,ψd-AF=(ψu=1,ψl=0,ψr=1,ψd=0) , (22)

as can be checked by inserting these values into eqs. (19) and (21). The spin reversal symmetry implies that is also a solution of the self-consistent equations. For the AF phase this corresponds to the excange of the two sublattices (i.e., the exchange of with ). These solutions exist for any value of the vertex weights and . Conversely, the stability of the solutions, that will be discussed in more detail in Sec. 4.3.4, depends on the fugacities. The numerical iteration of the self-consistent equations, eqs. (19) and (21), confirms that the fixed points given in eq. (22) are the only possible attractive solutions in the different regions of the phase diagram.

In order to calculate the free-energy, it is useful to consider the contributions to the partition function coming from a vertex, an horizontal edge and a vertical edge. These quantities are defined as follows:

 Zv[ψl,ψr,ψu,ψd]=a [ψlψuψrψd+(1−ψu)(1−ψl)(1−ψr)(1−ψd)]+ b [(1−ψl)ψu(1−ψr)ψd+ψl(1−ψu)ψr(1−ψd)]+ c [(1−ψu