Introduction—Continuous-variable (CV) quantum systems, particularly optical platforms, constitute a central architecture for modern quantum technologies due to their experimental accessibility, high-bandwidth operation, and intrinsic scalability [1]–[3]. In such systems, quantum states are naturally described by phase-space quasiprobability distributions [4]–[7]. Among them, the Wigner function plays a distinguished role, providing a complete phase-space representation of quantum states while retaining many features of a classical probability distribution [8]–[10].

The significance of the Wigner function can be further highlighted by its connection to the Gottesman-Knill theorem, which establishes that quantum circuits restricted to Clifford operations admit efficient classical simulation [11]. This efficient simulability, however, generally breaks down in the presence of Wigner negativity and is therefore widely regarded as a necessary resource for quantum computational advantage [12]. More broadly, states having negative Wigner-functions demonstrate quantum advantages in a variety of information tasks, including quantum computational speedup [13], quantum error correction [14], quantum state distillation [15]–[17], and many more [12], [18].

Driven by the fundamental role of Wigner negativity in quantum information processing, significant effort has been devoted to its detection and quantification [19]–[30]. Existing methods largely rely on quantum-state tomography and subsequent reconstruction of the Wigner function [31]–[34]. However, these procedures require a number of state copies that grows rapidly with system size [35], [36], rendering the detection of Wigner negativity increasingly demanding for large-scale and arbitrary quantum states. This motivates the search for methods that can detect and quantify Wigner negativity only from few copies, without full phase-space reconstruction.

In this Letter, we address this question affirmatively by developing an efficient framework for the detection and quantification of Wigner negativity. Our approach is based on moments of the Wigner function, a family of global phase-space quantities obtained by integrating powers of the distribution. We show that positivity of the Wigner function imposes nontrivial constraints on these moments and exploit this structure to derive complementary hierarchies of negativity detection criteria. Beyond detection, we obtain certifiable lower bounds on logarithmic Wigner negativity measure, directly from a finite set of moments and identify an optimal experimentally accessible quantifier requiring only the second and fourth moments.

Moreover, we establish an exact multicopy representation of arbitrary Wigner moments in terms of parity observables, thereby providing a direct operational interpretation and a scalable route to their experimental estimation. Combined with classical-shadow protocols, this enables detection and quantification of Wigner negativity directly from randomized measurements using only a small number of state copies, without prior knowledge of the state or full phase-space tomography. Our framework can be directly applied for entanglement detection.

Preliminaries—We consider an \(M\)-mode CV quantum system described by bosonic annihilation operators \(\vec{a}\mathrel{\vcenter{:}}=(a_1,\dots,a_M)\) satisfying the canonical commutation relations \([a_m,a_{m'}]=[a_m^\dagger,a_{m'}^\dagger]=0\) and \([a_m,a_{m'}^\dagger]=\delta_{m,m'}\mathbb{I}\). The associated Hilbert space is \(\mathcal{H}=\bigotimes_{m=1}^M \mathcal{H}_m\), where each \(\mathcal{H}_m\) is an infinite-dimensional space spanned by the Fock basis \(\{\ket{n}\}_{n\in\mathbb{N}}\), and physical states are represented by density operators \(\rho \in \mathcal{D}(\mathcal{H})\). Introducing the quadrature operators \(\hat{x}_m = ({a_m + a_m^\dagger})/{\sqrt{2}}, ~\hat{p}_m = ({a_m - a_m^\dagger})/{\sqrt{2}i}\) which satisfy \([\hat{x}_m,\hat{p}_{m'}]=i\delta_{m,m'}\mathbb{I}\) (throughout the letter, we consider \(\hbar=1\)), we obtain a \(2M\)-dimensional phase-space description with canonical coordinates \((\vec{x},\vec{p})\). In this representation, the quantum state \(\rho\) can be described by the Wigner function [37], [38] \[W_{\rho}(\vec{x},\vec{p}) = \frac{1}{(2\pi)^M} \int d^M y \; \bra{\vec{x}+\frac{\vec{y}}{2}} \rho \ket{\vec{x}-\frac{\vec{y}}{2}} \, e^{-i\vec{p}\cdot\vec{y}},\] which is real and normalized, i.e., \(\int d\vec{x}\, d\vec{p} \; W(\vec{x},\vec{p}) = 1\). Unlike classical probability distributions, however, the Wigner function can attain negative values, a hallmark of nonclassicality [9], that constitutes a key resource for CV quantum information processing. The relevant and essential features pertaining to the Wigner function are summarized in Supplemental Material (SM) A.

Moment-based detection of Wigner negativity—To detect Wigner negativity of arbitrary quantum states, we first develop a moment-based framework that provides an efficient and experimentally accessible alternative to full quantum state tomography (Basic features of application of moments in entanglement detection are presented in Supplemental Material (SM) B). To this end, we define the \(m\)th-order moment of the Wigner function (\(w_m\)) as follows:

The moments \(\{w_m\}\) encode global properties of the phase-space distribution. Importantly, when \(W\) is everywhere nonnegative, the sequence \(\{w_m\}\) is constrained by a hierarchy of inequalities inherited from positive measures. Any violation of these constraints therefore certifies Wigner negativity. Below, we derive three complementary families of such constraints.

1. Detection hierarchy using \(\mathcal{L}_p\) norms (LP). Our first hierarchy is rooted in fundamental \(\mathcal{L}_p\)-norm inequalities obeyed by positive phase-space distributions. Expressed in terms of Wigner moments, these relations yield the following theorem.

Any violation of the above inequality implies Wigner negativity. The proof follows from the properties of \(\mathcal{L}_p\) norms [39], [40] \(\|f\|_p = \left( \int |f|^p \, d\mu \right)^{1/p}\), where \(d\mu = W \, dx^N dp^N\), which is a probability measure whenever \(W(\vec{x}, \vec{p}) \ge 0\). The proof is completed using Hölder’s inequality in \(\mathcal{L}_p\) spaces. Full details are provided in SM C.0◻

While Theorem 1 establishes a systematic set of constraints relating consecutive moments, each inequality involves only two moments at a time. This motivates the search for stronger consistency conditions involving multiple moments simultaneously.

2. Detection hierarchy using log-convexity (LC) condition. Beyond the \(\mathcal{L}_p\)-norm hierarchy, positivity of the Wigner function also imposes log-convexity constraints on the moment sequence. This yields the following condition.

3. Moment-matrix approach via Hankel positivity. All the hierarchy proposed so far do not necessarily provide strictly stronger criteria for detecting Wigner negativity, as they are based on a restricted subset of moments rather than the complete moment sequence and hence fail to fully capture the information contained in the complete moment sequence. A stronger approach is to consider all moments up to a given order simultaneously through Hankel matrices (moment matrices), which play a central role in classical moment problems (see SM B for a brief discussion). This leads to a hierarchy of moment-matrix based criteria that strictly strengthen the previously obtained moment inequalities.

For the proof, we refer to SM D. Theorem 3 establishes that positivity of the Wigner function imposes an infinite hierarchy of constraints in the form of positive semidefinite Hankel matrices. These matrices are naturally nested, as each higher-order matrix incorporates all lower-order moments. This nested structure leads to a monotonic hierarchy of constraints, where higher-order conditions are strictly stronger than the lower-order ones. The following corollary formalizes this property.

Thus the Hankel-based criteria form a strictly increasing hierarchy of constraints, where higher-order matrices provide progressively stronger conditions for detecting Wigner negativity, i.e., \[H_1 \preceq H_2 \preceq H_3 \, \cdots.\]

The lowest-order members of the three hierarchies coincide, yielding the common criterion \(\Gamma_{3} =w_2^2-w_3>0\) for the detection of Wigner negativity. Notably, this criterion is consistent with the previously established third-moment witness of Ref. [41].

To assess the performance of the proposed criteria, we analyze mixed Fock (Eq. 6 ) and noisy cat states (Eq. 7 ) in Table 1 and Fig. 1, respectively. We compare the Wigner-negative regions detected by the LP, LC, and Hankel hierarchies for the above states respectively. In both cases, the Hankel criteria exhibit a systematic improvement with increasing order and detect substantially larger portions of the nonclassical region, eventually recovering almost the full negativity range. A general comparison of the strength of different criteria and their detection performance is provided in the End Matter (see SM E for detailed calculation).

Operational quantifier of Wigner negativity—While the criteria developed above provide practical and efficient way to identify Wigner negativity, they do not directly quantify its magnitude. This naturally motivates the problem of quantification: can one estimate the amount of Wigner negativity using only a finite number of experimentally accessible moments?

A standard quantifier of phase-space nonclassicality is the logarithmic Wigner negativity [42], \[L(W) =\log \,\!\left( \int |W(\vec{x},\vec{p})| \, d^N x\, d^N p \right), \label{eq:logWignerNeg}\tag{1}\] which measures the total volume of negativity contained in the Wigner distribution. By construction, \(L(W)=0\) whenever the Wigner function is nonnegative and \(L(W)>0\) whenever Wigner negativity is present. However, evaluating Eq. 1 generally requires knowledge of the full Wigner function, making direct experimental estimation challenging. To connect logarithmic Wigner negativity with the moment framework developed above, we introduce the absolute-valued Wigner moments \[\tilde{w}_r =\int |W(\vec{x},\vec{p})|^r \, d^N x\, d^N p, \qquad r\in\mathbb{N}. \label{eq:absolute95moments}\tag{2}\] These quantities contain information about the distribution of negativity across phase space and satisfy \(L(W)=\log(\tilde{w}_1)\). Now we define the family of moment-based quantities \[L_r(W) =\frac{1}{2-r} \Big[ \log(\tilde{w}_r) + (1-r)\log(\tilde{w}_2) \Big], \quad r>2. \label{eq:Lr}\tag{3}\] An important feature of this construction is that for every even integer \(r=2n\), \(\tilde{w}_{2n}=w_{2n}\). Consequently, the quantities \(L_{2n}(W)\) can be evaluated directly from the experimentally accessible Wigner moments introduced previously. Moreover, these quantities can be directly obtained from the same measurement data obtained while detecting Wigner negativity using Wigner moments. The significance of \(L_r(W)\) is established by the following theorem.

The proof is provided in SM F. The theorem shows that low-order Wigner moments alone are sufficient to construct experimentally accessible lower bounds to the logarithmic Wigner negativity. The quantities \(L_r(W)\) also retain a direct witnessing interpretation.

Thus every positive value of \(L_r(W)\) simultaneously certifies Wigner negativity and provides a nontrivial lower bound on the total negativity quantified by \(L(W)\). Finally, among all experimentally accessible members of the hierarchy, the lowest-order quantifier yields the strongest certified bound.

Therefore \(L_4(W)\) provides the tightest lower bound to logarithmic Wigner negativity obtainable from experimentally accessible even-order moments. The proofs of Corollary 2 and Proposition 1 along with further properties of \(L_r\) can be found in SM F.

Experimental feasibility—The Wigner moments \(\{ w_n \}\), introduced in Definition 1 are nonlinear functionals of the Wigner function and therefore appear, at first sight, to require full phase-space reconstruction. However, we show that this is not the case. Remarkably, every Wigner moment admits an exact representation as the expectation value of a universal observable acting on a finite number of copies of the quantum state. Starting from the phase-space representation of the Wigner function, we introduce the multi-copy operator \[A_n =2\int d^2\alpha\, \Delta(\alpha)^{\otimes n},\] where \(\Delta(\alpha)\) is the displaced parity operator and \(\alpha = (x+ip)/\sqrt{2}\), such that \(dx \,dp \equiv 2 d^2\alpha\). Consequently, the \(n\)-th Wigner moment can be written as \[w_n =\mathrm{Tr} \!\left[ \rho^{\otimes n}A_n \right]. \label{eq:wn95multicopy}\tag{4}\] The operators \(\{A_n\}\) are Hermitian by construction. The exact structure of \(A_n\) may be obtained by analyzing its position-space kernel and exploiting a decomposition into collective and relative degrees of freedom (see SM G). The resulting operator takes a remarkably simple and universal form, which we state in the following theorem.

The proof can be found in SM G. Eq. (?? ) reveals a simple geometric structure underlying all Wigner moments. In the collective-coordinate basis, one mode remains unchanged while all \((n-1)\) relative modes undergo a parity reflection. Consequently, the entire hierarchy of Wigner moments is encoded in the response of an \(n\)-copy state to parity operations acting only on the relative degrees of freedom. The complexity of the phase-space integral is thus replaced by a remarkably simple multi-copy observable. Beyond its conceptual significance, Theorem 5 immediately yields an operational measurement scheme. Instead of reconstructing the Wigner function, one prepares \(n\) identical copies of the state, applies the interferometric transformation \(F_n\), and measures the joint parity of the relative modes. Explicitly, \[w_n =\frac{1}{n\pi^{\,n-1}} \left\langle I\otimes \Pi^{\otimes(n-1)} \right\rangle_{F_n\rho^{\otimes n}F_n^\dagger}. \label{eq:operational95formula}\tag{5}\] Thus every Wigner moment acquires a direct operational interpretation as a measurable multi-copy parity expectation value. This establishes the fundamental link between phase-space nonclassicality and experimentally accessible observables that underlies the detection and quantification protocols developed in this Letter. Furthermore, since the moments \(w_n\) are expressed as multicopy nonlinear functions of the density operator, they can be estimated efficiently using classical-shadow techniques [43]. Numerical simulations for truncated Fock states in Fig 2 support that both the detection witness \(\Gamma_3\) and the quantifier \(L_4\) can be efficiently recovered from randomized measurements using only a modest number of classical shadows, thereby substantially reducing the number of state copies required compared to full quantum state tomography (see End Matter for details of the numerical experiment).

Applications—The moment hierarchies developed in this Letter are not limited to witnessing Wigner negativity. More generally, whenever negativity of a suitably constructed Wigner function certifies a quantum resource [12], [17], [44], the same framework yields efficient resource witnesses directly in terms of Wigner moments. As a representative example, we consider detection of quantum entanglement. Recent studies have established that negativity of suitably defined reduced or collective Wigner functions certifies bipartite NPT entanglement [45] or even, genuine multipartite entanglement [46], [47]. Consequently, the moment hierarchies developed here can be directly elevated from witnesses of Wigner negativity to experimentally accessible entanglement witnesses. In particular, violation of any of the criteria developed in this letter for the corresponding reduced Wigner function certify entanglement of the underlying state. Further details and explicit examples are presented in SM H.

Conclusions—In summary, we have developed a unified moment-based framework for the detection and quantification of Wigner negativity in arbitrary CV states. By exploiting fundamental constraints obeyed by positive phase-space distributions, we have established complementary hierarchies of negativity detection criteria based on \(\mathcal{L}_p\)-norm inequalities, log-convexity relations, and Hankel-matrix positivity. While the \(LP\) and \(LC\) hierarchies yield witnesses involving only a few consecutive moments, the Hankel hierarchy provides a monotonic sequence of increasingly powerful tests, substantially enhancing the detection of Wigner negativity, particularly in mixed states. The three hierarchies naturally support a progressive detection strategy: one may first test the \(LP\) and \(LC\) conditions, and invoke the more powerful Hankel criteria only when the first two witnesses are inconclusive. The additional moment information estimated for the Hankel hierarchy is rewarded by a corresponding increase in detection power, yielding a flexible trade-off between experimental cost and detection capability.

Beyond detection, we have obtained a family of certifiable lower bounds to the logarithmic Wigner negativity, with the optimal experimentally accessible bound given by \(L_4\), depending only on the second and fourth Wigner moments. Further, by establishing an exact multi-copy representation of arbitrary Wigner moments in terms of parity measurements on collective and relative degrees of freedom, we have provided an operational route to their experimental estimation. Combined with randomized-measurement and classical-shadow protocols, our framework provides a practical and scalable approach to certifying and quantifying nonclassical resources using only a modest number of state copies without prior state information or full phase-space reconstruction.

Finally, our results establish Wigner moments as a unifying link between phase-space nonclassicality, efficient quantitative resource characterization, and experimentally accessible observables. Beyond the results presented here, an intriguing open direction is to derive the optimal and complete set of moment conditions for certifying Wigner negativity of arbitrary quantum states and consequently to study whether the complete hierarchy of Wigner moments can be used to determine the geometry of the set of Wigner-positive states.

Acknowledgements.– BM acknowledges DST INSPIRE fellowship program for financial support.

End Matter↩︎

Appendix I: Wigner function for Fock and Cat states↩︎

For completeness, we collect here the explicit Wigner functions of the representative states used as examples throughout the main text.

1. Fock states. The Fock states of a harmonic oscillator are defined as \(|m\rangle \equiv \frac{(a^\dagger)^m}{\sqrt{m!}} |0\rangle\) with Wigner function, \[W_{|m\rangle}(x,p) = \frac{(-1)^m}{\pi} e^{-(x^2+p^2)} \mathcal{L}_m\!\left(2(x^2+p^2)\right), \nonumber\] where \(\mathcal{L}_m(z)\) denotes the \(m\)-th Laguerre polynomial. While the vacuum state \(|0\rangle\) is Gaussian and therefore Wigner positive, all states with \(m\ge1\) are Wigner negative. Consistent with this fact, we find that already the lowest-order Hankel condition detects nonclassicality of the Fock states, with \(\det[H_1] < 0\) for all \(m \ge 1\). Now we consider a mixed state, \[\label{fock95mixed} \rho = \lambda |0\rangle\langle 0| + (1-\lambda) |1\rangle\langle 1| ,\tag{6}\] with \(0 \le \lambda \le 1\). The associated Wigner function is \[W(x,p) = \frac{1}{2\pi} \Big[(1-\lambda)(x^2+p^2) + 2\lambda - 1\Big] e^{-\frac{1}{2}(x^2+p^2)} . \nonumber\] The state is Wigner negative for \(0<\lambda<0.5\). We compare the different moment-based criteria by determining the range of \(\lambda\) over which negativity is detected; the results are summarized in Table 1.

2. Cat states. Schrödinger cat states, defined as the superpositions of two coherent states, \(\ket{\mathrm{cat}_{\pm} (\alpha)} = \frac{1}{\mathcal{N}} \left( \ket{\alpha} \pm \ket{-\alpha} \right)\), with \(\mathcal{N} = \sqrt{2(1 \pm e^{-2\alpha^2})}\), possess the Wigner function [48] \[\begin{align} W_{\text{cat}_{\pm}} (x,p) = \frac{1}{2 \pi \mathcal{N}^2} &\Bigg( e^{ -\frac{1}{2} \left\{ (x-\alpha)^2 +p^2 \right\}} +e^{ -\frac{1}{2} \left\{ (x+\alpha)^2 +p^2 \right\}} \nonumber\\ & \pm \cos (\sqrt{2} \alpha x ) e^{ -\frac{1}{2} \left( x^2 +p^2 \right)} \Bigg) \nonumber \end{align}\] whose interference term gives rise to Wigner negativity, a clear signature of nonclassicality. In contrast, the classical mixture of the same coherent states, \(\rho_{\mathrm{cl}} = \frac{1}{2} \left( \ket{\alpha}\bra{\alpha} + \ket{-\alpha}\bra{-\alpha} \right)\), although non-Gaussian, possesses the positive Wigner function \[\begin{align} W_{\mathrm{cl}}(x,p) = \frac{1}{4 \pi} \Bigg( e^{ -\frac{1}{2} \left\{ (x-\alpha)^2 +p^2 \right\}} +e^{ -\frac{1}{2} \left\{ (x+\alpha)^2 +p^2 \right\}} \Bigg). \nonumber \end{align}\] We consider the noisy cat state, \[\label{cat95mixed} \rho_{\mathrm{mix}} = (1-\lambda)\,\ket{\mathrm{cat}_{\pm}}\bra{\mathrm{cat}_{\pm}} + \lambda\, \rho_{\mathrm{cl}},\tag{7}\] with \(0\leq\lambda\le1\) and corresponding Wigner function \[W_{\mathrm{mix}}(x,p) = (1-\lambda)\, W_{\mathrm{cat}}(x,p) + \lambda\, W_{\mathrm{cl}}(x,p). \nonumber\] We compare the moment-based criteria developed here by determining the range of \(\lambda\) for which Wigner negativity is detected for \(2\leq \alpha\leq 5\) in Fig. 1.

Appendix II: Details of numerical simulations: Shadow estimation of Wigner moments↩︎

The numerical implementation adapts the framework of Ref. [43] to CV phase-space quantities, where the fundamental objects are moments of the Wigner distribution. The simulations are performed in a truncated Fock basis containing the first \(d\) Fock states. For each copy of the quantum state, a unitary matrix \(U\in U(d)\) is sampled independently from the Haar measure. The use of Haar-random unitaries guarantees that every measurement basis is drawn uniformly from the unitary group, thereby implementing the standard global-measurement version of the classical-shadow protocol. After applying \(U\), a projective measurement is performed in the Fock basis, producing an outcome \(b\) with probability \[p(b|U) =\langle b|U\rho U^\dagger|b\rangle .\]

Each measurement outcome is converted into a classical shadow through the inverse measurement channel \[\hat{\rho} =(d+1)\, U^\dagger |b\rangle\!\langle b|U -I,\] which satisfies \(\mathbb{E}[\hat{\rho}] =\rho.\) Consequently, every shadow provides an unbiased estimator of the underlying quantum state. Repeating the procedure over \(N\) independent copies generates a collection of independent shadow snapshots \(\{\hat{\rho}_1,\hat{\rho}_2,\ldots,\hat{\rho}_N\}\). To connect the shadow data with phase-space quantities, we precompute the Wigner kernel operators \(\Delta(\alpha)\) on a discretized phase-space grid. For each shadow and each phase-space point, we evaluate \[\hat{W_j}(\alpha) =\mathrm{Tr}\!\left[\hat{\rho}_j\,\Delta(\alpha)\right].\] Because the shadow estimator is unbiased, one has \(\mathbb{E}\!\left[\hat{W}_j(\alpha)\right] = W(\alpha)\), so that the collection of shadow-derived quasiprobabilities provides direct access to phase-space moments. Estimating higher-order moments requires products of several independent copies of the Wigner function. Rather than multiplying shadow estimates obtained from the same measurement record, which would introduce finite-sample bias, we employ unbiased U-statistic estimators. For a moment of order \(k\), the estimator is constructed from all distinct \(k\)-tuples of shadows.

In Fig. 2, for each state considered in the numerical analysis, we generate \(N=500\) shadows and repeat the entire estimation procedure over \(100\) independent Monte-Carlo realizations. The reported values correspond to sample means, while the shaded regions represent \(95\%\) confidence intervals obtained from the observed fluctuations among realizations.

The numerical demonstrations presented in the main text focus on Fock states. This choice is made solely for clarity and benchmarking purposes. However, the method itself is more general. Any single-mode CV state admits a Fock-basis representation \(\rho =\sum_{m,n=0}^{\infty} \rho_{mn} |m\rangle\!\langle n|\), and finite-energy states can be approximated arbitrarily accurately by truncation to a sufficiently large Fock subspace. Once such a truncation is chosen, the classical-shadow protocol described above applies without modification. The estimation procedure is therefore not tied to Fock states and can be implemented for any arbitrary states. The extension to multimode systems is equally natural. For an \(M\)-mode state, the phase-space kernel becomes \(\Delta(\boldsymbol{\alpha}) =\bigotimes_{j=1}^{M} \Delta(\alpha_j)\), and the corresponding moments are defined through integrals over the full \(2M\)-dimensional phase space. Since classical shadows naturally estimate operators acting on tensor-product Hilbert spaces, the same measurement and reconstruction protocol immediately yields unbiased estimators for multimode Wigner moments. Consequently, the moment hierarchy admits a direct shadow-based implementation beyond the single-mode setting.

For a CV state truncated to a \(d\)-dimensional Fock subspace, quantum-state tomography requires \(O(d^2/\epsilon^2)\) copies for generic mixed states, where \(\epsilon\) denotes the desired estimation accuracy [35], [36]. By contrast, classical shadows estimate observables directly, with copy complexity \(O(\|A\|_{\rm shadow}^2/\epsilon^2)\) for an observable \(A\) [43]. Consequently, estimating a small number of Wigner moments may require substantially fewer copies than reconstructing the entire density matrix, while the same randomized measurements can be reused to estimate multiple moments simultaneously. The numerical results presented in the main text (see Fig. 2) show that sufficiently accurate estimation can already be achieved with a modest number of copies, indicating the practical feasibility of the proposed moment-based framework.

Appendix III: Performance study of different moment hierarchies↩︎

Here the main objective is to understand how efficiently the different criteria detect Wigner negativity when the maximum accessible moment order is fixed. To perform this analysis, we generate random quantum states in the truncated Fock space \(\mathcal{H}_N = \mathrm{span}\left\{|0\rangle,|1\rangle,\dots,|N\rangle \right\}\), with cutoff values ranging from \(N=1\) to \(N=10\). The random states are sampled according to the Hilbert-Schmidt ensemble. Explicitly, each density matrix is generated as \[\rho = \frac{GG^\dagger}{\text{Tr}(GG^\dagger)},\] where \(G\) is a complex Ginibre matrix with independently distributed Gaussian entries. By construction, the resulting states are positive semidefinite and properly normalized. For each random state, the Wigner function is computed numerically using the standard Fock-basis representation implemented in QuTiP. The Wigner moments \(\{ w_n\}\) are then evaluated numerically through direct integration over phase space. A state is classified as Wigner negative whenever the minimum value of the Wigner function attains negative values on the numerical phase-space grid, below certain threshold \(( -10^{-15})\). The fraction of randomly generated Wigner-negative states detected by each criterion is summarized in Fig. 3. Several important features emerge from the numerical analysis.

First, Wigner negativity rapidly becomes generic as the Fock cutoff increases. In particular, for \(N\ge4\), essentially all randomly generated states exhibit Wigner negativity. This behavior reflects the fact that positivity of the Wigner function represents a constrained class within the full CV state space.

Second, the numerical results clearly demonstrate the hierarchy of detection strengths (see SM E), \[\text{Hankel} \Longrightarrow \text{LC} \Longrightarrow \text{LP}.\] For every cutoff value, the Hankel hierarchy consistently detects the largest fraction of Wigner-negative states, followed by the LC hierarchy and finally the LP hierarchy.

Third, increasing the order of the scalar LP conditions does not improve the detection capability. In fact, the higher-order LP conditions rapidly become ineffective, with LP\((7)\) failing to detect any states in the present numerical analysis. This behavior indicates that isolated higher-order scalar inequalities capture only limited information about the underlying moment structure. A similar behavior is observed for the LC hierarchy. While LC\((4)\) substantially improves over LC\((2)\), the higher-order condition LC\((6)\) eventually loses detection power as the cutoff increases.

In contrast, the Hankel hierarchy becomes significantly stronger with increasing order. In particular, the higher-order Hankel conditions \(H_2\) and \(H_3\) detect almost all Wigner-negative states for sufficiently large cutoff values. This demonstrates that incorporating the global structure of the moment sequence through positivity of the full moment matrix is significantly more powerful than considering isolated scalar inequalities.

 

1 A. Quantum information in phase space: Wigner function↩︎

Here, we briefly review the phase-space formulation of quantum mechanics based on Wigner function, which provides a powerful description of continuous-variable (CV) quantum systems. Such systems describe bosonic degrees of freedom, such as modes of the electromagnetic field. An \(M\)-mode system is described by annihilation operators \(\vec{a} \mathrel{\vcenter{:}}= (a_1,\dots,a_M)\) satisfying the canonical commutation relations \[[a_m,a_{m'}]=[a_m^\dagger,a_{m'}^\dagger]=0, \qquad [a_m,a_{m'}^\dagger]=\delta_{m,m'}\mathbb{I}.\]

The associated Hilbert space is \(\mathcal{H}=\bigotimes_{m=1}^M \mathcal{H}_m\), where each \(\mathcal{H}_m\) is an infinite-dimensional space spanned by the Fock basis \(\{\ket{n}\}_{n\in\mathbb{N}}\), and physical states are represented by density operators \(\rho \in \mathcal{D}(\mathcal{H})\). A convenient representation of CV systems is obtained by introducing quadrature operators \[\hat{x}_m = \frac{a_m + a_m^\dagger}{\sqrt{2}}, \qquad \hat{p}_m = \frac{a_m - a_m^\dagger}{\sqrt{2}i},\] which satisfy \([\hat{x}_m,\hat{p}_{m'}]=i\delta_{m,m'}\mathbb{I}\) (we set \(\hbar=1\)). These operators define a \(2M\)-dimensional phase space with canonical coordinates \((\vec{x},\vec{p})\). Equivalently, one may use the complex variables \(\alpha_m = (x_m + i p_m)/\sqrt{2}\), which provide a compact description of phase space commonly used in quantum optics.

The phase-space formulation of quantum mechanics represents quantum states through quasiprobability distributions. Among these, the Wigner function plays a central role. For a state \(\rho\), it can be defined in the \((\vec{x},\vec{p})\) representation as \[W_{\rho}(\vec{x},\vec{p}) = \frac{1}{(2\pi)^M} \int d^M y \; \bra{\vec{x}+\frac{\vec{y}}{2}} \rho \ket{\vec{x}-\frac{\vec{y}}{2}} \, e^{-i\vec{p}\cdot\vec{y}},\] or, equivalently, in the complex phase-space picture via the displaced parity operator, \[W_\rho(\vec{\alpha}) = \mathrm{Tr}\!\left[\rho \, \Delta (\vec{\alpha})\right],\] where \(\Delta(\vec{\alpha}) = \left(\frac{1}{\pi}\right)^MD(\vec{\alpha}) \Pi D^\dagger(\vec{\alpha})\), with \(\Pi = \prod_m e^{-i\pi a_m^\dagger a_m}\) is the parity operator and \(D(\vec{\alpha}) = \prod_m e^{\alpha_m a_m^\dagger - \alpha_m^* a_m}\) is the displacement operator.

These two representations are completely equivalent and provide complementary perspectives: the \((x,p)\) formulation emphasizes the analogy with classical phase space, while the \(\alpha\)-representation highlights the role of displacement operations and is closely connected to experimental implementations.

The Wigner function is real and normalized, \[\int dx^M \, dp^M \; W(\vec{x},\vec{p}) = 1.\] Further, the Wigner function of a single mode system corresponding to a density matrix \(\hat{\rho}\) can be expressed as \[W(x,p)= \frac{1}{2 \pi} \int\limits_{ -\infty}^{\infty} \bra{x+\frac{y}{2}}\hat{\rho} \ket{x-\frac{y}{2}}e^{-ipy}dy.\] For a pure state \(\ket{\psi}\) with wave function \(\psi(x)\), this Wigner function simplifies to \[W(x,p)=\frac{1}{2\pi} \int\limits_{-\infty}^{\infty} \psi(x+\frac{y}{2})\psi^{*}(x-\frac{y}{2}) e^{-ipy} dy,\] and its marginals reproduce measurable probability distributions, \[\begin{align} \int W(x,p)\,dp = |\psi(x)|^2, \quad \int W(x,p)\,dx = |\phi(p)|^2, \end{align}\] where \(\psi(x)\) and \(\phi(p)\) denote the position and momentum wavefunctions, respectively. Despite the classical-like features, \(W(\vec{x},\vec{p})\) can take negative values, reflecting the nonclassical nature of quantum states, such negativities cannot be reproduced by any classical joint probability distribution and are therefore regarded as a signature of nonclassicality.

Another key aspect of the phase-space approach is the correspondence between operators and phase-space functions via the Weyl transform. For an operator \(\hat{A}\), one defines \[\tilde{A}(x,p) = \int dy \; e^{-ipy} \bra{x+\frac{y}{2}} \hat{A} \ket{x-\frac{y}{2}},\] such that expectation values take the classical-like form \[\langle \hat{A} \rangle = \int dx\,dp \; W(x,p)\,\tilde{A}(x,p).\] This mapping provides a direct bridge between operator algebra in Hilbert space and functions on phase space.

Within this framework, a fundamental distinction arises between Gaussian and non-Gaussian states. Gaussian states are characterized by Gaussian Wigner functions and therefore possess non-negative phase-space representations. Hudson’s theorem states that a pure state has a non-negative Wigner function if and only if it is Gaussian [49], [50]. This equivalence, however, does not extend to mixed states, for which non-Gaussian states may still exhibit positive Wigner functions [51], [52]. Consequently, Wigner-negative states constitute a distinguished subclass of non-Gaussian states, and Wigner negativity is widely regarded as a key indicator of nonclassicality and a valuable resource for continuous-variable quantum information processing.

2 B. The classical moment problem and its applications in quantum information theory↩︎

The classical moment problem concerns the characterization of sequences that arise as moments of a positive measure. Given a sequence of real numbers \(\{m_k\}_{k=0}^\infty\), the problem is to determine whether there exists a positive measure \(\mu\) such that \[m_k = \int x^k \, d\mu(x), \quad \forall k \in \mathbb{N}.\] Sequences admitting such a representation are referred to as moment sequences. A central result in this context is that positivity of the underlying measure imposes strong consistency conditions on the sequence \(\{m_k\}\). In particular, the associated Hankel matrices constructed from the moments must be positive semidefinite at all orders. More explicitly, defining the Hankel matrix \(H_n(\mathbf{m})\) with elements \[[H_n(\mathbf{m})]_{ij} = m_{i+j}, \quad i,j \in \{0,1,\dots,n\},\] a necessary and sufficient condition for \(\{m_k\}\) to correspond to a positive measure \(\mu\) is \[H_n(\mathbf{m}) \geq 0 \quad \forall n.\] Equivalently, this condition can be understood as the requirement that all quadratic forms generated by polynomials are nonnegative, \[\int \left( \sum_{i=0}^k c_i x^i \right)^2 d\mu(x) \ge 0,\] for any real coefficients \(\{c_i\}\). These constraints form a hierarchy of increasingly stringent conditions that any valid moment sequence must satisfy.

The framework of moment problems has found significant applications in quantum information theory, particularly in the characterization of quantum correlations and nonclassical features. One prominent example is the use of moments of the partially transposed density matrix, referred to as PT moments, for the detection of quantum entanglement [53]–[60]. Given a bipartite quantum state \(\rho_{AB}\), the PT moments are defined as \[p_n = \text{Tr}\left[(\rho_{AB}^{T_B})^n\right], \label{PTmoments}\tag{8}\] for integers \(n \geq 1\). These quantities provide indirect access to the spectrum of the partially transposed state \(\rho^{T_B}_{AB}\), whose eigenvalues \(\{\lambda_i\}\) determine entanglement properties. In particular, the characteristic polynomial \[\text{Det}(\rho^{T_B}_{AB} - \lambda I) = \sum_k a_k \lambda^k,\] has coefficients \(a_k\) that are functions of the moments \(\{p_n\}\). Since the transposition map is not physically implementable, these moments offer an experimentally accessible route to probe the spectrum.

A key advantage of moment-based approaches is their experimental accessibility. Techniques such as shadow tomography [43], [61], [62] allow efficient estimation of moments without requiring full state reconstruction. Since only a polylogarithmic number of copies is needed, these methods are particularly well-suited for high-dimensional and many-body systems. Moment-based methods have been further generalized to probe other forms of nonclassicality, including Kirkwood–Dirac nonpositivity [63], non-Markovianity [64] and quantum imaginarity [65].

3 C. Proof of Theorem 1↩︎

Here we aim to prove that for any positive normalized Wigner function \(W\) and every integer \(n\geq 2\), \((w_n)^{n-2} - (w_{n-1})^{n-1} \ge 0\). Since \(W \ge 0\) and \(\int W \, dx^N dp^N= 1\), we define a probability measure \[d\mu = W \, dx^N dp^N, \;\; \text{such that}\;\; \int d\mu =1.\] With respect to \(\mu\), the moments can be written as: \[w_n = \int W^{n-1} \, d\mu, \quad w_{n-1} = \int W^{n-2} \, d\mu.\] We now consider \(\mathcal{L}_p\) norms with respect to this probability measure, defined as \[\|f\|_p = \left( \int |f|^p \, d\mu \right)^{1/p}.\]

Hölder’s inequality in \(\mathcal{L}_p\) space: Let \(1 \leq p < \infty\) and \(1 < q < \infty\) be such that \(\frac{1}{p} + \frac{1}{q} = 1\). For functions \(f \in \mathcal{L}_p\) and \(g \in \mathcal{L}_q\), their product \(fg\) belongs to \(\mathcal{L}_1\). Furthermore, Hölder’s inequality states that \[\|fg\|_{1} \leq \|f\|_{p} \, \|g\|_{q}. \label{eq:32holder}\tag{9}\] As a special case, when \(p = q = 2\), this inequality reduces to the well-known Cauchy-Schwarz inequality. Here we choose, \[f = W^{n-2}, \quad g = 1,\] and exponents \[p = \frac{n-1}{n-2}, \quad q = n-1.\] Then, using Eq. 9 we get, \[\|W^{n-2}\|_1 \le \|W^{n-2}\|_p \|1\|_q.\] We now evaluate each term; \[\begin{align} \|W^{n-2}\|_1 &= \int |W^{n-2}| \, d\mu \nonumber \\ & = \int W^{n-2} \, d\mu \qquad \text{(since W\geq0)} \nonumber\\ &= w_{n-1},\\ \|W^{n-2}\|_p &= \left( \int |W^{n-2}|^p \, d\mu \right)^{1/p} \nonumber \\ &= \left( \int W^{(n-2) \left( \frac{n-1}{n-2} \right)} \, d\mu \right)^{1/p} \nonumber \\ &=\left( \int W^{n-1} \, d\mu \right)^{1/p} \nonumber \\ &=w_n^{1/p}. \end{align}\] Finally, \[\|1\|_q = \left( \int d\mu \right)^{1/q} = 1,\] since \(\mu\) is a probability measure. Therefore, \[w_{n-1} \le w_n^{1/p}.\] Substituting \(p = \frac{n-1}{n-2}\), we obtain \[w_{n-1} \le w_n^{\frac{n-2}{n-1}}.\] Raising both sides to the power \((n-1)\) gives \[\begin{align} &(w_{n-1})^{n-1} \le (w_n)^{n-2} \nonumber\\ \Rightarrow \;\;&(w_n)^{n-2} - (w_{n-1})^{n-1} \geq 0, \end{align}\] which completes the proof.

4 D. Proof of Theorem 3↩︎

We now aim to prove that for any positive normalized Wigner function \(W\), the Hankel matrix \(H_n(\boldsymbol{w})\) with elements, \([H_n(\boldsymbol{w})]_{ij} = w_{i+j+1}\), \(i,j \in \{0,1,\dots,n\}\) is positive semidefinite for all \(n \in \mathbb{N}\).

Let \(\boldsymbol{y} = (y_0, y_1, \dots, y_n) \in \mathbb{R}^{n+1}\) be an arbitrary real vector. Consider the quadratic form associated with the Hankel matrix \(H_n(\boldsymbol{w})\): \[\begin{align} \boldsymbol{y}^T H_n(\boldsymbol{w}) \boldsymbol{y} &= \sum_{i=0}^m \sum_{j=0}^m y_i y_j \, w_{i+j+1}. \end{align}\] Substituting the definition of the Wigner moments from Eq. ?? , we obtain \[\begin{align} \boldsymbol{y}^T H_n(\boldsymbol{w}) \boldsymbol{y} &= \sum_{i=0}^n \sum_{j=0}^n y_i y_j \int_{\mathbb{R}^{2N}} W^{\,i+j+1} \, d^N x \, d^N p . \end{align}\]

Since \(W \ge 0\) and the summations are finite, Fubini’s theorem allows us to interchange the order of summation and integration, yielding \[\begin{align} \boldsymbol{y}^T H_n(\boldsymbol{w}) \boldsymbol{y} &= \int_{\mathbb{R}^{2N}} \left( \sum_{i=0}^n \sum_{j=0}^n y_i y_j W^{\,i+j+1} \right) d^N x\, d^N p . \end{align}\]

We now rewrite the integrand by factoring out a single power of \(W\): \[\begin{align} \sum_{i,j=0}^n y_i y_j W^{i+j+1} &= W \sum_{i,j=0}^N y_i y_j W^i W^j . \end{align}\] The double sum appearing above can be recognized as a perfect square, \[\begin{align} \sum_{i,j=0}^n y_i y_j W^i W^j = \left( \sum_{i=0}^n y_i W^i \right)^2 . \end{align}\] Therefore, the quadratic form can be written as, \[\begin{align} \boldsymbol{y}^T H_n(\boldsymbol{w}) \boldsymbol{y} &= \int_{\mathbb{R}^{2N}} W \left( \sum_{i=0}^n y_i W^i \right)^2 d^N x\, d^N p . \end{align}\]

By assumption, \(W \ge 0\) in all phase points, and the squared term is manifestly nonnegative. Hence, the integrand is nonnegative pointwise over phase space, implying \[\boldsymbol{y}^T H_n(\boldsymbol{w}) \boldsymbol{y} \ge 0 \qquad \forall \;\boldsymbol{y} \in \mathbb{R}^{n+1}.\] This establishes that the Hankel matrix \(H_n(\boldsymbol{w})\) is positive semidefinite, completing the proof.

5 E. Comparison between the LP, LC and Hankel hierarchies↩︎

In the main text, we have introduced three different hierarchies of moment-based conditions for detecting Wigner negativity, namely: (i) the LP hierarchy, (ii) the LC hierarchy, and (iii) the Hankel-matrix hierarchy. While all these conditions follow from positivity of the Wigner function, they capture different structural aspects of the moment sequence. It is natural to ask how they are related to one another and which hierarchy provides the strongest detection power. In this section, we establish the precise analytical hierarchy among these conditions.

This implies that, Hankel matrix condition is stronger than the LC condition when the maximum available moment order is fixed. Consequently, violation of any such LC condition, guaranties the violation of the Hankel condition. However the opposite is not true, which is supported by the examples studied in this letter. Next, we compare the LC hierarchy with the LP hierarchy.

Combining Propositions 2, 3, we obtain the hierarchy for a fixed maximum accessible moment order, \[\text{Hankel} \Longrightarrow \text{LC} \Longrightarrow \text{LP}.\]

To justify the relative strength of the hierarchies, we emphasize that the LP conditions involve only two consecutive moments and therefore probe only limited information about the full moment sequence. In contrast, the LC conditions involve correlations among three neighboring moments and are therefore strictly stronger. Similarly, the LC conditions arise only from positivity of the \(2\times2\) principal minors of the Hankel matrices. Positivity of the full Hankel matrices imposes significantly stronger constraints, since it incorporates all moments up to a given order simultaneously. Consequently, there exist moment sequences satisfying all LC conditions while still violating positivity of the corresponding Hankel matrix.

Therefore, when the maximum accessible moment order is fixed, the Hankel hierarchy always provides the strongest detection criterion, followed by the LC hierarchy and finally the LP hierarchy. This explains the behavior observed in the examples studied in the main text, where the Hankel conditions consistently detect the largest set of Wigner-negative states.

6 F. Quantifying Wigner negativity using moments↩︎

A natural way to quantify Wigner negativity is through the absolute volume of the Wigner function. For a \(N\)-mode state with Wigner function \(W\), the negativity measure is given by \[\mathcal{N}(W) = \int |W(\vec{x},\vec{p})| \, d^N x \, d^N p .\] For all quantum states with a positive Wigner function, one has \(\mathcal{N}(W) = 1\), since \(W\) is normalized. In contrast, whenever \(W\) attains negative values, the positive and negative contributions no longer cancel under the absolute value and thus \(\mathcal{N}(W) > 1\). Hence, \(\mathcal{N}(W)\) serves as a faithful indicator of Wigner negativity. In analogy with logarithmic negativity in entanglement theory [66], the logarithmic measure of Wigner negativity can be defined as [42] \[L(W) = \log \left( \int |W| \, d^N x \, d^N p \right) = \log(\tilde{w}_1),\] By construction, \(L(W)=0\) for all Wigner-positive states and increases monotonically with the amount of negativity. However, evaluating this quantity requires pointwise knowledge of \(|W|\), making it experimentally demanding as it generally necessitates full state tomography. This motivates the development of alternative quantifiers that can be estimated efficiently from limited moment data, without requiring full reconstruction of the state.

To connect these quantifiers with the moment-based framework developed in this letter, it is natural to consider moments of the non-negative function \(|W|\). We therefore define the absolute-valued Wigner moments as \[\tilde{w}_m = \int |W(\vec{x},\vec{p})|^m \, d^N x \, d^N p \;, \qquad m \in \mathbb{N}.\]

Since \(|W|\) is non-negative, but generally not normalized whenever the Wigner function contains negative regions, it is natural to associate with it the normalized probability distribution \[f(\vec{x},\vec{p}) = \frac{|W(\vec{x},\vec{p})|^2}{\tilde{w}_2},\] which satisfies \(f\ge 0\) and \(\int f =1\).

We define the Rényi entropy of order \(r\) associated with this distribution as \[{R}_r(W) = \frac{1}{1-r} \log \left( \int f(\vec{x},\vec{p})^r \, d^N x \, d^N p \right).\] Using the definition of absolute-valued Wigner moments, this can be written as \[{R}_r(W) = \frac{1}{1-r} \log \left( \frac{\tilde{w}_{2r}}{\tilde{w}_2^r} \right).\] We now define the following quantity:

\[\begin{align} L_r(W) &= \frac{1}{2-r} \left[ \log(\tilde{w}_r) + (1 - r)\log(\tilde{w}_2) \right] \\ &= \frac{1}{2} \left( {R}_{r/2}(W) - \log(\tilde{w}_2) \right). \end{align}\]

We note that \(\tilde{w}_r = w_r\) whenever \(r\) is an even integer. This also clarifies the choice of \(\tilde{w}_2\) in the definition of \(f(\vec{x},\vec{p})\), as it is the lowest-order nontrivial absolute Wigner moment that can be directly estimated within the framework developed in this letter. As a result, \(L_r(W)\) can be efficiently evaluated for any even \(r\). Importantly, the required moments are extracted from the same measurement data used for negativity detection, enabling simultaneous detection and quantification without additional experimental overhead.

Before proving the results stated in the main text, we first establish a standard property of the Rényi entropies associated with the normalized distribution \(f(\vec{x},\vec{p})\). This observation will play a central role in identifying the optimal member of the moment-based quantifier hierarchy.

Proof of Theorem 4↩︎

Using the definitions of \(L_r(W)\) and \(L(W)\), \[L_r(W)-L(W) = \frac{1}{2-r} \log\!\left( \frac{\tilde{w}_r\,\tilde{w}_1^{\,r-2}}{\tilde{w}_2^{\,r-1}} \right).\] Since \(r>2\), we have \(2-r<0\). Therefore it suffices to prove \[\tilde{w}_2^{,r-1} \le \tilde{w}_1^{,r-2}\tilde{w}_r.\]

Applying the interpolation inequality for \(L_p\) norms, \[|W|_2 \le |W|_1^{\frac{r-2}{2(r-1)}} |W|_r^{\frac{r}{2(r-1)}},\] and raising both sides to the power \(2(r-1)\) yields \[|W|_2^{\,2(r-1)} \le |W|_1^{\,r-2} |W|_r^{\,r}.\] Using \[|W|_1=\tilde{w}_1,\qquad |W|_2^2=\tilde{w}_2,\qquad |W|_r^r=\tilde{w}_r,\] we obtain \[\tilde{w}_2^{\,r-1} \le \tilde{w}_1^{\,r-2}\tilde{w}_r,\] which proves \[L_r(W)\le L(W).\]

Proof of Corollary 2↩︎

If \(W(\vec{x},\vec{p})\ge0\), then \(\tilde{w}_1 = 1\) and therefore \(L(W)=0\). Then Theorem 4 immediately gives \[L_r(W)\le L(W)=0.\]

Conversely, if \(L_r(W)>0\), then Theorem 4 implies \(L(W)>0\), which can only occur when the Wigner function attains negative values. This completes the proof.

Proof of Proposition 1↩︎

From \(L_r(W) = \frac{1}{2} \left( R_{r/2}(W)-\log\tilde{w}_2 \right)\), together with Lemma 1, it follows immediately that \[L_r(W)\le L_s(W), \qquad r>s,\] because the term \(\tilde{w}_2\) is independent of \(r\). Hence the hierarchy \(\{L_r(W)\}_{r>2}\) is monotonically decreasing. Restricting to experimentally accessible even orders yields \[L_4(W)\ge L_6(W)\ge L_8(W)\ge \cdots.\] Since every \(L_r(W)\) provides a certifiable lower bound to the logarithmic Wigner negativity, the largest member of this hierarchy furnishes the tightest bound. Hence, \(L_4(W)\) is the optimal experimentally accessible moment-based quantifier. 0◻

Other properties of \(L_r(W)\)↩︎

We now discuss some other fundamental properties of the quantity \(L_r(W)\), which further clarify its role as a quantitative measure of Wigner negativity.

1. Invariance under Gaussian unitaries. For any symplectic transformation \(S\) in phase space, \[L_r(W \circ S) = L_r(W).\]

   Proof. Under a Gaussian unitary, the Wigner function transforms according to \[W(\vec{x},\vec{p}) \mapsto W(S^{-1}(\vec{x},\vec{p})),\] where \(S\) is a symplectic transformation preserving phase-space volume. Since the Lebesgue measure is invariant under such transformations, it follows that \[\tilde{w}_k = \int |W(\vec{x},\vec{p})|^k \, d^N x \, d^N p\] remains unchanged for all \(k\). As \(L_r(W)\) depends only on \(\tilde{w}_r\) and \(\tilde{w}_2\), it is therefore invariant under Gaussian unitaries. ◻

2. Additivity. For product states, \[L_r(W_1 \otimes W_2) = L_r(W_1) + L_r(W_2).\]

   Proof. For a product state, the Wigner function factorizes as \[W_{12}(\vec{x}_1,\vec{p}_1,\vec{x}_2,\vec{p}_2) = W_1(\vec{x}_1,\vec{p}_1)\, W_2(\vec{x}_2,\vec{p}_2).\] Consequently, \[\tilde{w}_k(W_1 \otimes W_2) = \int |W_1 W_2|^k = \tilde{w}_k(W_1)\,\tilde{w}_k(W_2),\] for all \(k\). Substituting this multiplicative relation into the definition of \(L_r(W)\) and using the additivity of the logarithm yields \[L_r(W_1 \otimes W_2) = L_r(W_1) + L_r(W_2).\] ◻

Illustrative examples↩︎

To demonstrate the effectiveness of the proposed moment-based quantifiers, we evaluate them for two representative classes of non-Gaussian states, as illustrated in Fig. 4. The first example we consider here is the \(N00N\) state, which is a two-mode quantum state where all \(N\) photons can be found in either one of the two modes \(a\) and \(b\), while the other remains vacant. It is a maximally path-entangled number state of the form \[|\psi, N\rangle = \frac{1}{\sqrt{2}}\left(|N\rangle_{a} |0\rangle_{b} + e^{i\phi}|0\rangle_{a} |N\rangle_{b}\right).\]

For the choice \(\phi = \pi\), the corresponding Wigner function in terms of the dimensionless quadratures \((x_1,p_1)\) and \((x_2,p_2)\) is given by [2] \[\begin{align} W_{N}(x_1,p_1,x_2,p_2) = \frac{1}{2\pi^{2}N!} e^{-(x_1^2+x_2^2+p_1^2+p_2^2)} &\times \Big[ -2^{N} \big( (x_1 + ip_1)^{N} (x_2 - ip_{2})^{N} + (x_1 - ip_1)^{N} (x_2 + ip_{2})^{N} \big) \\ &+(-1)^{N} N!\big( \mathcal{L}_{N}[2(x_1^2+p_1^2)] + \mathcal{L}_{N}[2(x_2^2+p_2^2)]\big)\Big], \end{align}\] where \(\mathcal{L}_N(x)\) denotes the Laguerre polynomial. This state exhibits Wigner negativity for all \(N \geq 1\).

The second example we consider here is the photon added coherent state, which is obtained by repeated application of the photon creation operator on a coherent state, which initially has properties like a classical field. The state is defined by

\[\begin{align} \ket{\alpha, m} = \frac{ {a^{\dagger}}^m \ket{\alpha} }{ (\langle\alpha|a^m {a^\dagger}^m | \alpha\rangle)^{1/2}}. \end{align}\] Here \(\ket{\alpha}\) is a coherent state, with \(\alpha \in\mathbb{C}\) and m is an integer. The Wigner function for this state is given by [67], \[\begin{align} W(z) = \frac{2 (-1)^m \mathcal{L}_m(|2z-\alpha|^2)) }{\pi \mathcal{L}_m(-|\alpha|^2)} \exp( -2|z-\alpha|^2 ), \;\; z\in \mathbb{C} \end{align}\]

For all of these families of states, the logarithmic Wigner negativity \(L(W)\) increases monotonically with the excitation number, indicating the enhancement of nonclassical features (see Fig. 4). The quantities \({L}_4(W)\) and \({L}_6(W)\) consistently act as lower bounds to \({L}_{N}(W)\), in accordance with our theoretical framework. Furthermore, the ordering \({L}_4(W) \geq {L}_6(W)\) is observed across all cases, with lower-order moments providing tighter bounds. These observations confirm that the proposed moment-based measures successfully capture both qualitative trends and quantitative aspects of Wigner negativity. At the same time, there exist Wigner-negative states for which \(L_r(W)\) does not yield a strictly positive value. This is because the quantities \(L_r(W)\) are constructed from only a finite set of low-order Wigner moments and therefore capture only partial information about the underlying phase-space distribution. The resulting loss of sensitivity is the natural trade-off for obtaining experimentally accessible and efficiently computable lower bounds on Wigner negativity.

7 G. Operator Representation and Accessibility of Wigner Moments↩︎

In this section, we derive an exact operator representation of the Wigner moments. We show that the \(n\)-th Wigner moment can be expressed as the expectation value of a universal multi-copy observable acting on \(n\) identical copies of the quantum state. This representation naturally reveals an experimentally feasible measurement scheme based on multi-mode interferometry and parity detection.

Multi-copy representation of Wigner moments↩︎

The Wigner function of a single-mode quantum state \(\rho\) can be expressed as \[W(\alpha) =\mathrm{Tr}\!\left[\rho\,\Delta(\alpha)\right], \label{eq:wigner95operator95representation}\tag{12}\] where \[\Delta(\alpha) =\frac{1}{\pi} D(\alpha)\Pi D^{\dagger}(\alpha). \label{eq:phase95point95operator}\tag{13}\] denotes the displaced parity operator. Here, \[D(\alpha) =\exp\!\left(\alpha a^{\dagger}-\alpha^{*}a\right)\] is the displacement operator and \[\Pi =(-1)^{a^{\dagger}a}\] is the photon-number parity operator. Throughout this work, we adopt the phase-space convention \(\alpha=({q+ip})/{\sqrt2}\), such that \(dq\,dp \equiv 2 d^2\alpha\). We freely switch between the notations \(W(\alpha)\) and \(W(q,p)\) using this convention. In this convention, the vacuum state has the Wigner function \[W_0(q,p) = \frac{1}{\pi} e^{-q^2-p^2},\] and the \(n\)-th Wigner moment is defined as \[w_n =2\int d^{2}\alpha\,W(\alpha)^n . \label{eq:wn95definition}\tag{14}\] Substituting Eq. (12 ) into Eq. (14 ) yields \[\begin{align} w_n &= 2\int d^{2}\alpha \left[ \mathrm{Tr}\!\left( \rho\Delta(\alpha) \right) \right]^n . \label{eq:wn95step1} \end{align}\tag{15}\] Using the identity \[\left[ \mathrm{Tr}(XY) \right]^n =\mathrm{Tr} \!\left[ X^{\otimes n} Y^{\otimes n} \right],\] valid for arbitrary operators \(X\) and \(Y\), Eq. (15 ) can be rewritten as \[\begin{align} w_n &= 2\int d^{2}\alpha\, \mathrm{Tr} \!\left[ \rho^{\otimes n} \Delta(\alpha)^{\otimes n} \right] \nonumber\\ &= \mathrm{Tr} \!\left[ \rho^{\otimes n} \int 2d^{2}\alpha\, \Delta(\alpha)^{\otimes n} \right]. \end{align}\] We therefore introduce the operator \[A_n =2\int d^{2}\alpha\, \Delta(\alpha)^{\otimes n}, \label{eq:An95definition}\tag{16}\] which allows the \(n\)-th Wigner moment to be expressed as \[w_n =\mathrm{Tr} \!\left[ \rho^{\otimes n}A_n \right]. \label{eq:wn95operator95form}\tag{17}\]

Equation (17 ) provides a multi-copy representation of the Wigner moments. Rather than viewing \(w_n\) as a nonlinear functional of the Wigner function, one may equivalently regard it as the expectation value of a state-independent observable \(A_n\) acting on \(n\) identical copies of the quantum state. Consequently, the problem of measuring Wigner moments reduces to determining the structure of the operator \(A_n\).

Position-space kernel of the moment operator↩︎

To determine the operator \(A_n\), we evaluate its kernel in the position representation. We begin with the position-space kernel of the displaced parity operator. In the convention adopted throughout this work, it is given by \[\langle x|\Delta(q,p)|y\rangle = \frac{1}{\pi} e^{ip(x-y)} \delta(x+y-2q). \label{eq:delta95kernel}\tag{18}\]

Using Eq. (18 ), the kernel of \(A_n\) can be written as \[\begin{align} \langle x_1,\ldots,x_n| A_n |y_1,\ldots,y_n\rangle &= \int dq\,dp \prod_{j=1}^{n} \langle x_j| \Delta(q,p) |y_j\rangle \nonumber\\ &= \frac{1}{\pi^n} \int dq\,dp\, e^{ip\sum_j(x_j-y_j)} \prod_j\delta(x_j+y_j-2q). \label{eq:kernel95step1} \end{align}\tag{19}\]

The integration over the momentum variable can be performed immediately. Using \[\int dp\,e^{ipS} = 2\pi\,\delta(S), \label{eq:p95integral}\tag{20}\] where \[S = \sum_{j=1}^{n}(x_j-y_j),\] we obtain \[\begin{align} & \langle x_1,\ldots,x_n| A_n |y_1,\ldots,y_n\rangle = \frac{2}{\pi^{n-1}} \delta\!\left( \sum_{j=1}^{n}(x_j-y_j) \right) \int dq \prod_{j=1}^{n} \delta(x_j+y_j-2q). \label{eq:kernel95step2} \end{align}\tag{21}\]

Using one of the delta functions to perform the \(q\)-integration yields \[\int dq \prod_{j=1}^{n} \delta(x_j+y_j-2q) = \frac{1}{2} \prod_{j=2}^{n} \delta\!\Big[ (x_j+y_j)-(x_1+y_1) \Big].\]

The remaining delta functions constrain all quantities \(x_j+y_j\) to be equal. Explicitly, \[x_1+y_1 = x_2+y_2= \cdots =x_n+y_n. \label{eq:equal95sum95constraint}\tag{22}\] Let the common value of these expressions be denoted by \(s\). Then \[x_j=s-y_j, \qquad j=1,\ldots,n. \label{eq:xj95relation1}\tag{23}\]

Substituting Eq. (23 ) into the remaining constraint appearing in Eq. (21 ) gives \[\begin{align} 0 &= \sum_{j=1}^{n}(x_j-y_j) \nonumber\\ &= \sum_{j=1}^{n}(s-2y_j) \nonumber\\ &= ns-2\sum_{j=1}^{n}y_j. \end{align}\] Hence, \[s = \frac{2}{n} \sum_{j=1}^{n}y_j. \label{eq:s95value}\tag{24}\]

Substituting Eq. (24 ) back into Eq. (23 ), we obtain \[x_j = \frac{2}{n} \sum_{k=1}^{n}y_k -y_j, \qquad j=1,\ldots,n. \label{eq:xj95final}\tag{25}\]

It is convenient to collect these relations into a matrix equation. Defining \[\mathbf{x} = \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{pmatrix}, \qquad \mathbf{y} = \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix},\] together with the matrix \[M_n = \frac{2}{n}J-I, \label{eq:Mn95definition}\tag{26}\] where \(J\) denotes the \(n\times n\) matrix whose entries are all equal to unity, Eq. (25 ) assumes the compact form \[\mathbf{x} = M_n\mathbf{y}. \label{eq:x95equals95My}\tag{27}\]

Equation (27 ) reveals the geometric structure of the operator \(A_n\). The kernel of \(A_n\) is supported only on configurations related by the linear transformation \(M_n\). In the following subsection, we show that this transformation is a reflection in the \((n-1)\)-dimensional relative-coordinate subspace and derive the exact kernel of \(A_n\), including the normalization factor arising from the associated Jacobian.

Reduction to a reflection operator↩︎

Equation (27 ) shows that the support of the kernel of \(A_n\) is restricted to points satisfying the linear relation \[\mathbf{x} =M_n\mathbf{y}, \qquad M_n = \frac{2}{n}J-I.\] To obtain the exact form of the kernel, it remains to determine the Jacobian associated with the transformation from the original set of delta-function constraints to the compact matrix constraint above.

Using the result of the previous subsection, the kernel of \(A_n\) can be written as \[\langle \mathbf{x}|A_n|\mathbf{y}\rangle = \frac{1}{\pi^{\,n-1}} \delta\!\left( \sum_j(x_j-y_j) \right) \int dq \prod_j \delta(x_j+y_j-2q). \label{eq:kernel95before95jacobian}\tag{28}\] The constraints appearing in Eq. (28 ) can be expressed in terms of the functions \[\begin{align} g_1(\mathbf{x},\mathbf{y}) &= \sum_{j=1}^{n}(x_j-y_j), \tag{29} \\ g_i(\mathbf{x},\mathbf{y}) &= x_i+y_i-x_1-y_1, \qquad i=2,\ldots,n. \tag{30} \end{align}\] The kernel therefore contains the product of delta functions \[\delta(g_1) \prod_{i=2}^{n} \delta(g_i). \label{eq:delta95constraints}\tag{31}\]

As shown in the previous subsection, the simultaneous solution of these constraints is precisely \[\mathbf{x} = M_n\mathbf{y}. \label{eq:constraint95solution}\tag{32}\] The transformation from the variables \((g_1,\ldots,g_n)\) to \((x_1,\ldots,x_n)\) introduces a Jacobian factor. To evaluate it, we compute the matrix \[J_g =\left( \frac{\partial g_i}{\partial x_j} \right).\] Using Eqs. (29 ) and (30 ), one finds \[J_g =\begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ -1 & 1 & 0 & \cdots & 0 \\ -1 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1 & 0 & 0 & \cdots & 1 \end{pmatrix}. \label{eq:jacobian95matrix}\tag{33}\]

The determinant of this matrix can be evaluated conveniently using the Schur complement. Writing Eq. (33 ) in block form, \[J_g =\begin{pmatrix} 1 & \mathbf{1}^{T} \\ -\mathbf{1} & I_{n-1} \end{pmatrix},\] where \(\mathbf{1}\) denotes the \((n-1)\)-dimensional column vector with all entries equal to unity and \(I_{n-1}\) is the \((n-1)\times(n-1)\) identity matrix, we obtain \[\begin{align} \det J_g &= \det(I_{n-1}) \, \det\!\left( 1-\mathbf{1}^{T}I_{n-1}^{-1}(-\mathbf{1}) \right) \nonumber\\ &= 1+\mathbf{1}^{T}\mathbf{1} \nonumber\\ &= 1+(n-1) \nonumber\\ &= n. \label{eq:jacobian95determinant} \end{align}\tag{34}\]

The multidimensional delta-function identity therefore gives \[\delta(g_1) \prod_{i=2}^{n} \delta(g_i) =\frac{1}{|\det J_g|} \, \delta(\mathbf{x}-M_n\mathbf{y}),\] which, together with Eq. (34 ), yields \[\delta(g_1) \prod_{i=2}^{n} \delta(g_i) =\frac{1}{n} \delta(\mathbf{x}-M_n\mathbf{y}). \label{eq:delta95reduction}\tag{35}\]

Substituting Eq. (35 ) into Eq. (28 ), we finally obtain \[\langle \mathbf{x}|A_n|\mathbf{y}\rangle = \frac{1}{n\pi^{\,n-1}} \delta(\mathbf{x}-M_n\mathbf{y}). \label{eq:kernel95final95result}\tag{36}\]

Equation (36 ) reveals that the action of \(A_n\) is entirely determined by the matrix \(M_n\). The operator maps a configuration \(\mathbf{y}\) to the transformed configuration \(M_n\mathbf{y}\) and is therefore generated by a linear orthogonal transformation in the \(n\)-dimensional configuration space. The structure of this transformation becomes particularly transparent after diagonalizing \(M_n\). As we show in the next subsection, \(M_n\) acts trivially on a single collective coordinate while reversing the sign of all relative coordinates. Consequently, \(A_n\) can be interpreted as a reflection operator in the \((n-1)\)-dimensional relative-coordinate subspace.

Exact operator form of \(A_n\): Proof of Theorem 5↩︎

The structure of the kernel obtained in Eq. (36 ) is entirely determined by the matrix \[M_n =\frac{2}{n}J-I. \label{eq:Mn95recalled}\tag{37}\] To understand the action of the operator \(A_n\), it is therefore necessary to analyze the spectral properties of \(M_n\).

We begin by recalling that the matrix \(J\) has rank one and can be written as \[J = \mathbf{u}\mathbf{u}^{T},\] where \[\mathbf{u} =(1,1,\ldots,1)^{T}.\] Consequently, \[J\mathbf{u} =n\mathbf{u},\] while every vector orthogonal to \(\mathbf{u}\) belongs to the kernel of \(J\). Thus, the spectrum of \(J\) consists of a single nonzero eigenvalue \(n\) and \(n-1\) vanishing eigenvalues, \[\mathrm{spec}(J) =\{n,0,\ldots,0\}.\]

Since \(M_n\) is related to \(J\) through Eq. (37 ), the corresponding eigenvalues are obtained immediately, \[\mathrm{spec}(M_n) =\left\{ \frac{2}{n}n-1, \frac{2}{n}\times 0-1, \ldots, \frac{2}{n}\times 0-1 \right\},\] which simplifies to \[\mathrm{spec}(M_n) =\{1,-1,\ldots,-1\}. \label{eq:spectrum95Mn}\tag{38}\]

Equation (38 ) shows that \(M_n\) possesses a unique eigendirection with eigenvalue \(+1\), while every direction orthogonal to it acquires eigenvalue \(-1\).

The distinguished eigenvector is \[\mathbf{e}_0 =\frac{1}{\sqrt n} (1,1,\ldots,1)^T,\] which naturally defines the collective coordinate \[Q_0 =\mathbf{e}_0^{T}\mathbf{x} =\frac{1}{\sqrt n} \sum_{j=1}^{n}x_j. \label{eq:collective95coordinate}\tag{39}\] The remaining \(n-1\) coordinates are chosen as an arbitrary orthonormal basis spanning the subspace orthogonal to \(\mathbf{e}_0\). Denoting these coordinates by \[Q_1,Q_2,\ldots,Q_{n-1},\] the transformation from the original coordinates to the new coordinates may be written as \[\mathbf{Q} =F_n\mathbf{x},\] where \(F_n\) is an orthogonal matrix whose first row is given by \(\mathbf{e}_0^{T}\).

In this basis, the matrix \(M_n\) becomes diagonal, \[F_nM_nF_n^{T} =\mathrm{diag}(1,-1,\ldots,-1). \label{eq:Mn95diagonal}\tag{40}\]

Equation (40 ) immediately reveals the geometric meaning of \(M_n\). The collective coordinate remains unchanged, \[Q_0 \longmapsto Q_0,\] whereas every relative coordinate changes sign, \[Q_j \longmapsto -Q_j, \qquad j=1,\ldots,n-1.\] Thus, \(M_n\) acts as the identity on the one-dimensional collective subspace and as a reflection on the \((n-1)\)-dimensional relative-coordinate subspace.

Using Eq. (40 ), the kernel in Eq. (36 ) can be expressed in the transformed coordinates. Since \(F_n\) is orthogonal, the multidimensional Dirac delta distribution remains unchanged under the coordinate transformation. Consequently, \[\begin{align} & \langle Q_0,\ldots,Q_{n-1} | A_n | Q_0',\ldots,Q_{n-1}' \rangle \nonumber\\ &= \frac{1}{n\pi^{n-1}} \delta(Q_0-Q_0') \prod_{j=1}^{n-1} \delta(Q_j+Q_j'). \label{eq:kernel95Q} \end{align}\tag{41}\]

The first factor in Eq. (41 ) is precisely the position-space kernel of the identity operator, \[\langle Q|I|Q'\rangle =\delta(Q-Q').\] The remaining factors can be identified with the position-space kernel of the parity operator, \[\langle Q|\Pi|Q'\rangle =\delta(Q+Q'). \label{eq:parity95kernel}\tag{42}\] Therefore, \[\delta(Q_0-Q_0') \prod_{j=1}^{n-1} \delta(Q_j+Q_j') =\left\langle \mathbf{Q} \left| I\otimes\Pi^{\otimes(n-1)} \right| \mathbf{Q}' \right\rangle .\]

Substituting this identity into Eq. (41 ) yields the operator representation \[A_n =\frac{1}{n\pi^{\,n-1}} F_n^{\dagger} \left( I\otimes \Pi^{\otimes(n-1)} \right) F_n, \label{eq:An95final}\tag{43}\] which proves Theorem 5.

Since \(I\otimes\Pi^{\otimes(n-1)}\) acts identically on all relative modes, the resulting operator \(A_n\) is independent of the specific orthonormal basis chosen for the relative-coordinate subspace. The theorem provides an exact operator representation of all Wigner moments. Remarkably, the structure of \(A_n\) is universal: irrespective of the value of \(n\), the operator consists of a parity measurement acting on the relative-coordinate subspace, while the collective mode contributes only an identity operator.

Experimental feasibility↩︎

The operator representation derived in Theorem 5 provides a direct experimental protocol for measuring arbitrary Wigner moments. In contrast to conventional approaches based on quantum state tomography, the present method allows the moments to be estimated directly from expectation values of a multi-copy observable. Consider the \(n\)-th Wigner moment, \[w_n =\frac{1}{n\pi^{\,n-1}} \mathrm{Tr} \!\left[ \rho^{\otimes n} F_n^{\dagger} \left( I\otimes \Pi^{\otimes(n-1)} \right) F_n \right] = \frac{1}{n\pi^{\,n-1}} \left\langle I\otimes \Pi^{\otimes(n-1)} \right\rangle_{F_n\rho^{\otimes n}F_n^{\dagger}}. \label{eq:wn95measurement}\tag{44}\] Equation (44 ) immediately suggests a three-step measurement protocol.

1. The first step consists of preparing \(n\) identical copies of the quantum state. The observable associated with the \(n\)-th Wigner moment acts on the tensor-product Hilbert space \(\mathcal{H}^{\otimes n}\), and therefore requires simultaneous access to \(n\) copies of the state \(\rho\).

2. The second step is the implementation of the orthogonal interferometric transformation \(F_n\). As shown in the previous subsection, the role of \(F_n\) is to separate the collective degree of freedom from the relative degrees of freedom. In the transformed basis, the first mode corresponds to the collective coordinate \(Q_0 =\frac{1}{\sqrt n} \sum_{j=1}^{n}x_j\), while the remaining modes span the relative-coordinate subspace. Physically, the transformation \(F_n\) can be realized using a passive Gaussian interferometric network implementing the corresponding orthogonal transformation of the quadratures. Such networks may be decomposed into beam splitters and phase shifters according to standard linear-optical constructions.

3. The final step is the measurement of parity on the relative modes. Since the operator acting in the transformed basis is \(I\otimes \Pi^{\otimes(n-1)}\), the collective mode contributes only an identity operator and therefore does not affect the measurement outcome. All relevant information is contained in the relative modes. The measurement protocol therefore requires estimation of the expectation value of the joint parity observable \[\Pi^{\otimes(n-1)} =\Pi_1\Pi_2\cdots\Pi_{n-1},\] acting on the relative-coordinate subspace. Combining these three steps yields an estimator for the Wigner moment, \[w_n =\frac{1}{n\pi^{\,n-1}} \left\langle I\otimes \Pi^{\otimes(n-1)} \right\rangle_{F_n\rho^{\otimes n}F_n^{\dagger}}.\]

Consequently, the determination of the Wigner moments requires neither reconstruction of the Wigner function nor full quantum state tomography. Instead, each moment is obtained directly from a single expectation value measured after an interferometric processing of \(n\) identical copies of the state. The physical interpretation of this result is particularly transparent. The transformation \(F_n\) isolates a single collective mode that remains invariant under the reflection generated by the operator \(A_n\). All nontrivial information about the Wigner moment is encoded in the relative-coordinate subspace, where the action of \(A_n\) reduces to a simultaneous parity measurement. The \(n\)-th Wigner moment therefore quantifies the response of the \(n\)-copy state to a reflection of all relative degrees of freedom while leaving the collective degree of freedom unchanged.

Theorem 5 thus establishes that the complete hierarchy of Wigner moments can, in principle, be accessed through multicopy interferometry followed by parity measurements on the relative modes. Since the protocol directly estimates the desired moment without reconstructing the underlying phase-space distribution, it provides an experimentally motivated alternative to full state tomography for characterizing phase-space properties of continuous-variable quantum states.

The experimental protocol described above requires the implementation of an orthogonal transformation \(F_n\) that separates the collective degree of freedom from the relative degrees of freedom of the \(n\)-copy state. In general, many such transformations are possible. The only requirement is that the first output mode coincide with the normalized collective coordinate \(Q_0\). while the remaining output modes form an orthonormal basis of the (n-1)-dimensional relative-coordinate subspace. For low-order moments, particularly simple choices of \(F_n\) may be constructed analytically. These examples illustrate how the collective and relative coordinates arise from the original position coordinates of the n-copies.

For \(n=2\), we choose \[F_2 \equiv\frac{1}{\sqrt2} \begin{pmatrix} 1 & 1\\ 1 & -1 \end{pmatrix}.\] The transformed coordinates are \[Q_0 =\frac{x_1+x_2}{\sqrt2}, \qquad Q_1 =\frac{x_1-x_2}{\sqrt2}.\] The coordinate \(Q_0\) describes the collective motion of the two copies, while \(Q_1\) measures their relative displacement. This transformation is precisely the one implemented by a balanced \(50:50\) beam splitter.

For \(n=3\), a convenient choice is \[F_3 \equiv \begin{pmatrix} \dfrac1{\sqrt3} & \dfrac1{\sqrt3} & \dfrac1{\sqrt3} \\[10pt] \dfrac1{\sqrt2} & -\dfrac1{\sqrt2} & 0 \\[10pt] \dfrac1{\sqrt6} & \dfrac1{\sqrt6} & -\sqrt{\dfrac23} \end{pmatrix}\] The corresponding coordinates are \[\begin{align} Q_0 = \frac{x_1+x_2+x_3}{\sqrt3}, \quad Q_1 = \frac{x_1-x_2}{\sqrt2}, \quad Q_2 = \frac{x_1+x_2-2x_3}{\sqrt6}. \end{align}\]

For \(n=4\), we may choose \[F_4 \equiv\; \begin{pmatrix} \dfrac12 & \dfrac12 & \dfrac12 & \dfrac12 \\[4mm] \dfrac1{\sqrt2} & -\dfrac1{\sqrt2} & 0 & 0 \\[4mm] \dfrac12 & \dfrac12 & -\dfrac12 & -\dfrac12 \\[4mm] 0 & 0 & \dfrac1{\sqrt2} & -\dfrac1{\sqrt2} \end{pmatrix}\] The transformed coordinates become \[\begin{align} Q_0 = \frac{x_1+x_2+x_3+x_4}{2}, \quad Q_1 = \frac{x_1-x_2}{\sqrt2}, \quad Q_2 = \frac{x_1+x_2-x_3-x_4}{2}, \quad Q_3 = \frac{x_3-x_4}{\sqrt2}. \end{align}\] In each case, the rows of \(F_n\) are orthonormal and satisfy \[F_nF_n^{T} =I.\] Furthermore, the transformations diagonalize the reflection matrix (\(M_n\)), \[F_nM_nF_n^{T} =\mathrm{diag}(1,-1,\ldots,-1).\] Consequently, the first output mode is left invariant by the action of (\(M_n\)), whereas every relative mode acquires a sign inversion. This is precisely the structure responsible for the appearance of the operator \[I\otimes\Pi^{\otimes(n-1)}\] in Theorem 5. The transformations \(F_2\), \(F_3\), and \(F_4\) provide explicit examples of collective-coordinate transformations entering the moment-measurement protocol.

8 H. Applications: Detection of Quantum Entanglement via Wigner Moments↩︎

Recent advances have uncovered direct links between entanglement and Wigner negativity [45], [46], shedding light on the structure of quantum correlations. In particular, negativities in appropriately constructed Wigner functions, such as reduced or collective phase-space distributions have been shown to imply the presence of bipartite entanglement or even genuine multipartite entanglement. We exploit this connection to lift our moment-based nonclassicality criteria to the detection of quantum entanglement.

8.1 Detection of bipartite entanglement↩︎

First we demonstrate how the moment-based criteria developed in the previous sections can be employed to detect bipartite entanglement in CV systems utilizing the recently established connections between Wigner negativity and the negativity of the partial transposed density matrix.

Let \(\rho\) be an \(N\)-mode CV state. We consider an equal bipartition of the modes into two subsystems, \[\label{eq:equal95bipartition} \vec{a}_A := \{a_m\}_{m=1}^{N/2}, \qquad \vec{a}_B := \{a_m\}_{m=N/2+1}^{N}.\tag{45}\] We denote the partial transpose of \(\rho\) over the modes \(\vec{a}_B\) as \(\rho^{T_B}\). Then entanglement across this bipartition can be characterized using the negativity of the partial transpose criterion: if the partially transposed state \(\rho^{T_B}\) is not positive semidefinite, then \(\rho\) is necessarily entangled.

A powerful connection arises from a phase-space formulation of the PPT condition. In particular, it has been shown that bipartite entanglement can be inferred from negativities in suitably defined reduced Wigner functions. To describe this construction, we introduce collective modes defined componentwise as \[\vec{a}_{\pm} = \frac{1}{\sqrt{2}}\left(\vec{a}_A \pm \vec{a}_B\right),\]

Given the state \(\rho\), we calculate the reduced state \(\mathrm{Tr}_{-}\rho\), obtained by tracing out all relative modes \(\vec{a}_-\). The Wigner function of this reduced state, denoted \(W_{\mathrm{Tr}_{-}\rho}\), is obtained from the full Wigner function \(W_\rho\) via marginalization over the phase-space variables associated with \(\vec{a}_-\). A key result, proven in Ref. [45], establishes that negativity of this reduced Wigner function implies violation of the PPT condition. More precisely, one has \[W_{\mathrm{Tr}_{-}\rho}(\vec{\alpha}) \not\geq 0 \;\;\Longrightarrow\;\; \rho^{T_B} \not\succeq 0,\] which certifies NPT entanglement across the chosen bipartition. This result provides a direct operational link between phase-space nonclassicality and bipartite entanglement.

Importantly, this approach reduces the problem of entanglement detection to identifying negativities in a reduced phase-space distribution. However, directly resolving such negativities typically requires high-resolution phase-space measurements. This motivates the use of moment-based criteria, which provide coarse-grained yet experimentally accessible signatures of Wigner negativity. Combining these observations, we obtain the following result.

The above theorem demonstrates that moment-based phase-space criteria can be systematically promoted to entanglement witnesses. A key advantage of this approach is that it requires only partial information about the Wigner function encoded in a finite set of moments, rather than a full phase-space reconstruction. This makes the resulting witnesses particularly suitable for experimental platforms where phase-space measurements are available but limited.

8.2 Detection of genuine multipartite entanglement↩︎

We now extend the above framework to the multipartite regime, where the structure of quantum correlations becomes significantly richer. In contrast to bipartite entanglement, multipartite systems admit different layers of correlations depending on how the system can be partitioned. The strongest form of such correlations is genuine multipartite entanglement (GME), which captures the presence of inseparability across all possible bipartitions.

Formally, an \(N\)-mode state \(\rho\) is said to be genuinely multipartite entangled if it cannot be written as a convex mixture of states that are separable with respect to any bipartition, i.e., \[\rho \neq \sum_{(A|\bar{A})} \sum_i p_{A}^{(i)} \, \rho_{A}^{(i)} \otimes \rho_{\bar{A}}^{(i)}, \qquad \text{where}\;\;\; p_{A}^{(i)} \geq 0.\] \(A (\bar{A})\) denote a subset of modes i.e., \(A= \{ m_n\}_{n=1}^{M}\) (complementary subset of modes, \(\bar{A} = \{m_n\}_{n=1}^{M} \setminus A\)) and runs over all bipartitions \((A|\bar{A})\). Such states exhibit correlations that cannot be reproduced by classical coordination across any grouping of subsystems, and play a central role in tasks such as quantum networks, distributed computation, and multipartite quantum communication.

Detecting GME is considerably more challenging than bipartite entanglement. Standard criteria, including PPT-based methods or covariance-matrix approaches, are often either insufficient or experimentally demanding, particularly for non-Gaussian states. This motivates the search for experimentally accessible criteria that can capture genuinely multipartite correlations using limited information.

Phase-space approach and centre-of-mass mode. Recent developments have shown that phase-space methods provide a natural route to GME detection by exploiting collective degrees of freedom of multimode systems. In particular, one introduces a centre-of-mass mode of the form \[a_+ = \frac{1}{\sqrt{M}} \sum_{m=1}^M \left( y_m a_m + z_m a_m^\dagger \right),\] where the coefficients \(\vec{y}, \vec{z} \in \mathbb{C}^M\) are choosen such that \(\vec{y}\circ\vec{y}^* - \vec{z}\circ\vec{z}^* = \vec{1}\), where \([\mathbf{A}\circ\mathbf{B}]_{m,n} = [\mathbf{A}]_{m,n}[\mathbf{B}]_{m,n}\) is the elementwise product and \(\vec{1} = (1,1,\dots,1)\) is a vector of ones. Given the full state \(\rho\), one finds the reduced state describing the centre-of-mass motion by tracing out all orthogonal (relative) modes, \(\mathrm{Tr}_{-}\rho\), which retains only the collective degree of freedom associated with \(a_+\). The corresponding Wigner function \(W_{\mathrm{Tr}_{-}\rho}(\alpha)\) can be obtained from the full Wigner function \(W_\rho(\vec{\beta})\) by marginalizing over all relative coordinates.

A key result, established in recent works (see Ref. [46]), is that negativity of a suitably processed Wigner function of the centre-of-mass mode provides a sufficient condition for GME. Since the reduced Wigner function may undergo smoothing due to marginalization, one considers a smoothed Wigner function defined via convolution, \[\widetilde{W}_{\mathrm{Tr}_{-}\rho}(\alpha \;; \mathcal{R}) = \int d^{2}{\beta} \;\; W_{\mathrm{Tr}_{-}\rho}(\beta) \, K(\alpha - \beta \,; \mathcal{R}),\] where \(K(\alpha - \beta \,; \mathcal{R})\) is an appropriate kernel determined by auxiliary states or filtering procedures.

The central implication can be stated as follows: if the negativities present in the centre-of-mass Wigner function are sufficiently robust to persist under such smoothing, then the underlying state cannot be decomposed into biseparable components. In particular, \[\exists \alpha \;:\; \widetilde{W}_{\mathrm{Tr}_{-}\rho}(\alpha \,; \mathcal{R}) < 0 \quad \Longrightarrow \quad \rho \;\text{is GME}.\]

This result establishes a direct and operational bridge between phase-space nonclassicality and genuine multipartite entanglement. Importantly, it shifts the problem of GME detection to identifying negativities in a single-mode quasiprobability distribution associated with collective degrees of freedom.

Moment-based detection of GME. While the above criterion is powerful, directly resolving negativities in \(\widetilde{W}_{\mathrm{Tr}_{-}\rho}\) may still require fine-grained phase-space measurements. This is precisely where the moment-based hierarchies developed in this work become useful. Since the \(L_p\)-norm, log-convexity, and Hankel-matrix criteria provide sufficient conditions for Wigner negativity, they can be directly applied to the smoothed centre-of-mass Wigner function to yield experimentally accessible witnesses of GME. We formalize this observation in the following theorem.

This result demonstrates that moment-based phase-space criteria provide a unified and experimentally viable framework for detecting multipartite entanglement. In contrast to conventional GME witnesses, which often require full tomography or access to high-order correlations, our approach relies only on low-order moments of a single effective mode. This significantly reduces the experimental overhead, particularly in platforms where phase-space measurements are naturally available.

8.3 Illustrative examples↩︎

We now illustrate the above results using representative families of multimode states. In particular, we consider parametrized non-Gaussian states for which the reduced Wigner function can be computed explicitly, and compare the detection thresholds obtained from moment-based criteria with those obtained from direct Wigner negativity and the PPT condition.

Example 1: Entangled single-photon state. We consider the bipartite single-photon state \[\ket{\theta} = \cos\theta\, \ket{1,0} + \sin\theta\, \ket{0,1},\;\; 0\leq \theta\leq2\pi\] where \(\ket{1,0} = a_1^\dagger \ket{0}\) and \(\ket{0,1} = a_2^\dagger \ket{0}\) denote single-photon excitations in modes \(1\) and \(2\), respectively. We introduce the collective modes \[a_\pm = \frac{a_1 \pm a_2}{\sqrt{2}},\] which can be inverted as \[a_1 = \frac{a_+ + a_-}{\sqrt{2}}, \qquad a_2 = \frac{a_+ - a_-}{\sqrt{2}}.\] The vacuum state is invariant under this change of modes, i.e., \(a_\pm \ket{0} = 0\). We now rewrite the basis states in terms of the new modes. Using the above relations, \[\begin{align} \ket{1,0} &= a_1^\dagger \ket{0}\\ &= \frac{1}{\sqrt{2}}(a_+^\dagger + a_-^\dagger)\ket{0}\\ &= \frac{1}{\sqrt{2}}\big(\ket{1}_+ \ket{0}_- + \ket{0}_+ \ket{1}_-\big), \\ \ket{0,1} &= a_2^\dagger \ket{0}\\ &= \frac{1}{\sqrt{2}}(a_+^\dagger - a_-^\dagger)\ket{0}\\ &= \frac{1}{\sqrt{2}}\big(\ket{1}_+ \ket{0}_- - \ket{0}_+ \ket{1}_-\big), \end{align}\] where we have used the definitions \[\begin{align} \ket{1}_+ = a_+^\dagger \ket{0}, \quad \ket{1}_- = a_-^\dagger \ket{0}. \end{align}\]

Substituting these expressions into \(\ket{\theta}\), we obtain \[\begin{align} \ket{\theta} &= \frac{1}{\sqrt{2}} \Big[ \cos\theta \big(\ket{1}_+ \ket{0}_- + \ket{0}_+ \ket{1}_-\big) \nonumber + \sin\theta \big(\ket{1}_+ \ket{0}_- - \ket{0}_+ \ket{1}_-\big) \Big] \nonumber \\ &= \frac{1}{\sqrt{2}} \Big[ (\cos\theta + \sin\theta)\, \ket{1}_+ \ket{0}_- \nonumber + (\cos\theta - \sin\theta)\, \ket{0}_+ \ket{1}_- \Big]. \end{align}\] Defining \[c_+ = \frac{\cos\theta + \sin\theta}{\sqrt{2}}, \qquad c_- = \frac{\cos\theta - \sin\theta}{\sqrt{2}},\] the state can be written compactly as \[\ket{\theta} = c_+ \ket{1}_+ \ket{0}_- + c_- \ket{0}_+ \ket{1}_-.\]

We now compute the reduced state obtained by tracing out the \(-\) mode: \[\rho_+ = \text{Tr}_- \big( \ket{\theta}\bra{\theta} \big).\] Expanding the density matrix, \[\begin{align} \ket{\theta}\bra{\theta} &= |c_+|^2 \ket{1}_+\bra{1} \otimes \ket{0}_-\bra{0} + |c_-|^2 \ket{0}_+\bra{0} \otimes \ket{1}_-\bra{1} \nonumber + c_+ c_-^* \ket{1}_+\bra{0} \otimes \ket{0}_-\bra{1} + c_+^* c_- \ket{0}_+\bra{1} \otimes \ket{1}_-\bra{0}. \end{align}\] Taking the partial trace over the \(-\) mode we get, \[\rho_+ = |c_+|^2 \ket{1}_+ \bra{1} + |c_-|^2 \ket{0}_+\bra{0}.\]

The Wigner function of the reduced state is therefore a convex combination of the vacuum and single-photon Wigner functions, \[W_{\text{Tr}_-\rho}(x_+,p_+) = |c_-|^2 W_0(x_+,p_+) + |c_+|^2 W_1(x_+,p_+),\] which is equivalent to the example of mixed Fock states presented in the main text, with \(\lambda = |c_-|^2\). Hence \(W_{\text{Tr}_-\rho}(x_+,p_+)\) is Wigner negative whenever, \[\begin{align} |c_-|^2<0.5 \; &\Rightarrow \sin (2\theta)>0 \nonumber\\ & \Rightarrow \theta \in \left(0,\frac{\pi}{2}\right) \cup \left(\pi,\frac{3\pi}{2}\right) \end{align}\] Hence, for a range of \(\theta\), the reduced Wigner function exhibits negativity, which, by Proposition 4, certifies the presence of bipartite entanglement in the original state, and that entire range of \(\theta\) can be detected by the moment based criteria (see Table. 1).

Example 2: Tripartite Dicke state. As an illustrative example of the detection of GME using Wigner moments, we consider the tripartite two-excitation Dicke state \[|D_2^3\rangle =\frac{1}{\sqrt3} \Bigl( |110\rangle + |101\rangle + |011\rangle \Bigr).\]

Following Ref. [47], we construct the reduced Wigner function using the kernel specified by \(\mathcal{R}={\ket{1}}\). The resulting reduced Wigner function is \[\widetilde{W}_{D_2^3}(\alpha;{\ket{1}}) =\frac{e^{-\frac{3}{2}|\alpha|^2}}{32\pi} \left( 81|\alpha|^6 -234|\alpha|^4 + 216|\alpha|^2 -16 \right).\]

First we calculate the Hankel matrix of order \(1\), whose determinant evaluates to \[\det(H_1) = 5.13\times10^{-3}>0.\] Therefore, the lowest-order moment condition is inconclusive for this state. Proceeding to the next level of the hierarchy, we find

\[\det(H_2) =-2.23\times10^{-7}<0.\] Hence, the GME of \(|D_2^3\rangle\) is therefore certified.

References↩︎