Static properties of 2D spinice as a sixteenvertex model
Abstract
We present a thorough study of the static properties of models of spinice type on the square lattice or, in other words, the sixteenvertex 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 meanfield approximation, by defining the model on carefully chosen trees. We employ the cavity (BethePeierls) method to derive selfconsistent equations, the fixed points of which yield the equilibrium properties of the model on the treelike graph. We compare meanfield and finite dimensional results. We discuss our findings in the context of experiments in artificial two dimensional spinice.
Contents
 1 Introduction
 2 Vertex models
 3 Numerical study of the sixteenvertex model

4 BethePeierls meanfield approximation: Vertex models on a tree
 4.1 The cavity method
 4.2 The trees
 4.3 The six and eightvertex models on a tree of vertices
 4.4 The six and eightvertex models on a tree of plaquettes
 4.5 Duality in the eightvertex model
 4.6 The sixteenvertex model on a tree of vertices
 4.7 The sixteenvertex model on the tree of plaquettes
 4.8 The sixteenvertex model: phase diagram and discussion
 4.9 Fluctuations in the symmetrybroken phases
 5 Summary and conclusions
 6 Acknowledgements
 A The CTMC algorithm
 B Equation for the single vertex
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. Antiferromagnets on a triangular lattice and spinice samples are materials of this kind.
In conventional spinice [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 noncolinearity of the crystal field and the effective ferromagnetic exchange, and from longrange 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 shortrange 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 twoin and twoout 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 (twoin  twoout) 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 groundstate entropy of the magnetic spinice 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 ‘icetype 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 twoin – twoout rule. This fact has been interpreted as a zerodivergence 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 powerlaw 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 spinliquid 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 icebreaking 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 threein and oneout spins contains a negative charge of magnitude two and the reversed configuration a positive charge also of magnitude two. In turn, the fourin cells carry a charge minus four and the fourout 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].
Spinice can be projected onto Kagome planes by applying specially chosen magnetic fields or pressure. Recently, the interest in spinice 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 microstates that can be directly visualized with microscopy. Following the same line of reasoning exposed in the previous paragraphs, such icetype systems should be modeled by the complete sixteenvertex model on a square lattice, where all kind of vertices are allowed.
The exact solution of the ice model (i.e., the sixvertex 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 eightvertex model vertices with four ingoing or four outgoing 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 longrange 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 sixvertex 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 manybody interactions, loop models, threecoloring 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 sixteenvertex model in two and three dimensions. At present, the experimental interest in classical frustrated magnets of spinice 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 sixteenvertex 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 (nonintegrable) complete sixteenvertex models, where all defects are allowed.
In this paper we describe the equilibrium properties of spinice systems with defects. We proceed in two directions. On the one hand, we study the static properties of the sixteenvertex 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 meanfield approximation, defined on a wellchosen tree of vertices, and we adapt the cavity (BethePeierls) method to derive selfconsistent equations on such trees, the fixed points of which yield the exact solution of the meanfield 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 meanfield 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 Isinglike 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 sublattices 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 sixvertex model
In the sixvertex model (i.e., spinice) arrows (or Ising spins) sit on the edges of a coordination four square lattice and they are constrained to satisfy the twoin twoout 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 outgoing minus the number of incoming arrows. Accordingly, the sixvertex 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 spinice 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 antiferromagnetic (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 antiferroelectrics, 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 spinice 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 sixvertex model is integrable. The freeenergy 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 freeenergy 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 shortrange statistical physics models. Very powerful analytic methods such as the YangBaxter 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 spin^{1}^{1}1A different order parameter can also be introduced in the eightvertex model. The latter can be reformulated as an Ising model with nearest, nexttonearest 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].
(1) 
with the horizontal and vertical fluctuating components given by
(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
(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 freeenergy per vertex is given by [27]
(5) 
At () the system experiences a frozentocritical 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 freeenergy is and the phase transition towards the SL phase is also of frozentocritical 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 longrange ordered phase with an infinite correlation length. At there is a KosterlitzThouless (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 spinice point for which . At this special point the ground state is macroscopically degenerate giving rise to the residual entropy [23]
(6) 
with the number of vertices in the sample.
The exact solution found by Baxter yields the freeenergy 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 freeenergy density can be approximated by
(7) 
with being the reduced distance from criticality, , and an exponent that plays the role of the one of the heatcapacity and takes the value here. The first derivative of with respect to the distance from the transition shows a step discontinuity at the SLFM transition as it would in a firstorder phase transition, even though the FM phase is frozen. This corresponds to a criticaltofrozen 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 freeenergy has an essential singularity at the critical temperature (towards the SL phase)
(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 SLFM and for the SLAF) 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 icemodel. 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 SLFM transition. The exponent is taken from the expansion of the freeenergy 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 eightvertex model
The eightvertex model is a generalization of the sixvertex model, first introduced to get rid of some of its very unconventional properties due to the hard icerule constraint (frozen FM state, quasi longrange 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 eightvertex 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 zerofield case (i.e., with symmetry) [46, 47]. In order to do that he introduced the celebrated StarTriangle relations (now called YangBaxter equations).
The phase diagram is characterized by the anisotropy parameter
(9) 
which becomes the sixvertex 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 PMFM ones and for the PMAF ones. The projection of the critical surfaces on the plane yields straight lines translated by with respect to the ones of the sixvertex model, in the direction of enlarging the PM phase, as shown by the dashed lines in Fig. 2^{2}^{2}2Beyond the parameter which is crucial in the classification of the phases of the eightvertex 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 frozentocritical 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 ^{3}^{3}3 In this work we will generically refer to the eightvertex model as the model where all fugacities, and are different from zero, in contrast with the sixvertex 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 freeenergy 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 FMPM transition is obtained from the principal regime by replacing the parameter by and by . The critical exponents for the PMFM transitions in this model are given in the second column of Table 1. The values of , and for the sixvertex model given in the first column are consistent with the eightvertex model results as the limit (when ).
sixvertex  eightvertex  

2  
0  
2  
1  
0  
1  
1/2 
2.4 The sixteenvertex model
The most general model obtained by removing the icerule is the sixteenvertex 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) threein oneout and threeout onein 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 spinreversed 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 fourvertex 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 spinreversed 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 sixteenvertex 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 nextnearest neighbor twobody 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 secondorder phase transition^{4}^{4}4The 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 secondorder 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 secondorder phase transition with the same properties of its AF analog discussed above. For and (i.e. and ) the system also exhibits a secondorder phase transition in the same universality class as the previous case.
Since the phase diagram of the generic sixteenvertex 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 eightvertex cases. Single spinflip 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 nonintersecting (sixvertex) or intersecting (eightvertex) loop configuration, stochastic nonlocal updates of the loops have been used to sample phase space [55, 56]. By imposing the correct probabilities all along the construction of nonlocal moves, cluster algorithms can be designed [56, 57]. Nontrivial issues such as the effect of boundary conditions have been explored in this way [56, 58].
Loopalgorithms, as usually presented in the context of Quantum Monte Carlo methods, exploit the worldline representation of the partition function of a given quantum lattice model [59]. It is well known that the six and eightvertex 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 loopalgorithms in the vertex models literature. A configuration in terms of bosonic world lines of the quantum spin chain in imaginary time can be onetoone 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 threein – oneout and onein – threeout defects for the study of spinice systems [61, 62]. However, the simultaneous inclusion of fourin and fourout 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 sixteenvertex model
In this Section we summarize the results obtained for the sixteenvertex 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 Nonequilibrium 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 singlespin 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: BortzKalosLebowitz [63], nfold 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 Nonequilibrium relaxation method
The fact that dynamic scaling applies during relaxation at a critical point [65] suggested to use shorttime 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 shorttime 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 nonmagnetized, , configurations can be used as starting conditions and the critical relaxation
(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 fourthorder reduced Binder’s cumulant
(11) 
where is the distance from the critical point, we extracted the critical exponent . From the maximum of the magnetic susceptibility
(12) 
we extracted , then as was already known. From the maximum of the specific heat
(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 hyperscaling, i.e. whether is satisfied by the exponent values obtained, that we summarize in Table 2 for the SL/PMFM transition and two choices of parameters.
3.2.1 The PMFM transition
In order to reduce the number of parameters in the problem we studied the PMFM transition for the special choice .
As the direct magnetization density defined in eqs. (1)(3) is the order parameter for the PMFM transition in the sixvertex model, we study this quantity to investigate the fate of the FM phase in the sixteenvertex 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 frozentocritical) in presence of defects. The data displayed with black points for the same parameters are the result of the meanfield analysis of the model, which will be discussed in Sec. 4. In [33] we suggested that the equilibrium phases of the sixteenvertex model with could be characterized by a generalization of the anisotropy parameter of the eightvertex model recalled in eq. (9):
(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 FMtransition hyperplanes onto the plane should then be parallel to the ones of the six and eightvertex 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 meanfield 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,
(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 eightvertex 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 fourthorder cumulant defined in eq. (11). Raw data on such Binder’s cumulant across the FMPM 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 eightvertex 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 powerlaw 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 sixvertex 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 sixvertex 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 sixvertex 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 eightvertex 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.
sixvertex  MC ()  MC ()  Ising  
7/4  
1/8  
0  
1/8  
7/4  
1  
?  yes  yes  yes  yes 
3.2.2 The AFPM transition
We now focus on the transition between the AF and PM phases. For the sixvertex model this is a KT transition while for the eightvertex 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 powerlaw relaxation, typical of the critical point, is clearly identifiable from the figure. We extract the critical value for the AFPM phase transition. Moreover, the data demonstrate that the PM phase is not of SLtype 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 AFPM transition is not as clean as for the FMPM 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 eightvertex model depend on the fugacities is known from the exact solution. In the sixteenvertex 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 AFPM 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 BethePeierls meanfield approximation: Vertex models on a tree
In this Section we study the properties of the six, eight and sixteenvertex 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 treelike graphs in the thermodynamic limit [68, 69]. The method, which is equivalent to the BethePeierls (BP) approximation, is based on the assumption that due to the treelike 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 abovementioned factorization property one can write relatively simple recursion equations for the marginal probabilities of the rooted sites. Such equations have to be solved selfconsistently, 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 shortscale () correlations and it is expected to give more accurate quantitative results than the standard fullyconnected meanfield approach. It is wellknown, for instance, that for the Ising model the BP method yields the meanfield critical exponents with a better estimate of than the one obtained on the fullyconnected 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 nonEuclidean 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 (leftright ) and vertical (updown ) 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 meanfield approximation for the pyrochlore spinice system in , based on the same treelike 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 singlevertex 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 spinice in . This is due to the absence of fluctuations in the frozen FM phase.^{5}^{5}5This is also due to the fact that coincides with the upper critical dimension of the problem in the absence of defects (i.e., the sixvertex model): the meanfield BP approximation is almost exact in apart from logarithmic corrections. However, such a singlevertex 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 icerule for the sixvertex model, or the “parity” rule for the eightvertex 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 treelike graphs, with techniques similar to those adopted here for the tree of single vertices [71].
4.3 The six and eightvertex models on a tree of vertices
In this Subsection we study the six and the eightvertex models on a tree of single vertices. We obtain the selfconsistent recursive equations for the marginal probabilities along with their fixed points. We perform the stability analysis of the solutions, compute the freeenergy 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
(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,
(17)  
and
(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 selfconsistent 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 selfconsistent equations for the probabilities are obtained:
(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 :
(20) 
The functions have been defined in Eq. (19). The first term in this sum corresponds to the unnormalized 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 freeenergy
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 sublattice and the other one for the sublattice , and to solving the following set of coupled equations:
(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
(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 selfconsistent 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 selfconsistent 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 freeenergy, 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: