Ising Fracton Spin Liquid on the Honeycomb Lattice

We study a classical Ising model on the honeycomb lattice with local two-body interactions and present strong evidence that at low temperature it realizes a higher-rank Coulomb liquid with fracton excitations. We show that the excitations are (type-I) fractons, appearing at the corners of membranes of spin flips. Because of the three-fold rotational symmetry of the honeycomb lattice, these membranes can be locally combined such that no excitations are created, giving rise to a set of ground states described as a liquid of membranes. We devise a cluster Monte-Carlo algorithm purposefully designed for this problem that moves pairs of defects, and use it to study the finite-temperature behavior of the model. We show evidence for a first order transition from a high-temperature paramagnet to a low-temperature phase whose correlations precisely match those predicted for a higher-rank Coulomb phase.

We study a classical Ising model on the honeycomb lattice with local two-body interactions and present strong evidence that at low temperature it realizes a higher-rank Coulomb liquid with fracton excitations.We show that the excitations are (type-I) fractons, appearing at the corners of membranes of spin flips.Because of the three-fold rotational symmetry of the honeycomb lattice, these membranes can be locally combined such that no excitations are created, giving rise to a set of ground states described as a liquid of membranes.We devise a cluster Monte-Carlo algorithm purposefully designed for this problem that moves pairs of defects, and use it to study the finite-temperature behavior of the model.We show evidence for a first order transition from a high-temperature paramagnet to a low-temperature phase whose correlations precisely match those predicted for a higher-rank Coulomb phase.
One of the central concepts in condensed matter physics is the notion of quasiparticles: the idea that low energy excitations of a system are weakly interacting particle-like objects.In general, these quasipaticles are capable of independent motion, and it is via this motion that energy inserted locally into the system can spread out, thus allowing equilibration.Fractons are quasiparticles outside this paradigm, being completely immobile when isolated [1][2][3][4][5][6][7][8][9][10][11].Fractons are intimately connected with the conservation laws of exotic gauge theories involving not only charge but higher moments (e.g.dipole moment) of the charge density [12][13][14][15][16].It is these conservation laws which render isolated fractons immobile.
Recent years have seen a concerted effort to establish theoretical models in which higher moment conservation laws and fracton physics appear [17][18][19][20][21][22][23][24][25].Various lattice models have been proposed to give rise to fractonic behavior, although these often require complicated multibody interactions.From the point of view of identifying routes to experimental realization, it is preferable to find models built-from short-ranged, two body interactions.
One setting in which the construction of such models has been successful is classical spin systems [26][27][28].Classical models have some advantages: a Hamiltonian can be readily be constructed to enforce a local constraint of choice in the low energy sector, and this constraint can be chosen in such a way as to reproduce the Gauss's law(s) of a given gauge theory.Once the model is constructed, it is in principle accessible to Monte Carlo simulation.
Thus far these models have been constructed from continuous degrees of freedom, but this has the drawback that one cannot isolate and study discrete fractons.Moreover, it is difficult to reintroduce quantum fluctuations in a controlled way.If fracton physics could be demonstrated in a classical Ising model, fractons will naturally be present as discrete excitations, and there exists a clear route to studying quantum effects by introducing off-diagonal terms, e.g.via a transverse field.
Here, we present an Ising model exhibiting a fractonic spin liquid regime.The model has an extensive degen- eracy of ground states, which may be thought of as the classical limit of a gapless fracton phase.This is distinct from fracton topological order, in which there is a subextensive degeneracy of locally indistinguishable quantum ground states.The excitations of the model are fractons appearing at the corners of membranes of flipped spins.Pairs of fractons can be bound into lineons, which can move along a certain lattice direction, but not along the perpendicular direction.We identify microscopic resonance processes which take the system between different ground states, consideration of which yields a lower bound on the observed entropy.In order to simulate the model's properties we devise -to our knowledge, the first -cluster algorithm for Monte Carlo simulations of frac-tons, which is purpose-built for the study of fractonic Ising models.This allows numerical access to relatively large systems and low temperatures.
Simulating the model, we find that it exhibits a first order phase transition at low temperatures, but that the low temperature state nevertheless lacks signs of conventional order.Instead, the ground state correlations exhibit four-fold pinch points in momentum space, which are characteristic of systems described by a gauge theory of tensor fields [29].
Fractons in the honeycomb model.We consider an Ising model on the honeycomb lattice with Hamiltonian which is a sum over constraints M h , defined on hexagons h and their exteriors ⟨h⟩ as illustrated in Fig. 1 (a).A Hamiltonian of this form was first considered for O(3) Heisenberg spins in Ref. 27.
There, it was shown that the system upon coarse graining can be described in terms of a suitably defined rank-2 tensor field m, subjected to a Gauss law with Tr[m] = 0.The structure of this Gauss's law implies conservation of the dipole moments of charge density, and hence fractonic excitations, as we discusss further below.
The field m relates to the microscopic degrees of freedom via a staggered, position dependent mapping.As a result, the relationship between the charge ρ and the microscopic constraint M h breaks lattice symmetry by hand.We choose a subset of hexagons such that each site is a member of exactly one of them.One possible choice is illustrated in Fig. 1 by a darker shade of some hexagons.Denoting this subset as + hexagons and the rest as −, the charge is then defined as Charges and moments of the higher-rank gauge theory can be determined explicitly from the microscopic model.The single spin flip, shown in Fig. 1 (b), preserves both total charge ρ and the dipole moments d ν = r ν ρ.The lowest order moment of the charge distribution which changes is the quadrupole moment q µν = r µ r ν ρ.
An excitation with the lowest (non-zero) energy of Eq. (1) involves four defects (hexagons with M h = ±1); it can be constructed from four spin flips, Fig. 1 (c).These defects are fractons: no single defect can be moved by any local combination of spin flips since that would change the total dipole moment.Pairs of (oppositely charged) fractons are lineons since they can be moved in the direction perpendicular to their dipole moment by flipping another four spins on the next-next hexagon, as indicated in Fig. 1 (c) in dashed-orange.Generally, fractons in our model appear at the corners of a "membrane" of flipped bonds.This, along with the ability to bind fractons into lineons implies this is a "Type-I" fracton model, constrasting with "Type-II" models where fractons appear on fractal structures without the possibility of mobile bound states [8].
While the mapping between the microscopic model and the coarse-grained field m µν in principle generalizes to the Ising model, it is an open question whether the system still realizes a higher-rank Coulomb phase or whether restricting the degrees of freedom to be discrete yields a set of ground states the average over which no longer corresponds to the deconfined phase of the gauge theory.For the case of the 'conventional' Coulomb liquid, cases are known where the hard-spin Ising and Heisenberg behaviors are (e.g.pyrochlore [30]), and are not (e.g.kagome [31][32][33][34]) , in accord with that of the soft-spin theory.In the following, we argue for the former case, that is that the Ising model [Eq.( 1)] does also realize a higher-rank Coulomb regime.To this end, we first show explicitly that the number of ground states grows exponentially with the number of sites N by identifying local "resonance processses" between different ground states.Second, using a novel cluster Monte-Carlo algorithm that moves pairs of defects we gain access to the thermodynamics of the model at large system size and low temperature.This allows us to both quantitatively extrapolate the residual entropy and compute the low-temperature correlations.
Extensive ground-state degeneracy.For periodic boundary conditions, the fact that pairs of fractons are lineons already implies a subextensive ground state entropy: we can create a pair of lineons, move one of them around the system in a nontrivial way and annihilate the pair again, reaching a different ground state.In addition to moving a lineon in a particular direction, we can also split it into two, as shown in Fig. 2  (a).Note that such a move is only possible because of sixfold rotation symmetry and would not be possible with cubic symmetry and more generally in the absence of at least three distinct orientations for the membranes.Crucially, in our case it can be done in two ways which we call forward and backwards split, shown on the left and right of Fig. 2 (c) respectively.The backwards split is particularly important since it allows us to close the worldline of lineons locally, resulting in a local move between different ground states.The minimal such move is shown in Fig. 2 (b) and is a combination of six membranes such that each corner defect is annihilated with exactly one other corner of opposite charge.This corresponds to flipping 24 spins simultaneously, as indicated in dashed-orange in panel (c) of the same figure.There, we also explicitly show a state with a finite density of such flippable motifs.It is constructed by starting form the Neél state (which has zero energy) and flipping spins along lines as indicated in gray, which also preserves the ground state constraint.The 24-spin move as shown in Fig. 2 (b,c) could be centered on any of the hexagons colored in green, and non-overlapping motifs can be flipped independently.This establishes a number of states exponential in the number of sites N and in particular implies a lower bound on the residual entropy Since this is a particular construction and there is a wealth of possibilities to combine membranes locally such that they create no defects, we do not expect this bound to be tight.A more quantitative estimate can be obtained from Monte-Carlo simulations (cf.Fig. 3) as discussed in the following.
Cluster Monte-Carlo algorithm.It is well established that in the presence of fracton excitations, any local algorithm will have a rapidly diverging relaxation time at low tempearture [1,2,5].Since single defects are immobile, local algorithms generally fail to anneal them out at low temperature and to gain access of the thermodynamics of Eq. ( 1), clearly a cluster algorithm is desirable.We have designed such an algorithm, which moves pairs of defects by effectively attempting to span one row of a membrane as shown in Fig. 2 (c).In the crucial step, the algorithm starts by flipping two bonds on a single hexagon and accepting this move with metropolis probability.In the case that the move is rejected, the algorithm attempts to flip another four spins [as indicated in dashed orange in Fig. 2 (c)], again trying to accept the FIG. 4. The structure factor in the low-temperature phase at the system sizes studied shows no sharp features except clearly visible four-fold pinch points.These pinch points are the signature of an emergent higher-rank Gauss law [29] and are of the same form as found for a higher-rank spin liquid in a related continuous spin model [27].
move with metropolis probability, multiplied by a factor to account for the probability of rejecting the first move.This additional factor is needed to ensure detailed balance [35].This is repeated until either the cluster is accepted or is ultimately rejected while spanning the full (linear) system size.If the move is accepted with zero energy cost, it either moves a distant pair of defects by one step in the direction perpendicular to the long side of the (single-row) membrane or it changes the ground state sector.While the relaxation time of this algorithm still scales significantly with system size at low temperature (we estimate τ ∼ L 7.4 ), it constitutes a major improvement over local dynamics.To demonstrate this, we performed a simulated annealing simulation using 100 intermediate temperatures between T = J and T = J/10 and 10 6 sweeps per temperature.When using a spin-flip metropolis algorithm at L = 18, out of 500 runs not a single one manages to anneal out all defects.In contrast, our cluster algorithm in the same configuration reaches zero energy in all cases.When augmented with the local 24 spin move shown in Fig. 2 (b-c), accepted with metropolis probability, and using feedback-optimized parallel tempering [36], we are able to equilibrate systems with up to N = 1152 (L = 24) spins in the ground state regime.To ensure equilibration, we compute the specific heat both directly from energy fluctuations and also as the derivative of internal energy with respect to temperature, and verify that these two estimators agree (cf.Fig. 3).This is considered a "stringent criterion" and used for example in computational studies of model glass formers [37].
Thermodynamic properties and low-temperature correlations.We use the setup described in the previous paragraph to study the thermodynamic properties of Eq.
(1) as a function of temperature for a range of system sizes as shown in Fig. 3.Both internal energy and specific heat are compatible with a first-order transition from a high-temperature paramagnetic phase to a lowtemperature phase corresponding to the ground state regime ⟨E⟩ ≈ 0. Integrating the specific heat from high temperature yields an estimate of the residual entropy per site, S 0 /N for each system size.In the figure, we show the result of integrating both the variance of the internal energy (blue triangles) as well as its derivative with respect to temperature (orange squares).Fitting the system size dependence using both a linear fit to the last four sizes and a quadratic fit to all system sizes yields an estimate of S 0 /N = 0.037 ± 0.010.This is quite a bit above the lower bound derived above [Eq.( 4)], but also comes with a large error bar due to significant finite size effects.The presence of these even at relatively large system sizes is not surprising given the size of the smallest local move as shown in Fig. 2 (b,c).As a rough quantification of the importance of finite-size effects, we note that taking into account periodic boundary conditions, the largest system size L = 24 has ⌊L/2⌋ = 12 independent degrees of freedom on the boundary, compared to ⌊L/6⌋ 2 ≃ 16 in the bulk.
Finally, we show the structure factor within the lowtemperature phase in Fig. 4. It is fully consistent with a low-temperature higher-rank Coulomb phase and shows four-fold pinch-points at the zone boundaries [29], which are sharp (that is exactly one pixel) for all system sizes considered.
The presence of a sharp transition as a function of temperature is somewhat surprising since in two dimensions, the higher-rank Coulomb liquid expected to describe the ground state regime is continuously connected to the paramagnet.It is however not inconsistent with a low-temperature liquid phase since a first-order transition is always possible also between continuously connected phases as demonstrated by the famous transition between gaseous and liquid water.An alternative possibility is that the low-temperature phase is a so called "fragmented liquid" [38], that is the set of ground states, although extensive, breaks some symmetry on average.A possible hint in this direction is that the maximum of the structure factor along the line cut shown in Fig. 4 scales roughly with linear system size.However, as discussed already above there are still significant finite size effects.
Ultimately, it is impossible to exclude this possibility of fragmentation in the absence of a more efficient algorithm and we leave this question open for future studies.
Conclusion.In summary, we have demonstrated the appearance of a fractonic spin liquid in the low temperature state of an Ising model on the honeycomb lattice.This low temperature state is separated from the high temperature paramagnet by a first order phase transition, and exhibits correlations matching those predicted for a Coulomb phase of rank-2 electric fields with scalar charges [29].Elementary excitations are Type-I fractons, appearing at the corners of membranes of flipped spins.
The discovery of a relatively simple Ising model, with finite-range, two-body interactions establishes a useful platform for the further exploration of fractonic physics.This could include the perturbative introduction of quantum effects via transverse fields or transverse exchange.This may be a better setting in which to study quantum effects on fractons than in the Heisenberg models suggested in Ref. 27, for which numerical calculations suggest that quantum fluctuations wash out the multifold pinch points [39].Here, the emergent Gauss's law is protected by a finite gap, so it may be more robust.Even if instanton effects drive the emergent gauge theory into a confined phase (as they do for the ordinary U(1) gauge theory in 2+1 D), the low temperature physics can still show interesting features related to the liquid phase.
Our purpose-built Monte Carlo algorithm provides a template for future numerical studies of Type-I fractonic models.Future work could address in more detail the topics of relaxation and disorder-free glassiness.The successful demonstration of an Ising fracton spin liquid, based on a Hamiltonian originally constructed for continuous spins [27], also raises the question of whether other classical spin liquids with higher-moment conservation laws [40,41] have Ising realizations, and what their properties may be.
Given the rapid progress in physical simulation platforms based on Rydberg atoms and superconducting qubits, the dynamics of fractons in our model may not be far from exploration in the laboratory.In Fig. S4, we show the structure factor as a function of temperature for a system with L = 24.Above the ground state regime, at T = J/4, the four fold pinch points are already beginning to form, but only become sharp after the first-order transition.

B. Higher-order correlations
As discussed in the main text, an open question in the absence of an exact solution is the degree to which the gorund state manifold breaks symmetry.There are no strong indicators for symmetry breaking in the two-point correlation function (see structure factor shown in Fig. 4 of the main text).However, in a system with such unusual constraints one can imagine also order showing up only -or at least more prominently -in higher-order correlation functions.While an exhaustive study of those is obviously not possible, a good intuition for order through entropy is often obtained from correlations corresponding to the local moves connecting different ground states.To this end, we show in Fig. S5 what we call the "gauge structure factor" in the low temperature regime for a system size of L = 24.This quantity if defined as the Fourier transform of the correlation function between flippable gauge motifs.That is for each hexagon, we define a "gauge flippability" as g(h) = 0 gauge flip centered at h costs nonzero energy ±1 else, sign depending on the current configuration (S2) and accordingly While g(h) quantity shows short range correlations as evidenced by the structure factor in Fig. S5, we find no FIG. 1. Ground states of the honeycomb model fulfill the constraint M h = 0, illustrated in (a).A single spin flip, as indicated by a red circle in (b) preserves total charge ρ as well as all dipole moments.The minimal-energy excitation, shown in (c) hence creates four distinct defects.These can be moved in pairs by flipping another four spins as indicated by orange dashed lines.
FIG. 2. Splitting of lineons and local resonance moves.(a) Pairs of fractons form lineons, mobile along a one-dimensional submanifold.These lineons can split in two, either the forwards or backward directions.(b) A local combination of six membranes, can create, split and recombine lineons, in such a way as to reach a new ground state configuration.This amounts to flipping 24 spins.Application of this move to the ground state in (c) does not create any excitation and could be centered around any of the green shaded hexagons.The move can be understood as a local combination of six membranes as shown in (b), where each corner overlaps with exactly one other corner of opposite charge.This is possible because as shown in (c), a single two-fracton lineon can split into two both in the forward (left) and in the backward direction (right).

FIG. 3 .
FIG.3.Thermodynamics of the fractonic Ising model.Energy E, Specific heat Cv, and residual entropy S0 are obtained from Monte-Carlo simulations.Energy/specific heat are consistent with a first-order transition from a high-temperature paramagnetic phase to a low-temperature phase where the constraint M h = 0 is satisfied.Finite-size extrapolation of the ground state entropy per site tends to a value significantly higher than the lower bound discussed in the main text, indicating a strongly fluctuating regime below the first order transition.

FIG. S4 :
FIG.S4: Structure factor as a function of temperature for system size L = 24.While the onset of the four-fold structure at the Brillouin-zone corners is already visible above the transition (T = J/4), the pinch points only sharpen below the transition (T = J/5)