Revealing nonclassicality beyond Gaussian states via a single marginal distribution

TL;DR

This paper introduces a novel method to detect nonclassicality in quantum states using only a single marginal distribution, applicable to both Gaussian and non-Gaussian states, with experimental validation on trapped ion motional states.

Contribution

The authors propose a framework that reveals nonclassicality from a single marginal distribution, providing a practical and axis-independent detection method for various quantum states.

Findings

Successfully demonstrated nonclassicality detection in trapped ion states

The method detects both Gaussian and non-Gaussian nonclassical states reliably

Provides a lower bound for entanglement potential using a single measurement axis

Abstract

A standard method to obtain information on a quantum state is to measure marginal distributions along many different axes in phase space, which forms a basis of quantum state tomography. We theoretically propose and experimentally demonstrate a general framework to manifest nonclassicality by observing a single marginal distribution only, which provides a novel insight into nonclassicality and a practical applicability to various quantum systems. Our approach maps the 1-dim marginal distribution into a factorized 2-dim distribution by multiplying the measured distribution or the vacuum-state distribution along an orthogonal axis. The resulting fictitious Wigner function becomes unphysical only for a nonclassical state, thus the negativity of the corresponding density operator provides an evidence of nonclassicality. Furthermore, the negativity measured this way yields a lower bound for entanglement potential---a measure of entanglement generated using a nonclassical state with a beam splitter setting that is a prototypical model to produce continuous-variable (CV) entangled states. Our approach detects both Gaussian and non-Gaussian nonclassical states in a reliable and efficient manner. Remarkably, it works regardless of measurement axis for all non-Gaussian states in finite-dimensional Fock space of any size, also extending to infinite-dimensional states of experimental relevance for CV quantum informatics. We experimentally illustrate the power of our criterion for motional states of a trapped ion confirming their nonclassicality in a measurement-axis independent manner. We also address an extension of our approach combined with phase-shift operations, which leads to a stronger test of nonclassicality, i.e. detection of genuine non-Gaussianity under a CV measurement.

Equations

$${\mathcal{D}}_1$$

$${\mathcal{D}}_2$$

$$W_{{\rho }_{\mathrm{cl}}}(x,y)=\int d^2\beta P({\beta }_x,{\beta }_y)W_{|\beta ⟩⟨\beta |}(x,y),$$

$${\mathcal{D}}_j[W_{{\rho }_{\mathrm{cl}}}(q,p)]=\int d^2\beta Q_1({\beta }_x)Q_j({\beta }_y)W_{|\beta ⟩⟨\beta |}(x,y)$$

$${\mathcal{D}}_1[W_{\sigma }(q,p)]$$

$${\mathcal{D}}_2[W_{\sigma }(q,p)]$$

$${\mathcal{N}}_{\mathrm{DM}}(\rho )\equiv \underset{\theta \in (0,\pi )}{\max }\frac{||{\rho }_{\mathrm{DM}}^{\theta }|{|}_1-1}{2},$$

$$W_{\rho }(q_1,p_1)W_{\overline{\rho }}(q_2,p_2)$$

$$\stackrel{\mathrm{BS}}{\to }$$

$$W_{\rho }(\frac{q_1+q_2}{\sqrt{2}},\frac{p_1-p_2}{\sqrt{2}})W_{\overline{\rho }}(\frac{q_1-q_2}{\sqrt{2}},\frac{p_1+p_2}{\sqrt{2}})$$

$$\stackrel{\mathrm{BS}}{\to }$$

$${\mathcal{P}}_{\mathrm{ent}}[\rho ]\equiv \frac{||[{\hat{U}}_{\mathrm{BS}}({\rho }_1\otimes |0⟩⟨0{|}_2){\hat{U}}_{\mathrm{BS}}^{†}{]}^{\mathrm{PT}}|{|}_1-1}{2},$$

$${\mathcal{N}}_{DM2}[\rho ]\le {\mathcal{P}}_{\mathrm{ent}}[\rho ].$$

$${\mathcal{M}}_{jk}=C({\xi }_j-{\xi }_k)e^{\frac{1}{2}({\xi }_j{\xi }_k^{\ast }-{\xi }_j^{\ast }{\xi }_k)},$$

$$:(\hat{a}+{\hat{a}}^{†}{)}^m:=\underset{\lambda \to 0}{\lim }\frac{{\partial }^m}{\partial {\lambda }^m}:e^{\lambda (\hat{a}+{\hat{a}}^{†})}:$$

$$=\underset{\lambda \to 0}{\lim }\frac{{\partial }^n}{\partial {\lambda }^m}{e}^{-{\lambda }^2/2}{e}^{\lambda (\hat{a}+{\hat{a}}^{†})}$$

$$=\underset{\lambda \to 0}{\lim }\sum _{k=0}^m\left(\genfrac{}{}{0ex}{}{m}{k}\right)\frac{{\partial }^k}{\partial {\lambda }^k}{e}^{-{\lambda }^2/2}\frac{{\partial }^{m-k}}{\partial {\lambda }^{m-k}}{e}^{\lambda (\hat{a}+{\hat{a}}^{†})}$$

$$=\underset{\lambda \to 0}{\lim }\sum _{k=0}^m\left(\genfrac{}{}{0ex}{}{m}{k}\right)2^{-k/2}{e}^{-{\lambda }^2/2}{H}_k(-\lambda /\sqrt{2})(\hat{a}+{\hat{a}}^{†}{)}^{m-k}{e}^{\lambda (\hat{a}+{\hat{a}}^{†})}$$

$$=2^{-m/2}\sum _{k=0}^m\left(\genfrac{}{}{0ex}{}{m}{k}\right)H_k(0)\{\sqrt{2}(\hat{a}+{\hat{a}}^{†}){\}}^{m-k}$$

$$\displaystyle=2^{-m/2}H_{m}\bigg{(}\frac{\hat{a}+\hat{a}^{{\dagger}}}{\sqrt{2}}\bigg{)},$$

$$\sigma =\hat{D}(\alpha )\hat{S}(r,\varphi ){\sigma }_{\mathrm{th}}(\bar{n}){\hat{S}}^{†}(r,\varphi ){\hat{D}}^{†}(\alpha ),$$

$$W_{\sigma }(q,p)=\frac{2\mu }{\pi }\exp [-2\mu (e^{2r}{q}^2+e^{-2r}{p}^2)],$$

$$M_{\sigma}(x)=\frac{1}{\sqrt{2\pi V_{\theta}}}\exp\bigg{(}-\frac{x^{2}}{2V_{\theta}}\bigg{)},$$

$${\eta }_{\mathrm{sq}}=\frac{|\theta {|}_{\mathrm{sq}}}{\frac{\pi }{2}}=\frac{1}{\pi }\arccos \frac{\cosh 2r-\mu }{\sinh 2r}.$$

$$W_{{\mathcal{D}}_2[\rho ]}(q,p)=M_{\rho }(q)\sqrt{\frac{2}{\pi }}{e}^{-2p^2},$$

$${\mathcal{D}}_2[\rho ]={\int }_{-\infty }^{\infty }dx{T}_{\rho }(x)|x⟩⟨x|.$$

$$⟨\psi |{\mathcal{D}}_2[\rho ]|\psi ⟩$$

$$=\sum _{n=0}^{\infty }\sum _{m=0}^{\infty }\frac{c_n^{\ast }}{\sqrt{n!}}\frac{c_m}{\sqrt{m!}}{\int }_{-\infty }^{\infty }dx{T}_{\rho }(x)e^{-x^2}{x}^{n+m}$$

$$<0.$$

$$\sum _{n=0}^N\sum _{m=0}^N{\beta }_n^{\ast }{\beta }_m{a}_{n+m}\ge 0,$$

$$M_{\rho }(x)=\frac{1}{N}\sum _{k=0}^{N-1}\frac{1}{\sqrt{2\pi v_k}}{e}^{-\frac{x^2}{2v_k}},$$

$$M_{\rho }(x)=\frac{1}{N}\sum _{k=0}^{N-1}\frac{1}{\sqrt{2\pi v_k}}{e}^{-\frac{x^2}{2v_k}}=\sum _{k^{\prime }}\frac{p_{k^{\prime }}}{\sqrt{\pi /2}}{e}^{-2(x-x_{k^{\prime }}{)}^2},$$

$$\frac{1}{N}+\frac{1}{N}\sum _{k\ne k_{\min }}{e}^{-\frac{1}{2}(v_k-v)y^2}=\sum _{k^{\prime }}{p}_{k^{\prime }}{e}^{-i{x}_{k^{\prime }}y}{e}^{-\frac{1}{2}(\frac{1}{4}-v)y^2},$$

$$W_{\sigma}^{\mathrm{DM2}}(x,y)=\frac{1}{\pi\sqrt{V}}\exp\bigg{(}-\frac{x^{2}}{2V}-2y^{2}\bigg{)},$$

$${\sigma }^{\prime }=\sum _{n=0}^{\infty }{\lambda }_n\hat{S}(r,\frac{\pi }{2})|n⟩⟨n|{\hat{S}}^{†}(r,\frac{\pi }{2}),$$

$$N_{\mathrm{DM2}}(V)=\max\bigg{(}\frac{1}{4\sqrt{V}}-\frac{1}{2},0\bigg{)}.$$

Full text

Nonclassicality is a fundamentally profound concept to identify quantum phenomena inaccessible from classical physics. It also provides a practically useful resource, e.g. entanglement, making possible a lot of applications in quantum information processing beyond classical counterparts ReviewCVQI ; Cerf ; ReviewGQI . A wide range of quantum systems, e.g. field amplitudes of light, collective spins of atomic ensembles, and motional modes of trapped ions, Bose-Einstein condensate and mechanical oscillators, can be employed for quantum information processing based on continuous variables (CVs) Cerf . It is of crucial importance to establish efficient and reliable criteria of nonclassicality for CV systems, desirably testable with less experimental resources, e.g. fewer measurement settings Richter2002 ; Richter2002' ; Mari2011 ; Park2015 ; ParkNha and with the capability of detecting a broad class of nonclassical states. In this paper, in view of the Glauber-Sudarshan P-function Glauber1963 ; Sudarshan1963 , those states that cannot be represented as a convex mixture of coherent states are referred to as nonclassical.

A standard method to obtain information on a CV quantum state is to measure marginal distributions along many different axes in phase space constituting quantum state tomography Lvovsky2009 . This tomographic reconstruction may reveal nonclassicality to some extent, e.g. negativity of Wigner function making only a subset of whole nonclassicality conditions. However it typically suffers from a legitimacy problem, i.e., the measured distributions do not yield a physical state when directly employed due to finite data and finite binning size D'Ariano1994 ; Lvovsky2009 . Much efforts was made to employ estimation methods finding a most probable quantum state closest to the obtained data Hradil1997 ; Banaszek1999 ; Rehacek2010 ; Teo2011 . There were also numerous studies to directly detect nonclassicality, e.g. an increasingly large number of hierarchical conditions Richter2002 requiring information on two or more marginal distributions or measurement of many higher-order moments Vogel ; Vogel' ; Bednorz . An exception would be the case of Gaussian states, with its nonclassical squeezing demonstrated by the variance of distribution along a squeezed axis.

Here we theoretically propose and experimentally demonstrate a simple, powerful, method to directly manifest nonclassicality by observing a single marginal distribution applicable to a wide range of nonclassical states. Our approach makes use of a phase-space map that transforms the marginal distribution (obtained from measurement) to a factorized Wigner distribution by multiplying the same distribution or the vacuum-state distribution along an orthogonal axis. We refer to those mathematical procedures as demarginalization maps (DMs), since a 1-dimensional marginal distribution is converted to a fictitious 2-dimensional Wigner function. The same method can be applied equally to the characteristic function as well as the Wigner function. We show that a classical state, i.e. a mixture of coherent states, must yield a physical state under our DMs. That is, the unphysicality emerging under DMs is a clear signature of nonclassicality. Remarkably, for all non-Gaussian states in finite dimensional space PN , our test works for an arbitrary single marginal distribution thus experimentally favorable. It also extends to non-Gaussian states in infinite dimension, particularly those without squeezing effect. We introduce a quantitative measure of nonclassicality using our DMs, which provides a lower bound of entanglement potential Asboth2005 —an entanglement measure under a beam-splitter setting versatile for CV entanglement generation Kim2002 ; Wang ; Asboth2005 ; Tahira2009 . Along this way, our method makes a rigorous connection between single-mode nonclassicality and negative partial transpose (NPT) entanglement Peres ; Horodecki ; Duan ; Simon ; MM ; MM' ; Walborn ; Walborn' ; Nha2008 ; Nha2012 , which bears on entanglement distillation EntanglementDistillation and nonlocality PeresConjecture ; PC1 ; PC2 ; PC3 ; PC4 ; PC5 ; PC6 ; PC7 .

As the measurement of a marginal distribution is highly efficient in various quantum systems, e.g. homodyne detection in quantum optics, our proposed approach can provide a practically useful and reliable tool in a wide range of investigations for CV quantum physics. We here experimentally illustrate the power of our approach by manifesting nonclassicality of motional states in a trapped-ion system. Specifically, we confirm the nonclassicality regardless of measured quadrature axis by introducing a simple faithful test using only a subset of data points, not requiring data manipulation under numerical methods unlike the case of state reconstruction. We also extend our approach combined with phase-randomization in order to obtain a criterion on genuine non-Gaussianity.

I Demarginalization maps and nonclassicality measure

I.1 Nonclassicality test via demarginalization maps

We first introduce our main tools, i.e. demarginalization maps (DMs),

[TABLE]

where $(x,y)^{T}=\mathcal{R}(\theta)(q,p)^{T}$ is a pair of orthogonal quadratures rotated from position $q$ and momentum $p$ with $\mathcal{R}(\theta)=\begin{pmatrix}\cos\theta&\sin\theta\\ -\sin\theta&\cos\theta\end{pmatrix}$. $M_{\rho}(x)=\int dyW_{\rho}(x,y)$ is a marginal distribution of the Wigner function $W_{\rho}(q,p)=\frac{2}{\pi}\mathrm{tr}[\rho\hat{D}(\alpha)(-1)^{\hat{n}}\hat{D}^{{\dagger}}(\alpha)]$, where $\hat{D}(\alpha)=e^{\alpha\hat{a}^{{\dagger}}-\alpha^{*}\hat{a}}$ is a displacement operator with $\alpha=q+ip$ Barnett ; Scully .

Our DM methods proceed as follows. Given a state with its Wigner function $W_{\rho}(q,p)$, we measure a marginal distribution $M_{\rho}(x)$ along a certain axis, $x=q\cos\theta+p\sin\theta$. We then construct a fictitious, factorized, Wigner function $W_{\rho}^{\rm DM}(x,y)$ either by replicating the obtained distribution as $M_{\rho}(x)M_{\rho}(y)$ (DM1) or by multiplying the marginal distribution of a vacuum state as $M_{\rho}(x)M_{|0\rangle\langle 0|}(y)$ (DM2), with $M_{|0\rangle\langle 0|}(y)=\sqrt{\frac{2}{\pi}}e^{-2y^{2}}$ (Fig. 1.). We test whether $W_{\rho}^{\rm DM}(x,y)$ is a legitimate Wigner function to represent a physical state.

Nonclassicality criteria—The constructed functions in Eqs. 1 and 2 are both in factorized forms, so judging their legitimacy is related to the problem what quantum states can possess a factorized Wigner function Note1 . Every coherent state $|\beta\rangle$ has a factorized Wigner function against all pairs of orthogonal quadratures, $W_{|\beta\rangle\langle\beta|}(x,y)=\frac{2}{\pi}e^{-2(x-\beta_{x})^{2}}e^{-2(y-\beta_{y})^{2}}$ Park2015 . Owing to this factorizability, the maps $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ transform a classical state into another classical one. A mixture of coherent states has a Wigner function

[TABLE]

with the probability density $P(\beta_{x},\beta_{y})$ for a coherent state $|\beta\rangle$ ($\beta=\beta_{x}+i\beta_{y}$). Applying each DM leads to

[TABLE]

$(j=1,2)$, where $Q_{1}(\beta_{x})=\int d\beta_{y}P(\beta_{x},\beta_{y})$ and $Q_{2}(\beta_{x})=\delta(\beta_{x})$ are nonnegative. The resulting distributions in Eq. 4 also represent a certain mixture of coherent states, hence a physical state. Therefore, if an unphysical Wigner function emerges under our DMs, the input state must be nonclassical.

Gaussian states—Let us first consider a Gaussian state $\sigma$ that has a squeezed quadrature $\hat{x}$ with $V_{x}\equiv\Delta^{2}\hat{x}<\frac{1}{4}$. Taking the squeezed marginal $M_{\sigma}(x)=\frac{1}{\sqrt{2\pi V_{x}}}e^{-\frac{(x-\langle\hat{x}\rangle)^{2}}{2V_{x}}}$ yields

[TABLE]

both of which violate the uncertainty relation $\Delta\hat{x}\Delta\hat{y}\geq\frac{1}{4}$. Thus, the squeezed state turns into an unphysical state under our DMs. This method, of course, succeeds only when the observed marginal distribution is along a squeezed axis that generally extends to a finite range of angles, if not the whole angles supple . We can further make the test successful regardless of quadrature axis by introducing a random phase rotation on a Gaussian state supple . Note that a mixture of phase-rotations, which transforms a Gaussian to a non-Gaussian state, does not create nonclassicality, so the nonclassicality detected after phase rotations is attributed to that of the original state.

Non-Gaussian states—More importantly, we now address non-Gaussian states. Every finite-dimensional state (FDS) in Fock basis, i.e. $\rho=\sum_{j,k=0}^{N}\rho_{jk}|j\rangle\langle k|$ is nonclassical, since all coherent states (except vacuum), and their mixtures as well, have an extension to infinite Fock states. It is nontrivial to demonstrate the nonclassicality of FDS when one has access to limited information, e.g., a noisy state $f|0\rangle\langle 0|+(1-f)|1\rangle\langle 1|$ for $f\geq\frac{1}{2}$ has no simple signatures of nonclassicality like squeezing and negativity of Wigner function. We prove that our DMs are able to detect all non-Gaussian states in finite dimension of any size, with details in Sec. S4 of supple . The essence of our proof is that there always exists a submatrix of the density operator corresponding to DMs, which is not positive-definite. Remarkably, this non-positivity emerges for a marginal distribution along an arbitrary direction, which means that the nonclassicality of FDS is confirmed regardless of the quadrature axis measured, just like the phase-randomized Gaussian states introduced in supple . This makes our DM test experimentally favorable, while the degree of negativity may well depend on the quadrature axis except rotationally-symmetric states. Our criteria can further be extended to non-Gaussian states in infinite dimension, particularly those without squeezing effect supple .

As an illustration, we show the case of a FDS $|\Psi\rangle=\frac{1}{\sqrt{2}}\left(|0\rangle+|2\rangle\right)$, whose original Wigner function and matrix elements are displayed in Fig. 1 (a). Our DM methods yield matrix elements as shown in Fig. 1 (b) and (c). The non-positivity of the density operator is then demonstrated by, e.g., $\langle 0|\rho|0\rangle\langle 8|\rho|8\rangle-|\langle 0|\rho|8\rangle|^{2}<0$ under DM1 and $\langle 0|\rho|0\rangle\langle 4|\rho|4\rangle-|\langle 0|\rho|4\rangle|^{2}<0$ under DM2, respectively.

I.2 Nonclassicality measure and entanglement potential

We may define a measure of nonclassicality using our DMs as

[TABLE]

where $||\cdot||_{1}$ is a trace norm and $\rho_{\mathrm{DM}}^{\theta}$ a density matrix under DM using a marginal distribution at angle $\theta$. Our DM negativity possesses the following properties appropriate as a nonclassicality measure, with details in supple . (i) $\mathcal{N}_{\mathrm{DM}}=0$ for a classical state, (ii) convex, i.e. non-increasing via mixing states, $\mathcal{N}_{\mathrm{DM}}(\sum_{j}p_{j}\rho_{j})\leq\sum_{j}p_{j}\mathcal{N}_{\mathrm{DM}}(\rho_{j})$, and (iii) invariant under a classicality-preserving unitary operation, $\mathcal{N}_{\mathrm{DM}}(\hat{U}_{c}\rho\hat{U}_{c}^{{\dagger}})=\mathcal{N}_{\mathrm{DM}}(\rho)$, where $\hat{U}_{c}$ refers to displacement or phase rotation. Combining (ii) and (iii), we also deduce the property that (iv) $\mathcal{N}_{\mathrm{DM}}$ does not increase under generic classicality preserving operations (mixture of unitary operations).

Our nonclassicality measure also makes a significant connection to entanglement potential as follows. A prototypical scheme to generate a CV entangled state is to inject a single-mode nonclassical state into a beam splitter (BS) Kim2002 ; Wang ; Asboth2005 ; Tahira2009 . It is important to know the property of those entangled states under PT, which bears on the distillibility of the output to achieve higher entanglement. Our formalism makes a connection between nonclassicality of single-mode resources and NPT of output entangled states. The effect of PT in phase space is to change the sign of momentum, $W_{\rho_{12}}(q_{1},p_{1},q_{2},p_{2})\rightarrow W_{\rho_{12}}(q_{1},p_{1},q_{2},-p_{2})$. If the resulting Wigner function is unphysical, the state $\rho_{12}$ is NPT. We first show that all nonclassical states detected under our DMs can generate NPT entanglement via a BS setting.

We inject a single-mode state $\rho$ and its rotated version $\overline{\rho}=e^{i\frac{\pi}{2}\hat{n}}\rho e^{-i\frac{\pi}{2}\hat{n}}$ into a 50:50 beam splitter (BS), described as

[TABLE]

Applying PT on mode 2 and injecting the state again into a 50:50 BS, we have

[TABLE]

Integrating over $q_{2}$ and $p_{2}$, the marginal Wigner function for mode 1 is given by $M_{\rho}(q_{1})M_{\rho}(p_{1})$, which is identical to DM1 of the state $\rho$ in Eq. 1. The other DM2 in Eq. 2 emerges when replacing the second input state $\overline{\rho}$ by a vacuum $\overline{\rho}=|0\rangle\langle 0|$. Therefore, if the original state $\rho$ is nonclassical under our DMs, the output entangled state via the BS scheme must be NPT.

In Ref. Asboth2005 , single-mode nonclassicality is characterized by entanglement potential via a BS setting, where a vacuum is used as an ancillary input to BS to generate entanglement. We may take negativity, instead of logarithmic negativity in Asboth2005 , as a measure of entanglement potential, i.e.,

[TABLE]

where $\hat{U}_{\mathrm{BS}}$ and $[\cdot]^{\mathrm{PT}}$ represent 50:50 beam-splitter operation and partial transpose on the mode 2, respectively. We then prove in supple that our DM2 measure provides a lower bound for the entanglement potential as

[TABLE]

Thus the nonclassicality measured under our framework indicates the degree of entanglement achievable via BS setting.

II Experiment

We experimentally illustrate the power of our approach by detecting nonclassicality of several motional states of a trapped 171Yb*+* ion. For the manipulation of motional state, the single phonon-mode $\hat{a}$ along $\rm X$-direction in 3-dimensional harmonic potential with trap frequencies $\left(\omega_{\rm X},\omega_{\rm Y},\omega_{\rm Z}\right)=2\pi~{}(2.8,3.2,0.6)$ MHz is coupled to two internal levels of the $S_{1/2}$ ground state manifold, $|F=1,m_{F}=0\rangle\equiv|\uparrow\rangle$ and $|F=0,m_{F}=0\rangle\equiv|\downarrow\rangle$ with transition frequency $\omega_{\rm HF}=(2\pi)~{}12.642821$ GHz. We implement the anti-Jaynes-Cumming interaction $H_{\rm aJC}=\frac{\eta\Omega}{2}\hat{a}^{\dagger}\hat{\sigma}_{+}+{\rm h.c.}$ and the Jaynes-Cumming interaction $H_{\rm JC}=\frac{\eta\Omega}{2}\hat{a}\hat{\sigma}_{+}+{\rm h.c.}$ with $\sigma_{+}=|\uparrow\rangle\langle\downarrow|$. $H_{\rm aJC}$ is realized by two counter-propagating laser beams with beat frequency near $\omega_{\rm HF}+\omega_{\rm X}$ and $H_{\rm JC}$ with frequency near $\omega_{\rm HF}-\omega_{\rm X}$ Park2015 . $\eta=\Delta k\sqrt{\hbar/2M\omega_{\rm X}}$ is the Lamb-Dicke parameter, $\Omega$ the Rabi frequency of internal transition, $\Delta k$ the net wave-vector of the Raman laser beams and $M$ the ion mass.

For our test, we generate the Fock states $|n=1\rangle$ and $|n=2\rangle$, together with the ground state $|n=0\rangle$. First, we prepare the ground state by applying the standard Doppler cooling and the Raman sideband cooling. Then we produce the Fock states by a successive application of $\pi$-pulse of $H_{\rm aJC}$ transferring the state $|\downarrow,n\rangle$ to $|\uparrow,n+1\rangle$, and the $\pi$-pulse for internal state transition $|\uparrow,n+1\rangle$ to $|\downarrow,n+1\rangle$. We also generate a superposition state $\frac{1}{\sqrt{2}}(|0\rangle+|2\rangle)$ by applying the $\pi/2$ pulse of $H_{\rm aJC}$ and then the $\pi$ pulse of $H_{\rm JC}$.

Nonclassicality test— We measure a characteristic function $C_{\rho}(k_{\theta})\equiv\langle e^{-2ik\hat{x}_{\theta}}\rangle$ with $\hat{x}_{\theta}=\hat{x}\cos\theta+\hat{p}\sin\theta$, by first making the evolution $\hat{U}=e^{-ik\hat{x}_{\theta}\hat{\sigma}_{x}}$ (simultaneously applying $H_{\rm aJC}$ and $H_{\rm JC}$ with proper phases) and then measuring internal state $\hat{\sigma}_{z}=|\uparrow\rangle\langle\uparrow|-|\downarrow\rangle\langle\downarrow|$ at times $t_{i}$ ($k=\eta\Omega t$) Roos ; Roos' . Using $\hat{U}^{\dagger}\hat{\sigma}_{z}\hat{U}=\cos(2k\hat{x}_{\theta})\hat{\sigma}_{z}+\sin(2k\hat{x}_{\theta})\hat{\sigma}_{y}$, we obtain $\langle\cos(2k\hat{x}_{\theta})\rangle$ and $\langle\sin(2k\hat{x}_{\theta})\rangle$, with the internal state initially prepared in the eigenstates $|+\rangle_{z}$ and $|+\rangle_{y}$ of $\sigma_{z}$ and $\sigma_{y}$, respectively. The Fourier transform of $C_{\rho}(k_{\theta})$ gives the marginal distribution of $\hat{x}_{\theta}$ Roos ; Roos' . In contrast, we directly use it without the Fourier transform, for which our DMs work equally well as for the Wigner function. We test $C^{\rm DM1}\equiv C_{\rho}(k_{x})C_{\rho}(k_{y})$ or $C^{\rm DM2}\equiv C_{\rho}(k_{x})C_{|0\rangle\langle 0|}(k_{y})=C_{\rho}(k_{x})e^{-\frac{1}{2}\lambda_{k_{y}}^{2}}$, with its density operator $\rho=\frac{1}{\pi}\int dk_{x}dk_{y}C(k_{x},k_{y}){\hat{D}}^{{\dagger}}(k_{x},k_{y})$ unphysical for a nonclassical state.

To set a benchmark (noise level) for classical states, we prepared the motional ground state $|n=0\rangle$ and obtained its marginal distributions along six axes with 1000 repetitions for each time $t_{i}$. It yielded the negativity $\mathcal{N}_{\mathrm{DM2}}=0.019\pm 0.02$ as represented by gray shading in Fig. 2 supple . On the other hand, the Fock states $|n=1\rangle$ and $|n=2\rangle$ clearly manifest nonclassicality for each marginal distribution taken at three different angles $\theta$ in Fig. 2 (a), at much higher negativity with error bars considering a finite data 1000. To further show that our method works regardless of measured axis, we also tested a superposition state $\frac{1}{\sqrt{2}}(|0\rangle+|2\rangle)$ not rotationally symmetric in phase space. As shown in Fig. 2 (b), its nonclassicality is well demonstrated for all measured angles individually while the degree of negativity varies with the measured axis.

Compared to our DM, one might look into nonclassicality directly via deconvolution, i.e. examining whether a marginal distribution $P(x)$ can be written as a sum of coherent-state distributions as $P(x)=\sqrt{\frac{2}{\pi}}\int d\tilde{x}\tilde{P}(\tilde{x})e^{-2(x-\tilde{x})^{2}}$, where $\tilde{P}(\tilde{x})$ must be positive-definite for classical states. $\tilde{P}(\tilde{x})$ is nothing but the marginal of Glauber-Sudarshan $P$-function, thus typically ill-behaved. One can test the positivity of $\tilde{P}(\tilde{x})$ alternatively using an $n\times n$ moment-matrix with elements $M_{ij}\equiv\langle\tilde{x}^{i+j}\rangle$ $(i,j=0,\dots,n-1)$Agarwal1993 . Figs. 2 (c) and (d) show the results under deconvolution using the same experimental data as we employed in Figs. 2 (a) and (b). To confirm nonclassicality, the degree of negativity must be large enough to beat that of the vacuum state including the statistical errors. Although those states produce negativity under deconvolution, their statistical errors are substantially overlapped with that of the vacuum state providing a much weaker evidence of nonclassicality than our DM. Full details are given in supple .

Instead of employing an entire characteristic function, we can also test our criterion by examining a subset of data using the Kastler-Loupias-Miracle-Sole (KLM) condition KLM ; KLM' ; KLM'' ; LegitimateW . This simple test provides a clear evidence of nonclassicality against experimental imperfections e.g. coarse-graining and finite data acquisition in other experimental platforms as well. The KLM condition states that the characteristic function $C_{\rho}(\xi)\equiv\mathrm{tr}[\rho\hat{D}(\xi)]$ for a legitimate quantum state must yield a $n\times n$ positive matrix $\mathcal{M}>0$ with matrix elements

[TABLE]

for an arbitrary set of complex variables $\{\xi_{1},\xi_{2},...,\xi_{n}\}$. In our case, we test the positivity of a matrix ($n=9$) constructed using $3\times 3$ points of rectangular lattice of size $d$ for the characteristic function under DM2 supple . As shown in Fig. 3 (a), the ground state $|0\rangle$ shows nonnegativity (thus the mixture of coherent states as well due to convexity of our method) for all values of $d$, whereas a nonclassical state $|\Psi\rangle=\frac{1}{\sqrt{2}}\left(|0\rangle+|2\rangle\right)$ manifests negativity in a certain range of $d$ confirming nonclassicality for each measured distribution at $\theta=\frac{\pi}{3},\frac{\pi}{2}$, and $\frac{2\pi}{3}$ (red solid curves) in Fig. 3 (b), (c) and (d), respectively. Furthermore, note that a mixture of the vacuum and the nonclassical state, $f|0\rangle\langle 0|+(1-f)|\Psi\rangle\langle\Psi|$ possesses a positive-definite Wigner function for $f\geq 0.66$, so even a full tomography may not directly show its nonclassicality via negativity. In contrast, our simple method manifests nonclassicality for $f=0.66$, as shown by blue dashed curves in Fig. 3 (b), (c) and (d). For Fock states, we consider the matrix test using $5\times 5$ lattice-points, which confirms negativity at the mixing $f=0.5$ with vacuum giving a non-negative Wigner function for both states $|1\rangle$ and $|2\rangle$ in Fig. 3 (e) and (f).

Genuine non-Gaussianity—We further extend our approach combined with phase-randomization to derive a criterion on genuine non-Gaussianity. Notably there exist quantum tasks that cannot be achieved by Gaussian resources, e.g. universal quantum computation Lloyd , CV nonlocality test NhaCV ; NhaCV' , entanglement distillation Eisert ; Fiurasek ; Giedke and error correction Niset . It is a topic of growing interest to detect genuine non-Gaussianity that cannot be addressed by a mixture of Gaussian states. Previous approaches all address particle nature like the photon-number distribution Filip2014 ; Jezek ; Straka2014 and the number parity in phase space Park2015 ; Genoni2013 ; Hughes2014 for this purpose. Here we propose a method to examine a genuinely CV characteristics of marginal distributions. Our criterion can be particularly useful to test a class of non-Gaussian states diagonal in the Fock basis, $\rho=\sum p_{n}|n\rangle\langle n|$, thus rotationally symmetric in phase space. For this class, one may detect nonclassicality using photon-number moments Vogel' , which can be experimentally addressed efficiently by phase-averaged quadrature measurements Munroe1995 ; Banaszek1997 . Lvovsky and Shapiro experimentally demonstrated the nonclassicality of a noisy single-photon state $f|0\rangle\langle 0|+(1-f)|1\rangle\langle 1|$ for an arbitrary $f$ Shapiro using the Vogel criterion Vogel'' . In contrast, we look into the genuine non-Gaussianity of non-Gaussian states as follows.

For a Gaussian state $\rho_{\rm G}$, the phase-randomization gives $\sigma\equiv\frac{1}{N}\sum_{k=0}^{N-1}e^{-i\theta_{k}{\hat{n}}}\rho_{\rm G}e^{i\theta_{k}{\hat{n}}}$ with $\theta_{k}\equiv\frac{k}{N}\pi$. As the number $N$ of phase rotation grows, the DM negativity of Gaussian states decreases. With $N\rightarrow\infty$ (full phase randomization), we obtain the Gaussian bound $\mathcal{B}_{\rm G}\approx 0.0887$ supple . Thus, if a state manifests a larger DM negativity as $\mathcal{N}>\mathcal{B}_{\rm G}$, it confirms genuine non-Gaussianity. We plot the Gaussian bounds for finite rotations $N=6$ and $N=12$ with $\mathcal{B}_{\rm G}\approx 0.0887$ against energy $n$ in Fig. 4. Our data for the state $|2\rangle$, which shows negativity insensitive to measured angles in Fig. 2, indicates genuine non-Gaussianity for the mixed states $f|0\rangle\langle 0|+(1-f)|2\rangle\langle 2|$ with $f=1-\frac{n}{2}$. For example, the $N=12$ case (brown dot-dashed) as well as the full phase randomization (black dashed horizontal) confirms quantum non-Gaussianity at $f=\frac{1}{2}$ corresponding to a positive Wigner-function.

III Conclusion and Remarks

Measuring marginal distributions along different axes in phase space forms a basis of quantum-state tomography with a wide range of applications. A marginal distribution is readily obtained in many different experimental platforms, e.g. by an efficient homodyne detection in quantum optical systems Lvovsky2009 ; Lvovsky2001 ; Ourjoumtsev2006 ; Huisman2009 ; Cooper2013 and by other quadrature measurements in trapped ion Roos ; Roos' ; ion , atomic ensembles as , optomechanics opto ; opto' , and circuit QED systems CQED ; CQED' . We here demonstrated that only a single marginal distribution can manifest nonclassicality by employing our demarginalization maps (DMs). Our DM methods are powerful to detect a wide range of nonclassical states, particularly non-Gaussian states. They provide a practical merit with less experimental efforts and make a stronger test of nonclassicality by analyzing data without numerical manipulation unlike state tomography.

Remarkably, nonclassicality can be demonstrated regardless of measured quadrature axis for all FDS states, which was also experimentally confirmed using a trapped ion system. We clearly showed that the proposed method provides a reliable nonclassicality test directly using a finite number of data, which can be further extended to other CV systems. In addition to the KLM test used here, we can manifest nonclassicality by looking into single marginal distributions under other forms, e.g. functional Nha2012 and entropic EUR ; EUR' inequalities. We also extended our approach to introduce a criterion on genuine non-Gaussianity employing marginal distributions combined with phase-randomization process. Our nonclassicality and non-Gaussianity tests were experimentally shown to successfully detect non-Gaussian states even with positive-definite Wigner functions whose nonclassicality is thus not immediately evident by the tomographic construction of Wigner function. As a note, for those nonclassical states with positive Wigner functions, one may use generalized quasi probability distributions like a filtered P-function Agarwal1970-1 ; Agarwal1970-2 ; Agarwal1970-3 . For example, the experiment in Kiesel2011 introduced a nonclassicality filter to construct a generalized P-function that yields a regularized distribution with negativity as a signature of nonclassicality for the case of photon-added thermal states. On the other hand, our DM method does not require a tomographic construction and provides a faithful test reliable against experimental imperfections like finite data and coarse graining.

Moreover, we established the connection between single-mode nonclassicality and NPT entanglement via BS setting—a prototypical model of producing CV entanglement. The negativity under our DM framework provides a quantitative measure of useful resource by identifying the minimum level of entanglement achievable in Eq. 10 NN . Nonclassicality and non-Gaussianity are important resources making a lot of quantum tasks possible far beyond classical counterparts. We thus hope our proposed method could provide a valuable experimental tool and a novel fundamental insight for future studies of CV quantum physics by critically addressing them.

Acknowledgements— M.S.Z. and H.N. were supported by National Priorities Research Program Grant 8-751-1-157 from the Qatar National Research Fund and K.K. by the National Key Research and Development Program of China under Grants 2016YFA0301900 and 2016YFA0301901 and the National Natural Science Foundation of China under Grants 11374178 and 11574002.

S1. Experiment: Analysis of fictitious characteristic function

We here explain how our experimental data are analyzed to yield a quantitative measure of DM negativity for each state in Fig. 2 and the lowest eigenvalue of KLM test in Fig. 3 of main text.

• •

DM2 negativity

First, to determine the DM negativity (noise level) of classical states, we prepared the motional ground state $|0\rangle$ and observed the marginal distributions $C(k_{x_{\theta}})\equiv\langle e^{-2ik\hat{x}_{\theta}}\rangle$, where $\hat{x}_{\theta}=\hat{x}\cos\theta+\hat{p}\sin\theta$, over six different angles $\theta$ with 1000 data at each time $t_{i}$ with $k=\eta\Omega t_{i}$. The resulting distribution including all data with error bars at each $k_{x}$ is drawn in Fig. 5 (a). A fictitious characteristic function is then constructed under DM2 as $C(k_{x},k_{y})=C(k_{x})e^{-\frac{1}{2}k_{y}^{2}}$, of which contour plot is given in Fig. 5 (c). We calculate the negativity of the corresponding density operator $\rho=\frac{1}{\pi}\int dk_{x}dk_{y}C(k_{x},k_{y}){\hat{D}}^{{\dagger}}(k_{x},k_{y})$ in the number-state basis. In Fig. 5 (e), we show the matrix elements $\rho_{mn}^{\rm DM2}$ of the density operator. $\frac{||\rho^{\mathrm{DM2}}||_{1}-1}{2}$ gives the result $\mathcal{N}_{\mathrm{DM2}}=0.019\pm 0.02$ in the main text.

On the other hand, for the case of a nonclassical state, e.g. $\frac{1}{\sqrt{2}}(|0\rangle+|2\rangle)$, the same procedures are taken for each measured axis, of which result, e.g. at $\theta=\frac{\pi}{2}$, is displayed in Fig. 5 (b,d,f).

• •

KLM test

For the purpose of KLM test mentioned in the main text, we choose a matrix using $n\times n$ lattice points in the space of characteristic function. For example, we illustrate the case of $3\times 3$ square lattice (black dots) in Figs.  5 (c) and (d). We may look at matrices by changing the lattice size $d$ according to data availability, all of which must give nonnegative eigenvalues for classical states. In other words, if there exists a negative $\lambda_{\rm min}$ for any $d$, it confirms nonclassicality.

S2. Comparison of DM and deconvolution

Given a marginal distribution, one can also look into nonclassicality via a deconvolution method instead of our DM. That is, one examines whether the obtained probability distribution $P(x)$ can be written as a convex sum of coherent-state distributions as $P(x)=\sqrt{\frac{2}{\pi}}\int d\tilde{x}\tilde{P}(\tilde{x})e^{-2(x-\tilde{x})^{2}}$, where $\tilde{P}(\tilde{x})$ must be positive-definite for classical states. In other words, if one finds negativity of $\tilde{P}(\tilde{x})$ under a Gaussian deconvolution of $P(x)$, it is a signature of nonclassicality. For their corresponding characteristic functions $C(k)=\int dxP(x)e^{-2ikx}$ and ${\tilde{C}}(k)=\int d\tilde{x}\tilde{P}(\tilde{x})e^{-2ik\tilde{x}}$, respectively, the deconvolution gives the relation ${\tilde{C}}(k)=C(k)e^{\frac{1}{2}k^{2}}$.

First, note that $\tilde{P}(\tilde{x})$ is nothing but the marginal distribution of Glauber-Sudarshan $P$-function, thus it is usually ill-behaved, e.g. delta-function for a coherent state and singular for Fock states. However, one can test the positivity of $\tilde{P}(\tilde{x})$ alternatively based on a moment test. For a classical state, a $n\times n$ matrix whose elements are given by $M_{ij}\equiv\langle\tilde{x}^{i+j}\rangle$ $(i,j=0,\dots,n-1)$ must be positive-definite at all levels of $n$. The moments $\langle\tilde{x}^{m}\rangle$ can be obtained via the moments $\langle x^{m}\rangle$ of the regular probability distribution $P(x)$ by way of $\langle\tilde{x}^{m}\rangle=2^{-3m/2}\langle H_{m}(\sqrt{2}x)\rangle$ with $H_{m}(x)$ the Hermite polynomial of order $m$. To derive it, let us use a method of operator ordering as follows. Due to the Baker–Campbell–Hausdorff formula, we have $e^{\lambda(\hat{a}+\hat{a}^{{\dagger}})}=e^{\lambda^{2}/2}e^{\lambda\hat{a}^{{\dagger}}}e^{\lambda\hat{a}}$, i.e. $:e^{\lambda(\hat{a}+\hat{a}^{{\dagger}})}:=e^{\lambda\hat{a}^{{\dagger}}}e^{\lambda\hat{a}}=e^{-\lambda^{2}/2}e^{\lambda(\hat{a}+\hat{a}^{{\dagger}})}$. We then obtain

[TABLE]

which yields the desired relation $:\hat{q}^{m}:=2^{-3m/2}H_{m}(\sqrt{2}\hat{q})$ with $\hat{q}=\frac{\hat{a}+\hat{a}^{{\dagger}}}{2}$. Note that we have used $H_{m}(x)=(-1)^{m}e^{x^{2}}\frac{d^{m}}{dx^{m}}e^{-x^{2}}$, $H_{m}(-x)=(-1)^{m}H_{m}(x)$ and $H_{m}(x+y)=\sum_{k=0}^{m}\binom{m}{k}H_{k}(x)(2y)^{m-k}$.

For a fair comparison with our DM, we here show the results of moment-matrix test under deconvolution for the states $|1\rangle$, $|2\rangle$ and $\frac{1}{\sqrt{2}}(|0\rangle+|2\rangle)$ together with the vacuum $|0\rangle$ using the same experimental data as we employed in Fig. 2 of the main text. To confirm nonclassicality, the degree of negativity must be large enough to beat that of the vacuum state including the statistical errors. As shown in Fig. 6, those nonclassical states produce negativity under deconvolution, however their statistical errors are substantially overlapped with that of the vacuum state making the interpretation weak (with the only exception at $\theta=\frac{\pi}{2}$ for the state $\frac{1}{\sqrt{2}}(|0\rangle+|2\rangle)$). This is clearly contrasted with the results under our DM in Fig. 2 of the main text.

One may try to optimize the matrix test by increasing the size $n$ of the matrix. This is theoretically valid, since the negativity not appearing in a low-dimension matrix can be found in a higher-dimension matrix as the latter encompasses the former. However, under practical situations with a finite number $N$ of data, the statistical error of the moment $\langle\tilde{x}^{m}\rangle$ given by $\Delta_{m}=\frac{2^{-3m/2}}{\sqrt{N}}\sqrt{\langle H_{m}^{2}(\sqrt{2}x)\rangle-\langle H_{m}(\sqrt{2}x)\rangle^{2}}$ exponentially increase with $m$. For instance, $\Delta_{m}=\sqrt{m!4^{-m}/N}$ for the vacuum state whereas it grows at a higher rate for typical nonclassical states as shown in Fig. 6 (g) and (h). It is thus harder with a larger $m$ to beat the reference level (vacuum) of nonclassicality with significant errors, and the matrix test becomes optimized at a moderate level of matrix dimension. As shown in Fig. 6, the test does not necessarily improve by increasing the matrix size. If one obtains much more data $N$, the effect of statistical error may be reduced, however, the fact that our DM method works well already with a low number of data is an evidence of superioity in manifesting nonclassicality. Note that our DM method employs well-behaved square-integrable functions by its construction unlike the deconvolution method. Furthermore, unlike the deconvolution method aiming only at detecting nonclassicality, our DM formalism also constitute a useful framework to connect nonclassicality and quantum entanglement in CV setting and to address a genuine non-Gaussianity of CV systems that has been of growing interest.

S3. Testing Gaussian states under DMs

As stated in the main text, a Gaussian nonclassical state, i.e. squeezed state, can be detected under our DMs if the measured marginal distribution is along a squeezed axis. For a given Gaussian state, we first determine a success range of measurement angles to manifest its nonclassicality. We later show that the nonclassicality can be demonstrated regardless of measurement axis by converting the Gaussian state to a non-Gaussian state under a finite number of phase rotations, which does not create nonclassicality.

A single-mode Gaussian state can generally be represented as a displaced squeezed thermal state,

[TABLE]

where $\hat{S}(r,\phi)=\exp[-\frac{r}{2}(e^{2i\phi}\hat{a}^{{\dagger}2}-e^{-2i\phi}\hat{a}^{2})]$ is a squeezing operator ($r$: squeezing strength, $\phi$: squeezing direction), and $\sigma_{\mathrm{th}}(\bar{n})=\sum_{n=0}^{\infty}\frac{\bar{n}^{n}}{(\bar{n}+1)^{n+1}}|n\rangle\langle n|$ is a thermal state with mean photon number $\bar{n}$. The displacement does not affect the physicality issue under our DM methods, so we may set $\alpha=0$ without loss of generality. Then a Wigner function with $\phi=0$ (squeezed along $q$-axis) is given by

[TABLE]

with $\mu=\frac{1}{1+2\bar{n}}$ the purity of the state. Its marginal distribution along the direction $\theta$ rotated from position $q$ is obtained as

[TABLE]

where $V_{\theta}=\frac{1}{4\mu}(\cosh 2r-\cos 2\theta\sinh 2r)$ is the variance. We identify the squeezing range $V_{\theta}<\frac{1}{4}$ to be $|\theta|<\frac{1}{2}\arccos\frac{\cosh 2r-\mu}{\sinh 2r}\equiv|\theta|_{\mathrm{sq}}$. We thus obtain a success probability $\eta_{\mathrm{sq}}$ that a randomly chosen angle $\theta\in[0,\frac{\pi}{2}]$ for a marginal distribution can detect nonclassicality as

[TABLE]

We plot $\eta_{\mathrm{sq}}$ as a function of the purtiy $\mu\equiv\mathrm{tr}\rho^{2}$ and squeezing strength $r$ in Fig. 7. In general, $\eta_{\mathrm{sq}}$ monotonically increases with the purity $\mu$. On the other hand, it decreases with $r$ as the range of squeezed quadratures decreases with the degree of squeezing. For instance, in the extreme squeezing $r\rightarrow\infty$, the range of squeezing direction becomes $|\theta|_{\mathrm{sq}}\rightarrow 0$.

• •

Detecting Gaussian states under a classicality-preserving operation

We now introduce a finite number of phase rotations, which does not create nonclassicality, and show that it enables us to detect the nonclassicality of a Gaussian state regardless of measurement axis. The phase-randomizing operation generally transforms a Gaussian state to a non-Gaussian state and it was employed, e.g., for the distillation of squeezing Franzen and the reconstruction of a Wigner function of rotationally-symmetric state (Fock states) Grangier . We here use our DM2 approach and give proof in two-steps: (i) DM2 detects every nonclassical marginal distribution. (ii) The marginal distribution obtained from a squeezed state under phase-randomization is always nonclassical.

(i) Connection between DM2 and the Hamburger moment problem

For a given state $\rho$, DM2 gives the Wigner function as

[TABLE]

which is equivalent to the Glauber-P representation as

[TABLE]

Here $\ket{x}$ represents a coherent state with real amplitude $x$ and the quasi-probability density $T_{\rho}(x)$ connects to the marginal distribution by $M_{\rho}(q)=\int_{-\infty}^{\infty}T_{\rho}(x)\sqrt{\frac{2}{\pi}}e^{-2(q-x)^{2}}$. (Note that $\mathcal{D}_{2}[\rho]$ is classical if and only if $T_{\rho}(x)$ is non-negative. We thus see that if the marginal distribution $M_{\rho}(q)$ cannot be expressed as a positive sum of normal distributions with $\sigma=\frac{1}{2}$ (vacuum fluctuation), the corresponding $\mathcal{D}_{2}[\rho]$ cannot represent a classical state.)

In addition, the density operator $\mathcal{D}_{2}[\rho]$ is unphysical if and only if there exists a pure state $\ket{\psi}=\sum_{n}c_{n}\ket{n}$ satisfying

[TABLE]

In measure theory, it is known that for the moments $a_{n}=\int_{-\infty}^{\infty}x^{n}d\rho(x)$, the measure $\rho$ is positive if and only if

[TABLE]

is satisfied for every non-negative integer $N$ and every set of complex variables $\{\beta_{0},...,\beta_{N}\}$ (Hamburger moment problem) ReedSimon . Comparing Eqs. (19) and (20), we see that $\mathcal{D}_{2}[\rho]$ is unphysical if and only if $T_{\rho}(x)$ fails to be positive. It therefore proves that $\mathcal{D}_{2}[\rho]$ becomes unphysical for every nonclassical marginal distribution.

(ii) Gaussian states under a finite number of phase rotations

For a general Gaussian state in Eq. (13), we apply $N$ (finite) phase rotations of angles $\{0,\frac{1}{N}\pi,...,\frac{N-1}{N}\pi\}$ with $N\geq\lceil\frac{1}{\eta_{\mathrm{sq}}}\rceil$, where $\eta_{\mathrm{sq}}$ is given in Eq. (16). An equal mixture of those rotations gives a marginal distribution

[TABLE]

with $v_{k}=\frac{1}{4\mu}[\cosh 2r-\cos 2(\theta+\frac{k}{N}\pi)\sinh 2r]$. Importantly, there always exists a set of $k$ satisfying $v_{k}<\frac{1}{4}$ (squeezing) regardless of $\theta$ (measured axis), whenever $N\geq\lceil\frac{1}{\eta_{\mathrm{sq}}}\rceil$ is satisfied. Below we show that the distribution in Eq. (21) cannot be represented as a mixture of Gaussian distributions all with vacuum noise. Together with the property (i), this proves that the resulting non-Gaussian state can be detected regardless of quadrature axis under DM2 method.

Suppose that $M_{\rho}(x)$ in Eq. (21) be written as a mixture of coherent-state distributions, i.e.

[TABLE]

then one encounters a contradiction. Let us take a Fourier transform, $\int dxe^{-ixy}M_{\rho}(x)$, and then multiply both sides by $e^{\frac{v}{2}y^{2}}$ with $v\equiv{\rm min}\{v_{k}\}<\frac{1}{4}$. It gives

[TABLE]

where the sum in the LHS excludes the terms with $v_{k}=v={\rm min}\{v_{k}\}$. It is readily seen that Eq. (23) cannot be satisfied at all $y$, e.g., LHS=$\frac{1}{N}$ and RHS=0 as $y\rightarrow\infty$.

• •

Gaussian bound under phase randomization

Let us now consider how the degree of nonclassicality (DM negativity) is bounded for all Gaussian states by a full phase-randomization. When the measured quadrature distribution is Gaussian with variance $V$, the Wigner function under DM2 is given by

[TABLE]

and the corresponding density matrix by

[TABLE]

where $r=-\frac{1}{2}\log 2\sqrt{V}$ and $\lambda_{n}=\frac{2}{2\sqrt{V}+1}(\frac{2\sqrt{V}-1}{2\sqrt{V}+1})^{n}$ represent the eigenvalues of $\sigma^{\prime}$. A direct calculation gives its DM negativity as

[TABLE]

Using the convexity of the DM negativity (proved in the next section), we deduce that the DM negativity of a squeezed state $\rho_{\rm G}$ under phase randomization, i.e. $\sigma\equiv\frac{1}{2\pi}\int_{0}^{2\pi}d\theta e^{-i\theta{\hat{n}}}\rho_{\rm G}e^{i\theta{\hat{n}}}$ must be bounded as

[TABLE]

where $\mathcal{F}(\phi,m)=\int_{\theta=0}^{\phi}(1-m\sin^{2}\theta)^{-1/2}d\theta$ represents the elliptic integral of the first kind, and $\theta_{c}=\frac{1}{2}\arccos(\tanh r)$ sets the boundary for squeezing $V_{\theta_{c}}=\frac{1}{4}$. While the DM negativity for a Gaussian state approaches infinity with squeezing, $\lim_{V\rightarrow 0}N_{\mathrm{DM2}}(V)=\infty$, the phase randomization reduces it to $\mathcal{B}_{\rm G}\approx 0.0887$, thus

[TABLE]

We plot the upper bound of Eq. (S3. Testing Gaussian states under DMs) (red solid curve) in Fig. 8, which takes maximum $\mathcal{B}_{\rm G}\approx 0.0887$. Due to the convexity of the DM negativity, we can further say that all Gaussian states and their mixtures also obey the same inequality. Therefore, if a state under a full phase-randomization violates the inequality as $\mathcal{N}_{\mathcal{D}[\rho]}>\mathcal{B}_{\rm G}$, it is a clear signature of genuine non-Gaussianity. For example, all Fock states achieve DM negativity above the Gaussian bound as shown in Fig. 9. Our method can detect genuine non-Gaussianity even for states having a non-negative Wigner function as shown in main text.

We can also numerically find the Gaussian bound for the case of finite $N$ phase-rotations, i.e. $\sigma^{\prime}\equiv\frac{1}{N}\sum_{k=0}^{N-1}e^{-i\theta_{k}{\hat{n}}}\rho_{\rm G}e^{i\theta_{k}{\hat{n}}}$ with $\theta_{k}\equiv\frac{k}{N}\pi$. In Fig. 8, the cases of $N=6$ and $N=12$ are plotted against energy $n$.

S4. Properties of DM negativity as a measure of nonclassicality

Let us first define our DM negativity for each measured distribution as

[TABLE]

where $||\cdot||_{1}$ represents a trace norm and $\rho_{\mathrm{DM}}^{\theta}$ the fictitious density matrix obtained by applying our DM method to the marginal distribution measured at angle $\theta$. It satisfies the following properties:

$\mathcal{N}_{\mathrm{DM}}^{\theta}(\rho)=0$ for classical states. 2. 2.

$\mathcal{N}_{\mathrm{DM}}^{\theta}(\sum_{j}p_{j}\rho_{j})\leq\sum_{j}p_{j}\mathcal{N}_{\mathrm{DM}}^{\theta}(\rho_{j})$. 3. 3.

$\mathcal{N}_{\mathrm{DM}}^{\theta}(\hat{D}(\alpha)\rho\hat{D}^{{\dagger}}(\alpha))=\mathcal{N}_{\mathrm{DM}}^{\theta}(\rho)$.

[P1] In the main text, we already showed that $\rho_{\mathrm{DM}}^{\theta}$ becomes a physical state when the initial state is classical, i.e. $||\rho_{\mathrm{DM}}^{\theta}||_{1}=1$ giving $\mathcal{N}_{\mathrm{DM}}^{\theta}=0$.

[P2] We can derive the convexity of DM negativity from the fact that the trace norm satisfies the homogeneity, i.e., $||cM||_{1}=|c|||M||_{1}$, and the triangle inequality, i.e., $||M+N||_{1}\leq||M||_{1}+||N||_{1}$ Vidal2002 ; Plenio2002 .

[TABLE]

[P3] A displacement operation only shifts the center of Wigner function and marginal distributions.

[TABLE]

DM thereby yields a displaced version $\overline{\rho}_{DM}^{\theta}=\hat{D}(\overline{\alpha})\rho_{DM}^{\theta}\hat{D}^{{\dagger}}(\overline{\alpha})$ for the displaced state $\overline{\rho}=\hat{D}(\alpha)\rho\hat{D}^{{\dagger}}(\alpha)$ with $\overline{\alpha}=\mathrm{Re}[\alpha e^{-i\theta}]$, which therefore does not affect the physicality of the output density operator. As the eigenspectrum of a Hermitian operator is invariant under a unitary operator, we conclude that $\mathcal{N}_{\mathrm{DM}}^{\theta}(\hat{D}(\alpha)\rho\hat{D}^{{\dagger}}(\alpha))=\mathcal{N}_{\mathrm{DM}}^{\theta}(\rho)$.

While a phase rotation is also a classicality preserving operation, i.e., $e^{i\hat{n}\phi}\ket{\alpha}=\ket{\alpha e^{i\phi}}$, it can increase the value of $\mathcal{N}_{\mathrm{DM}}^{\theta}[\rho]$ as the degree of DM negativity depends on the measured angle $\theta$. Thus, to make the measure invariant under phase rotations as well, we introduce an optimized version of DM negativity as

[TABLE]

The DM negativity then satisfies

$\mathcal{N}_{\mathrm{DM}}(\rho)=0$ for classical states. 2. 2.

$\mathcal{N}_{\mathrm{DM}}(\sum_{j}p_{j}\rho_{j})\leq\sum_{j}p_{j}\mathcal{N}_{\mathrm{DM}}(\rho_{j})$. 3. 3.

$\mathcal{N}_{\mathrm{DM}}(\hat{U}_{c}\rho\hat{U}_{c}^{{\dagger}})=\mathcal{N}_{\mathrm{DM}}(\rho)$ where $\hat{U}_{c}$ is a classicality preserving unitary operation, i.e., displacement and phase rotation.

[P1’] It is a direct consequence of [P1].

[P2’] We also obtain the convexity of optimized DM negativity from [P2].

[TABLE]

[P3’] The property [P3] implies that $\mathcal{N}_{\mathrm{DM}}(\hat{D}(\alpha)\rho\hat{D}^{{\dagger}}(\alpha))=\mathcal{N}_{\mathrm{DM}}(\rho)$. As a phase rotation only rotates the marginal distributions, i.e., $\mathcal{N}_{\mathrm{DM}}^{\theta}(e^{-i\hat{n}\phi}\rho e^{i\hat{n}\phi})=\mathcal{N}_{\mathrm{DM}}^{\theta+\phi}(\rho)$, the DM negativity optimized over the angles is not changed at all. Combining 2 and 3, we also deduce that

4. $\mathcal{N}_{\mathrm{DM}}$ does not increase under generic classicality preserving operations (mixture of unitary operations).

Let us here look into other nonclassicality measures in the literature considering the above properties. Note that nonclassical depth Lee1991 does not satisfy the convexity. As an example, a mixture of vacuum and a squeezed state, i.e., $\sigma=(1-p)|0\rangle\langle 0|+p|\psi\rangle\langle\psi|$, satisfies that $\tau(\sigma)=\tau(|\psi\rangle\langle\psi|)\neq 0$ for $0<p\leq 1$ and $\tau(\sigma)=0$ for $p=0$.

Gehrke et al. addressed a set of conditions for a proper nonclassicality measure in Gehrke2012 . They have only addressed two conditions: (1) the measure is zero only for classical state. (2) the measure is nonincreasing for any classical Kraus operator. They also proposed a degree of nonclassicality which quantifies the number of superposition between coherent states required for representing the given state. As an example, the degree of nonclassicality for a pure state in the form $\sum_{j=0}^{d}c_{j}\ket{\alpha_{j}}$ is $d$. Similar to the case of nonclassical depth, the degree of nonclassicalty does not satisfy the convexity. For example, a mixture of vacuum ($d=0$) and a pure state $\ket{\psi}$ with non-zero degree $k$, $\sigma=(1-p)|0\rangle\langle 0|+p|\psi\rangle\langle\psi|$, satisfies that $d(\sigma)=k$ for $0<p\leq 1$ and $d(\sigma)=0$ for $p=0$.

Connection to entanglement potential— With the entanglement potential defined in the main text as

[TABLE]

we here prove the relation

[TABLE]

To this aim, we slightly modify the procedures in the main text as follows. After the beam-splitting operation for entanglement generation, we perform local unitary operations $e^{-i\theta\hat{n}_{1}}\otimes e^{-i\theta\hat{n}_{2}}$, which has no effect on entanglement and simply rotates the output Wigner function in Eq. (7) of main text as

[TABLE]

where $q_{\pm}=\frac{1}{\sqrt{2}}(q_{1}\pm q_{2})$, $p_{\pm}=\frac{1}{\sqrt{2}}(p_{1}\pm p_{2})$, together with the fact that the Wigner function for vacuum is invariant under rotation. PT and BS operations change it to

[TABLE]

Integrating over $q_{2}$ and $p_{2}$, we obtain the marginal Wigner function for one mode as $M_{\rho}(x_{1})M_{|0\rangle\langle 0|}(y_{1})$, which represents the demarginalization of the single mode state $\rho$ in Eq. (2) along with the rotated axis $x_{1}=q_{1}\cos\theta+p_{1}\sin\theta$.

As a unitary operation preserves the eigenspectrum of an input state, we have

[TABLE]

where $\rho^{\prime}=[\hat{R}(\theta)\hat{U}_{\mathrm{BS}}(\rho_{1}\otimes|0\rangle\langle 0|_{2})\hat{U}_{\mathrm{BS}}^{{\dagger}}\hat{R}^{{\dagger}}(\theta)]^{\mathrm{PT}}$ and $\hat{R}(\theta)=e^{-i\theta\hat{n}_{1}}\otimes e^{-i\theta\hat{n}_{2}}$. Using the fact that trace norm is nonincreasing under partial trace Vidal2002 ; Plenio2002 , e.g., $||\rho_{AB}||_{1}\geq||\rho_{A}||_{1}$, we obtain

[TABLE]

which is satisfied for every rotation angle $\theta$. From the obtained inequality, we derive the relation in Eq. (33).

Furthermore, we can also show that two measures are identical for a Gaussian state, $\mathcal{N}_{\mathrm{DM2}}[\sigma]=\mathcal{P}_{\mathrm{ent}}[\sigma]$ by direct calculations.

S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.

• •

Derivation of Wigner function for FDS

The Wigner function of the FDS $\rho=\sum_{j,k=0}^{N}\rho_{jk}|j\rangle\langle k|$ is given by

[TABLE]

where $W_{|j\rangle\langle k|}(x,y)$ represents the operator $|j\rangle\langle k|$

[TABLE]

with $L_{n}^{(\alpha)}(z)$ a generalized Laguerre polynomial of order $n$, $\alpha_{\theta}=(x+iy)e^{-i\theta}$ and $W_{|j\rangle\langle k|}(x,y)=W_{|k\rangle\langle j|}^{*}(x,y)$ for $j<k$, as shown in Ref. UL . With $\alpha_{\theta=0}=q+ip$, we have a simple relation $W_{|j\rangle\langle k|}(x,y)=e^{i\theta(j-k)}W_{|j\rangle\langle k|}(q,p)$, thus we are only concerned with $W_{|j\rangle\langle k|}(q,p)$ letting $\theta=0$.

We first reexpress Eq. (39) using $\frac{d^{k}}{dx^{k}}L_{n}^{(\alpha)}(x)=(-1)^{k}L_{n-k}^{(\alpha+k)}(x)$,

[TABLE]

Employing the formulas in Ref. AS , i.e. $L_{n}^{(\alpha+\beta+1)}(x+y)=\sum_{j=0}^{n}L_{j}^{(\alpha)}(x)L_{n-j}^{(\beta)}(y)$ and $H_{2n}(x)=(-4)^{n}n!L_{n}^{(-1/2)}(x^{2})$, we obtain

[TABLE]

Combining Eqs. (S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.) and (S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.) with $\frac{d}{dx}H_{n}(x)=2nH_{n-1}(x)$ and $\frac{d^{n+1}}{dx^{n+1}}H_{n}(x)=0$, we obtain

[TABLE]

with $s=\max[0,2\ell-j-k]$ and $t=\min[j-k,2\ell]$, where we used $\frac{d^{n}}{dx^{n}}[f(x)g(x)]=\sum_{s=0}^{n}\frac{n!}{(n-s)!s!}\frac{d^{n-s}f}{dx^{n-s}}\frac{d^{s}g}{dx^{s}}$. Finally, we can recast Eq. (S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.) to

[TABLE]

where

[TABLE]

We also provide an alternative expression of $A_{|j\rangle\langle k|}(n)$ for even $j+k-n$ as

[TABLE]

This can be obtained by comparing two methods of deriving marginal distribution for $|j\rangle\langle k|$, that is, $M_{|j\rangle\langle k|}(q)=\int dpW_{|j\rangle\langle k|}(q,p)$ and $M_{|j\rangle\langle k|}(q)=\braket{q}{j}\braket{k}{q}$, where $\ket{q}$ is an eigenstate of the position operator, $\hat{q}\ket{q}=q\ket{q}$. Using the formula in Ref. Wang2008 , $H_{j}(x)H_{k}(x)=\sum_{n=0}^{j+k}a_{jk}(n)H_{n}(\sqrt{2}x)$ with

[TABLE]

and the expression in Ref. BR , $\braket{q}{n}=(\frac{2}{\pi})^{\frac{1}{4}}\sqrt{\frac{1}{2^{n}n!}}e^{-q^{2}}H_{n}(\sqrt{2}q)$, we obtain

[TABLE]

On the other hand, from Eq. (43), we obtain

[TABLE]

where we have used $\int_{-\infty}^{\infty}dxe^{-2x^{2}}H_{2m}(2x)=\sqrt{\frac{\pi}{2}}\frac{(2m)!}{m!}$ and $\int_{-\infty}^{\infty}dxe^{-2x^{2}}H_{2m+1}(2x)=0$ for an integer $m\geq 0$. Comparing Eqs. (S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.) and (S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.), we obtain Eq. (45).

• •

Unphysicality of all FDSs under DM1 and DM2 methods

Starting with a FDS state $\rho=\sum_{j,k=0}^{N}\rho_{jk}|j\rangle\langle k|$ with its highest excitation $N$, i.e. $\rho_{NN}\neq 0$, the density matrix elements for $W_{\rho}^{\mathrm{DM1}}(x,y)$ under DM1 are obtained by

[TABLE]

Using the orthogonality of Hermite polynomials, $\int dxe^{-4x^{2}}H_{a}(2x)H_{b}(2x)=\sqrt{\pi}2^{b-1}b!\delta_{a,b}$, we have

[TABLE]

Looking into Kronecker delta functions and $j+k+j^{\prime}+k^{\prime}\leq 4N$, we first have $\rho_{n_{1},n_{2}}^{\prime}=0$ for $n_{1}+n_{2}>4N$ that naturally imposes $\rho^{\prime}_{n,n}=0$ for $n>2N$. We now derive compact expressions for $\rho_{4N,0}^{\prime}$ and $\rho_{0,4N}^{\prime}$, which correspond to $j=k=j^{\prime}=k^{\prime}=N$ and $m=m^{\prime}=0$, as

[TABLE]

where we have used $A_{|N\rangle\langle N|}(2N)=\frac{1}{4^{N}N!}$ and $A_{|4N\rangle\langle 0|}(2N)=\frac{\sqrt{(4N)!}}{(-16)^{n}\{(2N)!\}^{2}}$ from Eqs. (S5. All FDS can be detected via a single marginal distribution regardless of measurement axis.) and (45). For the determinant of a sub-matrix, we find that

[TABLE]

which confirms the nonclassicality of the original FDS state $\rho=\sum_{j,k=0}^{N}\rho_{jk}|j\rangle\langle k|$.

Similarly, the density matrix elements for $W_{\rho}^{\mathrm{DM2}}(x,y)$ under DM2 are given by

[TABLE]

which manifests that $\rho_{n_{1},n_{2}}^{\prime\prime}=0$ for $n_{1}+n_{2}>2N$. In addition, we obtain compact expressions for $\rho_{2N,0}^{\prime\prime}$ and $\rho_{0,2N}^{\prime\prime}$, for which $j=k=N$ and $m=0$, as

[TABLE]

where $A_{|2N\rangle\langle 0|}(2N)=\frac{1}{4^{N}\sqrt{(2N)!}}$ from Eq. (45). This leads to

[TABLE]

again confirming the nonclassicality of the original state. Note that the above results under both of DM1 and DM2 hold regardless of $\theta$ (quadrature axis), which makes our criteria experimentally favorable. That is, the nonclassicality for an arbitrary FDS state can be verified by observing a single marginal distribution along any directions.

S6. Detection of non-Gaussian states in infinite dimension

We here consider the detection of nonclassicality for non-Gaussian states in infinite dimension, which has practical relevance to CV quantum information processing. In particular, we examine those non-Gaussian states without squeezing effect in order to demonstrate the merit of our formalism while employing a marginal distribution.

(i) photon-added coherent state $\hat{a}^{{\dagger}}\ket{\gamma}$— We note that a photon-added coherent state do not have squeezing effect for $|\gamma|<1$ Agarwal1991 , which restricts the detection of its nonclassicality under a variance test. In contrast, our criteria detect it for all $\gamma$. We can recast the state as a displaced FDS,

[TABLE]

where we have used $\hat{D}(\gamma)\hat{a}^{{\dagger}}\hat{D}^{{\dagger}}(\gamma)=\hat{a}^{{\dagger}}+\gamma^{*}$. As our DM criteria are able to detect every FDS and invariant under displacement operation, we can detect every photon added coherent state using our criteria. Remarkably, our DM methods successfully detect its nonclassicality regardless of quadrature axis as for the case of FDS.

(ii) photon-added thermal state $\rho=\hat{a}^{{\dagger}}\rho_{\mathrm{th}}(\bar{n})\hat{a}$— Its Wigner function is given by

[TABLE]

with its marginal distribution given by

[TABLE]

Applying DM method to the marginal distribution and obtaining relevant density matrix elements, we consider $\{\ket{0},\ket{4}\}$ subspace for DM1,

[TABLE]

and $\{\ket{0},\ket{2}\}$ subspace for DM2,

[TABLE]

then a negative eigenvalue can be found at low mean photon number $n\lesssim 0.40$ and $n\lesssim 0.45$, respectively. For a higher mean photon number, we need to investigate density matrix elements with higher Fock numbers.

(iii) Dephased odd cat states—$\rho_{\mathrm{odd}}\sim|\gamma\rangle\langle\gamma|+|-\gamma\rangle\langle-\gamma|-f(|\gamma\rangle\langle-\gamma|+|-\gamma\rangle\langle\gamma|)$ where $\ket{\gamma}$ is a coherent state with amplitude $\gamma$ and $1-f$ denotes the degree of dephasing. They show no squeezing for all $\gamma>0$ and $f>0$. As we obtain the density matrix elements for $\mathcal{D}_{j}[W_{\rho_{\mathrm{odd}}}](x,y)$ employing the marginal distribution for the rotated quadrature $q_{\theta}$,

[TABLE]

we find a non-positive submatrix for $\gamma>0$ and $f>0$. More explicitly, for the case of $\mathcal{D}_{1}[W_{\rho_{\mathrm{odd}}}]=W_{\rho_{\mathrm{odd}}^{\prime}}$, we find negativity in $\{\ket{0},\ket{4}\}$ subspace,

[TABLE]

particularly its negative determinant

[TABLE]

for all $\gamma>0$ and $f>0$ independent of the quadrature axis. Note that $9(e^{2\gamma^{2}}+f^{2})+12(e^{2\gamma^{2}}-f^{2})t^{2}+4(e^{\gamma^{2}}-f)^{2}(2t^{2}-1)$ is positive for every $t$.

For the case of $\mathcal{D}_{2}[W_{\rho_{\mathrm{odd}}}]=W_{\rho_{\mathrm{odd}}^{\prime\prime}}$, we obtain density matrix elements in $\{\ket{0},\ket{2}\}$ subspace as

[TABLE]

with its negative determinant

[TABLE]

for all $\gamma>0$ and $f>0$ independent of the quadrature axis.