# Non-Pauli topological stabilizer codes from twisted quantum doubles

Julio Carlos Magdalena de la Fuente<EMAIL_ADDRESS>JARA Institute for
Quantum Information, RWTH Aachen University, Aachen, Germany  Dahlem Center
for Complex Quantum Systems, Freie Universität Berlin, 14195 Berlin, Germany
Nicolas Tarantino  Dahlem Center for Complex Quantum Systems, Freie
Universität Berlin, 14195 Berlin, Germany  Jens Eisert  Dahlem Center for
Complex Quantum Systems, Freie Universität Berlin, 14195 Berlin, Germany

###### Abstract

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.

##  1 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 [25, 30, 8]. Since the no-
cloning theorem [28] prevents qubits from being copied, quantum error
correction cannot rely on simply copying the necessary information 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 [16, 31] 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 [16,
14, 31]. The _toric code_ [20], its associated planar embedding known as the
_surface code_ [15], and _color codes_ [3], 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_ [28] – and a lack
of self-correctability [6] will be common to any stabilizer approach in two
spatial dimensions, several techniques have already been identified to
circumvent these limitations, including _magic state distillation_ [26] and
just-in-time decoders [7, 5]. Moreover, codes 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$ [9] or
enhanced bit flip stability [2]. Other generalisations involving non-commuting
stabiliser sets [27] 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.

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 [21], 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 well-studied instances of this
are equivalent to multiple copies of the toric code phase [4]. 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 [24, 17, 11]. This omission is also significant given the
fact that, from a technological perspective [1], the identification of new
topological codes seems imperative.

In this work, we use a large class of topological orders hosting 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 Non-Pauli 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. [13], where the _double semion string-net model_ [24] 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. The main limitation of our
scheme is that it is not applicable to TQD models with Non-Abelian anyons.

This work is structured as follows: In Section 2 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
$\mathbb{Z}_{2}$ and a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ model and then
summarized for the general $\mathbb{Z}_{N}$ and $\mathbb{Z}_{N}^{2}$ cases in
Sec. 3, 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. 4, 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.

##  2 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 [20] introduced by Kitaev. They can be obtained by promoting the
global symmetry of a _symmetry protected topological (SPT) phase_ [11, 23, 10]
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. [17]
for more comprehensive perspective.

Figure 1: Local patch of a translation-invariant lattice on which we define
our model. It is a oriented triangulation of a compact surface. Around each
vertex, the edges directly adjacent to it are labelled from $l_{1}$ to $l_{6}$
and the other 6 edges that share a triangle with the vertex are labelled by
$l_{7}$ to $l_{12}$. Together, these edges constitute the _neighborhood_ of
the vertex and are marked in red above.

###  2.1 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
$\mathcal{H}_{l}$ is spanned by states labeled by elements of a finite group
$G$,

$\displaystyle\mathcal{H}_{l}=\operatorname{span}_{\mathbb{C}}\\{\ket{g},\;g\in
G\\}$ (1)

with $\innerproduct{g}{h}=\delta_{g,h}$. 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
$\mathcal{H}=\bigotimes_{\text{edges }l}\mathcal{H}_{l}$. While $G$ can be
chosen to be any finite group, we will only be treating Abelian cases in this
work. Additionally, the TQD model takes a cocycle function $\omega:G^{2}\to
U(1)$ as an input.111To be precise, the cocycle $\omega$ has to be normalized,
i.e. $\omega(0_{G},a,b)=\omega(a,0_{G},b)=\omega(a,b,0_{G})=1\;\forall a,b\in
G$. It will define the action of the vertex operators in the Hamiltonian.

The full TQD Hamiltonian is given by

$\displaystyle H_{TQD}=-\sum_{\text{plaquettes }p}B_{p}-\sum_{\text{vertices
}v}A_{v}.$ (2)

The first sum runs over all (triangular) plaquettes and the plaquette operator
acting on a triangular face is defined by

$\displaystyle\begin{split}B_{p}&\Bigg{|}\begin{minipage}[c]{36.98866pt}
\includegraphics[width=433.62pt]{triangle_labeled}
\end{minipage}\Bigg{\rangle}\\\ &=\delta_{(-l_{3})\oplus l_{2}\oplus
l_{1}}\Bigg{|}\begin{minipage}[c]{38.41139pt}
\includegraphics{triangle_labeled} \end{minipage}\Bigg{\rangle}\\\
&=\frac{1}{\absolutevalue{G}}\sum_{i\in C}\underbrace{\chi_{i}(-l_{3}\oplus
l_{2}\oplus l_{1})}_{\eqqcolon
B_{p}^{i}}\Bigg{|}\begin{minipage}[c]{36.98866pt}
\includegraphics[width=433.62pt]{triangle_labeled}
\end{minipage}\Bigg{\rangle},\end{split}$ (3)

where $\oplus$ denotes the group operation (modular addition for cyclic
groups) and $-l$ the inverse element of $l$ in $G$. The projector
$\delta_{g}=1$ for $g=0_{G}$ (identity element in $G$) and $\delta_{g}=0$
otherwise. The plaqeutte terms ensure that the ground space is flux free,
i.e., $B_{p}=1\;\forall p$. In the second step, we decomposed the projector
into a sum of group characters $\chi_{i}$ over conjugacy classes $C$, which
can be done for any finite group. While this decomposition is usually a
curiosity, this is useful when $G$ is an Abelian group. In this case, each
conjugacy class contains only a single element, and so the characters are in
one-to-one correspondence with the group elements. For cyclic groups, where
$g\in\\{0,\dots,N-1\\}$, this means that these character operators can be
thought of as the products _phase operators_ $Z_{l}$,

$\displaystyle B_{p}^{g}=$ $\displaystyle\prod_{i\sim
p}\sum_{n=0}^{N-1}\left(e^{\frac{2\pi
i}{N}n}\right)^{s(p,l)g}\outerproduct{n}{n}$ (4a) $\displaystyle=$
$\displaystyle\prod_{l\sim p}Z_{l}^{s(p,l)g},$ (4b)

with $s(p,i)=+1(-1)$ if edge $l$ is oriented in counterclockwise(clockwise)
direction around plaquette $p$. For $G=\mathbb{Z}_{2}$ for example, $Z_{i}$ is
the Pauli $Z$ operator acting on edge $i$ and one recovers the toric code
plaquette operator. Arbitrary finite abelian groups can always be decomposed
into cyclic factors, and thus their character operators can be factored into
products of phase operators as well.

The second sum in the Hamiltonian runs over all vertices and the summand
$A_{v}$, acting on the neighborhood of a vertex $v$, is defined by a similar
decomposition as $B_{p}$, i.e.,

$\displaystyle A_{v}=\frac{1}{\absolutevalue{G}}\sum_{g\in G}A_{v}^{g},$ (5)

where $\absolutevalue{G}$ is the order of the group $G$, and $A_{v}^{g}$ the
_vertex operator_ corresponding to the group element $g$ defined by its action
on a basis element,

$\displaystyle A_{v}^{g}\Bigg{|}\begin{minipage}[c]{62.59596pt}
\includegraphics[width=433.62pt]{basis_config}
\end{minipage}\dots\Bigg{\rangle}=$
$\displaystyle\frac{\omega(l_{9},l_{3}^{\prime},g)\omega(l_{3}^{\prime},g,l_{4})\omega(g,l_{4},l_{11})}{\omega(l_{8},l_{1}^{\prime},g)\omega(l_{1}^{\prime},g,l_{6})\omega(g,l_{6},l_{12})}\Bigg{|}\begin{minipage}[c]{62.59596pt}
\includegraphics[width=433.62pt]{basis_config_prime}
\end{minipage}\dots\Bigg{\rangle},$ (6)

where the label of the central vertex $v$ is left out for readability. The
vertex operator can be decomposed into two consecutive actions. First, it
changes the values of the edges adjacent to the vertex depending on their
orientation. On our lattice,

$\displaystyle l_{i}\stackrel{{\scriptstyle
A_{v}^{g}}}{{\longmapsto}}l_{i}^{\prime}=l_{i}\oplus(-g)$
$\displaystyle\text{for }i=1,2,3,$ (7a) $\displaystyle
l_{j}\stackrel{{\scriptstyle A_{v}^{g}}}{{\longmapsto}}l_{j}^{\prime}=g\oplus
l_{j}$ $\displaystyle\text{for }j=4,5,6,$ (7b)

where $\oplus$ denotes the group operation which we choose to be modular
addition for cyclic groups. Second, it scales the wavefunction by a phase
factor given by the product of $6$ cocycles $\omega:G^{3}\to U(1)$, 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 pre-factor is
determined by the orientation structure of the lattice. For a detailed
explanation of constructing the pre-factor for a general lattice, see Ref.
[17]. The cocycles encode the _topological data_ of the theory modelled by
$H_{TQD}$. Their defining property is the so-called cocycle condition

$\displaystyle\begin{split}\frac{\omega(g_{1},g_{2},g_{3})\omega(g_{0},g_{1}\oplus
g_{2},g_{3})\omega(g_{0},g_{1},g_{2})}{\omega(g_{0}\oplus
g_{1},g_{2},g_{3})\omega(g_{0},g_{1},g_{2}\oplus g_{3})}=1\\\\[6.0pt] \forall
g_{0},g_{1},g_{2},g_{3}\in G.\end{split}$ (8)

Obviously, $\omega(a,b,c)\equiv 1$ is always a solution and is called trivial.
If we use this trivial solution in Eq. (6) to define $H_{TQD}$, we obtain the
_quantum double Hamiltonian_ from [19]. Since – in general – there are non-
trivial solutions to this equation as well, 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 topologically ordered Hamiltonian.
However, not all solutions yield distinct orders but they are classified in
equivalence classes. A close investigation of the cocycle condition reveals
that if we have one solution $\omega$, we can always obtain another solution

$\displaystyle\tilde{\omega}(g_{1},g_{2},g_{3})=\omega(g_{1},g_{2},g_{3})\frac{\beta(g_{2},g_{3})\beta(g_{1},g_{2}\oplus
g_{3})}{\beta(g_{1}\oplus g_{2},g_{3})\beta(g_{1},g_{2})},$ (9)

where $\beta:G^{2}\to U(1)$ is an arbitrary function mapping two group
elements to a phase factor. If we have two TQD Hamiltonians defined by two
cocycles $\omega_{1}$ and $\omega_{2}$ in Eq. (6) so that they are in
different topological orders, we know and there exists no $\beta$ to map
$\omega_{1}$ onto $\omega_{2}$ by Eq. (9). Hence, inequivalent Hamiltonians
$H_{TQD}$ (in the sense of topological order) are classified by distinct
equivalence classes of functions $\omega$, which define elements of the third
group cohomology of $G$ over $U(1)$ with $\omega\sim\tilde{\omega}$ iff
$\exists\beta:G^{2}\to 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
$\mathbb{Z}_{2}$ and $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$. For an introduction
into group cohomology, see App. A.

###  2.2 Topological data

The Hamiltonian we have constructed on a triangulation of a compact surface
from a group $G$ together with a cocycle $\omega$ is indeed topologically
ordered. It has anyonic excitations and a robust ground state degeneracy
(GSD). For models with Abelian anyons only,
$\text{GSD}=\absolutevalue{G}^{2g}$ where $g$ is the genus of the surface on
which it is defined. 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.
[17].

###  2.3 Ground space and failure of commutativity

The ground space of $H_{TQD}$ can be found exactly and is defined implicitly
by the conditions

$\displaystyle\begin{split}B_{p}\ket{\psi}=A_{v}\ket{\psi}=\ket{\psi}\quad\forall&\text{
plaquettes }p\\\ &\text{ and vertices }v,\end{split}$ (10)

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$).222For a precise statement on how commutativity fails for
general Abelian groups, we refer to App. B. In the TQD model for
$G=\mathbb{Z}_{2}$ and the non-trivial cocycle [17]

$\displaystyle\omega_{1}(1,1,1)=-1,\mbox{\quad\quad}\omega(a,b,c)=1\quad\text{else,}$
(11)

a version of the double semion phase, the vertex operators do not commute when
acting on the following configuration:

$\displaystyle A_{2}^{1}A_{v}^{1}\ket{\begin{minipage}[c]{42.67912pt}
\includegraphics[width=433.62pt]{flux_config1}
\end{minipage}}=\ket{\begin{minipage}[c]{42.67912pt}
\includegraphics[width=433.62pt]{flux_config3} \end{minipage}}$ $\displaystyle
A_{v}^{1}A_{2}^{1}\ket{\begin{minipage}[c]{42.67912pt}
\includegraphics[width=433.62pt]{flux_config1}
\end{minipage}}=-\ket{\begin{minipage}[c]{42.67912pt}
\includegraphics[width=433.62pt]{flux_config3} \end{minipage}}.$

In the case of $G=\mathbb{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, $\ket{0}=\circ$
and $\ket{1}=\bullet$. The rest of the lattice is not explicitly shown, but we
assume that all other qubits are in $\ket{0}$. The vertex operators will be
defined in Sec. 3.1 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}=\prod_{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.

##  3 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,

$\displaystyle\tilde{A}_{v}^{g}=D_{v}^{g}A_{v}^{g},$ (12)

where $D_{v}^{g}$ is a full-rank operator that is diagonal in the edge basis,
with entries being roots of unity, and is equal to $\mathds{1}$ on the ground
space.

By imposing that the modified operators $\\{\tilde{A}_{v}^{g}\\}$ commute, we
obtain consistency equations for the modification phases $\\{\eta_{v}^{g}\\}$.
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 $\mathbb{Z}_{2}$ and
$\mathbb{Z}_{2}\times\mathbb{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
$\mathbb{Z}_{N}$ and $\mathbb{Z}_{N}^{2}$. 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 $\tilde{A}_{v}^{g}$ together with the plaquette
phase operators $B_{p}^{g}$ (see Eq. (3)) generate the Non-Pauli stabilizer
group

$\mathcal{S}_{TQD}=\langle\\{B_{p}^{g}\\},\\{\tilde{A}_{v}^{g}\\}\rangle,$
(13)

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_{p}^{g}$ and
$\tilde{A}_{v}^{g}$.

In the following, we will illustrate the construction for an arbitrary finite
group $G$ and a cocycle $\omega$ 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.

###  3.1 $\mathbb{Z}_{2}$ – double semion code

We first investigate the non-trivial $\mathbb{Z}_{2}$ model that is in the
same phase as the double semion string-net model [24]. We represent
$\mathbb{Z}_{2}$ as the set $\\{0,1\\}$ together with the group operation
being addition modulo two. With that,
$\mathcal{H}_{l}=\operatorname{span}_{\mathbb{C}}\\{\ket{0},\ket{1}\\}=\mathbb{C}^{2}$,
so it is a model of interacting qubits. $\mathbb{Z}_{2}$ has two inequivalent
cocycle classes, one trivial class $[\omega_{0}\equiv 1]$ and one non-trivial
class $[\omega_{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. (11). Inserting this cocycle into Eq. (6) yields
the Hamiltonian

$\displaystyle H_{\mathbb{Z}_{2}}=$ $\displaystyle-\sum_{\text{plaquettes
}p}\frac{1}{2}\Big{(}\mathds{1}+\prod_{l\sim
p}Z_{l}\Big{)}-\sum_{\text{vertices
}v}\frac{1}{2}\Big{(}\mathds{1}+A_{v}^{1}\Big{)}$ (14a)
$\displaystyle\mbox{\quad
with\quad}A_{v}^{1}=(-1)^{P_{9}^{-}P_{3}^{-}+P_{3}^{-}P_{4}^{+}+P_{4}^{+}P_{11}^{-}+P_{8}^{-}P_{1}^{-}+P_{1}^{-}P_{6}^{+}+P_{6}^{+}P_{12}^{-}}\prod_{l=1}^{6}X_{l},$
(14b)

where $Z_{l}$ and $X_{l}$ are the Pauli $Z$ and $X$ matrices on (the qubit
sitting on) edge $l$ and $P_{l}^{\pm}=\frac{1}{2}(1\pm Z_{l})$ is the
projector onto the space where edge $l$ caries the value $0$ or $1$,
respectively. Each term of the form $(-1)^{P^{-}_{i}P^{-}_{j}}$ corresponds to
a $CZ$ gate between qubit $i$ and $j$. These entangling gates are essential to
build up long range entanglement that is inequivalent to the one in the toric
code. The above Hamiltonian is in the so-called _double semion phase_.

By construction, $A_{v}^{1}$ always flips an even number of qubits adjacent to
a plaquette, and thus $\commutator{A_{v}^{1}}{B_{p}}=0\;\forall v,p$.333In
fact, $A_{v}^{g}$ is flux preserving in any twisted quantum double model due
to the orientation structure of the edges. Hence, the only obstruction for the
operators $B_{p}$ and $A_{v}^{1}$ to form a commuting set of operators that we
can use for stabilizer error correction comes from the vertex operators
$A_{v}^{1}$.

####  3.1.1 Obstruction in the original model

For the operators $A_{v}^{1}$ to generate a stabilizer group they have to
represent the group action of $\mathbb{Z}_{2}$ on site. In particular, any
representative of an element in $\mathbb{Z}_{2}$ should square to the
identity. Unfortunately, it turns out that

$\displaystyle\begin{split}\left(A_{v}^{1}\right)^{2}=&(-1)^{P_{1}^{-}+P_{3}^{-}+P_{4}^{-}+P_{6}^{-}+P_{8}^{-}+P_{9}^{-}+P_{11}^{-}+P_{12}^{-}}\mathds{1}\\\
\neq&\mathds{1},\end{split}$ (15)

where we have used $X_{l}^{2}=\mathds{1}$, and the decomposition of the
identity, $\mathds{1}=P_{l}^{+}+P_{l}^{-}$ for any edge $l$. In fact,
$\left(A_{v}^{1}\right)^{2}=-\mathds{1}$ exactly when an odd number of the
edges $\\{1,3,4,6,8,9,11,12\\}$ is in the state $\ket{1}$, which coincides
with

$\displaystyle\left(A_{v}^{1}\right)^{2}=(-1)^{B_{3,9,8,1}+B_{6,12,11,4}}\mathds{1},$
(16)

where $B_{i,j,k,l}=\frac{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_{v}^{g}$ must commute pairwise. Due to
the translational invariance of the model, we need only calculate the
commutation relation between $A_{v}^{1}$ with the three operators
$\\{A_{i}^{1},\;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

$\displaystyle A_{v_{1}}^{1}A_{v}^{1}=$ $\displaystyle
A_{v}^{1}A_{v_{1}}^{1},$ (17a) $\displaystyle A_{v_{3}}^{1}A_{v}^{1}=$
$\displaystyle A_{v}^{1}A_{v_{3}}^{1},$ (17b)
$\displaystyle\begin{split}A_{v_{2}}^{1}A_{v}^{1}=&(-1)^{P_{1}^{-}+P_{3}^{-}+P_{8}^{-}+P_{9}^{-}}A_{v}^{1}A_{v_{2}}^{1}\\\
=&(-1)^{B_{3,9,8,1}}A_{v}^{1}A_{v_{2}}^{1},\end{split}$ (17c)

using the same relations used to produce Eq. (15) and (16). 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 orientation444Choosing a different edge orientation in the first
place would only shift the commutativity obstruction to different edges, not
remove it. which determines which arguments enter in the cocycles in Eq. (6).

####  3.1.2 Modifying vertex operators by local phase

We have found that the vertex operators in the original TQD model fail to be
commuting stabilizers on the whole Hilbert space because, on one hand, they do
not implement the group action on site, i.e. $(A_{v}^{1})^{2}\neq\mathds{1}$,
and, on the other, fail to commute $A_{v_{1}}^{1}A_{2}^{1}\neq
A_{2}^{1}A_{v_{1}}^{1}$. 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 diagonal unitary as described in the beginning of Sec. 3,

$\displaystyle\tilde{A}_{v}^{1}=D_{v}A_{v}^{1}$ (18)

so that $\tilde{A}_{v}^{1}$ are stabilizers. The modified operators should
therefore satisfy

$\displaystyle\left(\tilde{A}_{v}^{1}\right)^{2}=$ $\displaystyle\mathds{1},$
(19a) $\displaystyle\tilde{A}_{v_{1}}^{1}\tilde{A}_{v_{2}}^{1}=$
$\displaystyle\tilde{A}_{v_{2}}^{1}\tilde{A}_{v_{1}}^{1},$ (19b)

for all vertices $v,v_{1},v_{2}$. Imposing Eq. (19a) and using Eq. (16) yields
the condition on the modification phase

$\displaystyle
1=D_{v}A_{v}^{1}D_{v}\left(A_{v}^{1}\right)^{-1}(-1)^{B_{3,9,8,1}+B_{6,12,11,4}}.$
(20)

In fact, this can be solved by the Ansatz

$\displaystyle D_{v}=i^{B_{3,9,8,1}+B_{6,12,11,4}}\bar{D}_{v}$ (21)

with $\bar{D}_{v}A_{v}^{1}\bar{D}_{v}\left(A_{v}^{1}\right)^{-1}=1$ since
$A_{v}^{1}$ leaves fluxes invariant. Now, inserting this family of solutions
into Eq. (19b) and using Eqs. (17) yields a second condition on $\bar{D}_{v}$,

$\displaystyle\begin{split}1=&(-1)^{B_{3,9,8,1}}\bar{D}_{2}A_{v}^{1}(\bar{D}_{2})^{-1}\left(A_{v}^{1}\right)^{-1}\\\
\times&(\bar{D}_{v})^{-1}A_{2}^{1}\bar{D}_{v}\left(A_{2}^{1}\right)^{-1},\end{split}$
(22)

which is solved by the one-parameter family

$\displaystyle\begin{split}\bar{D}_{v}(p)=&e^{\pi
i[pB_{3,9,8,1}L_{2}+(p-1)B_{6,12,11,4}L_{5}]}\\\ \times&e^{-\frac{\pi
i}{2}[(p+2)B_{3,9,8,1}+(p+1)B_{6,12,11,4}]},\end{split}$ (23)

with $p\in\mathbb{R}$. $L_{i}=\frac{1}{2}(1-Z_{i})$ is the operator that
measures the value of the edge $l_{i}$ in $\mathbb{Z}_{2}$. Note that
$\bar{D}$ is a periodic function in $p$ with periodicity $4$, i.e.
$\bar{D}(p+4)=\bar{D}(p)$. Also, the second factor containing $i$ only depends
on fluxes and ensures that
$\bar{D}_{v}A_{v}^{1}\bar{D}_{v}\left(A_{v}^{1}\right)^{-1}=1$ and does not
affect commutativity.

Figure 2: Geometric illustration of the support of $D_{v}^{1}(p)$ in the
neighborhood of vertex $v$. For any $\mathbb{Z}_{N}$ model (see App. C), the
obstructions (and therefore also the modification phases) only depend on the
horizontal edges, labelled by $l_{2}$ and $l_{5}$ and the two fluxes through
the regions around them, $B_{3,9,8,1}$ and $B_{6,12,11,4}$, depicted by a
dashed area in the figure above. The edges $l_{1},l_{8},\dots$ that are
written in gray only enter implicitly via the fluxes that are composed of
them. Note the orientation structure can be neglected in the $\mathbb{Z}_{2}$
case.

Combined, we obtain the modification phase

$\displaystyle\begin{split}D_{v}(p)=&e^{-\frac{\pi
i}{2}[(p+1)B_{3,9,8,1}+pB_{6,12,11,4})]}\\\ \times&e^{\pi
i[pB_{3,9,8,1}L_{2}+(p-1)B_{6,12,11,4}L_{5}]},\end{split}$ (24)

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.
$D_{v}(p+4)=D_{v}(p)$, $p\in[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\in[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

$\displaystyle D_{v}(1)=$
$\displaystyle(-i)^{B_{6,12,11,4}}(-1)^{B_{3,9,8,1}}(-1)^{B_{3,9,8,1}L_{2}},$
(25)

so that it does not depend on $l_{5}$.

Explicitly quantifying the obstructions of the group property (Eq. (16)) and
the commutativity of the operators acting on neighboring vertices (Eqs. (17)
in the original TQD model enabled us to remove them with a local phase
modification such that the modified operators
$\tilde{A}_{v}^{1}=D_{v}A_{v}^{1}$ faithfully represent the group
$\mathbb{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_{v}^{1}$ on the ground space, and thus the modified Hamiltonian
it still is in the double semion phase.

###  3.2 $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ – twisted color codes

In the previous section, we have 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=\mathbb{Z}_{2}\times\mathbb{Z}_{2}=\\{g=(g_{1},g_{2});\;g_{1},g_{2}\in\mathbb{Z}_{2}\\}$
as an input group for the TQD model, the local Hilbert space becomes
$\mathcal{H}_{l}=\operatorname{span}_{\mathbb{C}}\\{\ket{0,0},\ket{0,1},\ket{1,0},\ket{1,1}\\}\simeq\mathbb{C}^{2}\otimes\mathbb{C}^{2}$
which can be realized using two qubits per edge. As was shown in Ref. [22], an
untwisted $\mathbb{Z}_{2}\times\mathbb{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=\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ are classified by
$H^{3}(\mathbb{Z}_{2}\times\mathbb{Z}_{2},U(1))=\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}$,
and thus there are 8 different cocycle classes we can choose as input,
labelled by $(s_{1},s_{2},s_{3})\in\mathbb{Z}_{2}^{3}$. In an appropriate
gauge, a general cocycle $\omega\in
H^{3}(\mathbb{Z}_{2}\times\mathbb{Z}_{2},U(1))$ can be written as

$\displaystyle\begin{split}\omega(a,b,c)=&\omega_{1}^{s_{1}}(a_{1},b_{1},c_{1})\omega_{1}^{s_{2}}(a_{2},b_{2},c_{2})\\\
&\times\omega_{II}^{s_{3}}(a,b,c),\end{split}$ (26)

with $s_{i}=0,1$ and the group elements are represented by pairs of
$\mathbb{Z}_{2}$ variables, i.e., $a=(a_{1},a_{2})$. When $s_{3}=0$, only the
cocycles $\omega_{1}$ appear. They are the same as those seen for the
$\mathbb{Z}_{2}$ phases, defined in Eq. (11), 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 $\mathbb{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

$\displaystyle\omega_{II}(a,b,c)=$
$\displaystyle\omega_{1}(a_{1},b_{2},c_{2})$ (27a)
$\displaystyle\begin{split}=&\begin{cases}-1&a_{1}=b_{2}=c_{2}=1\\\
1&\text{else}\end{cases},\end{split}$ (27b)

which mixes the two tensor factors and therefore is unique to the
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ case.555Analogously, we could represent
this cocycle by $\omega_{1}(a_{2},b_{1},c_{1})$. However, it is gauge-
equivalent to the one we are using [29]. 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. (27) into
Eq. (6) reads

$\displaystyle\begin{split}H_{\mathbb{Z}_{2}\times\mathbb{Z}_{2}}=&-\sum_{\text{plaquettes
}p}\frac{1}{4}\Big{(}\mathds{1}+\prod_{l\sim
p}Z_{l}^{(1)}\Big{)}\Big{(}\mathds{1}+\prod_{l\sim p}Z_{l}^{(2)}\Big{)}\\\
&-\sum_{\text{vertices
}v}\frac{1}{4}\Big{(}\mathds{1}+A_{v}^{(1,0)}+A_{v}^{(0,1)}+A_{v}^{(1,1)}\Big{)}\end{split}$
(28a) $\displaystyle\mbox{\quad with\quad}A_{v}^{(1,0)}=$
$\displaystyle(-1)^{P_{4^{2}}^{-}P_{11^{2}}^{-}+P_{6^{2}}^{-}P_{12^{2}}^{-}}\prod_{l=1}^{6}X^{(1)}_{l},$
(28b) $\displaystyle A_{v}^{(0,1)}=$
$\displaystyle(-1)^{P_{9^{1}}^{-}P_{3^{2}}^{-}+P_{3^{1}}^{-}P_{4^{2}}^{+}+P_{8^{1}}^{-}P_{1^{2}}^{-}+P_{1^{1}}^{-}P_{6^{2}}^{+}}\prod_{l=1}^{6}X^{(2)}_{l}\mbox{\quad
and\quad}$ (28c) $\displaystyle A_{v}^{(1,1)}=$
$\displaystyle(-1)^{P_{9^{1}}^{-}P_{3^{2}}^{-}+P_{3^{1}}^{-}P_{4^{2}}^{+}+P_{4^{2}}^{+}P_{11^{2}}^{-}+P_{8^{1}}^{-}P_{1^{2}}^{-}+P_{1^{1}}^{-}P_{6^{2}}^{+}+P_{6^{2}}^{+}P_{12^{2}}^{-}}\prod_{l=1}^{6}X^{(1)}_{l}X^{(2)}_{l},$
(28d)

where $X^{(i)}_{l},Z^{(i)}_{l}$ are the qubit Pauli matrices acting on the
$i$th tensor factor. And $P_{l^{i}}^{\pm}=\frac{1}{2}(1\pm Z_{l}^{(i)})$ are
the projectors on th values of the $i$th tensor factor of edge $l$.

Again,
$\commutator*{B_{p}^{(1)}}{A_{v}^{g}}=\commutator*{B_{p}^{(2)}}{A_{v}^{g}}=0\;\forall
g\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ since each vertex operator flips an
even number of qubits around each plaquette. The only obstructions preventing
$\\{B_{p}^{(i)},A_{v}^{g};\;i=1,2;g\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\}$
from forming a pairwise commuting set come from the vertex operators. We will
quantify the obstructions below.

####  3.2.1 Obstructions in the original model

Each element in $\mathbb{Z}_{2}\times\mathbb{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 $\mathds{1}$.
Since the representative we chose for the type-II cocycle in Eq. (27) does not
depend on $a_{2}$, $b_{1}$ and $c_{1}$, we find

$\displaystyle\left(A_{v}^{(1,0)}\right)^{2}=$
$\displaystyle\mathds{1}\mbox{\quad and\quad}$ (29a) $\displaystyle
A_{v}^{(0,1)}A_{v}^{(1,0)}=$ $\displaystyle A_{v}^{(1,1)}.$ (29b)

For the other products of non-trivial group elements $(0,1),(1,0)$ and $(1,1)$
however, we obtain explicit obstructions

$\displaystyle\begin{split}\left(A_{v}^{(0,1)}\right)^{2}=&(-1)^{B_{3,9,8,1}^{(1)}}\mathds{1}\neq\mathds{1},\end{split}$
(30a)
$\displaystyle\begin{split}\left(A_{v}^{(1,1)}\right)^{2}=&(-1)^{B_{3,8,9,1}^{(1)}+B_{6,12,11,4}^{(2)}}\mathds{1}\neq\mathds{1},\end{split}$
(30b)
$\displaystyle\begin{split}A_{v}^{(1,0)}A_{v}^{(0,1)}=&(-1)^{B_{6,12,11,4}^{(2)}}A_{v}^{(1,1)}\neq
A_{v}^{(1,1)}\end{split}$ (30c)

where
$B_{j,k,l,m}^{(i)}=\frac{1}{2}(\mathds{1}-Z_{j}^{(i)}Z_{k}^{(i)}Z_{l}^{(i)}Z_{m}^{(i)})$
measures the flux in the $i$th tensor factor through the region enclosed by
the edges $\\{j,k,l,m\\}$. As in the $\mathbb{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_{v}^{g}$ and the three operators
$\\{A_{i}^{h},\;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

$\displaystyle A_{v_{2}}^{(0,1)}A_{v}^{(1,0)}=$
$\displaystyle(-1)^{B_{3,9,8,1}^{(1)}}A_{v}^{(1,0)}A_{v_{2}}^{(0,1)},$ (31a)
$\displaystyle A_{v_{2}}^{(1,0)}A_{v}^{(0,1)}=$
$\displaystyle(-1)^{B_{3,9,8,1}^{(2)}}A_{v}^{(0,1)}A_{v_{2}}^{(1,0)},$ (31b)

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).

####  3.2.2 Modifying vertex operators by a local phase

We have found that in the original $\mathbb{Z}_{2}\times\mathbb{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.
$A_{v}^{g}A_{v}^{h}\neq A_{v}^{g\oplus h}$ 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 $\mathbb{Z}_{2}$ TQD model,
namely factors of $-1$ that only depend on fluxes. To resolve the obstructions
for the three operators $A_{v}^{(1,0)}$, $A_{v}^{(0,1)}$ and $A_{v}^{(1,1)}$,
we modify them by a local phase $D_{v}^{g}$ that is the identity on the ground
space (Eq. (12)). For the modified operators to be stabilizers, we need them
to fulfill

$\displaystyle\left(\tilde{A}_{v}^{g}\right)^{2}=$
$\displaystyle\mathds{1}\quad\forall g\in\mathbb{Z}_{2}\times\mathbb{Z}_{2},$
(32a) $\displaystyle\tilde{A}_{v}^{(1,0)}\tilde{A}_{v}^{(0,1)}=$
$\displaystyle\tilde{A}_{v}^{(0,1)}\tilde{A}_{v}^{(1,0)}=\tilde{A}_{v}^{(1,1)},$
(32b)
$\displaystyle\commutator{\tilde{A}_{v}^{g}}{\tilde{A}_{v^{\prime}}^{h}}=$
$\displaystyle 0,\quad v\neq v^{\prime},\forall g,h.$ (32c)

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 (32a) gives us independent constraints on each of the generating
phases $D_{v}^{(0,1)}$ and $D_{v}^{(1,0)}$. Using Eq. (30a), we find that the
first phase must satisfy

$\displaystyle\begin{split}1=&D_{v}^{(0,1)}A_{v}^{(0,1)}D_{v}^{(0,1)}\left(A_{v}^{(0,1)}\right)^{-1}\\\
&\times(-1)^{B_{3,9,8,1}^{(1)}},\end{split}$ (33)

which is solved by any solution of the form

$\displaystyle D_{v}^{(0,1)}=i^{B_{3,9,8,1}^{(1)}}\bar{D}_{v}^{(0,1)}$ (34)

together with
$\bar{D}_{v}^{(0,1)}A_{v}^{(0,1)}\bar{D}_{v}^{(0,1)}\big{(}A_{v}^{(0,1)}\big{)}^{-1}=1$.
Since $A_{v}^{(1,0)}$ squares to $\mathds{1}$ already, the condition on
$D_{v}^{(1,0)}$ is simpler. Eq. (32a) imposes that
$1=D_{v}^{(1,0)}A_{v}^{(1,0)}D_{v}^{(1,0)}\big{(}A_{v}^{(1,0)}\big{)}^{-1}$.
Now, using Eq. (32b) to express $D_{v}^{(1,1)}$ in terms of $D_{v}^{(1,0)}$
and $D_{v}^{(0,1)}$ and inserting it into the condition
$\big{(}\tilde{A}_{v}^{(1,1)}\big{)}^{2}=\mathds{1}$ yields an additional
constraint on $\bar{D}_{v}^{(0,1)}$ and $D_{v}^{(1,0)}$,

$\displaystyle\begin{split}1=&A_{v}^{(0,1)}\bar{D}_{v}^{(0,1)}\left(A_{v}^{(0,1)}\right)^{-1}A_{v}^{(1,0)}\bar{D}_{v}^{(0,1)}\left(A_{v}^{(1,0)}\right)^{-1}\\\
\times&A_{v}^{(1,1)}D_{v}^{(1,0)}\left(A_{v}^{(1,1)}\right)^{-1}D_{v}^{(1,0)}(-1)^{B_{6,12,11,4}^{(2)}}.\end{split}$
(35)

A close inspection of this equation shows that this equation is satisfied by

$\displaystyle\bar{D}_{v}^{(0,1)}=$ $\displaystyle 1\mbox{\quad and\quad}$
(36a) $\displaystyle D_{v}^{(1,0)}=$
$\displaystyle(-1)^{L_{5}^{(2)}B_{6,12,11,4}^{(2)}},$ (36b)

where $L_{j}^{(i)}=\frac{1}{2}(1-Z_{5}^{(i)})$ measures the value of the $i$th
tensor factor on edge $l_{j}$ giving rise to the generating phases

$\displaystyle D_{v}^{(0,1)}=$ $\displaystyle i^{B_{3,9,8,1}^{(1)}},$ (37a)
$\displaystyle D_{v}^{(1,0)}=$
$\displaystyle(-1)^{L_{5}^{(2)}B_{6,12,11,4}^{(2)}}$ (37b)

The remaining modification phase can be calculated using Eq. (32b) to obtain

$\displaystyle D_{v}^{(1,1)}=$ $\displaystyle
i^{B_{3,9,8,1}^{(1)}}(-1)^{(L_{5}^{(2)}+1)B_{6,12,11,4}^{(2)}}.$ (38)

These solutions define the modified operators

$\displaystyle\tilde{A}_{v}^{(0,1)}=$ $\displaystyle
i^{B_{3,9,8,1}^{(1)}}A_{v}^{(0,1)},$ (39a)
$\displaystyle\tilde{A}_{v}^{(1,0)}=$
$\displaystyle(-1)^{L_{5}^{(2)}B_{6,12,11,4}^{(2)}}A_{v}^{(1,0)}\mbox{\quad
and\quad}$ (39b) $\displaystyle\tilde{A}_{v}^{(1,1)}=$ $\displaystyle
i^{B_{3,9,8,1}^{(1)}}(-1)^{(L_{5}^{(2)}+1)B_{6,12,11,4}^{(2)}}A_{v}^{(1,1)},$
(39c)

which form a faithful representation of $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$.
When inserting these operators into Eq. (32c), one finds that they also
commute for neighboring vertices which completes our construction of non-Pauli
stabilizers
$\\{\tilde{A}_{v}^{(1,0)},\tilde{A}_{v}^{(0,1)},\tilde{A}_{v}^{(1,1)}\\}$
based on the $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ TQD model constructed with a
type-II cocycle. Note that – unlike in the $\mathbb{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_{v}^{(0,1)}$
adds an $i$ whenever $B_{3,9,8,1}^{(1)}=1\mod 2$, which can be seen as an
$S$-gate on the flux. Similarly, one can see that the phase for
$A_{v}^{(1,0)}$, $(-1)^{B_{6,12,11,4}^{(2)}L_{5}^{(2)}}$, is a controlled
$Z$-gate between the flux $B_{6,12,11,4}$ and the edge $l_{5}$.

Figure 3: Graphical representation of the stabilizers
$\\{\tilde{A}_{v}^{(0,1)},\tilde{A}_{v}^{(1,0)},\tilde{A}_{v}^{(1,1)}\\}$
(left to right). They are composed of TQD vertex operators (Eqs. (28)) –
represented by the filled dots on the vertex – and modification phases. The
latter are only supported on a restricted neighborhood which is colored red
(blue) when acting on the first (second) tensor factor of the local
$\mathbb{C}^{2}\otimes\mathbb{C}^{2}$ Hilbert space. The edges in gray only
enter implicitly via the fluxes.

With these modifications, we have constructed a set of stabilizers whose code
space is the ground space of a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$
topological order that cannot be _factored_ into two (possibly twisted)
$\mathbb{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.

###  3.3 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.
The method can also be applied for the input groups
$\mathbb{Z}_{N}=\\{0,1,\dots,N-1\\}$ and $\mathbb{Z}_{N}\times\mathbb{Z}_{N}$
with type-I and type-II cocycles (and combinations thereof). 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 $\mathbb{Z}_{N}$ and $\mathbb{Z}_{N}\times\mathbb{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 $\mathbb{Z}_{N}$,

$\displaystyle\omega_{I}(a,b,c)=e^{\frac{2\pi i}{N^{2}}a(b+c-[b+c]_{N})},$
(40)

where $[a+b]_{N}=(a+b)\mod N$, allows us to explicitly solve the consistency
equations and obtain pairwise commuting vertex operators
$\\{\tilde{A}_{v}^{g}\\}$ that represent the group action of $\mathbb{Z}_{N}$.
We exploit the cyclicity of $\mathbb{Z}_{N}$ by imposing that every vertex
operator to the $N$th power should equal the identity, just as in Eq. (19a)
for $N=2$. In particular, this should hold for the generating vertex operator
$A_{v}^{1}$, which allows us to determine a suitable ansatz for the
corresponding modification phase $D_{v}^{1}$. From this, we find a family of
solutions $D_{v}^{1}(p)$ that ensure that the generating vertex operators not
only represent the group but also commute pairwise. The fact that
$\tilde{A}_{v}^{1}$ generates every other $\tilde{A}_{v}^{g}$ allows us to
compute every other modification phase $D_{v}^{g}$ iteratively. One of the
resulting modification phases for any $g\in\mathbb{Z}_{N}$ reads

$\displaystyle D_{v}^{g}=$ $\displaystyle e^{\frac{2\pi
i}{N^{2}}gB_{6,12,11,4}}e^{-\frac{2\pi
i}{N}g\sum_{n=0}^{-L_{2}-g}\left(P_{-2,8,1}^{(n)}-P_{-2,9,3}^{(n)}\right)}$
$\displaystyle\times$ $\displaystyle e^{-\frac{2\pi
i}{N}\sum_{n=0}^{g-1}n\left(P_{-2,8,1}^{(-L_{2}-n)}-P_{-2,9,3}^{(-L_{2}-n)}\right)}$
(41) $\displaystyle\times$ $\displaystyle e^{-\frac{2\pi
i}{N}\sum_{n=0}^{g-1}n\left(P_{6,12,-5}^{(-L_{5}+n)}-P_{4,11,-5}^{(-L_{5}+n)}\right)},$

where $P_{i,j,k}^{(n)}$ 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_{-2,8,1}^{(n)}$ 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 $\mathbb{Z}_{2}$ case. Unlike the $\mathbb{Z}_{2}$
case, we have to take the orientation of the edges into account by subtracting
$l_{11}$ and $l_{4}$ since the elements of $\mathbb{Z}_{N}$ are not self-
inverses. 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 $\mathbb{Z}_{2}$ solution, Eq. (25). The second term, which only
includes the edge value $l_{2}$ in the upper summation bound, reduces to the
second term in Eq. (25) 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. (25). For a detailed derivation of the modification phase
for $G=\mathbb{Z}_{N}$ we refer to App. B and C.

When considering the TQD model with gauge group $\mathbb{Z}_{N}^{2}$ and a
type-II cocycle (Eq. (27)), the construction follows a similar path. For each
tensor factor, the vertex operators must fulfill the same closure relation as
in $\mathbb{Z}_{N}$. This allows us to find suitable ansatzes for the
modification phases $D_{v}^{(0,1)}$ and $D_{v}^{(1,0)}$ for the two generating
vertex operators $\tilde{A}_{v}^{(0,1)}$ and $\tilde{A}_{v}^{(1,0)}$ so that
the modified operators $\\{\tilde{A}_{v}^{(0,1)},\tilde{A}_{v}^{(0,1)}\\}$
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 $\mathbb{Z}_{N}\times\mathbb{Z}_{N}$,
we can again iteratively construct the modification phase for any
$g=(g_{1},g_{2})\in\mathbb{Z}_{N}^{2}$,

$\displaystyle D_{v}^{(g_{1},g_{2})}$ $\displaystyle=e^{\frac{2\pi
i}{N^{2}}g_{2}B_{3,9,8,1}^{(1)}}$ (42) $\displaystyle\times e^{-\frac{2\pi
i}{N}g_{1}\sum_{i=0}^{-L_{5}^{(2)}-1+g_{2}}\big{(}P_{(6,12,-5)^{(2)}}^{(i)}-P_{(4,11,-5)^{(2)}}^{(i)}\big{)}},$

where we indicate the operators only supported on the $i$th tensor factor with
an upper index in brackets, e.g. $L_{5}^{(1)}$ acts like
$L_{5}\otimes\mathds{1}$ on the $\mathbb{C}^{N}\otimes\mathbb{C}^{N}$ Hilbert
space on edge 5. The flux
$B_{3,9,8,1}^{(i)}=(-L_{3}^{(i)}-L_{9}^{(i)}+L_{8}^{(i)}+L_{1}^{(i)})\mod N$
measures the $i$th tensor factor of the flux through the diamond left of the
vertex (see. Fig. 3) and the projectors $P_{(i,j,k)^{(l)}}^{(n)}$ are defined
as in the $\mathbb{Z}_{N}$ case but on the $l$th tensor factor. As with the
$\mathbb{Z}_{N}$ solution, we can identify the terms derived in the previous
subsection for $N=2$ (see Eqs. (38) and (37)). The first term only depends on
fluxes and is non-trivial when $g_{2}\neq 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_{5}^{(2)}$ enters only in the upper summation bound
of the second term and has a similar influence on the sum of projectors as
$g_{2}$. In contrast to the one-parameter family of the $\mathbb{Z}_{N}$ case,
we have found only a single solution. For a detailed derivation, see App. D.

Although we have only calculated the modification phase for $\mathbb{Z}_{N}$
and $\mathbb{Z}_{N}\times\mathbb{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 finite Abelian
group can be decomposed into $\mathbb{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[29]. Whilst 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.

###  3.4 Towards quantum error correction

In this subsection, we discuss the potential of the new topological codes
devised here for notions of fault tolerant quantum computing and quantum error
correction [31]. We have constructed pairwise commuting operators
$\tilde{A}_{v}^{g}$ that, together with $B_{p}^{g}$ (see Eq. (4)) generate the
Non-Pauli stabilizer group $\mathcal{S}_{TQD}$. Just as in other topological
codes, a measurement of all generators of this group gives rise to a unique
excitation pattern. This diagnostic information can be used in _decoders_ that
allow for an eventual correction by annihilating them with suitable string-
like operators. The same string operators, when closed around a non-trivial
loop, give rise to logical operators of the codes.

Let us illustrate the qualitative difference of the string operators in our
codes to the toric code-like string operators in Kitaev’s _quantum double
models_ [20] with the previously discussed TQD model with $G=\mathbb{Z}_{N}$
and a non-trivial cocycle from Eq. (40). The group action can be implemented
by a generalized Pauli operator $X:\ket{n}\mapsto\ket{n+1\mod N}$. With that,
the vertex operators are of the form

$\displaystyle A_{v}^{g}=\alpha_{v}^{g}\prod_{l\sim v}X^{gs(l,v)}_{l},$ (43)

where $\alpha_{v}^{g}$ is a (in the computational basis) diagonal unitary
operator containing the pre-factor of the original vertex operator made up of
cocycles (see Eq. (6)) and the correction phase $D_{v}^{g}$ (see Eq. (41)).
The sign $s(l,v)=(-1)+1$ for an (outgoing)incoming edge $l$ w.r.t. vertex $v$
(compare Eqs. (7)). The plaquette operators are the known plaquette operators
$B_{p}^{g}$ made up of generalized Pauli $Z$s around the plaquette. The
elementary string operators of the corresponding code have the property that
they are supported on paths on the (dual) lattice and commute everywhere with
the stabilizers except at their endpoints where they only commute up to a
phase $e^{\pm 2\pi i/N}$ with either the plaquette or vertex operators. Since
the plaquette operators are exactly the same as in the _( $N$-level) toric
code_, the string operators that create a pair of vertex-excitations (charges)
are strings of Pauli $Z$s along a path connecting the two vertices that carry
the excitations. The elementary string operators supported along a (oriented)
path $\mathcal{P}$ in the dual lattice creating fluxes at their end-
plaquettes, on the other hand, have to be of the form

$\displaystyle
S^{g}_{\mathcal{P}}=\beta_{\mathcal{P}}^{g}\prod_{l\in\mathcal{P}}X_{l}^{gs(l,\mathcal{P})},$
(44)

where $\beta_{\mathcal{P}}^{g}$ is a diagonal operator similar to
$\alpha_{v}^{g}$ but supported on the path $\mathcal{P}$ and its nearest
neighbors and $s(l,\mathcal{P})=(-1)1$ for an edge $l$ that crosses
$\mathcal{P}$ from the (left)right. The reason for that is that a simple $X$
string on the dual lattice does not commute with the vertex operators due to
the non-trivial pre-factor $\alpha_{v}^{g}$. For the double semion code, the
string operators are depicted in Fig. 4.

Figure 4: Example of two elementary string operators of the double semion
code. A $Z$-string along a path on the lattice creates vertex excitations at
its endpoints (red dots). Similarly, a pair of plaquette excitations (blue
dots) is created by a $X$-string along a path on the dual lattice which is
dressed with a phase factor $\beta_{\mathcal{P}}$. The latter is supported on
the neighborhood of the path, whose boundary is depicted in green. The logical
operators of the code are generated by these string operators acting on along
non-trivial loops on the surface.

To find the pre-factors for the string operators $\beta_{\mathcal{P}}^{g}$ it
suffices to look at the generator $g=1$ and then define
$S_{\mathcal{P}}^{g}=\left(S_{\mathcal{P}}^{1}\right)^{g}$. Moreover, strings
along longer paths can be formed by concatenation of elementary strings along
paths of length 1. This allows us to generate any flux-string operator for the
code (and with that the logical operators as well) from generating string
operators that create excitations on neighboring plaquettes and then use the
same operators to move one excitation away from the other. These generating
strings can be obtained as solutions to the consistency conditions

$\displaystyle\left(S^{1}_{\mathcal{P}}\right)^{N}=$
$\displaystyle\mathds{1},$ (45a)
$\displaystyle\commutator{S^{1}_{\mathcal{P}}}{A_{v}^{1}}_{G}=$
$\displaystyle\mathds{1}\mbox{\quad\quad}\forall v,$ (45b)

where $\commutator{A}{B}_{G}=ABA^{-1}B^{-1}$ denotes the group commutator and
$N$ is the order of the corresponding excitation. The first condition enforces
the correct fusion rules of the underlying excitations and the second one
ensures that $S_{\mathcal{P}}^{1}$ is indeed an elementary string operator of
the code as discussed above. The authors of Ref. [13] Ref. have encountered
similar conditions in their version of a semion code ($N=2$ in our case) and
have identified an algorithmic way to construct the phase $\beta_{v}^{1}$.
This method should be generalizable to higher local dimensions $N>2$ and
models from product groups like the _twisted color codes_ introduced above.
This – as well as identifying ways of making the best use of the diagnostic
information obtained in non-Pauli stabilizer measurements – will be the topic
of future work.

What all TQD codes from a non-trivial cocycles have in common is the exotic
structure of the vertex operators and the flux-string operators, as
illustrated above. This results in strings of (generalized) Pauli $X$ errors
no longer only creating flux-excitations at their endpoints but also creating
correlated vertex excitations in its interior. In particular, such a Pauli
$X_{\mathcal{P}}$ string along a path $\mathcal{P}$ will decompose into
elementary string operators as

$\displaystyle
X_{\mathcal{P}}=\sum_{\\{P_{z}\\}}c_{\mathcal{P}}(P_{z})P_{z}S_{\mathcal{P}}^{1},$
(46)

where $\\{P_{z}\\}$ is the set of possible products of (generalized) Pauli
$Z$s supported on the neighborhood of path $\mathcal{P}$ and
$\\{c_{\mathcal{P}}(P_{z})\\}$ are the (complex) expansion coefficients that
can be obtained by projecting $\beta_{\mathcal{P}}^{1}$ onto the Pauli words
$\\{P_{z}\\}$. This means that Pauli errors act on the TQD codes in a similar
way as coherent errors in Pauli stabilizer codes. Investigating TQD codes with
Pauli noise therefore could also bring insight into the study of coherent
errors in Pauli codes.

##  4 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 $\mathbb{Z}_{2}\times\mathbb{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 $\mathbb{Z}_{N}$ and
$\mathbb{Z}_{N}^{2}$ 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 string operators
for these modified models, which provide both the recovery operations and the
logical operations on the code space. We expect that the feature of our codes
that Pauli $X$ strings (and their higher dimensional analogues) not only
create plaquette excitations at their endpoints but also some vertex
excitations along their path can be used to design a tailored decoder with
potentially increased performance under Pauli noise. A first step in this
direction has been taken in Ref. [32] using a neural network decoder for their
_semion code_. Although they did not find substantial evidence of increased
performance compared to the _toric code_ , it is far from clear whether there
does not exist a better decoder tailored to the exotic string operators and
what can be achieved with Stabilizer codes implementing other topological
orders – like the _twisted color code_ introduced in Sec. 3.2. Moreover, a
more sophisticated analysis of the behaviour of our codes under Pauli noise
would not only yield a deeper understanding of Non-Pauli topological codes in
general but also connect to Pauli (topological) codes with coherent errors: In
fact, here one has to deal with similar a decomposition of errors into
elementary strings, but in reversed roles. As far as we know, mappings to
statistical mechanics models for such error models [12], for example, are yet
to be discovered.

On a more abstract level, our codes host fundamentally different anyons to
those in Pauli stabilizer codes resulting in a different algebra of the
logical operators, twist defects and domain walls [18] and with that allow for
a different fault-tolerant gate set. With our stabilizers at hand, it is now
possible to implement these twisted anyon theories in an error correcting code
inviting research into their potential for topological quantum computation.
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 the comprehensive investigation into their
properties and their potential for fault-tolerant quantum computing.

## Acknowledgements

This work has been supported by the DFG (CRC 183, EI 519/15-1), and the
Templeton Foundation. We thank B. J. Brown, M. Kesselring, C. Wille, J.
Conrad, and J. Haferkamp for fruitful discussions and the JARA Institute for
Quantum Information at RWTH Aachen for its hospitality during a large time
period of this research. Moreover, we thank F. Hassler for supporting this
collaboration.

## References

  * Acín et al. [2018] Antonio Acín, Immanuel Bloch, Harry Buhrman, Tommaso Calarco, Christopher Eichler, Jens Eisert, Daniel Esteve, Nicolas Gisin, Steffen J Glaser, Fedor Jelezko, and et al.  The quantum technologies roadmap: a european community view.  _New Journal of Physics_ , 20(8):080201, Aug 2018\.  ISSN 1367-2630.  doi: 10.1088/1367-2630/aad1ea.
  * Andrist et al. [2015] Ruben S. Andrist, James R. Wootton, and Helmut G. Katzgraber.  Error thresholds for abelian quantum double models: Increasing the bit-flip stability of topological quantum memory.  _Phys. Rev. A_ , 91:042331, Apr 2015.  doi: 10.1103/PhysRevA.91.042331.
  * Bombin and Martin-Delgado [2006] H. Bombin and M. A. Martin-Delgado.  Topological quantum distillation.  _Phys. Rev. Lett._ , 97:180501, Oct 2006.  doi: 10.1103/PhysRevLett.97.180501.
  * Bombin et al. [2012] H Bombin, Guillaume Duclos-Cianci, and David Poulin.  Universal topological phase of two-dimensional stabilizer codes.  _New Journal of Physics_ , 14(7):073048, jul 2012\.  doi: 10.1088/1367-2630/14/7/073048.
  * Bombin [2018] Hector Bombin.  Transversal gates and error propagation in 3d topological codes.  2018\.
  * Bravyi and Terhal [2009] Sergey Bravyi and Barbara Terhal.  A no-go theorem for a two-dimensional self-correcting quantum memory based on stabilizer codes.  _New Journal of Physics_ , 11(4):043029, apr 2009\.  doi: 10.1088/1367-2630/11/4/043029.
  * Brown [2020] Benjamin J. Brown.  A fault-tolerant non-clifford gate for the surface code in two dimensions.  6(21), 2020.  doi: 10.1126/sciadv.aay4929.
  * Calderbank and Shor [1996] A. R. Calderbank and Peter W. Shor.  Good quantum error-correcting codes exist.  _Phys. Rev. A_ , 54:1098–1105, Aug 1996.  doi: 10.1103/PhysRevA.54.1098.
  * Campbell [2014] Earl T. Campbell.  Enhanced fault-tolerant quantum computing in $d$-level systems.  _Phys. Rev. Lett._ , 113:230501, Dec 2014.  doi: 10.1103/PhysRevLett.113.230501.
  * Chen [2017] Xie Chen.  Symmetry fractionalization in two dimensional topological phases.  _Reviews in Physics_ , 2:3–18, 2017.  ISSN 2405-4283.  doi: https://doi.org/10.1016/j.revip.2017.02.002.
  * Chen et al. [2013] Xie Chen, Zheng-Cheng Gu, Zheng-Xin Liu, and Xiao-Gang Wen.  Symmetry protected topological orders and the group cohomology of their symmetry group.  _Phys. Rev. B_ , 87:155114, Apr 2013.  doi: 10.1103/PhysRevB.87.155114.
  * Chubb and Flammia [2019] Christopher T. Chubb and Steven T. Flammia.  Statistical mechanical models for quantum codes with correlated noise.  2019\.
  * Dauphinais et al. [2019] G Dauphinais, L Ortiz, S Varona, and M A Martin-Delgado.  Quantum error correction with the semion code.  _New Journal of Physics_ , 21(5):053035, may 2019\.  doi: 10.1088/1367-2630/ab1ed8.  URL https://doi.org/10.1088/1367-2630/ab1ed8.
  * Dennis et al. [2002] Eric Dennis, Alexei Kitaev, Andrew Landahl, and John Preskill.  Topological quantum memory.  _Journal of Mathematical Physics_ , 43(9):4452–4505, 2002.  doi: 10.1063/1.1499754.
  * Ghosh et al. [2012] Joydip Ghosh, Austin G. Fowler, and Michael R. Geller.  Surface code with decoherence: An analysis of three superconducting architectures.  _Phys. Rev. A_ , 86:062318, Dec 2012.  doi: 10.1103/PhysRevA.86.062318.
  * Gottesman [1997] D. Gottesman.  Stabilizer codes and quantum error correction.  _arXiv quant-ph/9705052_ , 1997.
  * Hu et al. [2013] Yuting Hu, Yidun Wan, and Yong-Shi Wu.  Twisted quantum double model of topological phases in two dimensions.  _Phys. Rev. B_ , 87:125114, Mar 2013.  doi: 10.1103/PhysRevB.87.125114.
  * Kesselring et al. [2018] Markus S. Kesselring, Fernando Pastawski, Jens Eisert, and Benjamin J. Brown.  The boundaries and twist defects of the color code and their applications to topological quantum computation.  _Quantum_ , 2:101, October 2018.  ISSN 2521-327X.  doi: 10.22331/q-2018-10-19-101.
  * Kitaev [2006] Alexei Kitaev.  Anyons in an exactly solved model and beyond.  _Annals of Physics_ , 321(1):2–111, 2006.  ISSN 0003-4916.  doi: 10.1016/j.aop.2005.10.005.  January Special Issue.
  * Kitaev [2003] A.Yu. Kitaev.  Fault-tolerant quantum computation by anyons.  _Annals of Physics_ , 303(1):2–30, 2003.  ISSN 0003-4916.  doi: 10.1016/S0003-4916(02)00018-0.
  * Knill and Laflamme [1997] Emanuel Knill and Raymond Laflamme.  Theory of quantum error-correcting codes.  _Phys. Rev. A_ , 55:900–911, Feb 1997.  doi: 10.1103/PhysRevA.55.900.
  * Kubica et al. [2015] Aleksander Kubica, Beni Yoshida, and Fernando Pastawski.  Unfolding the color code.  _New Journal of Physics_ , 17(8):083026, aug 2015\.  doi: 10.1088/1367-2630/17/8/083026.
  * Levin and Gu [2012] Michael Levin and Zheng-Cheng Gu.  Braiding statistics approach to symmetry-protected topological phases.  _Phys. Rev. B_ , 86:115109, Sep 2012.  doi: 10.1103/PhysRevB.86.115109.
  * Levin and Wen [2005] Michael A. Levin and Xiao-Gang Wen.  String-net condensation: A physical mechanism for topological phases.  _Phys. Rev. B_ , 71:045110, Jan 2005.  doi: 10.1103/PhysRevB.71.045110.
  * Lidar et al. [2009] Daniel A. Lidar, Todd A. Brun, and Todd Brun, editors.  _Quantum Error Correction_.  Cambridge University Press, 2009.  doi: 10.1017/cbo9781139034807.
  * Litinski [2019] Daniel Litinski.  Magic State Distillation: Not as Costly as You Think.  _Quantum_ , 3:205, December 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-12-02-205.
  * Ni et al. [2015] X. Ni, O. Buerschaper, and M. Van den Nest.  A non-commuting stabilizer formalism.  _J. Math. Phys._ , 56:052201, 2015.  doi: 10.1063/1.4920923.
  * Nielsen and Chuang [2010] M. A. Nielsen and I. Chuang.  _Quantum computation and quantum information_.  Cambridge University Press, 2010.
  * Propitius [1995] M. de Wild Propitius.  Topological interactions in broken gauge theories.  _arXiv hep-th/9511195_ , 1995.
  * Steane [1996] A. M. Steane.  Error correcting codes in quantum theory.  _Phys. Rev. Lett._ , 77:793–797, Jul 1996.  doi: 10.1103/PhysRevLett.77.793.
  * Terhal [2015] Barbara M. Terhal.  Quantum error correction for quantum memories.  _Rev. Mod. Phys._ , 87:307–346, Apr 2015.  doi: 10.1103/RevModPhys.87.307.
  * Varona and Martin-Delgado [2020] S. Varona and M. A. Martin-Delgado.  Determination of the semion code threshold using neural decoders.  _Phys. Rev. A_ , 102:032411, Sep 2020.  doi: 10.1103/PhysRevA.102.032411.

##  Appendix A Group cohomology with $U(1)$ coefficients

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. [11] 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.

In our context, we deal with the cohomology of groups over $U(1)$.666Formally,
here $U(1)$ is a $G$-module with trivial group action. To define the $n$th
cohomology group, we start by defining maps from multiple copies of $G$ onto
$U(1)$, $\omega_{n}:G^{n}\to U(1)$. Such a general map is called $n$-cochain
and we denote the set of all such functions by $\mathcal{C}^{n}(G,U(1))$.
Moreover, we define the so called coboundary operator mapping $n$-cochains to
$n+1$-cochains,
$\delta_{n}:\mathcal{C}^{n}(G,U(1))\to\mathcal{C}^{n+1}(G,U(1))$ with

$\displaystyle\begin{split}(\delta_{n}\omega_{n})(g_{0},g_{1},\dots,g_{n})=&\omega_{n}^{s(g_{0})}(g_{1},\dots,g_{n})\prod_{i=0}^{n-1}\omega_{n}^{(-1)^{i+1}}(g_{0},\dots,g_{i-1},g_{i}g_{i+1},\dots,g_{n})\\\
&\times\omega_{n}^{(-1)^{n+1}}(g_{0},\dots,g_{n-1}),\end{split}$ (47)

where $s(g)=-1$ if $g$ is antiunitary and 1 if it is unitary and the group
multiplication symbol on $G$ between $g_{i}$ and $g_{i+1}$ is implicit. A
short calculation shows that
$(\delta_{n+1}\circ\delta_{n})\omega_{n}=1\;\forall\omega_{n}\in\mathcal{C}^{n}(G,U(1))$
which is the defining property of any coboundary operator. For $n=2$ and $n=3$
for example the coboundary operator acts as

$\displaystyle(\delta_{2}\omega_{2})(g_{0},g_{1},g_{2})=$
$\displaystyle\frac{\omega_{2}(g_{1},g_{2})\omega_{2}(g_{0},g_{1}g_{2})}{\omega_{2}(g_{0}g_{1},g_{2})\omega_{2}(g_{0},g_{1})},$
(48a)
$\displaystyle\begin{split}(\delta_{3}\omega_{3})(g_{0},g_{1},g_{2},g_{3})=&\frac{\omega_{3}(g_{1},g_{2},g_{3})\omega_{3}(g_{0},g_{1}g_{2},g_{3})\omega_{3}(g_{0},g_{1},g_{2})}{\omega_{3}(g_{0}g_{1},g_{2},g_{3})\omega_{3}(g_{0},g_{1},g_{2}g_{3})}.\end{split}$
(48b)

The first thing we do with an algebraic map is to define its kernel and its
image. We call any element in the kernel of $\delta_{n}$ $n$-cocycle and
denote the set of $n$-cocycles as
$\mathcal{Z}^{n}(G,U(1))=\\{\omega\in\mathcal{C}^{n}(G,U(1));\;\delta_{n}\omega=1\\}$.
We call any element in the image of $\delta_{n-1}$ n-coboundary and denote the
set of $n$-coboundaries by
$\mathcal{B}^{n}(G,U(1))=\\{\omega\in\mathcal{C}^{n}(G,U(1));\;\omega=\delta_{n-1}\beta,\beta\in\mathcal{C}^{n-1}(G,U(1))\\}$.

Due to $(\delta_{n}\circ\delta_{n-1})\omega_{n-1}=1$ we can multiply any
$n$-cocycle with an $n$-coboundary and the result will still be a $n$-cocycle.
Group cohomology classifies all the inequivalent cocycles under such a gauge
freedom, i.e. the $n$th cohomology group of $G$ over $U(1)$ is defined by

$\displaystyle
H^{n}(G,U(1))=\faktor{\mathcal{Z}^{n}(G,U(1))}{\mathcal{B}^{n}(G,U(1))}.$ (49)

By comparison of Eqs. (8) and (9) with Eqs. (48a) and (48b), we see that the
classification of the topological orders of a twisted gauge theory with gauge
group $G$ in Eq. (LABEL:eq:formalclassification) is exactly what we defined
here as the third cohomology group of $G$ over $U(1)$.

##  Appendix B Obstruction in general TQD models from Abelian groups

Consider a finite Abelian group $G$ and a 3-cocycle $\omega\in H^{3}(G,U(1))$.
By construction, the vertex operators $\\{A_{v}^{g}\\}$ in the original TQD
model defined by $(G,\omega)$ (see Eq. (6)) form a representation of the group
on the flux-free Hilbert space, where $B_{p}=1\;\forall p$, i.e., they
implement the group action on-site and commute with any vertex operator acting
on a different vertex. In the following, we will investigate the how they act
on the total Hilbert space $\mathcal{H}=\bigotimes\mathcal{H}_{l}$. In
particular, we find two types of obstructions – one which capture the failure
of the group multiplication rule when implemented on site (the _on-site
obstruction_) and one which captures the failure of two neighboring operators
to commute (the _commutativity obstruction_).

###  B.1 On-site obstruction

We will first investigate the successive action of two vertex operators
$A_{v}^{g}$ and $A_{v}^{h}$ on an arbitrary basis element. Using the action
defined in Eq. (6), we obtain

$\displaystyle\begin{split}A_{v}^{g}A_{v}^{h}\Bigg{|}\begin{minipage}[c]{62.59596pt}
\includegraphics[width=433.62pt]{basis_config}
\end{minipage}\Bigg{\rangle}=&\frac{\omega(l_{9},l_{3}h^{-1},h)\omega(l_{3}h^{-1},h,l_{4})\omega(h,l_{4},l_{11})}{\omega(l_{8},l_{1}h^{-1},h)\omega(l_{1}h^{-1},h,l_{6})\omega(h,l_{6},l_{12})}\\\
&\times\frac{\omega(l_{9},l_{3}^{\prime\prime},g)\omega(l_{3}^{\prime\prime},g,hl_{4})\omega(g,hl_{4},l_{11})}{\omega(l_{8},l_{1}^{\prime\prime},g)\omega(l_{1}^{\prime\prime},g,hl_{6})\omega(g,hl_{6},l_{12})}\Bigg{|}\begin{minipage}[c]{62.59596pt}
\includegraphics[width=433.62pt]{basis_config_prime_prime}
\end{minipage}\Bigg{\rangle},\end{split}$ (50)

where $l_{i}^{\prime\prime}=l_{i}h^{-1}g^{-1}$ for $i=1,2,3$ and
$l_{j}^{\prime\prime}=ghl_{j}$ for $j=4,5,6$. Note that we have left the group
multiplication symbol implicit. Besides shifting the edge values of $l_{1}$ to
$l_{6}$, the two vertex operators multiply the basis vector with a phase
factor given by a product of 12 cocycles. We can simplify this large product
using the cocycle condition Eq. (8) so that it reduces to a product of 10
cocycles

$\displaystyle\begin{split}\frac{\omega(l_{9}l_{3}^{\prime\prime},g,h)\omega(g,h,l_{4}l_{11})\omega(l_{9},l_{3}^{\prime\prime},gh)\omega(l_{3}^{\prime\prime},gh,l_{4})\omega(gh,l_{4},l_{11})}{\omega(l_{8}l_{1}^{\prime\prime},g,h)\omega(g,h,l_{6}l_{12})\omega(l_{8},l_{1}^{\prime\prime},gh)\omega(l_{1}^{\prime\prime},gh,l_{6})\omega(gh,l_{6},l_{12})},\end{split}$
(51)

where we can identify the pre-factor of $A_{v}^{gh}$

$\displaystyle A_{v}^{g}A_{v}^{h}\Bigg{|}\begin{minipage}[c]{65.44142pt}
\includegraphics{basis_config}
\end{minipage}\Bigg{\rangle}=\frac{\omega(l_{9}l_{3}^{\prime\prime},g,h)\omega(g,h,l_{4}l_{11})}{\omega(l_{8}l_{1}^{\prime\prime},g,h)\omega(g,h,l_{6}l_{12})}A_{v}^{gh}\Bigg{|}\begin{minipage}[c]{65.44142pt}
\includegraphics{basis_config} \end{minipage}\Bigg{\rangle}.$ (52)

Since this equality holds for any basis element and any group elements $g,h$,
we have established a relation between the operators $A_{v}^{g}A_{v}^{h}$ and
$A_{v}^{gh}$. By rewriting the product of edges in terms of the fluxes
$b_{3,9,2}=l_{3}^{-1}l_{9}^{-1}l_{2}$, $b_{2,8,1}=l_{2}^{-1}l_{8}l_{9}$,
$b_{5,11,4}=l_{5}l_{11}^{-1}l_{4}^{-1}$ and
$b_{6,12,5}=l_{6}l_{12}l_{5}^{-1}$, this relation reads

$\displaystyle\begin{split}A_{v}^{g}A_{v}^{h}=&\frac{\omega(l_{2}(b_{3,9,2})^{-1},g,h)\omega(g,h,(b_{5,11,4})^{-1}l_{5}h^{-1}g^{-1})}{\omega(l_{2}b_{2,8,1},g,h)\omega(g,h,b_{6,12,5}l_{5}h^{-1}g^{-1})}A_{v}^{gh}\\\
\eqqcolon&\Omega_{\\{b\\}}^{(l_{2},l_{5})}(g,h)A_{v}^{gh},\end{split}$ (53)

where we have omitted the group multiplication symbol for clarity and
introduced the on-site obstruction phase
$\Omega_{\\{b\\}}^{(l_{2},l_{5})}(g,h)$, defined with a fixed flux
configuration $\\{b\\}$, that is only supported on the horizontal edges
$l_{2}$ and $l_{5}$ and depends on the two group elements $g,h$.

By introducing the fluxes as above, we immediately see that the on-site
obstruction phase equals 1 in the flux-free subspace. Hence,
$A_{v}^{g}A_{v}^{h}=A_{v}^{gh}$ holds in the absence of fluxes, as
anticipated. On the whole Hilbert space however, group multiplication is only
faithfully implemented up to the phase factor $\Omega$.777One might be eager
to see the vertex operators as a projective representation of $G$. However,
one has to be careful with this, since $\Omega$ is in fact an operator that
itself does not commute with the vertex operators. Hence, where one might
think only of projective representations characterized by group cohomology
over a module with trivial action as discussed in the previous section, one
really is dealing with a group cohomology over a module with non-trivial
action characterizing the obstruction phase, sometimes referred to as twisted
group cohomology.

###  B.2 Commutativity obstruction

To quantify the (non-)commutativity of the vertex operators in the original
TQD model, we compute the group commutator
$\commutator{A_{v^{\prime}}^{g}}{A_{v}^{h}}_{G}\coloneqq
A_{v^{\prime}}^{g}A_{v}^{h}(A_{v^{\prime}}^{g})^{-1}(A_{v}^{h})^{-1}$. It is
clear from the definition that two vertex operators that act on vertices
separated by two or more edges commute. We therefore only have to consider
neighboring vertices $v$ and $v^{\prime}$. Due to the translation symmetry of
our lattice, we only need to consider three cases: Whether $v$ and
$v^{\prime}$ are connected by $l_{1}$, $l_{2}$ or $l_{3}$ (see Fig. 1).

In the first case, adapting the labelling of the vertices in Fig. 1 and using
Eq. (6), the group commutator, acting on an arbitrary basis vector, reads

$\displaystyle\begin{split}\commutator{A_{1}^{g}}{A_{v}^{h}}_{G}\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond2}
\end{minipage}\Bigg{\rangle}=&\frac{\omega(g,g^{-1}l_{1},l_{6})\omega(l_{8},g,g^{-1}l_{1})}{\omega(l_{8}g,g^{-1}l_{1},h)\omega(g^{-1}l_{1},h,h^{-1}l_{6})}\\\
&\times\frac{\omega(l_{8},l_{1},h)\omega(l_{1},h,h^{-1}l_{6})}{\omega(l_{8},g,g^{-1}l_{1}h)\omega(g,g^{-1}l_{1}h,h^{-1}l_{6})}\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond2}
\end{minipage}\Bigg{\rangle},\end{split}$ (54)

where we have only explicitly shown the part of the lattice on which both
operators act non-trivially. Using the two cocycle conditions (Eq. (8))

$\displaystyle\frac{\omega(l_{8},l_{1},h)\omega(l_{8},g,g^{-1}l_{1})}{\omega(l_{8}g,g^{-1}l_{1},h)\omega(l_{8},g,g^{-1}l_{1}h)}\stackrel{{\scriptstyle\eqref{eq:cocycle_condition}}}{{=}}$
$\displaystyle\frac{1}{\omega(g,g^{-1}l_{1},h)}\mbox{\quad and\quad}$ (55a)
$\displaystyle\frac{\omega(g,g^{-1}l_{1},l_{6})\omega(l_{1},h,h^{-1}l_{6})}{\omega(g^{-1}l_{1},h,h^{-1}l_{6})\omega(g,g^{-1}l_{1}h,h^{-1}l_{6})}\stackrel{{\scriptstyle\eqref{eq:cocycle_condition}}}{{=}}$
$\displaystyle\omega(g,g^{-1}l_{1},h),$ (55b)

we see that this pre-factor is in fact equal to $1$ for any cocycle. Since
this holds for any basis vector,
$\commutator{A_{1}^{g}}{A_{v}^{h}}_{G}=\mathds{1}$ on the whole Hilbert space
and the vertex operators $A_{1}^{g}$ and $A_{v}^{h}$ commute in the original
TQD model. We observe the same for the vertex operators $A_{3}^{g}$ and
$A_{v}^{h}$, since

$\displaystyle\begin{split}\commutator{A_{3}^{g}}{A_{v}^{h}}_{G}\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond3}
\end{minipage}\Bigg{\rangle}=&\frac{\omega(g,g^{-1}l_{1},l_{6})\omega(l_{8},g,g^{-1}l_{1})\omega(l_{8},l_{1},h)\omega(l_{1},h,h^{-1}l_{6})}{\omega(l_{8}g,g^{-1}l_{1},h)\omega(g^{-1}l_{1},h,h^{-1}l_{6})\omega(l_{8},g,g^{-1}l_{1}h)\omega(g,g^{-1}l_{1}h,h^{-1}l_{6})}\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond3}
\end{minipage}\Bigg{\rangle}\\\
\stackrel{{\scriptstyle\eqref{eq:cocycle_condition}}}{{=}}&\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond3}
\end{minipage}\Bigg{\rangle},\end{split}$ (56)

where we have again used the cocycle condition to identify the pre-factor in
the first line with $1$.

This leaves only the second pair of vertices, $v$ and $2$, where we keep to
the same procedure. However, when we compute the commutator, Eq. (6),

$\displaystyle\begin{split}\commutator{A_{2}^{g}}{A_{v}^{h}}_{G}\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond1}
\end{minipage}\Bigg{\rangle}=\frac{\omega(g^{-1}l_{9},l_{3},h)\omega(g,g^{-1}l_{9},l_{3})\omega(g,g^{-1}l_{8},l_{1}h)\omega(l_{8},l_{1},h)}{\omega(g^{-1}l_{8},l_{1},h)\omega(g,g^{-1}l_{8},l_{1})\omega(g,g^{-1}l_{9},l_{3}h)\omega(l_{9},l_{3},h)}\Bigg{|}\begin{minipage}[c]{39.83368pt}\includegraphics[width=433.62pt]{diamond1}
\end{minipage}\Bigg{\rangle},\end{split}$ (57)

we find that – in contrast to the previous cases – no amount of cocycle
manipulation can remove the pre-factor, leaving a phase of the form

$\displaystyle\frac{\omega(g^{-1}l_{9},l_{3},h)\omega(g,g^{-1}l_{9},l_{3})\omega(g,g^{-1}l_{8},l_{1}h)\omega(l_{8},l_{1},h)}{\omega(g^{-1}l_{8},l_{1},h)\omega(g,g^{-1}l_{8},l_{1})\omega(g,g^{-1}l_{9},l_{3}h)\omega(l_{9},l_{3},h)}=\frac{\omega(g,g^{-1}l_{8}l_{1},h)}{\omega(g,g^{-1}l_{9}l_{3},h)}.$
(58)

Since Eq. (57) holds for any basis vector we established an identity for the
operator $\commutator{A_{2}^{g}}{A_{v}^{h}}_{G}$. Expressed in terms of the
fluxes $b_{2,8,1}$ and $b_{3,9,2}$ introduced in Eq. (53) it reads

$\displaystyle\begin{split}\commutator{A_{2}^{g}}{A_{v}^{h}}_{G}=&\frac{\omega(g,g^{-1}l_{2}b_{2,8,1},h)}{\omega(g,g^{-1}l_{2}(b_{3,9,2})^{-1},h)}\eqqcolon\Pi_{\\{b\\}}^{(l_{2})}(g,h),\end{split}$
(59)

where we introduced the commutativity obstruction phase
$\Pi_{\\{b\\}}^{(l_{2})}(g,h)$ similar to the on-site obstruction phase
$\Omega$ in the previous section. Note that, for a given flux configuration
$\\{b\\}$ and group elements $g,h$, the commutativity obstruction is only
supported on one edge, namely the edge connecting the two vertices $v$ and
$2$. When expressed in this form, we can directly see that the commutativity
obstruction phase is 1 in the flux-free subspace, as anticipated. The goal of
our stabilizer construction is to modify the vertex operators in such a way to
remove $\Omega$ and $\Pi$ simultaneously.

##  Appendix C Constructing stabilizers from input group $\mathbb{Z}_{N}$ and
a type-I cocycle

In this section, we will explicitly calculate the obstructions defined above
for $G=\mathbb{Z}_{N}$. We represent the group by the set
$\\{0,1,\dots,N-1\\}$ together with the group multiplication being addition
modulo $N$, i.e. $g_{1}\cdot g_{2}=(g_{1}+g_{2})\mod{N}\eqqcolon g_{1}\oplus
g_{2}$. We first introduce operators that generalize the Pauli $Z$ and $X$
matrices to $N$ level systems. These are defined y their action on
$\mathcal{H}_{l}=\operatorname{span}_{\mathbb{C}}\\{\ket{0},\ket{1},\dots,\ket{N-1}\\}$,

$\displaystyle Z=$
$\displaystyle\sum_{n=0}^{N-1}\lambda^{n}\outerproduct{n}{n}\mbox{\quad
with\quad}\lambda\coloneqq e^{2\pi i/N}\quad\mbox{\quad and\quad}\quad
X=\sum_{n=0}^{N-1}\outerproduct{n\oplus 1}{n}.$ (60)

Note that $X^{N}=Z^{N}=\mathds{1}$ and that they satisfy the commutation
relation

$\displaystyle XZ=\lambda^{-1}ZX.$ (61)

Having defined those operators, we can write $A_{v}^{g}$ of any
$\mathbb{Z}_{N}$ TQD model in terms of a phase factor (defined in Eq. (6) by
the chosen cocycle $\omega$) times the product $\prod_{l\sim
v}X_{l}^{s(v,l)g}$, where $s(v,l_{j})=(-)1$ for an edge $l_{j}$ pointing (away
from)towards vertex $v$. From that, it directly follows that the phase
operator $Z$ fulfills the following commutativity relation with $A_{v}^{g}$:

$\displaystyle A_{v}^{g}Z_{j}=\lambda^{-s(v,l_{j})g}Z_{j}A_{v}^{g},$ (62)

where $j$ labels the edge on which $Z_{j}$ acts and
$s(v,l_{j})=1\leavevmode\nobreak\ (-1)$ for an edge $l_{j}$ pointing towards
(away from) vertex $v$ (compare Eqs. (7)). Moreover, we introduce projectors
$P^{(n)}_{l_{1},\dots,l_{k}}\coloneqq\sum_{n_{1}\oplus\dots\oplus
n_{k}=n}\outerproduct{n_{1},\dots,n_{k}}{n_{1},\dots,n_{k}}$ projecting onto
the subspace on which values888To be precise, we say an edge has a value
$l\in\mathbb{Z}_{N}$ when it is in a state vector $\ket{l}$. Hence, we can sum
up the values according to the group multiplication on $\mathbb{Z}_{N}$. of
the $k$ edges $l_{1},\dots,l_{k}$ sum up to $n\mod N$. Since $n$ is understood
modulo $N$, we define $P^{(n+N)}_{j}=P^{(n)}_{j}$ for any (set of) edge(s)
$j$.

To evaluate the obstruction phases $\Omega$ and $\Pi$ that we have defined in
App. B, we first need to choose a (non-trivial) cocycle representative
$[\omega]\in H^{3}(\mathbb{Z}_{N},U(1))$ to insert into Eqs. (53) and (59).
For $\mathbb{Z}_{N}$, there are $N$ cocycle classes that can be labelled by
$p\in\mathbb{Z}_{N}$ are represented by

$\displaystyle\omega_{p}(a,b,c)=e^{\frac{2\pi i}{N^{2}}pa(b+c-[b\oplus
c])}=\begin{cases}\lambda^{ap}&b+c\geq N\\\ 1&\text{else}\end{cases},\quad
a,b,c\in\mathbb{Z}_{N}.$ (63)

Such cocycles, depending only on elements from the same group
$\mathbb{Z}_{N}$, are called type-I cocycles [29]. All representatives are
generated by $\omega_{1}$ and $(\omega_{1})^{0}=(\omega_{1})^{N}\equiv 1$
represents the trivial cocycle class. Investigating the TQD model defined by
the cocycle $\omega_{1}$ and lifting the obstructions for that particular
model with a generating modification phase is therefore enough to lift the
obstructions for all $\mathbb{Z}_{N}$ TQD models. The modification phase that
lifts the obstruction for a model defined with an $\omega_{p}$ cocycle is
given as $p$th power of the generating modification phase.

Inserting $\omega_{1}$ from Eq. (63) into Eq. (53) gives the on-site
obstruction

$\displaystyle\Omega_{\\{b\\}}^{(l_{2},l_{5})}(g,h)=\left(\lambda^{-B_{3,9,8,1}}\right)^{\Delta_{g+h,N}}\lambda^{g\sum_{i=g}^{g+h-1}\left(P_{4,11,-5}^{(i-l_{5})}-P_{6,12,-5}^{(i-l_{5})}\right)}=\Omega_{\\{b\\}}^{(l_{5})}(g,h),$
(64)

where $B_{3,9,8,1}=(-l_{3})\oplus(-l_{9})\oplus l_{8}\oplus l_{1}$ and
$\Delta_{g+h,N}=1$ for $g+h\geq 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
$\omega_{1}$ into Eq. (59) gives the commutativity obstruction

$\displaystyle\Pi_{\\{b\\}}^{(l_{2})}(g,h)=\lambda^{g\sum_{i=g-h}^{g-1}(P_{-2,8,1}^{(i-l_{2})}-P_{-2,9,3}^{(i-l_{2})})}.$
(65)

Note that these obstructions coincide with the ones calculated for $N=2$,
where $g=h=1$ is the only non-trivial case, in Sec. 3.1.

Having quantified the obstructions for the $\mathbb{Z}_{N}$ TQD model, we
introduce the modified vertex operators $\tilde{A}_{v}^{g}=D_{v}^{g}A_{v}^{g}$
with the phase modification $D_{v}^{g}$ being a (in the edge basis) diagonal
operator with all entries taking values being roots of unity. We impose that
$\evaluated{D_{v}^{g}}_{\text{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

$\displaystyle\tilde{A}_{v}^{g}\tilde{A}_{v}^{h}=\tilde{A}_{v}^{g\oplus
h},\quad\forall g,h.$ (66)

Inserting the definition of $\tilde{A}_{v}^{g}$ and Eq. (53) yields the on-
site consistency condition on the phases $\\{\eta_{v}^{g}\\}$,

$\displaystyle\eta_{v}^{g\oplus
h}=A_{v}^{g}\eta_{v}^{h}(A_{v}^{g})^{-1}\eta_{v}^{g}\Omega_{\\{b\\}}^{(l_{5})}(g,h).$
(67)

Since $\eta_{v}^{g}$ is diagonal in the edge basis, conjugation with
$A_{v}^{g}$ modifies it only by shuffling edge values. In particular,
$A_{v}^{g}l_{2}(A_{v}^{g})^{-1}=l_{2}\oplus g$ and
$A_{v}^{g}l_{5}(A_{v}^{g})^{-1}=l_{5}\oplus(-g)$. Using the cyclic property of
$\mathbb{Z}_{N}$, we can set $g=1$ and $h=N-1$ in Eq. (67) to obtain an
equation for $\eta^{0}$. Since $A_{v}^{0}=\mathds{1}$, $\eta_{v}^{0}=1$, and
we obtain

$\displaystyle 1=\eta_{v}^{0}=\eta_{v}^{(N-1)\oplus
1}\stackrel{{\scriptstyle\eqref{eq:app_groupmultiplicationconsistency}}}{{=}}A_{v}^{1}\eta_{v}^{N-1}(A_{v}^{1})^{-1}\eta_{v}^{1}\Omega_{\\{b\\}}^{(l_{5})}(1,N-1).$
(68)

We expand this further by recursively writing $N-1=(N-2)\oplus 1$,
$N-2=(N-3)\oplus 1,\dots,2=1\oplus 1$ and using Eq. (67) to rewrite
$\eta_{v}^{N-1},\eta_{v}^{N-2},\dots,\eta_{v}^{2}$, thereby obtaining the
closure relation

$\displaystyle\begin{split}1=&\prod_{n=0}^{N-1}A_{v}^{n}\Omega_{\\{b\\}}^{(l_{5})}(1,-(n+1))(A_{v}^{n})^{-1}\prod_{n=0}^{N-1}A_{v}^{n}\eta_{v}^{1}(A_{v}^{n})^{-1}\\\
=&\prod_{n=0}^{N-1}\Omega_{\\{b\\}}^{(l_{5}-n)}(1,-(n+1))\prod_{n=0}^{N-1}A_{v}^{n}\eta_{v}^{1}(A_{v}^{n})^{-1},\end{split}$
(69)

where we have used that a diagonal operator (such as $\Omega$ and $\eta$)
conjugated by $A_{v}^{g}$ is still a diagonal operator and therefore commutes
with any other diagonal operator. Eq. (69) makes is possible to find a
solution for the phase corresponding to the generator of $\mathbb{Z}_{N}$,
$\eta_{v}^{1}$. Moreover, the recursion process that led us to Eq. (69) can be
used iteratively to generate all other phases $\eta_{v}^{g}\;\forall
g\in\mathbb{Z}_{N}$ so that Eq. (66) is fulfilled. Once that is achieved, the
general commutativity problem reduces to restoring commutativity to the
generator $\\{\tilde{A}_{v}^{1}\\}$ since any other modified vertex operator
decomposes as $\tilde{A}_{v}^{g}=(\tilde{A}_{v}^{1})^{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
$\Omega_{\\{b\\}}^{(l_{5})}$ calculated before, the equation reads

$\displaystyle
1=\lambda^{B_{3,9,8,1}}\lambda^{\sum_{n=0}^{N-1}\sum_{i=1}^{N-1-n}\left(P^{(i-l_{5}+n)}_{4,11,-5}-P_{6,12,-5}^{(i-l_{5}+n)}\right)}\prod_{n=0}^{N-1}A_{v}^{n}\eta_{v}^{1}(A_{v}^{n})^{-1}.$
(70)

The double sum in the exponent can be simplified with some projector algebra.
It reads

$\displaystyle\begin{split}\sum_{n=0}^{N-1}\sum_{i=1}^{N-1-n}\left(P^{(i-l_{5}+n)}_{4,11,-5}-P_{6,12,-5}^{(i-l_{5}+n)}\right)=&\sum_{n=0}^{N-1}\sum_{i=1+n}^{N-1}\left(P^{(i-l_{5})}_{4,11,-5}-P_{6,12,-5}^{(i-l_{5})}\right)\\\
=&-\sum_{n=1}^{N-1}n\left(P^{(n-l_{5})}_{4,11,-5}-P_{6,12,-5}^{(n-l_{5})}\right)=-B_{6,12,11,4},\end{split}$
(71)

where we have used that $\sum_{n=0}^{N-1}P^{(n)}_{j}=\mathds{1}$ and noting
that the final expression is the operator measuring the flux through the
diamond to the right of the vertex, $B_{6,12,11,4}=l_{6}\oplus
l_{12}\oplus(-l_{11})\oplus(-l_{4})$. With that, the closure condition on
$\eta_{v}^{1}$ reads

$\displaystyle
1=\lambda^{-B_{3,9,8,1}-B_{6,12,11,4}}\prod_{n=0}^{N-1}A_{v}^{n}\eta_{v}^{1}(A_{v}^{n})^{-1}.$
(72)

Since the vertex operators do not change fluxes, this equation can be solved
by any expression of the form

$\displaystyle\eta_{v}^{1}=\lambda^{(B_{3,9,8,1}+B_{6,12,11,4})/N}\bar{\eta}_{v}^{1}\mbox{\quad
with\quad}\prod_{n=0}^{N-1}A_{v}^{n}\bar{\eta}_{v}^{1}(A_{v}^{n})^{-1}=1.$
(73)

We are left with a freedom $\bar{\eta}_{v}^{1}$ to solve an additional
consistency equation coming from the commutativity obstruction phase (65). To
be precise, imposing that the generating vertex operators $\\{A_{v}^{1}\\}$
commute, i.e., $\commutator*{\tilde{A}_{2}^{1}}{\tilde{A}_{v}^{1}}_{G}=1$,
yields the commutativity consistency condition on $\\{\eta_{v}^{1}\\}$

$\displaystyle\begin{split}\Pi_{\\{b\\}}^{(l_{2})}(1,1)=&\eta_{v}^{h}A_{2}^{1}(\eta_{v}^{h})^{-1}(A_{2}^{1})^{-1}A_{v}^{1}\eta_{2}^{1}(A_{v}^{1})^{-1}(\eta_{2}^{1})^{-1}\\\
=&\commutator{\eta_{v}^{h}}{A_{2}^{1}}_{G}\commutator{A_{v}^{1}}{\eta_{2}^{1}}_{G}.\end{split}$
(74)

Substituting in Eqs. (73) and (65), we obtain the two consistency conditions

$\displaystyle\lambda^{P_{-2,8,1}^{(-l_{2})}-P_{-2,9,3}^{(-l_{2})}}=\bar{\eta}_{v}^{1}A_{2}^{1}(\bar{\eta}_{v}^{1})^{-1}(A_{2}^{1})^{-1}A_{v}^{1}\bar{\eta}_{2}^{1}(A_{v}^{1})^{-1}(\bar{\eta}_{2}^{1})^{-1}\mbox{\quad
and\quad}\prod_{n=0}^{N-1}A_{v}^{n}\bar{\eta}_{v}^{1}(A_{v}^{n})^{-1}=1$ (75)

to lift the on-site and commutativity obstructions. We find a one-parameter
family of solutions

$\displaystyle\bar{\eta}_{v}^{1}(p)=\lambda^{-(pB_{3,9,8,1}+(1-p)B_{6,12,11,4})/N}\lambda^{-p\sum_{n=0}^{-l_{2}-1}\left(P_{-2,8,1}^{(n)}-P_{-2,9,3}^{(n)}\right)-(1-p)\sum_{m=0}^{-l_{5}}\left(P_{6,12,-5}^{(m)}-P_{4,11,-5}^{(m)}\right)},$
(76)

where we have used a similar manipulation as in Eq. (71) to show the second
condition for $\bar{\eta}_{v}^{1}$. Note that
$\eta_{v}^{1}(p+N^{2})=\eta_{v}^{1}(p)$, and thus all distinct solutions in
this family are labeled by $p\in[0,N^{2})$. Putting these together with Eq.
(73) yields the full modification phase

$\displaystyle\eta_{v}^{1}(p)=\lambda^{((1-p)B_{3,9,8,1}+pB_{6,12,11,4})/N}\lambda^{-p\sum_{n=0}^{-l_{2}-1}\left(P_{-2,8,1}^{(n)}-P_{-2,9,3}^{(n)}\right)-(1-p)\sum_{m=0}^{-l_{5}}\left(P_{6,12,-5}^{(m)}-P_{4,11,-5}^{(m)}\right)}.$
(77)

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. (67), we generate all other modification phases
$\\{\eta_{v}^{g}\\}$ iteratively so that the on-site condition is fulfilled.
The general modification phase then reads

$\displaystyle\begin{split}\eta_{v}^{g}(p)=&\lambda^{g((1-p)B_{3,9,8,1}+pB_{6,12,11,4})/N}\lambda^{-g\left[p\sum_{n=0}^{-l_{2}-g}\left(P_{-2,8,1}^{(n)}-P_{-2,9,3}^{(n)}\right)+(1-p)\sum_{n=0}^{-l_{5}}\left(P_{6,12,-5}^{(n)}-P_{4,11,-5}^{(n)}\right)\right]}\\\
&\times\lambda^{-\sum_{n=0}^{g-1}n\left[p\left(P_{-2,8,1}^{(-l_{2}-n)}-P_{-2,9,3}^{(-l_{2}-n)}\right)+(2-p)\left(P_{6,12,-5}^{(-l_{5}+n)}-P_{4,11,-5}^{(-l_{5}+n)}\right)\right]},\quad
p\in\mathbb{R}.\end{split}$ (78)

With that, we have obtained a solution for a very general case, namely all
$\mathbb{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.
3.1).

##  Appendix D Constructing stabilizers from input group
$\mathbb{Z}_{N}\times\mathbb{Z}_{N}$ and a type-II cocycle

In this section, we will calculate the obstruction phases derived in App. B
for $G=\mathbb{Z}_{N}^{2}$ and a type-II cocycle. Each element in the input
group can be written as a pair of $\mathbb{Z}_{N}$ elements, $g=(g_{1},g_{2})$
with group multiplication naturally carrying over from $\mathbb{Z}_{N}$.999We
use the same symbol “$\oplus$” for the group multiplication on
$\mathbb{Z}_{N}$ and $\mathbb{Z}_{N}^{2}$. 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[29]. The two type-I
generators are the same as those for $\mathbb{Z}_{N}$, depending on the data
from a single tensor factor, i.e.,

$\displaystyle\omega_{I,p}(a,b,c)=\omega_{p}(a_{i},b_{i},c_{i}),\quad i=1,2,$
(79)

where $\omega_{p}$ was defined in Eq. (63). Using such a cocycle to define a
TQD model will result in the same functional form of the obstructions $\Omega$
and $\Pi$ found for $\mathbb{Z}_{N}$. Hence, they also can be removed by the
same modification phase $\eta_{v}^{g}$ from Eq. (78) where all
$\mathbb{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

$\displaystyle\omega_{II,p}(a,b,c)=\omega_{p}(a_{1},b_{2},c_{2})=\begin{cases}\lambda^{a_{1}p}&b_{2}+c_{2}\geq
N\\\ 1&\text{else}\end{cases},$ (80)

where $\lambda=e^{2\pi 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 [29]. Eq. (80) shows that cocycles of type II mixes the two tensor
factors of the input group elements in a non-trivial way. Whereas the
$\mathbb{Z}_{N}^{2}$ TQD model with a type-I cocycle can be decomposed into
two (possibly inequivalent) $\mathbb{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 $\Omega$ and $\Pi$ with a
type-II cocycle and investigate how to lift these with appropriately chosen
phase modifications.

Inserting the chosen type-II cocycle representative $\omega_{II,1}$ from Eq.
(80) into the obstruction phases calculated in Eqs. (53) and (59), we find
that the vertex operators of the type-II $\mathbb{Z}_{N}^{2}$ TQD model fail
to implement the group action faithfully on site, generating the obstruction
phase

$\displaystyle\begin{split}\Omega_{\\{b\\}}^{(l_{5})}(g,h)=&\lambda^{-B_{3,9,8,1}^{(2)}\Delta_{g^{(1)}+h^{(1)},N}-B_{3,9,8,1}^{(1)}\Delta_{g^{(2)}+h^{(2)},N}}\\\
&\times\lambda^{g^{(2)}\sum_{i=g^{(1)}}^{g^{(1)}+h^{(1)}-1}\big{(}P_{4,11,-5^{(1)}}^{(i-l_{5}^{(1)})}-P_{6,12,-5^{(1)}}^{(i-l_{5}^{(1)})}\big{)}}\lambda^{g^{(1)}\sum_{i=g^{(2)}}^{g^{(2)}+h^{(2)}-1}\big{(}P_{4,11,-5^{(2)}}^{(i-l_{5}^{(2)})}-P_{6,12,-5^{(2)}}^{(i-l_{5}^{(2)})}\big{)}},\end{split}$
(81)

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_{4,11,-5^{(1)}}^{(n)}$ projects onto the subspace where
$l_{4}^{(1)}\oplus l_{-1}^{(1)}\oplus(-l_{5}^{(1)})=n$. The quantity
$\Delta_{\bullet,\bullet}$ is defined as in Eq. (64). Note that this
obstruction phase consists of similar terms as those found for
$\mathbb{Z}_{N}$ case, but with added mixing between tensor factors. The
analoguous calculation of the commutativity obstruction phase yields

$\displaystyle\Pi_{\\{b\\}}^{(l_{2})}(g,h)=\lambda^{g^{(2)}\sum_{i=g^{(1)}-h^{(1)}}^{g^{(1)}-1}\big{(}P_{-2,8,1^{(1)}}^{(i-l^{(1)}_{2})}-P_{-2,9,3^{(1)}}^{(i-l^{(1)}_{2})}\big{)}}\lambda^{g^{(1)}\sum_{i=g^{(2)}-h^{(2)}}^{g^{(2)}-1}\big{(}P_{-2,8,1^{(2)}}^{(i-l^{(2)}_{2})}-P_{-2,9,3^{(2)}}^{(i-l^{(2)}_{2})}\big{)}}.$
(82)

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 $\mathbb{Z}_{N}^{2}$. 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
$\tilde{A}_{v}^{g}=\eta_{v}^{g}A_{v}^{g}\;\forall g\in\mathbb{Z}_{N}^{2}$ and
imposing Eq. (66). In particular, this directly implies a general condition on
the modification phases, Eq. (67). This enables us to use cyclicity to
iteratively derive the three consistency conditions from the relations
$(N-1,0)\oplus(1,0)=(0,0)$, $(0,N-1)\oplus(0,1)=(0,0)$ and
$(N-1,N-1)\oplus(1,1)=(0,0)$, yielding

$\displaystyle 1=$
$\displaystyle\eta_{v}^{(0,0)}=\eta_{v}^{(1,0)\oplus\dots\oplus(1,0)}=\prod_{n=0}^{N-1}A_{v}^{(n,0)}\Omega_{\\{b\\}}^{(l_{5})}((1,0),(-(n+1),0))(A_{v}^{(n,0)})^{-1}\prod_{m=0}^{N-1}A_{v}^{(m,0)}\eta_{v}^{(1,0)}(A_{v}^{(m,0)})^{-1}$
$\displaystyle=$
$\displaystyle\lambda^{-B_{\text{3,9,8,1}}^{(2)}}\prod_{m=0}^{N-1}A_{v}^{(m,0)}\eta_{v}^{(1,0)}(A_{v}^{(m,0)})^{-1},$
(83a) $\displaystyle 1=$
$\displaystyle\eta_{v}^{(0,0)}=\eta_{v}^{(0,1)\oplus\dots\oplus(0,1)}=\prod_{n=0}^{N-1}A_{v}^{(0,n)}\Omega_{\\{b\\}}^{(l_{5})}((0,1),(0,-(n+1)))(A_{v}^{(0,n)})^{-1}\prod_{m=0}^{N-1}A_{v}^{(0,m)}\eta_{v}^{(0,1)}(A_{v}^{(0,m)})^{-1}$
$\displaystyle=$
$\displaystyle\lambda^{-B_{3,9,8,1}^{(1)}}\prod_{m=0}^{N-1}A_{v}^{(0,m)}\eta_{v}^{(0,1)}(A_{v}^{(0,m)})^{-1}\mbox{\quad
and\quad}$ (83b) $\displaystyle 1=$
$\displaystyle\eta_{v}^{(0,0)}=\eta_{v}^{(1,1)\oplus\dots\oplus(1,1)}=\prod_{n=0}^{N-1}A_{v}^{(n,n)}\Omega_{\\{b\\}}^{(l_{5})}((1,1),(-(n+1),-(n+1)))(A_{v}^{(n,n)})^{-1}\prod_{m=0}^{N-1}A_{v}^{(n,n)}\eta_{v}^{(1,1)}(A_{v}^{(n,n)})^{-1}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:app_groupmultiplicationconsistency}}}{{=}}$
$\displaystyle\prod_{n=0}^{N-1}A_{v}^{(n,n)}\Omega_{\\{b\\}}^{(l_{5})}((1,1),(-(n+1),-(n+1))(A_{v}^{(n,n)})^{-1}$
$\displaystyle\times\prod_{m=0}^{N-1}A_{v}^{(m,m)}A_{v}^{(1,0)}\eta_{v}^{(0,1)}(A_{v}^{(1,0)})^{-1}\eta_{v}^{(1,0)}\Omega((1,0),(0,1))(A_{v}^{(m,m)})^{-1}$
$\displaystyle=$
$\displaystyle\lambda^{-\sum_{i=1}^{2}\left(B_{3,9,8,1}^{(i)}+B_{6,12,11,4}^{(i)}\right)}\prod_{m=0}^{N-1}A_{v}^{(m+1,m)}\eta_{v}^{(0,1)}(A_{v}^{(m+1,m)})^{-1}A_{v}^{(m,m)}\eta_{v}^{(1,0)}(A_{v}^{(m,m)})^{-1}.$
(83c)

The final pre-factors from Eqs. (83a), (83b) and (83c), where the root of
unity $\lambda$ appears with fluxes in the exponent, are produced by the
obstruction phases using the same identities as in Eq. (71). The first two
conditions are analogous to the one found in Eq. (72), and so we begin by
solving the first two equations Eq. (83a) and (83b) in a similar fashion as
the closure relation in the $\mathbb{Z}_{N}$ case, Eq. (72). We find that they
are easily solved by

$\displaystyle\eta_{v}^{(1,0)}=$
$\displaystyle\lambda^{B_{3,9,8,1}^{(2)}/N}\bar{\eta}_{v}^{(1,0)}\mbox{\quad
and\quad}$ (84a) $\displaystyle\eta_{v}^{(0,1)}=$
$\displaystyle\lambda^{B_{3,9,8,1}^{(1)}/N}\bar{\eta}_{v}^{(0,1)},$ (84b)

where
$\prod_{m=0}^{N-1}A_{v}^{(m,0)}\bar{\eta}^{(1,0)}(A_{v}^{(m,0)})^{-1}=\prod_{m=0}^{N-1}A_{v}^{(0,m)}\bar{\eta}^{(0,1)}(A_{v}^{(0,m)})^{-1}=1$.
Note that these two constraints on $\bar{\eta}^{(1,0)}$ and
$\bar{\eta}^{(1,0)}$ are fulfilled by any term that is a $N$th root of unity
and only depends on the second respectively the first tensor factor of the
link variables. Inserting this ansatz into the third closure relation Eq.
(83c), we obtain a closure relation for $\bar{\eta}$,

$\displaystyle
1=\lambda^{-B_{6,12,11,4}^{(1)}-B_{6,12,11,4}^{(2)}}\prod_{m=0}^{N-1}A_{v}^{(m+1,m)}\bar{\eta}_{v}^{(0,1)}(A_{v}^{(m+1,m)})^{-1}A_{v}^{(m,m)}\bar{\eta}_{v}^{(1,0)}(A_{v}^{(m,m)})^{-1}.$

Given the large freedom available when fulfilling the first two closure
relations, we can construct $\bar{\eta}^{(1,0)}$ such that it cancels out
$\lambda^{-B_{6,12,11,4}^{(2)}}$ and $\bar{\eta}^{(0,1)}$ such that it cancels
out $\lambda^{-B_{6,12,11,4}^{(1)}}$. One solution is given by

$\displaystyle\bar{\eta}^{(1,0)}=$
$\displaystyle\lambda^{-\sum_{i=0}^{-l_{5}^{(2)}-1}\big{(}P_{6,12,-5^{(2)}}^{(i)}-P_{4,11,-5^{(2)}}^{(i)}\big{)}}\mbox{\quad
and\quad}$ (85a) $\displaystyle\bar{\eta}^{(0,1)}=$
$\displaystyle\lambda^{-\sum_{i=0}^{-l_{5}^{(1)}-1}\big{(}P_{6,12,-5^{(2)}}^{(i)}-P_{4,11,-5^{(2)}}^{(i)}\big{)}}$
(85b)

since the tensor factors of $l_{5}$ are shifted by $m$ and $m+1$,
respectively, by the conjugation with $A_{v}^{(m,m)}$ and $A_{v}^{(m+1,m)}$,
respectively. Then, summing over all $m$ in the exponent (due to the product
over all $m$) as in Eq. (71) gives exactly the desired flux.

The full modification phases for the generators therefore read

$\displaystyle\eta^{(1,0)}=$
$\displaystyle\lambda^{B_{3,9,8,1}^{(2)}/N}\lambda^{-\sum_{i=0}^{-l_{5}^{(2)}-1}\big{(}P_{6,12,-5^{(2)}}^{(i)}-P_{4,11,-5^{(2)}}^{(i)}\big{)}},$
(86a) $\displaystyle\eta^{(0,1)}=$
$\displaystyle\lambda^{B_{3,9,8,1}^{(1)}/N}\lambda^{-\sum_{i=0}^{-l_{5}^{(1)}-1}\big{(}P_{6,12,-5^{(2)}}^{(i)}-P_{4,11,-5^{(2)}}^{(i)}\big{)}}.$
(86b)

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. (67) to
produce a proper representation of the group action. We will give an explicit
expression the modification phases for any group element
$g\in\mathbb{Z}_{N}^{2}$ 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 $\tilde{A}_{v}^{(0,1)}$ and
$\tilde{A}_{v}^{(1,0)}$. Therefore, it is sufficient to resolve the
commutativity obstruction for those two operators while still fulfilling Eqs.
(83). Imposing
$\commutator*{\tilde{A}_{2}^{g}}{\tilde{A}_{v}^{h}}_{G}=\mathds{1}\;\forall\,g,h\in\\{(0,1);(1,0)\\}$
and evaluating the obstruction phases $\Pi_{\\{b\\}}^{(l_{2})}(g,h)$ for the
corresponding $g,h$ yields the conditions

$\displaystyle
A_{2}^{(1,0)}(\eta_{v}^{(1,0)})^{-1}(A_{2}^{(1,0)})^{-1}\eta_{v}^{(1,0)}A_{v}^{(0,1)}\eta_{2}^{(1,0)}(A_{v}^{(0,1)})^{-1}(\eta_{2}^{(1,0)})^{-1}=$
$\displaystyle\Pi_{\\{b\\}}^{(l_{2})}((1,0),(1,0))\stackrel{{\scriptstyle\eqref{eq:app_ZNZNcommutativityobstruction}}}{{=}}1,$
(87a) $\displaystyle
A_{2}^{(0,1)}(\eta_{v}^{(0,1)})^{-1}(A_{v}^{(0,1)})^{-1}\eta_{v}^{(0,1)}A_{v}^{(0,1)}\eta_{2}^{(0,1)}(A_{v}^{(0,1)})^{-1}(\eta_{2}^{(0,1)})^{-1}=$
$\displaystyle\Pi_{\\{b\\}}^{(l_{2})}((0,1),(0,1))\stackrel{{\scriptstyle\eqref{eq:app_ZNZNcommutativityobstruction}}}{{=}}1,$
(87b) $\displaystyle
A_{2}^{(1,0)}(\eta_{v}^{(0,1)})^{-1}(A_{2}^{(1,0)})^{-1}\eta_{v}^{(0,1)}A_{v}^{(0,1)}\eta_{2}^{(1,0)}(A_{v}^{(0,1)})^{-1}(\eta_{2}^{(1,0)})^{-1}=$
$\displaystyle\Pi_{\\{b\\}}^{(l_{2})}((1,0),(0,1))\stackrel{{\scriptstyle\eqref{eq:app_ZNZNcommutativityobstruction}}}{{=}}\lambda^{P_{-2,8,1^{(2)}}^{(-l_{2}^{(2)}-1)}-P_{-2,9,3^{(2)}}^{(-l_{2}^{(2)}-1)}},$
(87c) $\displaystyle
A_{2}^{(0,1)}(\eta_{v}^{(1,0)})^{-1}(A_{2}^{(0,1)})^{-1}\eta_{v}^{(1,0)}A_{v}^{(1,0)}\eta_{2}^{(0,1)}(A_{v}^{(1,0)})^{-1}(\eta_{2}^{(0,1)})^{-1}=$
$\displaystyle\Pi_{\\{b\\}}^{(l_{2})}((0,1),(1,0))\stackrel{{\scriptstyle\eqref{eq:app_ZNZNcommutativityobstruction}}}{{=}}\lambda^{P_{-2,8,1^{(1)}}^{(-l_{2}^{(1)}-1)}-P_{-2,9,3^{(1)}}^{(-l_{2}^{(1)}-1)}}.$
(87d)

Surprisingly, the generating modification phases derived from the on-site
condition (see Eqs. (86)) also fulfill these four equations. As mentioned, Eq.
(67) now allows us to iteratively generate all modification phases for any
group element $g=(g_{1},g_{2})$. It reads

$\displaystyle\eta_{v}^{g}=\lambda^{(g_{1}B_{3,9,8,1}^{(2)}+g_{2}B_{3,9,8,1}^{(1)})/N}\lambda^{-g_{1}\sum_{i=0}^{-l_{5}^{(2)}-1+g_{2}}\big{(}P_{6,12,-5^{(2)}}^{(i)}-P_{4,11,-5^{(2)}}^{(i)}\big{)}-g_{2}\sum_{i=0}^{-l_{5}^{(1)}-1+g_{1}}\big{(}P_{6,12,-5^{(1)}}^{(i)}-P_{4,11,-5^{(1)}}^{(i)}\big{)}}.$
(88)

With that, we have resolved both obstructions for the $\mathbb{Z}_{N}^{2}$ TQD
model with a type-II cocycle.