## Abstract

We present a numerical study of magnetic ordering in spin ice on kagome, a two-dimensional lattice of corner-sharing triangles. The magnet has six ground states and the ordering occurs in two stages, as one might expect for a six-state clock model. In spin ice with short-range interactions up to second neighbours, there is an intermediate critical phase separated from the paramagnetic and ordered phases by Kosterlitz–Thouless (KT) transitions. In dipolar spin ice, the intermediate phase has long-range order of staggered magnetic charges. The high- and low-temperature phase transitions are of the Ising and 3-state Potts universality classes, respectively. Freeze-out of defects in the charge order produce a very large spin correlation length in the intermediate phase. As a result of that, the lower-temperature transition appears to be of the KT type.

## 1. Introduction

Spin ice is a ferromagnet with peculiar properties [1]. Its discovery in the rare-earth pyrochlore Ho_{2}Ti_{2}O_{7} by Harris *et al.* [2] attracted the interest of physicists because geometric frustration results in an enormous number of ground states in spin ice. The degeneracy is similar to that of proton positions in water ice [3] and manifests itself in a large residual entropy at low temperatures measured by Ramirez *et al.* [4]. Later, it was shown that spin ice has another peculiar feature: even though it remains disordered as the temperature goes to zero, spin correlations decay with the distance algebraically rather than exponentially [5,6]. Such critical behaviour is characteristic of systems with constraints exemplified by lattice models with hardcore dimers. In the case of spin ice, at low temperatures, each tetrahedron of rare-earth ions is forced to have two spins pointing towards the centre of the tetrahedron and two away from it.

The current resurgence of interest in spin ice stems from an unusual character of magnetic excitations in this class of materials [7,8]. Reversing a single spin in a spin-ice state violates the above-mentioned constraint on the two tetrahedra sharing that spin. One of them has three spins pointing in and one pointing out, the other one spin pointing in and three out. By reversing additional spins, the two defects can be moved around independently from each other. These two defect tetrahedra carry equal and opposite magnetic charges defined as sources and sinks of magnetic field **H**. A defect with magnetic charge *Q* in an external magnetic field **H**_{ext} experiences a Zeeman force *μ*_{0}*Q***H**_{ext}. Two defects with magnetic charges *Q*_{1} and *Q*_{2} interact with one another via a Coulomb potential *μ*_{0}*Q*_{1}*Q*_{2}/(4*πr*). Various properties of spin ice can be most naturally described by focusing on the motion of these defects [9–12].

The concept of magnetic charge is not exactly novel. It can be found in earlier research articles [13] and even in textbooks [14]. Magnetic charges in spin ice are remarkable because they are mobile and represent a rare example of fractionalized excitations in three spatial dimensions: the underlying degrees of freedom, magnetic dipoles, must be split in half, so to speak, to create magnetic monopoles.

It is worth pointing out that magnetic charges in spin ice are fundamentally different from the magnetic monopoles introduced by Dirac [15] to explain quantization of electric charge. They are sources and sinks of magnetic field **H**, rather than of magnetic induction **B**. For that reason, there are no restrictions on the possible values of magnetic charges in spin ice.

In this paper, we illustrate the utility of the concept of magnetic charge for the problem of spin ice on a different lattice. A natural generalization of the pyrochlore network is a lattice consisting of corner-sharing simplexes [16]. In two dimensions, corner-sharing triangles form kagome, figure 1. A spin is shared by two triangles and points along the line connecting their centres. The easy axes of the three spins on a triangle make angles of 120^{°} with each other. Like on the pyrochlore lattice, interactions of the antiferromagnetic kind are not frustrating: an isolated triangle has two ground states with all spins pointing in or all pointing out. A kagome lattice with nearest-neighbour antiferromagnetic interactions also has two ground states, in which spins point into triangles of the same orientation. On the other hand, ferromagnetic interactions are frustrating and yield six ground states on a single triangle, with two spins pointing in (and one out) or two pointing out (and one in). The number of ground states grows exponentially with the number of spins when interactions are ferromagnetic and only include nearest neighbours. In that, spin ice on kagome is quite similar to its analogue on the pyrochlore lattice. However, the addition of longer-range interactions, especially dipolar ones, reveals some important differences between the two systems.

The main difference between the pyrochlore network and kagome concerns the allowed values of magnetic charge *Q* of a simplex defined as the flux of magnetization **M** into the simplex. For a tetrahedron, *Q* is even (in the appropriate units): zero in a two-in-two-out state, ±2 in three-in-one-out or three-out-one-in states and ±4 in the all-in or all-out states. The energy of spin ice penalizes states with a large (in absolute value) charge *Q* so that simplexes are neutral, *Q*=0, in low-energy states. The lack of magnetic charges on tetrahedra explains a puzzling observation by Gingras & den Hertog: spin-ice states remain very nearly degenerate, even in the presence of long-range dipolar interactions [17]. As Castelnovo *et al.* [8] showed, the energy of dipolar interactions is well approximated by the Coulomb energy of magnetic charges residing on simplexes. Because these charges vanish in two-in-two-out states, their Coulomb energy is the same. The residual interactions, neglected in the Coulomb approximation, force the system into a phase with long-range magnetic order at temperatures low compared with the strength of dipolar interactions [18].

By contrast, the magnetic charge of a triangle is ±1 in a spin-ice state and ±3 when its spins point all in or all out. Even if the interactions tend to suppress magnetic charges, the lowest (in magnitude) values of charge on a triangle are ±1. This has important consequences that we briefly outline below and consider in detail in the rest of the paper.

First, the presence of a charge degree of freedom on simplexes makes spin ice on kagome much more fluid than its pyrochlore counterpart. In pyrochlore spin ice, spin dynamics slows down considerably at low temperatures. Single-spin flips result in the creation of two magnetic charges and thus become forbidden when the temperature falls below the energy cost of a monopole. Moves within the low-energy sector require flipping entire chains of alternating spins with a minimal length of six [18]. Alternatively, spin fluctuations can be mediated by the motion of a few remaining magnetic monopoles. Experimental studies confirm a very slow relaxation in spin ice at liquid-helium temperatures [19,20], although the primary mechanism of low-*T* spin dynamics in the rare-earth pyrochlores remains to be clarified. No such impediment to spin fluctuations exists in kagome ice because a single-spin flip does not necessarily take the magnet out of the low-energy sector. Such a local move is allowed if the magnetic charges of the two simplexes sharing the spin change from +1 and −1 to −1 and +1. Although at present, there are no materials realizing two-dimensional kagome spin ice, several experimental groups have made artificial magnetic arrays with this geometry [21]. Qi *et al.* [22] demonstrated that their artificial spin ice on kagome stayed strictly within the low-energy sector in which the magnetic charges remain ±1. Ladak *et al.* [23] observed triple charges, but this is likely the result of strong quenched disorder in their samples [24].

Second, the presence of a charge degree of freedom opens a possibility of new kinds of ordered phases in spin ice. Pyrochlore spin ice is expected to have only two thermodynamics phases. As the temperature is lowered from , the paramagnetic state gradually crosses over to the spin-ice state, which still preserves all of the symmetries of the Hamiltonian and thus remains a paramagnet, albeit a correlated one, from the standpoint of the Landau theory. At a sufficiently low temperature, the residual spin interactions not captured by the Coulomb model [8] induce a phase transition into an ordered state that breaks the time-reversal and lattice symmetries [25]. On kagome ice, one has good reasons to expect an intermediate ordered phase with staggered magnetic charges that minimizes the Coulomb energy [26,27]. This phase is characterized by very strong spin fluctuations. At the lowest temperatures, the symmetry is broken further, and the system settles into a state with long-range magnetic order.

To test these heuristic considerations, we have performed extensive numerical studies of spin ice on kagome, both with short-range interactions (nearest and next-nearest neighbours) and with long-range dipolar interactions. The problem turned out to be subtle and the answers were found to be different for the models with short- and long-range interactions. In both cases, there is an intermediate phase. However, the nature of the intermediate phase is drastically different. The dipolar spin ice has the expected charge-ordered phase, whereas the short-range model has a critical intermediate phase without true long-range order. The determination of the universality classes of the phase transitions required rather large system sizes because some of the transitions were of the Kosterlitz–Thouless (KT) type. Simulating large systems is particularly difficult in the presence of long-range interactions.

## 2. Kagome spin ice: the model

### (a) The Hamiltonians

The model of spin ice on kagome was first introduced by Wills *et al.* [28]. It has Ising spins living on the vertices of corner-sharing triangles of a kagome lattice. The unit vector points along the line connecting the centres of triangles, while the Ising variable *σ*_{i} encodes the state of the spin. It is convenient to choose the unit vectors in such a way that they point into triangles of one orientation (say, △). The Hamiltonian of the model with nearest-neighbour (nn) exchange *J*_{1} can be written in two ways,
2.1
where we relied on the result for nearest neighbours. As usual, a ferromagnetic exchange *J*_{1}>0 for spins **S**_{i} translates into an antiferro-magnetic interaction for the Ising variables *σ*_{i}. It is well known that the Ising antiferromagnet on kagome remains paramagnetic down to zero temperature [29]. It retains a finite entropy density in the ground state, 0.502 per spin [30].

To understand this residual spin degeneracy, we rewrite *H*_{1} in terms of magnetic charges of triangles
2.2
Here, we have neglected an irrelevant constant. The magnetic charge is defined as
2.3
with the plus sign for one type of triangles and the minus for the other. As discussed earlier, the allowed values of magnetic charges are *Q*_{α}=±3 and ±1. Equation (2.2) shows that triangles with triple charges are energetically disfavoured. Consequently, as temperature is lowered, the magnet gradually enters the spin-ice phase consisting exclusively of triangles with *Q*_{α}=±1, corresponding to the two-in-one-out and one-in-two-out ice rules on kagome. The number of the spin-ice microstates grows exponentially with the number of spins. Invoking the famous Pauling estimation gives an entropy density per spin, which agrees very well with the numerical results.

Wills *et al.* [28] added interactions between next-nearest neighbours (nnn),
2.4
They considered both signs of *J*_{2}, but we are interested in the antiferromagnetic case *J*_{2}<0. Möller & Moessner [26] and Chern *et al.* [27] added long-range dipolar interactions,
2.5
where is the characteristic strength of dipolar coupling, **r**_{i} are spin locations, , *μ* is the magnetic dipolar moment of a spin and *r*_{nn} is the distance between nearest neighbours. For brevity, we shall refer to the model of Wills *et al.* [28], with the Hamiltonian *H*_{1}+*H*_{2} and antiferromagnetic *J*_{2}<0, as the *short-range kagome ice* and call the model with dipolar interactions, with the Hamiltonian *H*_{1}+*H*_{d}, the *long-range kagome ice*.

### (b) The ground states

In both models, interactions between Ising variables *σ*_{i} are antiferromagnetic for nearest neighbours and ferromagnetic for next-nearest neighbours. On the basis of that, one might expect the same type of magnetic order at low temperatures. The ground states for the Ising model on kagome with first- and second-neighbour interactions are known exactly [31,32]. They have an enlarged () magnetic unit cell with nine spins. The corresponding ground states of the ice model are shown in figure 2*a*. With the aid of a complex-order parameter,
2.6
where **Q**=(4*π*/3,0) and *N* is the number of spins, we can see that they are similar to those of the six-state clock model. The phase of the order parameter *M*=|*M*|e^{iϕ} takes on the six values *ϕ*=*nπ*/3, where *n* is an integer (figure 2*a*).

We were able to show that the long-range model with dipolar interactions has the same ground states. To do so, we treated the Hamiltonian *H*_{1}+*H*_{d}=*σ*_{i}*A*_{ij}*σ*_{j}/2 as a quadratic form in the Ising variables. We replaced the constraint *σ*_{i}=±1 with a less stringent normalization condition, *σ*_{i}*σ*_{i}=*N*, where *N* is the number of sites, and minimized the quadratic form by finding the lowest (most negative) eigenvalue of the matrix *A*_{ij}. Although this procedure minimizes the energy over an enlarged space of states, the method gives the right answer if the eigenstate with the lowest eigenvalue has *σ*_{i}=±1. That is indeed the case for the Hamiltonian *H*_{1}+*H*_{d}. The corresponding eigenstates are the six ground states of the Ising model on kagome.

### (c) Phase transitions: general considerations

Previous theoretical studies have revealed that the six-state clock model in two dimensions may order in a number of scenarios [33]. In all of them, there is a high-temperature disordered phase with 〈*M*〉=0 and a low-temperature ordered phase with 〈*M*〉≠0.

(i) Two KT transitions with an intermediate critical phase in which 〈

*M*〉=0 and spin correlations decay as a power of the inverse distance.(ii) From the paramagnetic phase, the system enters a partially ordered phase via an Ising transition. A second transition of the 3-state Potts universality class takes the system to the ordered state. The intermediate phase has an Ising-order parameter 〈

*M*^{3}〉≠0, whereas 〈*M*〉=0.(iii) Similar to (ii) but with the Ising and Potts transitions exchanged. The intermediate phase has a 3-state Potts-order parameter 〈

*M*^{2}〉≠0, whereas 〈*M*〉=0.(iv) A discontinuous phase transition between the paramagnetic and fully ordered phase.

Numerical studies [31,32] provide compelling evidence for scenario (i) with two KT transitions in short-range Ising models on kagome. Takagi & Mekata [31] mistakenly ascribed a non-zero-order parameter 〈*M*〉 to the intermediate phase, which is expected to have quasi-long-range order with power-law spin correlations and 〈*M*〉=0.

Wills *et al.* [28] also suggested an intermediate ordered phase with six possible states shown in figure 2*b*. They are equivalent to the intermediate states of Takagi & Mekata and have a non-zero-order parameter 〈*M*〉 with *ϕ*=(*n*+1/2)*π*/3. Such an intermediate state is clearly inconsistent with any of the four scenarios mentioned earlier.

The intermediate critical phase in scenario (i) is most simply understood in the continuous version of the clock model, namely the XY model with a sixfold anisotropy term in the free-energy functional. This term is irrelevant in the renormalization-group sense above the lower phase transition, so that the intermediate phase behaves as if the sixfold anisotropy were absent. In principle, there is a possibility that the sixfold anisotropy changes sign: at low temperatures, *a*_{6}<0 yields the states shown in figure 2*a*, whereas at higher temperatures, *a*_{6}>0 stabilizes the states shown in figure 2*b*. However, if the intermediate phase is indeed ordered and fully breaks the *Z*_{6} symmetry, then the transition between that phase and the fully symmetric paramagnet should follow one of the above scenarios, which requires yet another intermediate phase or a discontinuous phase transition. We see no evidence for either possibility.

Our proposal of the intermediate state with ordered magnetic charges corresponds to scenario (ii). The charge-ordering transition is of the Ising type. It takes the system from the paramagnetic phase into a state where the Ising-order parameter 〈*M*^{3}〉 has a non-zero expectation value. If 〈*M*^{3}〉=〈|*M*|^{3}e^{3iϕ}〉>0, the state can be thought of as a mixture of the clock states with *ϕ*=0, 2*π*/3 and −2*π*/3 (figure 2*a*). At a lower temperature, the system selects one of these states in a transition of the 3-state Potts universality class.

To test which of the scenarios are realized in kagome spin ice, we have performed large-scale numerical simulations of both the short- and long-range models with the Hamiltonians: *H*_{1}+*H*_{2} and *H*_{1}+*H*_{d}.

## 3. Numerical simulations

### (a) Kagome ice with short-range interactions

Here we provide numerical evidence for scenario (i), namely two successive KT transitions, in the short-range ice model on kagome. A standard Metropolis importance sampling was used to simulate the behaviour of the model on *L*×*L* unit cells with periodic boundary conditions; the number of spins *N*=3*L*^{2}. In a single Monte Carlo step, each spin in the lattice is updated once with a probability , where *ΔE* is the energy change of flipping a spin. It is worth noting that the single-spin update is sufficient for both the spin-ice regime and even the intermediate critical phase; the corresponding acceptance rates are roughly 30 per cent and 15 per cent, respectively. About 10 000 initial Monte Carlo sweeps were discarded to allow the system to equilibrate. We then retained results from 10^{5} to 10^{6} Monte Carlo sweeps for computing averages. Our most complete dataset was taken with parameters *J*_{2}=−*J*_{1}/3, for which we describe our analysis in detail in the following.

Figure 3 shows the temperature variation of specific heat, staggered charge order and magnetic-order parameters for different lattice sizes. While the appearance of two peaks in the specific-heat curve suggests two phase transitions, the locations and amplitudes of the peaks vary only slightly with lattice size, a feature characteristic of KT transitions. The intermediate critical phase, determined by finite-size scaling to be discussed below, is marked by the shaded region in the figures. The charge-order parameter is defined as the difference of magnetic charges on the two different types of triangles,
3.1
where (−1)^{α}=±1 for up and down triangles, respectively, and 2*L*^{2} is the number of triangles in the lattice. The charge order is proportional to the Ising-order parameter, 〈*Q*〉∝〈*M*^{3}〉. As can be seen from figure 3*b*, the staggered charge order decreases rapidly with increasing *L* above the low-*T* transition, indicating that magnetic charges remain disordered in the intermediate phase.

On the other hand, the magnetic-order parameter shows quite different finite-size behaviour across the three different regimes (figure 3*c*). Above the high-*T* peak, the rapid decrease of *M* with increasing *L* is consistent with a magnetically disordered phase, whereas a negligible finite-size effect implies an ordered state below the lower transition at *T*_{c1}. In the intermediate regime, the order parameter *M* falls off rather slowly with increasing system size. These observations are summarized in figure 4, which shows magnetic-order parameter *M* as a function of *L* at different temperatures. The magnetic-order parameter extrapolates to zero at high temperatures, whereas it levels out to a constant at low *T*. In the intermediate regime, the linear behaviour in the double-logarithmic plot indicates a power-law dependence
3.2
which is consistent with an infinite correlation length in a critical phase. The extracted value of critical exponent *η* as a function of temperature is shown in the inset of figure 4.

As discussed in §2*c*, the sixfold anisotropy in the Landau expansion of the free-energy functional is irrelevant in the intermediate critical phase. To confirm this, we show in figure 5 the distribution of instantaneous-order parameter *M*=|*M*| e^{iϕ} in the complex (Re*M*,Im*M*) plane at various temperatures. As can be seen in figure 5*c*, which corresponds to a temperature in the intermediate regime, the order parameter has a finite amplitude |*M*|≠0 owing to a finite system size (*L*=180). More importantly, a continuous O(2) symmetry emerges in this critical phase. In the ordered phase below *T*_{c1}, the anisotropy term reasserts itself and restores the *Z*_{6} symmetry, as demonstrated in figure 5*a,b*.

To further corroborate the proposed scenario of two KT transitions, we resort to the method of finite-size scaling [34]. For a KT transition, the correlation length near the critical temperature *T*_{c} diverges as [35]
3.3
where *a* is a non-universal constant and *t*=|*T*−*T*_{c}|/*T*_{c} is the reduced temperature. The order parameter and its susceptibility exhibit power-law behaviour, *M*∝*ξ*^{−η/2} and *χ*_{M}∝*ξ*^{2−η}, in a KT transition. For finite systems, the singular part of the free energy is expected to depend only on *ξ*/*L*. The order parameter and susceptibility are assumed to have the functional forms
3.4
where and are unknown universal functions. At the critical point, substituting into equation (3.4) gives the power-law relation (3.2). For an extended regime consisting of a line of critical points, the power law *M*∼*L*^{−η/2} should hold over the entire temperature range from *T*_{c1} to *T*_{c2}, which is indeed observed in our simulation, figure 4.

The finite-size scaling relations (3.4) can also be used to determine the upper and lower critical temperatures of the intermediate KT phase. For example, with an appropriately chosen set of parameters *a*, *b* and *T*_{c1}, the plot of *ML*^{b} versus should collapse on a universal curve for different system sizes. The same should hold for *χ*_{M}*L*^{−c} versus . The two constants, *b* and *c*, are related to the critical exponent: *b*=*η*/2 and *c*=2−*η*. As shown in figure 6, excellent data collapse was obtained using the following set of parameters: *a*=1.22, *b*=0.0615, *c*=1.746, *T*_{c1}=0.735*J*_{1} and *T*_{c2}=0.845*J*_{1}. From this, we estimated the critical exponent at the two transition temperatures: *η*(*T*_{c1})=2*b*=0.123 and *η*(*T*_{c2})=2−*c*=0.255. These values agree reasonably well with those extracted by linear fitting of the power-law relation (3.2) within the critical phase (see the inset of figure 4). Both values are close to the theoretical predictions and for the six-state clock model [36], and the slight overestimation could be due to the finite-size effects on KT transitions.

### (b) Kagome ice with dipolar interactions

Although the single-spin Metropolis algorithm is quite efficient for simulating short-range kagome ice, it suffers from a dynamical freezing in the charge-ordered phase. A similar problem arises in the low-temperature simulation of pyrochlore spin ice [18] in which a single-spin flip violates the ice rules, leading to a large energy cost and low acceptance rate of the updates. To overcome this problem, we added the non-local loop moves first introduced by Barkema & Newman [37] for square ice models. In a non-local update, a loop is first formed by randomly tracing a path through triangles satisfying the ice rules, alternating between spins pointing in and spins pointing out of the triangles. The process is completed when the path closes upon the starting spin or encounters any spin already included in the loop. Flipping all the spins in such a loop leaves the magnetic charges {*Q*_{α}} intact and conserves the nearest-neighbour energy *H*_{1}. The loop move results in a small gain or loss of the dipolar energy *ΔE*_{d}, and the update is accepted with the probability .

We used a combination of single-spin flips and loop moves to simulate the long-range ice model *H*_{1}+*H*_{d}. With periodic boundary conditions, dipolar interactions were summed over periodic copies up to a distance of 500*L*. Most of the results presented here are obtained with a dipolar coupling *D*=2*J*_{1}. Figure 7 shows the temperature dependence of specific-heat-, charge- and magnetic-order parameters for various system sizes. Similar to the case of short-range kagome ice, two peaks can be seen in the specific-heat curves, indicating two phase transitions. While the high-temperature peak becomes sharper with increasing *L*, the low-*T* transition barely shows any finite-size dependence. The behaviours of charge- and magnetic-order parameters shown in figure 7*b*,*c* are consistent with scenario (ii) discussed in §2. We discuss both transitions below.

The staggered charge order corresponds to the partially ordered phase in scenario (ii) that is characterized by an Ising-order parameter 〈*Q*〉∝〈*M*^{3}〉. The charge-ordering transition is thus expected to be in the universality class of the Ising model. To verify this conjecture, we performed finite-size scaling analysis and found excellent data collapse, figure 8, using the critical exponents of the two-dimensional Ising universality class. The critical temperature *T*_{c2}=0.29*D* is determined from the crossing of Binder’s fourth-order cumulant for different lattice sizes.

The origin of an intermediate phase with ordered charges can be understood by expressing the energy of the dipolar ice in terms of magnetic charges. This is achieved through the so-called dumbbell approximations [8], in which the dipoles are stretched into bar magnets of length such that their poles meet at the centres of triangles,
3.5
Here, *q*_{m}=*μ*/*a* is the magnetic charge of the dumbbell, *μ* is the moment of the spin. The self-energy for kagome ice. It can be shown that the dipolar energy *H*_{1}+*H*_{d} is well approximated by equation (3.5), with corrections that vanish with distance at least as fast as 1/*r*^{5} for each dipole pair.

The dumbbell approximation also illustrates an important difference between pyrochlore and kagome ices. Because the allowed charges are even on a tetrahedron and odd on a triangle, minimization of the self-energy term results in *Q*_{α}=0 on the pyrochlore lattice and ±1 on kagome. The conditions of minimal allowed charges on simplexes correspond to the ice rules on the respective lattices. For pyrochlore ice, the vanishing of magnetic charge also makes the Coulomb interaction, the second term in equation (3.5), strictly zero. This observation explains why all states satisfying the ice rules are essentially degenerate in pyrochlore spin ice over a wide range of temperatures [17]. The degeneracy is lifted by the residual interactions neglected in the Coulomb approximation and a long-range magnetic order with wavevector **Q**=(0,0,2*π*) sets in at a lower temperature compared with the dipolar energy scale [18].

The situation on kagome is quite different as non-zero charges generate a magnetic field. This results in substantial energy differences between states obeying the ice rule *Q*_{α}=±1. Consequently, interactions between uncompensated charges induce an Ising transition into a state with staggered charge order that minimizes the Coulomb energy. This partially ordered phase is closely related to the ice states of pyrochlore spin ice in a 〈111〉 magnetic field [38,39]. Triangles with positive (negative) charges on kagome correspond to three-in-one-out (three-out-one-in) tetrahedra of the pyrochlore lattice.

Spins remain disordered in the charge-ordered ice phase, as evidenced by the loop moves discussed at the beginning of this section. Such non-local loop updates flip spins on a closed path while maintaining the charge configuration. To quantify this residual degeneracy, we note that each triangle in a state with perfect charge order has two majority spins pointing into (or out of) the triangle and a minority spin pointing the other way. Such states can be mapped to dimer coverings on the honeycomb lattice by identifying the minority spins with the dimers. The number of dimer coverings on a honeycomb grows exponentially with the lattice size, giving rise to a residual entropy density *S*=0.108 per spin [38].

The remaining entropy of the charge-ordered phase is completely removed by a magnetic phase transition corresponding to the low-*T* peak of the specific-heat curve (figure 7*a*). The magnetic order is shown in figure 1*d* and is again characterized by order parameter *M*; it has an enlarged unit cell containing nine spins. Similar to pyrochlore spin ice, the magnetic ordering is induced by residual interactions beyond the dumbbell approximation equation (3.5).

As the *Z*_{6} symmetry of the spin-ice Hamiltonian is reduced to a threefold symmetry in the charge-ordered phase, the magnetic transition is expected to be in the universality class of 3-state Potts model. However, Monte Carlo simulations on systems up to *L*=36 fail to turn up any signature of the Potts criticality. Instead, the lack of a singularity in the specific heat, figure 9*a*, seems to be consistent with a KT transition.

### (c) Dimers on a honeycomb lattice

To shed light on this puzzling result, we consider a similar transition in the honeycomb dimer model. As discussed above, the charge-ordered states can be uniquely mapped to dimer coverings on the honeycomb. To induce the corresponding ordering of dimers, we introduce a next-nearest-neighbour attractive interaction between the dimers. Note that nearest-neighbour dimers are precluded by the hardcore constraints that each lattice site is attached to exactly one dimer. The partition function of the modified dimer model reads
3.6
where the sum is over all dimer coverings of the hexagonal lattice and counts the number of next-nearest-neighbour pairs in a given covering . An attractive interaction between next-nearest-neighbour dimers corresponds to *v*_{2}>0.

We used the direct-loop Monte Carlo algorithm [40] to simulate the dimer model (3.6). The non-local update is performed by initially breaking up an arbitrary dimer into a pair of monomers and then moving one monomer across the lattice by flipping a sequence of dimers along a path. The process is completed, and a new dimer configuration is created once the two monomers meet and recombine with each other. Because detailed balance is maintained at each step of the monomer’s movement, a large number of dimers can be efficiently updated without rejection. This method significantly reduces the autocorrelation time in the Monte Carlo process and allows us to simulate large lattices up to *L*=150. Figure 9*b* shows the specific heat versus temperature for various system sizes; also shown in the inset is the temperature dependence of the order parameter characterizing the dimer order. The rather weak finite-size dependence and a lack of singularity in specific heat again shows that the dimer ordering is characterized by a KT transition similar to the case of the dimer model with aligning interactions on the square lattice [41].

The appearance of a KT transition indicates the critical nature of the disordered hardcore dimer phase, or equivalently the charge-ordered ice phase. Indeed, by further mapping the dimer covering to a ‘height’ field [42], the dimer model is described by a sine–Gordon model in the coarse-grained approximation. It is well known that the height model undergoes a KT transition into a rough phase at high temperatures [43]. As the confining cosine potential term is irrelevant at high *T*, the rough phase is described by a Gaussian field theory with critical correlation functions in two dimensions.

The existence of such a critical phase relies on the absence of charge defects. At finite temperatures, however, thermally excited defects spoil the mapping to hardcore dimer coverings. There are two types of defects in the charge-ordered phase: triangles with ±3 charges (ice defects) and triangles that violate the charge order (charge defects). In the dimer model, the triply charged sites become monomers, whereas defects in charge order become sites with two dimers, figure 10. We ignore the triply charged defects in the intermediate phase on account of their high energy cost and focus on the charge defects. To that end, we modified the honeycomb model (3.6) by allowing dimer pairs with fugacity . The hardcore dimer model is recovered in the limit . Figure 11 shows the temperature variation of the order parameter, susceptibility and density of dimer pairs for *ε*=2, 10 and in units of *v*_{2}. The model with larger fugacity for dimer pairs shows a dramatically different behaviour from that of the hardcore dimer covering (figure 11).

For , a finite density of dimer pairs sets a bound on the dimer correlation length, thus altering the criticality of magnetic ordering to the universality class of 3-state Potts model. To confirm this, we conducted a finite-size scaling study on the dimer model with *ε*=2*v*_{2}. As can be seen from figure 12, one obtains excellent data collapse with the critical exponents of the two-dimensional 3-state Potts model. However, when the average separation between charge defects exceeds the lattice size, we are back to hardcore dimers with power-law spatial correlations for all distances. This explains the KT-like behaviour observed in spin ice at small lattice sizes. The critical behaviour characteristic of the 3-state Potts universality only reveals itself for sufficiently large systems [44].

## 4. Discussion

Magnetic charges in spin ice are an example of an emergent phenomenon. Although the fundamental degrees of freedom in these materials are magnetic dipoles, the behaviour of the system at low energies is often more conveniently expressed in the language of magnetic charges. In spin ice on kagome, magnetic charge of a simplex can take on odd values ±1 and ±3 in natural units. Even though magnetostatic interactions tend to minimize magnetic charge, low-energy spin-ice states have non-vanishing magnetic charges on simplexes. Such a magnet may have, in addition to the paramagnetic and fully ordered states, a distinct intermediate phase with ordered magnetic charges and disordered spins. This phase possesses considerable residual entropy (0.108 per spin) because a given pattern of magnetic charges can be realized by many configurations of magnetic dipoles.

Our numerical simulations provide evidence for an intermediate phase with staggered magnetic charges in spin ice with dipolar interactions on kagome. As the system is cooled down from the paramagnetic state, it first gradually enters the spin-ice regime in which magnetic charges of triangles are restricted to the smallest (in magnitude) values of ±1. A phase transition of the Ising universality class takes the magnet into the intermediate phase with staggered magnetic charges.

At a lower temperature, the system enters a magnetically ordered phase with six ground states shown in figure 2*a*. Because only three of them are accessible from each of the two states of the charge-ordered phase, this transition is expected to be in the universality class of the 3-state Potts ferromagnet. However, observing the corresponding critical behaviour proved challenging because the transition takes place in the background of nearly-perfect magnetic charge order. When the charge order has no defects at all, magnetic configurations can be mapped onto states of hardcore dimers with attraction. This mapping establishes that spins in the background of perfectly ordered magnetic charges are in a quasi-ordered phase with algebraic spatial correlations. The transition to the fully ordered phase is then of the KT type. Thermal fluctuations generate defects in charge order, thus violating the hardcore constraint for dimers. The spin correlation length is then set by the typical distance between isolated charge defects. As figure 7*b* shows, the charge-order parameter quickly approaches saturation upon cooling below *T*_{c2}=0.29*D*. At *T*=0.24*D*, i.e. considerably above the onset of magnetic order at *T*_{c1}, the concentration of charge defects is about 0.01 per triangle. Furthermore, most of these defects come in pairs with zero net charge and minimal separation (resulting from single-spin flips); so they just renormalize the spin correlation function without shortening the correlation length. Not surprisingly, the intermediate phase exhibits power-law spin correlations up to very large distances and the transition to the magnetically ordered phase appears to be of the KT type. Observing the true critical behaviour of the 3-state Potts universality class would require extremely large system sizes inaccessible to us.

## Acknowledgements

The authors thank Paula Mellado, Gunnar Möller and Roderich Moessner for useful discussions. G.W.C. was supported by the US National Science Foundation grant no. DMR-0844115. O.T. was supported in part by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under award no. DE-FG02-08ER46544.

## Footnotes

One contribution of 8 to a Theo Murphy Meeting Issue ‘Emergent magnetic monopoles in frustrated magnetic systems’.

- This journal is © 2012 The Royal Society