Generalized Toric Codes Coupled to Thermal Baths

We have studied the dynamics of a generalized toric code based on qudits at finite temperature by finding the master equation coupling the code's degrees of freedom to a thermal bath. As a consequence, we find that for qutrits new types of anyons and thermal processes appear that are forbidden for qubits. These include creation, annihilation and diffusion throughout the system code. It is possible to solve the master equation in a short-time regime and find expressions for the decay rates as a function of the dimension $d$ of the qudits. Although we provide an explicit proof that the system relax to the Gibbs state for arbitrary qudits, we also prove that above a certain crossing temperature, qutrits initial decay rate is smaller than the original case for qubits. Surprisingly this behavior only happens with qutrits and not with other qudits with $d>3$.


Introduction
It is known that the fragility of quantum states in the presence of interaction with an environment represents the main challenge for the large-scale implementation of quantum information devices in quantum computation and communication. Quantum error correction is the theoretical method that was devised to protect a quantum memory or communication channel from external noise [1][2][3][4][5][6][7][8]. In these quantum error correction schemes, to improve the stability of quantum information processing, the logical qubits should be implemented in manyparticle systems, typically N physical spins per logical qubit. This is the quantum version of the classical method based on encoding information by repetition or redundancy of logical bits in terms of physical bits [9,10]. The logical qubits should be stable objects with efficient methods of state preparation, measurements and application of gates. By efficiency we mean a certain scaling behavior, e.g. the lifetime of a logical qubit should grow with N .
In order to implement fault-tolerant methods for quantum information processing, we need to find a physical system with good enough properties to accomplish this protection from a noisy environment and decoherence. One promising candidate is topological orders in strongly correlated systems. Here, the ground state is a degenerate manifold of states whose degeneracy depends on the topological properties of a certain lattice of qubits embedded into a surface with nontrivial topology [11]. Many-body interacting terms in a Hamiltonian are responsible for the existence of this topological degeneracy. The logical qubits are stored in global properties of the system represented by nontrivial homological cycles of the surface. In these topological codes, 4 enough, while with qubits in 2D lattice models the topological protection is lost under the action of thermal fluctuations [57], it is however possible to set up a fully fledged topological quantum computation using certain types of TCCs in higher-dimensional lattices [58]. Under these conditions, it is possible to prove that self-correcting quantum computation, including state preparation, quantum gates and measurement can be carried out in the presence of the disturbing thermal noise. Additionally, note that thermal noise does not always turn out to be detrimental in quantum information, even for systems without topological order [59,60].
In this work, we extend those results regarding the thermal effects on generalized toric codes constructed out of qudits. Here in we summarize briefly some of our main results.
1. We formulate the dynamics of a generalized toric code based on qudits at finite temperature. To this end, we find the master equation coupling the qudits of the system code to a thermal bath. 2. We study and classify the different types of thermal processes that may occur when the anyonic excitations are created, annihilated or diffused throughout the system. In particular, we find that for qutrits new types of anyons and thermal processes appear that are forbidden for qubits. 3. The master equation is too involved so as to yield an explicit expression for the decay rate of the topological order initially present in the code. However, in a short-time regime it is possible to solve it and find expressions for the decay rates as a function of the dimension d of the qudits. Interestingly enough, we find that the decay rate for qutrits presents a crossover temperature T c that is absent for any other qudits. 4. We can give an explicit proof that for long enough times, the non-local order parameter representing the topological order in the system decays to zero.
This paper is organized as follows. In section 2, we review the formulation of the master equation of the 2D Kitaev code for qubits in order to establish the notation and the necessary tools to study thermal effects in more general toric codes. We also introduce a non-local order parameter and study the fate of topological orders for two different regimes: the shorttime regime and the long-time regime. In section 3, we find the master equation describing topological qutrits coupled to a thermal bath. This allows us to see new energy processes for the anyonic excitations that are not present when the toric code is made up of qubits. Likewise, the short-time regime has a different behavior that can be seen in the initial decay rate of the topological order. In particular, we can define a crossover temperature for qutrits where the decay rate is better than that with other qudits. Section 4 presents the conclusions. See appendix A for the evolution of the order parameter for qutrits and appendix B for a proof of the irreducibility of the computational representation of the d-Pauli group needed to study the master equation in the long-time regime.

Thermal stability of the Kitaev two-dimensional (2D) model
We shall not dwell upon the details of Kitaev's toric code [4]; however, we will introduce the basic ideas to understand how to apply a thermal stability analysis to it, as well as to establish the notation and methods. We will consider a k × k square lattice embedded in a 2-torus. Let us attach a qubit, such as a spin 1/2, to each edge of the lattice. So we have N = 2k 2 qubits. For each vertex s and each face p, we denote the stabilizer operators in the following form: where X j and Z j are the Pauli matrices applied to the qubit on site j. A s and B p commute among each other for they have either 0 or 2 common edges. They are also Hermitian and have eigenvalues 1 and −1 (see figure 1). Therefore, they constitute an Abelian subgroup of the Pauli group of n qubits that is a stabilizer group. Let H be the Hilbert space of all n = 2k 2 qubits and define the topological quantum code or protected subspace C ⊆ H as follows: This construction defines a quantum code called the toric code. The operators A s and B p are the stabilizer operators of this code, i.e. operators that leave trivially invariant the code space. As we want to analyze the physical properties of this code, in particular the thermal properties of the topological order, it is convenient to define its associated Hamiltonian in the form Complete diagonalization of this Hamiltonian is possible since operators A s , B p commute. In particular, the ground state coincides with the protected subspace of the code C; it is fourfold degenerate (see figure 2). All excited states are separated by an energy gap E 4. This is due to the fact that the difference between the eigenvalues of A s (B p ) is equal to 2. Excitations come in pairs since they correspond to violations of the plaquette and/or vertex stabilizer operators and these must comply with the overall constraints s A s = 1 and p B p = 1. Thus, excitations are represented as open strings in the direct or the dual lattice of the original square lattice. An essential feature of this Hamiltonian is its locality in terms of four-body interactions, very useful for practical purposes. Another key property is that this Hamiltonian model is gapped, which led to the initial expectation that all types of 'errors', i.e. noise-induced excitations, will be removed automatically by some relaxation processes. Of course, this requires cooling, i.e. some coupling to a thermal bath with low temperature (in addition to the Hamiltonian (3)), as we shall describe later on. It can be shown that this Hamiltonian is robust under local quantum perturbations at zero temperature [57]: there would be a level splitting that will vanish as exp(−ak), where k is the length of the lattice [4].
Due to this unavoidable coupling to a thermal bath, our system is subject to thermal errors as well. These can be seen as violations on the plaquette and vertex conditions A s | = | , B p | = | . Moreover, A s and B p are unitary, and also Hermitian in the case of qubits. Therefore, violations on the plaquette and/or vertex condition are given by for a certain number of sites s and/or plaquettes p. These violations cost energy to our system, thereby becoming excitations. As long as they always come in pairs (to satisfy the conditions s A s = 1 and p B p = 1), they can be seen (pictorially) as string operators with plaquette or vertex violations at the ends. Errors on the system can be expressed in terms of operators σ x , σ z or products among themselves. These operators act on each edge j where the physical qubits are placed. We use the notation σ x for a Pauli operator of type X when it refers to an error, i.e. a bump operator acting due to the coupling to the thermal bath. Likewise with σ z . It is just a matter of notation to distinguish when we have an operator that defines our stabilizer operators in A s , B p and when we have an error acting on the system. To see what effect they produce, we will see how the ground state changes by applying these σ x,z . We will see that this corresponds to the creation, annihilation and movement of a pair of excitations, which from now on we shall refer to as anyons. These are called anyons since their wave function picks up a different phase than fermions or bosons when we exchange the end-particles of string operators of x-type with ztype. According to this notation, when we apply a bump operator from the thermal bath, it will 7 act on the ground state of the system as follows: where | is the ground state of the system where our information is encoded. This means that the physical qubit at the edge j has been bumped. The energy cost will be E = 4 in energy units of the system corresponding to the definition of H sys . As a first step, one is interested in designing a stable quantum memory, i.e. an N -particle system that can support at least a single encoded logical qubit for a long time, preferably with this time growing exponentially with N . This is the notion of stability we shall refer to from now on. In the paper by Alicki et al [50], they provide a rigorous method to prove thermal instability of the 2D Kitaev model and obtain a master equation that describes the dynamics of the system weakly coupled to a thermal environment. We will study the problem of thermal instability within the framework of topological orders obtaining complementary and interesting results.

Davies' formalism
Let us consider a small and finite system that is coupled to one or more heat baths at the same inverse temperature β = (k B T ) −1 , leading to the total Hamiltonian Here H sys represents the Hamiltonian of the system where the quantum information is encoded and which we want to protect from the external thermal noise. H bath is the bath Hamiltonian, i.e. it describes the internal dynamics of the bath that is out of our control. Finally, V represents the coupling between the system and the thermal bath. S α and f α are operators that act on the system and bath, respectively. Both the coupling operators S α and f α are assumed to be Hermitian (without loss of generality [41]). In the weak-coupling regime that we shall assume throughout this work, the Fourier transformĝ α of the auto-correlation function of f α plays an important role, as it describes the rate at which the coupling is able to transfer energy between the bath and the system [39][40][41][42]. Often a minimal coupling to the bath is chosen, minimal in the sense that the interaction part of the Hamiltonian is as simple as possible but still addresses all energy levels of the system Hamiltonian in order to have an ergodic reduced dynamics. This last condition is ensured if [41,[43][44][45][46] i.e. no system operator apart from those proportional to the identity commutes with all the S α and H sys . The weak-coupling limit [39][40][41][42] results in a Markovian evolution for the system given in the Heisenberg picture by the master equation The generator of the evolution G(X ) is a sum of two terms: the first is a usual Liouville-von Neumann term as in the quantum mechanics of closed systems, while the second is a particular type of the Kossakowski-Lindblad generator: Here the S α (ω) are the Fourier components of S α as it evolves under the system Hamiltonian where the ωs are the Bohr frequencies of the system Hamiltonian (hω = E 1 − E 2 , for two energy levels E 1 and E 2 ). In addition, the temperature of the environment appears in (11) through β, and this generator is the so-called Davies' generator [39] or the Born-Markov generator in the quantum optics literature.

Master equation for the 2D Kitaev model with qubits
Given the simplicity of Kitaev's model, we can apply Davies' theory for studying its stability in the presence of thermal noise. This is represented pictorially in figure 3.

9
The interaction Hamiltonian is assumed to be local and associated with σ x and σ z errors: where f x j and f z j are associated with two different baths. Thus, first of all, we need to compute the Fourier transform of the system operators e it H sys σ x j e −it H sys and e it H sys σ z j e −it H sys in order to define the dynamical operators of the system. Here H sys : Thus, stabilizers A s only play a role in the Fourier transform of σ x j and B p only in σ z j . By computing this Fourier transform, we obtain the dynamical operators of the system due to the coupling to the thermal bath. With = 4 denoting the gap of the Toric code Hamiltonian, then the expressions for these operators S α (ω) that appear in equation (11) are as follows [50].
and the projectors: . These operators have a nice interpretation in terms of anyonic properties of the system.
, where R(4) and R(0) are the exchange rates between the system and the bath associated with each Bohr frequency, namely ω = 0, 4, assuming units of J = 1.

Topological order
We shall study the evolution of the expectation value GS| X c |GS as a simple order parameter, where X c is the tensor product of σ x Pauli operators along one non-contractible loop on the surface of the torus and |GS denotes a generic ground state of the system Hamiltonian.
This ground state is a superposition of the degenerate states in the ground state manifold of H sys , namely C. This gives us a sufficient measure of the topological order of the system [49]. If this quantity falls to zero during the time evolution for every element of C, there is not a global and self-protected way to encode quantum information. The evolution of the operator X c is given by equation (47), In order to simplify the computation, we remove the free evolution by performing the transformationX Since the dissipator is invariant under this transformation, we obtain Interestingly, for the expectation value we obtain GS| X c (t) |GS = GS|X c (t) |GS , as |GS is an eigenstate of H sys . Taking into account expressions (14) and (15), the action of the dissipators on X c can be simplified to and where we have used the fact that [P j ±,0 , X c ] = 0 for every j, as these projectors are only functions of vertex operators. However, the same assertion is not true for R j ±,0 in general. If j ∈ c, i.e. j does not belong to the path where X c is acting on, every element commutes with each other and their contribution is zero. On the other hand, if j ∈ c, as σ z j σ x j σ z j = −σ x j , the string operator yields σ z j X c σ z j = −X c . Therefore, simplifying we obtain where |c| is the number of points in the path c.

Short-time regime
The solution to the master equation (18) is formally written asX c (t) = e L(t) X c . However, this expression is too involved to be computed analytically except for short and long times to be specified hereby. In the first case, at lowest order we havẽ 12 The evolution of GS| X c (t) |GS is given by To arrive at this equation, we have used the fact that for all j: whereas for L z , we have Finally, as GS| X c (t) |GS = GS|X c (t) |GS , the desired equation valid at short times is It is important to remark that R(0) does not appear in the initial decay rate, as long as short times are concerned. The diffusion of anyons is a second-order process in time as it requires first the creation of a pair of anyons with R( ), and later free diffusion with R(0).

Long-time regime
On the other hand, in order to analyze the thermal properties for long times, we write the Davies generator in the Schrödinger picture through the relation Tr[L † (ρ)X ] = Tr[ρL(X )] for any X and ρ. It is a well-known result [39][40][41][42] that the Gibbs state is a stationary state for L † , where ρ β = e −β H sys /Z , β is the same to the inverse temperature as the surrounding bath, and Z := Tr(e −β H sys ) is the system partition function. To guarantee that any initial state of the system relaxes to ρ β , we can resort to condition (7). In our case this follows from Schur's lemma as S α = σ x j , σ z j and {1, σ x , σ z , σ x σ z } form an irreducible representation of the Pauli group.
Thus GS| X c (t) |GS Tr[X c ρ β ] for large t, and we have Tr[X c ρ β ] = 0. This is simply due to the fact that ρ β is diagonal in any of the possible eigenbases of H sys , and it is not difficult to choose one such that X c vanish on diagonal elements, for some eigenbases {|ψ i } of H 0 , Kitaev's Hamiltonian.
In conclusion, whatever the initial value of the order parameter X c (0) , it decays to zero during the time evolution of the system, provided that the temperature is finite. The decay rate at short times is equal to 2 |c|R( )e − β . Note the detrimental effect of the factor |c|: the larger size of the system, the higher the decay rate. In order to keep the order parameter above a certain finite value such that X c (0) = 0, this decay rate must decrease, which is not the case when increasing the system size.

The Kitaev 2D model for qudits
In this section, we consider again a 2D toric code, but instead of assuming that we have a twolevel system on each site, we will consider that particles arranged on the torus have d accessible levels. We will first derive a general theory for qudits and then consider the case d = 3 (qutrits). A qutrit can be represented, for instance, as a particle of spin 1 or a three-level system in an atom, etc.
This problem is very interesting since qutrits have certain advantages with respect to qubits.
1. Qutrits have a larger capacity for information storage.
2. Quantum channels are more robust for qutrits. For example, Bell inequalities are proved with more accurate bounds. This is relevant for quantum key distribution. 3. Entanglement quantum distillation is more efficient with qutrits than with qubits [61]. 4. Qutrit logic gates [62] are also capable of providing universal quantum computation, i.e. the necessary computational power to construct all possible logic gates [8].
To build a system like that, we will try to choose the Hamiltonian and the operators acting on the system in the same way as before. Previously, for two-level systems, we have considered the Pauli matrix algebra to be the basis of operators in our system. Now, we have to use a proper generalization for dimension d. As iX Z = Y gives the second Pauli matrix, it is enough to consider X and Z in this generalization to quantum states with d multilevels. However, the generalization of Pauli matrices to dimension d is not unique 2 . Thus, we shall select the most important properties of Pauli matrices of dimension 2 for our purpose of quantum error correction.
In d = 2, we defined a basis: |0 , |1 in the Fock space of each particle. They are defined as the eingenstates of the Z Pauli matrix. And the X Pauli matrix takes |0 to |1 and vice versa.
The key important properties of these matrices for doing error correction are the following.
• They satisfy a cyclic condition (i.e. applying twice Z or X Pauli matrices is the identity), i.e. they are unitary.
• They anticommute, which means X Z = −Z X .
Those are the properties that are generalized to the d-dimensional case. Hermiticity is not taken into account as a basic ingredient, as we can always add the Hermitian conjugate obtaining a Hermitian operator, e.g.Z = Z + Z † , thenZ is Hermitian. Now we consider a basis for the particle Fock space: |0 , |1 , . . . , |d − 1 , which will be the eingenvectors of the generalized Z matrix with a certain eingenvalue. We define X as the operator which takes the state |0 to |1 , then |1 to |2 and so on. We will also ask for a cyclic condition as in the previous case: All these requirements can be cast on to the following defining relations: Looking at equation (29) we can deduce the meaning of operators X and Z . X is the displacement operator in the computational basis (i.e. in the Fock space basis of the physical qudits). Z is the dual operator of X under a discrete Fourier transform. In other words, Z is diagonal in the computational basis and its eingenvalues are the weights of the Fourier transform. Thus, X plays the role of the displacement operator and Z is the dual operator on a system with discrete states of qudits [8].
Due to the cyclic condition (28) of Z (Z d = 1), we have the relation ω d = 1 where, in general, ω is a complex number. This implies that ω is a primitive d-root of unity, Additionally, we can easily verify that Z X = ωX Z , as follows from equation (29).
We have already the algebra of operators that we are going to use in order to built the stabilizer operators on this qudit toric code. The problem is that if we construct the vertex and plaquette operators as before, namely, then [A s , B p ] = 0 for all s and p. They commute with each other provided that they do not share any common edge, but that is not the case if they share two. This happens because in this case the operators X and Z are no longer Hermitian. As shown in figure 7, we have Hermitian operators. We need to think of another way to define our operators to have the same commutation rules as before, and this leads to defining an orientation on the lattice. This is shown in figure 7. Defining an orientation on the lattice is a direct consequence of the non-Hermiticity of operators X and Z . Using the orientation of the lattice, we define the stabilizer operators in the following way. To build the vertex operators A s we assign an operator X or X −1 depending on the arrows of the edges of the lattice. If an arrow is pointing towards the vertex j, we will use X −1 j to build A s , and if the arrow is pointing out another vertex, we use X j . For plaquette operators B p , Z k is taken if the arrow is pointing clockwise and Z −1 k for anti-clockwise, as shown in figure 7. To see now that we obtain the correct commutation rule, we look again at figure 7 and check, Then, the Hamiltonian could be written as follows: Although, according to the definition of A s and B p , this operator is unitary, it is important to note that the operators A s and B p are no longer Hermitian, so H aux is no longer Hermitian. However, we may redefine the Hamiltonian in the following way: where H sys is Hermitian now. The effect that H † aux has in the system is a redefinition of the orientation on the lattice. So we have a superposition of a lattice orientated in the way of figure 7 (arrows up and right) and another with arrows down and left. Nevertheless, one can always think in terms of H aux for the pictorial image and then use H sys to compute energies and derive equations.

The anyon model
The theory developed above was done for the general case of qudits. From now on and to be concrete concerning thermal effects, we will focus on the case where d = 3 (qutrits). Later on we will be able to extract conclusions for qudits as well. There are still many important aspects to be studied about this model and its coupling to a thermal bath. We need to compute the energy gap of the Hamiltonian, i.e. the energy difference between the ground state where the code lies and the excited states which represent the errors. It is also important to calculate the anyon statistics, as long as they are associated with the excitations of a topological system with qutrits.
At d = 3, the phase factors are ω = e i(2π/3) , ω 2 = e i(4π/3) , ω 3 = 1. We will see, for this particular case, how excitations can be created, moved and annihilated. This will give us the properties of the anyon model that is going to be associated with the group Z 3 .
As before, we use a notation in which σ z j = Z j and σ x j = X j , except that we use the symbol σ to denote errors acting on the system, i.e. bump operators acting because of the coupling to the thermal bath, whereas we shall use X, Z for the Hamiltonian interactions defined by the vertex and plaquette operators of H sys . Errors on the system can be expressed in terms of operators σ x , σ z or products containing them, and acting on each edge j where the qutrits are placed. And the same goes for σ z . To see what effect these errors have on the system, we will see how the ground state changes by applying σ x,z . We will see that this corresponds to processes in which anyons are created, annihilated or moved throughout the torus.
Let us see what happens when we bump a qutrit in a position j from the outside and then act with the Hamiltonian H aux , H aux σ z j |ψ . Note that every operator of the Hamiltonian commutes with this σ z j except two A s operators which share a leg with this qubit j. But, contrary to the case of d = 2, there is an orientation defined on the lattice. So, for instance, if an error (σ z j ) occurs in a certain vertical edge, one of these A s (the one below) is defined with an X j , thus: but the A s above the edge is defined with X −1 j , then: Hence, we have two violations of the vertex condition, one with charge ω and the other with ω 2 . This is one of the two types of anyons that we will have in this system, and we shall denote it as an ω 2 -ω anyon. It is important to point out that these are only labels to classify the excitations based on the violations of the operator A s (and B p ). In principle, we could classify anyons based on the violation of stabilizers A −1 s (and B −1 p ) that appears in H † aux . It is just a matter of labeling; the physics is the same. Now we can act with σ z j again and obtain the other anyon type called ω-ω 2 . Actually, they could be considered as the same anyon type as before but with opposite orientation. However, it is convenient to define them as two types of anyons as they will have different braiding properties. Moving anyons of the same type around each other will be different from the case of having anyons of different types. Likewise, it will be necessary to have anyons of different types in order to have fusion of anyons without annihilation. We shall explain this in the next subsection in more detail.
Note that acting twice with σ z j is equivalent to acting with (σ z j ) −1 . Thus, although every error can be expressed in terms of X and Z operators, it will be useful to think sometimes as if we act with either X, Z or X −1 , Z −1 . All these arguments are exactly the same in the case of B p operators and σ x errors. Therefore, we have four types of anyons, two of plaquette type and two of vertex type.
Let us study now the braiding of the anyons. We will consider two chains of different types: plaquette anyon and vertex anyon (as in figure 8). In this case, we get something remarkably different from the d = 2 case. Now it is not the same to let one anyon remain still and move the other around it as it is to do it the other way around. Thus, let us move particles around each other. For example, let us move an x-type particle around a z-type particle (see figure 9). Then, because S x (c) and S z (t) cross each other on just one qutrit satisfying the relation and S x (c)|ψ x (q) = |ψ x (q) . We see that the global wave function, i.e. the state of the entire system, acquires the phase factor ω 2 . Nonetheless, if the operation is the opposite, i.e. if we move a z-type particle around a x-type particle, then since S x (q) and S z (c) cross each other just on one qutrit again satisfying the relation Z X = ωX Z (38) and S z (c)|ψ z (t) = |ψ z (t) . We see that the global wave function acquires now the phase factor ω. Figure 9. Anyons of type Z (red) on the direct lattice attached to a string t.
Anyons of type x (green) on the dual lattice associated with a string q. The x-type particle moves around a z-type particle on a closed string c.
Thus, we arrive at a very important novelty for qutrits that is different from the case when we dealt with qubits in two aspects.
1. The phase that the anyon picks up is different from −1.
2. The phase depends on the orientation in which the braiding close path is traversed.

New anyon energy processes
First of all, let us look at the gap of the Hamiltonian. We will reach our first excited state by applying a σ z or σ x operator to the ground state. Let us see what is the energy difference between the ground state and the first excited state. Remember that 2H sys = H aux + H † aux . We denote by P and S the number of plaquette and vertex operators, respectively, with P + S = N being the number of qutrits in the lattice, and {l, l } are the adjacent vertices of the site of a qutrit j: Thus, the energy difference is The action of σ x produces the same energy increment but we have to do the commutation with the operators B p .
This calculation can be easily extended to the case of qudits with arbitrary d, obtaining the gap equation Note that there is a reduction of the energy gap for d = 3 in comparison with the case of qubits, where it was 4. It is also important to point out that if we act again on the same bond of the lattice with (σ z ) −1 , there would be an energy reduction of the same amount of energy. Moreover, if at the endpoint of an anyon ωω 2 we act with σ z , we obtain the same pair of anyons again, and the same energy, but longer (see figure 10. (2)). In this process, the energy is preserved, E = 0. This means that there is no energy exchange between the thermal bath and the system. We can understand the process as a diffusion of the anyon with no energy cost. In analogy to the case d = 2, this is what is called moving an anyon. It is also important to remark that for qutrits, all processes that involve moving a simple pair of anyons still have no energy cost.
Until now, there is a complete analogy with the case of d = 2. But we are going to see now a process that only occurs at d > 2. Imagine that there have been two excitations on the system, and two anyons of opposite orientation have been created. Moreover, they are separated by just one vertex operator. The situation is plotted in figure 10. (1).
Imagine that we act now with a σ z on the bond, which is error-free, that links the anyons ω 2 -ω and ω-ω 2 (opposite orientation). Let us analyze the energy process.
so the energy difference is What has occurred is that two anyons have been tied together, but not annihilated. This process lowers the energy of the system by a smaller amount than the process of annihilation. If in this situation we would act with a (σ z j ) −1 on the point where the two pairs of anyons are tied together, the two anyons would split apart, and this process would cost energy E = 3/2. This could be analyzed exactly the same way with σ x errors and B p operators.
It is remarkable that this phenomenon cannot happen at d = 2, as at d = 2 the product ωω = (−1)(−1) = 1. Therefore, d = 3 is the first nontrivial case to have processes like these in a toric code with qudits.

Master equation for topological qutrits
As we have seen, all these processes are generated by the action of operators σ z , (σ z ) 2 and σ x , (σ x ) 2 ; as in this case, the square of the Pauli operators is their Hermitian conjugate. Nevertheless, the energy exchange depends on the situation of the system when we bump it with the thermal bath from outside. Before writing the master equation that describes the dynamics of the system, it will be useful to distinguish between these situations by local projectors. The answer to the question whether this is possible or not in this case is not trivial. However, we show that it is possible to classify into groups of processes that have the same energy gain from the bath. Furthermore, they could be distinguished by certain projection operators that only involve two adjacent vertex or plaquette operators.
We arrive at the following classification: Figure 11. (1) Initial state 1-1 ⇒ Final state ω 2 -ω. (2) Initial state ω 2 -ω ⇒ Final state ω-ω 2 . This is an example of what happens to the topological charges when there is a bump from the thermal bath outside. The first one the energy gain is E = 3. The second one E = 0.
In this table, we have represented all combinations of two adjacent topological charges. In the first column, we depict a representation of the different types of anyons, with two topological charges attached at their ends and linked by a dash. Correspondingly, all these anyons have an intrinsic orientation. At the left side of the dash there is the eigenvalue of the operator A s and at the right side, the eigenvalue of the adjacent operator A s . A physical qutrit j would be in the middle of the dash (see the example in figure 11). In the second column, we write the projector that gives 1 for that situation and 0 for the others.
Here we have defined the following operators in order to simplify the notation: α=0,+1,−1 (s, s ) : , where s and s are the two vertexes surrounding the qutrit j. The index α takes values on the exponent of the phases ω that appear from the braiding processes. These projectors tell us which charges of the system surround a certain qutrit. That is why they are local projectors. Moreover, it is easy to verify that they form a set of orthogonal projectors: (P j α ) 2 = P j α . As we have already explained, we classify the situation of the system in terms of the charges according to the eingenvalues of the operators A s associated with the part of the Hamiltonian H aux . One could do the same thing for A −1 s , but the situation of the system will be the same independently of the label we assign them. So these projectors can discriminate perfectly between eigenstates of the Hamiltonian H sys . Now, given a certain state of the system |ψ , by applying these projectors we can figure out which situation we have. This means that if an operator σ z or σ x (or their Hermitian conjugate) is going to act on our system, we will know which energy process is bound to happen. Based on this, and studying the different situations that we can encounter, one can define a set of operators that tells us whether an anyon has been moved, created, annihilated or fused when we apply the generalized Pauli operators (as we did in figure 10). This is done by analyzing the initial and the final state after the action of a bump operator and seeing which would be the energy after and before the process, as shown in figure 11. Therefore, we have: Here the upper indices of operators a j are related to the energy cost of the process.
• a (1) † j creates a pair of anyons of z-type and a (1) j annihilates it. The energy cost is E = 3. • a (2) † j and a (2) j are related to the process of fusion or separation, respectively, of anyons as in figure 10.(1) and also to the process of creation (and annihilation) of a pair of anyons tied to a previous pair. The energy cost is E = 3 2 . • a 0 j moves anyons and also it can invert the orientation of a pair of anyons (as in figure 11.(2). There is no energy cost in these processes.
For the plaquette operators B p we proceed in the same way, obtaining a similar result. The corresponding local projectors that we denote as R j are built analogously just by changing A s for B p , where p and p are the adjacent plaquettes to the qutrit j. Then the operators that describe the analogous process for x-type anyons are Some of these operators are associated with more than one projector, unlike for qubits. That is because for three-level systems, the possibilities for different excitation scenarios have grown significantly.
As we have seen in the previous section, these operators arise naturally as the Fourier transform of the interaction Hamiltonian when a thermal bath is weakly coupled with our system, In this case, the interaction Hamiltonian will be of the form and it is quite important to remark that there are only three Bohr frequencies this time, ω = 0, ± 3 2 , ±3. We can check that the dynamical operators obtained are indeed compatible with this interaction potential as α S α = α S α (ω). In our case, it is trivial to check: with n = 0, 1, 2, using equations (44) and (45). Moreover, [H, a n ] ∝ a n , based on the fact that H sys is made of stabilizers, which at most introduces a phase when they are applied to states a i |φ . Thus, A s (B p )a i |φ ∝ a i |φ and a i A s (B p ) |φ ∝ a i |φ ; therefore [H, a i ] ∝ a i , ∀a i a dynamical operator of our system. With this proviso, the Davies generator turns out to be given by with

Topological order
Similarly to the case of qubits, we will study the evolution of the expectation value GS| X c |GS , where X c is the tensor product of σ x generalized Pauli operators (d = 3) along a noncontractible loop, and |GS denotes a certain ground state in the stabilizer subspace; namely, a superposition of the degenerate states in the ground state manifold of H sys . In the weak-coupling limit, the master equation that describes the dynamics of this quantity is dX c (t) In order to simplify the calculation, we remove the free evolution part of the equatioñ both the dissipator L and the mean value GS| X c |GS being invariant under this transformation.

Short-time regime
In the short-time regime, we can approximateX c (t) (1 + tL)X c ; here we denote X c := X c (0). Thus, the evolution of GS| X c (t) |GS is We need to calculate GS| L(X c ) |GS , with L(X c ) = L x (X c ) + L z (X c ). This calculation is done in appendix A, obtaining Hence, we can define := 2 R( )e − β |c| as the initial decay rate of the system. For qutrits, = 3, while for qubits (see equation (25)) we have obtained an analogous expression but with = 4 instead.

25
This result can be generalized for the case of qudits with arbitrary d. We have already seen that, at short times, only the creation of anyons contributes to the decay of topological order. The free diffusion of anyons and the fusion processes among them will not appear as they are second order processes in time. However, as we increase d there are more types of anyons with different energies. Moreover, a pair of anyons should always be compatible with the conditions s A s = 1 and p B p = 1. That means that the possible types of anyons with different energies are of the form ω n − ω d−n with n = 1, . . . , d 2 , and respective energies n = 2(1 − cos 2π n d ). Note that n = 1 refers to the lowest energy pair of anyons, i.e. the energy gap of the Hamiltonian. Thus, the initial decay rate has to be the sum of all these contributions: It is important to point out that in the case of qudits, an analogous expression for the interaction with the environment to (13) All nontrivial powers of σ x and σ z are included to allow for excitations of physical qudits from one level to another, at first order in time.
Using equation (53) it will be possible to establish a crossover temperature T c as the limit for which the initial decay rate will be larger for qubits than for qudits. For the sake of comparison, we take R( n ) the same for qubits and qudits. This is reasonable since n are of the same order, and R( n ) are the Fourier transforms of the bath coupling that induces the excitations on the physical qudits. Thus, we set up the condition d (T c ) := 2 (T c ). Using equation (53) we arrive at the following expression: For other values of d, the initial decay rate for qudits will always be larger than for qubits. This happens as n n increases almost linearly with d, and d = 3 is the only case when this quantity is smaller than 4, i.e. the gap in the case of qubits. Let us now compute T c for qutrits: with E 0 being the natural energy unit of the system. This leads to the following crossover temperature: The meaning of this temperature is the following. Above this temperature T c , the initial decay rate for qutrits is smaller than that for qubits, something that makes qutrits better in this comparison. For E 0 ∼ 100 kHz used in the proposal of a Rydberg quantum simulator [63] for the operators of the 2D toric code, we obtain an estimate of T c ∼ 20 µK.
In addition, it could be computed a T c comparing systems with d odd and (d − 1) even. There is always a temperature above which the system of qudits with d odd has a lower initial decay rate than the previous (d − 1) even.
It is also important to point out that is only the initial decay rate. It is possible that the dynamics of anyons, with free diffusion, etc, play an important role in the loss of topological order. Beyond short times, our conjecture is that the new processes that appear in the case of qutrits, i.e. fusion of anyons that end tied up, will be an obstacle to the free diffusion of anyons. This would represent an improvement for the stability of the generalized toric code in some intermediate time regime for this is the cause of the loss of topological order in the system.

Long-time regime
Now we want to study the master equation (49) in the opposite time regime. We are interested in the fate of the non-local order parameter we are using to describe the topological order in a system of qudits in a generalized toric code. We conjecture that the final state will be given by a thermal Gibbs state. To show that our observable for the order parameter X c approaches the expectation value of X c in the Gibbs state for times long enough, we resort again to the condition (7). In the generalized case, it reads as The generalization of the toric code leads to new physics. Indeed, we have specialized for the case of qutrits and obtained very interesting results. First of all, new Abelian anyons have arisen with novel braiding properties, i.e. new statistics by exchange of particles. For instance, let us move a pair of anyons around another pair that stays still. We would pick up a different phase, letting the first pair remain still and moving the other one around. Furthermore, new energy processes appear which are forbidden for qubits, d = 3 being the first nontrivial system where these new processes can be observed. Moreover, we present a master equation that describes the dynamics of any observable of the system coupled to a thermal bath, giving a complete description of the problem.
We have proposed a new way to study thermal stability regarding the loss of topological order in the system. At short times, the system starts losing its order with a certain decay rate that we are able to compute explicitly. We have checked that the system relaxes to the thermal state for any value of d, as expected. However, we have proved that above a certain crossover temperature, the initial decay rate for qutrits is smaller than in the original case for qubits. Surprisingly, this behavior only happens with qutrits and not with other qudits with d > 3.
vanishes. On the other hand, because σ z σ x = ωσ x σ z and the cyclic property of the trace, we conclude that the character of every element of the form σ m x σ k z is zero for any representation. The rest of the terms are proportional to the identity ω n 1, and so χ(ω n 1) = ω n d.
The irreducibility criterion asserts [9,100] that a representation of a group G is irreducible if and only if the scalar product of characters is the identity, that is thus, the representation is irreducible.