Parafermions in a Kagome lattice of qubits for topological quantum computation

Engineering complex non-Abelian anyon models with simple physical systems is crucial for topological quantum computation. Unfortunately, the simplest systems are typically restricted to Majorana zero modes (Ising anyons). Here we go beyond this barrier, showing that the $\mathbb{Z}_4$ parafermion model of non-Abelian anyons can be realized on a qubit lattice. Our system additionally contains the Abelian $D(\mathbb{Z}_4)$ anyons as low-energetic excitations. We show that braiding of these parafermions with each other and with the $D(\mathbb{Z}_4)$ anyons allows the entire $d=4$ Clifford group to be generated. The error correction problem for our model is also studied in detail, guaranteeing fault-tolerance of the topological operations. Crucially, since the non-Abelian anyons are engineered through defect lines rather than as excitations, non-Abelian error correction is not required. Instead the error correction problem is performed on the underlying Abelian model, allowing high noise thresholds to be realized.


I. INTRODUCTION
Non-Abelian anyons exhibit exotic physics that would make them an ideal basis for topological quantum computation [1][2][3]. It has recently become apparent that truly scalable quantum computation with non-Abelian anyons can only be achieved when invoking active error correction, despite the protection provided by a finite anyon gap [4][5][6]. The development of practical systems in which non-Abelian anyons may be created, manipulated, and detected is therefore highly important. Systems in which non-Abelian anyons arise typically suffer from one of two drawbacks: either they are experimentally extremely challenging to realize (as is the case for quantum double [1] or string-net models [7]), or it is not clear how they can be made compatible with the active error correction required for fault-tolerance (as is the case for FQH systems).
A particularly attractive approach for building a faulttolerant quantum computer is to use a system of physical qubits (spin-1 2 particles). A number of technologies allow for precise qubit control, such as superconducting qubits [8], trapped atomic ions [9], spin qubits [10], or cold atoms or polar molecules in optical lattices [11]. A qubit lattice with two-body nearest neighbour interactions would therefore be an ideal system to realize non-Abelian anyons. While the model we propose involves three-body interactions, we discuss how these can be obtained from two-body interactions.
Non-Abelian anyons supported by a qubit system typically are Majorana zero modes, also known as Ising anyons [12][13][14][15]. A variety of proposals for experimental realization of Majorana zero modes in solid state systems have also been developed [16]. These anyons can be used to perform universal quantum computation when enhanced by non-topological operations [17,18]. However, these additional operations are highly resource intensive. Anyon models with a richer set of topological operations would therefore be much more practical for the realization of topological quantum computation. Here we solve this by introducing a model composed of two-qubit Hamiltonian interactions that can realize a more complex model of non-Abelian anyons, known as Z 4 parafermions. The error correction problem for these is studied in detail.
Extrinsic defects in Abelian topological states can behave like non-Abelian anyons [35]. The idea of non-Abelian anyons at the ends of defect lines, first introduced for FQH states [36], has been adapted to the D(Z d ) quantum double models in Refs. [37,38]. These anyons are Majorana zero-modes for d = 2 and more powerful parafermions for d > 2. Unfortunately, the generalized Pauli operators appearing in the D(Z d ) quantum double models models coincide with the physically relevant spin-operators only for d = 2. Otherwise, their structure makes them highly difficult to realize experimentally. The case d = 4, however, allows us to combine the best of both worlds. The joint Hilbert space of two qubits allows the 4-dimensional generalized Pauli operators to be expressed in terms of two-qubit operators. Using this, we show how Z 4 parafermions can emerge in a lattice of qubits with nearest-neighbor interactions only. This allows the computational power of Z 4 parafermions to be harnessed in a qubit system.
The fact that our system is built on top of a system supporting Abelian anyons (the D(Z 4 ) quantum double model) proves very useful. The non-Abelian parafermion modes can not only be braided with each other, but also with Abelian excitations of the quantum double model, allowing us to generate the entire Clifford group for d = 4 by quasi-particle braiding. This extends beyond the limited set of gates found using the same parafermions in previous work [21]. Furthermore, we do not have to perform non-Abelian error correction (a still poorly under-stood problem [4,5,[39][40][41][42]) to guarantee fault-tolerance, but can correct the underlying Abelian model. This Abelian error correction problem is nevertheless more involved than the well-studied error correction problem for the standard D(Z d ) models, and we study it in detail.
The rest of this paper is organized as follows. In Sec. II we show how Z 4 parafermion operators can be expressed in terms of qubit operators. Sec. III introduces a qubit Hamiltonian whose low-energetic excitations correspond to the D(Z 4 ) quantum double model. In Sec. IV we discuss how Z 4 parafermion modes appear at the ends of defect strings in our model. We demonstrate in Sec. V how the non-Abelian braiding statistics of these modes can be used to perform logical gates. Appendix D contains a proof that the set of gates which can be performed this way generates the entire Clifford group C 4 , which may be of independent interest. In Sec. VI we study the error correction problem of our model in detail and conclude in Sec. VII.

II. Z4 PARAFERMION OPERATORS IN TERMS OF QUBIT OPERATORS
We consider d-dimensional generalizations of the Pauli matrices X and Z. These are unitary operators satisfying X d = Z d = 1 1 and ZX = ωXZ, where ω = e 2πi/d with integer d > 1. If we define Y = ω (d+1)/2 X † Z † , we also have Y d = 1 1, XY = ωY X, and Y Z = ωZY . Operators X i and Z i act on qudit i and hence [X i , These operators are related to those of parafermions. Given a total ordering on the qudits {i}, one can obtain parafermion operators via a non-local transformation [43] These satisfy the Z d parafermion relations, The operators X, Y , and Z can be represented as ddimensional matrices. It is thus natural to seek a representation of these operators for the case d = 4 on the Hilbert space of two qubits (spins- 1 2 ). Indeed, given two qubits 1 and 2, one easily verifies that the operators are 4-dimensional generalized Pauli operators, and Z 4 parafermions can be obtained from these via Eq. (1). We also note that X 2 = σ x 1 σ x 2 , Y 2 = σ y 1 σ y 2 , and Z 2 = σ z 1 σ z 2 .
FIG. 1. Two qubits are located at each vertex of a Kagome lattice. Each pair of qubits hosts two Z4 parafermions. To unlock their potential for non-Abelian braiding, two such parafermions need to become unpaired, which is achieved by adding a defect line to the lattice. These are strings of strong local operators acting on qubit pairs (encircled). They create unpaired parafermion modes located at their ends (light pentagon-shaped regions consisting of a hexagon and a triangle).

III. MODEL
We consider a two-dimensional Kagome (trihexagonal) lattice as in Fig. 1. Each vertex of the lattice hosts one 4-dimensional qudit (one pair of Z 4 parafermions) or, in other words, two qubits. The Hamiltonian of our model is given by Here, the first term is a sum of equivalent terms for each triangle in the Kagome lattice. We label the vertices around one triangle a, b, and c, and the two qubits which are present at vertex a are called a 1 and a 2 , etc. The triangle terms in the Hamiltonian are then given by The second sum i in Eq. (4) runs over all vertices in the lattice. The two qubits located at vertex i are called i1 and i2. This second sum thus represents a uniform magnetic field in x-direction. Our Hamiltonian involves three-qubit terms of the form σ z a σ z b σ z c . It is in principle straightforward to generate these from one-body terms and two-body interactions by use of perturbative gadgets [44][45][46]. Consider a "mediator qubit" u coupled to qubits a, b, and c. Starting from a Hamiltonian and consider the perturbative regime ∆ |α|, |β|. In this regime, it is possible to integrate out qubit u. Taking up to third-order terms into account, one finds an effective Hamiltonian Choosing δ = −β and γ = 2 α 2 ∆ produces the desired three-qubit term without any undesired one-or two-qubit terms.
The generation of three-body interactions in optical lattices has been discussed in detail in Refs. [47,48]. These proposals would make the perturbative gadgets unnecessary. A "toolbox" for generating spin-lattice models such as ours in optical lattices has also been developed [49]. Generating non-Abelian anyons other than Majorana zero modes by use of perturbative gadgets from two-body interactions has previously been discussed in Refs. [50,51].
The spin-Hamiltonian in Eq. (4) can be exactly rewritten as Here again the first sum runs over all triangles in the lattice and the corners of a triangle are labeled a, b, and c. The second sum runs again over all vertices of the lattice. We now consider the perturbative limit h J and regard the second sum in Eq. (8) as a perturbation to the first term. Note that all terms in the first sum in Eq. (8) commute, so the unperturbed Hamiltonian is trivially solved. The lowest-order non-vanishing terms appear in sixth-order perturbation theory. We find an effective Hamiltonian where the second sum runs over all hexagons in the Kagome lattice and r, s, t, u, v, w label the six vertices around each hexagon. The effective Hamiltonian in Eq. (9) is derived in Appendix A. We note that all summands in H eff commute, so the system is exactly solvable. The excitations of this system are Abelian anyons corresponding to the D(Z 4 ) quantum double model. The topological degeneracy of the model can be made manifest by studying non-local loop degrees of freedom that commute with all stabilizers Z a Z b Z c , X r X † s X t X † u X v X † w , and their Hermitian conjugates, and fullfil themselves Z 4 relations. A possible choice of operators is illustrated in Fig. 2.
In passing, we note that the Z 2 version of Eq. (8), in which the Z 4 operators X and Z are replaced by Pauli operators σ x and σ z , leads to an effective Hamiltonian analogous to Eq. (9) and thus provides a very simple model with topological order. While this model requires three-body operators σ z σ z σ z as opposed to Kitaev's honeycomb Hamiltonian [12] which involves two-body interactions only, all of these interactions connect the same spin-component, which may provide a significant practical simplification over the honeycomb model.

IV. PARAFERMION MODES AND DEFECT LINES
The model is constructed from the cyclic qudit operators Z and X, which are related to parafermion operators. It is therefore natural to seek an interpretation of the model in terms of parafermionic modes.
To do this we must first fix the exact form of the stabilizers, which define the anyonic charge carried by each excitation. Let us use E p (M p ) to denote the stabilizer for a hexagonal (triangular) plaquette, p. For hexagonal plaquettes we use the convention that where r refers to the top-right corner and the other corners are labelled in an anti-clockwise fashion. For triangular plaquettes we use M p = Z a Z b Z c for all triangles of the form and M p = Z † a Z † b Z † c for all triangles of the form . The stabilizer operators E p and M p are unitary operators with eigenvalues ω k , k ∈ {0, 1, 2, 3}, where here and in the following ω = i for d = 4. An eigenvalue ω g of the E P corresponds to a charge anyon of the form e g , while M P similarly detects flux anyons m h . Fusion of charge anyons forms a representation of Z 4 , as does that of fluxes. The convention for the stabilizer operators chosen before ensures that the anyonic charge of both charge and flux type anyons is independently conserved (modulo 4). A full clockwise monodromy of an e g around an m h , or vice versa, yields a phase ω gh , see Fig. 3 for illustration.
Just as Majorana modes (Ising anyons) in the qubit toric code [14,38], parafermions appear in our system at the ends of defect strings. For the interpretation in terms of parafermions, it will be useful to introduce a new set of composite anyons defined as ψ g = e g × m g . These also obey Z 4 fusion with each other, and their braiding behavior can be inferred from the behavior of the constituent charge and flux particles. The particles {ψ 0 , ψ 1 , ψ 2 , ψ 3 } form a chiral Abelian anyon model with Chern number ν = 2 [12].
Note that We now perform a local transformation from the set of stabilizer generators {E p , M p }, detecting the charges on the left-hand-side of Eq. (10), to a new set {S p , R p } which detects the two charges on the right-hand-side. Let H denote the set of hexagonal plaquettes and T denote the set of triangular plaquettes. Note that |T | = 2|H|. Consider an injective map ϕ : H → T , which to each hexagonal operator E p assigns one of the six adjacent triangular operators M ϕ(p) . Typically, we choose M ϕ(p) to be the top-right neighbor of E p , while other choices become necessary next to defect lines. The transformation from the old to the new set of stabilizers reads S p = E p for p ∈ H and for p ∈ T . Here, Im(ϕ) denotes the image of the map ϕ.
Since p∈H S p = p∈T R p = 1 1, the charges detected by the new stabilizers are separately conserved (modulo 4). Just like the ψ g anyons detected by the S p stabilizers, the R g charges detected by the R p stabilizers also form an anyon model obeying Z 4 fusion. However, while these two anyon models have the same fusion rules, they are not equivalent, as they exhibit different braiding behavior. A full clockwise monodromy of a ψ g around a ψ h gives a phase of ω 2gh , a monodromy of a r g around an r h gives a phase of 1, and a monodromy of a ψ g around a r h gives a phase of ω gh , see again where the three particle models other than {ψ 0 , ψ 1 , ψ 2 , ψ 3 } correspond to the simple Z 4 model. The stabilizer operators S p detect the presence of ψ g anyon which are pinned to a pentagon-shaped double plaquette, made up of a neighbouring pair of triangular and hexagonal plaquettes. These anyons can be regarded as generalizations of Dirac fermions to the group Z 4 (rather than Z 2 ). Just as Dirac modes can be decomposed into two Majorana modes, so too can the ψ modes be decomposed into two parafermion modes. Two parafermion modes, P a and P b , are therefore associated with each double plaquette, P . These are described using parafermion operators satisfying Eq. (2). The parity operator for the ψ mode associated with a pair (j, k) is defined ω (d+1)/2 γ j γ † k for j < k, and so S P = ω 5/2 γ Pa γ † P b . For a stabilizer state, the system is within a definite eigenstate of all S P . The parafermion modes are therefore all paired, with the pairs corresponding to the two within each double plaquette. In order to use the parafermion modes as non-Abelian anyons, some must be allowed to become unpaired. The creation and transport of unpaired parafermion modes can be done by adapting the method of Ref. [38] to the Kagome lattice. The method can be interpreted in terms of anyonic state teleportation [52,53], as explained for the Majorana case in Ref. [54].
The method introduces unpaired parafermion modes at the endpoints of defect lines. These are lines on which additional single qudit terms are added to the Hamiltonian, of one of the two following forms , Specific examples are shown in Fig. 4. The single qudit terms added along defect lines are much stronger than any other interactions, and thus effectively remove the qudits on which they act from the code. This means that the E P and M P operators for the double plaquettes along these lines no longer commute with the Hamiltonian, and so can no longer be used as stabilizer generators. Their pentagon-shaped product, R P , is used instead. The pentagons in Fig. 4 show how next to a defect line the mapping ϕ needs to pick the bottom-left triangular-shaped stabilizer of a hexagonshaped stabilizer to ensure that their product still commutes with the Hamiltonian.
This change of the stabilizer generators of the code has a drastic effect. Consider an e g anyon moved towards a point along a defect line from one direction, and an m g moved towards the same point from the other direction. Both of these are detected by R P type stabilizers. When they meet on the same double plaquette, they will fuse to form a ψ g , and so not be detected by the R P stabilizers anymore. In fact, since the S P stabilizer is removed for double plaquettes along a defect line, they will not be detected by any stabilizer operator. The ψ g occupancy of a defect line corresponds to an increased groundstate degeneracy of the system, referred to as a synthetic topological degeneracy [38].
In the following section, the {ψ g , r h } decomposition of the D(Z 4 ) model will prove more convenient than the {e g , m h } decomposition. A process in which a defect line converts an m g into an e −g can equivalently be described as one in which a r g passes a defect line which emits a ψ −g .

V. PARAFERMIONS AS NON-ABELIAN ANYONS
Since unpaired parafermion modes reside at the endpoints of defect strings, it is natural to use them to explain the properties of the modified stabilizer. Parafermion modes are described by a non-Abelian anyon model with particle species {ψ 0 , ψ 1 , ψ 2 , ψ 3 , σ}. Here ψ 0 ≡ 1 corresponds to the anyonic vacuum and σ is an unpaired parafermion mode. The fusion rules of this anyon model are where ⊕ denotes addition modulo 4. A pair of parafermions (or the defect line between them) may therefore collectively hold any of the four types of ψ anyon.
The anyon model with the fusion rules given in Eq. (14) does not allow for a non-trivial solution of the pentagon and hexagon equations. As such, it obeys only projective non-Abelian statistics. The computational power of braiding parafermions has recently been studied in Ref. [55]. In our setup, we cannot only braid the parafermions with each other, but can also braid the Abelian e-and m-particles around them, which provides the possibility to perform additional gates. In the following, we want to study the gate set that can be generated this way.
As in the Majorana/Ising case, we use four parafermion modes (two defect strings) for which the total fusion sector is vacuum to store one logical qudit. The natural logical operators are parity operators for the pairs of parafermions. An eigenvalue ω g corresponds to a ψ g occupancy for the pair, and so the specific result σ ×σ = ψ g if they would be fused. We associate the Z basis of the logical qudit with the ψ occupancy of vertical pairs (connected by defect lines).
Specific choices of logical operator are illustrated in Fig. 4. TheZ L corresponds to a clockwise loop of an e 1 around a defect line. The braiding of this e 1 around the ψ g held in the pair yields the required phase of ω g . Thẽ X L corresponds to clockwise loop of an e −1 anyon which is converted to an m 1 through one defect line and back to an e −1 through the other. Equivalently, we can describe it as a clockwise loop of a r 1 and a transfer of a ψ 1 from the right to the left defect line.
Let us denote a state in which the left defect line holds a mode ψ g and the right defect line holds a mode ψ h by |ψ g , ψ h . Two defect lines create a 4 × 4-fold synthetic topological degeneracy. For computational purposes, we restrict to the 4-dimensional subspace of states of the form |g L ≡ |ψ g , ψ −g . This is the set of states which can locally be created from the anyonic vacuum. The effect of the logical operators on these states isX L |g L = |g ⊕ 1 L andZ L |g L = ω g |g L . In addition to the logical operatorsX L andZ L , which can be performed in our model by braiding the Abelian D(Z 4 ) anyons around the parafermion modes (ends of defect strings), we can perform further topologically protected single-qudit and two-qudit gates by braiding the parafermion modes themselves. Defect lines used for braiding are shown in Appendix B. Crucially, braiding parafermions allows one to perform an entangling gate by topological means, which is in contrast to Majorana fermions [21]. What is more, exploiting the fact that our non-Abelian system is built on top of an Abelian D(Z 4 ) system allows us to generate the entire 4-level Clifford group by braiding quasi-particles, as we discuss in the following.
For the rest of this section, X and Z refer to the logical operators calledX L andZ L before, respectively. The first column in Fig. 5 illustrates how braiding of D(Z 4 ) charges and fluxes can be used to perform logical X and Z gates. Whether an e 1 or an m 1 anyon is used to perform the logical Z gate is irrelevant.
Consider two ?parafermion modes storing a ψ g particle. A full clockwise monodromy of one parafermion around the other can be understood as a monodromy of the constituent e g around the m g , yielding an ω g 2 phase. We can thus expect a single exchange of the two parafermion modes storing a ψ g to yield a square root of this phase, such as ω g 2 /2 . This is demonstrated directly by studying the necessary microscopic operations in App. C.
For a logical qudit stored in four parafermion modes, let S denote a clockwise exchange of a vertical pair of parafermion modes, and T an exchange of a horizontal pair, see Fig. 6. As discussed, we have S = g ω g 2 /2 |g g|, while T is diagonal in the logical X ba-sis. In the logical Z basis, T reads (for d = 4) FIG. 6. Generators S and T of all gates that can be performed on a qudit stored in the fusion space of four parafermions by braiding them.
Again, in contrast to Majorana fermions, parafermions support an entangling gate between two logical qudits by braiding operations [21]. The controlled phase-gate Λ is defined by its action on a logical two-qudit basis state, Λ |g, h = ω gh |g, h . In our parafermion scheme, an entangling gate can be performed by braiding of a pair of parafermions from one qudit with a pair from the other. Let us consider, for example, the braiding of the left vertical pair for both qudits. For an initial logical product state |g, h , the process corresponds to braiding a ψ g clockwise around a ψ h , which yields a phase of ω 2gh . The resulting operation is therefore the squared controlled phase-gate Λ 2 . For d = 2, corresponding to the Ising/Majorana case, Λ 2 = 1 1, and so this operation is trivial. For d > 2, however, it is a non-trivial entangling gate, akin to the one proposed in Ref. [21].
Clearly a more powerful entangling gate would be Λ itself. This can be achieved for Z d parafermions for odd d by taking the (d + 1)/2-th power of Λ 2 . However these do not admit the simple decomposition into qubits that we have used in defining the model. Fortunately, we can make use of the underlying charge and flux anyons to realize Λ despite the even qudit dimension.
The defect line may be interpreted as a hole for ψ type anyons: an area in which they may be placed such that their state becomes delocalized along the line and they are no longer detected by the stabilizers [56][57][58]. Similar holes can also be engineered for the constituent charge and flux anyons. A defect line is therefore a special case of the combination of a charge and flux hole, in which only ψ g = e g × m g type anyons may reside rather than general e g × m h anyons. Nevertheless, we can consider a process in which a defect line is transformed into a separate charge and flux hole. Details on these holes and the transformations between them can be found in Appendix E.
When only the charge hole of one qubit is braided around the defect line of another, the process for an initial state |g, h corresponds to braiding an e g around a ψ h , which would yield the phase ω gh . The charge and flux holes can then be recombined into a defect line. The net effect of the entire process is to apply the controlled phase gate Λ. Such a process is illustrated in Fig. 7.
One process which could split the defect line in this way is simply to intersect it with two others. One would be a line along which charge anyons are hopped by highstrength terms. The other would similarly hop flux anyons. The stabilizers that detect charges and fluxes, respectively, along these lines would then be suppressed. By adiabatically removing the defect line which delocalizes ψ modes, its ψ g anyon occupation would be transferred to these two lines. The recombination of the defect line would be done by the reverse process.
For a tensor product of d-level systems, the Clifford group C d consists of gates that map tensor products of d-level Pauli operators to other such tensor products under conjugation. In Appendix D, we prove the following theorem.
Theorem. The single-qudit gates S, T , and Z, and nearest-neighbor controlled phase-gates Λ generate the entire Clifford group C 4 .
As an example,H = ST S = T ST satisfiesHXH † = Z andHZH † = X † , so it can be identified with the logical Hadamard gate, up to a phase. Indeed, using the standard definition one verifies that √ ωH =H. One possible implementation of H (up to a phase) is a cyclic permutation of the four parafermion modes, as in Fig. 5. This can be pictorially understood as follows.
An X corresponds to a transfer of a ψ 1 from the right to the left defect line, accompanied by a clockwise loop of a r 1 around a horizontal pair. A Z corresponds to a clockwise loop of a r 1 around a vertical pair. A π/2 rotation as performed by H thus maps these two operations onto each other, up to the fact that we do not perform a vertical ψ 1 transfer, as the ψ occupancy of the vertical pair is delocalized along the defect line.

VI. ERROR CORRECTION
For any system with a finite energy gap at finite temperature, excitations will appear with a finite density. This corresponds to finite length scale on which quantum computation can be performed before errors are almost certain to appear. This length scale can be increased by increasing the gap or lowering the temperature. However, neither of these methods is truly scalable. Error correction is therefore required if scalable quantum computation is to be performed.
For non-Abelian systems, the first studies of the corresponding error correction problem have recently appeared [4,5,[39][40][41][42]. Error correction for non-Abelian anyons is still poorly understood and its feasibility has not been demonstrated for the (realistic) time-continuous case. It comes thus very welcome that while our system provides the computational power of non-Abelian parafermions, its physical excitations still are Abelian D(Z 4 ) anyons, and the error correction problem for D(Z n ) quantum double models (including the timecontinuous case) is well-studied [39,40,[59][60][61]. However, when correcting these D(Z 4 ) anyons, we face a number of difficulties not considered in previous studies [39,[59][60][61]: (i) Our stabilizer operators are products of Z 4 qudit operators X, X † , Z, Z † , while an error model is realistically expressed in terms of single-qubit operators σ x , σ y , σ z . These do not map eigenstates of the stabilizer operators to other eigenstates and one single-qubit operator can produce a product of up to three qudit operators (see below).
(ii) We consider quantum information stored in a synthetic topological degeneracy, which involves a defect line allowing anyons to change from one sublattice to the other (stars to hexagons and vice versa). We thus cannot decode each sublattice separately, as usually done for the toric code and other D(Z d ) quantum double models, but have to correct both of them simultaneously while taking the possibility of transferring anyons from one to the other into account.
(iii) Besides simplistic i.i.d. error models (such as depolarizing noise), we are particularly interested in Hamiltonian protection of a quantum state subject to thermal errors. (iv) We do not consider a square lattice, but a trihexagonal one, which makes moving anyons and defining their distance more involved.

A. Error model
Since our 4-level qudits are composed of two qubits, it is natural to consider an error model in terms of singlequbit operations σ x , σ y , and σ z . For a qudit hosted in two qubits 1 and 2, single-qubit Pauli operators can be expressed in terms of Z 4 operators by inverting Eq. 3. We find If we start from an eigenstate of all stabilizer operators, applying single-qubit Pauli operators will generate a superposition of states corresponding to different syndrome outcomes. By measuring all stabilizer operators, we can project again into a subspace with definite syndrome values. Each single-qubit Pauli operator thereby translates into a product of up to three qudit operators. Table I summarizes (up to irrelevant phases) into which qudit operators a certain single-qubit Pauli operator will translate with equal probability. As a first simple error model, which does not involve a notion of Hamiltonian protection, we consider depolarizing noise. That is, for each qubit of the code we apply a Pauli operator with some probability p (the depolarization rate), where each of the three Pauli operators is chosen with equal probability.
More interesting from a physical perspective is a thermal error model. We consider a quantum state stored in the degenerate groundstates of the Hamiltonian given in Eq. (9), and assume that the system is weakly coupled to a heat bath at some temperature T . Following e.g. Ref. [5], we assume that evolving the system according to the Metropolis algorithm provides a reasonable approximation of the thermalization process, since the evolution obtained by means of the Metropolis algorithm is local, Markovian, and has the thermal state as its unique fixed point.
During our simulation, we proceed as follows. We first pick one of the spins-1 2 of the system at random, then pick one of the three single-qubit opertors acting on that qubit at random, and convert that to a 4-dimensional generalized Pauli operator according to Table I. We then calculate the energy cost ∆ tot of applying that generalized Pauli operator (or products thereof). This energy cost is of the form where ∆ and ∆ are the energy costs of creating a single triangle/hexagon-type anyon with charge 1 or 3 in Eq. (9), respectively. (That is, ∆ = 2J and ∆ = 2 63 8 h 6 (2J) 5 .) Creating an anyon with charge 2 will have an energy cost 2∆ or 2∆ . The coefficients m and n are elements of {0, ±2, ±4}, depending on the change in anyonic charge. The proposed error is then accepted with probability min{1, e −∆tot/k B T }. The noise model we apply is thus the standard classical Metropolis algorithm that maps eigenstates of Eq. (9) to other eigenstates. At any given time during our simulation, the system is "classical" in the sense that it does not involve superpositions of different anyon configurations.
If the proposed error is accepted, we copy the current state of the system and try to correct it. If correction is successful, we continue our simulation with the uncorrected version of the system. If correction fails (for at least one logical operator), we interpret this as the quantum information having survived for a time which is given by the number of Metropolis steps divided by the number of spins in the system. The thermal error model has three relevant energy scales k B T , ∆ , and ∆ . Since ∆ appears in higherorder perturbation theory than ∆ , we expect ∆ > ∆ . Furthermore, effective protection requires k B T < ∆ , ∆ . We introduce a parameter λ which quantifies the separation of these three energy scales, i.e., ∆ = λk B T and ∆ = λ 2 k B T . Very high values of λ are uninteresting, since they exponentially suppress errors from occurring.

B. Without defects
If there are no defect lines present, the anyonic charge of both types of anyons is conserved (modulo 4), and they can be corrected separately. Various techniques have been developed for correcting general D(Z n ) quantum double models [39,[59][60][61]. However, correcting the D(Z 4 ) case is particularly easy, since we can exploit the relation Z 4 /Z 2 Z 2 . Specifically, we can first fuse all oddly-charged anyons in pairs. In a second round, the  Fig. 2 as a function of the qubit depolarization rate p for code sizes L = 20, 28, 36, 44, 52. Each data point represents 10 4 logical errors, such that error bars are negligible. We recognize a threshold error rate pc ≈ 24% forX1 and pc ≈ 10% forZ1.
remaining anyons, which are all of charge 2, are fused in pairs. In order to find these pairings, we use the library Blossom V [62], which is the latest implementation of the efficient minimum-weight perfect matching algorithm due to Edmonds [63]. The weight between two equal-type anyons is thereby defined as the minimal number of generalized Pauli operators that need to be applied to create a pair of anyons at the two given locations. Fig. 8 shows our results for the depolarizing noise model, i.e., the logical error rates of the the logical op-eratorsX 1 andZ 1 illustrated in Fig. 2 as a function of the depolarization rate p. One clearly recognizes threshold error rates p c ≈ 24% and p c ≈ 10%, respectively. The equivalent figures for the logical operatorsZ 2 and X 2 look very similar and yield equivalent threshold error rates p c .
These thresholds are best compared with those for an equivalent code based on Z 2 anyons, and so with only a single qubit on each vertex. For independent bit and phase flips, the thresholds forX 1 andZ 1 are p c ≈ 16.4% and p c ≈ 6.7%, respectively [64,65]. When the hexagonal and triangular plaquettes are decoded separately, these correspond to thresholds of p c ≈ 24.6% and p c ≈ 10.5% for depolarizing noise. The similarity of these Z 2 values with those of Z 4 is remarkable. This qudit code is therefore just as adept at suppressing qubit noise as its qubit counterpart. It is well-known that the finite-temperature lifetime of a two-dimensional quantum memory with local interactions only is upper-bounded by a constant independent of the system size, see e.g. Ref. [66]. Fig. 9 shows the lifetime of a logical qudit with logical operators X 1 and Z 1 subject to the thermal error model. We notice lifetimes that decrease to an asymptotic value for large L and considerable finite-size tails. These tails correspond to the regime in which the breakdown of error correction is not due to the density of anyons becoming so high that pairing them becomes ambiguous, but where the breakdown is caused by one of the first pairs wandering along a topologically non-trival path around the torus. The smaller the system, the longer it takes to produce an anyon pair, leading to the observed tails for small enough L and T (large enough λ).

C. With defects
When defect lines as in Fig. 4 are present, the error correction problem becomes more involved. It is no longer possible to correct the two anyon types (hexagons and triangles in our case) separately, as is usually done for the D(Z n ) models [39,[59][60][61]. Instead, error correction needs to take the possibility of converting between different anyon types into account. We thus pair all oddlycharge anyons of both types in a first round and all remaining charge 2 anyons of both types in a second round. Pairings can involve anyons which are of equal or of different type. The weight for connecting two anyons is defined as the minimal number of generalized Pauli operators needed to create a pair of anyons at their respective positions from the vacuum. For equal-type anyons, this will be an error string that crosses an even number of defect lines, while for different-type anyons this will be an error string that crosses an odd number of defect lines. This can mean, for instance, that connecting two equal-type anyons can have a large weight despite them being geometrically nearby, if there is a defect line between them.
For a code of linear size L in both dimensions, with periodic boundary conditions and L even, we choose defect lines involving L/2 + 1 qudits, as shown in Fig. 4 for L = 20. The logical operatorsX L andZ L then have a distance L + 2 and L/2 + 4, respectively.
For the depolarizing error model, we find the threshold error rates p c for both of the logical operatorsX L and Z L given in Fig. 4. The results are given in Fig. 10. For the defect operatorX L , we find a threshold error rate p c ≈ 24%, as for the operatorsX 1 andX 2 in the defectfree case (Figs. 2 and 8), while for the defect operator Z L we find a threshold error rate p c ≈ 10%, as for the operatorsZ 1 andZ 2 in the defect-free case.
The fact that these values coincide with the defect free case is not unexpected. The introduction of the defects essentially corresponds to a change in the boundary conditions. However, the vast majority of errors have large support within the bulk. The value of the threshold is therefore dominated by bulk effects rather than boundary effects. Fig. 11 shows the average lifetime of the qudit stored in the synthetic topological degeneracy in Fig. 4 for the thermal error model. We note that for a given parameter λ, the asymptotic lifetimes (L → ∞) are close to those in the defect-free case given in Fig. 9.

VII. CONCLUSIONS
We have proposed a system which, on the physical level, involves only nearest-neighbor two-qubit interactions, allows one to perform all Clifford gates through quasi-particle braiding, and has a well-understood error correction problem.
We have greatly benefitted from the fact that our non-Abelian system is built on top of a system whose excitations correspond to an Abelian anyon model. This allows us to perform the logical operators X, Z, and Λ through quasi-particle braiding. It also makes our error correction problem manageable, despite some subtleties such as the fact that single-spin Pauli operators generate superpositions between different syndrome outcomes and the ability to convert between different anyon species during error correction.
Universal quantum computation requires the ability to perform non-Clifford gates, such as "small-angle" unitaries. While it is not difficult to perform a non-Clifford operation by non-topological means in our system, this abandons fault-tolerance. The technique of magic statedistillation [67] is typically used to restore fault-tolerance. While research on magic state distillation has so far focused on prime qudit dimensions d [68,69], qudit codes with the right transversality properties to perform magic state distillation in non-prime dimensions, including d = 4, also exist [70]. Unfortunately, for non-prime d it is not known whether Clifford gates plus an arbitrary non-Clifford gate are sufficient to achieve universality [71]. It is our hope that our work fuels interest in the d = 4 case, being a power of 2 and thus allowing to employ qubits, as demonstrated in our work, while being the smallest power of 2 that allows one to go beyond the Ising/Majorana case.
Alternatively, one could imagine energetically penalizing one of the degrees of freedom of a two-qubit Hilbert space to obtain a synthetic qutrit (d = 3). Magic state distillation for qutrits is well-studied [72], potentially allowing to perform fault-tolerant universal quantum computation with Z 3 parafermions in a qubit system [73].
The authors thank M. Barkeshli for elaborations on the development of the idea of generating non-Abelian defects in topological systems. This work was supported by the SNF, NCCR QSIT, and IARPA.
Appendix A: Sixth-order degenerate perturbation theory For our perturbation theory, we employ a Schrieffer-Wolff transformation [74], as formalized in Ref. [75].
Consider an unperturbed Hamiltonian H 0 whose spectrum can be separated into a low-and a high-energy subspace, which are energetically separated by a gap. Given a perturbation V , we want to find an effective Hamiltonian H eff describing the "effective" physics on the lowenergy subspace. The effective Hamiltonian can be developed in a perturbative series in powers of some small expansion parameter.
Let P denote the projector onto the low-energy subspace and Q = 1 1 − P the projector onto the highenergy subspace. We define V d = P V P + QV Q and V od = V − V d = P V Q + QV P . For some operator A, we define the superoperatorÂ viaÂ(O) = [A, O]. Let H 0 = i E i |i i| be the spectral decomposition of H 0 and define the superoperator L via We employ the convention that unless indicated otherwise by use of brackets, a superoperator L acts on all operators to its right. For the sixth-order effective Hamiltonian, one derives from Ref. [75] the expression where In our case, the low-energy subspace onto which P projects is given by the space in which all triangle operators in Eq. (8) have minimal energy, i.e., Z a Z b Z c ≡ 1 for all triangles (a, b, c). This subspace is fully degenerate. The lowest-energetic excitations change the eigenvalue of a stabilizer Z a Z b Z c from 1 to ±i. Since the eigenvalue of −(Z a Z b Z c + H.c.) is thereby changed from −2 to 0, this has an energy cost ∆ = 2J. Note, however, that stabilizer eigenvalues can only be changed in pairs, such that the gap between the low-energetic (groundstate) subspace and the space of excited states is in fact given by 2∆.
A crucial property of our Hamiltonian is that there is no lower-than-sixth-order perturbation that acts within the groundstate space. Therefore, we are only interested in terms of the form P V od (V d ) 4 V od P , which allows to greatly simplify the effective Hamiltonian. Namely, only the first summand in all expressions in Eqs. (A3) and (A4) is relevant in our case. We find Using now that in our case V d P = 0, this can be further simplified to There are 6! = 720 possibilities for applying the six factors X r X † s X t X † u X v X † w around one hexagon which leads the system back to the groundstate. Table II lists all possible routes the excitation energy above the groundstate can take, together with their numbers of possibilities.
In conclusion, we find the sixth-order effective Hamiltonian is given by the multiplicities in Table II, divided by the product of all excitation energies (in multiples of ∆) along the virtual process.  A selection of double plaquettes used in a braiding operation. Each double plaquette corresponds to two parafermion modes. For a double plaquette P , the parafermion to the right is labelled P1, and that to the left is P2.

Appendix B: Moving unpaired parafermion modes
To consider the creation and braiding of unpaired parafermionic modes, we must first decide on the double plaquettes with which we will work. Let us consider those of Fig. 12. To visualize the two parafermion modes within each double plaquette we use light blue circles. The one to the right of a double plaquette P is labelled P 1 , and that to the left is P 2 .
Parity operators for ψ modes are defined on pairs of parafermion modes. We are primarily concerned with two types of pairing: those of the two modes within the same double plaquette, and those of two modes from neighbouring double plaquettes. Relevant examples of the latter type are shown in Fig. 12 by red, orange and green lines connecting the corresponding modes.
For the two modes within each double plaquette, the parity operator ω 5/2 γ P1 γ † P2 corresponds to the stabilizer S P . The orange and red lines connecting modes P 2 to (P + 1) 1 denote the parity operators ω 5/2 γ P2 γ † (P +1)1 . For orange lines, these correspond to the operator Y on the vertex through which the line passes. For red lines they correspond to the operator X † Z.
Consider a state initially within the stabilizer space. The parity operators for the pairs of parafermion modes within each double plaquette are therefore part of the stabilizer. The state therefore corresponds to this definite pairing of the modes.
Let us now consider the removal of the operator S A from the set of stabilizer generators (while R A remains). The corresponding parafermion modes are now, in some sense, unpaired. This contributes a factor of four to the ground space degeneracy [76]. However, due to the fact that the 'unpaired' parafermions are not well separated, it is not difficult for local perturbations to lift the degeneracy of this space. To become truly unpaired, and benefit from topological protection, they must be separated.
To do this, we can add a term K(Y + Y † ) to the Hamiltonian, which corresponds to the parity operator ω 5/2 γ A2 γ † B1 . This acts on the vertex through which the orange line connecting these modes passes.
For K J, this new term will overwhelm the S B term. The pairing of B 1 and B 2 will then be broken, and B 1 will become paired with A 2 instead. The unpaired mode originally at A 2 is therefore effectively moved to B 2 . If the new term is introduced adiabatically, the degenerate subspace associated with the unpaired parafermion modes will remain in the same state during this process.
Similar processes can be used to move the unpaired modes further. The Hamiltonian term K(X † Z + X † Z) corresponding to ω 5/2 γ B2 γ † A1 can then be used to move the parafermion at B 2 to C 2 , for example. Unpaired parafermion modes can therefore be separated by arbitrary distances, at the endpoints of lines on which single qudit terms are added to the Hamiltonian. In terms of qubits, these correspond to two-body interactions between qubits in the same site.
In order to unlock the potential of parafermions for quantum computation, it must be possible to braid the parafermion modes. Let us consider a specific example of this, using the system of Fig. 12. Consider an initial state within the stabilizer space of all S P except A and L. At these two points, we have the unpaired parafermion modes A 1 , A 2 , L 1 and L 2 . Let us now consider operations such as those described above to move A 1 and L 2 away, beyond the bottom of the figure. All four parafermion modes are then well separated, and so the ground state degeneracy is topological protected.
We will now consider the exchange of A 2 with L 1 . We do this by first moving A 2 to K 2 , then L 1 to A 2 , and finally K 2 to L 1 . The two modes have then swapped places. An exchange of opposite chirality would correspond to first moving L 1 to K 1 , and so on. Note that all modes are kept well separated during the exchange, and so topological protection is always maintained. During the exchange, the movement of the modes is mostly achieved using Hamiltonian terms that correspond to the red and orange pairings in Fig. 12. These are all single qudit terms. At the junction, however, terms corresponding to the green pairings are used. We must therefore consider these in detail.
For the pairing show by the light green line, the parity operator is ω 5/2 γ P2 γ † O1 . This has the effect of creating an ψ g , ψ −g pair on the double plaquettes O and P . This requires the two qudit operator X † 3 X 4 Z † 4 . For the dark green pairing, the ω 5/2 γ P1 γ † G2 parity operator similarly requires the three qudit operator ω 5/2 X 1 X † 2 Z 2 Z 3 . These terms correspond to four-and six-body quasi-local interactions on the corresponding qubits, respectively. They can be realized by standard methods of perturbative gadgets. However, note that they need only be implemented while an exchange is in progress.

Appendix C: Exchange of two parafermion modes
To determine the effects of a single clockwise exchange we consider the smallest possible implementation. This involves the double plaquettes labelled F , G and P in Fig. 12, which are shown in more detail in Fig. 13. We consider a state in which the modes at P 1 and F 2 are unpaired, and those at G 1 and G 2 are paired by the Hamiltonian term S G . Using this, we determine the effects of exchanging the unpaired modes. The method used to exchange the two modes is similar to previous methods proposed in order to perform anyon braiding [52,53].
The results of the exchange are most easily understood in terms of the ψ mode formed by this pair. The parity operator, Γ for this mode is an operator that creates a ψ 1 , ψ −1 pair and places them in double plaquettes P and F , respectively. Also it must have eigenvalues of the form ω g , and so Γ 4 = 1. These conditions are satisfied by We similarly require operations that can move parafermions between the relevant plaquettes. These correspond to the green line between G and P and the red line between F and G. These are respectively. In these relations a, b and c are all elements of Z 4 . We have freedom in choosing the values a, b, c ∈ Z 4 for these relations. The corresponding freedom also exists for all operators used to move parafermion modes, as well as the logical operators. The values used do not simply correspond to differences in a global phase. Instead they determine which eigenspace of these operators has eigenvalue ω 0 = 1, and so which one corresponds to the vacuum occupancy ψ 0 of the ψ mode. These phases therefore cannot be chosen entirely arbitrarily, since the overall conservation constraint of ψ modes must be maintained. However, since here we do not explicitly consider the operations that placed unpaired parafermion modes at P 1 and F 2 , we can assume that their phases are chosen in a way that maintains this conservation. We will therefore consider a free choice of a, b, and c.
The first step in exchanging the parafermions is to move the one at P 1 to G 2 . This is done by adiabatically changing the Hamitonian to one in which the term Π + Π † is present and stronger than S G . This causes the modes at G 1 and P 1 to pair, moving the mode once at P 1 to G 2 . The mode at F 2 is then moved to that at P 1 by adiabatically changing the Hamitonian to one in which the Φ + Φ † term is present and stronger than S G , and the Π + Π † term is removed, pairing F 2 with G 1 . The mode at G 2 is then moved to F 2 by adiabatically removing the Π + Π † term and so allowing S G to become dominant and G 1 and G 2 to pair. This process then results in the clockwise exchange of the modes.
The first step of this transformation takes a state that is initially in the ω 0 eigenspace of S G and projects it to one in the ω 0 eigenspace of Π. The next step projects the state into the ω 0 eigenspace of Φ. The final step projects back into the ω 0 eigenspace of S G . The end effect is then P G P Φ P Π P G . Here P G is the projector onto the ω 0 eigenspace of S G , etc. The rightmost P G simply reflects the fact that the initial state lies within this eigenspace.
The parity operator Γ for the pair of unpaired modes commutes with S G , and so can be mutually diagonalized with the above operator. We can therefore interpret its effects in terms of the phase factor assigned to each of the possible ψ g eigenspaces of the ψ mode of the pair.
When doing this, different values of a, b, and c will result in different operations. This may seem to contradict the standard notion of a topological protected operation. However, these differences can be most easily understood by considering movement of parafermion modes implemented by measurement rather than adiabatic Hamiltonian manipulation. This method forces pairing of parafermion modes by measuring the occupancy of their corresponding ψ mode, and so forcing it to have a definite value. Ideally, this measurement will give the vacuum result ψ 0 . The effect is then the same as the adiabatic manipulation. If a different ψ g results, it must be removed by fusing it with the unpaired parafermion mode being moved. The different values used for the phases when moving unpaired modes, such as a, b and c here, determine how the measurement results are interpreted in terms of ψ anyons, and so determine the net ψ g fused with the modes being moved. As such, differences in the conventions used for an exchange will change the resulting operation only by a factor of Γ g , for some value of g that depends on a, b, and c.
For the standard convention used throughout this paper, with a = b = c = 2, the phase assigned to a ψ g occupation by the exchange is ω g 2 /2 ω g(g+1) . These are indeed all square roots of ω g 2 , as predicted in the main text. However, they are not of the elegant form ω g 2 /2 that would be more conducive for the proofs of Appendix D.
For an exchange operation that does have the required form, consider a = 1 and b = c = 2. The phase assigned to ψ g is ω −g 2 /2 in this case, up to a global phase of ω 1/2 . The required phases ω g 2 /2 would therefore be obtained from an anticlockwise exchange.
The fact that we obtain the phase ω −g 2 /2 for a clockwise exchange in this case means we would get ω −g 2 for a full clockwise monodromy. This would seem to contradict the arguments of the main text, which predict a phase of ω g 2 . However, note that these phases differ only by a factor of ω 2g 2 = ω 2g , and so are equivalent up to a factor of Γ 2 . Since such factors are to be expected for different choices of a, b and c, this monodromy does not contradict our expectations. The ambiguity in these phases reflects the fact that the anyon model described by the fusion rules in Eq. (14) obeys only projective non-Abelian statistics.
Note that the arguments above do not assume anything about the initial state of the exchanged parafermions. Only their initial positions and the operations used to move them are required. The effect of the braiding is expressed in terms of Γ, the parity operator for their shared ψ mode, which can be defined for any pair of parafermions. It therefore does not matter what state the fusion space of the parafermions was initially in, and it does not matter whether or not they exist at the end of the same defect line. The effect of the braiding is the same in all cases.

Appendix D: Generators of the Clifford group
For a tensor product of n d-level systems (C d ) ⊗n , the Pauli group P d is defined as the group generated by the generalized Pauli operators X i and Z i , and the Clifford group C d is defined as the normalizer of P d in the unitary group on (C d ) ⊗n . That is, elements in C d map tensor products of d-level Pauli operators to other such tensor products under conjugation.
We start with some general remarks on the action of C ⊗n Since the number of distinct eigenvalues is preserved under conjugation, this implies that the Pauli group P ⊗n d decays into distinct orbits when C ⊗n d acts on it by conjugation. This is in stark contrast to the case where d is an odd prime, which is studied in Ref. [77], where there is only one non-trivial orbit. The orbit containing the elements X 1 , Z 1 , . . ., X n , Z n consists of elements of the form where k, a 1 , b 1 , . . . , a n , b n ∈ Z 4 , at least one of the exponents a 1 , b 1 , . . ., a n , b n is odd, and p = n i=1 a i b i determines whether integer or half-integer powers of ω appear as phases (we sometimes write ω for i to avoid confusion with indices).
The following proof is an adaption of the proof in Appendix A of Ref. [77], where it is shown that a certain set of gates generate the Clifford group C ⊗n d for the case where d is an odd prime. The general structure of our proof is identical to the one in Ref. [77], while the generating set and individual lemmas and their proofs are different. After completion of this work, we became aware of the more general proof in Ref. [78].
Let us define the single-qudit unitaries and Lemma 1. The gates S † , T † , Z † , X, X † , and √ iH can all be generated from S, T , and Z.
Proof. As S 8 = T 8 = 1 1, we have S † = S 7 and T † = If for some Clifford gate U ∈ C ⊗n d we have with |α| = 1, we write M (U )(a 1 , b 1 , . . . , a n , b n ) T = (a 1 , b 1 , . . . , a We define the controlled Pauli-operators and Note that C Z has been called Λ in the main part of this work. We write A → U B as a shorthand for U AU † = B. We have and showing that C X , C Z ∈ C ⊗2 4 . Lemma 3. The gate C X can be generated from S, T , and C Z .
Proof. We note that Thus which together with Lemma 1 completes the proof.
Let us define a more general controlled operator as It acts by conjugation as Up to the phase ω st/2 , this action is identical to the one of the conditional Pauli gate C X s Z t studied in Ref. [77]. We point out again that such a phase is unavoidable for Z 4 , as there is, for instance, no unitary U such that X 1 → U X 1 X 2 Z 2 , since these two operators are not isospectral. Let us define the SWAP gate S via S |j |k = |k |j . Evidently, it acts as The gate S thus allows to generate non-local entangling gates from nearest-neighbor ones.
Lemma 4. The gate iS can be generated from S, T , and C Z .
Note that the gates S and iS act identically by conjugation.
Proof. Let One verifies that which together with Lemmas 1 and 3 completes the proof. Let All arithmetics involving the exponents a j , b j , c j , and d j that follow are to be understood modulo 4. It follows from the commutation relation ZX = ωXZ that P Q = ω (P,Q) QP where Lemma 5. Given P, Q ∈ P ⊗n 4 with (P, Q) = 1, we can generate U ∈ C ⊗n 4 from S, T , and nearest-neighbor C Z such that with |α P | = |α Q | = 1, and there exists j ∈ {1, . . . , n} such that a j d j − b j c j = 1.
Proof. Let P and Q be as in Eq. (D18). Since by assumption, there exists j such that If r j = 1, we are done. If r j = −1, then there is k = j with r k = 2 or r k = −1. From Lemma 2, we know that from S i and T i we can generate single-qudit unitaries that change ai ci bi di in arbitrary ways as long as r i = a i d i −b i c i is preserved. So up to gates that can be generated from S j and T j , we can assume that a j = 0, b j = c j = 1, and d j = 0. If r k = 2 then, up to gates that can be generated from S k and T k , we can assume that a k = 1, b k = c k = 0, and d k = 2. Finally, if r k = −1 then, up to gates that can be generated from S k and T k , we can assume that a k = 1, b k = 1, c k = −1, and d k = 2. We note that application of a phase-gate C Z(j,k) changes r j to r j = r j + (b k d j − b j d k ), and recall that non-local phase gates C Z(j,k) can be generated from nearest-neighbor ones and SWAP gates, which we can generate according to Lemma 4. In both cases (r k = 2 and r k = −1), we find that application of a phase-gate C Z(j,k) gives r j = 1. Lemma 6. Given P, Q ∈ P ⊗n 4 with (P, Q) = 1, we can generate U ∈ C ⊗n 4 from S, T , and nearest-neighbor C Z such that with P , Q ∈ P ⊗n−1
Proof. Let P and Q be as in Eq. (D18). By Lemma 5, we can assume that there is j with a j d j − b j c j = 1. Employing Lemma 4, we can perform a SWAP between qudits 1 and j. Finally, we perform a single-qudit unitary L on qudit 1 which is such that M (L) = dj −cj −bj aj . As det M (L) = 1 (mod 4), such a unitary L can be constructed from S and T according to Lemma 2. We note that with |α X | = |α Z | = 1, which completes the proof.
Lemma 7. For any V ∈ C ⊗n 4 we can construct U from S, T , Z, and nearest-neighbor Proof. Clearly, so by Lemma 6, we can assume that V X 1 V † = X ⊗ P and V Z 1 V † = Z⊗Q , up to gates that can be constructed from S, T , and nearest-neighbor C Z . Now let We define with i ∈ {2, . . . , n}. The gate with can be constructed from S, T , and nearest-neighbor C Z according to Lemmas 1, 3 and 4. Using Eqs. (D11) and (D14), we find the sequences of mappings and up to phases. Using again that XZX † =ωZ and ZXZ † = ωX, and that X can ge generated according to Lemma 1, allows us to generate arbitrary phases compatible with Eq. (D1). Since Theorem 1. Any Clifford gate V ∈ C ⊗n 4 can be constructed from S, T , Z, and nearest-neighbor C Z .
Proof. The proof is done by induction over n. The case n = 1 is given by Lemma 2. For n > 1, let U be as in Lemma 7. Since U † V commutes with X 1 and Z 1 , we have U † V = 1 1 ⊗ V , where V ∈ C ⊗n−1 4 acts on qudits {2, . . . , n}. Assuming that the induction hypothesis holds for n − 1, V and hence V can be constructed from S, T , and nearest-neighbor C Z . this basis will correspond to measurement (by the environment) in the Z basis. Any unitary that introduces unwanted relative phases can be expressed as a sum of powers of Z. As such, these processes must have support on all sites on which Z has support, which are all sites around the hole. Since the shrinking (and expansion) of holes does not have such support, it cannot cause any decoherence.
Corresponding processes can also be applied to hexagonal plaquettes. In that case, a logical qudit can be stored in the m g occupations of plaquettes. Such qudits can be topologically protected by using lines of Z j + Z † j terms to combine neighbouring horizontal plaquettes.
Using the processes of extending and shrinking holes, it is possible to move them. Gates can then be implemented through braiding. Braiding an e-type hole of triangular plaquettes in state |g around an m-type hole of hexagonal ones in state |h corresponds to braiding the e g anyon held by the former around the m h of the latter, yielding a phase ω gh . This is a qudit generalization of the controlled phase gate.

Fusing holes into defect lines
By considering the alternative stabilizer generators S and R discussed in Sec. IV, holes can also be created which hold ψ g anyons. These are formed by similarly combining the double plaquettes. Indeed, these are exactly the defect lines considered in the bulk of this paper. These have the property that anyons crossing the defect line undergo an automorphism that preserves the structure of the underlying Abelian state: it maps e anyons to their dual, the m anyons, and vice versa. This property is not shared by the e-and m-type holes. In these cases, the lines form a boundary along which one type of anyon can condense, but the other cannot cross. It is this difference that gives the ψ-holes additional properties, namely the localized parafermion modes at their endpoints, that the e-and m-holes do not possess. The topological degeneracy and protection, however, is a property shared by all three.
An e-type hole and an m-type hole together correspond to a two-qudit space. However, let us consider the subspace spanned by states |g, g . These are such that the e-type hole carries an e g anyon whenever the m-type hole carries an m g . A single qudit can be stored in this subspace Since ψ g = e g × m g , two holes as described above hold a net ψ g . Similar fusion can also be applied to holes, as we will now show. Specifically an e-type and an m-type hole can be combined into a ψ-type hole, and a ψ-type hole can be split into an e-type and m-type one.
These processes are in fact a simple generalization of the hole extension and shrinking processes described above. Suppose we have a defect line along double plaquettes P , P , P , . . . This stores a qudit whose basis states |g are ω g eigenstates of W P W P W P . . . Let us now extend this line. However, rather than adding another double plaquette, we instead add a triangular plaquette p. This can be done by adiabatically introducing the strong term X j + X † j to the Hamiltonian on a site shared by p and the triangular part of P . This term does not commute with either E p or W P . However, it does commute with their product. The final state will then have the qudit basis states defined by the operator E p W P W P W P . . .. Further such processes can be used to extend the e-type part of the defect line. Corresponding processes on the hexagonal plaquettes can be used to grow an m-type part of the defect line. An illustration is given in Fig. 14. FIG. 14. Red circles denote defects as in Eq. (13), with a pair of parafermions (purple) emerging at the ends of the defect line. Blue circles correspond to terms of the form Xj + X † j , while green circles correspond to terms of the form Zj + Z † j . These grow the e-and m-part of the defect line, respectively. Once both the blue and green defect lines have been added, the red defect line can be removed. The ψg particle initially stored in the red defect line has then been split into its eg and mg components.
When both e-type and m-type parts have been added, the original ψ-type defect line can be removed. This is done simply by removing the parity operator terms along its length and allowing the S P terms to again dominate. The end result is that the ψ g originally stored in the defect line now resides in the e-type hole as an e g and the m-type hole as an m g , corresponding to the state |g, g of their individual qudits. As for the shrinking of holes, this process does not have sufficient support to distinguish between different basis states. The process therefore does not decohere any superpositions of these states, nor does it assign any relative phases.
To recombine the two holes into a single defect line, the process is simply reversed. This will be straightforward if the two holes are in a state of the form |g, g , since the state of the defect line will simply become |g . However some processes, such as mistakes during error correction, could result in holes whose states are not of this form. This will introduce frustration that will not allow all of the involved triangular and hexagonal plaquettes to return to their ground state after recombination.
As an example, consider a state of the form |g, h . This corresponds to an e g and an m h , and could arise from an initial state |g, g if an m h−g were added in error to the m-type hole, or from |h, h with an e g−h error on the e-type hole.
The anyons e g and an m h can combine either to a ψ g and m h−g , or a ψ h and e g−h . The former will be energetically favourable, due to the weaker strength of the M p terms. The adiabatic process will therefore result in ψ g being stored on the defect line and an m h−g anyon present as an excitation on one of the triangular plaquettes that was once part of the holes. Syndrome measurement will then detect the m g−h . However, since the value g − h gives information only about the error that occurred, and not the value of g or h, this does not extract any information about the stored qudit.