Pauli stabilizer models of twisted quantum doubles

We construct a Pauli stabilizer model for every two-dimensional Abelian topological order that admits a gapped boundary. Our primary example is a Pauli stabilizer model on four-dimensional qudits that belongs to the double semion (DS) phase of matter. The DS stabilizer Hamiltonian is constructed by condensing an emergent boson in a $\mathbb{Z}_4$ toric code, where the condensation is implemented by making certain two-body measurements. We rigorously verify the topological order of the DS stabilizer model by identifying an explicit finite-depth quantum circuit (with ancillary qubits) that maps its ground state subspace to that of a DS string-net model. We show that the construction of the DS stabilizer Hamiltonian generalizes to all twisted quantum doubles (TQDs) with Abelian anyons. This yields a Pauli stabilizer code on composite-dimensional qudits for each such TQD, implying that the classification of topological Pauli stabilizer codes extends well beyond stacks of toric codes - in fact, exhausting all Abelian anyon theories that admit a gapped boundary. We also demonstrate that symmetry-protected topological phases of matter characterized by type I and type II cocycles can be modeled by Pauli stabilizer Hamiltonians by gauging certain 1-form symmetries of the TQD stabilizer models.

Quantum error correcting codes are essential for protecting quantum information from environmental noise and faulty operations [1][2][3]. One of the most prominent classes of quantum error correcting codes are the two-dimensional topological quantum codes [4], in which quantum information is encoded in degenerate eigenstates of topologically ordered systems [5][6][7][8][9]. Topological quantum codes offer protection from arbitrary local errors, due to the local indistinguishability of the degenerate eigenstates [10], and admit fault-tolerant operations with low-overhead, derived from the properties of the anyons in the underlying topological order [4,[11][12][13]. Despite the merits of two-dimensional topological quantum codes, it can be challenging to assess the error correcting properties of a general topologically ordered system.
A class of topological quantum codes with particularly transparent error correcting properties are the topological Pauli stabilizer codes [4,10,14]. These combine the robust error correction of topological quantum codes with the simple algebraic structures of Pauli stabilizer codes [15]. In general, topological Pauli stabilizer codes are defined by Pauli stabilizer models -i.e., Hamiltonians whose terms are mutually commuting local products of Pauli operators [4]. By convention, the ground state subspace of the Pauli stabilizer model defines the logical subspace of the error correcting code. Topological Pauli stabilizer codes are exemplified by the Z 2 toric code (TC) and the generalization to N -dimensional qudits referred to as the Z N TC (or Z N quantum double), both introduced in Ref. [4].
The characteristic properties of a topological Pauli stabi-arXiv:2112.11394v4 [quant-ph] 15 Dec 2022 lizer code can be understood in terms of the topological order of the underlying Pauli stabilizer model. The Z N TC, for example, is based on the topological order of a Z N gauge theory [16], and the properties of the gauge charges and fluxes are central to fault-tolerant quantum computations using the Z N TC [10,17,18]. This motivates developing a classification of the topological phases of matter captured by Pauli stabilizer models. Such a classification informs us about the universal properties that can be exploited for quantum computation using a topological Pauli stabilizer code. This leads to the natural question: what topological phases of matter can be described by Pauli stabilizer models?
This question was partially addressed in the works of Refs. [19][20][21]. Under a technical assumption [22], it was shown in Refs. [19] and [20] that all (translation invariant) Pauli stabilizer models on qubits belong to the same phase as the Z 2 TC or decoupled copies of the Z 2 TC. Later, this was generalized to prime-dimensional qudits and made rigorous by Ref. [21]. Ref. [21] proved that, for any prime p, every (translation invariant) Pauli stabilizer model on p-dimensional qudits [23] has the same topological order as some number of decoupled copies of the Z p TC [24]. This completes the classification of (translation invariant) Pauli stabilizer models on prime-dimensional qudits, showing that the universal properties of Pauli stabilizer models on prime-dimensional qudits are fully captured by TCs. However, the works of Refs. [19][20][21] leave the classification of Pauli stabilizer models on composite-dimensional qudits (i.e. products of primes) unaddressed. One can ask: are there Pauli stabilizer models on composite-dimensional qudits that capture more exotic topological phases of matter -beyond that of decoupled layers of the Z N TC?
In this work, we answer this question in the affirmative. We construct Pauli stabilizer models defined on compositedimensional qudits that realize all Abelian anyon theories that admit a gapped boundary. These models define topological Pauli stabilizer codes whose underlying topological order is beyond that of the Z N TC or decoupled copies of the Z N TC. This represents a substantial step towards a full classification of topological Pauli stabilizer codes, which is evidently much richer than suggested by the classification for primedimensional qudits. Furthermore, based on the expectation that commuting projector Hamiltonians in two spatial dimensions must have gapped boundaries [25,26], we expect that the Pauli stabilizer models presented here yield a complete classification of topological Pauli stabilizer codes.
We emphasize that the Pauli stabilizer models introduced in this work can be distinguished from copies of the Z N TC by the properties of their anyonic excitations. For example, we define a Pauli stabilizer model that belongs to the double semion (DS) topological phase of matter. Similar to the Z 2 TC, the DS topological order has four types of anyons. However, in marked contrast, the DS topological order features semionic excitations -which produce a statistical phase of i upon interchange. As there is no such excitation for a Z 2 TC, the DS stabilizer model belongs to a distinct phase of matter. Instead, it can be interpreted as a twisted gauge theory, more formally known as a twisted quantum double (TQD), where the 'twist' manifests in the exotic exchange statistics of the gauge fluxes [27][28][29].
More generally, the Pauli stabilizer models presented here capture all TQDs with Abelian anyons, which according to Refs. [30] and [31], account for every Abelian topological order with gapped boundaries [32]. Therefore, we lay the groundwork for using universal properties beyond those of the TC for quantum computing using topological Pauli stabilizer codes. We also define Pauli stabilizer models of symmetryprotected topological phases of matter. While the ground state subspaces of such models are nondegenerate on a closed manifold -and hence, cannot be used to encode quantum information -they are of interest as they evade the no-go theorems of Ref. [33].
The Pauli stabilizer models described in this work are constructed from copies of the Z N TC by condensing (i.e., proliferating) certain bosonic anyons [34][35][36]. This exploits the fact that all TQDs with Abelian anyons can be obtained from decoupled copies of the Z N TC through condensation, as pointed out in Refs. [34] and [37]. We show further that the condensation can be implemented within the stabilizer formalism using few-body Pauli measurements. As a result, condensation maps Pauli stabilizer models to Pauli stabilizer models. We note that this builds off of earlier works, in which boson condensation was used to construct tensor network representations of TQD ground states [37,38] and symmetryprotected topological states [39].
It is worth noting that exactly solvable models of TQDs have been introduced in previous works. In particular, commuting projector Hamiltonians of TQDs were first described in Refs. [28,29,40]. These models were then leveraged to build non-Pauli stabilizer models of TQDs in Refs. [41,42]. While these models indeed capture the characteristic properties of TQDs, their potential for fault-tolerant quantum computation is relatively opaque (although, see Ref. [43]). The TQD models presented here have the key feature that they can be fully understood within the familiar Pauli stabilizer formalism. This opens up the possibility for novel applications of TQDs to quantum computation and provides algebraically simple models for simulating topological phases of matter.
The outline of this paper is as follows. In Section II, we introduce the paradigmatic example of our construction -the DS stabilizer model. We then derive the DS stabilizer Hamiltonian by condensing a boson in the Z 4 TC, and establish its topological order by mapping its ground state subspace to that of the DS string-net model using a finite-depth quantum circuit (with ancilla). Next, to generalize the DS stabilizer model to Pauli stabilizer models of TQDs, we provide a primer on Abelian anyon theories in Section III. Subsequently, in Section IV, we construct Pauli stabilizer models of TQDs, starting with a review of the characteristic data of TQDs with Abelian anyons. The associated Pauli stabilizer models are constructed by condensing anyons in decoupled layers of TCs. In Section V, we derive Pauli stabilizer models of SPT phases by condensing gauge charges in the Pauli stabilizer models of TQDs. Appendix A gives further details on the mapping of the DS stabilizer Hamiltonian to the DS string-net model, and Appendices B and C provide examples of TQDs and a general formula for the group structure of the anyons in TQDs.

II. DOUBLE SEMION STABILIZER MODEL
Before describing the general construction of topological Pauli stabilizer models of twisted quantum doubles with Abelian anyons, we provide a concrete example. Our example is a Pauli stabilizer model defined on four-dimensional qudits that belongs to the double semion (DS) phase of matter. In Section II A, we define the Pauli stabilizer Hamiltonian and show explicitly that its anyonic excitations exhibit the characteristic properties of the anyons in the DS topological order. In the subsequent section, Section II B, we construct the DS stabilizer model from a Z 4 toric code (TC) by condensing certain excitations. We then conclude the section with Section II C, where we use a finite-depth quantum circuit to map the ground state subspace of the DS stabilizer Hamiltonian to that of the DS string-net model of Refs. [44] and [45].
We start by recalling the characteristic properties of the anyons in the DS phase, which we label by {1, s,s, ss}. Similar to the anyons in the Z 2 TC phase, these form a Z 2 × Z 2 group under fusion, with the fusion rules given by: The exchange statistics and braiding of the anyons are determined by the function θ : {1, s,s, ss} → U (1), defined by: Physically, θ(a) is the phase obtained upon exchanging two identical a anyons. Eq. (2) tells us that s is a semion,s is an anti-semion, and ss is a boson, i.e., the wave function incurs a phase of i, −i, and 1, respectively. The statistics of the anyons are an invariant of topological order -they cannot be changed by any local perturbations of the Hamiltonian that preserve the energy gap. Therefore, the statistics in Eq. (2) distinguish the DS topological order from the Z 2 TC phase. The braiding relations of the anyons can then be deduced from the exchange statistics θ using the expression in Eq. (68), introduced later in the text. Here, we simply state that the semions and anti-semions have the following braiding relations: where B θ (a, a ) is the phase accrued from a full braid of the anyons a and a . We show that the anyonic excitations of the DS stabilizer model, described in the next section, obey Eqs. (1) and (2), which indicates that it belongs to the DS phase.
A. Definition of the Pauli stabilizer model Here, we describe a Pauli stabilizer Hamiltonian belonging to the DS phase. The Hamiltonian is defined on a square lattice with a four-dimensional qudit at each edge, as shown in where we have labeled the computational basis states by j ∈ Z 4 . These operators are sometimes referred to as the shift and clock operators, respectively. They satisfy the relations: and for any pair of edges e and e , they obey the commutation relations: Z e X e = iX e Z e e = e X e Z e e = e .
The DS stabilizer Hamiltonian then takes the form: where the sums are over all vertices v, plaquettes p, and edges e, respectively. We find it convenient to represent the Hamiltonian terms graphically, as shown below: Note that the definition of C e depends on whether the edge e is vertical or horizontal. Using the commutation relations in Eq. (6), it can be checked that the terms are mutually commuting.
Given that the Hamiltonian terms commute with one another, they define the stabilizer group S DS : where the angled bracket notation denotes that S DS is generated by the terms of H DS . By construction, the logical subspace of the stabilizer group S DS coincides with the ground state subspace of H DS . Explicitly, the logical subspace H L is the mutual +1 eigenspace of the stabilizers: On a torus, the DS Hamiltonian has a four-fold ground state degeneracy. This can be seen by counting the number of independent constraints imposed on the ground state subspace by the stabilizers. Letting N v denote the number of vertices, there are N v vertex terms, N v plaquette terms, and 2N v edge terms. These constraints are not entirely independent, however. In particular, the vertex term squares to a product of edge terms and plaquette terms. Therefore, although A v has order four, it only contributes N v order two constraints. There are also two global relations among the vertex terms and plaquette terms: Consequently, there are only N v − 1 independent vertex terms and N v − 1 independent plaquette terms. All together, we find that there are 4N v − 2 independent order two constraints. Given that there are two four-dimensional qudits per vertex, yielding a total Hilbert space of dimension 4 2Nv , the ground state subspace is four-dimensional: We now characterize the topological order of H DS by considering its anyonic excitations. The anyonic excitations of H DS are created by string operators such as those depicted in Fig. 2. There are three types of string operators: W s γ , Ws γ , and W ss γ . The first two are defined along an oriented pathγ in the dual lattice, while the third is defined along an un-oriented path γ in the direct lattice. To make the string operators explicit, we decompose them into products of short string operators: where the dashed red lines denote the orientation of the path on the dual lattice. Note that the expressions for W s γ and Ws γ in Eq. (13) are ambiguous up to a phase, since we have not specified an ordering of the short string operators. Nonetheless, any choice of ordering yields a string operator. Importantly, the operators in Eq. (13) commute with the Hamiltonian terms along the length of the path, and only fail to commute with vertex terms and plaquette terms at the endpoints. Notice also that the A v and B p terms of the Hamiltonian are small loops of the string operators, i.e.: whereγ v is the counter-clockwise oriented path through the FIG. 3. The exchange statistics of the anyon s can be computed using the formula in Eq. (18).γ1,γ2, andγ3 are oriented paths in the dual lattice incident upon the same plaquette p. The string operators (W s γ 2 ) † and W s γ 3 fail to commute, giving the statistics θ(s) = i. edges connected to v and γ p is the path formed by the edges bordering p.
We have suggestively labeled the string operators by the anyons s,s, and ss of the DS phase. Indeed, the string operators create anyonic excitations of the DS phase, as verified below. The first property to check is the fusion rules of the excitations. The fusion rules are obtained by multiplying string operators that share the same endpoints. We find that fusing two s string operators along the same pathγ in the dual lattice gives [46]: Up to stabilizers alongγ, this is equivalent to [47]: Therefore, the string operators labeled by s multiply to local operators at the endpoints. Since the trivial anyon 1 represents local excitations created by local operators, this means that s×s = 1. Similarly, it can be checked that the string operators Ws γ and W ss γ satisfy the fusion rules of the anyons in the DS phase up to stabilizers and local operators at the endpoints. Thus, the anyonic excitations of H DS satisfy the fusion rules in Eq. (1).
The next property to check is that the anyonic excitations created by the string operators in Eq. (13) have the same statistics as the anyons of the DS topological order. This is accomplished by following the prescription described in Refs. [48][49][50]. Letγ 1 ,γ 2 , andγ 3 be paths sharing a common endpoint p and ordered counter-clockwise around p, as in Fig. 3. Then, the exchange statistics θ(a) of the anyon a is computed by the expression: The computation of θ(s) is pictured in Fig. 3. We find that θ(s), θ(s), and θ(ss) are: FIG. 4. The logical operators of the DS stabilizer code on a torus can be represented by string operators along non-contractible loops such as α (blue) and β (red).
which match the statistics of the DS anyons in Eq. (2). Therefore, the anyonic excitations of H DS have the same fusion and statistics as the anyons in the DS topological order. This implies that H DS indeed belongs to the DS phase. The four-fold degeneracy of H DS on a torus can now be understood in terms of the DS anyons. Let α and β be generators of the non-contractible cycles on the torus, pictured in Fig. 4. Then, the long string operators W s α , W s β , Ws α and Ws β form Pauli X and Pauli Z on the logical subspace H L . To make this explicit, we define: In the logical subspace, the long string operators satisfy: where ∼ emphasizes that the relations hold only in the logical subspace. The only nontrivial commutation relations between the long string operators are: Thus, they have the same relations as Pauli X and Pauli Z operators on a pair of qubits. The relations in Eq. (21) follow from the fusion rules of s ands, while the commutation relations in Eq. (22) follow from the braiding relations in Eq. (3). Indeed, in this model, the braiding relations of s ands can be computed from the commutation relations of string operators that intersect at a single point (see Ref. [50] for a more general statement).

B. Construction of the Pauli stabilizer model
The DS stabilizer model can be derived from a Z 4 TC by condensing a certain bosonic anyon. To demonstrate this, we first describe the construction of the DS stabilizer model at the level of the anyons. We then implement the construction at the lattice level to arrive at the DS stabilizer model.

Anyon-level construction
Recall that the anyons of the Z 4 TC form a Z 4 × Z 4 group, generated by the anyons e and m. The sixteen anyons of the Z 4 TC are shown in table below: The statistics of an arbitrary anyon e p m q (with p, q ∈ Z 4 ) is given by: Note that the e and m anyons are often interpreted as the gauge charge and gauge flux, respectively, of a Z 4 gauge theory. According to Eq. (23), e and m are bosons and have the Aharanov-Bohm phase i, as expected: The construction of the DS stabilizer model from the Z 4 TC is motivated by the fact that the Z 4 TC has both semionic excitations and anti-semionic excitations. In particular, em is a semion and em 3 is an anti-semion. Further, their fusion product e 2 is a boson, similar to the fusion of the semion s and the anti-semions of the DS anyon theory. The only obstacle to identifying em and em 3 with s ands, respectively, is that em and em 3 have order four under fusion, while s ands have order two. Indeed, both em and em 3 square to the boson e 2 m 2 : The resolution is to condense the boson e 2 m 2 . As described below, condensing e 2 m 2 ensures that em and em 3 represent anyons with the same fusion rules as s ands. Moreover, the condensation of e 2 m 2 leaves us with precisely the anyon theory of the DS phase. The condensation of e 2 m 2 has two effects. The first is that anyons that braid nontrivially with e 2 m 2 become confined. As demonstrated later in this section, after condensation, the string operators of the confined anyons create excitations along the length of the string. Therefore, the string operators create nonlocal excitations, and the confined anyons do not correspond to anyonic excitations in the condensed theory. We call the remaining anyons, i.e., those that braid trivially with e 2 m 2 , deconfined anyons.
The second effect of condensation is that deconfined anyons related by fusion with e 2 m 2 become identified. This means, in particular, that e 2 m 2 is identified with the trivial anyon 1.
To emphasize the identification of anyons after condensation, we use square brackets. Thus, in the condensed theory, we have: More generally, for any deconfined anyon a in the Z 4 TC, we write: Note that the confined anyons are inconsistent with the identification of e 2 m 2 with the trivial anyon, since braiding with e 2 m 2 differs from braiding with 1.
The effects of condensing e 2 m 2 are shown in Table I. The anyons that remain after condensation can be labeled by: Given Eq. (27), the anyons of the condensed theory satisfy the fusion rules: The exchange statistics above are well-defined, as can be checked using Eq. (68) in Section III and the fact that em and

Lattice-level construction
With this, we are prepared to construct the DS stabilizer model at a lattice level, starting with the lattice model for the Z 4 TC. The Z 4 TC is defined on a square lattice with a fourdimensional qudit at each edge (Fig. 1). The Hamiltonian is the following sum of vertex terms and plaquette terms: where A TC v and B TC p are given by: We define the stabilizer group S TC as the group generated by the vertex terms and plaquette terms of H TC : The anyonic excitations of the Z 4 TC are generated by products of the e and m string operators, illustrated in Fig. 5. Fig. 5 also shows a string operator for a pair of e 2 m 2 excitations. The e 2 m 2 string operators can be generated by the short string operators: Note that there is ambiguity in how the e 2 and m 2 excitations are bound together to form e 2 m 2 . In Eq. (34), we have made an arbitrary choice for the short string operators. This choice ultimately determines the form of the vertex terms and the string operators of the DS stabilizer model. We point out that the short string operators W e 2 m 2 e chosen in Eq. (34) are precisely the C e terms from Eq. (8), and throughout the rest of this section, we use the notation C e to agree with Section II A. The next step in the construction of the DS stabilizer model is to condense the e 2 m 2 excitations. To this end, we define S C to be the group of e 2 m 2 string operators. This is generated by the set of short e 2 m 2 string operators {C e }: In particular, S C includes open string operators, which create, annihilate, and more generally proliferate the e 2 m 2 excitations. Therefore, after condensing e 2 m 2 , the states |ψ in the ground state subspace should satisfy: In other words, the stabilizer group for the condensed theory should contain S C as a subgroup. This is certainly the case for the DS stabilizer group in Eq. (9). Only certain elements of the stabilizer group S TC are compatible with the condition in Eq. (36). These are, in particular, the elements of S TC that commute with every element of S C . We define S C TC to be the subgroup formed by the stabilizers of S TC that commute with the elements of S C : To understand the elements of S C TC , it is useful to note that the stabilizer group S TC is generated by small loops of e and m string operators, corresponding to B TC p and A TC v , respectively. This means that the elements of S TC are, in general, products of closed string operators of the Z 4 TC anyons, such as loops of em string operators or e 2 string operators. Therefore, the elements of S C TC are loops of string operators that commute with the open e 2 m 2 string operators. These correspond to closed string operators of anyons that braid trivially with e 2 m 2 . As a result, S C TC can be generated by small loops of em string operators and small loops of e 2 string operators. These are exactly the A v and B p terms of the DS Hamiltonian in Eq. (8): S C TC is designed to exclude the closed string operators of confined anyons.
Finally, the stabilizer group of the condensed theory is generated by S C TC and S C , i.e., it is given by: The justification for taking S C TC , S C to be the stabilizer group of the condensed theory comes from considering the corresponding Hamiltonian, which is precisely the DS stabilizer Hamiltonian H DS . We show that H DS exhibits the two effects of condensation described below Eq. (25). Namely, the string operators of the confined Z 4 TC anyons create confined excitations, and the string operators of deconfined Z 4 TC anyons create identical anyonic excitations, if they differ by an e 2 m 2 string operator.
To give an example of a confined excitation in the DS stabilizer model, let us consider the e string operators of the Z 4 TC. An open e string operator W e γ , such as the one depicted in Fig. 5, fails to commute with the vertex terms A v of the DS stabilizer Hamiltonian at the endpoints of the path γ. It also fails to commute with the edge terms C e along the length of γ. Therefore, separating the vertex excitations at the endpoints of γ comes with an energetic penalty that grows linearly with the length of the string operator. This is indicative of confinement. On an infinite plane, it is not possible to create a single vertex excitation (i.e., change the eigenvalue of A v by i) without creating either an extensive number of edge excitations or a second vertex excitation. We note that confined excitations, such as those created by the e string operator in the DS stabilizer model, are only possible in Pauli stabilizer models on composite-dimensional qudits. For translationally invariant Pauli stabilizer models built from prime-dimensional qudits on an infinite plane, it is always possible to choose the local stabilizer terms of the Hamiltonian so that one can violate any given single stabilizer term using some (possibly nonlocal) product of Pauli operators [21]. In addition to confinement, we see that anyons created by string operators of the deconfined Z 4 TC anyons become identified, if the corresponding Z 4 TC anyons differ by e 2 m 2 . As a concrete example, the string operators for em and e 3 m 3 differ by an e 2 m 2 string operator. Since the e 2 m 2 string operators are stabilizers in the DS stabilizer model, the em and e 3 m 3 string operators create the same anyons when applied to a ground state of H DS . Therefore, the deconfined anyons em and e 3 m 3 have become identified, and H DS indeed describes the theory obtained after condensing e 2 m 2 in a Z 4 TC.
Interestingly, the construction of the DS stabilizer code described above can be understood in terms of the effects of Pauli measurements on stabilizer groups. Recall that, given a stabilizer group S, the measurement of a product of Pauli operators P results in a modified stabilizer group. The modified stabilizer group is generated by ±P (depending on the measurement outcome [51]) and the elements of S that commute with P . Starting with the stabilizer group S TC , measurements of the short string operators {C e } produce the stabilizer group: By post-selecting for +1 measurement outcomes or by error correction, we obtain the stabilizer group S DS . Hence, the condensation of e 2 m 2 can be implemented by simply measuring the set of operators {C e }. This yields an efficient construction of the DS stabilizer code from a Z 4 TC, requiring only two-body Pauli measurements.

C. Relation to the string-net model
We now identify a finite-depth quantum circuit U (with ancillary degrees of freedom) that maps the ground state subspace of the DS stabilizer Hamiltonian H DS to the ground state subspace of the DS string-net model H s-n DS . This establishes that H DS belongs to the same topological phase of matter as H s-n DS , by the arguments of Ref. [52]. Our strategy is to map H DS to H s-n DS using the following operations: • adding and removing ancillary degrees of freedom, • conjugating H DS by layers of U, • making ground space-preserving changes to the Hamiltonian.
This transformation shows that the ground state subspace of the Hamiltonian UH DS U † is equivalent to that of H s-n DS (up to adding ancillary degrees of freedom to H DS and H s-n DS ). We refer to Appendix A for a more explicit mapping of the ground states using notation from simplicial cohomology.
Let us start by briefly recalling the DS string-net model of Refs. [44] and [45]. The DS string-net model is defined on a hexagonal lattice with a qubit at each edge (see Fig. 6d). The Hamiltonian H s-n DS is a sum of two types of terms and is given by: The vertex terms A s-n v and plaquette terms B s-n p are represented as: where the Pauli X and Pauli Z operators in Eq. (42) are defined on qubits and S is the phase gate with diagonal elements diag(1, i). H s-n DS is not a stabilizer model, since the B s-n p terms only commute with one another up to products of A s-n v terms. Nonetheless, the ground state subspace is the mutual +1 eigenspace of the Hamiltonian terms.
To relate the DS stabilizer Hamiltonian to H s-n DS , the first step is to map H DS to a DS stabilizer model on a triangular lattice. The motivation is that the triangular lattice is dual to the hexagonal lattice, on which the DS string-net model is defined. This mapping is accomplished by introducing an ancillary four-dimensional qudit to the center of each square plaquette (see Fig. 6a). For convenience, we rotate the lattice by 45 • , so that the ancillary qudit lies on the horizontal edge of a regular triangular lattice, as shown in Fig. 6b. We label the vertices, edges, and faces of the triangular lattice by v, e, and f , respectively. We take the qudits on the horizontal edges e to be stabilized by the group X 2 e , Z 2 e . The Hamiltonian on the system with ancillary qudits is then: where the sums are over the set of horizontal edges E h . We now couple H DS to the ancillary qudits to create a DS stabilizer Hamiltonian on the triangular lattice. We do so by conjugating H DS by a finite-depth quantum circuit U CX composed of control-X gates. We denote a control-X gate with control qudit at edge e and target qudit at edge e by CX ee . The gate CX ee is defined by the mapping of operators: To specify the finite-depth quantum circuit U CX , we label the vertices of the upward pointing triangles according to: .
With this, U CX is the product of control-X gates: where the product is over the set of all upwards pointing triangles F up , and the control and target qudits are specified by pairs of vertices according to the labeling in Eq. (45). U CX forms the first two layers of U [53]. We construct the DS stabilizer Hamiltonian H I DS on a triangular lattice by conjugating the Hamiltonian in Eq. (43) by U CX . Here, we have started to label the intermediate Hamiltonians with Roman numerals. Up to redefining the generators of the stabilizer group (which preserves the ground state subspace), conjugation by U CX yields the Hamiltonian: The terms A I v , B I f , and C I e are given graphically as: Notice that the X 2 e terms in Eq. (43) have been absorbed into the definition of C I e , and the face term B I f associated to the upwards pointing triangle comes from conjugating the Z 2 e horizontal edge term. We point out that the construction of the DS stabilizer model in Section II B carries over to the triangular lattice and produces precisely the Hamiltonian in Eq. (47).
Since the DS string-net model is defined on a system of qubits, the next step is to map each four-dimensional qudit of the DS stabilizer Hamiltonian to a pair of qubits. This can be accomplished by introducing a pair of qubits to each edge. To make this explicit, we label a pair of qubits at an edge e by A and B (see Fig 6c) and denote the Pauli X and Pauli Z operators by X A e , Z A e , X B e , Z B e . The system on four-dimensional qudits can then be mapped to qubits by conjugating with a finite-depth quantum circuit U 2,4 . The unitary operator U 2,4 is defined by the mapping of operators [54]: The operators on the left-hand sides of Eq. (49) are defined on the four-dimensional qudit, while the operators on the righthand sides act on the pair of qubits. Here, CX AB e is the control-X B gate with the A qubit as the control qubit. The mapping of operators in Eq. (49) maps the Pauli stabilizer model H I DS to a non-Pauli stabilizer model H II DS . Note that the qudit degrees of freedom become decoupled from the system, and may be disregarded [55].
For simplicity, let us focus on the effects of the transformation in Eq. (49) on the C I e terms. In terms of the operator algebra on the A and B qubits, C II e = U 2,4 C I e U † 2,4 is: The form of C II e is noteworthy, because, as argued below, the B qubits can be disentangled from the ground states by applying control-Z gates. That is, the control-Z gates are designed to map each C II e in Eq. (50) to a single site X B e operator. The ground state subspace of the DS string-net model can now be prepared from the ground state subspace of H II DS by applying a finite-depth quantum circuit U CZ composed of control-Z gates. U CZ consists of two layers and takes the form: The first layer U AB is needed to decouple the B qubits, leaving us with one qubit per edge. The second layer U AA ensures a12 a13 a23 CZ12,23|a12a13a23 S † 12 S13S † 23 |a12a13a23 0 0 0 1 1 that the ground state wave functions have the correct amplitudes, i.e., the amplitudes match those of the ground states of H s-n DS defined on the hexagonal dual lattice. Let us define the layers of U CZ more carefully, starting with U AB . U AB is the product of control-Z gates: The product above is over all faces 123 , and CZ AB 12,23 is the control-Z gate between the A site on the edge 12 and the B site on the edge 23 with the ordering of vertices: Given an arbitrary face 123 , conjugation by U AB maps the Pauli X operators X A 12 and X B 23 according to: where X B vv denotes the B site Pauli X operator on the edge vv . Importantly, this maps each edge operator C II e to a single site operator X B e : The next layer, U AA , is a product of control-Z gates between the A sites. We let CZ AA 12,23 denote the control-Z operator between the A site on edge 12 and the A site on edge 23 using the ordering of vertices in Eq. (53). U AA is then defined to be the finite-depth quantum circuit: where the product is over the set of all upward pointing triangles F up . The Pauli X operators X A 12 and X A 23 defined on the edges of the upwards pointing triangle 123 are mapped as: Notice that conjugation by U AA does not affect the edge terms in Eq. (55).
After conjugating H II DS by U CZ , the ground states are in a product state on the B sites. This follows from the fact that the ground states are +1 eigenstates of the U CZ C II e U † CZ = X B e terms. This means that, without affecting the ground state subspace, we can freely replace all B site operators with the identity. In what follows, we suppress the the A site label and consider systems with a single qubit at each edge. Composing the transformations in Eqs. (46), (49), (52), and (56), we find that H DS is mapped by the finite-depth quantum circuit: to the Hamiltonian: where we have suppressed the ancillary degrees of freedom.
The vertex terms A III v and face terms B III f are pictured below: which holds in the B III f = 1 subspace, as shown in Table II [56]. Note that after replacing all of the control-Z gates using the identities above, the Hamiltonian terms are no longer mutually commuting. Nevertheless, they are unfrustrated and the ground state subspace is the mutual +1 eigenspace of the Hamiltonian terms.
Using these transformations, H III DS can be mapped to: where A IV v and B IV f are graphically represented by: When viewed as a model on the hexagonal dual lattice, H IV DS is precisely the DS string-net model in Eq. (41). By construction, H IV DS has the same ground state subspace as UH DS U † , up to adding ancillary degrees of freedom to H DS and H IV DS . Therefore, we have found a finite-depth quantum circuit with ancillary degrees of freedom that maps the ground state subspace of H DS to the ground state subspace of the DS string-net model.

III. PRIMER ON ABELIAN ANYON THEORIES
In this section, we review the characteristic data of Abelian anyon theories with the aim of generalizing the DS stabilizer model to a wider class of Abelian topological orders. We also introduce notation for boson condensation, as it is essential to the construction of the topological stabilizer models defined in the next section. Further details on Abelian anyon theories can be found in Ref. [57], and we refer to Ref. [58] for background on boson condensation.
An Abelian anyon theory A is specified by a pair (A, θ), where A is a finite Abelian group and θ is a function from A to U (1): that satisfies certain constraints. Intuitively, the group A corresponds to the Abelian group formed by the anyons under fusion, and θ encodes the statistics and braiding of the anyons. Given this interpretation, we refer to the elements of A as anyons. The DS anyon theory, for example, is specified by A = Z 2 × Z 2 with θ given by: where the elements of Z 2 × Z 2 have been labeled by the anyons {1, s,s, ss}. We call an anyon b a bosonic anyon if it has trivial exchange statistics: θ(b) = 1.
Based on physical arguments, only certain choices of the function θ define a consistent anyon theory. The first requirement is that θ satisfies: for any anyon a ∈ A and integer n. In words, this means that exchanging two a n anyons is equivalent to exchanging n 2 pairs of a anyons. This is shown graphically in Fig. 8a. The second constraint on θ is most naturally stated in terms of the function: defined by [25]: Physically, B θ (a, a ) captures the braiding relations between the anyons a and a in A. This can be seen by considering the exchange statistics of aa , as depicted in Fig. 8b. The second constraint on θ is then given by the fact that braiding a n around a is (tautologically) the same as braiding n copies of a around a . That is, the function B θ satisfies: for all a, a ∈ A and any integer n.
In more formal language, the function θ corresponds to a quadratic form. Specifically, one can define the function q : A → [0, 1) by the expression: According to Eq. (66), we have: and by Eq. (69), the function: is bilinear. Eqs. (71) and (72) are the defining properties of a quadratic form. Consequently, Abelian anyon theories are classified by quadratic forms [57]. Furthermore, in this work, we consider anyon theories with the property that every anyon braids nontrivially with at least one anyon in A. That is, for every a ∈ A, there exists an a ∈ A such that B θ (a, a ) = 1.
Such anyon theories are referred to as modular anyon theories and are classified, in particular, by nondegenerate quadratic forms. We note that, more generally, anyon theories are specified by the following three pieces of data (i) a set of anyon labels, (ii) the fusion rules of the anyons, and (iii) the socalled F -and R-symbols. In Ref. [59], it is shown that the pair (A, θ) determines the data (i)-(iii) for Abelian anyon theories.
For later purposes, let us introduce the concept of stacking Abelian anyon theories. Given the anyon theories A 1 = (A 1 , θ 1 ) and A 2 = (A 2 , θ 2 ), we can form an anyon theory A 1 A 2 , which we refer to as a stack of A 1 and A 2 . Physically, A 1 A 2 is the anyon theory of two decoupled topological orders that realize the anyon theories A 1 and A 2 , respectively. The stack of A 1 and A 2 is defined by the pair for any anyons a 1 ∈ A 1 and a 2 ∈ A 2 . Using this language, Refs. [19] and [21] showed that the anyon theories of translationally-invariant Pauli stabilizer models on primedimensional qudits are necessarily stacks of toric codes. An important tool for the construction of the Pauli stabilizer models in the next section is boson condensation. Thus, we next discuss the effects of boson condensation in Abelian anyon theories. At a physical level, the condensation of a boson b amounts to proliferating the b quasiparticles. As a consequence, pairs of bosons can be freely created and destroyed in the ground state. To make the effects on the anyon theory explicit, we let {b i } be a set of bosons with trivial mutual braiding relations. The two important effects of condensation are then as follows.
(i) The first effect is that the anyons that braid nontrivially with at least one b i in {b i } become confined. A single confined anyon is not a local excitation in the condensed theory, since the energetic cost of separating a pair of confined anyons grows linearly with the separation (see Section II B). The remaining anyons are referred to as deconfined anyons. A deconfined anyon a satisfies: It is important that the bosons {b i } have trivial mutual braiding relations, so that they do not confine one another upon condensation. (ii) The second effect is that anyons that are related by products of the condensed bosons become identified. More precisely, after condensation, the anyons can be organized into equivalence classes, where two anyons a and a are equivalent if they satisfy: for some set of integers {p i }. We denote the equivalence class containing the anyon a by [a]. In particular, the bosons b i are identified with the trivial anyon, i.e., We note that the statistics and braiding of the deconfined anyons are unaffected by the condensation. For any [a], [a ], we have: These are well-defined given the condition in Eq. (74). Lastly, we define the notion of a Lagrangian subgroup. A Lagrangian subgroup L is a subgroup of A composed of bosonic anyons with the following two properties: (i) The elements of L have trivial braiding relations with each other, i.e., for every b, b ∈ L, B θ (b, b ) = 1.
(ii) For every anyon a ∈ A − L, there exists a b ∈ L such that a braids nontrivially with b, i.e., B θ (a, b) = 1.
For example, the ss anyon of the DS anyon theory generates a Lagrangian subgroup. It trivially satisfies condition (i), and satisfies condition (ii) because it braids nontrivially with s and s. The existence of a Lagrangian subgroup signals the potential for a gapped boundary in a topologically ordered system. Specifically, a Hamiltonian whose excitations are described by the anyon theory A admits a gapped boundary if and only if A has a Lagrangian subgroup [30,60]. Moreover, Refs. [30] and [31] proved that every Abelian topological order that admits a gapped boundary can be described by a twisted quantum double (TQD). Here, we use TQD to refer to any model obtained by gauging the symmetries of a Hamiltonian belonging to a symmetry-protected topological (SPT) phase. The argument is that, if the bosons of the Lagrangian subgroup are condensed, then all of the anyons become confined (or identified with the trivial anyon). This is due to condition (ii) of a Lagrangian subgroup. Furthermore, the condensation can be implemented by gauging the 1-form symmetries associated with the closed string operators of the bosons in L. Since all of the anyons become confined, gauging the 1-form symmetry produces a model for an SPT phase protected by a 0-form symmetry. By subsequently gauging the 0-form symmetry of the SPT model, we recover the initial anyon theory. Thus, TQDs, including the Pauli stabilizer models constructed in the next section, must exhaust all Abelian anyon theories that can be realized in systems with gapped boundaries.

IV. TWISTED QUANTUM DOUBLE STABILIZER MODELS
We now construct Pauli stabilizer models corresponding to twisted quantum doubles (TQDs). The construction is a generalization of the construction of the DS stabilizer model in Section II. In particular, we start with a toric code (TC) (or stack of TCs) defined on composite-dimensional qudits and condense certain bosonic anyons so that the remaining deconfined excitations match those of a TQD with Abelian anyons. We refer to the resulting Pauli stabilizer models as TQD stabilizer models. Before describing the construction in detail, we recall some of the characteristic properties of TQDs. Note that TQD Hamiltonians were first presented in Refs. [28] and [61] based on the spacetime formulation in Ref. [27]. These models exhibit the characteristic properties of TQDs described below, although the details of the models are not essential to follow the discussion.

A. Anyon theories of twisted quantum doubles
In this work, "TQD" refers to any two-dimensional topologically ordered system that can be obtained by gauging the symmetry of a model with symmetry-protected topological (SPT) order. Given the classification of two-dimensional SPT phases [62], TQDs are characterized by a group G and a group cocycle ω belonging to H 3 [G, U (1)] -the same data that specifies the SPT phase. Moving forwards, we assume familiarity with group cohomology and refer to Refs. [62] and [63] for the necessary background. For certain choices of G and ω, the corresponding TQDs host non-Abelian anyonssee, for example, the models of Refs. [61,63,64]. Stabilizer models, however, are unable to model topological orders with non-Abelian anyons [65]. Therefore, we restrict our attention to TQDs whose anyons are Abelian. We refer to such TQDs as Abelian TQDs.
A TQD is an Abelian TQD if and only if G is a finite Abelian group and ω ∈ H 3 [G, U (1)] is cohomologous to a product of type I and type II cocycles (defined below) [66] [67]. To make the form for ω explicit, we write G as a general finite Abelian group: and introduce the set of integers I: Here, n i is an element of Z Ni , and n ij belongs to Z Nij , where N ij denotes the greatest common divisor of N i and N j . Furthermore, n ij satisfies n ij = n ji and n ii = 2n i . The set of integers I then specifies the following group cocycle ω I ∈ H 3 [G, U (1)] of an Abelian TQD [63]: Here, g, h, k are elements of G, β is an arbitrary 2-cochain, and ω i and ω ij are so-called type I and type II cocycles. Explicitly, ω i and ω ij are given by [63]: where g i , h i , k i ∈ Z Ni are the i th components of g, h, k, and [· · · ] Ni denotes addition modulo N i . In words, Eq. (79) tells us that ω I is cohomologous to a product of type I and type II cocycles determined by the integers I. Abelian TQDs are, therefore, conveniently parameterized by an Abelian group G and the set of integers I in Eq. (78). We now describe the characteristic properties of the anyonic excitation of Abelian TQDs. We rely on the fact that Abelian TQDs are gauge theories, derived by gauging the symmetry of SPT models. As such, a subgroup of the anyons in an Abelian TQD can be labeled as gauge charges [16,68]. The gauge charges are uniquely determined by the gauging procedure and reproduce the group G = M i=1 Z Ni under fusion. We use c i to denote the unit of gauge charge that generates the subgroup Z Ni . The gauge charges are always bosons with trivial mutual braiding, i.e.: All other anyonic excitations of Abelian TQDs carry nonzero gauge flux, implying that they braid nontrivially with at least one of the gauge charges. We call an anyon ϕ i an elementary flux if it carries a single unit of gauge flux, i.e., its braiding with gauge charges satisfies: where N ij denotes the least common multiple of N i and N j .
Using the properties of θ in Eqs. (66) and (69), it can be checked that Θ i and Θ ij are well-defined -they do not change if the elementary fluxes are modified by attachment of gauge charges.
In what follows, it useful to express Θ i and Θ ij in terms of the group cohomological data that specifies an Abelian TQD. The results of Ref. [63] imply that Θ i and Θ ij can be expressed in terms of the parameters in I as: We see that a nonzero n i implies that the elementary flux ϕ i has nontrivial statistics, and a nonzero n ij implies that the elementary fluxes ϕ i and ϕ j have nontrivial mutual braiding. To provide intuition for the quantities in Eq. (85), let us consider the TQDs corresponding to the group G = Z 2 . The group cohomology H 3 [Z 2 , U (1)] is equivalent to Z 2 , and the nontrivial element is cohomologous to a type I cocycle. This implies that there are two distinct Abelian TQDs for G = Z 2 , labeled by n 1 = 0 and n 1 = 1. The TQD corresponding to n 1 = 0 is in the same phase as the Z 2 TC and has anyons {1, e, m, ψ}. The gauge charge is the boson e, and the elementary flux can be chosen to be ϕ 1 = m or ϕ 1 = ψ. In either case, Θ 1 is: The TQD corresponding to n 1 = 1, on the other hand, gives a model belonging to the DS phase with anyons {1, s,s, ss}.
Here, ss is the gauge charge and the elementary flux is either ϕ 1 = s or ϕ 1 =s. In this case, Θ 1 takes the value: In Appendix B, we give further illustrative examples of Abelian TQDs using the K-matrix formalism. An important relation between elementary fluxes and gauge charges can be derived from Eq. (85). To this end, we consider the product of N i elementary fluxes ϕ i . The resulting anyon ϕ Ni i carries vanishing gauge flux, as can be seen by braiding a charge c j around ϕ Ni i . Therefore, ϕ Ni i must be generated by gauge charges. To probe the gauge charges carried by ϕ Ni i , we consider braiding an elementary flux ϕ j around ϕ Ni i . There are two cases: Since θ satisfies Eq. (66), this reduces to: Finally, we employ Eq. (85) to find: Eqs. (83) and (90) tell us that ϕ Ni i carries 2n i copies of the gauge charge c i .
(ii) If j = i, we use the property of B θ in Eq. (69): This can then be written in terms of Θ ij using Eq. (84): Substituting the expression for Θ ij in Eq. (85) into the formula above, we have: where we have used that N ij N ij = N i N j . Thus, ϕ Ni i carries n ij copies of the gauge charge c j .
Combining (i) and (ii), we find that the gauge charges and elementary fluxes satisfy [69]: This relation is independent of the choice of ϕ i due to the fact that Θ i and Θ ij , used in the derivation, do not depend on the choice of elementary fluxes. We also point out that Eq. (94) implies that the group formed by the anyons of an Abelian TQD can differ from G × G, if the parameters n i and n ij are nonzero. We refer to Appendix C for an aside on the group structure of anyons in Abelian TQDs.
In summary, the anyons of an Abelian TQD corresponding to G and I form a group A TQD with the presentation: where In the next section, our goal is to construct a Pauli stabilizer model for each choice of G = M i=1 Z Ni and parameters I such that its anyons satisfy Eqs. (82), (83), and (85) and form the group in Eq. (95). This implies that the anyonic excitations of the Pauli stabilizer Hamiltonian form the same anyon theory A TQD as an Abelian TQD specified by G and I.

B. Construction of the Pauli stabilizer models
Having described the anyonic excitations of TQDs, we are ready to construct the TQD stabilizer models. The construction proceeds in two steps. In the first step, we identify a TC (or stack of TCs) on composite-dimensional qudits with the property that its anyonic excitations contain the anyons of a TQD as a quotient group. In the second step, we condense certain emergent bosons so that the deconfined excitations have the same properties as the anyons of a TQD. We describe the construction at the level of the anyons before providing explicit lattice models for the Pauli stabilizer models.

Anyon-level construction
To capture the statistics and braiding of anyons in an Abelian TQD corresponding to a group G = M i=1 Z Ni , we start with a stack of M TCs, where the i th TC is a Z N 2 i TC. We note that, in some cases, it is possible to use a smaller group than Z N 2 i . However, Z N 2 i is sufficient for the general construction of TQD stabilizer models, as the anyons in stacks of Z N 2 i TCs are able to reproduce the statistics and braiding needed for the characteristic Θ i and Θ ij of the Abelian TQD corresponding to G.
We label the e and m excitations of the i th TC by e i and m i , respectively. Recall that the e and m excitations are bosons with the braiding relations: Anyons in the stack of TCs can be written as products of e and m excitations of the form: M i=1 e pi i m qi i , for integers p i , q i ∈ Z N 2 i . The statistics and braiding of the anyons follow from the properties of the e and m excitations and are given by the general formulas: Our first objective is to identify anyons in this system that exhibit the same statistics and braiding as the anyons of an Abelian TQD corresponding to G = M i=1 Z Ni and I, where I specifies a product of type I and type II cocycles, as described in the previous section.
Let us first consider the excitations {e Ni i }. These excitations generate the group G under fusion and are bosons with trivial mutual braiding: Therefore, they satisfy the same conditions as in Eq. (82) for gauge charges of an Abelian TQD. We suggestively define: To describe excitations that mimic the properties of elementary fluxes, we set a convention for the ordering of the factors Z Ni of G. We first assume (without loss of generality) that each N i is a power of a prime. We then (implicitly) order the factors of G so that N i < N i+1 , for all i. For example, we order Given this ordering on the factors of G, we consider the excitations {φ i } for 1 ≤ i ≤ M , whereφ i is defined as: Here, the product is over all j with j > i. Note that due to our choice of ordering,

Nj
Ni n ij is an integer. Indeed, if N i does not divide N j , then the greatest common divisor of N i and N j is N ij = 1, implying that n ij must be zero, since it is an element of Z Nij . From the general formula in Eq. (98), the braiding of ϕ i and the excitationc i is: This agrees with the braiding of an elementary flux and a gauge charge, as in Eq. (83). Furthermore, the statistics and mutual braiding ofφ i tell us, for all i, j: which matches the values for Θ i and Θ ij in Eq. (85).
Thus far, we have found excitationsc i andφ i that exhibit the same statistics and braiding as the gauge charges and elementary fluxes of an Abelian TQD. Now, we need to ensure that they also have the same fusion rules. However, using the definition ofφ i , the product of N i copies ofφ i is equivalent to:φ This implies thatc i andφ i fail to satisfy the relations of gauge charges c i and elementary fluxes ϕ i in Eq. (94), rewritten here: This discrepancy can be resolved by condensing a particular set of bosonic anyons {b i }, with 1 ≤ i ≤ M . The bosons are chosen so that the anyons of the resulting theory obey the relations in Eq. (105). In particular, we construct an anyon theory C by condensing the collection of bosons {b i }, defined by: where the product is over j with j < i. The first thing to note is that the bosons have trivial mutual braiding relations: Therefore, they can be condensed without confining one another. Second, the bosons have trivial braiding relations with the anyonsc i andφ i : This means that the equivalence classes [c i ] and [φ i ] represent deconfined anyons in C. Here, we have used the notation [ · ] described in Section III to represent equivalence classes of excitations related by fusion with condensed bosons. The condensation of the bosons produces relations between [c i ] and [φ i ] analogous to the relations in Eq. (105). To see this, we express b i in terms ofφ i andc i . A straightforward calculation shows that b i is equivalent to: Therefore, after condensation, we have: Given that the equivalence relation is well-defined under fusion, we obtain: Eq ]} with the same statistics, braiding, and fusion rules as the gauge charges and elementary fluxes of an Abelian TQD corresponding to G and the parameters I. Hence, they generate an anyon theory that is equivalent to the anyon theory A TQD of the Abelian TQD. This means that the condensed theory C contains A TQD as a subtheory.
The last step of the construction is to argue that C For the sake of contradiction, suppose that there exists an anyon a that remains deconfined and nontrivial after condensing the sets of bosons {c i } and {b i }. Since a is deconfined after the condensation of thec i anyons, it must braid trivially with eachc i = e Ni i . This means that a takes the form: with p i ∈ Z N 2 i and s i ∈ Z Ni . After condensing the bosons in {c i }, each b i can be identified with a product of m i excitations: for r i ∈ Z Ni . The expression in Eq. (114) conflicts with the assumption that a is nontrivial after the condensation of the two sets of bosons. Indeed, after condensing the bosonsc i , we have: The observation above tells us that C is modular, i.e., for every anyon [a] in C there exists an anyon [a ] that braids nontrivially with [a]. This is because, if an anyon braided trivially with all other anyons in C, then it would be deconfined after condensing {[c i ]} [70]. Since C is modular, we can apply the following result from Ref. [71]: any modular anyon theory A, that contains a modular subtheory A 1 , factorizes as: for some modular theory A 2 . The product , defined in Section III, implies that the anyons in A 1 braid trivially with the anyons in A 2 . A TQD forms a modular subtheory of C. Therefore, C factorizes as: for some modular anyon theory B. The anyon theory B must be trivial, because, after condensing the [c i ] bosons, there are no deconfined anyons. If B was nontrivial, then the anyons in B would remain deconfined after condensing {[c i ]}. We conclude that C is equal to A TQD .

Lattice-level construction
Having described the construction at the level of the anyons, we now turn to an explicit lattice construction of the TQD stabilizer models. The stabilizer model for the TQD specified by G = M i=1 Z Ni and the set of integers I [Eq. (78)] is defined on a square lattice, where each edge hosts M qudits, and the i th qudit has dimension N 2 i . The Pauli X and Pauli Z operators for the i th qudit at the edge e are given by: When it is clear from context, we omit the edge label e.
The construction of the TQD stabilizer model begins with a stack of Z N 2 i TCs. The Hamiltonian for the decoupled Z N 2 i TCs is: with A TC v,i and B TC p,i represented pictorially by: The associated stabilizer group for the Hamiltonian H TC is generated by the full set of vertex terms and plaquette terms of the Hamiltonian H TC : The TQD stabilizer model is then derived by condensing the bosons b i , defined in Eq. (106) as: We take the short string operators for the bosons b i to be: whereẐ i is shorthand for the product of Pauli Z operators: To build longer string operators, one would need to account for the orientation of a path along the dual lattice, but for our purposes the operators C e,i suffice. The short string operators C e,i generate the stabilizer group S C : We now follow the same logic as used in Section II B to condense the set of bosons {b i }. Let us define S C TC to be the group of stabilizers in S TC that commute with the elements of S C : With some foresight, this is generated by the set of operators , defined graphically as: Here,Ž i is notation for the operator: The stabilizer group of the condensed theory is then: The corresponding TQD stabilizer Hamiltonian H TQD is defined as: The unit gauge charge c i and an elementary flux ϕ i can be created by the short string operators below: whereZ i is given by the product of Pauli Z operators: (134) These agree with the formulas forc i andφ i in Eqs. (100) and (101). Longer string operators, along oriented paths γ in the direct lattice andγ in the dual lattice, are given by: Lastly, as a consistency check, we compute the ground state degeneracy of H TQD on a torus. Similar to the calculation in Section II A, we count the number of independent constraints on the states in the ground state subspace, defined by: We expect the ground state degeneracy to have order |G| 2 , since the ground state degeneracy of a topological order on a torus is equivalent to the number of anyons in the theory [25].
To count the number of constraints on the ground state subspace, we let N v denote the number of vertices in the lattice. Then, for each layer i, there are N v vertex terms, N v plaquette terms, and 2N v edge terms. Each of these yields an order N i constraint. The vertex term A v,i gives an order N i constraint, since A Ni v,i is a product of plaquette and edge terms. On a torus, there are also global relations among the vertex stabilizers and plaquette stabilizers that need to be considered. In particular, for every i, the vertex terms and plaquette terms satisfy: Therefore, we only have N v −1 independent vertex stabilizers and plaquette stabilizers. This gives us a total of 4N v − 2 constraints with order N i . Since there are two N 2 i -dimensional qudits per vertex, the dimension of the Hilbert space H i in the layer i is: meaning that the dimension of the total Hilbert space H is: Given that there are 4N v −2 constraints with order N i for each layer i, the dimension of the ground state subspace is equal to: This is equivalent to |G| 2 , the order of the fusion group for the TQD.
Note that the calculation above is sufficient to show that the set of operators TC . This is because, in the absence of the global relations (i.e., on a simply connected manifold), the ground state is nondegenerate. Consequently, any element of S TC that commutes with the elements of S TQD , must already belong to S TQD (see Appendix C of Ref. [33]).

V. PAULI STABILIZER MODELS OF SPT PHASES
As established in Sections III and IV, TQDs are derived by gauging the symmetry of models belonging to symmetryprotected topological (SPT) phases. Specifically, TQDs corresponding to a group G and cocycle ω ∈ H 3 [G, U (1)], are obtained by gauging the G symmetry of the associated SPT models. Conversely, models of SPT phases can be constructed from TQDs by gauging 1-form symmetries [72]. In this section, we demonstrate that gauging certain 1-form symmetries of the TQD stabilizer models produces Pauli stabilizer models of SPT phases. To make the discussion concrete, we focus on a stabilizer model for a Z 2 SPT phase, obtained by gauging a 1-form symmetry of the DS stabilizer Hamiltonian. Similar arguments can be used to construct Pauli stabilizer models for SPT phase characterized by group cocycles of the form in Eq. (79).
Indeed, TQDs, and topologically ordered systems more generally, possess 1-form symmetries generated by anyon string operators along closed paths. Models belonging to the DS phase, for example, have a 1-form Z 2 ×Z 2 symmetry with generators given by loops of semion s and boson ss string operators. We point out that 1-form symmetries in a topologically ordered system need not correspond to anyonic excitations, however. The DS stabilizer Hamiltonian in Section II has a 1-form Z 4 × Z 2 symmetry generated by the vertex terms A v and plaquette terms B p . Only the Z 2 × Z 2 subgroup corresponds to anyon string operators.
Given the connection between 1-form symmetries and anyons, we can make a more precise statement about the construction of SPT models from TQDs. Namely, models for SPT phases arise from gauging the 1-form symmetries associated to the gauge charges of TQDs. After gauging the 1-form symmetry associated to the gauge charges, the system (i) gains a 0-form symmetry and (ii) has trivial topological order (in the absence of symmetries) -both of which are required for SPT phases. The theory obtained by gauging the symmetry has trivial topological order, since gauging 1-form symmetries has the same effect as condensing the corresponding anyons (elaborated upon for the DS stabilizer model below). The gauge charges of the TQD form a Lagrangian subgroup (Section III) -thus, condensing the gauge charges confines all of the anyons.
We now focus on the DS stabilizer model. For the DS stabilizer model, in particular, the gauge charge is the boson ss with Z 2 fusion rules. Therefore, we obtain a model for a Z 2 SPT phase by gauging the 1-form symmetry corresponding to ss. For convenience, we re-write the DS stabilizer Hamiltonian H DS : with A v , B p , and C e defined by: We also recall that the string operators for the anyon ss are generated by the short string operators: Explicitly, the 1-form symmetry corresponding to ss is generated by operators of the form: where λ denotes a loop in the direct lattice. At this point, one could gauge the 1-form symmetry following the usual minimal coupling prescription [45,73,74]. However, we find it illuminating to take a different approach. For Pauli stabilizer Hamiltonians, gauging a 1-form G symmetry corresponding to a boson b is equivalent to condensing b anyons bound to 0-form G charges. More concretely, for the DS stabilizer model, the 1-form symmetry associated to ss can be gauged by condensing ss anyons bound to 0-form Z 2 charges, as described below.
To make this explicit, we add a qubit to each vertex of the square lattice, as in Fig. 9. We use X v and Z v to denote the Pauli X and Pauli Z operators at the vertex v. In what follows, we graphically represent X v and Z v with red operators. We define a modified HamiltonianĤ DS with ancillary qubits to FIG. 9. The Pauli stabilizer model of the Z2 SPT phase is defined on a square lattice with a four-dimensional qudit (blue) at each edge and a qubit (red) at each vertex. be:Ĥ which corresponds to the stabilizer group: By construction, this Hamiltonian has an additional 0-form Z 2 symmetry generated by the product of Pauli X operators v X v . The Pauli Z operator Z v creates a 0-form Z 2 charge at the vertex v. In terms of the minimal coupling prescription, the qubits on the vertices are the gauge fields.
Next, we introduce the short string operators: which proliferate ss anyons bound to 0-form Z 2 charges. The claim is that gauging the 1-form symmetry is equivalent to condensing the excitations created by the short string operators D e . From the perspective of minimal coupling, the operators in Eq. (147) can be understood as the 1-form Gauss's law.
In accordance, products of D e along a closed path recover a 1-form symmetry operator.
To condense the ss anyons bound to Z 2 charges, we follow the logic described in Sections II B and IV B. We start by defining the stabilizer group S D generated by the set of short string operators {D e }: We then define the stabilizer groupŜ D DS formed by the elements ofŜ DS that commute with the operators in S D . Formally, this is the stabilizer group: The only generators ofŜ DS that fail to commute with the elements of S D are the vertex terms A v and X v . Their products of the form A v X v , however, commute with every element of S D . We claim that the stabilizer group in Eq. (149) is generated by the following elements: Note that the plaquette terms B p are not included as generators, since they can be generated by products of the edge terms. The products A v X v can be interpreted as coupling the vertex terms A v to the gauge fields. This ensures that the Hamiltonian terms are gauge invariant, i.e., commute with the D e terms. Finally, the stabilizer group of the condensed theory is: and the corresponding stabilizer Hamiltonian is: The Hamiltonian terms are pictured below: Since the condensation procedure does not affect the Z 2 symmetry ofĤ DS , the Hamiltonian H SPT has a 0-form Z 2 symmetry generated by v X v .
In what follows, we confirm that H SPT has no anyonic excitations by computing the ground state degeneracy of H SPT on a torus. We then use the methods of Ref. [75] to diagnose its SPT order. We find that the Hamiltonian in Eq. (152) describes the Z 2 SPT phase characterized by the nontrivial el- Similar to the previous calculations of the ground state degeneracy in Sections II A and IV B, we count the number of constraints in terms of the number of vertices N v . Given that A v X v squares to a product of C e and D e terms, the Hamiltonian terms give 5N v order two constraints. Unlike the case of the TQD stabilizer models, there are no global relations amongst these constraints. Taking into account the qubit at each vertex and the four-dimensional qudit at each edge, the dimension of the total Hilbert space is 2 Nv 4 2Nv . After imposing the 5N v order two constraints, we find that the ground state subspace H L has dimension: dim(H L ) = 2 Nv 4 2Nv /2 5Nv = 1.
FIG. 10. The action ofPR(1) is shown above. It is equivalent to PR(1) for sites that are both inside the region R (shaded grey) and away from the boundary of R.
Therefore, H SPT has a unique ground state on a torus (in fact on any manifold without boundary) and does not admit anyonic excitations. We now deduce the SPT phase described by H SPT by following the methods of Ref. [75] -i.e., by considering an effective boundary symmetry action, we compute a group cocycle ω ∈ H 3 [Z 2 , U (1)]. For simplicity, we assume H SPT is defined on an infinite plane and use R to denote the lower half plane. We then study the effective boundary symmetry action along the boundary of R. To this end, we define P R (g) to be the symmetry action restricted to R, i.e.: where g is a {0, 1}-valued element of Z 2 . We also definẽ P R (g) to be:P with the product over all A v X v whose support is contained within R.P R (1) is portrayed in Fig 10. Notice that away from the boundary of R, the action of P R (1) matches that ofP R (1). As argued in Ref. [33], their difference gives the effective boundary symmetry action P(g): P(1) acts as the identity away from the boundary of R and near the boundary of R can be graphically represented as: We emphasize that P(g) is a tensor product of Pauli operators. Conventional wisdom says that, if the effective boundary symmetry action is a tensor product, then the SPT model belongs to the trivial SPT phase [76]. However, the more accurate statement is that the SPT model belongs to the trivial SPT phase if the effective boundary symmetry action is a tensor product of linear representations of the symmetry. P(g) is not a tensor product of linear representations of Z 2 . P(g) only satisfies the Z 2 group laws in the boundary Hilbert space H B , given by the set of states: Here, we have used supp(S) to denote the support of the stabilizer S. (See Refs. [33,73,75] for more details on the boundary Hilbert space of SPT models.) Indeed the effective boundary symmetry action squares to: which is a product of stabilizers whose support is contained within R. The next step in diagnosing the SPT order is to truncate the effective boundary symmetry action. We truncate the effective boundary symmetry action in Eq. (158) to an interval along the boundary of R with endpoints v L and v R . According to the arguments in Refs. [33] and [75] the ambiguities at the endpoints of the truncation do not affect the characterization of the SPT phase. We choose the truncation of P(1) to take the form: Furthermore, we take P (0) to be the identity. In the boundary Hilbert space, P (g) only satisfies the Z 2 group laws up to operators at the endpoints: where ∼ denotes that the relation holds in H B . We define Ω(g, h), for g, h ∈ Z 2 to be the right hand side of Eq. (162) if g = h = 1 and the identity otherwise. Then the truncated effective boundary symmetry action satisfies: Ω(g, h) can be decomposed as: where Ω v L (g, h) and Ω v R (g, h) are localized near v L and v R , respectively. We take Ω v L (g, h) to be: The last step is to consider the associativity of the truncated effective symmetry actions for g, h, k ∈ Z 2 :

P (g)[P (h)P (k)] = [P (g)P (h)]P (k).
(166) Substituting Eq. (163) into Eq. (166), we obtain: Given the decomposition of Ω(g, h) in Eq. (164), the condition above is satisfied up to a phase at each endpoint, i.e., at the endpoint v L , we have: Here, ω is the U (1)-valued group cocycle in H 3 [Z 2 , U (1)] that characterizes the SPT phase. From Eq. (168), we see that ω is given by: This represents the nontrivial element of H 3 [Z 2 , U (1)], implying that H SPT describes the nontrivial Z 2 SPT phase. We would also like to point out that the calculation of ω above, is analogous to the calculation of the F -symbol [50,77] for the semion s in the DS stabilizer model. This is because the effective boundary symmetry action only differs from the s string operator by X v operators, which do not affect the calculation.
We have now shown that the nontrivial two-dimensional Z 2 SPT phase can be modeled by the Pauli stabilizer Hamiltonian H SPT . We emphasize that this does not conflict with the results of Ref. [33], since the argument in Ref. [33] assumes that the Z 2 SPT model is defined on qubits, as opposed to fourdimensional qudits. Furthermore, the construction of Pauli stabilizer models for SPT phases presented here is restricted to SPT phases characterized by group cocycles that are products of type I and type II cocycles. We are unable to construct Pauli stabilizer models for SPT phases characterized by type III cocycles, in agreement with Ref. [33]. This is also consistent with the fact that type III cocycles correspond to TQDs with non-Abelian anyons, which cannot be modeled by Pauli stabilizer Hamiltonians.
We note that the Z 2 SPT model belongs to a trivial SPT phase, when considered as an SPT model protected by a Z 4 symmetry. This is because the effective boundary symmetry action in Eq. (158) is a tensor product of linear representations of Z 4 . In other words, the SPT model can be trivialized by extending the Z 2 symmetry to a Z 4 symmetry [78,79]. More generally, an SPT model protected by a i Z Ni symmetry, which is characterized by a product of type I and type II cocycles, can be trivialized by extending the symmetry to i Z N 2 i . This yields the commutative diagram shown in Fig. 11, which further motivates the construction of the i Z Ni TQD stabilizer model from a i Z N 2 i TC.
11. An SPT model protected by a i ZN i symmetry (and characterized by a product of type I and type II cocycles) is trivialized by extending the symmetry to i Z N 2 i . Subsequently gauging the 0-form symmetry produces a model in the same phase as a i Z N 2 i TC. One can then condense bosons according to the prescription in Section IV to obtain a i ZN i TQD stabilizer model and gauge the 1-form symmetries corresponding to the gauge charges to arrive at a stabilizer model belonging to the i ZN i SPT phase.

VI. DISCUSSION
We have constructed a Pauli stabilizer model for every twisted quantum double (TQD) with Abelian anyons. Our strategy was to condense bosonic anyons in decoupled toric codes (TCs) so that the remaining deconfined excitations are those of the Abelian TQD. As an example, we constructed a Pauli stabilizer model of the double semion (DS) topological order in Section II. We explicitly verified that the DS stabilizer model belongs to the same phase as the DS string-net model by finding a finite-depth quantum circuit (with ancillary degrees of freedom) that maps between the ground state subspaces. In addition, we described how Pauli stabilizer models of SPT phases can be obtained from the TQD stabilizer models. We made the construction explicit for the nontrivial Z 2 SPT phase of Ref. [28].
Our work builds on the classification of topological Pauli stabilizer codes initiated by Refs. [19][20][21]. We conjecture that our models give a complete classification of topological Pauli stabilizer codes up to finite-depth Clifford circuits with ancilla. For four-dimensional qudits, for example, this would imply that every Pauli stabilizer model is locally equivalent to decoupled copies of Z 4 TCs, Z 2 TCs, DS stabilizer codes, and six-semion stabilizer codes (defined in Appendix B). However, the proof of such a classification is challenging, as it involves working with polynomial rings over finite rings as opposed to polynomial rings over finite fields in the case of prime-dimensional qudits. More specifically, the current approach to classification requires a rigorous statement about the existence of a bosonic anyon. The argument in Ref. [49] uses Hasse-Minkowski theorem to prove that there is a gapped one-dimensional boundary, implying that there is a nontrivial Lagrangian subgroup for a topologically ordered system (see Section III). This does not directly apply to Z 4 and hence, new techniques need to be developed for this case. Furthermore, the current classification requires a proof that all local excitations of a Pauli stabilizer Hamiltonian are fully mobile in two dimensions. We expect that a general-ization of Hilbert's syzygy theorem is needed for systems of composite-dimensional qudits. We leave the details to forthcoming works.
Another important future objective is to understand the quantum error correcting properties of the TQD stabilizer codes. As a first step, it would be interesting to compute the optimal error thresholds using the statistical mechanical mappings of Refs. [10,[80][81][82]. In such mappings, a spin is associated with each local stabilizer generator and coupling strengths are determined by the so-called Nishimori conditions [83]. This yields a disordered statistical mechanical Hamiltonian, where the disorder realization depends on the configuration of Pauli errors on the code space. Subsequently, the phase diagram for the model can be studied using Monte Carlo methods, and the optimal error threshold of a topological Pauli stabilizer code can be read off from the critical point along the Nishimori line. One could also consider other methods for estimating the error thresholds, such as those of Refs. [84][85][86], which were used to compute the thresholds of Z p TCs with prime p. Beyond studying the error correcting properties, it would also be interesting to explore whether there are fault-tolerant operations that are more natural to implement on TQD stabilizer codes as compared to TCs.
The Pauli stabilizer models presented here may be of interest beyond their potential for quantum error correction. In particular, they are a substantial simplification from their stringnet counterparts and may be amenable to simulation on manybody qudit platforms. Moreover, as mentioned in Section II B, the ground state of the DS stabilizer Hamiltonian can be prepared by making two-body measurements of a Z 4 TC ground state. This gives an efficient construction of the DS ground state starting from the Z 4 TC ground state. One could make further measurements to obtain the ground state of the Z 2 SPT model of Section V. Anyon condensation in TCs also allows for the construction of symmetry-enriched topological phases, so similar methods as Section V could be used to create simulable models of these phases. It may be insightful to compare the condensation approach employed here to other means of preparing topological states, such as in Ref. [87].
Finally, while TQDs capture all Abelian anyon theories that admit gapped boundaries [30,31], there are of course, twodimensional Abelian anyon theories that possess obstructions to gapped boundaries (e.g., a nonzero chiral central charge). It is expected that such Abelian anyon theories cannot be described by Pauli stabilizer models in two spatial dimensions [25,26]. However, Abelian anyon theories without gapped boundaries have been realized on the boundary of commuting projector Hamiltonians in three spatial dimensions [88]. Refs. [89] and [90] have taken this one step further and constructed three-dimensional Pauli stabilizer models that host chiral Abelian anyon theories (and ungappable Abelian anyon theories, more generally) on their two-dimensional surface. Up to stacking with two-dimensional TQD stabilizer models and condensing bosons, these provide a three-dimensional Pauli stabilizer model for every two-dimensional Abelian anyon theory, wherein the anyon theory is realized on the boundary. We also note that it may be interesting to consider fault-tolerant quantum computation in these models along the lines of Refs. [91] and [92].
Abelian anyon theories without gapped boundaries have also been identified in the context of topological subsystem codes, such as the three-fermion subsystem code in Refs. [80,[92][93][94]. Roughly speaking, topological subsystem codes correspond to a parameter space of frustrated Hamiltonians with common conserved quantities, which are taken to be the stabilizers of the subsystem code [80]. In the honeycomb model of Ref. [25], for example, loops of fermionic string operator are preserved throughout the phase diagram and define the stabilizer group [95]. With this, one can assign an anyon theory to a topological subsystem code based on the conserved quantities. The basic idea is that the stabilizers of the subsystem code generate a 1-form symmetry associated to an Abelian anyon theory. Preliminary work suggests that we can leverage our TQD stabilizer codes to build a topological subsystem code for every two-dimensional Abelian topological order, regardless of whether the theory admits gapped boundaries. Such models may provide natural candidates for systems that host non-Abelian anyons, similar to the honeycomb model of Ref. [25]. In what follows, we consider two-dimensional triangulated manifolds M equipped with a branching structure, i.e., an assignment of an orientation to each edge of the lattice with the property that there are no cycles around any face (see Fig. 12a, for example). The branching structure determines an ordering of the vertices of a triangle, as shown in Fig. 12b. We refer to the vertices v, edges e, and faces f as 0-simplices, 1simplices, and 2-simplices, respectively. We label a p-simplex σ p by its vertices 1 . . . p + 1 , where the vertices are ordered according to the branching structure. For example, a face can be labeled as 123 with edges 23 , 13 , and 12 and vertices 1 , 2 , and 3 .
A p-chain on M is a formal sum of p-simplices with coefficients in an Abelian group A. An arbitrary p-chain s p takes the form: Explicitly, arbitrary 1-chains s 1 and 2-chains s 2 can be written in the form: The p-chains over A form a group denoted by C p [M, A]. For two-dimensional manifolds, the group C p [M, A] is taken to be the trivial group for p < 0 and p > 2.
The boundary operator ∂ is a linear map from the group of (p + 1)-chains to the group of p-chains: The action of the boundary operator on the arbitrary chains in Eq. (A2) is given by: α 123 ( 23 − 13 + 12 ).

(A4)
While p-chains and the boundary operator ∂ do not appear in the main text, they are essential to defining cochains and the coboundary operator, described below. A p-cochain a is a linear map from the group of p-chains to A: The Abelian group formed by p-cochains is denoted by C p [M, A]. Letting A = Z N , for some N , we use v and e to represent the cochains: Any 0-cochain a or 1-cochain c can then be expressed as a linear combination of the cochains v and e: for some α v , α e ∈ Z N . The coboundary operator δ is a linear map from the group of p-cochains to the group of (p + 1)-cochains: Given a p-cochain a, the coboundary of a is defined by: where s is an arbitrary (p + 1)-chain. For example, taking A to be Z 2 , the coboundary of v satisfies: This says that δv(e) is 1 for every edge connected to v. More generally, let b be the for some β v ∈ Z 2 . Then the coboundary of b is equal to: which, given Eq. (A10), evaluates to 1 on edges along closed paths in the dual lattice (see Fig. 13). This shows that the states |δb in Appendix A 2 correspond to configurations of loops in the dual lattice.
Lastly, we introduce the cup product ∪. The cup product maps a p-cochain and a q-cochain to a (p + q)-cochain: The cup product of the p-cochain a p and the q-cochain a q evaluated on arbitrary (p+q)-simplex σ p+q = 1 . . . p+q+1 is: Here, 1 . . . p + 1 is the p-simplex formed by the first p + 1 vertices of σ p+q and p + 1 . . . p + q + 1 is the q-simplex formed by the last q + 1 vertices of σ p+q .
As an example of the cup product, we consider the Z 2valued cochain e ∪e from Eq. (A19). e and e are 1-cochains, so the cup product is a 2-cochain. Evaluated on a face 123 , we have: e ∪ e( 123 ) = e ( 12 )e ( 23 ).
(A15) This is nonzero if and only if e = 12 and e = 23 . Another interesting example is the cup product v ∪ δv, which appears in Eq.
for any edge 12 . This is nonzero if and only if v = 1 and 12 is connected to v. In other words, v ∪ δv evaluates to 1 on an edge e if and only if e is oriented outwards from v.

Ground state mapping
Having introduced notation from simplicial cohomology, we are now ready to show explicitly that the ground states of the DS stabilizer model can be mapped to the ground states of the DS string-net model using a finite-depth quantum circuit. To simplify the discussion, we start with a DS stabilizer model defined on a triangular lattice equipped with a branching structure (Fig. 12a). Although the argument below assumes that the manifold is simply-connected (so that the ground state of the DS stabilizer model is unique) it can be generalized to arbitrary two-dimensional orientable manifolds straightforwardly.
We consider a DS stabilizer Hamiltonian of the form: with the terms given graphically as: , .
Note that the form of C e depends on the orientation of the edge e. It is convenient to write C e uniformly, using a cup product. In particular, C e can be written as: for any edge e. Here, e denotes the Z 2 -valued 1-cochain that evaluates to 1 on the edge e and zero otherwise. The ground state |ψ DS of the DS stabilizer Hamiltonian can be obtained by projecting the computational zero state to the ground state subspace. Letting N e denote the number of edges and |0 ⊗Ne be the computational zero state, the ground state is [97]: For simplicity, here and throughout this section, we ignore the normalization of the state. Note that it is not necessary to include projectors for the face terms B f , since |0 ⊗Ne is already in the +1 eigenspace of the face terms. After expanding the C e and A v projectors, they may be written as: where the sums on the right-hand side run over cochains c ∈ C 1 [M, Z 2 ] and a ∈ C 0 [M, Z 4 ]. The A v projector can be further simplified by acting with the Pauli Z operators on the computational zero state. More specifically, for any vertex v, A v applied to |0 ⊗Ne is equal to: where δ is the coboundary operator and v is the Z 2 -valued 0-cochain that evaluates to 1 on v and 0 otherwise [98]. Pictorially, the operator on the right-hand side of Eq. (A23) is: We now have that the ground state of the DS stabilizer Hamiltonian is equal to: where we have substituted C e in Eq. (A19) into the expression for the projector in Eq. (A21).
To compare |ψ DS to the ground state of the DS string-net model, we map |ψ DS to a system of qubits. We introduce a pair of qubits to each edge e and label the qubits by A and B. As in Section II C, the operator algebra for the pair of qubits at edge e is generated by the Pauli X and Pauli Z operators: The state on four-dimensional qudits is mapped to the system of qubits by applying the finite-depth quantum circuit U 2,4 , defined by: This unitary maps the computational zero state on qudits to the computational zero state on qubits: where the two entries of the state on the right-hand side correspond to the A and B sites, and the qubits on the left-hand side and the qudits on the right-hand side have been suppressed. The ground state in Eq. (A25) is mapped to qubits by conjugating the projectors by U 2,4 . This gives us: Note that the factors of CX AB e act trivially on the computational zero state |0, 0 ⊗Ne . Since a only appears in the exponent of order two operators, it can be replaced by a Z 2 -valued The ground state is then equivalent to: To evaluate the expression in Eq. (A29) further, we introduce cohomological notation for the computational basis states. Let We define the Z 2 -valued 1-cochains a A and a B by the conditions: The computational basis state |{a A e }, {a B e } can then be labeled by the corresponding 1-cochains, i.e.: Returning to the expression for U 2,4 |ψ DS in Eq. (A29), we apply the Pauli operators to the computational zero state. Then, using the cohomological notation for the basis states, we find: This can be reduced to a more familiar form by shifting the sum over c, so that: Making this substitution, the state U 2,4 |ψ DS becomes: The sign (−1) δb∪c(f ) in Eq. (A35) can be produced by applying a control-Z gate between the A site on an edge 12 and the B site on an edge 23 of the face f = 123 , where the ordering of the vertices is determined by the branching structure (see Appendix A 1). Letting CZ AB 12,23 denote this control-Z gate, we define U AB to be the following finite-depth quantum circuit: Applying U AB to U 2,4 |ψ DS yields: With this, the A and B sites have been disentangled. Ignoring the product state on the B sites, we are left with: The state in Eq. (A38) is mapped to the ground state of the DS string-net model by applying another finite-depth quantum circuit of control-Z gates. We let CZ AA 12,23 denote the control-Z operator between the A site on edge 12 and the A site on edge 23 . Then, we define U AA to be the finite-depth quantum circuit: where the product is over all upward pointing triangles. Acting with U AA on U AB U 2,4 |ψ DS gives the state: where Ψ(b) is the amplitude: In the next section and Refs. [96] and [99], it is argued that Ψ(b) is equivalent to: Here, N loops (δb) is the number of loops (on the dual lattice) in the configuration corresponding to δb (see Fig. 13). In other words, U AA U AB U 2,4 |ψ DS is equal to: which is precisely the ground state of the DS string-net model on the hexagonal dual lattice. Since U AA U AB U 2,4 is a finitedepth quantum circuit, the DS stabilizer Hamiltonian must belong to the same phase as the DS string-net model [52]. In summary, U 2,4 maps from four-dimensional qudits to pairs of qubits, U AB disentangles the B site qubits from the A site qubits, and U AA fixes the amplitude to be that of the DS stringnet model ground state.

String-net model ground state amplitudes
In the argument above, we claimed that the amplitudes of the state U AA U AB U 2,4 |ψ DS are precisely those of the ground state of the DS string-net model. Let us complete the argument that U AA U AB U 2,4 |ψ DS is equivalent to: To see this, we consider mapping the state to a symmetryprotected topological (SPT) state by gauging the 1-form symmetry. This can be implemented by an operator duality, (the two-dimensional Kramers-Wannier duality) which maps between a system of qubits on the edges of the lattice and a system of qubits on the vertices [28,73,74]. For a set of Z 2 values {b v }, we label the computational basis state |{b v } by the 0-cochain b. We let X v and Z v denote the Pauli X and Pauli Z operators on the vertex v. With this, the explicit mapping of operators is: where the product is over edges connected to v, and vv is an arbitrary edge. Under the duality, the state |δb is mapped to |b . Furthermore, it can be checked that U AA U AB U 2,4 |ψ DS is mapped to: Given the identification of 0-cochains b with computational basis states |{b v } , this can alternatively be written as: We now simplify the state |ψ SPT by describing a geometric interpretation for the amplitude Ψ({b v }), following Ref. [100]. On the triangular lattice with the branching structure shown in Fig. 6b, Ψ({b v }) is equivalent to: To arrive at the expression above, we used the explicit formula for the cup product in Appendix A 1. The amplitude in Eq. (A48) can be further reduced to the form: where V 1 (b), E 1 (b), and F 1 (b) are the number of vertices, edges, and faces, respectively, contained entirely within the domains formed by the b v = 1 vertices in the configuration {b v }. We let Σ 1 (b) denote the domains formed by the vertices with b v = 1. Then, Ψ(b) can be expressed in terms of the Euler characteristic χ(Σ 1 (b)): The Euler characteristic of an orientable surface Σ is: where g(Σ) is the genus of Σ and n(Σ) is the number of boundary components of Σ. Therefore, letting N dw (b) be the number of domain walls in the configuration associated with b, the amplitude Ψ(b) is equivalent to: This gives us: The final step of the argument is to gauge the 0-form symmetry of the SPT state. This should return us to the state U AA U AB U 2,4 |ψ DS . The 0-form symmetry is gauged by applying the duality in Eq. (A45), leaving us with: Since N dw (b) is the same as the number of loops N loops (δb) in the configuration corresponding to δb, we have the desired result: Appendix B: K-matrix formulation of twisted quantum doubles For completeness, we describe Abelian TQDs and the construction of the associated stabilizer models in terms of the K-matrix formalism. The K-matrix formalism is based on an integer symmetric matrix K, whose matrix elements encode the couplings between U (1) Chern-Simons theories. We refer to Refs. [101][102][103] for further details on the connection to Chern-Simons theories. What is important for our discussion is that the universal properties of the Abelian anyon theory can be deduced entirely from the matrix K.
We restrict our focus to K-matrices for Abelian TQDs. To this end, we take G to be the finite Abelian group of the form G = M i=1 Z Ni and recall that the cocycle of an Abelian TQD can be labeled by a set of integers I = {n i } M i=1 ∪ {n ij } M i,j=1 , as described in Section IV A. We define the K matrix associated to the TQD to be the 2M × 2M matrix: where N is the M × M diagonal matrix: and S I is the M × M symmetric matrix: 2n 1 n 12 n 13 · · · n 1M n 12 2n 2 n 23 · · · n 2M n 13 n 23 2n 3 · · · n 3M . . . . . . . . . . . . . . .
We note that the ij matrix element of S I is precisely n ij . This choice of K TQD has the property that the gauge charges c i and (some choice of) elementary fluxes ϕ i correspond to unit vectors l ci and l ϕi of length 2M . In particular, l ci and l ϕi are the unit vectors with a 1 in the (M + i) th or i th entry, respectively. Further, since every anyon a in the TQD can be generated by c i and ϕ i : we can assign an integer vector l a to each anyon a such that the i th entry is p i and the (M + i) th entry is q i : l T a = (p 1 p 2 · · · p M q 1 q 2 · · · q M ).
For any two anyons a and a , the vectors l a and l a satisfy: There is some ambiguity in the assignment of integer vectors to anyons within the K-matrix formalism. To see this, notice that the columns of K encode the relations of the gauge charges and elementary fluxes. For example, the first column corresponds to the relation c N1 1 = 1, while the (M + 1) th column tells us: Therefore, we may freely redefine the vectors l a by integer multiples of the columns of K, i.e.: l a ∼ l a + j m j col j (K) where m j are integers and col j (K) is the j th column of K.
The exchange statistics and braiding relations for anyons a and a can be extracted from the inverse of K using the formulas: θ(a) = e πil T a K −1 la , B θ (a, a ) = e 2πil T a K −1 l a . (B9) In terms of the quadratic form q and associated bilinear function b q (both defined in Section III), we have: The relations in Eq. (B10) can be checked by computing the inverse of the K TQD , which we find to be From which, we obtain: These agree with the statistics and braiding relations of the gauge charges and elementary fluxes in Eq. (85).
We now see that all of the characteristic properties of the TQD can be determined from the matrix K TQD . In particular, the fusion rules of the anyons are given by Eqs. (B6) and (B8), and the statistics and braiding can be computed from Eq. (B10).
Before turning to examples, we clarify some freedom in the definition of a K-matrix K. The freedom comes from choosing a different set of generators for the anyons. For example, for K TQD , the generators are implicitly chosen to be the gauge charges and a set of elementary fluxes. A different set of generators corresponds to transforming the vectors l a as W l a , where W is an integer matrix belonging to GL(2M, Z). To preserve the statistics and braiding, K −1 must transform as: which means that K is mapped according to: Therefore, any two K-matrices related by a transformation as in Eq. (B14) describe the same anyon theory.

Examples
We are now prepared to consider examples, starting with the Z N TC. In this case, K is the 2 × 2 matrix: The anyons of the Z N TC are generated by e and m, which correspond to the unit vectors: Let us now explore examples where the fusion group of the anyons differs from G × G. First, we consider the TQD with G = Z 3 and n 1 = 1. This corresponds to the matrix: Interestingly, the anyons in this theory have Z 9 fusion rules. The generating anyon of order nine can be chosen to be the elementary flux l T ϕ1 = (1 0), with q(ϕ 1 ) = 1/9. The excitation ϕ 3 1 is a boson and corresponds to the charge c 1 of order 3. This can be seen by noticing that: Next, we consider the TQD with G = Z 2 ×Z 2 and n 12 = 1. Notably, this TQD is characterized by a type II cocycle. Given the general form for K, we find: Therefore, the anyons have Z 4 ×Z 4 fusion rules. Furthermore, the elementary fluxes have nontrivial mutual braiding: In fact, this TQD is equivalent to the Z 4 TC by the identification ϕ 1 = e and ϕ 2 = m. To see this explicitly, we can perform a basis transformation to relabel the anyons as those of the Z 4 TC. To do so, we define the matrix W , given by: Thus, the theory describes a decoupled Z 4 TC and Z 1 TC, where the latter factor of Z 1 can be ignored, as it has no anyonic excitations. With this, we have mapped the anyon content of the TQD to that of the Z 4 TC. Similarly, if n 12 = 1 and either n 1 = 1 or n 2 = 1, but not both, then the TQD is again equivalent to the Z 4 TC. As a final example, let us consider the TQD with G = Z 2 × Z 2 and n 1 = n 2 = n 12 = 1. In this case, the K matrix is: (B28) Similar to the previous example, this anyon theory has Z 4 ×Z 4 fusion rules. However, in contrast, the elementary fluxes are semions. All together, this theory has six semions, six antisemions, and four bosons. This differs from other Z 2 × Z 2 TQDs, which all have at least one fermionic excitation. We refer to this theory as the six-semion theory, since it is closely related to the three-fermion theory, according to the classification of Abelian anyon theories in Ref. [57]. Specifically, they both belong to the family of anyon theories denoted by F 2 r in Ref. [57], with r = 1 and r = 2 corresponding to the three-fermion theory and six-semion theory, respectively. For the anyon theories in F 2 r , the fusion group is Z 2 r × Z 2 r and the quadratic form is given by: q((a 1 , a 2 )) = a 2 1 + a 2 2 + a 1 a 2 2 r , where (a 1 , a 2 ) is an element of Z 2 r × Z 2 r . The chiral central charge is 4 mod 8 if r is odd and 0 mod 8 if r is even.

Construction of TQDs from boson condensation
The general construction of TQDs starting from TCs, as described in Sec. IV B, can be concisely stated in terms of K-matrices. To see this, we begin with a stack of TCs, given by the following 2M × 2M K-matrix: The next step in the construction is to condense the set of bosons in Eq.
We can verify that the b i are bosons with trivial mutual braiding by noting that: q(b i ) = 1 2 l T bi K −1 T C l bi = 0 mod 1, ∀i, b q (b i , b j ) = l T bi K −1 T C l bj = 0 mod 1, ∀i, j. (B32) The vectors of the bosons can be compiled into a 2M × M matrix Q, where the vectors l bi form the columns of Q: Here, U I is the M × M upper triangular matrix: n 1 n 12 n 13 · · · n 1M n 2 n 23 · · · n 2M n 3 · · · n 3M . . . . . .
Indeed, Eq. (B32) is satisfied by noting that The deconfined anyons, i.e., those that braid trivially with the bosons b i , are generated by the columns of the following matrix: where I M ×M is the M × M identity matrix. Note that in Eq. (B36), the first M columns correspond to the elementary fluxes of Eq. (101), and the following M columns correspond to the gauge charges of Eq. (100). We see that confirming that the deconfined anyons braid trivially with the condensed bosons. Note that the exchange statistics and braiding relations of the deconfined anyons are captured by the matrix L T K −1 T C L. Therefore, after condensation, the corresponding K-matrix is given by L −1 K T C (L −1 ) T . The matrix L −1 K T C (L −1 ) T is precisely the K-matrix K T QD , given in Eq. (B1).
the group structure by considering the fusion of elementary fluxes.
To start, we choose an elementary flux ϕ i for every i in {1, . . . , M }. After choosing the elementary fluxes, each anyon can be expressed as a unique product of elementary fluxes and gauge charges. That is, the anyons can be written in the form: