Optimal Thresholds for Fracton Codes and Random Spin Models with Subsystem Symmetry

Fracton models provide examples of novel gapped quantum phases of matter that host intrinsically immobile excitations and therefore lie beyond the conventional notion of topological order. Here, we calculate optimal error thresholds for quantum error correcting codes based on fracton models. By mapping the error-correction process for bit-flip and phase-flip noises into novel statistical models with Ising variables and random multi-body couplings, we obtain models that exhibit an unconventional subsystem symmetry instead of a more usual global symmetry. We perform large-scale parallel tempering Monte Carlo simulations to obtain disorder-temperature phase diagrams, which are then used to predict optimal error thresholds for the corresponding fracton code. Remarkably, we found that the X-cube fracton code displays a minimum error threshold ($7.5\%$) that is much higher than 3D topological codes such as the toric code ($3.3\%$), or the color code ($1.9\%$). This result, together with the predicted absence of glass order at the Nishimori line, shows great potential for fracton phases to be used as quantum memory platforms.

Introduction. The study of quantum phases constitutes a cornerstone of quantum many-body physics and can potentially enable technological advances. Among its modern applications is the possible realization of a fully fledged quantum computer by means of fault-tolerant methods for processing quantum information [1][2][3]. Topological codes [4,5] stand among the best options to performing fault-tolerant quantum computation due to their high thresholds and linear scaling of the system qubit resources [6,7]. Nevertheless, 2D topological stabilizer codes like the most studied surface code [4,6,8] permit only topological implementations of Clifford gates [9], while non-Clifford gates are necessary for realizing the desired quantum advantages [10], motivating the quest for new 3D codes. Fracton models [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] represent a generalization to 3D topological orders and provide alternatives to quantum memories beyond the standard paradigm of topological computing. These models host intrinsically immobile point-like excitations called fractons [15] which make a key difference from conventional topological orders and have potential beneficial applications. While a few decoders [28,29] and several experimental platforms [30][31][32] are proposed, the theoretical limit on error thresholds of fracton codes is unexplored, which is nevertheless crucial for devising new decoders and for justifying the practical relevance of fracton codes.
The goal of this Letter is to investigate how a fracton model behaves as an active error correcting code against stochastic Pauli errors which are widely used for benchmarking quantum memories. We have defined error corrections in the presence of a subextensive ground state degeneracy and computed the optimal thresholds for one of the most representative fracton models in three dimensions-the Xcube model [16]. We address the problem by a combina-tion of theoretical analyses and numerical simulations. Using a statistical-mechanical mapping method [6] that has previously produced error thresholds for codes beyond those for which it was initially conceived [33][34][35][36][37][38][39][40], we derive two statistical models related to Pauli errors of the X-cube model, in the formulation of classical spin variables that are suited for simulations. The numerical simulation of statistical models in three dimensions with randomness is generally challenging, and the required resources are even higher for our models as they possess lower-dimensional subsystem symmetries rather than a more conventional global symmetry. Only by utilizing state-of-the-art parallel tempering Monte Carlo methods and performing large-scale simulations have we been able to compute various many-body correlation functions and determine the phase diagrams for the two statistical models up to moderate error rates.
We estimate the optimal error thresholds of the X-cube code against bit-flip (X) and phase-flip errors (Z) as p X c 0.152 (4) and p Z c 0.075 (2), respectively. The minimum of these thresholds is remarkably higher than what was found in conventional 3D topological codes such as the toric code (0.033) [41,42] and the color code (0.019) [38], which signals the potential of the X-cube model as a fault-tolerant quantum memory. This is further confirmed by the analytical result that the Nishimori line is free of fracton glass order through which the resilience of the quantum code may be lost (see SM [43]). In addition, our results represent the first study of spin models with both subsystem symmetries and quenched random disorder in three dimensions, hence are also of interest for the statistical mechanics community. X-cube model as quantum memory. Consider a cubic lattice L with periodic boundary condition (PBC). We introduce Figure 1. Stabilizer generators, excitations, and logical operators of the X-cube model. A physical qubit is assigned to each edge ( ) of the lattice. A Pauli X (Z) operator is presented as a cyan (red) edge. Stabilizer generators A c and B µ v are defined at unit cubes (c) and vertices (v), respectively. The ground states satisfy A c = 1 and B µ v = 1 for all c, v and µ ∈ {x, y, z}. Cyan strings (S) and gray membranes (M) represent string and membrane operators of the form X S ∈S X and Z M ∩M ∅ Z , respectively. Fracton excitations A c = −1 are created by membrane operators at their corners and colored as cyan cubes. String operators create lineon excitations at their turns and ends. The red ellipsoids represent lineons with where µ is indicated by the elongated directions. The extended strings and membranes S 1 , S 2 , M 1 , and M 2 represent two pairs of logical operators creating no excitations due to PBC. a qubit to every edge and define stabilizer generators A c and B µ v at each unit cube c and vertex v of the lattice. Specifically, A c is defined to be the tensor product of Pauli X operators on the 12 edges of a cube, and B µ v is a tensor product of Pauli Z operators on the four edges adjacent to a vertex and perpendicular to the spacial direction µ = x, y, z. Namely, as visualized in Fig. 1. For convenience, we label the set of qubits as Q and those of stabilizer generators as A {A c } and B {B x v , B y v }, respectively. The X-cube model is a paradigmatic fracton model constructed by summing over these stabilizer generators, As A c and B µ v commute, this Hamiltonian is exactly solvable [16], and its ground states satisfy A c = 1 and B µ v = 1 for all c, v, µ. The elementary excitations are two types of gapped topological defects. A unit cube with A c = −1 is referred to as a fracton (solid cyan cube in Fig. 1), which is an intrinsi- corresponds to an excitation termed a lineon (red ellipsoid in Fig. 1), which can move along the µ direction but is immobile in the two directions ν µ. The three possible lineons at each vertex are subject to a constraint B x v B y v B z v = 1. On a lattice of size L 3 with PBC, the X-cube model has 2 6L−3 degenerated ground states, scaling subextensively with system size [16,20]. These ground states are indistinguishable by local operations, hence provide a fault-tolerant code Hilbert space. We can view them as 6L − 3 logical quibits by introducing 6L − 3 pairs of non-local operators (X S j , Z M j ).
Here, X S j ∈S j X and Z M j ∩M j ∅ Z are defined on extended strings and membranes winding around the lattice (see Fig. 1).
Error correction. Fractons and lineons can be used to diagnose errors in the stabilizer code. For example, as illustrated in Fig. 1, a phase-flip Z error on a single qubit will cause four fractons at each of its adjacent cubes. Similarly, a single bitflip X error will create two lineons at the vertices sharing the edge. Therefore, an ensemble of fractons or lineons can act as an A-or B-syndrome reflecting Z or X errors in the system.
For simplicity, we consider a situation where each qubit is affected by phase-flip and bit-flip errors independently and assume perfect measurements for all stabilizer generators. Moreover, since the X-cube model is a Calderbank-Shor-Steane (CSS) code [44], i.e., the type-A and type-B stabilizer generators involve either purely X or Z operators, we can correct the bit-flip and phase errors separately.
The error-correction process can be described by introducing a (co)chain complex, where Z A 2 , Z Q 2 and Z B 2 denote the Z 2 = {0, 1} vector spaces for labeling configurations of type-A stabilizer generators (A), physical qubits (Q), and type-B stabilizer generators (B), respectively. The boundary maps ∂ A and ∂ B are linear and specify the qubits involved in every type-A and type-B stabilizer generator. Correspondingly, the transpose operator ∂ † A (∂ † B ) maps an error configuration η ∈ Z A(B) 2 to an ensemble of frac- for all CSS codes. Only certain error configurations are compatible with a given syndrome. Among those, error configurations are equivalent if and only if they can be connected by type-A and type-B stabilizer generators. Namely, provided η − η ∈ im ∂ A or im ∂ B , two errors η and η will have the same effect on the encoded quantum state, where im ∂ denotes the image of the boundary map. Thus, the spaces of X η and Z η can be divided into equivalence classes by the quotients Z Q 2 /im ∂ A and Z Q 2 /im ∂ B . We denote the equivalence classes as η X η + im ∂ A and η Z η + im ∂ B , respectively.
For a possible B-or A-syndrome σ with probability Pr(σ), the total probability of those equivalence classes compatible with σ satisfies λ pr η σ + λ = Pr(σ), with λ ∈ Z Q 2 labeling inequivalent logical operators. The correction can be real- ized successfully if, for typical syndromes, there exists a most probable equivalence class η * σ such that pr η * σ → Pr(σ) in the large system limit [6]. However, this is only possible when the rate p for local X and Z errors lie below some optimal threshold values p X c and p Z c . For p > p c , the error class cannot be unambiguously identified, and the code becomes ineffective. Finding the optimal error thresholds is therefore crucial to any quantum code.
Mapping to statistical-mechanical models. An elegant and numerically preferable way to determine p X c and p Z c of the X-cube code is utilizing a statistical mapping method [6] which maps bit-and phase-flip errors to suitably chosen statistical-mechanical models.
Suppose both X and Z errors occur independently at each qubit at rate p. Then the probability the system is affected by an X or Z error configuration η ∈ Z Q 2 is where η( ) = 1 or 0 on edges with or without an error.
This probability can be interpreted as a Boltzmann weight by introducing an effective temperature T satisfying Eq. (5) defines the so-called Nishimori line and allows us to control the rate of random qubit errors through the auxiliary temperature (see SM [43]). Accordingly, the total probability of a bit-flip error equivalence class η X η+im ∂ A is mapped to the partition function of an interacting spin Hamiltonian H A η , where β = 1/T is the inverse temperature, f ≡ { f (c)} c∈A ∈ Z A 2 represents a configuration of type-A stabilizer generators, ∂ A f ( ) = c∈∂ † A f (c) labels the edges of f , and S c = (−1) f (c) ∈ {±1} denotes effective Ising variables on the center of cubes.
The form of H A η realizes a 3D random plaquette Ising (RPI) model on a dual lattice with quenched disorder, where c ∈ ∂ † A specifies the four type-A stabilizer generators sharing an edge (see Fig. 2), while the probabilities of their coupling to be anti-ferromagnetic and ferromagnetic are p and 1 − p, respectively. Moreover, aside from a global Z 2 symmetry, H A η is invariant under a subsystem symmetry flipping spins in individual planes and may be viewed as a novel random spin model.
Analogously, the modeling of a phase-flip error equivalence class η Z η + im ∂ B leads to a 3D random anisotropically coupled Ashkin-Teller (RACAT) model with quenched disorder on the original cubic lattice, where 2 denotes the indicator vectors for type-B stabilizer generators. The spins are coupled only along with the unitμ direction (Fig. 2). In contrast to the usual 3D Ashkin-Teller model [45], the RACAT model Eq. (8) has the planar symmetries of flipping all S µ and S ν spins in an arbitrary µ-ν plane besides a global Z 2 × Z 2 symmetry.
The disorder-free limits (p = 0) of H A η and H B η are dual to each other [46], as for general fracton and topological CSS codes [43]. There is no exact duality in the presence of disorder, nevertheless, our results suggest an approximate duality relation between the error thresholds p X c and p Z c [43]. On the side of the statistical-mechanical models, the relative probability between two X (or Z) error equivalence classes under the error rate p is given by the difference between their free energies, where λ ∈ Z Q 2 represents logical X (or Z) operators of the Xcube model and can flip a sequence of coupling coefficients in H A η (H B η ). The condition of existing pr η * σ X(Z) → Pr(σ) for the most probable equivalence class η * σ X(Z) requires that the free energy to introduce a nontrivial string (membrane) defect (λ 0) diverges in the thermodynamical limit, namely, δF A(B) η,λ → ∞ (see SM [43]). This is only possible when H A η and H B η are in their ordered phases. Hence p X c and p Z c can be determined from the order-disorder phase transitions of the two random spin models.   (2). The error thresholds are determined by the largest p values exhibiting an order-disorder phase transition (see SM [43]). The effective temperature T is an auxiliary quantity and related to p by the Nishimori line (dotted line). A correctable X-cube code corresponds to the part of the Nishimori line inside the ordered phases.  . For a second-order phase transition the curves for different system sizes intersect near the critical temperature T c . The shaded areas show a conservative estimate of T c . Such intersections cannot be recognized for p > p X c or p Z c where we conclude no phase transition.
Error thresholds and phase diagrams. The phase diagrams of the RPI and RACAT models are shown in Fig. 3, obtained by large-scale parallel tempering Monte Carlo simulations. We locate the phase transitions by cross-checking the energy histogram, specific heat, order parameter, its susceptibility, and the correlation length [43].
To construct the appropriate order parameters, the planar symmetries of H A η and H B η have to be taken into account as they can lead to trivial cancellation of local orders. For the RPI model, we define with . and [.] denoting the thermal and disorder average, respectively. The inner sum in Eq. (10) involves a subextensive number (∝ L 2 ) of spins, while the norm enforces the planar-flip invariance [47]. Thus, Q A defines a long-range order which is sub-dimensional and made of plane-like objects. The order parameter for the RACAT model is constructed similarly which describes an order for extended line objects. Here, the S z v spins are taken for simplicity, as the three Ising variables in Eq. (8) are permutable.
At low X and Z error rates, the energy histograms reveal a first-order phase transition for both H A η and H B η [43]. This agrees with previous studies on the disorder-free (p = 0) limit of the two models [46,47]. Hence the transition temperatures can be estimated in a relatively straightforward way [43].
For larger p, the phase transitions are softened to continuous ones in line with the Imry-Ma scenario [48,49]. We can then locate the transitions by studying the second-moment correlation length whereG (k) r G (r) e −ik·r is a Fourier transform of the spatial correlator G (r), and k min denotes any smallest nonzero wave vector [50].
The spatial correlators related to the order parameters Q A and Q B are given by so that G A(B) (r) → Q A(B) 2 in the limit of |r| → ∞.
Provided a continuous phase transition exists, ξ L /L is expected to be scaled as g L 1/ν (T − T c ) near the critical point, and the curves for different system sizes should intersect at (T c , g (0)), where g is a universal scaling function and ν denotes the critical exponent of ξ. We can use this property to search the largest error rates where H A η and H B η continue showing a continuous phase transition, which in turn implies the error thresholds of the X-cube code.
Our simulations show that the ξ L (T )/L curves exhibit clear intersections up to an error rate p X c 0.152(4) for the RPI model and p Z c 0.075(2) for the RACAT model (see Fig. 4 and also SM [43]). Thereafter, a clear intersection cannot be recognized [43], indicating lack of an order-disorder transition. Namely, for error rates larger than p X c and p Z c , H A η and H B η host no long-range order, and the X-cube code hence becomes uncorrectable.
Conclusions. The X-cube model is the archetypal stabilizer code exhibiting the fascinating quantum physics of fracton topological orders. In this work we investigated its capability as a quantum memory through a combination of theoretical analyses and detailed numerical simulations. We estimated its optimal error thresholds as p X c 15.2(4)% against bit-flip noise and p Z c 7.5(2)% against phase-flip noise, featuring a remarkably higher minimum error rate (7.5%) compared to the 3D toric code (3.3%) [41,42] and color code (1.9%) [38]. Our work establishes the general connection between the fault tolerance of fracton codes and statisticalmechanical models with subsystem symmetries. The Pauli error thresholds in any CSS code with zero-encoding rate [6,33,34,38] obey the inequality H p X c + H p Z c ≤ 1 imposed by the quantum Gilbert-Varshamov bound [51][52][53], where H (p) −p log 2 (p) − (1 − p) log 2 (1 − p) is the Shannon entropy. Our results give H p X c + H p Z c 1.00 (2), which not only satisfies this constraint but is close to its upper bound, similar to the situations found in conventional topological codes [33,38,54]. We formulate this near saturation via an approximate duality in the SM [43] and conjecture it for general fracton and topological CSS codes. Our work can guide further studies of fracton models, and the approximate duality predicts even higher thresholds for the self-dual checkerboard [16] and Haah's [13] codes.  Here we show that the 3D random plaquette Ising (RPI) model H A η and the 3D random anisotropically coupled Ashkin-Teller (RACAT) model H B η can be related by a Kramers-Wannier duality. Our discussion is not restricted to specific models and can apply to the analysis of general fracton and topological Calderbank-Shor-Steane (CSS) codes.

A. Exact Kramers-Wannier duality for disorder-free models
The duality between H A η and H B η is exact in the disorder-free (p = 0) limit. In this limit, the relevant error equivalence classes are the trivial ones (η ≡ 0), and the partition function Z A η modeling X-errors reduces to where f ≡ { f (c)} c∈A ∈ Z A 2 labels the configuration of type-A stabilizer generators, ∂ A f ∈ Z Q 2 specifies the corresponding qubit configuration, and W β (ξ) e β(−1) ξ( ) denotes the Boltzmann weight for a general qubit configuration ξ ∈ Z Q 2 . The Kramers-Wannier duality can be viewed as a Fourier transform [59]. The dual of W β can be expressed as in terms of a dual inverse temperatureβ specified by where ζ ∈ Z Q 2 is the conjugate variable of ξ, and ζ, ξ ∈Q ζ ( ) ξ ( ) denotes an inner product. The Fourier transform rewrites Z A 0 as Using the identities ζ, one finds that only those ζ ∈ ker ∂ † A contribute. Moreover, as in the thermodynamical limit the free energy density is independent of the boundary conditions, we can choose an open boundary condition such that ker with g labelling configurations of type-B stabilizer generators, and ∂ B g = ζ ∈ im ∂ B . Therefore, the Kramers-Wannier duality for CSS codes in the disorder-free limit can be established as As the X-cube model is a CSS code, the duality between the disorder-free limit of the RPI model and the RACAT model, namely, H A η=0 and H B η=0 , then follows immediately.
B. Approximate duality between the optimal bit-flip and phase-flip error thresholds In the presence of disorder (p > 0), there is no exact duality between the pair of error models H A η and H B η . Nevertheless, the current work for error models with subsystem symmetries and the previous studies for models with global or local symmetries [33,38,54] suggest that the optimal error thresholds p X c and p Z c satisfy an approximate duality relation where H (p) −p log 2 (p) − (1 − p) log 2 (1 − p) is the Shannon entropy. Below, we show that this approximate duality may be understood by a replica analysis. Consider a disorder average of the partition function Z A η over n replicas, where f j ∈ Z A 2 labels configurations of type-A stabilizer generators for the j-th replica, and f . At the error rate p, the coupling coefficient (−1) η( ) equals to ±1 with the probability 1 − p and p, respectively. Thus, the disorderaveraged Boltzmann factor associated with each edge can be expressed as and |h| counts the number of the h j = 1 components. As in the disorder-free case, we can define a dual Boltzmann factor with t (t 1 , t 2 , · · · , t n ) ∈ Z n 2 being the conjugate variable of h, and |t| denoting the number of t j = 1. Analogous to Eq. (S5), the disorder-averaged partition functions Z A n and Z B n can be related as where Z A n and Z B n are viewed as functions of n + 1 Boltzmann factors labelled by |h| or |t| = 0, 1, · · · , n. With the principle factors w P (p, β) max w p,β (|h|) = w p,β (0) and w P (p, β) max w p,β (|t|) = w p,β (0) factored out, Eq. (S11) becomes where u(p, β) 1 w P (p,β) (w p,β (0), w p,β (1), · · · , w p,β (n)) and u(p, β) 1 w P (p,β) ( w p,β (0), w p,β (1), · · · , w p,β (n)) are the normalized Boltzmann factors.
Analogously, we also have Assume that, for both H A η and H B η , there are only one ordered phase and one disordered phase separated by a single phase transition at p X c , β p X c and p Z c , β p Z c , respectively, on the Nishimori line β = β(p), where The phase transition of Z A η (Z B η ) will then occur along the path of normalized Boltzmann factors u(p, β(p)) at p X c (p Z c ), and also along the dual path u(p, β(p)) at p Z c (p X c ), respectively, using the relations in Eqs. (S13) and (S12). As the phase transition is unique, one may expect u p X c , β p X c ≈ u p Z c , β p Z c and u p Z c , β p Z c ≈ u p X c , β p X c . Hence, by multiplying Eq. (S12) and Eq. (S13), and using the identity 2 |A|+|B|−|Q| = |ker ∂ A | |ker ∂ B |, we obtain the following relation between the principle Boltzmann factors w P and w P , Moreover, given the expressions of the Boltzmann factors in Eq. (S9) and Eq. (S10), one has ln w P (p, β (p)) w P (p, β (p)) = ln w p,β(p) (0) w p,β(p) (0) = ln 2 n 2 (1 − p) n+1 + p n+1 = n ln 2 1 2 − H (p) + O n 2 (S16) in the replica limit n → 0 and along the Nishimori line Eq. (S14). Therefore, the approximate duality relation Eq. (S7) can be established, namely, (S17)

S.II. ERROR PROBABILITY AND CORRECTION
To be able to recover encoded quantum information from qubit errors, one needs to identify the error equivalence class unambiguously. For clarity, we take the detection of bit-flip (X) errors as an example, whereas phase-flips (Z) errors can be analyzed in the same way.
Let ρ = |ψ ψ| be the density matrix of a generic state in the X-cube code space and pr(η) denote the probability of an X error configuration η. The error affected state ρ 1 is given by where in the second line we have grouped error configurations into equivalence class η X η + im ∂ A . We relabel error equivalence classes by syndromes and logical operators. This has a practical relevance: Error syndromes are local measurements of stabilizer generators and cannot distinguish topological operators, such as X logical operators which are strings winding around the lattice. Namely, a set of distinct X error equivalence classes η σ + λ X η σ + λ + im ∂ A will lead to the identical syndrome σ, where λ labels inequivalent logical operators. Such a surjective relation between error equivalence classes and error syndromes roots in the nature of a topological code. Hence, ρ 1 can be expressed as Without loss of generality, we use η * σ X with λ = 0 to represent the most probable error class, such that pr( η * σ X ) ≥ pr( η * σ + λ X ), ∀λ. Then, if we clean the error syndrome by X η * σ , the quantum state becomes To ensure the recovery of the initial state, namely, ρ 2 → ρ, we shall have pr η * σ + λ X /pr η * σ X → 0 for λ 0 and large enough code blocks. This is indeed possible if the code is in its ordered phase and, given our established mapping between the X-cube code and statistical-mechanical models, it can be demonstrated by a procedure originally developed for surface codes [6].
It is intuitive to rewrite the Nishimori line in Eq. (S14) as The probability of each error configuration η compatible with σ then resembles a Boltzmann weight, Accordingly, the relative probability between two X error equivalence classes relates to the partition function of the RPI model as (reproduced from the main text for convenience) Therefore, the relative probability can be estimated by the free energy cost δF A η * σ ,λ 0 of a topological defect (generalization of domain walls) in H A η . In the ordered phase of H A η , such defects are suppressed exponentially. In the disordered phase, topological defects costs zero energy, hence all error classes η σ + λ X have equal probability and a most probable error class does not exist.

S.III. ABSENCE OF GLASS ORDER ALONG NISHIMORI LINE
Whether there exists a spin glass order can crucially affect the error thresholds of a CSS code. In particular, in establishing the approximate duality relation Eq. (S7), the assumption of a unique phase transition also implies spin glass phases are irrelevant. Below, we prove that there is indeed no spin glass order along Nishimori line neither for the RPI model nor the RACAT model, by generalizing Nishimori's argument which was originally conceived for the Edwards-Anderson model [60].
Moreover, as p ∝ e β(p) and 1 − p ∝ e −β(p) in terms of the auxiliary inverse temperature β(p) − 1 2 ln p 1−p , one has σ pr (Γ (σ) τ; p) σ ξ σ pr (Γ (σ) τ; p) By Eq. (S29) and Eq. (S30), the disorder average [ S ξ τ,β ] p becomes Further, using Eq. (S27) again, and noticing σ 2 ξ = 1, and reindexing Γ(σ)τ → τ, one realizes an identity In particular, on the Nishimori line β = β(p), we have which indicates an equivalence between a normal order paramter and a spin-glass (SG) order parameter. Specifically, if we choose S ξ to be the RPI correlation function, the identity Eq. (S33) becomes G A (r) = G A SG (r), with In the ordered phase of the RPI model, both G A (r) and G A SG (r) will develop a finite value in the limit |r| → ∞, while they vanish in the disordered phase. A spin glass phase would require lim |r|→∞ G A (r) = 0 but lim |r|→∞ G A SG (r) 0, which is nevertheless precluded by Eq. (S33).
Similarly, for the RACAT correlator, we can also derive an identity G B (r) = G B SG (r), where Hence, a spin glass order is also absent in the RACAT model along the Nishimori line.

S.IV. DETECTING PHASE TRANSITIONS
The stability of the X-cube code against local X and Z errors can be related to the order-disorder phase transitions of the RPI and RACAT model, respectively. Here we describe in detail how these phase transitions and the optimal error thresholds p X c and p Z c are determined.  Figure S2. The histograms P(E) for the 3D RACAT model at different error rates. As in the case of the RPI model, the system experiences discontinuous phase transitions in the low p regime as shown in (a) and (b). The quenched disorder weakens the discontinuity as the distance ∆E between the two P(E) peaks shrinks upon increasing p. The non-divergent peaks in (c) suggest finite-size effects of a second-order phase transition.

A. Observables and order of phase transitions
For both models and each p value, we compute the energy histogram P(E), the specific heat C v , the susceptibility χ, and the second-moment correlation length ξ L , ξ L (β, L) = 1 2 sin (|k min |/2) Here E is the energy density of a state, Q measures the value of the order parameters Q A or Q B , and . and [.] represent the thermal and disorder average, respectively. The definition of ξ L is reproduced from the main text for convenience. We find the order of the phase transitions by examining the behaviors of the energy histogram. For continuous phase transitions, P(E) typically shows a single peak for all temperatures. When a phase transition is discontinuous, P(E) features a double-peak structure on large enough system sizes and at temperatures near the transition, reflecting the phase separation of two metastable phases [61]. In addition, with increasing system sizes, the double peaks shall grow sharper and evolve towards separated δ functions. This is what we observed in the low p regions of the RPI model H A η and the RACAT model H B η , as shown in Figure S1 and Figure S2.
The distance ∆E between the double P(E) peaks reflects the latent heat [62], which shrinks upon increasing p as the quenched disorder weakens first-order phase transitions. Close to the error thresholds p X c 0.152(4) and p Z c 0.075(2), although P(E) remains showing two peaks [ Figure S1(c) and Figure S2(c)], the weight of the valley between them does not evolve towards zero when increasing L. This implies the double peaks will not evolve into separate δ-functions in the infinite size limit, in contrast to the case of a first-order phase transition. Such a non-diverging behavior is also observed as a finite-size effect in simulations of the 2D 4-state Potts model [63], where the phase transition is analytically known to be continuous [64].
Therefore, these data suggest that, in the finite but low p regions, the phase transitions of H A η and H B η are discontinuous similar as in the disorder-free limits [47,65]. Nevertheless, in both cases, the first-order transition line terminates at a critical point, where the p value coincides with, or is situated slightly before, the corresponding threshold error rate p c .

B. Non-standard first-order scaling
For a typical first-order phase transition, such as that in the 2D q-state Potts model with q ≥ 5, the finite-size transition temperature T c (L) is shifted by an amount L −d from T c = T c (∞) [62,66], where d is the dimension of the system. However, it has been understood recently that, in the presence of planar spin-flip symmetries such as in the disorder-free models H A η=0 and H B η=0 , the leading order finite-size correction to T c is modified to L −(d−1) due to the sub-extensive degeneracy (∼ 2 3L ) [47,65]. As the quenched disorder does not affect the subsystem symmetries of H A η and H B η , and their first-order phase transitions are robust against finite p (see Fig. S1 and Fig. S2), we expect that the modified scaling relation still holds, at least for the low p regions, and we extrapolate T c by fitting With increasing p, the first-order phase transitions gradually soften and the correlation length will eventually exceed our simulation system sizes. Nevertheless, as shown in Fig. S3, for all p values showing a first-order phase transition, the estimated T c is sufficiently above the Nishimori line which defines the relevant (p, T ) values for an error code. Therefore, while the precision on the transition temperatures is expected to be improved from larger system sizes (which are beyond reach however), this should only have minimal consequences in determining the error thresholds p X c and p Z c (see S.IV C).

C. Optimal error-threshold values
In the modeling of local errors, temperature is an auxiliary variable introduced by the constraint Eq. (S22). Only (p, T ) values on the Nishimori line are relevant for quantum error correction. Namely, a correctable X-cube code corresponds to the part of the Nishimori line inside the ordered phase of the RPI or RACAT model. Clearly, the high-temperature phase of both models is trivially disordered. We can thereby estimate the optimal error-threshold values p X c and p Z c by the largest error rates exhibiting an order-disorder phase transition.
As the phase transitions are continuous for large enough p values, we use ξ L as an estimator. Following the scale invariance, in the vicinity of a critical point ξ L scales as and ξ L /L for different sizes is expected to intersect near T c , where g is a universal scaling function and ν is the critical exponent of correlation length. We show in Figures S4 the curves of ξ L /L in the vicinity of the estimated optimal error thresholds. A clear intersection can be observed at p X c 0.152(4) for the RPI model and at p Z c 0.075(2) for the RACAT model, indicating the existence of a second-order phase transition. For p > p X c or p Z c , the curves do not meet or their intersection becomes very ambiguous, and we conclude no order-disorder phase transition.

S.V. DETAILS OF NUMERICAL SIMULATIONS
Simulating systems with quenched disorder and first-order phase transitions is generally a challenging task. In the simulations of the 3D RPI and RACAT model, we employ parallel tempering (PT) jointly with the heat bath and over-relaxation algorithms to equilibrate the systems [50]. The distribution of temperatures is carefully chosen and tested to ensure the acceptance ratio of PT updates [67] for each of the disorder strengths p and system sizes L. Large numbers (N d ) of random coupling configurations are considered, with N d = 200 in the low p regimes and N d = 800, 1600 for p values near the error thresholds. Statistical error bars are estimated by the bootstrap method [68]. Equilibration is tested by a binning analysis. For each simulation, the measurements are averaged over the Monte Carlo time interval [2 τ−1 , 2 τ −1], where τ labels the bins [34]. The system is considered equilibrated when, at least, the last three bins agree within statistical uncertainty. Simulation parameters are summarized in Table S1 and  Table S2. The simulations ran at the CoolMUC-2 cluster and the KCS cluster at Leibniz-Rechenzentrum (LRZ) and used over three million CPU hours without taking into account the intensive tests for optimizing temperature distributions in parallel tempering.  Table S1. Simulation parameters for the 3D RPI model. L is the linear size of the system. N d denotes the number of random coupling configurations at the error rate p. Larger N d is considered near the X error threshold p X c 0.152(4). 2 τmax gives the number of Monte Carlo sweeps in a simulation. N T temperatures between T min and T max are simulated in parallel. The same conventions are used for the 3D RACAT model in Table S2.  Table S2. Simulation parameters for the 3D RACAT model.