Thermalization, Error-Correction, and Memory Lifetime for Ising Anyon Systems

We consider two-dimensional lattice models that support Ising anyonic excitations and are coupled to a thermal bath. We propose a phenomenological model for the resulting short-time dynamics that includes pair-creation, hopping, braiding, and fusion of anyons. By explicitly constructing topological quantum error-correcting codes for this class of system, we use our thermalization model to estimate the lifetime of the quantum information stored in the encoded spaces. To decode and correct errors in these codes, we adapt several existing topological decoders to the non-Abelian setting. We perform large-scale numerical simulations of these two-dimensional Ising anyon systems and find that the thresholds of these models range between 13% to 25%. To our knowledge, these are the first numerical threshold estimates for quantum codes without explicit additive structure.


I. INTRODUCTION
One of the most interesting features of two-dimensional quantum systems with anyonic excitations [1] is their application to quantum information science, where the topological nature of the anyons offers some intrinsic protection of quantum coherence [2].These topologically ordered systems are insensitive to local perturbations [3][4][5], have a ground-space degeneracy that is a function of the topology of the system [6,7] and can realize quantum gates by braiding anyons, which are naturally robust to error due to their topological nature [2,8].Many common anyon models such as the toric code [2] and the color codes [9] have an Abelian structure that simplifies their analysis greatly.However, Abelian anyon models are very restricted when compared to the general class of anyonic systems.Notably, Abelian anyons lack the capacity to perform quantum computation by braiding of anyons alone.By contrast, non-Abelian anyon systems have excitations that can be used to implement topological quantum computation [10,11], and so they are of interest for quantum information processing schemes as well as quantum error correction.It is also known that the dynamics of Abelian anyons can significantly differ from those of non-Abelian anyons in some circumstances [12][13][14][15].
A particular class of non-Abelian anyonic excitations called Ising anyons are thought to be among the most promising candidates for experimental realization of topological quantum computing.Ising anyons are experimentally motivated [16] as the expected excitations of the ν = 5 2 fractional quantum Hall states [17,18], and also appear as the excitations of several lattice spin models [19][20][21].In particular, Ising anyons can robustly perform Clifford gates via braiding and have a high noise threshold for universal quantum computation when coupled with magic state distillation [22,23].
These general features of anyon models, and of Ising anyons in particular, have attracted a lot of attention during the past decade.This picture of an intrinsically robust topological quantum computer is spoiled at finite temperature, however.Indeed, if left unattended, thermal excitations diffusing in the system can corrupt the quantum information stored in the ground space.While these excitations are suppressed by a mass gap (ideally) much larger than the temperature, they nonetheless appear in a finite density at any non-zero temperature, and for scalability issues it is essential to devise a scheme that can cope with their existence.One approach is to engineer a mass gap that grows with the system size, leading to a self-correcting behavior of the memory.However, there is accumulating evidence that this cannot be realized in two spatial dimensions [24][25][26][27].Additionally, even a self-correcting system in contact with a thermal bath will accumulate a finite density of excitations over time, and it is not clear that the quantum information stored in the ground state can be manipulated or readout without first cooling back to the correct ground state [28].
Here we pursue a different approach that can be seen as a viable alternative to self-correction as well as a reliable method for readout of noisy stored quantum information.It consists of actively monitoring the presence of thermal anyons and correcting the unwanted operations resulting from their presence in the system.The Ising anyon model is an obvious candidate to benchmark such a proposal because it is numerically tractable via a mapping to free Majorana fermions [22].
The motivation behind our work stems from two distinct potential applications.On the one hand, as outlined above, anyons can arise as localized gapped excitations in two-dimensional systems, e.g., [17,18].In this setting, our method is needed to eliminate thermal defects in the systems without spoiling the topologically encoded information.On the other hand, our results can also be understood from a purely coding-theoretic perspective.Some topologically ordered systems can arise in lattice models with local commuting Hamiltonians [2,20,29].The local terms of the Hamiltonian can be thought of as check operators defining a local commuting projector code, a class of codes going beyond the usual stabilizer formalism that has been the object of recent study [24][25][26][27].Such codes allow for the implementation of quantum error-correction schemes in a geometrically local setting [2,30].The local check operators can be measured using standard circuitry [31] on an ordinary, circuit-model quantum computer.Thus, this second scenario we consider is not contingent on the existence of physical systems that naturally present topological quantum order; it could in principle be used in any quantum computing architecture with nearest-neighbor interactions on a two-dimensional lattice.In this setting, our results provide an efficient decoding algorithm for the corresponding quantum error-correcting code.To our knowledge, this provides the first efficient decoding schemes for a family of non-additive quantum codes, i.e. codes without explicit Pauli-matrix tensor product structure.

A. Summary of main results
We begin by proposing a phenomenological model of dynamics in an Ising anyon system and several noise models.Importantly, these models encompass all of the major physical error processes associated to non-Abelian anyons, including pair-creation, hopping, braiding, and fusion of anyons.We define two distinct families of noise models, following the two physical applications of our method described in the previous paragraph.To simulate the noise affecting a topologically-ordered system at finite temperature, we use a Metropolis procedure to choose among the possible physical error processes.Because the Gibbs thermal equilibrium state is the fixed point of this local noise process, we believe that it captures the essential properties of a real thermalization process.To simulate the noise affecting a quantum computer that uses an Ising error-correcting code, we use a more generic "white noise" model where each type of physical error process occurs at a predetermined fixed rate.This model is in essence an infinite temperature limit of the previous noise model, and so it probably also accurately captures the short-time dynamics of a topologically-ordered system at finite temperature.Note that we do not necessarily expect the physical noise occurring in a circuitbased quantum computer to be described by these physical processes.Nonetheless, in absence of detailed microscopic descriptions of the device, they provide a good starting point to benchmark the code (in the same sense that Pauli errors are used to benchmark stabilizer codes).
We explore the effect of these noise channels on two different encoding schemes (degenerate ground spaces).The first encoding we examine stores the information in the degenerate ground space of a system with Ising anyonic excitations embedded in a torus while the second stores information in a degenerate space associated with the fusion space of several well-separated Ising anyons.Both of these encodings lend themselves equally well to the two physical scenarios described above (although embedding an actual topologically ordered system on a torus probably poses an additional challenge).These two encodings are motivated by topological quantum memories [30] and topological quantum computing [22,32] schemes respectively.
Following this, we adapt existing techniques from topological quantum error correction to the non-Abelian setting.The non-Abelian nature of the model significantly changes the decoding problem and presents a main challenge of our study.In an Abelian model, decoding is realized by grouping nearby excitations that fuse to the vacuum.In contrast, for a non-Abelian model, it is not possible to predict with certainty the outcome of a fusion process.Thus, decoding needs to be an iterative process: in a first round particles are grouped in a certain way; if all groups do not fuse to the vacuum, then a second round groups the remaining particles, and so on.In particular, we adapt the clustering renormalization group decoder of Bravyi and Haah [33] and the perfect matching algorithm decoders [30,34] to our non-Abelian topological codes to determine their error-correction thresholds and memory lifetimes under various noise channels.Our decoders provably run in polynomial time in the size of the lattice.
Using these decoders, we explicitly calculate numerical error-correction threshold estimates for these codes and for various noise regimes, including those where hopping and braiding of anyons dominates pair-creation.Our results seem to indicate that the perfect matching decoder outperforms the clustering decoder, and that the thresholds are broadly comparable for both types of codes, and for various types of noise.We find numerical threshold estimates that vary between 13% and 25%, where the percentages are in terms of a total noise rate density which we define and discuss in detail in Sections II and VI.These thresholds are obtained by large-scale two-dimensional simulations of the dynamics of a dense "gas" of non-abelian anyons on lattice sizes of up to 48 × 48 nodes.All previous results have focused on quasi-one-dimensional models or single-particle quantum walks [12][13][14][15].
The paper is organized as follows.In Sec.II we introduce the Ising anyons, as well as our phenomenological models for their dynamics.Next, Sec.III goes into some detail about the efficient simulation of these dynamics on a classical computer.Following this in Sec.IV we give explicit constructions for the topological quantum codes that we study.We then detail our decoding algorithms in Sec.V and present numerical threshold estimates in Sec.VI.We discuss the significance of these numerical results in Sec.VII, and offer two appendices, the first of which summarizes the definition of the Ising anyons, and the second discusses the regime of validity of one of the noise models we analyse.

II. PHENOMENOLOGICAL MODEL FOR ISING ANYON DYNAMICS
An anyon model is specified by its (gapped) excitation types as well as several different ways these excitations can interact, specifically by braiding, fusion, and splitting.These three processes in fact completely characterize the universal topological properties of a system with anyonic excitations.Other processes will depend on microscopic details of the system and will be largely irrelevant for the behavior of the global topological degrees of freedom of the system that we will be interested in.For this reason they can safely be ignored in this work.
We will now describe a phenomenological model of these topological dynamical processes that allows us to simulate the dynamics of Ising anyons.We will use this to discuss several important classes of noise models and quantum codes and study the effects of these noise models on the codes we define.
The Ising anyon model consists of two distinct non-trivial excitation types (or charges) conventionally labelled ψ and σ, and for convenience we label the vacuum (no particle) by I. Since this is a non-Abelian anyon model, there are multiple ways in which two anyonic charges can combine by fusion.These are specified by the fusion rules [35] The full data specifying Ising anyon dynamics are given in Appendix A and describe the braiding, fusion, and splitting of Ising anyons.
We consider dynamics on a discretized space defined by a graph G = G(V, E).Associated with each site i ∈ V of the graph are a set of occupation variables n i q (taking nonnegative integer values) which specify the number of particles with charge q ∈ {ψ, σ} located at i.The collection of all occupation variables specifies the charge configuration on the lattice.Although we define charge configurations with multiple particles at each site of the graph, we will mostly regard all particles at a given site as physically being fused together, though we denote them separately for convenience.
In order to complete the specification of the total state of the system, we also require a Hilbert space, which we call the fusion space, associated with fusion outcomes of any σ particles on the lattice (ψ particles have unique fusion outcomes, and so will not contribute to this space).For a charge configuration with 2m σ particles (recall that σ particles can only be created or destroyed in pairs), the fusion space will be 2 m dimensional.
In order to model the most general allowed processes on a topological system, processes of braiding, fusion, and splitting are most easily recast as the following elementary operations on the graph G.
1. Creation of particle-antiparticle pairs of the form q × q for some q ∈ ψ, σ on sites incident to edge (i, j) ∈ E; 2. Hopping of all charges from one site i to a neighboring site j; 3. Exchange of all the charges from one site i with those on a neighboring site j either clockwise or anticlockwise; 4. Decoherence of the total charge of a site i.
Here we use the notation q for the antiparticle of q, i.e., the particle type such that q × q = I + . ... For Ising anyons, each particle type is its own antiparticle.The decoherence process is the only one Charges at i, j Process I, I I, q q, I q, q Pair Creation Hopping Exchange TABLE I. Allowed elementary noise processes for our phenomenological noise model of anyon dynamics.The check marks denote allowed processes on a directed edge i → j of the underlying graph G, and q,q denote non-vacuum anyon charges.
of these four which may require some explanation: this process corresponds to projecting into a particular fusion outcome of the set of charges at a site.In general, our model allows superpositions of different fusion outcomes at a site, and depending on the physical system we are attempting to model these may tend to decohere slowly or rapidly.
Our aim is to simulate the dynamics of a system with Ising anyon excitations in contact with a thermal bath for a short period of time.In order to do this, we have two complimentary methods for sampling from the above processes.

A. Fixed Rate Sampling
The first sampling mechanism is the most naive and the most general.It proceeds by taking a set of rates for each fundamental anyonic noise process: γ q c for pair creation of charge q, γ h for hopping and γ e for exchange.Since decoherence is a slightly different type of noise process, we will treat it separately for convenience, and associate a decoherence probability p d to decoherence events.Since we will be interested in calculating memory lifetimes for a fixed strength noise channel, any rescaling of the sum of the γ noise rates will simply correspond to a rescaling of the memory lifetime.For this reason, the rates γ should simply be regarded as relative rates of the different noise processes, and we will always normalize the total rate to 1.
We simulate the dynamics of the system for T timesteps with T drawn from a Poisson distribution of mean T 0 .At each step a single operation from the set {pair creation, hopping, exchange} will be performed.Specifically, this proceeds by randomly choosing a directed edge e = (i, j) uniformly from the graph, and then selecting from the nontrivial processes allowed on that edge by their relative rates.We call a process trivial if it leaves the state invariant.The states on which each process may act nontrivially are summarized in the Table I, where we use q, q to denote any non-vacuum anyonic charge.
Given the set of nontrivial processes that are allowed on the selected edge, an operation to perform is chosen according to the relative rates γ of these nontrivial processes.Following each such operation, the decoherence process is applied to every site with probability p d .After T timesteps, this will approximately correspond to a simulation of the system running for time t sim ∝ T 0 |E| for a graph with |E| edges.
This fixed rate sampling mechanism is a natural noise model in a quantum computing architecture which directly implements our error correcting codes, and may also be a good approximation to short-timescale thermalization processes in Ising-anyon quantum systems for appropriately chosen rates.However, the timescale on which it is a good approximation in this latter case may be very small.For this reason, we also consider an alternative sampling mechanism.

B. Metropolis Sampling
In order to more faithfully capture the thermalization process in our model, we consider a second sampling mechanism based on the Metropolis sampling method.In order to specify a Metropolis procedure we give a mass to each particle type m ψ and m σ , and then define the Hamiltonian for P q i the projector to total charge q on site i.The Metropolis method proceeds by proposing an elementary noise operation uniformly at random, and then calculating the energy difference between the state before and after this operation ∆ E .The operation is then accepted with probability min{1, e −β∆ E } for a system at inverse temperature β = 1/k B T .This procedure guarantees convergence to a thermal state obeying detailed balance.
Unfortunately, the hopping and pair-creation processes we wish to consider will not always take energy eigenstates to energy eigenstates, and so it will be convenient to consider a restricted set of noise operations that will ensure that our system remains in an energy eigenstate at each timestep of its evolution.To achieve this, after each operation is applied we decohere every site of the lattice, projecting into an energy eigenstate of the system.This will allow us to consistently define Metropolis acceptance probabilities at the cost of taking the extreme decoherence limit of our noise model.Because it is a local noise process that converges to the thermal state with detailed balanced condition, we believe that this sampling method captures the essence of real thermalization processes.
We can consider this a semiclassical noise model, since we effectively disallow superpositions of different total charges at each site.We discuss the limitations of this noise model in Appendix B.

III. SIMULATING ISING ANYON PHENOMENOLOGY
It is known that the dynamics of Ising anyons can be efficiently simulated [22].By this we mean that the operations induced on the fusion space and charge configuration space by pair-creation, braiding, and fusion of n Ising anyons can be simulated, modulo global phases, by a classical computer in time polynomial in n.We will make use of this result to simulate the dynamics of our phenomenological model under various noise channels.This will allow us to determine an errorcorrection threshold for our system by observing when it is able to preserve quantum information.
Recall that in the Ising anyon model, the ψ particles do not contribute to the fusion space as they have unique fusion outcomes.Thus we can associate the fusion space with the set of σ particles.In order to explicitly quantify the effects of braiding operations on this fusion space, it is convenient to assign a Majorana fermion mode to each σ particle present in our system.This allows us to define a set of Majorana operators ĉα for each σ particle α with ĉα ĉβ = 2δ αβ − ĉβ ĉα (3) For 2m σ particles, this defines a 2 m dimensional space which will serve as the fusion space for these particles.A simulation of Ising anyon dynamics requires the simulation of both the charge configurations and their associated fusion spaces.

A. Planar graph simulations
For simplicity, let us first describe a simulation where G is a planar graph.With this constraint, all physical operations will preserve global charge.That is, any physical observable O satisfies where Q is the total charge operator The eigenvalues of Q correspond to the possible total charges of all σ particles, and after fixing this global charge to be vacuum, the remaining fusion space is 2 m−1 dimensional.
Braiding processes between anyons act as unitary transformations on the fusion space.In order to specify the transformation corresponding to a braid operation between anyons on adjacent sites of G, we must first specify a basis for this space by defining an ordering of the σ particles in the system.It is convenient to choose the ordering of the σ particles relative to a natural ordering of the sites in the graph G.This is given by a bijective function over sites of G (i.e.simply a fixed labeling of sites).We can then choose an ordering function # σ for the σ particles in our system in a way that is consistent with the # G order.By consistent, we mean that for σ particles α and β at sites i and j respectively with # G (i) < # G (j), we require that # σ (α) < # σ (β).For clarity, we will conventionally use Roman letters to refer to sites of G, while using Greek letters to denote σ particles in the # σ ordering.
Braiding operations involving ψ particles act only on the fusion space by accumulating global phases (see Appendix A for details).For this reason, we only need explicitly consider the action on the fusion space of braiding operations between σ particles.The most convenient generating set of operators for braiding and fusion of σ particles is that given by operations between # σ -neighboring particles.Note that we will conventionally take the subscript of the Majorana operators ĉα to refer to this ordering of the corresponding σ particle.When performing a clockwise braid between two # σ -adjacent particles σ α and σ α+1 , we enact the following operation on the fusion space (up to global phases) [22] The anticlockwise exchange of these particles corresponds to the operation These operators can be seen to map the Majorana operators as follows: We track the state of our system in the fusion space in the Heisenberg representation; that is, we will represent the state by a group of (commuting) stabilizer operators whose common +1 eigenspace defines the desired state.For each pair of σ particles in our system, we will include an additional generator for the stabilizer group S. The stabilizer group will always contain the total charge operator Q.
As well as the operation on the fusion space corresponding to braiding particles, we must also specify an analogous operator corresponding to fusion of neighboring particles.The fusion product (corresponding to the combined charge) of a pair of # σ -neighboring particles is represented by the operator The +1 eigenvalue corresponds to combined vacuum charge, while the −1 eigenvalue corresponds to combined ψ charge.When two σ particles are created from vacuum, two new modes are created and initialized in the +1 eigenstate of a corresponding M (2) operator.
Similarly, the fusion product of 2l neighboring σ anyons is given by Decoherence or measurement of the total charge of a set of σ particles corresponds to a measurement of the relevant M operator.We are able to neglect global phases acquired during braiding operations as long as we guarantee that our system remains in a common charge eigenstate at each stage of the simulation.That is, as long as the occupation variables n ψ and n σ can be treated classically.Our phenomenological model involves only processes which take charge eigenstates to charge eigenstates, and so the (unitary) braiding and (projective) measurement operations above suffice to generate all required dynamics of our system.
Although the braid and measurement operations introduced above are enough in principle to specify all the dynamics of our system, there is no sense in which they correspond to the elementary processes of our phenomenological model of dynamics.The phenomenological processes are operations on the graph G, while the braid moves described by ( 9)-( 12) corresponds to operations on the line specified by the # σ ordering.In order to translate between these two pictures, we introduce the notion of linearized braid sequences.

B. Linearized braid sequences
In order to explicitly compute the effect on the fusion space of an operation involving Gneighboring sites (e.g.exchange), we must map it to a sequence of # σ -neigboring operations.We will call this sequence the "linearized" braid sequence corresponding to a graph edge e = (i, j).This will give us a prescription to take some subset of the charges at i to site j, where another operation may be performed such as exchange or hopping.
The linearized braid sequence is most easily understood graphically, and is shown in Fig. 1.To define it, we embed G in the plane, and overlay a planar path connecting each node of G in the order # G .We can then smoothly deform this graph until the # G path is a straight line, and the image of an edge e = (i, j) ∈ E under this deformation provides us with the prescription to translate a G-neighboring operation on e into a sequence of # G -and # σ -neighboring operations.
FIG. 1.An example choice of G (shown with gray edges) and # G (shown with solid black edges).We also highlight several edges of G to demonstrate the map between G-neighbors and paths over # G .After embedding G in the plane, we overlay a linear graph corresponding to the # G ordering.Following this, we can map any desired lattice edge to a braid sequence on the line smoothly.
To move the charges at site i to site j along edge e, we follow the deformed edge as it moves along the # G path from i.If the edge moves over the top of a node k in # G , then moving right along that edge corresponds to a sequence of clockwise exchanges between the σ particles from i with those at k. Moving left corresponds to a sequence of anticlockwise exchanges.If the edge moves under a node in # G , the braid sense is reversed.The fact that the # G and # σ orderings are consistent guarantees that these operations are sequences of # σ -neighboring exchanges.These exchanges are implemented mathematically by conjugating the stabilizer group by the relevant braid operator B α or B † α .We successively exchange the charges from i along the path until they arrive at j, where the desired physical process can be performed.
Using these linearized braid sequences, we can simulate the phenomenological processes defined in Sec.II as follows: 1. To pair-create q-type charges along edge e = (i, j), we create two charges of type q at site i, initialized in the vacuum fusion channel.Following this, we perform a linearized braid sequence to take one of these charges from i to j, where it is added to the total charge at j; 2. To hop charges along edge e from i to j, we perform a linearized braid sequence taking charges from i to j, where they are then simply added to the existing charges at j; 3. To exchange charges along edge e = (i, j) either clockwise or anticlockwise, we perform the linearized braid sequence taking charges at site i to j, where they can be exchanged with the j charges in the relevant sense, and then return charges from site j to site i by the reverse linearized braid sequence; 4. To decohere the total charge at site i, we have no need to linearize the process as it acts only on a single site.We simply perform a projective measurement of the total charge at i.If n i σ is odd, the total charge will necessarily be σ, otherwise the total charge can be calculated using Eq. ( 14) (also taking into account any existing ψ particles at i).
In this way we can simulate the phenomenology of an Ising anyon system in contact with a thermal bath.
In order to perform an error-correction routine on this system, we will require one further algorithm which simulates fusion of the total charge within a region of the graph.We achieve this by using linearized braid sequences to hop all charges in the region to a common site, after which we successively projectively fuse all pairs of σ particles and calculate the total charge on this site.

C. Non-planar graph simulations
As well as planar graphs, we will also consider anyon dynamics on non-planar graphs.Specifically, we will treat a square lattice graph with periodic boundary conditions.We will restrict discussion to this case, though discretizations of general 2-manifolds follow in a similar manner.
The behavior of anyons braiding on a torus will differ from braiding on a plane insofar as the torus admits braiding processes around non-contractible loops on the manifold.Mathematically, the braid group on a nontrivial manifold requires additional generators to describe these processes compared to the braid group on the plane, as is discussed in Ref. [36].
In contrast to a sphere, anyon models on a torus define a degenerate Hilbert space even in the absence of anyonic charges.A basis for this space can be defined by fluxes corresponding to each anyonic charge running through a nontrivial loop of the torus1 .An anyon braiding around this nontrivial loop will behave as if it had performed a monodromy of the corresponding enclosed charge.
Although there are two nontrivial loops on a torus, they do not each carry an independent associated flux.The fluxes through these two loops are related by topological S transformations, given for the Ising anyon model in Appendix A. The behavior of anyons braiding around the nontrivial loops of the torus can also be analyzed by the extended diagrammatic calculus of Pfeifer et al. [37].
We will now describe the simulation of Ising anyon phenomenology on a torus, concentrating on the ways in which it differs from the planar case.Consider the graph in Fig. 2.Although not planar, this graph has a natural embedding in the torus.As in the planar graph case, it is convenient to choose a linear ordering of sites of G as shown.Again, this allows us to find linearized braid sequences for a given graph edge from its image under a smooth deformation of the graph to the line.An Ising anyon system on a torus has an associated 3-dimensional "topological" Hilbert space in addition to the fusion space and the charge configuration space.This represents the flux running through a nontrivial loop of the manifold.There are two natural orthonormal bases for this space corresponding to the two nontrivial loops.Calling these loops h and v as shown in Fig. 3, we denote these bases as {|I h , |ψ h , |σ h } and {|I v , |ψ v , |σ v }.Given the charge configuration corresponding to vacuum at all sites of G except charge q at a single site (all charge configurations can be transformed into a state of this form by braiding and fusion operations that do not cross the seams of the torus; see Fig. 2), we can transform from the h basis to the v basis by topological S transformation S q (see Appendix A for details).For Ising anyons, these have the form Note that on a nontrivial manifold, the total charge of the system is no longer a conserved quantity, though for the Ising anyons we can guarantee that the total charge will never be σ and so we need not consider S σ .
Each time an anyon moves around a nontrivial loop of the torus, it effectively braids around the flux through that loop [37].For example, if an anyon braids along the green path in Fig. 3 while the topological state is |q h , this anyon acts precisely as if it had performed a double exchange with charge q.In contrast, the blue path corresponds to performing a monodromy of the flux in the v basis.In order to determine the effect of the relevant operations, we must use an S transformation to change between the h and v bases in between these processes.The S-matrix required will depend on the total charge on the lattice at any given point in time.
We can derive the rules for anyons braiding around the nontrivial loops of the torus by using the Pfiefer et al. diagrammatic calculus [37] or by simply considering the S transformations and the braiding rules of the Ising anyons.Consider a process labeled L h q , by which we create a pair of q × q particles from vacuum, braid one around an h loop, then fuse them together again.We find that the action of this operator on the various basis states is as follows: where |I bulk has I charge at every site of G, and |ψ bulk denotes a state with a single ψ charge at one site, with all other sites taking charge I (note that all of these states have trivial fusion spaces).
We also consider the analogous operators corresponding to braiding around the v loop, L v q : Note that these relations make explicit the connection between the h and v loops of the torus, insofar as the identity (L v q ) † = SL h q S † holds, where S acts as S I or S ψ depending on the bulk charge.We have made a number of convention choices here, most significantly choosing a preferred braiding sense around each of the nontrivial loops.When an anyon braids in the preferred sense, it corresponds to performing the relevant L operation, while braiding in the reverse direction corresponds to the adjoint operation L † .
The definitions of L v and L h operators allow us to not only study the specific processes by which they were defined, but also general braiding processes on the torus.A general braiding operation on the graph in Fig. 2 will perform some action on the fusion space by braiding with other anyons, as well as applying a "topological correction" operation on the ground space according to the relevant L operator.This topological correction must be performed each time an anyon moves along a "seam" edge in Fig. 2 (as well as the standard linearized braid operations that are required).This prescription will ensure that operations on the torus satisfy appropriate Fox braid group relations for this nontrivial manifold [36,38].
Although we could directly simulate the processes corresponding to ( 16)- (23), in practice there is a more straightforward method to simulate these processes up to global phases that takes advantage of the structure of the Ising anyons [8,39].This corresponds to tracking the evolution of the (commuting) operators L h ψ and L v ψ .Since they are self-inverse, these operators each have eigenvalues ±1.The corresponding four eigenstates span the space and similarly under exchange of h and v), this is sufficient to track the ground state of the system under topological correction operations.Finally, when a σ particle crosses a seam for which the corresponding L ψ operator has eigenvalue −1, the change in bulk charge that occurs causes the σ particle to acquire an additional ψ charge (mathematically, its Majorana operator's eigenvalue changes sign).

IV. QUANTUM ERROR CORRECTING CODES IN ISING ANYON SYSTEMS
In order to discuss error correction in an Ising anyon system, we must define an encoding of a logical Hilbert space in our model.In anyonic models, the construction of this encoding generally follows one of two paths.The first is motivated by topological quantum computation schemes, where quantum information is encoded in the fusion space associated with a fixed number of non-Abelian anyons.This kind of scheme is well developed for Ising anyons [22].The second is motivated by topological codes such as the toric code, where the degenerate ground space on a nontrivial manifold is used as a codespace [30].We present constructions for both types of code and compare their properties.
During this discussion, we will abstract away much of the explicit formalism introduced earlier for simulating Ising anyon dynamics.As before, we regard our system as being defined over a graph G with charges located at each site of G.Here we represent the charge configuration space by 3-dimensional Hilbert spaces for each site s ∈ V with basis {|I s , |ψ s , |σ s } representing the charge present at s.We will represent the state within the fusion space (if one exists) by kets of the form |q f αβ for the state where the σ particles labelled α and β fuse to outcome q ∈ {I, ψ}.

A. Ising Fusion Code
The first code we consider is constructed in the fusion space of several anyons placed at preferred sites of a graph G.We call the resulting code the Ising fusion code (IFC).
Consider a system of L 2 sites arranged in a square lattice as shown in Fig. 4. The boundaries are chosen to give G a spherical topology.There are four preferred sites on the lattice chosen to be as far apart as possible in graph distance (which will maximize the code distance), and we call these preferred sites "code" sites.We refer to all non-code sites as "bulk" sites.
Since G has spherical boundary conditions, it follows that all physical actions on our system will conserve global charge.We choose to work within the total vacuum charge sector of our model.The IFC is defined by vacuum charge at each bulk site, and a σ charge at each of the code sites.Although the code space corresponds to a fixed charge configuration, a 2-dimensional fusion space is also associated with the four σ particles in the system.We can define a basis for this space by considering the fusion product of a preferred pair of σ particles, which we choose to be those at code sites c t and c r .These basis states have the form: where we denote the σ particles at each code site by t, r, l, and b.
Logical operators in this code can be taken to correspond physically to braiding the code site σ anyons around one another.Equivalent operations can also be achieved e.g. by tunneling processes between code sites.These operations have been well studied, and all Clifford group operations on one qubit can be achieved in this way [40,41].As an example, the logical Pauli X operator (note that the codespace is a qubit) can be implemented by performing a monodromy between the r and b σ particles, while a logical Pauli Z operator corresponds to a monodromy between the t and r anyons.
We remark that we could consider a version of the code we have described here to be implemented as the ground space of a noninteracting Hamiltonian with the form Such a Hamiltonian could emerge as a low-energy effective Hamiltonian in a spin-based model, for example [19][20][21].

B. Ising Topological Code
We also consider a code defined as the ground space of an Ising anyon system embedded in a nontrivial manifold.Specifically, we choose G to be an L × L square lattice with periodic boundary conditions as in Fig. 2, and imagine a local gapped Hamiltonian whose ground space corresponds to vacuum charge at every site of G.An arbitrary anyon system of this form embedded in a torus has a ground space dimension equal to the number of distinct anyon species (including vacuum), and so the ground space of our Ising anyon system will be three dimensional.We call the resulting code the Ising topological code (ITC).
Basis states for this ground space can be labelled by fluxes taking values {I, ψ, σ} running through one of the nontrivial loops (labelled h and v) of the torus.We denote these basis states by for µ ∈ {h, v}.As noted in Sec.II, the basis transformation corresponding to exchanging the two nontrivial loops on the torus (exchanging the h and v bases) is given by a topological S transformation.
The logical operators of this code correspond to braiding of anyons around these nontrivial loops of the torus.Explicitly, they correspond to the operators L q defined in Sec.III C.Although L h σ and L v σ are not strictly logical operators (as they do not preserve the codespace), they correspond to an irrecoverable decoding failure for this code, as with any of the rest of these operations.The only distinction is that errors corresponding to L h σ and L v σ may be detectable.One might think that if this error is detected by a measurement of the total bulk charge we could correct it by performing the inverse braid, but by detecting it we have performed a partial measurement on the state of the system and so corrupted the encoded information.Note however that (L h σ ) 2 and (L v σ ) 2 are valid logical operators.
We could consider a version of the code we have described here to be implemented as the ground space of the Hamiltonian or any alternative local Hamiltonian that energetically penalizes any non-vacuum charge in the bulk.

C. Codes with boundaries
When discussing topological codes, often planar codes are considered as they may offer a more feasible implementation than a code which requires a nontrivial manifold (such as the ITC).The most well known of these is the surface code, the analogue of Kitaev's toric code with boundary [30,42].However, these codes require the ability to construct a local gapped boundary at which some subset of the charges in the anyon model may condense.The construction of these boundaries is not always possible for an arbitrary anyon model, and the exact boundaries that can be constructed may even depend on whether the microscopic Hamiltonian is based on bosonic or fermionic degrees of freedom.
Abelian anyon boundaries are discussed in [43,44], and although the non-Abelian case (including Ising anyons) is not treated it seems likely that an analysis of this type could be extended to show that no gapped boundaries are possible for an Ising anyon system based on an underlying bosonic Hamiltonian.However, it may be possible to construct a gapped boundary which absorbs ψ charges in an Ising anyon system based on a fermionic Hamiltonian.This would allow for the implementation of an Ising anyon code defined on an annulus (topologically equivalent to the cylinder) with a 2-dimensional codespace, or similarly a family of codes with many internal boundaries.We do not treat this case explicitly here, though the construction and simulation of these codes would follow straightforwardly from standard topological code techniques and other techniques developed here.

D. Determining code thresholds
A threshold for a code is a critical value of some parameter of the noise channel such that in the thermodynamic limit (i.e. as L → ∞), we can guarantee that the recovery operation calculated by the decoder will exactly reverse the effects of the noise channel (that is it will restore the encoded quantum state with unit probability) if the parameter is below the critical value.Of course, this threshold will depend on the decoder and noise channel in general.
To the best of our knowledge, all previous estimates of thresholds (numerical or otherwise) for quantum codes have only been made for additive codes (i.e.codes described by the Pauli stabilizer formalism [45] or the higher dimensional generalizations thereof [46][47][48]).In these codes, it is straightforward to calculate the precise effect of the noise and recovery operations on the codespace in generality.In our codes, this is not so trivial, and so we introduce a new tool to the study of quantum codes to calculate thresholds for the IFC and the ITC.
In order to estimate the threshold for the codes and decoders we have defined, we first benchmark our error-correction protocol by calculating quasi-thresholds2 below which our recovery procedure will perfectly preserve a set of particular code states (as opposed to preservation of the whole codespace) in the thermodynamic limit.This can easily be done by simulating the initialization of the codespace in this state, and performing a measurement in an appropriate basis after the recovery operation has been performed.
Next we will see how we can guarantee preservation of the entire code space by simply looking at a few quasi-thresholds for specific states, thus avoiding the problem of tracking the full logical subspace.The structure of the space that is perfectly preserved by the noise and recovery operations in the infinite limit must form a closed matrix algebra [49,50].For a d-dimensional codespace, if we choose d non-orthogonal pure states that span C d and guarantee that they are each preserved, it follows that the entire Hilbert space must be preserved.This is a consequence of the fact that the density matrices corresponding to these d states generate all density matrices over this space as a matrix algebra.Thus, the threshold for the code will be the minimum of the d quasi-thresholds arising from these states (below which, we can guarantee that each of these states are preserved in the thermodynamic limit).This cannot be achieved with fewer than d pure states.
To prove these claims, define d non-orthogonal states ρ i = |φ i φ i |, and note that ρ i ρ j ∝ |φ i φ j | for any φ i |φ j = 0. Since the {|φ i } span the space, we can construct any density matrix as a linear combination of |φ i φ j | (since we can take linear combinations of |φ i to form any pure state, we can take linear combinations of |φ i φ j | to form any density matrix).The fact that this cannot be achieved with k < d pure states can be seen by noting that the set {ρ i ρ j } is closed under multiplication (up to proportionality) and contains k 2 unique elements for 1 ≤ i, j ≤ k.The Hermitian part of the span of these matrices is also spanned by {(ρ i ρ j + ρ j ρ i )} (by linear independence of the ρ i ρ j ), of which there are k(k+1) 2 unique elements.Since an arbitrary density matrix is specified by d(d+1) 2 − 1 parameters, the smallest k to generate all density matrices as linear combinations of {ρ i ρ j } is d.
Note that mixed states can in general reduce this upper bound below d to as little as 3, but they introduce other challenges in verifying that they are preserved by the joint noise-recovery operation, so we eschew them here.
In terms of our simulations this means that since the IFC codespace is 2-dimensional, if we calculate a quasi-thresholds for a computational basis state |0 and a conjugate basis state |+ then we can guarantee that the threshold for the code will be the minimum of these two quasi-thresholds.Similarly for the ITC, the codespace is 3-dimensional, and so we can obtain the code threshold by taking the minimum of three quasi-thresholds corresponding to independent non-orthogonal states.

V. DECODING ALGORITHMS
In order to calculate a memory lifetime for the codes described in Sec.IV under a noise channel described by the phenomenology of Sec.II, we initialize in a preferred codestate |φ 0 C , before sampling noise processes by either a fixed rate or Metropolis mechanism for T timesteps.For a lattice of linear size L (where L = |E| 2 for our models), we choose T according to a Poisson distribution with mean T 0 = t sim |E| for some t sim > 0. This t sim represents the physical time over which we have coupled the system to a thermal bath.If, after applying a decoding algorithm, we arrive back in the |φ 0 C state, then the decoding succeeded; otherwise, we say there was a logical error.Loosely speaking, the threshold is the critical time t * such that for all t sim < t * the logical error rate decreases asymptotically to zero as a function of the lattice size L for all |φ 0 C in the code space.As noted in Sec.IV D, for a d-dimensional codespace, the preservation of d non-orthogonal linearly independent states is sufficient to guarantee the preservation of the entire space.
The critical lifetime t * of the code will depend not only on the rates γ, but also on the decoder used to readout the quantum information.
We adapt two algorithms used in decoding Abelian anyon systems, a clustering renormalization group (RG) decoder due to Bravyi and Haah [33], and the perfect matching decoder [30].These decoders rely on the metric structure of the surface in order to make decisions about which (nearby) anyons to fuse together as part of the recovery process.Compared to the Abelian decoders, several additional complications arise from the presence of non-Abelian anyons in our system.The most obvious is the fact that in the Abelian setting, it is always possible to predict fusion outcomes of any pair of anyons with certainty.In the non-Abelian setting, this is no longer possible in general, and so decoding must proceed iteratively as the decoder proposes fusion operations and queries the system to learn the outcome of these fusions.Another complication involves the fact that the precise path taken by two anyons being fused (not just the homology class of the path) plays a role in the outcomes of the process.For Abelian anyons, braiding between two anyons simply results in a global phase that we can disregard.When dealing with non-Abelian anyons, two paths that differ by a braid around another anyon may give rise to differing fusion outcomes, and so we must be more precise when the decoder specifies the recovery operation to be performed.

A. The Simple Clustering Decoder
The first decoder we consider uses an agglomerative clustering algorithm inspired by the clustering RG decoder of Bravyi and Haah [33].Related decoders have also recently been proposed [51,52].The Bravyi-Haah RG decoder proceeds by placing every charge in its own cluster and then iteratively: growing each cluster, merging the overlapping clusters, and then fusing all charges within each cluster.
When applied to Abelian anyon systems, the cluster decoder only needs to measure the system once to learn the initial charge configuration, as this implies the charge configurations at all points during the decoding and recovery routine.For the same reason, the fusion of charges within each cluster need not be applied in sequence, as they can be simulated on a classical computer and applied to the system altogether in one final step.
In contrast, when applied to a system of non-Abelian Ising anyons, the clustering decoder must measure the charge configuration and physically perform the fusions before each round of clustering.In our simulation this is done by the decoherence/measurement and fusion protocols defined in Sec.II.Thus, our simple adaptation of the clustering decoder can be summarized as follows: function Cluster Simple Measure charges on each site Initialize a cluster for each site with non-trivial charge while more than one cluster remains do Fuse all charges within each cluster Eliminate any cluster with vacuum charge Grow each cluster by a constant length Merge overlapping clusters end while end function Clusters are grown according to the graph distance metric.In order to merge overlapping clusters, we construct a spanning tree for each cluster and then fuse along the branches of the tree.The charge resulting from the fusion of all anyons is placed at the root of the tree.
When applied to the IFC we must also slightly modify this algorithm to ensure that σ particles are returned to the code sites if they have been moved away by thermal processes.This can be achieved simply by initializing clusters on the code sites and returning any fusion results from clusters containing the code sites to these sites (also in this case the termination condition for the algorithm must be modified).
The original clustering RG decoder [33] has a provable threshold and we expect that their analysis can be extended to the non-Abelian case, as it relies primarily only on the separation of defect clusters (though we note that the proof of threshold for the clustering RG decoder in Ref. [33] crucially depends on an exponential schedule of cluster size increases, as compared to the linear schedule of our algorithm).We also note that our adapted clustering decoder inherits its polynomial runtime guarantee from the arguments in Ref. [33].

B. The Fusion-Aware Clustering Decoder
The simple clustering algorithm defined above does not take advantage of all the structure of the anyonic excitations present in our system.In particular, we can attempt to optimize this decoder by making use of details of the anyon fusion algebra.Instead of clustering together all nearby charges on the lattice, we can choose to cluster together only particular types of charge at each iteration.Although we will only describe this process explicitly for our Ising anyon system, the ideas extend naturally to arbitrary anyon models.
The algorithm is based on the following subroutine.
function Cluster Typed(set of charge types Q i ) Measure charges on each site Initialize cluster for each site with non-trivial charge in Q i while more than one cluster remains do Fuse all charges in Q i within each cluster Eliminate any cluster without charges in Q i Grow each cluster by a constant length

Merge overlapping clusters end while end function
The decoder then proceeds by performing Cluster Typed on several sets of charges Q i for i = 1 to χ. Explicitly, we define the following function.
function Cluster Aware for i = 1 to χ do Cluster Typed(Q i ) end for end function We can see that the Cluster Simple decoding routine simply corresponds to choosing χ = 1 and Q 1 = {I, ψ, σ}.We customize this process for the Ising anyons by choosing χ = 2 and Q 1 = {I, σ}, Q 2 = {I, ψ}.This amounts to fusing all pairs of σ particles into either vacuum or ψ, and then fusing all remaining ψ particles separately.It generally produces slightly higher thresholds than the simple clustering decoder in the regimes of interest.

C. The Perfect Matching Decoder
The second kind of decoder we use is based on the perfect matching algorithm (PMA) of Edmonds [53] and its descendants [54].The PMA decoder is typically used for Abelian anyon systems whose charges are self-inverse (i.e.q = q for all charges q) such as the toric code [2].It proceeds by calculating the minimum-weight matching between anyons that can fuse to vacuum.
The naive generalization of the PMA decoder to non-self-inverse anyon models would calculate minimum weight matchings between hyperedges of three or more particles; deciding the existence of such hypergraph perfect matchings is NP-complete in general [55] as is the minimization problem.It is therefore very unlikely to be efficient (though heuristic methods have been used to approximate the solution to this problem [56]).
Fortunately, the Ising anyon model has only self-inverse charges, and so we need not consider the hypergraph matching problem here, though similar methods could be applied to arbitrary anyon models if a heuristic hypergraph matching algorithm were used in place of the PMA.
Our adaptation is algorithmically similar to the standard implementation of PMA-based decoders for Abelian anyon systems, but there are qualitative distinctions in its implementation due to the indeterminacy of the fusion outcomes of non-Abelian anyons.The PMA decoder is based on the following subroutine: function PMA Single(charge q, configuration of q-type anyons n s q ) Construct a complete graph K q on the set of q-type anyons, with edge weights given by the minimum length path between anyons Use the PMA to compute a minimum-weight perfect matching of this graph Fuse pairs of anyons that are matched in K q end function In our implementation of this algorithm, the path along which pairs are fused is taken to be any shortest length path between the two anyons.
In the case of the IFC, this algorithm must again be modified slightly to take into account the fact that the codespace is not the vacuum state.This can be achieved by considering the code sites as boundaries and calculating a matching over a modified graph as in Ref. [34].Note that when applied to the ITC, the fact that the global charge is not conserved may lead to a graph without a perfect matching.This corresponds to decoding failure.
When applied to Ising anyons, we make use of the structure of the fusion algebra as in the fusion-aware cluster decoding algorithm Cluster Aware.A more general algorithm could be devised for arbitrary self-inverse anyon models (or arbitrary anyon models if a heuristic hypergraph matching algorithm replaces PMA), but we restrict discussion to the specific Ising anyon decoder.This has the form function PMA Ising Measure charges on each site s PMA Single(σ, n s σ ) Measure charges on each site s PMA Single(ψ, n s ψ ) end function The PMA decoder has complexity O(N 3 ) for a N site system [57] (again assuming constant-time fusion and measurement protocols), though related decoders [58] achieve average complexity O(N ) and are parallelizable to average time O(1).The thresholds given by the PMA decoder are typically higher than those obtained from the clustering decoders.

VI. NUMERICAL RESULTS
To estimate error threshold values for our codes, we use Monte Carlo sampling with various noise models and decoders.The four main questions we pose are: 1. What are the thresholds/memory lifetimes of Ising anyon codes? 2. How do these values vary between different Ising anyon codes? 3. How do different decoding schemes perform on these codes? 4. How do different noise models affect the thresholds/lifetimes of these codes?
While the first question is the most natural, it is the least clearly posed and so we will postpone discussion of it until Sec.VI E. We will first discuss questions 2-4 in turn.
In each of our threshold plots, the x-axis is given by t sim , or equivalently the average number of errors per edge T 0 |E| .This is the natural measure of noise strength in our system and is equivalent to the time spent in contact with the thermal bath before decoding.We discuss this measure further in Sec.VI E. The quasi-thresholds and thresholds we extract from our data will thus be given as critical values t * , where of course the specific value of t * depends on the choice of code, noise model, and decoder.

A. Comparing code performance
In order to compare the performance of the IFC and the ITC, we restrict to a particular simple noise channel.Specifically, we consider the case where γ ψ c = γ σ c = 0, all other γ = 0, and p d = 0.Although we present explicit data only for this choice of parameters, further numerics suggest that the results we find are fairly insensitive to variation of the noise channel.In particular, under no noise model does the existence of a threshold seem to change.For simplicity, we will also exclusively restrict to use of the PMA decoder in this subsection.The implementation of PMA we use is Blossom V [59].

Ising Fusion Code
As described in Sec.IV D, in order to determine a threshold for the IFC we first determine quasithresholds below which we can guarantee the preservation of the |0 and |+ states respectively in the thermodynamic limit.The minimum of these two quasi-thresholds is the threshold for this code, below which we can guarantee the preservation of the entire logical qubit in the thermodynamic limit.
Preservation probabilities for the |0 and |+ states are shown in Fig. 5.We find that the two quasi-thresholds for this system are indistinguishable up to numerical uncertainty, and thus find the threshold for the IFC as t * IFC ≈ 0.24.

Ising Topological Code
Our simulations of the system embedded on a torus Fig. 6 yield a threshold of t * ITC ≈ 0.245, very close to the threshold obtained on the sphere with four localized anyons.We note however that the decoding failure rate is systematically lower on the torus, a phenomenon that we attribute to the fact that the code's minimum distance is essentially doubled on the torus.

B. Comparing decoder performance
In order to compare the different decoders described in Sec.V, we will again restrict to the simple noise model described by γ ψ c = γ σ c = 0, all other γ = 0, and p d = 0.As in the previous section, the results we find seem insensitive to this We will the IFC threshold as a benchmark to compare the relative performance of our alternative decoding algorithms.We will present only plots for preservation of the |0 state, as plots for the |+ state are quite similar and give identical threshold estimates.
The threshold estimate obtained from the PMA decoder was determined in the previous section to be t * IFC ≈ 0.24.This compares to the numerical results for the simple clustering decoder and fusion aware clustering decoder shown in Figure 7.These can be seen to demonstrate threshold values of t * SC ≈ 0.14 for the simple clustering decoder and t * FAC ≈ 0.15 for the fusion aware clustering decoder.
Note that as well as giving a much higher threshold than the clustering decoders, the PMA decoder also produces numerical data in which the threshold is much more easily identifiable and the uncertainty in the threshold value is clearly lower.This gives a justification for our use of the PMA decoder in the comparison of the different codes and noise models above and in subsequent sections.

C. Variation with relative (fixed) rates
Here we consider the fixed rate sampling mechanism described in Sec.II.By varying the parameters of our noise model, we can vary the relative size of effects from the non-Abelian nature of the Ising anyons.We can interpolate between a purely Abelian anyon model (γ σ c = 0) and one with non-Abelian charges, as well as varying the relative strengths of processes that contribute to the non-Abelian nature of the anyon dynamics.In particular, setting γ h very high would allow braiding processes to easily perform non-local unitary gates on the fusion space (the hallmark of a non-Abelian anyon system), while setting p d to be very high suppresses some effects of having using the fusion-aware clustering decoder.These plots both use the Ising Fusion Code.We see that the performance of the clustering-type decoders is poorer than the PMA-type decoder for the same noise models, though they still exhibit a threshold.
a non-trivial fusion space (as discussed in Appendix B).Thus by studying the variation of the threshold value for our system under different noise channels, we can infer the relative effects (if any) that the non-Abelian nature of our anyons has on this quantity.
To this end, we consider the following three distinct scenarios to compare to that studied in Sec.VI A 1. We use the IFC and the PMA decoder to study these noise channels.

ψ-only pair creation
We first consider the case where γ ψ c = 0, all other γ = 0, and p d = 0.This situation disallows the presence of σ anyons, and so never allows a non-trivial fusion space to arise.The resulting model is completely analogous to the toric code under the bit-flip channel, and this correspondence will be made explicit in Sec.VI E.
The data for this noise model is shown in Fig. 8a.We estimate the threshold for this system as t * ψ ≈ 0.13.This is around half of the value obtained when both the σ and ψ anyons are created, as might be expected from similar study of e.g. the toric code under the bit-flip channel compared to the combined bit and phase-flip channel.The fact that the sum of the γ rates are always normalized to one means that by introducing two almost independent types of errors, we can preserve our quantum code space for around twice as many time steps.

Pair creation and decoherence
We now study the system in a regime where decoherence is the dominant mechanism.Explicitly, we have γ ψ c = γ σ c = 0, all other γ = 0, and p d = 1.The results are shown in Fig. 8b.We find a threshold of t * dec ≈ 0.25, consistent with values where no decoherence occurs.The fact that local charge decoherence has little effect on information storage supports the idea that Metropolis sam- pling accurately models thermalization despite the fact that it imposes complete charge decoherence at every time step, which may not be true in real thermalization processes.

Pair creation dominated by hopping
Finally, we study a regime where instead of pair-creation, hopping is the dominant mechanism for transporting charge around the lattice.We take γ h γ ψ c = γ σ c = 0, all other γ = 0, and p d = 0. We find a threshold of t * h ≈ 0.25, again consistent with the situation presented in Fig. 5.We conclude that the threshold values are largely insensitive to variation of the γ rate parameters, except in predictable ways such as the threshold decreasing by up to a factor of 1  2 as γ ψ c γ σ c varies away from 1.We do not present data for non-zero values of γ e , as any braiding process can also be viewed as a higher-order process consisting of many hopping operations, and we again find the threshold values to be insensitive to variation in this parameter.

D. Metropolis sampling
Here we consider the Metropolis sampling mechanism described in Sec.II B for the Ising Topological Code.In this case, the rates of the various noise mechanism are not fixed, but instead are given by a rate equation that depends on the state of the system and a temperature.We have set the mass of both the fermion and anyon to m = 1, which fixes the energy units.We estimate the threshold t * at a fixed temperature T by varying the lattice size L and finding the intersection point as in the previous sections.By repeating for different temperatures, we obtained the memory phase diagram of Fig. 9.
We find that the memory lifetime scales exponentially with the inverse temperature.We can explain this observation by assuming that the memory corruption is dominated by pair creation, and that other noise processes are largely irrelevant.This assumption is justified by the fact that the other noise processes (exchange, hopping, fusion) can only occur on an edge neighbouring an anyon, and that the anyon density is relatively low.Under this assumption, the Metropolis rule will simply cause the creation process to be rejected with probability e −2m/k B T since it involves a 2m increase of the energy.This has the effect of reducing the error rate, or equivalently increasing the memory lifetime by e 2m/k B T .Since the infinite-temperature Metropolis noise model is equivalent to a fixed-rate model, we can use the lifetime t * ITC ≈ 0.245 obtained in the fixed-rate model to predict the critical time t * as a function of temperature to be t * = 0.245e 2m/k B T , in excellent agreement with the numerical results Fig. 9.We expect this behaviour to saturate at sufficiently low temperature where processes that so not change the energy will be greatly favoured.These processes tend to delocalize defects and we believe that these will lead to more failures.

E. Threshold analysis
In comparing our results to existing topological code thresholds (e.g. with the toric code [33,34,57,60,61]), one should note that we use a slightly different scale on the x-axis of our plots than is conventional.In previous analyses, the x-axis has typically been represented by some form of iid noise strength.Since the non-Abelian nature of our noise processes means that the noise model cannot obviously be recast as an iid error model, we instead use the related measure of average  error operations per edge.As noted previously, our fixed rate sampling (see Sec. II A) simulations proceed by applying a number T of error operations where T is drawn from a Poisson distribution of mean T 0 , and we define the average error operations per edge, or simulation time, as t sim ≡ T 0 |E| = t as our measure of noise strength.
We can directly compare this value to the equivalent iid noise parameter in some cases where both are equally valid.Specifically, consider the toric code under the iid bitflip error channel with strength p iid .Comparing this with our noise model for any choice of rates such that γ d = 0, we can directly compute a mapping between the Poisson and iid binomial noise models to obtain an equivalence of the parameters given by Thus we have for t sim 1 that p iid ≈ t sim .This error model coincides with that of our Ising anyons in the case where σ anyons are excluded.Since ψ anyons are self-inverse Abelian anyons, they share all the features which are expected to affect the threshold with one species of toric code anyons.For this reason, we can directly compare the thresholds we find for our model in the ψ-only case to those for the toric code (note that without σ anyons present, the decoherence mechanism is trivial, and the pair-creation and hopping operations are identical).
As is seen in Fig. 8a, we find a threshold for the IFC of around t * ψ ≈ 0.13.This corresponds to an iid threshold of p ψ iid ≈ 0.11.This may seem surprising as the theoretical maximum achievable threshold for the toric code under this error model is estimated as 0.1094 [30].We believe that the discrepancy in this case lies in an additional finite-size effect which is not present in the toric code, as well as statistical errors that contribute to the uncertainty in our threshold estimates.The lowest-order logical operators for the IFC correspond to tunneling processes between code sites.These operators have fixed end points, in comparison to the logical operators for the toric code which consist of any closed loops.It is reasonable to suspect that this fact may affect the likelihood of a logical operator occurring, particularly on finite size lattices.Indeed, we obtained initial numerical estimates of these effects though a simple percolation model, and found evidence that they indeed contribute to an overestimate of the threshold.A precise and rigorous accounting of all systematic and statistical error sources would allow for a more precise determination of the thresholds for our codes and decoders, but is beyond the scope of this paper.
Note that these extra sources of finite-size scaling effects do not apply to the Ising Topological Code, and so we expect that the threshold values quoted above for the ITC apply in the thermodynamic limit.

VII. DISCUSSION
Our results demonstrate that error correction is possible in non-Abelian anyon systems and indicate that the non-Abelian nature of the excitations do not have a significant effect on the memory lifetime of stored quantum information.In fact, the various values of the error correction thresholds were remarkably insensitive to the details of the noise processes, including the details of the noise rate equation.
Our results have application to both topological quantum memories and topological quantum computing schemes based on systems with Ising anyon excitations, but more broadly they demonstrate the feasibility in more general 2-dimensional topologically ordered systems.Any such scheme of sufficient size would require some form of error correction protocol in order to succeed.We show how an appropriate error correction scheme may be implemented and how it may behave under a range of noise models.
In so doing, we have constructed several topological codes of interest.These codes are to our knowledge the first non-additive codes to be numerically simulated, and as such we develop a number of new tools to the study of quantum error-correcting codes and in particular topological quantum codes.Most notably, the decoding algorithms we describe may be more broadly applicable to other similar decoding problems.
Based on our numerics, whether or not the system rapidly decoheres in the charge basis plays little role in the memory lifetime of these systems.This can be understood as follows: although the correspondence between the decoding problem and percolation in topological codes is not exact [62], there is an intuitive correlation between percolation events and decoding failures in these systems.That is, for logical errors to occur, anyonic charge must be transported along a non-local path.Decoherence processes do not affect this intuitive picture at all, and so we might naively expect that they will not significantly alter the nature of the logical noise channel.
The fact that decoherence does not significantly affect the nature of the error-correction problem in our Ising anyon system indicates that the semiclassical noise model used in our Metropolis sampling simulations is justified, and the restriction of the noise model in this case is not expected to significantly affect the results.
Our work lends itself to extension in several directions.Firstly, it would be interesting to apply a similar analysis to other numerically tractable anyon models and/or to concrete microscopic Hamiltonians with anyonic excitations.During preparation of this manuscript, some work in this direction has been done [63].It is also an open problem as to how our error correction scheme may be leveraged into a demonstration of fault tolerant error correction, which must be the true goal of any theoretical error correction protocol.
At higher temperatures or longer timescales, this approximation may produce deviations from the true physics.One example of a process whose outcomes will be simulated incorrectly by the semiclassical model is as follows: If charge measurements were made at times t 1 , each fusion product has equal probability of being observed as I or ψ.If this measurement is made, then made again at time t 2 , the fusion outcomes are again equally likely to be observed as I or ψ.
However, if the measurement is not made at time t 1 , and the superposition of charges is allowed to persist, then at time t 2 the fusion outcomes measured will deterministically be ψ.
Thus we see that using a semiclassical simulation that collapses local charge superpositions at each timestep may misclassify some phenomena that would otherwise occur.A similar situation can be constructed for other anyon models, such as that studied in [63].For this reason, we find it necessary to perform a full simulation of the fusion space to ensure that the true behavior of braiding anyons is captured.Fortunately, the results of Sec.VI demonstrate that the appearance and value of an error-correction threshold are relatively insensitive to these details, and so a semiclassical simulation may give results that are close to the true threshold values under reasonable circumstances, a situation which is not obvious a priori.

FIG. 2 .
FIG.2.A non-planar graph G (shown with gray edges) and # G (shown with solid black edges).The edges along the "seams" of the torus are shown as dashed lines.As in the planar case of Fig.1, we also show the deformation of several edges into linearized braid sequences.

FIG. 4 .
FIG.4.The Ising fusion code uses a square graph with spherical boundary conditions.Four preferred code sites are highlighted.These sites were chosen to maximize the distance between them (and hence to maximize the code distance).We label them c t , c r , c b , and c l as shown.
FIG. 5. a) Logical |0 state failure probability versus simulation time (average errors per edge) for γ ψ c = γ σ c using the PMA decoder on the Ising Fusion Code.b) Logical |+ state failure probability for the same decoder/code/noise combination.Both quasi-thresholds are equal to within sampling and finite-size errors.

FIG. 6 .
FIG.6.Failure probability versus simulation time (average errors per edge) for γ ψ c = γ σ c using the PMA decoder on the Ising Topological Code.
FIG. 7. a) Logical |0 state failure probability versus simulation time (average errors per edge) for γ ψ c = γ σ c using the simple clustering decoder.b) Logical |0 state failure probability versus simulation time for γ ψ c = γ σ c

FIG. 9 .
FIG. 9. Memory lifetime t * scales exponentially with the inverse temperature.The linear fit is t men = 0.245e 2m/k B T .Inset: Thermal simulations used to determine the memory lifetime t * ITC ≈ 0.245 of the system at temperature T = 50m/k B .Similar simulations (with L = 16 & 12) were preformed for different values of T = 50, 20, 10, 5, 1, 0.75, and 0.5 to determine the memory lifetime as a function of temperature.