(Four) Dual Plaquette 3D Ising Models

A characteristic feature of the 3d plaquette Ising model is its planar subsystem symmetry. The quantum version of this model has been shown to be related via a duality to the X-Cube model, which has been paradigmatic in the new and rapidly developing field of fractons. The relation between the 3d plaquette Ising and the X-Cube model is similar to that between the 2d quantum transverse spin Ising model and the Toric Code. Gauging the global symmetry in the case of the 2d Ising model and considering the gauge invariant sector of the high temperature phase leads to the Toric Code, whereas gauging the subsystem symmetry of the 3d quantum transverse spin plaquette Ising model leads to the X-Cube model. A non-standard dual formulation of the 3d plaquette Ising model which utilises three flavours of spins has recently been discussed in the context of dualising the fracton-free sector of the X-Cube model. In this paper we investigate the classical spin version of this non-standard dual Hamiltonian and discuss its properties in relation to the more familiar Ashkin–Teller-like dual and further related dual formulations involving both link and vertex spins and non-Ising spins.


Introduction
The Kramers-Wannier dual [1] of the classical Ising Hamiltonian with nearest neighbour ij couplings on a 3d cubic lattice is the 3d Ising gauge theory where the sum runs over plaquettes and the gauge spins U live on the edges of the plaquettes. The coupling β in the partition function Z(β) = ∑ {σ} exp(−βH Ising ) and its dual β * in Z(β * ) = ∑ {U} exp(−β * H Gauge ) are related by β * = − 1 2 log[tanh(β)]. We use un-superscripted variables, e.g., U, σ i , τ i , µ i , to denote spins in classical Hamiltonians and superscripted variables, e.g., σ x,z i , τ x,z i , µ x,z i , to denote the Pauli matrices appearing in quantum Hamiltonians. The positional subscript indices i, j, k . . . are occasionally omitted for brevity.
In this paper we will investigate the relation between four (apparently) different formulations of the dual to the 3d plaquette Ising model, which has also been dubbed the gonihedric Ising model [2][3][4][5] H κ=0 = − ∑ σ i σ j σ k σ l .
(3) The degeneracy affects the finite size scaling behaviour at the first order transition [25][26][27][28][29], changing the universal 1/L 3 finite size scaling shift in estimates of a first order transition point on an L 3 lattice (with periodic boundary conditions) [30,31] to 1/L 2 . For non-zero κ the planar subsystem symmetry appears to be broken at finite temperature [32,33] and the transition becomes second order. The Kramers-Wannier dual to H κ=0 takes the form of an anisotropic Ashkin-Teller model [34]. It still possesses a planar subsystem symmetry and degenerate low temperature phase, so the modified finite size scaling at the first order transition is observed there also [25][26][27][28][29].
The subsystem symmetry in the quantum spin version of the 3d plaquette Ising model has recently been shown to be closely linked to the properties of the X-Cube model [35], which has become a paradigmatic model for the new and rapidly developing field of fractons, which are quasiparticles with restricted mobility in isolation. Some recent reviews of what is now a burgeoning fracton literature can be found in [36,37]. To see the role played by the subsystem symmetry in constructing the X-Cube model, first consider gauging the global Z 2 symmetry in the case of the 2d quantum transverse spin Ising model This can be done by introducing τ z on the links and an additional plaquette flux term to endow the link spins with dynamics, which gives a gauge-invariant Ising (or Z 2 gauge-Higgs [38][39][40]) model where we have dropped the link indices on the τ z for conciseness. The gauge-invariant sector of the high temperature phase, β → 0, of this model, where σ x i ∏ k∈+,i τ x k = 1 and k labels the four edges (+) incident to vertex i, gives Kitaev's Toric Code model [41,42] We use the gauge invariance to trade σ x i for ∏ k∈+,i τ x k , leaving the mutually commuting terms and customarily set h = β p = 1. The Toric code displays topological order and has anyonic quasiparticle excitations. On the other hand, gauging the subsystem symmetry of the 3d plaquette Ising model in a similar manner leads to the X-Cube model [35]. In this case, when we start with the quantum transverse spin 3d plaquette Ising model Hamiltonian gauging the Z 2 subsystem symmetry requires inserting a τ z which lives on the plaquettes The equivalent of the plaquette flux term in the Toric Code derivation is now a set of three "X" terms as shown in Figure 2, one in each lattice plane B xy,yz,xz i = ∏ j∈+,i τ z j . If we again consider the gauge invariant sector σ x i ∏ k∈ ,i τ x k = 1, where the τ x k live on the twelve incident plaquettes impacted by flipping the single central spin at site i, the high temperature limit β → 0 produces the X-Cube Hamiltonian where it is simpler to think of the τ x , τ z 's residing on the links of the dual lattice. The A term is a product of all the τ x around a cube and the B terms are the three "crosses" of τ z 's shown in Figure 2.  Figure 2. The terms contributing to the X-Cube Hamiltonian. The cube A term is a product of the twelve τ x spins on the edges of the cube and the three B "X" terms composed of τ z spins lie in each of the three lattice planes as shown on the corner.
The remaining couplings have again been set to one. The quasiparticles arising from defects in the A terms are fractons and cannot move in isolation, whereas the defects in the B terms give lineons, which can only move in straight lines. The order in the X-Cube model is not topological. It has an exponential, but sub-extensive, ground state degeneracy inherited from the plaquette Ising model as a consequence of the subsystem symmetry.
It was observed recently in [43] that the Hamiltonian for the fracton-free subsector (where all the A cube terms are +1) of the X-Cube model in a transverse field could be written in terms of a dual Hamiltonian (at the risk of causing confusion we have kept the notation of [35] for the A and B terms rather than [43], which swaps A and B, though we denote the edge Pauli matrices by τ rather than σ in distinction to both [35,43]) with three flavours of Ising spins σ i , τ i , µ i living on the vertices of the cubic lattice rather than the links where the nearest neighbour sums in the four spin terms each run along one of the orthogonal axes, with ij, ik and jk representing the z, y and x axes respectively. The constraint on the A terms is automatically resolved by these spins. In this paper we discuss the properties of the classical spin version of this Hamiltonian, dubbed H dual2 for reasons to be explained in the next section. We shall see that it is closely related via a gauge-fixing to the Ashkin-Teller-like [44] Hamiltonian, H dual1 , constructed using the classical Kramers-Wannier duality from the 3d plaquette Ising model, as well as though a decoration transformation to a third Hamiltonian, H dual3 , which mixes edge and vertex spins. We find that the characteristic planar subsystem symmetry of the 3d plaquette Ising model is still present in H dual1,2,3 and also that the interesting, possibly glassy, dynamical properties of the 3d plaquette model are also apparent in the duals. The Hamiltonians H dual2,3 are already implicit in the discussion by Savvidy and Wegner in [45] in the context of the general framework for dualities [46] in spin models.

Duals Galore
The Kramers-Wannier dual to H κ=0 was initially constructed by Savvidy et al. [34] by considering the high temperature expansion of the plaquette Hamiltonian which can be written as on an L 3 cubic lattice, where the sum runs over closed surfaces with an even number of plaquettes at any vertex. In the summation n(S) is the number of plaquettes in a given surface. The low temperature expansion, i.e., high temperature in the dual variable of the following anisotropic Hamiltonian produced the requisite diagrams. In H dual0 the sums are one-dimensional and run along the orthogonal axes, with ij, ik and jk again representing the z, y and x axes respectively using our conventions. The spins are non-Ising and live in the fourth order Abelian group, since the geometric constraints on having an even number of plaquettes at each vertex mean that with e being the identity element. They can be thought of as representing differently oriented matchbox surfaces such as that shown in Figure 3, which are combined by facewise multiplication. The shaded faces carry a negative sign and the associated spin variable lives at the centre of the matchbox. Any spin cluster boundary in the model can be constructed from such matchboxes while still satisfying the local constraint on the number of incident plaquettes.
The spins may also be taken to be Ising (±1) variables if we set η i = σ i τ i , which is more convenient for simulations. This modifies H dual0 to an anisotropically coupled Ashkin-Teller Hamiltonian [44] This formulation of the dual model was first investigated numerically in [47] and it was found that it displayed a first order phase transition and a similar planar subsystem symmetry to that of H κ=0 . The continued presence of the subsystem symmetry was a consequence of the anisotropic couplings, which allowed a greater freedom in transforming the spin variables than in the isotropically coupled version of Equation (19), which is just the Ashkin-Teller model at its four-state Potts point. It is possible to construct H dual1 and its higher dimensional equivalents [45] using the general framework for duality in Ising lattice spin models that was first formulated by Wegner in [46]. Suprisingly, there are two further possible ways to write the dual to H κ=0 in three dimensions with this machinery, using either the general formula for the dual of codimension one surfaces or the formula for the dual of two dimensional surfaces in d dimensions. If we temporarily use the notation of [45], the dual Hamiltonian for a codimension one surface in d dimensions is given there by where the Λ spins live on each of the (d − 3) dimensional (hyper)vertices situated at the vertices r of the hypercubic lattice and the indices α, β, γ run from 1 to d. The unit vectors e γ point along the lattice axes. On the other hand, the dual Hamiltonian for a two-dimensional gonihedric surface embedded in d dimensions is of the form where we now have Γ spins on each (hyper)edge in addition to the Λ spins at each vertex.
If we specialize to two dimensional plaquette surfaces embedded in a cubic lattice in three dimensions, which is the case for the dual of H κ=0 , either formulation may be employed since this is both a codimension one surface and a two-dimensional surface embedded in three dimensions. Returning to our own notation, the codimension one Hamiltonian of Equation (20) in three dimensions may be written as which is just the Hamiltonian of Equation (14) that appeared as the classical spin limit of the dual to the fracton-free subspace of the X-Cube model. The three flavours of spins living at each vertex display a local Ising gauge symmetry σ i , τ i , µ i → γ i σ i , γ i τ i , γ i µ i in addition to the planar subsystem symmetry shared with H κ=0 and H dual1 , as we shall see presently.
Still within the general approach of Savvidy and Wegner [45,46], in three dimensions the Hamiltonian of Equation (21) for the two-dimensional surface variant also contains three flavours of vertex spins σ i , τ i , µ i , but in addition there are gauge-like spin variables U 1,2,3 ij living on the lattice edges which couple in an anisotropic manner to the vertex spins We thus have four different Hamiltonian formulations for the dual of the plaquette Hamiltonian H κ=0 in three dimensions: • H dual0 in Equation (17) (14) containing purely four spin interactions.
• H dual3 in Equation (22) containing both vertex spins and gauge-like edge spins.
We have already seen that setting η i = σ i τ i in H dual0 in Equation (17), with σ i , τ i being Ising spins, keeps the algebra of Equation (18) intact and gives the Ashkin-Teller Hamiltonian of H dual1 in Equation (19).
In the next section we discuss the relation between the four-spin Hamiltonian H dual2 of Equation (14) and the gauge-spin Hamiltonian H dual3 of Equation (22), and thereafter that between H dual2 and H dual1 .

Decoration
The equivalence between H dual3 and H dual2 is a consequence of a variation of the classical decoration transformation [48]. In the standard transformation an edge with spins σ 1 , σ 2 at each vertex is decorated with a link spin s as in Figure 4.  If the coupling between s and σ 1 and σ 2 isβ, summing over the central spin s gives rise to a new effective coupling β between the primary vertex spins σ 1 , σ 2 ∑ s exp β s(σ 1 + σ 2 ) = A exp(βσ 1 σ 2 ).
Both the prefactor A and the coupling β may be expressed in terms ofβ by enumerating possible spin configurations in Equation (23). This gives We can repeat this procedure with the U spins on each edge in H dual3 . In this case each direction has two flavours of vertex spin and performing the sum generates the four-spin couplings of H dual2 , for example where A and the relation between β,β are the same as in Equation (24). The sum over U may be carried out globally over every edge which immediately demonstrates equivalence of the partition functions for H dual3 and H dual2 The overall factor B coming from a product of A's on the individual links is irrelevant for calculating physical quantities and the two couplings are again related by the decoration relation,

Gauge Fixing and Subsystem Symmetry
The equivalence between H dual2 and H dual1 , on the other hand, is a consequence of a gauge symmetry which is present in H dual2 [49] We are at liberty to choose the Ising spin gauge transformation parameter γ i to be equal to one of the spin values, say µ i , at each site so the gauge transformation then becomes which, using the fact that the sum over the remaining spin variables σ i , τ i is invariant under the transformation, relates the partition functions for the two Hamiltonians as The coupling β is not transformed in this case and we can, of course, choose to eliminate any one of the three spins, which simply amounts to relabelling the axes. From this perspective H dual1 is simply a gauge-fixed version of H dual2 . This can be confirmed by Monte-Carlo simulations which measure the same energies (and energy distributions) and transition points for the observed first order phase transitions [49]. The equivalence between H dual3 and H dual2 described in the preceding section via the decoration transformation also sheds light on the presence of this gauge symmetry in H dual2 . All the terms in H dual3 are of the gauge-matter coupling form σ i U ij σ j , so this action possesses a similar, standard gauge invariance to that seen in other gauge-matter systems such as the gauge-Ising model of Equation (7), namely When the U spins are summed over to give H dual2 , the gauge symmetry of the σ, τ and µ spins in Equation (27) remains as an echo of this symmetry. In both cases if we look at a single site transformation all three spins σ i , τ i and µ i must be transformed. In H dual3 this is a consequence of the way in which the three edge spins U 1,2,3 ij couple to the vertex spins. A characteristic feature of H dual1 is the planar subsystem symmetry intermediate between a gauge and a global symmetry, just as with the 3d plaquette Ising Hamiltonian. For H dual2 the anisotropic couplings mean that it is still possible to flip planes of one of the spins (the one which is "missing" from the interactions in the direction perpendicular to the planes) at zero energy cost as shown in Figure 5. It is also possible to flip two or three orthogonal faces on the cube, so tiling the entire lattice with such combinations we can see that in addition to the purely ferromagnetic ground state we may have arbitrary (and possibly intersecting) flipped planes of spins.
The ground state structure, and the mechanism of anisotropic couplings which allows the plane spin flips, is thus identical to that in H dual1 , whose possible ground states on a single cube we recall for comparison in Figure 6. Flipping the third spin µ in H dual2 , which is absent in H dual1 , is replaced by flipping both the σ and τ spins in H dual1 , consistent with the gauge transformation relating the two Hamiltonians. In summary, the ground state structure of H dual2 shows an interesting interplay between the gauge symmetry of Equation (27) and the subsystem symmetry. The local gauge symmetry allows one to reduce the effective number of degrees of freedom and recover the ground state structure of H dual1 . H dual3 is a similar case since the geometrical arrangement of the couplings is similar in spite of the presence of the additional spins U on the links. Each of the σ, τ, µ spins couple in two directions, which define the lattice planes in which they can be flipped without affecting the energy.

Indicative Monte-Carlo
Low precision Monte-Carlo simulations using simple Metropolis updates found a first order phase transition in both H dual2 and H dual1 in the region β 1.3 − 1.4 [47,49]. Much higher precision simulations using multicanonical methods were later carried out for the original 3d plaquette Ising model H κ=0 and the Ashkin-Teller dual H dual1 in order to accurately determine the transition point (β ∞ = 1.31328 (12) in the case of the dual model) and confirm the non-standard finite size scaling that is a consequence of the exponential degeneracy of the low temperature phase [25][26][27][28][29]. Even with the modest statistics and the use of a Metropolis update in the simulations of [47,49] a sharp drop in the energy, as would be expected for a first order transition, is clearly visible in the region of the transition point. A plot of the energy is shown for various lattice sizes in Figure 7 for H dual2 and the values for H dual1 are essentially identical.  The first order nature of the transition for H dual2 and H dual1 can be further confirmed by observing a dual peak structure in the energy histogram P(E) near the transition point and a non-trivial value of Binder's energy cumulant as a consequence of the shape of P(E) [25,49]. Based on these observations, and allowing for a factor of 1/2 in our definitions of H dual1 and H dual2 in [47,49], we would expect to see a transition in H dual3 at the the value ofβ found by inverting the decoration transformation, namely 1 2 cosh −1 (exp(1.3 − 1.4)) = 0.99 − 1.04 in the thermodynamic limit. To confirm this expectation, we carried out Monte-Carlo simulations using 10 3 , 12 3 , 16 3 and 18 3 lattices with periodic boundary conditions for all spins at various temperatures, again with a simple Metropolis update. After an appropriate number of thermalization sweeps, 10 7 measurement sweeps were carried out at each lattice size for each temperature.
Looking at measurements of the energy from our simulations of H dual3 in Figure 8 we can see that a similar sharp drop in the energy consistent with a first order transition is still present. The relatively low statistics and the use of a Metropolis update for the data presented here for H dual2 and H dual3 mean that a high accuracy extrapolation using the correct 1/L 2 finite size scaling for the transition point is not feasible, which would require more extensive multicanonical simulations along the lines of [25,26,29]. Nonetheless, the agreement of the suitably transformed finite size lattice transition points in the Monte-Carlo simulations confirm that H dual3 and H dual2 are related by the decoration transformation and the first order nature of the transition for H dual3 is clear from the dual peak form of P(E) near the transition point in Figure 9.

Dynamics
Another interesting feature of the original plaquette Hamiltonian H κ=0 is its dynamical behaviour. It possesses a region of strong metastability around the first order phase transition and displays glassy characteristics at lower temperatures [18][19][20][21][22][23][24] with non-trivial ageing properties. We found that the Ashkin-Teller dual Hamiltonian H dual1 also shares these characteristics since it failed to relax to the equilibrium minimum energy of E = −1.5 when cooled quickly from a hot start [47]. Note that in this case, since we are exploring the real time dynamics of the system, simulations with a Metropolis update are preferable to more sophisticated algorithms.
H dual2 displays identical behaviour under cooling to H dual1 . We consider 20 3 , 60 3 and 80 3 lattices which are first equilibrated in the high temperature phase at T = 3.0 and then cooled at different rates to zero temperature. The energy time series is recorded during this process. In Figure 10 we can see that with a slow cooling rate of δT = 0.00001 per sweep, the model still relaxes to a ground state with E = −1.5 for all the lattice sizes. However, as can be seen in Figure 11 with a faster cooling rate of δT = 0.001 per sweep the model no longer relaxes to the ground state energy of E = −1.5, but is trapped at a higher value, which is around −1.415 for the larger two (60 3 and 80 3 ) lattices. Whether the observed behaviour under cooling is a sign of genuine glassiness in H κ=0 or not remains a matter of debate and similar considerations would apply to the dual Hamiltonians discussed here.

Discussion
Motivated by recent work on the quantum X-Cube Hamiltonian and related dual models [43] we revisit various formulations of classical spin Hamiltonians dual to the 3d plaquette Ising model. We describe the following chain of relations between these models where we have indicated the operations relating the various Hamiltonians. A variant of the decoration transformation in which edge spins are summed out relates H dual3 to H dual2 . In transforming H dual3 to H dual2 the coupling is therefore transformed as β = 1 2 ln[cosh(2β)]. The gauge-invariant nature of H dual3 due to the presence of both edge and vertex spins leaves an echo in the vertex spin gauge symmetry of H dual2 , which in turn ensures the equivalence of H dual2 and H dual1 via a gauge-fixing. Allowing non-Ising spins gives a final equivalence between the dual models H dual1 and H dual0 and a standard Kramers-Wannier duality transformation then takes us back to the 3d plaquette Ising Hamiltonian of H κ=0 where the story began.
The planar subsystem symmetry of this 3d plaquette Ising Hamiltonian H κ=0 remains a feature of the various dual Hamiltonians and affects the finite size scaling properties at the first order transition displayed by these models, just as with H κ=0 . The nature of the order parameter for the various duals, and indeed H κ=0 itself, remains to be satisfactorily clarified. An attempt at this has been made for H κ=0 in [50] and indeed appeared to give sensible numerical results in [25][26][27][28][29]. This was based on the observation by Suzuki et al. [51][52][53] that an anisotropic version of the 3d plaquette Ising model with open boundary conditions could be transformed to an uncoupled stack of 2d Ising models and suggested that a two spin correlator summed perpendicular to lattice planes might still serve as an order parameter in the isotropically coupled case of H κ=0 . It would be more satisfactory to have a less heuristic approach to an order parameter based on a clearer understanding of the nature of the low temperature order in the 3d plaquette Ising Hamiltonian. In this respect a study of the order parameters in the various dual models might be helpful.
In a similar vein, the principal interest in [43] was actually investigating "odd" variants of fracton models, in which the signs of some of the terms in the Hamiltonians were reversed, leading to frustrated models. The geometrical nature of the order in such frustrated spin models would be of interest in the classical case too. Finally, the dynamics of the various classical Hamiltonians discussed here display glassy features. The question of whether the glassy dynamics of the quasiparticle excitations in the quantum models [54,55] offers any insights into this behaviour may be worth pursuing.