Entropy as a function of magnetization for a 2D spin-ice model exhibiting a Kasteleyn transition

We present a combined analytical and numerical study of the entropy as a function of magnetization for an orientable 2D spin-ice model that exhibits a Kasteleyn transition. The model that we use is related to the well known six-vertex model but, as we show, our representation of it is more convenient for constructing approximate expressions for the entropy at fixed magnetization. We also discuss directions for further work, including the possibility of deforming our model into one exhibiting a quantum Kasteleyn transition.

We present an Ising model on a two-dimensional lattice with mixed ferro-and antiferromagnetic nearest-neighbor and next-nearest-neighbor exchange couplings. We show that this model exhibits 'spin ice' behavior, retaining a non-zero entropy density even at zero temperature. We further demonstrate that, in the presence of an applied longitudinal magnetic field, it undergoes a Kasteleyn (line-defect-driven) transition from a fully magnetized to a 'stringy' phase upon heating. Finally, we discuss directions for further work, including the possibility of deforming our model into one exhibiting a quantum Kasteleyn transition.

I. INTRODUCTION
A popular version of the third law of thermodynamics is that the entropy density of a physical system tends to zero in the T → 0 limit 1 . However, there is a class of theoretical models that violate this law [2][3][4][5][6][7][8] : models in this class exhibit a ground-state degeneracy which grows exponentially with the system size, leading to a non-zero entropy density even at T = 0. Nor can these be easily dismissed as theorists' abstractions, since one also sees ample evidence in experiment 9-12 that there are systems in which the entropy plateaus at a non-zero value over a large range of temperature. In many such cases it is suspected that it eventually falls to zero at a much lower temperature scale, though recent theoretical work on skyrmion magnets suggests that this intuition may not always be reliable 13 .
Whatever the ultimate low-temperature fate of these materials, it is clear that over a broad range of temperatures they exhibit physics which is well captured by models with a non-zero residual entropy density. One important class of these are so-called ice models, in which the ground-state manifold consists of all configurations which satisfy a certain local 'ice rule' constraint [14][15][16] . The first such model was Pauling's model for the residual configurational entropy of water ice 3 . Here the local constraint is that two of the four hydrogens neighboring any given oxygen should be chemically bonded to it to form a water molecule. Similar models were subsequently discovered to apply to the orientations of spins along local Ising axes in certain rare-earth pyrochlores 7,14 , which by analogy were dubbed 'spin ice' compounds. Such models develop power-law spin-spin correlations at low temperatures, with characteristic 'pinch points' in the momentum representation of the spin-spin correlation function 17 , but they do not order. Their low-temperature state is often referred to as a 'co-operative paramagnet' 18 .
One interesting feature of such co-operative paramagnets is their response to an applied magnetic field. The configurations that make up the ice-rule manifold usually have different magnetizations; thus an applied field, depending on its direction, may either reduce 12,19,20 or entirely eliminate 21,22 the degeneracy. In the latter case, further interesting physics may arise when the system is heated, especially if the ice-rule constraints do not permit the thermal excitation of individual flipped spins. In such cases the lowest-free-energy excitation may be a string of flipped spins extending from one side of the system to the other. A demagnetization transition mediated by such excitations is known as a Kasteleyn transition [22][23][24] .
In spin ice research to date, insight has often been gained from the study of simplified models where the dimensionality is reduced or the geometry simplified while retaining the essential physics [25][26][27] . In that spirit, we present in this paper a two-dimensional ice model which exhibits a Kasteleyn transition in an applied magnetic field. The model is especially interesting since, unlike its three-dimensional counterparts, it has the same Ising quantization axis for every spin. This raises the possibility that it could be extended to include a transverse magnetic field, thereby allowing the exploration of quantum Kasteleyn physics.
The remainder of this paper is structured as follows. In section II, we present our spin ice model, along with some analytical and numerical results on its thermodynamic properties in the absence of an applied magnetic field. In section III, we analyse the model in the presence of a magnetic field: we show that it has a Kasteleyn transition, and we characterize it. In section IV, we use an alternative representation of the ice-rule states -the 'string representation' -to determine the model's entropy as a function of its magnetization. Finally, in section V, we summarize our findings and discuss possible future lines of work.

II. THE MODEL
The model that we shall consider has the following Hamiltonian: Here i and j label the sites of a two-dimensional square lattice, σ i = ±1 is an Ising variable on lattice site i, and h is an externally applied (longitudinal) magnetic field.  Fig. 1. The first six configurations listed are those that, in the absence of an external magnetic field, constitute the sixfold-degenerate ground-state (or 'ice rule') manifold.
The exchange interaction J ij is given by: and r j = r i +x +ŷ; −J r i = nx + mŷ (n + m even) and r j = r i −x +ŷ; 0 otherwise, where r i is the position vector of site i,x andŷ are the unit vectors of a Cartesian system in the two-dimensional plane, and J is a positive constant. In this paper, we shall always work in the limit J ≫ |h|, k B T . Furthermore, where necessary we shall take the number of sites in the lattice to be N , always assuming N to be large enough that edge effects can be neglected. When we refer to the density of something (e.g. the entropy density), we shall always mean that quantity divided by the number of spins -not, for example, by the number of plaquettes. The lattice described by (2) is shown in the upperleft inset of Fig. 1, with ferromagnetic bonds represented by solid lines and antiferromagnetic bonds represented by dotted lines. One may view this lattice as made of corner-sharing plaquettes, one of which is shown in the lower-right inset of Fig. 1. It is easy to see that the bonds on this plaquette cannnot all be satisfied at once: the model (1) is therefore magnetically frustrated.
The sixteen spin configurations of the elementary plaquette, together with their energies, are shown in Table  I. When h = 0, i.e. in the absence of an external magnetic field, there are six degenerate ground-state configurations. They are shown in the left-hand inset of Fig. 3: we shall call them the 'ice-rule configurations,' and the manifold spanned by them the 'ice-rule manifold. ' Since this Ising model is magnetically frustrated, we do not expect it to show an ordering transition as the temperature is reduced. Rather, we expect a smooth crossover into a co-operative paramagnetic state in which every plaquette is in one of the ice-rule configurations.
The density of defects (a measure of how many plaquettes are not in an ice-rule configuration) should vanish smoothly as the temperature tends to zero, and the specific heat will show a corresponding Schottky-like peak at temperatures T ∼ J/k B but no sharp features.
Because the ground-state degeneracy is exponential in the system size, the model will have a non-zero entropy density even at zero temperature. A naïve estimate would suggest a value of k B ln 6 per plaquette, i.e. 1 2 k B ln 6 ≈ 0.896 k B per spin, due to the six-fold groundstate degeneracy. This estimate, however, is too naïve, since it ignores the important constraint that the icerule configurations chosen for two neighboring plaquettes must agree on the orientation of the spin at their shared corner.
We may easily improve our estimate of the zerotemperature entropy density by taking this constraint into account at a local level. Imagine 'growing' a spin configuration of the lattice from top to bottom. Each time a new row is added, the orientations of spins 1 and 2 of each plaquette of the row being added (j) will be fixed by the (already chosen) configuration of the row above (j − 1). The ice rules for this model do not favor any particular direction for any single site on the plaquette; hence the probabilities of the four configurations of this pair of spins are simply P ↑↑ = P ↑↓ = P ↓↑ = P ↓↓ = 1/4. The number of ice-rule configurations consistent with these constraints is (see Fig. 3 Thus half the plaquettes in the new row have no choice of configuration, while the other half may choose between two. This gives an average entropy per plaquette of 1 2 k B ln 2, which corresponds to an entropy density of 1 4 k B ln 2 ≈ 0.173 k B per spin. This estimate is still rather crude, since it neglects correlations between the configurations of neighboring plaquettes in row j − 1, which will be induced by their connections to a common plaquette in row j − 2. However, it was shown by Lieb 5 that such correlation corrections may be resummed to yield an exact result for the ground-state entropy density of such 'square ice' models: We shall call this value the 'Lieb entropy density,' and denote it s First, we demonstrate the increasing predominance of ice-rule configurations as the temperature is lowered. For this it is useful to define the number of defects on a plaquette as the number of single spin-flips by which the spin configuration deviates from the closest ice-rule configuration. By this measure, the states in the top line of Table I have zero defects, those in the second line have one, and those in the third line have two. Fig. 1 shows the density of defects as a function of temperature.
The asymptotic high-temperature value of this quantity can be easily calculated. In the infinite-temperature limit all configurations of a plaquette are equally probable, i.e. each has a probability 1 16 . From Table I, we see that there are six configurations with no defects, eight configurations with one, and two configurations with two. Hence the average number of defects per plaquette at infinite temperature is 0 × 6 16 + 1 × 8 16 + 2 × 2 16 = 3 4 . Since there are twice as many spins as plaquettes, the defect density is simply half of this, i.e. ρ defects → 3 8 = 0.375 as k B T /J → ∞.
Second, we calculate the entropy density of the system as a function of temperature, using the Wang-Landau method 28 . The results are shown in Fig. 2. At high temperatures the entropy density tends to k B ln 2, the Ising paramagnetic value. At low temperatures it tends to a non-zero constant value which is in good agreement with the Lieb entropy density s Lieb 0 given above. In between there are no sharp features, confirming that the model exhibits only a crossover from high-temperature paramagnetic to low-temperature cooperative-paramagnetic behavior.
Third, we obtain the specific heat capacity as a function of temperature, also using the Wang-Landau method. The results are shown in Fig. 3. In keeping with our results for the entropy density in Fig. 2, we see that although there is a broad Schottky-like peak at temperatures of order J/k B there are no sharp features, supporting our expectation that this model would not exhibit a phase transition.

III. KASTELEYN TRANSITION
Our real interest in this model, however, is in its unusual response to an externally applied longitudinal magnetic field. We call this field h, and in the following we shall take it to be positive. As shown in the first line of Table I, the degeneracy between the six ice-rule configurations is lifted as soon as the field h is applied.
Indeed, for any non-zero h (and remembering that we always work in the h ≪ J limit) the ground state of a plaquette is the unique 'all up' configuration. It follows that, at T = 0, the entire lattice simply has σ i = +1 for all sites i. Now let us consider what happens to this fully magnetized state as the temperature is increased. One might expect the appearance of a dilute set of 'down' spins. However, a feature of this model is that a single spinflip takes the system out of the ice-rule manifold, and at h, k B T ≪ J this will not occur. To understand what will happen instead, let us introduce a representation of the states in the ice-rule manifold in terms of strings.
We begin with a single plaquette. If we take as our reference state the one in which all the spins are up, we may represent the six ice-rule configurations in terms of lines joining the spins that are down. This is shown in the right-hand inset of Fig. 3. Representing the 'all down' configuration as two vertical lines rather than two horizontal ones is in principle arbitrary, but it has the advantage of yielding a model in which these lines of down spins can neither cross each other nor form closed loops.
To make an ice-rule-obeying configuration of the entire lattice, we must put these plaquettes together in such a way that any string that leaves one plaquette enters its neighbor. Thus there is a one-to-one mapping between ice-rule-obeying configurations of the spins σ i and configurations of these strings. Each string must extend all the way across the lattice.
To proceed further, let us suppose that the lattice consists of L x sites in the horizontal direction and L y sites in the vertical direction, so that N = L x L y . Each string, irrespective of its configuration, contains precisely L y spins, so that a configuration with N s strings has N s L y down spins and thus an energy of 2hN s L y relative to the fully magnetized state (or 'string vacuum'). Such a string is the minimal demagnetizing excitation of the system that is consistent with the ice rule.
Since a single string has an energy cost proportional to the linear size of the system, it might appear that such strings cannot be thermally excited. This is not true, however, because a single string also has two choices about which way to go every time it enters a new plaquette, meaning that its entropy of k B L y ln 2 is also proportional to L y . Thus the free-energy cost of introducing a single string into the fully magnetized state is When the temperature reaches the critical value T c = 2h/(k B ln 2), this free-energy cost flips sign, and the system becomes unstable to the proliferation of strings. (This is somewhat similar to what happens in a Berezinskii-Kosterlitz-Thouless transition 29,30 , except that in our model we do not have 'positive' and 'negative' strings, so the physics of screening plays no rôle.) In fact the increase in the string density from zero for T > T c -which corresponds directly to the decrease in the magnetization from its saturated value -is continuous. This is because the above argument applies strictly only to a single string introduced into the fully magnetized state. Once a finite density of strings has been created the entropy associated with new ones is reduced, and thus the temperature at which it becomes free-energetically favorable to create them goes up. This kind of transition, in which the elementary thermal excitations are system-spanning strings, is called a Kasteleyn transition. It was first described by Kasteleyn in the context of dimer models 23 .
The above predictions are again borne out by our Monte Carlo simulations, the results of which are shown in Figs. 4-6. Fig. 4 shows a three-dimensional plot of the equilibrium value of the magnetization, M , as a function of the temperature and the applied magnetic field. At all temperatures below T c (h) the magnetization takes its saturated value; above T c (h) it decreases smoothly with increasing temperature, tending to zero only as T → ∞. This may be understood in the string representation of the problem. As more and more strings are introduced, the entropy density of each new one decreases; in the limit where half the lattice sites are populated by strings it tends to zero, meaning that this will occur only in the infinite-temperature limit. Fig. 5 shows the magnetic susceptibility, determined at three different values of the applied field. In each case, one sees at T = T c (h) the asymmetric peak characteristic of a Kasteleyn transition. This highlights an intriguing consequence of the physics of the Kasteleyn strings: below T c (h) the linear susceptibility is strictly zero, while as T c (h) is approached from above the susceptibility diverges. For a two-dimensional Kasteleyn transition one expects to find β = 1/2 on the high-temperature side 31,32 , that is, where µ ≡ (M sat −M )/M sat is the reduced magnetization and t ≡ (T − T c )/T c is the reduced temperature. This is indeed the case in our simulations: the inset of Fig. 5 is a logarithmic plot of µ as a function of t, calculated for a system of 8192 spins and with an applied field of h/J = 0.017 (grey filled circles), compared with the expected t 1/2 behavior (solid red line). Similar behavior is found for all simulated fields below 0.1h/J. In Fig. 6 we collect our data into a phase diagram. The filled red circles show the temperature of the Kasteleyn transition, determined from the data in Fig. 4 as the temperature at which the magnetization departs from its saturated value. The thick black line is the prediction T c (h) = 2h/(k B ln 2) derived above. The departure of the red points from this line at larger fields and temperatures is due to the violation of the condition h, k B T ≪ J.
In the pink region the thermal excitations are not full strings, but instead string fragments extending from one ice-rule-violating plaquette to another. The physics of such string fragments, and their signatures in neutron scattering, were discussed by Wan and Tchernyshyov 27 . The red dots show the Kasteleyn temperature as determined from the magnetization curves, i.e. the temperature at which the magnetization first departs from its saturated value. The black line is the theoretical prediction Tc(h) = 2h/(kB ln 2). As expected, the simulation results depart from the theoretical prediction at temperatures where the condition that the spin configuration remain strictly in the ice-rule manifold, h, kBT ≪ J, is no longer fulfilled (pink area).

IV. ENTROPY AS A FUNCTION OF MAGNETIZATION
Finally, let us demonstrate the usefulness of the string representation by using it to calculate the entropy density of the system, s, at a fixed value of the magnetization density, m ≡ M/M sat . Clearly s(m) is an even function of m, so we may restrict our calculation to the case m 0. The magnetization density may equivalently be expressed as the density of strings, η s , via the formula η s = (1 − m)/2.
To determine the entropy density corresponding to a given value of η s , consider propagating the string configuration downwards from the top of the lattice. We shall assume that this propagation has reached a certain row j, and concentrate on a single string in that row. As it enters a new plaquette in row j + 1, this string has in principle two choices: to continue vertically downwards, or to cross the plaquette diagonally. However, if another string is entering the same plaquette, it has only one choice, since the strings cannot cross (see Fig. 3).
The probability that a second string enters the same plaquette in row j + 1 as the first is simply η s . Thus the average number of choices available to the first string upon entering the new plaquette is η s × 1 + (1 − η s ) × 2 = 2 − η s . This means that each string has a total entropy S s ≈ k B L y ln (2 − η s ); with a total number of strings η s L x , it follows that the total entropy is S ≈ In Fig. 7 we compare this approximation with numerical results for the entropy density obtained using the Wang-Landau method. The filled black circles are the numerical results, while the dashed red curve is our analytical approximation (5). It is clear that these were never going to coincide, since the m → 0 limit ofs 0 (m) is the Pauling entropy density, 1 2 k B ln 3 2 , while the m → 0 limit of the actual entropy density is the Lieb entropy density, 3 4 k B ln 4 3 . The origin of the difference between Lieb's exact result and Pauling's approximation lies in positive correlation of closed loops 4,5 , which increases by a small factor the number of possible configurations obeying the ice rule. If one makes the crude assumption that this factor is independent of m, this results in a constant additive change to the logarithm in (5): If we choose the constant α to match the known result at m = 0, the resulting curve (shown in blue) gives a very reasonable fit to the numerical data points over the whole range 0 m 1.

V. SUMMARY AND FUTURE WORK
In this paper, we have presented a new spin-ice model defined on a two-dimensional lattice of mixed ferro-and antiferromagnetic bonds. We have analysed its thermodynamic properties in zero applied magnetic field, and we have also characterized the Kasteleyn transition that it exhibits when a magnetic field is applied. Finally, we have shown that its entropy when the magnetization is non-zero is well captured by the string representation.
One appealing feature of this model is that, unlike full three-dimensional spin ices, the Ising quantization axis is the same on each lattice site. This makes it natural to consider adding to the model a spatially uniform transverse magnetic field, Γ. The results of this should be particularly interesting in the h, Γ, k B T ≪ J regime, where the applied field is expected to stabilize the string phase at low temperatures, leading to a line of quantum Kasteleyn transitions in the zero-temperature (h, Γ) plane. This extension of the model (1) is the subject of a forthcoming work 33 .