Non-Pauli topological stabilizer codes from twisted quantum doubles

It has long been known that long-ranged entangled topological phases can be exploited to protect quantum information against unwanted local errors. Indeed, conditions for intrinsic topological order are reminiscent of criteria for faithful quantum error correction. At the same time, the promise of using general topological orders for practical error correction remains largely unfulfilled to date. In this work, we significantly contribute to establishing such a connection by showing that Abelian twisted quantum double models can be used for quantum error correction. By exploiting the group cohomological data sitting at the heart of these lattice models, we transmute the terms of these Hamiltonians into full-rank, pairwise commuting operators, defining commuting stabilizers. The resulting codes are defined by non-Pauli commuting stabilizers, with local systems that can either be qubits or higher dimensional quantum systems. Thus, this work establishes a new connection between condensed matter physics and quantum information theory, and constructs tools to systematically devise new topological quantum error correcting codes beyond toric or surface code models.

It has long been known that long-ranged entangled topological phases can be exploited to protect quantum information against unwanted local errors. Indeed, conditions for intrinsic topological order are reminiscent of criteria for faithful quantum error correction. At the same time, the promise of using general topological orders for practical error correction remains largely unfulfilled to date. In this work, we significantly contribute to establishing such a connection by showing that Abelian twisted quantum double models can be used for quantum error correction. By exploiting the group cohomological data sitting at the heart of these lattice models, we transmute the terms of these Hamiltonians into full-rank, pairwise commuting operators, defining commuting stabilizers. The resulting codes are defined by commuting non-Pauli stabilizers, with local systems that can either be qubits or higher dimensional quantum systems. Thus, this work establishes a new connection between condensed matter physics and quantum information theory, and constructs tools to systematically devise new topological quantum error correcting codes beyond toric or surface code models.

I. INTRODUCTION
Any architecture proposed for information storage must be equipped with an error correction strategy to avoid the corruption of the data encoded, whether the information is classical or quantum in nature [1][2][3]. Since the no-cloning theorem [4] prevents qubits from being copied, quantum error correction cannot rely on data redundancy at any point. Thankfully, the fact that errors are usually local, i.e., they affect a small number of qubits, has lead to fruitful alternative strategies. By distributing the relevant data over a whole system, it is possible to detect the errors without ever needing to copy the original state.
Building from this insight, stabilizer codes [5] have taken a particularly prominent role in the search for encoding strategies for scalable and fault-tolerant quantum computing. In stabilizer codes, the subspace in which the quantum information is stored is the joint eigenspace of pairwise commuting operators, called stabilizers. Among these are a class of codes -so called topological codes -where error detection can be performed with the measurement of local stabilizers. These measurement outcomes, repackaged into syndromes, determine the errors that have occurred. By construction of such codes, the measurement does not destroy the stored quantum information and makes it possible to restore it with a suitable error correction scheme [5,6]. The toric code [7], its associated planar embedding known as the surface code [8], and color codes [9], are by far the most studied codes, and have emerged as the gold standard of this class of error-correction protocols. Their simple construction -with stabilisers built out of Pauli words -means that they collectively provide a wide range of easily understood schemes. That said, there are strong reasons to seek for new codes beyond these Pauli stabilizer models. While the lack of a universal and fault-tolerant gate set -by virtue of the Eastin-Knill theorem [4] -and a lack of self-correctability [10] will be common to any stabilizer * julio.carlos.magdalena@rwth-aachen.de approach in two spatial dimensions, several techniques have already been identified to circumvent these limitations. Both magic state distillation [11] and just-in-time decoders [12,13] give rise to universal computational power at the cost of some overhead, depending strongly on the specific code architecture chosen [14]. Moreover, codes with a transversal Clifford gate built out of d-level systems have been found to have superior error correction capabilities compared to qubit-based codes, with an increasing performance with increasing d [15] or enhanced bit flip stability [16]. Other generalisations involving non-commuting stabiliser sets [17] have demonstrated the ability to produce gate sets which, while not universal, have enhanced computation power. Taken together, these findings strongly motivate the quest for new topological quantum error correction codes with stabilizers outside the Pauli group that may be better-suited to practical implementations and for which the overheads are more manageable.
In light of this search, we present a wealth of new topological codes. To do so, we have taken inspiration from the closely related field of topological phases of matter. The conditions for quantum error correction, the Knill-Laflamme criteria [18], are highly reminiscent of conditions for the topological order in quantum many-body theory. However, this connection is rarely made explicit beyond the toric code, which can be seen as defining a gapped Hamiltonian with 4 anyon types and a topological ground state degeneracy. While it is true that all topological error correcting codes can ultimately be understood as defining a system containing anyonic excitations and therefore being in a topological phase, all wellstudied instances of this are equivalent to multiple copies of the toric code phase [19]. What is sorely lacking in this picture is a way of reverse-engineering topological quantum error correcting codes from the wealth of topological phases of matter. This seems a remarkable omission in the light of the powerful and highly developed classification of such phases from the perspective of condensed matter and mathematical physics [20][21][22]. This omission is also significant given the fact that, from a technological perspective [23], the identification of new topological codes seems imperative.
In this work, we use a large class of topological orders host-ing Abelian anyons to construct new topological error correcting codes. In particular, we modify existing lattice models for topological orders -twisted quantum double models -so that they give rise to stabilizers. In their original form, the local terms of these Hamiltonians do not commute in a particular excited subspace of the Hilbert space, which makes them -on first glance -unsuitable for stabilizer error correction. Practically speaking, commutativity is a highly desirable property in the context of quantum error correction, in that it allows for error correction schemes based on independent local measurements of such stabilizers without perturbing the stored quantum information. We restore commutativity by first deriving the quantities that obstruct this property from the group cohomology data of twisted gauge theories. In most cases -namely, for Abelian twisted quantum doubles -these obstructions can be lifted completely by carefully modifying the offending terms in the Hamiltonian, yielding a true stabilizer code, consisting of commuting non-Pauli operators. A first step in this direction was taken in Ref. [24], where the double semion string-net model [20] was modified with a local phase factor to overcome the same commutativity problem. However, this approach lacks a systematic and quantitative understanding of the failure of commutativity, and as such it cannot be generalized to other lattice models for more exotic topological orders. Our results go a significant step further, providing a robust framework for deriving quantum error correcting codes from not only a Hamiltonian in the double semion phase, but from a huge family of Abelian phases as well. This paper is structured as follows: In Section II we give a comprehensive introduction to twisted quantum double models for general (Abelian) groups. We have kept the mathematical details to a minimum while still presenting our results in a self-contained manner. Our construction is done explicitly for a Z 2 and a Z 2 × Z 2 model and then summarized for the general Z N and Z 2 N cases in Sec. III, with full details in Apps. B-D. Moreover, we give a brief overview on the properties of the newly constructed codes and how they relate to known schemes for topological quantum error correction. In Sec. IV, we conclude our work and give an outlook on future directions and potential use cases of the codes and discuss potential applications in topological quantum information processing.

II. INTRODUCTION INTO TWISTED QUANTUM DOUBLE MODELS
Twisted quantum double (TQD) models are lattice models for topological order in 2+1 dimensions which can be viewed as a generalization of the quantum double model [7] introduced by Kitaev. They can be obtained by promoting the global symmetry of a symmetry protected topological (SPT) phase [22,25,26] to a local gauge symmetry via minimal coupling to the original "spins" of the SPT. We will restrict the discussion of the model only to the aspects necessary for our construction. An interested reader is referred to Ref. [21] for more comprehensive perspective. Around each vertex, the edges directly adjacent to it are labelled from l1 to l6 and the other 6 edges that share a triangle with the vertex are labelled by l7 to l12. Together, these edges constitute the neighborhood of the vertex and are marked in red above.

A. The Hamiltonian
We define our model on a translation-invariant, oriented triangulation of a general compact surface shown in Fig. 1. We label the edges in the neighborhood of a vertex v from l 1 to l 12 . Each edge l i carries a degree of freedom (gauge field) whose local Hilbert space H l is spanned by states labeled by elements of a finite group G, Its local dimension is |G|, so a group with |G| = 2 will be a qubit model , |G| = 3 a qutrit model and so on. The total Hilbert space is then simply given by H = edges l H l . While G can be chosen to be any finite group, we will only be treating Abelian cases in this work Now that we have a Hilbert space, we can define a Hamiltonian on it. By keeping to the basis used in Eq. (1), the Hamiltonian terms have straightforward descriptions in terms of their matrix elements. This Hamiltonian is given by The first sum runs over all (triangular) faces and the plaquette operator acting on a triangular face is defined by where δ g = 1 for g = 1 G (identity element in G) and δ g = 0 otherwise, with ". . . " standing in for the rest of the graph on which B p acts trivially. It projects out G-fluxes through the face p and ensures that the ground space is flux free, i.e., B p = 1 ∀p. For any finite group, this projector can be expressed in terms of group characters χ i , where i labels the conjugacy Second, it scales the wavefunction by a phase factor given by the product of 6 special functions ω : G 3 → U (1) called cocycles, with one per triangle adjacent to v. The order of the arguments in the cocycles and whether they appear in the numerator or denominator of this prefactor is determined by the orientation structure of the lattice. For a detailed explanation of constructing the prefactor for a general lattice, see Ref. [21]. The cocycles encode the topological data of the theory modeled by H T QD . Their defining property is the socalled cocycle condition ω(g 1 , g 2 , g 3 )ω(g 0 , g 1 · g 2 , g 3 )ω(g 0 , g 1 , g 2 ) ω(g 0 · g 1 , g 2 , g 3 )ω(g 0 , g 1 , g 2 · g 3 ) = 1 Obviously, ω(a, b, c) = 1 is always a solution and is called trivial. If we use this trivial solution in Eq. (6) to define H T QD , we obtain the quantum double Hamiltonian from [28].
Since -in general -there are non-trivial solutions to this equa-tion, the TQD model covers a much broader class of Hamiltonians than the pure quantum double Hamiltonians. In principle, one can choose any function satisfying condition (8), insert it into Eq. (6) and obtain a consistent model with topological order. However, not all solutions yield distinct orders. A close investigation of the cocycle condition reveals that if we have one solution ω, we can always obtain another solutioñ where β : G 2 → U (1) is an arbitrary function mapping two group elements to a phase factor. If we have two TQD Hamiltonians defined by two cocycles ω 1 and ω 2 in Eq. (6) so that they are in different topological orders, we know and there exists no β to map ω 1 onto ω 2 by Eq. (9). Hence, inequivalent Hamiltonians H T QD are classified by distinct equivalence classes of functions ω, which define elements of the third group cohomology of G over U (1) with ω ∼ω iff ∃β : G 2 → U (1) such that they are related by Eq. (9). In the next section, we will see examples of such functions for simple groups such as Z 2 and Z 2 × Z 2 . For an introduction into group cohomology, see App. A.

B. Topological data
The Hamiltonian we have constructed on a triangulation of a compact surface from a group G together with a cocycle ω is indeed topologically ordered. It has anyonic excitations and a robust ground state degeneracy (GSD). For Abelian models, GSD = |G| 2g where g is the genus of the surface on which it is defined. This causes a non-zero topological entanglement entropy, a characteristic feature of topologically ordered systems. The topological quantum numbers (topological spin and S-matrix) of the excitations are uniquely defined by the input group G and the cocycle one chooses to define A v . Moreover, they are gauge invariant in the sense that one can choose any cocycle in the same equivalence class to define A v and still obtain the same topological data. This corresponds to a transformation of the cocycles like in Eq. (9). For the derivation and explicit expressions of those quantities in terms of the input data see Ref. [21].

C. Ground space and failure of commutativity
The ground space of H T QD can be found exactly and is defined implicitly by the conditions so that it is the simultaneous eigenspace of the plaquette and vertex operators. As stated before, this subspace has dimension larger than 1 on a surface with non-trivial topology and therefore we hope to use that space as a code space of an error correction stabilizer code. Unfortunately, the Hamiltonian is not exactly solvable on the whole Hilbert space, i.e. one cannot simultaneously diagonalize all vertex and plaquette operators. In particular, the vertex operators fail to commute in the presence of certain fluxes (B p = 0). 1 In the TQD model for G = Z 2 and the non-trivial cocycle [21] ω 1 (1, 1, 1) = −1, ω(a, b, c) = 1 else, a version of the double semion phase, the vertex operators do not commute when acting on the following configuration: In the case of G = Z 2 , the local Hilbert space corresponds to the one of a qubit which we represented by a circle in the state vectors above. We have labelled the state vectors of the qubits with circles, |0 = • and |1 = •. The rest of the lattice is not explicitly shown, but we assume that all other qubits are in |0 . The vertex operators will be defined in Sec. III A where we discuss the double semion phase in detail.
From an error correction perspective, there are many (single qudit) errors that create such fluxes and for the vertex operators to be proper stabilizers they also have to commute in that sector of the Hilbert space. In principle, one can make the Hamiltonian exactly solvable by multiplying it with a projector on the flux-free subspace, P B=1 = p B p . However, we will loose information about the excited sector in doing so in the sense that it makes it impossible to identify an excitation uniquely by measuring those local operators. Luckily, it turns out that a slight modification of the vertex operators that does not change the topological order resolves the obstruction of commutativity entirely.

III. CONSTRUCTION OF FULLY COMMUTING MODELS
In order to successfully overcome the obstacle to commutativity, the vertex operators have to be modified. This modification should not alter the topological phase of the Hamiltonian. This means that the vertex operators must be altered in such a way so that the ground space is left unchanged and that the spectral gap is preserved. At the same time, the modification should be minimal. The latter means that the modification should be local -in the sense of not increasing the support at all -and leave the large sections of the Hilbert space where the Hamiltonian is solvable undisturbed. In particular, only for certain flux configurations do the operators need to be minimally altered.
It constitutes the main result of this work that, for large classes of twisted quantum doubles, the above desiderata can be achieved with a modification of the form, where η g v is a full-rank operator that is diagonal in the edge basis, with entries being roots of unity, and is equal to 1 on the ground space.
By imposing that the modified operators {Ã g v } commute, we obtain consistency equations for the modification phases The derivation of these consistency equations is only possible with the TQD model and exploits the machinery of group cohomology in its construction. We will illustrate the procedure of solving these equations by means of two simple examples where the input group is Z 2 and Z 2 × Z 2 . In the second case, we investigate a topological order that is entirely new to the context of quantum error correction, but is still in principle realizable with a qubit architecture. We apply the same procedure to more exotic models derived from the groups Z N and Z 2 N . With those completed, we have resolved the commutativity issue for every Abelian topological order that can be obtained from a twisted quantum double model. These results are described in the last part of this section. For the general formalism and the calculation see App. B, C, and D.
The newly obtained operatorsÃ g v together with the plaquette phase operators B g p (see Eq. (4)) generate the Non-Pauli stabilizer group which defines the corresponding code. Although stabilizer groups are normally defined as a subgroup of the Pauli group, this restriction is not necessary for the general concepts of stabilizer error correction. These require only that there exist an invariant subspace under the action of those generalized stabilizers, namely the common +1 eigenspace of all B g p andÃ g v . In the following, we will illustrate the construction for an arbitrary finite group G and a cocycle ω for two simple examples that already go beyond the toric code. We will see that we arrive at a wealth of new topological quantum error correcting codes produced directly from twisted quantum double models for topological order.

A. Z2 -double semion code
We first investigate the non-trivial Z 2 model that is in the same phase as the double semion string-net model [20]. We represent Z 2 as the set {0, 1} together with the group operation being addition modulo two. With that, H l = span C {|0 , |1 } = C 2 , so it is a model of interacting qubits. Z 2 has two inequivalent cocycle classes, one trivial class [ω 0 ≡ 1] and one non-trivial class [ω 1 ]. The TQD model with the trivial cocycle would yield a Hamiltonian in the Toric Code phase. The canonical representative of the non-trivial class is given by Eq. (12). Inserting this cocycle into Eq. (6) yields the Hamiltonian where Z l and X l are the Pauli z and x matrices on (the qubit sitting on) edge l and P ± l = 1 2 (1 ± Z l ) is the projector onto the space where edge l caries the value 0 or 1, respectively. This Hamiltonian is in the so-called double semion phase.
By construction, A 1 v always flips an even number of qubits adjacent to a plaquette, and thus A 1 v , B p = 0 ∀v, p. 2 Hence, the only obstruction for the operators B p and A 1 v to form a commuting set of operators that we can use for stabilizer error correction comes from the vertex operators A 1 v .

Obstruction in the original model
For the operators A 1 v to generate a proper stabilizer group they have to represent the group action of Z 2 on site. In particular, any representative of an element in Z 2 should square to the identity. Unfortunately, it turns out that where we have used X 2 l = 1, and the decomposition of the identity, 1 = P + l + P − l for any edge l. In fact, where B i,j,k,l = 1 2 (1 − Z i Z j Z k Z l ) measures the flux through the region enclosed by the edges {i, j, k, l}. This shows that the operators fail to represent the group action on the part of the Hilbert space where B 3,9,8,1 + B 6,12,11,4 = 1 mod 2. In particular, on the ground space (in which no flux is present) the group action is implemented correctly. This is not the only obstruction in the original TQD model. To form stabilizers, the vertex operators A g v must commute pairwise. Due to the translational invariance of the model, we need only calculate the commutation relation between A 1 v with the three operators {A 1 i , i = 1, 2, 3} acting on the three vertices {1, 2, 3} connected to v by the edges {l 1 , l 2 , l 3 } (see Fig. 1 for the labelling) to confirm this. It turns out that using the same relations used to produce Eq. (16) and (17). Again, we find that they commute in the zero-flux sector of the Hilbert space. Interestingly, vertex operators on neighboring vertices only fail to commute in the last case, when they are connected by a horizontal edge (labeled by l 2 in Fig. 1) which is neighboring a nontrivial flux B 3,9,8,1 = 1 mod 2.
This particular "locality" of the commutativity obstruction is a consequence of our chosen edge orientation 3 which determines which arguments enter in the cocycles in Eq. (6).

Modifying vertex operators by local phase
We have found that the vertex operators in the original TQD model fail to be proper stabilizers on the whole Hilbert space because, on one hand, they do not implement the group action on site, i.e. (A 1 v ) 2 = 1, and, on the other, fail to commute However, we were able to quantify the obstructions and found that they take a very particular form and only depend on fluxes. Because of this, there are no obstructions in the ground space. To remove the obstructions on the whole Hilbert space, we want to modify the vertex operators by a local phase as described in the beginning of Sec. III, so thatÃ 1 v are proper stabilizers. The modified operators should therefore satisfy for all v, v 1 , 2 ∈ V (Γ). Imposing Eq. (20a) and using Eq. (17) yields the condition on the modification phase In fact, this can be solved by which is solved bȳ is the operator that measures the value of the edge l i in Z 2 . Note thatη is a periodic function in p with periodicity 4, i.e.η(p + 4) =η(p). Also, the second factor containing i only depends on fluxes and ensures that Combined, we obtain the modification phase where the first factor ensures that the group property is fulfilled on site and the second factor fixes the commutativity on the whole Hilbert space. Since η(p + 4) = η(p), p ∈ [0, 4) parametrizes all the distinct modification phases in this family of solutions. The geometric support structure for general p is depicted in Fig. 2. Note that the parameter p sets the dependence on l 2 and l 5 in the second factor. The freedom to choose p ∈ [0, 4) may be useful in an actual error correction scheme since different modification phases yield different stabilizers that in turn could have different properties in the decoding process. For p = 1, for example, we obtain 9,8,1+B6,12,11,4 (−1) B3, 9,8,1l2 =i B6,12,11,4 (−1) B3,9,8,1 (−1) B3, 9,8,1l2 , so that it does not depend on l 5 . Explicitly quantifying the obstructions of the group property (Eq. (17)) and the commutativity of the operators acting on neighboring vertices (Eqs. (18) in the original TQD model enabled us to remove them with a local phase modification such that the modified operatorsÃ 1 v = η v A 1 v faithfully represent the group Z 2 on the whole Hilbert space. The constructed operators are outside the Pauli group and can be used as stabilizers in the context of quantum error correction. Moreover, the modification does not change the action of A 1 v on the ground space, and thus the modified Hamiltonian it still is in the double semion phase.
In the previous section, we constructed a set of stabilizers defined on a lattice of qubits such that its code space corresponds to a double semion ground space. However, the double semion phase is not the only twisted gauge theory one can implement with a qubit architecture. By taking G = Z 2 × Z 2 = {g = (g 1 , g 2 ); g 1 , g 2 ∈ Z 2 } as an input group for the TQD model, the local Hilbert space becomes H l = span C {|0, 0 , |0, 1 , |1, 0 , |1, 1 } C 2 ⊗ C 2 which can be realized using two qubits per edge. As was shown in Ref. [29], an untwisted Z 2 × Z 2 quantum double model -otherwise known as two copies of the toric code -is equivalent to the color code. In this section we investigate the twisted versions thereof, which we call twisted color codes.
The possible topological orders of a TQD model with G = and thus there are 8 different cocycle classes we can choose as input, labelled by (s 1 , s 2 , s 3 ) ∈ Z 3 2 . In an appropriate gauge, a general cocycle ω ∈ H 3 (Z 2 × Z 2 , U (1)) can be written as with s i = 0, 1 and the group elements are represented by pairs of Z 2 variables, i.e., a = (a 1 , a 2 ). When s 3 = 0, only the cocycles ω 1 appear. They are the same as those seen for the Z 2 phases, defined in Eq (12), only now they depend explicitly on a particular tensor factor, and are referred to as type-I cocycles. Cocycles of that type yield TQD models describing a topological order that is decomposable into Z 2 phases. For example, choosing (s 1 , s 2 , s 3 ) = (1, 1, 0) produces a Hamiltonian describing a product of two double semion phases. In this case, one can make the vertex operators from each copy fully commuting using the same phase modification derived in the previous section. When s 3 = 1, we have a cocycle that can be represented by which mixes two tensor factors and therefore is unique to the Z 2 × Z 2 case. 4 To distinguish it from the previously studied cocycles, it is referred to as a type-II cocycle. A TQD model with such a cocycle as input requires a different modification of the Hamiltonian, which we construct in this section. The TQD Hamiltonian built by inserting the type-II cocycle from Eq. (28) into Eq. (6) reads

and (29c)
where X p , A g v ] = 0 ∀g ∈ Z 2 × Z 2 since each vertex operator flips an even number of qubits around each plaquette. The only obstructions preventing from forming a pairwise commuting set come from the vertex operators. We will quantify the obstructions below.

Obstructions in the original model
Each element in Z 2 × Z 2 is its own inverse. For the vertex operators to generate a proper stabilizer group, they must represent the group action on site and therefore also square to 1. Since the representative we chose for the type-II cocycle in Eq. (28) does not depend on a 2 , b 1 and c 1 , we find and (30a) For the other products of non-trivial group elements (0, 1), (1, 0) and (1, 1) however, we obtain explicit obstruc-tions 6,12,11,4 1 = 1, (31b) m ) measures the flux in the ith tensor factor through the region enclosed by the edges {j, k, l, m}. As in the Z 2 case, the operators square to the identity in the flux-free subspace which includes the ground space.
In addition to the on-site obstructions, the vertex operators fail to commute for neighboring vertices. Due to translation invariance, we need only calculate the commutation relation between A g v and the three operators {A h i , i = 1, 2, 3} connected to v by the edges {l 1 , l 2 , l 3 } (see Fig. 1 for the labelling) for any pair (g, h). We obtain the commutativity relations with the remaining pairs either commuting or not giving independent obstruction phases(see App. D). In particular, we find that only the vertex operators acting on different tensor factors fail to commute and they fail precisely when the vertices on which they act are connected by a horizontal edge neighboring a flux. In fact, this is a general property of our model, and can be traced back to the original choice of edge orientations (see App. B).

Modifying vertex operators by a local phase
We have found that in the original Z 2 × Z 2 TQD model the vertex operators fail to be proper stabilizers because, on one hand, the group action is not represented correctly on-site, i.e.
in general, and, on the other, some of them fail to commute. We were able to quantify the obstructions and found that they have a similar structure as in the Z 2 TQD model, namely factors of −1 that only depend on fluxes. To resolve the obstructions for the three operators A , we modify them by a local phase η g v that is the identity on the ground space (Eq. (13)). For the modified operators to be stabilizers, we need them to fulfill The first two conditions are on-site conditions reflecting that the vertex operators should form a representation of the Abelian input group and the third condition is the commutativity condition necessary for the vertex operators to be stabilizers.
Condition (33a) gives us independent constraints on each of the generating phases η (0,1) v and η (1,0) v . Using Eq. (31a), we find that the first phase must satisfy which is solved by any solution of the form Eq. (33b) to express η A close inspection of this equation shows that this equation is satisfied byη where l The remaining modification phase can be calculated using Eq. (33b) to obtain 6,12,11,4 .
These solutions define the modified operators which form a faithful representation of Z 2 × Z 2 . When inserting these operators into Eq. (33c), one finds that they also commute for neighboring vertices which completes our construction of non-Pauli stabilizers {Ã (1,1) v } based on the Z 2 × Z 2 TQD model constructed with a type-II cocycle. Note that -unlike in the Z 2 case -we have not found a one-parameter family of solutions, though we do not claim that our solution is unique. Thankfully, the modification phases we derived here are already quite simple in form, depending only on a restricted neighborhood of the vertex, illustrated in Fig. 3. The modification phase of A With these modifications, we have constructed a set of stabilizers whose code space is the ground space of a Z 2 × Z 2 topological order that cannot be factored into two (possibly twisted) Z 2 phases. It is just one example of how our analysis of the on-site and commutativity obstructions in a TQD model allows us to obtain stabilizers from various topological orders, since the techniques used here can be extended to more general models.

C. General Abelian topological order
In the previous subsections, we explicitly calculated and corrected the obstructions to the construction of stabilizers from qubit-based TQD models. However, these are just two exemplary cases of our result for any local dimension, i.e. for the input groups Z N = {0, 1, . . . , N − 1} and Z N × Z N with type-I and type-II cocycles respectively. For those topological orders, the construction of the modification phases follows a similar line as that of the qubit based models. In the following, we will sketch the general construction and state the resulting modification phases for Z N and Z N × Z N . The detailed calculation can be found in App. C and D.
The action of the vertex operator in terms of cocycles (Eq. (6)) allows us to quantify the on-site and the commutativity obstructions for a generic TQD model. Moreover, we can derive consistency equations for the modification phase by imposing that the vertex operators should represent the group action on-site and commute for neighboring vertices (see App. B). Using the canonical representative of a type-I cocycle for Z N , where [a + b] N = (a + b) mod N , allows us to explicitly solve the consistency equations and obtain pairwise commuting vertex operators {Ã g v } that represent the group action of Z N . We exploit the cyclicity of Z N by imposing that every vertex operator to the N th power should equal the identity, just as in Eq. (20a) for N = 2. In particular, this should hold for the generating vertex operator A 1 v , which allows us to determine a suitable ansatz for the corresponding modification phase η 1 v . From this, we find a family of solutions η 1 v (p) that ensure that the generating vertex operators not only represent the group but also commute pairwise. The fact thatÃ 1 v generates every otherÃ g v allows us to compute every other modification phase η g v iteratively. The resulting modification phase for any g ∈ Z N reads i,j,k is the projector onto the space in which the sum of edge values l i + l j + l k = n mod N and a minus sign in front of an index states that the inverse element enters in the sum. For example, the projector P (n) −2,8,1 projects onto the space in which −l 2 + l 8 + l 1 = n mod N . The flux B 6,12,11,4 = (l 6 + l 12 − l 11 − l 4 ) mod N is defined in a similar fashion as in the Z 2 case. Unlike the Z 2 case, we have to take the orientation of the edges into account by subtracting l 11 and l 4 since Z N is not an involutory group. We note that the modification phase above is 1 in the zero-flux subspace, where B 6,12,11,4 = 0 and all the projectors appearing cancel pairwise in the above expression. The first term is a N 2 th root of unity and only depends on fluxes. It is the higher dimensional analogue of the first factor in the Z 2 solution, Eq. (26). The second term, which only includes the edge value l 2 in the upper summation bound, reduces to the second term in Eq. (26) when N = 2. The only term where l 2 enters in the argument of the projectors is the third term. For N = 2, this term reduces to the last term in Eq. (26). For a detailed derivation of the modification phase for G = Z N we refer to App. B and C.
When considering the TQD model with gauge group Z 2 N and a type-II cocycle (Eq. (28)), the construction follows a similar path. For each tensor factor, the vertex operators must fulfill the same closure relation as in Z N . This allows us to find suitable ansatzes for the modification phases η } represent the group action of the two generators (0, 1) and (1, 0) on-site in a consistent fashion on both tensor factors and commute pairwise. Since they correspond to the two generators of Z N × Z N , we can again iteratively construct the modification phase for any g = (g 1 , g 2 ) ∈ Z 2 N , where we indicate the operators only supported on the ith tensor factor with an upper index in brackets, e.g. l  (i,j,k) (l) are defined as in the Z N case but on the lth tensor factor. As with the Z N solution, we can identify the terms derived in the previous subsection for N = 2 (see Eqs. (39) and (38)). The first term only depends on fluxes and is non-trivial when g 2 = 0. In contrast, the second term depends on both tensor factors of g, where g 1 defines which N th root of unity is appended and g 2 selects which projectors appear in the sum. The explicit edge value l for Z N and Z N × Z N , these results readily generalize to any Abelian group yielding an Abelian anyon theory. Since, by the fundamental theorem of finitely generated Abelian groups, any finitely generated Abelian group can be decomposed into Z N factors. As a consequence of this, one can also decompose the cocycles of an Abelian group into cocycles of of this cyclic decomposition [30]. While there exist cocycles beyond the type-I and type-II classes considered here, their inclusion in a TQD model produces topological order containing non-Abelian anyons, and are therefore not suitable for stabilizer error correction based on commuting syndrome measurements. In our construction, we have explicitly constructed the modification phase for type-I and type-II cocycles. Since the cocycle defining an arbitrary Abelian TQD model will be a product of type-I and type-II cocycles, the modification phases can be computed as above for each factor in that product. The resulting phases can then be multiplied together as they are all diagonal operators and therefore commute. This finalizes the argument that our construction carries over to any Abelian topological order derived from a TQD model.

D. Towards Quantum Error Correction
We have constructed pairwise commuting operatorsÃ g v that, together with B g p (see Eq. (4)) generate the Non-Pauli stabilizer group S T QD . Just as in other topological codes, a measurement of all generators of this group gives rise to a unique excitation pattern which can then be corrected by annihilating them with suitable string-like operators.
The fact thatÃ g v not only consists of X-like operators but also a diagonal operator means that they -in general -do not commute with an X error. Thus, the vertex stabilizers are sensitive to bit-flip as well as phase-flip errors, in contrast to the toric code where vertex operators detect solely the phaseflip errors. Moreover, a string of such bit-flip errors not only gives rise to an elementary excitation at the end points of the string (as in the toric code) but also along its length. In fact, the number of excitations along an X-string will scale with its size. This extra sensitivity is a manifestation of the fact that the our TQD-derived codes are not of CSS type, and may reduce the need for long recovery strings, possibly increasing the performance of standard decoder approaches based on minimum-weight perfect matching.

IV. CONCLUSION AND OUTLOOK
In this work, we have exploited the deep connections between topological phases of matter and topological error correction to construct a new class of stabilizer codes built from twisted quantum double models hosting Abelian anyons. To do so, we have established a systematic and quantitative understanding of how the vertex operators of twisted quantum double models fail to commute outside of the ground space and therefore precluding their use as stabilizers without further modification. We began with the relatively straightforward task of deriving the obstructions for the fixed-point Hamiltonian of the double semion phase -the twisted version of the toric code Hamiltonian -and of a twisted Z 2 ×Z 2 phase. By appropriately modifying the vertex operators, we have obtained commuting stabilizers from both models that, in principle, can be implemented with a two dimensional qubit architecture. This approach readily generalizes to other twisted quantum double models with higher local dimensions. We have explicitly constructed commuting stabilizers from the twisted Z N and Z 2 N models, making it possible to construct stabilizers for every Abelian TQD model.
Our findings invite further research into explicit error correction schemes based on these novel stabilizers so that their potential may be fully explored. For any decoder proposed, one has to find suitable the string operators for these modified models, which provide both the recovery operations and the logical operations on the code space. We expect that the computational power of these codes will differ significantly from that of previously studied ones, as the anyons of the topological order used to construct the code fix the algebra of the logical operators, the twist defects and domain walls [31], and thereby determining the fault-tolerant gate set. To the best of our knowledge, these have yet to be determined for TQD-derived codes, and would represent a first step in determining the performance of these when used in conjunction with known universal computation strategies such as magic state distillation or the recently-proposed just-in-time decoders [12,13]. Though these codes in and of themselves are a novel take on the idea of stabilizer-based error correction, it is our hope that the many questions raised by this class of codes spurs further investigation in to their properties.
In this section, we will give an algebraic definition of group cohomology with U (1) coefficients. In some sense it can be thought of a condensed version of the appendix in Ref. [22] in which we give all the necessary background to understand our general framework described in the next appendices. Besides that, it is an interesting subject on its own and pops up in different fields of physics and mathematics.
To evaluate the obstruction phases Ω and Π that we have defined in App. B, we first need to choose a (non-trivial) cocycle representative [ω] ∈ H 3 (Z N , U (1)) to insert into Eqs. (B4) and (B10). For Z N , there are N cocycle classes that can be labelled by p ∈ Z N are represented by Such cocycles, depending only on elements from the same group Z N , are called type-I cocycles [30]. All representatives are generated by ω 1 and (ω 1 ) 0 = (ω 1 ) N ≡ 1 represents the trivial cocycle class. Investigating the TQD model defined by the cocycle ω 1 and lifting the obstructions for that particular model with a generating modification phase is therefore enough to lift the obstructions for all Z N TQD models. The modification phase that lifts the obstruction for a model defined with an ω p cocycle is given as pth power of the generating modification phase.
Inserting = Ω (l5) where B 3,9,8,1 = (−l 3 ) ⊕ (−l 9 ) ⊕ l 8 ⊕ l 1 and ∆ g+h,N = 1 for g + h ≥ N and 0 otherwise. Interestingly, the l 2 dependence drops out such that, for a fixed flux configuration, the modification phase is only supported on l 5 . Similarly, inserting ω 1 into Eq. (B10) gives the commutativity obstruction Note that these obstructions coincide with the ones calculated for N = 2, where g = h = 1 is the only non-trivial case, in Sec. III A.
Having quantified the obstructions for the Z N TQD model, we introduce the modified vertex operatorsÃ g v = η g v A g v with the phase modification η g v being a (in the edge basis) diagonal operator with all entries taking values being roots of unity. We impose that η g v flux-free = 1 so that the ground space properties -and with them the topological data of the model -remain unchanged. Let us first consider two vertex operators acting on the same vertex. The modified operators should fulfill Inserting the definition ofÃ g v and Eq. (B4) yields the on-site consistency condition on the phases {η g v }, Since η g v is diagonal in the edge basis, conjugation with A g v modifies it only by shuffling edge values. In particular, Using the cyclic property of Z N , we can set g = 1 and h = N − 1 in Eq. (C8) to obtain an equation for η 0 . Since A 0 v = 1, η 0 v = 1, and we obtain We expand this further by recursively writing , . . . , η 2 v , thereby obtaining the closure relation where we have used that a diagonal operator (such as Ω and η) conjugated by A g v is still a diagonal operator and therefore commutes with any other diagonal operator. Eq. (C10) makes is possible to find a solution for the phase corresponding to the generator of Z N , η 1 v . Moreover, the recursion process that led us to Eq. (C10) can be used iteratively to generate all other phases η g v ∀g ∈ Z N so that Eq. (C7) is fulfilled. Once that is achieved, the general commutativity problem reduces to restoring commutativity to the generator {Ã 1 v } since any other modified vertex operator decomposes asÃ g v = (Ã 1 v ) g . We will therefore first solve the on-site consistency condition and then derive a second consistency condition for the commutativity of the modified vertex operators.
Inserting the explicit form of the obstruction phase Ω The double sum in the exponent can be simplified with some projector algebra. It reads where we have used that Since the vertex operators do not change fluxes, this equation can be solved by any expression of the form 9,8,1+B6,12,11,4) We are left with a freedomη 1 v to solve an additional consistency equation coming from the commutativity obstruction phase (C6). To be precise, imposing that the generating vertex operators Substituting in Eqs. (C14) and (C6), we obtain the two consistency conditions to lift the on-site and commutativity obstructions. We find a one-parameter family of solutions where we have used a similar manipulation as in Eq. (C12) to show the second condition forη 1 v . Note that η 1 v (p + N 2 ) = η 1 v (p), and thus all distinct solutions in this family are labeled by p ∈ [0, N 2 ). Putting these together with Eq. (C14) yields the full modification phase 9,8,1+pB6,12,11,4) This expression consists of two parts. The rightmost one, built from two sums of projectors and having an explicit l 2 and l 5 dependence ensures commutativity, whereas the other part, depending only on fluxes and therefore not altering the commutativity properties, is there to fulfill the on-site condition. Using the recurrence relation from the derivation of the closure relation Eq. (C8), we generate all other modification phases {η g v } iteratively so that the on-site condition is fulfilled. The general modification phase then reads (C19) With that, we have obtained a solution for a very general case, namely all Z N TQD models. On a first glance this expression seems complex, but a little inspections shows that, for specific choices of p and certain (small) local dimensions, it reduces to a manageable expression (see Sec. III A). In this section, we will calculate the obstruction phases derived in App. B for G = Z 2 N and a type-II cocycle. Each element in the input group can be written as a pair of Z N elements, g = (g 1 , g 2 ) with group multiplication naturally carrying over from Z N . 8 It is clear that this group is generated by two elements, namely (0, 1) and (1,0). When looking at the cohomology of this group, one finds that the resulting cocycle classes are generated by three elements split into two types [30]. The two type-I generators are the same as those for Z N , depending on the data from a single tensor factor, i.e.
where ω p was defined in Eq. (C4). Using such a cocycle to define a TQD model will result in the same functional form of the obstructions Ω and Π found for Z N . Hence, they also can be removed by the same modification phase η g v from Eq. (C19) where all Z N variables now carry a tensor factor index i.
Besides these type-I cocycles, there are type-II cocyles that depend on both tensor factors simultaneously and can be represented by where λ = e 2πi/N . One could also define the cocycle with indices 1 and 2 interchanged, but this is known to be gauge equivalent to the above definition [30]. Eq. (D2) shows that cocycles of type II mixes the two tensor factors of the input group elements in a non-trivial way. Whereas the Z 2 N TQD model with a type-I cocycle can be decomposed into two (possibly inequivalent) Z N TQD models, a type-II cocyle gives rise to a different topological order that cannot be factored in this way. In the following, we calculate the obstructions Ω and Π with a type-II cocycle and investigate how to lift these with appropriately chosen phase modifications.
Inserting the chosen type-II cocycle representative ω II,1 from Eq. (D2) into the obstruction phases calculated in Eqs. (B4) and (B10), we find that the vertex operators of the type-II Z 2 N TQD model fail to implement the group action faithfully on site, generating the obstruction phase where the fluxes and projectors are defined as in the previous section but with every group element and edge value carrying an additional (upper) index (i) labelling the tensor factor of the corresponding variable. To avoid notation clutter in the projectors, we only write the tensor factor index once. For example, P 5 ) = n. The quantity ∆ •,• is defined as in Eq. (C5). Note that this obstruction phase consists of similar terms as those found for Z N case, but with added mixing between tensor factors. The analoguous calculation of the commutativity obstruction phase yields The procedure to lift these obstructions begins identically to that in the previous section. Again, we derive closure relations from group multiplication in every cyclic sub-factor of Z 2 N . In addition to these, one also finds extra constraint equations coming from group multiplication of elements from different sub-factors when the model includes a non-trivial type-II cocycle. Just as in the previous section, once the group multiplication is implemented consistently on-site, we only need to make the vertex operators for the generators (0, 1) and (1, 0) commute and can then iteratively construct all other modification phases so that both the on-site and commutativity obstructions are removed.
We start by introducing the modified vertex operatorsÃ g v = η g v A g v ∀g ∈ Z 2 N and imposing Eq. (C7). In particular, this directly implies a general condition on the modification phases, Eq. (C8). This enables us to use cyclicity to iteratively derive the three consistency conditions from the relations (N − 1, 0 Given the large freedom available when fulfilling the first two closure relations, we can constructη (1,0) such that it cancels out λ −B (2) 6,12,11,4 andη (0,1) such that it cancels out λ −B (1) 6,12,11,4 . One solution is given bȳ Interestingly, the phase for the first generator only depends on the second tensor factor of fluxes and edges and vice versa. This is a direct consequence of the way the type-II cocycle couples the two tensor factors. The other modification phases can be calculated iteratively using Eq. (C8) to produce a proper representation of the group action. We will give an explicit expression the modification phases for any group element g ∈ Z 2 N at the end of this section after having discussed the commutativity obstruction. Once the vertex operators faithfully represent the group action on-site, every modified vertex operator can be decomposed in terms of the generating vertex operatorsÃ (0,1) v andÃ (1,0) v . Therefore, it is sufficient to resolve the commutativity obstruction for those two operators while still fulfilling Eqs. (D5). Imposing [Ã g 2 ,Ã h v ] G = 1 ∀ g, h ∈ {(0, 1); (1, 0)} and evaluating the obstruction phases Π Surprisingly, the generating modification phases derived from the on-site condition (see Eqs. (D8)) also fulfill these four equations. As mentioned, Eq. (C8) now allows us to iteratively generate all modification phases for any group element g = (g 1 , g 2 ).