Thermodynamics of the classical spin-ice model with nearest neighbour interactions using the Wang-Landau algorithm

In this article we study the classical nearest-neighbour spin-ice model (nnSI) by means of Monte Carlo simulations, using the Wang-Landau algorithm. The nnSI describes several of the salient features of the spin-ice materials. Despite its simplicity it exhibits a remarkably rich behaviour. The model has been studied using a variety of techniques, thus it serves as an ideal benchmark to test the capabilities of the Wang Landau algorithm in magnetically frustrated systems. We study in detail the residual entropy of the nnSI and, by introducing an applied magnetic field in two different crystallographic directions ([111] and [100]), we explore the physics of the kagome-ice phase, the transition to full polarisation, and the three dimensional Kasteleyn transition. In the latter case, we discuss how additional constraints can be added to the Hamiltonian, by taking into account a selective choice of states in the partition function and, then, show how this choice leads to the realization of the ideal Kasteleyn transition in the system.


Introduction
In magnetism, when a spin cannot fully minimise its interactions with its neighbours, the system is called "frustrated". This situation can arise under a number of different circumstances such as bond disorder, further neighbour interactions and lattice geometry. Geometrically frustrated magnets are the cleanest and better controlled in experimental systems. Frustration precludes the formation of simple ordered ground states. Rather, it typically leads to a degenerate manifold of ground states which are unstable to small perturbations. Frustrated systems exhibit a rich variety of behaviour, including order by disorder, fractionalisation and magnetic analogues of solids, liquids, glasses, ice, quantum liquids and Bose condensation. They represent ideal model systems for the study of generic concepts relevant to collective phenomena, where simple classical Hamiltonians can give rise to a wealth of different phenomena [1][2][3]. All this make geometrically frustrated systems the focus of attention of both theoretical and experimental research, making for ideal test-grounds for numerical techniques.
Among the developments in classical simulation techniques, the Wang-Landau algorithm (WLA) has stood out a e-mail: sag2@st-and.ac.uk as a powerful tool for the determination of phase diagrams [4,5]. The aim of this technique is to give a very accurate estimate of the density of states (DoS) of the system, from which the thermodynamic behaviour can be calculated using canonical ensemble statistical mechanics. There are two particular points of this technique that make it attractive for the study of frustrated systems: first, by simply providing an estimate of the DoS it gives direct information about the degeneracy of the ground state and its low temperature excitations; second, in its algorithmic construction it has a built in mechanism to avoid the trapping of the system into local minima in the energy landscape. The existence of many local minima in frustrated magnets is one of the main reasons of the under-performance of the usual Monte Carlo techniques. Recently, the WLA has been successfully applied to specific problems in frustrated systems and dimer models, such as the formation of magnetisation plateaus in the Ising Shastry-Sutherland model [6] or to determine the order of a transition in the Heisenberg stacked triangular antiferromagnet [7] and in dimer models with next-nearest-neighbour interactions [8]. In the present work we explore thoroughly a simple classical frustrated model, the nearest-neighbour spin-ice model (nnSI), by means of the WLA. This model, despite its apparent simplicity -Ising spins on a pyrochlore lattice interacting ferromagnetically at nearest neighbours -exhibits an extremely rich behaviour. The degeneracy of the zero field ground state grows exponentially with the system size and has construction rules analogous to those of protons in water ice. By applying an external magnetic field it is possible to find a metamagnetic transition or to tune into an effective two dimensional model (the kagome-ice, also with extensive residual entropy) and to find two-and three-dimensional Kasteleyn transitions. The physics of this model system has been well studied using a variety of techniques [9][10][11][12][13], and thus serves as an ideal benchmark to test the capabilities of the WL algorithm in this type of systems. Beyond that, the WLA provides new thermodynamic results that have not been obtained by other methods, such as the free energy as a function of the order parameter, which, to our knowledge, has been calculated for a Kasteleyn transition for a first time.
The structure of this article is as follows. In the next section we give a brief discussion of the model and the simulation technique for completeness. In Section 3, we study the model at zero applied field, calculate the residual entropy and compare it with different existing estimates using the WLA technique. In the same section, we study the system under field applied along [111], and explore the entropy of the kagome-ice state, as well as the entropy peak that rises when this state is destroyed by increasing the temperature. Finally, we study the case of field applied along [100], and do a characterisation of the three dimensional Kasteleyn transition into the fully polarised state.

Wang-Landau algorithm
In 2001, Wang and Landau introduced an algorithm to estimate the density of states of a system by performing a random walk in energy space [4,5]. This method is closely related to "umbrella sampling" techniques [14] and multicanonical Monte Carlo [15]. The algorithm provides a very good estimate of the DoS of the system over a bounded region of the energy spectrum. The DoS can be calculated as a function of the energy, if one works in the canonical ensemble, but also as a function of other variables like pressure or magnetisation if one is interested in other ensembles such as the isothermal-isobaric or its magnetic equivalent (as we will use in Sects. 3.2 and 3.3). In this section we describe the procedure for one variable, which is then straightforwardly extended to the case of several variables.
The algorithm requires the knowledge of the Hamiltonian of the system and a method to sample configurations − in our case it is simply random single spin flips. The starting point is an arbitrary initial configuration with energy E 0 . Initially, the DoS is taken as homogeneous: g(E) = 1 for all E. One step of the calculation consists then in choosing new random configuration, calculating its energy E 1 , and accepting it or discarding it with probability At each step, the histogram H(E) and the DoS g(E) of the final configuration are modified according to g(E) → g(E)f and H(E) → H(E) + 1. Initially the modification factor f for g(E) is taken from a larger value (usually, f 0 = e 1 ) which is then reduced as the algorithm progresses. The rule by which the modification factor f i is reduced is an important choice that conditions the accuracy of the DoS and speed at which the process converges. In principle, one can choose any function that tends monotonically to one, and stop the process once f reaches a given value. In the original work, Wang and Landau propose the rule f k+1 = √ f k , but other choices are possible which give faster convergence. In particular, Belardinelli and Pereyra (BP) [16] proposed an alternative method with improved convergence times, in which the modification factor is eventually scaled as the inverse of the Monte Carlo time, t. For this work we have adopted the BP algorithm.
In our case, within one Monte Carlo step this procedure was repeated with single spin flips until the spin of each site had been chosen to change at least once on average. In practice, the quantity ln g(E) → ln g(E) + ln f is used for simplicity (hence the choice f 0 = e). Notice that this procedure guarantees that the random exploration in phase space will not jam at local minima: each time the transition to the new state is not accepted, the DoS of the initial state is increased and thus the probability of accepting any future transition, which will be proportional to g(E 0 ), is also increased. The detailed balance condition is satisfied to within ln f accuracy.
During this random walk through phase space, the histogram is accumulated and checked periodically. When H(E) becomes flat, the histogram is reset, and the next random walk begins with a finer modification factor f 1 . The criterion of flatness varies according the size and complexity of the system. Usually, the criterion used for flatness is that every entry in H(E) is not smaller than a percentage of the average histogram for all E.
The final result is a relative DoS. To calculate the absolute values one needs some additional information, e.g. a known point in the DoS -in our case the high temperature value of the DoS tends to 2, or the knowledge of the integral of the DoS -in our case 2 N . Once the DoS of the system is known, it is straightforward to calculate the thermodynamic quantities in the canonical ensemble: for example, from the free energy can be written as: the entropy and the specific heat, using the usual linear fluctuation relation,

Spin ice
The nnSI Hamiltonian under an external magnetic field reads: where the S i are Ising spins situated at the corners of a pyrochlore lattice (see inset of Fig. 1), H is the external magnetic field, g the gyromagnetic ratio and J eff is the effective exchange interaction which is taken to be positive. The model is applicable to a certain class of materials, of which the most notable examples are Dy 2 Ti 2 O 7 and Ho 2 Ti 2 O 7 . In these materials, the magnetic ions sit at the corners of a pyrochlore lattice and are constrained by the crystalline field to point along the local 111 quantization axes (i.e. to or from the centre of the tetrahedra). The nnSI model provides a very good description of these materials between 0.2 K and 10 K [17]. The main difference arises from the fact that in the materials, in addition to the exchange interaction -which is antiferromagneticthere is a large long range dipolar interaction. If the latter is truncated beyond the nearest-neighbour spins, J eff is effectively ferromagnetic [18].
The ground state of the nnSI scales exponentially with system size and obeys the local construction rule that within any tetrahedra, two of its spins should point inwards, and two outwards [17,18]. This rule is called the "ice-rule" given its analogy to the Bernal and Fowler rules for protons in water-ice [19] and is the origin of the epithet "spin-ice" given to these models. The exponential degeneracy of this ground state leads to an extensive residual entropy at zero temperature.
We analysed the nnSI model using the WLA with both the original implementation and the one proposed by Belardinelli and Pereyra. In order to reach larger lattices we also performed simulations dividing the energy range in multiple regions. To normalize the DoS, we generally used the condition that there are only two states with maximum energy (the "all-in" and "all-out" configurations). The comparison between several runs performed with different seeds, or different normalisation, allowed us to estimate the errors on the residual entropy for different sizes. We explored the configuration space through random single spin-flip moves. In addition, we used a conventional cubic cell for the pyrochlore lattice, which contains 16 spins, and simulated systems with L×L×L cells, with L ranging from 1 to 8 (16 to 8192 spins). The DoS was estimated as a function of energy with the modification factor changing from f 0 = e 1 to f final = exp(10 −9 ).
For the results of Section 3.1 we used the singlevariable WLA, while for Sections 3.2 and 3.3 we accumulated the DoS as a function of both the energy and the magnetisation in the direction of the applied magnetic field. In this latter case, the thermodynamic quantities are calculated using a partition function which is summed in the energy and in a variable conjugate to the magnetic field, in this way the derivatives of Z provide information about the average magnetisation. We have chosen J = 1.11 K in order to match the effective nn exchange constant for Dy 2 Ti 2 O 7 .

Spin-ice with no applied external field: residual entropy
As a first step we analysed the nnSI model with no applied external field (H = 0), in particular, we were interested in the behaviour of the specific heat and the entropy. In Figure 1 it is shown the logarithm of the DoS of the nnSI model for a system size L = 4. Its shape, in contrast with the highly symmetrical DoS of the ferromagnetic Ising model (see e.g. Ref. [5]), is characteristically asymmetric, starting from a high value at the lowest energy, a direct measure of the high degeneracy of the ground state.
In the case of the specific heat, the expectation is that as the system is cooled down, there will be an onset of correlations as the temperature T approaches J/k, that will be evidenced by a Schottky-like peak [17]. The specific heat calculated using the WLA shows a peak close to T = 1 K (green crosses in Fig. 2), in accordance with those calculated by other techniques [9] and with the high temperature experimental results in spin-ice materials (T > 0.6 K) [20]. This same figure shows the temperature dependence of the entropy per mole of the system. The high temperature value (not shown in the figure) is R ln 2 ≈ 5.76 J/mol K, indicating that the system is behaving as an uncorrelated paramagnet. As the temperature is lowered the entropy is reduced until it reaches a residual value S 0 close to 1.7 J/mol K.
The residual entropy is a characteristic feature of ice models; in real spin-ice materials, such as Dy 2 Ti 2 O 7 or Ho 2 Ti 2 O 7 it is expected that the degeneracy of the spinice manifold will be lifted by additional interactionschiefly the dipolar interaction -and that the system at T = 0 will be ordered [21,22]. In the nnSI model, however, the construction rules are strictly those of Bernal and Fowler, and one expects to find an exponentially degenerate state with the same extensive residual entropy of the three dimensional ice models. The determination of the value of this entropy has a long history, starting with Linus Pauling's famous estimate in 1935 [23]. S 0 can be written as: where k is Boltzmann's constant and W N = W NT is the number of microstates that form the ground state of the system, with N T the number of tetrahedra. Pauling's estimate gives W = 3/2, which translated to the entropy per mole gives S 0 = R/2 ln 3/2 ≈ 1.68 J/molK, in reasonable agreement with the results of Figure 2.
Pauling's estimate neglects correlation effects (in the form of loops) and it can be shown that it is a lower bound on the true S 0 [24]. While an exact solution exist for two dimensional ice models [25], this is not true for three dimensions. Currently the best estimate of the entropy for three dimensional ice is that due to Nagle [26], who, building on work by Di Marzio and Stillinger [27], used a series expansion method to derive the estimate W Nagle = 1.50685 (15).
(4) Fig. 3. The number of microstates of the ground state, W , as a function of the inverse of the number of tetrahedra 1/NT . As expected, the residual entropy decreases as the system size is increased. Shown for reference are the thermodynamic value of W from Pauling's estimate (blue square) and from Nagle's calculation (red circle). The line is a fit according to equation (5).
The inset shows the full range of the fit (from L = 1 to L = 8).
Corrections to Pauling's entropy have been studied in the specific case of the pyrochlore lattice, such as the Numerical Linked Cluster expansion of reference [28]. The most accurate calculation of S 0 for the nnSI model in the literature comes from the integration of energy and magnetisation data obtained by loop Monte Carlo simulations [10] and gives W I = 1.5071 (3), very close to Nagle's result.
The WLA provides a direct determination of the entropy without the need of integrating the specific heat and specifying additional constants, and is thus ideally suited for an accurate determination of the residual entropy of the nnSI model. A variant of the WLA has been successfully applied to determine the residual entropy of two simple nearest neighbours ice models in three dimensional hexagonal lattices: the six-state H 2 O molecule model and the two-state H-bond model [29]. A similar method is applied to nnSI. We have calculated the DoS for lattice sizes L = 1 to 8 which correspond to N T = 8 to 4096 tetrahedra, and from those determined W for the ground state. Figure 3 shows W as a function of the inverse of the number of tetrahedra 1/N T . We have used two different criteria to normalise the DoS: in one case we used a known density for a given state (the highest energy state) and in the other we used the sum of all states (2 N in our system). The differences in W given by the different criteria of normalisation, or that obtained in independent runs using a different set of random numbers, are smaller than the size of the symbols in the figure. Figure 3 shows that, as expected, the residual entropy decreases as the system size is increased. In order to obtain the thermodynamic value of W we perform a least-squares fit this data to the form From this fit we obtain W ∞ = 1.50682 (9), with a 1 = 1.557(9) and θ = 0.883 (3). The sub-linear value of θ is an indication of bond correlations in the ground-state manifold [29,30]. Our value for W in the thermodynamic limit is perfectly consistent with the results by Nagle (see Eq. (4)) and by previous calculations [10,29].

Spin-ice: field along [111]
A remarkable feature of the nnSI model is that an external magnetic field can tune the system into regions of different physics. Two notable cases happen when the field is oriented (i) along the crystallographic [111] direction, which will be described in this subsection, and (ii) along [100], which will be described in the following subsection.
To describe the effects of an applied external field in the WLA, on one hand the Zeeman term must be included in the Hamiltonian (the second term in Eq. (2)) and, on the other hand, it is more convenient to calculate the DoS as a function of two indices: energy and magnetisation M , the latter being the quantity conjugate to the field. In this way is it possible to work in the magnetic equivalent of the isothermal-isobaric ensemble and obtain directly from the doubly summed partition function the Gibbs free energy G(H, T ) and the mean value of M . The disadvantage of this approach is that it is too expensive in calculation resources and, therefore, the calculations are usually constrained on smaller system sizes. Figure 4 shows the logarithm of the DoS as a function of E and M along [111] calculated for a system of size L = 3. As expected, the DoS is still asymmetric along the E direction (cf. Fig. 1) and symmetric in the M axis. The figure also shows the projection of the DoS over the M -E plane, which takes the shape of a pentagon. In this projection one can see that while the highest energy state corresponds to M = 0 the ground state manifold comprises a range of magnetisation values; for H [111] this range goes from −3.33 to 3.33 μ B which is the highest value of the magnetisation that can be obtained along [111] without breaking the ice rules.
In references [10,11] the case of the external field along [111] has been studied for the nnSI model both analyti- cally and through conventional Monte Carlo simulations. Along this crystallographic direction the system can be thought of as a stack of alternating kagome and triangular planes. In the triangular planes the projection of the spins in the [111] direction is 1, while in the kagome planes it is 1/3. The evolution of the system under field can be obtained by analysing the behaviour of the magnetisation at low temperatures. In Figure 5 we show M as a function of the applied magnetic field as calculated using the WLA with two variables (energy and magnetisation) for L = 3 and for the parameters of Dy 2 Ti 2 O 7 , which is in perfect coincidence with all the features found in reference [10]. At very low fields the magnetisation rises linearly with a slope which is proportional to 1/T . In this regime the field is merely selecting a subset with non-zero magnetisation along [111] from the ground-state manifold. The slope is given by the competition between the gain in Zeeman energy and the entropy, because the number of states at a given magnetisation decreases sharply as M increases. This is shown in the inset of Figure 5 where the logarithm of the DoS is plotted as a function of the magnetisation at a fixed energy (k 0.1 K above the ground state). This shows how the logarithm of the number of available states is halved as the magnetisation is raised towards the 3.33μ B /Dy ion. When the field is increased to a value high enough to overcome the entropic effects, but low enough that the ice-rules are still obeyed, the spins in the triangular planes (with projection 1) all align with the field, the kagome planes decouple, and the system becomes effectively two-dimensional termed as "kagome-ice" (KI). The spin-ice rules in the KI still allow for an exponentially degenerate number of configurations, and it still possesses an extensive residual entropy (as shown below). The magnetisation reaches a plateau in this field range at m = 3.33μ B per Dysprosium ion, which is easily calculated by taking into account the ice-rules and the different projections of the spins in the kagome and triangular planes 1 . If the field is increased further, it eventually overcomes the exchange interaction and the system goes through a metamagnetic increase of the magnetisation (at around 1 T in Fig. 5) to a fully polarised state where all spins maximise their projection with the magnetic field (m = 5μ B /Dy ion). In the nnSI model this sudden increase is merely a crossover and its width is strongly temperature dependent. If the dipolar interactions are added to the model, the metamagnetism becomes a first order transition [31].
As mentioned before, the ice-rules in the KI phase do not fully constrain the system and allow for an exponentially large number of possible configurations. In this case it is possible to obtain an exact solution for the residual entropy due to a mapping from the KI into dimers in a honeycomb lattice [11,32]. This mapping allows the study of the effects of slightly misaligned fields, which lead to a two-dimensional Kasteleyn transition (see Ref. [11]), as well as the transition to the fully polarised state, which can be interpreted as a dimer-monomer transition and is expected to be accompanied with a peak in the entropy as a function of field [10].
In conventional MC simulations, the calculation of the field dependent entropy usually includes the integration of the specific heat as a function of temperature, with the integration constant for each field point determined by the value at an appropriate fixed point. In the case of the WLA this calculation is straightforward from the known g(E, M ) and requires no further input.
We have calculated the entropy as a function of field along [111] for different temperatures (see Fig. 6) from the DoS for L = 3. The green curve shows the T = 0 result: here the entropy jumps from its zero field value (the 3D ice-entropy) to the KI value (close to 0.7 J/mol K) and at higher fields (around 1 T) jumps down to zero as the system becomes fully polarised. The red line is the exact value for the KI entropy (S 0 = 0.6715 J/mol K) as calculated in references [11,32]. The slight discrepancy between the calculated S 0 and the exact result is due to finite size effects. It is worth pointing out that despite its name and the fact that the coordination number of the lattice is 4, the KI is not sensu stricto an ice-type model 2 and thus its residual entropy is smaller: W KI = 1.175, compared 1 In the KI configuration, the apical spin is fully aligned with the field, thus contributing μ to the magnetisation. The spins in the kagome plane have a 1/3 projection into [111]; of them, two point with this component towards the field and one points opposite, the net contribution is thus 1/3 μ. The total magnetisation per spin is then: (1 + 1/3)/4 μ ≈ 3.33μB .
2 Ice-models are defined on lattices of coordination number 4 and due to the ice-rules have a degeneracy of six states for a single independent plaquette. In the kagome-ice, there are only three possible states per plaquette. Even its extension, the kagome spin-ice (see Ref. [33]) is only ice-like, in the sense that the vertices (the centres of the triangles) form a lattice of coordination number 3. The T = 0 curve shows how the entropy is reduced to zero through two steps: a first one into the KI state, and a second one (at about 1 T) when the system becomes fully polarised. At higher temperatures the most noticeable feature is a giant entropy peak at the polarisation transition (see text).
to W 2D = (4/3) 3/2 = 1.539 . . . At high temperatures, the most notable feature is the appearance of a giant peak in the entropy, which at low temperatures (below 0.7 K) is larger than the residual entropy, and becomes a peak of infinitesimal width at T = 0 (not seen in the figure due to the discretisation of the data). It might seem counterintuitive that the application of an ordering field might result in an increase of the entropy. This can be explained simply, with the aid of the dimer mapping, as coming from the crossing of an extensive number of energy levels, corresponding to different numbers of monomer defects, which have macroscopic entropies (see Ref. [10]). In real materials though, this feature of the nnSI model, which has potential applications for magnetocaloric manipulations, is almost completely suppressed by additional interactions.

Spin-ice: field along [100]
A case of particular interest arises when an external field is applied in the [100] direction. In this case, contrary to the previous one, the fully saturated state belongs to the zero field ground-state manifold, that is, it satisfies the icerules. At low temperatures (kT J eff ) and for any value of the magnetic field, there are no excitations in the form of local violations to the ice-rules. In references [12,13] it was shown that in this regime the competition between entropy gain and loss of Zeeman energy gives rise to a three-dimensional Kasteleyn transition [34] where strings of negative magnetisation proliferate and span the whole length of the sample. This transition takes place at a field dependent critical temperature given by [12]: below which no string is present. This critical temperature comes from equating in the free energy the loss in Zeeman energy per segment of a string of negative magnetisation (2/ √ 3h, due to the spin projection) with the term arising from the entropic gain per segment (T ln 2). Since the line spans the whole sample, an equal number of in-pointing and out-pointing spins are flipped and the ice-rule is preserved in the whole sample, that is to say, the Kasteleyn transition occurs between different spin-ice states. In the case of our simulation (where we have imposed periodic boundary conditions) these lines take the shape of noncontractible closed loops in the torus. Notice that at low temperatures this process does not involve any violation of the ice-rules. This means that its characteristic energy is completely independent of the value of J eff and that it will thus be present even in the limit J eff /kT → ∞.
The main characteristic of a Kasteleyn transition is its asymmetry: excitations are only possible at the disordered side of the transition. This is shown in Figure 7 where the magnetic susceptibility is plotted as a function of temperature at a series of fixed fields as obtained from our WLA simulation for a system size of L = 3 (circles). There, it is clearly seen that at low temperatures, while the susceptibility tends to diverge when T K is approached from the disordered side -resembling a second-order phase transition -it is flat on the other side -a behaviour more akin to a first-order phase transition. This transition was initially termed as "3/2-order" [35,36] and later as "K-type" [37].
A similar asymmetric peak should be expected in the specific heat, C H . Figure 8 shows C H as a function of temperature for a series of fixed fields as extracted from our simulations. At low fields, it is easy to distinguish two . For low fields, it is easy to distinguish the high temperature Schottky-type peak characteristic of the onset of spin-ice correlations from the low temperature asymmetric peak due to the K-type transition. As the field is raised, the transition moves to higher T and it is gradually affected by the presence of additional local excitations. The lines show the specific heat calculated using only configurations obeying the ice-rule, consequently, the Schottky-type peak is absent.
clearly defined peaks, the first one, at high temperatures corresponds to the Schottky-like peak already mentioned in the H = 0 section corresponding to the onset of spinice correlations in the system. The low temperature peak corresponds to the Kasteleyn transition, and shows the expected sharp edge at low temperatures and gradual rise on the high temperature side. As the field is increased, the K-type transition is moved towards higher temperatures. In our simulations J eff /k = 1.11 K so the condition of J eff kT is no longer satisfied at T K (H) and point like excitations are seen at both sides of the transition, gradually changing the peak into a more symmetrical shape. As point excitations become more important, the simple argument sketched above for the field dependence of T K is no longer valid, and the transition widens and deviates from a linear dependence (see the blue triangles in the inset of Fig. 7).
In the WLA, it is possible to impose additional constraints on the Hamiltonian by selectively choosing the states used to construct the partition function. This means that it is possible to study simultaneously the cases with and without the presence of point defects without the need to introduce any additional mechanism than single-spin flips. In our case in particular, it is very simple for finite size lattices to identify the states that strictly obey the icerule by their magnetisation value. In this way, it is possible to calculate the different thermodynamic quantities for the ideal Kasteleyn transition, isolating the effect of pointlike defects. The results for the susceptibility and specific heat are shown as solid lines in Figures 7 and 8, respectively. They coincide at low temperatures, when kT J eff and the unconstrained curves show a transition at a lower temperature. The most noticeable change is seen in the specific heat, where the Schottky-like peak is absent from the ice-rule obeying curves. Furthermore, if we repeat the analysis to extract T K (H) from this latter set of data, the curve (red triangles in the inset to Fig. 7) follows very closely the theoretical prediction (solid line).
One of the unique characteristics of the WLA is that it makes it simple to calculate the dependence of the free energy as a function of a chosen order parameter -in our case the magnetisation. This provides valuable information regarding the nature of a phase transition, and is particularly interesting for an unusual case such as the K-type transition.
As we have mentioned before, the Kasteleyn transition takes place when J eff /kT is small enough that excitations that break the ice-rule are extremely improbable. In this case, the energy of the system is a constant, the free energy to a purely entropic term and the Gibbs potential is given by G = −T S − M H. This resembles the case of a simple paramagnet, however, as discussed by Jaubert and collaborators (see Refs. [38,39]) a crucial difference arises from the ice-rule constraint: if, contrary to the case of the paramagnet, this constraint brings the entropy to zero at a finite H/kT , this is sufficient to drive a Kasteleyn transition in the system. This ad-hoc supposition can be put to test using the WLA. The inset of Figure 9 shows the behaviour of the entropy per spin, s as a function of the energy in the neighbourhood of s = 0 for a fixed field of 0.05 T. As seen in the figure, the slope at which the entropy vanishes is indeed finite and, furthermore, it is given precisely by 1/T K , with T K the Kasteleyn temperature for this field determined from χ and C.
The main panel of the figure shows the Gibbs potential as a function of the magnetisation, G(M ), at T = 0.2 K for three different fields along [100]: H 1 > H K , H 2 = H K and H 3 < H K as determined using the WLA for a system of L = 4 using only ice-rule configurations. This figure captures the characteristic features expected for a K-type transition. The low field curve resembles that of a paramagnet, with a wide minimum at a non-zero magnetisation. As the field is raised towards H K this minimum becomes wider and moves to higher values of M ,while fluctuations increase. At the critical field H K , the system becomes singular, the minimum sits at M sat and the curve becomes flat (dG/dM = 0) in its neighbourhood. For H > H K , the absolute minimum sits at M sat , the neighbourhood to the minimum is linear, with dG/dM finite and negative, showing the complete absence of fluctuations in the ordered state.

Conclusions
In this work we have explored by means of the WLA the nearest-neighbour spin-ice model, an example of a simple classical frustrated model. We have determined the value of the residual entropy S 0 by doing a finite size analysis of S 0 (L) which can be calculated directly from the DoS determined by the WLA for samples of size L 3 . By including the magnetisation as one of the parameters of the DoS, we were able to calculate the thermodynamic properties of the nnSI model as a function of field. In particular, we have investigated the cases where the field was applied in the [111] and the [100] directions. In the first case, we show how the magnetisation as a function of field calculated with the WLA has all the features characteristic to this direction, namely, a 1/T rise to the kagome-ice plateau at 3.33μ B /Dy ion and a sudden jump at H ≈ 1 T to the fully saturated state. The method allowed for a direct calculation of the field dependent entropy, S(H), without the need of any additional fixed parameter. As expected, S(H) has a plateau at the KI phase with a value that tends to that determined by the mapping of the system into dimer configurations on a honeycomb, and shows a marked peak at the polarisation transition. In a similar fashion, we showed that when a field is applied along [100], a Kasteleyn transition takes place between a S = 0 state and one where line-like excitations pierce the system. Additionally, the WLA provides information not obtained through other methods, such as the free energy as a function of the order parameter, which, to our knowledge, has been calculated for the first time for a Kasteleyn transition. By this mean, we were also able to prove that the ad-hoc assumption that the local constraint of the spin-ice rule, brings the entropy to zero at a finite H/kT , is perfectly valid. We have, furthermore, showed that the WLA allows the computation of the thermodynamic properties of the system when additional constraints -not present in the Hamiltonian -are put into the system. In particular, we showed that by selecting only the ice-rule states, we could calculate the behaviour of the ideal Kasteleyn transition, that is, the one that takes place in the absence of point defects. In summary, we show that the WLA is a very useful tool for simulations of frustrated magnetic systems and provides accurate information for the thermodynamic properties in equilibrium. As a result, quantities such as the entropy and the free energy, that is cumbersome to obtain through other methods, can be easily computed. In addition, the algorithm can be used to study the system under the influence of additional constraints.