Fast computation of spherical phase-space functions of quantum many-body states

Quantum devices are preparing increasingly more complex entangled quantum states. How can one effectively study these states in light of their increasing dimensions? Phase spaces such as Wigner functions provide a suitable framework. We focus on phase spaces for finite-dimensional quantum states of single qudits or permutationally symmetric states of multiple qubits. We present methods to efficiently compute the corresponding phase-space functions which are at least an order of magnitude faster than traditional methods. Quantum many-body states in much larger dimensions can now be effectively studied by experimentalists and theorists using these phase-space techniques.


I. INTRODUCTION
Current (and near-term) quantum devices are expected to prepare increasingly more complex entangled quantum states [1][2][3][4]. How can one effectively illustrate and analyze these states in light of their increasing dimensions? Phase spaces [5][6][7][8] such as Wigner functions have been widely used to meet this challenge. We will focus in this work on representing (finite-dimensional) quantum states of single qudits or permutationally symmetric states of multiple qubits using spherical phase spaces [9,10].
In this work, we consider spherical phase spaces of finite-dimensional quantum states and we develop a novel approach to efficiently compute these phase-space representations. For up to which dimensions can phase spaces be practically utilized? Our approach has a significant advantage in this regard as it allows for much larger dimensions to be addressed in a reasonable time frame. Therefore, phase-space descriptions of quantum manybody states are now feasible for dimensions which were beyond the reach of prior approaches. In summary, our results will enable practitioners and experimentalistbut also theorists-to visualise and study complex quantum states in considerably larger dimensions. This is accomplished by applying an efficiently computable Fourier series expansion and a fast Fourier transform (FFT) [55]. In particular, Fig. 1(a) compares the the run time of our Method C (as detailed in Sec. V) to the traditional Methods A and B (see Sec. III) and, indeed, our Method C is at least an order of magnitude faster. Moreover, Fig. 1(b) highlights that the rootmean-square error of certain test cases is comparable to machine precision for the considered dimensions and this suggests that our approach is numerically stable. We provide implementations in various programming environments (see Sec. V D and [56]), including C [57], MAT-LAB [58], Mathematica [59], and Python [60].
Our work has the following structure: We first discuss our motivation and highlight applications in Sec. II. Prior computational approaches to determine phase-space representations of finite-dimensional quantum systems are considered in Sec. III. In order to set the stage, we shortly recall the parity-operator description of spherical phase spaces which we have developed in [9]. Section V constitutes the main part of our manuscript where we develop our novel approach to efficiently compute spherical phase-space representations up to arbitrarily fine resolutions. We continue with a discussion of our results and further applications in Sec. VI, before we conclude. Important details are explained in appendices.

II. MOTIVATION AND APPLICATIONS
Various quantum-technology efforts (such as quantum computing or metrology) aim at creating large entangled multi-qubit states. Here, we focus in particular on the important class of states that are symmetric under permutations of qubits. These states include important families such as GHZ or squeezed states which are central in, e.g., quantum metrology [11] or entanglement veri-   [46][47][48][49]. Our Method C (Sec. V) combines spherical sampling techniques [50,51], explicit descriptions of rotation operators [52,53], Fourier series expansions, and fast Fourier transforms (FFT). The run times depend only on d and not the quantum state. (b) Root mean square (RMS) errors for certain quantum states [54] relative to their analytically known formula. Method C shows a high numerical precision comparable to machine precision.
fication [2,3]. They are also typically illustrated and analyzed in their phase-space representation (see, e.g., [3,11]) which can be naturally plotted on the surface of a sphere. This is reflected by the inherent symmetries and reduced degrees of freedom as compared to general multiqubit states. Before starting the technical discussion in Sec. III, we will now motivate our topic and highlight applications.
We first recall that permutationally symmetric states with N = 2J qubits can be mapped to states of a single spin J (or qudit with d = 2J+1) where J denotes a positiv integer or half-integer [10,[61][62][63][64][65]. Permutation symmetry appears in various applications including probe states in quantum metrology for optimal sensing, e.g., magnetic fields [11][12][13][14]. Permutationally symmetric qubit states can be efficiently reconstructed and are used for entanglement verification [2,3,11,[63][64][65]. We will illustrate a few practically relevant, high-dimensional examples for which traditional methods (see Sec. III) take an impractically large amount of time in order to determine the desired phase-space function. Further discussions and applications are deferred to Sec. VI.
The first example considers and highlights the Greenberger-Horne-Zeilinger (GHZ) state ( 0⟩ ⊗N + 1⟩ ⊗N ) √ 2 as the superposition of the all-zero and all-one state for N qubits which can be interpreted as the spin-up and spin-down states of a single qudit. Their high degree of entanglement supports the ultimate quantum precision in metrology, which is known as the Heisenberg limit [11]. GHZ states have been successfully created in numerous experiments with, e.g., trapped ions [28], superconducting qubits [3], and Rydberg atoms [2] for up to 20 qubits. Although phase-space functions of GHZ states can be analytically approximated for large dimensions [9,61], we are interested in computing them exactly within numerical precision and without relying on approximations. Figure 2(a) shows Wigner functions of GHZ states for an increasing number of qubits with N ∈ {8, 16, 32, 64}. Already the case N = 32 is currently beyond the exper-imental state of the art [2,3], but near-term quantum hardware are expected to deliver GHZ states of larger dimensions via, e.g., linear-depth quantum circuits [4].
We also consider so-called symmetric Dicke states [62] which are defined [10,63] as a superposition of all permutations of computational basis states with a fixed number of zeros and ones in a multi-qubit system. In particular, where the sum runs over all p = N n distinct permutations P k of the N qubits. These states are isomorphic to the single-qudit states Jm⟩ by mapping N to J 2 and m to (N 2 − n). We plot the Dicke state Jm⟩ with d = 2J+1 = 129 and m = 0 in Fig. 2(b). This corresponds to a highly entangled quantum state of 128 indistinguishable qubits where 64 qubits are in the 0⟩ state and 64 qubits are in the 1⟩ state (refer to Eq. (1)). One observes an axial symmetry (i.e. invariance under global Z rotations) and strong entanglement results in heavily oscillating Wigner functions in Fig. 2 Finally, squeezed states ξ⟩ ∶= exp[−iξ I 2 x ] 0⟩ ⊗N are obtained from the spin-up state of a single qudit or, equivalently, the all-zero state of N qubits under the influence of a squeezing interaction Hamiltonian I 2 x . The corresponding evolution time ξ is known as the squeezing angle [66] and I x is the x component of the total angular momentum operator, i.e., proportional to the sum of all Pauli σ x operators that act on different qubits. These states have been created in various experiments including Bose-Einstein condensates [17][18][19][20][21][22][23][24][25]67] for up to thousands of atoms. In such experiments, these finitedimensional squeezed states correspond to internal degrees of freedom (which we treat as an effective qudit) of fundamentally indistinguishable atoms. We plot their Wigner functions for the case of d = N + 1 = 500 and an increasing squeezing angle ξ in Fig. 2(c). For such large dimensions squeezed states with small squeezing x ] 0⟩ ⊗N with varying squeezing angle ξ and fixed d = N +1 = 500 in a plane for a small spherical subset with θ ≤ 0.05π (run time ≈ 1 min): Gaussian for ξ = 0 (left); squeezed Gaussians for ξ < 0.05. Larger ξ ≥ 0.05 lead to non-trivial and rapidly oscillating shapes which are nicely recovered, while analytical approximations fail in this regime. Red (dark gray) and green (light gray) for positive and negative values, respectively. The brightness indicates the absolute value of the function relative to its global maximum.
angles can be approximated well using the techniques described in [9,61]. In particular, the spin-up state 0⟩ ⊗N for ξ = 0 in Fig. 2(c) is a Gaussian-like function because the sphere can be approximated locally as a plane. For small squeezing angles, these states can be analytically approximated using star products [61,68]. Their phasespace representations are squeezed Gaussian functions which are very similar to the ones known in quantum optics [66,69]. This is illustrated in Fig. 2(c) where the aforementioned approximations apply to the cases ξ = 0, ξ = 0.003125, and ξ = 0.0125. For larger squeezing angles, Wigner functions will, however, deviate strongly from simple squeezed Gaussian states and non-trivial, heavily oscillating contributions become dominant as is shown in Fig. 2(c) for ξ = 0.05 and ξ = 0.2. This motivates our numerical approach to exactly determine phase-space functions for large spin-like systems (and permutationally symmetric multi-qubit states) where analytical approximations do usually fail.
(2). The expansion coefficients c jm ∶= Tr [ρ T † jm ] are computed from the density matrix ρ and the tensoroperator coefficients T jm [73][74][75][76]. The matrix elements Equation (2) describes the standard approach for numerically computing spherical phase-space functions. In a first step, it relies on efficient approaches to calculate Clebsch-Gordan coefficients. The calculation of the expansion coefficients c jm is, however, computationally expensive for large dimensions d = 2J+1 ≫ 1. In particular, one needs to determine O(d 2 ) distinct tensor-operators T jm and their matrix entries. Appendix A clarifies that O(d 3 ) Clebsch-Gordan coefficients have to be calculated which dominates the run time for computing all of the O(d 2 ) expansion coefficients c jm in Eq. (2).
Two different approaches to calculate Clebsch-Gordan coefficients result in two different methods (Method A and B) to the coefficients c jm . Method A uses the builtin Mathematica [59] function that performs arbitraryprecision integer arithmetic. In Method B, the run time can be significantly reduced by numerically computing Clebsch-Gordan coefficients using a FORTRAN [82] implementation [49] of a recursive algorithm [46][47][48]. Methods A and B are compared in Fig. (3). For Method A (B), all tensor operators for certain dimensions d ≤ 300 (d ≤ 500) have been determined and we estimate a complexity O(d 4 ) in this range.
After the expansion coefficients c jm have been obtained, the phase-space function F ρ (θ, φ, s) is spherically sampled in a second step by applying a fast spherical harmonics transform which might rely on equiangular samples or Gauss-Legendre grids. The second step requires a practically and asymptotically negligible time of O(d 3 ) when compared to the first step. Spherical harmonics transforms are widely used in various scientific contexts and efficient implementations are available [50,[83][84][85][86].

IV. THE PARITY-OPERATOR DESCRIPTION OF SPHERICAL PHASE SPACES
We recall the parity-operator description of spherical phase spaces developed in [9] in order to develop faster methods to compute spherical phase-space functions in Sec. V. We keep the notation introduced in Sec. III and specify the rotation operator as R(θ, φ) ∶= e iφJz e iθJy , where J z and J y are components of the angular momentum operator [88]. Building on [71,78,79,89,90], the s-parametrized phase-space functions are defined in [9] as expectation values of rotated parity operators M s by This extends work [91][92][93][94][95] on rotated parity operators to all s-parametrized phase spaces. The parity operator is defined by its expansion into diagonal tensor operators T j0 of order zero. The corresponding matrix elements are given by for j ∈ N∪{0} and m, m ′ ∈ {−J, . . . , J}. Equation (2) could be recovered by applying the rotation operators to the tensor operators in Eq.
. For an increasing spin number J, spherical phase spaces converge to their infinite-dimensional counterparts while rotations transform into translations along the tangent of a sphere [9,10,61,96]. While we focus here on single qudits (and permutationally symmetric quantum states of multiple qubits), generalizations of the parityoperator approach to arbitrary coupled quantum states are also available [93,94,[97][98][99].

V. EFFICIENT COMPUTATION OF SPHERICAL PHASE-SPACE FUNCTIONS
We develop now our main results on efficiently computing spherical phase-space functions. Section V A presents a first approach using parity operators (see Sec. IV), an explicit form for rotation operators, and a spherical sampling strategy. This does-by itself-not lead to an effective approach. But it provides the necessary ingredients to specify spherical phase-space functions as a finite Fourier series in Sec. V B which includes our efficient algorithm for the corresponding Fourier coefficients. A fast Fourier transform is then applied as detailed in Sec. V C to recover an equiangular spherical sampling of the phase-space function. Finally, we discuss implementations of our efficient algorithms in Sec. V D.

A. A first approach via parity operators, matrix entries of rotations, and spherical sampling
Equation (4) can be directly applied to calculate phasespace functions as expectation values of rotated parity operators. The parity operators are determined by Eq. (5) and the matrix entries of the rotation operator [80] are analytically given as Wigner-D functions (which are widely available in software environments such as Mathematica). We also use results of [52,53] to compute the matrix entries of the rotation operator using fast Fourier transforms (see Appendix B). The phase-space function is then computed as the trace of the matrix product of the operators in (4).
to their spherical harmonics decompositions, we can apply spherical sampling schemes with a discretized grid of spherical angles (θ k , φ ). One can uniquely represent a phase-space function by sampling on an equiangular grid with n 2 ≥ (4J+2) 2 = (2d) 2 rotation angles [50,51]. One then evaluates Eq. (4) at all angles in Eq. (6) to obtain a equiangular spherical sampling of the phase-space function. However, this first approach requires matrix multiplications for each of the O(d 2 ) spherical angles. This leads to inefficiencies and an overall run time of O(d m ), where 4.2 ⪅ m ≤ 5 depending on the efficiency of the matrix-multiplication algorithm (and m = 5 corresponds to a naive implementation) [100]. More effective methods are presented in Sec. V B. The presented approach can be combined with the algorithm of [50,51] to recover the spherical-harmonics expansion coefficients c jm in Eq. (2).

B. Efficient algorithms for the Fourier coefficients
We now expand on the approach in Sec. V A by exploiting the structure of the rotated parity operators and by analytically evaluating the matrix products in Eq. (4). This facilitates a novel computational scheme for computing the Fourier expansion of spherical phase-space functions which significantly differs from the methods in [52,53]. We begin by computing the Fourier expansion coefficients of the rotation operators R(θ, φ). Recall that any (unitary) matrix can be written in terms of its spec-tral resolution which also holds for As detailed in Appendix B, A and B m are projection operators that project onto the eigenvectors of the spin operators J y and J z , respectively. The dependence on the rotation angles has been completely absorbed into the Fourier components e i θ e imφ .
We can now analytically evaluate the trace of matrix products in Eq. (4) and we prove in Appendix C that the phase-space function can be decomposed into a finite, band-limited Fourier series. The Fourier expansion coefficients F m implicitly depend on the density matrix ρ and the parity operator M s (as well as s) and they can be obtained from ρ via a linear transformation: where −2J ≤ , m ≤ 2J and ρ m1,m2 ∶= ⟨Jm 1 ρ Jm 2 ⟩ are the density-matrix entries in the standard qudit basis.
A proof of Result 1 is given in Appendix C. The transformation matrices K ∈ C d×d implicitly depend on the parity operator M s (and s). They can be efficiently calculated as a finite sum (see Appendix D) Here,M s denotes the parity operator M s transformed into the eigenbasis of the operator J y , and U ν ⟩ are the eigenvectors of J y , such that J y U ν ⟩ ∶= ν U ν ⟩. The matrix entries ofM s are therefore given as [M s ] ab = ⟨U a M s U b ⟩. Result 1 leads to two different algorithms to compute the Fourier coefficients in Eq. (8) (as detailed in Appendix D). These algorithms are then combined with a fast Fourier transform (which has a much smaller run time) in order to effectively compute an equiangular spherical sampling of the spherical phase-space function (as discussed in Sec. V C). The first algorithm to compute the Fourier coefficients is denoted as Method C: The transformation matrix K is computed for a fixed via Eq. (10)  The run time of a C implementation of Method C is compared in Fig. 3 to the traditional Methods A and B from Sec. III. We empirically observe an asymptotic scaling of O(d 4 ) for all three methods and d ≤ 500, which is visible as near-parallel lines in the log-log plot of Fig. 3. However, Method C is evidently much faster. Figure 1 (a) shows the relative runtimes of Methods A and B compared to Method C highlighting that Method C is at least an order of magnitude faster. Consequently, Method C can be used for much larger dimensions.
The second algorithm to compute the Fourier coefficients in Eq. (8)  This results in a significantly faster implementation (see Fig. 3) which also suggests a better asymptotic scaling (with a smaller slope in Fig. 3). The disk storage and RAM requirements for Methods C and D are detailed in Table I while assuming double precision. Method D is preferable (at least) for dimensions d ≤ 500 as it significantly reduces the run time with a reasonable amount of disk storage. For larger dimensions, one has to balance speed with storage requirements.

C. Spherical sampling of the phase-space function via a fast Fourier transform
We now utilize the Fourier series from Sec. V B to obtain an equiangular spherical sampling of a phasespace function by applying a fast Fourier transform. We start with the (4J+1) × (4J+1) Fourier coefficients F m   (d)) time complexity and results in a grid with (4J+1) × (4J+1) spherical samples of the phase-space function. But this is only the coarsest grid possible for a complete reconstruction (refer to Eq. 6) and finer girds can correct for non-uniformities and lead to smoother spherical representations.
In order to obtain a finer grid, it is preferable to add zero padding to the Fourier coefficients which results in a n × 2n coefficient array with additional zeros where n ≥ 4J+2. Many FFT implementations are optimized for n being a power of two. After applying the FFT, one essentially obtains two copies of the phase-space function as θ varies over 0 ≤ θ < 2π in the result (while the phasespace function is only defined for 0 ≤ θ < π). However, by straightforwardly discarding the redundant half one recovers the desired n × n sampling of the phase-space function.
Note that this equiangular sampling is compatible with (equiangular) spherical harmonics transforms (see Sec. (6) and, e.g., [50,51,84]) that could be used to compute the coefficients c jm in Eq. (2). We also remark that performing fast Fourier transforms is usually preferable to fast spherical transforms (which are used in Methods A and B). This is particularly relevant when one aims at sampling phase-space functions for a fixed dimension d to an arbitrarily high resolution n. The two-dimensional FFT takes O(n 2 log 2 (n)) time. Practical spherical harmonics transforms have, however, a time complexity between O(n 5 2 log(n)) and O(n 3 ) depending on the implementation [50,[83][84][85][86] and asymptotically faster implementations might introduce numerical errors and only become superior for very fine resolutions [85].

D. Implementations of our algorithms
We have made implementations of our algorithms for computing spherical samplings of phase-space functions freely available [56]. The algorithm for precomputing the coefficients K in Eq. (10) for a fixed dimension d has been implemented in C without any external dependencies. For convenience, we provide a program (with external dependencies as LAPACK [101]) to precompute the parity operators [M s ] ξξ and eigenvectors U ν ⟩ (Sec. B 2), even though their computation time and storage requirements are negligible (see Table I). We currently interface with the precomputed data for d ≤ 500. Using the precomputed data, implementations of Method D with suitable zero padding (Sec. V C) are available for C, MAT-LAB, Mathematica, and Python [102].

VI. DISCUSSION
Traditional approaches to efficiently compute spherical phase-space functions rely heavily on expensive evaluations of Clebsch-Gordan coefficients and use spherical harmonics transformations (see Sec. III). We provide much faster algorithms by going beyond these techniques and by applying a suitable Fourier expansion and a fast Fourier transform. This leads to the two variants (Method C and D) which involve different time-memory tradeoffs. Method C calculates the transformation matrices K on-the-fly and they are then employed to spherically sample the phase-space function in O(d 4 ) time. Method D precomputes the transformation matrices K and stores them using O(d 3 ) disk space. The stored transformation matrices enable us to spherically sample the phase-space functions in O(d 3 ) time. We have implemented our algorithms in various programming environments such as C, MATLAB, Mathematica, and Python [56].
We also remark that our C implementation can be further optimized, e.g., with regard to memory handling and loops. The overall run time of the discussed algorithms could be reduced by truncating spherical-harmonics or Fourier coefficients which could be motivated by prior knowledge or symmetry considerations. In addition, the disk storage of Method D can be optimized to O(d) if the summation in Eq. (8) can be restricted to Fourier coefficients F m with , m ≤ t for some suitable constant t. But this might not be a good approximation for general quantum states and we are focussing on computing phase-space function exactly up to numerical precision.
We finally discuss how our results could be applied to compute analytical derivatives with respect to spherical rotation angles. Following Sec. V and Result 1, one obtains the Fourier coefficients F m and this representation helps us to compute derivatives analytically by multiplying the coefficients F m with i × (or i × m): These derivatives are particularly relevant for the computation of star products of phase-space functions (see [61]). This can be extended to analytical gradients which enables us to search for local extrema of phasespace functions (e.g., minima of locally negative regions) via gradient descent optimizations.

VII. CONCLUSION
In this work, we have considered spherical phase spaces of large quantum states and have provided effective computational methods for them. Our methods allow now for much larger dimensions than before. Going beyond approaches using tensor-operator decompositions and spherical-harmonics transforms, we can directly harness the efficiency of the fast Fourier transform applied to an efficiently computable Fourier series expansion. Our C implementation [56] is at least an order of magnitude faster than prior implementations when compared for up to dimension 500 (or up to 499 qubits in permutationally symmetric states). Our data also suggest an asymptotic speed-up by utilizing suitable precomputations.
The presented computational methods for phase spaces of single-qudit and permutation-symmetric multi-qubit states have applications to many-body physics, quantum metrology, and entanglement validation. We have illustrated many-body examples in Sec. II some of which are pursued in current quantum hardware. Our results will enable both theoreticians and experimentalists to more effectively work with phase-space representations in order to study high-dimensional quantum effects. This will help to guide future experimental advancements in generating complex quantum states of high fidelities [1][2][3][4].  should be directly available in memory and the overall computation time of this approach is therefore dominated by evaluating the Clebsch-Gordan coefficients. We expect that computing a single one of them requires O(d n ) time with n > 0 and based on our numerical computations in Fig. 3 we speculate that n ≈ 1.

Appendix B: Fourier series representation of the rotation operator
We now establish how the rotation operator in Eq. (4) can be decomposed into a Fourier series. This step is crucial for deriving our Result 1, which finally allows us to efficiently decompose a phase-space function into Fourier components.
Recall that the rotation operator defined in Eq. (4) is parametrized in terms of Euler angles as R(θ, φ) = e iφJz e iθJy via the spin operators J y and J z . These spin operators are defined via their commutation relations [J j , J k ] = i ∑ jk J for j, k, = x, y, z and jk is the Levi-Civita symbol, refer to, e.g., [88,103]. For an N -qubit system these are proportional to sums of Pauli acting on individual qubits and j ∈ {x, y, z}. These operators are unitarily equivalent and have the eigenvalues m ∈ {−J, −J + 1, . . . , J} due to the eigenvalue equation Note that in an N qubit system 2J = N . Here we denote eigenvectors of the J y operator as U m ⟩ and recall the orthogonality condition ⟨U m U n ⟩ = ⟨Jm Jn⟩ = δ mn . The spectral resolution of these spin operators is obtained in terms of the rank-1 projectors U m ⟩⟨U m =∶ A m and Jm⟩⟨Jm =∶ B m as It immediately follows that rotation operators decompose into the following sum of rank-one projectors Note that the dependency on the rotation angles θ and φ is now completely absorbed by the Fourier components e iθm and e iφm . The rank-1 matrices A m and B m are projections onto the eigenvectors of the spin operator J y from Eq. (B2) and we define their matrix elements as The explicit form of the Fourier coefficients [A m ] m1m2 was derived analytically in [52] as with summation bounds a = max (0, m 2 − m 1 ) and b = min (J − m 1 , J + m 2 ). The explicit form of the coefficients appearing in the above summation are with summation bounds c = max (0, −J + m + λ) and here d = min (λ, J + m) and (. . . )! denotes the factorial function while λ are the binomial coefficients.

Numerical computation of the eigenvectors
A simple and efficient way for numerically evaluating the coefficients [A m ] m1m2 in Eq. (B5) was proposed in [53]. This approach first computes the eigenvectors U m ⟩ from Eq. (B1) by numerically diagonalizing the spin operator J y . One then obtains the numerical representation of the eigenvectors U m ⟩ that define the rank-1 projector A m = U m ⟩⟨U m . Its matrix elements can then be obtained straightforwardly as products of vector entries of eigenvectors of J y from Eq. (B1) and here [. . . ] * denotes complex conjugation. The matrix J y can be diagonalised to numerical precision (it is tridiagonal and Hermitian) which provides a highprecision numerical representations of [A m ] m1m2 . This has been demonstrated in [53] using the ZHBEV diagonalisation routine of the software package LAPACK [101]. We use this approach in this work for numerically computing eigenvectors.

Appendix C: Derivation of Result 1
Substituting the expansion of rotation operators from Eq. (B3) into our definition of phase spaces in Eq. (4) and using that the rank-one projectors A m and B m are self adjoint we obtain This is a Fourier series decomposition of the phase-space functions. It is our aim now to express its Fourier coefficients explicitly. In particular, one can rearrange the terms in the trace and obtain where the first term in the trace is simply a projection of the density matrix onto a single matrix element in the z basis as B λ ρ B κ = Jλ⟩⟨Jκ ρ λκ . Here, matrix elements of the density operator are denoted as ρ λκ ∶= ⟨Jλ ρ Jκ⟩ assuming the standard z basis.
We now explicitly express this phase-space function as a Fourier series and denote its expansion coefficients as F m via  Here we have introduced the set of matrices K which simply multiply the density matrix element-wise and we define their explicit form as a summation over the matrix products Note that the Fourier coefficients F m depend both on the density operator ρ and on the parity operator M s , and implicitly on the eigenvectors of J y . We have introduced the matrices K , which completely determine the dependence on the parity operator and on the eigenvectors of J y . These matrices can be precomputed and stored or computed on-the-fly. The Fourier coefficients can then be completely determined via the efficient summation We now fix and evaluate Eq. (D1) for this fixed . We compute the matrix K element-wise as [K ] ab using the explicit expression [ U ν ⟩⟨U ν+ ] ab = [U ] νa ([U ] ν+ ,b ]) * , where * denotes complex conjugation. Computing such a matrix K in Eq. (D1) requires O(d 3 ) time for a fixed . We therefore conclude that computing every coefficient matrix K with ∈ {−2J, . . . 2J} requires O(d 4 ) time.
After computing K for a fixed , one can proceed according to two distinct strategies, which we refer to as Method C and D in the main text. In case of Method D, we store the matrix K and repeat this procedure for each ∈ {−2J, . . . 2J}. This requires O(d 3 ) disk storage space. These precomputed matrices can be used later in Result 1 for computing phase spaces in O(d 3 ) time which requires only O(d 2 ) memory, i.e., for ρ, U andM s , and one only reads in a single matrix K at a time. In case of Method C, we compute K for a fixed , and use it immediately for evaluating the summation in Result 1 for a fixed . We can then repeat this procedure for each ∈ {−2J, . . . 2J}. Therefore, Method C does not require disk storage space for the matrices K , but allows for calculating phase-spaces via Result 1 in O(d 4 ) time and similarly using O(d 2 ) memory.