Qutrit magic state distillation

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 that is 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 of protecting the coherent quantum state and allowing 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 FT 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 the |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 nonlinear 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 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 that 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 straightaway 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 figure 1.

Zoom In Zoom Out Reset image size

Despite these obstacles, Bravyi and Kitaev [5] showed that some mixed non-stabilizer states can allow 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. Although their methods are extremely powerful, they only considered the qubit case. For non-generic noise models, it is also possible to distil certain three-qubit states that can be used to implement the Toffoli gate [8]. Here we tackle the problem of distilling pure three-dimensional (3D), 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 three 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 d-dimensional 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 an 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 a resource state.

MSD for qubits has an elegant geometrical visualization 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 call these as Hadamard-type or H-type rotations. The second type is the set of 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 seven-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, non-iterative 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 FT schemes [18]. Moreover, an experimental implementation of MSD has been recently demonstrated in an NMR system for the 5-qubit code [19].

Surprisingly, very few studies have been conducted on whether magic state distillation can be generalized to higher dimensions. There is 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 to understand 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 that in higher dimensions, where it has been intensively studied in the context of symmetric, informational complete, positive operator valued measure (SIC-POVMs) [22], mutually unbiased bases (MUBs) [23, 24] and discrete Wigner functions (DWF) [25].

Recently, expanding on their work in [26], van Dam and Howard [27] studied noise thresholds with d-dimensional quantum systems. They have calculated the set of robust qudit states that are most resilient to depolarizing noise and found that 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 have 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 show both similarities and differences with previously known MSD protocols for qubits. We have chosen to focus on the 3D qutrit space for this first study since it has the benefit of prime dimension and computational tractability. Nevertheless, many of the features that we uncover are likely to be generic.

We start in section 2 by defining the notation and discussing the state space and the Clifford group in higher dimensions. In section 3, we present a generalized 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 five-qubit code covered in section 4 and appendix B. Here we find magic state distillation occurring for two 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 contains four eigenstates of the qutrit Hadamard-squared 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 five-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 five-qutrit code can be converted into the desirable phase states. We overview the whole process in figure 2, which shows that the promotion of the Clifford group can be achieved by following five-qutrit distillation by parity checking and equatorialization, which we describe in section 5.

Zoom In Zoom Out Reset image size

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 will assume that d is a prime dimension.

2.1. Higher dimensions and stabilizer codes

The d-dimensional single-qudit X and Z Pauli operators are defined as [28]

$$\begin{equation} X=\sum^{d-1}_{j=0}|{(j+1)\ {\rm mod}\ d}\rangle \langle{j}|, \qquad Z=\displaystyle\sum^{d-1}_{j=0}\omega^{j} \ |{j}\rangle \langle {j}|, \end{equation} \tag{ 1 }$$

where ω = e^2πi/d is the dth root of unity. From this definition we see that X and Z are traceless non-Hermitian unitaries. They obey the commutation relation XZ = ω⁻¹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:

$$\begin{equation} \{\sigma_{j,k}=\omega^{cjk}X^{j}Z^{k}|\forall (j,k)\in\mathbb{Z}_{d}^{2}\}, \end{equation} \tag{ 2 }$$

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 2.2. 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

$$\begin{equation} \sigma_{j,k}\sigma_{j',k'}= \omega^{j'k-c(jk'+j'k)}\sigma_{(j+j'),(k+k')}. \end{equation} \tag{ 3 }$$

In the case of a composite system of n qudits, we use the symplectic notation to represent the n-fold tensor products of Pauli operators:

$$\begin{eqnarray} \sigma_{j_{1},k_{1}}\otimes\sigma_{j_{2},k_{2}}\otimes\dots\otimes\sigma_{j_{n},k_{n}}&\equiv&\sigma_{j_{1}j_{2}\dots j_{n},k_{1}k_{2}\dots k_{n}}\nonumber\\ &\equiv&\sigma_{\boldsymbol{j},\boldsymbol{k}}, \end{eqnarray} \tag{ 4 }$$

where j and k are vectors in $\mathbb {Z}_{d}^{n}$. The Pauli operators satisfy a generalized commutation relation,

$$\begin{equation} \sigma_{\boldsymbol{j},\boldsymbol{k}}\sigma_{\boldsymbol{j}',\boldsymbol{k}'}=\omega^{\boldsymbol{k}.\boldsymbol{j}'-\boldsymbol{j}.\boldsymbol{k}'}\sigma_{\boldsymbol{j}',\boldsymbol{k}'}\sigma_{\boldsymbol{j},\boldsymbol{k}}, \end{equation} \noindent \tag{ 5 }$$

where kċj' − jċk' is the symplectic inner product. Based on the above definitions, the n-qudit Pauli group $\mathcal {P}_n$ is

$$\begin{equation} \mathcal{P}_{n}=\{\omega^{l}\sigma_{\boldsymbol{j},\boldsymbol{k}} \ | \ l \in \mathbb{Z}_{d}\}. \end{equation} \tag{ 6 }$$

All elements of $\mathcal {P}_{n}$ have eigenvalues of the form ω^m for $m\in \mathbb {Z}_d$.

We shall 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 $\mathcal {S}_{n}\in \mathcal {P}_n$ called the stabilizer group. The subgroup $\mathcal {S}_{n}$ is generated by n independent and mutually commuting generators $\left < g_{1},g_{2},\dots ,g_{n-k}\right >$ with $|\mathcal {S}_{n}|=d^{n}$ and $\omega ^{m}I\notin \mathcal {S} _{n}$ for all non-zero $m \in \mathbb{Z}_{d}$. The stabilizer group $\mathcal {S}_{n}$ forms a code, called the stabilizer code, with the stabilizer states being the code words in the code space. Stabilizer codes encode k logical qudits (d^k-dimensional Hilbert space $\mathcal {H}_{k}$) into a larger Hilbert space $\mathcal {H}_{n}$ of n physical qudits. We denote such codes by [[n,k,δ]]_d, where δ is the distance of the code and 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 $\sigma _{j,k}^{L} \in \mathcal {P}_{n}\backslash \mathcal {S}_{n}$ and commutes with all elements of $\mathcal {S}_n$.

2.2. Qudit space and Pauli group orbits

When we discuss the general structure of a magic state distillation in section 5.1, 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 follows:

$$\begin{equation} \rho(\boldsymbol\alpha)=\frac{1}{d}\displaystyle\sum_{(j,k)}\alpha_{j,k}\sigma_{j,k}, \end{equation} \tag{ 7 }$$

where the summation is over all pair elements of $\mathbb {Z}_{d}^{2}$ and we have assumed that α_0,0 = 1. We will use σ_0,0 and the conventional $\mathbbm{1}$ interchangeably. We denote α_j,k as the Bloch components, as a generalization of the qubit convention, and α as the Bloch vector that 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 ρ = ρ^†:

$$\begin{equation} \displaystyle\sum_{(j,k)}\alpha_{j,k}\sigma_{j,k} = \displaystyle\sum_{(j,k)}\alpha_{j,k}^{*}\sigma_{j,k}^{\dagger}. \end{equation} \tag{ 8 }$$

The importance of the extra phase factor of ω^cjk in our definition in equation (2) is to ensure that σ^†_j,k = σ_−j,−k, which can be easily verified. Using this fact and after relabelling, equation (8) reduces to

$$\begin{equation} \alpha_{j,k}^{*}=\alpha_{-j,-k}. \end{equation} \tag{ 9 }$$

A direct implication of equation (9) is that only half of the Bloch components, or (d² − 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² − 1) real parameters defining the density operator. A Bloch component α_j,k can be evaluated using the following relation:

$$\begin{equation} \alpha_{j,k}={\mathrm{tr}}(\rho\sigma_{j,k}^{\dagger})={\mathrm{tr}}(\rho\sigma_{-j,-k}). \end{equation} \noindent \tag{ 10 }$$

Beyond the qubit case, it is not possible to visualize 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 $\boldsymbol {\alpha } \in \mathbb {C}^{(d^2-1)/2}$. This complex vector space has a natural inner product, $\langle \boldsymbol {\alpha }, \boldsymbol {\beta } \rangle = \sum _{j}\alpha _{j}^{*}\beta _j{}$, a norm $| \boldsymbol {\alpha } | = \sqrt { \langle \boldsymbol {\alpha }, \boldsymbol {\alpha } \rangle }$ and distance D(α,β) = |α − β|. These geometric concepts are related to the density matrix representation via

$$\begin{equation} {\mathrm{tr}} ( \rho_{\boldsymbol{\alpha}}^{\dagger} \rho_{\boldsymbol{\beta}} ) = (1+2 \langle \boldsymbol{\alpha}, \boldsymbol{\beta} \rangle ) / d. \end{equation} \noindent \tag{ 11 }$$

Since all pure state satisfy tr(ρ²_α) ⩽ 1, this entails

$$\begin{equation} |\boldsymbol{\alpha }|^{2} \leqslant \frac{d-1}{2}, \end{equation} \noindent \tag{ 12 }$$

and all physical states are within a Bloch-like ball of radius $\sqrt {(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 correspond to positive physical states. The additional condition required for ensuring a positive pure state, as shown in [34], is tr(ρ²) = tr(ρ³) = 1.

The geometry vectors α corresponding to physical states ρ_α are 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

$$\begin{equation} \langle \boldsymbol{\beta}_{j} , \boldsymbol{\beta}_{k} \rangle = \textstyle \frac{1}{2}( d \delta_{j,k}-1) . \end{equation} \noindent \tag{ 13 }$$

We consider the complex hyperplane spanned by {β_j}, such that $\boldsymbol {\alpha } = \sum _{j} b_{j} \boldsymbol {\beta }_{j}$. It follows that an operator ρ_α is positive if and only if $\sum _{j} b_{j} \leqslant 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 $\mathcal {P}_{n}$ when acting on a general qudit state ρ(α) with conjugation being the group action. The orbits of a general state ρ(α), denoted by Orb(ρ(α)), are defined as

$$\begin{equation} {\mathrm{Orb}}(\rho(\boldsymbol{\alpha}))=\{\rho(\boldsymbol{\alpha '})=\sigma_{\boldsymbol{j '},\boldsymbol{k '}} \ \rho(\boldsymbol{\alpha}) \ \sigma_{\boldsymbol{j '},\boldsymbol{k '}}^{\dagger}\quad \forall \ \sigma_{\boldsymbol{j'},\boldsymbol{k'}}\in \mathcal{P}_{n}\}. \end{equation} \noindent \tag{ 14 }$$

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\in \mathbb {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 given by

$$\begin{equation} \rho(\boldsymbol{\alpha'})=\displaystyle\sum_{(\boldsymbol{j},\boldsymbol{k})}\alpha_{\boldsymbol{j},\boldsymbol{k}}\sigma_{\boldsymbol{j'},\boldsymbol{k'}}\sigma_{\boldsymbol{j},\boldsymbol{k}}\sigma_{\boldsymbol{-j'},\boldsymbol{-k'}} \end{equation} \tag{ 15 }$$

$$\begin{equation} \hspace{2.5pc}=\displaystyle\sum_{(\boldsymbol{j},\boldsymbol{k})}\omega^{\boldsymbol{j}\boldsymbol{k'}-\boldsymbol{j'}\boldsymbol{k}}\alpha_{\boldsymbol{j},\boldsymbol{k}}\sigma_{\boldsymbol{j},\boldsymbol{k}} , \end{equation} \tag{ 16 }$$

where σ^†_j',k' = σ_−j',−k', the commutation and composition relations were used to derive equation (16).

2.3. 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 the above results explicitly for the d = 3 case. Our definition of the qutrit Pauli basis set in equation (2) is σ_j,k = ω^−jkX^jZ^k, where ω = e^2πi/3 and c = −1. The explicit qutrit ρ(α) state is

$$\begin{eqnarray} &&\fl\rho(\boldsymbol\alpha)=\textstyle\frac{1}{3}(\sigma_{0,0}+\alpha_{1,0}\sigma_{1,0}+\alpha_{1,0}^*\sigma_{2,0}+\alpha_{0,1}\sigma_{0,1}\nonumber\\&&\ms \fl\hphantom{\rho(\boldsymbol\alpha)=}+\alpha_{0,1}^*\sigma_{0,2}+\alpha_{1,1}\sigma_{1,1}+\alpha_{1,1}^*\sigma_{2,2}+\alpha_{1,2}\sigma_{1,2}+\alpha_{1,2}^*\sigma_{2,1}). \end{eqnarray}\noindent \tag{ 17 }$$

As we can see, completely specifying a qutrit state would only require four complex independent parameters. In terms of the Bloch components, the purity condition tr(ρ²) = 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 equation (16). We are interested in knowing how the phases of the four independent Bloch components change when an element from the nine-qutrit σ_j,k operators is conjugated with the general qutrit state. The result is summarized in table 1. We refer to these phases as the orbital Bloch phases. These represent the phases that generate the set of Pauli states 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 1 into the Bloch components of the magic states would also give a valid magic state with the same distillation properties.

Table 1. The qutrit orbital Bloch phases.

σj',k'	σj',k'ρ(α1,0,α0,1,α1,1,α1,2)σ†j',k'
σ0,0	ρ(α1,0,α0,1,α1,1,α1,2)
σ±1,0	ρ(α1,0,ω∓1α0,1,ω∓1α1,1,ω±1α1,2)
σ0,±1	ρ(ω±1α1,0,α0,1,ω±1α1,1,ω±1α1,2)
σ±1,±1	ρ(ω±1α1,0,ω∓1α0,1,α1,1,ω∓1α1,2)
σ±1,∓1	ρ(ω∓1α1,0,ω∓1α0,1,α1,1ω±1,α1,2)

2.4. 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 an n-qudit stabilizer state |ψ_s〉 (which would require the complete description of the map on dⁿ 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 operations are called the Clifford unitaries. More formally, the Clifford group is the normalizer of the stabilizer group, defined as

$$\begin{equation} \mathcal{C}_{n}=\{C | CP_{i}C^{\dagger}=P_j \hspace{2mm}\forall \hspace{1mm}P_{i},P_{j}\in \mathcal{P}_{n}\}. \end{equation} \tag{ 18 }$$

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

$$\begin{equation} H|{j}\rangle=\frac{1}{\sqrt{d}}\displaystyle\sum_{k\in \mathbb{Z_{d}}}\omega^{jk}|{k}\rangle. \end{equation} \tag{ 19 }$$

Under conjugation, the Hadamard gate transforms Pauli operators as follows:

$$\begin{equation} H\sigma_{j,k}H^{\dagger}=\sigma_{-k,j}. \end{equation} \tag{ 20 }$$

For completeness, the d-dimensional S and Λ(X) gates are given by

$$\begin{equation} S|{j}\rangle= \omega^{\frac{j}{2}\left(j-1\right)}|{j}\rangle, \end{equation} \tag{ 21 }$$

$$\begin{equation} \Lambda(X)|{j}\rangle|{k}\rangle=|{j}\rangle |{\left(j+k\right)\,{\rm mod} \ d}\rangle. \end{equation} \tag{ 22 }$$

In addition to unitary dynamics, we allow for the so-called Pauli measurements. Although the qudit Pauli operators are not Hermitian, they are unitary operators with non-degenerate orthogonal eigenstates, and therefore can be treated as observables, where the measured quantity corresponds to a complex eigenvalue instead of a real eigenvalue.

Using the definitions and notations we have developed in the previous section, we will show how the three steps of an 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:

$$\begin{equation} \rho^{\otimes n}=\frac{1}{d^n}\displaystyle\sum_{(\boldsymbol j,\boldsymbol k)\in\mathbb{Z}_{d}^{n}}\alpha_{j_{1}\cdots j_{n},k_{1}\cdots k_{n}}\sigma_{j_{1}\cdots j_{n},k_{1}\cdots k_{n}}. \end{equation} \tag{ 23 }$$

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 input to the protocol, output a state with higher fidelity with respect to a pure non-stabilizer state, and ultimately distilling this non-stabilizer pure state. If the search is done systematically, one can in principle identify the entire region of resource states.

Stabilizer 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\in \mathbb {Z}_{\mathrm {d}}$), then the protocol is aborted, and the procedure is repeated with a fresh state ρ^⊗n. Also, it is important to note that the error correction code is not being used for the usual purpose of correcting errors since the syndrome measurements are made on the product state ρ^⊗n. If successful, the measurement of the stabilizers simply projects the state to the code's subspace. The projector operator describing this measurement procedure can be put in a convenient form to use as follows:

$$\begin{equation} \Pi=d^{k-n}\displaystyle\sum_{\boldsymbol m \in \mathbb{Z}_{d}^{n-k}}g_{1}^{m_{1}}g_{2}^{m_{2}}\cdots g_{n-k}^{m_{n-k}}. \end{equation}\noindent \tag{ 24 }$$

After the measurements, the n copies of ρ will be projected into the code's subspace, and the following mapping will be performed:

$$\begin{equation} \rho^{\otimes n}\mapsto\frac{\Pi\rho^{\otimes n}\Pi^{\dagger}}{{\mathrm{tr}}\left(\rho^{\otimes n}\Pi\right)}. \end{equation}\noindent \tag{ 25 }$$

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 correspond to the components of the logical operators prior to decoding σ^L_j,k. After one round of distillation, these can be evaluated as follows:

$$\begin{equation} \alpha_{j,k}^{\mathrm{out}}=\frac{{\mathrm{tr}}(\Pi\rho^{\otimes n}\Pi^{\dagger}(\sigma_{j,k}^{L})^{\dagger})}{{\mathrm{tr}}(\rho^{\otimes n}\Pi)}. \end{equation}\noindent \tag{ 26 }$$

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 behaviour and the fixed points to high accuracy.

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 take the same form in all dimensions. It is usually presented in terms of the conventional generalized Pauli operators of equation (1), as shown in table 2.

Based on these stabilizers we have studied the distillation map of equation (B.5) 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 is one of an equivalent 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 2 is the canonical one, although we found that the behaviour was simplified by following each iterate with an additional Clifford unitary given by:

$$\begin{equation} R = \left( \begin{array}{@{}ccc@{}} 1 & \omega & \omega \\ \omega^{2} & \omega & \omega^{2} \\ \omega^{2} & \omega^{2} & \omega \\ \end{array} \right) , \end{equation} \tag{ 27 }$$

where this maps a Hermitian operator ρ_α such that

$$\begin{equation} \{ \alpha_{1,0} , \alpha_{0,1} , \alpha_{1,1} , \alpha_{1,2} \} _{R} \mapsto\{ \alpha_{1,2}^{*} , \alpha_{1,1} , \alpha_{0,1}^{*} , \alpha_{1,0} \}. \end{equation} \tag{ 28 }$$

Without this corrective Clifford one observes cycling behaviour throughout the distillation process; see appendix 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.

4.1. 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]]₂.⁴ Since the exact generalized form of the Hadamard gate is defined in equation (19), a good starting point would be to investigate whether the qutrit Hadamard eigenstates can be distilled by [[5,1,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

$$\begin{equation} H=\frac{1}{\sqrt{3}}\left(\begin{array}{@{}ccc@{}} 1 & 1 & 1 \\ 1 & \omega & \omega^{2}\\ 1 & \omega^{2} & \omega\\ \end{array}\right), \end{equation} \tag{ 29 }$$

where ω = e^2πi/3 and it has the eigenvalues (+1,−1,i). We label the corresponding three eigenstates as (|H₊〉,|H₋〉,|H_i〉). The density operators of the eigenstates have the form

$$\begin{equation} |{H_{+1}}\rangle \langle{H_{+1}}|\equiv \rho\left(a,a,b,b\right), \end{equation} \tag{ 30 }$$

$$\begin{equation} |{H_{-1}}\rangle \langle{H_{-1}}|\equiv \rho\left(b,b,a,a\right), \end{equation} \tag{ 31 }$$

$$\begin{equation} |{H_{\rm i}}\rangle \langle{H_{\rm i}}|\equiv \rho\left(c,c,c,c\right), \end{equation} \tag{ 32 }$$

where $a= \frac {1}{4}\left (1+\sqrt {3}\right )$, $b= \frac {1}{4}\left (1-\sqrt {3}\right )$ and $c=-\frac {1}{2}$ are real parameters. This basis of pure states all lie in the hyperplane of the 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 ρ:

$$\begin{equation} \rho\mapsto\displaystyle\sum_{j=1}^{4}\frac{1}{4}H^{j}\rho H^{j^{\dagger}}. \end{equation} \tag{ 33 }$$

This twirling operation maps the general Bloch components as follows:

$$\begin{equation} \alpha_{1,0} \ \mbox{and} \ \ \alpha_{0,1}\mapsto\frac{\mbox{Re}(\alpha_{1,0}+\alpha_{0,1})}{2}, \end{equation} \tag{ 34 }$$

$$\begin{equation} \alpha_{1,1} \ \mbox{and} \ \ \alpha_{1,2}\mapsto\frac{\mbox{Re}(\alpha_{1,1}+\alpha_{1,2})}{2}. \end{equation}\noindent \tag{ 35 }$$

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 figure 3(i)), but the results of Veitch et al [14] prove undistillability of all qutrit states with a positive Wigner function. The numerically calculated positive region is shown in figure 3(i) by the yellow area.

Zoom In Zoom Out Reset image size

Running the [[5,1,3]]₃ distillation for all the remaining points as input states we have discovered that both the |H₊〉 and |H₋〉 states are distillable, but that |H_i〉 is not an attractor. The distillable regions are enclosed by blue dashed triangles. The states |H₊〉 and |H₋〉 are equally valuable as magic states because R|H_±〉∝|H_∓〉, where R is a Clifford unitary (see equation (27)), which perhaps also explains the symmetry in their distillation regions.

The path of the distillation takes the form shown in figure 3(ii), where we have chosen the |H₊〉 blue triangle as an example. The small black points are an example of input states to the distillation, and the black lines represent the distillation paths towards the |H₊〉 state. Note 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 following form:

$$\begin{equation} \rho_{\triangle}=(1-\epsilon_{1}-\epsilon_{2})|{H_{+}}\rangle\langle{H_{+}}|+\epsilon_{1}|{H_{-}}\rangle\langle{H_{-}}|+\epsilon_{2}|{H_{\rm i}}\rangle\langle{H_{\rm i}}|, \end{equation}\noindent \tag{ 36 }$$

with ₁ + ₂ ⩽ 1. For clarity, we write the Bloch components of the above state as ρ_△(A_△,B_△,C_△,D_△). They can be calculated explicitly using equation (10):

$$\begin{eqnarray} A_{\triangle}=B_{\triangle}&=&\textstyle\frac{1}{4}(1+\sqrt{3} - 2\sqrt{3}\epsilon_1 - (3+ \sqrt{3}) \epsilon_{2}),\nonumber\\ C_{\triangle}=D_{\triangle}&=&\textstyle\frac{1}{4}(1-\sqrt{3} + 2\sqrt{3}\epsilon_1 - (3- \sqrt{3}) \epsilon_{2})\!. \end{eqnarray} \tag{ 37 }$$

Since we know the general distillation map for any set of Bloch components (see appendix B), we can simply substitute the above expressions into equations (B.11)–(B.14) 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₁ and ^out₂ to the first-order term as follows:

$$\begin{equation} \epsilon_{1}^{\mathrm{out}}(\epsilon_1 ,\epsilon_2)=\langle{H_{-}}|\rho_{\triangle}^{\mathrm{out}}|{H_{-}}\rangle\approx\left(0.38+0.09\epsilon_2\right)\epsilon_1, \end{equation} \tag{ 38 }$$

$$\begin{equation} \epsilon_{2}^{\mathrm{out}}(\epsilon_1 ,\epsilon_2)=\langle{H_{\rm i}}|\rho_{\triangle}^{\mathrm{out}}|{H_{\rm i}}\rangle\approx{\left(0.77+3.55\epsilon_1\right)\epsilon_2}. \end{equation} \tag{ 39 }$$

The above expressions show an asymmetric error suppression in the ₁ (along the |H₊〉 − |H₋〉 line) and ₂ (along the |H + 〉 − |H_i〉 line) directions. The particular distillation paths of figure 3(ii) can be explained by observing the difference in the coefficients of ^out(₁,0) and ^out(0,₂), where we see that in the distillation region of the |H₊〉 state there is a stronger attraction towards the |H_i〉 state compared with the |H₋〉 state.

The above analysis shows that the performance of the [[5,1,3]]₃ code in distilling the qutrit Hadamard states is not as good as in the qubit case where the 15-qubit code of [6] has an output error probability of ^out ≈ 35³. This is as expected given the similar performance of the five-qubit code (see footnote ⁴) in distilling the H-type qubit magic states.

The state |H_i〉 is not distillable by [[5,1,3]]₃. In fact, this state belong to the family of states with a maximally non-positive 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 depolarizing noise with an error threshold of 75% (i.e. d/(d + 1) for d = 3). Whether such a state is distillable is still an open question.

To improve the size of the distillation region we have investigated a qutrit version of the seven-qubit code [[7,1,3]]₂ distillation proposed in [13]. We start with the stabilizer generators of the [[7,1,3]]₂ code and by adding the (−1) power to the appropriate X- and Z-Pauli operators, we constructed a set of generalized seven-qudit commuting stabilizer generators as shown in table 3. 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 distils states towards the non-stabilizer segments of the line joining the |H₊〉 and |H₋〉 states with the distillation region enclosed by the green triangle in figure 3. In other words, the seven-qutrit code distils not pure, but mixed states. Regardless of this, the protocol may be useful in bringing states into the region distillable by the five-qutrit code. The distillation path for the |H₊〉 state is shown in figure 3(iii). This code increases the distillation region as shown by the solid blue curve in figure 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]]₃ code is used to distil the |H_±〉 states.

Table 3. The stabilizer generators of the seven-qudit code.

g1=	I	I	I	X−1	X	X	X−1
g2=	X	I	X−1	I	X−1	I	X
g3=	I	X	X−1	I	I	X−1	X
g4=	I	I	I	Z	Z	Z	Z
g5=	Z	I	Z	I	Z	I	Z
g6=	I	Z	Z	I	I	Z	Z
XL=	X	X−1	X	X	X−1	X	X−1
ZL=	Z	Z	Z	Z	Z	Z	Z

4.2. Hadamard-squared subspace

In this section, we introduce a second class of magic states distilled by the [[5,1,3]]₃ code. The state in question is an eigenstate of the H² 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

$$\begin{equation} |{\varphi}\rangle = a |{0}\rangle + b |{1}\rangle + b |{2}\rangle, \end{equation} \tag{ 40 }$$

where up to four decimal places

$$\begin{equation} a = -0.1203-0.0272 {\rm i}, \end{equation} \tag{ 41 }$$

$$\begin{equation} b = 0.7017. \end{equation} \tag{ 42 }$$

The equality of the |1〉 and |2〉 components follows from the H² symmetry. In the Bloch representation, the α vector is

$$\begin{equation} \boldsymbol{\alpha} = \{ 0.3236, -0.4772, 0.5438, 0.6098 \} , \end{equation} \tag{ 43 }$$

with the feature that all components are real being again related to the H² symmetry. Our numerical analysis shows that these states are clearly attractor fixed points of the distillation protocol, but we do not have 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. Let us start with a |φ〉 state undergoing depolarizing noise:

$$\begin{equation} \rho_{\mathrm{dep}}=(1-\epsilon)|{\varphi}\rangle \langle{\varphi}|+\epsilon \mathbbm{1}/3. \end{equation} \tag{ 44 }$$

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:

$$\begin{equation} \epsilon^{\mathrm{out}}=1-\langle{\varphi}|\rho_{\mathrm{dep}}^{\mathrm{out}}|{\varphi}\rangle. \end{equation} \tag{ 45 }$$

In general, we expect that ^out ≈ ⁿ for very small . Therefore, the power n can be evaluated as the gradient of a log–log plot of ^out versus , as shown in figure 4(i), and we found that n ≈ 1.

Zoom In Zoom Out Reset image size

For completeness, we include the success probability p_s of the syndrome measurements. Successful syndrome measurements, where all outputs of the stabilizer measurements are +1, are described by the projector Π given in equation (B.15). Hence, the probability of this measurement is simply

$$\begin{equation} p_{\mathrm{s}}=\textrm{tr}(\Pi\rho^{\otimes 5 }\Pi)=\textrm{tr}(\rho^{\otimes 5 }\Pi), \end{equation} \tag{ 46 }$$

which is given in equation (B.16) for all sets of Bloch components. We have computed p_s for both |φ〉 and |H_±〉 undergoing depolarizing noise as input states to the distillation. A plot of p_s is given in figure 4(i).

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 sub-protocols can be used to prepare the phase states that are useful for gate injection. The phase states are called so because they only hold phase information, having the following form:

$$\begin{equation} |{\Phi_{\theta, \phi}}\rangle = |{0}\rangle + \mathrm{e}^{\mathrm{i} \theta} |{1}\rangle + \mathrm{e}^{\mathrm{i} \phi} |{2}\rangle . \end{equation} \noindent \tag{ 47 }$$

In section 5.3, we show how to use these states for gate injection. However, first we describe, in section 5.1, the parity-checking protocol, which is used to convert both |H₊〉 and |φ〉 into the plus-state |Ψ⁺〉 = |0〉 + |1〉. These plus-states are then input into the equatorialization procedure, in section 5.2, which finally outputs a desired phase state. For an overview of how these protocols fit together, see figure 2.

5.1. 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:

• 1.  
  preparation 1: uniformly randomly choose from the set of unitaries {1,H²} and apply;
• 2.  
  preparation 2: apply X^†;
• 3.  
  preparation 3: uniformly randomly choose from the set of unitaries {1,S,S²} and apply;

where S = |0〉〈0| + |1〉〈1| + ω|2〉〈2|. For eigenstates of H² 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

$$\begin{equation} \rho(\delta_{0}, \eta) = (1-\eta_{0} - \delta_{0}) {\vert \Psi^{+} \rangle \langle \Psi^{+} \vert} + \delta_{0} {\vert \Psi^{-} \rangle \langle \Psi^{-} \vert} + \eta_{0} {\vert 2 \rangle \langle 2 \vert}, \end{equation} \noindent \tag{ 48 }$$

where $|{\Psi ^{\pm }}\rangle =( |{0}\rangle \pm |{1}\rangle )/\sqrt {2}$. Note that |Ψ⁻〉 is Clifford equivalent to the Hadamard eigenstate |H_i〉. For imperfect |H₊〉 and |φ〉 states, with depolarizing noise , we have

$$\begin{equation} \eta_{0} = c(|{\psi}\rangle) + {\epsilon} / 3, \end{equation} \tag{ 49 }$$

$$\begin{equation} \delta_{0} = {\epsilon} / 3, \end{equation} \tag{ 50 }$$

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 η₀, whereas δ₀ will linearly increase. However, this is not problematic as δ₀ can be made arbitrarily small via distillation by the five-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₁Z^†₂ 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 that the iterative relations are

$$\begin{equation} \eta_{n+1} = \eta_{n}^{2} / p_{n}, \end{equation} \tag{ 51 }$$

$$\begin{equation} \delta_{n+1} = \delta_{n} (1- \eta_{n}-\delta_{n}) / p_{n}, \end{equation} \tag{ 52 }$$

where p_n is the success probability

$$\begin{equation} p_{n} = (1+\eta_{n}(3 \eta_{n} - 2))/2. \end{equation} \tag{ 53 }$$

When the small noise component is zero, so δ₀ = 0, and the large noise is not too large, η₀ < 1/3, then η_n vanishes exponentially quickly such that η_n ∼ (2η₀)ⁿ. Allowing for non-zero δ₀, the protocol can be iterated for approximately n ∼ log(δ₀) rounds before the δ noise becomes problematic.

Let us consider a concrete example. If we have a |φ〉 magic state with depolarization noise  = 10⁻⁸, this is first prepared into a noisy plus state with η₀ ∼ 0.0152 and δ₀ = 10⁻⁸/3. The total noise, η_n + δ_n, will decease for the first three rounds of parity checking. After the fourth round we have a plus state with a total error of only 2.707 × 10⁻⁸. This illustrates that high-fidelity plus states can be prepared from high-fidelity |φ〉 states in a small number of rounds.

5.2. 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 plus state, |Ψ⁺〉. We measure a two-qutrit stabilizer operator and postselect such that we project onto the subspace spanned by

$$\begin{eqnarray*} |{0_{L}}\rangle & = & ( |{0,0}\rangle + \omega |{1,2}\rangle + \omega^{2} |{2,1}\rangle ) , \\ |{1_{L}}\rangle&=& X_{1}X_{2}|{0_{L}}\rangle = ( |{1,1}\rangle + \omega |{2,0}\rangle + \omega^{2} |{0,2}\rangle ) , \\ \nonumber |{2_{L}}\rangle&=&X_{1}^{2}X_{2}^{2}|{0_{L}}\rangle = ( |{2,2}\rangle + \omega |{0,1}\rangle + \omega^{2} |{1,0}\rangle), \end{eqnarray*}$$

and then decode onto a single qutrit. When successful this produces the following transformation:

$$\begin{equation} |{\Psi^{+}}\rangle^{\otimes 2} \rightarrow ( |{0}\rangle + |{1}\rangle + (\omega+\omega^{2}) |{2}\rangle) /\sqrt{3}. \end{equation} \tag{ 54 }$$

Noting that ω + ω² = −1, we find that the output is a phase state

$$\begin{equation} |{\Phi_{0, \pi}}\rangle = ( |{0}\rangle + |{1}\rangle - |{2}\rangle) /\sqrt{3}. \end{equation} \tag{ 55 }$$

Finally, we have a magic state of the desired form.

5.3. Gate injection

Let us begin by considering a general phase state |Φ_θ,ϕ〉. Given such a magic state, and a second qutrit state in any state, |ψ〉, we perform gate injection by measuring Z₁Z^†₂. Given measurement outcome ω^k, we perform the decoding |x,y〉 → |y + k〉. The ω^k outcome effects a unitary, U_k,θ,ϕ, on |ψ〉, that is diagonal in the computational basis with eigenvalues

$$\begin{eqnarray} U_{0,\theta, \phi} & = & ( 1,\,\mathrm{e}^{\mathrm{i} \theta},\, \mathrm{e}^{\mathrm{i} \phi}) , \nonumber\\ U_{1,\theta, \phi} & = & ( \mathrm{e}^{\mathrm{i} \theta}, 1, \,\mathrm{e}^{\mathrm{i} \theta}) , \nonumber\\ U_{2,\theta, \phi} & = & ( \mathrm{e}^{\mathrm{i} \phi}, \,\mathrm{e}^{\mathrm{i} \theta}, 1 ) . \end{eqnarray} \noindent \tag{ 56 }$$

Each unitary occurs with equal probability because the phase state 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 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 $\mathcal {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 appendix C, where we show that the group $\mathcal {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 currently do not have a complete proof.

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]]₃ 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 $\frac {1}{3}(3-\sqrt {6})\approx 18.3\%$ (see note 4). We have introduced a seven-qutrit code and shown how this code can improve the distillation resource region of [[5,1,3]]₃. Interestingly, in the case of the seven-qubit code, the distillation threshold of the Hadamard eigenstates is tight. And although the seven-qutrit code does not distil the Hadamard eigenstates, the numerical results of section 1 suggest that this tightness is retained in the qutrit case, where all the mixed states along the non-stabilizer segments of the convex line of |H_±〉 are distillable. An important difference with respect to the qubit case is the existence of one Hadamard eigenstate which is not an attractor to the distillation protocols studied, but which forms an unstable fixed point. This state is also the Hadamard eigenstate furthest from the stabilizer region.

After searching the state space for distillable regions, we found another set of distillable states outside the Hadamard plane. These states are eigenstates of a non-degenerate Clifford operator, and it is an open question whether they are distillable by any other stabilizer code or are unique to the five-qutrit code. Under depolarizing noise, these states are distillable up to a noise threshold of 34.4%.

We have also shown how to convert the outputs of the five-qutrit 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 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² 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 realize all or part of the qudit Clifford group?

We thank Matty Hoban and Joe Fitzsimons for useful discussions and Chris Ferrie for his valuable comments on the manuscript. We acknowledge financial support from the EPSRC, the Leverhulme Trust and the EU (QESSENCE).

We have chosen a particular decoding for which the distillation protocol has the simplest behaviour. If instead, the canonical decoding were used, without the addition of an R rotation, then purification would still occur, albeit between each iterate the output would cycle between different states. In the Hadamard plane, we would observe an oscillation between |H_±〉, which is obvious since R|H_±〉 = |H_∓〉. In contrast, for the |φ〉 state there is a more complex four-cycle, illustrated in figure (A.1), such that the distillation map for one iterate, $\mathcal {E}$, performs $\mathcal {E}(|{\varphi _{j}}\rangle )=|{\varphi _{j+1}}\rangle$ and |φ〉 = |φ₁〉 = |φ₅〉. The four cycling states are related by

$$\begin{eqnarray} |{\varphi_{2}}\rangle & = & R^{\dagger} |{\varphi_{1}}\rangle, \nonumber\\ |{\varphi_{3}}\rangle & = & H |{\varphi_{1}}\rangle, \\|{\varphi_{4}}\rangle & = & R^{\dagger}R^{\dagger} |{\varphi_{1}}\rangle.\nonumber \end{eqnarray} \noindent \tag{ A.1 }$$

However, by considering $\mathcal {E'}(\rho )= R \mathcal {E'}(\rho ) R^{\dagger }$, this cycling behaviour vanishes. Note also that this cycling behaviour is seen not only for the pure states but also for depolarized states and so all these states are distilled by the five-qutrit code.

Zoom In Zoom Out Reset image size

Based on the generalized magic state protocol constructed in section 5.1, we will compute the exact distillation map on the Bloch components after a single round of distillation of the five-qutrit code.

We start by expressing the stabilizer generators and the logical operators of the [[5,1,3]]₃ in terms of the σ_j,k operators as shown in table B.1. The input to the distillation protocol is five identical copies of a qutrit state ρ. Using equation (23), the input state is expressed as

$$\begin{equation} \rho^{\otimes 5}=\frac{1}{3^5}\displaystyle\sum_{(\boldsymbol j,\boldsymbol k)\in\mathbb{Z}_{3}^{5}}\alpha_{j_{1}\dots j_{5},k_{1}\dots k_{5}}\sigma_{j_{1}\dots j_{5},k_{1}\dots k_{5}}. \end{equation} \tag{ B.1 }$$

Using equation (24), a successful measurement of the four stabilizer generators with outcome +1 corresponds to the following projector:

$$\begin{equation} \Pi=\frac{1}{3^{4}}\displaystyle\sum_{\boldsymbol q\in\mathbb{Z}_{3}^{4}}g_{1}^{q_{1}}g_{2}^{q_{2}}g_{3}^{q_{3}}g_{4}^{q_{4}}. \end{equation} \tag{ B.2 }$$

Substituting the stabilizer generators in table B.1 into the above expression and using the composition law in equation (3), the projector becomes

$$\begin{eqnarray}&&\fl \Pi=\frac{1}{81}\displaystyle\sum_{\boldsymbol q\in\mathbb{Z}_{3}^{4}}(\sigma_{(q_{1}-q_{3}),(-q_{4})}\otimes\sigma_{(q_{2}-q_{4}),(q_{1})}\otimes \sigma_{(q_{3}),(-q_{1}+q_{2})}\otimes\sigma_{(-q_{1}+q_{4}),(-q_{2}+q_{3})}\otimes\sigma_{(-q_{2}),(-q_{3}+q_{4})}).\nonumber\\ \end{eqnarray} \noindent \tag{ B.3 }$$

We can simplify the notation of the above expression by using equation (4), which gets rid of the tensor product sign, as shown in equation (B.15).

The distillation map in equation (26) can be put into a simpler form as follows:

$$\begin{equation} \alpha_{j,k}^{\mathrm{out}}=\frac{{\mathrm{tr}}(\Pi\rho^{\otimes n}\Pi^{\dagger}(\sigma_{j,k}^{L})^{\dagger})}{{\mathrm{tr}}(\rho^{\otimes n}\Pi)}, \end{equation} \noindent \tag{ B.4 }$$

$$\begin{equation} \hspace{1.7pc}=\frac{{\mathrm{tr}}(\rho^{\otimes n}\Pi\sigma_{-j,-k}^{L})}{{\mathrm{tr}}(\rho^{\otimes n}\Pi)}. \end{equation} \noindent \tag{ B.5 }$$

where in the last step equation (10), Π = Π^†, [σ^L_j,k,Π] = 0 and the cyclic property of the trace were used.

The remaining task is to substitute equations (B.1) and (B.3) into  (B.5) to calculate the distillation map on the Bloch components.

Let us start by evaluating tr(ρ^⊗5Π). Recall that all the σ_j,k operators are traceless except for the identity operator σ_0,0. Therefore, the only terms that will survive in tr(ρ^⊗5Π) are the coefficients of the identity operator. We get the identity operator in ρ^⊗5Π when the σ_j,k operators in ρ^⊗5 and the σ_j',k' in Π have the opposite subscripts (i.e. σ_j,kσ_j',k' = σ_0,0 iff j = j' and k = k'). As a result, tr(ρ^⊗5Π) will be the sum of all the Bloch components that are the coefficients of the σ_j,k operators such that the subscripts (j,k) are the negative of the subscripts in equation (B.3). In fact, since the summation is over all the elements of the ring $\mathbb {Z}^{4}_{3}$ it is possible to multiply all the subscripts by (−1) without changing the actual value of the summation. Hence, tr(ρ^⊗5Π) can be compactly expressed as shown in equation (B.16). In a similar way, we can express tr(ρΠσ^L_−j,−k) for all four logical operators. For example, in the case of evaluating the output Bloch component α^out_1,0, equation (B.5) becomes

$$\begin{equation} \alpha_{1,0}^{\small{\mathrm{out}}}=\frac{{\mathrm{tr}}(\rho^{\otimes 5}\Pi\sigma_{-1,0}^{L})}{{\mathrm{tr}}(\rho^{\otimes 5}\Pi)}, \end{equation} \noindent \tag{ B.6 }$$

with tr(ρ^⊗5Πσ^L_−1,0) given in equation (B.17).

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 equation  (B.18), where the subscript r represents the number of 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 be compactly expressed in terms of a single function. This function is given in equation (B.19), and based on this function, the distillation map can be expressed as

$$\begin{equation} A_{r+1}=\mathcal{F}(A,B,C,D), \end{equation} \tag{ B.7 }$$

$$\begin{equation} B_{r+1}=\mathcal{F}(B^*, A, D, C^*), \end{equation} \tag{ B.8 }$$

$$\begin{equation} C_{r+1}=\mathcal{F}(A^*, C, B, D), \end{equation} \tag{ B.9 }$$

$$\begin{equation} D_{r+1}=\mathcal{F}(B^*, D^*, A^*, C). \end{equation} \tag{ B.10 }$$

These four expressions represent 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 section 4), and will result in the cycling behaviour in appendix A. We need to ensure that the map in equation (28) is applied after every iteration. This can be easily achieved in our formalism by appropriate relabelling as follows:

$$\begin{equation} A_{r+1}=\mathcal{F}(D^{*},C,B^{*},A), \end{equation} \tag{ B.11 }$$

$$\begin{equation} B_{r+1}=\mathcal{F}(C^*, D^{*}, A, B), \end{equation} \tag{ B.12 }$$

$$\begin{equation} C_{r+1}=\mathcal{F}(D, B^*, C, A), \end{equation} \tag{ B.13 }$$

$$\begin{equation} D_{r+1}=\mathcal{F}(C^*, A^*, D, B^*), \end{equation} \tag{ B.14 }$$

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 multi-variable polynomials of the order of five. It is known from the famous Abel–Ruffini theorem that there is no algebraic solution for a general polynomial of the order of five or above. Therefore, the best way to discover the fixed points of the distillation is by 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 towards a fixed point.

$$\begin{equation} \fl \Pi=\frac{1}{81}\displaystyle\sum_{\boldsymbol{q}\in\mathbb{Z}_3^4}\sigma_{(q_{1}-q_{3})_{1}(q_{2}-q_{4})_{2}(q_{3})_{3}(-q_{1}+q_{4})_{4}(-q_{2})_{5},(-q_{4})_{1}(q_{1})_{2}(-q_{1}+q_{2})_{3}(-q_{2}+q_{3})_{4}(-q_{3}+q_{4})_{5}} , \end{equation} \tag{ B.15 }$$

$$\begin{equation} \fl {\mathrm{tr}}(\rho^{\otimes 5}\Pi) = \frac{1}{81}\displaystyle\sum_{\boldsymbol{q}\in\mathbb{Z}_3^4}\alpha_{(q_{1}-q_{3})_{1}(q_{2}-q_{4})_{2}(q_{3})_{3}(-q_{1}+q_{4})_{4}(-q_{2})_{5},(-q_{4})_{1}(q_{1})_{2}(-q_{1}+q_{2})_{3}(-q_{2}+q_{3})_{4}(-q_{3}+q_{4})_{5}} , \end{equation} \noindent \tag{ B.16 }$$

$$\begin{eqnarray} &&\fl {\mathrm{tr}}(\rho^{\otimes 5}\Pi\sigma_{-1,0}^{L})\!=\!\frac{1}{81}\displaystyle\!\sum_{\boldsymbol{q}\in\mathbb{Z}_3^4}\alpha_{(q_{1}-q_{3})_{1}(q_{2}-q_{4})_{2}(q_{3})_{3}(-q_{1}+q_{4})_{4}(-q_{2})_{5},(-q_{4}+1)_{1}(q_{1}+1)_{2}(-q_{1}+q_{2}+1)_{3}(-q_{2}+q_{3}+1)_{4}(-q_{3}+q_{4}+1)_{5}},\nonumber\\ \end{eqnarray} \noindent \tag{ B.17 }$$

$$\begin{eqnarray} &&\fl {\mathrm{tr}}(\rho^{\otimes 5}\Pi)\!=\!\frac{1}{81}(1+ 10 (\left|A_r\right|{}^2+\left|D_r\right|{}^2) (\left|B_r\right|{}^2+\left|C_r\right|{}^2)+5(B_r^2 A_r^* C_r^{*2}+ D_r^2 A_r^{*2} B_r^* \nonumber\\ &&+ D_r (A_r^2 D_r C_r^*+B_r^2 C_r^2)+ B_r^{*2} (A_r C_r^2+C_r^{*2} D_r^{*})+D_r^{*2} (A_r^2 B_r+C_r A_r^{*2}) )), \end{eqnarray} \noindent \tag{ B.18 }$$

$$\begin{eqnarray} &&\fl \mathcal{F}(A,B,C,D)\!=\!\frac{1}{81}(B_r^5 +10 B_r (D_r A_r^*+B_r^*) (A_r C_r^*+ C_r D_r^*)+ 5 (A_r C_r^2 \left|A_r\right|^2+D_r^2 (A_r B_{\mathrm{r}}^{*2}+C_r)\nonumber\\ &&+A_r^{*2} (B_r^{*2}D_r^*+C_r^*)+D_r^{*2} (A_r^2+D_r C_r^{*2})+|C_r|^4 B_r^*))/{\mathrm{tr}}(\rho^{\otimes 5}\Pi). \end{eqnarray} \tag{ B.19 }$$

Here we show that $\mathcal {C}_{+N}$ contains an infinite number of unitaries. We do so by considering a particular unitary,

$$\begin{equation} K = H.N.H.N.H , \end{equation} \tag{ C.1 }$$

which is an alternating sequence of Hadamards and our non-Clifford gate. The eigenvalues e^iπλ_n of K can be analytically calculated by using Mathematica, such that

$$\begin{equation} \lambda_{1} = \mathrm{arctan} ( \sqrt{2}) / \pi , \end{equation} \tag{ C.2 }$$

$$\begin{equation} \lambda_{2} = -\mathrm{arctan} ( \sqrt{2}) / \pi , \end{equation} \tag{ C.3 }$$

$$\begin{equation} \lambda_{3} = 1 / 2. \end{equation} \tag{ C.4 }$$

Since both λ₁ and λ₂ are irrational numbers, every power of K gives a different unitary, and so the single-qutrit group, $\mathcal {C}_{+N}$, is infinite in order. The irrationality of λ₁, and hence also λ₂, can be shown by expressing it as $\lambda _{1}= \mathrm {arccos}(1/\sqrt{3})/\pi$ and then using a theorem from [37].