Probing nonclassicality with matrices of phase-space distributions

We devise a method to certify nonclassical features via correlations of phase-space distributions by unifying the notions of quasiprobabilities and matrices of correlation functions. Our approach complements and extends recent results that were based on Chebyshev's inequality [Phys. Rev. Lett. 124, 133601 (2020)]. The method developed here correlates arbitrary phase-space functions at arbitrary points in phase space, including multimode scenarios and higher-order correlations. Furthermore, our approach provides necessary and sufficient nonclassicality criteria, applies to phase-space functions beyond $s$-parametrized ones, and is accessible in experiments. To demonstrate the power of our technique, the quantum characteristics of discrete- and continuous-variable, single- and multimode, as well as pure and mixed states are certified only employing second-order correlations and Husimi functions, which always resemble a classical probability distribution. Moreover, nonlinear generalizations of our approach are studied. Therefore, a versatile and broadly applicable framework is devised to uncover quantum properties in terms of matrices of phase-space distributions.


Introduction
Telling classical and quantum features of a physical system apart is a key challenge in quantum physics. Besides its fundamental importance, the notion of (quantum-optical) nonclassicality provides the basis for many applications in photonic quantum technology and quantum information [1,2,3,4,5]. Nonclassicality is, for example, a resource in quantum networks [6], quantum metrology [7], boson sampling [8], or distributed quantum computing [9]. The corresponding free (i.e., classical) operations are passive linear optical transformations and measurement. By exceeding such operations, protocols which utilize nonclassical states can be realized. Furthermore, nonclassicality Martin Bohmann: martin.bohmann@oeaw.ac.at is closely related to entanglement. Each entangled state is nonclassical, and single-mode nonclassicality can be converted into two-and multi-mode entanglement [10,11,12].
Consequently, a plethora of techniques to detect nonclassical properties have been developed, each coming with its own operational meanings for applications. For example, quantumness criteria which are based on correlation functions and phase-space representations have been extensively studied in the context of nonclassical light [13,14].
The description of physical systems using the phase-space formalism is one of the cornerstones of modern physics [15,16,17]. Beginning with ideas introduced by Wigner and others [18,19,20,21], the notion of a phase-space distribution for quantum systems generalizes principles from classical statistical theories (including statistical mechanics, chaos theory, and thermodynamics) to the quantum domain. However, the nonnegativity condition of classical probabilities does not translate well to the quantum domain. Rather, the notion of quasiprobabilitiesi.e., normalised distributions that do not satisfy all axioms of probability distributions and particularly can attain negative values-was established and found to be the eminent feature that separates classical concepts from genuine quantum effects. (See Refs. [22,14] of thorough introductions to quasiproabilities.) In particular, research in quantum optics significantly benefited from the concept of phase-space quasiprobability distributions, including prominent examples, such as the Wigner function [19], the Glauber-Sudarshan P function [23,24], and the Husimi Q function [25]. In fact, the very definition of nonclassicality-the impossibility of describing the behaviour of quantum light with classical statistical optics-is directly connected to negativities in such quasiprobabilities, more specifically, the Glauber-Sudarshan P function [26,27]. Because of the general success of quasiprobabilities, other phase-space distributions for light have been conceived [28,29,30], each coming with its own advantages and drawbacks. For example, squeezed states are represented by nonnegative (i.e., classical) Wigner functions although they form the basis for continuous-variable quantum information science and technology [31,32,33], also having a paramount role for quantum metrology [34,35].
Another way of revealing nonclassical effects is by using correlation constraints which, when violated, witness nonclassicality. Typically, such conditions are formulated in terms of inequalities involving expectation values of different observables. Examples in optics are photon anti-bunching [36,37,38] and sub-Poissonian photon-number distributions [39,40], using intensity correlations, as well as various squeezing criteria, being based on field-operator correlations [41,42,43,44]. They can follow, for instance, from applying Cauchy-Schwartz inequalities [45] and uncertainty relations [46], as well as from other violations of classical bounds [47,48,49]. Remarkably, many of these different criteria can be jointly described via so-called matrix of moments approaches [50,51,52,53,54]. However, each of the mentioned kinds of nonclassicality, such as squeezed and sub-Poissonian light, requires a different (sub-)matrix of moments, a hurdle we aim at overcoming.
Over the last two decades, there had been many attempts to unify matrix-of-moment-based criteria with quasiprobabilities. For example, the Fourier transform of the P function can be used, together with Bochner's theorem, to correlate such transformed phase-space distributions through determinants of a matrix [55,56], being readily available in experimental applications [57,58,59,60], and further extending to the Laplace transformation [61]. Furthermore, a joint description of field-operator moments and transformed phase-space functions has been investigated as well [62]. Rather than considering matrices of phasespace quasiproabilities, concepts like a matrix-valued distributions enable us to analyzed nonclassical hybrid systems [63,64]. Very recently, a first successful strategy that truly unifies correlation functions and phase-space functions has been conceived [65]. However, these first demonstrations of combining phasespace distributions and matrices of moments are still restricted to rather specific scenarios.
In this contribution, we formulate a general framework for uncovering quantum features through correlations in phase-space matrices which unifies these two fundamental approaches to characterizing quantum systems. By combining matrix of moments and quasiprobabilities, this method enables us to probe nonclassical characteristics in different points in phase space, even using different phase-space distributions at the same time. We specifically study implications from the resulting second-and higherorder phase-space distribution matrices for single-and multimode quantum light. Furthermore, a direct measurement scheme is proposed and non-Gaussian phase-space distributions are analyzed. To benchmark our method, we consider a manifold of examples, representing vastly different types of quantum features. In particular, we show that our matrixbased approach can certify nonclassicality even if the underlying phase-space distribution is nonnegative. In summary, our approach renders it possible test for nonclassicality by providing easily accessible nonclassicality conditions. While previously derived phase-space-correlation conditions [65] were restricted to single-mode scenarios, the present approach straightforwardly extends to multimode cases. In addition, our phase-space matrix technique includes nonclassicality-certification approaches based on phase-space distributions and matrices of moments as special cases, resulting in an overarching structure that combines both previously separated techniques.
The paper is structured as follows. Some initial remarks are provided in Sec. 2. Our method is rigorously derived and thoroughly discussed in Sec. 3. Section 4 concerns several generalizations and potential implementations of our toolbox. Various examples are analyzed in Sec. 5. Finally, we conclude in Sec. 6.

Preliminaries
In their seminal papers [23,24], Glauber and Sudarshan showed that all quantum states of light can be represented diagonally in a coherent-state basis through the Glauber-Sudarshan P distribution. Specifically, a single-mode quantum state can be expanded asρ where |α denotes a coherent state with a complex amplitude α. Then, classical states are identified as statistical (i.e., incoherent) mixtures of pure coherent states, which resemble the behavior of a classical harmonic oscillator most closely [66,67]. For this diagonal representation to exist for nonclassical states as well, the Glauber-Sudarshan distribution has to exceed the class of classical probability distributions [26,27], particularly violating the nonnegativity constraint, P 0. This classification into states which have a classical correspondence and those which are genuinely quantum is the common basis for certifying nonclassical light.
As laid out in the introduction, nonclassicality is a vital resource for utilizing quantum phenomena, ranging from fundamental to applied [6,7,8,9]. In this context, it is worth adding that, contrasting other notions of quantumness, nonclassicality is based on a classical wave theory. That is, it is essential to discern nonclassical coherence phenomena from those which are accessible with classical statistical optics, as formalized through Eq. (1) with P ≥ 0. See, e.g., Ref. [68] for a recent experiment that separates classical and quantum interference effects in such a manner. For instance, free operations, i.e., those maps which preserve classical states, do include beam splitter transformations, resulting in the generation of entanglement from single-mode nonclassical states via such a free operation [10,11,12] that is vital for many quantum protocols.

Phase-space distributions
Since the Glauber-Sudarshan distribution can be a highly singular distribution (see, e.g., Ref. [69]), generalized phase-space functions have been devised. Within the wide range of quantum-optical phasespace representations, the family of s-parametrized distributions [29,30] is of particular interest. Such distributions can be expressed as where colons indicate normal ordering [70] is the displaced photon-number operator, written in terms of bosonic annihilation and creation operators,â andâ † , respectively. It is worth recalling that the normal ordering acts on the expression surrounded by the colons in such a way that creation operators are arranged to the left of annihilation operators whilst ignoring commutation relations. Note that, for convenience, we parametrize distributions via the width parameter σ, rather than using s. Both are related via From this relation, we can identify the Husimi function, Q(α) = P (α; 1), for s = −1 and σ = 1; the Wigner function, W (α) = P (α; 2), for s = 0 and σ = 2; and the Glauber-Sudarshan function, P (α) = P (α; ∞), for s = 1 and σ = ∞. Whenever a phase-space distribution contains a negative contribution, i.e., P (α; σ) < 0 for at least one pair (α; σ), the underlying quantum state is nonclassical [26,27]. In such a case, the distribution P (α; σ) refers to as a quasiprobability distribution which is incompatible with classical probability theory. Nonetheless, for any σ ≥ 0 and any state, this function represents a real-valued distribution which is normalized, P (α; σ) = P (α; σ) * and d 2 α P (α; σ) = 1. In addition, it is worth mentioning that the normalization of the state is guaranteed through the limit

Matrix of moments approach
Besides phase-space distributions, a second family of nonclassicality criteria is based on correlation functions; see, e.g., Refs. [71,72] for introductions. For this purpose, we can consider an operator function f = f (â,â † ). Then, holds true for all P ≥ 0. Now, one can expandf in terms of a given set of operators, e.g.,f = i c iÔi , resulting in :f †f : = i,j c * i c j :Ô † iÔ j : . Furthermore, this expression is nonnegative [cf. Eq. (5)] iff the matrix ( :Ô † iÔ j : ) i,j is positive semidefinite. This constraint can, for example, be probed using Sylvester's criterion [73] which states that a Hermitian matrix is positive-definite if and only if all its leading principal minors have a positive determinant. It is worth mentioning that Eq. (5) defines the notion of a nonclassicality witness, where :f †f : < 0 certifies nonclassicality.
The above observations form the basis for many experimentally accessible nonclassicality criteria, such as using basis operators which are powers of quadrature operators [50], photon-number operators [44], and general creation and annihilation operators [71,72]. See Refs. [13] for an overview of momentbased inequalities. In the following, we are going to combine the phase-space distribution technique with the method of matrices of moments to arrive at the sought-after unifying approach of both techniques.

Matrix of phase-space distributions
Both phase-space distributions and matrices of moments exhibit a rather dissimilar structure when it comes to formulating constraints for classical light. Consequently, a full unification of both approaches is missing to date, excluding the few attempts mentioned in Sec. 1. In this section, we bridge this gap and derive a matrix of phase-space distributions which leads to previously unknown nonclassicality criteria, also overcoming the limitations of earlier methods.

Derivation
For the purpose of deriving our criteria, we consider an operator functionf = i c i exp[−σ ini (α i )]. Then, the normally ordered expectation value off †f can be expanded as Based on the above relation, we may define two matrices, one for classical amplitudes, and one for the quantum-optical expectation values, which can be expressed in terms of phase-space distributions using Eq. (2). Specifically, M (q) corresponds to a matrix of phase-space distributions. Moreover, the fact that the normally ordered expectation value off †f is nonnegative for classical light [Eq. (5)] is then identical to the entry-wise product (i.e., the Hadamard product •) of both matrices being positive semidefinite, defining our phase-space matrix M . For classical light, all principal minors of M have to be nonnegative according to Sylvester's criterion. Conversely, the violation of this constraint certifies a nonclassical state, where M is defined through arbitrary small or large sets of parameters σ i and σ j and coherent amplitudes α i and α j . Therefore, inequality (10) enables us to formulate various nonclassicality conditions which correlate distinct phase-space distributions as it typically only done for matrix-of-moments-based techniques when using different kinds of observables. We finally remark that the expression in Eq. (10) resembles a nonlinear nonclassicality witnessing approach. As a first example, we may explore the first-order criterion, i.e., a 1 × 1 matrix of quasiprobabilities. Selecting arbitrary σ-parameters and coherent amplitudes, i.e., (α 1 ; σ 1 ) = (α; σ), we find the following restriction for classical states [cf. Eq. (10)]: This inequality reflects the fact that finding negativities in a parametrized phase-space distribution P (α; 2σ) is sufficient to certify nonclassicality. Also recall that we retrieve the Glauber-Sudarshan distribution in the limit σ → ∞. Since the nonnegativity of this distribution defines the very notion of a nonclassical state [26,27], we can conclude from this examples that our approach is necessary and sufficient for certifying nonclassicality. However, the Glauber-Sudarshan distribution has the disadvantage of being a highly singular for many relevant nonclassical states of light and, thus, hard to reconstruct from experimental data. Consequently, it is of practical importance (see Secs. 4 and 5) to consider higher-order criteria beyond this trivial one.

Second-order criteria
We begin our consideration with an interesting second-order case. We chose (α 1 ; σ 1 ) = (0; 0) and (α 2 ; σ 2 ) = (α; σ). This yields the 2 × 2 phase-space matrix M = 1 : exp(−σn(α)): : exp(−σn(α)): : exp(−2σn(α)): . (12) Up to a positive scaling, the determinant of this matrix results in the following nonclassicality criterion: In particular, we can set σ = 1 to relate this condition to the Wigner and Husimi functions, leading to W (α) − 2πQ(α) 2 < 0. This special case of our general approach has been recently derived using a very different approach, using Chebyshev's integral inequality [65]. There it was shown that, by applying the inequality (13) for σ = 1, it is possible to certify nonclassicality even if the Wigner function of the state under study is nonnegative. In this context, remember that the Husimi function, Q(α) = α|ρ|α /π, is always nonnegative, regardless of the stateρ. Beyond this scenario, we now study a more general 2 × 2 phase-space matrix M . For an efficient description, it is convenient to redefine transformed parameters as Note that these parameters relate to the two-body problem. That is, the quantities in Eq. (14a) define the relative position and barycenter in phase-space, respectively; and the two quantities in Eq. (14b) resemble the reduced and total mass in a mechanical system, respectively. In this alternative parametrization, the two matrices, giving the total phase- 2 1 , and (15a) If this determinant is negative for the state of light under study, its nonclassicality is proven. In terms of phase-space distributions, this condition can be also recast into the form Interestingly, this nonclassicality criterion correlates different points in phase space for different distributions, P (α 1 ; 2σ 1 ) and P (α 2 ; 2σ 2 ), with a phase-space distribution with the total width Σ at the barycenter A of coherent amplitudes, P (A; Σ).

Higher-order cases
The next natural extension concerns the analysis of higher-order correlations. Clearly, one can obtain an increasingly large set of nonclassicality tests with an increasing dimensionality of M , determined by the number of pairs (α i ; σ i ). In order to exemplify this potential, let us focus on one specific 3 × 3 scenario and more general scenarios for specific choices of parameters.
Another higher-order matrix scenario corresponds to having identical coherent amplitudes, i.e., α i = α for all i. In this case, we find that the two Hadamardproduct components of the matrix M simplify to thus resulting in M = M (q) . Therefore, we can formulate nonclassicality criteria which correlate an arbitrary number of different phase-space distributions, defined via σ i , at the same point in phase-space, α. Analogously, one can consider a scenario in which all σ parameters are identical, σ i = σ. Then, we get Consequently, we obtain nonclassicality criteria which correlate an arbitrary number of different points in phase space, α i , for a single phase-space distribution, parametrized by σ.

Comparison with Chebyshev's integral inequality approach
As mentioned previously, a related method based on Chebyshev's integral inequality has been introduced recently [65]. It also provides inequality conditions for different phase-space distributions. The nonclassicality conditions based on Chebyshev's integral inequality take the form To compare both approaches, let us discuss their similarities and differences.
In its simplest form, involving only σ 1 and σ 2 , the condition in Eq. (22) resembles the tests based on the 2 × 2 matrix in Eq. (12). In particular, for the case σ 1 = σ 2 = σ both methods yield the exact same conditions. For σ 1 = σ 2 such an agreement of both methods cannot be found because of the inherent symmetry of the phase-space matrix approach, M = M † , which stems from its construction via a quadratic form; cf. Eq. (6). Also, for more general, higher-order conditions, i.e. D > 2, such similarities cannot be found either. Conditions of the form in Eq. (22) consist of only two summands. The first term is a single phasespace function with the width parameter Σ which is associated to the highest σ parameter involved in the inequality. The second term is a product of D phasespace distributions, each individual distribution has a width parameter σ i , together bound by the condition Σ = D i=1 σ i . By comparison, our phase-space matrix approach yields, in general, a richer and more complex set of higher-order nonclassicality tests, such as demonstrated in Sec. 3.3.
Let us point out further differences between the two approaches. Firstly, we observe that the inequalities based on Chebyshev's integral inequality only apply to one single point in phase space. In contrast, the phase-space matrix method devised here includes conditions that combine different points in phase space; cf. Eq. (6). Secondly, Chebyshev's integral inequality approach cannot be extended to multimode settings. Such a limitation does not exist for the matrix approach either, as we show in the following Sec. 4.1. We conclude that both the technique in Ref. [65] and our phase-space matrix approach for obtaining phase-space inequalities yield similar secondorder conditions but, in general, give rise to rather different nonclassicality criteria. In particular, the phase-space matrix framework offers a broader range of variables-be it coherent amplitudes or widthsthat lead to a richer set of nonclassicality conditions.

Extended relations to nonclassicality criteria
To finalize our first discussions we now focus on the relation to matrices of moments. Previously, we have shown that, already in the first order, our criteria are necessary and sufficient to verify the nonclassicality, and we discussed our method in relation to Chebychev's integral inequality. Furthermore, indirect techniques using transformed phase-space functions, such as the characteristic function [62] and the two-sided Laplace transform [61], have been previously related to moments. Thus, the question arises what the relation of our direct technique to such matrices of moments is.
For showing that our framework includes the matrix of moments technique, we may remind ourselves that derivatives can be understood as a linear combination, specifically as a limit of a differential quotient, . This enables us to write [75] a †mân = σ −(m+n) ∂ m α ∂ n α * e σ|α| 2 :e −σn(α) : α=0 and σ=0 , (23) expressing arbitrary momentsâ †mân via linear combinations of the normally ordered operators that represent σ-parametrized phase-space distributions. Thus, in the corresponding limits, we can identify the operatorf in Eq.
In conclusion, we find that our necessary and sufficient methodology not only includes nonclassicality criteria based on phase-space functions themselves [cf. Eq. (11)], but it also includes the technique of matrices of moments as a special case. In a hierarchical picture, this means that our family of nonclassicality criteria, including arbitrary orders of σ-parametrized phase-space functions, encompasses both negativities of phase-space functions and matrices of moments. Because of the above relation, the order of moments that is required to certify nonclassicality also sets an upper bound to the size of the matrix of phase-space distributions so that it certifies nonclassicality. Therefore, our approach unifies and subsumes both earlier types of nonclassicality conditions.

Generalizations and implementation
In this section, we generalize our approach to arbitrary multimode nonclassical light and propose a measurement scheme to experimentally access the matrix of phase-space distributions. In addition, we show that our approach applies to phase-space distributions which are no longer limited to σ parametrizations and relate these findings to the response of nonlinear detection devices.

Multimode case
After our in-depth analysis of single-mode phasespace matrices, the multimode case follows almost straightforwardly. For the purpose of such a generalization, we consider N optical modes, represented via the annihilation operatorsâ m for m = 1, . . . , N and extending to the displaced photon-number operatorsn m (α (m) ) = (â m − α (m) ) † (â m − α (m) ). Now, σ-parametrized multimode phase-space functions can be expressed as where we allow for different s parameters for each mode, with s (m) = 1 − 2/σ (m) [Eq. ≥ 0 (25) holds true for classical light and for any dimension (or order) of the multimode matrix M and any sigma and alpha values. Conversely, det(M ) < 0 is a nonlinear witness of multimode nonclassicality. Similarly to the single-mode case, an increasingly large matrix M with increasingly dense sets of parameters for the various alpha and sigma values then enables one to probe the nonclassicality of arbitrary multimode states.
Since we have already exemplified various scenarios for single-mode phase-space correlations, in the following, we restrict ourselves to a particular multimode case. Specifically, we focus on two optical modes and a 3 × 3 phase-space matrix M is, where quasiprobabilities as a function of single-mode parameters indicate marginal phase-space distributions. Adopting a notation of pairs of coherent amplitudes and widths, M is thus defined via the following two-mode parameters: (α 1 , α 3 ) = (0, α (2) ; 0, σ). In particular, we can express the nonclassicality constraint from the determinant of M [74] for σ = 1 via joint and marginal Wigner and Husimi functions, Violating this inequality verifies the nonclassicality of the two-mode state under study.

Direct measurement scheme
The reconstruction of phase-space distributions can be a challenging task [76]. For this reason, we are going to devise a directly accessible setup to infer the phase-space matrix. See Fig. 1 for an outline which is based on the approaches in Refs. [77,72,78]. For convenience, we restrict ourselves to a single optical mode; the extension to multiple modes follows straightforwardly. That is, each of the multiple modes can be detected individually by a correlationmeasurement setup as depicted in Fig. 1. Furthermore, it is noteworthy that our phase-space matrix approach is not limited to this specific measurement scheme proposed here and generally applies to any detection scenario which allows for a reconstruction of quasiprobability distributions. For the setup in Fig. 1, we begin our considerations with a coherent state |β , representing our signalρ = |β β|. Firstly, we split this signal equally into 2 modes, resulting in a two-mode coherent state |β/ √ 2, β/ √ 2 . In addition, local oscillator states are prepared, |β i and |β j for each mode. Each of the two signals is then mixed with its local oscillator on a |t| 2 :|r| 2 beam splitter, where |t| 2 + |r| 2 = 1. One output of each beam splitter is discarded, namely the lower and upper one for the top and bottom path in Fig. 1, respectively. This results in the input-output relation which is then detected as follows. Each of the resulting modes is measured with a detector or detection scheme based on photon absorption, thus being described by a positive operatorvalued measure (POVM) which is diagonal in the photon-number representation [79]. Consequently, one or a combination of detector outcomes (e.g., in a generating-function-type combination [80]) corresponds to a POVM element of the form Π(n) = :e −Γ(n) :.
Using |m m| = :e −nnm /m!: for an m-photon projector, this means that we identify ∞ m=0 π m |m m| = :e −n m=0 π mn m /m!: = Π(n) = :e −Γ(n) :, where the eigenvalues π m corresponds to the Taylor expansion coefficients of the function z → exp[z − Γ(z)]. Accordingly, the function Γ(n) models the detector response [79,70]. Finally, the correlation measurement of this response for our coherent signal states takes the form Now it is convenient to defineΓ(n) = Γ(|t| 2n /2) and for all LO choices i and, similarly, for j. Furthermore, we generalize this treatment to arbitrary states, ρ = d 2 β P (β)|β β|, using the Glauber-Sudarshan representation [Eq.
(1)]. Therefore, the correlations measured as described above [Eq. (28)] obey M i,j = :e −Γ(n(αi)) e −Γ(n(αj )) : , which corresponds to a directly measured phase-space matrix element, e.g., for a linear detector responsẽ Γ(n) = σn. The other way around, we can choosê f = i c i exp(−Γ(n(α i ))) for the general classicality constraint in Eq. (5), even for nonlinear detector responses. Then, the matrix of phase-space distribution approach applies, regardless of a linear or nonlinear detection model. (See also Refs. [81,80] in this context.) As an example, we consider a case with two single on-off click detectors (represented by Π(n) in Fig. 1) with a non-unit quantum efficiency η det and a nonvanishing dark-count rate δ [82], which represents realistic detectors in experiments. In addition, we can introduce neutral density (ND) filters to attenuate the light that impinges on each detector. The POVM element for the no-click event in combination with the ND filters then readsΠ(n) = : exp(−(ηn + δ)):, where 0 ≤ η ≤ η det is a controllable efficiency. The measured correlation for this scenario takes the form Therein, the adjustable efficiency η i plays role of σ i . Also, the positive factor that includes the dark counts is irrelevant because it does not change the sign of the determinant of M , i.e., the verified nonclassicality.
In summary, the measurement layout in Fig. 1 enables us to directly measure the entries of our phasespace matrix M . As an experimental setup, this scheme also underlines the strong connection between correlations and their measurements and phase-space quasiprobabilities and their reconstruction. We may emphasize that all experimental techniques and components that are used in the proposed setup are readily available; see, e.g., the related quantum state reconstruction experiments reported in Refs. [83,80]

Generalized phase-space functions
The σ-parametrized phase-space distributions we considered so far are related to each other via convolutions with Gaussian distributions [28,29,30]. However, there are additional means to represent a state without relying on Gaussian convolutions only.
Such generalized phase-space function can be obtained from the Glauber-Sudarshan P function via P Ω (α) = d 2α P (α) Ω(α;α,α * ) = :Ω(α;â,â † ): (32) for a kernel Ω ≥ 0 [30,84]. The construction of this socalled filter or regularizing function Ω can be done so that the resulting distribution P Ω is regular (i.e., without the singular behavior known from the P function) and is positive semidefinite for all classical states [84]. For instance, a non-Gaussian filter Ω has been used to experimentally characterize squeezed states via regular distributions which exhibit negativities in phase space [85]; this cannot be done with s-parametrized quasiprobability distributions, which are either nonnegative or highly singular for squeezed states.

Examples and benchmarking
In the following, we apply our method of phase-space matrices to various examples and benchmark its performance. For the latter benchmark, we could consider different phase-space functions. Using the P function would be impractical as it is often a highly singular distribution. The Wigner function is regular and can exhibit negativities. But error estimations from measured data can turn out to be rather difficult because it requires diverging pattern functions [86,87] (see Ref. [88] for an in-depth analysis). Beyond those practical hurdles, we focus on the Q function here because, already in theory, it is always nonnegative. Thus, it is hard to verify nonclassical features based on this particular phase-space distribution. Additionally, the Q function is easily accessible in experiments and can be directly measured via the widely-used double-homodyne (aka, eight-port homodyne) detection scheme [70]. Nonetheless, we are going to demonstrate that, with our method, it is already sufficient for many examples to consider second-order correlations of Q functions. For this purpose, we use the condition in Eq. (17), which follows from the 2 × 2 matrix condition with σ 1 = σ 2 = 1/2. This special case of that condition then reads as Meaning that, when the correlations from Q functions at different points in phase space fall below the classical limit zero, nonclassical light is certified with the nonnegative family of Q distributions.
Moreover, since Q functions are nonnegative, the second term in Eq. (35) is subtractive in nature. Thus, it is sufficient to find a point α 1 in phase space for which Q(α 1 ) = 0 holds true-together with an α 2 with Q(α 2 /2) > 0, which has to exist because of normalization-in order to certify nonclassicality through Eq. (35). Setting α 1 = α, this leads to the simple nonclassicality condition Q(α) = 0, which applies to arbitrary quantum states. In Ref. [89], this specific condition has been independently verified as a nonclassical signature of non-Gaussian states. Here, we see that this nonclassical signature is indeed a corollary of our general approach. Furthermore, we remark that this condition only holds if the Q function is exactly zero. In experimental scenarios, in which errors have to be accounted for, it is infeasible to get this exact value. Therefore, the condition Eq. (35) is more practical as it allows us to certify nonclassicality through a finite negative value. Furthermore, this condition is applicable even if Q(α) = 0 does not hold true.

Discrete-variable states
We start our analysis of nonclassicality by considering discrete-variable states for a single mode. In the case of quantized harmonic oscillators, such as electromagnetic fields, a family of discrete-variable states that are of particular importance are number states |n . They represent an n-fold excitation of the underlying quantum field and show the particle nature of said fields, thus being nonclassical when compared to classical electromagnetic wave phenomena. However, photon-number states require Glauber-Sudarshan P distributions that are highly singular because they involve up to 2nth-order derivatives of delta distributions [70]. On the other hand, the Q function of photon-number states, is an accessible and smooth, but nonnegative function. Thus, by itself, it cannot behave as a quasiprobability which includes negative contributions that uncover nonclassicality. Except for vacuum, the Q |n function is zero for α = 0 and positive for all other arguments α [Eq. Note that this family of discrete-variable number states is rotationally invariant, rendering the phase of α 2 irrelevant. In Fig. 2, we visualize the results of our analysis. For all number states, we observe a successful verification of nonclassicality in terms of inequality Eq. (35). The singlephoton state shows the largest violation for this specific nonclassicality test, and the negativity of det(M ) decreases with the number of photons. A possible explanation for this behavior is that this condition is most sensitive towards the particle nature of the quantum states, being most prominent in the single excitation of the quantized radiation field. Again, let us emphasize that we verified nonclassciality via a matrix M of classical (i.e., nonnegative) phase-space functions.

Continuous-variable states
After studying essential examples of discretevariable quantum states, we now divert our attention to typical examples of continuousvariable states.
For this reason, we consider squeezed vacuum states which are defined as |ξ = (cosh r) −1/2 ∞ n=0 (−e iϕ tanh[r]/2) n (2n)!|2n /n!, for a squeezing parameter r = |ξ| and a phase ϕ = arg(ξ). Without a loss of generality, we set ϕ = 0. Squeezed states are widely used in quantum optical experiments and provide the basis of continuous-variable quantum information processing [31]. Their parametrized phase-space distributions are known to be either highly singular or nonnegative Gaussian functions (see, e.g., Refs. [32,69]). For example, the Q function of the states under study can be written as In the context of earlier discussions, note that this Q function is not zero for α = 0, or anywhere else. In Fig. 3, the left-hand side of inequality Eq. (35) is shown for the Q |ξ function of a squeezing parameter r. The points in phase space are determined by choosing α 1 = 0 and minimizing det(M ), being solved for α 2 = [(2/λ) ln[(1 + λ)/(1 + λ/2)]] 1/2 , where λ = tanh(r). We observe negative values as a direct signature of the nonclassicality of squeezed states. Remarkably, this is achieved using the same criterion that applies to photon-number states, typically vastly different correlation functions are required (using either photon numbers [39] or quadratures [41]). While inequality Eq. (35) is violated for any squeezing parameter r > 0, we see that there exists an optimal region of squeezing values around r = 0.6 (likewise, 5 dB of squeezing) for which the considered criterion is optimal. In particular, this shows that this condition works optimally in a range of moderate squeezing values and, thus, is compatible with typical experiments. We also want to recall that the Q |ξ are a Gaussian distributions which do not have any zeros in the phase space. Thus, criteria based on the zeros of the Husimi (a) (b) Figure 4: In plot (a), the two-mode Q function in Eq. (39) for the mixed and weakly correlated stateρ is depicted for |λ| 2 = 1/2 and phase-space points with Im(α (1) ) = Im(α (2) ) = 0. Part (b) visualizes the application of the nonclassicality inequality (40) to this state for N = 2 modes and for the parameter pairs (α 2 ) = (0, β). Nonclassicality is verified because of det(M ) < 0, and maximized for |α| and |β| around one.
Q function [89] cannot detect nonclassicality in this scenario. In contrast, our inequality condition can even certify this Gaussian nonclassicality, hence providing a more sensitive approach to detecting quantum light.

Mixed two-mode states
To further challenge our approach, we now consider a bipartite mixed state.
This state presents a particular challenge for nonclassicality verification because it shows only weak nonclassicality and quantum correlations. Namely, this state is not entangled, has zero quantum discord, and has classical marginal single-mode states (i.e., the partial traces tr 1 (ρ) = tr 2 (ρ) yield thermal states) [90]. However, it shows nonclassical photon-photon correlations [91,90,92]. The state's two-mode Q function can be computed using Gaussian functions and the phase averaging in Eq. (38), which gives where I 0 denotes the zeroth modified Bessel function of the first kind. See also Fig. 4(a) in this context. To apply our approach, whilst using Q functions only, we can directly generalize our criterion in Eq.
In Fig. 4(b), we apply the case N = 2 of this inequality to identify the nonclassicality ofρ for |λ| 2 = 1/2. Again, the same approach as used in both singlemode scenarios enables us yet again to uncover the nonclassical behavior of this bipartite state for all nonzero choices of parameters |α| and |β|. Note in this context that the phase of these parameters does not contribute because of the fully phase-randomized structure of the mixed state in Eq. (38).

Multimode superposition states
To further exceed the previous, bipartite state, we consider an N -mode state in this part. Specifically, we focus on a multimode superposition of coherent states [93], which consists of two N -fold tensor products of polar opposite coherent states, |±γ . Specifically, the skewsymmetric state |Ψ γ,N is of interest because it yields a GHZ state for |γ| → ∞ and W state for |γ| → 0, combining in an asymptotic manner two inequivalent forms of multipartite entanglement [94,14].
The Q functions for the states in Eq. (41) can be straightforwardly computed; they read To apply our criteria in Eq. (40), and for simplicity, we set α (m) j = α j for all mode numbers m and points in phase space, α j . In Fig. 5, we exemplify the certification of nonclassicality for the state |Ψ   41)] as a function of α 1 = Re(α 1 ) and for a fixed α 2 = 0. We remark that, for other mode numbers N , the plot looks quite similar. Most pronounced are nonclassical features for γ close to zero, relating to a W state in which a single photon is uniformly distributed over three modes. For large γ values, relating to a GHZ state, the negativities decrease, but det(M ) remains below zero. We reiterate that our relatively simple, second-order correlations of Q functions render it possible to certify the nonclassical properties of multimode, non-Gaussian states.

Generalized phase-space representations and nonlinear detection model
For demonstrating how our phase-space matrix approach functions beyond s-parametrized distributions, we consider an on-off detector that is based on two-photon absorption [95]. In this case, the POVM where :(ηn) 2n e −ηn /(2n)!: describes a measurement operator for 2n-photon states with a linear quantum efficiency η. In this context, it is worth mentioning that χ [eη 2 ]/[4n] has to be satisfied to ensure that the approximated POVM element correctly applies for photon numbers up to 2n [96]. The parameter χ relates to the nonlinear absorption efficiency.
In Fig. 6, we apply this approach and consider the single-mode even coherent state |Ψ (+) γ,1 [cf. Eq. (41) for N = 1], which is a non-Gaussian state, because we focused on the odd coherent state in the previous example. It is worth emphasizing that other methods to infer nonclassical light (e.g., the Chebyshev approach from Ref. [65]) are incapable to detect this state's quantum features. Here, we can directly certify nonclassicality of this non-Gaussian state despite the challenge of also having a non-Gaussian detection model.

Conclusion
In summary, we devised a generally applicable method that unifies nonclassicality criteria from correlation functions with quasiprobability distributions. Thereby, we created an advanced toolbox of nonclassicality tests which exploit the capabilities of both phase-space distributions and matrices of moments to probe for nonclassical effects. Furthermore, our framework is applicable to an arbitrary number of modes, arbitrary orders of correlation, and even phase-space functions perturbed through convolutions with non-Gaussian kernels. A measurement scheme was proposed to directly determine the elements of the phase-space matrix, the underlying key quantities of our method. In addition, we showed and discussed in detail that our treatment includes previous findings as special cases, is experimentally accessible even if other methods are not, and overcomes challenges of previous techniques when identifying nonclassicality.
The phase-space-matrix approach incorporates nonclassicality tests based on negativities of the phase-space distributions, including the Glauber-Sudarshan P function, and the matrix-of-moments approach as special cases. Thus, we were able to unify two major techniques for certification of nonclassicality. As the P function and the matrix of moments themselves are already necessary and sufficient conditions for the detection of nonclassicality, the introduced phase-space-matrix approach obeys the same universal feature. In other words, for any nonclassical state there exists a phase-space matrix condition which certifies its nonclassicality.
By applying our nonclassicality criteria to a diverse set of examples, we further demonstrated the power and versatility of our method. These examples covered discrete-and continuous-variable, singleand multimode, Gaussian and non-Gaussian, as well as pure and mixed quantum states of light. Remarkably, we used for all these states only the family of second-order correlations and phase-space distributions which are always nonnegative. Nevertheless, these basic criteria were already sufficient to certify distinct nonclassical effects on one common ground, further demonstrating the strength of our method. When compared to matrices of moments, the kinds of nonclassicality under study would require very different moments for determining the states' distinct quantum properties. Finally, we put forward an ex-perimental scheme, only relying on readily available optical components, to directly measure the quantities required to apply our method. This scheme applies even for imperfect detectors with a nonlinear response. Furthermore, we want to add that the practicality and strength of the matrix of phase-space distributions in certifying nonclassicality of lossy and noisy quantum state can be experimentally demonstrated [97].
Here, we focused on nonclassical effects of light, owing to their relevance for photonic quantum computation and optical quantum communication. The introduced approach may be further developed for the certification of other quantum features, such as non-Gaussianity. Currently, our method detects nonclassicality for Gaussian and non-Gaussian states equally, which could be further developed for a more finegrained quantumness analysis. However, instead of applying normal ordering, the construction of linear and nonlinear witnesses has to be adapted for this purpose. Furthermore, other kinds of quantum effects, such as entanglement, can be interpreted in terms of quasiprobabilities [22] and are similarly witnessed through correlations [98]. Thus, an extension to entanglement might be feasible as well. Therefore, our findings may provide the starting point for uncovering quantum characteristics through matrices of quasiprobabilities in other physical systems. Additionally, the derived framework can be utilized in the context of quantum information theory, such as in recently formulated resource theories of nonclassicality [6,7] and other measures of nonclassicality [99], which employ the phase-space formalism [22], thus potentially benefiting from our phase-space correlation conditions for future applications.
Note added. After finalizing this work, we have been made aware of a related work in preparation by J. Park, J. Lee, and H. Nha [100].