Measuring non-Gaussianity with Correlation

TL;DR

This paper introduces a correlation-based framework for quantifying quantum non-Gaussianity, linking it to operational measures and providing an experimentally feasible protocol for estimation, advancing both theory and practice.

Contribution

It presents a novel correlation-based approach to measure quantum non-Gaussianity, including an experimental protocol and theoretical bounds, unifying theory and practical estimation.

Findings

Correlation generation characterizes non-Gaussianity.

Operational measures are monotonic under Gaussian channels.

Sample-efficient protocol estimates non-Gaussianity in experiments.

Abstract

Quantum non-Gaussianity is a key resource for quantum advantage in continuous-variable systems. We introduce a general framework to quantify non-Gaussianity based on correlation generation: two copies of a state become correlated at a $50{:}50$ beam splitter if and only if the state is non-Gaussian, with correlations reducing to entanglement in the pure-state case. This connection enables operational measures of non-Gaussianity, defined through correlation quantifiers such as R\'enyi-$\alpha$ entropy for pure states and R\'enyi-$\alpha$ mutual information for mixed states. We prove that all such measures are monotonic under Gaussian channels. Building on this framework, we propose a sample-efficient experimental protocol to estimate non-Gaussianity using standard optical components, even in the state agnostic setting. Finally, we establish a lower bound on the sample complexity of estimating Wigner negativity, allowing a direct comparison with our protocol. Our results provide both a unifying theoretical framework for non-Gaussianity and a practical route toward its experimental quantification.

Figures

Equations

$$\text{quantity}[\hat{\mathbf{r}},{\hat{\mathbf{r}}}^T]=\mathrm{i}\Omega$$

$${\chi }_{{\rho }_G}(\mathbf{r})={\mathrm{e}}^{\frac{1}{4}{\mathbf{r}}^T{\Omega }^{T}Y\Omega \mathbf{r}-\mathrm{i}\mathbf{\mu }\Omega \mathbf{r}}.$$

$${\hat{U}}_{BS}=\exp \text{quantity}(\mathrm{i}\frac{\pi }{4}\sum _{i=1}^m\left[{\hat{p}}_{i,A}{\hat{q}}_{i,B}-{\hat{q}}_{i,A}{\hat{p}}_{i,B}\right]).$$

$$S_{BS}=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}\text{mathds}{1}_m& \text{mathds}{1}_m\\ -\text{mathds}{1}_m& \text{mathds}{1}_m\end{array}\right)$$

$${\mathcal{E}}_G(\rho )={\text{Tr}}_C\text{quantity}[{\hat{U}}_G\rho \otimes \text{outerproduct}00^{\otimes k}{\hat{U}}_G],$$

$${\chi }_{{\hat{U}}_{BS}\rho \otimes \rho {\hat{U}}_{BS}^{†}}\text{quantity}(\mathbf{r}=\text{quantity}({\mathbf{r}}_1,{\mathbf{r}}_2{)}^T)={\chi }_{\rho }\text{quantity}(\frac{{\mathbf{r}}_1-{\mathbf{r}}_2}{\sqrt{2}}){\chi }_{\rho }\text{quantity}(\frac{{\mathbf{r}}_2+{\mathbf{r}}_1}{\sqrt{2}}).$$

$$E_{\alpha }(\rho )=\frac{1}{1-\alpha }{\log }_2\text{quantity}(\text{Tr}\text{quantity}[{\text{Tr}}_B\text{quantity}(\rho {)}^{\alpha }])$$

$$\text{Tr}\text{quantity}[\text{quantity}({\text{Tr}}_B\text{quantity}[{\hat{U}}_{BS}\rho \otimes \rho {\hat{U}}_{BS}^{†}]{)}^2]=\frac{1}{(2\pi {)}^m}\int \text{differential}\mathbf{r}\text{absolutevalue}{{\chi }_{\rho }\text{quantity}(\frac{\mathbf{r}}{\sqrt{2}})}^4.$$

$$I_{\alpha }({\rho }_{AB})=\underset{{\sigma }_B}{\min }{D}_{\alpha }\text{quantity}({\rho }_{AB}\Vert {\rho }_A\otimes {\sigma }_B),$$

$$I_1(A;B{)}_{\rho }$$

$$\text{absolutevalue}\text{Tr}\text{quantity}[\rho \sigma ]-\text{Tr}\text{quantity}[{\text{SWAP}}_{2M}\rho \otimes \sigma ]\le 1-q_{2M}\le 1-q_M^{\rho }{q}_M^{\sigma }$$

$$[\hat{q},\hat{p}]=\mathrm{i}.$$

$$\text{quantity}[\hat{\mathbf{r}},{\hat{\mathbf{r}}}^T]=\mathrm{i}\Omega$$

$${\hat{U}}_G$$

$$\hat{H}$$

$$S={\mathrm{e}}^{\Omega \text{H}}$$

$${\hat{U}}_G\hat{\mathbf{r}}{\hat{U}}_G^{†}=S\hat{\mathbf{r}}.$$

$$S\Omega S^T=\Omega .$$

$$\hat{D}(\bar{\mathbf{r}})={\mathrm{e}}^{\mathrm{i}{\bar{\mathbf{r}}}^T\Omega \hat{\mathbf{r}}}$$

$$\hat{D}(\bar{\mathbf{r}})\hat{\mathbf{r}}\hat{D}(\bar{\mathbf{r}}{)}^{†}=\hat{\mathbf{r}}+\bar{\mathbf{r}}.$$

$${\hat{U}}_G=\hat{U}\hat{S}(s)\hat{D}(\mathbf{r})\hat{V},$$

$$\hat{S}(s_i)={\mathrm{e}}^{\mathrm{i}\frac{s_i}{2}\text{quantity}(\hat{q}\hat{p}+\hat{p}\hat{q})}.$$

$${\hat{U}}_R\text{quantity}(\varphi )={\mathrm{e}}^{\mathrm{i}\frac{\varphi }{2}\text{quantity}({\hat{q}}^2+{\hat{p}}^2)}$$

$${\hat{U}}_{BS}(\theta )=\exp \text{quantity}(\mathrm{i}\theta \sum _{i=1}^m\left[{\hat{p}}_{i,A}{\hat{q}}_{i,B}-{\hat{q}}_{i,A}{\hat{p}}_{i,B}\right])$$

$${\rho }_G=\frac{e^{-\beta \hat{H}}}{\text{Tr}[e^{-\beta \hat{H}}]}$$

$$\bar{\mathbf{r}}=\text{Tr}[\hat{\mathbf{r}}\rho ]$$

$$\sigma =\text{Tr}[\{(\hat{\mathbf{r}}-\bar{\mathbf{r}}),(\hat{\mathbf{r}}-\bar{\mathbf{r}}{)}^T\}\rho ]$$

$${\chi }_{\rho }(\mathbf{r})=\text{Tr}\text{quantity}[\hat{D}(\mathbf{r})\rho ].$$

$${\chi }_{{\rho }_G}(\mathbf{r})={\mathrm{e}}^{\frac{1}{4}{\mathbf{r}}^T{\Omega }^T\sigma \Omega \mathbf{r}-\mathrm{i}\bar{\mathbf{r}}\Omega \mathbf{r}}.$$

$${\chi }_{\Phi (\rho )}(\mathbf{r})={\mathrm{e}}^{-\frac{1}{4}{\mathbf{r}}^T{\Omega }^{T}Y\Omega \mathbf{r}-\mathrm{i}\bar{\mathbf{r}}\Omega \mathbf{r}}{\chi }_{\rho }\text{quantity}({\Omega }^T{X}^T\Omega \mathbf{r})$$

$$Y+i\Omega \ge iX\Omega X^T.$$

Full text

Quantum information processing exploits inherently quantum effects to achieve advantages over classical counterparts. Among the various frameworks for implementing quantum protocols, continuous-variable (CV) systems, which are described by infinite-dimensional Hilbert spaces associated with observables possessing a continuous spectrum, are particularly promising. This paradigm has been realized in various experimental platforms, including quantum optical systems [1], microwave cavities coupled to superconducting qubits [2, 3], and trapped ions [4, 5].

Gaussian states and operations play a central role in continuous-variable quantum information processing. Despite the infinite dimensionality of the underlying Hilbert space, Gaussian states can be fully characterized by their first and second moments — namely, the mean vector and the covariance matrix. Moreover, their evolution under Gaussian operations is efficiently described at the level of these moments, enabling powerful analytical techniques. Beyond their mathematical tractability, Gaussian states and operations are readily accessible in experimental platforms, particularly in quantum optics. Their experimental availability has enabled the implementation of various foundational quantum protocols, including quantum teleportation [6, 7, 8] and quantum key distribution [9].

Restricting to Gaussian states and Gaussian operations imposes significant limitations. Within Gaussian quantum information, several no-go theorems prohibit key tasks such as entanglement and resource distillation [10, 11, 12, 13], quantum error correction [14], and the violation of Bell inequalities and contextuality [15]. Additionally, purely Gaussian processes can be efficiently simulated on classical computers [16], which precludes any computational advantage. Therefore, non-Gaussianity, whether in the form of states or operations, is essential for overcoming these limitations and enabling genuinely quantum protocols — though realizing such resources experimentally remains a major challenge.

The operational distinction between Gaussian and non-Gaussian has motivated the development of a resource theory of non-Gaussianity. Resource theories aim to quantify the amount of a given resource — in this case, non-Gaussianity — required to implement specific tasks or protocols. Several quantifiers have been proposed, including the relative entropy of non-Gaussianity [17, 18, 19], robustness [20, 21, 22], stellar rank [23, 24], Gaussian rank or extent [24, 25] and Wigner negativity [26, 27, 28].

These measures not only serve as formal tools within the resource-theoretic framework but also allow one to quantify how non-Gaussian—and thus how genuinely quantum—a given experiment is. However, a major limitation of many existing measures is that their evaluation typically requires full quantum state tomography, which becomes impractical as the number of modes increases [29].

In this work, we introduce a novel framework for defining measures of non-Gaussianity, which allow in some cases efficient experimental access, based on the generation of correlation [30]. Specifically, when two identical copies of a state are sent through a $50{:}50$ beam splitter, the output state is correlated if and only if the input state is non-Gaussian [31] (see Figure 1). This observation enables us to use measures of correlation applied to the output of the beam splitter as operational quantifiers of non-Gaussianity. We establish general monotonicity properties for this class of measures, ensuring their validity within a resource-theoretic framework.

We focus in particular on two such measures: the Rényi-$\alpha$ entropy and Rényi-$\alpha$ mutual information. The Rényi-$2$ entropy, in particular, is experimentally attractive, as it can be efficiently estimated using a SWAP test, even without prior knowledge of the quantum state — making the protocol state-agnostic. Additionally, we derive a lower bound on the sample complexity required to estimate the Wigner negativity and compare it to the sampling cost of estimating the Rényi-2 entropy. Our analysis reveals that the latter can be significantly more efficient.

Overall, our results provide a new operational interpretation of non-Gaussianity through its ability to generate correlation. This connection not only provides an operational interpretation of non-Gaussian resources but also offers practical tools for experimental quantification. Moreover, our approach can be viewed as a generalization of the Hong–Ou–Mandel experiment to arbitrary quantum states and correlation measures.

Background. —

We consider an infinite-dimensional separable Hilbert space consisting of $m$ modes. The canonical operators or quadratures are $\hat{\bm{r}}=(\hat{q}_{1},\hat{p}_{1},...,\hat{q}_{m},\hat{p}_{m})^{T}$ with commutation relation

[TABLE]

with $\Omega=\bigoplus^{m}_{j=1}\begin{pmatrix}0&1\\ -1&0\end{pmatrix}$ being the symplectic form. Note that $\hat{\bm{r}}\hat{\bm{r}}^{T}$ is the outer product. It is often convenient to represent quantum states using the characteristic function $\chi_{\rho}(\bm{r})=\Tr\quantity[\hat{D}(\bm{r})\rho]$, where $\hat{D}(\bm{r})=\mathrm{e}^{\mathrm{i}\bm{r}\Omega\hat{\bm{r}}}$ is the displacement operator. This allows us to define the set of Gaussian states in an instructive manner. We define the set of Gaussian states $\rho_{G}\in\mathcal{G}$ as the states that cen be represented as

[TABLE]

Thus the characteristic function of a Gaussian state is represented by a multivariate Gaussian function. Gaussian unitaries map Gaussian pure states to Gaussian pure states and are generated by Hamiltonians that are at most quadratic in the quadrature operators. The Gaussian unitary that is central to this work is the $m-$mode $50:50$ beam splitter

[TABLE]

Every Gaussian unitary $\hat{U}_{G}$ has an associated symplectic matrix $S_{G}$ that transforms the vector of quadratures as $\hat{U}_{G}\hat{\bm{r}}\hat{U}_{G}^{\dagger}=S_{G}\hat{\bm{r}}$. The symplectic matrix of the beam splitting unitary defined above is

[TABLE]

where the ordering is $\bm{r}=[q_{1},p_{1},...,q_{m},p_{m}]^{T}$. For more details, see Appendix A.

Gaussian channels are completely-positive and trace-preserving maps that map any Gaussian states to another Gaussian states. We can write any Gaussian channel from $n$ modes to $m$ modes using a Stinespring dilation [32] using vacuum states and Gaussian unitaries as

[TABLE]

where the subsystem $C$ contains the last $n+k-m$ modes.

Correlation through beam splitting. —

We start with the central construction of this manuscript. We consider two copies of the same state and act on these states with an $m-$mode $50{:}50$- beam splitter

[TABLE]

Importantly, the output state of the beam splitter is a product state if and only if the input state is Gaussian [31]. This implies that the state $\rho_{AB}=\hat{U}_{BS}\rho\otimes\rho\hat{U}_{BS}^{\dagger}$ is correlated if and only if $\rho$ is a non-Gaussian state. Moreover, in the case of pure states, $\Psi_{AB}=\hat{U}_{BS}\psi\otimes\psi\hat{U}_{BS}^{\dagger}$ is entangled if and only if $\psi$ is a pure non-Gaussian state.

Having the direct connection between non-Gaussianity and correlation generation established, we can define our framework to define measures of non-Gaussianity. Let $C$ be an arbitrary measure of correlation and define $N_{C}(\rho)\coloneqq C(\hat{U}_{BS}\rho\otimes\rho\hat{U}_{BS}^{\dagger})$. By identifying the class of operations on the input state copies that cannot increase the resulting correlations, we can derive general structural properties of our construction. The resulting monotone $N_{C}$ is then a faithfull monotone under Gaussian channels.

The proof can be found in Appendix B.

Another general property is that if the input state is a product state, $\rho=\rho_{A}\otimes\rho_{B}$, then the output of the beam splitter remains factorized., i.e. $\hat{U}_{BS}\rho\otimes\rho\hat{U}_{BS}^{\dagger}=\quantity[\hat{U}_{BS}\rho_{A}\otimes\rho_{A}\hat{U}_{BS}^{\dagger}]\otimes\quantity[\hat{U}_{BS}\rho_{B}\otimes\rho_{B}\hat{U}_{BS}^{\dagger}]$. This is a useful property if the measure $E$ is multiplicative or additive. Otherwise, this property can be useful to derive bounds.

As mentioned above $N_{C}$ is a faithful measures of non-Gaussianity. Considering the fact that preparing the probabilistic mixture of Gaussian states is operationally easily accessible, another concept known as genuine non-Gaussianity has been discussed [27, 28]. Since the set of Gaussian states is non-convex, a probabilistic mixture of two Gaussian states can be non-Gaussian. A state is called genuinely non-Gaussian if it cannot be written as a convex combination of Gaussian states. The natural question in the context of this work is whether our measure serves as a measure of genuine non-Gaussianity. However, as stated before the output is a product state if and only if the input states are Gaussians, which implies that the output is generally correlated. Thus it is not the case by showing that our measure generally takes a non-zero value for a convex mixture of Gaussian states, clarifying that our measure is a valid measure for the resource theory of non-Gaussianity that takes the set of (non-convex) Gaussian states as its free states as in Refs. [17, 18].

With the general construction and their properties in place, we can now focus on specific measures of correlation to quantify the non-Gaussianity of the input states. For pure states the correlation of the state and the entanglement of the state are equivalent. Therefore we consider the entanglement Rényi-$\alpha$ entropy [33]

[TABLE]

for the pure state case. For $\alpha=1$, this reduces to the von Neumann entropy of the reduced state, and for $\alpha=2$ it is directly related to the purity of the reduced state. The Rényi-$\alpha$ entropy serves as an entanglement measure for pure states and is particularly straightforward to compute for $\alpha=2$, as it requires only the reduced-state purity — a quantity that, as we will show later, can be accessed efficiently in experiments.

The direct connection between purity of the reduced state and the non-Gaussianity allows us to find a requirement on the characteristic function. Note that purity is a quantity involving the square of the characteristic function, i.e. $\Tr\quantity[\rho^{2}]=\frac{1}{(2\pi)^{m}}\int\differential\bm{r}\absolutevalue{\chi_{\rho}(\bm{r})}^{2}$, while non-Gaussianity involves the quartic power of the characteristic function

[TABLE]

The measures discussed so far are valid only for pure states, where correlations and entanglement are equivalent. For mixed states, however, one must consider all correlations, not just entanglement, since there exist non-Gaussian states that remain separable but are no longer product after the beam splitter. Further details are provided in Appendix C. In order to treat mixed states as well, we consider the Rényi-$\alpha$ mutual information [34]

[TABLE]

where $D_{\alpha}\quantity(\cdot\lVert\cdot)$ is the minimal Rényi divergence with $\alpha\in[\frac{1}{2},\infty)$ and $\rho_{A}=\Tr_{B}[\rho_{AB}]$. For $\alpha=1$ one recovers the known von Neumann mutual information

[TABLE]

where $D\quantity(\cdot\lVert\cdot)$ is the Umegaki’s relative entropy. Evaluating the Rényi-$\alpha$ mutual information for arbitrary $\alpha$ is challenging as it involves a minimization over states. However, there are upper and lower bounds using Rényi-$\alpha$ entropy and Rényi-$\alpha$ conditional entropy [34], which could admit direct estimation. Consult Ref. [35] for more information on Rényi-$\alpha$ (conditional) entropy.

In Appendix D, we present analytical examples for a variety of states, evaluating the Rényi-$\alpha$ entropy. In Figure 2 we present the entanglement Rényi-$2$ entropy for various states with respect to the mean photon number. Interestingly, we observe that different classes of standard non-Gaussian states take different values, where Fock states achieve the largest. This contrasts with the case of Wigner negativity, where these states were found to have nearly the same values for fixed mean photon number [27]. This means that our measure can provide further restrictions on feasible transformation between different classes of non-Gaussian resource states that could not be detected by the Wigner negativity.

Experimental access. —

Having established the theoretical framework, we now turn to the practical question of how these quantities can be measured in an experimental setting.

As noted previously, Gaussian states are amenable to many analytical tools because they are fully characterized by their covariance matrix and mean vector. However, every quantum state—including non-Gaussian ones—has an associated covariance matrix, even though it does not fully determine the state in the non-Gaussian case. A natural first step experimentally is to measure the covariance matrix and apply established methods to infer partial information about correlations. However, the correlations generated in our setting cannot be captured by the covariance matrix alone.

Consult Appendix E for the proof.

Since methods relying solely on the covariance matrix fail to detect the correlation generated in our approach, we turn to a promising alternative. As demonstrated earlier, the Rényi-$\alpha$ entropy is an effective tool to detect and quantify non-Gaussianity for pure states as then entanglement and correlation are equivalent. Equivalently, for $\alpha=2$ this corresponds to measuring the purity of one of the beam splitter outputs. This approach is closely related to convolution operations in continuous-variable systems [36] and finds analogues in recent works on quantifying non-Gaussianity in fermionic systems [37, 38] and qudits [39, 40, 41].

It is well known that the purity of an unknown quantum state can be determined experimentally via a SWAP test [42], a ubiquitous tool in quantum information processing. The standard implementation involves a controlled-SWAP gate acting on two copies of the state, conditioned on the state of an auxiliary qubit. The purity is then obtained by repeated measurement on the auxiliary qubit. In consequence, measuring non-Gaussianity involves generating two copies of the state after the beam splitter and determining the purity of one subsystem. This procedure requires four copies of the original state we want to test.

A SWAP test using a controlled-SWAP gate has been experimentally realized in a microwave architecture through coupling with a transmon qubit [43]. Building on this, a recent proposal aims to improving the error resilience of the experimental protocol [44]. Rather than implementing a controlled-SWAP gate directly, one can apply beam splitters to the two CV states, combined with a controlled phase-shift between one of the states and a coupled qubit [45]. The purity of the state can then be extracted through repeated measurements of the coupled qubit. All of these implementations of a SWAP-test of CV-modes require coupling to an auxiliary qubit. Coupling between qubits and continuous-variable modes is common in microwave architectures, making this approach particularly well-suited for such platforms. The circuit of a standard SWAP test is depicted in Figure  3a.).

The protocols for SWAP-test discussed so far are not easy to implement in quantum optical experiments, as coupling between CV nodes and qubits remains challenging. However, in some specific cases such as restricting to Fock states controlled-SWAP gates have been implemented using photonic circuits [47]. More generally, one can implement a destructive SWAP-test using photon number resolving (PNR) detectors [48]. A destructive SWAP-test using this approach and Fock states as the input is equivalent to the seminal Hong-Ou-Mandel experiment [48]. More recently a generalization of the Hong-Ou-Mandel experiment has been implemented experimentally using parity measurements to obtain the purity between different Fock states  [49].

While destructive SWAP tests discussed so far have been restricted to Fock states, it has been shown that a destructive SWAP test can be implemented for arbitrary states using a balanced beam splitter and two photon-number-resolving detectors [46]. This approach also works in the multi-mode scenario and was demonstrated experimentally [46]. A circuit of the destructive SWAP test can be found in Figure  3b.). Realistic PNR detectors saturate beyond a finite photon number. We model this limitation by assuming that the detectors can resolve photon numbers up to a total of $m+n\leq 2M$, where $n,m$ are the measured photons of each detector respectively. Without this restriction, only a single PNR detector would be required. Under this assumption, the estimator for the SWAP test restricted to the $2M$-photon Fock subspace takes the form $\Tr\quantity[\text{SWAP}_{2M}\rho\otimes\sigma]=\frac{1}{S}\sum_{s=1}^{S}(-1)^{n(s)}\Theta\quantity(2M-n(s)-m(s))$, where $\Theta$ is the Heaviside function and $n(s),m(s)$ the measurement result of the PNR detectors at round $s$. The systematic error given by the finite detection capability of the PNR detectors can be bounded by

[TABLE]

with $q_{M}^{\rho/\sigma}=\Tr\quantity[\rho/\sigma\sum_{n=0}^{M}\outerproduct{n}{n}]$and $q_{2M}=\Tr\quantity[Q_{2M}\rho\otimes\sigma]$ while $Q_{2M}$ is the projector onto the $2M$ Fock state subspace on $2$ modes.

In addition to this systematic error incurred through the PNR detector, we have a statistical error given by finite measurement statistic. We can estimate the expectation value with standard deviation $\delta$ by measuring $\quantity(\frac{1}{\delta^{2}})$ number of times. The total error is upper bounded by the sum of both errors. Thus in order to bound the total error by $\epsilon$, we require the total measurable photon number $2M$ has to fulfill $1-q_{2M}\leq\epsilon$. This is a requirement on the quality of the PNR detectors. We then have to run the algorithm $\quantity(\frac{1}{\delta^{2}})$ times where $0<\delta\leq\epsilon-\quantity(1-q_{2M})$. The number of samples to obtain error $\epsilon$ therefore depends indirectly on the energy of the states one considers as the error depends how the tail in the Fock state decomposition decays.

In all the cases discussed so far, the number of measurements required —while ignoring systematic errors— to estimate the purity with error $\epsilon$ scales as $\mathcal{O}\quantity(\frac{1}{\epsilon^{2}})$, similar to the cost of estimating a single observable. In the standard SWAP test using qubit measurements, there are only two possible measurement outcomes. In contrast, PNR detectors theoretically yield a countably infinite set of outcomes, which might suggest an energy-dependent scaling. However, the key advantage of using PNR detectors lies in the fact that the purity depends solely on the parity of the detector’s outcome, simplifying the estimation considerably.

Overall, all methods proposed here to measure the non-Gaussianity of an arbitrary and unknown pure state requires four copies at the time, constant depth and $\mathcal{O}\quantity(\frac{1}{\epsilon^{2}})$ rounds. This protocol allows us to measure non-Gaussianity efficiently and state-agnostically in an experiment. The advantage becomes clear when contrasted with full state tomography, which for pure states scales as $\mathcal{O}(\frac{E^{m}}{\epsilon^{2m}})$ [29], where $E$ is the energy of the state.

The mixed-state case presents additional challenges. In particular, an experimental determination of the Rényi-$\alpha$ mutual information is demanding, and we leave the design of a concrete protocol to future work. Nevertheless, a few general remarks are in order. Even though Rényi-$\alpha$ mutual information includes a minimization we can use the bounds provided in [34] to find bounds on the resource content. Then, one can use the expression for the conditional Rényi-$\alpha$ entropies provided in [35] to evaluate the bounds, which could be measured using, e.g., the quantum singular value transformation [50], while thorough investigation on a specific protocol is left for further work. See Appendix F for more details.

To conclude this paragraph, we note that if one cannot have access to several copies simultaneously, one can use classical shadow techniques [51]. See Appendix G for details.

Lower bound on estimating Wigner negativity. —

In the previous section, we have seen that the amount of non-Gaussianity for pure states can be measured efficiently in experiments — in principle independent of the number of modes or the energy. To assess the practical relevance of our proposal, we compare its experimental attainability with that of estimating the Wigner negativity.

Wigner negativity [26, 27, 28] is one of the mostly studied measures for the resource theory of non-Gaussianity and is widely used in experiments to show non-classicality of experimental states [52, 53]. Wigner negativity is a standard benchmark for assessing the quality of experimentally generated states. In practice, it is typically inferred via tomography adapted to the experimental platform, with the negativity then computed from the reconstructed state. However, this is an unfair comparison, as a full state reconstruction contains far more information than a single Wigner negativity value. Consequently, directly estimating the Wigner negativity—without reconstructing the entire state—could in principle be more efficient than state tomography. To level the playing field, we establish a lower bound on the sample complexity required to estimate the Wigner negativity. Consult Appendix H for details. The basic idea of the lower bound is to connect the difference in Wigner negativity of two states with a state discrimination task. Aside from the general connection between estimating Wigner negativity and state discrimination, we use a family of cubic phase states to quantitatively evaluate the lower bound. Using a lower bound on the Wigner negativity of the cubic phase state, we show that the sample complexity diverges if the energy is not taken into account. Restricting to finite-mean photon number states renders the sample lower bound finite. In Figure 4, we illustrate the scaling of the lower bound on estimating the Wigner negativity as a function of the mean photon number. Estimating the Wigner negativity therefore requires a sample complexity that scales at least as the third root of the mean photon number. Comparing this with the lower bound for state tomography, $\mathcal{O}\left(\frac{E}{\epsilon^{2}}\right)$ for a single mode, we see a gap in which a protocol that directly estimates the Wigner negativity could offer an advantage. Nonetheless, this shows the efficiency of our approach to detect and measure non-Gaussianity in an experiment, since the number of samples required to measure non-Gaussianity using our protocols requires constant number of samples independent of the energy.

Discussions. —

In summary, we have introduced a novel framework to detect and quantify non-Gaussianity in continuous-variable systems by linking it to the correlation generated when two copies of a state are mixed on a $50{:}50$ beam splitter. This fundamental connection allows us to define non-Gaussianity measures through the correlation produced at the beam splitter, extending earlier ideas applied to non-classical states [54, 55]. We have established general properties of this construction such as monotonicty under Gaussian channels and illustrated it using the Rényi-$\alpha$ entropy and the Rényi-$\alpha$ mutual information. Experimentally, our framework enables the use of SWAP tests to measure non-Gaussianity with constant sample complexity representing a generalization of the seminal Hong–Ou–Mandel experiment. We also derived a lower bound on the number of sample required to estimate the Wigner negativity, showing the advantage of our proposed protocol.

Looking ahead, it will be valuable to develop continuous-variable–specific methods for directly measuring mutual information in experiments. Another promising direction is to explore the topology of non-Gaussian states within our framework, identifying which states are most resourceful for a given energy or which non-Gaussian states lead to genuinely entangled states and which to separable states. Finally, our sample lower bound on estimating Wigner negativity motivates the design of experimental protocols to estimate it directly—potentially without full tomography—or to establish a formal equivalence between direct estimation and full state reconstruction.

Note added.

During the preparation of this manuscript, we became aware of the independent work of Bu and Li [56].

Appendix A Background: Continuous variables

We start the supplemental material by introducing the necessary background for continuous-variable quantum systems [32] . In this work, we will use the canonical operators $\hat{q},\hat{p}$ that fulfill the following commutation relations

[TABLE]

The commutation relations can be conveniently written for a vector of quadratures of $m-$modes $\hat{\bm{r}}=(\hat{q}_{1},\hat{p}_{1},...,\hat{q}_{m},\hat{p}_{m})^{T}$ as

[TABLE]

with $\Omega=\bigoplus^{m}_{j=1}\begin{pmatrix}0&1\\ -1&0\end{pmatrix}$ being the symplectic form. Note that $\hat{\bm{r}}\hat{\bm{r}}^{T}$ is the outer product.

Gaussian states and Gaussian unitaries play a central role for continuous variables. Gaussian unitaries transform the canonical operators linearly are defined as

[TABLE]

where we used the notation $\bm{\bar{r}}=(r_{q_{1}},r_{p_{1}},...r_{q_{m}},r_{p_{m}})^{T}$ a vector of real numbers. The matrix H is a symmetric matrix, while $\hat{H}$ denotes the operator. Every quadratic Gaussian unitary operator has an associated symplectic matrix

[TABLE]

fullfilling

[TABLE]

Symplectic matrices $S\in\text{Sp}_{2m}(\mathds{R})$ are defined as the matrices that keep the symplectic form invariant

[TABLE]

The operators generated by the Hamiltonian of linear order are the so-called Heisenberg-Weyl operator or Displacement operator

[TABLE]

that transform the canonical operators like

[TABLE]

Any (multimode) Gaussian unitary can, as a consequence of the Euler (sometimes called Bloch–Messiah) decomposition, be decomposed as

[TABLE]

where $\hat{U}$ and $\hat{V}$ are passive Gaussian unitary operators. The unitary $\hat{S}(s)=\bigotimes_{i=1}^{m}\hat{S}(s_{i})$ is called squeezing

[TABLE]

Passive Gaussian unitaries are energy-preserving and can in turn be decomposed into a product of phase rotations

[TABLE]

and beam splitter

[TABLE]

respectively. The unitary implements for $\theta=\frac{\pi}{4}$ a $50{:}50$-beam splitter. In the following we will denote $\hat{U}_{BS}$ being the $50{:}50$-beam splitter.

Analogously, we define Gaussian states as the states that are generated by quadratic Hamiltonian

[TABLE]

including the case $\beta\rightarrow\infty$ being pure states. Gaussian states have the important property that they are fully determined by their mean

[TABLE]

and their covariance matrix

[TABLE]

where $\{,\}$ is the anti-commutator. Note that although Gaussian states are fully determined by their first two statistical moments, higher-order moments generally do not vanish. Non-Gaussian states also possess a well-defined mean and covariance matrix, but they are not fully characterized by these quantities.

Gaussian states can be conveniently be represented the characteristic function

[TABLE]

The characteristic function $\chi_{\rho}:\mathds{R}^{2m}\rightarrow\mathds{C}$ is the symplectic Fourier transform of the Wigner function. The characteristic function fulfills the requirement for a state $\Tr\rho$ that $\chi_{\rho}(0)=1$ and for $\rho^{\dagger}=\rho$ that $\chi_{\rho}(-\bm{r})=\chi_{\rho}(\bm{r})^{*}$. The characteristic function of a Gaussian state with mean $\bm{\bar{r}}$ and covariance matrix $\sigma$ is given by

[TABLE]

Using the characteristic function one can represent the action of a general Gaussian channel $\Phi_{G}$— the maps that map Gaussian states to Gaussian states—

[TABLE]

with the condition.

[TABLE]

Appendix B Monotonicity

In this section, we present the fundamental properties of our construction. These properties can then be used to establish monotonicity, once a measure is chosen. Specifically, we demonstrate which operations are equivalent to operations that act locally after the beam splitter.

B.1 $n$ mode to $n$ mode Gaussian channels

In this section we show monotonicity of the measures using our construction under $n$ mode to $n$ mode Gaussian channels. In later sections we will show monotonicity under partial trace and Gaussian state preparation, which then directly implies monotonicity under Gaussian channels. We first consider a $n$ mode to $n$ mode Gaussian channel without displacements

[TABLE]

with the short-hand notation $\bar{Y}=\Omega^{T}Y\Omega$, $\bar{X}^{T}=\Omega^{T}X^{T}\Omega$. We deal with displacements separately later in the section.

If we apply the Gaussian channel to both output systems of the beam splitter, we have

[TABLE]

If we now instead apply the channel before the beam splitter we have

[TABLE]

If we now apply the beam splitter we get

[TABLE]

It holds that

[TABLE]

and for the exponential term

[TABLE]

We can thus multiply out the characteristic function and get

[TABLE]

which is exactly the same as applying the channel after the beam splitter. We can thus commute the Gaussian channel through the beam splitter. We are left to show how we can commute displacement operators through the beam splitter. It holds that

[TABLE]

Thus we have shown that Gaussian $n$ mode to $n$ mode channels can be moved through the beam splitter to act locally on both outputs. This cannot increase the total correlation. Note that if $Y=0$, the matrix $X$ is a symplectic matrix and the channel describes Gaussian unitary evolution.

B.2 Tracing out subsystem

In this section we will show monotonicity under partial trace. We start by considering as the input a state out of which a sub system $C$ is traced out. The input state is then given as $\sigma=\Tr_{C}\quantity[\rho]$, while the output of the beam splitter is

[TABLE]

This means that the output state with the input $(\Tr_{C}[\rho])^{\otimes 2}$ coincides with the state obtained by tracing out $C$ systems of the output with input $\rho^{\otimes 2}$, whose characteristic function is

[TABLE]

Thus tracing out before the beam splitter is equivalent to tracing out the outputs locally, the correlations cannot increase.

B.3 Adding Gaussian subsystem

We note that the output of tensor product states in our setup is a tensor product of the separate subsystems. As Gaussian states are product states after the beam splitters, we can add Gaussian subsystems without increasing the correlation of the output state.

B.4 Multiplicative output

In this section we show that if the input states are product states then the output will be multiplicative along the same axis. We take two states of the form $\rho\otimes\sigma$ as the input of the beam splitter. The characteristic function before the beam splitter is

[TABLE]

and after the beam splitter

[TABLE]

At the same time, we have

[TABLE]

which is equivalent to the case where we used the states $\rho\otimes\sigma$ as the input. The output states are thus the same as the output of two separate beam splitters.

This can be a very useful property depending on the measure of correlation employed. Some measures such as have the multiplicative property. Other measures do not have this property, but the fact the states factorizes can be used to derive useful bounds.

Appendix C Correlation without entanglement

In this work, we measure the amount of non-Gaussianity in two copies of the state by quantifying the correlations of the state after a $50{:}50$ beam splitter. For pure states, the correlations are equivalent to the entanglement in the output state. This is not true for mixed states. If one wants a faithful measure one has to consider all types of correlations of the output state. In the following we will provide a non-Gaussian state whose state psot beam splitter is a seperable state. More specifically, we consider states of the form

[TABLE]

where the Gaussian state $\rho_{G}$ have zero mean. Then it holds that

[TABLE]

In the second line we used that $\hat{U}_{BS}\rho_{G}\otimes\rho_{G}\hat{U}_{BS}^{\dagger}=\rho_{G}\otimes\rho_{G}$ if $\rho_{G}$ has zero mean, i.e. $\mu$. Indeed we obtain a separable state.

Appendix D Examples

In this section we provide examples for the amount of non-Gaussianity measured by the correlation after the beam splitter. As we are considering pure state we will use the entanglement monotones the von Neumann entropy and the 2-Rényi entropy [33].

The Schmidt decomposition of a pure state $\psi$ is given by

[TABLE]

The entanglement Rényi-$\alpha$ entropy [33] is then

[TABLE]

with the special case $\alpha=1$

[TABLE]

and

[TABLE]

If we have two copies of pure state $\psi$ in the input of the beam splitter then the output is $\chi_{\psi}(\frac{\bm{r_{2}}+\bm{r_{1}}}{\sqrt{2}})\chi_{\psi}(\frac{\bm{r_{2}}-\bm{r_{1}}}{\sqrt{2}})$. The reduced state is then $\chi_{\rho}=\chi_{\psi}(\frac{\bm{r_{1}}}{\sqrt{2}})\chi_{\psi}(\frac{-\bm{r_{1}}}{\sqrt{2}})$ and thus

[TABLE]

In Figure 2 we show the entanglement $2$-Rényi entropy for Fock states, $0N$ states and cubic phase states. The $0N$ state we are considering is defined as a equal superposition between the vacuum and a Fock state $\ket{N}$

[TABLE]

In Appendix H, we discuss additional properties of the cubic phase state; for convenience, we repeat its definition here:

[TABLE]

D.1 Fock states

Here, we derive analytical expressions quantifying the non-Gaussianity of Fock states. The output state for two identical Fock states after the beam splitter is [54]

[TABLE]

with

[TABLE]

Thus the non-Gaussianity measured using the entanglement entropy is

[TABLE]

and the 2-Rényi entropy

[TABLE]

D.2 Cat states

In this section we are providing analytical expression for squeezed cat states of the form $\ket{\psi}=\sum_{i=0}^{1}c_{i}\hat{S}(s)\ket{\alpha_{i}}$, where $\hat{S}(s)$ is the squeezing operator and a coherent state $\ket{\alpha_{i}}$ with real $\alpha_{i}$. See Appendix A for the precise definitions. The output after the beam splitter is given by

[TABLE]

with $\bm{k}=\quantity(i,j)$, $\bar{c}_{\bm{k}}=c_{i}c_{j}$, $\bar{\alpha}_{\bm{k}}=\frac{1}{\sqrt{2}}(\alpha_{j}+\alpha_{i})$, and $\bar{\alpha}_{\bm{k}}^{\prime}=\frac{1}{2}(\alpha_{j}-\alpha_{i})$, where in the third line we used that the squeezing can be commuted through the beam splitter.

As long as $\absolutevalue{\bar{\alpha}_{\bm{k}}-\bar{\alpha}_{\bm{k}^{\prime}}}\gg 1$ or $\absolutevalue{\bar{\alpha}_{\bm{k}}^{\prime}-\bar{\alpha}_{\bm{k}^{\prime}}^{\prime}}\gg 1$ all $\ket{\bar{\alpha}_{\bm{k}},\bar{\alpha}_{\bm{k}}^{\prime}}$ are approximately orthogonal. This holds for standard cat states with large enough $\alpha$. The non-Gaussianity measured through the entanglement entropy and 2-Rényi entropy can then be approximated as

[TABLE]

[TABLE]

Appendix E The correlation cannot be detected by covariance matrix

In this section, we show that the correlation generated by beam splitter cannot be detected solely through the covariance matrix. While every quantum state admits an associated covariance matrix, only Gaussian states are fully characterized by their first and second moments—that is, by their mean vector and covariance matrix.

Recall that the symplectic transformation corresponding to the beam splitter unitary is given by:

[TABLE]

The covariance matrix of the two copies of the input state of the beam splitter is

[TABLE]

for a general—not necessarily Gaussian—state $\rho$. $\sigma_{\rho}=\Tr\quantity[\{\hat{\bm{r}}-\mu,\quantity(\hat{\bm{r}}-\mu)^{T}\}\rho]$ is the covariance matrix of the state $\rho$, where $\mu$ is the mean of the state $\rho$ and $\{,\}$ the anti-commutator.

We then track how the covariance matrix is transformed through the beam splitting interaction. The covariance matrix of the output is

[TABLE]

We can see that the covariance matrix before and after the beam splitter is identical. Thus, we cannot measure or detect the correlation by only considering the covariance matrix.

Appendix F Measuring mutual information

Experimentally accessing the Rényi-$\alpha$ mutual information remains an open problem. Nevertheless, we provide a few comments on possible avenues for its measurement. Although its definition involves a minimization, the bounds established in [34] can be used to obtain estimates of the resource content. Theorem 1 in Ref. [34] states that for $\alpha,\beta>0$ and $\gamma\geq\frac{1}{2}$ satisfying $\frac{\alpha}{\alpha-1}=\frac{\beta}{\beta-1}+\frac{\gamma}{\gamma-1}$ that for $(\alpha-1)(\beta-1)(\gamma-1)>0$

[TABLE]

and for $(\alpha-1)(\beta-1)(\gamma-1)<0$

[TABLE]

We can use the expression provided in [35] to evaluate the bounds:

[TABLE]

for all $\alpha\in(0,1)\cup(1,\infty)$. The conditional entropy $H^{\downarrow}$ has a convenient form for $\alpha=2$

[TABLE]

It could be possible to estimate these entropies using the quantum singular value transform [50]. We leave a specific protocol for future work.

Appendix G Classical shadow

Here, we summarize the use of classical shadows in our framework. Classical shadows are a powerful framework for quantum state estimation with a wide range of applications [57]. Beyond state tomography, classical shadows can be used to estimate quantities such as the von Neumann entropy and the purity of subsystems [57].

Unlike the techniques discussed earlier in this supplemental material—which require simultaneous access to multiple copies of the quantum state—classical shadows only require sequential access to individual copies. This makes them particularly appealing for experimental implementations. This section follows the approach introduced in Ref. [51].

The shadows are obtained by using homodyne detection along a random axis. Assume we want to obtain classical shadows for an $m-$mode quantum state $\rho$. We sample Haar random $m$ matrices $\text{Sp}(2)\cap\text{SO}(2)\cong\text{U}(1)$ $S_{i}$. To this end, we sample random angles $\theta_{1},\dots,\theta_{m}$ to get

[TABLE]

Then we use a copy of the state $\rho$ and apply the rotation $\hat{U}_{S}\rho\hat{U}_{S}^{\dagger}$ with $\hat{U}_{S}=\hat{U}_{S_{1}}\otimes...\otimes\hat{U}_{S_{m}}$ and thus one random rotation or one sampled angle per mode. Afterwards we perform a homodyne measurement along the position axis and obtain $\bm{x}=\quantity(\bm{x_{1}},...,\bm{x_{m}})^{T}\in\mathds{R}^{2m}$. The effect of averaging over all random matrices $S$ and all measurement outcomes $\bm{x}$ and undoing the transformation of $S$ in the post-measurement state $\tilde{\rho}_{x}$, acts like the quantum channel

[TABLE]

For homodyne detection, $\mathcal{E}$ is a simple linear bosonic channel, whose action can be represented as a random displacement channel. The channel is injective and invertible on its range.

We then define the classical shadow as

[TABLE]

$\bar{\rho}_{S,x}$ is an unbiased estimation of $\rho$. The information of $\rho$ is contained in the probability distribution of obtaining measurement outcome $x$.

We only need to perform the initial rotation and the measurement of the quantum state using a quantum device. The rest of the operation, such as the counter-rotation and the inverse of the channel, can be done in an efficient postprocessing, as everything involved is Gaussian.

In the case of homodyne detection, the post-measurement “state” of subsystem $A$ including the inverse of the channel is then given by

[TABLE]

for any $\bm{r_{A}}=\{\bm{r_{j}}\}_{j\in A}\in\mathds{R}^{2\absolutevalue{A}}$ and $(S_{j}\bm{r_{j}})_{2}=\sum_{k}(S_{j})_{2k}(\bm{r_{j}})_{k}$. This characteristic function is not associated to a quantum state. We can however use it to construct a random matrix of size $(M+1)^{\absolutevalue{A}}$ by

[TABLE]

where $\chi_{\ket{n_{1}}\bra{n_{2}}}$ are the matrix elements of the displacement operator in Fock basis.

By repeating this $N$ times and compute the average we obtain the classical shadow

[TABLE]

G.1 Non-linear functionals

Our main interest is to use classical shadows to obtain non-linear functionals of the state $\rho$. The examples we are interested are purity of a subsystem $P(\rho_{A})=\Tr\quantity[\rho_{A}^{2}]$ and the von Neumann entropy of a subsystems $S(\rho_{A})=-\Tr\quantity[\rho_{A}\log\rho_{A}]$. We consider only the reduced state $\rho_{A}$ where $A$ denotes the set of modes containing $\absolutevalue{A}=r$ modes with $r<m$. We assume that the state $\rho$ has local finite energy $\Tr\quantity[(\mathds{1}+n^{j})\rho_{j}]=E_{j}<\infty$ on mode $j$, which is true for any experiment.

von Neumann entropy:

We start with the von Neumann entropy. The derivation can be found in [51]. Given the shadow $\sigma_{A}^{(N)}(M)$ obtained from the scheme described in the previous section and a constant $d_{p}\in\mathds{N}$, we use the following matrix functional

[TABLE]

where $P_{M}$ projects into the subspace of dimensions $M+1$ that contains up to Fock number $M$ in order to approximate the entropy. $M$ thus describes the highest Fock number one takes into account. It holds then that for any $\epsilon>0$ and $M\in\mathds{N}$ that

[TABLE]

with $\epsilon^{\prime}=\frac{\epsilon^{2}}{14(M+1)^{r}\mathrm{e}}2^{-4(M+1)^{r}/\epsilon}$,

[TABLE]

where $\bm{n_{1}},\bm{n_{2}}\in\{0,..,M\}^{r}$, $M_{j}=\max\{n_{1}(j),n_{2}(j)\}$, $m_{j}=\max\{n_{1}(j),n_{2}(j)\}$ and $L^{(k)}_{j}$ being Laguerre polynomials.

Purity:

The other non-linear functional we are interests in is the purity. We can use Theorem 3 [51] to upper bound the number of samples. If the energy of the state is bounded like $E_{r}=\max_{\absolutevalue{A}\leq r}\Tr\quantity[\rho_{A}H_{r}]\leq\infty$, the number of samples is

[TABLE]

and

[TABLE]

then it holds for $L$ observables $O_{j}$ with $\norm{O_{j}}_{\infty}\leq 1$ on regions $A_{j}$ of size at most $r$ that

[TABLE]

with probability at least $1-\delta$.

A quantum state $\rho_{A}$ is an example of such an observable as it holds that $\norm{\rho_{A}}_{\infty}\leq 1$ and then $L=1$.

Appendix H Sampling cost lower bound for estimating Wigner negativity

We are interested in comparing the experimental difficulty of estimating non-Gaussianity using our protocol with that of estimating the Wigner negativity. To this end, we compute a lower bound on the number of samples required for Wigner negativity estimation. The Wigner negativity of the state $\rho$ is defined by

[TABLE]

The standard approach to determining the Wigner negativity of a quantum state in an experiment involves full state or Wigner tomography, followed by reconstruction of the Wigner function. However, such complete reconstructions may be unnecessary when the goal is to estimate the Wigner negativity itself, which contains significantly less information than the full state.

H.1 State discrimination

We prove the lower bound on the number of samples required to obtain the Wigner negativity of a quantum state by connecting the task of estimating the Wigner negativity to a state discrimination task. The optimal probability to succeed in discriminating between the states $\rho$ and $\sigma$ using $n$ copies is given by

[TABLE]

Suppose we can estimate $\mathcal{W}(\rho)$ $\forall\rho$ with additive error $\epsilon$ and failure probability $1-\delta$ with samples $N(\epsilon,\delta)$. Suppose

[TABLE]

then we can discriminate $\rho$ and $\sigma$ with probability $1-\delta$ and $N(\epsilon,\delta)$ samples. This implies

[TABLE]

In particular, if $\rho$ and $\sigma$ are pure states $\psi$ and $\phi$, we have $\norm{\psi^{\otimes n}-\phi^{\otimes n}}_{1}=2\sqrt{1-F(\psi^{\otimes n},\phi^{\otimes n})}=2\sqrt{1-F(\psi,\phi)^{n}}$, which gives

[TABLE]

Then

[TABLE]

which provides a sampling lower bound for estimating the Wigner negativity as

[TABLE]

In the following, we aim to put a general lower bound for the number of samples required to estimate the Wigner negativity of an arbitrary unknown state. We achieve this by putting a lower bound for $N(\epsilon,\delta)$ via (98) for a specific choice of $\psi$ and $\phi$ satisfying (94), which serves as a lower bound for the sampling cost required for estimating arbitrary unknown states. To this end, we consider a class of states known as cubic phase states, which we discuss in the next subsection.

H.2 Cubic phase state

The general bound derived above can now be made explicit by considering concrete examples. In what follows, we focus on the case of cubic phase state. This section provides a summary of the most important properties of the cubic phase state as well as a lower bound on their Wigner negativity.

The cubic phase state,

[TABLE]

is a paradigmatic example of a highly non-Gaussian quantum state. It plays a central role in continuous-variable quantum computation, where it enables universal gate sets. In particular, the cubic phase state allows for the implementation of a non-Clifford $T$-gate on the Gottesman-Kitaev-Preskill (GKP) code space via gate teleportation [58]. More generally, it also features in the Braunstein–Lloyd scheme [59], where teleporting a (possibly noisy or Gaussian-diluted) cubic phase gate provides the necessary nonlinearity for universal quantum computation.

The fidelity between two cubic phase states with same $\gamma$ is given as [27]

[TABLE]

We see that the fidelity between cubic phase states with same $\gamma$ only depends on the squeezing difference. This implies we can increase the difference in Wigner negativity by increasing $c$ and $r$ while keeping the fidelity and thus $r-r^{\prime}$ constant.

The Wigner function of the cubic phase state $\ket{\gamma,r}$ is given by [27, 60]

[TABLE]

It has to be noted that there are many cubic phase states that will have the same Wigner negativity. As is noted in [28] it holds that

[TABLE]

The Wigner negativity of a state $\rho$ is defined by

[TABLE]

Since the Wigner negativity is invariant under Gaussian unitaries, Wigner negativity is constant on lines $x=\gamma\mathrm{e}^{3r}$

[TABLE]

In Figure 5 we show the Wigner negativity for increasing $x=\mathrm{e}^{3r}\gamma$. We see that $\mathcal{W}(\mathrm{e}^{3r}\gamma)$ can be lower bounded by the function $f(x)=x^{\frac{1}{3}}$. Importantly, for large $x$, we observe that the true Wigner negativity grows faster than $x^{\frac{1}{3}}$. This is especially useful for values of $x$ that we cannot compute numerically any longer. We thus want to emphasize that this heuristic lower bound is applicable for large $x$.

As we noted above, there are many cubic phase states with the same Wigner negativity. We are interested in which of those cubic phase states with constant $x$ has the lowest photon number. The mean photon number of a cubic phase state is given by [27]

[TABLE]

where we allowed for a displacement of $P$. Displacing the cubic phase state does not change the Wigner negativity and thus we choose $P=-6\gamma\mathrm{e}^{2r}$. We will write $x=\gamma\mathrm{e}^{3r}$ and thus

[TABLE]

We set $x$ constant and see which $r$ minimizes the mean photon number. The first derivative is

[TABLE]

The minimal mean photon number for a given value of $x$ is then

[TABLE]

Lower bound for Wigner negativity

While an analytic expression for the Wigner negativity is not known, and numerical evaluation is limited for large $x$, asymptotic estimates provide meaningful insight. In this section we provide a lower bound on the Wigner negativity in the limit of large $x$. More specifically, we want to find a lower bound of

[TABLE]

For $x$ large enough, the term $\frac{1}{\sqrt{\frac{1}{2}}^{4}(6x)^{\frac{4}{3}}}$ does not contribute in $\text{Ai}\quantity(\frac{p+3xq^{2}}{\frac{1}{2}(6x)^{\frac{1}{3}}}+\frac{1}{\sqrt{\frac{1}{2}}^{4}(6x)^{\frac{4}{3}}})$. Then

[TABLE]

The Airy function decays faster than exponential for positive argument. For a positive argument and large $x$, the corresponding integral is independent of $x$. We then bound the integral by first neglecting the integration over $q$, since $q$ only appears quadratic in the argument of the Airy function

[TABLE]

For large $x$ we use the asymptotic for negative arguments $\text{Ai}(x)\simeq\frac{\sin\quantity(\frac{2}{3}x^{\frac{3}{2}}+\frac{\pi}{4})}{\sqrt{\pi}x^{\frac{1}{4}}}$. Since we are only interested in lower bounding the Wigner negativity asymptotically for large $x$, we neglect constant factors and obtain

[TABLE]

For large $x$ the term $\sin$ oscillates very fast compared to the exponential decay. We can then approximate it with a constant $c$. In the end we obtain

[TABLE]

We can therefore see that we obtain a crude lower bound for large $x$ as $\norm{W_{\psi}}_{1}(x)\geq x^{\frac{1}{6}}$. In any case, we observe that the Wigner negativity for large $x$ follows a power law. As mentioned before in the section, we see numerically that the function $x^{\frac{1}{3}}$ is a better bound.

H.3 Sampling lower bound using the cubic phase state

Having established a general sample lower bound for estimating Wigner negativity, we now employ the cubic phase state to obtain concrete bounds.

Assume now we have for the Wigner negativity of both cubic phase states that

[TABLE]

so that $\phi$ and $\psi$ satisfy (94) and $\Delta r$ is a small difference in the squeezing of the cubic phase states. As we have seen in the previous section the fidelity between two cubic phase states only depends on the difference in squeezing if the $\gamma$ is the same in both states.

To summarize, the fidelity depends only on the squeezing difference $\Delta r$, whereas the difference in Wigner negativity, and thus in $\epsilon$, is governed by its scaling with $x=e^{3r}\gamma$. As we have shown in the previous section, the Wigner negativity scales faster than a power law with some $\alpha$. Then $\epsilon$ is

[TABLE]

and therefore

[TABLE]

We see that if we have a lower bound on the scaling of the Wigner negativity with $x$ then we have an upper bound on $\Delta r$ and a lower bound on the fidelity $F(\psi,\phi)$. Thus in the limit $x=\mathrm{e}^{3r}\gamma\rightarrow\infty$ we obtain $\mathrm{e}^{3\alpha r}\gamma^{\alpha}\rightarrow\infty$ and in consequence $\Delta r\rightarrow 0$. This means that $F(\psi,\phi)\rightarrow 1$ and the number of samples $N(\epsilon,\delta)$ diverge for constant $\epsilon$ and $\delta$. Since the Wigner negativity of the cubic phase state grows faster than the lower bound we derived, maintaining a constant $\epsilon$ requires $\Delta r$ to decrease more rapidly than in the case of the bound. This means that generally the sampling lower bound diverges, even when only considering one mode.

Therefore, to make a meaningful statement, we must take the mean photon number into account. While keeping the difference in Wigner negativity constant, $\Delta r$ decreases as the Wigner negativity increases.

Let us put all the ingredients together and investigate how the number of samples scale with the mean photon number of the cubic phase state. In general, for fixed accuracy $\epsilon$, $\Delta r$ depends on the Wigner negativity $\mathcal{W}$, which itself is a function of the parameter $x$. This dependence allows us to relate $\mathcal{W}$ to the mean photon number of the state. In order to compute the lower bound we fix $\epsilon,\delta$ and then compute Wigner negativity for increasing $x$ such that the difference is $2\epsilon$. Using these $x$ we compute the mean photon number. In Figure 6, we plot this lower bound on the required number of samples for estimating the Wigner negativity.

To recap, the fidelity depends on the difference in squeezing in two cubic phase states, whose difference in Wigner negativity is constant $2\epsilon$. A lower bound on the scaling of the Wigner negativity yields an upper bound on the squeezing difference, and consequently a lower bound on the bound we obtain by directly evaluating the winger negativity exactly. To reveal how the number of required samples scales with $x$, we will evaluate the lower bound assuming the Winger negativity follows a power law scaling. In the previous section we showed that the Wigner negativity scales for large $x$ at least with $x^{\frac{1}{6}}$, while we see numerically that it grows faster than $x^{\frac{1}{3}}$. Then

[TABLE]

We directly see how $\Delta r$ decreases by increasing the Wigner negativity $\mathcal{W}$. The sampling lower bound is then given as

[TABLE]

We see that a lower bound of the Wigner negativity, especially a lower bound on the scaling with $x$, leads to a lower bound of the number of samples.