Qutrit Magic State Distillation

Magic state distillation (MSD) is a purification protocol that plays a central role in fault tolerant quantum computation. Repeated iteration of the steps of a MSD protocol, generates pure single non-stabilizer states, or magic states, from multiple copies of a mixed resource state using stabilizer operations only. Thus mixed resource states promote the stabilizer operations to full universality. Magic state distillation was introduced for qubit-based quantum computation, but little has been known concerning MSD in higher dimensional qudit-based computation. Here, we describe a general approach for studying MSD in higher dimensions. We use it to investigate the features of a qutrit MSD protocol based on the 5-qutrit stabilizer code. We show that this protocol distills non-stabilizer magic states, and identify two types of states, that are attractors of this iteration map. Finally, we show how these states may be converted, via stabilizer circuits alone, into a state suitable for state injected implementation of a non-Clifford phase gate, enabling non-Clifford unitary computation.


I. INTRODUCTION
Quantum computers hold the promise of solving certain computational tasks at an exponentially faster rate than is currently believed to be possible with classical computers [1]. One of the main obstacles making the task of building a quantum computer difficult is due to quantum decoherence, where errors can, if not corrected, accumulate and spread rapidly. The theory of quantum fault-tolerance (FT) [2] provides a means to protect the coherent quantum state and allow a reliable quantum computation. This can be achieved provided the physical error rate is below a certain threshold value. The exact value of the threshold error depends on the fault-tolerance scheme adopted for the computation.
All FT schemes have a limited set of operations allowing direct implementation on the encoded quantum information [3], and in most known schemes these are the stabilizer operations. Stabilizer operations consist of a family of unitary circuits known as the Clifford group, preparation of 0⟩ state and measurement in the computation basis. Although these operations can produce highly entangled states the Gottesman-Knill theorem [4] tells us that a computation consisting of the stabilizer operations alone can be efficiently simulated by a classical computer. The stabilizer state operations are not quantum universal; in fact, Clifford circuit operations alone cannot even implement non-linear classical logic gates, such as the Toffoli gate. Nevertheless, in most quantum FT schemes the stabilizer operations are those which can be most readily achieved fault tolerantly. It is natural then to ask, what additional resources are needed to promote the stabilizer operations to universality?
One answer to this question is given by the theory of magic state distillation (MSD) [5,6]. Almost any single qubit gate is approximately universal for SU (2) [7], and similarly, if the stabilizer operations are augmented by a supply of many copies of almost any ancillary pure state, the states can be used * Electronic address: hussain.anwar.09@ucl.ac.uk as a resource for the implementation of non-Clifford gates via "state-injection" methods [6]. A supply of qubits in a particular state may suffice, provided it is not a stabilizer state.
To compute fault tolerantly, we require a scheme that tolerates preparation noise for these ancillary resources. Hence, we are interested in which mixed quantum states provide a suitable resource. One can straight away identify a set of states which will not be useful in this regard. These are the single qubit states that can be prepared via the stabilizer operations, together with classical randomness. These states are the convex hull of the Pauli eigenstates, and we shall call such states the stabilizer states. In the Bloch picture they occupy an octahedron, as shown in Fig. (1).
Despite these obstacles, Bravyi and Kitaev [5] showed that some mixed non-stabilizer states can enable universal quantum computing. They proposed a process of "magic state distillation" whereby suitable resource states are efficiently converted, using only stabilizer operations, into a smaller number of purer non-stabilizer states, the so-called magic states. While their methods are extremely powerful, they only considered the qubit case. For non-generic noise models, it is also possible to distil certain 3-qubit states that can be used to implement the Toffoli gate [8]. Here we tackle the problem of distilling pure 3-dimensional, or qutrit, non-stabilizer states, and herein refer to all pure non-stabilizer states as magic states.
Known magic state distillation schemes have an iterative structure [6,9], with each iterate having 3 steps: 1. Initialization, prepare n−copies of the qudit resource state ρ r ; 2. Projection, measure Pauli-observables that stabilize a d-dimensional subspace and postselect on the +1 outcome; 3. Decoding, perform a Clifford unitary that maps the ddimensional subspace onto a single physical qudit ρ out r .
When successful, the output state ρ out r is used as one of the inputs on the next level of iteration. All known magic state distillation protocols achieve higher purities via iteration, although non-iterative protocols are an interesting possibility; for example, compare with the hashing protocol [10,11] and quantum polar codes [12] used in analogous context of entanglement distillation. If there exists a protocol that iteratively, or by other means, converts a resource state ρ r into a magic state of arbitrarily high purity, then ρ r is said to be resource state.
Magic state distillation for qubits has an elegant geometrical visualisation in terms of the Bloch sphere. In this representation, the stabilizer states form an octahedron, whose vertices are the Pauli eigenstates. Clifford Group unitaries coincide with symmetries of this octahedron. The magic states distilled by the protocol of [6] correspond to the states invariant under certain rotational symmetries. The first of these is the set of 180 ○ rotations around the edges of the octahedron, which contains the Hadamard gate, and we shall call Hadamard-type, or H-type rotations. The second type are 120 ○ rotations around the faces of the octahedron, known as T-type rotations.
There has been a considerable amount of work to improve the original qubit schemes in [6] and to modify their noise model. Reichardt [13] showed how all the states above the edges of the octahedron can be distilled by Steane's 7 qubit code. In contrast, two interesting no-go theorems for qubit and odd-dimensional systems show that not all non-stabilizer states are useful resources. Recently, Veitch et. al. [14] showed that for odd-dimensional systems there exist bound states that cannot contribute any enhancement to the computational power of the stabilizer operations. The qubit problem is more subtle. Campbell and Browne [15] showed that for any iterative protocol there will always exist undistillable qubit states above the faces of the octahedron. However, noniterative qubit protocols are still poorly understood. Furthermore, Campbell [16] introduced an activation protocol that can activate qubit states from above the octahedron face to the distillable regions. The magic state distillation protocol allows us to upper bound the error tolerance threshold [17] and is at the heart of state-of-the-art fault tolerance schemes [18]. Moreover, an experimental implementation of MSD has been FIG. 2: An overview of how different protocols in this paper are related. We consider distillation of two different classes of states using the 5-qutrit code. The output of either can be converted by 2 subsequent sub-protocols, called parity-checking and equatorialization, to produce different magic states. The output of equatorialization can then be exploited to perform a non-Clifford unitary. The relevant section is noted in the top-right corner of each protocol box.
recently demonstrated in an NMR system for the 5−qubit code [19]. Surprisingly, very little study has been made of whether magic state distillation can be generalized to higher dimensions. There are a number of reasons why such a study would be important. In particular, topological quantum systems have anyons with braiding statistics and a dimensionality that is fixed by the physical system. Research into understanding these systems is ongoing [21,25], but it is entirely plausible that the most promising systems provide braiding statistics corresponding to Clifford operations of dimensions greater than 2. Furthermore, the original magic state schemes have some surprising features (hence the name "magic states"), and it is unclear how much these depend on special features of the qubit Clifford group, or whether they are generic. There are therefore many open questions. Can the region of states which converge under distillation be identified? A study of higher dimensional MSD may give new insights into the structural differences between the Clifford group in two-dimensions, which plays a central role in quantum computation theory, and in higher-dimensions, where it has been intensively studied in the context of SIC-POVMS [22], MUBs [23,24] and DWF (discrete Wigner functions) [25].
Recently, expanding on their work in [26], van Dam and Howard [27] have studied noise thresholds with ddimensional quantum systems. They have calculated the set of robust qudit states that are most resilient to depolarizing noise and found the degree of noise needed to map such states to the set of stabilizer states increases with dimension scaling with d (d + 1). Thus, the higher dimensional states have the potential to offer higher magic state distillation thresholds. Until now, however, no magic state distillation schemes for d > 2 had been presented.
In this paper, we present the first generalization of magic state distillation to higher dimensional systems. We demonstrate that magic state distillation can be achieved in a qutrit system, and find both similarities and differences with previously known MSD protocols for qubits. We have chosen to focus on the three-dimensional qutrit space for this first study since it has the benefit of prime dimension, and computational tractability. Nevertheless, many of the features which we uncover are likely to be generic.
We start in Sec. II by defining notation and discussing the state space and the Clifford group in higher dimensions. In Sec. III we present a generalised approach to study the distillation properties of any stabilizer code of any prime dimension. The main protocol we consider is a qutrit generalization of the 5-qubit code covered in Sec. IV and App. B. Here we find magic state distillation occurring for 2 distinct families of states. The first family contains a pair of eigenstates of the qutrit Hadamard operator. Under depolarizing noise, they are distilled up to an error threshold of 23.3%. This family is generic in that all quantum states can be mapped into the Hadamard plane by random application of the Hadamard unitaries, a process known as Hadamard-twirling. The second family contain four eigenstates of the qutrit Hadamardsquared operator, but lies within a degenerate eigenspace of this operator, and so is not uniquely defined by it. Under depolarizing noise, they are distilled up to an error threshold of 34.5%.
Not all magic states, such as those output by the 5-qutrit code, are known to be useful for state injected non-Clifford gates. However, certain special magic states, which we call phase-states, do have the capacity to implement non-Clifford unitaries. We show that, by executing additional protocols, the outputs from the 5-qutrit code can be converted into the desirable phase-states. We overview the whole process in Fig. (2), which shows that the promotion of the Clifford group can be achieved by following 5-qutrit distillation by parity-checking and equatorialization, which we describe in Sec. V.

II. DEFINITIONS AND NOTATION
A quantum dit, or a qudit, is a d−level quantum system where d ≥ 2. The corresponding quantum state ρ can be represented by a positive semidefinite d × d matrix of unit trace. In the language of linear algebra, such a matrix can be decomposed as a linear sum of a basis set. For our purposes we have adopted the Weyl basis (also known as the Heisenberg-Weyl basis). This is because the Weyl basis set is a natural generalization of the conventional qubit Pauli operators which makes them very convenient in studying MSD. Herein we assume that d is a prime dimension.

A. Higher Dimensions and Stabilizer Codes
The d−dimensional single-qudit X and Z Pauli operators are defined as [28]: (1) where ω = e 2πi d is the d-th root of unity. From this definition we see that X and Z are traceless non-Hermitian unitaries. They obey the commutation relation XZ = ω −1 ZX and are cyclic in the power of d (i.e X d = Z d = I). We define a slightly different form of the single-qudit Pauli operators as follows: where c = (1−d) 2. The main reason for choosing this definition and the significance of the extra phase ω cjk will become apparent in section II B. But we shall first outline some of the properties of the Pauli operators based on this definition. For a single qudit, the composition of two Pauli operators can easily be verified to be For the case of a composite system of n−qudits we use the symplectic notation to represents the n−fold tensor products of Pauli operators where j and k are vectors in Z n d . The Pauli operators satisfy a generalized commutation relation where k.j ′ − j.k ′ is the symplectic inner product. Based on the above definitions, the n−qudit Pauli group P n is All element of P n have eigenvalues of the form ω m for m ∈ Z d . We will now briefly review stabilizer formalism and the basic properties of stabilizer codes [29]. A stabilizer state ψ⟩ is a simultaneous +1 eigenvector of an Abelian Pauli subgroup S n ∈ P n called the stabilizer group. The subgroup S n is generated by n independent and mutually commuting generators ⟨g 1 , g 2 , . . . , g n−k ⟩ with S n = d n and ω m I ∉ S n for all non-zero m ∈ Z d . The stabilizer group S n forms a code, called the stabilizer code, with the stabilizer states being the codewords in the code-space. Stabilizer codes encodes k logical qudits (d k −dimensional Hilbert space H k ) into a larger Hilbert space H n of n physical qudits. We denote such code by [[n, k, δ]] d , where δ is the distance of the code and the subscript d is the dimension of the Hilbert space. Finally, we represent the logical operators on the code subspace by σ L j,k , such that σ L j,k ∈ P n S n and commutes with all elements of S n .

B. Qudit Space and Pauli Group Orbits
When we discuss the general structure of a magic state distillation in Sec. V A we will see that it is most convenient to use σ j,k as the basis set to represent a qudit state ρ. But recall that σ j,k is a non-Hermitian unitary operator. To guarantee the Hermiticity of ρ we must therefore impose the Hermiticity condition ρ = ρ † . We begin by expressing ρ as where the summation is over all pair elements of Z 2 d and we have assumed that α 0,0 = 1. We will use σ 0,0 and the conventional 1 1 interchangeably. We refer to α j,k as the Bloch components, as a generalization of the qubit convention, and α as the Bloch vector which has the Bloch components as its elements. Observe that for d = 2, the set of σ j,k is Hermitian and the Bloch components will be real, but in the general qudit case the Bloch components are complex, constrained by a Hermiticity relation in order for ρ = ρ † to hold. To work out the Bloch components' relation we start by explicitly writing ρ = ρ † : The importance of the extra phase factor of ω cjk in our definition in Eq. (2) is to ensure that σ † j,k = σ −j,−k , which can be easily verified. Using this fact, and after relabelling, Eq. (8) reduces to A direct implication of Eq. (9) is that only half the Bloch components, or (d 2 − 1) 2, are independent, as the other half are simply the complex conjugates. Hence, only half the Bloch components are needed to define the density operator. Of course, the independent Bloch components are complex and we still have (d 2 − 1) real parameters defining the density operator. A Bloch component α j,k can be evaluated using the following relation: Beyond the qubit case, it is not possible to visualise the entire state space with a geometrical picture similar to the Bloch sphere. However, there have been some attempts to study the geometry of the state space in higher dimensions [31], and for a Bloch type representation for qutrits [32,33]. Having fixed normalization and Hermiticity, a mixed qutrit state is described by a complex vector α ∈ C (d 2 −1) 2 . This complex vector space has a natural inner product, ⟨α, β⟩ = ∑ j α * j β j , a norm α = ⟨α, α⟩, and distance D(α, β) = α − β . These geometric concepts are related to the density matrix representation via Since all pure state satisfy tr(ρ 2 α ) ≤ 1, this entails and all physical states are within a Bloch-like ball of radius (d − 1) 2 about the origin, with pure states on the surface of the Bloch ball. The qubit state space is of course a special case in which all the points on the surface of the Bloch sphere corresponds to positive physical states. The additional condition required to ensure a positive pure state, as shown in [34], is tr(ρ 2 ) = tr(ρ 3 ) = 1.
The geometry vectors α corresponding to physical states ρ α is quite intricate. However, within certain hyperplanes the substructure is very simple. Consider hyperplanes defined by a set of d positive orthonormal operators, {ρ β j }, such that tr(ρ β j ρ β k ) = δ j,k . Note that in the geometric picture, this entails We consider the complex hyperplane spanned by {β j }, such that α = ∑ j b j β j . It follows that an operator ρ α is positive if and only if ∑ j b j ≤ 1 and b j = b j for all j. Hence, the physical α lie within the convex polytope with β j as vertices.
We have d vertices all equally separated from each other and residing within a real d − 1 dimensional hyperplane. Hence, the polytope has the structure of a standard simplex. For d = 3, and an appropriate plane, the physical states reside within an equilateral triangle. Finally, we shall discuss the orbits of the Pauli group P n when acting on a general qudit state ρ(α) with conjugation being the group action. The singular orbit of a general state ρ(α), denoted by Orb(ρ(α)), is defined as In our Bloch representation we are using Pauli group elements as a basis set for the states, thus conjugation by Pauli operators will not transform the basis elements, but will add a phase of the form ω l for some l ∈ Z d . Therefore, the overall effect of this conjugation is to add certain phases to the Bloch components. The exact form of the phases is exactly given by: where σ † j ′ ,k ′ = σ −j ′ ,−k ′ , the commutation and composition relations where used in the last step.

C. Qutrits
So far the discussion has been for all prime dimensions, but in the remainder of this paper we will discuss the qutrit case only. Therefore, we shall outline some of above results  (2) is σ j,k = ω −jk X j Z k , where ω = e 2πi 3 and c = −1. The explicit qutrit ρ(α) state is: As we can see, completely specifying a qutrit state would only require 4 complex independent parameters. In terms of the Bloch components, the purity condition tr(ρ 2 ) = 1 for a general qutrit state can be shown to be α ≤ 1.
The Pauli group orbits for a single qutrit state can be evaluated using Eq. (16). We are interested in knowing how the phases of the four independent Bloch components change when an element from the 9 qutrit σ j,k operators is conjugated with the general qutrit state. The result is summarised in table I. We refer to these phases as the orbital Bloch phases. These represent the phases which generate the set of states Pauli equivalent to any state. The magic states that we will find are unique up to a Bloch orbital phase. In other words, inserting one of the phases from the set in table I into the Bloch components of the magic states would also give a valid magic state with the same distillation properties.

D. Stabilizer Dynamics and the Clifford Group
One the most important advantages of the stabilizer formalism is the simplicity it provides when studying the dynamics of the stabilizer states. Instead of studying the action of a unitary U on a n−qudit stabilizer state ψ s ⟩ (which would require the complete description of the map on d n parameters), the problem can be reduced to studying the evolution of the stabilizer operators by U in the Heisenberg picture (which is linear in n) [35]. An important class of unitaries are those that map the elements of the stabilizer group back to the stabilizer group under conjugation. These operation are called the Clifford unitaries. More formally, the Clifford group is the normalizer of the stabilizer group, defined as For any prime dimension, the Clifford group has been shown to be generated by three gates [28,36]. These are the Hadamard gate H, the Phase gate S and the Controlled-NOT gate Λ(X). The d−dimensional Hadamard gate is given by: Under conjugation, the Hadamard gate transforms Pauli operators as For completeness, the d−dimensional S and Λ(X) gates are given by: In addition to unitary dynamics we allow for so-called Pauli measurements. For qubit systems the Pauli group contains Hermitian operators where for qudit they are unitary but not Hermitian, with the exception of the identity. Rather, when we say we measure a non-trivial Pauli P this is taken to mean that we perform a POVM measurement that has the eigenvectors of P as POVM elements.

III. MAGIC STATE DISTILLATION
Using the definitions and notations we have developed in the previous section, we will show how the three steps of a MSD protocol described in the introduction can be formulated to study the distillation properties of any stabilizer code of any prime dimension.
Resource state preparation: The computational model considered when studying MSD consists of perfect stabilizer operations and the ability to prepare n identical copies of a noisy resource state ρ r . It should be noted that the resource state is noisy due to the imperfect preparation process. By repeating the preparation procedure n times, the state ρ ⊗n r will be prepared. As an input to the MSD protocol we consider a general state By performing the remaining steps of the iteration on the above general form, we will determine the map on the Bloch components of the initial general state ρ. Then, by searching the state space for different initial states, we can identify the resource states as those that when used as an input to the protocol the output state has a higher fidelity with respect to a pure non-stabilizer state, and ultimately distilling this nonstabilizer pure state. If the search is done systematically, one can in principle identify the entire region of resource states ρ r .
Stabilizers measurement and Decoding: The (n − k) stabilizer generators of a stabilizer code [[n, k, δ]] d are measured successively along with postselecting on the +1 outcome of each measurement. That is, if one of the outcomes is ω k (for some non-zero k ∈ Z d ) then the protocol is aborted, and the procedure is repeated with a fresh state ρ ⊗n . Also, it is important to notice that the error correction code is not being used for the usual purpose of correcting errors since the syndrome measurements are performed on the product state ρ ⊗n . If successful, the measurement of the stabilizers simply project the state to the code's subspace. The projector operator describing this measurement procedure can be put into a convenient form to us as After the measurements, the n copies of ρ will be projected into the code's subspace, and the following map will be performed: The state is decoded via a Clifford operator [9]. In a Heisenberg picture, the decoding operation maps logical operators on the code-space to unencoded operators acting on a single qudit. The output Bloch components after decoding α out j,k therefore corresponding to the components of the logical operators prior to decoding σ L j,k . After one round of the distillation, these can be evaluated as follows: The resultant expressions for the output Bloch components will be multi-variable complex polynomials of order n. For the qutrit codes we consider here, we have not found analytic solutions for the fixed points of the map. However, the problem is tractable by using numerical methods to study the distillation behaviours and the fixed points to a high accuracy.

IV. THE 5-QUTRIT CODE
Using the generalized formulation of MSD in the previous section we have studied the distillation properties of the five qutrit code. The stabilizer generators of the general five qudit code [ [5,1,3]] d takes the same form in all dimensions. It is usually presented in terms of the conventional generalized Pauli operators of Eq. (1), as shown in table II.  [30].
Based on these stabilizers we have studied the distillation map of Eq. (B5) for the four Bloch components of a general input qutrit state. The exact distillation calculation and the expressions of the distillation map are given in appendix B.
The decoding of a stabilizer code is not unique, but one of an equivalence class of unitaries, a coset of the Clifford group, which are all equally valid choices. The choice of decoding will affect the iterative distillation behaviour. The decoding specified by the logical operators in table II is the canonical one, though we found behaviour was simplified by following each iterate with the following additional Clifford unitary: where this maps a Hermitian operator ρ α such that: Without this corrective Clifford one observes a cycling behaviour throughout the distillation process, see App A. We identify two qualitatively different families of states. Firstly those in the Hadamard plane, which satisfy HρH † = ρ. Secondly, we investigate distillation of an interesting set of states outside the Hadamard plane.

A. Hadamard-like Distillation
In the qubit case, the eigenstates of the Hadamard gate are known to be magic states, distillable by the five qubit code [[5, 1, 3]] 2 [37]. Since the exact generalized form of the Hadamard gate is defined in Eq. (19), a good starting point would be to investigate whether the qutrit Hadamard eigenstates can be distilled by [[5, 1, 3]] 3 . We begin by outlining some of the structural properties of the qutrit Hadamard eigenspace. In the matrix representation the qutrit Hadamard is given by where ω 2πi 3 and it has the eigenvalues (+1, −1, i). We label the corresponding three eigenstates as ( The density operators of the eigenstates have the form: where a = 1 4 1 + √ 3 , b = 1 4 1 − √ 3 and c = − 1 2 , are real parameters. This basis of pure states all lie on the hyperplane of operators of the ρ(x, x, y, y). Probabilistic mixtures of these states form an equilateral triangle and as reviewed earlier, all points outside this triangle correspond to non-physical operators. Furthermore, any qutrit state can be projected onto the Hadamard plane by applying the following twirling operation to each copy of the input state ρ: This twirling operation maps the general Bloch components as α 1,1 and α 1,2 ↦ Re(α 1,1 + α 1,2 ) 2 .
As such, we are interested in studying the distillable regions within the Hadamard plane, i.e. the states corresponding to the points inside the triangle.
In studying the distillable region in the Hadamard plane, it is also informative to chart out regions for which distillation is impossible by any protocol. Clearly, all stabilizer states are undistillable (red region in Fig. (3i)), but the results of Veitch et al [14] prove undistillability of all qutrit states with a positive Wigner function. The numerically calculated positiveregion is shown in Fig. (3i)  The path of the distillation takes the form shown in Fig.  (3ii) where we have chosen the H + ⟩ blue triangle as an example. The small black points are few examples of input states to the distillation, and the black lines represent the distillation paths toward the H + ⟩ state. Notice how the distillation does not follow a straight line (along the magic axes) as in the qubit case. In fact, in analogy to the qubit case, the above plane is the Hadamard magic plane. The curved distillation path can be understood by studying how the noise of a resource state in the Hadamard plane is suppressed in different directions by the distillation. We start by considering a general state ρ △ inside the triangle of the form with 1 + 2 ≤ 1. For clarity, we write the Bloch components of the above state as ρ △ (A △ , B △ , C △ , D △ ). They can be calculated explicitly using Eq. (10) as Since we know the general distillation map for any set of Bloch components (see App. B), we can simply substitute the above expressions into Eqs. (B11-B14) to evaluate the output Bloch components α out = (A out △ , B out △ , C out △ , D out △ ). The output state is then ρ out △ = ρ(α out ). We have numerically calculated the output out 1 and out 2 to the first order term as follows: The above expressions shows an asymmetric error suppression in the 1 (along the H + ⟩ − H − ⟩ line) and 2 (along the H+⟩− H i ⟩ line) directions. The particular distillation paths of Fig. (3ii) can be explained by observing the difference in the coefficients of out ( 1 , 0) and out (0, 2 ), where we see that in the distillation region of the H + ⟩ state there is a stronger attraction toward the H i ⟩ state compared to the H − ⟩ state. The above analysis shows that the performance of the [[5, 1, 3]] 3 code in distilling the qutrit Hadamard states is not as good as the qubit case where the 15 qubit code by [6] has an output error probability of out ≈ 35 3 . This is to be expected given the similar performance of the five qubit code [37] in distilling the H-type qubit magic states.
The state H i ⟩ is not distillable by [ [5,1,3]] 3 . In fact, this state belong to the family of states with maximally nonpositive Wigner function [27]. As we can see this state is the furthest away from the stabilizer region in the Hadamard plane and to bring it to the stabilizer region would require a depo- To improve the size of the distillation region we have investigated a qutrit version of the seven qubit code [[7, 1, 3]] 2 distillation proposed by [13]. We start with the stabilizer generators of [[7, 1, 3]] 2 code and by adding the (−1) power to the appropriate X and Z Pauli operators, we constructed a set of generalized 7−qudit commuting stabilizer generators as shown in table III. We repeated the distillation procedure for this set of generators for the case d = 3 (exact calculations are omitted here) and we have investigated its distillation capability in the Hadamard plane. We found that this code attracts towards the non-stabilizer segments of the line joining the H + ⟩ and H − ⟩ states with the distillation region enclosed by the green triangle in Fig. (3). In other words, the 7−qutrit code distils not pure, but mixed states. Regardless, the protocol may be useful for bringing states into the region distillable by the 5-qutrit code. The distillation path for the H + ⟩ state is shown in Fig. (3iii). This code increase the distillation region as shown by the solid blue curve in Fig. (3). For example, a state between the solid blue line and the dashed blue triangle is first distilled by the seven qutrit code to a state within the dashed blue triangle, after which the [[5, 1, 3]] 3 code is used to distil the H ± ⟩ states.

B. Hadamard-squared subspace
In this section, we introduce a second class of magic states distilled by the [[5, 1, 3]] 3 code. The state in question is an eigenstate of the H 2 operator, but lies within a degenerate eigenspace for this operator, and so is not uniquely defined by it. The magic state considered here has the form: where up-to 4 decimal places The equality of the 1⟩ and 2⟩ components follows from the H 2 symmetry. In the Bloch representation, the α vector is, ii) The success probability for the trivial syndrome measurements for the case where the magic states ϕ⟩ and H ± ⟩ are undergoing depolarizing noise.
with the feature that all components are real being again related to the H 2 symmetry. Our numerical analysis shows that these states are clearly attractor fixed points of the distillation protocol, but we do not have a closed form analytic expressions for them. As a consequence, we will not be able to determine analytically how the error is suppressed as we did in the previous section for the Hadamard magic states. Nevertheless, we can still gain a numerical indication of how the error of ϕ⟩ states is suppressed. Lets start with a ϕ⟩ state undergoing depolarizing noise: For a sufficiently small the distilled state ρ out dep will also be of the above form (i.e on the depolarizing axis). We can then calculate the output error probability as follows: In general, we expect that out ≈ n for very small . Therefore, the power n can be evaluated as the gradient of a log − log plot of out versus , as shown in Fig. (4i), and we found that n ≈ 1.
For completeness, we include the success probability p s of the syndrome measurements. Successful syndrome measurements, where all outputs of the stabilizer measurements is +1, are described by the projector Π given in Eq. (B15). Hence, the probability of this measurement is simply: which is given in Eq. (B16) for all sets of Bloch components. We have computed p s for both ϕ⟩ and H ± ⟩ undergoing depolarizing noise as an input states to the distillation. A plot of p s is given in Fig. (4ii).

V. PROMOTING THE CLIFFORD GROUP
An interesting practical problem is how to use the magic states, H + ⟩ and ϕ⟩, or their Clifford equivalent states, to perform injection of a non-Clifford gate. While their utility is not immediately apparent in their current form, additional subprotocols can be used to prepare the phase-states that are useful for gate injection. The phase-states are so-called because they only hold phase information, having the following form: In Sec. V C, we show how to use these states for gate injection. However, first we describe, in Sec V A, the parity-checking protocol, which is used to convert both H + ⟩ and ϕ⟩ into the plus-state Ψ + ⟩ = 0⟩ + 1⟩. These plus-states are then input to the equatorialization procedure, in Sec. V B, which finally outputs a desired phase-state. For an overview of how these protocols fit together the reader may refer back to Fig. (2).

A. The parity-checker protocol
Here we introduce a simple qutrit distillation protocol that is very efficient against a specific type of noise, but vulnerable against another type of noise. However, both H + ⟩ and ϕ⟩ have zero overlap with the "bad" noise term and so the protocol can be efficiently used to convert these states into a plus-state. Before beginning the iterative protocol some manipulation of the input states H + ⟩ and ϕ⟩ is required: where S = 0⟩⟨0 + 1⟩⟨1 + ω 2⟩⟨2 . For eigenstates of H 2 the first step is not strictly required, but it is listed to increase the generality of the protocol. In particular, this preparation procedure maps all quantum states to where Ψ ± ⟩ = ( 0⟩± 1⟩) √ 2. Note that Ψ − ⟩ is Clifford equivalent to the Hadamard eigenstate H i ⟩. For imperfect H + ⟩ and ϕ⟩ states, with depolarizing noise , we have where for the two magic states of interest c( H + ⟩) = 0.2113 and c( ϕ⟩) = 0.0152. The parity-checker protocol will exponentially suppress the value of η 0 , whereas δ 0 will linearly increase. However, this is not problematic as δ 0 can be made arbitrarily small via distillation by the 5-qutrit protocol. The iterative parity-checker is now fairly simple, on the (n+ 1)th round we have 1. Take two copies of ρ(δ n , η n ); 2. Measure the observable Z 1 Z † 2 and postselect on +1; 3. Decode the state such that j, j⟩ → j⟩; 4. Use the output state ρ(δ n+1 , η n+1 ) as an input in the next iterate.
It is straightforward to verify the iterative relations are where p n is the success probability When the small noise component is zero, so δ 0 = 0, and the large noise is not too large, η 0 < 1 3, then η n vanishes exponentially quickly such that η n ∼ (2η 0 ) n . Allowing for non-zero δ 0 , the protocol can be iterated for approximately n ∼ log(δ 0 ) rounds before the δ noise becomes problematic. Let us consider a concrete example. If we have a ϕ⟩ magic state with depolarization noise = 10 −8 , this is first prepared into a noisy plus-state with η 0 ∼ 0.0152 and δ 0 = 10 −8 3. The total noise, η n +δ n , will decease for the first 3 rounds of parity checking. After the 4th round we have a plus-state with a total error of only 2.707 × 10 −8 . This illustrates that high-fidelity plus-states can be prepared from high fidelity ϕ⟩ states in a small number of rounds.

B. Equatorialization
Here we describe a simple magic state protocol that converts the plus-states to the phase-states. The phase-states lie on a generalization of the qubit Bloch sphere equator, hence the term Equatorialization. This protocol is probabilistic but not iterative. We take two highly purified copies of a plusstate, Ψ + ⟩. We measure a 2-qutrit stabilizer operator and postselect such that we project onto the subspace spanned by: and then decode onto a single qutrit. When successful this produces the following transformation: Noticing that ω + ω 2 = −1, we find that the output is a phasestate Finally, we have a magic state of the desired form.
Each unitary occurs with equal probability because the phasestate contains no variation in amplitudes. Although random, a desired unitary can eventually be reached by repeated attempts, with each attempt corresponding to a step of a random walk on a manifold of phase gates with a toroidal topology. For the Φ 0,π ⟩ state, the corresponding unitaries take a simple form. The closure of U k,0,π gives a group of order 4, up to to a global phase, which is composed of U k,0,π and the identity. For such a small group any desired unitary will be reached, with high probability, within a small number of attempts.
For brevity, let us herein denote N = U 0,0,π . Clearly, N is non-Clifford and so we label the group C +N as the single qutrit group generated by N and the Clifford group. The Gottesman-Knill theorem no longer applies and so we know no method of simulating computations using these unitaries. However, does this provide a dense cover of SU(3)? We discuss this question in App. C, where we show that the group C +N is of infinite order. The size of the group, and that it contains basis changing gates of the Clifford group, make it highly plausible that it provides a dense cover of SU(3). However, we presently do not have a complete proof.

VI. SUMMARY
In this paper we have demonstrated that magic state distillation protocols do exist for higher dimensional systems. We have provided a generic formulation that can be used to study the distillation capabilities of any stabilizer code for any prime dimension. We have shown that the five qutrit code [ [5,1,3]] 3 is capable of distilling the qutrit Hadamard eigenstates, implying that, in analogy to the qubit case, there exist H-type qutrit magic states. Under depolarizing noise the H⟩ ± are distillable up to a noise threshold of 23.3%. This is a higher threshold than the threshold achieved if the five qubit code were used to distil the qubit Hadamard eigenstates, which is We have also shown how to convert the outputs of the 5qutrit code into other magic states. In turn, these have been shown to provide a unitary gate that promotes the Clifford group. We have offered good, but inconclusive, evidence that this promoted Clifford group enables universal quantum computing.
Finally, we list some few open questions that merit future investigation. Do analogous magic state distillation protocols work similarly in arbitrary dimensions (or maybe just prime dimensions)? Are the magic states listed here exhaustive (could there be other attractive fixed points)? Can we find a closed form for the H 2 eigenstate coefficients? Can we better understand the relationship between a code and the magic states associated with the code? Is the set of unitaries (Clifford + (1,1,-1)) approximately universal? Are there topological quantum computing models which realise all or part of the qudit Clifford group?
in Fig. (5), such that for the distillation map for one iterate, E, performs E( ϕ j ⟩) = ϕ j+1 ⟩ and ϕ⟩ = ϕ 1 ⟩ = ϕ 5 ⟩. The 4-cycling states are related by; However, by considering E ′ (ρ) = RE ′ (ρ)R † , this cycling behaviour vanishes. Note also, that this cycling behaviour is not only seen for the pure states but for depolarized states, and so all of these states are distilled by the 5-qutrit code. the case of evaluating the output Bloch component α out 1,0 , Eq. (B5) becomes α out 1,0 = tr(ρ ⊗5 Πσ L −1,0 ) tr(ρ ⊗5 Π) , with tr(ρ ⊗5 Πσ L −1,0 ) given in Eq. (B17). We have evaluated the expressions for the four output Bloch components. However, writing them out in terms of α j,k notation is cumbersome. Therefore, for clarity, we will relabel the four qutrit Bloch components as follows (α 1,0 , α 0,1 , α 1,1 , α 1,2 ) ≡ (A, B, C, D). For example, tr(ρ ⊗5 Π) is given in Eq. (B18), where the subscript r represent the number of the distillation rounds with r = 0 corresponding to the initial input state. Furthermore, it can be shown that the resultant expressions for the four output Bloch components can compactly be expressed in terms of a single function. This function is given in Eq. (B19), and based on this function the distillation map can be expressed as (B10) These four expression represents the complete distillation map, as the remaining four components are simply the complex conjugates of these expressions. However, the above expressions do not incorporate the additional corrective Clifford (see Sec. IV), and will result to the cycling behaviour in App. A. We need to ensure that the map in Eq. (28) is applied after every iteration. This can easily be achieved in our formalism by the appropriate relabelling as follows: A r+1 =F(D * , C, B * , A), (B11) B r+1 =F(C * , D * , A, B), (B12) C r+1 =F(D, B * , C, A), (B13) D r+1 =F(C * , A * , D, B * ), which is the corrected distillation map. We can see that in order to calculate the fixed points of this map analytically, one would have to solve the above simultaneous complex multivariable polynomials of order 5. It is known from the famous Abel-Ruffini theorem that there is no algebraic solution for a general polynomial of order five or above. Therefore, the best way to discover the fixed points of the distillation is through numerical means. We started with initial states ρ(A, B, C, D) for certain Bloch components and computed the above expressions for a number of iterations, and observed whether there is a convergence toward a fixed point.