A Hierarchy of Anyon Models Realised by Twists in Stacked Surface Codes

Braiding defects in topological stabiliser codes can be used to fault-tolerantly implement logical operations. Twists are defects corresponding to the end-points of domain walls and are associated with symmetries of the anyon model of the code. We consider twists in multiple copies of the 2d surface code and identify necessary and sufficient conditions for considering these twists as anyons: namely that they must be self-inverse and that all charges which can be localised by the twist must be invariant under its associated symmetry. If both of these conditions are satisfied the twist and its set of localisable anyonic charges reproduce the behaviour of an anyonic model belonging to a hierarchy which generalises the Ising anyons. We show that the braiding of these twists results in either (tensor products of) the S gate or (tensor products of) the CZ gate. We also show that for any number of copies of the 2d surface code the application of H gates within a copy and CNOT gates between copies is sufficient to generate all possible twists.


Introduction
The question of which quantum gates can be performed fault-tolerantly in a particular quantum error correcting code is of vital importance if we wish to use the code in quantum computation. The conventional method of achieving fault-tolerant operations is via the application of transversal gates [1]. In topological codes these have been generalised to locality preserving logical operators [2] [3] which map local errors to local errors. However, these codes also allow us to perform fault-tolerant gates using another method: T. R. Scruby: thomas.scruby.17@ucl.ac.uk braiding of topological defects. Examples of topological defects include punctures (produced by the removal of stabilisers) and twists (the endpoints of domain walls in the code). Braiding of these defects can allow us access to gates which are difficult or impossible to implement transversally. For example, the S gate cannot be implemented transversally in the 2d surface code (outside of folded surface code constructions such as [4]) but can be implemented via braiding of twist defects [5].
A recent paper by M. Kesselring et al [6] fully categorises the twists of the 2d colour code, sorting them into nine conjugacy classes. In light of this result it seems natural to ask what gates we can implement via the braiding of these twists. In this work we attempt to answer this question for at least some of the colour code twists. In [5] the fact that braiding twists produces an S gate is shown by considering the action of this braid on the logical operators of the code but the same result can be obtained by considering the twists as (Ising) anyons and analysing their braiding relations, as in [7]. Formally the twists in a topological code are described by G-crossed braided tensor categories [8] and cannot in general be considered as anyons. We will discuss below the cases in which neglecting the full G-crossed category treatment of these defects is permissable and we will see that three of the nine conjugacy classes of colour code twist are examples of cases where twists can be analysed using anyonic models. These models are members of a hierarchy of anyonic models generalising the standard Ising anyon model used to study twists in the surface code.
This paper is organised as follows: in section 2 we present a short overview of anyon fusion and braiding relations and twists in topological codes, touching briefly on the G-crossed braided tensor category formalism and the occasions when it is acceptable to disregard it. We then define a hi-erarchy of "extended Ising models" in section 3 and discuss the general fusion and braiding relations for these models in sections 4 and 5 respectively. In section 6 we clarify the correspondence between these relations and the possible logical operations that can be performed using these anyons. Finally, in section 7 we discuss how general models in this hierarchy can be realised in stacks of 2d surface codes with special attention given to the case of the 2d colour code.

Anyons and Twists
In this section we briefly review the theoretical background necessary for the rest of the paper. We assume that readers are familiar with topological stabiliser codes and so our focus is on providing an outline of anyon fusion and braiding relations in sections 2.1 and 2.2. We refer readers who are not familiar with these codes to references [9,10,11,12,13]. In section 2.3 we present a similar discussion regarding twists in topological codes and briefly touch on the category theory formalism that describes these objects.

Fusion and Braiding
There exist a wide variety of ways to describe anyon models. They can be described in terms of topological charges [7], unitary braided tensor categories [8] and through the lens of conformal field theory [14] [15]. For our purposes the topological charge description is largely sufficient, although we will very briefly use the categorytheoretic approach when discussing twists in section 2.3. This means that we have some finite number of anyonic species, each possessing a unique label and topological charge. The anyons obey a set of fusion and braiding relations. Fusion relations are generally written in the form [16] a × b = c N c ab c (1) where N c ab is an integer counting the number of ways anyons a and b can fuse into c. For a given charge a if for any charge b N c ab is non-zero for at most a single charge c, in other words the fusion of a with any other anyon has only one possible result, then we say that a is Abelian. Otherwise it is non-Abelian. All anyon models must contain a unique vacuum charge 1 such that a × 1 = a. Additionally, each charge a in the model must have a unique inverseā with which it can fuse to the vacuum in a unique way (N 1 aā = 1). The total anyonic charge within a given region is a topological invariant and so cannot be altered by operations within this region. However, through alterations to anyon fusion order and intermediate fusion outcomes we can arrive at this same total charge via different paths, called fusion channels. This gives rise to the notion of an anyonic fusion space with dimension equal to the number of possible fusion channels. The quantum dimension of an anyon is defined to be meaning that d a = 1 for all Abelian anyons and d a > 1 for non-Abelian anyons. The dimension of the fusion space of N anyons a grows asymptoticlly as (d a ) N in the limit of large N . Clearly if we wish to use anyons in quantum computation then only those models which contain non-Abelian anyons are of interest to us. Changes of basis in the fusion space can be described using F-moves where u (v) can be any of the fusion outcomes of a and b (b and c). More generally we should include additional indices describing the precise fusion channel by which a and b (b and c) fuse to u (v) but in what follows we will only consider fusion rules with N c ab equal to either 1 or 0 and so such generality is unnecessary. The Fmatrices associated with an anyon model can be found from the fusion rules by solving the pentagon equation [16] which can be understood diagramatically in Fig.  1.
Logical operations on the fusion space can be performed via the braiding of anyon pairs. This is an operation which fuse to c. R c ab will be a phase if a and b have only a single possible fusion outcome or a diagonal matrix indexed by c if there are multiple possible outcomes. The charges a and b and the outcome of their fusion c are left unchanged by the braiding in accordance with the fact that anyonic charges cannot be modified through local operations. However, the fusion outcome of a or b with a third anyon may be modified by this braid. If the F-matrices for a model are known then we can find the R-matrices for that model using the hexagon equation As with the pentagon equation, the hexagon equation can more easily be understood when presented diagramatically as in Fig. 2.

Examples
Two anyon models of central importance in our work are the quantum double of Z 2 [7] and the Ising anyons [16]. The former is an Abelian model with charges 1, e, m and , fusion rules and braiding relations This model describes the excitations that arise in the toric code, with e and m anyons corresponding to X and Z errors and corresponding to a combination of the two (i.e. a Y error). It also possesses a symmetry: we can exchange the e and m charge labels without affecting any of the fusion or braiding relations.
In contrast, the Ising anyon model is non-Abelian. It contains three charges: 1, ψ and σ. The Ising anyon fusion rules are and the F matrix for the fusion of three σs is

Twists in Topological Codes
We noted above that the quantum double of Z 2 is symmetric under exchange of e and m charges. We can consider a domain wall which applies precisely this symmetry, achievable in the toric code via a line of modified stabilisers each containing two Z and two X operators as shown in Fig. 3 [5] [7]. An X error moved across such a domain wall will be transformed to a Z error and vice versa. The end points of domain walls such as this are called twists and are formally described by G-crossed braided tensor categories [8]. We now give a very brief outline of some of the basic ideas of this formalism. This will be limited to the minimum details required for drawing a connection between twists and anyons and readers interested in a rigourous mathematical description of this formalism should refer to the sources cited. An anyon model can be described by a unitary braided tensor category C 0 which has charges a 0 and a (possibly trivial) symmetry group G. The elements of G are labelled g and correspond to the symmetries of the anyon model. The identity a) element of this group is labelled 0. The action of g on C 0 is an invertible map from C 0 to itself. In a physical realisation of this anyon model each g will correspond to a twist and braiding an anyon around this twist will apply the symmetry g to that anyon, with this action denoted as g a 0 . The topological charge of a twist can be measured by braiding it with an anyon a 0 which is invariant under the symmetry g as in Fig. 4 a). For each symmetry g we have a new category C g which has charges a g . The number of distinct a g is equal to the number of g-invariant charges in C 0 . When fusing charges a g and b h we must have that a g × b h = c g·h where g·h is the composition of elements of G. This is called G-graded fusion.
Since g · 0 = g we have that a g × b 0 = a g . If the g-invariant charge(s) in C 0 cannot distinguish b 0 from the vacuum then a g = a g and we say that anyon b 0 is "localised" by twist a g . Such charges and twists have fusion/splitting rules such that •ā g also localises b 0 .
• b 0 is one of the possible fusion outcomes of a g ×ā g .
• All charges in the orbit of b 0 under the action of g are also localised by a g andā g .
Additionally we note that the set of localisable charges for a particular twist must be closed under fusion since if a 0 and b 0 both braid trivially with the g-invariant charges in C 0 then so must the result of their fusion.
Braiding of charges a g and b h involves the action of the relevant symmetries such that charges can be modified by these braids (this is the Gcrossed braiding part of the formalism). However, we will not require this part of the theory for reasons that will become clear below.
Finally, we note that fusion rules in this formalism must still satisfy the pentagon equation. Braiding rules are generalised to follow a "heptagon equation" which accounts for the fact that braiding with twists can alter charge labels.
We return now to the previously discussed case of twists in the toric code. C 0 in this case is the quantum double of Z 2 which has only one symmetry so G has only two elements, 0 and g, where 0 is the identity. C 0 contains two g-invariant charges (1 and ) so C g has two charges which can be distinguished by braiding with . More specifically there is one charge corresponding to the fusion of the twist with either 1 or and one charge corresponding to fusion with e or m. This twist possesses two significant features: (1) it is selfinverse and (2) its associated invariant charges are also its localisable charges. This means that the subset of charges consisting of this twist and its localisable charges is closed under fusion. Furthermore, none of the charges in this subset can be altered by braiding with any of the others. In other words this subset functions as an anyon model -specifically the Ising anyon model. This is precisely what was noticed by H. Bombin in [7], although the argument in that paper was formulated in terms of topological string operators.
In the 2d colour code the situation is not quite so simple. In [6] the 72 twists of the colour code are identified and arranged into nine conjugacy classes. The authors of [6] point out that the action of these twists can best be understood by considering the nine (non-trivial) bosonic anyons of the colour code arranged as in Table. 1. These anyons are all self-inverse and Abelian. Two anyons in a row or column fuse to the third and braid trivially with each other. Two anyons which do not share a row or column fuse to a fermion and aquire a phase of -1 under full exchange (monodromy).
The symmetries of the colour code anyon model are the permutations of this table which preserve the rows and columns. These permutations are the column permutations (6 options), row permutations (6 options) and the transpose (2 options) giving 6 × 6 × 2 = 72 possible symmetries. Twists belonging to three of the nine conjugacy rx gx bx ry gy by rz gz bz Table 1: The nine non-trivial bosonic anyons of the 2d colour code arranged as in [6]. rx implies an X error on a red plaquette and so on.
classes possess the same properties as the surface code twist described above: they are self-inverse and their associated sets of invariant and localisable charges are equivalent. One of the other six classes is trivial (it contains only the identity twist) and twists in the other five classes possess neither of these properties.
In what follows we consider the general case of the anyon model associated with these self-inverse twists and their localisable charges. In section 7 we will show that twists in any number of stacked surface codes can only have an invariant set of localisable charges if they are self-inverse.

A Hierarchy of Models
Recall the standard Ising model in which we have a single non-Abelian anyon σ and two Abelian anyons 1 and ψ such that σ × σ = 1 + ψ. We can extend this model by including additional Abelian anyons and modifying the outcome of σ × σ to include these anyons. Such extended models already exist in the literature in the form of parafermions, for example in [17], but the Abelian anyons in these models are not generally self-inverse and so cannot be the natural anyons of the colour code. If we write the Abelian anyons of an extended Ising model as α i (where α 0 = 1) and the non-Abelian anyon as β and require that each α i must be its own antiparticle then we obtain the following fusion relations . These are exactly the fusion relations we observe for self-inverse twists and their localisable anyons in the colour code. For example, the colour code twist B is self-inverse and can localise the anyons bx, by, bz from table 1 (as well as the vacuum anyon). The set {1, bx, by, bz} is closed under fusion and so if we write 1 = α 0 , bx = α 1 , by = α 2 , bz = α 3 and B = β then fusion relations of these five charges exactly match (11).
Only specific values of n yield valid extended Ising models. For example, the model {α 0 , α 1 , α 2 , β} is not valid because it is not closed under fusion (α 1 ×α 2 must have a fusion outcome α k where k = 0, 1, 2). Given a valid extended Ising model containing m αs we can find the next valid model with n > m by adding a single new charge α m+1 to the model, fusing α m+1 with all existing charges, and adding all fusion outcomes to the model. If we write these fusion outcomes as α i × α m+1 = α i+m+1 then we can see that the resulting model must be closed under fusion since in another α k with k ≤ m since the initial model was closed under fusion.
• Fusion of α m+1 with any anyon in the model results in another anyon in the model due to the above procedure.
• Fusion of any α i , α j with i ≤ m and j ≥ m + 2 can be written as α i≤m × α j≥m+2 = α i≤m ×α k≤m ×α m+1 by definition of α j≥m+2 , and this is in the model due to the above two points.
Thus we can inductively define all extended Ising models beginning from the standard Ising model: {α 0 , α 1 , β}. We can label these models by I k where k = 1 is the standard Ising model. The number of α i in a given I k is n k = 2n k−1 = 2 k−1 n 1 and n 1 = 2 so n k = 2 k . The β anyon for each I k can be written as β k , and it has quantum dimension √ 2 k . Note that we can equivalently define these models from multiple copies of Z 2 . The Abelian charges of the standard Ising anyon model form a group (with composition of group elements described by the model's fusion rules) that is isomorphic to Z 2 . Similarly the Abelian charges of I 2 form a group that is isomorphic to Z 2 × Z 2 and so on. In general the group of Abelian charges of I k will be isomorphic to k copies of Z 2 . A set consisting of a finite group A and a single additional element β with composition of these elements defined as in (11) is called a Tambara-Yamagami (TY) category [18]. Thus our "extended Ising hierarchy" can be described more formally as a hierarchy of TY categories with categories in the k th level of the hierarchy having base group (Z 2 ) k .

F Matrices
In this section we show the possible F matrices F β βββ for general β k . The full derviation of these matrices can be found in appendix A. These matrices (up to symmetry-preserving row and column permutations) are found to be a subset of the Hadamard matrices called Sylvester or Walsh matrices [19] multiplied by a constant. The implications of this fact on the trace of these matrices is examined as this will be relevant in section 5.
A general F-matrix for extended Ising anyons can be found from the pentagon equation (3). Every F d abc is equal to zero (if it is disallowed by the fusion rules) or a phase except for F β βββ which involves only non-Abelian anyons. In appendix A we use the pentagon equation to show that up to a choice of gauge this matrix has the form where φ is a symmetric Hadamard matrix with elements given by Hadamard matrices are n × n matrices where all entries are ±1 and M M T = nI [20]. Additionally, φ as we have defined it has the property that all entries in the first row and column are +1 (such Hadamard matrices are called "normalised" but we will avoid using this term to prevent confusion with its more common usage in quantum physics). We can further restrict φ to a subset of Hadamard matrices by using the group theory and TY category connections discussed previously. The F matrices of TY categories are related to symmetric non-degenerate bicharacters where "character" refers to multiplicative character, i.e. a group homomorphism from a group A to the multiplicative group of a field F × [21]. A bicharacter is then a function χ : A × A → F × which satisfies [22] χ(a 1 a 2 , a 3 ) = χ(a 1 , a 3 )χ(a 2 , a 3 ) and χ(a 1 , a 2 a 3 ) = χ(a 1 , a 2 )χ(a 1 , a 3 ).

(14)
We can see that φ satisfies this by considering the pentagon equation with 1 = 3 = β, 2 = α i , 4 = α j and 5 = α k which yields a constraint (15) which can be rewritten as using (13) and the fact that Thus φ is a bicharacter on (Z 2 ) k . The above definition means that if we fix one argument of a bicharacter on a group A then the bicharacter as a function of the other argument defines a character on A. In other words, each row and column of φ is a character on (Z 2 ) k . Because they are homomorphisms the action of these characters on (Z 2 ) k is defined by their action on each of the k copies of Z 2 which make it up. There are two valid ±1 valued characters on Z 2 which are given by the rows/columns of the 2×2 Hadamard matrix and coincide with the irreducible representations of Z 2 . Thus for (Z 2 ) k there are 2 k possible characters corresponding to the rows/columns of H ⊗k 1 . These matrices are a subset of the Hadamard matrices called Sylvester or Walsh matrices [19]. The possible bicharacters for (Z 2 ) k are then the Sylvester matrix H ⊗k 1 and any other matrices which can be obtained from this matrix via symmetry-preserving row/column permutations.
H 1 is the unique bicharacter for k = 1. For k = 2 there are four possible bicharacters which we can write in matrix form as One of these matrices has trace 4 while the other three have trace 0. These three all correspond to the same anyon model up to relabelling of charges. In appendix B we show that symmetry-preserving permutations of columns of a symmetric Hadamard matrix must alter the trace of the matrix by either 0 or ±2 k , with the latter only being possible for even k. The Sylvester matrices all have trace 0 and since one row/column contains only +1s we cannot have Tr(φ) = −2 k so the only possible traces are 0 for odd k and 0, 2 k for even k. In comparison general symmetric 2 k × 2 k Hadamard matrices must also have trace 0 for odd k but the trace can take any value 2 k − 2 k/2+1 ≥ 0 for even k [23].

R Matrices
We now discuss the possible R matrices for extended Ising models which are found from the hexagon equation (4). This equation was solved for the case of TY categories in [24]. Once again a similar solution in the language of section 2 can be found in appendix C.
We are concerned only with R ββ which are diagonal matrices. From the hexagon equation we obtain the constraint with R α 0 ββ having possible values {±1, ±i, ±e iπ/4 , ±e −iπ/4 } for even k and {±e iπ/8 , ±e −iπ/8 , ±ie iπ/8 , ±ie −iπ/8 } for odd k. This equation tells us that (up to a global phase) the non-zero elements of R are equal to ±1 and ±i, with the number of real (imaginary) elements equal to the number of +1s (−1s) in the diagonal of φ.
The hexagon equation also gives a constraint on the trace of R, which tells us about the relative numbers of positive and negative diagonal elements. We can use (20) to rewrite this as where a and b are integers. In order for |R α 0 ββ | 2 = 1 we need a 2 + b 2 = 2 k . For even k the solutions k Tr(φ) Number of ±1 Number of ±i to this equation are a = ±2 k/2 , b = 0 and a = 0, b = ±2 k/2 and for odd k they are a = ±b = ±2 (k−1)/2 (see appendix C). Combined with (20) this completely determines the number of +1, −1, +i and −i elements on the diagonal of R ββ for a given k and Tr(φ) up to global phase, with the number of each shown for each case in table 2.
We also show in appendix C that the elements φ ii describe the exchange statistics of the anyons α i , with φ ii = 1 indicating that α i is bosonic and φ ii = −1 telling us that α i is fermionic. Thus we expect that for both odd and even k we can obtain models containing 2 k−1 bosonic and 2 k−1 fermionic Abelian charges and a single non-Abelian charge. We additionally expect that for even k we can obtain models containing 2 k bosonic Abelian charges, no fermionic charges and a single non-Abelian charge. The models with 2 k−1 bosonic and 2 k−1 fermionic charges can be viewed as "sub-models" of k copies of the standard Ising model (e.g. in some kind of multi-layer system) containing all Abelian charges and only a single non-Abelian charge (namely, the charge corresponding to σ 1 ⊗ σ 2 ... ⊗ σ k ). We note also that this "sub-model" simply corresponds to the case where we neglect some of the charges in the original model and does not mean that these charges are no longer present. It is therefore different from procedures such as that of Bais and Slingerland [14] in which an actual change to the model is made.
The models containing 2 k bosonic charges cannot be produced from copies of the standard Ising model and instead correspond to copies of a different anyonic model with four bosonic Abelian charges and a single non-Abelian charge.

Logical Gates
In this section we will rephrase the findings from the past two sections in terms of the possible logical gates which we can perform using these anyons.
Consider a specific anyon model containing 2 k Abelian charges and a single non-Abelian charge. This model has an F matrix F β βββ associated with changing the fusion order of three of the non-Abelian anyons and an R matrix R ββ associated with braiding two of the non-Abelian anyons. A system containing four such anyons has a 2 k dimensional fusion space for which the 2 k Abelian anyons of the model form a canonical basis. The F and R matrices provide us with two logical operations which can be performed on this space. The F matrix is a mapping between the canonical basis and a basis of equal superpositions of the canonical basis vectors. The R matrix selectively applies one of the phases {+1, −1, +i, −i} to each vector of the canonical basis with the total number of +1s, −1s, +is and −is applied consistent with table 2.
Both of these operations may be interpreted in terms of Clifford gates on k qubits: the F matrix has the same rows and columns as a tensor product of k Hadamard gates (up to a global phase) and the same is true for the R matrices and tensor products of either k S gates or k/2 CZ gates. Notice that our canonical basis vectors are currently labelled only by anyonic charges and we have not yet defined an encoding for qubits in this space. We can always choose this encoding such that the ordering of the diagonal elements of R in the computational basis is consistent with the respective tensor product of Clifford gates. The same is also true for the trace-0 F matrices but not for the trace-√ 2 k F matrices since the trace is independent of our choice of encoding. These matrices cannot be decomposed into a tensor product of single qubit gates and instead correspond to tensor products of the trace-4 matrix in (19) multiplied by 1/2. This matrix is equivalent to SWAP·(H ⊗ H) and so is also Clifford. Thus, up to global phases and a choice of encoding, all F and R matrices implement Clifford operations on our Hilbert space.

Stacked Surface Codes
Part of our motivation for obtaining the results presented in the previous sections was to examine the braiding relations of twists in the 2d colour code. In this section we will see that we can obtain twists belonging to the first and second levels of the hierarchy in this code and in general we can obtain twists belonging to the k th level of the hierarchy in a stack of k 2d surface codes. When visualising such a stack we might imagine multiple copies of the surface code placed one above the other, but we allow operations between any two layers regardless of how "far apart" in the stack they are and thus there is no notion of distance in a third dimension and the entire stack is embeddable in 2d (although this will have an impact on the locality of the stabilisers).
The fact that models from the second level of the hierarchy can be realised in the 2d colour code is readily apparent from [6]. For example the twist B which exchanges the r and g columns of Table 1 has four nontrivial localisable charges: In general if we consider equivalence up to global phases and choose our qubit encoding as discussed in the previous section then there are two R-matrices for k = 2: Recall that the first of these matrices belongs to a model with four bosonic Abelian charges while the second belongs to a model with two bosonic and two fermionic charges. In the notation of Kesselring et al [6] the models containing four bosons are those associated with a twist from conjugacy class B while those containing two bosons and two fermions are associated with twists from conjugacy class C. The twists in conjugacy class G have only two localisable charges and correspond to models from the first level of the hierarchy.
So far we have seen that we can realise the k = 1 level of the extended Ising hierarchy in the surface code, and the k = 2 level in the colour code, which is equivalent to two copies of the sur-face code [25]. We now consider the general case of a stack of k copies of the surface code.
Consider the anyon model of k stacked surface codes (in the absence of twists). The topological charges in this model are the elements of a finitely-generated free group, whose generating set can be written {e 1 , m 1 , e 2 , m 2 , ..., e k , m k } where the subscript shows the layer in the stack which the charge belongs to. Twists in the code stack correspond to symmetries of the anyon model. These symmetries can be formally defined as the elements of the automorphism group of the anyon model which preserve braiding relations. The action of these symmetries can be described via a set of orbits, each of which can be written as with the "trivial orbit" defined as We first show that only self-inverse symmetries g of this anyon model can have a g-invariant set of localisable charges. Consider a twist t g and a charge b 0 = a 0 × g a 0 . b 0 can be localised by t g because we can split it into a 0 and g a 0 , then braid a 0 around t g and fuse it to the vacuum with g a 0 . If we braid b 0 around t g we obtain g b 0 = g a 0 × g 2 a 0 and so in order for b 0 to be g-invariant we must have g 2 a 0 = a 0 , implying that g is self-inverse.
From this we can see that if a twist in a stack of surface codes (together with its set of localisable charges) can be considered as an anyon model this model will belong to the hierarchy of extended Ising models defined in section 3.
All non-trivial orbits associated with a selfinverse symmetry have the form a → b → a which can also be written a ↔ b.
The full automorphism group of a finitely generated free group with ordered basis [x 1 , ..., x n ] can be generated by the elementary Neilsen transformations [26]: The second transformation is equal to the identity transformation in our case because all charges in our model are their own inverse. We thus consider only the first and third transformations, but not all applications of these transformations are valid because we must also preserve braiding relations. In order to do this we require that if we map x i → x j then we must also map x i → x j and if we map x i → x i x j then we must map x j → x i x j , where x i can be either e i or m i and x i x i = i . We also cannont map e i → e i m i within a layer because this exchanges a boson with a fermion. In other words all symmetries of the model can be generated by the transformations e i ↔ m i (26) e i ↔ e j and m i ↔ m j (27) e i ↔ e i e j and m j ↔ m i m j (28) which are simply the generators of all colour code symmetries generalised to act on a stack of more than two surface codes [6]. A simple way to obtain a twist corresponding to a β k anyon is simply to combine twists associated with symmetry (26) on k different levels. The domain walls produced by these symmetries in the code correspond respectively to lines of H, SWAP and CNOT gates applied in the code stack. Since SWAP can be generated from CNOTs the generating set of symmetries can be reduced to just (26) and (28). This is consistent with the set of generating symmetries identified in [3] although we arrive at this result by a different method. Any product of these symmetries thus corresponds to a product of Clifford gates in the code stack and so the code containing the twists will also be a 2d stabiliser code. Braiding operations are performed using predefined sets of single-qubit Pauli measurements and additional modifications to stabilisers by the Clifford gates listed above [5] and such operations in a 2d stabiliser code should not result in a logical non-Clifford gate. Thus all twists produced by composition of these symmetries should have Clifford braiding relations. This result is valid for more than just selfinverse twists since (26)(27)(28) are the generators of all symmetries of the anyon model. Thus the restriction to Clifford braiding operations is valid for all twists in stacked surface codes. This is in agreement with recent results regarding the power of defect braiding in topological codes [27] [28].
Finally we comment briefly on the faulttolerance of such braiding procedures. As mentioned above, braiding operations with twists can be performed using the standard code deformation techniques of (1) measurement of modified stabilisers and (2) single-qubit Pauli measurements to remove physical qubits from the code and provide information for decoding [5]. In a code with local stabilisers these operations will also be local so we expect that braids with these generalised twists should remain fault-tolerant under existing decoding procedures.

Summary
We have constructed a hierarchy of anyonic models extending the standard Ising anyon model and identified the significant properties of the F and R matrices for these models in the general case. These models are specific cases of Tambara-Yamagami categories and have F and R matrices belonging to the Clifford group. These anyon models can be realised using Abelian anyons and twists in stacked surface codes: given a stack of k surface codes we can realise models belonging to level k of the hierarchy. The restriction to Clifford group braiding relations of twists extends even to those twists which do not reproduce the behaviour of anyons. This result is consistent with other recent results in this area.
A number of possible future research directions exist following the results outlined above. Although we have shown that twists in stacked surface codes will always have Clifford braiding relations we have only characterised these relations for a small subset of these twists. Finding the braiding relations for the remaining twists will likely require the use of the full G-crossed braided tensor category formalism described in [8].
In topological codes of dimension greater than 2 the braiding relations of defects and excitations are very poorly understood. Braiding in general is a more complicated concept in these higher dimensions and must be performed with non-pointlike objects to be non-trivial. Recent work shows that defects can be constructed in higher dimensions that reproduce the braiding relations of Ising anyons [28] and it is also known that using domain walls and puncture encodings we can perform a braided version of any singlequbit transversal gate within a code. A rigorous examination of such braiding schemes should be possible using the full G-crossed category formalism since puncture encodings utilise the fact that some of the natural anyons of the code can condense at the puncture's boundary. Passing these punctures through a domain wall should therefore be equivalent to braiding these natural anyons with a twist and modifying the anyon's charge in the process. It remains to be seen whether or not it is possible to implement a braided non-Clifford gate that is not also a transversal gate of the code.
The results of this work could also be extended to qudit codes, provided the anyon models of these codes possess symmetries. The twists that arise in such cases may again allow us to obtain non-Abelian anyons in codes where all the natural anyons are Abelian, but in cases where the natural anyons are non-Abelian but non-universal there is the possibility that additional non-Abelian charges may be obtained which make these models universal. These models may therefore be a fruitful topic for future research.

Appendices A Derivation of F Matrices
Since we will be using the pentagon equation (3) extensively in this section we quote it again here (29) β is non-Abelian so (F β βββ ) is a matrix which we will henceforth refer to as F . If we let 1 = α 0 then the elements of F are F ij = (F β βββ ) α j α i . All other (F s pqr ) j i are phases. We note some important properties of these phases and their inverses: 1. Every phase (F s pqr ) j i has an inverse (F s rqp ) i j since changes of fusion order are always reversible.
2. (F s pqr ) j i where p,q or r = 1 are trivial reorderings and correspond to a phase of 1.
3. Phases of the type (F s pqp ) i i are self-inverse and so have value either 1 or −1.
We begin with the case where 1, 2, 3, 4 = β, 5 = 1, b = d = β, a = α i and c = α j . We obtain constraints In the case that i = j the LHS of this equation is equal to 1 by property 1 as listed above. If i = j the LHS is equal to 0 since α i × α j = 1 is only possible when i = j. If we write (F 1 βαxβ ) β β = θ x and sum over repeated indices then we can rewrite (30) as where each element F ij = θ i · F ij (we do not sum over the repeated index here). The LHS of (31) is the matrix F and its inverse. F is unitary so we must have that By property 3 θ i = ±1 and by property 2 θ 0 = 1. Thus and and we can combine these to show Thus if θ x = −1 we must have F 0x = F x0 = 0 for all x ≥ 0.
The next configurations of the pentagon equation that we consider are 1 = α x , 2, 3, 4 = β and 2 = α x , 1, 3, 4 = β. The first of these yields constraints The second configuration yields From (36) we have that and from (37) Setting x = i and x = j we can combine these equations to obtain 40) Thus all elements F ij have magnitude equivalent to that of F 00 and so the only way for θ x = −1 is to have F 00 = 0 but this contradicts the constraint that F ix F xi = 1 since all the terms in this sum would be 0. Thus all θ x = 1 and F must be Hermitian. Additionally, the magnitude of all elements in the matrix must be 1/ √ 2 k and since the diagonal elements must be real F 00 = ±1/ √ 2 k . Setting we can rewrite (40) as By property 3 all φ ij = ±1 (with φ i0 = φ 0j = 1 by property 2) while f i0 can be any phase. We also have that f 0j = (f j0 ) * which can be verified by considering the five-anyon fusion tree Equating these and eliminating trivial terms using property 2 we get which can be rewritten as using property 1.
We can then write F as the Hadamard product of a matrix of f s (f ) and a matrix of φs (φ) multiplied by 1/ f is Hermitian and so φ must also be Hermitian. The multiplication of F by itself gives where we have used the fact that f 0x f x0 = 1. The final equivalence implies where I is the 2 k × 2 k identity matrix. The diagonal elements of f are all 1 (since f 0i = (f i0 ) * ) so we must have that φ 2 = 2 k I. Thus φ is a symmetric Hadamard matrix [20]. φ has only a finite number of solutions while f has an infinite number. The obvious interpretation of these two matrices is that different φ correspond to different anyon models (with the discreteness of these solutions consistent with Ocneanu rigidity [29]), while the different f correspond to a choice of gauge. We can see that we cannot transform between solutions of φ by changes to f by observing that f is completely characterised by the values of f i0 , whereas φ i0 are always 1, so changes to f are always reflected in the first row of F β βββ while changes to φ are not. We also show in appendix C that the braiding matrices of the models depend only on φ and not on f . Thus we can make the gauge choice that all f = 1 so Note that instead of (38) and (39) we could have obtained from (36) and from (37). Setting x = i and x = j and combining these as before gives (51) Using property 1 we can see that the first three terms are equal to (f j0 f 0i ) −1 which is equal to 1 due to our choice of gauge. Thus we have that

B Column Permutations of Symmetric Hadamard Matrices
In this appendix we show that column permutations that preserve the symmetry of a symmetric 2 k × 2 k Hadamard matrix must alter the trace by either 0 or ±2 k . Consider swapping the second and third columns of the following matrix from (19). This matrix corresponds to H 1 ⊗ H 1 and this column swap will result in the matrix from (19) with trace 4.
We have divided this pair of columns into three sections. Any changes to the top or bottom (offdiagonal) sections will result in a non-symmetric matrix but the central (diagonal) section can be modified while preserving symmetry. Because these two columns are identical in the off-diagonal sections and the diagonal section is symmetric both before and after the exchange the symmetry of the overall matrix is preserved. However, for k > 2 such exchanges cannot be symmetry preserving because the diagonal section will always contain 2 elements from each column while the off-diagonal sections will contain the other 2 k − 4. In order for the columns of the matrix to be orthogonal each pair of columns must match in exactly half their entries so for matrices larger than 4 × 4 it is impossible to exchange two columns in such a way that the off-diagonal sections are unchanged. Instead, we must exchange sets of columns. Consider a general 2 k × 2 k matrix broken into 2 k−2 × 2 k−2 blocks as follows where we have once again marked diagonal and off-diagonal sections.
Lemma 1: The exchange of EACG and F BDH can only preserve the symmetry of the matrix if the trace of the diagonal section is negated by the exchange We can only exchange EACG and F BDH if E = F and G = H. Additionally we must have that A = −B and C = −D or the columns of the matrix would not be orthogonal. Thus the trace of the diagonal section is negated by this exchange.
Lemma 2: The exchange of EACG and F BDH can only preserve the symmetry of the matrix if A is a symmetric Hadamard matrix For the overall matrix to be symmetric both before and after the exchange we require that B T = C and A T = D and therefore A = −B = −C T = D T and the entire diagonal section is determined by A.
A is must be symmetric for the overall matrix to be symmetric.
To show that the columns of A are orthogonal we consider two columns, c 1 and c 2 , of our full 2 k × 2 k matrix such that both columns pass through A. We note that if this pair of columns have m matched elements within A then they must have 2m matched elements within the diagonal section and 2 k−1 − 2m matched elements in the off-diagonal sections (since exactly half the elements of each column must match). We now consider the pair c 1 and −c 2 where −c 2 is the column passing through B that is equal to −1 × c 2 within the diagonal section and identical to c 2 outside of this section (i.e. the column we wish to exchange c 2 with). c 1 and −c 2 have 2 k−1 − 2m matched elements within the diagonal section and so must have 2m matched elements in the off-diagonal sections. Since c 2 and −c 2 are identical in the off-diagonal sections we have that 2m = 2 k−1 − 2m and m = 2 k−3 . Thus for each pair of columns in A half of the elements are matched and the other half are unmatched and the columns of A are orthogonal.
Since A is a symmetric Hadamard matrix of order 2 k−2 it can be partitioned into blocks such that A is a symmetric Hadamard matrix of order 2 k−4 . This process can be repeated recursively, eventually terminating when we arrive at a matrix of order 1 (for odd k) or order 2 (for even k). The possible symmetric Hadamard matrices with these orders are (up to a possible global phase of -1) (18) and (19). These matrices are restricted to have trace=0 and trace=0 or 4 respectively, and going back up the chain of recursion we see that Tr(A) = 0 for odd k and Tr(A) = 0, ±2 k−1 for even k. The trace of the central section is then 2Tr(A) and by Lemma 1 a symmetry preserving exchange of columns must negate this trace, changing the trace of the overall matrix by either 0 for odd k or 0, 2 k for even k.
So far we have only considered swapping columns in a very specific arrangement, but any desired set of column swaps can be rewritten in this form such that it is apparent that the same constraints apply. This is achieved by applying column permutations such that the columns we wish to swap are correctly arranged at the centre of the matrix and then applying a matching set of row permutations (since matching column and row permutations preserve the symmetry of the matrix). Following the exchange of column blocks as described above we apply the same set of row and column permutations again. By considering how permutations transform under conjugation and the fact that row and column permutations commute we can see that this operation is equivalent to exchanging the columns without first moving them to the center.

C Derivation of R Matrices
The R-matrix is defined from the hexagon equation When we set 1, 2, 3, 4 = β, a, c = α i this gives: and from (42) Setting 1 = α i and 4 = α j we instead obtain So we see that the braiding relations are dependent on φ but not on f as we would expect.
Consider the case in (56) where i = j. We have that The RHS is independent of i so Additionally, if we set i = 0 and sum over j all terms on the RHS cancel except for those where b = 0 which sum to 2 k R β βα 0 (since all rows/columns of φ except for the first contain an equal number of 1s and −1s), so since R β βα 0 just describes braiding with the vacuum and is therefore trivial. Using (59) we can rewrite this as The ± in the sum do not all need to be the same, but they must be chosen such that |R α 0 ββ | 2 = 1 or the matrix R ββ would not be unitary.
The number of +1 (−1) terms in the diagonal of φ tells us the number of ±1 (±i) terms in the sum of (62) which we can rewrite as In order for |R α 0 ββ | 2 = 1 we require that a 2 + b 2 = 2 k . For even k this means that either a = ±2 k/2 , b = 0 or a = 0, b = ±2 k/2 . For odd k we have that a = ±b = ±2 (k−1)/2 . The proofs are as follows: Consider a right-angled triangle with sides a ≤ b < c opposed by angles A, B, C.
• Even: √ 2 k = 2 k/2 is an integer so from a 2 + b 2 = 2 k we assume the existence of a Pythagorean triple (|a|, |b|, 2 k/2 ). a and b must have the same parity for the sum of their squares to be even. If they are both odd then this triple is primitive since |a| and |b| have no even factors and 2 k/2 has no odd factors. If they are both even then there must be some associated primitive triple (|a|/2 x , |b|/2 x , 2 k/2−x ) where the first two elements are both odd. All primitive triples can be constructed using Euclid's formula a = m 2 − n 2 , b = 2mn, c = m 2 + n 2 (64) where m and n are a pair of coprime integers, one of which is even. However, this means that c is odd, giving a contradiction. Thus this primitive triple does not exist and neither does the triple (|a|, |b|, 2 k/2 ). The only remaining solutions to the equation a 2 + b 2 = 2 k are a = ±2 k/2 , b = 0 and a = 0, b = ±2 k/2 • Odd: √ 2 k = 2 k/2 = 2 (k−1)/2 √ 2 where 2 (k−1)/2 is an integer power of two. Given c = 2 (k−1)/2 √ 2 we have that a = 2 (k−1)/2 √ 2sin(A) and b = 2 (k−1)/2 √ 2sin(B). We require that both a and b are integers and thus sin(A) and sin(B) must both be integer multiples of 1/ for odd k.
Finally, we note that by setting i = j in (57) we can show that where R α i α i 1 = 1 since braiding with the vacuum is trivial. Using this and instead setting j = 0 we find In other words the elements φ ii tell us the selfexchange statistics of the charges α i .