Exotic Symmetry Breaking Properties of Self-Dual Fracton Spin Models

Fracton codes host unconventional topological states of matter and are promising for fault-tolerant quantum computation due to their large coding space and strong resilience against decoherence and noise. In this work, we investigate the ground-state properties and phase transitions of two prototypical self-dual fracton spin models -- the tetrahedral Ising model and the fractal Ising model -- which correspond to error-correction procedures for the representative fracton codes of type-I and type-II, the checkerboard code and the Haah's code, respectively, in the error-free limit. They are endowed with exotic symmetry-breaking properties that contrast sharply with the spontaneous breaking of global symmetries and deconfinement transition of gauge theories. To show these unconventional behaviors, which are associated with sub-dimensional symmetries, we construct and analyze the order parameters, correlators, and symmetry generators for both models. Notably, the tetrahedral Ising model acquires an extended semi-local ordering moment, while the fractal Ising model fits into a polynomial ring representation and leads to a fractal order parameter. Numerical studies combined with analytical tools show that both models experience a strong first-order phase transition with an anomalous $L^{-(D-1)}$ scaling, despite the fractal symmetry of the latter. Our work provides new understanding of sub-dimensional symmetry breaking and makes an important step for studying quantum-error-correction properties of the checkerboard and Haah's codes.


I. INTRODUCTION
A major task in the theoretical study of topological quantum computation is the search for the most resilient quantum error correction (QEC) codes.Although the two-dimensional (2D) toric code [1][2][3] and color code [4][5][6] have a code capacity p c ≃ 10.9% against random qubit errors and are ideal for topological quantum memory, they permit only transversal implementations of logical Clifford gates as any 2D topological stabilizer code [7,8].Such gates can be efficiently simulated in polynomial time on classical computers [9].Therefore, despite the striking experimental advances of these codes on trapped-ion [10][11][12] and superconducting-qubit [13][14][15] platforms, in the long term, we have to go beyond 2D to realize the desired quantum advantages with topological protection.However, in three dimensions, thresholds of standard topological stabilizer codes are significantly reduced by the multi-spin interactions and gauge symmetries in their error modeling.For instance, the effective spin models determining the stability of 3D toric codes and color codes induce random Z 2 and Z 2 × Z 2 gauge theories and exhibit optimal minimum thresholds ≃ 3.3% [3,16] and ≃ 1.9% [17] against random qubit errors, respectively.This hence motivates the quest for novel codes in three and higher dimensions.
Fracton codes represent a novel class of 3D topological stabilizer codes [18][19][20][21].Such codes support subextensive ground-state degeneracies (GSDs) and can provide larger cod-ing spaces than standard topological codes with constant degeneracy.Moreover, their gapped excitations conform to intrinsic mobility constraints that can suppress error propagations and maintain the correctability of the codes.Strong error resilience of fracton codes because of this latter feature is indeed confirmed in our recent work [22].There, we determined the error thresholds of the representative X-cube code with high accuracy and found that its code capacity, p c ≃ 7.5%, is remarkably higher than other 3D codes [3,16], including 3D toric and color codes [17].
This work aims to carry out a necessary move toward understanding the QEC properties of the type-I checkerboard code and type-II Haah's code.In Ref. [22], we conjectured that these two self-dual fracton codes may reach the optimal error threshold p c ∼ 11% as a quantum memory.Accurate determination of their thresholds requires constructing a disordertemperature phase diagram for each model and demands exceeding computational resources owing to the interplay between fracton and spin-glass physics.As a first step, we focus in the current work on the error-free limit of the two codes and investigate phase transitions and construct order parameters of their dual spin models.Understanding the associated sub-  3 and face centers of an FCC lattice, respectively.The red and green tetrahedra show an example of the a s v+a and a s v−a interaction terms, respectively, where a ∈ {(0, 0, 0), (1, 1, 0), (1, 0, 1), (0, 1, 1)} labels the four FCC sublattices.The entire lattice can be intuitively visualized with small cubes: Each shaded (empty) cube contains a single red (green) tetrahedron interaction.(b) A unit cell of the fractal Ising model on a simple cubic lattice.Vertices are again denoted by blue circles, but here v = (x, y, z) ∈ Z 3 .The red and green tetrahedra represent the nearest-neighboring and nextnearest-neighboring interaction terms a 1 s v+a 1 and a 2 s v+a 2 , respectively, with a 1 ∈ {(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)} and a 2 ∈ {(0, 0, 0), (1, 1, 0), (1, 0, 1), (0, 1, 1)}.
dimensional symmetry breaking (SDSB) and non-local order parameters can guide future studies at finite error rates.
At the error-free limit, duals of the checkerboard and Haah's codes lead to the disorder-free tetrahedral Ising model (TIM) and fractal Ising model (FIM), respectively.Subextensive topological GSDs of the fracton codes are reflected by distinct sub-dimensional symmetries of the dual-spin models.Such symmetries act on lower-dimensional manifolds or fractals of the system [58] and are intermediate between global and gauge symmetries.Spontaneous symmetry breaking remains possible but produces long-range orders without a local order parameter.The associated phase transitions fit neither the Landau paradigm nor the standard confinementdeconfinement scenario.Their finite-size scalings also exhibit anomalous behaviors, despite both spin models investigated here experiencing a strong first-order transition.
The manuscript is organized as follows.In Sec.II, we define the two fracton spin models, discuss their symmetry properties, and construct appropriate correlators and order parameters.Sec.III is devoted to the fractal symmetry generators and the consequent degeneracy of the fractal Ising model, utilizing a polynomial ring formalism.Sec.IV discusses the phase transitions and scaling behaviors of the two models.We conclude at Sec. V with an outlook.App.A and App.B include details of fractal generators and multi-canonical simulations.

II. FRACTON SPINS MODELS
The tetrahedral and fractal Ising models relate to the checkerboard and Haah's codes in two ways.One way is using a generalized Wegner's duality [59,60] where the two fracton codes are obtained by gauging these two spin models [18].In this correspondence, immobile excitations in a fracton code map to domain-wall corners of its dual spin model.Furthermore, a fracton topological order is dual to an SDSB phase on the spin-model side.
The other way is through a statistical-mechanical mapping that describes the error-correction procedure of a topological code.This approach is initially developed for standard gauge codes, such as toric and color codes [2,6,17], and has been recently extended to fracton codes [22].Under this mapping, a fracton spin model corresponds to the error-free limit of a fracton code, while qubit errors are reflected by including quenched disorders.
Here, we focus on the disorder-free limit, but will study the more involved disorder-full case in a future work.
The Hamiltonian is given by which consists of four-body interactions living on elementary tetrahedra, and each site is shared by four J + and four J − tetrahedra.We set the interaction to be isotropic as H TIM is invariant under flipping all spins, including face centers, of an arbitrary xy, yz, or zx-layer.This is easy to verify as either no or two spins in the interaction terms of Eq. ( 1) are affected.There are in total 3L such plane-flip symmetries for a lattice of size L × L × L. Specifically, we use g x i to denote the symmetry of flipping spins in the i th xy-layer, located at x = i.Similarly, g y j and g z k are used to denote the symmetries at y = j and z = k, respectively.These symmetries are subject to three global constraints, Thus, only 3L−3 of them are independent, leading to a subextensive ground state degeneracy with log 2 GSD = 3L − 3.
As two degenerate ground states differ by flipping at least 1 2 L 2 spins, i.e., an entire plane of the lattice, no finite-order perturbations can connect them at the thermodynamical limit.Therefore, a long-range order spontaneously breaking the plane-flip symmetry remains permitted.However, characterizations of such ordering are fundamentally distinguished from breaking a global symmetry.As a necessary condition, the magnitude of a physical correlator or an order parameter needs to be preserved under arbitrary plane flips, which excludes any local order parameter.Without losing generality, we consider the following minimal correlator where the hat symbol indicates a unit vector and r ∈ 2Z + 1. G z TIM (r) spans an irregular tetrahedron whose first and last pairs of spins belong to two different xy-layers separated by a distance r.At the r → ∞ limit, it relates to an order parameter , with Here, q z TIM can be viewed as an extended or semi-local ordering moment in an xy-layer, but itself is a linewise object that represents the total moments of local correlators s v s v+ x+ŷ and s v+ŷ+ẑ s v+ x+ẑ along an entire z-line of FCC unit cells.In view of dimensional reduction, it has a characteristic dimension dim(q z TIM ) = 1 and a codimension codim(q z TIM ) = D − dim(q z TIM ) = 2.We hence refer to Q z TIM as a sub-dimensional order parameter to distinguish it from local order parameters and Wilson-loop order parameters.
One can similarly construct two equivalent order parameters Q x TIM and Q y TIM .It is sufficient to use any of them.
The Hamiltonian consists of two four-body interactions on the tetrahedra specified by {a 1 } and {a 2 }, as depicted in Fig. 1, where one s v participates four J 1 tetrahedra and four J 2 tetrahedra.For simplicity, we take J 1 = J 2 = 1.
In correspondence of Haah's code, H FIM may have a fractal symmetry under PBCs.There is no simple algebraic expression of its fractal symmetry generators and ground states [61].However, they can be constructed systematically by the approach provided in Sec.III.
The fractal symmetry leads to more exotic order parameters.We can define invariant correlators by isotropically scaling the two interaction terms in H FIM , and G ′ FIM (r) measure the long-range correlations of spins at the four corners of tetrahedra.
We will discuss in Sec.III that locations of spins transformed by a symmetry operation depend on choices of fractal generators.As a consequence, neither local ordering nor semi-local ordering like Q z TIM can be well defined.This property also relates to the lack of string operators in Haah's code [20,62] whose locally creatable excitation patterns correspond to the interaction terms of H FIM [63].
Nonetheless, G FIM 0 and G ′ FIM 0 at r → ∞ indicates a long-range order, and G FIM = G ′ FIM = 1 is only possible for ground states.Hence, both G FIM and G ′ FIM can serve as order parameters to describe the fractal symmetry breaking of H FIM .We refer to these non-local correlators as fractal order parameters to distinguish them from sub-dimensional order parameters that admit a semi-local or dimension-reduction interpretation.
The non-local fractal order parameter is also essentially distinct from Wilson loop correlators in gauge theory.Wilson loops generally vanish; deconfined and confined phases are distinguished by their decaying behaviors, namely the perimeter law and the area law, instead of the expectation values [64].In contrast, in the fractal symmetry breaking phase, G FIM and G ′ FIM can robustly maintain a finite value at an arbitrarily large distance.

C. Self-duality
From the viewpoint of quantum error correction, selfduality of the tetrahedral and fractal Ising models reflects the equivalent role of Pauli X and Z errors in the checkerboard and Haah's codes.At the level of spin models, this self-duality can also be understood directly through the Kramers-Wannier duality [65,66].
The dual of H TIM is given by placing a dual Ising spin s v ⋆ +a at the center of an orignal J + tetrahedron a s v+a .It is convenient to work with the shaded and empty unit cubes in Fig. 1, which are exclusively occupied by J + and J − tetrahedra, respectively.The dual lattice remains an FCC lattice with a uniform shift of 1 2 ( x, ŷ, ẑ); namely, the dual vertices and face centers are labeled by ).An original spin s v+a has four neighboring shaded cubes, whose centers span a dual J ⋆ − tetrahedral a s v ⋆ −a .Correspondingly, centers of the four empty neighboring cubes of s v+a span a dual J ⋆ + tetrahedral a s v ⋆ +a .Thus, the dual Hamiltonian preserves the form of H TIM , while the shaded and empty cubes are swapped on the dual side.
The dual of H FIM can be analyzed analogously.We again place a dual Ising spin s v ⋆ at the center of an original cube.The dual lattice is a simple cubic lattice with vertices given by v ) .The two tetrahedron interactions in Eq. ( 5) are dual to a 1 s v ⋆ −a 1 and a 2 s v ⋆ −a 2 , respectively, where the reverse sign of a 1 and a 2 merely swaps the incidence of spins in each unit cube.Hence, the fractal Ising model is also self-dual.
In virtue of this self-duality, we immediately have Here, i denotes the (dual) Boltzmann factor of a single interaction term h (⋆)  i , and i is a shorthand label distinguishing those individual tetrahedra in (dual) H TIM or k ω ⋆ β ⋆ ,k e −inkπ , namely, Given the self-duality and provided they only have a single phase transition, which is the case for both of H TIM and H FIM (Sec.IV), Z (β) and Z β ⋆ shall experience the same singularity.Hence, at the transition point whose solution gives the self-dual point

III. FRACTAL SYMMETRY & POLYNOMIAL RING REPRESENTATION
We now discuss the symmetry of the fractal Ising model H FIM utilizing a polynomial ring formulation.This formalism is introduced in Ref. [63] for describing translational invariant spin Hamiltonians and offers a convenient way to characterize fractal symmetries.Before proceeding to compute the fractal symmetry and GSD of H FIM , we include a brief review of this formalism.

A. Polynomial ring formalism
The polynomial ring formalism is set up through a group ring R = Z 2 [Λ], with Λ = {x i y j z k |i, j, k ∈ Z L } denoting the lattice translation group.Namely, the coordinate of each vertex in the lattice is represented in a multiplicative notation, and PBCs are imposed by identifications x L = y L = z L = 1.The ring R can be viewed as a set consisting of all polynomials of the form The coefficients a i jk = 0, 1 ∈ Z 2 are integers modulo 2. Intuitively, a polynomial f ∈ R represents a collection of vertices or lattice vectors whose coefficients are a i jk = 1.
For example, the sets {a 1 } and {a 2 } in Eq. ( 5) holding the four-spin interactions correspond to the following two polynomials Furthermore, we use X ( f ) to denote the operation of flipping spins on the set of vertices specified by f .In the group ring formalism, the excitation pattern resulting from X( f ) acting on a ground state can be specified by an R-linear map from R to R 2 , which takes the form Here, standard polynomial multiplications are employed, and are the spatial inversion of ε 1 and ε 2 , respectively.We use the notations x ≡ x −1 , y ≡ y −1 , and z ≡ z −1 for brevity.The Rlinearity in Eq. ( 13) is a consequence of translation symmetry.Physically, the polynomial pair (ε 1 f, ε 2 f ) describes the excitation pattern of the two interaction terms in H FIM .This can be understood by considering a simplest polynomial f = 1.In this example, φ(1) = (ε 1 , ε 2 ) labels the locations of the excited J 1 and J 2 tetrahedra in H FIM due to a single spin flip at x 0 y 0 z 0 .

B. Polynomial representation of spin-flip symmetries
Operation X ( f ), which flips the spins on the set of vertices specified by f , is a symmetry if and only if it creates no excitations, namely, φ ( f ) = 0. Therefore, the kernel of φ, represents the group of spin flip symmetries of H FIM .Explicitly, ker φ is specified by finding all the solutions f ∈ R that satisfy ε 1 f = ε 2 f = 0.As we show in App.A, the first equation ε 1 f = 0 implies that f can be expressed as Here, t ≡ 1 + x + y for brevity.The factor b ∈ Z 2 x, y /(x L − 1, y L − 1) is a polynomial in x and y, subject to the constraint The form of b is further solved from ε 2 f = 0.For convenience, we introduce Clearly, α f = 0 is equivalent to ε 2 f = 0, provided ε 1 f = 0 is satisfied.Nevertheless, α is easier to handle as in Eq. ( 18) the dependence on z has been eliminated.We next treat b as a polynomial in y with coefficients in Z 2 x /(x L − 1).In order to find all the solutions of α f = 0, we temporarily lift the y L = 1 PBC.Then, α f = 0 reduces to αb = (c 1 y + c 0 ) (y L − 1).Thus, b can formally be expressed as a polynomial long division and the division should produce zero reminder.Here, the coefficients c 0 and c 1 ∈ Z 2 x /(x L − 1) are polynomials in x; two terms are needed because the degree of α in y is 2.
The polynomial long division in Eq. ( 19), where PBC y L = 1 is lifted, allows for the following physical interpretation.It ensures that X(b) (namely, the operation that flips the set of spins specified by b in the xy-plane) generates α-type excitations only at the two open boundaries y 0 and y L .Further, to ensure that excitations at both the y 0 end and the y L end get canceled when we reconnect the y-PBC, the polynomial long division has to produce zero remainder.The excitation patterns at both ends are described by c 1 y + c 0 .It uniquely determines the form of b by Eq. (19).
There is no simple explicit formula of b [61], and, in practice, it is more convenient to compute it on the fly.Specifically, one can pick up two independent polynomials in x as the c 0 and c 1 coefficients and examine if c 1 y+c 0 is divisible by α and if Eq. ( 17) is satisfied.In Sec.III C we discuss that the choices of c 0 and c 1 are straightforward in certain situations.

C. Fractal symmetry & ground state degeneracy
A most distinctive feature of the fractal Ising model is that its symmetry group strongly depends on the system size L. In particular, only specific Ls support fractal symmetries under PBC.As a consequence, the ground state degeneracy, which is determined by the number of spin-flip symmetries, i.e., GSD = |ker φ|, is also size-dependent.
The simplest choice of Ls leading to a fractal symmetry is L = 2 n , with n ∈ Z + .In such a case, Eq. ( 17) is automatically fulfilled since coefficients of f ∈ R are Z 2 valued, hence t L = (1 + x + y) L ≡ 1.The choices of c 0 and c 1 are then only subject to the requirement that (c 1 y + c 0 ) (y L −1) is divisible by α.We demonstrate in App.A that, for L = 2 n , this remaining constraint also reduces to a much simpler version: c 0 + c 1 containing an even number of terms, denoted |c 0 + c 1 | = 0 mod 2.
As an example, we can choose c 0 = c 1 = 1, which leads to a fractal pattern depicted in Fig. 2. Specifically, c 0 = c 1 = x 0 serves as a seed for Eq.(19).The resultant b represents the pattern of flipped spins in the xy-layer at z = 0.And other layers are specified accordingly by those higher degree terms in Eq. (16).
By taking other combinations of c 0 and c 1 , we can readily generate different fractal patterns and also global symmetries.For instance, c 0 = c 1 = L−1 i=0 x i leads to the obvious global symmetry which flips all spins in the entire lattice.
To determine the GSD, it is sufficient to know the number of spin-flip symmetries.For L = 2 n , it is simply the number of choices for c 0 and c 1 , under the constraint that |c 0 + c 1 | is even.Thus, the GSD for this class of system sizes is |ker φ| = 2 2L−1 .
When L 2 n , the fractal Ising Hamilton realizes different symmetry groups and degeneracies.It is not yet clear whether a simple expression exists to describe the ground state degeneracy for all values of L. Nevertheless, we are able to enumerate several situations that cover typical systems sizes We notice that the GSD in Eq. ( 20) equals the square root of the GSD for Haah's code [63].This is because the fractal Ising model is associated with one type of qubit error (either Pauli X or Z) in Haah's code.Hence, two copies of H FIM make up the degeneracy of Haah's code.In fact, a generalized version of this relation can be proved for general topological Calderbank-Shor-Steane (CSS) codes, which we will present in a separate work.

IV. PHASE TRANSITIONS AND SCALING
Having established the symmetries, order parameters, and degeneracies of the two fracton spin models, we next discuss their phase transitions and scaling behaviors.Numerical simulations are crucial as self-duality alone tells neither the number of phases nor the nature of their phase transitions.
We perform large-scale Monte Carlo simulations by jointly utilizing heat-bath updates and a multi-canonical method [67], and PBCs are considered.Both H TIM and H FIM turn out to experience a very strong first-order phase transition.The transi-  tion point in such situations is notoriously hard to locate precisely due to large hysteresis [68].We overcome this problem with multi-canonical simulations by learning a nearly flat histogram.Equilibrations are then ensured by the convergence of physical observables within statistical error bars.See App.B for algorithm details and simulation parameters.Phase transitions are determined by cross-checking the behaviors of the energy histogram P(E), specific heat C V , susceptibility χ, Binder cumulant B, where Here, P(E) measures the distribution of total energy E at a given temperature, O ∈ {Q z TIM , G FIM } denotes the subdimensional and fractal order parameters of the two models, and N is the number of spins in the system.The first-order nature of the phase transitions is concluded from diverging double peaks in histograms and negative dips in Binder cumulants (App.B).Finite-size transition points β c (L) are located by the extrema of C V , χ and B and equal-weighted peaks of P(E).
In a standard first-order phase transition, β c (L) is known to satisfy a scaling relation  This relation can be derived from a two-phase model by assuming the system has an Ω-degenerate ordered phase and a non-degenerate disordered phase, with an energy jump E d − E o ∼ L D between the two phases [69,70].At a firstorder phase transition, the ordered (W o ) and disordered (W d ) phases are expected to have the same weight, namely, Eq. ( 25) is obtained by a series expansion of ln W o W d around β ∞ c .However, in the current problem, the above scaling is subject to modification of subextensive degeneracies.It is first demonstrated in studies of a plaquette Ising model (PIM) and an anisotropically coupled Ashkin-Teller model (ACATM).Up to the leading order, the modification amounts to factoring out the size dependence in GSD [71,72].The consequent scaling thus turns into an anomalous one where b is a non-universal constant prefactor, and d reflects the power of log 2 GSD ∼ L d .Eq. ( 27) also indicates stronger finite-size effects in SDSB phase transitions, as the leading-order correction converges slower than in an usual first-order transition.Nevertheless, for H TIM and H FIM , once we know by numerics there is only a single phase transition, the thermochemical transition point β ∞ c = 1 2 ln 1 + √ 2 is fixed by self-duality.We are then left with a one-parameter fitting at the leading order.
We expect H TIM to display an 1 L 2 scaling as in the cases of PIM and ACATM [71], as these 3D models all have a planeflip symmetry.This is indeed confirmed by our simulations by simulating large systems up to L = 28.As shown in Fig. 3, the transition temperatures computed using distinct estimates all firmly collapse onto a 1 L 2 -line.It may be less obvious at first thought for H FIM .However, considering the degeneracies given in Eq. (20), this fractal symmetric model should also exhibit a 1  L 2 scaling in the sequence of system sizes L = 2 n (or L = 4 n − 1).Because of the constraints, only three lattice sizes L ∈ {8, 16, ∞} are available in the L = 2 n sequence; others are either too small for a meaningful fit or too large for practical simulations.In Fig. 4, we compare the 1  L 2 fit with a standard 1 L 3 fit and a reference 1 L fit.It is remarkable that such a minimal three-point fitting neatly identifies the expected non-standard scaling.
In the other non-trivial sequence L = 4 n − 1, the only accessible size is L = 15 as the next relevant one is already L = 63.Nonetheless, given the agreement between theory and simulations so far, we presume that this class also falls into a 1 L 2 scaling but converges to β ∞ c with a different slope due to the non-universal prefactor in Eq. (27).
The results of H FIM and H TIM indicate that, at the leading order, the anomalous first-order scaling depends only on the sub-extensive part in GSD, not specifically on their symmetry generators and global constraints.In fact, SDSB phase transitions without a fracton relevance can also show the same scaling behavior as in the case of a hybrid symmetry breaking [38].However, order parameters in these situations are intrinsically different as discussed in Sec.II.
We conclude this section by briefly commenting on the system size sequences L = 2 m−1 (2 n +1) and L = 2 m−1 (2 2n−1 −1) in Eq. (20).In these cases, if we keep m fixed, the fractal Ising model has distinct global-symmetry generators and constant degeneracies at different system sizes.A strong first-order phase transition is detected for all cases, but fitting the scaling is again very resource-consuming, as one has to group them according to their degeneracy classes and resort to large system sizes.However, given their constant degeneracies, one naturally expects them to fall into a conventional 1 L 3 scaling but converge to the thermodynamic limit with different slopes.

V. SUMMARY
Fracton systems provide novel schemes of fault-tolerant quantum computation and unconventional states of matter and call for new phenomenologies.In this work, we carried out a comprehensive study on two representative self-dual fracton spin models: the tetrahedral Ising model (TIM) and the fractal Ising model (FIM) that are the ungauging correspondence of checkerboard code and Haah's code, respectively.We constructed their order parameters and analyzed their ground-state degeneracies and phase transitions.
For the tetrahedral Ising model, its planar-flip symmetries lead to order parameters built from sub-dimensional linewise moments, as in Eq. ( 4).The long-range order of the system emerges from correlations of these extended objects, instead of local correlators.It represents a new type of order parameter, distinguished from those of Landau-types and Wilson loops.The construction may be generalized to general type-I fracton models that admit a foliation structure [25,26] and also non-fracton sub-dimensional symmetric models that can be partitioned into coupled lines or layers [38].
For the fractal Ising model, its fractal symmetry excludes the presence of local and semi-local ordering moments.Instead, the long-range symmetric correlator serves as an effective order parameter, as in Eq. ( 6).Namely, the system develops a long-range symmetry-broken order without ordering moments.This sharply contrasts with the lack of local order parameters in symmetry-unbroken topological phases and may be viewed as a characteristic property of fractal symmetry breaking.Furthermore, as there are no simple expressions of fractal generators, the utilization of new algebraic tools, such as polynomial rings, becomes necessary for investigating such systems.
The phase transitions of both models belong to the regime of sub-dimensional symmetry breaking (SDSB) that features a sub-extensive number of degenerate ground states and longrange orders with non-local order parameters.Such phase transitions appear to be commonly first-order in three dimensions.Aside from the two models studied here, other examples include the dual spin models of the X-cube code [71,72] and two sub-dimensional symmetric spin models without a fracton correspondence [38].They are also in line with the first-order quantum phase transitions of perturbed Xcube [34], checkerboard [73], and Haah's codes [35]; namely, all three representative fracton codes when their excitations are condensed in the simplest manner.Therefore, understanding the origin of these first-order transitions or finding exceptions will be an interesting exploration.
Despite being first-order, these phase transitions display an anomalous finite-size scaling, including the fractal Ising model (Sec.IV).The exponent of the scaling reflects the size dependence of the subextensive GSDs.Thus, this anomalous scaling represents a common feature for fracton spin models and other sub-dimensional symmetric models.
It would also be interesting to investigate the Z N forms of these fracton spin models.Such Z N generalizations are a canonical topic in studying models with global symmetries and local symmetries.In particular, N's value may affect the structure of the underlying phase diagram and the nature of the associate phase transitions.For instance, the N-state clock models [74,75] and Z N lattice gauge theories [76,77] in three dimensions have a single continuous phase transition for any finite N, which extrapolate to the 3D XY model and a confined U(1) gauge theory when N → ∞, respectively [64,78].However, the N-state clock models in two dimensions [79][80][81] and the Z N lattice gauge theories four dimensions [82,83] have an intermediate phase with algebraically decayed correlations for N ≳ 5, and consequently two separate phase transitions.The appearance of such intermediate phases is deeply related to the self-duality of 2D spin models and 4D gauge theories [80,83].One can analogously ask whether self-dual 3D fracton spin models, particularly the Z N generalizations of the tetrahedral and fractal Ising models, could also support an intermediate phase at some larger N.
The breaking of a sub-dimensional symmetry constitutes a different type of phase transition than breaking a global symmetry or a gauge symmetry.Our work advances the development of a complete understanding of SDSB phase transitions.Moreover, the order parameters and scaling established here can further incorporate effects of quenched disorders and provide a crucial guide to studying the error resilience of the checkerboard code and Haah's code.The project/research is part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus.
The data used in this work are available in Ref. [84].where in the last line we used the PBC z = z L−1 .
Notice the following identities Hence, for L = 2 n , Eq. (A5) reduces to which holds if and only if c 1 + c 0 is divisible by p, hence by x + 1 using the original variable x.
For simplification, the condition can be equivalently stated as requiring c 0 + c 1 to contain an even number of terms.This equivalence can be demonstrated by observing the dichotomy c 1 + c 0 ≡ 0 or 1 (mod x + 1) and its correspondence with whether c 1 + c 0 contains an even or odd number of terms.

Appendix B: Multi-canonical simulations
As shown in Fig. 5, both fracton spin models H TIM and H FIM showcase a strong first-order phase transition with a significant energy barrier ∆E = E d − E o between the disordered and ordered phase.Such an energy barrier strongly suppresses intermediate states and makes tunneling events between the two phases unlikely.As a result, canonical simulations will be trapped in one of the metastable states.To overcome this issue, we performed multicanonical (MC) Monte Carlo simulations that modify the canonical probability distribution and promote the exploration of intermediate states.We now summarize the main steps of the multicanonical algorithm and refer to Ref. [67] for a detailed introduction.
The multicanonical algorithm aims to learn a flat energy distribution in a sufficiently large temperature interval covering E d and E o .To do so, we define the following multicanonical partition function where ρ(E) is the density of states at energy E, and e −g(β,E) is an unknown weighting factor to be learnt.Here, H(E) = N R j=1 h j (E)/N R is the averaging histogram at a given energy E, and ⟨. . .⟩ E ′ denotes a thermal average over the energy interval of interest.
After each iteration, the new g(β, E) will suppress those frequently visited states but promote the less visited states in the previous Monte Carlo runs.Such processes are repeated until H(E) satisfies a flatness condition In practice, we consider H(E) is flat if |H(E) − ⟨H(E ′ )⟩ E ′ | /⟨H(E ′ )⟩ E ′ ≲ 0.1 for three consecutive iterations.
At the end of the above weight-learning procedure, we obtain g(β, E) that gives us a nearly flat energy distribution H MC (E) over [E o , E d ], as shown in Fig. 6.We can use the learned weight to adjust the Boltzmann factor at any nearby temperature β ′ by setting g(β ′ , E) = g(β, E) + (β − β ′ )E and perform new Monte Carlo simulations according to Eq. (B1) for measurement runs.In general, g(β  The global normalization factor is Z Z MC = E ′ H MC (E ′ ) e g(β,E ′ ) .In Table I

Figure 2 .
Figure 2. Fractal symmetry X ( f ) that corresponds to (c 0 , c 1 ) = (1, 1) in a periodic system of linear size L = 128.The pattern of flipped spins (black pixels), encoded by f , is illustrated by several layers in the z-direction.(a, b) 3D and 2D representations for the pattern of spin flips in the xy-layer at z = 0. (c-h) 2D representations for patterns of flipped spins in some other xy-layers; the presented xy-layers are located at z = 1 32 L, 1 8 L, 1 4 L, • • • , 7 8 L.
a binary function, we can rewrite the corresponding Boltzmann factors as ω β,n = e β cos nπ and ω ⋆ β ⋆ ,k = e β ⋆ cos kπ , with n, k = 0, 1.The original and dual Boltzmann factors are related by a discrete Fourier transform, ω β,n = 1 √ 2 with integers n, m ≥ 1.The first two classes in Eq. (20) are subextensive degeneracies due to different fractal symmetries; both have an exponent linear in L. The latter two classes offer various lattice sizes (by changing n) that support a constant degeneracy for each fixed m.For instance, L = 3, 5, 7, ... (m = 1) have GSD = 2 and only the obvious global symmetry.Meanwhile, L = 6, 10, 18, ... (m = 2) give GSD = 8, indicating the presence of sophisticated global symmetries.

Figure 3 .
Figure 3. Anomalous scaling for the tetrahedral Ising model.Simulations are performed on lattices with L = 16, 18, ..., 28 (N = 1 2 L 3 ).Finite-size transition temperatures β c (L) are determined from the extrema of specific heat C v , susceptibility χ, and Binder cumulant B and the equal-weight peaks of energy histogram P(E).The thermodynamical β ∞ c is fixed by self-duality.All fits collapse onto the expected 1 L 2 scaling.

Figure 4 .
Figure 4. Anomalous scaling for the fractal Ising model, with L = 8, 16, ∞.The minimal three-point fit clearly prefers the expected 1 L 2 scaling (middle), and excludes the a conventional 1 L 3 fit (left) and a reference 1 L fit (right).
Acknowledgments.G.C., K.L., and L.P. acknowledge support from FP7/ERC Consolidator Grant QSIMCORR, No. 771891, and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868.M.A.M.-D.acknowledges support from grants MINECO/FEDER Projects, PID2021-122547NB-I00 FIS2021, the "MADQuantumCM" project funded by Comunidad de Madrid.M.A.M.-D.has been financially supported by the Ministry of Economic Affairs and Digital Transformation of the Spanish Government through the QUANTUM ENIA project call-Quantum Spain project, and by the European Union through the Recovery, Transformation and Resilience Plan-NextGenerationEU within the framework of the Digital Spain 2026 Agenda.M.A.M.-D.has also been partially supported by the U.S.Army Research Office through Grant No. W911NF-14-1-0103.H.S. acknowledges support from the National Natural Science Foundation of China (Grant No. 12047503).K.L. acknowledges support from the New Cornerstone Science Foundation through the XPLORER PRIZE, Anhui Initiative in Quantum Information Technologies, and Shanghai Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01).

Figure 5 .
Figure 5. Energy histograms (left column) and Binder cumulants (right column) for TIM (first row) and FIM (second computed with multicanonical simulations in the proximity of the phase transitions.The behavior of the dips in the Binder cumulants and the double peaks in the energy histograms for growing lattice sizes confirm the presence of a first-order phase transition in both models.

Figure 6 .
Figure 6.Energy the estimated transition point obtained from the multicanonical simulations and the resulting canonical reweighted distribution for the Tetrahedral Ising model at L=20. ) It is not hard to see that Z MC (β) gives a flat energy distribution in the interested interval [E o , E d ] only if the weighting factors satisfying e g(β,E) = ρ(E)e −βE for all energies therein.This is achieved by initializing g(β, E) = 0 ∀ E and updating them iteratively by repeating the following steps.1. Run a set of Monte Carlo simulations in parallel at a given temperature β near the phase transition with N R different replicas (independent initializations).We use the standard heat-bath algorithm for Monte Carlo updates, but the Boltzmann weights are now evaluated according to Eq. (B1).Compute the energy histogram h j (E) for the j-th replica; namely, keep tracking the frequency of each energy E ∈ [E o , E d ].
2. Once a previous set of Monte Carlo runs finishes, update g(β, E) with an inversion rule g(β, E) → g(β, E) + ln H(E) − ln H(E ′ ) E ′ .(B2) ′ , E) does not have a perfect flat energy distribution for β ′ β, but it is enough to efficiently explore the intermediate states in the interested interval [E o , E d ].Physical observables are defined by canonical expectation values.Given the learned reweighting factors e −g(β,E) , we can

Table I .
, we summarize the simulation parameters.It is worth noting that, as multicanonical simulations explore a larger configuration space than a canonical one, larger amounts of Monte Carlo updates are required to obtain accurate estimates of physical expectation values.Simulation parameters for the weight-learning process and the subsequent multicanonical simulations.N R denotes the number of independent simulations from which histogram averages are taken in the weight learning process.N WL is the number of sweeps used for each replica in the last three iterations.Multicanonical simulations then run at N T different temperatures for N S number of Monte Carlo sweeps at each temperature.