Verifying single-mode nonclassicality beyond negativity in phase space

While negativity in phase space is a well-known signature of nonclassicality, a wide variety of nonclassical states require their characterization beyond negativity. We establish a framework of nonclassicality in phase space that addresses nonclassical states comprehensively with a direct experimental evidence. This includes the negativity of phase-space distribution as a special case and further analyzes quantum states with positive distributions effectively. We prove that it detects all nonclassical Gaussian states and all non-Gaussian states of arbitrary dimension remarkably by examining three phase-space points only. Our formalism also provides an experimentally accessible lower bound for a nonclassicality measure based on trace distance. Importantly, this foundational approach can be further adapted to constitute practical tests in two directions looking into particle and wave nature of bosonic systems, via an array of nonideal on-off detectors and coarse-grained homodyne measurement, respectively. All these tests are practically powerful in characterizing nonclasssical states reliably against noise, making a versatile tool for a broad range of quantum systems in quantum technologies.


I. INTRODUCTION
Describing a quantum state of light or matter in phase space, e.g., Wigner function [1], is profoundly important to study quantum dynamics. It is a crucial tool to delineate the boundary between classical and quantum physics widely used in quantum optics [2], continuous variable (CV) quantum informatics [3,4] and other fields of quantum science [5]. As a classical phase-space distribution takes non-negative values like a probability distribution, a negativity emerging in quantum distribution is regarded as a signature of nonclassicality. However, negativity is just one aspect of multifaceted nonclassicality characterizing only a subset of nonclassical states. There exist quantum states with positive distributions that can nevertheless be classified as nonclassical, e.g., a squeezed state of light that is a key resource for CV quantum informatics [4] and single photons under a high-loss channel that are elementary information carriers for quantum informatics [6,7].
Nonclassical states are essential resources broadly for quantum informatics generating entangled states [8][9][10][11][12], providing advantage for quantum metrology [13][14][15] and quantum computation [16,17], etc.. A recent resource theory identified all quantum non-Gaussian states, even with positive Wigner functions, as a resource for quantum tasks, e.g., subchannel discrimination [18], which was further generalized to all CV nonclassical states [19]. It is thus crucial to establish a framework that can widely analyze nonclassical states beyond negativity. Specific properties were often used to characterize nonclassical states such as squeezing and photon-number statistics * Electronic address: jiyong.park@hanbat.ac.kr † Electronic address: hyunchul.nha@qatar.tamu.edu (sub-Poissonian) [20][21][22][23] extended also to multi-mode cases [24,25]. Distillation of nonclassicality can also be used to verify the nonclassicality of an initial state, however, requiring multiple copies of the same nonclasscial state and postselection [26][27][28]. More broadly, a quantum state tomography may be used to obtain complete information on a state thereby confirming nonclassicality [29]. However, it requires extensive measurements for sufficient data, and more seriously, an optimization process to find a physical state closest to obtained data. The data itself does not directly represent a legitimate quantum state rendering its significance weaker. It is necessary to characterize nonclassicality by examining phase-space in a faithful and resource-efficient way.
Adhering to negativity as a nonclassical feature, some works proposed to display negativity by modifying phasespace distributions, e.g., a regularized P -function under filtering process [30,31]. Other distributions closely related to the so-called s-parametrized functions [2] were also studied in view of photon statistics from on-off detectors [32,33]. Phase-space inequalities were also obtained by combining different s-parametrized functions useful to some extent [34]. Nevertheless, it is worth asking if the original Wigner function contains substantial information on nonclassicality beyond negativity. In this respect, Banaszek and Wódkiewicz proposed a Bell test examining four phase-space points to manifest nonlocality of two-mode states with positive Wigner functions [35], which were extended to generalized quasiprobability distributions [36] and genuine multipartite nonlocality [37][38][39]. The works in [40,41] demonstrated Bell-like tests also for single-mode nonclassicality and quantum non-Gaussianity. While conceptually remarkable and practically useful, these methods do not address a broad range of nonclassical states, e.g., squeezed states with purity < 0.86 are out of reach.
In this article, we propose a hierarchy of nonclassicality criteria in phase space that yields an efficient and broadly applicable test for CV systems. Our formalism addresses the Wigner function at n(n+1) 2 phase-space points progressively (n = 1, 2, . . . ). It includes the negativity of Wigner function at n = 1. Remarkably, it can detect all nonclassical Gaussian states and all non-Gaussian states of arbitrary dimension at the next level n = 2, i.e., looking into three phase-space points only. This opens a new possibility for a faithful and efficient test. We show that our foundational approach can constitute two practical tests characterizing nonclassical states reliably and efficiently from a particle and a wave point of views, respectively. It thus makes our method a versatile tool for a wide range of CV systems in quantum physics and technologies. We illustrate the practical power of our approach by examples. Our proposed approach is fruitful also in other aspects. It provides an experimentally accessible lower bound for nonclassical distance defined via trace norm [42], which is hard to obtain even theoretically. It can also be further extended to identify quantum non-Gaussianity [43][44][45][46][47][48][49][50][51][52][53][54][55][56][57] under energy constraint.

II. CRITERIA
Let us start with a general condition on classicality. A classical state, i.e., a mixture of coherent states, must satisfy for an arbitrary f (α) since its Sudarshan-Glauber-P function P ρc (α) is positive definite [58,59]. Our aim is to establish criteria that deal with the Wigner function at discrete phase-space points by choosing f (α) properly. Not only providing a fundamental insight, the Wignerfunction approach also leads to two general practical tests broadly applicable for CV systems, as shown later. To our aim, invoking the convolution between the Pfunction and the s-parametrized function [2] For the classicality to hold for arbitrary c i 's, we deduce the following theorem.

III. HIERARCHY
By its construction, M (n+1) 0 implies M (n) 0 since the matrix M (n+1) includes M (n) as its submatrix. That is, there naturally occurs a hierarchy of criteria with n increasing. If nonclassicality is confirmed at the level of n, it must be so at the next levels of n + 1, etc., but the converse is not always true.
Our formulation includes the negativity of Wigner function at the lowest n = 1, M (1) = π 2 W (β) 0. Then, it is fundamentally interesting, and practically important, to know how many phase-space points are required to verify noclassicality for states with positive Wigner functions. We prove below that our method can detect nonclassical states comprehensively using only three points {β 1 , β1+β2 2 , β 2 } on a line, i.e., M (2) 0 with

IV. GEOMETRIC INTERPRETATION
Before demonstrating its usefulness, let us briefly discuss the meaning of the classicality condition M (2) ≥ 0. One readily finds that all coherent states satisfy ρi , then makes a general classicality condition M (2) ≥ 0 for a mixture of coherent states. For a classical state, we thus see that the Wigner function at midpoint β1+β2 2 must be bounded by the geometric mean of the Wigner functions at end points β 1 and β 2 , importantly with a scaling factor e − 1 2 |β1−β2| 2 . In fact, this factor results from the commutator [â,â † ] = 1 representing the size of vacuum fluctuation.
Gaussian states. Every single-mode Gaussian state can be expressed as a displaced squeezed thermal state is a squeezing operator with strength r and angle φ of squeezing axis. ρ th (n) = ∞ n=0n n (n+1) n+1 |n n| is a thermal state with mean numbern. We can readily show M (2) 0 by taking three points along a squeezed axis [ Fig. 1(a)], with two end points at a distance 2d and the middle point at the origin. Our test turns out to be successful for a wide range of d as shown in Fig. 1(c). Without loss of generality, we consider an x-squeezed thermal state (γ, φ = 0), whose Wigner function is given by W σ (q, p) = 2µ π e −2e 2(r−rc ) q 2 e −2e −2(r+rc ) p 2 , with purity µ = (1 + 2n) −1 and critical squeezing r c = − 1 2 log µ. Section S1 of the Supplemental Material (SM) [60] gives its Wρ(d) 2 exp(−4d 2 ) for a Fock state |n under 80 % loss channel for n = 1 (red solid), n = 2 (gray dashed), n = 3 (black dot-dashed). R < 1 (shaded region) confirms nonclassicality for a wide range of displacement d.
Non-Gaussian states: More importantly, the threepoints test M (2) 0 can detect a broad range of non-Gaussian states. We first demonstrate its success for all non-Gaussian states of arbitrary truncation in Fock space. This includes as examples all noisy Fock states having positive Wigner functions. In Sec. S3 of SM [60], we further demonstrate that it can be extended to states of practical relevance having infinite Fock-state components.
The Wigner function of an arbitrary Fock-space truncated state (FSTS), ρ = N i,j=0 ρ jk |j k|, takes a form W ρ (α) = N i,j=0 ρ jk W |j k| (α), with W |j k| (α) given in Sec. S2 of SM [60]. As the case of negative Wigner functions is already treated at n = 1, we focus on the case of positive Wigner functions. Choosing β 1 = 2de iϕ and β 2 = 0 gives det ρ (de iϕ )e −4d 2 whose value less than 1 verifies nonclassicality. R(d) is a continuous function of d satisfying R(0) = 1. For the FSTS, we always find lim d→∞ R(d) = 0 with details in Sec. S2 of SM [60]. Therefore, there must be a finite d satisfying R(d) < 1 confirming nonclassicality. Remarkably, it works regardless of ϕ, i.e. insensitive to the axis of three points.
As an illustration, in Fig. 1, we plot the ratio R for (c) squeezed states and (d) Fock states under a 80%-loss channel. We confirm nonclassicality, R < 1, for a broad range of displacement d.

V. NONCLASSICALITY DISTANCE
It is also a topic of great interest to quantify the degree of nonclassicality for a given state ρ. A typical approach is to measure a distance between ρ and its closest classical state ρ c as N d (ρ) ≡ 1 2 min ρc∈C ||ρ − ρ c || 1 , with || · || 1 the trace norm and C the set of classical states. This is, however, very hard to obtain even if the state is completely known. Our formalism provides a lower bound for this nonclassical distance [42,61] enabling its practical estimation. With details in Sec. S4 of SM [60], we obtain where λ min is the least eigenvalue of M (n) at the level n. At n = 1, Eq. (7) shows that a negative value in phase space directly provides a reliable estimate for nonclassical distance. We can further estimate the nonclassical distance of a state with a postive Wigner function by using M (2) . For instance, for a general Gaussian state σ, using Eq. (6), which is beyond the results in Refs. [42,61] addressing only pure Gaussian states. We also establish connection between our approach and nonclassical depth [62] in S9 of SM [60].
With details in Sec. S5 of SM [60], if the least eigenvalue of M (2) ρ for a state ρ with energy E satisfies it confirms QNG. In Fig. 2, we plot ∆λ min = B(E) − λ min for a non- As seen from Fig. 2(a), our criterion detects QNG with f < 1 2 (positive Wigner function) for a squeezing r 0.237. Note that the squeezing operation does not create QNG as it is a Gaussian operation. In this context, the result also represents the QNG of f |2 2|+(1−f )|0 0| without squeezing. A recent ion-trap experiment realized a measurement in squeezed Fock basis, {Ŝ(r)|n : n = 0, 1, · · · } [69]. This can be adopted to verify QNG of stateŝ S(r)ρ nGŜ † (r) without performing squeezing on a non-Gaussian state ρ nG enhancing the range of QNG detection. Fig. 2(b) gives another example of a positive Wigner function with its QNG verified.

VII. PRACTICAL TESTS
The Wigner function corresponds to the number parity after displacement, i.e., W ρ (α) = 2 π tr[D † (α)ρD(α)(−1)n]. It is routinley measured in various systems, e.g., ion-trap [40,71] and circuit-QED [72], for which our proposed test M (2) can thus directly characterize nonclassicality. On the other hand, we can also derive alternative, practical, schemes out of Wigner-function framework, which can test nonclassicality reliably and efficiently against experimental imperfections. First, we present a generalized formalism to use on-off detectors registering photons without photon-number resolving (PNR). Second, we also (red solid) against n with binning size σ = 0.1. n: bin number for quadrature q = nσ. Grey shades represent the size of error due to finite data (a) ∼ 10 5 and (b) ∼ 10 6 .
present a marginal version of Wigner-function test, i.e., using M (q) = dpW ρ (q, p), which can be measured by homodyne detection well established for a wide variety of quantum systems including quantum optics [29], trapped ion [73], atomic ensemble [74], circuit cavity QED [75,76], and optomechanics [77,78]. Both of our proposed tests are powerful against noise with wide applicability.
On-off detector array. When an input light is equally divided via beam-splitting to impinge on N on-off detectors, the probability of k-detectors clicking is [79] p k [ρ] = tr ρ : with η detector efficiency and :Ô : normal-ordering.
The counting statistics p k above can also be obtained via the time-multiplexing approach using a single detector [33]. We first generalize our criterion to sparametrized function W ρ (α; s) [2] to use the counting Our classicality condition is readily generalized to M (s,n) 0 for an arbitrary s using elements in Eq. (11). We find the connection between the s-parametrized functions and the counting statistics p k in S6 of SM [60] as with each s m = 1 − 2N (N −m)η (m = 0, · · · , N − 1). Eq. (12) means that N different s-parametrized distributions W ρ (α; s m ) are determined by the counting statistics {p 0 , · · · , p N } obtained for a displaced stateD † (α)ρD(α). Furthermore, we prove in Sec. S2 of SM [60] that M (s,n=2) (three points test) can detect all nonclassical Gaussian and non-Gaussian states (FSTSs), importantly for an arbitrary s. This broader applicability beyond Wigner function makes our test robust against noise.
Let us illustrate the case of testing a Fock state |n under 50% loss channel mixed with a thermal photonn = 0.05 by using only N = 1 on-off detector of efficiency η = 0.7 [80]. We further consider an error due to finite data acquisition ∼ 10 5 (Sec. S8 of SM [60]). Our 3-points test Fig. 3(a), there exists a range of displacement to detect nonclassicality substantially beating the error. For instance, we have the signal to noise ratio as 1−R ∆R = 2.48 at d = 1.1. We also demonstrate the successful detection for other noisy Fock states with error analysis in Sec. S8 of SM [60].
Homodyne test. We next present a test using a marginal distribution M (q) = dpW ρ (q, p). Homodyne detection to measure M (q) is highly efficient, but requires a careful analysis. It is because that the actual homodyne data is coarse-grained due to finite binning, which may lead to a false detection of nonclassical effects [81][82][83]. Let σ be the binning size of homodyne data. Then all data in the range [−σ/2, σ/2] belong to the same bin yielding a coarse-grained distribution M σ ρ [n] ≡ σ/2 −σ/2 dδM ρ (nσ + δ) (n: bin number representing mean quadrature q = nσ). A classicality condition M (H) ≥ 0 then emerges with its elements where k can be either 0 or 1, with details in Sec. S7 of SM [60]. We prove in SM [60] that this marginal test even with a coarse-grained information detects all nonclassical Gaussian and non-Gaussian states (FSTSs). In Fig. 3(b), we show the result for a squeezed state (r = 0.3) un- 2∆ 2 e inφ ρe −inφ leaving no squeezing at ∆ = 1.2. Our homodyne test under coarse-graining (σ = 0.1) clearly manifests nonclassicality over 7 standard deviation, 1−R ∆R = 7.11. We also illustrate other cases in Sec. S8 of SM [60].

VIII. CONCLUSION
A phase-space approach usually provides us with a valuable insight into quantum physics [84]. While nega-tivity is one manifestation of nonclassicality, recent studies made it clear that all nonclassical states even without negativity are valuable resources for quantum information science [8][9][10][11]18]. It is thus critically important to establish a comprehensive framework of addressing nonclassical states with and without negativity covering a wide range of quantum systems. We have introduced a hierarchy of nonclassicality conditions that can address nonclassicality beyond negativity effectively and efficiently. Our approach makes it possible to analyze all nonclassical Gaussian states and non-Gaussian states using three phase-space points. Our formalism further provides a lower bound for nonclassical distance and a criterion to detect quantum non-Gaussianity with positive Wigner functions. Remarkably, our foundational approach also constitutes two practical tests looking into particle nature (number parity) and wave nature (marginal distribution), making a versatile tool for CV systems broadly. We illustrated the practical power of our tests adopting nonideal on-off detectors without resolving photon numbers and coarse-grained homodyne detection, respectively.
We hope our work could further stimulate works related to nonclassical effects from both a fundamental and a practical perspective. Our approach here clearly indicates that the information on nonclassicality is sufficiently imbedded in phase space even at a few points. Our geometric interpretation on classicality has stipulated the relation among the values of Wigner function, which is fundamentally associated with quantum fluctuation represented by a commutation relation or uncertainty principle. This seems worthwhile to further purse in studying nonclassicality for quantum multipartite systems as well. In the near term, we anticipate our framework can be useful for both theoretical and experimental analysis of quantum systems. In particular, our proposed tests can address all different CV systems including quantum optics, nano-or opto-mechanics, atomic ensemble, and circuit cavity QED, and so on.
Note added. We recently became aware of a closely related work by Bohmann, Agudelo and Sperling [85]. We note that our main idea and some results were earlier presented at an international conference [86].

S1. Optimal phase-space test for a Gaussian state
We first consider a 2 × 2 matrix A whose matrix elements are given by where F (x, y) = e −ax 2 −by 2 is a Gaussian function with a > c > b > 0. Then, we show that the minimum lowest eigenvalue of the matrix A is given by For a given function F (x, y) = e −ax 2 −by 2 , we can show that the points (x 1 , y 1 ) and (x 2 , y 2 ) minimizing the lowest eigenvalue of A must be (1) on x-axis, i.e., y 1 = y 2 = 0 and (2) symmetric with respect to the origin, i.e., x 1 + x 2 = 0. First, the lowest eigenvalue of A is given by For fixed values of A 11 = u and A 22 = v, the set of the points (x 1 , y 1 ) and (x 2 , y 2 ) satisfying F (x 1 , y 1 ) = u and F (x 2 , y 2 ) = v form two ellipses having the same center, directrix and major axis. Under the conditions a − c > 0 and b − c < 0, the ratio R is maximized at We now set y 1 = y 2 = 0 and rewrite Eq. (S3) as Now that λ in Eq. (S6) is a function of two variables A 11 = u and A 22 = v, we further fix its product A 11 A 22 = uv ≡ w. This also fixes the value of x 2 1 + x 2 2 = − 1 a log w. Note that X − √ X 2 + Y with Y > 0 decreases as X decreases or Y increases. Using the relation between the arithmetic mean and the geometric mean, we observe that A 11 + A 22 and A 11 A 22 (R − 1) are minimized and maximized, respectively, when |x 1 | = |x 2 | and u = v for a given uv = w.
All things considered together with |x 1 | = |x 2 | = x, we now need to optimize Examining its first derivative, we obtain the minimum lowest eigenvalue of A in Eq. (S2) at The Wigner function of the x-squeezed thermal state is given by W σ (α = q + ip) = 2µ π e −2e 2(r−rc ) q 2 e −2e −2(r+rc ) p 2 , with purity µ = (1 + 2n) −1 and critical squeezing r c = − 1 2 log µ. In view of Eq. (S1), its parameters read a = 2e 2(r−rc) , b = 2e −2(r+rc) and c = 2 satisfying the condition a > c > b > 0 for a nonclassical state r > r c . Therefore, its least eigenvalue of M  We here show that all nonclassical Gaussian states and non-Gaussian states of arbitrary Fock-space truncation can be detected via the matrix M (s,2) criterion for any s. This naturally includes the case of Wigner function test, s = 0.
Gaussian states-The s-parametrized quasiprobability function of a x-squeezed thermal state is given by where with r and r c representing the squeezing strength and the critical squeezing strength for nonclassicality, respectively. Using the minimum lowest eigenvalue of M (s,2) σ in Eq. (S2) of S1 with the parameters a and b for the Gaussian states, we see that it becomes negative for all nonclassical states r > r c regardless of s < 1.
Non-Gaussian states-The s-parametrized quasiprobability function of an arbitrary Fock-space truncated state (FSTS), ρ = N j=0 N k=0 ρ jk |j k|, is written as where W |j k| (α; s) for j ≥ k [1] is given by with a generalized Laguerre polynomial L We thus investigate a function of d 1−s (S16) whose value less than 1 confirms nonclassicality. R(d) is a continuous function of d satisfying R(0) = 1. For an FSTS, we always find lim d→∞ R(d) = 0 with a dominant contribution given by Therefore, there must be a finite d satisfying R(d) < 1, i.e. det M (2) < 0 for all s < 1.
Remarkably, the above proof works regardless of ϕ, i.e. insensitive to the axis of three points. We illustrate it by an example of non-rotationally symmetric state in phase space, that is, a superposition state 1 √ 2 (|0 + |1 ) under a 50 % loss. In Fig. S1, we show the results of testing it by choosing three points along x-axis (red solid), 1 √ 2 (x + p)axis (gray dashed) and p-axis (black dot-dashed), respectively. We can clearly see that the test is successful for a wide range of displacement d whatever axis is taken. Wρ(d) 2 exp(−4d 2 ) , where d represents the distance from the origin along x-axis (red solid), 1 √ 2 (x + p)-axis (gray dashed) and p-axis (black dot-dashed), respectively. R < 1 confirms the detection of nonclassicality for a wide range of d in each case.

S3. Non-Gaussian states with infinite Fock-state components
We here illustrate the detection of non-Gaussian states having infinite Fock-state components that are of current practical interest for CV quantum informatics.
Setting β 1 = 0 and β 2 = iπ 4γ , we observe that which manifests that our method can detect all nonclassical dephased cat states regardless of γ and f .

B. Photon added coherent states
Another infinite-dimensional quantum state of current interest is the photon-added coherent state, |Ψ ∼ a † |γ . As |Ψ is a pure non-Gaussian state with negativity, let us treat a realistic situation under a loss channel, i.e. ρ = Tr E U BS |Ψ Ψ| ⊗ |0 0| E U † BS , where U BS represents a beam splitting operation. We first note that a † |γ = a †D (γ)|0 =D(γ)(a † + γ * )|0 =D(γ)(|1 + γ * |0 ). That is, |Ψ is nothing but a displaced FSTS. As the local displacement followed by the beam splitting can be expressed in different order as the beam splitting followed by another displacement, we finally see that the mixed state ρ is just a lossy FSTS followed by a displacement. We have already proved that any FSTS can be detected under our formalism. As the final displacement can be incorporated to the choice of three phase-space points accordingly, it proves that all photon added states under a loss channel can be detected as well. The same idea is readily extended also to the multiple-times photon-added coherent states, |Ψ ∼ a †m |γ =D(γ)(a † + γ * ) m |0 , under a noisy channel.

S4. Nonclassicality distance
The nonclassical distance is defined as N d (ρ) ≡ 1 2 min ρc∈C ||ρ − ρ c || 1 , with || · || 1 the trace norm and C the set of classical states. If N d (ρ) = D, it allows a decomposition ρ = ρ c +D(ρ + −ρ − ) where ρ c is its nearest classical state under the trace measure. Dρ + and −Dρ − represent the mixtures of eigenstates for ρ − ρ c with positive and negative eigenvalues, respectively (trρ + = trρ − = 1). For a general operatorÔ, we have where Q represents the set of all quantum states. We setÔ = n i,j=1 u * i u jMij with its connection to M ij of Eq. (3) in main text as M ij = tr[ρM ij ]. Here u = {u 1 , u 2 , ..., u n } is specifically taken to be the eigenvector for the lowest eigenvalue of M. From tr[ρ cÔ ] ≥ 0, we obtain We have above used |tr[σÔ]| ≤ n for any state σ ∈ Q from |tr[σM ij ]| ≤ 1 and n i=1 |u i | ≤ √ n.

S5. Quantum non-Gaussianity test
Our formalism can further identify quantum non-Gaussianity (QNG), i.e. those states beyond a mixture of Gaussian states. Eq. (S9) is the lowest eigenvalue of M (2) for a Gaussian state, which can be rewritten in terms of energy E = tr[ρn]. We verify QNG if the least eigenvalue of a given state is smaller than those of Gaussian states with the same E We here derive the above Gaussian bound B(E) in the following steps. To begin with, note that the phasespace matrix M (n) is linear with respect to states, i.e. M In this case, the smallest eigenvalue λ of M (2) ρ cannot be less than the weighted sum of the smallest eigenvalues λ i of M (2) σi , i.e. λ ≥ i p i λ i . This can be readily seen by considering the eigenvectors corresponding to the least eigenvalues as |λ and |λ i , i.e. λ|M (2) |λ = λ and λ i |M (2) σi |λ i = λ i , respectively. Then, from M (2) = p i M (2) σi , we obtain λ = λ|M (2) |λ = p i λ|M (2) σi |λ ≥ p i λ i due to the condition λ|M (2) σi |λ ≥ λ i |M (2) σi |λ i = λ i .
(ii) A pure Gaussian state is a displaced squeezed state, σ =D(γ)Ŝ(r, φ)|0 0|Ŝ † (r, φ)D † (γ), which has the energy E = sinh 2 r + |α| 2 . On the other hand, as we explained in the main text, the displacement does not affect the least eigenvalue of M (2) σ , which is given by λ min,σ = −2e −r coth r sinh r in Eq. (S9). This means that the squeezed state without displacement is energyefficient to achieve the same level of least eigenvalue. We thus consider only the squeezed states with energy E = sinh 2 r, which gives the expression for its least eigenvalue.
(iii) Importantly, we note that λ E min,σ above is a convex function of E (Fig. 1). Thus, among all Gaussian states with the same average energy E, the pure squeezed state achieves the smallest eigenvalue. If a state is a mixture of pure squeezed states with E = i p i E i , its least eigenvalue cannot be less than that of a single pure squeezed states with energy E due to the convexity of the function in Eq. (S21). Therefore, λ E min,σ represents the smallest possible eigenvalue of the matrix M (2) σ among all Gaussian states σ with the energy E.
Combining (i), (ii) and (iii), we obtain a QNG criterion. That is, if the least eigenvalue of M (2) ρ for a state ρ with energy E is less than B(E) ≡ λ E min,σ , it confirms QNG. The state cannot be a mixture of Gaussian states.
In Fig. S2, we plot ∆λ min = B(E) − λ min for a non- For a fraction f > 1 2 , the state has the negativity of Wigner function, which is a trivial evidence of QNG. As can be seen from Fig. S2, our criterion can detect QNG with f < 1 2 for a squeezing r 0.237. Note that the squeezing operation does not change QNG as it is a Gaussian operation. In this respect, the result in Fig. S2 also represents the QNG of f |2 2|+(1− f )|0 0| without squeezing.
A recent ion-trap experiment realized a measurement in squeezed Fock basis, {Ŝ(r)|n : n = 0, 1, · · · } by controlling the interaction between the spin and the motional states [2]. This measurement can be adopted to verify QNG of those statesŜ(r)ρ nGŜ † (r) without performing squeezing on a non-Gaussian state ρ nG . Fig. S2 (c) gives another example of a positive Wigner function with its QNG verified. The QNG of a mixed cat state f |C C| + (1 − f )|0 0|, with a cat state |C ∼ |γ + | − γ , is confirmed up to the vacuum fraction 1 − f for each γ. Black solid represents the value of 1 − f at which the Wigner function becomes positive, above which our criterion detects QNG for the mixed cat state (red dashed). The yellow-filled region thus represents the detection of QNG for a non-Gaussian state with a positive Wigner function.

S6. Connection between the s-parametrized functions and the counting statistics from on-off detectors
We start with the counting statistics Eq. (10) in main text and a relation between the expectation value of a normally ordered operator for a quantum state ρ with its Glauber-P function, i.e., tr[ρ : f (n) :] = d 2 βP ρ (β)f (|β| 2 ) for an arbitrary well-defined function f [3]. We first obtain Using PD † (α)ρD(α) (β) = P ρ (α + β), we get Employing the convolution relation in Eq.
(2) of main text and d 2 βP ρ (β) = 1, we obtain where s m = 1 − 2N (N −m)η and T km is defined as We thus see that the counting statistics p k for a displaced state ρ α ≡D † (α)ρD(α) is composed of quasiprobability functions via a triangular matrix T . The relation can be represented as which leads to using the inverse matrix T −1 . A triangular matrix is invertible if and only if all elements on its principal diagonal are non-zero [4]. We find that T ii = 0 for all i, which means that we can always obtain N different s-parametrized quasiprobability distributions from the photocounting statistics via N on-off detectors using Eq. (S29).

S7. Power of criteria using marginal distributions
Marginal test: Choosing β i = q i + ip with the same p for all i in Eq. (3) of main text and integrating over Under a coarse-graining with binning size σ, we choose test points as q i = (2m i + k)σ + δ with integer m i , k = 0 or 1, and δ ∈ [− σ 2 , σ 2 ]. Integrating over δ, we have is the coarse-grained marginal distribution with n the bin number. Our whole procedures corresponds to the classicality condition as ij ≥ 0 thus constituting a matrix test with elements In the above, k = 0 (1) corresponds to the case of choosing all even (odd)-numbered bins, making two different tests. If one chooses even-numbered (n 1 = 2m 1 ) and odd-numbered (n 2 = 2m 2 + 1) bins together, the middle bin ( n1+n2 2 ) is not well-defined. Therefore, the matrix test in Eq. (S30) has been established according to the quadrature values q i = (2m i + k)σ + δ for a fixed k (0 or 1), while m i may vary, from the beginning.

Marginal distribution of Gaussian state
Without loss of generality, we again deal with the case of a x-squeezed thermal state. Its marginal distribution is given by with µ its purity, r the squeezing strength and r c = − 1 2 log µ the critical squeezing strength for nonclassicality. Considering the ratio R = which indicates that all nonclassical Gaussian states, r > r c , can be detected (R < 1) with an arbitrary non-zero d. More practically, we now examine a coarse-grained marginal distribution of Gaussian states.

Coarse-grained marginal distribution of Gaussian state
The coarse-grained marginal distribution of the Gaussian state is given by Here, erf(z) = 2 √ π z 0 e −t 2 dt is the error function that can be expanded as which indicates that erf(z) = e −z 2 z √ π {1+O(z −2 )} for z ≫ 1. Setting z 1 = √ 2e r−rc (n− 1 2 )σ and z 2 = √ 2e r−rc (n+ 1 2 )σ, we have which yields erf(z2) erf(z1) ≪ 1 for r > r c and |n| ≫ 1. If we investigate the ratio with z = √ 2e r−rc (2n − 1 2 )σ. For a sufficiently large n, we thus see R → 0 for r > r c . That is, R[n] < 1 in a certain range of large n. It indicates that all nonclassical Gaussian states can be witnessed by using coarse-grained marginal distribution with an arbitrary binning size σ.

Marginal distribution of FSTS
The marginal distribution of an arbitrary FSTS, ρ = N j=0 N k=0 ρ jk |j k|, is given by where the marginal distribution for the operator |j k| is given by Looking into the ratio R = , we obtain R ∝ d −2N for a large d ≫ 1. Thus, R → 0 for a very large d. This means that the nonclassicality of every FSTS can be verified, R < 1, in a certain range of large d. More practically, we now examine a coarse-grained marginal distribution of FSTS.

Coarse-grained marginal distribution of FSTS
The coarse-grained marginal distribution of FSTS is given by using M |j k| (q) in Eq. (S38). We are going to deal with M σ ρ [n] for a large n, so with Eq. (S40) in mind, we first look at the integration Here Γ(a, x) = ∞ x e −t t a−1 dt is an incomplete Gamma function [5], which can be expanded as which yields for n ≫ 1. Using Eqs. (S40), (S41) and (S45), we thus obtain for n ≫ 1. Considering the ratio , we thus see R ∝ (nσ) 1−2N for n ≫ 1, which indicates that the nonclassicality of every FSTS can be verified by using the coarse-grained marginal distribution with an arbitrary finite binning size σ.

S8. Error analysis and practical examples
The data acquisition procedure can be thought of as picking up outcomes randomly from a multinomial distribution. Let us assume that the total number of possible outcomes is k and the probability to obtain i-th outcome is p i with k i=1 p i = 1. If we pick up outcomes N s times and the number of observation for i-th outcome is X i , the variance of the probability estimator Xi Ns and the covariance between the probability estimators Xi Ns and Xj Ns are given by pi(1−pi) Ns and − pipj Ns , respectively [6]. For the test based on the on-off detectors, we obtain s-parametrized quasiprobability distributions by a linear transformation of the measured clickcounting probability as described in Sec. S5. For a linear combination Z = i a i X i , the variance of Z is given by (∆Z) 2 = i a 2 i (∆X i ) 2 + i i =j a i a j Cov(X i , X j ) where ∆x and Cov(x, y) represent the standard deviation of x and the covariance between x and y, respectively. Using the propagation of the uncertainty, we can estimate the statistical uncertainty of the obtained s-parametrized quasiprobability distibutions.
Furthermore, the statistical uncertainty of the ratio R = AB C 2 can be estimated as We can also estimate the statistical uncertainty of the ratio for the test based on the coarse-grained homodyne detection by a similar error analysis. We now illustrate the power of two practical tests using the on-off detectors and the homodyne detection, respectively, including the analysis of error due to a finite data acquisition.
On-off detector test-We first show the results for the case of Fock state |n under a 50% loss channel, which makes all Wigner functions positive definite. In Fig. S3, we demonstrate the detection of nonclassicality using N = 2 on-off detectors with a data number N s ∼ 10 5 . The red curve represents R = as a function of displacement d, with the grey shades representing the error size ∆R for each d. We see that there always exists a range of d for each state in which the value R is well below 1 with the error ∆R included, that is, 1 > R + ∆R. For instance, we have 1−R ∆R = 2.21 at d = 1.87 for the noisy Fock state |3 , so the signal beats the error over 2 standard deviation. Interestingly, the signal to noise ratio increases with n, as seen from figures, as 1 Adopting this strategy, we analyze the case of lower detector efficiency η = 0.75 in Fig. S4. We again confirm the successful detection of nonclassicality reliably against noise with a lower detector efficiency for a broad range of displacements.
Homodyne test-We now demonstrate the detection of nonclassicality using a coarse-grained marginal distribution with a binning size σ = 0.1 and a data number N s ∼ 10 6 . In Fig. S5, we show the results for the same noisy Fock states as in Fig. S3. Black circles repre- at each value of m with m 1 = 10, which corresponds to the choice of three quadratures q 1 = 2m 1 σ, q 2 = (m 1 + m)σ, and q 3 = 2mσ for test. The grey shades represent the error size ∆R due to finite data N s ∼ 10 6 . We again see that there always exists a range of m for each state to achieve 1 > R+ ∆R. For instance, we have 1−R ∆R = 21.34 at m = 4 for the noisy Fock state |3 , so the signal beats the error over 21 standard deviation. For other Fock states |4 , |5 , and |6 , the signal to noise ratio is 1−R ∆R =17.73, 11.84, and 6.48 at m=5, 6, and 7, respectively.
As seen from the above examples, the corase-grained homodyne test is very robust against experimental imperfections. We further demonstrate its practical power m designates a bin number corresponding to quadrature q = mσ. Grey shades represent the size of error due to finite data ∼ 10 6 . R < 1 confirms nonclassicality. The results clearly manifest nonclassicality beating error in the orange region, where the error bar is negligible so not appreciable in the plot.
with other examples. Fig. S6 (a) confirms the case of a noisy single-photon state f |1 1|+(1−f )|0 0| at f = 0.3, and Fig. S6 (b) for a single photon state under a 50% loss channel mixed with a thermal photonn = 0.1. For comparison, there have been some protocols proposed to distill squeezing by which one may confirm nonclassicality upon the postselected (distilled) copies. For example, R. Filip remarkably proposed a distillation of squeezing for a noisy single-photon state [8]. Its single-copy distilla- n designates a bin number corresponding to quadrature q = nσ. Grey shades represent the size of error due to finite data ∼ 10 6 . R < 1 confirms nonclassicality. The results clearly manifest nonclassicality beating error in the orange region. tion ( Fig. 1 of [8]) entails a low probability of success, e.g. < 10 −2 for the case of f = 0.3, with distilled squeezing low to detect. In constrast, our method does not require postselection of data and manifests nonclassicality robust against practical errors. Next we move on to consider another practical noise, i.e. phase-diffusion of an optical signal. A phase diffusion can be described as where ∆ is the standard deviation of random phase-shift. The phase diffusion affects the marginal distribution of a quantum state ρ as where M ρ (x θ ) represents the marginal distribution of the quantum state ρ for the quadraturex θ =â e −iθ +â † e iθ 2 .
As shown in Fig. S7 (a), the phase-diffusion destroys squeezing at ∆ > 1.074, ∆ > 0.901 and ∆ > 0.696 for squeezed states with squeezing parameter r = 0.1, r = 0.2 and r = 0.4, respectively. The squeezing can be distilled after phase-diffusion, as shown in Fig. S7 (b), e.g. by applying the protocol in [9] that generally requires multi copies of the same nonclassical states and postselection. Our homodyne test directly confirms nonclassicality for these phase-diffused squeezed states. Fig.  S8 shows a strong result from our test for the case of ∆ = 1.2, at which all considered states completely lose squeezing, against coarse-graining of homodyne data and finite data acqusition.
Comparison with moment test-The above examples illustrate that our homodyne test verifies nonclassicality when the usual squeezing test by homodyne detection fails. We may further compare our test with moment-based homodyne tests. Suppose that one obtains a marginal distribution M (q) = dpW ρ (q, p) from a homodyne measurement, where W ρ (q, p) is the Wigner function. The distribution M (q) can be related to the marginal distributionM (q) = dpP ρ (q, p) of the Sudarshan-Glauber P -function P ρ (q, p). That is, M (q) = 2 π dqM (q)e −2(q−q) 2 . If the state is classical, the dis-tributionM (q) must be positive-definite. This can be checked in terms of moments q k ≡ dqM (q)q k . We can readily see that an n × n matrix Q, whose elements are given by Q ij = q i+j , must be positive definite. For example, the violation of the condition Q ≥ 0 at n = 2 represents the usual quadrature squeezing.
We have performed the analysis on the phase-diffused squeezed states using Q matrix test. This test fails to detect nonclassicality at n = 2 (squeezing) and n = 3 for a severely decohered state, so we move to the level n = 4. Fig. S9 shows the results for the squeezed state with r = 0.4 and r = 0.5 against the phase diffusion ∆. For a fair comparison, we include the effects of finite data ∼ 10 6 and the coarse-graining σ = 0.1 for the error analysis. To confirm nonclassicality, the least eigenvalue (red solid) of the matrix Q must be negative. Consdering the error level (black dashed), we see that the test fails at a large phase-diffusion, ∆ > 1.76 (r = 0.4 case) and ∆ > 1 (r = 0.5 case). One may try to enhance the test by further going up to the higher level of n, which however is not much favoarable in the presence of errors due to finite data and coarse-graining. In contrast, as we show in Fig. S10, our homodyne test clearly manifests nonclassicality even at the large phase-diffusion ∆ = 2 for the same squeezed states.

S9. Estimating the nonclassical depth
Our formulation also provides an estimate for another nonclassicality measure, i.e. nonclassical depth [10], as follows.
If the nonclassical depth of a quantum state ρ is τ , the