Discrete Wigner Function Derivation of the Aaronson-Gottesman Tableau Algorithm

The Gottesman-Knill theorem established that stabilizer states and operations can be efficiently simulated classically. For qudits with dimension three and greater, stabilizer states and Clifford operations have been found to correspond to positive discrete Wigner functions and dynamics. We present a discrete Wigner function-based simulation algorithm for odd-$d$ qudits that has the same time and space complexity as the Aaronson-Gottesman algorithm. We show that the efficiency of both algorithms is due to the harmonic evolution in the symplectic structure of discrete phase space. The differences between the Wigner function algorithm and Aaronson-Gottesman are likely due only to the fact that the Weyl-Heisenberg group is not in $SU(d)$ for $d=2$ and that qubits have state-independent contextuality. This may provide a guide for extending the discrete Wigner function approach to qubits.


I. INTRODUCTION
A general set of unitary quantum operators on n-qudit states requires a number of operators that grows exponentially with n. An important exception to this rule involves the set of Clifford operators acting on stabilizer states. These states play an important role in quantum error correction [1] and are closed under action by Clifford gates. Efficient simulation of such systems was demonstrated with the tableau algorithm of Aaronson and Gottesman [1,2] for qubits (d = 2). Finding the underlying reason for why such an efficient algorithm is possible for Clifford circuit simulation has since been the subject of much study [3][4][5].
Recent progress has been the result of work by Wootters [6], Eisert [4], Gross [7], and Emerson [5], who have formulated a new perspective based on the discrete phase spaces of states and operators in finite Hilbert spaces using discrete Wigner functions. They have shown that stabilizer states have positive-definite discrete Wigner functions and that Clifford operators are positive-definite maps. This implies that Clifford circuits are efficiently simulatable on classical computers. In odd-dimensional systems, stabilizer states have been shown to be the discrete analogue to Gaussian states in continuous systems [7] and Clifford group gates have been shown to have underlying harmonic Hamiltonians that preserve the discrete Weyl phase space points [8]. This means Clifford circuits are expressible by path integrals truncated at order 0 and are thus manifestly classical [8,9].
This poses the question of what the relationship is between past efficient algorithms for Clifford circuits and the propagation of discrete Wigner functions of stabilizer states under Clifford operators. Here, we show that the original Aaronson-Gottesman tableau algorithm is actually equivalent to such a discrete Wigner function propagation and that the tableau matrix coincides with the discrete Wigner function of a stabilizer state. We accomplish this by first developing a Wigner function-based algorithm that classically simulates stabilizer state evolution under Clifford gates and measurements in theẐ Pauli basis for odd d. We then show its equivalence to the well-known Aaronson-Gottesman tableau algorithm [2]. Both algorithms require O(n 2 ) qudits to represent n stabilizer states, O(n) operations per Clifford operator, and deterministic and random measurements require O(n 2 ) operations.
The Aaronson-Gottesman tableau algorithm makes use of the Heisenberg representation. This means that evolution is accomplished by updating an associated tableau or matrix representation of the Clifford operators instead of the stabilizer states themselves. Our algorithm is framed in the Schrödinger picture, and involves evolving the Wigner function of stabilizer states. Since the two algorithms are equivalent, the formulation of Clifford simulation in the Heisenberg picture is a choice and not a necessity for its efficient simulation. Furthermore, we reveal the purely classical basis of both algorithms and the physically intuitive phase space structures and symplectic properties on which they rely.

II. DISCRETE WIGNER FUNCTION OF THE STABILIZER FORMALISM
Before we discuss the discrete Wigner function, we introduce a basic framework that defines how a phase space behaves for odd d-dimensional Hilbert spaces. To begin, we associate the computational basis with the position basis, such that the PauliẐ j operator on the jth qubit for n qubits acts as a "boost" operator: The discrete Fourier transform operator is defined by: d kj lj |k 1 , . . . , k j , . . . , k n l 1 , . . . , l j , . . . , l n | , This is equivalent to the Hadamard gate and allows us to define the PauliX j operator as follows: WhileẐ j is a boost,X j is a shift operator becausê X δq j |k 1 , . . . , k j , . . . , k n ≡ |k 1 , . . . , k j ⊕ δq, . . . , k n , (4) where ⊕ denotes integer addition mod d.
We can reexpress the boostẐ j and shiftX j operators in terms of their generators, which are the conjugateq j andp j operators respectively: Thus, we can refer to theX j basis as the momentum (p j ) basis, which is equivalent to the Fourier transform of the q j basis:p These bases form the discrete Weyl phase space (p, q). The Wigner function Ψ x (p, q) of a pure state |Ψ is defined on this discrete Weyl phase space: We will shortly be interested in the discrete Wigner function of stabilizer states. But first we introduce the effect that the Clifford gates have in this discrete Weyl phase space.

A. Clifford Gates
A Clifford group gateV is related to a symplectic transformation on the discrete Weyl phase space, governed by a symplectic matrix MV and vector αV [8]: Wigner functions Ψ x (x) of states evolve under Clifford operatorsV by where x ≡ (p, q). When considering Clifford gate propagation, we can restrict to the generators of the group. One such set of generators is made up of the phase-shift gateP i , the Hadamard gateF i , and the controlled-notĈ ij (which act on the ith and jth qudits).
The phase shiftP i is a one-qudit gate with the underlying Hamiltonian HP i = − d+1 2 q 2 i + d+1 2 q i [8]. Without loss of generality, we will instead consider which we will refer to as the phase-shift gate in this paper. We note that the usual phase-shift can be obtained from the new one within the Clifford group: where [P i ,Ẑ i ] = [P i ,Ẑ i ] = 0. Hence,P i is an adequate replacement generator forP i , and we will use it instead ofP i from now on. Since its Hamiltonian has no linear term (HP i = −q 2 i ), this leads to an easier presentation ahead since αP The corresponding equations of motion forP i areṗ i = 2q i andq i = 0. Hence, for ∆t = 1, The Hadamard gateF i is a one-qudit gate and has the underlying Hamiltonian HF i = − π 4 (p 2 i + q 2 i ) [8]. The corresponding equations of motion areṗ i = π 2 q i andq i = − π 2 p i . Hence, for ∆t = 1, +δ i,j δ n+i,k − δ n+i,j δ i,k , and αF i = 0. Finally, the two-qudit CNOTĈ ij on control qudit i and second qudit j has the corresponding Hamiltonian HĈ ij = p i q j [8]. The corresponding equations of motion are (ṗ i ,ṗ j ) = −(0, p i ) and (q i ,q j ) = (q j , 0). Hence, for ∆t = 1, and αĈ ij = 0.

B. Wigner Functions of Stabilizer States
A discrete Wigner function associated with the boost and shift operators defined in Eqs. 5 and 6 is given by the following theorem [8]: Theorem 1. The discrete Wigner function Ψ x (x) of a stabilizer state Ψ for any odd d and n qudits is δ Φ×x,r for 2n × 2n matrix Φ and 2n vector r with entries in Z/dZ.

III. WIGNER STABILIZER ALGORITHM
With the discrete Wigner function of a stabilizer state defined in Theorem 1 and the effect of the Clifford group generators on discrete Wigner functions defined in Eq. 10, we can now examine the effect Clifford operators have on stabilizer states.

A. Stabilizer Representation
From Theorem 1, propagation of the stabilizer state Ψ can be represented by considering the state in the Wigner representation Ψ x (x) = δ Φt·x,rt . In this way, Φ t and r t specify a linear system of equations in terms of p t and q t . The first n rows of Φ t are the coefficients of (p t , q t ) T in p 0 (p t , q t ) and the last n rows of Φ t are the coefficients of (p t , q t ) T in q 0 (p t , q t ): The Kronecker delta function sets this linear system of equations equal to r t . In this way, an affine map-a linear transformation displaced from the origin by r t -is defined. This system of equations must be updated after every unitary propagation and measurement. Since the Wigner functions Ψ x (x) of stabilizer states propagate under M as Ψ x (Mx), it follows that The importance of vector r t and when it must be updated will become evident when we consider random measurements. Hence, after n operations M 1 , M 2 , . . ., M n , The matrices are ordered chronologically left-to-right instead of right-to-left.
and I n is the n × n identity matrix. Thus the the stability matrices M forF i ,P i andĈ ij given in Eqs. 13-15 differ from their inverses only by sign changes in their off-diagonal elements: −δ i,j δ n+i,k + δ n+i,j δ i,k , and M −1 We assume the quantum state is initialized in the computational basis state Ψ 0 = |0 ⊗ · · · ⊗ |0 n and so initially we should set Φ 0 = 0 0 0 I n and r 0 = 0. The initial stabilizer state is Ψ x0 = δ q t ,0 . However, it will become clear when we discuss measurements that it is practically useful to instead set thereby setting Ψ x0 = δ (p t ,q t ),(0,0) -not a true Wigner function. This is equivalent to the last matrix if the first n rows in Φ t x and r t are ignored-the same as ignoring p 0 (p t , q t ). In fact, we have two Wigner functions here: one defined by the first n rows and another by the last n rows. We proceed in this manner, ignoring the first n rows, until their usefulness becomes apparent to us. For n qudits unitary propagation requires O(n 2 ) dits of storage to track Φ t and r t . More precisely, since Φ t is a 2n × 2n matrix and r t is an 2n-vector, 2n(2n + 1) dits of storage are necessary.

B. Unitary Propagation
Φ t contains the coefficients of the linear equations relating x 0 to x t . Each row is one equation relating q 0i or p 0i to x t . When manipulating rows of Φ t we shall refer to the linear equations that these rows define.
Examining Eqs. 20, 21, and 22, we see that the inverse stability matrices of the generator gatesF i ,P i and C ij are the sum of an identity matrix and a matrix with a finite number of non-zero off-diagonal elements. The number of these off-diagonal elements is independent of the number of qudits, n. Hence, multiplying Φ t with a new stability matrix in Eq. 17 and evaluating the matrix multiplication is equivalent to performing a finite number of n-vector dot products and so requires O(n) operations. Therefore, keeping track of propagation of stabilizer states by Clifford gates can be simulated with O(n) operations.
Let us examine these unitary operations more closely. Defining ⊕ and to be mod d addition and subtraction respectively, we find: Phase gate on qudit i (P i ). For all j ∈ {1, . . . , 2n}, set Φ j,n+i → Φ j,n+i 2Φ j,n .
Hadamard gate on qudit i (F i ). Negate Φ j,i mod d, and then swap Φ j,i and Φ n+i,j .
This confirms that unitary propagation in this scheme requires O(n) operations.

C. Measurement
The outcome of a measurementẐ i on a stabilizer state can be either random or deterministic. As described above, the bottom half of Φ t defines q 0j for j ∈ {1, . . . , n}, each of which is a linear combination of q ti and p ti . The entries in the (n + j)th row of Φ t give the coefficient of p ti and q ti in q 0j for j ∈ {1, . . . , n}. If the coefficient of p ti in any q 0j is non-zero then the mea-surementẐ i will be random. If all coefficients of p ti are zero for q 0j ∀j, then the measurement ofẐ i will be deterministic. This can be seen from the fact that if our stabilizer state |Ψ is an eigenstate ofẐ i thenẐ i |Ψ = e iφ |Ψ for some φ ∈ R and (discrete) Wigner functions do not change under a global phase. Thus, measuringẐ i leaves the Wigner function of |Ψ invariant if the measurement is deterministic. SinceẐ i is a boost operator that increments the momentum of a state by one, its effect on the linear system of equations specified by the Wigner function is: (24) Thus, if the lower half of the ith column of Φ t is zero thenẐ i leaves the Wigner function invariant (and so the measurement is deterministic). Verifying that these coefficients are all zero takes O(n) operations for eachẐ i .
In other words, to see if a given measurement ofẐ i is random or deterministic, a search must be performed for non-zero Φ tn+j,i elements. If such a non-zero element exists then the measurement is random since it means that the final momentum of qudit i affects the state of the stabilizer and so its position must be undetermined (by Heisenberg's uncertainty principle). If no such finite Φ tn+j,i element exists, then the measurementẐ i is deterministic.
We now describe the algorithm in detail for these two cases: Case 1: Random Measurement Let the (n + j)th row in the bottom half of Φ t have a non-zero entry in the ith column, Φ tn+j,i . Since the random measurementẐ i will project qudit i onto a position state, we will replace the (n + j)th row with q 0i = q i (the uniformly random outcome of this measurement). After this projection onto a position state, none of the other qudit's positions should depend on qudit i's momentum, p ti . To accomplish this, before we replace row (n+j), we solve its equation for p ti and substitute every instance of p ti in the linear system of equations with this solution. As a result, every equation will no longer depend on p ti and we can go ahead and replace the (n + j)th row with q 0i = q i .
There is one more thing to do, which will be important for deterministic measurements: replace the jth row with the old (n + j)th row. This sets p 0i = q 0j (p t , q t ), which becomes the only remaining equation explicitly dependent on p ti . In other words, p 0i ∝ p ti , similar to the beginning when we set p 0i = p ti by setting Φ = I 2n .
However, now we also conserve any dependence p 0i has on the other qudits incurred during unitary propagation. In other words, we conserve p ti 's dependence upon the other qudits, but only in the Wigner function specified by the top n rows, which we ignore otherwise.
After replacing the equation specified by row (n + j) of Φ t and r t with a randomly chosen measurement outcome q i (i.e. q 0i = q i ), the identification of rows (n + i) and (n + j) are exchanged, so that the former now specifies q 0j (p t , q t ) while the latter specifies q 0i (p t , q t ), respectively. p 0i has also been updated by replacing the jth row in the first half of Φ t , with the (n + j)th row we just changed. Again, this row now describes p 0i (p t , q t ) while the ith row now specifies p 0j (p t , q t ). Overall, this takes O(n 2 ) operations since we are replacing O(n) rows with O(n) entries.
Case 2: Deterministic Measurement Since the measurement is deterministic, Φ t and r t do not change. The n equations specified by the bottom half of Φ t x t = r t can be used to solve for q ti -the deterministic measurement outcome. In general, this can also be done by inverting Ψ t and evaluating x t = Φ −1 t · r t for q i . Aaronson and Gottesman themselves noted that such a matrix inversion is possible, but practically takes O(n 3 ) operations [10].
Fortunately, there is another method that scales as O(n 2 ) and requires use of the n equations represented by the top n rows of Φ t , which were included in our description by setting Φ 0 = I 2n . The linear system of n equations represented by Φ t x t = r t can be written as where we are interested in linear combinations of the bottom half, q 0 (p t , q t ), to solve for the measurement outcome q ti : n j=1 c ij q 0j = q ti .
where c ij ∈ Z/dZ. Lemma 1. The coefficient in front of p ti in the row of Φ t that specifies p 0j (p t , q t ), Φ tji , is equal to the coefficient c ij in front of q 0j that makes up q ti in Eq. 27. Or, equivalently, Proof. Under evolution under the Clifford group operators, This means that we can express the matrix inversion as follows: Therefore M −1 t 11i,j = (M t ) 22i,j , and so This property can also be seen in the drawing of phase space shown in Fig. 1. There, initial perpendicular p 0j and q 0j manifolds are drawn along with harmonically evolved p ti and q ti manifolds, which remain perpendicular to each other and make an angle α to the first p 0j and q 0j manifolds, respectively. The projection of q ti (p 0 , q 0 ) onto q 0j can be represented as the length b of a right triangle's adjacent side to the angle α, with an opposite side set to some length a. The projection of p 0j (p t , q t ) onto p ti is similarly represented by the length b of a right triangle's adjacent side to the angle α, with an opposite side also set to length a. It follows that the third angle β in both triangles must be the same, and so by the law of sines Therefore b = b and so these two projections are equal to one another.
In the discrete Weyl phase space such manifolds must lie along grid phase points and obey the periodicity in x p and x q , but the premise is the same.
Overall, the procedure outlined in Lemma 1 for deterministic measurements takes O(n 2 ) operations since Eq. 28 is a sum of O(n) vectors made up of O(n) components. So, the overall measurement protocol takes O(n 2 ) operations. Note that this formulation of the algorithm shows that it is the symplectic structure on phase space and the linear transformation under harmonic evolution that allows the inversion (Eq. 33) to be performed efficiently.

IV. AARONSON-GOTTESMAN TABLEAU ALGORITHM
The Aaronson-Gottesman tableau algorithm was originally defined for qubits (d = 2) [2]. Like the algorithm we presented in the previous section, it only requires overall O(n 2 ) operations for propagationan and measurement for n qubits. The algorithm has been proven to be extendable to d > 2 [11] and similar algorithms have been formulated in d > 2 [12]. Alternatives have also been developed to the tableau formalism, though they prove to be equally efficient in worst-case scenarios [13]. However, we are not aware of any direct extension of the Aaronson-Gottesman tableau algorithm to dimensions greater than two. In this and the next section, we will show that the Wigner algorithm presented in Section III is equivalent to the Aaronson-Gottesman tableau algorithm.

A. Stabilizer Representation
The Aaronson-Gottesman algorithm is defined the stabilizer formalism. It keeps track of the evolution of a stabilizer state by updating the generators of the stabilizer group, elements of which are defined as follows: Definition 1. A set of operators that satisfies S = {ĝ ∈ P such thatĝ |ψ = |ψ } are called the stabilizers of state |ψ , where P is the set of Pauli operators, each of which has the form e πi 2 αP 1 · · ·P n where α ∈ {0, 1, 2, 3} for n qubits withP i ∈ {Î i ,Ẑ i ,X i ,Ŷ i }.
For the sake of completeness, we present here a summary of the Aaronson-Gottesman algorithm, in order to compare it to our algorithm. For more details, see [1,2].
Each n-qubit stabilizer state is uniquely determined by 2 n Pauli operators. There are only n generators of this Abelian group of 2 n operators. Therefore, an nqubit stabilizer state is defined by the n generators of its stabilizer state. Every element in this set of generators, {ĝ 1 ,ĝ 2 , . . . ,ĝ n }, is in the Pauli group, and each generator has the form:ĝ Any unitary propagation by Clifford operators or measurement of the stabilizer state changes at least some of theP ij elements of the n generators of the state's stabilizer. This includes the ±1 phase in Eq. 36, which must also be kept track of in Aaronson-Gottesman's algorithm.

B. Unitary Propagation
For each Clifford operation, Aaronson and Gottesman showed that only O(n) operations are necessary to update all generators [2]. Specifically, according to the update rules in Table I, each generator can be updated with a constant number of operators for every single Clifford gate, therefore O(n) in total. However, it is a little more complicated to update the generators after each measurement. To do this efficiently, Aaronson introduced "destabilizers": Definition 2. Destabilizers {ĝ 1 , . . . ,ĝ n } are the operators that generate the full Pauli group with the stabilizers {ĝ 1 , . . . ,ĝ n }. They have the following properties: (i)ĝ 1 ,ĝ 2 , . . .,ĝ n commute.
(ii) Each destabilizerĝ h anti-commutes with the corresponding stabilizerĝ h , and commutes with all other stabilizers.
To incorporate the destabilizers, a tableau becomes useful to see how they play a role in updating the stabilizer generators during measurement [2].
Aaronson-Gottesman defined such a 2n × (2n + 1) binary tableau matrix as:  This matrix contains 2n rows. The first n rows denote the destabilizersĝ 1 toĝ n while rows (n + 1) to 2n represent the stabilizersĝ 1 toĝ n . The (n + 1)th bit in each row denotes the phase (−1) ri for each generator. We encode the jth Pauli operator in the ith row as shown in Table II. We can update the stabilizers and destabilizers as follows: Hadamard gate on qubit i For all j ∈ {1, 2, ..., 2n}, r j → r j ⊕ x ji z ji , then swap x ji with z ji .
CNOT gate on control qubit i and target qubit j For all k ∈ {1, 2, ..., 2n}, r k → r k ⊕ x ki z kj (x kj ⊕ z ki ⊕ 1), These actions correspond to those given in Table I. Notice the striking similarity of these tableau transformation rules under unitary propagation to the Φ transformation rules in Section IV. The most glaring difference is that the Aaronson-Gottesman algorithm involves updates of the vector r. We will discuss this and its connection to the dimension d = 2 of the system in Section V. It is clear that these transformations also take O(n) operations each.

C. Measurement
To describe the measurement part of the algorithm, we need to first define a rowsum operation in the tableau which corresponds to multiplying two Pauli operators together.
As defined in [2]: Rowsum: To sum row i and j, first update the bits that represent operators by x ik ⊕ x jk and z ik ⊕ z jk for k = 1, . . . , n. To calculate the resultant phase, Aaronson and Gottesman first defined the following function: (37) Since each stabilizer generator is the tensor product of n single qubit Pauli operators (see Eq. 36), they must be multiplied together to obtain the phase: (38) Having defined the rowsum function, let us now consider a measurement ofẐ i on qubit i. For d = 2, Pauli group operators can only commute or anti-commute with each other. IfẐ i anti-commutes with one or more of the generators then the measurement is random. IfẐ i commutes with all of the generators then the measurement is deterministic. We consider these two cases: Case 1: Random Measurement Z i anti-commutes with one or more of the generators. If there is more than one, we can always pick a single anti-commuting generator,ĝ j , and update the rest by replacing them with their product withĝ j (i.e. taking the rowsum of their corresponding rows) such that they commute withẐ i . These updates take O(n 2 ) operations. Finally, we only need to replaceĝ j byẐ i .
In other words, with respect to the tableau, there should exist at least one j ∈ {n + 1, n + 2, ..., 2n} such that x ji = 1. Replacing all rows where x ki = 1 for k = j with the sum of the jth and kth row (using the rowsum function) sets all x ki = 0 for k = j.
Finally, we replace the (j − n)th row with the jth row and update the jth row by setting z ji = 1 and all other x jk s and z jk s to 0 for all k. We output r j = 0 or r j = 1 with equal probability for the measurement result. This procedure takes O(n 2 ) operations because each rowsum operation takes O(n) operations and up to n−1 rowsums may be necessary.
Case 2: Deterministic Measurement Z i commutes with all generators. In this case, there is no j ∈ {n+1, n+2, ..., 2n} such that x ji = 1 and we don't need to update any of the generators. But we do need to do some work to retrieve the measurement outcome.
MeasurementẐ i commutes with all of the stabilizers, therefore either +Ẑ i or −Ẑ i is a stabilizer of the state. So it must be generated by the generators. The sign ±1 is the measurement outcome we are looking for. This means that where c j = 1 or 0. For those destabilizers g k that satisfy c k = 1. Otherwise c k = 0. This can be seen from where we used part (ii) of Definition 2 of the destabilizers and Eq. 40. The last equality requires c k = 1. Therefore, to find the deterministic measurement outcome, the stabilizers whose corresponding destabilizer anti-commutes with the measurement operationẐ i must be multiplied together. Every row (n + j) in the bottom half of the tableau, such that x ji = 1 (for j ∈ {1, . . . , n}), can be added up together and stored in a temporary register. The resultant phase ±1 of this sum is the measurement result we are looking for.
Checking if each destabilizer commutes or anticommutes withẐ i takes a constant number of operations. One multiplication takes O(n) operations, and there are O(n) multiplications needed. Therefore, a measurement takes O(n 2 ) operations overall.

V. DISCUSSION
As we made clear throughout Section IV, the scaling of the number of required operations with respect to number of qubits n is exactly the same in the Aaronson-Gottesman algorithm in the Wigner algorithm presented in Section III. The two algorithms also require the same number of dits of temporary storage for performing the deterministic measurement. Moreover, there is a correspondence between the tableau employed by Aaronson-Gottesman and the matrix Φ t and vector r t we use. In particular, the tableau is equal to Φ t r t : and This can be seen through the following equation: wherex ≡ (p,q). Multiplying the right hand side of the first equation and the second equation by exp − 2πi d r ti , it follows that (45) In other words, r ti specifies the phase exp − 2πi d r ti of the ith stabilizer, which is itself specified by Φ tn+i,j for j ∈ {0, . . . , 2n}. These are the same roles for r and the tableau in the Aaronson-Gottesman tableau algorithm [2].
Indeed, both algorithms check the bottom half of their matrices for finite elements of Φ n+j,i to determine if a measurement on the ith qudit will be random or not. They also use a very similar protocal to determine the outcome of deterministic measurements. The Wignerbased algorithm motivates these manipulations in terms of the symplectic structure of Weyl phase space and the relationship between the two Wigner functions specified by the top and bottom of Φ, providing a strong physical intuition for their effects. Aaronson and Gottesman motivate these manipulations using the anti-commutation relations between the stabilizer and destabilizer generators. Also, the latter half of both the Wigner function's r t and Aaronson-Gottesman's r are used to determine measurement outcomes. The only fundamental algorithmic difference between the approaches is that the Wigner-based algorithm does not require updates of r t during unitary propagation. The reason for this lies in the fact that Aaronson-Gottesman's algorithm deals with systems with d = 2 while the Wigner-based algorithm is restricted to odd d.
In particular, for the one-qubit Clifford group gate operatorÂ = {P i ,F i } ∀i = {1, . . . , n}, the Aaronson-Gottesman algorithm specifies that for a q-or p-state, its Wigner function evolves by: But for |r = 1 √ 2 (|0 ± i |1 ), a Y -state which is diagonal to q and p, its Wigner function must first be translated: where the translation β can be (1, 0) or (0, 1) equivalently. There is a similar state-dependence for the twoqudit CNOT gateĈ ij . This demonstrates that the Aaronson-Gottesman algorithm is state-dependent on the qubit stabilizer state it is acting on. On the other hand, our algorithm on qudit stabilizer states is state-independent. This likely is a consequence of the fact that the Weyl-Heisenberg group, which is made up of the boost and shift operators defined in Eqs. 5 and 6 that underlie the discrete Wigner formulation, are a subgroup of U (d) instead of SU (d) for d = 2 [14]. Furthermore, qubits exhibit stateindependent contextuality while odd d qudits do not [15].

VI. EXAMPLE OF STABILIZER EVOLUTION
FIG. 2. A decomposition of the two qutrit Wigner function into nine 3 × 3 girds, where each 3 × 3 grid denotes the value of the Wigner function at all pt 1 and qt 1 for a fixed value of pt 2 and qt 2 denoted by the external axes. This organization is used in Fig. 3 below.
As a demonstration of what stabilizer state propagation looks like in the Wigner formalism, we proceed to go through an example of Bell state preparation and measurement starting from the state |0 ⊗ |0 . The prepared Wigner function is illustrated in Fig. 3 with the color black, and the Wigner function represented by setting  (3) Qudit 1 is then measured producing the random outcome 1, which collapses qudit 2 into the same state, so that d) |1 ⊗ |1 results. The black color indicates the Wigner function specified by the lowest n rows of δ Φ t x,r t , and the gray color indicates the Wigner function specified by the highest n rows (q 0 (p t , q t ) and p 0 (p t , q t ), respectively). The evolution and algorithmic implementation are explained in the text.
denoting an initially prepared state of |0 ⊗ |0 . This is clear in Figure 3a by the black band that lies along all Weyl phase space points with q t1 = 0 and q t2 = 0. On the other hand, the gray manifold is perpendicular to the black one, and lies along Weyl phase space points with p t1 = 0 and p t2 = 0.
Acting on this state withF 1 produces 3 2×0 |2 ⊗ |0 . Applying the algorithm specified at the end of Section III B, we find: Thus, the momentum of qudit 1 is now determined and is 0 while the second qudit is unchanged. This can be seen in Fig. 3b, where the q t2 values of the non-zero Weyl phase space points are the same, while the state has rotated by −π/2 in (p t1 , q t1 )-space. A similar transforma-tion has occured for the perpendicular gray manifold.
Acting next withĈ 12 produces the Bell state 1 √ 3 (|00 + |11 + |22 ), which is represented by the following Wigner function: The entanglement between the two qudits is evident in both of their dependence on each other's momenta and positions, p t1 = −p t2 and q t1 = q t2 , specified by the last two rows. Fig. 3c shows that the state is still representable as lines in Weyl phase space, except they now traverse through the different planes of (q t1 , p t1 ) associated with each value of (q t2 , p t2 ). However, if you consider the left column in Fig. 3c corresponding to q t2 = 0, you can see that the only black Weyl phase points are at q t1 = 0. Similarly, the middle column corresponding to q t2 = 1 shows that q t1 = 1, and the right column corresponding to q t2 = 2 shows that q t1 = 2 too, confirming that |P hi = 1 √ 3 (|00 + |11 + |22 ). Thus, the entanglement of the two qudits positions is clearly evident in this figure of the Wigner function.
We then proceed to measure qudit 1. Since the lower two equations involve p t1 , we know that this is a random measurement. Let us pick the outcome to be 1 and set the third row as such, replacing the first row with the old third row. This collapses qudit 2 into the same state: The lower two rows show that now q t1 = 1, as we chose, and q t2 = q t1 = 1. The collapse of qudit 2 into |1 can also been seen in Fig. 3c by the fact that q t1 = 1 only in the 3 × 3 grids that correspond to q t2 = 1 too.
Finally, the fact that the measurement of q t2 is deterministic can be seen in that p t2 is not present in the last two rows of Φ t . Furthermore, it is clear since the first row has a coefficient of 1 in front of p t1 , that the corresponding third row must be added with weight 1 to the fourth row to obtain this deterministic measurement outcome of q t2 = 1. This can also be seen in Fig. 3 by finding the projection of p 01 onto p t2 , which are shown by the gray manifolds in panels a) and d) respectively. They are collinear and so the projection is equal to 1. (Perpendicular manifolds corresponds to a projection of 0, and those that lie π/4 diagonally with respect to each other have a projection equal to 2 in this discrete geometry.)

VII. CONCLUSION
In summary, we introduced an algorithm that efficiently simulates stabilizer state evolution under Clifford gates and measurements in theẐ Pauli basis for odd d qudits. We accomplished this by relying on the phase-space perspective of stabilizer states as discrete Gaussians and Clifford operators as having underlying harmonic Hamiltonians. We showed the equivalence of our algorithm, through Eqs. 44 and 45, to the wellknown Aaronson-Gottesman tableau algorithm [2], revealing that Aaronson-Gottesman's tableau corresponds to a discrete Wigner function. As a consequence, we revealed the physically intuitive phase space perspective of Aaronson-Gottesman's algorithm, as well as its extension to higher odd d.
This work illustrates that no efficiency advantage is gained by using the Heisenberg representation for stabi-lizer propagation. Eq. 45 indicates that the Heisenberg representation is equivalent to the Schrödinger representation in this context; evolving the operators is just as efficient as evolving the states, as perhaps expected.
Lastly, the correspondence between the Wigner-based algorithm and the Aaronson-Gottesman tableau algorithm may point the direction on how to resolve the long-standing issue of describing the Wigner-Weyl-Moyal and center-chord formalism for d = 2 systems. We have shown that the Aaronson-Gottesman algorithm is essentially a d = 2 treatment of the Wigner approach. The salient ingredient appears to be the state-dependence of this evolution, and likely is related to the fact that the Weyl-Heisenberg group is not in SU (d) for d = 2 and that qubits exhibit state-independent contextuality, which odd d qudits do not. Exploring the details of this state-dependence is a promising subject of future study.