Discrete Wigner Function Derivation of the Aaronson–Gottesman Tableau Algorithm

Abstract

The Gottesman–Knill theorem established that stabilizer states and Clifford operations can be efficiently simulated classically. For qudits with odd 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 for qubits. We show that the efficiency of both algorithms is due to harmonic evolution in the symplectic structure of discrete phase space. The differences between the Wigner function algorithm for odd-d and the Aaronson–Gottesman algorithm for qubits are likely due only to the fact that the Weyl–Heisenberg group is not in $SU\left(d\right)$ for $d=2$ and that qubits exhibit state-independent contextuality. This may provide a guide for extending the discrete Wigner function approach to qubits.

1. Introduction

The cost of brute-force classical simulation of the time evolution of n-qubit states grows exponentially with n. An important exception to this involves the set of Clifford operators acting on stabilizer states. This set of states plays an important role in quantum error correction [1] and is 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], Gross [7], Veitch et al. [8,9], Mari et al. [4], and Howard et al. [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. In odd-dimensional systems, 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 non-contextual and 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 [10]. This means Clifford circuits are expressible by path integrals truncated at order ${\hslash }^0$ and are thus manifestly classical [10,11].

This poses the question: what is the relationship between past efficient algorithms for Clifford circuits and the propagation of discrete Wigner functions of stabilizer states under Clifford operators? In the present paper, we show that the original Aaronson–Gottesman tableau algorithm for qubit stabilizer states 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 $\widehat{Z}$ Pauli basis for odd d. We then show its equivalence to the well-known Aaronson–Gottesman tableau algorithm [2] for qubits ($d=2$). Both algorithms require $\mathcal{O}\left(n^2\right)$ dits to represent n stabilizer states, $\mathcal{O}\left(n\right)$ operations per Clifford operator, and both deterministic and random measurements require $\mathcal{O}\left(n^2\right)$ operations.

The Aaronson–Gottesman tableau algorithm makes use of the Heisenberg representation. This means that time evolution is accomplished by updating an associated tableau or matrix representation of the Clifford operators instead of the stabilizer states themselves. The algorithm we present is framed in the Schrödinger picture and involves evolving the Wigner function of stabilizer states. By demonstrating that the two algorithms are equivalent, we show that the formulation of Clifford simulation in the Heisenberg picture is a choice and not a necessity for its efficient simulation. Furthermore, by instead working in the Schrödinger picture we are able to more easily reveal the purely classical basis of both algorithms and the physically intuitive phase space structures and symplectic properties on which they rely.

2. Discrete Wigner Function for Odd d Qudits

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 ${\widehat{Z}}_j$ operator on the jth qudit for n qudits acts as a “boost” operator:

$${\widehat{Z}}_j\left|k_1,\dots ,k_j,\dots ,k_n\right.〉=e^{\frac{2\pi i}{d}{k}_j}\left|k_1,\dots ,k_j,\dots ,k_n\right.〉,$$

(1)

where $k_j\in \mathbb{Z}/d\mathbb{Z}$ for $1\le j\le n$.

The discrete Fourier transform operator is defined by:

$$\begin{array}{c}{\widehat{F}}_j=\frac{1}{\sqrt{d}}\sum _{\begin{array}{c}{k}_j,l_j\in \\ \mathbb{Z}/d\mathbb{Z}\end{array}}{e}^{-\frac{2\pi i}{d}{k}_j{l}_j}\left|k_1,\dots ,k_j,\dots ,k_n\right.〉\left.〈l_1,\dots ,l_j,\dots ,l_n\right|.\hfill \end{array}$$

This is the d-dimensional equivalent of the Hadamard gate and allows us to define the Pauli ${\widehat{X}}_j$ operator as follows:

$${\widehat{X}}_j\equiv {\widehat{F}}_j{\widehat{Z}}_j{\widehat{F}}_j^{†}.$$

(2)

While ${\widehat{Z}}_j$ is a boost, ${\widehat{X}}_j$ is a shift operator because

$${\widehat{X}}_j^{\delta q}\left|k_1,\dots ,k_j,\dots ,k_n\right.〉\equiv \left|k_1,\dots ,k_j\oplus \delta q,\dots ,k_n\right.〉,$$

(3)

where ⊕ denotes integer addition mod d.

We can reexpress the boost ${\widehat{Z}}_j$ and shift ${\widehat{X}}_j$ operators in terms of their generators, which are the conjugate ${\widehat{q}}_j$ and ${\widehat{p}}_j$ operators, respectively:

$${\widehat{Z}}_j=e^{\frac{2\pi i}{d}{\widehat{q}}_j}$$

(4)

and

$${\widehat{X}}_j=e^{-\frac{2\pi i}{d}{\widehat{p}}_j}.$$

(5)

Thus, we can refer to the ${\widehat{X}}_j$ basis as the momentum ($p_j$) basis, which is equivalent to the Fourier transform of the $q_j$ basis:

$${\widehat{p}}_j={\widehat{F}}_j{\widehat{q}}_j{\widehat{F}}_j^{†}.$$

(6)

These bases form the discrete Weyl phase space $(\mathit{p},\mathit{q})$.

The Wigner function $W_{\mathsf{\Psi }}(\mathit{p},\mathit{q})$ of a pure state $\left|\mathsf{\Psi }\right.〉$ is defined on this discrete Weyl phase space:

$$\begin{array}{ccc}\hfill W_{\mathsf{\Psi }}(\mathit{p},\mathit{q})& =& d^{-n}\sum _{\begin{array}{c}{\mathit{\xi }}_q\in \\ {(\mathbb{Z}/d\mathbb{Z})}^n\end{array}}{e}^{-\frac{2\pi i}{d}{\mathit{\xi }}_q·\mathit{p}}\mathsf{\Psi }\left(\mathit{q}+\frac{(d+1){\mathit{\xi }}_q}{2}\right){\mathsf{\Psi }}^{*}\left(\mathit{q}-\frac{(d+1){\mathit{\xi }}_q}{2}\right).\hfill \end{array}$$

(7)

This is equivalent to the discrete Wigner function introduced by Gross [7]. We will shortly be interested in the discrete Wigner function of stabilizer states. However, first, we introduce the effect that the Clifford gates have in this discrete Weyl phase space.

2.1. Clifford Gates

A Clifford group gate $\widehat{V}$ is related to a symplectic transformation on the discrete Weyl phase space, governed by a symplectic matrix ${\mathcal{M}}_{\widehat{V}}$ and vector ${\mathit{\alpha }}_{\widehat{V}}$ [7]:

$$\left(\begin{array}{c}{\mathit{p}}^{\prime }\\ {\mathit{q}}^{\prime }\end{array}\right)={\mathcal{M}}_{\widehat{V}}\left[\left(\begin{array}{c}\mathit{p}\\ \mathit{q}\end{array}\right)+\frac{1}{2}{\mathit{\alpha }}_{\widehat{V}}\right]+\frac{1}{2}{\mathit{\alpha }}_{\widehat{V}}.$$

(8)

Wigner functions $W_{\mathsf{\Psi }}\left(\mathit{x}\right)$ of states evolve under Clifford operators $\widehat{V}$ by

$$W_{\mathsf{\Psi }}\left({\mathcal{M}}_{\widehat{V}}\left(\mathit{x}+{\mathit{\alpha }}_{\widehat{V}}/2\right)+{\mathit{\alpha }}_{\widehat{V}}/2\right),$$

(9)

where $\mathit{x}\equiv (\mathit{p},\mathit{q})$. When considering Clifford gate propagation, we can restrict to a set of gates which are generators of the Clifford group. One such set of generators is made up of the phase-shift gate ${\widehat{P}}_i$, the Hadamard gate ${\widehat{F}}_i$, and the controlled-not (CNOT) ${\widehat{C}}_{ij}$ (which act on the ith and jth qudits).

The phase shift ${\widehat{P}}_i$ is a one-qudit gate with the underlying Hamiltonian $H_{{\widehat{P}}_i}=-\frac{d+1}{2}{q}_i^2+\frac{d+1}{2}{q}_i$ [10]. Without loss of generality, we will instead consider

$${\widehat{P}}_i^{\prime }={\widehat{P}}_i{\widehat{P}}_i{\widehat{Z}}_i,$$

(10)

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:

$${\widehat{P}}_i={\widehat{P}}_i^{\prime }{\widehat{P}}_i^{\prime }{\widehat{Z}}_i,$$

(11)

where $[{\widehat{P}}_i,{\widehat{Z}}_i]=[{\widehat{P}}_i^{\prime },{\widehat{Z}}_i]=0$. Hence, ${\widehat{P}}_i^{\prime }$ is an adequate replacement generator for ${\widehat{P}}_i$, and we will use it instead of ${\widehat{P}}_i$ from now on. Since its Hamiltonian has no linear term ($H_{{\widehat{P}}_i^{\prime }}=-q_i^2$), this leads to an easier presentation ahead since ${\mathit{\alpha }}_{{\widehat{P}}_i^{\prime }}=\mathbf{0}$. The corresponding equations of motion for ${\widehat{P}}_i^{\prime }$ are ${\dot{p}}_i=2q_i$ and ${\dot{q}}_i=0$. Hence, for $\Delta t=1$,

$${\left({\mathcal{M}}_{{\widehat{P}}_i^{\prime }}\right)}_{j,k}={\delta }_{j,k}+2{\delta }_{i,j}{\delta }_{n+i,k}.$$

(12)

The Hadamard gate ${\widehat{F}}_i$ is a one-qudit gate and has the underlying Hamiltonian $H_{{\widehat{F}}_i}=-\frac{\pi }{4}(p_i^2+q_i^2)$ [10]. The corresponding equations of motion are ${\dot{p}}_i=\frac{\pi }{2}{q}_i$ and ${\dot{q}}_i=-\frac{\pi }{2}{p}_i$. Hence, for $\Delta t=1$,

$$\begin{array}{c}\hfill {\left({\mathcal{M}}_{{\widehat{F}}_i}\right)}_{j,k}={\delta }_{j,k}-{\delta }_{i,j}{\delta }_{i,k}-{\delta }_{n+i,j}{\delta }_{n+i,k}\\ \hfill +{\delta }_{i,j}{\delta }_{n+i,k}-{\delta }_{n+i,j}{\delta }_{i,k},\end{array}$$

(13)

and ${\mathit{\alpha }}_{{\widehat{F}}_i}=\mathbf{0}$.

Finally, the two-qudit CNOT ${\widehat{C}}_{ij}$ on control qudit i and second qudit j has the corresponding Hamiltonian $H_{{\widehat{C}}_{ij}}=p_i{q}_j$ [10]. The corresponding equations of motion are $({\dot{p}}_i,{\dot{p}}_j)=-(0,p_i)$ and $({\dot{q}}_i,{\dot{q}}_j)=(q_j,0)$. Hence, for $\Delta t=1$,

$${\left({\mathcal{M}}_{{\widehat{C}}_{ij}}\right)}_{k,l}={\delta }_{k,l}-{\delta }_{i,k}{\delta }_{j,l}+{\delta }_{n+j,k}{\delta }_{n+i,l},$$

(14)

and ${\mathit{\alpha }}_{{\widehat{C}}_{ij}}=\mathbf{0}$.

2.2. Wigner Functions of Stabilizer States

A discrete Wigner function for stabilizer states associated with the boost and shift operators defined in Equations (4) and (5) is given by the following theorem [10]:

Theorem 1.

The discrete Wigner function $W_{\mathsf{\Psi }}\left(\mathit{x}\right)$ of a stabilizer state Ψ for any odd d and n qudits is ${\delta }_{\mathsf{\Phi }\times \mathit{x},\mathit{r}}$ for $2n\times 2n$ matrix Φ and $2n$ vector $\mathit{r}$ with entries in $\mathbb{Z}/d\mathbb{Z}$.

An equivalent form was proven by Gross [7] who also showed that these discrete Wigner functions of stabilizer states are non-negative. In particular, if we begin with a stabilizer state defined as $\left|{\mathsf{\Psi }}_0\right.〉=\left|{\mathit{q}}_0\right.〉$, then $W_{{\mathsf{\Psi }}_0}\left(\mathit{x}\right)={\delta }_{{\mathsf{\Phi }}_0\times \mathit{x},{\mathit{r}}_0}$, where ${\mathsf{\Phi }}_0=\left(\begin{array}{cc}0& 0\\ 0& {\mathbb{I}}_n\end{array}\right)$ for ${\mathbb{I}}_n$ the $n\times n$ identity matrix, and ${\mathit{r}}_0=(\mathbf{0},{\mathit{q}}_0)$.

3. Wigner Stabilizer Algorithm for Odd d Qudits

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 Equation (9), we can now examine the effect Clifford operators have on stabilizer states. We note that since the discrete Wigner functions of stabilizer states are non-negative and Clifford operations take stabilizer states to stabilizer states, it follows that Clifford operations (if associated positive-operator valued measures (POVMs) also have non-negative Wigner functions) can always be efficiently classically simulated by sampling from these Wigner functions as probability distributions [4]. However, here we pursue a description that is not dependent on classical sampling.

3.1. Stabilizer Representation

From Theorem 1, propagation of the stabilizer state $\mathsf{\Psi }$ can be represented by considering the state’s Wigner function: $W_{\mathsf{\Psi }}\left(\mathit{x}\right)={\delta }_{{\mathsf{\Phi }}_t·\mathit{x},{\mathit{r}}_t}$. In this way, ${\mathsf{\Phi }}_t$ and ${\mathit{r}}_t$ specify a linear system of equations in terms of ${\mathit{p}}_t$ and ${\mathit{q}}_t$. The first n rows of ${\mathsf{\Phi }}_t$ are the coefficients of ${({\mathit{p}}_t,{\mathit{q}}_t)}^T$ in ${\mathit{p}}_0({\mathit{p}}_t,{\mathit{q}}_t)$ and the last n rows of ${\mathsf{\Phi }}_t$ are the coefficients of ${({\mathit{p}}_t,{\mathit{q}}_t)}^T$ in ${\mathit{q}}_0({\mathit{p}}_t,{\mathit{q}}_t)$:

$$\left(\begin{array}{c}{\mathit{p}}_0\\ {\mathit{q}}_0\end{array}\right)={\mathsf{\Phi }}_t\left(\begin{array}{c}{\mathit{p}}_t\\ {\mathit{q}}_t\end{array}\right).$$

(15)

The Kronecker delta function sets this linear system of equations equal to ${\mathit{r}}_t$. In this way, an affine map—a linear transformation displaced from the origin by ${\mathit{r}}_t$—is defined. This system of equations must be updated after every unitary propagation and measurement.

Since the Wigner functions $W_{\mathsf{\Psi }}\left(\mathit{x}\right)$ of stabilizer states propagate under $\mathcal{M}$ as $W_{\mathsf{\Psi }}\left(\mathcal{M}\mathit{x}\right)$, it follows that

$${\mathsf{\Phi }}_t\to {\mathsf{\Phi }}_t{\mathcal{M}}_t^{-1}.$$

(16)

(The importance of vector ${\mathit{r}}_t$ and when it must be updated will become evident when we consider random measurements.) Hence, after n operations ${\mathcal{M}}_1$, ${\mathcal{M}}_2$, …, ${\mathcal{M}}_n$,

$${\mathcal{M}}_t^{-1}={\mathcal{M}}_1^{-1}{\mathcal{M}}_2^{-1}\dots {\mathcal{M}}_n^{-1}.$$

(17)

The matrices are ordered chronologically left-to-right instead of right-to-left.

Since $\mathcal{M}$ is symplectic, ${\mathcal{M}}_t^{-1}=-\mathcal{J}{\mathcal{M}}_t^T\mathcal{J}$ where

$$\mathcal{J}=\left(\begin{array}{cc}0& -{\mathbb{I}}_n\\ {\mathbb{I}}_n& 0\end{array}\right).$$

(18)

Thus, the the stability matrices $\mathcal{M}$ for ${\widehat{F}}_i$, ${\widehat{P}}_i^{\prime }$ and ${\widehat{C}}_{ij}$ given in Equations (12)–(14) differ from their inverses only by sign changes in their off-diagonal elements:

$${\left({\mathcal{M}}_{{\widehat{P}}_i^{\prime }}^{-1}\right)}_{j,k}={\delta }_{j,k}-2{\delta }_{i,j}{\delta }_{n+i,k},$$

(19)

$$\begin{array}{c}\hfill {\left({\mathcal{M}}_{{\widehat{F}}_i}^{-1}\right)}_{j,k}={\delta }_{j,k}-{\delta }_{i,j}{\delta }_{i,k}-{\delta }_{n+i,j}{\delta }_{n+i,k}\\ \hfill -{\delta }_{i,j}{\delta }_{n+i,k}+{\delta }_{n+i,j}{\delta }_{i,k},\end{array}$$

(20)

and

$${\left({\mathcal{M}}_{{\widehat{C}}_{ij}}^{-1}\right)}_{k,l}={\delta }_{k,l}+{\delta }_{i,k}{\delta }_{j,l}-{\delta }_{n+j,k}{\delta }_{n+i,l}.$$

(21)

We assume the quantum state is initialized in the computational basis state ${\mathsf{\Psi }}_0={\underbrace{\left|0\right.〉\otimes \cdots \otimes \left|0\right.〉}}_n$ and so initially we should set ${\mathsf{\Phi }}_0=\left(\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbb{I}}_n\end{array}\right)$ and ${\mathit{r}}_0=\mathbf{0}$. The initial stabilizer state is $W_{{\mathsf{\Psi }}_0}={\delta }_{{\mathit{q}}_t,\mathbf{0}}$. However, it will become clear when we discuss measurements that it is practically useful to instead set

$${\mathsf{\Phi }}_0=\left(\begin{array}{cc}{\mathbb{I}}_n& \mathbf{0}\\ \mathbf{0}& {\mathbb{I}}_n\end{array}\right),$$

(22)

thereby setting $W_{{\mathsf{\Psi }}_0}={\delta }_{({\mathit{p}}_t,{\mathit{q}}_t),(\mathbf{0},\mathbf{0})}$—not a true Wigner function. This new matrix ${\mathsf{\Phi }}_0$ is equivalent to the last matrix if the first n rows in ${\mathsf{\Phi }}_t\mathit{x}$ and ${\mathit{r}}_t$ are ignored—the same as ignoring ${\mathit{p}}_0({\mathit{p}}_t,{\mathit{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 $\mathcal{O}\left(n^2\right)$ dits of storage to track ${\mathsf{\Phi }}_t$ and ${\mathit{r}}_t$. More precisely, since ${\mathsf{\Phi }}_t$ is a $2n\times 2n$ matrix and ${\mathit{r}}_t$ is an $2n$-vector, $2n(2n+1)$ dits of storage are necessary.

3.2. Unitary Propagation

${\mathsf{\Phi }}_t$ contains the coefficients of the linear equations relating ${\mathit{x}}_0$ to ${\mathit{x}}_t$. Each row is one equation relating ${q_0}_i$ or ${p_0}_i$ to ${\mathit{x}}_t$. When manipulating rows of ${\mathsf{\Phi }}_t$ we shall refer to the linear equations that these rows define.

Examining Equations (19)–(21), we see that the inverse stability matrices of the generator gates ${\widehat{F}}_i$, ${\widehat{P}}_i$ and ${\widehat{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 ${\mathsf{\Phi }}_t$ with a new stability matrix in Equation (16) and evaluating the matrix multiplication is equivalent to performing a finite number of n-vector dot products and so requires $\mathcal{O}\left(n\right)$ operations. Therefore, keeping track of propagation of stabilizer states by Clifford gates can be simulated with $\mathcal{O}\left(n\right)$ 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 (${\widehat{P}}_i^{\prime }$). For all $j\in \{1,\dots ,2n\}$, set ${\mathsf{\Phi }}_{j,n+i}\mapsto {\mathsf{\Phi }}_{j,n+i}\ominus 2{\mathsf{\Phi }}_{j,n}$.

Hadamard gate on qudit i (${\widehat{F}}_i$). For all $j\in \{1,\dots ,2n\}$, negate ${\mathsf{\Phi }}_{j,i}$ mod d, and then swap $\ominus {\mathsf{\Phi }}_{j,i}$ and ${\mathsf{\Phi }}_{j,n+i}$.

CNOT from control i to target j (${\widehat{C}}_{ij}$). For all $j\in \{1,\dots ,2n\}$, set ${\mathsf{\Phi }}_{k,j}\mapsto {\mathsf{\Phi }}_{k,j}\oplus {\mathsf{\Phi }}_{k,i}$ and ${\mathsf{\Phi }}_{k,n+i}\mapsto {\mathsf{\Phi }}_{k,n+i}\ominus {\mathsf{\Phi }}_{k,n+j}$.

This confirms that unitary propagation in this scheme requires $\mathcal{O}\left(n\right)$ operations.

3.3. Measurement

The outcome of a measurement ${\widehat{Z}}_i$ on a stabilizer state can be either random or deterministic. As described above, the bottom half of ${\mathsf{\Phi }}_t$ defines ${q_0}_j$ for $j\in \{1,\dots ,n\}$, each of which is a linear combination of ${q_t}_i$ and ${p_t}_i$. The entries in the $(n+j)$th row of ${\mathsf{\Phi }}_t$ give the coefficient of ${p_t}_i$ and ${q_t}_i$ in ${q_0}_j$ for $j\in \{1,\dots ,n\}$. If the coefficient of ${p_t}_i$ in any ${q_0}_j$ is non-zero then the measurement ${\widehat{Z}}_i$ will be random. If all coefficients of ${p_t}_i$ are zero for ${q_0}_j$ $\forall j$, then the measurement of ${\widehat{Z}}_i$ will be deterministic. This can be seen from the fact that if our stabilizer state $\left|\mathsf{\Psi }\right.〉$ is an eigenstate of ${\widehat{Z}}_i$, then ${\widehat{Z}}_i\left|\mathsf{\Psi }\right.〉=e^{i\varphi }\left|\mathsf{\Psi }\right.〉$ for some $\varphi \in \mathbb{R}$ and (discrete) Wigner functions do not change under a global phase. Thus, measuring ${\widehat{Z}}_i$ leaves the Wigner function of $\left|\mathsf{\Psi }\right.〉$ invariant if the measurement is deterministic. Since ${\widehat{Z}}_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:

$${\mathsf{\Phi }}_t\left(\begin{array}{c}{p_t}_1\\ ⋮\\ {p_t}_i\\ ⋮\\ {p_t}_n\\ {\mathit{q}}_t\end{array}\right)=\left(\begin{array}{c}{\mathit{r}}_{\mathbf{t}}p\\ {\mathit{r}}_{\mathbf{t}}q\end{array}\right)\underset{{\widehat{Z}}_i}{\mapsto }{\mathsf{\Phi }}_t\left(\begin{array}{c}{p_t}_1\\ ⋮\\ {p_t}_i+1\\ ⋮\\ {p_t}_n\\ {\mathit{q}}_t\end{array}\right)=\left(\begin{array}{c}{\mathit{r}}_{\mathbf{t}}p\\ {\mathit{r}}_{\mathbf{t}}q\end{array}\right).$$

(23)

Thus, if the lower half of the ith column of ${\mathsf{\Phi }}_t$ is zero, then ${\widehat{Z}}_i$ leaves the Wigner function invariant (and so the measurement is deterministic). Verifying that these coefficients are all zero takes $\mathcal{O}\left(n\right)$ operations for each ${\widehat{Z}}_i$.

In other words, to see if a given measurement of ${\widehat{Z}}_i$ is random or deterministic, a search must be performed for non-zero ${{\mathsf{\Phi }}_t}_{n+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 ${{\mathsf{\Phi }}_t}_{n+j,i}$ element exists, then the measurement ${\widehat{Z}}_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 ${\mathsf{\Phi }}_t$ have a non-zero entry in the ith column, ${{\mathsf{\Phi }}_t}_{n+j,i}\ne 0$. Since the random measurement ${\widehat{Z}}_i$ will project qudit i onto a position state, we will replace the $(n+j)$th row with ${q_0}_i=q_i^{\prime }$ (the uniformly random outcome of this measurement). After this projection onto a position state, none of the other qudits’ positions should depend on qudit i’s momentum, ${p_t}_i$. To accomplish this, before we replace row $(n+j)$, we solve its equation for ${p_t}_i$ and substitute every instance of ${p_t}_i$ in the linear system of equations with this solution. As a result, every equation will no longer depend on ${p_t}_i$ and we can go ahead and replace the $(n+j)$th row with ${q_0}_i=q_i^{\prime }$.

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_0}_i={q_0}_j({\mathit{p}}_t,{\mathit{q}}_t)$, which becomes the only remaining equation explicitly dependent on ${p_t}_i$. In other words, ${p_0}_i\propto {p_t}_i$, similar to the beginning when we set ${p_0}_i={p_t}_i$ by setting $\mathsf{\Phi }={\mathbb{I}}_{2n}$. However, now we also preserve any dependence ${p_0}_i$ has on the other qudits incurred during unitary propagation. In other words, we preserve ${p_t}_i$’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 ${\mathsf{\Phi }}_t$ and ${\mathit{r}}_t$ with a randomly chosen measurement outcome $q_i^{\prime }$ (i.e., ${q_0}_i=q_i^{\prime }$), the identification of rows $(n+i)$ and $(n+j)$ are exchanged, so that the former now specifies ${q_0}_j({\mathit{p}}_t,{\mathit{q}}_t)$ while the latter specifies ${q_0}_i({\mathit{p}}_t,{\mathit{q}}_t)$. ${p_0}_i$ has also been updated by replacing the jth row in the first half of ${\mathsf{\Phi }}_t$, with the $(n+j)$th row we just changed. Again, this row now describes ${p_0}_i({\mathit{p}}_t,{\mathit{q}}_t)$ while the ith row now specifies ${p_0}_j({\mathit{p}}_t,{\mathit{q}}_t)$. Overall, this takes $\mathcal{O}\left(n^2\right)$ operations since we are replacing $\mathcal{O}\left(n\right)$ rows with $\mathcal{O}\left(n\right)$ entries.

Case 2: Deterministic Measurement

Since the measurement is deterministic, ${\mathsf{\Phi }}_t$ and ${\mathit{r}}_t$ do not change. The n equations specified by the bottom half of ${\mathsf{\Phi }}_t{\mathit{x}}_t={\mathit{r}}_t$ can be used to solve for ${q_t}_i$—the deterministic measurement outcome. In general, this can also be done by inverting ${\mathsf{\Psi }}_t$ and evaluating ${\mathit{x}}_t={\mathsf{\Phi }}_t^{-1}·{\mathit{r}}_t$ for $q_i$. Aaronson and Gottesman themselves noted that such a matrix inversion is possible, but practically takes $\mathcal{O}\left(n^3\right)$ operations. (However, we are not certain if Aaronson and Gottesman were referring to the ${\mathsf{\Phi }}_t$ matrix corresponding to the $2n\times 2n$ part of their tableau when they discuss matrix inversion in [2].)

Fortunately, there is another method that scales as $\mathcal{O}\left(n^2\right)$ and requires use of the n equations represented by the top n rows of ${\mathsf{\Phi }}_t$, which were included in our description by setting ${\mathsf{\Phi }}_0={\mathbb{I}}_{2n}$. The linear system of n equations represented by ${\mathsf{\Phi }}_t{\mathit{x}}_t={\mathit{r}}_t$ can be written as

$$\begin{array}{ccc}\hfill {\mathsf{\Phi }}_t{\mathit{x}}_t& =& {\mathit{r}}_t,\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{\mathit{p}}_0({\mathit{p}}_t,{\mathit{q}}_t)\\ {\mathit{q}}_0({\mathit{p}}_t,{\mathit{q}}_t)\end{array}\right)& =& \left(\begin{array}{c}{\mathit{r}}_{\mathbf{t}}p\\ {\mathit{r}}_{\mathbf{t}}q\end{array}\right),\hfill \end{array}$$

(25)

where we are interested in linear combinations of the bottom half, ${\mathit{q}}_0({\mathit{p}}_t,{\mathit{q}}_t)$, to solve for the measurement outcome ${q_t}_i$:

$$\sum _{j=1}^n{c}_{ij}{q_0}_j={q_t}_i,$$

(26)

where $c_{ij}\in \mathbb{Z}/d\mathbb{Z}$.

Lemma 1.

The coefficient in front of ${p_t}_i$ in the row of ${\mathsf{\Phi }}_t$ that specifies ${p_0}_j({\mathit{p}}_t,{\mathit{q}}_t)$, ${{\mathsf{\Phi }}_t}_{ji}$, is equal to the coefficient $c_{ij}$ in front of ${q_0}_j$ that makes up ${q_t}_i$ in Equation (26). Equivalently,

$$c_{ij}={q_0}_j·{q_t}_i({\mathit{p}}_0,{\mathit{q}}_0)={p_0}_j({\mathit{p}}_t,{\mathit{q}}_t)·{p_t}_i={{\mathsf{\Phi }}_t}_{ji}.$$

(27)

Proof. 

Under evolution under the Clifford group operators,

$$\left(\begin{array}{c}{\mathit{p}}_t\\ {\mathit{q}}_t\end{array}\right)={\mathcal{M}}_t\left(\begin{array}{c}{\mathit{p}}_0\\ {\mathit{q}}_0\end{array}\right).$$

(28)

${\mathcal{M}}_t^{-1}=-\mathcal{J}{\mathcal{M}}_t^T\mathcal{J}$ since ${\mathcal{M}}_t$ is symplectic. This means that we can express the matrix inversion as follows:

$$\begin{array}{}\mathrm{(29)}& \hfill \left(\begin{array}{c}{\mathit{p}}_0\\ {\mathit{q}}_0\end{array}\right)& =& {\mathcal{M}}_t^{-1}\left(\begin{array}{c}{\mathit{p}}_t\\ {\mathit{q}}_t\end{array}\right)\hfill \mathrm{(30)}& & =& -\mathcal{J}{\mathcal{M}}_t^T\mathcal{J}\left(\begin{array}{c}{\mathit{p}}_t\\ {\mathit{q}}_t\end{array}\right)\hfill \mathrm{(31)}& & =& -\mathcal{J}{\left(\begin{array}{cc}{\left({\mathcal{M}}_t\right)}_{11}& {\left({\mathcal{M}}_t\right)}_{12}\\ {\left({\mathcal{M}}_t\right)}_{21}& {\left({\mathcal{M}}_t\right)}_{22}\end{array}\right)}^T\mathcal{J}\left(\begin{array}{c}{\mathit{p}}_t\\ {\mathit{q}}_t\end{array}\right)\hfill \mathrm{(32)}& & =& \left(\begin{array}{cc}{\left({\mathcal{M}}_t\right)}_{22}& {\left(-{\mathcal{M}}_t\right)}_{12}\\ {\left(-{\mathcal{M}}_t\right)}_{21}& {\left({\mathcal{M}}_t\right)}_{11}\end{array}\right)\left(\begin{array}{c}{\mathit{p}}_t\\ {\mathit{q}}_t\end{array}\right).\hfill \end{array}$$

Therefore, ${{\left({\mathcal{M}}_t^{-1}\right)}_{11}}_{i,j}={{\left({\mathcal{M}}_t\right)}_{22}}_{i,j}$, and so

$$c_{ij}={q_0}_j·{q_t}_i({\mathit{p}}_0,{\mathit{q}}_0)={p_0}_j({\mathit{p}}_t,{\mathit{q}}_t)·{p_t}_i={{\mathsf{\Phi }}_t}_{ji}.$$

(33)

This property can also be seen in the drawing of phase space shown in Figure 1. There, initial perpendicular ${p_0}_j$ and ${q_0}_j$ manifolds are drawn along with harmonically evolved ${p_t}_i$ and ${q_t}_i$ manifolds, which remain perpendicular to each other and make an angle $\alpha$ to the first ${p_0}_j$ and ${q_0}_j$ manifolds, respectively. The projection of ${q_t}_i({\mathit{p}}_0,{\mathit{q}}_0)$ onto ${q_0}_j$ can be represented as the length b of a right triangle’s adjacent side to the angle $\alpha$, with an opposite side set to some length a. The projection of ${p_0}_j({\mathit{p}}_t,{\mathit{q}}_t)$ onto ${p_t}_i$ is similarly represented by the length $b^{\prime }$ of a right triangle’s adjacent side to the angle $\alpha$, with an opposite side also set to length a. It follows that the third angle $\beta$ in both triangles must be the same, and so by the law of sines

$$\frac{a}{sin\alpha }=\frac{b}{sin\beta }=\frac{b^{\prime }}{sin\beta }.$$

(34)

Therefore, $b=b^{\prime }$ 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 $\mathcal{O}\left(n^2\right)$ operations since Equation (27) is a sum of $\mathcal{O}\left(n\right)$ vectors made up of $\mathcal{O}\left(n\right)$ components. Therefore, the overall measurement protocol takes $\mathcal{O}\left(n^2\right)$ 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 (Equation (32)) to be performed efficiently.

4. Aaronson–Gottesman Tableau Algorithm for Qubits ($d=2$)

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 $\mathcal{O}\left(n^2\right)$ operations for propagation and measurement for n qubits. The algorithm has been proven to be extendable to $d>2$ [12] and similar algorithms have been formulated in $d>2$ [13]. Alternatives have also been developed to the tableau formalism, though they prove to be equally efficient in worst-case scenarios [14]. 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 3 is equivalent to the Aaronson–Gottesman tableau algorithm extended to odd d.

4.1. Stabilizer Representation

The Aaronson–Gottesman algorithm is defined in 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=\{\widehat{g}\in \mathcal{P}$ such that $\widehat{g}\left|\psi \right.〉=\left|\psi \right.〉\}$ are called the stabilizers of state $\left|\psi \right.〉$, where $\mathcal{P}$ is the set of Pauli operators, each of which has the form $e^{\frac{\pi i}{2}\alpha }{\widehat{P}}_1\otimes \cdots \otimes {\widehat{P}}_n$ where $\alpha \in \{0,1,2,3\}$ for n qubits with ${\widehat{P}}_i\in \{{\widehat{I}}_i,{\widehat{Z}}_i,{\widehat{X}}_i,{\widehat{Y}}_i\}$.

For the sake of completeness, we present here a summary of the qubit Aaronson–Gottesman algorithm, in order to compare it to our odd d qudit 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 n-qubit stabilizer state is defined by the n generators of its stabilizer state. Every element in this set of generators, $\{{\widehat{g}}_1,{\widehat{g}}_2,\dots ,{\widehat{g}}_n\}$, is in the Pauli group, and each generator has the form:

$${\widehat{g}}_i=\pm {\widehat{P}}_{i1}\dots {\widehat{P}}_{in}.$$

(35)

Any unitary propagation by Clifford operators or measurement of the stabilizer state changes at least some of the ${\widehat{P}}_{ij}$ elements of the n generators of the state’s stabilizer. This includes the $\pm 1$ phase in Equation (35), which must also be kept track of in Aaronson–Gottesman’s algorithm.

4.2. Unitary Propagation

For each Clifford operation, Aaronson and Gottesman showed that only $\mathcal{O}\left(n\right)$ operations are necessary to update all generators [2]. Specifically, according to the update rules in

Table 1. Transformation of stabilizer generators under Clifford operations.

Gates	Input	Output
Hadamard	$\widehat{X}$	$\widehat{Z}$
	$\widehat{Z}$	$\widehat{X}$
phase	$\widehat{X}$	$\widehat{Y}$
	$\widehat{Z}$	$\widehat{Z}$
CNOT	$\widehat{X}\otimes \widehat{I}$	$\widehat{X}\otimes \widehat{X}$
	$\widehat{I}\otimes \widehat{X}$	$\widehat{I}\otimes \widehat{X}$
	$\widehat{Z}\otimes \widehat{I}$	$\widehat{Z}\otimes \widehat{I}$
	$\widehat{I}\otimes \widehat{Z}$	$\widehat{Z}\otimes \widehat{Z}$

, each generator can be updated with a constant number of operators for every single Clifford gate, therefore $O\left(n\right)$ 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 $\{{\widehat{g}}_1^{\prime },\dots ,{\widehat{g}}_n^{\prime }\}$ are the operators that generate the full Pauli group with the stabilizers $\{{\widehat{g}}_1,\dots ,{\widehat{g}}_n\}$. They have the following properties:

(i) 

${\widehat{g}}_1^{\prime }$, ${\widehat{g}}_2^{\prime }$, …, ${\widehat{g}}_n^{\prime }$ commute.

(ii) 

Each destabilizer ${\widehat{g}}_h^{\prime }$ anti-commutes with the corresponding stabilizer ${\widehat{g}}_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\times (2n+1)$ binary tableau matrix as:

$$\left(\begin{array}{ccccccc}{x}_{11}& \dots & x_{1n}& z_{11}& \dots & z_{1n}& r_1\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮\\ x_{n1}& \dots & x_{nn}& z_{n1}& \dots & z_{nn}& r_n\\ x_{(n+1)1}& \dots & x_{(n+1)n}& z_{(n+1)1}& \dots & z_{(n+1)n}& r_{n+1}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮\\ x_{\left(2n\right)1}& \dots & x_{\left(2n\right)n}& z_{\left(2n\right)1}& \dots & z_{\left(2n\right)n}& r_{2n}\end{array}\right).$$

This matrix contains $2n$ rows. The first n rows denote the destabilizers ${\widehat{g}}_1^{\prime }$ to ${\widehat{g}}_n^{\prime }$ while rows $(n+1)$ to $2n$ represent the stabilizers ${\widehat{g}}_1$ to ${\widehat{g}}_n$. The $(n+1)$th bit in each row denotes the phase ${(-1)}^{r_i}$ for each generator. We encode the jth Pauli operator in the ith row as shown in

Table 2. Binary representation of the Pauli operators and the Pauli group phase used in their tableau representation.

${\mathit{x}}_{\mathit{ij}}$	${\mathit{z}}_{\mathit{ij}}$	${\widehat{\mathit{P}}}_{\mathit{j}}$
0	0	${\widehat{I}}_j$
0	1	${\widehat{Z}}_j$
1	0	${\widehat{X}}_j$
1	1	${\widehat{Y}}_j$
${\mathit{r}}_{\mathit{i}}$	Phase
0	$+1$
1	$-1$

.

We can update the stabilizers and destabilizers as follows:

Hadamard gate on qubit i For all $j\in \{1,2,\dots ,2n\}$, $r_j\mapsto r_j\oplus x_{ji}{z}_{ji}$, then swap $x_{ji}$ with $z_{ji}$.

Phase gate on qubit i For all $j\in \{1,2,\dots ,2n\}$, $r_j\mapsto r_j\oplus x_{ji}{z}_{ji}$, $z_{ji}\mapsto z_{ji}\oplus x_{ji}$.

CNOT gate on control qubit i and target qubit j For all $k\in \{1,2,\dots ,2n\}$, $r_k\mapsto r_k\oplus x_{ki}{z}_{kj}(x_{kj}\oplus z_{ki}\oplus 1)$, $x_{kj}\mapsto x_{kj}\oplus x_{ki}$, $z_{ki}\mapsto z_{ki}\oplus z_{kj}$.

These actions correspond to those given in Table 1.

Notice the striking similarity of these tableau transformation rules under unitary propagation to the $\mathsf{\Phi }$ transformation rules in Section 4. The most notable difference is that the Aaronson–Gottesman algorithm involves updates of the vector $\mathit{r}$. We will discuss this and its connection to the dimension $d=2$ of the system in Section 5. It is clear that these transformations also take $O\left(n\right)$ operations each.

4.3. Measurement

To describe the measurement part of the algorithm, we need to first define a rowsum operation in the tableau that 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}\oplus x_{jk}$ and $z_{ik}\oplus z_{jk}$ for $k=1,\dots ,n$. To calculate the resultant phase, Aaronson and Gottesman first defined the following function:

$$f(x_{ik},x_{jk},z_{ik},z_{jk})=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{x}_{ik}=z_{ik}=0,\hfill \\ z_{jk}-x_{jk}\hfill & \mathrm{if}{x}_{ik}=z_{ik}=1,\hfill \\ z_{jk}(2x_{jk}-1)\hfill & \mathrm{if}{x}_{ik}=1,z_{ik}=0,\hfill \\ x_{jk}(1-2z_{jk})\hfill & \mathrm{if}{x}_{ik}=0,z_{ik}=1.\hfill \end{array}\right.$$

(36)

Since each stabilizer generator is the tensor product of n single qubit Pauli operators (see Equation (35)), they must be multiplied together to obtain the phase:

$$\left\{\begin{array}{cc}0\hfill & \mathrm{if}{r}_i+r_j+\sum _{k=1}^{n}f(x_{ik},x_{jk},z_{ik},z_{jk})\equiv 0\left(\mathrm{mod}4\right),\hfill \\ 1\hfill & \mathrm{if}{r}_i+r_j+\sum _{k=1}^{n}f(x_{ik},x_{jk},z_{ik},z_{jk})\equiv 2\left(\mathrm{mod}4\right).\hfill \end{array}\right.$$

(37)

Having defined the rowsum function, let us now consider a measurement of ${\widehat{Z}}_i$ on qubit i. For $d=2$, Pauli group operators can only commute or anti-commute with each other. If ${\widehat{Z}}_i$ anti-commutes with one or more of the generators, then the measurement is random. If ${\widehat{Z}}_i$ commutes with all of the generators, then the measurement is deterministic. We consider these two cases:

Case 1: Random Measurement

${\widehat{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, ${\widehat{g}}_j$, and update the rest by replacing them with their product with ${\widehat{g}}_j$ (i.e., taking the rowsum of their corresponding rows) such that they commute with ${\widehat{Z}}_i$. These updates take $\mathcal{O}\left(n^2\right)$ operations. Finally, we only need to replace ${\widehat{g}}_j$ by ${\widehat{Z}}_i$.

In other words, with respect to the tableau, there should exist at least one $j\in \{n+1,n+2,...,2n\}$ such that $x_{ji}=1$. Replacing all rows where $x_{ki}=1$ for $k\ne j$ with the sum of the jth and kth row (using the rowsum function) sets all $x_{ki}=0$ for $k\ne 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\left(n^2\right)$ operations because each rowsum operation takes $O\left(n\right)$ operations and up to $n-1$ rowsums may be necessary.

Case 2: Deterministic Measurement

${\widehat{Z}}_i$ commutes with all generators. In this case, there is no $j\in \{n+1,n+2,...,2n\}$ such that $x_{ji}=1$ and we don’t need to update any of the generators. However, we do need to do some work to retrieve the measurement outcome.

Measurement ${\widehat{Z}}_i$ commutes with all of the stabilizers; therefore, either $+{\widehat{Z}}_i$ or $-{\widehat{Z}}_i$ is a stabilizer of the state. Therefore, it must be generated by the generators. The sign $\pm 1$ is the measurement outcome we are looking for. This means that

$$\prod _{j=1}^n{\widehat{g}}_j^{c_j}=\pm {\widehat{Z}}_i,$$

(38)

where $c_j=1$ or 0.

For those destabilizers $g_k^{\prime }$ that satisfy

$$\{{\widehat{g}}_k^{\prime },\pm {\widehat{Z}}_i\}=0,$$

(39)

$c_k=1$. Otherwise, $c_k=0$. This can be seen from

$$\{{\widehat{g}}_k^{\prime },\pm {\widehat{Z}}_i\}=\{{\widehat{g}}_k^{\prime },\prod _{j=1}^n{\widehat{g}}_j^{c_j}\}=\prod _{\begin{array}{c}j=1\\ j\ne k\end{array}}^n{\widehat{g}}_j^{c_j}\{{\widehat{g}}_k^{\prime },{\widehat{g}}_k^{c_k}\}=0,$$

(40)

where we used part (ii) of Definition 2 of the destabilizers and Equation (39). 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 ${\widehat{Z}}_i$ must be multiplied together. Every row $(n+j)$ in the bottom half of the tableau, such that $x_{ji}=1$ (for $j\in \{1,\dots ,n\}$), can be added up together and stored in a temporary register. The resultant phase $\pm 1$ of this sum is the measurement result we are looking for.

Checking if each destabilizer commutes or anti-commutes with ${\widehat{Z}}_i$ takes a constant number of operations. One multiplication takes $O\left(n\right)$ operations, and there are $O\left(n\right)$ multiplications needed. Therefore, a measurement takes $O\left(n^2\right)$ operations overall.

5. Discussion

As we made clear throughout Section 4, the scaling of the number of required operations with respect to number of qudits n is exactly the same in the ($d=2$) Aaronson–Gottesman algorithm as in the (odd d) Wigner algorithm presented in Section 3. 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 ${\mathsf{\Phi }}_t$ and vector ${\mathit{r}}_t$ we use. In particular, the tableau is equal to $\left(\begin{array}{cc}{\mathsf{\Phi }}_t& {\mathit{r}}_t\end{array}\right)$:

$${\mathsf{\Phi }}_t=\left(\begin{array}{cc}\frac{\partial {\mathit{p}}_0}{\partial {\mathit{p}}_t}& \frac{\partial {\mathit{p}}_0}{\partial {\mathit{q}}_t}\\ \frac{\partial {\mathit{q}}_0}{\partial {\mathit{p}}_t}& \frac{\partial {\mathit{q}}_0}{\partial {\mathit{q}}_t}\end{array}\right)\equiv \left(\begin{array}{cccccc}{x}_{11}& \dots & x_{1n}& z_{11}& \dots & z_{1n}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ x_{n1}& \dots & x_{nn}& z_{n1}& \dots & z_{nn}\\ x_{(n+1)1}& \dots & x_{(n+1)n}& z_{(n+1)1}& \dots & z_{(n+1)n}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ x_{\left(2n\right)1}& \dots & x_{\left(2n\right)n}& z_{\left(2n\right)1}& \dots & z_{\left(2n\right)n}\end{array}\right)$$

(41)

and

$${\mathit{r}}_t=\left(\begin{array}{c}{\mathit{r}}_p\\ {\mathit{r}}_q\end{array}\right)\equiv \left(\begin{array}{c}{r}_1\\ ⋮\\ r_n\\ r_{n+1}\\ ⋮\\ r_{2n}\end{array}\right).$$

(42)

This can be seen through the following equation:

$$\begin{array}{ccc}\hfill \exp \left(\frac{2\pi i}{d}\sum _{j=1}^{2n}{{\mathsf{\Phi }}_t}_{n+i,j}{\widehat{x}}_j\right)\left|{\mathsf{\Psi }}_t\right.〉& =& \prod _{j=1}^{2n}exp\left(\frac{2\pi i}{d}{{\mathsf{\Phi }}_t}_{n+i,j}{\widehat{x}}_j\right)\left|{\mathsf{\Psi }}_t\right.〉\hfill \\ & =& exp\left(\frac{2\pi i}{d}{r_t}_i\right)\left|{\mathsf{\Psi }}_t\right.〉,\hfill \end{array}$$

(43)

where $\widehat{\mathit{x}}\equiv \left(\widehat{\mathit{p}},\widehat{\mathit{q}}\right)$. Multiplying the right-hand side of the first equation and the second equation by $exp\left(-\frac{2\pi i}{d}{r_t}_i\right)$, it follows that

$$exp\left(-\frac{2\pi i}{d}{r_t}_i\right)\prod _{j=1}^{2n}exp\left(\frac{2\pi i}{d}{{\mathsf{\Phi }}_t}_{n+i,j}{\widehat{x}}_j\right){\mathsf{\Psi }}_t={\widehat{g}}_i\left|{\mathsf{\Psi }}_t\right.〉=\left|{\mathsf{\Psi }}_t\right.〉.$$

(44)

In other words, ${r_t}_i$ specifies the phase $exp\left(-\frac{2\pi i}{d}{r_t}_i\right)$ of the ith stabilizer, which is itself specified by ${{\mathsf{\Phi }}_t}_{n+i,j}$ for $j\in \{0,\dots ,2n\}$. These are the same roles for $\mathit{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 ${\mathsf{\Phi }}_{n+j,i}$ to determine if a measurement on the ith qudit will be random or not. They also use a very similar protocol to determine the outcome of deterministic measurements. The Wigner-based 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 $\mathsf{\Phi }$, 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. In addition, the latter half of both the Wigner function’s ${\mathit{r}}_t$ and Aaronson–Gottesman’s $\mathit{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 ${\mathit{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 $\widehat{A}=\{{\widehat{P}}_i,{\widehat{F}}_i\}\forall i=\{1,\dots ,n\}$, the Aaronson–Gottesman algorithm specifies that for a q- or p-state, its Wigner function evolves by:

$$W_{\mathsf{\Psi }}({\mathcal{M}}_{\widehat{A}}\mathit{x}).$$

(45)

However, for $\left|r\right.〉=\frac{1}{\sqrt{2}}(\left|0\right.〉\pm i\left|1\right.〉)$, a Y-state which is diagonal in the $pq$ plane, its Wigner function must first be translated:

$$W_{\mathsf{\Psi }}\left({\mathcal{M}}_{\widehat{A}}\mathit{x}+\beta \right),$$

(46)

where the translation $\beta$ can be $(1,0)$ or $(0,1)$ equivalently. There is a similar state-dependence for the two-qubit CNOT gate ${\widehat{C}}_{ij}$.

This demonstrates that the Aaronson–Gottesman algorithm is state-dependent on the qubit stabilizer state it is acting on. On the other hand, the Wigner function algorithm on odd d 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 Equations (4) and (5) that underlie the discrete Wigner formulation, are a subgroup of $U\left(d\right)$ instead of $SU\left(d\right)$ for $d=2$ [15]. Furthermore, qubits exhibit state-independent contextuality while odd d qudits do not [16]. Recent progress on this subject relating non-contextuality to classical simulatability for qubits can be found here [17,18].

6. Example of Stabilizer Evolution

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 $\left|0\right.〉\otimes \left|0\right.〉$. To illustrate this process we decompose the two qutrit Wigner function of this state into nine $3\times 3$ grids, as shown in Figure 2. The prepared Wigner function is denoted in Figure 3 with the color black, and the Wigner function represented by setting ${\mathsf{\Phi }}_0=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)$ (i.e., considering the top n rows of $\mathsf{\Phi }$ to be a separate Wigner function, as discussed at the end of Section 3.1) is denoted with the color gray.

We begin with

$$\begin{array}{c}{W}_{\mathsf{\Psi }}\left(\mathit{x}\right)=\hfill \\ {\delta }_{\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}{p_t}_1\\ {p_t}_2\\ {q_t}_1\\ {q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)}={\delta }_{\left(\begin{array}{c}{p_t}_1\\ {p_t}_2\\ {q_t}_1\\ {q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)},\hfill \end{array}$$

(47)

denoting an initially prepared state of $\left|0\right.〉\otimes \left|0\right.〉$. This is clear in Figure 3a by the black band that lies along all Weyl phase space points with ${q_t}_1=0$ and ${q_t}_2=0$. On the other hand, the gray manifold is perpendicular to the black one, and lies along Weyl phase space points with ${p_t}_1=0$ and ${p_t}_2=0$.

Acting on this state with ${\widehat{F}}_1$ produces $\frac{1}{\sqrt{3}}\left(e^{\frac{2\pi i}{3}0\times 0}\left|0\right.〉+e^{\frac{2\pi i}{3}1\times 0}\left|1\right.〉+e^{\frac{2\pi i}{3}2\times 0}\left|2\right.〉\right)\otimes \left|0\right.〉$. Applying the algorithm specified at the end of Section 3.2, we find:

$$\begin{array}{c}{W}_{\mathsf{\Psi }}\left(\mathit{x}\right)=\hfill \\ {\delta }_{\left(\begin{array}{cccc}0& 0& -1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}{p_t}_1\\ {p_t}_2\\ {q_t}_1\\ {q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)}={\delta }_{\left(\begin{array}{c}-{q_t}_1\\ {p_t}_2\\ {p_t}_1\\ {q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right).}\hfill \end{array}$$

(48)

Thus, the momentum of qutrit 1 is now determined and is 0 while the second qutrit is unchanged. This can be seen in Figure 3b, where the ${q_t}_2$ values of the non-zero Weyl phase space points are the same, while the state has rotated by $-\pi /2$ in $({p_t}_1,{q_t}_1)$-space. A similar transformation has occurred for the perpendicular gray manifold.

Acting next with ${\widehat{C}}_{12}$ produces the Bell state $\frac{1}{\sqrt{3}}\left(\left|00\right.〉+\left|11\right.〉+\left|22\right.〉\right)$, which is represented by the following Wigner function:

$$\begin{array}{c}{W}_{\mathsf{\Psi }}\left(\mathit{x}\right)=\hfill \\ {\delta }_{\left(\begin{array}{cccc}0& 0& -1& 0\\ 0& 1& 0& 0\\ 1& 1& 0& 0\\ 0& 0& -1& 1\end{array}\right)\left(\begin{array}{c}{p_t}_1\\ {p_t}_2\\ {q_t}_1\\ {q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right)}={\delta }_{\left(\begin{array}{c}-{q_t}_1\\ {p_t}_2\\ {p_t}_1+{p_t}_2\\ -{q_t}_1+{q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right).}\hfill \end{array}$$

(49)

The entanglement between the two qutrits is evident in both of their dependence on each other’s momenta and positions, ${p_t}_1=-{p_t}_2$ and ${q_t}_1={q_t}_2$, specified by the last two rows. Figure 3c shows that the state is still representable as lines in Weyl phase space, except they now traverse through the different planes of $({q_t}_1,{p_t}_1)$ associated with each value of $({q_t}_2,{p_t}_2)$. However, if you consider the left column in Figure 3c corresponding to ${q_t}_2=0$, you can see that the only black Weyl phase points are at ${q_t}_1=0$. Similarly, the middle column corresponding to ${q_t}_2=1$ shows that ${q_t}_1=1$, and the right column corresponding to ${q_t}_2=2$ shows that ${q_t}_1=2$ too, confirming that $\left|\mathsf{\Phi }\right.〉=\frac{1}{\sqrt{3}}\left(\left|00\right.〉+\left|11\right.〉+\left|22\right.〉\right)$. Thus, the entanglement of the two qutrits’ positions is clearly evident in this Figure of the Wigner function.

We then proceed to measure qutrit 1. Since the lower two equations involve ${p_t}_1$, 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 qutrit 2 into the same state:

$$\begin{array}{c}{W}_{\mathsf{\Psi }}\left(\mathit{x}\right)=\hfill \\ {\delta }_{\left(\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& -1& 1\end{array}\right)\left(\begin{array}{c}{p_t}_1\\ {p_t}_2\\ {q_t}_1\\ {q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right)}={\delta }_{\left(\begin{array}{c}{p_t}_1+{p_t}_2\\ {p_t}_2\\ {q_t}_1\\ -{q_t}_1+{q_t}_2\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right).}\hfill \end{array}$$

(50)

The lower two rows show that now ${q_t}_1=1$, as we chose, and ${q_t}_2={q_t}_1=1$. The collapse of qutrit 2 into $\left|1\right.〉$ can also been seen in Figure 3c by the fact that ${q_t}_1=1$ only in the $3\times 3$ grids that correspond to ${q_t}_2=1$ too.

Finally, the fact that a measurement of ${q_t}_2$ would be deterministic at this point can be seen in the fact that ${p_t}_2$ is not present in the last two rows of ${\mathsf{\Phi }}_t$. Furthermore, it is clear, since the first row has a coefficient of 1 in front of ${p_t}_1$, that the corresponding third row must be added with weight 1 to the fourth row to obtain this deterministic measurement outcome of ${q_t}_2=1$. This can also be seen in Figure 3 by finding the projection of ${p_0}_1$ onto ${p_t}_2$, 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 $\pi /4$ diagonally with respect to each other have a projection equal to 2 in this discrete geometry.)

7. Conclusions

In summary, we introduced an algorithm that efficiently simulates stabilizer state evolution under Clifford gates and measurements in the $\widehat{Z}$ 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 Equations (43) and (44), to the well-known Aaronson–Gottesman tableau algorithm [2] for qubits, 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 stabilizer propagation. Equation (44) 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 difference appears to be the state-dependence of this evolution, and likely is related to the state-independent contextuality that qubits exhibit, which odd d qudits do not. Exploring the details of this state-dependence is a promising subject of future study.