Lifting topological codes: Three-dimensional subsystem codes from two-dimensional anyon models

Topological subsystem codes in three spatial dimensions allow for quantum error correction with no time overhead, even in the presence of measurement noise. The physical origins of this single-shot property remain elusive, in part due to the scarcity of known models. To address this challenge, we provide a systematic construction of a class of topological subsystem codes in three dimensions built from abelian quantum double models in one fewer dimension. Our construction not only generalizes the recently introduced subsystem toric code [Kubica and Vasmer, Nat. Commun. 13, 6272 (2022)] but also provides a new perspective on several aspects of the original model, including the origin of the Gauss law for gauge flux, and boundary conditions for the code family. We then numerically study the performance of the first few codes in this class against phenomenological noise to verify their single-shot property. Lastly, we discuss Hamiltonians naturally associated with these codes, and argue that they may be gapless.


I. INTRODUCTION
In order to reliably store and process quantum information, one needs to protect it against the deleterious effects of noise.Ideally, one would like to construct a self-correcting quantum memory [3][4][5], which is a quantum many-body system capable of passively protecting encoded information, analogously to magnetic storage for classical information.Such protection can be achieved by carefully engineering the Hamiltonian describing the system and the coupling of the system to a thermal bath [6][7][8].Unfortunately, it is not known how to realize a self-correcting quantum memory in three or fewer spatial dimensions without making strong assumptions about symmetries of the system and noise [9][10][11].
A typical solution to the problem of protecting quantum information invokes techniques of quantum error correction (QEC) [12,13], whereby one encodes information into some QEC code.By measuring certain operators, commonly referred to as parity checks, one can then detect and, subsequently, correct errors afflicting the encoded information.Topological codes constitute a particularly appealing class of QEC codes, which can be realized by arranging qubits on some lattice and measuring parity checks that are geometrically local.
Most prominent examples of topological codes include the toric code [3,14] and color code [15], which can be readily implemented in two spatial dimensions with currently pursued quantum hardware, such as superconducting circuits [16,17] or Rydberg atoms [18,19].
In the presence of unavoidable measurement errors, performing reliable QEC becomes more intricate and resource-intensive.With topological codes in two or three spatial dimensions, one usually resorts to performing repeated rounds of measurements in order to gain confidence in their outcomes [3,20].In particular, the number of measurement rounds needs to scale with the system size, unless one is willing to tackle a serious engineering challenge of measuring non-local operators of high weight [21][22][23].As a result, the time overhead of QEC increases and QEC itself becomes susceptible to time-correlated noise.
One can reduce the time overhead and achieve resilience to time-correlated noise by considering unorthodox topological codes, such as the recently-introduced three-dimensional subsystem toric code [2] and gauge color code [24][25][26], that allow for single-shot QEC [27].Roughly speaking, the single-shot QEC property asserts that, even in the presence of measurement errors, one can perform reliable QEC with a topological code by only performing a constant number of measurement rounds of geometrically-local operators.However, the origins of the single-shot QEC property remain elusive, with only few topological codes demonstrating it.
In this article, we provide a systematic construction of a novel class of QEC codes, which we refer to as subsystem abelian quantum double (SAQD) codes.This class can be viewed as a generalization of the subsystem toric code and is also closely related to the gauge color code.Roughly speaking, our construction allows to take any abelian quantum double model (that is natively defined in two spatial dimensions) and realize a corresponding topological subsystem code in three spatial dimensions.In addition, we investigate topological subsystem codes from the perspective of reliable quantum memories capable of single-shot QEC.In particular, we describe computationally-efficient decoding algorithms and benchmark their performance against phenomenological noise.We report memory thresholds around 1%, which are competitive with that of the surface code [28,29] given the advantage of single-shot QEC.Lastly, we argue that the Hamiltonians that naturally arise from topological subsystem codes are probably at a phase transition, providing an indication that these models will not serve as a self-correcting quantum memory.
We remark that the existence of the qudit generalization of the subsystem toric code is perhaps not surprising.However, our construction offers a radically different perspective on topological subsystem codes as constructed from (and thus closely related to) models natively defined in fewer spatial dimensions.Such a perspective provides a new way to understand possible boundary conditions and sheds light on the single-shot QEC property that topological subsystem codes exhibit.It also opens up possibilities of further generalizations, perhaps leading to concrete non-abelian models in three spatial dimensions.
The remainder of this article is structured as follows.In Section II, we discuss preliminaries, including the subsystem toric code and quantum double models.Then, in Section III, based on the intuition from the subsystem toric code, we introduce a class of subsystem abelian quantum double (SAQD) codes.In Section IV, we shift gears and numerically study QEC properties of SAQD codes, in particular we estimate memory thresholds for phenomenological noise.Finally, we discuss thermal stability of SAQD codes in Section V and provide concluding remarks in Section VI.We also enclose a number of appendices.In Appendix A, we define Pauli subsystem codes, their logical operators and code parameters.To ensure clarity, we provide a complete, and explicit, list of gauge generators for the SAQD in Appendix B. Finally, in Appendix C, we discuss the spectral gap of the SAQD Hamiltonian.

II. PRELIMINARIES
In this section, we provide some preliminaries required for the remainder of the article.We begin by briefly describing the subsystem toric code, followed by reviewing two-dimensional abelian quantum double models.

A. Intuition from the subsystem toric code
Stabilizer [30] and subsystem [31,32] codes constitute two important and widely-studied classes of QEC codes.A stabilizer code is defined by a stabilizer group S, which is an abelian subgroup of the Pauli group P and does not contain −I.By definition, logical qubits (or qudits) are encoded in the code space that is the +1 eigenspace of all the operators in S. A subsystem code is a generalization of a stabilizer code where logical information is encoded into only a subset of the logical qubits.The remaining (i.e., unused) logical qubits are commonly called gauge qubits.Subsystem codes can be specified by a gauge group G, which can be any subgroup of P. By definition, any operator in G has a trivial action on the encoded logical information.The stabilizer group of a subsystem code is the center of the gauge group G (up to phases).For a detailed discussion, we refer the reader to Appendix A. The three-dimensional subsystem toric code is a generalization of the stabilizer toric code in an analogous way as the gauge color code is a generalization of the stabilizer color code.A conventional way to understand both the subsystem toric code and the gauge color code is to relate them to the corresponding stabilizer codes in the same spatial dimension.Any state in the code space of the stabilizer code can be viewed as a code state of the subsystem code [33], and a code state of the subsystem code can be transformed into a code state of the stabilizer code by fixing the gauge qubits to be in a specific state.More formally, we have the inclusions where S stab and S denote the stabilizer groups of the stabilizer and subsystem toric codes respectively, and G denotes the gauge group of the subsystem toric code.Any state in the +1 eigenspace of S stab will be in the +1 eigenspace of S. Since S is the center of G, given any operator σ ∈ S stab such that σ / ∈ S, there must be an operator τ ∈ G that anticommutes with it.Consequently, any state in the +1 eigenspace of S can be transformed into a state in the +1 eigenspace of S stab by applying an operator (namely τ ) in G, which, by definition, has no effect on the encoded information.
It is non-trivial to construct subsystem codes from the corresponding stabilizer codes, as evidenced by the amount of time between the introduction of the three-dimensional stabilizer toric code [3] and the threedimensional subsystem toric code [2].Although the stabilizer and subsystem codes have a similar construction, they have strikingly different properties.For example, the aforementioned subsystem codes have the singleshot QEC property, whereas the corresponding stabilizer codes do not.And despite these subsystem code models arguably providing the most transparent illustration of the single-shot QEC property, the origins of this property, and its exact connection to self-correction, are unclear.
The subsystem toric code can be defined on the cubic lattice with open boundary conditions by placing qubits on the edges of the lattice, as depicted in Fig. 1a.By definition, each volume of the lattice supports one stabilizer generator and eight gauge generators, whose type is determined by the color associated with the volume.The product of gauge generators around each volume is the identity operator; similarly for gauge generators surrounding each vertex.These relations are key to the single-shot QEC property of the subsystem toric code, as they ensure that -1 measurement outcomes of gauge operators form closed loops.If measurement errors cause some of the loops to be broken, then we can close them up, repairing the damage and obtaining an updated estimate for the stabilizer syndrome used in a standard QEC decoder.
Note that the gauge generators of the subsystem toric code surrounding each vertex can be viewed as the stabilizer generators of the toric code supported in a twodimensional sphere around that vertex (Fig. 1b).Moreover, the whole gauge group can split into two commuting sets of operators, each of which is the collection of disjoint, two-dimensional toric codes.This simple observation is the starting point of our construction of threedimensional subsystem codes based on arbitrary abelian quantum double models which we present in Section III.

B. Two-dimensional quantum double models
To any finite abelian group there is an associated topological code, called the abelian quantum double model (AQD).The most familiar example is the toric code, which is the special case associated to Z 2 .For the remainder of our discussion, we specialize to the case G ≃ Z d as all other cases can be obtained via stacking layers, with one layer for each factor in Eq. ( 2).Let Λ be a directed graph embedded on an orientable 2-manifold.For simplicity, we will initially assume that Λ is closed.Degrees of freedom of dimension d are located on the edges of Λ.In the group basis the generalized Pauli operators act as where ω := exp(2πi/d).The commutation relations are The Hamiltonian defining the model has the familiar form For a fixed vertex v, the A v term consists of Z operators applied to each incoming edge of Λ, and Z † to each outgoing edge.For a fixed face f , the B f term consists of X operators applied to each edge whose direction matches the ambient orientation, and X † applied to all edges opposite to the orientation.For example, with analogous action on vertices and faces of arbitrary degree.These rules ensure that the Hamiltonian terms pairwise commute.Ground states of the model are simultaneous +1 eigenstates of all Hamiltonian terms.Interpreting a state |g⟩ on an edge as the presence of a string of color g, the vertex terms ensure strings cannot terminate or change color (without branching), while the face terms cause a given ground state to be a superposition of all homotopic, closed loops [34].On a 2-sphere, the ground state is unique, while on higher genus surfaces, the ground space is degenerate.There are two kinds of excitations above the ground state, denoted e g and m g , corresponding to violating vertex and face terms respectively, which are created by local operators.When brought close together, these particles fuse according to the group product, namely e g × e h = e g+h mod d and m g × m h = m g+h mod d .In general, any configuration of excitations on a closed manifold must fuse to vacuum e 0 m 0 since all local operators preserve the total charge, and the ground state is neutral.Unlike the toric code, for which G ≃ Z 2 , typical excitation configurations 'branch'.For example, whenever G ̸ ≃ Z 2 we can have the following configuration of excitations 1. Boundaries In addition to closed manifolds, AQD models can be constructed on manifolds with boundary.In turn, we will make use of these boundaries to construct threedimensional subsystem codes with boundary.For a general abelian group G, there are many possible gapped boundaries that can be used to modify Eq. (6), see for example Refs.[35][36][37].For our code construction, we will make use of the only two types of boundary that occur irrespective of the choice of G.These boundaries are commonly referred to as rough and smooth due to their lattice implementation.
Rough boundaries are characterized by all vertices of Λ that lie on the boundary having degree one.The bulk face terms of the Hamiltonian are replaced by terms of the form where the dashed line is the physical boundary.At smooth boundaries, all vertices have degree larger than one.The bulk vertex terms are replaced by terms of the form Strings of X operators can terminate at rough boundaries without creating excitations, while strings of Z can terminate at smooth boundaries.Physically, this ability to remove an excitation onto the boundary is commonly called condensation.In that language, the rough boundary condenses e type excitations, while the smooth boundary condenses m particles.

Interfaces
Boundaries to vacuum can be understood as a particular example of an interface between two (possibly identical) AQDs.The most familiar example is the symmetry of the toric code which exchanges the e and m excitations [38], however, even in AQD codes, they can be far richer [39].In two-dimensional codes, the possible self-interfaces (domain walls) classify the fault-tolerant logical operators [40][41][42][43].
Beyond domain walls, one can also consider interfaces between domain walls [38,44,45], which can be used to further boost the computational power of twodimensional topological codes.We have not explored their impact on the subsystem codes we introduce here.

III. SUBSYSTEM ABELIAN QUANTUM DOUBLE CODES IN THREE DIMENSIONS
We now proceed to define subsystem abelian quantum double (SAQD) codes in three spatial dimensions.We then explain how boundaries and interfaces can be obtained naturally from our construction.We finish by discussing the gauge flux in the SAQD codes and their relation to the qudit toric code.For simplicity, we begin with a concrete lattice implementation, however the code can be defined more generally (Section III 6).

Bulk
Our construction of the gauge group follows the following intuition: 1. Tile three-dimensional space with topological codes defined on 2-spheres.
3. Identify the degrees of freedom at the intersections.
Step 1 produces a stabilizer code which encodes no qudits.Steps 2 and 3 convert this stabilizer code into a subsystem code.
Let L be the cubic lattice, with volumes colored red and blue in a checkerboard pattern, and vertices colored green and yellow in a similar pattern (see Figs. 1a  and 2a).As mentioned in Section II B, defining an AQD requires a directed graph.The volume coloring is an auxiliary tool enabling consistent definition of these graphs throughout the code as follows.Connect the centers of nearest neighbor blue volumes of L with directed edges, as depicted in Fig. 2a, with the direction of each edge chosen arbitrarily.Qudits of the final code ultimately reside on the added, blue edges [46].With this auxiliary lattice in place, we are in a position to define the SAQD gauge group.
At each vertex of L, we place a microscopic (i.e., not scaling with code size) sphere, which we imagine to be hosting a copy of the appropriate AQD (Step 1).As in Fig. 2a, the auxiliary blue lattice induces a directed triangulation on these spheres, allowing an explicit lattice model.In addition to the directed graph, it is necessary to provide an orientation to the 2-spheres, which can be achieved by defining a vector from green vertices of L to yellow vertices.We now inflate these spheres until their triangulations intersect with the auxiliary lattice (Step 2).Each sphere is split by six directed edges into four triangular regions and thus is homeomorphic to the boundary of a tetrahedron.Finally, the qudits on the common edges are identified (Step 3), leading to the geometry depicted in Fig. 2b.Qudits reside on the dark black edges which form the triangulation in Fig. 2b.
With the geometry established, we aim to use the topological features of the AQD code to ensure the resulting code has the single-shot QEC property.The gauge group of the resulting code is given by where A v , B f are defined in Eq. ( 7).Gauge operators of different colors do not commute, so this fails to be a stabilizer code.Although the model can no longer be thought of as a collection of unbroken loops as discussed in Section II B, it remains true that the total charge on a sphere must be trivial.
In the remainder of this section, we discuss the stabilizers of the code, and how boundaries can be introduced.As in the gauge color code and subsystem toric code, the latter are necessary to ensure a non-trivial codespace.We show numerical evidence that the code has the single-shot property in Section IV.
A full accounting of the gauge generators is provided in Appendix B for completeness.

Stabilizers
Since the underlying topological codes are commuting, and distinct spheres of the same color are nonintersecting, the stabilizer group corresponding to the gauge group Eq. ( 11) is given by any operator that can be created from generators of either color.One can verify that locally generated stabilizers are associated with the volumes of the cubic lattice.Recall that we began with as bicolored cubic lattice L. Red volumes only contain faces of tetrahedra, such as that in Fig. 3c, and consequently support X-type stabilizer operator.Conversely, blue volumes contain vertices of tetrahedra such as those depicted in Fig. 3d, and support Z-type stabilizers.
In addition to the local stabilizers, the model on the 3-torus (T 3 ) supports sheet-like stabilizers in each coordinate plane as depicted in Eq. (B4).By introducing boundaries, these stabilizer operators can be promoted to logical operators.

Boundaries
Due to the sheet-like stabilizer operators on a closed manifold, boundaries are required to have a non-trivial codespace This is also seen in the subsystem toric and gauge color codes.This suggests that their nature, although not truly topological, is closely related to models in one fewer spatial dimension.
We will discuss two classes of boundaries for the SAQD, both inherited from the underlying AQD: I. Boundaries given by macroscopic copies of the underlying AQD (as opposed to the microscopic 2spheres at the vertices of L).
II. Boundaries of the subsystem code induced by the one-dimensional boundaries of the AQD.
The type of boundary (Class I) we consider is obtained by defining the code on T 2 × I, where I := [0, 1], as shown in Fig. 4a.Gauge generators associated to the boundaries T 2 ×{0} and T 2 ×{1} (referred to as bottom and top respectively) are those of the underlying AQD model.Unlike the spheres in the bulk, these surfaces are macroscopic (Fig. 4b).To encode qudits, it is necessary to color the boundaries top and bottom oppositely.As in Ref. [2], we need to introduce additional qudits as discussed in Appendix B to make this possible.
Orientation of the top and bottom surfaces is defined in the same way as the bulk spheres (Section III 1).The orientation vectors point from green to yellow.As an example, the top surface, colored green, has the gauge generators shown in Fig. 4c.

Open boundaries
The other type (Class II) of boundaries we consider is obtained by introducing boundaries on the 2-spheres of AQD, as introduced in Section II B 1.This allows definition of the code on more general manifolds (with boundaries), such as the cube.For the remainder of this section, the top and bottom boundaries will be macroscopic AQD models.
As discussed in Section II B 1, all AQD models can be equipped with rough and smooth boundaries.These boundaries of the two-dimensional models can be used to define boundaries of the associated three-dimensional SAQD code.When a green or yellow sphere intersects with a boundary, we can choose to equip the code on that sphere with rough boundaries, as shown in Fig. 5a.
We refer to such a boundary as a rough boundary of the SAQD.Conversely, we could choose smooth boundaries for all spheres that intersect the boundary, as in Fig. 5b, defining a smooth boundary for the SAQD [47].
At the y-edges of the cube, where different boundaries meet, the spheres carry mixed boundary conditions.Where the rough/smooth boundaries intersect the top and bottom boundaries, the macroscopic AQD inherits the appropriate class of boundary.

Three-body model
By construction, all bulk gauge generators, and those associated to rough and smooth boundaries have nontrivial support on at most three qudits.The macroscopic AQD models on the top and bottom boundary introduces four body terms, however this is not essential.It may be advantageous from an error correction perspective to reduce these to weight three [3].In Appendix B 4, following Ref.[2], we show how these boundaries can be modified to reduce the weight.

Logical operators
When defined on T 3 , the SAQD code has both local and macroscopic (sheet-like) stabilizer operators.Recall that the purpose of introducing boundaries was to promote these sheet-like stabilizers to logical operators.Sheets that do not intersect the boundary remain as stabilizers.For example, on T 2 × I, the sheet-like operators in the xz-plane (Fig. 4d) are in G Z d , so remain stabilizers.Conversely, sheet-like operators that do intersect the boundary can no longer be constructed from gauge generators, but remain in the centralizer of the gauge group.On T 2 × I, this gives sheet-like logical operators in both the xy and yz-planes, meaning the code supports two logical qudits.
In the bulk, these bare logical operators can be dressed with either green or yellow gauge generators to reduce their weight.Although these dressed logical operators no longer lie in the centralizer of the gauge group, they do commute with the stabilizer, and so remain good logical operators.By dressing the sheet-like logical operators with green generators, the sheet is reduced to a string-like operator on the green boundary, and vice versa.These strings correspond to the logical operators of the boundary AQD as shown in Fig. 4e.For a more explicit discussion of the logical operators we refer the reader toAppendix B.

Code parameters
On a closed manifold, SAQDs do not encode any qudits due to the macroscopic stabilizer generators.When placed on a manifold M × I (possibly where M has boundaries), with AQD models at the top and bottom, the codes have a number of logical qudits equal to the number for that AQD model on M .
For the manifolds discussed above, the code parameters are summarized in Table I, and obtained in Appendix B.

Interfaces
As discussed in Section II B 2, locality-preserving gates on topological codes in two spatial dimensions are in one-to-one correspondence with codimension-1 interfaces.Such gates can be lifted to the associated SAQD code, simply by applying the appropriate circuit to each microscopic or macroscopic copy of the AQD.Such an action will locally preserve the full gauge group, meaning its action can be restricted straightforwardly to the boundaries.Since the action on the dressed logical operators is identical to that in the input code, the logical action is also identical.Conversely, any non-Clifford gate that can be enacted on these SAQD codes cannot be a symmetry of the full gauge group in this way [48,49].In this work, we have not investigated the interpretation of higher codimension defects, such as point-like defects in two dimensions.This class includes both anyonic excitations, and more exotic objects such as twist defects [38].

General lattice implementation
In this section, we have restricted our discussion to the code defined on the cubic lattice, possibly with boundaries.As in Ref. [2], the definition can be extended to any colorable octahedral lattice (see Ref. [2] for the definition).Given such a lattice, we can connect the blue vertices as in Fig. 2 and inflate the green/yellow vertices to construct the code.

A. Physics of the gauge flux in the SAQD codes
We now briefly discuss the behavior of the gauge flux in the three-dimensional SAQD codes (which is qualitatively similar to that of the gauge color code [50,51]).For simplicity, we restrict our attention to the X-type gauge operators B f , in turn, this means we will only discuss red (X-type) stabilizer operators.We will refer to an outcome ω g when measuring B f , as a gauge flux of g flowing into the associated sphere.A flux of g −1 flowing in can be interpreted physically as g flowing out.
There are two colors of flux lines, green flows into green vertices and yellow flux flows into yellow vertices.Just like in the gauge color code, flux lines can split and fuse according to the product on the group G.

No errors
Since the gauge group is non-abelian, the outcome of any B f measurement appears random.This means that on any measurement round, the gauge flux configuration will be random, however there are relations: 1. Due to the topological order on the 2-spheres at each vertex, there can be no net charge on the sphere, and so the total flux in to any vertex must be 1 G .
2. Since every stabilizer operator can be formed from either only green B f or only yellow B f , the total green flux leaving any red volume must equal the yellow flux leaving the same volume.
3. By definition, the code space has all stabilizer operators equal to +1, so the total flux of either color leaving each red volume must be 1 G .
The first two relations, referred to as Gauss laws, are required to hold by construction, with the first being due to the two-dimensional topology and the second due to the three-dimensional geometry.

Physical errors
When there are physical errors acting on the code, the first two relations must still hold.Errors can be thought of as creating G charge at the red volumes, so can source flux, leading to violation of the third relation (Fig. 6a).This is then used to detect, and with the aid of a decoder to correct errors.

Measurement errors
When measurement errors occur, the resulting gauge flux configurations may violate any of the relations, as in Fig. 6b.Since relations [1] and [2] are necessarily true in the physical system, observed charge at the vertices, or mismatch of the yellow and green charge at a volume can be used to correct the defective measurements.
We discuss how the gauge flux can be used to correct both types of error in Section IV.

B. Relation to qudit toric codes
We conclude this section by briefly discussing the relationship between the SAQD codes and the qudit toric code in three spatial dimensions.For simplicity, we restrict our attention to the 3-torus.
Recall from Fig. 3 that the stabilizer group associated to the SAQD code is generated by operators associated to the volumes of the cubic lattice.
We can define another abelian subgroup of G via gauge fixing [24,52,53].This technique allows new stabilizer codes to be constructed from subsystem codes by choosing a maximal subset of commuting, local gauge generators to be the stabilizer group.Consider the subgroup of G Z d as the group generated by all locally generated stabilizers (associated to volumes) of Z-type and all X-type gauge operators.On the blue lattice (Fig. 2a), this gives Z-type stabilizers associated to vertices, and X-type stabilizers associated to faces.The code defined in this way deserves to be called the Z d toric code, since it reduces to the standard construction when d = 2.To the best of our knowledge, such a definition has not appeared explicitly in the literature.

IV. ERROR THRESHOLDS AND DECODING
In this section, we use Monte Carlo simulations to estimate the error threshold of the SAQD codes defined for Z d and supported on I 3 (with rough/smooth boundaries as discussed above) as a function of d.The simulation code is available on GitHub [54].

A. Setup and error model
The SAQD codes for Z d are CSS [55,56], meaning that we can choose a generating set for the gauge group where each generator consists of either Pauli X-type operators or Pauli Z-type operators only.Consequently, we can decode Pauli X-and Z-type errors independently and we therefore restrict our attention to Z-type errors.Specifically, we consider single-qudit error channels of the following form Error channels of this form have been used previously in studies of qudit topological codes [57][58][59][60][61].
We define a QEC cycle to consist of qudit Z-type errors, noisy X-type gauge flux measurements and the application of a recovery operator produced by a decoding algorithm (that only uses measurement outcomes from that QEC cycle).We simulate the above with the following sequence of steps.
1. Update the residual error by applying a Z-type gauge operator chosen uniformly at random from the full gauge group (not just from a set of generators).
2. Update the residual error by applying the error channel in Eq. ( 12) to each qudit in the model.
3. Calculate the X-type gauge flux.
4. For each flux measurement, with probability p return the outcome chosen uniformly at random from Z d \ {x}, where x ∈ Z d is the correct outcome.
5. Update the residual error with the qudit recovery operator output by the decoder described in Section IV B.
We apply the random gauge operator in step 1 above as we are implicitly assuming alternating rounds of X-type and Z-type gauge operator measurements.We begin the simulation by initializing the residual error to be trivial and we then simulate t QEC cycles.After the final QEC cycle we model readout by a round of ideal (no measurement error) decoding.

B. Clustering decoder
We break the decoding problem into two stages: (a) flux validation and (b) qudit correction.In the first stage we attempt to correct the measurement error using the local relations between gauge operators.And in the second stage we attempt to find a (qudit) recovery operator using the stabilizer eigenvalues, which we construct from the corrected flux measurements produced in the first stage.Both decoding stages can be captured by the following decoding problem.Let H ∈ F m×n d be a parity-check matrix, where each row is a parity check and the each column corresponds to a variable.We are given an error syndrome σ ∈ F m d and our task is to output a recovery operator y ∈ F n d such that Hy = σ (and we may further require that y has minimal Hamming weight).In the first decoding stage the variables are the flux measurement outcomes and the checks are the local relations, while in the second stage the variables are the qudits and the checks are the stabilizers.
In the Z 2 case, both decoding problems can be solved efficiently and exactly by the minimum-weight perfect matching algorithm [2], as in this case each variable is in the support of at most two parity checks.However, for qudits of local dimension d > 2 the generalization of matching is inefficient and we therefore use a clustering decoder [57,62].We note that two-stage decoders such as the one we have described have previously been used to decode a variety of quantum codes with the singleshot property [2,63,64].
To explain how the clustering decoder works, we review the construction of a syndrome graph G from a parity-check matrix H with the property that each column contains at most two nonzero entries (the paritycheck matrices for both decoding stages have this property).For each row of H we create a node in G and we connect two vertices if their corresponding rows in H have overlapping support.We also create an auxiliary node (representing the physical boundary of the model) and for any column with only one nonzero entry we connect the corresponding node the auxiliary node.
A given error syndrome σ can be interpreted as specifying the charge present at each node in G.We define a cluster to be a subgraph of G with the structure of a tree.The charge of a cluster is the product (in G) of the charges of its nodes.We define a cluster to be neutral if and only if either its charge is 1 G or it contains the auxiliary node (representing the physical boundary).The first step of the clustering decoder is to initialize a list of clusters, where there is a cluster for each of the nodes in G with non-trivial charge.The decoder then iterates the following three subroutines until the list of clusters is empty.
1. Grow: For each non-neutral cluster, iterate over its leaves as follows.For each neighbor of the leaf, check if the neighbor is in the cluster.If not, add it to the cluster and: (i) If the neighbor is the auxiliary node then stop growing the cluster and mark it as neutral.(ii) If the neighbor belongs to a different non-neutral cluster, then mark the clusters to be merged.If the combined charge of the clusters to be merged is neutral then stop growing the current cluster.
2. Merge: For each pair of clusters marked for merging, attach the smaller cluster as a subtree of the larger cluster.
3. Neutralize: For each neutral cluster compute a correction by moving all the nonzero charges to the root (if the charge of the cluster is trivial) or otherwise by moving all the nonzero charges to the auxiliary node.After neutralizing a cluster, remove it from the list of clusters.
We remark that the performance of our clustering decoder can be further improved by optimizing the "Growth" and "Merge" steps.In particular, we could adapt them to resemble the corresponding steps of the union-find decoder [65].Given that the path from each charge to the root of its cluster can be updated during the "Growth" and "Merge" steps, the runtime analysis for the union-find decoder is still applicable, therefore implying that our optimized decoder can be implemented in almost linear time.We emphasize however that the standard version of the union-find decoder is only applicable in our setting for the Z 2 model.

C. Simulation results
We investigate the behavior of the error threshold of the model as a function of the local dimension d and the number of QEC cycles t.We find that the error threshold is insensitive to the number of QEC cycles, but increases with the local dimension d, in agreement with previous results [57][58][59][60][61].The results are shown in Fig. 7.
We remark that the error channel in Eq. ( 12) may bias the results in favor of large d, as we are squeezing an ever larger number of error events into the same probability interval.As an alternative, we could consider a model, where p increases with d.For example, we can view m qubits as an effective qudit of dimension d = 2 m .According to Eq. ( 12) the probability of no error is 1 − p(d), where we have now included a dependence on d.On the other hand, the probability of no error for m qubits is (1 − p(2)) m .Equating these two values gives In Fig. 7d, we rescale the error thresholds according to the above equation, which dramatically changes the 0.2% 0.6% 1% 1.4%  Although the logical error rates increase, the error thresholds are essentially the same.(c) A plot of the error thresholds, p th , as a function of the number of QEC cycles for different values of d.The error thresholds are estimated from data such as those shown in (a) and (b) using finite-size scaling analysis [67], and the error bars are derived using bootstrapping with 100 estimates.(d) A plot of the rescaled error thresholds, p * th , as a function of the number of QEC cycles, where the error thresholds from (c) are rescaled according to Eq. (13).In this case the advantage for larger d disappears.
results.There is a clear degradation of performance for large d, giving the opposite conclusion to that of Fig. 7c.This suggests that the increased performance of qudit topological codes for large d reported in previous works may in fact be an artifact of the error model used in simulations.
The error threshold we find for the Z 2 case is very close to the error threshold found for the gauge color code in Ref. [63] using a similar decoder.However in that case the error threshold decays from a higher starting value whereas we observe no meaningful change in the error threshold as we vary the number of QEC cycles (we discuss a possible reason for this effect in Section IV D).

D. Decoder comparison
For the Z 2 model, we compare the performance of four different decoders, where for each of the two decoding stages we use either the clustering decoder described in Ref. Section IV B or a matching decoder [2].Using different decoders for each stage allows us to investigate which stage limits the overall performance.For each decoder we estimated the error threshold for t = 4 QEC cycles, the results are shown in the following The results suggest that the correction stage is the bottleneck of the decoder, as changing the validation decoder has a comparatively small impact on the performance.This makes intuitive sense as the flux validation stage uses measurement outcomes that need to satisfy local symmetries captured by the Gauss law, which are not present at the qudit correction stage.Moreover, this may also explain the fact that the error threshold does not seem to depend on the number of QEC cycles in Fig. 7c.Lastly, we emphasize that our implementation of the decoders for the SAQD codes is relatively naive.For instance, our results could likely be improved by optimizing the clustering decoder using a similar approach as previous work on decoding the two-dimensional qudit toric code [57][58][59][60][61].It may also be possible to improve the performance by adapting the "single-stage" decoder described in [68] to the case of subsystem codes.Another possibility of improving the performance is to view the SAQD codes through the lens of Floquet codes [69].In particular, we made an implicit assumption on the measurement schedule to comprise measuring all Pauli X-type operators, followed by measuring all Pauli Ztype operators.Considering other (interleaved) measurement schedules is likely to boost error thresholds, as numerically demonstrated with two-dimensional topological codes [70].

V. PHASES AND TRANSITIONS OF AN SAQD HAMILTONIAN MODEL
The question as to whether topological subsystem codes in three spatial dimensions are thermally stable is an important open problem [4,24,27,[71][72][73].Given a subsystem code, we can associate to it a Hamiltonian, and then a natural question related to its thermal stability is whether the associated Hamiltonian is gapped or gapless.In particular, the Hamiltonian being gapless would seem to preclude it from being thermally stable.
In this section, we argue that the most natural Hamiltonian associated to the SAQD codes is likely to be at a phase transition.Our argument follows that of Kramers and Wannier in Ref. [74].In particular, we provide a self-duality on the Hamiltonian, and show that there are at least two phases.If there is a single transition, it must occur at the point most naturally associated to the code.
Unlike the case of stabilizer codes, associating a Hamiltonian to a subsystem code is not unique.In this work, we have taken the perspective that the green and yellow AQD codes should form the building blocks of the code.It is therefore natural to define the code Hamiltonian associated to a group G to be corresponding to Eq. ( 6) on the green and yellow spheres with weights J g and J y respectively.Since, from the error correction perspective, there is no preference of either color, the code is associated to H(1, 1).From this perspective, the stabilizer operators in Fig. 3 are (1-form) symmetries of the Hamiltonian.In Appendix C, we show that Eq. ( 14) is local-unitary equivalent to a Hamiltonian of the form where are reminiscent of the Raussendorf-Bravyi-Harrington (RBH) model [75].
Eq. ( 15) has a self-duality map (see Appendix C) exchanging J g ↔ J y .When (J g , J y ) = (1, 0), the Hamiltonian is paramagnetic.Conversely, at the point (0, 1), there is non-trivial 1-form (Z 2 × Z 2 )-protected topological order, where the symmetries are generated by the images of the stabilizers.
From this self-duality, we can deduce that if the model has a single phase transition, then it must occur where J g = J y .If this transition is continuous, it follows that H(1, 1) is gapless.It could also be that the transition is first order, as in the symmetric sector [76], in which case the gap would remain open.
It has previously been noted that thermal stability of the RBH model with 1-form symmetries may be closely related to the single-shot properties of associated subsystem codes [2,77].Our perspective suggests that the relevant symmetries are the stabilizer operators of the subsystem code.

VI. DISCUSSION
In this article, we provided a systematic construction of the three-dimensional SAQD codes, which provide a natural generalization of the subsystem toric code.Although we focused on the cubic lattice, the SAQD codes can be realized on colorable octahedral lattices [2].It is currently unclear which aspects of these lattices are necessary and which are merely sufficient.
Our construction required inflating the vertices of the lattice to 2-spheres on which the two-dimensional topologically-ordered model was placed.The choice of spheres does not seem to be essential, rather one could consider inflating to more interesting 2-manifolds.In such a setting, the constraints on the gauge flux described in Section III A should remain in place, however other code properties may be more exotic.In order to ensure that the SAQD codes exhibit the single-shot QEC property, we had to carefully choose the geometry used to combine these 2-spheres.Manifolds with nontrivial genus would allow exploration of more interesting geometries, such as linking.
Central to our construction were topologicallyordered models in two spatial dimensions, which we assumed to be untwisted abelian quantum double models.It is natural to ask whether our construction can be extended to all abelian anyon models, or general Levin-Wen models.From a physical viewpoint, such a generalization may lead to novel phases, possibly without a gap.This raises the question of classifying the SAQD codes and other topological subsystem codes in three spatial dimensions, and whether there is a connection with the classification of two-dimensional topological orders.To properly formulate this question, we require a notion of equivalence between (topological) subsystem codes, which generalizes local unitary equivalence of topological subspace codes [78].Beyond pure-state properties, the authors of Ref. [79] have defined a notion of mixedstate topological quantum field theory (TQFT), and established a connection to the two-dimensional Levin-Wen models based on modular fusion categories.It would be interesting to explore this connection, and understand whether the SAQD codes give more examples of mixed-state TQFTs.
The gauge color code appears to fit into our perspective, but not directly into our construction.Around each volume of the gauge color code, the gauge generators form a copy of the stabilizer color code model in two spatial dimensions.This is reminiscent of the AQD models around vertices in our construction, and would suggest the gauge color code should be equivalent to the SAQD associated to Z 2 × Z 2 .This appears to be in slight conflict with the stabilizer case.In three spatial dimensions, the stabilizer color code (without boundaries) is locally-unitarily equivalent to three decoupled copies of the stabilizer toric code [80].We suspect that resolving this disparity would shed light on both the appropriate notion of equivalence and the necessary threedimensional geometrical properties.
The SAQD codes constitute an excellent illustration of the single-shot QEC property and physics behind it.Moreover, they can provide illuminating insights into a connection between single-shot QEC and self-correcting memories.We have argued that the most naive Hamiltonians associated to the SAQD codes are probably gapless, but are closely related to models with symmetryprotected self-correction [9][10][11].Given this close connection, one might ask whether the subsystem code Hamiltonian could be modified to energetically enforce the symmetry, possibly leading to a self-correcting quantum memory.Colors are reversed at the bottom boundary.

T 2 ×[0, 1]
We now consider partially open boundary conditions.In the bulk, nothing is changed from Appendix B 1. At the boundaries T 2 ×{0} and T 2 ×{1}, we place 2-dimensional AQD models of opposite colors.To achieve this, we need to modify the lattice at the boundary, adding qudits on 'half-edges' as indicated on the top surface of Fig. B1.In that figure, new spheres have been attached to the top (all yellow) and bottom (all green).Unlike those in the bulk, these are only in contact with two lattice edges, so are only decorated with two oriented edges.We choose the following everywhere for simplicity where we remark that the orientation bottom edge of the yellow spheres is dictated by the choice on bulk spheres, and similarly for the top edge of the green spheres.Inclusion of these additional spheres, and their AQD models, leads to the additional gauge generators With these extra qudits in place, we can assign operators to the top and bottom faces.These are colored with the opposite color to the spheres incident on the face completing the gauge generators for the model on T 2 × I. Compared to T 3 , there are L 2 additional qudits contributed by not identifying the top and bottom edges.Additionally, there are L 2 /2 extra qudits introduced on both the bottom and top faces.Therefore, there are n = 3L 3 + 2L 2  (B11) qudits on the lattice.
which act on the Z d Pauli operators as and reduces to the CNOT gate when d = 2.Additional qudits are placed on half of the plaquettes on the top and the bottom faces, for example on the top face where the gray dots indicate the added qudits, which we refer to as plaquette qudits.For each of the L 2 plaquette qudits, we add two gauge generators, X and Z.The transformed code therefore has The circuit is translationally invariant, acting on a unit cell as where the gates are applied in the order indicated, with all '1' operators applied first, followed by all '2' operators, and so on.Note that the gates act only on the qudits sticking out of the surface, and therefore any gauge operator that does not have support on these sites is left unchanged.The altered gauge generators on the top are

FIG. 1 :
FIG. 1: (a)The three-dimensional subsystem toric code can be defined on the cubic lattice.Qubits are placed on edges.Pauli X-and Z-type gauge generators are depicted as triangular red and blue faces, respectively, whereas stabilizer generators are associated with volumes.(b) Gauge generators can be partitioned into two commuting sets, associated to green and yellow vertices respectively.Around each vertex, the gauge generators form a 2-dimensional sphere supporting a toric code.

FIG. 2 :
FIG. 2: Given a cubic lattice with checkerboard red and blue volumes, and bicolored (green and yellow) vertices, the code lattice is defined by: (a) Joining the centers of all blue volumes by directed edges.Code qudits are associated to the blue edges.(b) Inflating the vertices to form 2-spheres until they touch, with directed triangulation induced by the blue edges.An interactive version of this figure is available [1].

3 :
FIG. 3: Stabilizer operators are those that can be recovered from gauge generators of either color.They are associated to the volumes of the cubic lattice.(a) and (b) show the faces participating in a given X-type stabilizer.In (c) and (d), we show how to recover the stabilizers by multiplying certain gauge generators in two different ways.An interactive version of this figure is available [1].

FIG. 4 :
FIG. 4: T 2 × I with class I boundaries.2-dimensional AQD models are on the top and bottom boundaries, with lattice periodic in other directions.(b) shows the view from above.All spheres incident on the top boundary are colored yellow, and the AQD model on the surface is colored green.Colors are reversed at the bottom boundary.(c) Gauge generators assigned to the top surface, realizing a class I boundary.To aid discussion, we fix the coordinate axes in (d) for the remainder of the manuscript.(e) Logical operators on the SAQD are sheet-like operators that intersect the boundary.They can be dressed to string-like operators supported on the boundary, corresponding to the logical operators of the boundary AQD.

FIG. 5 :
FIG. 5: Class II boundaries.At the rough (a) and smooth (b) boundaries, the half spheres are decorated with a quantum double model with the appropriate type of boundary.An interactive version of this figure is available [1].

FIG. 6 :
FIG. 6: Gauge flux configurations.Unlike the subsystem toric code, but similarly to the gauge color code, gauge flux can split for general SAQDs.In addition, the flux lines can carry g ∈ G units of flux (indicated by shading).(a) Typical gauge flux configuration in the absence of measurement errors.Flux obeys a Gauss law defined by the input AQD model.(b) When measurement errors occur, flux lines do not obey the required Gauss law.An interactive version of this figure is available [1].

FIG. 7 :
FIG. 7: (a)A plot of the logical error rate, p fail , as a function of the qudit and measurement error rate, p, over t = 2 QEC cycles for various Z d SAQDs.The error bars show the Agresti-Coull 95% confidence intervals[66].For each value of d, the error threshold is the point where the curves for different L intersect.(b) A plot analogous to (a) but for t = 16 QEC cycles.Although the logical error rates increase, the error thresholds are essentially the same.(c) A plot of the error thresholds, p th , as a function of the number of QEC cycles for different values of d.The error thresholds are estimated from data such as those shown in (a) and (b) using finite-size scaling analysis[67], and the error bars are derived using bootstrapping with 100 estimates.(d) A plot of the rescaled error thresholds, p * th , as a function of the number of QEC cycles, where the error thresholds from (c) are rescaled according to Eq.(13).In this case the advantage for larger d disappears.
FIG.B1: T 2 × 2-dimensional AQD models are on the top and bottom boundaries, with lattice periodic in other directions.We regard the 'dangling' edges as having length 1/2, meaning this lattice has linear size 6 × 6 × 6.(b) shows the view from above.All spheres incident on the top boundary are colored yellow, and the AQD model on the surface is colored green.Colors are reversed at the bottom boundary.

TABLE I :
Code parameters for the SAQDs on various manifolds of linear size L, possibly with boundary.The cube (I 3 ) is equipped with rough/smooth boundaries, while (T 2 × I) ′ is the code with only 3-body generators as discussed in Section III 2. The number of physical qudits is denoted n, while the number of logical qudits is k.The bare (dressed) distance is denoted D bare (D dr ).