Non-Uniform Noise Rates and Griffiths Phases in Topological Quantum Error Correction
Adithya Sriram Nicholas O'Dea Yaodong Li Tibor Rakovszky Vedika Khemani


APS/123-QED




Department of Physics, Stanford University, Stanford, CA 94305



§ ABSTRACT 

The performance of quantum error correcting (QEC) codes are often studied under the assumption of spatio-temporally uniform error rates. On the other hand, experimental implementations almost always produce heterogeneous error rates, in either space or time, as a result of effects such as imperfect fabrication and/or cosmic rays.
It is therefore important to understand if and how their presence can affect the performance of QEC in qualitative ways.
In this work, we study effects of non-uniform error rates in the representative examples of the 1D repetition code and the 2D toric code, focusing on when they have extended spatio-temporal correlations; these may arise, for instance, from rare events (such as cosmic rays) that temporarily elevate error rates over the entire code patch.
These effects can be described in the corresponding statistical mechanics models for decoding, where long-range correlations in the error rates lead to extended rare regions of weaker coupling.
For the 1D repetition code where the rare regions are linear,  
we find two distinct decodable phases: a conventional ordered phase in which logical failure rates decay exponentially with the code distance, and a rare-region dominated Griffiths phase in which failure rates are parametrically larger and decay as a stretched exponential. 
In particular, the latter phase is present when the error rates in the rare regions are above the bulk threshold.
For the 2D toric code where the rare regions are planar, we find no decodable Griffiths phase: rare events which boost error rates above the bulk threshold lead to an asymptotic loss of threshold and failure to decode. 
Unpacking the failure mechanism implies that techniques for suppressing extended sequences of repeated rare events (which, without intervention, will be statistically present with high probability) will be crucial for QEC with the toric code. 










Non-Uniform Noise Rates and Griffiths Phases
 in Topological Quantum Error Correction
    Adithya Sriram, Nicholas O'Dea, Yaodong Li, Tibor Rakovszky, Vedika Khemani
    September 5, 2024
=====================================================================================



 

§ INTRODUCTION


A quantum error correcting (QEC) code aims to protect quantum information against the effects of environmental decoherence by robustly and redundantly encoding logical qubits into entangled states of many physical qubits.
The threshold theorem of fault-tolerant quantum computing allows for arbitrarily long computations to be performed on the logical qubits with arbitrarily high accuracy; i.e., the logical error rate goes to zero as the number of physical qubits is scaled up, provided that noise rates are below the error threshold, whose value depends on both the code being used and the types of errors that occur. <cit.>.

Accurately estimating the threshold and the logical failure rates is of great importance for any hardware implementation.
While this often requires taking into consideration various details (including a detailed description of the error model), it is believed that much of the qualitative aspects of QEC can be understood within simplified toy models, involving only phenomenological noise.
Within this setting, mappings to statistical mechanics models of decoding <cit.> have been particularly useful, in informing both the theoretical upper bound on the threshold, and the asymptotic scaling of the sub-threshold logical failure rate.

Phenomenological noise is often studied under the simplifying assumption of uniform error rates on all qubits for all times, while realistic implementations almost always display some sort of heterogeneity, resulting in e.g. non-uniform noise rates.
It is conceivable that when the fluctuations are short-range correlated, the simplifying assumption is still valid.
On the other hand, with extended spatio-temporal correlations — which naturally occur due to effects ranging from fabrication errors <cit.> to stochastic events such as cosmic rays striking the quantum hardware <cit.> — 
it is unclear if and how things will be qualitatively different.
Studies of such effects are therefore necessitated for a qualitative understanding of code performance in practice <cit.>. 

In this work, we investigate said effects of non-uniform spatio-temporally correlated noise rates 
which map to extended, sub-dimensional rare regions in the corresponding stat mech models.
We show that such rare regions, despite being sub-dimensional and rare, can often have a disproportionately large effect and dominate the logical failure of the code, a phenomenon known as a Griffiths effect.
Griffiths effects have been extensively studied in disordered systems, and correspond to rare but large disorder fluctuations which dramatically change the properties of phase transitions (and the proximate phases) <cit.>.


We mainly focus on rare events (such as cosmic rays) that temporally increase the error rate from p_ bulk to p_ rare > p_ bulk over the entire code patch.
We study this for the representative examples of the 1D repetition code and the 2D toric code, in an phenomenological error model with both measurement and qubit errors.
They are described by the 2D random bond Ising model (RBIM) and the 3D random plaquette Ising gauge theory (RPGT), respectively, both with non-uniform coupling strengths, see Sec. <ref>

In the case of the the 1D repetition code, we leverage known results from disorder physics to argue that the rare regions lead to a new Griffiths phase above the conventional decodable phase, where the logical failure rate decays as a stretched exponential of the code distance, instead of the exponential scaling that occurs in the conventional decodable phase (when rare regions are absent).
Therefore, while the rare regions do not lead to a failure of error correction (there is still a finite threshold), they can qualitatively change the code's performance by parametrically increasing the logical failure rate.  
Our results are described in detail in Sec. <ref>.

In contrast, for the 2D toric code, rare events can be catastrophic: we find that as soon as p_ rare is above the bulk threshold, the entire code loses its threshold (i.e. the probability of a logical error no longer vanishes in the thermodynamic limit)[We often have in mind a situtation where p_ bulk can be tuned by improving the qubit quality, whereas p_ rare is an external parameter that we have no control over.
Therefore, we say the effects are catastrophic if there is a large enough  which makes the code undecodable even when p_ bulk→ 0.
However, it might be possible to gain some control over p_ rare in the hardware, as pointed out in Ref. <cit.>.]. Notably, this happens independently of the bulk error rate, i.e. even as p_ bulk→ 0. We note that similar error models were recently studied numerically by Tan et. al. <cit.>, where it was assumed that  rare events such as cosmic rays last for at most a finite duration. Within their setup, the rare events are found to be “benign”, in the sense that when the rare region error rate p_ rare is slightly above the bulk threshold, the code can remain in the decodable phase by correspondingly decreasing the bulk error rate p_ bulk. The key difference in our work is that we consider a more realistic setup where the rare events occur at a  finite rate in time. This, in turn, implies that the longest period with elevated error rates is typically unbounded from above: an L× L toric code system requires O(L) rounds of repeated measurements for decoding, which typically produces a largest rare-event sequence of size O(log(L)).  We argue that these largest rare regions dominate the logical failure rate, and lead to loss of threshold. Both our setup and that of Ref. <cit.>  can be understood within our theoretical framework, as we detail in Sec. <ref>.
Our results imply that techniques for suppressing lasting rare events <cit.> will be crucial to QEC with the toric code.

We also provide a physical interpretation of the difference between the 1D repetition code and the 2D toric code as follows.
Both models have excitations that are pointlike, and logical errors that are one-dimensional.
It is convenient to consider their dual models under Kramers-Wannier duality, after the random sign disorder in the stat mech models are neglected.
We obtain 2D and 3D Ising models with non-uniform couplings, respectively, where logical failure of the code can be related to dual Ising correlation functions in both cases.
The difference between the two phase diagrams can be attributed to whether the rare regions can order by themselves in isolation, in this dual picture. 

Finally, we note that our results for the 1D repetition code are obtained by leaning on known results for Griffiths physics in the celebrated McCoy-Wu model, which is a disordered 2D Ising model with correlated columnar disorder <cit.>. In contrast, the toric code problem yields a 3D RPGT with correlated disorder, which has not been studied much before in the context of Griffiths physics. Thus our analysis of this problem should also be of independent interest as a stat mech problem.






§ BACKGROUND AND MODELS



We will focus on two familiar examples of error correcting codes, the 1D repetition code and the 2D toric / surface code. Both are characterized by a set of stabilizer parity checks. Errors are characterized by their syndrome, the set of checks they violate. In both examples, these are point-like excitations (domain walls and anyons, respectively). The goal of error correction is to pair up these excitations in a way that undoes the effect of the error. One also has to deal with errors in the measurement of the stabilizers, which is usually done by combining information from many rounds of measurements; this introduces a time coordinate, making the two problems 1+1 and 2+1 dimensional, respectively. 

For any given error correcting algorithm, the logical failure rate refers to the probability that the algorithm fails to correctly recover the encoded logical information. It is related to the probability of large error chains whose syndromes are incorrectly paired up by the algorithm. The error threshold of the code is the maximal error rate (of both physical and measurement errors) such that the failure rate goes to zero as the number of (qu)bits is taken to infinity. Of particular interest is the maximum likelihood decoder, which corresponds to the theoretically optimal decoding and thus yields the largest threshold. 

While usually both the physical and the measurement error rates are taken to be constants, which are the same everywhere in the system, and for all the different measurement rounds, in realistic scenarios, they would be different for different qubits / stabilizers and can also fluctuate in time. We will be interested in the effect of such fluctuations, particularly those with long-range correlations in space or time (focusing on the former case). 












 §.§ Statistical mechanical models


The threshold and failure rate of stabilizer codes can be understood in terms of appropriate disordered statistical mechanics models <cit.>. 
Detailed derivations of these models are of secondary importance for our purposes, and we refer the reader to <cit.> for a more thorough explanation.
Throughout this paper, we focus on the illustrative examples of the 1D repetition code and the 2D toric code.
As both are CSS codes, we look solely at the Z stabilizers of the stabilizer code with stochastic single-qubit X noise acting on the system.

For the 1D repetition code, the stat mech model associated with decoding is known to be a 2D random bond Ising model (RBIM) on a square lattice,

    Z_ RBIM = ∑_{σ}exp(∑_⟨ j,k ⟩ K_jkτ_jkσ_j σ_k ),

and for the 2D toric code, the model is a 3D random plaquette (Ising) gauge theory (RPGT),

    Z_ RPGT = ∑_{σ}exp(∑_ K_τ_∏_⟨ j, k⟩∈σ_jk).

In Eq. (<ref>), the sum is over edges jk containing nearest-neighbor spins σ_j, σ_k.
In Eq. (<ref>), the sum is over plaquettes and each term is a product over spins which reside on the edges jk of the plaquette.
In both models, there is one Ising spin per constant time plane for each physical qubit.
Couplings that extend into the time direction (timelike) correspond to bitflip errors, and in Fig. <ref> they are represented as either vertical bonds (for the repetition code, see Fig. <ref>(a)) or plaquettes with normals pointing in spatial directions (for the toric code, see Fig. <ref>(b)).
Couplings within constant time planes (spacelike) take the form of Z stabilizers and correspond to measurement errors.
In Fig. <ref> they are represented by horizontal bonds (for the repetition code, see Fig. <ref>(a)) or plaquettes with normals pointing in the time z direction (for the toric code, see Fig. <ref>(b)).

For the ease of discussion, in the following we use α as a label for couplings, to refer to either a bond ⟨ j, k ⟩ in the RBIM, or a plaquette  in the RPGT. We denote the local error rate by p_α 
(i.e. it is a measurement error rate for a spacelike coupling, and a physical error rate for a timelike coupling). The local coupling strengths K_α are related to the local error rates via the Nishimori condition

    e^-2 K_α = p_α/1-p_α.

Here, this equation is understood to hold for each and every α, in both Eqs. (<ref>,<ref>). 
The τ_α variables incorporate the error history and introduce quenched random sign disorder into the system.
They take the value τ_α=-1 when an error happens on α i.e. with probability p_α, and take the value τ_α=+1 with probability 1-p_α.  




As we are interested in error models with spatio-temporal heterogeneity, we introduce a second type of randomness, and allow the local error rates p_α to differ for different α, while maintaining the Nishimori condition Eq. (<ref>) everywhere.
Correspondingly, the coupling constants K_α also vary in spacetime. Thus, we will study stat mech models with two types of disorder, namely (i) random sign disorder τ_α, and (ii) spatio-temporally varying K_α.
Physically, (i) can be understood as coming from different stochastic error realizations within a fixed error model, whereas (ii) is due to randomness in the error model itself, defined by the error rates at all locations {p_α}.





 §.§ Defects and failure rates




Condition (<ref>) ensures that the stat mech models correctly encode the success of maximal-likelihood decoders that have full knowledge of the error model.
In particular <cit.>, the failure and success probabilities of decoding are related to the free energy cost of the topological defect in the stat mech model which corresponds to a logical error.
For the RBIM, this is a domain wall; and for the RPGT, this is a magnetic flux tube (see Fig. <ref>), both extending perpendicular to the temporal direction.
As the defects in both models are linelike, they have a free energy cost which to leading order scales as 

    Δ F ∝σ(ℓ) ·ℓ,

where ℓ is the length of the defect and σ(ℓ) is its tension. In the ordered phases of the stat mech model (ferromagnetic for RBIM, deconfined for RPGT), defects are costly with diverging free-energy cost. This maps to the decodable phase with success probability approaching 1 as the system size is scaled up. A given defect may or may not correspond to a logical error. For instance, in the RBIM mapping of the 1+1D repetition code, only the defect in the horizontal direction is a logical operator (see Fig.  <ref>).[Throughout this work, we denote the line tension of defects associated with logical operators as σ_∥, as they are parallel to the rare regions. We denote the line tension associated with the transversal defect going in the temporal direction as σ_⊥, for it is perpendicular to the rare regions.]
In codes with uncorrelated disorder, these quantities are expected to scale the same in all phases. However, as we will see below, this need not be the case for long-range correlated disorder.

For a given error model {p_α}, which thereby fixes the set { K_α}, the relative success probability for this error model is given by <cit.>

    ({K_α}) = 
        [  e^- Δ F/1+e^- Δ F]_{τ_α},

where Δ F is the free energy cost of the defect within a particular realization of {τ_α} and [ …]_{τ_α} denotes the quenched average over {τ_α}.





Numerically, we use the minimium weight perfect matching (MWPM) decoder as implemented through the PyMatching package <cit.> to calculate .
Though many of our physical arguments are derived from considerations of the maximum likelihood decoder, in App. <ref> we provide arguments for why they are also expected to hold for the MWPM decoder.

We calculate  for each error model { p_α}.
This process yields a distribution of .
We do not expect  to be self-averaging, as it will be dominated by rare error models where the LRR are exceptionally large, see App. <ref> for further discussions.
Instead of the mean , we will mostly use the median  in our numerics.
From this we may numerically define the defect cost Δ F as

    Δ F ≡ -ln[ ()_ med/1- ()_ med],

in analogy with Eq. (<ref>).[Given the average over the sign disorder {τ_a} in Eq. (<ref>), Δ F in Eq. (<ref>) is not the median over error rate realizations of the sign-disorder-averaged free energy. This does not matter for our purposes.]





 §.§ Rare regions 
 

Our discussions so far have been fairly general.
Now, we turn to our primary focus, that is long range correlations in the distribution of non-uniform error rates, and the effect of these on the free-energy cost of defects. 
For concreteness, we focus on cases with uniform measurement error rates p_ meas.
They correspond to uniform spacelike couplings in the stat mech model.
In contrast, we set bit-flip error rates p_ bf to be spatially uniform on each time slice, but temporally varying between time slices, and are drawn from a Bernoulli distribution:

    p_ bf(t) =  
          p_bulk   with probability   1 - γ, 
    
          p_rare   with probability γ .
 
These produce varying timelike couplings in the stat mech model via Eq. (<ref>), see Fig. <ref>.
This is a minimal model of stochastic rare events  which globally affect the system, inspired by phenomena such as cosmic rays. 

Let us denote by p_0 the threshold of the code when the error rate is uniform, i.e. when = p_ bulk.
(This number is set by p_ meas, but we omit this dependence here.) In App. <ref>, we calculate p_0 for the models use in this work.
We are interested in the following regime,

    p_ bulk < p_0 < ,

which models an experiment that would be below threshold if not for rare events that temporarily elevates the error rate to .
Parts of the stat mech model where p_ bf = are henceforth referred to as rare regions.

Whenever γ > 0, 
we have rare regions whose typical temporal extent L_⊥ is finite, which should be benign from an error correction point of view.
However, 
we argue below that logical failure events will be dominated by the largest rare region (LRR), and with γ > 0 the largest L_⊥ will typically be unbounded from above in an infinitely long experiment.[Throughout the paper, we will take the time for the experiment (hence also the “height” of the stat mech models) to be proportional to L.]
In Sec. <ref> and Sec. <ref> below, we present this rare region analysis for the 1D repetition code and for the 2D toric code, respectively, using both analytic arguments and numerics.




§ LINEAR RARE REGIONS FOR THE 1+1D REPETITION CODE 




In this section, we analytically predict the effects of the correlated disorder in bitflip rates Eq. (<ref>) in the 1D repetition code using the model described in Sec. <ref>, and we numerically test these predictions.
The repetition code is a classical code that can serve as a simple theoretical model and as a experimental benchmark <cit.>.
In particular, here, the repetition code allows an analogy with the McCoy-Wu model, where the asymptotic scaling of the defect free energy can be predicted analytically <cit.>.


As we noted in the previous section, the stat mech model describing decoding the 1D repetition code of length L is an L × T RBIM for a time duration of T; we take T=L in the following.
Logical failure rates are controlled by defect free energy costs, and the defects in the RBIM are domain walls.
The correlated disorder in bitflip rates introduces rows of weak vertical bonds in the RBIM that also have a higher likelihood of sign errors; we call regions of consecutive weak bonds in time rare regions.





 §.§ Central assumption about largest rare region 


A central simplifying assumption we make in order to derive a qualitative phase diagram is that the logical failure of the code is dominated by the largest rare region (LRR) of the model[By “dominate," we do not require the LRR to strictly set the logical failure rates; rather, we only assume that studying the LRR will give correct predictions for the asymptotic behavior in the phases on either side of the threshold, even if the LRR by itself does not necessarily predict nonuniversal quantities like phase boundaries in the full L × L model.].
Correspondingly, in the stat mech model, we assume that the cost of the spatial defect Δ F is controlled by the LRR of weak couplings.
We will see that much of the numerical results can be understood qualitatively by focusing on the LRR and treating it as an isolated system. 

Our assumption can be understood from the intuitive picture that logical errors will most likely occur when the error rates are above threshold for the longest duration. One might also anticipate the assumption directly from the stat mech model by noting that the LRR has the most “room” for the domain wall to move (relative to smaller rare regions), and therefore the defect gains the most entropy (note that, since the bulk couplings are larger, it is energetically favorable for the defect to stay within rare regions). Both of these pictures neglect the interaction between the LRR and the bulk, as well as between LRR and other rare regions. We will find that much of the phase diagram is qualitatively described by LRR, although we note that there are discrepancies at large  sufficiently far above threshold.

We briefly summarize the leading order scaling of defect free energy in clean models where neither the random sign disorder nor the spatio-temporal randomness in {K_α} are present.
These clean models are a convenient toy picture for the LRR, and provide a basis for comparison with our numerical results of Δ F.


For the LRR in the 1D repetition code, we obtain the 2D Ising model living in a thin strip of dimensions L × L_⊥.
In the paramagnetic phase with K_ rare < K_c, a high-temperature expansion can be used to show that the line tension decays to zero with increasing L_⊥ as

    σ_∥(L) ∝ e^-α L_⊥,

where α depends continuously on K_ rare. We discuss this scaling in more detail in Appendices <ref> and <ref>.

For a rate of weak couplings γ > 0 and for a O(L) height of the stat mech model, we have that the typical L_⊥ for the LRR grows as

    L_⊥∝log L,

with a nonzero constant of proportionality. L_⊥ is distributed identically to the longest run of heads in a sequence of L biased coin flips, which is known to have mean ∝log(L) up to O(1) corrections <cit.>. This random variable also has just O(1) variance, meaning that its mean and median behave similarly. The mean ∝log(L) can be understood heuristically as striking a balance between L opportunities to have a length L_⊥ string of rare regions and the exponentially decaying probability of such a consecutive string: the expected number of rare regions of height L_⊥ is asymptotically proportional to L γ^L_⊥, which is one for L_⊥∼log L.

Together, Eqs. <ref> and (<ref>) reproduce the same scaling as in Eq. (<ref>). The value of z depends on several factors. It depends on the height of the rare region (and hence γ, which sets the proportionality constant in Eq. (<ref>)). It also depends on how deep into the paramagnetic phase the rare region is (and hence on K_rare and hence p_rare, which sets α in Eq. (<ref>)). 



 §.§ Comparison to McCoy-Wu

A useful comparison can be made between the stat mech model for the 1D repetition code and the McCoy-Wu model <cit.>.
The latter can be obtained from Eq. (<ref>) by removing the sign disorder (setting τ_jk = +1 everywhere) but still keeping the correlated disorder in the couplings {K_α}.
By analogy with Eq. (<ref>), we are interested in the regime, 

    K_ bulk > K_0^2D > K_ rare,

where K_0^2D is the critical coupling strength of the 2D uniform Ising model without any weak bonds.
The rare regions of the model are locally in the paramagnetic phase, whereas the bulk is still in the ferromagnetic phase, and the McCoy-Wu model is said to be in the Griffiths phase.
Refs. <cit.> showed that a defect parallel to the rare region direction has a line tension of the form (compare Eq. (<ref>))

    σ_∥(L) ∝ L^-z,

where L is the linear size of the system and z is a dynamical exponent which depends continuously on K_ rare <cit.>. z=0 when K_rare = K_0^2D. 

Unlike the LRR approximation, z tends to infinity at a finite value of K_rare = K_c^MW. This point coincides with the loss of spontaneous magnetization in the McCoy-Wu model, and is typically identified as the critical point of the McCoy-Wu model. By spontaneous magnetization, we mean the magnetization in the thermodynamic limit on taking a bulk longitudinal field h → 0^+. Note that other common metrics of order may show transitions that do not necessarily coincide with K_c^MW; vertical correlation functions that cut perpendicular to the rare regions can still decay to zero for some choices of K>K_c^MW. 

For K_rare < K_c^MW, the domain wall tension σ_∥ decays faster than any power law. This is not captured in the LRR approximation. 





 §.§ Predicted phase diagram

Using intuition from the McCoy-Wu model, we can now piece together the phase diagram shown in  Fig. <ref>(c) as a function of  at fixed p_ bulk < p_0^2D. Here, p_0^2D is the critical bitflip error rate of the L × L repetition code. This should be distinguished from the threshold p_D and the critical p_SM.

When we write  in this section, we consider the median .

First, when p_rare < p_0^2D, the entire system is in the decodable phase and the code is necessarily decodable. The RBIM exhibits ordinary ferromagnetic behavior such as a finite domain wall tension. This corresponds to a failure rate  decaying exponentially in system size.

For p_rare between p_0^2D and p_ SM, the rare regions are in the paramagnetic phase, but the stat mech model remains ordered. The stat mech model is in the ferromagnetic Griffiths phase with σ_∥ described by Eq. (<ref>), giving a domain wall cost of the form Δ F ≈ L σ_∥(L) ∝ L^1-z. 
We expect z to vary continuously with p_ rare. In particular, z=0 when = p_0^2D. This phase is denoted by the hatched region in Fig. <ref>(c).

The Griffiths phase further splits into two distinct regimes depending on whether z<1 or z>1. We denote the  at which z=1 as p_D (where p_D < p_SM is a distinct error rate from the critical p_SM), at which point  is a constant greater than zero that is asymptotically independent of system size.[A constant  at threshold assumes that the power law does not have subleading multiplicative corrections. For example, σ_∥∼log(L)/L^z would cause  at threshold to decay to 0 with L. We do not expect any such multiplicative corrections.]
When p_0^2D < p < p_D, we have 0<z<1, and the code is still decodable. 
Here, ℙ_fail tends to zero as L →∞, since the cost L^1-z diverges and  decays as a stretched exponential. 

On the other hand, p > p_D means z>1, and the domain wall cost L^1-z asymptotically vanishes, making → 1/2 and the code nondecodable. Nevertheless, this regime is still in ordered phase of the stat mech model, due to the fact that transversal domain walls (which are not related to ℙ_fail) continue to have a line tension σ_⊥∼𝒪(1), yielding a growing free energy cost; see Fig. <ref> below.



When p_rare is very large and above the critical disorder strength p_ SM of the stat mech model, the system fails to order and the domain wall tension decays faster than any power law, giving a vanishing domain wall cost and ≈ 1/2. This is captured by McCoy-Wu but is beyond the LRR approximation. We emphasize that the LRR nevertheless qualitatively captures the part of the phase diagram surrounding p_D, and only fails substantially above p_D. Note also that this discrepancy is small; the discrepancy is only in the qualitative description of how rapidly  approaches 1/2 at sufficiently large . 



 §.§ Numerical results






We now compare this predicted phase diagram to numerical simulations of the repetition code. For the numerics on the 1+1D repetition code, we take a length L system with L measurements.
We choose a uniform measurement error rate p_ meas = 0.11, which sets p_0^2D = 0.10.
We fix the bulk bitflip rate to be p_ bulk = 0.02 < p_0^2D, and vary the value of p_ rare.
We sample the rare regions according to the Bernoulli distribution Eq. (<ref>), where we take γ = 1/3 in both cases.

Our numerical results supporting the phase diagram are shown in Fig. <ref> and Fig. <ref>. In Fig. <ref>, we compute ℙ_fail as a function of p_rare. We identify the decodable-to-nondecodable transition from the crossing of curves with different L and find that it is clearly distinct from the threshold of the bulk system p_0^2D≈ 0.1[We also observe an empirical separation between the mean and typical (median) failure rate, not pictured in the main text, which we comment on in the Appendix <ref>.]. In Fig. <ref>, we depict the scaling of domain wall cost with system size. We see that there is an extended regime where the domain wall cost scales as a power-law with continuously varying exponent 1-z, consistent with our prediction for  within the Griffiths phase (p_0^2D, p_ SM). 

The inset of Fig. <ref> shows the numerically extracted z obtained from the linear fits in the main panel. We indeed find that it changes continuously with p_rare, dropping below 1 at p_D; however, we also observe that z remains finite even when p_rare is made smaller than the bulk threshold, contrary to our analytical argument. However, we note that our theoretical predictions neglected subleading contributions to the domain wall cost. In general, one might expect 

    Δ F ∼ L^1-z + O(L^α)

where L^α with α < 1-z is a subleading term dependent on . We do not attempt to compute such subleading corrections, but we note that these may cause significant finite-size effects, particularly in estimating z numerically.
We view the persistence of z>0, even for < p_0^2D where we expect z=0, as an artifact of such subleading terms.[In Appendix <ref>, we consider an L × L RBIM without correlated line disorder for comparison. There, we again fit z in the same manner as in the inset of Fig. <ref>, but we still find a z that fails to reach 0 even at small p_ bulk.]
To summarize, while our numerics on z in Fig. <ref> are suggestive, they are not yet conclusive.



To further investigate our theoretical predictions, we directly investigate an L ×log_2 L system with uniform bitflip error rates p_ bulk =  = p. A uniform system of this size mimicks the largest rare region of an L by L heterogeneous system.
With numerical results shown in Fig. <ref> for log_2 L ≤ 10, we confirm that a scale invariant point p_ D > p_0^2D exists in the latter systems, and that failure rates below this point decay as stretched exponentials in L.
We note however that the numerical value of p_ D is non-universal, so a direct comparison with the L × L heterogenous system cannot be made.

When  is between p_ D and p_ SM, the code is no longer decodable.
However, within this region, the overall phase of the RBIM is still ferromagnetic and the cost of the defect in the time direction continues to grow.
This defect can be thought of as a “temporal” logical operator, relevant for fault-tolerant gates involving moving logical information in space <cit.>.
Our numerical results in Fig. <ref> show that this logical 
has a higher threshold than when the logical was parallel to the correlation direction, confirming the phase diagram in Fig. <ref>.
As the order-disorder transition occurs when defects in all directions cease to grow in free energy cost with system size, it is at p_ SM where the system truly transitions to being paramagnetic.

Finally, we close this section by noting that p_D will additionally depend upon the value of γ used. Predictably, the code will be more tolerant to rare events if they are rarer and we show this explicitly in Appendix <ref>.





  §.§.§ Additional numerical details for the repetition code


Here we provide additional numerical details for Figs. <ref>, <ref>, <ref> and <ref>.
In Fig. <ref>, each point on each curve is determined from the median of 10^3 error model realizations { p_α}. For each error model realization,  is determined from averaging over 10^4 physical error configurations.
In Fig. <ref>, z is determined by a linear fit to logΔ F (determined from Eq. (<ref>)) vs log L. The extracted value of z should be viewed qualitatively. Due to potential subleading corrections which are present even in a system with uniform error rates (see Appendix  <ref>), we do not see z = 0 for < p_0^2D as predicted. Furthermore, at low error rates and large system sizes, the median  is often 0; which is why the smallest values of  are restricted to smaller system sizes. 
We note that techniques such as importance sampling (e.g. chapter 9 of Ref. <cit.>) may give greater access to this difficult regime of small , and leave detailed exploration to future work.

In Fig. <ref>, each point on each curve was determined from 10^6 physical error configurations. 

In Fig. <ref>, each point on each curve is determined from the median of 10^3 error model realizations.  for each error model realization is determined from 10^4 physical error configurations. To realize the transversal defect, we make the measurement error rate a random variable in the spatial direction, and set the bit flip error rate to 0.11. 




§ PLANAR RARE REGIONS FOR THE 2+1D TORIC CODE 






In this section we turn to time-dependent bitflip error rates in the 2D toric code.
They lead to rare regions in the random-plaquette gauge theory (RPGT) Eq. (<ref>) that are planar and have infinite extent in the two spatial dimensions, see Fig. <ref>(b).
Similarly to Sec. <ref>, we draw bitflip error rates on different time steps from the Bernoulli distribution in Eq. (<ref>), and choose the time duration to be T = L.

We proceed as in Sec. <ref>, first detailing our analytic predictions of the logical defect scaling for a quasi-2D system of dimensions L × L × L_⊥ with uniform error rates.
Treating this isolated system as a model of LRR, we then turn to the L × L × L heterogeneous system with non-uniform error rates, where quantitative agreement between numerical results and analytic predictions are found. Our key finding is the absence of a decodable Griffiths phase: as soon as p_rare exceeds the bulk threshold, and the rare regions are in the “wrong” phase, they make the decoder fail. As we will see, this is due to the 
2D nature
of the rare regions, which allows them to order by themselves.



 §.§ High temperature expansion within a clean planar LRR 


As we have seen in Sec. <ref>, clean models (where neither the random sign disorder or the spatial-temporal randomness in {K_α} is present) are useful for gaining intuition.

In parallel with Sec. <ref>, a high temperature expansion within a thin slab of dimensions L × L × L_⊥  yields

    σ_∥(L) ∝ e^-α L_⊥· L  ⇒   -logσ_∥(L) ∝ L · L_⊥,

see details in Appendix <ref>.
Once again, the LRR has width L_⊥∝log L, as in Eq. (<ref>).

As in Sec. <ref> and Eq. (<ref>), we perform this expansion for K < K_c^3D, where K_c^3D is the critical coupling strength of the 3D uniform bulk system[Strictly speaking, the expansion is only valid for K < K_c^2D(L_⊥), where K_c^2D(L_⊥) < K_c^3D is the L_⊥-dependent critical coupling strength of the slab with height L_⊥.
Our statement here can be justified by noting that K_c^2D(L_⊥) approaches K_c^3D as L_⊥ increases, see Sec. <ref>, <ref> for detailed discussions.].
This result suggest that σ_∥(L) vanishes exponentially with L in the limit L→∞.
The defect cost is therefore Δ F ∝ L^1-z L for some positive constant z, which vanishes in the thermodynamic limit when K < K_c^3D.
Translating this back to the decoding problem, an L × L toric code within time duration T = L_⊥∝log L becomes immediately undecodable whenever p > p_0^3D, where p_0^3D is the threshold of the 3D RPGT with uniform error rates.

Comparing Eq. (<ref>) with Eq. (<ref>), there is an extra multiplicative factor L in -logσ_∥(L) due to the increased dimension of the rare region.
The strong vanishing of σ_∥(L) in  Eq. (<ref>) suggests qualitative differences between linear and planar rare regions. 








 §.§ Dual picture 


Such qualitative differences can be more clearly appreciated in a dual picture. We summarize this here and explain this picture in depth in Appendix <ref>. Under a Kramers-Wannier duality, the linear and planar LRRs are described by quasi-1D and quasi-2D Ising models, respectively.
The defect line tension maps to the inverse dual correlation length;
the vanishing of the line tension (and therefore proliferation of line defects) corresponds to the ordering of the dual Ising models.
The difference in the scaling of σ can be therefore attributed to whether the dimension of the rare region is below or above the lower critical dimension of the Ising model. 

In particular, for the planar rare regions, the critical dual coupling strength K_c^2D, ∗(L_⊥) of the L × L × L_⊥ dual Ising model is finite, and can be arbitrarily close to K_c^3D, ∗ of a 3D bulk Ising model, as L_⊥∝log L can be arbitrarily large[In particular, under Kramers-Wannier duality, we have K_c^2D, ∗(L_⊥) > K_c^3D, ∗ for all L_⊥ < ∞.
As L_⊥→∞, K_c^2D, ∗(L_⊥) approaches K_c^3D, ∗ from above.
For any K^∗ > K_c^3D,∗, we will also have K^∗ > K_c^2D, ∗(L_⊥) for sufficiently large L_⊥.].
Therefore, for any K^∗ > K_c^3D,∗, the LRR will become ordered for sufficiently large L_⊥∝log L.

We can also understand these statements from the perspective of quantum models in one lower dimensions, related to the stat mech models via quantum-classical duality. We can relate the clean, L × L × L_⊥ classical problem to a quantum system on a L × L_⊥ lattice[Note that in this case, the “time” direction corresponds to one of the spatial dimensions of the toric code, rather than to physical time.]. In this quantum model, the calculation of σ(L) amounts to asking about the energy cost of inserting a point-like flux defect in the trivial paramagnetic phase (since we are dealing with rare regions). Using the quantum version of Kramers-Wannier duality, this is the same as the energy gap between the two symmetry broken ground states of a two-dimensional quantum Ising model in its ordered phase, which is exponentially small in the volume of the system, yielding the scaling in Eq. (<ref>). The same argument could be used to deduce the scaling (<ref>) in the case of 1D repetition code.



 §.§ Crossover scaling 


The picture above offered by the Kramers-Wannier duality for the clean model suggest that when analyzing numerical data of  at p_ rare≳ p_0^3D of the toric code problem, we must take into account the proximity of p_0^3D to the 2D critical point of the planar LRR.
We denote this L_⊥-dependent critical error rate as p_0^2D(L_⊥), and we have lim_L_⊥→∞ p_0^2D(L_⊥) = p_0^3D.
As we increase L_⊥, we expect a crossover from the 2D RBIM transition to the 3D RPGT bulk transition.



To capture the finite L_⊥ crossover, we propose the following phenomenological two-parameter scaling function

    (, L, L_⊥) ≈Φ[ (-p_0^3D) · L^1/ν_3, L_⊥ / L ].

We expect the scaling function to be descriptive of  when both L_⊥ and L are much larger than 1.
Here, ν_3 is the correlation length exponent for the 3D RPGT transition.
For any finite L_⊥ / L > 0, this reduces to the conventional scaling function of the 3D RPGT with a nonzero aspect ratio.
The planar limit is L_⊥ / L → 0, where we expect to see critical scaling near a 2D critical point.
With these considerations, we find that 

    p_0^2D(L_⊥) - p_0^3D = z_0 · L_⊥^-1/ν_3,

where z_0 is a non-universal positive constant. See Appendix <ref> for more detail.
We also find that at small L_⊥ / L the two-parameter scaling function Φ reduces to the following asymptotic form with a single parameter

    L_⊥ / L → 0=Φ^2D[ (-p_0^2D(L_⊥)) · L_⊥^1/ν_3·( L_⊥ / L)^-1/ν_2] 
       =
        Φ^2D[ ( (-p_0^3D) · L_⊥^1/ν_3 - z_0 ) ·( L_⊥ / L)^-1/ν_2].

Here, ν_2 is the correlation length exponent for the 2D RBIM transition, and Φ^2D is the corresponding universal scaling function.
Note that for any finite L_⊥, this function reduces to the standard 2D scaling form with the argument ∝ (-p_0^2D(L_⊥)) · L^1/ν_2.




 §.§ Predicted phase diagram



With these, and the assumptions that the logical failure of the code is dominated by the LRR (see Sec. <ref>),
we now discuss the predicted phase digram in Fig. <ref>(d).
Recall that we choose p_ bulk < p_0^3D, and vary  across p_0^3D.

When < p_0^3D, the entire system has error rate below p_0^3D, and the code is in its decodable phase.
Correspondingly, the RPGT is in its deconfined phase, where the conventional scaling - log∝ L holds.

When > p_0^3D, we predicted in Sec. <ref>, <ref> that the code becomes undecodable in the thermodynamic limit, i.e. lim_L→∞→ 3/4 for any p_ bulk.
In Fig. <ref>(d), we highlight the critical point p_ D(L) ≡ p_0^2D(L_⊥∝log L) of the planar LRR, which approaches p_0^3D from above as L →∞.
For  between p_0^3D and p_ D(L), the system is in the “mesoscopic crossover regime” (dashed green region), whose extent shrinks with increasing L.
Within this regime, we expect  to be suppressed with increasing L for L below a crossover length scale ξ_ crossover∝ ( - p_0^3D)^-ν_3.
Above the crossover length scale,  grows with L, and eventually saturates to 3/4 in the infinite L limit.
Therefore, a non-monotonic dependence of  on L is expected.
All these phenomenology should be captured by the scaling function Eq. (<ref>) for sufficiently large L, L_⊥, and for  sufficiently close to p_0^3D.

Similar to the 1+1D repetition code, there is again a range of values of  between p_D(L) and p_ SM, where the rare regions are disordered but the line tension σ_⊥ of the transveral defect (perpendicular to the planar rare regions) is still finite. Due to the presence of the disordered rare regions, we call this regime a Griffiths phase but unlike the 1+1D repetition code, it is non-decodable.





 §.§ Numerical results









We calculate  for both an L× L × L heterogeneous system with rare regions, and an L× L ×log_2 L isolated system with uniform error rates.
In both cases, we choose a uniform measurement error rate p_ meas = 0.01, which sets p_0^3D = 0.045.
For the L × L × L system, we choose p_ bulk = 0.01 < p_0^3D, and generate the rare regions according to Eq. (<ref>) with γ = 1/3.

Our numerical results are shown in Fig. <ref> and <ref>.
In Fig. <ref> for the L× L × L system, we find no scale-invariant point p_ D for the median .
Instead, we find a crossover regime with a non-monotonic dependence of  on L.
This is consistent with our expectations from the discussions above.
In Fig. <ref> for the isolated L× L × L_⊥ system, these predictions are more clearly borne out.
When setting L_⊥ = log_2 L, we observe a clear downward drift of p_D(L) towards p_0^3D with increasing L, as well as the non-monotonic dependence of  on L.
We also extract p_0^2D(L_⊥) at fixed finite values of L_⊥.
The results can be fitted well to Eq. (<ref>), see inset of Fig. <ref>.


Furthermore, the single-parameter scaling function Eq. (<ref>) is relevant as L_⊥≪ L, and can be directly tested.
In Fig. <ref>(triangles) we show the data collapse of the heterogeneous L× L× L system, and find good agreement with Eq. (<ref>) when setting L_⊥∝log L in its parameter.
In Fig. <ref>(squares), collapse of data from the L× L × L_⊥ system is shown, where agreements are again found.









Finally, the functional form of Φ^2D in Eq. (<ref>) is universal and thus can be extracted independently by decoding the 2D toric code with perfect measurements.
These results (Fig. <ref>(circles)), when plotted together with the numerical scaling collapses from the previous two numerical experiments, lie on top of each other (up to a constant rescaling of their parameters), see Fig. <ref>.
This strongly supports the crossover scaling picture provided by Eqs. (<ref>,<ref>).
In particular, the data collapse is strong evidence for our theoretical prediction that the 2D scaling function, modulo a rescaling of its argument, describes the failure probability close to threshold even in the (2+1)D planar-disordered toric code described by a 3D stat mech model. 







To summarize, in the toric code, whenever p_ rare > p_0^3D, planar rare regions render decoding entirely impossible in the thermodynamic limit[As a consistency check, this should follow directly from Eq. (<ref>). Taking L_⊥∝log L, we indeed find lim_L→∞ = 3/4 for any p > p_0^3D.]. However, a transverse defect perpendicular to the planar rare regions may still have a diverging Δ F. Indeed, the defect cost will diverge whenever the transvserse correlation length in the dual Ising model is finite. It is known <cit.> that one can have a spontaneous magnetization without a diverging transverse correlation length in 3d Ising models with planar disorder, corresponding in the primal gauge theory to a non-decodable region that has a diverging transvserse defect cost. These dualities are complicated by the presence of quenched random sign disorder, which we have largely neglected in our discussion of dualities. We do find such a region via numerical results on small systems L ≤ 18 in Appendix <ref>, and we find a transition in the transverse defect cost at a p_ SM > p_0^3D. We leave a more thorough investigation of this transition to future work.








  §.§.§ Additional Numerical Details


Here we provide additional numerical details for Fig. <ref> and Fig. <ref>. In Fig <ref>, each point on each curve is determined from the median of 10^3 error model realizations and  for each error model realization is determined from 10^4 physical error configurations. For Fig. <ref>, for L = 4,8,16,32, each point is determined from 10^5 physical error configurations. For L = 64,128,256 we use 10^4 physical error configurations. In the inset, for each value of L_⊥, we obtain the critical error rate from system sizes of L = 24 to L = 64. For each system size, we obtain  from 10^5 physical error configurations. 

In Fig. <ref>, we plot the same data from Fig. <ref> along with additional data collected for a toric code with no measurement errors.  In all data collapses, we take ν_3 = 1 and ν_2 = 1.5.
We use z_0 = 0.08 for the planar-disordered 2+1D toric code and z_0 = 0.09 for the log_2 L 2+1D toric code. 
z_0 is not universal, and it is reasonable to expect different z_0 values for these two systems.
For the 2D toric code with no meausrement errors, we use w = (p - p_0^2D)· L^1/ν_2 for the argument of the scaling form, where p_0^2D≈ 0.103 <cit.>.
A slight x axis rescaling was also necessary for good collapse. For the log_2 L toric code, we rescale the x axis by 1.1 and for the disordered toric code, we rescale by 1.2.





§ CONCLUSIONS



In this work we analyzed the performance of topological quantum codes in the presence of non-uniform error rates that are long range correlated.
We find that rare events with increased error rates (above the bulk threshold) can have dramatic effects on the code performance.
We point out crucial differences between linear and planar rare regions: the former lead to a new decodable phase for the 1D repetition code where the logical failure rate scales as a stretched exponential in the code distance, while the latter make decoding immediately impossible for the 2D toric code, as soon as the error rate in the rare region exceeds the global decoding threshold[Our analysis can be immediately extended to the case with fabrication errors, leading to parts of the code patch which have a larger than bulk error rate. They will result in linear rare regions parallel to the temporal direction, for both the 1D repetition code and the 2D toric code. For this reason, our analysis implies that that a decodable stretched exponential phase is present in both cases.].
We expect these analyses to be more broadly applicable to other topological codes with point-like excitations.

While we focused on toy error models for simplicity, the asymptotic scaling thus predicted are expected to be universal, and can hold for more realistic error models. 
For instance, fabrication errors may lead to qubits with different fixed error rates and this may be better modeled by a noise distribution other than Bernoulli <cit.>.
For cosmic ray events, an error model incorporating the relaxation time of the qubits can be considered.


While current quantum error correction (QEC) experiments are restricted to small code distances, we expect our results to be descriptive of future experiments when they approach the scaling limit.
More broadly, QEC experiments provide new platforms and new motivations for exploring disorder physics.
This work provides first explorations in this direction, and in future work it would be interesting to consider other models, including those in (generalized) gauge theories which are themselves under-motivated in solid state systems.

 



§ ACKNOWLEDGEMENTS


We acknowledge helpful discussions with Aditya Mahadevan, Arpit Dua, David Huse, Chaitanya Murthy and especially Akshat Pandey.
We thank Hengyun Zhou for helpful discussions and for bringing Ref. <cit.> to our attention.

A.S. acknowledges support from the National Science Foundation Graduate Research Fellowship. N.O.D. acknowledges support
from the ARCS Foundation for ARCS Scholar funding. Y.L. was supported in part by the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grant GBMF8686. Y.L. and T.R. were supported in part by the Stanford Q-FARM Bloch Postdoctoral Fellowship in Quantum Science and Engineering. 
V.K. acknowledges support from the Office of Naval Research Young Investigator Program (ONR YIP) under Award Number N00014-24-1-2098, the Alfred P. Sloan Foundation through a Sloan Research Fellowship and the Packard Foundation through a Packard Fellowship in Science and Engineering.
Numerical simulations were performed on Stanford Research Computing Center's Sherlock cluster. 











§ MINIMUM WEIGHT PERFECT MATCHING


Throughout this work, we use minimum weight perfect matching (MWPM) for decoding.
MWPM finds the minimum weight error consistent with the observed syndrome.
As the weight of an error is associated with an energy, MWPM neglects entropic contributions and can be thought of as a zero temperature decoder. 

MWPM works as follows.
First, a decoding graph is defined such that every node of the graph corresponds to a check of the code, and every edge in the graph represents a local error that connects all checks triggered by this error.
We assume that each edge connects at most two checks (i.e. no hyperedges), which is the case for all the model considered in this paper.
Second, a syndrome measurement marks which of the checks were triggered and those nodes are highlighted on the decoding graph.
Next, the decoder works to pair up each of the highlighted nodes on the check graph with the minimum number of edges.
From such a minimum weight matching, a proposed correction operation can be obtained.
The pairing can be carried out efficiently in O(n^3) time, where n is the number of nodes on the decoding graph <cit.>.

To account for non-uniform error rates, weights can be assigned to edges based upon the rate of the corresponding error.
If error a has error rate p_a, then its respective edge on the decoding graph has weight ln1-p_a/p_a. 
Throughout this work, we supply the decoder with the information about the non-uniform error rates.



 §.§ Repetition Code and Toric Code with Uniform Error Rates



Here we provide baseline numerics using MWPM for the repetition code and the toric code. We utilize a uniform error rate in both systems in order to obtain p_0^2D and p_0^3D. This data is shown in Fig. <ref>(a) and Fig. <ref> respectively. The values we obtain are consistent with those found in the literature <cit.>.

For the repetition code, we also calculate the defect cost similar to the procedure used in Fig. <ref>. This is shown in Fig. <ref>(b). The defect cost is expected to scale purely exponentially (i.e. z = 0). However, due to subleading corrections in the defect cost scaling, the empirical value of z is slightly greater than 0. Furthermore, the error rates shown here are close enough to the threshold that we cannot rule out contributions to the scaling form from critical scaling. This is to be contrasted with the case of the rare-event repetition code in the main text where the nonzero empirical value of z for > p_0^2D is theoretically understood to be coming from a sublinear leading term, and we see its effects reasonably far from the threshold.







 §.§ Griffiths physics for MWPM


We used free energies of defects to inform our theoretical treatment of success probabilities. This picture holds for the maximum likelihood (ML) decoder, where success probabilities can be reframed in terms of free energies of defects in statistical mechanical models with couplings K=β J set in terms of the error rates by “Nishimori conditions."
For MWPM,  is actually set by the probability that introducing a defect lowers the energy of the system. Quantities that depend on details of energies probe zero temperature physics, which can be in a different phase from the finite-temperature phase probed by the ML decoder.
However, we argue that the relevant Griffiths physics at play does not change much from the picture for the ML decoder.

For specificity, consider the case of the RBIM. Here, the defect is a domain wall formed by flipping the sign of a column of couplings. For sufficiently large disorder strengths at zero temperature, there is a phase transition from the ferromagnet to a spin glass phase. This spin glass phase will control the physics of the “wrong phase" spatial regions in the Griffiths phase.
Assume that  is largely set by the largest rare region living in the wrong phase. In an L by L system, this largest rare region has size L_⊥ by L with L_⊥ scaling as ln(L). The corresponding  for MWPM on a model with dimensions of the rare region will be bounded below by that of the maximal likelihood decoder on the same region. 

However, the maximal likelihood decoder has its failure probability controlled by free energy differences, and the free energy difference in the paramagnetic phase is exponentially suppressed in system size as a L e^-b L_⊥ for some a and b that are functions of error rates but are independent of system size. This gives ∼ 1/(1+e^a L e^-b L_⊥) for the ML decoder whenever the error rate is sufficient to enter the paramagnetic phase. Note again that  for the MWPM decoder is also bounded below by this quantity; though this is not a proof, we believe it plausible that  for the MWPM decoder in the spin glass behaves in a similar manner to that of the ML decoder in the paramagnet. That is, we believe ∼ 1/(1+e^a' L e^-b' L_⊥) for some new functions a',b' of the error rates. Furthermore, we expect this functional form to hold through the whole spin glass phase, not only where the error rate is sufficient for the better ML decoder to leave the ferromagnetic phase. Under this assumption, the phases of  in Fig. <ref>c) remain unchanged for the MWPM decoder, even if the phase boundaries shift relative to the ML decoder. 

An analogous argument holds for MWPM in the 3D gauge theory.



 §.§ Mean-Median Separation


In the main text, we noted that our numerical results yielded a separation between the mean and median logical failure rates as calculated from the MWPM decoder. This phenomenon is shown in Fig. <ref>. When there is uncorrelated or short range correlated disorder, this  distribution is expected to be approximately Gaussian. We find that this is emphatically not the case for long range correlations. There are several sources of randomness that could lead to asymptotically different scaling forms for the two quantities. For instance, the largest rare region in a given realization could be larger than its typical size by a big 𝒪(1) constant times log(L), which could enhance the failure rate to be 𝒪(1) instead of decaying. Those rare disorder realizations, which occur with probability ∼ 1/poly(L), would skew the mean failure probability to be at least ∼ 1/poly(L) rather than stretched exponential. Such uncommon realizations will cause the mean to scale differently than the median. Though we have theoretical reason to suspect that the mean and median should scale differently (power-law versus stretched exponential), we do not resolve this asymptotically different scaling in the numerics.





 §.§ Distribution Dependence
 

In Fig. <ref>,  we present data on how p_D for the repetition code changes with γ. We predictably find that p_D increases with decreasing γ.






§ DERIVATION OF DEFECT TENSIONS


Here, we will derive the defect tension for both the Ising model and the plaquette gauge theory using low and high temperature expansions. These derivations will be asymptotic in nature. Though the models discussed in the main text are disordered and contain couplings with both sign and magnitude disorder, we expect that the bulk phase properties to have the same asymptotic scaling behavior as the clean models.  



 §.§ 2D Ising Model


In a ferromagnet, the defect cost scales as L with a 𝒪(1) line tension. We can directly calculate this from the low temperature expansion. At lowest order in e^-K, the partition function of an L by L_⊥ clean Ising system is Z ∼ 2 e^2K L L_⊥. On introducing a domain wall via antiperiodic boundary conditions along the length L direction, the lowest order contribution is Z_DW∼ 2L_⊥ e^2 K L L_⊥ -2KL. This contribution comes from the L_⊥ locations that the shortest domain wall of length L can be placed. The definition of the domain wall free energy (up to a multiplicative factor of β) is Δ F = log(Z) - log(Z_DW), giving Δ F ∼ 2 K L - log(L_⊥) ∼ 2KL and a domain wall tension of 2K. Higher-order contributions renormalize the domain wall tension of 2K via a series in e^-K.

In a paramagnet, the defect cost should scale as L e^-zL_⊥. We may calculate this from the high temperature expansion, which sums over weighted loops on the real lattice. Recall the high temperature expansion of the Ising model:

    Z = (cosh K)^2N∑_{σ_i }∏_i,ν( 1 + (tanh K) σ_i σ_i+ν)

where ν is a lattice translation. Since σ_i^2 = 1 and ∑_σ_i = ± 1σ_i = 0, only terms which contain no factors of σ_i can contribute to the partition function. These terms correspond to closed loops built out of bonds. The first few terms of the partition function series are
 
    Z = (cosh K)^2N(1 + L L_⊥ (tanh K)^4 + 2 LL_⊥ (tanh K)^6 + …)

The effect of anti-periodic boundary conditions (i.e. flipping the sign of the couplings on some column) is to put a - sign on some of the loops. Note that any contractible loop must contain an even number of links on the antiferromagnetic column.
Thus, the first nontrivial contribution to

    Δ F = log(Z) - log(Z_DW)

will be from the shortest-length non-contractible loops that run perpendicular to the antiferromagnetic column. These loops have length L_⊥, and there are L of them. 
Since such loops carry opposite signs between Z and Z_DW, we have

    Δ F ∼ 2 L (tanh K)^L_⊥

with manifest exponential decay in L_⊥.




 §.§ 3D Gauge Theory


Magnetic flux tubes are the defects of the ℤ_2 lattice gauge theory stat mech model corresponding to the 2+1D toric code. In the low temperature phase of the clean gauge theory, we may obtain the cost of this defect through a low temperature expansion. This cost comes from the free energy difference of flipping the sign of the coupling of a column of plaquettes. The defect cost is

    Δ F    = F_ flipped column - F_ no flipped column
       = -logZ_ flipped column/Z_ no flipped column.


Similar to the 2D Ising model, at lowest order in e^-K, the partition function of a clean L_x × L_y × L_z system is Z_ no flipped column∼ 4e^3K L_x L_y L_z. Upon fixing a flux tube by flipping a the sign of row of plaquettes transverse to the xy plane, the lowest order contribution to Z_ flipped column is ∼ 4L_x L_y e^K(3L_x L_y L_z - 4L_z. This form is the result of the column of antiferromagnetic bonds lowering the energy of the ground state configuration. There are L_x L_y places to put this column. The flux tube cost becomes Δ F = log (Z_ no flipped column) - log (Z_ flipped column) again up to a multiplicative constant of β. This reduces to Δ F ∼ 4KL_z - log L_x L_y ∼ 4KL_z. The flux tube tension is an 𝒪(1) number, 4K and again, higher order contributions renormalize the tension in powers of e^-K.

Next, we calculate how the defect cost should scale in the high temperature phase. The first non-unity term in the high temperature expansion are closed surfaces of area 6. However, any closed surface which is contractible will always contain even numbers of plaquettes from the tube, similar to the 2D case. Therefore, the first term which is different between the Z_ flipped column and Z_ no flipped column will be the product of all the plaquettes within a plane pierced by the tube. This occurs at order L_xL_y. We obtain the defect cost in the high temperature phase to be

    Δ F ∼ 2 L_z (tanhK)^L_xL_y

with manifest exponential decay.
We see that the defect cost decreases exponentially in the area of the plane within the high temperature phase.




§ DISCUSSION OF THE PHASE DIAGRAMS WITH KRAMERS-WANNIER DUALITY


In this Appendix, we provide additional discussion of the phase diagrams in Fig. <ref>. We emphasize how differences in these phase diagrams arise from the dimensions of the rare regions. 

To build intuition, we consider success probabilities controlled by free energies in the context of maximum likelihood decoding. We neglect the random sign disorder in these stat mech models for simplicity, as we do not believe the sign disorder affects the relevant properties of the phases. At zero temperature (appropriate for describing the MWPM numerics in the main text), the phase transitions are in fact driven entirely by the random sign disorder; however, see the discussion in Appendix <ref>.

We make use of Kramers-Wannier dualities where defect free energy costs are mapped to spin correlation functions in the dual Ising models. Within this formulation, we attribute the difference in the phase diagrams to whether the dimension of the rare regions are below or above the lower critical dimension of the dual Ising models. In particular, the fact that a two dimensional Ising model can order is responsible for the absence of a stretched exponential phase of the toric code in (2+1)D. 






 §.§ Quasi-1D rare regions for repetition code in (1+1)D spacetime


As we explain in the main text, the stat mech model relevant to the repetition code in (1+1)D with non-uniform noise rates is the McCoy-Wu model (after neglecting weak random sign disorder). The z direction corresponds to time, where faulty measurements are performed up to a time t=L_z, and perfect measurements (simulating readout via single-site measurements) are performed at t=L_z.

In the stat mech model, 

    -  = Z_++-Z_+-/Z_+++Z_+-

where Z_++ refers to fixing both the top and bottom boundary to be all σ = 1 and Z_+- refers to fixing the top and bottom boundaries to have opposite orientation. Note that this quantity can be interpreted in two equivalent ways in this primal model. It is asymptotically 1 whenever the free-energy cost of inserting a domain wall grows unboundedly in system size. It additionally can be rewritten as -  = ⟨σ_i,z=0σ_j,z=L_z⟩_∞ in a model without fixed boundary conditions but with infinite-strength horizontal couplings K_x = ∞ on the top and bottom boundaries. 

When thinking in terms of free-energy costs, it's useful to consider just

    = Z_+-/Z_+++Z_+-

and in particular the ratio of Z_+-/Z_++≤ 1. Note that when Z_+-/Z_++ is small,  is ∼Z_+-/Z_++, and when Z_+-/Z_++∼ 1, ∼1/2.


We can rewrite Z_+-/Z_++ in terms of quantities in the Kramers-Wannier dual models, and we detail this construction in Fig. <ref> and Fig. <ref> for the respective cases of a periodic and open repetition code. For simplicity, in the figures we do not visually distinguish between rare and bulk regions; later, we will view these as uniform L by L_⊥ systems that model the largest rare regions of height L_⊥∼ln(L).





Under Kramers-Wannier duality, the spins of the dual model (which we will label as μ to distinguish them from the spins σ of the primal model) live on the dual lattice, and local couplings K = β J are in one-to-one correspondence with the dual couplings K → K^* = -1/2log(tanh(K)).
The rare regions have K_x, rare < K_c^2D < K_x, bulk, see Fig. <ref>(a).
Therefore, in the dual Ising model, the local couplings in the rare regions and in the bulk satisfy K_x,  rare^∗ > K_c^2D > K_x,  bulk^∗. In the main text, K_z is always set to K_c^2D in all regions; more generally, K_x,c^2D would then be a function of K_z.

For the repetition code in periodic boundary conditions, an additional η “domain wall" degree of freedom is required in the dual model. This η variable allows the low temperature expansion of the dual model to match the high temperature expansion of the primal model; the latter includes lines that cross the bulk and terminate at the bottom and top boundaries. Such lines correspond to a single vertical domain wall in the dual model. When the horizontal boundary conditions are periodic, domain walls naively come in pairs in the dual model. The introduction of η allows for single domain walls in the dual model by converting a column of couplings in the dual model to be antiferromagnetic instead of ferromagnetic. In particular, as noted in Fig. <ref>, the couplings on such a column are weighted by η. Summing over η=± 1 thus sums over all configurations with an even and odd number of vertical domain walls in the dual model. 

In the high temperature expansion of Z_++, configurations with even and odd numbers of vertical lines crossing the bulk come with the same sign. However, such configurations come with opposite signs for Z_+-. By weighting the Boltzmann weights of the dual model by η, we reproduce the relative sign differences. 

That is,

    Z_+-/Z_++ = ⟨η⟩_*

where * denotes an expectation in the dual model.

In our numerics, we always take periodic boundary conditions for the repetition code. However, the dual of the free energy cost of a domain wall can be massaged into a conceptually simpler form in open boundary conditions. We summarize the dual model for open boundary conditions in Fig. <ref>. 

In Fig. <ref>(b), we show the corresponding dual model. However, unlike periodic boundary conditions, this expression can be massaged into a more natural correlation function between dual spins on opposite edges.

By sending μ→ -μ to the right or left of the η-weighted bonds, we can convert configurations with η=-1 into configurations with η=1 but oppositely polarized left and right boundary conditions. We can view these fully polarized boundary conditions as coming from free boundary conditions with infinite couplings. Furthermore, the relevant signs of configurations can be found by weighting with μ_L μ_R instead of η; here μ_L is any dual spin on the left boundary and μ_R is any dual spin on the right boundary; the infinite couplings on the boundary make the choice immaterial.

Explicitly,

    Z_+-/Z_++ = ⟨μ_Lμ_ R⟩_*,∞     (OBC)

for any choice of μ_L on the leftmost boundary and any choice of μ_R on the rightmost boundary. Here ∞ denotes that the left and right boundaries have infinite vertical couplings; see Fig. <ref>(c).

From this relation, we can identify the line tension of the domain wall with the inverse correlation length of the dual Ising model, namely σ = (ξ^∗)^-1.

Below, we focus on the largest rare regions (LRR). We assume that - is controlled by configurations where the domain wall is entirely contained within largest paramagnetic rare region, where the domain wall gains the most entropy and where the line tension is the smallest. For conceptual simplicity, we use an L by L_⊥ system as a model of the LRR. The horizontal couplings are all (K_x, rare)^*. 

The height of the largest rare region L_⊥ goes as ln(L), but it is worthwhile to consider L_⊥ fixed and finite first. Importantly, the dual Ising model is (quasi) one-dimensional and does not order at any nonzero temperature (K_x, rare)^* < ∞, and always has a finite correlation length ξ^∗.
Correspondingly, the line tension in the primal model σ = (ξ^∗)^-1 is nonzero at any non-infinite temeperature K_x, rare > 0.

The infinite critical temperature for the vanishing of this domain wall line tension can be compared with the corresponding decoding problem at L_⊥ = O(1).
When introducing the random sign disorder back into the primal Ising model according to the Nishimori conditions, the primal model describes the decoding problem of a repetition code with length L and run for L_⊥ time steps.
In this case, perfect state initialization and perfect final syndrome measurement are assumed, so that the domain wall can only fluctuate between y ∈ [0, L_⊥].
It is easy to show that this code has threshold p_ th = 0.5 for any L_⊥ = O(1).

For both finite L_⊥ and L_⊥ = ln(L), the μ-μ correlation function can be obtained from a low temperature expansion within the dual model, yielding (ξ^∗)^-1∝ e^-2K^∗· L_⊥.
This expansion is similar to that leading to Eq. (<ref>).
Alternatively, viewing the dual partition function as a path integral with the x direction as a Euclidean time, the inverse correlation length is given by the splitting between ground state energies of the transfer matrix, i.e. a transverse field Ising model of length L_⊥.
This is again exponentially suppressed by L_⊥ when K_x,  rare^∗ > K_c^2D. 

This points to a difference between strictly finite L_⊥ and L_⊥ growing unboundedly with L; the latter can have an asymptotically vanishing inverse correlation length. For L_⊥∼log(L) and K_x,  rare^∗ > K_c^2D, the inverse correlation length in the dual model vanishes as a power law in L. Correspondingly, the domain wall tension in the primal model decays algebraically 1/L^z, and the domain wall cost goes as L^1-z, where z is a function of the coupling strength. 



 §.§ Quasi-2D rare regions for toric code in (2+1)D spacetime


After neglecting weak random sign disorder, the toric code in (2+1)D spacetime is described by a ℤ_2 lattice gauge theory in three dimensions, where the coupling strengths are uniform within each plane but may vary from plane to plane, see Fig. <ref>.
We again focus on the LRR, and treat it as a quasi-2D system. We similarly assume that a homologically nontrivial flux loop along the x or y directions receives the most contribution from configurations where the flux loop is completely contained within the LRR. More precisely, we assume that the behavior of the phases can be understood in terms of the properties of the LRR, even if some of the phase boundaries change when considering the additional effects of smaller rare regions.

Under Kramers-Wannier duality, the LRR is described by a quasi-2D Ising model with dimensions L_x × L_y × L_⊥, where L_x, L_y →∞. The height of the largest rare region L_⊥∼ln(L_z) will therefore slowly diverge with system size L, but it is again worthwhile to consider L_⊥ fixed and finite first.

Within the LRR we have T_ rare^∗ < T_c^3D, *, where T_c^3D, * is the critical temperature of the dual 3D Ising model.
Under appropriate boundary conditions, we also have the relation Eq. (<ref>) and can identify σ = (ξ^∗)^-1. The dual Ising model in the LRR can develop long range order at finite temperature, with a critical temperature denoted T_c^2D, ∗(L_⊥).  T_c^2D, ∗(L_⊥) increases monotonically with L_⊥, and approaches the 3d transition temperature T_c^3D, * as L_⊥→∞.
Thus, for any T_ rare^∗ > T_c^3D, *, we have T_ rare^∗ > T_c^2D, ∗(L_⊥) for L_⊥ larger than a T_ rare^∗-dependent constant. 
Correspondingly, whenever T_ rare^∗ > T_c^3D, *, the inverse correlation length and the domain wall line tension asymptotically vanish as L_⊥ increases.

In the toric code decoding problem with L_x = L_y = L_z = L, L_⊥∝ln(L) and hence L_⊥ is unbounded above as L grows. The above discussion about the approach of T_c^2D, ∗(L_⊥) to T_c^3D, * implies that whenever > p_c^3D, logical errors will then proliferate and the code is not decodable.
The contrast with the repetition code case can  therefore be attributed to the dimension of the rare regions.


When T_ rare^∗ < T_c^2D, ∗(L_⊥), the correlation length can either be obtained via an expansion along the lines of Appendix <ref>, or via the ground state splitting of a transverse field Ising model of size L_y × L_z, both yielding lnσ∝ - L_y · L_z for a defect in the x direction.










 §.§ Transversal defects for toric code in (2+1)D spacetime 



In our time dependent error rate toric code, we obtained a gauge theory where the plaquette coupling magnitudes were completely correlated in the spatial direction. By duality, this is related to a 3D Ising model with bond coupling magnitudes which are correlated in planes. The 3D Ising model with planar defects is believed to have a “smeared" phase transition <cit.>.
As each of the rare regions have infinite extent in two spatial dimensions, they are able to undergo phase transitions independently of the bulk. Then, as one lowers the temperature, different parts of the system order independently at different temperatures and the global order parameter develops smoothly from 0, when the smallest rare region orders. At the temperature where the smallest rare region orders, there is an essential singularity in the free energy. 

In Figure <ref>, we show numerics for the threshold of the transverse defect for the toric code in (2+1)D spacetime.
This defect goes in the temporal direction, perpendicular to the planar rare regions. In the repetition code, we were able to interpret the vanishing of the transverse defect as the point at which the underlying statistical mechanics model undergoes a phase transition. It is not clear whether that interpretation is true for the gauge theory. In future work, it would be interesting to further investigate this and whether symptoms of the smearing could be observed in the code. Furthermore, it would also be interesting to study what happens in the classical statistical mechanics model with planar defects when there is also uncorrelated sign disorder. 







§ CROSSOVER SCALING FUNCTION


Here we discuss in some detail the crossover scaling function for an L× L × L_⊥ system with uniform error rate p.
We posit that the mean failure probability of the code near p_0^3D is captured by the following phenomenological crossover scaling function (see Eq. (<ref>))

    (p, L, L_⊥) ≈Φ[ (p-p_0^3D) · L^1/ν_3, L_⊥ / L ] ≡Φ[x,y].

As usual for finite size scaling, we require both L and L_⊥ to be large for the scaling function to be descriptive.
We are particularly interested in the behavior of this function when y → 0, as the largest rare regions will have typically (L_⊥)_ max∝ln L.

For future convenience, we define single-parameter scaling functions Φ^3D_y[z] for each y,

    Φ^3D_y[z ≡ y^1/ν_3· x] ≡Φ[x,y].

They are therefore cross sections of Φ at constant values of y (up to a rescaling of x).
For any finite y > 0, Φ^3D_y is analytic, and describes the 3D bulk RPGT transition with aspect ratio y.


We extract the functional form of Φ^3D_y → 0[z] in two steps.

(i)
We first take L →∞, while keeping (p-p_0^3D) and L_⊥ finite.[In this limit, we have y → 0, x →∞, but z = y^1/ν_3· x = (p-p_0^3D) · L_⊥^1/ν_3 remains finite. This justifies our choice of the parameter z for Φ_y^3D, see Eq. (<ref>).]
In this case, the system is an infinite 2D slab with height L_⊥, which has an L_⊥-dependent critical error rate, denoted p_0^2D(L_⊥).
Therefore, as we take L →∞, we expect that

    = 
    3/4·Θ(p-p_0^2D(L_⊥)).

On the other hand, by definition of Φ^3D_y we have

    = Φ^3D_y=0[z = (p-p_0^3D) · L_⊥^1/ν_3].

Comparing Eqs. (<ref>, <ref>), we conclude that Φ^3D_y=0 is necessarily singular, with a step singularity at z = z_0.
Matching the location of the singularity, we have

    z_0 = (p_0^2D(L_⊥) - p_0^3D) · L_⊥^1/ν_3
       ⇒   p_0^2D(L_⊥) - p_0^3D = z_0 · L_⊥^-1/ν_3,

see also Eq. (<ref>).
Again, this relation holds when L_⊥ is sufficiently large.
Therefore, from the crossover scaling function we can infer how p_0^2D(L_⊥) approaches p_0^3D with increasing L_⊥.

(ii)
Next, to see how Φ^3D_y becomes singular as y → 0, we continue to keep (p-p_0^3D) and L_⊥ finite, and consider large but finite L, so that y = L_⊥ / L ≪ 1.
In this case, we expect to recover the scaling function near the 2D RBIM transition, namely

    = Φ^3D_y → 0[z] = Φ^2D [w = λ(L_⊥) · (p-p_0^2D(L_⊥)) · L^1/ν_2],

where λ(L_⊥) is an L_⊥-dependent multiplicative factor.
We emphasize that Φ^2D is an analytic universal scaling function, and can be extracted e.g. from logical failure rates of the 2D toric code with perfect measurements.
Noticing that

    (p- p_0^2D(L_⊥))· L^1/ν_2   = ((p- p_0^3D) - (p_0^2D(L_⊥) - p_0^3D )) · L^1/ν_2   = ((p- p_0^3D) - z_0 · L_⊥^-1/ν_3) · L^1/ν_2   = ((p- p_0^3D) · L_⊥^1/ν_3 - z_0) · L_⊥^-1/ν_3· L^1/ν_2   = (z - z_0) · y^-1/ν_2· L_⊥^1/ν_2-1/ν_3,

we may therefore define as our scaling variable

    w ≡ (z-z_0) · y^-1/ν_2 = λ(L_⊥) · (p- p_0^2D(L_⊥))· L^1/ν_2

where λ(L_⊥) = L_⊥^1/ν_3-1/ν_2,
and conclude that

    L_⊥ / L → 0=Φ^3D_y → 0[z = y^1/ν_3· x] 
              = Φ^2D[w = (z-z_0)· y^-1/ν_2],

see also Eq. (<ref>).
Based on Eq. (<ref>), we include in Fig. <ref> a schematic visualization of the function Φ^3D_y for small values of y as y → 0.







 §.§ Fluctuations in Size of Largest Rare Region


In the main text we found that this scaling function Eq. (<ref>) collapsed the data for the planar disordered L × L × L system, by setting L_⊥≡ln L.
We claimed that this worked because the physics of the disordered model was entirely controlled by the largest rare region which was of height (L_⊥)_ max∝ln L. However, there are expected to be subleading 𝒪(1) fluctuations in the height of the largest rare region. It is worth then asking how these fluctuations affect the mean and median .

In particular, we argue that scaling collapse of the mean is asymptotically destroyed by these small fluctuations in L_⊥ as L →∞; however, this only noticeably occurs in systems with at least hundreds of millions of spins. On the other hand, the median is more resilient with scaling collapse maintained at all sizes. The only caveat is that the form of the scaling variable w as a function of L at similarly large system sizes will need to change slightly to include subleading-in-L corrections to L_⊥.

For a fixed system size L, we may study the effect of a random L_⊥ through the scaling function Eq. (<ref>), where the analytic scaling function Φ^2D now has a random parameter w(L_⊥, L).
(For notational simplicity, we write Ψ for Φ^2D henceforth.)
We expand w around the mean of L_⊥ to the first order,

    w(L_⊥, L) = w(𝔼[L_⊥], L) + ϵ·∂ w/∂ L_⊥|_L_⊥ = 𝔼(L_⊥) + O(ϵ^2),

where ϵ≡ L_⊥ - 𝔼[L_⊥].
We define w_0 = w(𝔼[L_⊥], L), and write w' ≡ (∂ w / ∂ L_⊥)|_L_⊥ = 𝔼(L_⊥).
With these, we may write the mean failure probability as

    𝔼[ Ψ(w)] 
       = 𝔼[ Ψ(w_0 + ϵ· w') ] 
       =  𝔼[ Ψ(w_0) ] + 𝔼[ ϵ· w' ·Ψ'(w_0) ] +𝔼[ ϵ^2 · (w')^2 /2Ψ”(w_0) ] + …
       = Ψ(w_0) +  𝔼[ϵ] · w'·Ψ'(w_0) +𝔼[ ϵ^2 ] ·(w')^2/2·Ψ”(w_0) + …
       = Ψ(w_0) +𝔼[ ϵ^2 ]  (w')^2 /2Ψ”(w_0) + …
 
Here, we used 𝔼[ ϵ]=0, and we expect from standard extreme value statistics that 𝔼[ ϵ^2 ] = 𝒪(1); by 𝒪(1), we mean that it is asymptotically independent of system size. 
We expect Ψ”(w_0) to also be 𝒪(1), as Ψ is analytic and bounded.

Note in particular that 𝔼[ ϵ^2 ]  (w')^2 /2Ψ”(w_0) (and the higher-order corrections) will generically depend on parameters like L and error rates differently than w_0 depends on such parameters. This makes the expectation of the scaling function a sum of scaling functions with different scaling parameters, which will generically destroy single-parameter scaling if the corrections are not small.

Furthermore, we have that

    w' = A L^1 / ν_2 L_⊥^-1 / ν_2 - 1 + B L^1 / ν_2 L_⊥^-1/ν_2 + 1/ν_3 - 1.

Here, A and B are constants which depend on p - p_0^3D, z_0, ν_2 and ν_3.
Notably, w' diverges with L even when the scaling argument w is fixed and small.
Therefore, 𝒪(1) fluctuations in L_⊥ result in unbounded fluctuations in w, which will generically destroy the scaling collapse for the mean  when L →∞.

However, for L_⊥ = ln(L), w'<1 for L less than about 5000, so this destruction of scaling collapse only occurs at quite large sizes. We still see scaling collapse of the mean of  at the sizes of L that we probe (not pictured).

Because we do not expect single-parameter collapse for the mean in the thermodynamic limit, our plots of collapse are all for the median. We do expect good single-parameter collapse of the median of  at all L.

In particular, for a monotonic function like Ψ, 

    median(Ψ(w(L,L_⊥))) = Ψ(median(w(L,L_⊥)))

so collapse is maintained so long as the non-random median(w(L,L_⊥)) is used as the scaling variable.

As a last technical note, the form of median(w(L,L_⊥)) will slightly differ at sufficiently large sizes from what we use for collapse. We use w = (p-p_c^2d(L_⊥)) L^1/ν_2 L_⊥^1/ν_3 - 1/ν_2|_L_⊥ = c ln(L), taking L_⊥ directly proportional to ln(L) and neglecting the subleading 𝒪(1) corrections. At large sizes L, the multiplicative factor of L^1/ν_2 nevertheless makes the effect of these subleading corrections non-negligible; at sufficiently large L, single-parameter scaling collapse occurs for the slightly more complicated parameter (p-p_c^2d(L_⊥)) L^1/ν_2 L_⊥^1/ν_3 - 1/ν_2|_L_⊥ = c ln(L) + 𝒪(1) with explicit 𝒪(1) corrections included in L_⊥. However, this minor distinction only matters once the spins number in the hundreds of millions.





