Stationary Phase Method in Discrete Wigner Functions and Classical Simulation of Quantum Circuits

One of the lowest-order corrections to Gaussian quantum mechanics in infinite-dimensional Hilbert spaces are Airy functions: a uniformization of the stationary phase method applied in the path integral perspective. We introduce a"periodized stationary phase method"to discrete Wigner functions of systems with odd prime dimension and show that the $\frac{\pi}{8}$ gate is the discrete analog of the Airy function. We then establish a relationship between the stabilizer rank of states and the number of quadratic Gauss sums necessary in the periodized stationary phase method. This allows us to develop a classical strong simulation of a single qutrit marginal on $t$ qutrit $\frac{\pi}{8}$ gates that are followed by Clifford evolution, and show that this only requires $3^{\frac{t}{2}+1}$ quadratic Gauss sums. This outperforms the best alternative qutrit algorithm (based on Wigner negativity and scaling as $\sim\hspace{-3pt} 3^{0.8 t}$ for $10^{-2}$ precision) for any number of $\frac{\pi}{8}$ gates to full precision.


Introduction
This paper aims to bridge the gap between efficient strong simulation qubit methods, which use low stabilizer rank representations of magic states, and oddprime-dimensional qudit methods, which have relied on Wigner negativity in past studies. We introduce a strong simulation odd-dimensional qudit algorithm that is more efficient than the current state-of-the-art algorithm based on Wigner negativity.
A strong simulation refers to the computation of a quantum probability. This differs from weak simulation, which refers to the computation of sample outcomes of a quantum probability distribution. Strong simulation of quantum universal circuits up to multiplicative error is #P -hard [1][2][3]; the cost of classical strong simulation is expected to not only scale exponentially but to also be harder than tasks that can be performed efficiently on a quantum computer. However, this does not mean that there do not exist efficient methods in particular cases, such as those frequently covered in introductory textbooks on quantum mechanics.
For instance, by representing quantum circuits through the path integral perspective, efficient approximation schemes can be formed by only including lowestorder terms in . It is known that quantum circuits require terms at orders greater than 0 for them to attain quantum universality [4][5][6]. In the continuous-variable setting, this corresponds to requiring non-Gaussianity as a resource. The "minimal" non-Gaussianity that can be added takes the functional form of Airy functions, and these can be treated within the path integral formalism with just one more order of correction in terms of : O( 1 ) [7]. Perhaps a similar systematic expansion in orders of is also possible in the discrete case, and can result in a more efficient treatment of some universal gates.
In the discrete case, non-Gaussianity has been studied in other ways not explicitly dependent on orders of . Such past work on efficient classical simulation has utilized the negativity of Wigner functions or contextuality [8,9] as a resource of non-Gaussianity. More formally, the resource theory of contextuality in finitedimensional systems has also been explicitly considered recently [10,11].
For quantum statesρ, the Wigner function [12,13] ρ(x) is a non-negative function for stabilizer states. Since Clifford gates take stabilizer states to stabilizer states, classical simulation of stabilizer states under Clifford gate evolution using this discrete Wigner formulation can be accomplished in polynomial time and has been shown to be non-contextual [4].
However, Clifford gates and stabilizer states do not allow for universal quantum computation. The extension of this set by π 8 gates or magic states to allow for universal quantum computation introduces contextuality into ρ(x), which means that they can also take negative values [14][15][16][17]. Since the contextual part of a classical simulation is solely responsible for its transition from polynomial to exponential scaling, it is useful to treat contextuality as an operational resource [18] that can be added to a polynomially-efficient backbone in the appropriate amount to allow for a quantum process to be simulated. The goal of such a contextuality injection scheme is a classical computation cost that scales as efficiently as possible with the contextuality present.
Previous attempts at developing a framework for classical algorithms for simulation of qudit quantum circuits that leverage contextuality have done so indirectly by using Wigner function negativity [19]. Recently, there have also been some more direct approaches [11]. In continuous systems, a historically successful approach of leveraging contextuality efficiently is to use semiclassical techniques that rely on expansion in orders of . The relationship between contextuality and higher orders of in the Wigner-Weyl-Moyal (WWM) formalism has been recently established [4][5][6]. Previous derivation of WWM in odd dimensions [20][21][22] were not able to accomplish this because they began by using results from the continuous case of the WWM formalism. In the continuous regime, the stationary phase method is applies in its traditional form, and must then be "discretized" and "periodized". Following that approach does not allow for the derivation of higher order corrections because anharmonic trajectories cannot be obtained by "periodizing" to a discrete Weyl phase space grid [22].
Here we are able to develop a method that includes such higher order terms for discrete odd primedimensional systems. Without appealing to the continuous case, we use a periodized version of Taylor's theorem to develop the stationary phase method in the discrete setting. We explain how this "periodized" form of the stationary phase method differs from the continuous version of the stationary phase method.
This stationary phase method allows us to derive the WWM formalism with asymptotically decreasing order corrections to a subset of unitaries that include the T -gate and π 8 gates. These serve as contextual corrections to non-contextual Clifford operations. We obtain a new method for classical simulation of quantum circuits involving Clifford+ π 8 gates for qutrits. To demonstrate the effectiveness of our resultant method, we calculate its computational cost compared to a recent study that introduced a method to sample the same odd-dimensional Wigner distribution using Monte Carlo methods [19]. The stationary phase method scales as 3 t 2 +1 for calculations to full precision while the best alternative qutrit algorithm scales as ∼3 0.8t for calculations to 10 −2 precision. This allows the stationary phase method to be more useful for intermediate-sized circuits Related approaches include recent proposals to use non-Gaussianity as a resource in continuous (infinite dimensional) quantum mechanics such as optics [23][24][25], as well as the well-established field dealing with semiclassical propagators in continuous systems [7,26,27]. This paper is organized as follows: in Section 2 we review the stationary phase method in continuous systems. In Section 3 we briefly introduce aspects of the WWM formalism and develop its quantum channel representation of gates. We introduce closed-form results for quadratic Gauss sums in Section 4 and this motivates the introduction of the periodized stationary phase method for discrete odd prime-dimensional systems in Section 5. Section 6 compares and characterizes this discrete stationary phase method with its continuous analog and uniformization. The periodized stationary phase method is then used to evaluate the qutrit π/8 gate magic state in Section 7. We then establish a relationship between stabilizer rank, Wigner phase space dimension and number of critical points in Section 8. These results allow us to produce an expression with 3 t 2 critical points for t states in Section 9. We discuss future directions of study in Section 10 and close the paper with some concluding remarks in Section 11.

Stationary Phase in Continuous Systems
The stationary phase method is an exact method of reexpressing the integral of an exponential of a function in terms of a sum of contributions from critical points: (1) where S[x j ] is a functional over x j , and the sum is indexed by critical "trajectories" j defined by and we have chosen to exercise our freedom to factor out the term i from S, for later clarity. These contributions from critical points correspond to Gaussian integrals, or Fresnel integrals in the general case, which can be analytically evaluated, or higher-order uniformizations corresponding to two or more critical points. In general, the higher order contributions asymptotically decrease in significance. When S is the classical action of a particle, Eq. 2 defines the critical points as classical trajectories that satisfy the Hamiltonian associated with S. (In this case, S is functional of x because it is the integral of the Lagrangian over time, the latter of which is a function of x andẋ). Terminating Eq. 1's expansion of S[x] at second order corresponds to making a first order approximation in : where the sum is over all classical paths that satisfy the boundary conditions ∂S ∂x = ∂S ∂x = 0. The subscript j in S j (x, x ) indexes the generally infinite number of classical trajectories that satisfy the initial x and final x conditions (i.e. x and x are not generally sufficient to define the functional S[x(t)] ≡ S(x, x )).
Here we will be interested in the double-ended propagator, or quantum channel-the amplitude from one state to another state. Thus, we will discuss the stationary phase approximation in reference to two integrals. We will also be interested in the case that the action S can be expressed as a polynomial in its arguments and thus becomes a bona fide function. Given functions ψ(x ) and ψ (x ) that represent these states, and an action S that represents the system of interest (usually this is an integral of the system's Lagrangian and so is related to its Hamiltonian), the propagator between them can be written as where the action S is assumed to be a polynomial in C[x 1 , . . . , x n , x 1 , . . . , x n ] (a complex-valued function of x i and x n ), S j corresponds to the action evaluated at the jth zero of ∂S ∂(x ,x ) -the jth "critical point" of S-and S 2 j is a function of the second derivative of S evaluated at the jth critical point.
As an example, if one chooses to represent the propagator in terms of initial and final position, ψ| = q| and |ψ = |q respectively, then the above equation becomes where the sum is over the critical points q j and q j of the polynomial action S function. The first term is often called the Van Vleck propagator or, more precisely, the Van Vleck-Morette-Gutzwiller propagator [28][29][30].
However, we have purposely been general with our choice of representation ψ and ψ because instead of the position representation, we will employ the "centerchord" representation. This "center-chord" representation is more suitable for discrete Hilbert spaces [22] and produces the Wigner-Weyl-Moyal (WWM) formalism. Using this representation with the stationary phase method means that stationary phase expansions must be made around the centers x and x . We briefly develop aspects of this formalism that are pertinent to gate concatenation in the next section.

Discrete Wigner Functions: Wiger-Weyl-Moyal Formalism
The discrete WWM formalism of quantum mechanics is equivalent to the matrix and path integral representations [21,31]. A complete description of the odd-dimensional WWM formalism can be found elsewhere [20][21][22]. Here we introduce the basic formalism of Weyl symbols of products of operators or gates, which will be necessary in order to develop a propagator-like treatment.

WWM Formalism Basics
We consider the case of d odd. We will later restrict this further to d = p h for p odd prime and h > 0, not to be confused with the reduced Planck's constant. We set = d 2π . We label the computational basis for our system by j ∈ 0, 1, . . . , d − 1, for d odd. We identify the discrete position basis with the computational basis and define the "boost" operator as diagonal in this basis: where ω ≡ ω(d) = e 2πi/d . We define the normalized discrete Fourier transform operator to be equivalent to the Hadamard gate: This allows us to define the Fourier transform ofẐ as the "shift" operator:X which acts asX δq |j ≡ |j ⊕ δq , where ⊕ denotes mod-d integer addition. The Weyl relation holds forX andẐ: Consistent with [20][21][22], we define the symplectic matrix for I n the n-dimensional identity. The Weyl symbol (Wigner function) of an operatorρ can be written is a 2n-vector corresponding to a conjugate pair of n-dimensional center momenta and positions while λ ≡ (λ p , λ q ) corresponds to its dual chord momenta and positions. The Weyl symbol of a unitary gate is in C and any function in C can always be written in the form where S(x) ∈ C[x] and x U (x) = 1. Given x ∈ (Z/dZ) n we can make this more precise by restricting S(x) to be a polynomial in the ring C [(Z/dZ) n ] where deg(S) < d (without loss of generality). We call S(x) the center-generating action or just simply the action due to its semiclassical role [20].

Weyl Symbols of Quantum Circuits
We are interested in the Weyl symbol of products of operators because quantum circuits evolve by products of operators corresponding to elementary gates. The Weyl symbol of a product of an even 2N number of operatorsÂ 2N , . . . ,Â 1 has a phase determined by this symplectic area for the odd number 2N + 1 of centers: The generalized symplectic matrix for N degrees of freedom can be defined through the introduction of the matrix J N and H N [20]: The product of J N and the inverse of H N is The symplectic area for an odd number of centers is defined to be Notice here the convention of associating a factor of d −n to every sum over x i ≡ (x pi , x qi ). This is a convention that we will continue to use throughout this paper. For instance, we expect for any quantum gateÂ =Û or density matrix statê A =ρ.
For example, the Wigner function of two operatorŝ A 2Â1 can be written: is the symplectic area associated with three centers x 1 , Proceeding with a propagator-like treatment, here we will be interested in propagating from some stateρ toρ and so will be interested in extending these definitions to capture four operators-ρ ÛρÛ † . By Eq. 14, this requires the symplectic area for five degrees of freedom: More generally, we will be interested in an even numbers of operators:ρ Û 1 · · ·Û m ρÛ † m · · ·Û 1 . To accomplish this, we can extend the definition for the Weyl symbol of the product of 2N operators given in Eq. 14 to an odd number (2N − 1) of products by considering the Weyl symbols of 2N products of operators and settingÂ 2N =Î, which has the corresponding Weyl symbol I(x) = 1. This leaves the x 2N variable alone in the argument of the exponential, free to be summed over to produce a Kronecker delta function. The remaining exponentiated terms are defined to be Finally, consider the Weyl symbol of the product of the four operatorsÎÂ 3Â2Â1 : Notice that the prefactor fell by d 2n compared to Eq. 14 since that is the dimension of the vector inside the Kronecker delta function. This agrees with our established standard of a factor of d −n for every sum over x i since after the Kronecker delta function is summed over, the resulting expression involves only a double sum over phase space and so we expect two factors of d −n .
With the WWM formalism for the Weyl symbols of products of operators thus established, we move on to consider quantum channels, which are double-ended propagators.
In general, past treatment of the quantum propagation of a stateρ by a unitaryÛ in the WWM formalism have been content with Eq. 23, settingÂ 2 =ρ , and 20,21]. However, such a treatment intimately ties the propagator to the initial state,ρ , and to act on a final state necessarily involves using Eq. 20 settingÂ 2 =ρ, the final state, andÂ 1 =ÎÛρ Û † . This leaves an intermediary phases ∆ 3 to sum over and as a result does not produce a self-contained double-ended propagator (a quantum channel) that acts on states by summation without any additional functions or phases, as in Eq. 3. A double-ended propagator form is often both more familiar and more useful for generality.
Using the additional general formulae defined in this Section, notably Eq. 14, we can develop such a doubleended propagator or quantum channel.

Weyl Symbols for Quantum Channels
We define U U * (x, x ) to be the Weyl symbol of a quantum channel that can take any initial state ρ to a final states ρ by summing over their phase space variables: (25) To obtain U U * (x, x ), we use Eq. 14 for 2N = 4. We set A 4 (x 4 ) = 1 and sum x 4 away, set A 2 (x 2 ) = ρ but then discard it along with its sum over x 2 , and set . Relabelled the variables of summation, this produces: (26) where we also used the identity ∆ 5 (x, . This is equivalent to using Eq. 23 for A 2 (x 2 ) = ρ (x 2 ) and discarding the sum over x 2 . Notice in Eq. 26 the inclusion of the full phase ∆ 5 in the expression. As a result, acting on ρ(x) or ρ (x ) by U U * (x, x ) as in Eq. 26 simply involves summing over x or x respectively; there is no intermediate phase ∆ 3 as in Eq. 20 that must also be summed over.
Our choice of normalization in Eq. 26 allows us to continue the standard of including a factor of d −n to every sum over a pair of conjugate phase space degrees of freedom in Eq. 25. This means that d −2n x, for the Weyl symbol of this propagator is perhaps clunky, since it can be mistaken for the Weyl symbol ofÛÛ † , U U † (x). But as Eq. 26 denotes, it is the Weyl symbol of a double-ended propagator and so has two arguments, x and x 1 1 To avoid any confusion, we point out that U U * (x, x ) = U U † (x) = 1; U U † (x) is the Weyl symbol ofÛÛ † =Î and so is a function of only one variable, x; U U * (x, x ) can most appropriately be associated with the Weyl symbol of the superoperator Û •Û † , where and • denote the operators the superoperator acts on, and so is a function of two variables, (cont'd on pg. 7) For Clifford gatesÛ , the associated Weyl symbol U (x) is a non-negative map; U (x) takes non-negative states to non-negative states. This is clearer for U U * (x, x ), which becomes a non-negative real function: U U * (x, x ) ≥ 0. This parallels the non-negativity of Wigner functions ρ(x) of stabilizer statesρ. The extension of this set by π 8 gates or magic states to allow for universal quantum computation introduces contextuality into U U * (x, x ) and ρ(x), respectively, which means that they can now also have (real) negative values.
In the next section we will explore how the nonnegative real function U U * (x, x ) produced by a Clifford gate can be described by a single Gauss sum.

Gauss Sums
We now restrict to d = p h for p odd prime and h ∈ Z.
Consider, a simplified definition of the Gauss sum from [32,33]: The result of evaluating this sum can be simplified and summarized in the following list [32]: There is a more broad application of this result to padic number theory, a topic that is briefly introduced in the Appendix A, but is beyond the scope of this work.
A simple case of this result produces one of the most often used equations in this paper: Proof A = 0 and so by Proposition 1(b), the sum is equal to 0 if v = 0 since there does not exist any u ∈ Z/p n Z such that Au = 0. On the other hand, if v = 0, then Proposition 1(b) states that the sum is equal to This is the definition of the Kronecker delta function.
We now give an example of the application of this Proposition to Clifford gates.

Gauss Sums on Clifford Gates
As we saw in Eq. 13, the Weyl symbol of a unitary gate such as a Clifford gate is [22]. We can write the action for a Clifford gate as: where B is a 2n × 2n symmetric matrix with elements in Z/dZ. This simpler form makes Clifford gates act on states in a correspondingly simpler manner. Namely, Clifford gates transform Wigner functions by a linear symplectic transformation [22]. One way to see this is by deriving functions that act as Hamiltonians from the above action, which can be transformed into the computational basis by Legendre transform [22]. These Hamiltonians are similarly harmonic and so their effect on states can be described that the action of harmonic Hamiltonians acting on continuous degrees of freedom can be: by rigid symplectic transformations. In the discrete case, this translated to taking Wigner arguments (x) to themselves while preserving area (a symplectic linear transformation). We proceed to see how this is manifested from the perspective of the double-ended propagator Let us evaluate U U * (x, x ) from this perspective: where U (x) = exp 2πi Applying Proposition 1 produces where We can define a symplectic matrix M by relating it to the symmetric matrix B through the Cayley parametrization: (32) By Corollary 1, since A = 0, the sum over x 3 is nonzero if and only if v = 0 (see Appendix B): Or, in other words, we are presented with the usual plane wave sum identity for a Kronecker delta function.
For values of x and x that satisfy There is only one solution x given x , and so d −2n . This result affirms that the Clifford gate symplectically transforms Wigner functions by point-to-point symplectic (area-preserving) permutation as discussed in the beginning of this Section.
Note that this approach is no longer suitable if S(x) has powers higher than quadratic, or, if even though it is quadratic, B and α do not have elements in Z. This is because the Gauss sum results from Proposition 1 cannot apply in these cases. To handle a more general form for S(x), the the Gauss sum results of this section must be extended. We explore such a generalization in Section 5 where we develop the discrete and periodic stationary phase method.

Stationary Phase
We will now proceed to momentarily consider the more general case of a sum over an exponentiated polynomial As added incentive, U U * (x, x ) is naively a simpler function to deal with than the Weyl symbol of a unitary operator, U (x). U U * (x, x ) is real-valued just like the Weyl symbols (ρ(x)) of density functions (Wigner functions), while U (x) is generally complex-valued. Also, U U * (x, x ) resembles traditional propagators more closely in that it can be said to take states from x to x, whereas U (x) requires pairing with U * (x) and summation over intermediate values with the appropriate phase as in Eq. 23, to act on a state. but such that ∂S(x) ∂x still has coefficients in Z p . We do this so as to make use of a result from this more general treatment since, in this case, we can no longer derive a simple Gauss sum result. However, if we split the sum up, we can obtain a closed form in terms of several Gauss sums or sums of higher order with arguments in To accomplish this we express the sum, I, in terms of local terms Ix: where j < m. So far this is simply a reorganization of the order of terms in the summation. However, since these local terms are over coarser periodic domains, their polynomial power can be reduced in much the same way that polynomial powers can be lowered in proximal Taylor series expansions over a continuous domain. In this case, since we are dealing with a periodic discrete domain, such a simplification relies on the "polynomial version" of Taylor's theorem [34], where the sum is over the components of 0 ≤ α i ∈ Z n such that |α| ≡ i (α i !). Using the polynomial version of Taylor's theorem, Fisher proved the following Lemma [33]: , and a point a ∈ Z n p , let S (<k) a (x) denote the Taylor polynomial of degree k − 1 for S(a + x). Assume that p t S(x) and the partial derivatives have coefficients in Z p , for some integer t, and let = k = min{t, log k/ log p }.
We express the quadratic Taylor polynomial of S(x) at a as S <3 is the Hessian matrix of f at a. We proceed as follows: (c) Assume that det(H a ) = 0, and choose an integer Then the following are equivalent: ∇S(a) = 0 ∈ Z n p and α ≡ a mod p j , and Iā = p nm/2 e 2πiS(α)/p m G m (H α , 0).
As an example consider S(x) = x 3 + 2x 2 over Z/3 2 Z: We reexpress this as a sum of local terms Ix: where Since S(x) ∈ Z p it follows that t = 0 and so = 0. 0 mod 3 = ∇S(x) = 3x 2 + 4x has solution x = 0. Therefore, by Theorem 1(a), Therefore, due to Theorem 1, we are able to simplify the sum, I, to just one term, I0.

Diagonal Unitary Example
An important subset of non-Clifford gates include diagonal gates with rational eigenvalues. These can always be written asÛ = exp 2πi . As an example, we consider an action that encompasses the generalized π 8 -gates for qutrits [35,36]: where C, B ∈ Z and α ∈ Z n . We consider S 9 as a polynomial over Z/3 2 Z so that its Weyl symbol is . We note that though the coefficients of this polynomial are still in Z as the actions were in Section 4, it is cubic and thus must be treated by the method of stationary phase instead of as a single quadratic Gauss sum.
It is easy enough to verify that this corresponds to a unitary operator (see Appendix C): So we now consider the double-ended propagator corresponding to S 9 : where S(x 2q ) is a cubic polynomial shown in Appendix C. For some values of A, B, and α, some phase space points (x, x ) produce three Gauss sums while others produce none. This can be found by seeing which values of (x, x ) cause to equal zero. Examples of this can be found in Appendix C. This illustrates the usefulness of the periodized stationary phase approximation in not only formalizing the number of Gauss sums that are necessary, but also in providing an easy way to find this number (i.e. finding the zeros of For general diagonal gates, the coefficients in their action polynomials with fall in Q instead of Z. We will develop a technique in Section 7 that will show how to change the domain of summation from Z/p h Z to Z/p h Z, where h > h, such that the resultant equivalent action has coefficients in Z thereby allowing for its treatment by the method of periodized stationary phase developed here. But first, let us characterize and discuss this periodized stationary phase method.

Uniformization
The discrete "periodized" stationary phase method introduced in the Section 5 has many similarities to the stationary phase approximation in the continuous setting. This is discussed in more detail in Appendix D.
One of the properties of the continuous stationary phase approximation is that the number of Gaussian integrals it produces can be decreased by increasing the expansion in the exponentiated action to a higher order than quadratic. This produces a smaller number of higher order integrals, such as Airy functions (corresponding to a cubic expansion) and the Pearcey integral (corresponding to a quartic expansion). This procedure is called uniformization. Though this produces fewer terms, these integrals are correspondingly more difficult to evaluate compared to Gaussian integrals.
A similar uniformization procedure can be applied with the discrete "periodized" stationary phase method. As we shall see here, this reduces the number of quadratic Gauss sums that the method produces into a fewer number of higher-order Gauss sums.
According to Theorem 1, the discrete stationary phase method drops the size of summation by a factor of p m− m/3 and changes the summand to a sum of quadratic Gauss sums with closed-form solutions multiplied by phases. However, the stationary phase method can also decrease the domain by a larger multiple of p, but the resultant sum is over terms that are themselves sums with higher than quadratic order; they are no longer quadratic Gauss sums.
To prove part (b) of Theorem 1, Fisher appealed to Lemma 1 to show that if m ≤ kj − k = 3j − 3 and ∇S(a) ≡ 0( mod p min{j,m−j} ), then for any x ∈ which can be rewritten in terms of the Gauss sum notation as above in Theorem 1 (b). This trivially generalizes to higher order polynomials: . . , ∂S ∂xn ) has coefficients in Z p . Define k as in Eq. 36. Given integers m ≥ j ≥ 1 and a point a ∈ Z n p with reduction a ∈ (Z/p j Z) n , define Iā = Iā(f ; Z/p m Z) as in Eq. 34. Then: If m ≤ kj − k and ∇S(a) ≡ 0 mod p min(j,m−j) , then Proof The same proof as Fisher's can be used here.
In particular, for k = 4, Completing the cube for n = 1 by setting x = ∇ 3 S 1 3 x + p −j ∇ 2 S/∇ 3 S this can be rewritten as: for This is a discrete sum analog to the Airy function [37].
Similarly, for k = 5, setting for n = 1 allows the equation to be rewritten as: for This is a discrete sum analog to the Pearcey integral [38].
Higher order instances can be similarly developed leading to discrete sum analogs to integrals familiar in catastrophe theory [39].
For n degrees of freedom with a domain of summation of Z/p m Z (as might be produced by, for instance, n p m -dimensional qudits), the naive numerical summation involves a sum over p nm terms. Using the Gauss sum k = 3 simplification or appropriate uniformization level for k > 3, this sum can be reduced to a sum over p n terms involving Gauss, Airy, Pearcey etc. sums, which number j p m−j and can be pre-computed and stored for use during the summation. These terms can perhaps also be approximated numerically instead of tabulated, thereby eschewing the exponential cost in storage. This is the current approach with Airy functions and Pearcey integrals in the continuous case in computation, for instance [40].
We note that including non-Clifford gates within the WWM formalism was not possible under the previous derivation of the formalism in terms of powers of [20][21][22]. This is because earlier derivations of the discrete WWM formalism begain with the continuous case, where the stationary phase method is able to be applied in its traditional form, and then "discretized" and "periodized" the final propagator. As a result, the propagator could only be written to order , which just includes Clifford propagation, and higher order corrections could not be derived because anharmonic trajectories cannot be obtained by "periodizing" to a discrete Weyl phase space grid. In this way, this paper finally accomplishes this extension to higher orders of by instead using this different "periodized" stationary phase method, thereby solving this old problem for the first time. Just as in the old derivation, it finds that Clifford propagation is captured at order 0 , but unlike the old derivation it is able to formally ascribe a power of required to include (diagonal and their Clifford transformations) non-Clifford gates in a formal Taylor series, as discussed in Section D.

Application to universal gatesets
The T -gate is a non-Clifford single qudit gate. It extends the Clifford gateset to allow for universal quantum computation [35]. Generalizations of the T -gate to qudits are frequently called π/8 gates.
Here we will show how the discrete "periodic" stationary phase method can be applied to the Wigner function of the qutrit π/8 gate magic state to obtain a summation over quadratic Gauss sums: π/8 gates differ for dimension p = 3 and p > 3 [35,36]. For p = 3, with ζ = e 2πi/9 , where v = (0, 6z + 2γ + 3 , 6z + γ + 6 ) mod 9. (60) Despite these differences, these π/8 gates are all diagonal with entries that are rational powers of e 2πi/p . This means that their Weyl symbol actions Here we will show how to perform this change in the domain of summation for a gate from Z/p h Z to Z/p h Z, where h > h, such that the resultant equivalent action will go from having coefficients in Q to Z. This is a useful technique for general diagonal gates, since the coefficients in their action polynomials will generally fall in Q instead of Z. Reexpressing their Weyl symbols in this way will allow for their evaluation by the method of stationary phase. Though we will only show this for the qutrit π/8 gate as an example, this technique can be equivalently applied for all diagonal gates with rational coefficients, as well as their Clifford transformations, which only symplectically permute the action.
To demonstrate this we examine a particular π/8 gate for qutrits and its corresponding magic state.
Let z = 1, γ = 2 and = 0 so that The corresponding Weyl symbol is , when the domain is increased from 3 to 3 2 ; we want to consider the sum over x q of U π/8 (x) = U π/8 (x q ) and reexpress it in terms of an exponential such that it is a sum over the larger space Z/3 2 Z of a polynomial with integer coefficients (see E): Notice that the first line consists of a sum of an exponentiated polynomial with coefficients in Q while the second line is a sum over a larger space, but this time of an exponentiated polynomial with coefficients in Z. Thus, we have accomplished what we set out to do.
Since we will not be performing the sum above in isolation but with other exponentiated terms, so it is important to establish a slightly more general result: Applying this identity to U π/8 U * π/8 (x, x ) transforms it into an expression with exponentiated polynomials over a larger domain of summation and whose coefficients are in Z. This result is now amenable for the discrete "periodized" stationary phase method result given by Theorem 1. This produces (see Appendix E): , and G 0 is defined in Eq. 27.
We are interested in the magic state of this gate, which corresponds to it acting onĤ |0 = |p = 0 . The Wigner function of |p = 0 is ρ (x) ≡ 1 3 δ xp,0 . Hence, (where, again, details are shown in Appendix E).
If we let S(x p , x q ) be the phase, we note that ∂S ∂xp = ∂S ∂xq = 0 mod 3 and ∇ 2 S ≡ H = 0 mod 3 0 ∀ x p , x q . Hence, evaluation at any phase space point, or linear combination thereof, requires summation over the three reduced phase space pointsx 2q ∈ Z/3Z since they are all critical points.

Application
We consider as an example, computation of the qutrit circuit outcome: This corresponds to the probability of the outcome |0 0| after k qutrit magic states are acted on by a random Clifford circuit U C on n ≥ k total qutrits. A recent study introduced a method to sample this distribution using Monte Carlo methods on Wigner functions [19]. Results were demonstrated for one to ten magic states, in a system of 100 qutrits, that required 10 5 to > 10 8 samples, respectively, to attain precision (P sampled − P ) < 10 −2 with 95% confidence [19].
Before the Clifford gates are applied, in the Wigner picture the circuit can be described by k products of Eq. 66 and (n − k) products of δ xq j ,0 . The random Clifford gate can be described as an affine transformation by a symplectic matrix M and vector α [4,41], as described in Section 4.1, such that The form of the M and α for the Clifford gates are given in [4]. The final contraction to |0 0| corresponds to a sum over all n of the x pj and x qj except for the first degree of freedom, which is only considered at x q1 = 0 (and all x p1 ). We note that this is a strong simulation algorithm. Therefore, it is interesting to compare the number of terms produced by the stationary phase method that must be summed over in this strong simulation with the number of samples that must be taken in the prior strong simulation by Pashayan et al. However, a direct comparison of the two on an equal footing is a bit blurred by the fact that the latter only reports results based on a Monte Carlo sampling of their terms, producing a result correct only up to a precision = 0.01. Their explicit evaluation would produce a result correct up to precision = 0.0 and would likely scale far worse (we estimate at least 3 2t for a naive evaluation of the 2t dimensional Wigner function, see Fig. 1). Monte Carlo calculation of finite sums reduce the number of terms that must be evaluated for intermediate to large sums, and we assume that is the case with Pashayan et al.'s results. Thus we will proceed to compare the two methods with the caveat that Monte Carlo sampling of the sums produced by the stationary phase method will likely also improve its performance for intermediate to large values of t.
For k magic states in a 100-qutrit system, we have k one-dimensional sums (integrals) over x 2qj and 100 sums over all the x pj and 99 sums over all the x qj except x q1 , resulting in a total of (k + 199) sums. This can be concisely written by letting x ≡ (x p1 , . . . , x p100 , x q1 , . . . , x q100 ) ≡ (x p , x q ) and for There are 199 sums in Eq. 69 and the ρ π/8 contain k more, for a total of (k + 199).
The restriction on the sum over x to be in D, defined in Eq. 70, can be treated as an additional Kronecker delta function modulo 3 2 , δ(M −1 (x − α 2 ) − α 2 ) 101 ), multiplying the full sum over x .
We can begin reducing the number of sums by first summing out all δ(x qj ) in Eq. 69, which will replace all x qj in this additional Kronecker delta function for j ∈ {t + 1, . . . , 100}. We can then proceed to sum away all x pj for j ∈ {t + 1, . . . , 100} since they are not present anywhere in the full summand and put in the appropriate factors of 3. The additional delta function can now only include terms x pi and x qi for i ∈ {1, . . . , t}. Any such x pi term can now be chosen to be summed away, which replaces it on the corresponding ρ π/8 state. The remaining (t − 1) x pi terms can only be linear terms in the ρ π/8 exponential functions. Thus, they can be summed away to produce (t − 1) Kronecker delta functions that are linearly independent and so can be used to sum away (t − 1) x qi . There will remain left over one x qi variable along with the t x 2qi variables-a remaining total of (t + 1) sums. This process is illustrated schematically in Table 1.
If the additional Kronecker delta function representing the restricted sum never had any remaining x pi terms then it could simply be used to sum away and replace a x qi term with other x qi s, and one could then proceed as outlined before by summing away all x pi to produce Kronecker delta functions and summing away all x qi . This would result in only t remaining sums over the x 2qi .
Therefore, in the worst case, the total number of variables that will need to be summed in the t π/8 gate magic states will be (t + 1) consisting of Gauss sums and leading phases. This is worst-case scaling of 3 t+1 terms in a sum. We show how this scaling compares to the alternative algorithm in Fig. 1.
ActingÛ C on the ket instead produces M(x + α 2 ) + α 2 in the arguments for the initial stabilizers and ρ π/8 functions with a sum over x restricting x q1 = 0. The preceding argument then can be made again since the rows in the top half and bottom half of M(x + α 2 ) + α 2 correspond to linearly independent equations and so can be treated independently [4]. This results in a worstcase of (t + 1) sums again. We have thus found that a π/8 gate (and magic state) can be captured with only p critical points. This is commensurate with its role as the "minimally" non-Clifford gate [42] and that its direct products have the smallest stabilizer rank [43].
We proceed to improve this scaling further by developing a relationship between state stabilizer rank and the number of critical points necessary to represent that state and then leveraging this result with the π/8 gate magic state's optimal stabilizer rank for two qutrits.

Stabilizer Rank and Stationary Phase Critical Points
We will now see that a relationship can be established between the stabilizer rank of this state and the number of these critical points. This is in contrast to the relationship between the amount of negativity present magic states and contextuality, which appear to be inversely related [44]. This suggests that negativity is

Number of Terms in Sum
Number of π/8 Gates not the most efficient way to introduce "magic" or noncontextuality in a practical algorithm, and indeed we find this to be the case here.
In the stationary phase method applied to the infinite dimensional (continuous) case, the critical points correspond to intersections between Gaussian manifolds, the continuous generalization of stabilizer states [22]. However, in the "periodized" stationary phase method examined here, we cannot expect the same relationship to hold. Nevertheless, a relationship can still be estabilished between the number of stabilizer states that a state can be expressed in terms of-its stabilizer rankand the number of critical points necessary to represent that state:

Theorem 3 (Stabilizer Rank and Critical Points) Let
|Ψ be a n-qudit odd-prime-p-dimensional state that can be written as an equiprobable linear combination of the p n logical basis states. If |Ψ can also be expressed as an equiprobable linear combination of p m (m ≤ n) orthogonal stabilizer states, which can be written, after some Clifford transformationÛ C , as products of m orthogonal single-qudit stabilizer states {|φ ij } and (n−m) single-qudit stabilizer states {|ψ ik } that are Step #

Description
Resultant Equation Replace a x p i term using remaining delta function.
Sum away remaining (t − 1) x pi terms. Table 1: Step-by-step schematic illustration of the worst-case simplification of Eq. 69 to (t+1) sums described in the text. For clarity, in the equations in the third column, , and the rest, ρ rest π/8 (x q i , x2q i ), which is not dependent on x p i . Furthermore, Kronecker delta functions are mostly written in shorthand in terms of just their comma-separated arguments. It should be possible from the text to determine the domains of the sums and products as well as the form of the functions listed. not logical states, where c i ∈ C and |c i | 2 = |c j | 2 ∀ i, j, then the Wigner function ofρ = |Ψ Ψ|, ρ(x), can be expressed in terms of ≤ p m quadratic Gauss sums of dimension p n−m (i.e. ρ(x) can be written in terms of ≤ p m critical points of dimension p n−m ).
Proof For a fixed j, there must exist a one-qubit Clifford transformation that takes |φ ij to |0 . It follows after this transformation, for all other i, |φ ij are also orthogonal one-qudit stabilizer states and so must be the other logical states on that qudit. Proceeding in this manner on the remaining one-qudit {φ ij } produces a Clifford transformationÛ C such that where |i is the base p representation of a product of single-qudit logical (stabilizer) states. We can rewrite the (n − m) remaining qudits in the base p representation: where again |j is written in base p representation. |c i,j | 2 = |c k,l | 2 ∀ i, j, k, l since we are given that |c i | 2 = |c j | 2 ∀ i, j and {|ψ ik } are single-qudit stabilizer states that are not logical states, and so are equiprobable in terms of logical states.
It follows that where S(i, j) ∈ R is some polynomial function of i and j when they are broken up into their m and (n − m) base p constituents, respectively. Therefore, the Weyl symbol of the diagonal gatê U that takesρ = c i,j |i |j = |i |ψ i , and is thus a Clifford gate since it takes a stabilizer state to a stabilizer state. Hence, for fixed i, S(i, j) must be a quadratic polynomial in j.
Therefore, each i in the sum indexes a term made up of a p n−m -dimensional quadratic Gauss sum. This implies that the Wigner function has ≤ p m quadratic Gauss sums of dimension p n−m , or equivalently, critical points.

Two Qutrit π/8 gate magic States
A single qubit T -gate magic state is not a stabilizer state and has a stabilizer rank of two. Two copies of the qubit magic state are somewhat remarkable in that they can also be written in terms of only two stabilizer states: where each line consists of a stabilizer state.
In the qutrit case, a similar pattern holds. A single qutrit T -gate magic state consists of at least three stabilizer states. Moreover, two copies of the qutrit magic state can also be written in terms of only three stabilizer states: where again each line consists of a stabilizer state (F is the Hadamard gate or discrete Fourier transform). Qubit strong simulation algorithms based on stabilizer rank have been able to leverage this fact to halve their exponential scaling of terms to O(2 0.5t ) [43,45]. The performance can be improved to O(2 ∼0.468t ) by using the fact that six qubit T -gate magic states can be written in terms of only seven stabilizer states [46].
The dependence of the number of critical points on the intermediate values of x 2q for the π/8 gate magic state, means that generally the number of critical points scales exponentially with number of π/8 gates. However, if this stabilizer rank for two qutrit π/8 gate magic states can be used, the exponent is reduced by a factor of 2.
Eq. 76 can be transformed into a more presentable form by the action of the Clifford controlled-not gate, C 12 , which adds the value of the low qutrit to the high qutrit mod 3 and produces: We recognize this as a pointer state, in which the high qutrit tells us the state of the low qubit. These states are stabilizer states, and we may write: This form satisfies the description of Theorem 3 and so implies that the Wigner function ofÛ ⊗2 π/8F ⊗2 |00 can be represented by three one-dimensional quadratic Gauss sums, just like forÛ π/8F |0 -the single π/8 gate magic state.
To see how this is true in the Wigner representation, we consider two copies of the magic state given by Eq. 66: We are allowed to transform x 2q1 and x 2q2 by a symplectic transformation since the sum over these variables is invariant to symplectic transformations; the sum over a finite field of an exponential polynomial function that is linearly transformed merely changes the order of summation. However, to calculate the full trace as in Eq. 69 involves summing over x as well, which involves terms that are also exponential polynomial functions but over a subset of the full domain D given by Eq. 70. A symplectic transformation does not leave this restricted sum invariant, it additionally transforms the restriction (Eq. 70) by some matrix M and vector α. Fortunately, we consider all such transformations of the restriction in our analysis above due to already considering any Clifford transformation on the π/8 gates. Therefore, we are free to further symplectically transform the x variables within our analysis.
The two-qutritĈ 12 controlled-not gate transforms (x p1 , x p2 , x q1 , x q2 ) to (x p1 , x p2 − x p1 mod 3, x q1 + x q2 mod 3, x q2 ) [22]. Acting on both x and x 2q1 and x 2q2 with (Ĉ 12 ) 2 produces: where We can see that for each of the three values that x 2q1 takes in the sum, the exponent is a quadratic polynomial for x 2q2 for fixed x with coefficients in Z/3Z and so is a quadratic Gauss sum. We further see that the x q1 are similarly cubic and that for each value of x q1 the exponent is a quadratic polynomial for x q2 for fixed x 2q1 and x 2q2 . x p1 and x p2 are both linear terms just as before, and therefore the same arguments for the scaling of the marginal trace hold here.
Thus, the end result is a sum of three one-dimensional Gauss sums, exactly the same number of Gauss sums as we found for a single π/8 gate in Section 7.
Thus the number of critical points is reduced to 3 t/2+1 so that the cost of simulation of the example in Section 7 is O(3 t/2 ). This can be seen by the lowest black curve in Figure 1.
This performance is better than the Monte Carlo algorithm of Pashayan et al. [19] and is achieved for an algorithm with no Monte Carlo sampling error. Use of Monte Carlo would further improve this scaling. It is interesting to note that weak simulation of qutrits by the method of [47] scales as 3 0.32t and so it is perhaps possible that Monte Carlo could improve this strong simulation algorithm to be more efficient than weak simulation.

Future Directions
One of the central pillars of the stationary phase method used here is that non-contextual operations are efficiently classically simulable. This means that the Weyl symbols of Clifford gates reduce to a Gauss sum and affect simple symplectic transformations in Weyl phase space. Furthermore, stabilizer states have non-negative Wigner functions that characterize functions of affine subspaces of the discrete phase space.
For a similar approach to work for qubits, the same operations must be non-contextual in order to be able to be free resources. This requires the WWM formalism to be extended from two generators, p and q, to three generators that become Grassmann elements [5]. The resultant Grassmann algebra cannot be treated over disjoint states in phase space [48], such as we have done here, and so a Grassmann calculus must be used. It would be interesting to see if such the stationary phase method can be applied to the Grassmann algebra and produce a similar treatment of qubit π/8 gates.
Another interesting direction for future study regards the mathematical relationship between π/8 gate magic state stabilizer rank and the number of critical points in exponentiated multidimensional polynomials. The improvement in scaling found in Section 9 found by a symplectic transformation of the x 2q1 -x 2q2 plane is really due to the reduction in the cubic power of a polynomial with respect to one degree when it is reexpressed as an inseparable polynomial with the second degree of freedom. Theorem 3 strongly suggests that a similar simplification holds for higher numbers of π/8 gate magic states that are known to have lower stabilizer ranks [43]. The form of this relationship may be helpful in finding such lower stabilizer ranks for even higher numbers of π/8 gate magic states than are currently know, as well as establishing concrete bounds.
A related interesting question regards approximate stabilizer rank of magic states instead of their exact stabilizer rank. This is sufficient for weak simulation and frequently leads to a more efficient algorithm since we can get away with introducing error in the probability distribution we are sampling that is supposed to represent the exact probability distribution. A weak simulation result for qudits has recently been developed [47]. However, it would be interesting to see if there is a similar result to Theorem 3 that deals with approximate stabilizer rank and if there is some such "approximate" analog to the discrete stationary phase method that is useful for weak simulation.

Conclusion
This paper established the stationary phase method as a way to understand the order 0 non-contextual Clifford subtheory in the WWM formalism, producing single Gauss sums, and higher order contextual extensions of the subtheory, in terms of uniformizationshigher order sums-that can be reexpressed in terms of a sum over critical points or Gauss sums. This firmly tracks with the same relationships that exist in the continuous infinite-dimensional Hilbert space treatment of Gaussianity and non-Gaussianity, even though the stationary phase method introduced here for qudit systems differs in that it is a "periodized" stationary phase.
We discussed these differences and similarities between the continuous and discrete case. This involved comparing this measure of non-contextuality to negativity. We found that the usage of higher order uniformizations through the stationary phase method is more efficient than using negativity for π/8 gate magic states, at least in the manner that has so far been tried. This also seems to have been noticed in the discrete community, which has turned to favor stabilizer state decomposition of magic states [43,45]. By relating the stabilizer rank of magic states to the number of critical points necessary to treat them, this paper falls in line with this latter approach in treating non-Clifford gates.
We found that we are able to calculate a single qutrit marginal from a system consisting of t π/8 gate magic states that are then evolved under Clifford gates, with a sum consisting of 3 t+1 critical points corresponding to closed-form Gauss sums when the magic states are kept separable. This scaling improves to 3 t 2 +1 when pairs of magic states are rotated into each other in accord with the optimal two-qutrit π/8 gate magic state stabilizer rank. We showed that the latter scaling improves upon the current state-of-the-art.
All of this taken together establishes the usefulness of contextuality for practical application of classical simulation of qudit quantum algorithms through the venue of semiclassical higher order corrections in accomplished by the stationary phase method.

A p-adic Numbers
For prime p, we define the p-adic order or p-adic valuation, ν p (n), of a non-zero integer n to be the highest exponent ν such that p ν divides n and set ν p (0) = ∞.
p-adic numbers Q p can be written as: for c i ∈ {0, 1, 2, . . . , (p − 1)} and k is an integer. p-adic integers, Z p , are p-adic numbers where c i = 0 for all i < 0. Let us define the absolute value of a p-adic number n, |n| p , to be the inverse of p taken to the power of its valuation: |n| p = 1 p νp(n) . Hence, we can define the metric |n − m| p to denote the distance between two p-adic numbers m and n. Notice that this means that m and n are "close" together if their distance is a large power of p, which is the opposite expected from the Euclidean metric. As a result, this p-adic formalism presents an alternative way to complete the rational numbers to the real numbers R; completing the rationals with respect to the p-adic metric, (Q, |•| p ), produces the p-adic numbers.
The results presented in this paper come from p-adic number theory. However, it should be clear from this brief presentation that positive integers and positive rational numbers with terminating base p expansions will have terminating p-adic expressions (Eq. 82) that are identical to their base p expansions. Since these are primarily the cases we will be concerned with in this paper, a more thorough understanding of p-adic theory is not really strictly necessary.

B Gauss Sums on Clifford Gates
We define the Hessian and the scalar so that the sum can be rewritten to make use of Proposition 1b: where we assume det B = 0. We drop the arguments on u, v, and c and such subsequent terms for conciseness.
Since B can be expressed with entries in Z/pZ it follows that (B − J B −1 J ) −1 ∈ Z p and so u ∈ Z p . Thus, G 1 (A, 0) = 0 ∀ x, x , x 3 ∈ (Z/pZ) 2N by Proposition 1b.
In preparation of the final sum over x 3 , we expand out the exponent's phase: Gathering powers of x 3 , we define the Hessian where we made use of the property that B, B − and we define the scalar for β ≡ − 1 2 J T (B − J B −1 J ) −1 . In summary: the vector and

This allows us to rewrite
By Proposition 1b, since A = 0, the sum over x 3 is non-zero if and only if v = 0: Thus we find that the sum over x 3 is non-zero i.f.f. x = M x + α 2 + α 2 and for these values the phase 2c is equal to zero too. We can see this last fact by examining the quadratic, linear and constant parts of 2c w.r.t. x separately.
Quadratic terms in x : Using the results from the same simplification from the work on the quadratic part to simplify the first two terms in the linear part, we find: Constant terms: In particular, for these values of x and x , U * U (x, x ) = d 2n since the phase 2c is equal to zero and G 1 (0) = d 2n for a 4n dimensional summation variable (x 1 , x 2 ).

C First Non-Trivial Example
We verify that the operator with the action S 9 corresponds to a unitary operator: We replace the x 3q sum and sum away x 3p : and c = −(2x q −α q )(−x p +x p +α p +(2x q −α q )(−B+C(2x q −α q )).
(110) where in the last step we could have replaced x 2q instead of x 1q with similar results.
This leaves a single sum over x 2q of an exponentiated cubic polynomial S(x 2q ) = C x 3 2q + B x 2 2q + α x 2q + c with coefficients in Z. We proceed to find the exponential argument's zeros {x 2q }, the stationary phase points: where all the arithmetic operations are understood to be taken mod 3 2 .
Therefore, by Theorem 1(a) and (b), and so U 9 U * 9 ((0, 0), (2, 4)) = 0. This example illustrates the usefulness of this periodized stationary phase method for evaluating non-Clifford gate Weyl symbols and its potential for simplifying or reducing the full sum.

D Comparison between the Discrete Periodized Stationary Phase Method and the Stationary Phase Approximation in the Continuous Case
The stationary phase method in the continuous case can be described as representing a integral over a continuous domain by a set of discrete points {x i }. These are the stationary points of the phase of the integrand. The phase is expanded by Taylor's theorem at these discrete points {x i }, and if the expansion is truncated at quadratic order, then Gaussian integrals result.
The notion of the stationary phase method differs when there is no longer a proper Euclidean metric to define continuous distances or areas for the integrand's domain. Similarly, the traditional version of Taylor's theorem does not hold. Instead, we must use the "periodized" version of Taylor's theorem [34] in a "periodized" stationary phase approximation. In this version, the points {x i } where the phase's derivative is zero correspond to critical points where the "periodized" phase is stationary, i.e., where the phase, taken every x i points, slows down so that it ceases cancelling its opposing contributions around the unit circle with equal weights.
In the discrete case we consequently look for stationarity along equally-spaced periodic intervals that span the whole summand, instead of at particular points of the integrand; stationary phase points correspond to "reduced" points representing periodic intervals of the summand at which an expansion to second order in the action is made producing quadratic Gauss sum contributions.

E Qutrit π/8 Gate
when the domain is increased from 3 to 3 2 ; we want to consider the sum over x q of U π/8 (x) = U π/8 (x q ) and reexpress it in terms of an exponential such that it is a sum over the larger space Z/3 2 Z of a polynomial with integer coefficients. Hence, we want to find α, β ∈ Z s.t.
Substituting in x q = 1 and x q = 3 2 − 1 = −1( mod 3 2 ) we find the system of equations 4 × (−2) mod 3 2 = α + β mod 3 2 , and 14 × (−2) mod 3 2 = α − β mod 3 2 , which has the solution α = 0 and β = 1. Therefore, Notice that the first line consists of a sum of an exponentiated polynomial with coefficients in Q while the second line is a sum over a larger space, but this time of an exponentiated polynomial with coefficients in Z. Thus, we have accomplished what we set out to do.
Since we will not be performing the sum above in isolation but with other exponentiated terms, so it is important to establish a slightly more general result: We evaluate the x 3q sum and, as sum away x 3p : We are interested in the magic state of this gate, which corresponds to it acting onĤ |0 = |p = 0 . The Wigner function of |p = 0 is ρ (x) ≡ 1 3 δ xp,0 . Hence, If we let S(x p , x q ) be the phase, we note that ∂S ∂xp = ∂S ∂xq = 0 mod 3 and ∇ 2 S ≡ H = 0 mod 3 0 ∀ x p , x q . Hence, evaluation at any phase space point, or linear combination thereof, requires summation over all three reduced phase space pointsx 2q ∈ Z/3Z since they are all critical points. Again, we note that Eq. 134 is 3-periodic in all of its arguments, and so its simplification by Theorem 1 is again a particularly simple special case where its sums can simply be restricted to be over Z/3Z with the appropriate power of 3 added to compensate.