Discrete Wigner Function Derivation of the Aaronson-Gottesman Tableau Algorithm

TL;DR

This paper introduces a discrete Wigner function-based simulation algorithm for odd-dimensional qudits that matches the efficiency of the Aaronson-Gottesman algorithm, leveraging symplectic structure in phase space.

Contribution

It derives a new simulation algorithm using discrete Wigner functions for odd qudits, aligning with Aaronson-Gottesman in efficiency and offering insights for qubit extension.

Findings

The Wigner function algorithm has the same complexity as Aaronson-Gottesman.

Efficiency is due to harmonic evolution in symplectic phase space.

Differences are linked to group structure and contextuality in qubits.

Abstract

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.

Tables

Table 1. Table 1: Transformation of stabilizer generators under Clifford operations.

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

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

$x_{ij}$	$z_{ij}$	${\hat{P}}_j$
$0$	$0$	${\hat{I}}_j$
$0$	$1$	${\hat{Z}}_j$
$1$	$0$	${\hat{X}}_j$
$1$	$1$	${\hat{Y}}_j$

Equations

$${\hat{Z}}_j|k_1,\dots ,k_j,\dots ,k_n⟩=e^{\frac{2\pi i}{d}{k}_j}|k_1,\dots ,k_j,\dots ,k_n⟩.$$

$${\hat{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}|k_1,\dots ,k_j,\dots ,k_n⟩⟨l_1,\dots ,l_j,\dots ,{\text{l}}_n|,$$

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

$${\hat{X}}_j^{\delta q}|k_1,\dots ,k_j,\dots ,k_n⟩\equiv |k_1,\dots ,k_j\oplus \delta q,\dots ,k_n⟩,$$

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

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

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

$${\Psi }_x(\mathbf{p},\mathbf{q})$$

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

$${\Psi }_x\left({\mathcal{M}}_{\hat{\mathbf{V}}}\left(\mathbf{x}+{\mathbf{\alpha }}_{\hat{V}}/2\right)+{\mathbf{\alpha }}_{\hat{V}}/2\right),$$

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

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

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

$${\left({\mathcal{M}}_{{\hat{F}}_i}\right)}_{j,k}={\delta }_{j,k}-{\delta }_{i,j}{\delta }_{i,k}-{\delta }_{n+i,j}{\delta }_{n+i,k}$$

$$+{\delta }_{i,j}{\delta }_{n+i,k}-{\delta }_{n+i,j}{\delta }_{i,k},$$

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

$$\left(\begin{array}[]{c}\boldsymbol{p}_{0}\\ \boldsymbol{q}_{0}\end{array}\right)=\boldsymbol{\Phi}_{t}\left(\begin{array}[]{c}\boldsymbol{p}_{t}\\ \boldsymbol{q}_{t}\end{array}\right),$$

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

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

$$\boldsymbol{\mathcal{J}}=\left(\begin{array}[]{cc}0&-\mathbb{I}_{n}\\ \mathbb{I}_{n}&0\end{array}\right),$$

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

$${\left({\mathcal{M}}_{{\hat{F}}_i}^{-1}\right)}_{j,k}={\delta }_{j,k}-{\delta }_{i,j}{\delta }_{i,k}-{\delta }_{n+i,j}{\delta }_{n+i,k}$$

$$-{\delta }_{i,j}{\delta }_{n+i,k}+{\delta }_{n+i,j}{\delta }_{i,k},$$

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

$$\boldsymbol{\Phi}_{0}=\left(\begin{array}[]{cc}\mathbb{I}_{n}&\boldsymbol{0}\\ \boldsymbol{0}&\mathbb{I}_{n}\end{array}\right),$$

$$\boldsymbol{\Phi}_{t}\left(\begin{array}[]{c}{p_{t}}_{1}\\ \vdots\\ {p_{t}}_{i}\\ \vdots\\ {p_{t}}_{n}\\ \boldsymbol{q}_{t}\end{array}\right)=\left(\begin{array}[]{c}\boldsymbol{r_{t}}_{p}\\ \boldsymbol{r_{t}}_{q}\end{array}\right)\underset{\hat{Z}_{i}}{\mapsto}\boldsymbol{\Phi}_{t}\left(\begin{array}[]{c}{p_{t}}_{1}\\ \vdots\\ {p_{t}}_{i}+1\\ \vdots\\ {p_{t}}_{n}\\ \boldsymbol{q}_{t}\end{array}\right)=\left(\begin{array}[]{c}\boldsymbol{r_{t}}_{p}\\ \boldsymbol{r_{t}}_{q}\end{array}\right).$$

$${\mathbf{\Phi }}_t{\mathbf{x}}_t$$

$$\displaystyle\left(\begin{array}[]{c}\boldsymbol{p}_{0}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})\\ \boldsymbol{q}_{0}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})\end{array}\right)$$

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

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

$$\left(\begin{array}[]{c}\boldsymbol{p}_{t}\\ \boldsymbol{q}_{t}\end{array}\right)=\boldsymbol{\mathcal{M}}_{t}\left(\begin{array}[]{c}\boldsymbol{p}_{0}\\ \boldsymbol{q}_{0}\end{array}\right).$$

Full text

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 Gottesman98 and are closed under action by Clifford gates. Efficient simulation of such systems was demonstrated with the tableau algorithm of Aaronson and Gottesman Aaronson04 ; Gottesman98 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 Gottesman99 ; Mari12 ; Howard14 .

Recent progress has been the result of work by Wootters Wootters87 , Eisert Mari12 , Gross Gross06 , and Emerson Howard14 , 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 Gross06 and Clifford group gates have been shown to have underlying harmonic Hamiltonians that preserve the discrete Weyl phase space points Kocia16 . This means Clifford circuits are expressible by path integrals truncated at order $\hbar^{0}$ and are thus manifestly classical Kocia16 ; Dax17 .

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 $\hat{Z}$ Pauli basis for odd $d$. We then show its equivalence to the well-known Aaronson-Gottesman tableau algorithm Aaronson04 . Both algorithms require $\mathcal{O}(n^{2})$ qudits to represent $n$ stabilizer states, $\mathcal{O}(n)$ operations per Clifford operator, and deterministic and random measurements require $\mathcal{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 $\hat{Z}_{j}$ operator on the $j$th qubit for $n$ qubits acts as a “boost” operator:

[TABLE]

The discrete Fourier transform operator is defined by:

[TABLE]

This is equivalent to the Hadamard gate and allows us to define the Pauli $\hat{X}_{j}$ operator as follows:

[TABLE]

While $\hat{Z}_{j}$ is a boost, $\hat{X}_{j}$ is a shift operator because

[TABLE]

where $\oplus$ denotes integer addition mod $d$.

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

[TABLE]

and

[TABLE]

Thus, we can refer to the $\hat{X}_{j}$ basis as the momentum ($p_{j}$) basis, which is equivalent to the Fourier transform of the $q_{j}$ basis:

[TABLE]

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

The Wigner function $\Psi_{x}(\boldsymbol{p},\boldsymbol{q})$ of a pure state $\left|\Psi\right\rangle$ is defined on this discrete Weyl phase space:

[TABLE]

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.

II.1 Clifford Gates

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

[TABLE]

Wigner functions $\Psi_{x}(\boldsymbol{x})$ of states evolve under Clifford operators $\hat{V}$ by

[TABLE]

where $\boldsymbol{x}\equiv(\boldsymbol{p},\boldsymbol{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 gate $\hat{P}_{i}$, the Hadamard gate $\hat{F}_{i}$, and the controlled-not $\hat{C}_{ij}$ (which act on the $i$th and $j$th qudits).

The phase shift $\hat{P}_{i}$ is a one-qudit gate with the underlying Hamiltonian $H_{\hat{P}_{i}}=-\frac{d+1}{2}q_{i}^{2}+\frac{d+1}{2}q_{i}$ Kocia16 . Without loss of generality, we will instead consider

[TABLE]

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:

[TABLE]

where $[\hat{P}_{i},\hat{Z}_{i}]=[\hat{P}^{\prime}_{i},\hat{Z}_{i}]=0$. Hence, $\hat{P}^{\prime}_{i}$ is an adequate replacement generator for $\hat{P}_{i}$, and we will use it instead of $\hat{P}_{i}$ from now on. Since its Hamiltonian has no linear term ($H_{\hat{P}^{\prime}_{i}}=-q_{i}^{2}$), this leads to an easier presentation ahead since $\boldsymbol{\alpha}_{\hat{P}^{\prime}_{i}}=\boldsymbol{0}$.

The corresponding equations of motion for $\hat{P}^{\prime}_{i}$ are $\dot{p}_{i}=2q_{i}$ and $\dot{q}_{i}=0$. Hence, for $\Delta t=1$,

[TABLE]

The Hadamard gate $\hat{F}_{i}$ is a one-qudit gate and has the underlying Hamiltonian $H_{\hat{F}_{i}}=-\frac{\pi}{4}(p_{i}^{2}+q_{i}^{2})$ Kocia16 . 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$,

[TABLE]

and $\boldsymbol{\alpha}_{\hat{F}_{i}}=\boldsymbol{0}$.

Finally, the two-qudit CNOT $\hat{C}_{ij}$ on control qudit $i$ and second qudit $j$ has the corresponding Hamiltonian $H_{\hat{C}_{ij}}=p_{i}q_{j}$ Kocia16 . 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$,

[TABLE]

and $\boldsymbol{\alpha}_{\hat{C}_{ij}}=\boldsymbol{0}$.

II.2 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 Kocia16 :

Theorem 1**.**

The discrete Wigner function $\Psi_{x}(\boldsymbol{x})$ of a stabilizer state $\Psi$ for any odd $d$ and $n$ qudits is $\delta_{\boldsymbol{\Phi}\times\boldsymbol{x},\boldsymbol{r}}$ for $2n\times 2n$ matrix $\boldsymbol{\Phi}$ and $2n$ vector $\boldsymbol{r}$ with entries in $\mathbb{Z}/d\mathbb{Z}$.

In particular, if we begin with a stabilizer state defined as $\left|\Psi\right\rangle_{0}=\left|\boldsymbol{q}_{0}\right\rangle$, then it was shown in Kocia16 that ${\Psi_{x}}(\boldsymbol{x})=\delta_{\boldsymbol{\Phi}_{0}\times\boldsymbol{x},\boldsymbol{r}_{0}}$, where \boldsymbol{\Phi}_{0}=\left(\begin{array}[]{cc}0&0\\ 0&\mathbb{I}_{n}\end{array}\right), and $\boldsymbol{r}_{0}=(\boldsymbol{0},\boldsymbol{q}_{0})$.

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.

III.1 Stabilizer Representation

From Theorem 1, propagation of the stabilizer state $\Psi$ can be represented by considering the state in the Wigner representation $\Psi_{x}(\boldsymbol{x})=\delta_{\boldsymbol{\Phi}_{t}\cdot\boldsymbol{x},\boldsymbol{r}_{t}}$. In this way, $\boldsymbol{\Phi}_{t}$ and $\boldsymbol{r}_{t}$ specify a linear system of equations in terms of $\boldsymbol{p}_{t}$ and $\boldsymbol{q}_{t}$. The first $n$ rows of $\boldsymbol{\Phi}_{t}$ are the coefficients of $(\boldsymbol{p}_{t},\boldsymbol{q}_{t})^{T}$ in $\boldsymbol{p}_{0}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$ and the last $n$ rows of $\boldsymbol{\Phi}_{t}$ are the coefficients of $(\boldsymbol{p}_{t},\boldsymbol{q}_{t})^{T}$ in $\boldsymbol{q}_{0}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$:

[TABLE]

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

Since the Wigner functions $\Psi_{x}(\boldsymbol{x})$ of stabilizer states propagate under $\boldsymbol{\mathcal{M}}$ as $\Psi_{x}(\boldsymbol{\mathcal{M}\boldsymbol{x}})$, it follows that

[TABLE]

The importance of vector $\boldsymbol{r}_{t}$ and when it must be updated will become evident when we consider random measurements. Hence, after $n$ operations $\boldsymbol{\mathcal{M}}_{1}$, $\boldsymbol{\mathcal{M}}_{2}$, $\ldots$, $\boldsymbol{\mathcal{M}}_{n}$,

[TABLE]

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

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

[TABLE]

and $\mathbb{I}_{n}$ is the $n\times n$ identity matrix.

Thus the the stability matrices $\boldsymbol{\mathcal{M}}$ for $\hat{F}_{i}$, $\hat{P}^{\prime}_{i}$ and $\hat{C}_{ij}$ given in Eqs. 13-15 differ from their inverses only by sign changes in their off-diagonal elements:

[TABLE]

[TABLE]

and

[TABLE]

We assume the quantum state is initialized in the computational basis state $\Psi_{0}=\underbrace{\left|0\right\rangle\otimes\cdots\otimes\left|0\right\rangle}_{n}$ and so initially we should set \boldsymbol{\Phi}_{0}=\left(\begin{array}[]{cc}\boldsymbol{0}&\boldsymbol{0}\\ \boldsymbol{0}&\mathbb{I}_{n}\end{array}\right) and $\boldsymbol{r}_{0}=\boldsymbol{0}$. The initial stabilizer state is $\boldsymbol{\Psi_{x}}_{0}=\delta_{\boldsymbol{q}_{t},0}$. However, it will become clear when we discuss measurements that it is practically useful to instead set

[TABLE]

thereby setting ${\Psi_{x}}_{0}=\delta_{(\boldsymbol{p}_{t},\boldsymbol{q}_{t}),(\boldsymbol{0},\boldsymbol{0})}$—not a true Wigner function. This is equivalent to the last matrix if the first $n$ rows in $\boldsymbol{\Phi}_{t}\boldsymbol{x}$ and $\boldsymbol{r}_{t}$ are ignored—the same as ignoring $\boldsymbol{p}_{0}(\boldsymbol{p}_{t},\boldsymbol{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}(n^{2})$ dits of storage to track $\boldsymbol{\Phi}_{t}$ and $\boldsymbol{r}_{t}$. More precisely, since $\boldsymbol{\Phi}_{t}$ is a $2n\times 2n$ matrix and $\boldsymbol{r}_{t}$ is an $2n$-vector, $2n(2n+1)$ dits of storage are necessary.

III.2 Unitary Propagation

$\boldsymbol{\Phi}_{t}$ contains the coefficients of the linear equations relating $\boldsymbol{x}_{0}$ to $\boldsymbol{x}_{t}$. Each row is one equation relating ${q_{0}}_{i}$ or ${p_{0}}_{i}$ to $\boldsymbol{x}_{t}$. When manipulating rows of $\boldsymbol{\Phi}_{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 gates $\hat{F}_{i}$, $\hat{P}_{i}$ and $\hat{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 $\boldsymbol{\Phi}_{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 $\mathcal{O}(n)$ operations. Therefore, keeping track of propagation of stabilizer states by Clifford gates can be simulated with $\mathcal{O}(n)$ operations.

Let us examine these unitary operations more closely. Defining $\oplus$ and $\ominus$ to be mod $d$ addition and subtraction respectively, we find:

Phase gate on qudit $i$ ($\hat{P}^{\prime}_{i}$). For all $j\in\{1,\,\ldots,\,2n\}$, set $\boldsymbol{\Phi}_{j,n+i}\mapsto\boldsymbol{\Phi}_{j,n+i}\ominus 2\boldsymbol{\Phi}_{j,n}$.

Hadamard gate on qudit $i$ ($\hat{F}_{i}$). Negate $\boldsymbol{\Phi}_{j,i}$ mod $d$, and then swap $\ominus\boldsymbol{\Phi}_{j,i}$ and $\boldsymbol{\Phi}_{n+i,j}$.

CNOT from control $i$ to target $j$ ($\hat{C}_{ij}$). Set $\boldsymbol{\Phi}_{k,j}\mapsto\boldsymbol{\Phi}_{k,j}\oplus\boldsymbol{\Phi}_{k,i}$ and $\boldsymbol{\Phi}_{k,n+i}\mapsto\boldsymbol{\Phi}_{k,n+i}\ominus\boldsymbol{\Phi}_{k,n+j}$.

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

III.3 Measurement

The outcome of a measurement $\hat{Z}_{i}$ on a stabilizer state can be either random or deterministic. As described above, the bottom half of $\boldsymbol{\Phi}_{t}$ defines ${q_{0}}_{j}$ for $j\in\{1,\,\ldots,\,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 $\boldsymbol{\Phi}_{t}$ give the coefficient of ${p_{t}}_{i}$ and ${q_{t}}_{i}$ in ${q_{0}}_{j}$ for $j\in\{1,\,\ldots,\,n\}$. If the coefficient of ${p_{t}}_{i}$ in any ${q_{0}}_{j}$ is non-zero then the measurement $\hat{Z}_{i}$ will be random. If all coefficients of ${p_{t}}_{i}$ are zero for ${q_{0}}_{j}$ $\forall j$, then the measurement of $\hat{Z}_{i}$ will be deterministic. This can be seen from the fact that if our stabilizer state $\left|\Psi\right\rangle$ is an eigenstate of $\hat{Z}_{i}$ then $\hat{Z}_{i}\left|\Psi\right\rangle=e^{i\phi}\left|\Psi\right\rangle$ for some $\phi\in\mathbb{R}$ and (discrete) Wigner functions do not change under a global phase. Thus, measuring $\hat{Z}_{i}$ leaves the Wigner function of $\left|\Psi\right\rangle$ invariant if the measurement is deterministic. Since $\hat{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:

[TABLE]

Thus, if the lower half of the $i$th column of $\boldsymbol{\Phi}_{t}$ is zero then $\hat{Z}_{i}$ leaves the Wigner function invariant (and so the measurement is deterministic). Verifying that these coefficients are all zero takes $\mathcal{O}(n)$ operations for each $\hat{Z}_{i}$.

In other words, to see if a given measurement of $\hat{Z}_{i}$ is random or deterministic, a search must be performed for non-zero ${\boldsymbol{\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 ${\boldsymbol{\Phi}_{t}}_{n+j,i}$ element exists, then the measurement $\hat{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 $\boldsymbol{\Phi}_{t}$ have a non-zero entry in the $i$th column, ${\boldsymbol{\Phi}_{t}}_{n+j,i}$. Since the random measurement $\hat{Z}_{i}$ will project qudit $i$ onto a position state, we will replace the $(n+j)$th row with ${q_{0}}_{i}=q^{\prime}_{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_{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^{\prime}_{i}$.

There is one more thing to do, which will be important for deterministic measurements: replace the $j$th row with the old $(n+j)$th row. This sets ${p_{0}}_{i}={q_{0}}_{j}(\boldsymbol{p}_{t},\boldsymbol{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 $\boldsymbol{\Phi}=\mathbb{I}_{2n}$. However, now we also conserve any dependence ${p_{0}}_{i}$ has on the other qudits incurred during unitary propagation. In other words, we conserve ${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 $\boldsymbol{\Phi}_{t}$ and $\boldsymbol{r}_{t}$ with a randomly chosen measurement outcome $q^{\prime}_{i}$ (i.e. ${q_{0}}_{i}=q^{\prime}_{i}$), the identification of rows $(n+i)$ and $(n+j)$ are exchanged, so that the former now specifies ${q_{0}}_{j}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$ while the latter specifies ${q_{0}}_{i}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$, respectively. ${p_{0}}_{i}$ has also been updated by replacing the $j$th row in the first half of $\boldsymbol{\Phi}_{t}$, with the $(n+j)$th row we just changed. Again, this row now describes ${p_{0}}_{i}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$ while the $i$th row now specifies ${p_{0}}_{j}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$. Overall, this takes $\mathcal{O}(n^{2})$ operations since we are replacing $\mathcal{O}(n)$ rows with $\mathcal{O}(n)$ entries.

Case 2: Deterministic Measurement

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

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

[TABLE]

where we are interested in linear combinations of the bottom half, $\boldsymbol{q}_{0}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$, to solve for the measurement outcome ${q_{t}}_{i}$:

[TABLE]

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

Lemma 1**.**

The coefficient in front of ${p_{t}}_{i}$ in the row of $\boldsymbol{\Phi}_{t}$ that specifies ${p_{0}}_{j}(\boldsymbol{p}_{t},\boldsymbol{q}_{t})$, ${\boldsymbol{\Phi}_{t}}_{ji}$, is equal to the coefficient $c_{ij}$ in front of ${q_{0}}_{j}$ that makes up ${q_{t}}_{i}$ in Eq. 31. Or, equivalently,

[TABLE]

Proof.

Under evolution under the Clifford group operators,

[TABLE]

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

[TABLE]

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

[TABLE]

This property can also be seen in the drawing of phase space shown in Fig. 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}(\boldsymbol{p}_{0},\boldsymbol{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}(\boldsymbol{p}_{t},\boldsymbol{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

[TABLE]

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}(n^{2})$ operations since Eq. 32 is a sum of $\mathcal{O}(n)$ vectors made up of $\mathcal{O}(n)$ components. So, the overall measurement protocol takes $\mathcal{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. 51) to be performed efficiently.

IV Aaronson-Gottesman Tableau Algorithm

The Aaronson-Gottesman tableau algorithm was originally defined for qubits ($d=2$) Aaronson04 . Like the algorithm we presented in the previous section, it only requires overall $\mathcal{O}(n^{2})$ operations for propagationan and measurement for $n$ qubits. The algorithm has been proven to be extendable to $d>2$ deBeaudrap11 and similar algorithms have been formulated in $d>2$ Yoder12 . Alternatives have also been developed to the tableau formalism, though they prove to be equally efficient in worst-case scenarios Anders06 . 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.

IV.1 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=\{\hat{g}\in\mathcal{P}$ such that $\hat{g}\left|\psi\right\rangle=\left|\psi\right\rangle\}$ are called the stabilizers of state $\left|\psi\right\rangle$, where $\mathcal{P}$ is the set of Pauli operators, each of which has the form $e^{\frac{\pi i}{2}\alpha}\hat{P}_{1}\cdots\hat{P}_{n}$ where $\alpha\in\{0,1,2,3\}$ for $n$ qubits with $\hat{P}_{i}\in\{\hat{I}_{i},\hat{Z}_{i},\hat{X}_{i},\hat{Y}_{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 Gottesman98 ; Aaronson04 .

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, $\{\hat{g}_{1},\hat{g}_{2},\ldots,\hat{g}_{n}\}$, is in the Pauli group, and each generator has the form:

[TABLE]

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

IV.2 Unitary Propagation

For each Clifford operation, Aaronson and Gottesman showed that only $\mathcal{O}(n)$ operations are necessary to update all generators Aaronson04 . Specifically, according to the update rules in Table 1, 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 $\{\hat{g}^{\prime}_{1},\ldots,\hat{g}^{\prime}_{n}\}$ are the operators that generate the full Pauli group with the stabilizers $\{\hat{g}_{1},\ldots,\hat{g}_{n}\}$. They have the following properties:

• (i)

$\hat{g}^{\prime}_{1}$, $\hat{g}^{\prime}_{2}$, $\ldots$, $\hat{g}^{\prime}_{n}$ commute.

• (ii)

Each destabilizer $\hat{g}^{\prime}_{h}$ anti-commutes with the corresponding stabilizer $\hat{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 Aaronson04 .

Aaronson-Gottesman defined such a $2n\times(2n+1)$ binary tableau matrix as:

[TABLE]

This matrix contains $2n$ rows. The first n rows denote the destabilizers $\hat{g}^{\prime}_{1}$ to $\hat{g}^{\prime}_{n}$ while rows $(n+1)$ to $2n$ represent the stabilizers $\hat{g}_{1}$ to $\hat{g}_{n}$. The $(n+1)$th bit in each row denotes the phase $(-1)^{r_{i}}$ for each generator. We encode the $j$th Pauli operator in the $i$th row as shown in Table 2.

We can update the stabilizers and destabilizers as follows:

Hadamard gate on qubit $i$ For all $j\in\{1,2,...,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,...,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,...,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 $\boldsymbol{\Phi}$ transformation rules in Section IV. The most glaring difference is that the Aaronson-Gottesman algorithm involves updates of the vector $\boldsymbol{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.

IV.3 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 Aaronson04 :

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,\,\ldots,\,n$. To calculate the resultant phase, Aaronson and Gottesman first defined the following function:

[TABLE]

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

[TABLE]

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

Case 1: Random Measurement

$\hat{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, $\hat{g}_{j}$, and update the rest by replacing them with their product with $\hat{g}_{j}$ (i.e. taking the rowsum of their corresponding rows) such that they commute with $\hat{Z}_{i}$. These updates take $\mathcal{O}(n^{2})$ operations. Finally, we only need to replace $\hat{g}_{j}$ by $\hat{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\neq j$ with the sum of the $j$th and $k$th row (using the rowsum function) sets all $x_{ki}=0$ for $k\neq j$.

Finally, we replace the $(j-n)$th row with the $j$th row and update the $j$th row by setting $z_{ji}=1$ and all other $x_{jk}$s and $z_{jk}$s to [math] 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

$\hat{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. But we do need to do some work to retrieve the measurement outcome.

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

[TABLE]

where $c_{j}=1$ or [math].

For those destabilizers $g^{\prime}_{k}$ that satisfy

[TABLE]

$c_{k}=1$. Otherwise $c_{k}=0$. This can be seen from

[TABLE]

where we used part (ii) of Definition 2 of the destabilizers and Eq. 58. 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 $\hat{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,\,\ldots,\,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 $\hat{Z}_{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 $\boldsymbol{\Phi}_{t}$ and vector $\boldsymbol{r}_{t}$ we use. In particular, the tableau is equal to \left(\begin{array}[]{c|c}\boldsymbol{\Phi}_{t}&\boldsymbol{r}_{t}\end{array}\right):

[TABLE]

and

[TABLE]

This can be seen through the following equation:

[TABLE]

where $\hat{\boldsymbol{x}}\equiv\left(\hat{\boldsymbol{p}},\hat{\boldsymbol{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

[TABLE]

In other words, ${r_{t}}_{i}$ specifies the phase $\exp\left(-\frac{2\pi i}{d}{r_{t}}_{i}\right)$ of the $i$th stabilizer, which is itself specified by ${\Phi_{t}}_{n+i,j}$ for $j\in\{0,\ldots,2n\}$. These are the same roles for $\boldsymbol{r}$ and the tableau in the Aaronson-Gottesman tableau algorithm Aaronson04 .

Indeed, both algorithms check the bottom half of their matrices for finite elements of $\boldsymbol{\Phi}_{n+j,i}$ to determine if a measurement on the $i$th qudit will be random or not. They also use a very similar protocal 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 $\boldsymbol{\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. Also, the latter half of both the Wigner function’s $\boldsymbol{r}_{t}$ and Aaronson-Gottesman’s $\boldsymbol{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 $\boldsymbol{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 $\hat{A}=\{\hat{P}_{i},\hat{F}_{i}\}\,\forall i=\{1,\ldots,n\}$, the Aaronson-Gottesman algorithm specifies that for a $q$- or $p$-state, its Wigner function evolves by:

[TABLE]

But for $\left|r\right\rangle=\frac{1}{\sqrt{2}}(\left|0\right\rangle\pm i\left|1\right\rangle)$, a $Y$-state which is diagonal to $q$ and $p$, its Wigner function must first be translated:

[TABLE]

where the translation $\boldsymbol{\beta}$ can be $(1,0)$ or $(0,1)$ equivalently. There is a similar state-dependence for the two-qudit CNOT gate $\hat{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, 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$ Bengtsson17 . Furthermore, qubits exhibit state-independent contextuality while odd $d$ qudits do not Mermin93 .

VI 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\rangle\otimes\left|0\right\rangle$. The prepared Wigner function is illustrated in Fig. 3 with the color black, and the Wigner function represented by setting \boldsymbol{\Phi}_{0}=\left(\begin{array}[]{cc}1&0\\ 0&0\end{array}\right) (i.e. considering the top $n$ rows of our function as a separate Wigner function) is illustrated with the color gray.

We begin with

[TABLE]

denoting an initially prepared state of $\left|0\right\rangle\otimes\left|0\right\rangle$. 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 $\hat{F}_{1}$ produces $\frac{1}{\sqrt{3}}\left(e^{\frac{2\pi i}{3}0\times 0}\left|0\right\rangle+e^{\frac{2\pi i}{3}1\times 0}\left|1\right\rangle+e^{\frac{2\pi i}{3}2\times 0}\left|2\right\rangle\right)\otimes\left|0\right\rangle$. Applying the algorithm specified at the end of Section III.2, we find:

[TABLE]

Thus, the momentum of qudit $1$ is now determined and is [math] while the second qudit is unchanged. This can be seen in Fig. 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 occured for the perpendicular gray manifold.

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

[TABLE]

The entanglement between the two qudits 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. 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_{t}}_{1},{p_{t}}_{1})$ associated with each value of $({q_{t}}_{2},{p_{t}}_{2})$. However, if you consider the left column in Fig. 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|P\right\rangle hi=\frac{1}{\sqrt{3}}\left(\left|00\right\rangle+\left|11\right\rangle+\left|22\right\rangle\right)$. 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_{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 qudit $2$ into the same state:

[TABLE]

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 qudit $2$ into $\left|1\right\rangle$ can also been seen in Fig. 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 the measurement of ${q_{t}}_{2}$ is deterministic can be seen in that ${p_{t}}_{2}$ is not present in the last two rows of $\boldsymbol{\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 Fig. 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 [math], and those that lie $\pi/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 $\hat{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 Eqs. 62 and 63, to the well-known Aaronson-Gottesman tableau algorithm Aaronson04 , 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. Eq. 63 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.