Overview of the phase space formulation of quantum mechanics with application to quantum technologies

The phase-space formulation of quantum mechanics has recently seen increased use in testing quantum technologies, including metho ds of tomography for state verification and device validation. Here, an overview of quantum mechanics in phase space is presented. The formulation to generate a generalized phase-space function for any arbitrary quantum system is shown, such as the Wigner and Weyl functions along with the asso ciated Q and P functions. Examples of how these different formulations have b een used in quantum technologies are provided, with a focus on discrete quantum systems, qubits in particular. Also provided are some results that, to the authors' knowledge, have not been published elsewhere. These results provide insight into the relation between different representations of phase space and how the phase-space representation is a powerful tool in understanding quantum information and quantum technologies.

states in QED experiments [41] , and the rotation of qubit before taking a generalized parity measurement in Reference [42]. More examples of such experiments are given within this review.
In this review, we begin by presenting the formulation of quasiporbability distribution functions for systems of position and momentum, where we will then generalize this formulation for use in any arbitrary system in Section 2. This will provide an understanding of the general framework needed to understand how to map a given operator to its representation in phase space, and how such a mapping is bijective and informationally complete. Following this in Section 3, we will use this generalized formulation to introduce the phase-space formulation of various finite-dimensional quantum systems. This will include an in-depth discussion on different representations of qubits in phase space and how different representations are related; we will also present multiple ways one can represent multi-qubit states in phase space dependent on different group structures. We will then describe how these methods can be used in practice for use in quantum technologies in Section 4; where we will translate some important figures of merit for quantum computation into the language of phase space, as well discussing the importance of negative values in the Wigner function. This will be followed by examples of experiments performed on qubits where the phase space distributions were directly measured.

Formulation of phase space
In this section we will lay the groundwork for how one can generate a phase-space distribution for any arbitrary operator, where we will provide a general framework for any group structure in Section 2.3. We will begin with the original formulation in position-momentum space, within the Heisenberg-Weyl group (HW) where we will provide the general formula to create any quasiprobability distribution function and any characteristic function. This will be followed by looking into the dynamics of such states and how representing quantum mechanics in phase space allows one to see to what extent one can treat dynamics classically. This section will end with a general framework and a summery of the results in Reference [38].

Systems of position and momentum
Wigner's original formulation was generated in terms of a wavefunction, ψ(q), and is in essence the Fourier transform of a correlation function, where for a function defined in terms of position and momentum where we will assume that = 1. We note that in the original construction, the Wigner function was presented to describe the position-momentum phase space for many particles, here we will restrict the discussion to one concomitant pair of degrees of freedom where we will discuss the composite case later in this section. This formulation of the Wigner function was then soon generalized to apply to any arbitrary operator, A, as an expectation value of a kernel, Π(α), such that [3,4] where we now define the Wigner function in terms of the complex variable α = (q + ip) / √ 2 to fit inline with the quantum optics literature.
The kernel in Equation (2) is fully constructed as a displaced parity operator which consists of the parity operator Π = 2e ia † a , The displacement operator displaces any arbitrary state around phase space and is central to the definition of a coherent state, which is generated by displacing the vacuum state, |0 , such that where |α is an arbitrary coherent state. Note that there is a factor of 2 in Equation (4), and subsequently in Equation (3), this is chosen for two reasons. First, it aligns with the formulation of the Wigner function that was introduced by Glauber when writing it in a displaced parity formalism. But more importantly, it lines up more nicely with the definition of the Wigner function for discrete systems. Later, we will introduce a condition that requires the trace of the kernel to be one, which requires Equation (4) to have the factor of 2 out the front to hold.
The Wigner function is just one of a whole class of potential quasiprobability distribution functions for quantum systems, what makes it a preferable choice is that the kernel to transform to the Wigner function is the same as the kernel used to transform back to the operator, which is not always the case, such that showing that the Wigner function is an informationally complete transformation of any arbitrary operator. Another important advantage is that the marginals of the Wigner function are the probability distribution of the wavefunction. That is for the Wigner function of a state ψ(q) integrating over all the values of p to produce the probability for ψ(q), likewise for ψ(p) when integrating over q.
To represent quantum mechanics in phase space, some properties one would usually desire from a probability distribution must be sacrificed -hence the prefix 'quasi' in quasiprobability distribution functions. For the Wigner function this sacrifice is on positive-definiteness, a more in-depth conversation on the negativity of the Wigner function will be given in Section 4.2. If however one doesn't want to sacrifice a non-negative probability distribution function, the alternative is to consider the Husimi Q function [6] , where instead the sacrifice comes in the form of losing the correct marginals.
The Q function is generated in terms of coherent states, where where the last form of the Q function can be thought of as an expectation value of a general coherent state, where instead of the kernel being a displaced parity operator it is a displacement of the projection of the vacuum state. Another drawback to the Q function for some is that transforming back to the operator will require a different kernel, this is the kernel that transforms an arbitrary operator to the Glauber-Sudarshan P function [9,8] . Likewise, the kernel to transform back from the P function to the operator is the kernel for the Q function -this is the usual definition of the P function where P A (α) is the P function for an arbitrary operator, A. The P function seems to sacrifice more, as it doesn't have the correct marginals and it is not non-negative, in fact coherent states in the P functions are singular, which may be an attractive analogy when considering coherent states as so-called 'classical stats', although we must note that coherent states are still fundamentally quantum. Despite the apparent drawbacks, there have been a set of key results by analysis of the P function with Gaussian smoothing, see for example Reference [43] -this will also be discussed further in Section 4 Apart from the Wigner, Q and P functions, there are further choices of quasi-probability distribution function for quantum systems. In order to completely describe the full range of possible phase-space functions, including the appropriate kernels to transform both to and from the given function, it is first necessary to introduce the quantum mechanical characteristic functions. The characteristic function is a useful tool in mathematical analysis, and is the Fourier transform of the probability distribution function. In the same way, the characteristic functions in quantum mechanics are Fourier transforms of the quasi-probability distribution functions. The Fourier transform of the Wigner function is also known as the Weyl function and is defined such that relate the Wigner and Weyl functions.
The superscript in brackets (0) in Equation (11) is a specific case of the characteristic function, χ Note that when s = 0, D 0 (α) is the standard displacement operator from Equation (5). When s = 0 we say that the displacement operator is symmetically, or Weyl, ordered. Alternatively, when s = −1 or s = 1 the corresponding displacement operators, and therefore characteristic functions, are anti-normal and normal respectively; this can be realized by showing the operators explicitly D −1 (α) = e −α * a e αa † and D 1 (α) = e αa † e −α * a , demonstrating the opposite order of the arguments in the exponents. Note that, like the quasi-probability distribution functions, the characteristic functions are informationally complete, such that where the reverse transform kernel is the Hermitian conjugate of the generalized displacement operator with minus the s value of the function. The Q and P functions are then the Fourier transforms of the anti-normally ordered and normally ordered characteristic functions respectively, where (α) exp (αα * −α * α) dα, and, P (α) = 1 π ∞ −∞ χ (1) A (α) exp (αα * −α * α) dα. (16) More generally we can say that they are the normally and anti-normally quasi-probability distribution, where we can define a general s-values quasi-probability distribution function F (s) which hold for any value of s, generating an array of quasi-probability distribution functions, however only s = −1, 0, 1 will be considered here. Further, by using the definition of the characteristic function in Equation (13) and substituting it into Equation (17), we can see that by bringing the integral inside the trace operation we can define the generalized kernel for any quasi-probability distribution function as the Fourier transform of the displacement operator Where a general function in HW phase space is defined as For the reverse transform, to take a phase-space distribution to an operator, one can then perform providing the transforms between an arbitrary density operator and an arbitrary s-valued quasiprobability distribution function. To complete the picture, all that is needed is a discussion of the dynamics on quantum systems in phase space.

Evolution in phase space
Dynamics of quantum systems are particularly appealing in the phase space formulation. To understand why this is the case we return to classical Hamiltonian physics and understand the phase space approach in terms of the Wigner function from this perspective. The classical dynamics of any distribution A for a given Hamiltonian H is given by where the Poisson bracket defines the way that a quantities flow in phase space; as the Poisson bracket is an example of a Lie bracket, the fact that this represents a specific kind of infinitesimal transform for classical dynamics is natural. The Poisson bracket thus specifies the way in which the state of the system transforms incrementally in time for a given Hamiltonian as a flow in phase space. For example, in a time increment dt the coordinates and we can see how the Hamiltonian defines a "flow" in phase space (the above picture illustrates this in one-dimension). If we change the definition of the Poisson bracket this would result in different physics and the phase space formulation of quantum mechanics follows exactly this line of reasoning. In this way we will see that in redefining the Poisson bracket Quantum mechanics can be seen as a continuous deformation of this underlying "flow" algebra. Before proceeding we first note that if we set A equal to the probability density function ρ(t) then we observe two important corollaries. The first is that conservation of probability requires dρ dt = 0 (so it is a constant, but not necessarily an integral, of the motion) and the second is that {ρ, H } can be thought of as the divergence of the probability current. This is expressed as in the Liouville "conservation of probability" equation (which has obvious parallels with the von-Neumann equation for density operators): When Dirac motivated the Heisenberg picture of quantum mechanics he did so by arguing that equations of the above form and the Lie algebraic structure of the Poisson bracket be both retained, if quantum dynamics is to be defined in terms of some infinitesimal transform this is a sensible requirement. Dirac's argument then was to introduce other mathematical structures in which to achieve this end -operators in a vector space. The alternative we take here is the Deformation Quantization approach that leaves all quantities in exactly the same form as in classical Hamiltonian physics but to introduce quantum phenomena by deforming the phase space directly. In this view, function multiplication is replaced by a 'star' product, such that lim →0 f g = f g and that we introduce an alternative to the Poisson bracket -termed the Moyal bracket that satisfies the property This last expression guarantees that any such formulation of 'quantum mechanics' will smoothly reproduce classical dynamics in the limit → 0 and for this reason can be viewed as a deformation of classical physics. For this reason we briefly reintroduce in this section for comparison with classical mechanics, we will then return to = 1 in the next section. Understanding the Wigner function as the exact analog of the probability density function -and as a quantum constant of the motion -means that conservation of probability requires dW dt = 0, we also note that {{W, H }} can be thought of as the divergence of this generalized probability current. The Louiville equation is then replaced by the more general equation: We see that stationary states are then those for which {{H , W }} = 0 which lead to left and rightgenvalue equations that are the phase space analog of the time-independent Schrödinger equation (whose solutions will be integrals of the motion as the phase-space analogue of energy eigenstates). Although this is not a foundations paper it is worth noting: (i) the quantisation of phase space resolves issues around operator ordering in convectional quantum mechanics (see Groenewalds theorem for example [3] ); (ii) that this approach may be of particular value when looking to unify quantum mechanics and gravity as time is treated on the same footing as all other coordinates and the geometric view of introducing curvature through the metric tensor is perfectly natural in a phase space formulation. The star-product is usually written, for canonical position and momentum, in its differential form as Importantly for our discussion the star product can also be written in a convolution-integral form: f (q, p) g(q, p) = 1 π 2 2 dq dp dq dp f (q + q , p + p ) g(q + q , p + p ) exp 2i (q p − q p ) .
which has the advantage that it both highlights that quantum effects are non-local and, from a computational perspective, numerical integration can be less prone to error than differentiation (although more costly). We now can use a specific example to a result by Brif and Mann [33] to connect this form to a more general expression f (q, p) g(q, p) = dq dp dq dp f (q , p ) g(q , p ) Tr [Π(q, p)Π(q , p )Π(q , p )] .
where Π(q, p) is the usual displaced parity operator.
To summarise the quantum dynamics of any system or composite of systems is given by where we see this as an equivalent view to classical mechanics with a deformation of function multiplication that in the limit → 0 will, if the limit exits, reduce these dynamics to classical physics. This can be viewed as moving from a non-local to local theory as W (Ω) H (Ω) → W (Ω)H (Ω) as → 0. As a final note the classical dynamics of particles should even be recoverable in this limit as → 0 the Wigner function of, e.g., coherent states tends to Dirac delta functions centred at one point in phase space ρ(q i , p i ) = i δ(q i − Q i , p i − P i ) and the Liouville equation for point particles will recover that classical point trajectories in an analogous way to solutions of the Kontsevich equation used in plasma physics. [44] As Moyal brackets tend to Poisson brackets as → 0 it is natural to see that the Kontsevich equation is recoverable in this limit.

Quantum mechanics as a statistical theory
We will now generalize the discussion of phase space for use for any group structure. This will include us generalizing important equations from Section 2.1 and providing the generalization for the evolution in Equation (34).
To generalize the formulation of phase-space functions, a suitable candidate for the given group structure needs to be chosen to replace the Heisenberg-Weyl displacement operator. [38] Such a generalized displacement operator, D(Ω) then displaces the corresponding choice of vacuum state for the given system, |vac , creating the generalized coherent state in that system [31] |Ω = D(Ω) |vac .
For examples of such operators and the generation of generalized coherent states and the generation of generalized displacement operators for different systems, see for example References [31,45,46,47,36,48,49]. Note that when generating quasiprobability distribution functions, any suitable generalized displacement operator can produce an informationally complete function, however when considering just the Weyl characteristic function much more care needs to be taken in considering the correct operator, [38] since it is crucial to determine a generalized displacement operator that generates an informationally complete function, where the inverse transform to the operator exists, such that where dΩ is the volume-normalized differential element or Haar measure for the given system. Once a generalized displacement operator has been chosen, the generation of the Q and P functions simply requires the generalized coherent state from Equation (35), where however the construction of the Wigner function is somewhat more involved.
To generate the Wigner function we now need to construct the corresponding generalized parity operator, Π, to create the kernel Π(Ω), yielding the Wigner function S-W. 2 Reality: W A (Ω) is always real valued when A is Hermitian, therefore Π(Ω) must be Hermitian. From the structure of the kernel in Equation (38), this also means that the generalized parity operator, Π, must also be Hermitian.
is integrated over the full volume Ω, and dΩ is the volume normalized differential element, see Appendices in Reference [36] or Reference [38] for more details. Also note that, given the kernel for a system, the generalised parity can be recovered by taking this is a helpful result in future calculations. [38,12,13] We can therefore say in general that any quasiprobability distribution function is defined where is the generalized kernel for any s-valued phase-space function. As is demonstrated from Equation (39), the Wigner function is an informationally complete representation of the state, i. e. it contains all the same information as the density operator as one can transform between the two representations without loss of information. The same can be said of the Q and P functions, where the operator can be returned from the Q function by integrating over the kernel for the P function, likewise the operator can be returned by integrating the P function over the Q function kernel, explicitly which is the generalization of Equation (20). It is now natural to ask whether there is a simple comparison between the different phase-space functions without first going to the density operator. Indeed this can be done by way of a generalized Fourier transform. [33,28,38,50] We can simply bypass the calculation of the density operator by substituting Equation (43) for one value of s into Equation (41) for a different s value, yielding A (Ω) F(Ω , s 2 ; Ω, s 1 ) dΩ (44) where F(Ω , s 2 ; Ω, is the generalized Fourier operator. This can further be generalized to transform between any phase-space function, including the Weyl function, by substituting in Equation (36) instead. We can further extend the idea behind this Fourier kernel to perform a convolution between any two phase-space functions to create a third function, [15,34,33,28,51,50] such that where is the convolution kernel. Stratonovich showed in Reference [5] that the kernels can be combined for the individual systems to create a kernel for a composite system. The procedure to do so is straight-forward, it simply requires one to take the Kronecker product of the individual systems for n subsystems, where Ω = {Ω 1 , Ω 2 , ..., Ω n } and Ω i are the degrees of freedom for the individual systems. Similarly Π i (Ω i ) is the kernel for each individual system. The full composite Wigner function is then generated by using Equation (48) in Equation (38). Transforming back to the operator from the composite Wigner function requires a simple modification of Equation (39), where the new volume normalized differential operator is more details on this construction can be found in Reference [38]. We now have all the ingredients necessary to generate a quasi-distribution functions for any system, or any composition of individual systems, to represent quantum technologies in phase space.
3 Phase-space formulation of finite quantum systems We will now turn our attention to the generation of phase space functions for finite quantum systems. It will be shown how the general framework from Section 2.3 can be applied to specific finite-dimensional quantum systems to represent them in phase space. Where consideration of such structures is important in considering the use of phase-space methods in quantum technologies. As such, we will show the important case of representing qubits in phase space, providing a link between the two most well-known formulations of qubits in phase space. This will be followed by looking at how more general qudits and collections of qubits can be represented in phase-space in different situations. Given the phase-space representation of all these systems, along with the continuous-variable construction for optical quantum systems, we will then be able to consider many approaches to quantum technologies in phase space.

Qubits in phase space
The most important system to introduce now for use in quantum technologies is the construction of phasespace functions for qubits. The qubit Wigner function has two main formulations, there is the original construction by Stratonovich that considers the qubit over the Bloch sphere, with continuous degrees of freedom. The second formulation was generalized by Wootters, and considers the qubit states on a discrete toroidal lattice, with discrete degrees of freedom. Both formulations will be considered here, as each has had important results in understanding the properties of quantum mechanics in a quantum technology setting. These two formulations are also specific examples of a general formulation of qudits in phase space. It is also interesting to note that the discretelattice qubit Wigner function from Wootters can be considered a sub-quasiprobability distribution of the continuous-spherical Wigner function. This relation will be shown later, but this relation does not hold for general qudit states.
We will also consider the generation of the Q and P functions for qubits, which are already well-known in the Stratonovich formulation, however to our knowledge does not exist in the discrete case.

The Stratonovich kernel
Following the steps from Section 2, a suitable generalized displacement operator and generalized parity operator is necessary for the construction of a Wigner function for qubits. We will start by considering an appropriate generalized displacement operator. Since we are now concerned with distributions on the sphere, the generalized displacement operator is a rotation operator. An example of the change in geometry and in displacement can be seen in Figure 1 (b). A complicating factor in choosing an operator is that rotation around a sphere isn't a unique procedure, there are many ways in which one can perform such an operation. Out of the options, we highlight two choices that are heavily used in the appropriate literature and give equivalent results upon rotation of the generalized parity operator. These are which are the standard Euler angles, where we introduce the general notation U M N that in this case is set to SU(N ), M = 2j, for the quantum number j. The second operator R(ξ), where ξ = θ exp(−iφ)/2 and σ ± = (σ x ± iσ y )/2, is a generalization of the Heisenberg-Weyl displacement for SU (2), and can be found in much of the spin-coherent state literature -for example see References [31,29] -and can be related to the Euler angles by [38] R Although it is much more tempting to have a preference towards this rotation operator, due to the similarity to the Heisenberg-Weyl case, this rotation operator generates a Weyl function that is not informationally We therefore need to construct an informationally complete Weyl function by including all three Euler angles, [38] and so we prefer to keep this Φ in the definition of the rotation operator. Further, keeping all three Euler angles allows us to directly relate to operations on a quantum computer, where it is usual to include an arbitrary rotation in terms of Euler angles -or at least some operation with three degrees of freedom. This is also the reason we prefer to use U 1 2 (φ, θ, Φ) over R(ξ), despite the latter having the satisfying connection to D(α).
The reason for the requirement of all three angles can be thought of by considering the Weyl function as giving the expectation value of the rotation operator used. When considering the rotation operators R(ξ), the rotation is on CP 1 rather than the full SU(2). This means we have rotations restricted to a 2-sphere. In turn, the rotation operator cannot produce the state σ z , which is necessary for discerning its two eigenstates |↑ and |↓ . This is the reason for the result in Equation (52).
It is interesting to note that when our goal is to use the Weyl function to consider the expectation value of different operators, we can instead consider the Weyl function on a torus, where we set Φ = 0. It turns out that this restriction is informationally complete and the last degree of freedom is not strictly necessary in the qubit case, as long as the range is extended to −π <φ,θ ≤ π. However, we will generally choose keep the Φ degree of freedom to allow the rotation operator to also define the coherent state representation on CP 1 , and more importantly to keep this formulation informationally complete for larger Hilbert spaces.
We now need to use our choice of rotation operator to rotate the corresponding generalized parity operator. For the Wigner function, this is explicitly where σ z is the standard Pauli z operator. The appearance of the √ 3 in Equation (54) may seem somewhat strange, as is generally assumed that the generalized parity should simply be σ z , without the identity or the coefficients. Choosing σ z as a parity operator can be found in some works, and has proven useful in verification protocols, such as in Reference [52], this choice, however, leads to a function that is not informationally complete, i. e. Equation (39) does not hold. Note also that S-W. 3 implies that Tr[Π] = 1, which does not hold for σ z . To generate a generalized parity that is informationally complete, it is therefore necessary to follow the Stratonovich-Weyl correspondence S-W. 1 -5. For more detail of how this can be done, see Reference [5] or in the Appendix Reference [53]. This yields the Stratonovich qubit Wigner function kernel where the Φ degree of freedom gets cancelled out due to the generalized parity being a diagonal operator, this is the reason behind the equivalence of the rotation operators in Equation (50) when rotating the generalized parity operator. We can then simply return the operator by applying this case to Equation (43) If one, on the other hand, insisted on using just the σ z operator, creating a phase-space function W Z ρ (Ω) = Tr ρ U 1 2 (Ω)σ z U 1 2 (Ω) † , to transform back to the density operator, the identity needs to be reinserted however, we therefore have to assume that the function is of a trace-1 operator, and this doesn't hold for an arbitrary operator. This result is clearly shown by considering either Reference [5] or Reference [53] and setting the coefficient of the identity operator to 0. Although this can be useful for examining certain properties of a quantum state, it is important to bear in mind that it doesn't satisfy the Stratonovich-Weyl correspondence. Further, using the σ z operator as a parity creates a function similar to the Weyl function -that is, it is the expectation value of some arbitrary operator in phase space. The only problem is the identity operator is not represented in this distribution like the σ z operator in not represented in the case of the CP 1 rotation operator.
Examples of the single-qubit Wigner function for the eigenstates of the Pauli operators can be seen in Figure 1 (d), where |↑ and |↓ are the eigenstates of σ z with eigenvalues ±1 respectively. Likewise, |→ and |← are the eigenstates of σ x with eigenvalues ±1 and |+ and |− are the eigenstates of σ y with eigenvalues ±1. To the right of each function that is generated by the Stratonovich kernel is the same state generated by the Wootters kernel, that will be shown in the following section. More comparison between the two is therefore to follow. It is immediately clear that when generated with the Stratonovich kernel on a spherical geometry, the Wigner function is similar to the more familiar Bloch sphere, where the highest value of the probability distribution is where the Bloch vector would point. It is also worth noting that, unlike the Heisenberg-Weyl functions and further the discrete Wigner functions, every spin coherent state has a negative region. This is a result of the spherical geometry and the role of quantum correlations within phase space, this will be discussed in more detail in Section 4.2.
Alternatively, one may be interested in constructing the Q or P function for the qubit, where the generation of the suitable kernel is just as simple as in the qubit case. Like the Wigner function generalized parity operator, these too are diagonal operators. The two can be calculated explicitly as for the Q and P function respectively, where the ±1 in the bracket signifies the value of s for the generalized phase-space function. These both result in kernels where the Φ degree of freedom also cancels out when taking a similarity transform with the rotation operator. The cancellation of the Φ term can therefore mean that the last Euler angle is not necessary for the calculation of quasi-probability distribution functions, and one will occasionally see the kernel generated ignoring this term, for good reason to exclude redundant information. This also lines up nicely with the toroidal definition of the Weyl function, where we only care about the first two degrees of freedom in the rotation operator. As we will soon see, a toroidal geometry is a natural alternative to the Bloch sphere in considering qubit states. where the we have the real and imaginary values of the continuous Weyl function respectively. This is followed by the discrete Weyl function.

The Wootters kernel
The alternative formulation to consider qubits in phase space can be seen in the framework set out by Wootters. [15] In the same year, Feynman also gave an identical construction of a Wigner function for qubits. [14] Feynman's construction treated only the qubit case, whereas in the work by Wootters the formalism was extended to any prime-powered dimensional phase space. More importantly, in Reference [15] Wootters presented the discrete Wigner function in the Moyal formalism -the expectation value of a kernel, which Wootters intuitively names the 'phase-point operator'. Despite this, we will give the formulation of the discrete Wigner function presented by Feynman, as we feel it provides a different insight into the problem, and then subsequently relate it to Wootters's phase-point operator. Following this, we will show how this formulation can be considered a sub-quasiprobability distribution of the Stratonovich formulation.
The idea behind Feynman's formulation was to write down probability distributions which are the joint probabilities, p ±± , of finding the qubit having the spin aligned with the ±z and ±x axes respectively. The choices of operator and expectation values in Equation (59) were presented by Feynman without any derivation or motivation; although, as we will see later Wootters uses the same basis, this choice is not unique, in fact there is no unique choice of a probability distribution of this form, see References [54,55,56] for examples of generalizations of these results. 1 Following the results of Wootters, we can take the expressions in Equation (59) and turn them into a phase-point operator. We first introduce two degrees of freedom to span the discrete phase space, these are z, x ∈ {0, 1} which cover all four phase points. The phase-point operator is then defined resulting in the discrete Wigner function for any arbitrary operator where the subscript is dimension d = 2 for A d (z, x). We can then relate the results in Equation (59) by noting that W(0, 0) = p ++ , W(0, 1) = p +− , W(1, 0) = p −+ , and W(1, 1) = p −− . Like the Stratonovich kernel, we can define a generalized displacement operator to both generate a Weyl function and further to be used to construct the Wootters kernel in terms of a generalized displaced parity operator. The displacement operator for discrete systems such as this is which, like the SU(2) Wigner function, can be used to generate a discrete Weyl characteristic function The four values of the discrete displacement operator are simply the Pauli and identity operators. As such, there is a direct relationship between the Euler-angle parameterization and toroidal lattice approach. Note first that by setting Φ = 0 in the SU(2) case, the rotation operator defines translations around a torus. This then allows us to describe the values of the toroidal lattice as points in the continuous torus generates by the Euler angles. We also need to note that there can be a difference in phase between elements in U 1 2 (φ, θ, 0) and D 2 (z,x), therefore, a factor of i may be necessary to equate the two. However, the complex behaviour of characteristic functions is well-known, and it is instead standard practice to define the characteristic as the absolute value or absolute value squared, where such practice is commonplace in the signal processing community, for example when analyzing the analog of the Weyl function -the ambiguity function. [57,58] This allows us to relate the two approaches to a Weyl function by taking the absolute value squared of each or by adding an extra phase, where it's also interesting to see the Weyl function as a generalized autocorrelation function, where for a pure state We can think of the rotation operator as acting on the ket state and then taking the inner product. It then makes sense to take the absolute value squared to get a real-valued autocorrelation function.
Note that in Equation (64), we relate the continuous Weyl function to the discrete case by setting Φ = 0, this is the simplest way to reduce to two degrees of freedom for easy calculation and equivalent results. If, on the other hand, we chose to set Φ = −φ, like we do when equating with the Heisenberg-Weyl case, there would be no way to produce the σ z operator, hence the lack of informational completeness in this slice for the eigenstates of σ z in Equation (52).
The informational completeness for Equation (63) can be seen by observing which demonstrates the power of the use of the Weyl function as a tool for practical purposes in state reconstruction. For example, the discrete form of the Weyl function has proven to be useful in verification protocols, [59] where the Weyl function can be used as a probability distribution for choice of measurement bases. In the case of stabilizer states, this results in the Weyl function having values of 0 over the majority of the distribution and 1 when the Weyl kernel is the stabilizer for that state. This is simply seen by noting that a stabilizer state is an eigenstate of a product of Pauli operators. For a single qubit, this result is simply seen by noting the same naturally holds when extending to multiple qubits. Note, however, that in Reference [59] the Weyl function is normalized by a factor of 1/ √ 2 n , for n qubits. This is done so the factor of 1/2 is not needed in Equation (65) and z,x X (z,x) 2 = 1.
Having defined an appropriate displacement operator that generates an informationally complete Weyl characteristic function, to complete the discrete formulation in terms of a generalized displaced parity operator, it is necessary to generate the generalized parity operator. The generalized parity is, from Equation (40), simply and is the kernel written in displaced parity form. Interestingly, the eigenvalues for A 2 are the same as those for Π 1 2 , namely (1 ± √ 3)/2. This means that the non-diagonal Wootters generalized parity can be diagonalized as where, In fact, like the Weyl function, the full Wootters kernel can be written in terms of the Stratonovich kernel, where therefore, the Wootters kernel can be considered a subset of the Stratonovich kernel, and the Wootters qubit Wigner function is a sub-quasidistribution function of the Stratonovich qubit Wigner function. Note that the four points on the Stratonovich qubit Wigner function where the Wootters function lies form a tetrahedron within the Bloch sphere, as can be seen in Figure 1, coinciding with the formulation of symmetric informationally complete projective measurements. [60,61] We can transform between the Wigner and Weyl functions here by performing the discrete Fourier transform and The Fourier transform between the Wigner and Weyl function can therefore be seen as the transformation between Pauli basis measurements and symmetric informationally complete projective measurements. From Equation (69) and Equation (70), where we consider the Wootters kernel as a subset of the full Stratonovich kernel, we can define Q and P functions in discrete phase space. By using the parity operators from Equation (58) we can define the general kernel that yields For the Q function, this results in a set of symmetric informationally complete POVMs. Considering this in terms of frame theory in quantum measurement, [62,22] consider this discrete Q function as a frame, where the dual frame is then the P function kernel, such that More generally we can transform between any of these kernels by By understanding these phase-space functions in these terms, they can prove to be powerful tools for verification of quantum states, for example Reference [63] which we will discuss further in Section 4. The relation between the Stratonovich and Wootters formulation of phase-space functions for qubits is a special case that doesn't generalize as the dimension increases. This is because of the difference in geometry between the representations. Following this discussion, we will see two variations in how a qudit can be represented on a continuous phase space. One representation keeps the SU(2) group structure and is the natural extension of the Stratonovich kernel, that was also introduced in his seminal paper. [5] This method creates a discrete Wigner function on a sphere, where the radius is proportional to the dimension of the Hilbert space. The second continuous method uses the structure of SU(d), where d is the dimension of the Hilbert space. Geometrically, this is an oblate spheroid.
The Wootters function, on the other hand, is defined on Z(d) × Z(d) discrete toroidal lattice, where d is prime. There have been many formulations to extend this form of the Wigner function to higher dimensions, the original Wootters formulation was defined with respect to fields and so the dimensions size was restricted to primes. Wootter then extended this formulation with Gibbons et al to work for dimensions that are powers of primes. [18] Later work also considered the extension to odd numbered dimensions, [16,19,21] , and even numbered dimensions. [17] Given the different choices, we refer the reader to any of these articles for specifics, or to Reference [21,22,64] for more detailed overview and comparison. Since this construction above the qubit case isn't important to this review, the explicit forms will not be given.

SU(2), spin-j
We will now show the generalization of continuous phase space for larger finite systems. Here we begin by staying within the SU(2) Lie group structure before going on to the SU(N ) formulation in Section 3.3. Both cases can be used to consider similar situations, but also have some subtle differences. The SU(2) formulation requires another variable M = 2j, two times the azimuthal quantum number, that allows us to represent a spin-j qudit where the dimension of the Hilbert space is d = M + 1.
It can alternatively be used to model an M -qubit state in the symmetric subspace, and can be used over composite methods for certain processes, such as single-axis twisting or Tavis-Cummings interactions. [65,66,67,68] These methods can further be used on states that aren't fully symmetric, by taking the approach from Zeier, Glaser, and colleagues [69,70] a multi-qubit state can be represented by a collection of SU(2) Wigner functions, n ≥ M ≥ 0 for n qubits, to show the Wigner function for every M -valued symmetric-subspace. This leads to many different Wigner functions as the number of qubits increases.
We begin by considering the general SU(2) rotation operator. Like with the qubit case, there are many ways of formulating rotations on a sphere. The extension to the two ways given in Equation (50) are simply where are the the angular momentum operators, the higher-dimensional extension of the Pauli operators in SU (2), and J ± j = J M 2 (1) ± iJ M 2 (2) are the angular momentum raising and lowering operators. Here we also introduce the generalized SU(N ) operator notation J M N (i). The intuition from looking at R j (ξ) is that as j → ∞, the Heisenberg-Weyl displacement operator from Equation (5) is recovered. This result was shown by Arecchi et al. in Reference [29]. The procedure to show this was also demonstrated in Reference [12], where first we need a contraction of the operators, where we take such that, when c → 0 The rotation operator therefore becomes by setting we yield lim A similar treatment can be applied to the Euler angles by noting that J M . This leads to the result that, when taking Φ = −φ the limit on the right goes to 1, yielding an equivalence to D(α). Further, when applying this to the vacuum state, we yield equivalent to the coherent state generated by the standard displacement operator. Where the φ and Φ degrees of freedom act as a global phase and cancel out when we set Φ = −φ. Note also in this result that the limit of the highest weighted spin state is the vacuum state. By analogy, the definition of a generalized spin coherent state is just the highest weight state rotated by a suitable operator. Either rotation operator will work, however they may differ by a phase factor. For simplicity, we simply define the spin coherent state by This results in the simple definition of the Q and P functions with respect to the spin coherent states where it has been shown in Reference [51] that these two functions also reach the Heisenberg-Weyl Q and P functions for j → ∞.
The general idea for this infinite limit of these functions is to imagine a plane tangent to the north pole of the qubit's sphere. As the value of j increases, the radius of the sphere correspondingly increases, coming closer and closer to the tangential plane at the north pole. Where at the limit, the sphere completely touches the plane.
This visualization exercise also helps to understand the difference in Weyl function between SU(2) and the Heisenberg-Weyl group, where the general SU(2) Weyl function is As mentioned in the single qubit case, the general SU(2) Weyl function requires all three Euler angles. This can be seen in a variety of ways, most importantly it is because when using R j (θ, φ) as the kernel, the resulting function for |j, j and |j, −j are identical. Furthermore, any mixed state on the line inside the Bloch sphere between these two states are also identical. To resolve this, a third degree of freedom is needed, resulting in the preference for the Euler angles. This is not a problem for the Heisenberg-Weyl case as any state on the equator is mapped to infinity, and only the equivalent of the top hemisphere is considered in the phase space. It is also worth noting that in the qubit case, the Weyl function is still informationally complete when Φ = 0, this is not the case with a general qudit. As we saw when considering the Wootters kernel, there's a special connection between the qubit on the sphere and the qubit on a torus, this means that for the Weyl function, setting Φ = 0, we are modelling the qubit state on a continuous torus. However, due to how the su(2) algebra scales as as M increases, the state can no longer be considered on a torus and the Φ degree of freedom is necessary to generate an informationally complete function. Next we need to show how the Wigner function can be generated. Since the original formulation of the general SU(2) Wigner function, there have been many variations in how it is presented. From a harmonic series approach to presenting it as a kernel with which a group action can be taken. Since all of these representations are equivalent, we will present how the kernel can be generated through a multipole expansion.
For this, we need a new set of generators to generalize the Pauli operators from the qubit case. These are known as the tensor or Fano multipole operators, [71] T j lm = 2l + 1 2j + 1 j m ,n=−j C jn jm ,lm |j, n j, m | , where C jn jm ,lm are Clebsch-Gordan coefficients that couple two representations of spin j and l to a total spin j. The kernel can then be generated from there operators by where Y * lm (θ, φ) are the conjugated spherical harmonics. By following the same procedure from Equation (40), the SU(2) generalized parity can be derived where we note again that j = M/2. Note that when taking the limit j → ∞, this generalized parity operator becomes the parity operator defined in Equation (4). The derivation and proof of this can be found in Reference [12]. Resulting in the generalized displaced parity operator where, since the Φ degrees of freedom cancel out with the diagonal generalized parity, the full kernel is equivalent to the Heisenberg-Weyl kernel from Equation (3) as j → ∞. [12] More generally, the Q and P functions from Equation (91), and any s-parameterized quasi-probability distribution function, can be generated by a simple extension of Equation (94) and Equation (95), where the general kernel can be calculated [33,50] where the extra value (C jj jj,l0 ) −s is 1 when s = 0, returning Equation ( which for s = −1 yields Π Where, as mentioned earlier, these functions also become their infinite dimensional counterparts as j → ∞. [51] From Equation (44), we can then relate and transform between these different quasi-distribution functions by way of a generalized Fourier transform. However, there is an alternative method by way of spherical convolution [50] . By using the highest-weight state in the Wigner representation and performing a convolution using Equation (47) we can transform an s-valued function to the corresponding (s − 1)-valued function for the same state. For example, we can transform from the P function to the Wigner function, and the Wigner function to the Q function by using this method. Similarly, we can transform from the P function of a state to the corresponding Q function by using the convolution with the Q function highest weighted state.

SU(N)
We now consider an alternative route to treat a general d-level or multi-qubit quantum system. Like the SU(2) case, an SU(N ) representation can also be used for a general qudit state, where in this case d = N and we only consider M = 1. In order to represent an n-qubit state, the construction of an SU 2 d kernel is required -it is also possible to take the same approach as with SU(2) and consider only the symmetric subspace with an SU(d) representation. Alternatively, the generalized displacements and generalized parities can be mixed and matched to some extent as was shown in Reference [37], this case will also be considered here.
Although we provide a formalism to construct these Wigner functions, the details to actually construct the kernels can prove to be complicated, there are therefore some open questions for generalized SU(N ) kernels. the first open problem is generalizing SU(N ) for any value of M , effectively creating the analog of the symmetric subspace for M N -level qudits. However, the procedure to create the generalized displacement operators for this case can be found in [36].
The second open problem is the consideration of the hyperbolic space of SU(Q, P ). A construction of an SU(1, 1) Wigner function has recently been provided in References [72,73], however the general SU(Q, P ) case remains to be solved. Insight into these cases may lead to some interesting results in ideas relating to AdS/CFT in a phase space formalism. Note that if one is also interested in some of the more fantastical routes to the unification of quantum mechanics and gravity, in theory a construction of E 8 should be possible, however unwieldy the number of degrees of freedom prove to be. Here we will settle ourselves with providing the results from the SU(1, 1) phase space construction to close this section.
The kernel for an SU(N ) phase space now requires a specific generalized displacement and generalized parity. The first step is to construct the generators of SU(N ), which can be represented by N 2 − 1 matrices of size N × N . The method to construct these can be found in Reference [74]. This results in a set of In fact, only a subset of the generators are needed to construct the SU(N ) Euler angles, these are generalizations of the σ y and σ z operators, which are where 2 ≤ a ≤ N . Note that this is a subset of the full set of generator of su(N ); N 2 − 1 elements are needed to generate su(N ) however we just 2(N − 1) elements of the algebra are needed to generate the generalized Euler angles for SU(N ). We also note that there has been an alternative construction of phase space for SU (3) where they instead used three dipole and five quadrupole operators as the generators of the Lie algebra [75] . Following this, there is a procedure to generate the generalization of the Euler angles for SU(N ), that can be found in References [46,47,74]. This results in the generalized Euler angles where Finally yielding the SU(N ) rotation operator Just from this generalization of the Euler angles, many of the results from the SU(2) case can be extended to SU(N ). Most straight-forward is an informationally complete Weyl characteristic function that takes the form Note that a discrete alternative can be yielded by constructing all N 2 − 1 generators for SU(N ), resulting in where 0 ≤ k ≤ N 2 − 1, J 1 N (0) = 1 N and the other J 1 N (k) matrices are the generators of the SU(N ) algebra. Also from the rotation operator in Equation (105) one can generate generalized coherent states. [74] Like earlier examples of coherent states, SU(N ) coherent states are generated by applying the rotation operator to a suitable analog to the vacuum state. In a similar procedure SU(N ) coherent states can be generated by the |0 1 N state, where we define this state as being the eigenstate of λ N (N 2 − 1) with the lowest eigenvalue, −[2(N − 1)/N ] 1 2 , making it the lowest-weighted state, resulting in ignoring an overall global phase. We note that it is also possible to set the analog to the vacuum state as the eigenstates of J 1 N (3) with eigenvalue 1, as is the more standard choice for the use of SU(2) for qubits. However, given the structure of the generators for SU(N ), taking the highest-weight state results in much tidier calculations.
From the construction of a coherent state, the generations of the Q and P functions are simply where Ω are the degrees of freedom on SU(N ) and dΩ in the volume normalized differential element for SU(N ), the construction of which can be found in Reference [46]. Similarly, the kernel for the SU(N ) Wigner function can be generated with respect to coherent states, where [36] By setting the variable to 0, the generalized parity operator for SU(N ) is yielded, where producing an operator constructed by the identity operator and the diagonal operator from Equation (101), where a = N . Like the construction of an SU(2) generalized parity operator, the SU(N ) parity is a diagonal operator and commutes with theŜ 1 N (Φ) part of the rotation operator, resulting in the degrees of freedom cancelling out. The number of degrees of freedom are further reduced to allow the SU(N ) Wigner function to live on the complex projective space CP N −1 with 2(N − 1) degrees of freedom. By considering the operators, this can be realized by noting that each SU(N ) generator acts on N SU(2) subspaces of the full SU(N ), where each subspace is acted on by J N y (1, a). Since N − 1 elements of the generalized parity operator are the same, with just the last element differing, it in is effect a weighted identity operator for the rotations that do not act on the last SU(2) subspace. And so, the elements of the rotation operator that have as the generator J N y (1, a) for a < N cancel out, until we reach the element that rotates around J N y (1, N ). Note that if we stay with the fundamental representation of SU(N ) there is a simple transformation between the different quasiprobability distribution functions. This can be shown by noting that the generalization of Equation (110) is [36] Π 1(s) resulting in the general function We can then simply transform between functions by taking providing a general formula to consider qudits in SU(N ) similar to Equation (78) for qubits. Note in the same way as the qubit case, one can define a set of symmetric informationally complete measurements from the Q function, this can then be transformed into its dual by which corresponds to the transform between an N -dimensional frame and its dual frame. [22,62] This is all provides a structure for considering qudits in phase space, alternatively by considering the SU(2) subgroup structure of SU(N ), we can provide a link between the different representations of multiqubit states. It also means that we can alternatively represent multi-qubit states with a hybrid approach of SU(2) and SU(N ) that was introduced in Reference [37] and applied in Reference [42]. From applying the procedure from Equation (48) to generate a kernel for composite systems, for multi-qubit states we get where U n (φ, θ, Φ) is the n-qubit rotation operator. We can take this new n-qubit rotation operator and apply it to the SU(N ) generalized parity operator to create an alternative representation of a multi-qubit state in phase space that also obeys the Stratonovich-Weyl correspondence S-W. 1 -5. Each choice of representing qubit states has its own benefits and drawbacks. We can immediately see that the benefit of using the multi-qubit SU(2) rotation operator with the SU(N ) generalized parity operator is the reduced number of degrees of freedom. The multi-qubit rotation operator results in 2n degrees of freedom, whereas the SU(N ) operator will result in 2(2 n − 1) degrees of freedom for n qubits. While on the topic of the SU(N ) representation of quantum states, it's also worth noting a slightly different formulation, that is SU(P, Q) that has the same number of degrees of freedom as SU(P + Q). However, the geometry of SU(P, Q) differs significantly from the compact SU(N ) geometry, resulting from a difference in sign between the first P elements and the remaining Q elements in the metric, yielding a non-compact, hyperbolic geometry. The case of a generalized SU(P, Q) representation in phase space is still an open question. However the specific case of SU(1, 1) has been addressed in a recent paper. [72] The SU(1, 1) representation of certain quantum systems has proven to be useful in simplifying some particular problems in quantum mechanics. In particular Hamiltonians involving squeezing. This is because the SU(1, 1) displacement operator is a squeezing operator, resulting in the SU(1, 1) coherent states to actually be squeezed states. Since the squeezing of quantum states leads to many important results in quantum technologies, such as metrology in particular, the use of these phase-space methods for squeezing can prove to be a useful tool.
The generators of algebra of SU(1, 1) are similar to SU(2) with a difference in the commutation relations, where The raising and lowering operators are then constructed that, along with K M z , constitutes an alternative basis, where the generalized displacement operator is much like Equation (79), where ζ is the complex squeezing parameter. The central result from Reference [72] the generalized displaced parity operator for SU(1, 1) is constructed from Equation (120) and the SU(1, 1) generalized parity (−1) K M z , yielding as the kernel for SU(1, 1). Some of the authors from Reference [72] also looked at the equivalence between SU(2) and SO(3) in phase space . Such analysis could also be used to go between SU(1, 1) and SO(1, 2) in phase space, or even generalized to SO(1, 3) or higher, to consider relativistic quantum mechanics in phase space. There have been works that have looked at generalizing the Wigner function for relativistic particles, for instance see References [76,77,78,79]. Although the results are interesting and potentially useful, it is beyond the scope of this paper to discuss them here.

Quantum technologies in phase space
Given the mathematical framework, we will now look into how the representation of quantum systems in phase space can be used to understand processes in quantum technologies. We will give some examples of experiments that have taken place to directly measure the phase-space distribution of quantum states. This will include both discrete-variable and continuous-variable systems, as well as the hybridization of the two. We will also mention ideas that haven't yet been put into practice, but are worth discussion.
Before the practical results, we will consider how important metrics already present in the quantum information community can be described within the phase-space formulation. One metric that is unique to phase-space functions is the presence of 'negative probabilities'. Due to the uniqueness and unintuitive nature of such a phenomena, this will receive a full discussion shortly, where we will discuss what negative values do and do not say about a quantum state. Before looking at negative volume of the Wigner function as a metric, we will consider other figures of merit that can be calculated through phase-space methods

Measures of quality
There are measures that are widely used when considering the density matrix or state vector of a quantum state that can be expressed within the phase-space representation of quantum mechanics. In particular, considering the Wigner function and a key property of traciality from S-W. 4, where many similarities naturally fall out. The most natural starting point is to consider the expectation value of some observable A with respect to some state ρ such that In the case of a qudit in SU(2), we note that the expectation value of the angular momentum operators can be related to the centre of mass for the corresponding Wigner function, where where [80] f where the minus signs and normalization come from our formulation of the Wigner function kernel, which may differ in sign is following a different convention.
There are further important properties in quantum information that utilize the trace of two operators, as such their construction in the phase-space representation are straightforward. A simple example of this is the purity of a state, which for density matrices is P(ρ) = Tr ρ 2 , can be calculated in phase space by Next we consider the fidelity of a state ρ 2 with reference to a target state ρ 1 , where We are often interested in comparing an experimentally generated state ρ 2 with a pure target state ρ 1 , in which case the fidelity reduces to F(ρ 1 , ρ 2 ) = Tr [ρ 1 ρ 2 ], which when applied to Equation (122) we simply yield The phase-space calculation of fidelity leads to an interesting method of fidelity estimation of measured quantum states. We will provide the procedure how one can do that here, following the work on Flammia and Liu. [59] First, it is important to note that the fidelity can be calculated using any valid distribution F (Ω), including the Weyl function and the Q and P functions, so long as we have the corresponding dual functioñ F (Ω) resulting in where in this sense the Wigner function is self-dual, and the P function is the dual to the Q function, and vice versa. One drawback from the structure of Equation (129) is that we need to integrate over the full phase space, however we know from Section 3 that only four points in the phase space are required to produce an informationally complete function. Therefore for n qubits we require a sum over 2 n points in phase space. We can then calculate fidelity through the sum Note that for this to hold, and for later calculations, we also require the functions to be normalized such that k (F ρ1 (k)) 2 = k (F ρ1 (k)) 2 = 1. Following Reference [59], we then need to create an estimator for the fidelity. For this, select any measurement k with probability Pr(k) = (F ρ1 (k)) 2 . The estimator is then where E[X] is the expected value of X. When measuring in practice, calculating the function F ρ2 (k) perfectly is not possible, and so we need to make many copies of ρ 2 and perform sufficient measurements to assure a given confidence of fidelity. The details of which can be found inReference [59].
Here we would like to note that besides being applicable to any number of qubits, such a procedure can also be used for any number of qudits, where we can use the su(N ) algebra to generate a qudit Weyl function or a set of informationally complete measurements on the manifold of pure states in SU(N ).
Fidelity is a measure of how close a measured state is to some target state. We may conversely be interested in a distance measure, that measures how far away (in Hilbert space) a given state is from some target state. Such a distance measure is known as the trace distance, where For a qubit, Equation (132) can instead be expressed in terms of the Wigner function, where This can be seen by noting the two states can be expressed as resulting in the trace distance being Similarly the Wigner function can be expressed from Equation (110) along with the properties in Equation (125). This results in Such a procedure can then be generalized to multi-qubit states. Further, similarly to the calculation of fidelity, we can exchange the integral over the continuous degrees of freedom by performing a sum over discrete degrees of freedom. Where the procedure can be applied to not only the Wigner function, but the all distributions, albeit with a different normalization out the front. We also note that alternatively, one can find the Monge distance with respect to the Q function. [81] This effectively calculated the distance between the centers of mass of two distributions. Given the non-negative property of the Q function, it can be useful in translating procedures performed on classical probability distributions into quantum mechanics. As such, it can also be used in calculating the Rényi entropy. The Rényi entropy is defined with respect to the Q function as [82,83] This result was extended to Wigner functions in Reference [84], however this is limited to the Heisenberg-Weyl case. We note here that in the case of SU(N ), one can simply use Equation (114) to write Equation (139) in terms of the Wigner function.
The above set of equations provide a framework that show how the phase-space representation of quantum mechanics can prove a viable alternative to the more traditional density operator approach. All of the above can also be extended for any choice of phase-space function, such as the Q, P , and Weyl functions. And where the equations were specifically states for qudit states, general expressions can be found to calculate these values for SU(N ) structures of even continuous-variable systems. Further to these measures, negativity in the Wigner function plays an important role, where the manifestation of these negative values provide nuanced details for each system. A closer look into these measures will be considered next.

Negative probabilities
Negative probabilities have been a subject of interest since the early days of quantum mechanics. From the observation of negative values in the Wigner function, Dirac later discussed their presence along with negative energies in 1942 [85] , where he stated they 'should not be considered as nonsense. They are welldefined concepts mathematically, like a negative of money.' This growing acceptance to the concept of negative probabilities then lead to a number of other people to take the concept seriously. A more rigorous exploration was undertaken two years later by Bartlett [86] . Following this, in his PhD thesis, Groenewold argued that to satisfy other desirable properties of a probability distribution function in phase space for quantum mechanics, it is necessary for negative probabilities to exist. In the '80s, Feynman then gave a simple argument for the existence of negative probabilities in constructing a discrete Wigner function. [14] . Put simply, his argument was that negative probabilities are no more nonsense than a negative quantity of anything (in his example he gave apples), and that the negativity is just one step in a total sum. The negativity is never considered in isolation.
Any probability over phase space needs to be averaged over a quadrature, due to the inherent uncertainty of quantum states. When a measurement is made over a quadrature, the resulting probability is always nonnegative. The Wigner function is unique among different phase-space distributions as it results in the correct marginals, see Equation (8). One might be more inclined to choose the Q-function as it is always non-negative everywhere, however it doesn't satisfy the same property as Equation (8).
But, what does the negativity actually tell us when it appears in the Wigner function? It is generally accepted that it is a result of quantum correlations, or quantum interference. By putting two coherent states in a superposition, the coherences from the two states will interfere creating oscillating positive a negative interference in between, that we interpret as quantum correlations. An important thing to note here is that the volume of the positive and negative regions in these correlations are equal, and so cancel out when integrated over. Meaning that both the superposition and statistical mixture of two coherent states result in the same overall integral of 1.
The full story is that this negativity is actually a consequence of non-Gaussianity, a result that is known as Hudson's theorem. [87] Consequently negativity doesn't capture the important quantum correlations that arise from squeezing, in particular two-mode squeezing. Which, since the state is a squeezed Gaussian, has no negative quasi-probabilities. So some care needs to be taken in using negative volume as a measurement of quantumness, non-classicality, entanglement and so on. That being said, it doesn't mean negativity shouldn't be used, and there have been many results showing its use. [88,89,90,91] It just may be more important to look at the context in some cases, and look for how correlations manifest in particular ways.
In the same way, the interpretation of negative values in discrete quantum systems needs some care. Any coherent, single-qubit state generated by the Stratonovich kernel will exhibit negative values. The maximum positive value is in the direction that the Bloch vector is pointing, as we deviate from that point over phase space, the probability slowly decreases to zero and then negative, where the maximal negative value is at the point orthogonal to the Bloch vector. Because of the term 'classical states' for coherent states, this tends to lead to the confusion that coherent qubits states, such as the standard choice of basis in quantum information |↑ and |↓ , are classical. This term 'classical' is just a misnomer, a qubit state is fundamentally non-classical, likewise with any quantum coherent state. By definition, any coherent qubit state is quantum; the only non-quantum state that a qubit can exist in is the fully mixed state. Because of this quantumness and that the phase space is represented on the sphere, the center of the Gaussian on the sphere creates interference with itself on the opposite side of the sphere, resulting in the negative values. Negativity can therefore be used as a sign of quantumness, but it's not necessarily a sign of entanglement.
One could also make the argument that this negativity can be considered a sign of superposition, as |0 is a superposition of (|0 +|1 )/ √ 2 and (|0 −|1 )/ √ 2, and in fact any pure state for a qubit is the superposition of two antipodal pure states. In this way, the negativity is a result of this and the overall symmetry of qubit states; this really shows the beauty of the Stratonovich kernel, how it presents the symmetry in coherent qudit states, where every N -level coherent state is just the SU(N ) rotation of the chosen ground state.
If one decides that this symmetry is not so important and require that specific spin coherent states are non-negative, then they can consider the Wootters kernel instead. The main benefit to the Wootters kernel is that the eigenstates of the Pauli matrices all produce non-negative Wigner functions. In fact, if an octahedron is drawn within the Bloch sphere, where the vertices are the eigenstates of the Pauli matrices, the discrete Wigner function for every state that lies on and within this octahedron is non-negative. These are the so-called magic states, [92] if a pure state rotates between any of these magic states, then negativity arises -a result that is made clear by considering the Wootters kernel as a subset of the Stratonovich kernel. This gives rise to the literature on the equivalence between Wigner function negativity and contextuality, and further work on magic state distillation, proving to be an advantage to choose the Wootters kernel over the Stratonovich kernel. [93,27,25,26] The negativity that arises from using Wootters kernel can further be summarized by considering the extension of Hudson's theorem for finite-dimension systems produced by Gross. [21] Alternatively, when using the Stratonovich kernel, one could consider an analog of Hudson's theorem to apply to the decreasing volume of negativity in the highest weighted state, |j, j , as j → ∞. At the infinite limit the state becomes Gaussian and the Heisenberg-Weyl group is reached -the highest weighted state is simply the non-negative vacuum state. In the same way, by choosing the Dicke state with one energy level lower |j, j − 1 , as j → ∞ this slowly transforms into the one-photon Fock state. And so on for |j, j − n producing the discrete analog of the n-photon Fock state. This relation can even be seen in the qubit Wigner function using the Stratonovich kernel, where the spin coherent state appears to look like the vacuum state from one side, and then the one-photon Fock state from the orthogonal direction. This interesting relation between the representations for the two different systems further leads to ways to identify signatures of quantum correlations in systems that have heterogenous degrees of freedom, such as spin-orbit coupling [94] or light-matter interaction. [95] In fact, there is even correlation with the negative volume of the hybrid Wigner function and quantumness. [96] Where the negative volume of a Wigner function with continuous degrees of freedom is which can be seen as an additional measure of quality derived from the Wigner function.

Methods for optical systems
Phase-space methods have gained much attention in the optics community, where there is a large collection of books on quantum optics that include in-depth treatment of phase-space methods and the application thereof. This will therefore not be a thorough account of the advancements made in quantum optics using these methods, we will instead just point out some results we find particularly interesting and will prove useful later in this section. For anyone who wants to read more deeply into the application of phase-space methods in quantum optics we suggest the following heavily truncated list of book on quantum optics. [97,98,99,100,101] There have been further, in-depth reviews of the use of the Wigner function in optical systems [102,103] and in the use of optical quantum computing [104,105] . When generating the Wigner function experimentally, there are many approaches one can take. The two we are interested in here are first the generation of the density operator and then following Equation (2) to produce the Wigner function. The second is to take direct measurements of points in phase space, by passing the need to calculate the full density operator. The real power of measuring quantum systems in phase space comes in realizing the displaced parity formalism of the kernel. By realizing the invariance of cyclic permutations of elements within a trace operation, it is clear that the Wigner function generation can be seen as the parity measurement of a displaced state: where D(−α) = D † (α), and ρ(α) = D(α)ρD(α) † is the rotated density operator.
This insight amounts to an experimental procedure of creating the quantum state, then displacing the state by some value α and then measuring this state. At this point the parity operator can be applied through classical computation of diagonal elements, or by following the procedure in Reference [106], one can perform a parity measurement by coupling the optical state to a two-level atom, where the two levels act as the ±1 parity values.
Direct measurement of the Weyl function, on the other hand, hasn't gained so much attention. Similarly to the measurement of the Wigner function, this can be executed by coupling to a qubit, where instead of taking parity measurements, the qubit is measured in the x and y Pauli basis. Following the procedure laid out in References [107,108], we start by coupling the optical state, ρ f with a qubit in the state |+ = (|↑ +|↓ )/ √ 2, such that We then perform the unitary transformation note that this is the Pauli z operator with the displacement operator with a factor of a half. Once the state has been prepared in the appropriate point in phase space, α, the qubit is then measured in either the x or y basis, resulting in where σ x = Tr [σ x ρ tot (α)] and σ y = Tr [σ y ρ tot (α)]. We note that a similar procedure can be performed by coupling a qubit to a qudit to generate the Weyl characteristic function for qudit states.
Once we have generated the Weyl function, we can use Equation (17) to perform a simple Fourier transformation to generate any quasi-probability distribution function. An important part of the analysis of optical states is the investigation such transform along with Gaussian smoothing. One such case where this approach is insightful is when considering two-mode squeezed states. Such states are Gaussian and so have no negative values, despite being highly entangled. Generating a characteristic function and then filtering to a quasiprobability can result in negative values. [109,110] Allowing a phase-space analysis of such states. [111,112] A similar approach has also been taken for analyzing discrete states. [113] An experimental way to achieve these results should be possible in combination with an extension of Equation (144), we will now consider how else one can directly measure the phase space of a discrete quantum state.

Measuring qubits in phase space
By building on the techniques used to measure phase space for continuous variable systems, it is possible to visualize discrete variable systems in phase space. Some early work sought to first recreate the density operator for the state, before applying the appropriate kernel. One such example of this method was to show the creation of an atomic Schrödinger cat state in an experiment using a Rydberg atom with angular momentum j = 25. [114] This work takes the atom as a large angular momentum space, with a Hilbert space of dimension 51.
Alternatively, one can measure a multi-qubit state by considering the collection of symmetric subspaces that the full Hilbert space is made up of, in a tomography method known as DROPS (or discrete representation of operators for spin systems). [69,70] Other methods would rely on taking the j = n/2 symmetric subspace of an n-qubit system, the power of using the DROPS method is that you get much more information about the state -including information about antisymmetric states such as the two-qubit singlet state, which cannot be represented by the j = 1 symmetric subspace. Another benefit of this method is that the authors have made vast libraries of results available, including code to generate these functions. [?] The methods already mentioned rely on the full construction of the density matrix for the state, requiring full tomography before a phase-space representation can be constructed. As the number of qubits increases, this can get exponentially difficult, we therefore desire methods to represent a quantum state in phase space that don't require full density matrix reconstruction.
Our main focus from now on is the Wigner function and to present results that utilize the displaced parity structure of the kernel, where we can imitate the Heisenberg-Weyl group result from Equation (141). By using the framework outlined in Section 3 these results can be adapted into whichever phase-space distribution is preferred, as they are all informationally complete representations of a quantum state. First we will consider the simplest case of phase-space functions generated by a single SU(2) kernel construction from Equation (96). This will then be followed up with utilizing Equation (116) to show how a composition of single-qubit kernels can be used to measure the phase space of quantum states directly.
Adapting Equation (141) for use in general SU(2) systems results in The procedure is to rotate the generated state ρ to a desired point in phase space, creating the state ρ(φ, θ, Φ), followed by taking the appropriate generalized parity measurement for the group structure.
When experimentally measuring phase space, how Equation (145) is practically implemented depends on the system that you're measuring and how the measurement results are collected. This means that practically, the rotation and generalized parity measurement needs to be put in the right context. Two examples of this when measuring a single qubit can be found in recent works, Reference [115] and Reference [80], where in Reference [115] the authors measured a two-level caesium atom, and in Reference [80] authors measure a state produced in a nitrogen vacancy center in diamond. In both cases, they first produced a state in a superposition of the excited and the ground state.
Using Equation (145) by applying pulses to rotate this state by values of θ and φ to given point in phase space. Projective measurements are then taken at each point in phase space, resulting in the spin projection probabilities p m (θ, φ) for m = ±1/2, where we can think of these projective measurements as where the P ± 1 2 are the projectors of the eigenstates of σ z /2 with eigenvalues ±1/2 respectively. These can also be thought of as the diagonal entries of the density matrix ρ(φ, θ, 0) from Equation (145). The generalized parity is then applied to the spin projection probabilities by calculating where Π M 2 mm is the m th diagonal entry of the generalized parity operator. Note that here we're using the same normalization of the Wigner function used throughout the text, where as the results in Reference [115] and Reference [80] differ by a factor of 2π.
Results from the two experiments in are shown in Figure [115]. In Reference [80], they measured 1200 points in phase space, for θ values with step size π/60 and step size π/10 for values of φ. They also demonstrated the consequence of decoherence on the Wigner function for this state, in Figure 2  Like with the Wigner function for optical systems, phase-space methods prove useful for showing the coherence in qubits. However, for systems of qubits to be useful for the advancement of quantum technology, it is important to consider multi-qubit systems. This is where the power of phase-space methods becomes more apparent. For a single qubit, full tomography isn't that difficult a task, and requires fewer measurements than measuring all of phase space. However, when the number of qubits increases, the size of the Hilbert space increases exponentially, and so measuring points in phase space becomes more of a viable alternative.
As we have already mentioned, there are two main approaches to represent a multi-qubit state. We start with the symmetric subspace, where we are limited to considering symmetric states of composite qubits. In permission. [115] Copyright (2018) by the American Physical Society. many situations, symmetric states are all that are desired from the output for use in computational settings or metrology. Examples of such symmetric states are coherent states [29,30,31] ; GHZ states [116] , or in this case atomic Schrödinger cat states [65] ; W states, or Dicke states [67] ; and spin-squeezed states. [65,117,118] It is therefore important to understand how we can directly measure the phase space of these larger Hilbert spaces, presenting an alternative to traditional tomographical methods. It has already been noted that the Wootters-kernel Wigner function can be seen as being a sub-distribution of the Wigner function generated with the Stratonovich kernel. Since both functions are informationally complete, this means that by taking four carefully chosen points in phase space, the full state can be reconstructed. In fact, these points don't need to be the four that map the Stratonovich kernel to the Wootters kernel, they just need to be appropriately distributed around the phase space.
Similarly, it is possible to reconstruct the full quantum state, and therefore the full phase-space distribution, by taking a subset of phase-point measurements for any value of j. It was shown in Reference [50] that one can similarly reconstruct the full state by measuring a finite number of points. This can be realized by generalizing Equation (147) to an SU(2) Hilbert space of any arbitrary dimension, where note here we've restricted discussion to the Wigner function, but this can easily be generalized to any phasespace function, see Reference [50] for more details. By using the result from Equation (94) that the kernel can be generated from tensor operators and spherical harmonics, points in the phase space can be redefined in terms of a spherical harmonic expansion, with coefficients c lm that correspond to the spherical harmonic Y lm (θ, φ). The full phase-space function can then be reconstructed by calculating c lm from the Wigner function at (4j + 2) 2 points in phase space, W (θ a , φ b ), where θ a = aπ/(4j + 2) and φ b = 2bπ/(4j + 2). Note that the 4j + 2 points are put as a lower bound, and more points may be necessary given the intricacies of the state and phase space being measured. However, it is also noted in Reference [50] that this is just the case by equally distributing the angles over phase space resulting in bunching at the poles, using other methods to distribute points over the sphere, such as taking a Lebedev grid [119,120,121] may reduce the number of points in phase space needed for reconstruction. For example, taking the case where j = 1/2 would result in 16 points in phase space, where we know from Equation (71) the full state can be reconstructed from 4 points in phase space, distributed as a tetrahedron over the surface of a sphere. It stands to reason that this reduction in points will scale as we increase j. A method to reduce these points for specific states was demonstrated in Reference [122]. We add that since these (4j + 2) 2 points in phase space constitute an informationally complete set of measurements, one could then adapt Equation (130) to calculate the fidelity from these points in phase space.
Alternatively, to measure systems of multiple qubits, we can consider a kernel that is the composition of multiple SU(2) kernels from Equation (116). As was shown in Equation (116), measurement of the phase space for such a composite state can be performed by local rotations on each qubit and then taking a generalized parity measurement of the resulting state. Each qubit then has two degrees of freedom from the local rotations, this then results in a function that has 2n degrees of freedom.
In order to make sense out of this many degrees of freedom, it is necessary to choose certain slices of this high-dimensional phase space. The most natural choice is to take the equal-angle slice, this is where the number of degrees of freedom is reduced to two, (θ, φ), where we set every θ equal, such that θ 1 = θ 2 = ... = θ n = θ. Likewise, we set all the φ degrees of freedom to equal each other. The reason this is a natural choice is because the resulting phase-space function shares similarities with the symmetric subspace variant. For example a j = n/2 two-component atomic Schrödinger cat state is similar to the equal angle slice of an n-qubit GHZ state, where in both cases the Wigner function has two coherent states orthogonal to one-another, say at the north and south pole, with n oscillations around the equator between the two coherent states. These oscillations are the signature of quantum correlations arising from superposition in Schrödinger cat state and entanglement in GHZ states. In fact, there are similarities between any symmetric n-qubit state in the equal-angle slice and the corresponding j = n/2 function. It's also worth noting that for an n-qubit Q function for a symmetric state, the equal-angle slice is indistinguishable to the j = n/2 Q function for the comparable state.
The equal-angle slice for the Wigner function was also shown to be useful in considering the graph-isomorphism problem. [123] It was shown in Reference [123] that if you produce two graph states [124,125] from isomorphic graphs, then their equal-angle slices are equivalent, since taking the equal-angle slice can be thought of as an analogy of disregarding qubit, and in this case node, order. It was further shown in Reference [123] that non-isomorphic graphs have distinguishable equal-angle slices. However the results get more difficult to analyze as the number of qubits increase, as the detail in the Wigner function gets more and more intricate as the size of the Hilbert space increases. Apart from the equal-angle slice, there are many choices one may prefer to use, depending on whatever characteristic correlations may arise within the system of interest. And some of these different choices, along with the equal-angle slice, will be considered in more detail here.
Starting with the simplest case of a composite system, we will consider bipartite entanglement between two qubits. Two qubits famously exhibit entanglement in the well-known maximally entangled Bell states, explicitly where |↑↑ is shorthand for |↑ ⊗ |↑ . Apart from the equal-angle slice, there are other slices one can take to gain an understanding of the quantum correlations present within the system. One such example was shown in Reference [42], where the slice φ 1 = φ 2 = 0 was taken and the resulting θ 1 and θ 2 degrees of freedom were plotted against each other. An example of this slice for |Φ − and |Ψ + from Reference [42] is shown in Figure 3 (a). Two-qubit Bell-type entanglement in the Wigner function manifests in this slice as oscillations between negative and positive quasi-probabilities in phase space. Figure 3 (a) shows theoretical results for what should be measured at 81 points in phase space for these two states, along with what was actually measured using IBM's qx.
There are further ways to view phase space for bipartite entanglement, for instance the approach in works such as Reference [126] consider bipartite entanglement between two states that have larger Hilbert space, demonstrating the entanglement between two Bose-Einstein condensates. In Reference [126], the authors considered the Wigner function of these states by taking a marginal Wigner function and a conditional Wigner function. The marginal Wigner function is simply the Wigner function for one of the states, which can be calculated for either the reduced density matrix or by integrating out one set of degrees of freedom where ρ 1 = Tr 2 [ρ] is the reduced density matrix for the first subsystem. Note that for a maximally entangled Bell states, taking the marginal Wigner function results in a completely mixed state that has a value of 1/2 everywhere, see the video in the supplementary material of Reference [42] for an example of this. The conditional Wigner function is calculated by projecting a state onto the basis states of one of the subsystems, more details can be found in Reference [126]. As the number of qudits increases, the choices of slice becomes more difficult, due to the increasing number of degrees of freedom. In this case, the n-qudits can be represented by n marginal distributions, however if there is any entanglement present these become less useful. As a default, the equal-angle slice is the simplest starting point to getting an idea of the correlations that manifest in phase space.
Examples of the equal-angle slice have been shown to be useful in characterizing multi-qubit states in many experimental settings, where there have been examples ranging from three-qubit states up to 20-qubit multi-component Schrödinger cat states. In Reference [127], the authors created entangled states of three qubits, generated by two photonic qubits and a further path-encoded qubit. Through this method they generated a cluster GHZ state and a W state.
As we move up to a five-qubit state in Figure 3 (b) and Figure 3 (c), we can see there are now five oscillations. When plotting on a sphere, the variation in the φ degrees of freedom for θ = π/2 lies on the equator. This gives rise to the name 'equatorial slice'. The equatorial slice on it's own has been plotted in Figure 3 (b) with both theoretical and experimental results plotted, where there has been a least-squares best-fit curve plotted for the experimental data. By looking for these oscillations around the equator, this test serves as a verification of the generation of GHZ-type entanglement. Although there is clearly a lot of noise in the machine, the oscillating behaviour shows that a considerable about of GHZ-type entanglement remains in the state.  Figure 3 (b) considers this equatorial slice taking measurement with the Wigner function parity operator. However, this may not always be the most desirable choice to look for these oscillations, since as the number of qubits increases, the amplitude of the oscillations steadily decreases. Alternatively, one can try different choices of generalized parity operator, such as in Reference [52] where σ z was taken, as can be seen in Figure 3 (f). Experimental results for taking the equatorial slice using this parity is shown for 10, 12, 14, 16, 17, and 18 qubits in Figure 3 (f), showing good agreement between theory and experiment.
We note that the use of the σ z operator as a parity operator does not satisfy the Stratonovich-Weyl correspondence, requiring some extra assumptions for informational completeness, see Equation (57) and the discussion around it. However, if we're measuring the state of a quantum system it is safe to make the relevant assumptions -that the operator is trace 1. The kernel σ z (θ, φ) = U 1 2 (θ, φ)σ z U 1 2 (θ, φ) † can then be considered as a generalized observable, creating a characteristic function that has a Hermitian kernel, where the Pauli operators σ x , σ y , and, σ z are yielded by setting (θ, φ) = (−π/2, 0), (θ, φ) = (−π/2, −π/2), and (θ, φ) = (0, 0) respectively. This allows a direct comparison to verification and fidelity estimation protocols, for example References [128,63] and the fidelity estimation results in Section 4.1.
The experiments from Reference [52] were performed using a 20-qubit machine that generated entangled states by employing a single-axis twisting Hamiltonian. [65] The single-axis twisting Hamiltonian also creates different multi-component Schrödinger cat states throughout evolution before reaching the bipartite GHZ state. Figure 3

Hybrid continuous-discrete-variable quantum systems
We now look into ways to measure the phase-space distribution of states that are hybridizations of continuousand discrete-variable systems. Many experiments and processes that consider the state of qubits also implicitly involve the presence of quantum light. Likewise, many optical experiments also involve the interaction with some kind of atom, whether or not artificial. This interaction is the foundation of quantum electrodynamics and circuit quantum electrodynamics. Where the advancements made in this area has directly contributed to the current state of quantum technologies. In fact, the experiments on qubit states mentioned in the last section are performed by controlling the atomic state with quantum light pulses. Whether individually for each qubit or using a bus model to control all the qubits at once.
Conversely, the continuous-variable Wigner function can serve as the position and momentum of an atom, so creating a composite function of discrete a continuous variables may be interpreted as the atomic state depending on its position and momentum. This was shown useful to consider in Reference [94] when looking at the full Wigner function of an atom. By taking this composite approach there is much insight to be gained, from properties of light-matter interaction to other phenomena such as spin-orbit coupling. In this section we will focus on the former application as it more directly applies to use in quantum technologies.
There are many ways phase-space methods can be used in discerning correlations that arise in the interaction between continuous-and discrete-variable systems. Like in Equation (150) for two qudits, it is usual to simply consider the marginal functions for each system individually, tracing out either the discretevariable system or the continuous-variable system. This is also implicitly the case in measuring the Wigner function in qubit experiments when using quantum states of light to manipulate qubits. In theory, the qubit state should not have any correlations with the optical state, however there will always be some leakage, interpreted as decoherence in the qubit.
On the other side, it is typical to just consider the optical state in a light-matter interaction. An important example is the generation of an optical Schrödinger cat state. In such experiments, it is normal to entangle light with a two-level atom, by doing this one can then split the light into the superposition of two coherent states, generating a Schrödinger cat state. [41] . At the point of measurement of the optical state, the goal is to produce separable atomic and optical states, and therefore minimizing decoherence in the Wigner function. Like in the previous case, there may be some correlations that exist between the qubit and the field, therefore it may be more desirable to have an idea of the correlations that exist between the two systems.
Different examples of how one can plot such hybrid states are shown in Figure 4. We present some theoretical methods displayed in Figure 4. We will then discuss how some of these methods have been implemented experimentally and how the future possibilities of measurement of hybrid quantum systems.
We will first give examples of the results when taking the marginal Wigner functions. For the separable states, |0 |↑ and |0 |↓ , shown in Figure 4 (a) and (b), the marginal functions are as expected, in both cases the HW Wigner function is the vacuum state and the qubit is in the |↑ and |↓ state respectively. Figure 4 (c) show the marginal Wigner functions for the state (|0 |↓ + |1 |↑ )/ √ 2 -a hybrid analog of a Bell state. As a result, the quantum correlations are non-local and the marginal Wigner functions are mixed states Where the qubit state is maximally mixed, which can be seen in Figure 4 (c)(ii) with a constant value over the whole distribution.
Note that we have chosen an alternative projection to display the qubit state. We have used the Lambert azimuthal equal-area projection, [129] in order to display the qubit completely. It is similar to a stereographic projection with the north pole in the center, however the south pole is now projected to a unit circle, rather than to infinity. The function is then distributed around the disk to preserve the area -which is an important property when considering probability distributions. this results in the equator of the sphere existing on the concentric circle with radius 1/ √ 2. Such an approach was taken in Reference [68] by considering the Tavis-Cummings interaction [66] between an optical cavity and an atomic j = 5/2 state initially in an atomic Schrödinger cat state. By considering the Wigner functions of both systems, it is then insightful to visualize the transfer of the Schrödinger cat states between the two systems. Demonstrating how one can overcome decoherence in quantum systems by coupling a large atomic state to an optical mode. However, this lacks the consideration of the non-local correlations between the two systems, and any entangled state will result in Wigner functions similar to Figure 4 (c)(i) and (ii).
To get a sense of these non-local correlations a hybridization approach is needed. One approach has been to hybridize a density matrix for a qubit with the Wigner function for the field mode. This results in a 2 × 2 grid, representing the entries of the qubit density matrix, where on each entry of the density matrix the corresponding Wigner function for the field mode is shown. We can label this Wigner function w(i, j), where i, j ∈ {0, 1}, such that w(i, j) = Tr [(|i j| ⊗ Π(α)) ρ], where ρ is the state for the full system. Such an approach reveals much more about the non-local correlations between the two systems, and can be found in works such as Reference [130,131]. Examples of this approach are shown in Figure 4(iii). Note that the real value has been taken, which doesn't make a difference for the first two states, however (|0 |↓ + |1 |↑ )/ √ 2 has an imaginary component on the off-diagonal entries that has been ignored. Since we're concerned with eigenstates of σ z , the local terms are shown in diagonal entries and the non-local correlations can be seen in the off-diagonal entries, which is clearly seen in Figure 4 (c)(iii).
Alternatively one could take Pauli measurements of the qubit, and then calculate the HW Wigner function in that basis, the resulting distribution is shown in Figure 4(iv) and can be considered as a hybridization of the discrete Weyl function with the HW Wigner function. Similarly to the density matrix approach this produces a 2 × 2 grid, where this time the elements correspond to Wigner function, yielding a joint quasiprobability distribution. For each state, we take an identity measurement, this is equivalent to taking the marginal HW Wigner function, where the equivalence can be considered with Figure 4 (c)(i).
Such a method has gained much interested in practice, where these exact states were considered in Reference [132], where the authors took measurements in the Pauli basis at different points in position-momentum phase space. By building up measurement statistics each point in phase space was then represented by the expectation value of the given Pauli operator at this point in phase space. Similar results that consider entanglement between a qubit and a Schrödinger cat states have also shown how this method is useful for characterizing entanglement. [133,134] Alternatively we can consider a different method to hybridize Wigner functions for the two systems. One could hybridize the Wootters Wigner function with the HW Wigner function in a similar way to coupling the density matrix, as was done in Reference [135]. However here we want to consider hybridizing the Stratonovich-Kernel Wigner function with the the HW Wigner function. [94,95] At each pointα we can plot the Wigner function for the qubit, producing a lattice of Wigner functions where W (α, θ, φ) = Tr Π(α) ⊗ Π 1 2 (θ, φ) ρ . 3)/2. We can see that for the two separable states in (a) and (b) that the two marginal functions are pure states, conversely we see in (c) that the entangled state produces mixed states in the marginal functions, where the qubit has a uniform value over the distribution. Note that the qubit has been plotted in a Lambert azimuthal area preserving map, see text for details. In (iii) we have plotted the density matrix for the qubit, where in each entry we show the HW Wigner function. Similarly in (iv) we have a composition of the discrete Weyl function where in each entry the HW Wigner function has been plotted, note that this corresponds to taking Pauli measurements of the qubit and then plotting the HW Wigner function, as a result the marginal HW Wigner function is appears as the identity measurement of the qubit. In both (iii) and (iv) the colorbar values are η = 2. Finally in (v) we show a hybrid representation developed in Reference [94,95] which have colorbar values η = 1 + √ 3 We can then add transparency to each of the spheres according to the maximum value of the Wigner function at that pointα, max θ,φ W (α, θ, φ), resulting in an envelope over the distribution. Examples of this are given in Figure 4(a) -(c) (v), where for the separable states we can see the envelope of the vacuum state. At every discrete point in position and momentum phase space we then see the |↑ state, in Figure 4(a)(v), and the |↓ state in Figure 4(b)(v). Like with the Wootters Wigner function, when there are non-local correlations between the two systems, we see a dependence between the degrees of freedom. Where in Figure 4(c)(v) we can see a radial dependence on the state of the qubit. All these methods provide different snapshots into such coupled systems and in combination build up an understanding into the of these types of non-local correlations. We also see certain overlap between different methods of plotting such systems, which further allows us to get a sense of the overall picture of coupled quantum systems in phase space.

Conclusions and Outlook
In this review we have explored the phase-space formulation of quantum mechanics, considering how different systems can be represented in phase space. Our focus was providing the framework to represent discrete-variable systems in phase space, in particular the different ways single-and multi-qubit states can be considered.
We provided the general framework to generate an informationally complete phase-space distribution from any arbitrary operator, as well as how to transform between any of the phase-space representations and how the dynamics of quantum systems can be modelled in phase space. This was followed by examples of discrete systems in phase space, with a focus on the construction on qubit and multi-qubit states. We have provided a link between between the Stratonovich kernel and the Wootters kernel for a qubit state. Given the resulting structure of a simplex inside a sphere, this provides a path into alternative constructions of discrete Wigner function for qudits and symmetric informationally complete states, where we can consider an N 2 -dimensional simplex inside the SU(N ) N 2 − 1-dimensional surface.
From the construction on Wigner functions for discrete quantum states, we then explored how such distributions can be useful for applications to quantum technologies. Exploring how well-known metrics in quantum information can be written in the language of phase space, as well as methods to use them in experimental settings. There are still more metrics than can be translated into phase space, and perhaps even more that are peculiar to the phase-space representation, such as the negative volume of the Wigner function.
We then considered direct measurement of quantum systems in phase space, beginning with the phasespace reconstruction with single-qubit states before treating more complicated composite systems. Where we first looked at multi-qubit states, where the dimension of the Hilbert space increases as 2 N , whereas the degrees of freedom for a phase-space distribution only scale as 2N for N qubits. Therefore such methods can provide computational simplifications in larger systems -especially if we are only interested in finding certain types of correlations within a larger state.
We then moved onto the current work of hybrid quantum systems. Although the distributions considered here and simple single-qubit states these provide an outline into how come can adapt phase-space methods for composite systems with heterogeneous degrees of freedom. We note that there is also some exploration into how this can be extended to consider multi-qubit states coupled to a field mode via a Tavis-Cummings coupling in Reference [53]. Conversely we can consider three continuous-variable distributions to each qubit, acting as the position and momentum in three dimensions to each atom, such an approach was considered in Reference [94] when treating the states of atoms.