framework. Indeed, in the absence of correlations the combination of Eqs. (222) and (189) implies that $W(\vec{x}_g | \vec{x}_f) = W_g(\vec{x}_g)$. As a consequence, we see from Eq. (192) for the conditional expectation value $\langle \hat{A} \rangle_{g|\vec{x}_f} = \langle \hat{A} \rangle$, and thus from Eq. (190) that $W_{f|\hat{A}}(\vec{x}_f) = W_f(\vec{x}_f)$. In other words, the conditional operation has no effect whatsoever and gives the same result as tracing out the modes in $\mathbb{R}^{2\prime}$.

A closer look at the explicit expressions

$$W_{f|\hat{A}}(\vec{x}_f) = \frac{\langle \hat{A} \rangle_{g|\vec{x}_f} W_f(\vec{x}_f)}{\langle \hat{A} \rangle}$$

and

$$\langle \hat{A} \rangle_{g|\vec{x}_f} := (4\pi)^{\prime} \int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W(\vec{x}_g | \vec{x}_f),$$

shows that whenever there are correlations, and thus $\langle \hat{A} \rangle_{g|\vec{x}_f} \neq \langle \hat{A} \rangle$, the conditional Wigner function is a priori non-Gaussian. When we use explicitly that the initial state is Gaussian and thus that $W(\vec{x}_g | \vec{x}_f)$ is given by Eq. (193), this condition can be translated to the existence of nonzero components in $V_{gf}$ in Eq. (185). The precise properties of the resulting non-Gaussian quantum state depend on the conditional expectation value $\langle \hat{A} \rangle_{g|\vec{x}_f}$.

In the literature, some attention has been devoted to proposing different types of measurements for such heralding procedures. One may think of using on-off detectors [202], photon-number-resolving detectors [156], parity detectors [203], and more exotic multimode setups

[204,205]. However, these works usually assume that the initial quantum state is a pure Gaussian state obtained by an idealized source of multimode squeezed vacuum states. As we saw in Sec. 2, for pure-state correlations automatically imply entanglement, and it even turns out that all correlated pure states violate a Bell inequality [206]. In other words, for pure states all correlations are quantum correlations and all these quantum correlations are of the strongest type. When we no longer make such assumptions on the initial multimode Gaussian state, we see that $\langle \hat{A} \rangle_{g|\vec{x}_f}$ will not only depend on the chosen POVM $\hat{A}$, but also on the properties of $W(\vec{x}_g | \vec{x}_f)$. In Secs. 1 and 2, we explain that certain types of non-Gaussian features can only be achieved through certain types of quantum correlations in the initial Gaussian state. An overview of the results of this section is provided in Fig. 5.

## 1. Quantum non-Gaussianity and entanglement

To understand the role of quantum entanglement in a conditional preparation scheme, we contrast it to a system with only classical correlations. In that regard, let us suppose that the initial quantum state is separable such that its Wigner function can be cast in the form Eq. (222). By inserting this form in Eq. (188), we find that

$$W_{f|\hat{A}}(\vec{x}_f) = \int d\gamma p(\gamma) \frac{\int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W_g^{(\gamma)}(\vec{x}_g)}{\int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W_g(\vec{x}_g)} W_f^{(\gamma)}(\vec{x}_f). \quad (240)$$

<figure><img src="image_7.png" /><figcaption>FIG. 5. Different types of quantum correlations are required to be present in the initial Gaussian state $W(\vec{x})$ to create conditional states $W_{f|\hat{A}}(\vec{x}_f)$, as described in Eq. (190), that belong to a certain class. We thus show how the typical hierarchy of quantum correlations (right) can be connected to the structure of the CV state space that was introduced previously in Fig. 1. Throughout Sec. B, we prove that these different types of quantum correlations are necessary resources to achieve different types of states.</figcaption></figure>
As a next step, we define

$$ \tilde{p}_A(\gamma) := p(\gamma) \frac{\int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W_g^{(\gamma)}(\vec{x}_g)}{\int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W_g(\vec{x}_g)}, \quad (241) $$

and show that $\tilde{p}_A(\gamma)$ is a well-defined probability distribution. First, we use the definition of the reduced state

$$ W_g(\vec{x}_g) = \int_{\mathbb{R}^{2\prime}} d\vec{x}_f W(\vec{x}_f \oplus \vec{x}_g) \quad (242) $$

$$ = \int d\gamma p(\gamma) W_g^{(\gamma)}(\vec{x}_g), \quad (243) $$

and thus we immediately find that $\int d\gamma \tilde{p}_A(\gamma) = 1$. Furthermore, we note that $W_g^{(\gamma)}(\vec{x}_g)$ is the Wigner function of a well-defined quantum state $\hat{\rho}_g^{(\gamma)}$ and thus

$$ \int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W_g^{(\gamma)}(\vec{x}_g) = \text{tr}[\hat{\rho}_g^{(\gamma)} \hat{A}] \ge 0. \quad (244) $$

The final inequality follows from the fact that $\hat{A}$ is a positive semidefinite operator. As a consequence, we find that $\tilde{p}_A(\gamma) \ge 0$ for every possible $\gamma$. Thus, we find that for a separable initial state

$$ W_{f|A}(\vec{x}_f) = \int d\gamma \tilde{p}_A(\gamma) W_f^{(\gamma)}(\vec{x}_f). \quad (245) $$

Up to this point, we assumed only that the initial state is separable. As we saw in Sec. C, a mixed quantum state with a positive Wigner function cannot necessarily be decomposed in states with positive Wigner functions. Therefore, we can generally not infer much about the properties of the Wigner function $W_f^{(\gamma)}(\vec{x}_f)$ in Eq. (245).

As a next step, we use the fact that the initial state is also a Gaussian state. Recall from Eq. (103) that any mixed Gaussian state can be decomposed as a mixture of pure Gaussian states. A priori, however, it is not trivial that this decomposition is consistent with decomposition in separable states, Eq. (222). Thus, it remains to show that for Gaussian separable states the Wigner functions $W_f^{(\gamma)}(\vec{x}_f)$ and $W_g^{(\gamma)}(\vec{x}_g)$ in Eq. (222) are also Gaussian.

We start from a crucial observation on covariance matrices that was made in Ref. [207]: whenever an m-mode state with covariance matrix V is separable, there are covariance matrices $V'_f$ and $V'_g$ such that

$$ V \ge \begin{pmatrix} V'_f & 0 \\ 0 & V'_g \end{pmatrix} = V'_f \oplus V'_g. \quad (246) $$

Note that $V'_f$ and $V'_g$ are covariance matrices on the phase spaces $\mathbb{R}^{2\prime}$ and $\mathbb{R}^{2\prime}$, respectively. Nevertheless, $V'_f$ and $V'_g$ are generally not the same as the covariance matrices

$V_f$ and $V_g$ of Eq. (185) that describe the marginal distributions. We should emphasize that the Williamson, Eq. (92), and Bloch-Messiah, Eq. (95), decompositions offer the necessary tools to explicitly construct $V'_f$ and $V'_g$ (we come back to this point in Sec. C). This allows us to use similar techniques as in Eq. (103). Let us first define

$$ W'_f(\vec{x}_f) := \frac{e^{-\frac{1}{2} \vec{x}_f^T V'_f^{-1} \vec{x}_f}}{(2\pi)^m \sqrt{\det V'_f}}, \quad (247) $$

$$ W'_g(\vec{x}_g) := \frac{e^{-\frac{1}{2} \vec{x}_g^T V'_g^{-1} \vec{x}_g}}{(2\pi)^m \sqrt{\det V'_g}}. \quad (248) $$

We can then use Eq. (246) to define a positive definite matrix $V_c := V - V'_f \oplus V'_g$, such that a decomposition of the type Eq. (103) gives us

$$ W(\vec{x}_f \oplus \vec{x}_g) = \int_{\mathbb{R}^{2m}} d\vec{y} \, W'_f(\vec{x}_f - \vec{y}_f) W'_g(\vec{x}_g - \vec{y}_g) \\ \times \frac{e^{-\frac{1}{2} (\vec{y} - \vec{\xi})^T V_c^{-1} (\vec{y} - \vec{\xi})}}{(2\pi)^m \sqrt{\det V_c}}, \quad (249) $$

where we again impose the structure of the bipartition on $\vec{y} = \vec{y}_f \oplus \vec{y}_g$, with $\vec{y}_f \in \mathbb{R}^{2\prime}$ and $\vec{y}_g \in \mathbb{R}^{2\prime}$. Furthermore, recall that $\vec{\xi}$ is the mean field of the initial Gaussian state $W(\vec{x}_f \oplus \vec{x}_g)$. The structure we obtain in Eq. (249) exactly corresponds to Eq. (222), where $\vec{y}$ now labels the states and thus plays the role of the abstract variable $\gamma$.

We can then use the structure Eq. (249) in the derivation Eq. (245) and then we find that

$$ W_{f|A}(\vec{x}_f) = \int_{\mathbb{R}^{2\prime}} d\vec{y}_f W'_f(\vec{x}_f - \vec{y}_f) \tilde{p}_A(\vec{y}_f). \quad (250) $$

In any concrete choice of $\hat{A}$, one can use Eq. (241) to derive an explicit expression for $\tilde{p}_A(\vec{y}_f)$, which will generally be a non-Gaussian probability distribution, such that $W_{f|A}(\vec{x}_f)$ describes a non-Gaussian state. However, the resulting conditional state, Eq. (250), is clearly a statistical mixture of Gaussian states and thus lies in the convex hull of Gaussian states. In the language of Sec. C this means that the conditional state is non-Gaussian but not quantum non-Gaussian and has a stellar rank 0.

In summary, we have assumed that our initial state with Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ is a separable Gaussian state. Without making any assumptions on the POVM element $\hat{A}$ of the conditional operation, we retrieve that the conditional state always is a convex combination of Gaussian states, given by Eq. (250). Thus, when the initial state is Gaussian, entanglement is a necessary resource to produce quantum non-Gaussian states via conditional operations.
## 2. Wigner negativity and Einstein-Podolsky-Rosen steering

In Sec. E, we explained that Wigner negativity is a “stronger” non-Gaussian feature than quantum non-Gaussianity. Here, we show that also stronger types of quantum correlations are required to conditionally create Wigner negativity. To understand how Wigner negativity can be achieved through a conditional preparation scheme, it suffices to understand when the conditional expectation value $\langle \hat{A} \rangle_{g|\vec{x}_f}$ in Eq. (192) reaches negative values.

Regardless of the chosen POVM, $W_A(\vec{x}_g)$ is the Wigner function of a positive semidefinite operator $\hat{A}$ as defined by Eq. (42). Thus, whenever there is a quantum state $\hat{\rho}'$ that has $W(\vec{x}_g | \vec{x}_f)$ as associated Wigner function, Eq. (43) implies that $\langle \hat{A} \rangle_{g|\vec{x}_f} = \text{tr}[\hat{\rho}' \hat{A}] \ge 0$. Hence, to conditionally create a nonpositive Wigner function (190) the conditional probability distribution $W(\vec{x}_g | \vec{x}_f)$ cannot be a well-defined Wigner function. This observation holds whenever the initial state has a positive Wigner function.

When in addition we assume that the initial state is Gaussian, we find that $W(\vec{x}_g | \vec{x}_f)$ is a Gaussian distribution (193). Whether the conditional probability distribution Eq. $W(\vec{x}_g | \vec{x}_f)$ describes a Gaussian quantum state depends entirely in the properties of its covariance matrix, i.e., the Schur complement $V_{g|\vec{x}_f}$. Indeed, $W(\vec{x}_g | \vec{x}_f)$ describes a quantum state if and only if $V_{g|\vec{x}_f}$ satisfies the Heisenberg inequality. Because $V_{g|\vec{x}_f}$ does not depend on the choice $\vec{x}_f \in \mathbb{R}^{2l}$, it follows that $W(\vec{x}_g | \vec{x}_f)$ corresponds to a quantum state either for all $\vec{x}_f \in \mathbb{R}^{2l}$ [if the Schur complement Eq. (194) satisfies the Heisenberg inequality] or for none of the $\vec{x}_f \in \mathbb{R}^{2l}$ [if the Schur complement Eq. (194) violates the Heisenberg inequality].

If $V_{g|\vec{x}_f}$ satisfies the Heisenberg inequality, the conditional state's Wigner function $W_{f|A}(\vec{x}_f)$ must thus be positive. To better understand the physical resources required to conditionally create Wigner negativity, one must comprehend what it means for $V_{g|\vec{x}_f}$ to violate Heisenberg's inequality in terms of quantum correlations. It turns out that this condition is closely related to the original argument of the EPR paper [172]. The violation of Heisenberg's inequality by the Schur complement $V_{g|\vec{x}_f}$ corresponds to Gaussian quantum steering in the state $W(\vec{x}_g \oplus \vec{x}_f)$.

To understand the connection between the conditional covariance matrix $V_{g|\vec{x}_f}$ and quantum steering, we first express the Wigner function obtained by conditioning on a Gaussian measurement, such that the associated POVM element has a Wigner function $W_G(\vec{x}_f)$:

$$W_{g|G}(\vec{x}_g) = \frac{\int_{\mathbb{R}^{2l}} d\vec{x}_f W_G(\vec{x}_f) W(\vec{x}_f | \vec{x}_g)}{\int_{\mathbb{R}^{2l}} d\vec{x}_f W_G(\vec{x}_f) W_f(\vec{x}_f)} W_G(\vec{x}_g). \quad (251)$$

In a very similar way, we can also show that

$$W_{g|G}(\vec{x}_g) = \int_{\mathbb{R}^{2l}} d\vec{x}_f \frac{W_G(\vec{x}_f) W_f(\vec{x}_f)}{\int_{\mathbb{R}^{2l}} d\vec{x}_f W_G(\vec{x}_f) W_f(\vec{x}_f)} W(\vec{x}_g | \vec{x}_f). \quad (252)$$

Hence, when $W(\vec{x}_g | \vec{x}_f)$ is a bona fide Wigner function for every $\vec{x}_f$ this expression is an explicit manifestation of the local hidden variable model Eq. (230). In other words, whenever $W(\vec{x}_g | \vec{x}_f)$ describes a quantum state, the modes in $g$ cannot be steered by Gaussian measurements on the modes $f$. Note that we can generalize Gaussian measurements to any measurement with a positive Wigner function.

The remarkable feature of EPR steering is that the inverse statement also holds: when $W(\vec{x}_g | \vec{x}_f)$ is not a bona fide Wigner function Gaussian measurements can steer the state. Let us assume that Eq. (230) holds for Gaussian measurements. It then follows that a well-defined covariance matrix $U$ exists such that the covariance matrix $V_{g|G}$ of the conditional state $W_{g|G}(\vec{x}_g)$ satisfies $V_{g|G} \ge U$ for all Gaussian measurements. Furthermore, $U$ is physical and satisfies the Heisenberg inequality. Reference [176] shows that the existence of such a $U$ implies that the full covariance matrix of the system satisfies $V + 0_f \oplus i\Omega_g \ge 0$, which in turn implies that $V_{g|\vec{x}_f}$, the Schur complement of $V$, satisfies the Heisenberg inequality.

This shows that we can only generate Wigner negativity through Eq. (190) if the initial state can be steered by Gaussian measurements on the subsystem associated with phase space $\mathbb{R}^{2l}$. Note that the creation of Wigner negativity occurs in the opposite direction to the steering: We can produce Wigner negativity in the modes $f$ by performing a measurement on the modes $g$ if the modes $g$ can be steered by performing Gaussian measurements on the modes $f$. Somewhat counterintuitively, it turns out that the created Wigner negativity volume, Eq. (156), is not directly proportional to the strength of EPR steering [208].

As a final remark, we note that, in a multimode context, EPR steering is constrained by monogamy relations [209–211]. Notably, this implies that when a single mode $g$ can be steered by a single other mode $f$, it is impossible for any other mode to also steer $g$. This naturally has profound consequences for the conditional generation of Wigner negativity that we discussed in this section. The monogamy relations for quantum steering can be used to derive similar monogamy relations [208] for the created Wigner negativity volume (156).

## C. Quantum correlations through non-Gaussianity

In Sec. B, we extensively considered the use of quantum correlations as a resource to create non-Gaussian effects. In this subsection, we focus on the opposite idea where non-Gaussian operations increase or even create quantum
correlations. The subject of entanglement in non-Gaussian states is generally difficult to study, for some states it may be sufficient to evaluate lower-order moments [212] and when the density matrix in the Fock representation is available one can apply DV approaches to characterize entanglement [213]. However, these methods cannot always be applied and there are no universally applicable entanglement criteria that are practical to evaluate for arbitrary CV quantum states.

## 1. Entanglement measures on phase space

In Sec. 2, we argued that any pure state that manifests correlations between subsystems contains entanglement. Measuring entanglement in this case becomes equivalent to measuring the amount of correlation within the pure state. In particular, for pure states, one finds a wide range of entanglement measures in the literature [161]. In the case of CV systems, some measures are more appropriate than others, and here we focus on one particularly intuitive measure that is based on purity.

When we consider an arbitrary bipartite pure quantum state with Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ (with $\vec{x}_f \in \mathbb{R}^{2l}$ and $\vec{x}_g \in \mathbb{R}^{2l''}$), we find that its purity is $\mu = 1$ by definition. However, this is not necessarily true for the subsystems **f** and **g**. We can use Eq. (94) to evaluate the purity of any state based on its Wigner function, and we define

$$ \mu_f = \int_{\mathbb{R}^{2l}} d\vec{x}_f [W_f(\vec{x}_f)]^2, \quad \text{and} \quad \mu_g = \int_{\mathbb{R}^{2l''}} d\vec{x}_g [W_g(\vec{x}_g)]^2, \qquad (253) $$

where we again use the definitions (182), (183). Because the global state with Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ is pure, we always find that $\mu_f = \mu_g$ (this is a general consequence of the existence of a Schmidt decomposition for pure states). Furthermore, if the pure state is separable, we find $W(\vec{x}_f \oplus \vec{x}_g) = W_f(\vec{x}_f)W_g(\vec{x}_g)$ and as a consequence we obtain that $\mu_f = \mu_g = 1$. However, when $\mu_f = \mu_g < 1$ there must be correlations between the subsystems **f** and **g** and the smaller the purity of the subsystems, the stronger these correlations are. Without delving into the details, we stress that the opposite notion also holds: when there is a correlation between the subsystems, the purity of the subsystems is smaller than one.

To convert this quantity into an entanglement measure [214], it is useful to define the Rényi-2 entropy for subsystem **f**

$$ S_R := -\log \mu_f. \qquad (254) $$

We then find that $S_R \ge 0$ and $S_R = 0$ if and only if the state is separable. Furthermore, it should be clear that $S_R$ cannot be increased by local unitary operations on the subsystems **f** and **g**. We can thus define an entanglement measure

for the pure state $|\Psi\rangle$ with Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ by setting

$$ \mathcal{E}_R(|\Psi\rangle) := S_R. \qquad (255) $$

This constitutes a well-defined entanglement measure for any chosen bipartition and any pure state on the phase space.

To extend this measure to mixed states, we follow Ref. [214] and construct a convex roof. Any mixed state $\hat{\rho}$ can be decomposed in pure states as $\hat{\rho} = \int d\gamma p(\gamma) |\Psi^{(\gamma)}\rangle\langle\Psi^{(\gamma)}|$, we abbreviate this decomposition as the ensemble $\{p(\gamma), |\Psi^{(\gamma)}\rangle\}$. For each pure state in this ensemble, we can evaluate the entanglement $\mathcal{E}_R(|\Psi^{(\gamma)}\rangle)$ and subsequently average all of these values according to $p(\gamma)$. However, the decomposition of $\hat{\rho}$ in pure states is far from unique and different ensembles $\{p(\gamma), |\Psi^{(\gamma)}\rangle\}$ generally lead to a different value of entanglement even though they are all constrained to produce the same state $\hat{\rho}$. Therefore, it is common to define

$$ \mathcal{E}_R(\hat{\rho}) := \inf_{\{p(\gamma), |\Psi^{(\gamma)}\rangle\}} \int d\gamma p(\gamma) \mathcal{E}_R(|\Psi^{(\gamma)}\rangle) \qquad (256) $$

as the general “Rényi-2 entanglement” of the state $\hat{\rho}$.

Formally, this is an elegant definition that can in principle be calculated directly from the Wigner function. However, in practice it is nearly impossible to actually identify all possible decompositions $\{p(\gamma), |\Psi^{(\gamma)}\rangle\}$, which makes this measure notoriously hard to evaluate for mixed states. This has sparked some alternative definitions of entanglement measures for Gaussian states, where any Gaussian state can be decomposed in an ensemble of Gaussian states, Eq. (103). Thus, one can define “Gaussian Rényi-2 entanglement” by restricting Eq. (256) to only Gaussian decompositions [215]. In this sense, Gaussian Rényi-2 entanglement is by construction an upper bound to the general Rényi-2 entanglement.

As an alternative to entanglement measures, it is common to use entanglement witnesses. These have been particularly successful for Gaussian states [146,216–220], where one commonly applies methods based on the covariance matrix of the state. Due to the extremality of Gaussian states [12] these results also provide witnesses for entanglement if the state is non-Gaussian. However, there are several examples of non-Gaussian entangled states for which no entanglement can be detected from the covariance matrix. Notable progress was made by developing entanglement witnesses for non-Gaussian states with specific structure in their Wigner function [221].

It is noteworthy to emphasize that the positive-partial transpose (PPT) criterion of Ref. [216] can in principle be implemented on the level of Wigner functions. To make this apparent, let us first define the transposition operator $T$ that implements $\hat{\rho} \mapsto \hat{\rho}^T$. When $W(\vec{x})$ with $\vec{x} \in \mathbb{R}^{2m}$
denotes the Wigner function of the state $\hat{\rho}$, we can write the Wigner function of $\hat{\rho}^T$ as $W(T\vec{x})$. The matrix $T$ can be written as

$$T = \bigoplus_{i=0}^{m} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \qquad (257)$$

which can be derived from the definition of the Wigner function [222]. The concept of partial transposition in entanglement theory relies on the simple idea that one can apply a transpose only on one of the two subsystems in the bipartition. In our context, this means that the Wigner function changes as $W(\vec{x}_f \oplus \vec{x}_g) \mapsto W(\vec{x}_f \oplus T\vec{x}_g)$ (where $T$ is now taken only on the $l'$ modes of subsystem g). The PPT criterion is based on the idea that, in absence of entanglement, the function $W(\vec{x}_f \oplus T\vec{x}_g)$ still gives a well-defined Wigner function of a quantum state. However, there are entangled states for which this is no longer true and $W(\vec{x}_f \oplus T\vec{x}_g)$ becomes unphysical. This lack of physicality is expressed by the fact there exist positive semidefinite operators $\hat{A}$ for which

$$(4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x} \, W_A(\vec{x}) W(\vec{x}_f \oplus T\vec{x}_g) < 0. \qquad (258)$$

Finding such observables $\hat{A} \ge 0$ for a non-Gaussian state $W(\vec{x}_f \oplus \vec{x}_g)$ is generally a very hard task. For Gaussian states, on the other hand, the physicality of $W(\vec{x}_f \oplus T\vec{x}_g)$ is simply checked through Heisenberg's inequality. For more general non-Gaussian states, this is insufficient and one should check a full hierarchy of inequalities instead [212]. Nevertheless, one may yet uncover more direct methods to check the properties of $W(\vec{x}_f \oplus T\vec{x}_g)$.

## 2. Entanglement increase

One of the most well-known protocols for increasing entanglement is entanglement distillation. In this protocol, one acts with local operations on a large number of mixed entangled states that are shared by two parties and concentrates the entanglement in a smaller number of maximally entangled pairs [223]. When the initial states are pure and the local operation serves only to increase the entanglement and not the purity, we speak of entanglement concentration [224]. Conditional operations play an important role in these protocols, and we can alternatively think of entanglement distillation as the idea that a conditional operation can increase the entanglement of a state. For Gaussian quantum states, there is a notoriouse no-go theorem that states that Gaussian measurements (or Gaussian operations in general) cannot increase bipartite entanglement [225–227]. It was quickly realized that these no-go results can be circumvented by even the most basics non-Gaussian states: those created through a non-Gaussian noise process [228,229]. On the other hand, if

one wants to distill entanglement in a CV system starting from initial Gaussian states one really requires non-Gaussian operations. One such example is given in Refs. [230,231], where the authors propose to use a Kerr non-linearity to distill entanglement for mixed Gaussian states. In contrast, conditional schemes have also been proposed [147,232,233], avoiding the need for optical nonlinearities. In those protocols, one first uses conditional operations to create non-Gaussian states and subsequently uses Gaussianization to obtain states with higher entanglement. A narrowly related protocol [234] relies on the implementation of noiseless linear amplification [235], where the non-Gaussian element is injected in the form of auxiliary Fock states.

The realization that photon subtraction and addition can be used to increase the entanglement of a Gaussian input state was developed reasonably early [236–238] and was further formalized in works such as Refs. [239–241]. Remarkably, all of these works explicitly assume that the initial state under consideration is a two-mode squeezed state and the approach strongly relies on the structure of this type of state in the Fock basis. Beyond the two mode setting, the class of CV graph states has also been studied in the context of entanglement increase [242,243]. Here we provide an alternative approach, based on phase-space representations to understand entanglement increase due to the subtraction of a single photon.

Our approach relies on the fact that we can easily apply the entanglement measure (256) when the global state is pure. This means that we are focusing on a context of entanglement concentration. Furthermore, when we perform photon subtraction on a pure Gaussian state, the resulting photon-subtracted state is also pure, as we saw in Sec. 2. The starting point is the Wigner function of the photon-subtracted state, Eq. (214), which we rewrite as

$$W^{-}(\vec{x}) = \frac{W(\vec{x})}{\mathrm{tr}(V_{\mathrm{b}} - 1) + \|\vec{\xi}_{\mathrm{b}}\|^{2}} \left( \|B^{T}(1 - V^{-1})(\vec{x} - \vec{\xi}) + \vec{\xi}_{\mathrm{b}}\|^{2} + \mathrm{tr}[1 - B^{T}V^{-1}B] \right). \qquad (259)$$

The state $W^{-}(\vec{x})$ is thus obtained by subtracting a photon from the Gaussian state $W(\vec{x})$. As we consider a pure two-mode state we assume that the state has a $4 \times 4$ covariance matrix of the form $V = S^T S$, where $S$ is a symplectic matrix. We assume that the photon is locally subtracted in one of the modes of the mode basis, such that

$$B = \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{pmatrix}. \qquad (260)$$

However, to assess the entanglement in the system, we must obtain the Wigner function for the reduced state associated to either of the two modes. When we focus on mode
b where the photon is subtracted, we can simply obtain the reduced photon subtracted state $W_b^-(\vec{x}_b)$ by subtracting a photon from the reduced Gaussian state $W_b(\vec{x}_b)$. As such, we obtain

$$W_b^-(\vec{x}_b) = \frac{\|(1 - V_b^{-1})\vec{x}_b + V_b^{-1}\vec{\xi}_b\|^2 + \text{tr}[1 - V_b^{-1}]}{\text{tr}(V_b - 1) + \|\vec{\xi}_b\|^2} \\ W_b(\vec{x}_b). \qquad (261)$$

This is now a single-mode photon subtracted state, but it is no longer pure. This lack of purity is notably reflected by $V_b$, which is no longer symplectic. Nevertheless, we can use the Williamson decomposition, Eq. (92), and write

$$V_b = \nu \begin{pmatrix} r & 0 \\ 0 & r^{-1} \end{pmatrix}, \qquad (262)$$

where we set the phase such that the squeezing coincides with one of the axes of phase space. What remains is for us to calculate the purity

$$\mu_b^- = 4\pi \int_{\mathbb{R}^2} d\vec{x}_b [W_b^-(\vec{x}_b)]^2. \qquad (263)$$

The final expression for the purity is not very insightful. When on top we use that the purity $\mu_b$ of the Gaussian state $W_b(\vec{x}_b)$ is given by $\mu_b = 1/\nu$, an explicit calculation of $\mu_b^-$ makes it possible to prove (the motivated reader can use a combination of patience and software for symbolic algebra to do so) that

$$\frac{\mu_b^-}{\mu_b} \le \frac{1}{2}. \qquad (264)$$

In other words, photon subtraction reduces the purity at most by a factor of 2.

When we use Eq. (256) to define the entanglement of the two-mode photon-subtracted state, Eq. (259), we find that it is given by

$$\mathcal{E}_R(|\Psi^-\rangle) = -\log \mu_b^-, \qquad (265)$$

because the two-mode state is pure. The entanglement of the initial Gaussian state is given by $\mathcal{E}_R(|\Psi_G\rangle) = -\log \mu_b$, such that we can use Eq. (264) to find that

$$\Delta\mathcal{E}_R := \mathcal{E}_R(|\Psi^-\rangle) - \mathcal{E}_R(|\Psi_G\rangle) \le \log 2. \qquad (266)$$

In other words, photon subtraction can increase the Rényi-2 entanglement of an arbitrary Gaussian state, but at most by an amount $\log 2$. It turns out that this result can be generalized to all bipartitions of Gaussian pure states of an arbitrary number of modes [244]. Furthermore, the same work shows that when the entanglement measure $\mathcal{E}_R(|\Psi_G\rangle)$ is replaced with the Gaussian Rényi-2 entropy of Ref.

[215], the result holds for all bipartitions of all Gaussian states (including mixed ones).

For the particular case of a two-mode pure Gaussian state, we can directly evaluate $\Delta\mathcal{E}_R$ for some important examples. Say, for example, that we consider the EPR state that is obtained by mixing two oppositely squeezed vacuum states on a balanced beam splitter. In this case $\vec{\xi} = 0$ and $V$ is given by

$$V = \frac{1}{2} \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \end{pmatrix}^T \begin{pmatrix} s & & & \\ & s^{-1} & & \\ & & s^{-1} & \\ & & & s \end{pmatrix} \\ \times \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \end{pmatrix} \\ = \frac{1}{2s} \begin{pmatrix} s^2+1 & 0 & s^2-1 & 0 \\ 0 & s^2+1 & 0 & 1-s^2 \\ s^2-1 & 0 & s^2+1 & 0 \\ 0 & 1-s^2 & 0 & s^2+1 \end{pmatrix}. \qquad (267)$$

We then extract directly that

$$V_b = \frac{s^2 + 1}{2s} 1, \qquad (268)$$

such that we find that the parameters in Eq. (262) are set to $r = 1$ and $\nu = (s^2 + 1)/(2s)$. And thus we directly obtain

$$\Delta\mathcal{E}_R = \log(2) - \log\left(\frac{s^4 + 6s^2 + 1}{(s^2 + 1)^2}\right). \qquad (269)$$

We clearly see that the entanglement increase vanishes in absence of squeezing, whereas we achieve the $\log(2)$ limit for $s \to \infty$. Adding a mean field with $\vec{\xi}_b \ne 0$ immediately complicates the problem. As can be seen in Fig. 6, where we plot the case $\vec{\xi}_b = (0, 1)^T$, the presence of a mean field in the mode of photon subtraction lowers the entanglement increase $\Delta\mathcal{E}_R$. Nevertheless, in the limit $s \to \infty$ we reach the limit $\log(2)$ regardless of the displacement.

This example clearly shows that photon subtraction can be used as a tool to increase entanglement. The setting corresponds to the case that is typically studied in most works on CV entanglement distillation such as Refs. [239–241]. It turns out that one can further increase entanglement in such systems by subtracting more photons. Furthermore, photon addition and the combination of addition and subtraction on both modes have also been considered. The methods we use in this Tutorial are not easily generalized to the subtraction and addition of many photons, but in return they can be applied to a much wider class of initial Gaussian states.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Entanglement increase, Eq. (266), through photon subtraction in one mode of a pair of entangled modes. The initial Gaussian states are obtained by mixing either two equally squeezed modes (yellow curves) or one squeezed mode and one vacuum mode (red curves) on a beam splitter (see also sketches on the right). We show how a variation of squeezing (in dB compared to shot noise level) in these initial squeezed vacuum states influences the entanglement increase due to photon subtraction. We consider cases without mean field (solid curves) and with a mean field $\tilde{\xi} = (0, 0, 0, 1)^T$ (dashed curves).</figcaption></figure>

As a second example, we consider a single-mode squeezed state that is split in two on a balanced beam splitter. This means that the Gaussian state is given by

$$ V = \frac{1}{2} \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \end{pmatrix}^T \begin{pmatrix} s & & & \\ & s^{-1} & & \\ & & 1 & \\ & & & 1 \end{pmatrix} \\ \times \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \end{pmatrix} \\ = \frac{1}{2} \begin{pmatrix} s+1 & 0 & s-1 & 0 \\ 0 & \frac{s+1}{s} & 0 & \frac{1}{s}-1 \\ s-1 & 0 & s+1 & 0 \\ 0 & \frac{1}{s}-1 & 0 & \frac{s+1}{s} \end{pmatrix}, \quad (270) $$

such that we get

$$ V_b = \frac{1}{2} \begin{pmatrix} s+1 & 0 \\ 0 & \frac{s+1}{s} \end{pmatrix}, \quad (271) $$

such that we find that we identify the parameters of Eq. (262) as $\nu = \sqrt{2+s+s^{-1}}/2$ and $r = (1+s)/\sqrt{2+s+s^{-1}}$. In absence of any mean field, i.e., with $\tilde{\xi} = 0$, we then find an entanglement gain given by

$$ \Delta\mathcal{E}_R = \log(2) - \log\left(\frac{3+2s+3s^2}{2(s+1)^2}\right). \quad (272) $$

Interestingly, in this case we reach the maximal entanglement gain for vanishing squeezing $s \to 1$, where we reach

$\Delta\mathcal{E}_R \to \log(2)$. This case may seem somewhat counterintuitive, but it should be emphasized that the success probability of photon subtraction also vanishes in this case. Yet, our conditional approach assumes that we are in the scenario where a photon was subtracted and the negligible fraction of the state that is not in vacuum is enhanced. In the limit of vanishing squeezing, the photon subtracted state converges to the Bell state $(|1, 0\rangle + |0, 1\rangle)/\sqrt{2}$. On the other hand, in the limit where squeezing is high we still find a finite entanglement increase as $\Delta\mathcal{E}_R \to \log(4/3)$.

When we add a mean field given by $\tilde{\xi}_b \neq 0$, there is an importance of the phase because our state locally has some remaining asymmetry (which can be seen from $r \neq 1$). In Fig. 6 we particularly show the case where $\tilde{\xi}_b = (0, 1)^T$ such that the direction of the displacement coincides with the quadrature where the noise is minimal. In this case we observe that for some values of initial squeezing, the entanglement decreases due to photon subtraction. Note that this quite remarkably implies that in some cases photon subtraction can actually be used to increase the purity of a state.

We thus showed that photon subtraction is a useful non-Gaussian operation to increase entanglement. However, in the presence of a mean field in the subtraction mode, it is also possible to decrease entanglement. Even though this subject has been studied for nearly two decades, for arbitrary Gaussian input states, there are still many open questions. Notably, there has not been much work on the effect of photon subtraction on multipartite entanglement, nor on stronger types of quantum correlations. Our discussion in Sec. 2 suggests an important interplay between EPR steering and Wigner negativity, and thus it is intriguing to wonder whether well-chosen non-Gaussian operations can increase quantum steering. Since all steerable states are also entangled, it is a reasonable conjecture that some of the protocols that can increase quantum entanglement should also increase quantum steering.

We have followed the terminology found in the literature and referred to this process as entanglement distillation, because our conditional operation has only a finite success probability. This implies that we can use a large number of Gaussian entangled states and use photon subtraction to obtain a much smaller number of more entangled states. Yet, it must be stressed that there is a more subtle process happening: the entanglement is increased by adding non-Gaussian entanglement on top of the existing Gaussian entanglement. To get a better grasp of this non-Gaussian entanglement, it is useful to go to a setting where no other type of entanglement is present as we do in Sec. 3.

### 3. Purely non-Gaussian quantum entanglement

In this subsection, we explore an idea that is in many ways complementary to the previous subsection: rather
than using a local non-Gaussian operation to increase already existing entanglement, we now use a nonlocal non-Gaussian operation to create entanglement between unentangled modes.

Let us again assume that our state is initially Gaussian as described by Eq. (85), and we induce the non-Gaussian features through the conditional methods of Sec. B. The mean field of the initial state is given by $\tilde{\xi} = \tilde{\xi}_f \oplus \tilde{\xi}_g$, and

$$V = \begin{pmatrix} V_f & V_{fg} \\ V_{gf} & V_g \end{pmatrix}, \quad \text{with} \quad V_f = V_{f1} \oplus V_{f2}. \quad (273)$$

Here, we have introduced the modes of interest, labeled by **f** and a set of auxiliary modes **g** upon which a measurement will be performed to induce non-Gaussian features in the modes **f**. In the initial state, we consider a bipartition in the modes **f** without any direct correlations, hence $V_f = V_{f1} \oplus V_{f2}$. In other words, the modes in **f**1 are completely uncorrelated to the modes in **f**2.

To induce non-Gaussian effects, we resort to the conditional framework by acting with a POVM element $\hat{A}$ upon the auxiliary modes **g**, and we rewrite Eq. (190) as

$$W_{f|A}(\vec{x}_{f1} \oplus \vec{x}_{f2}) = \frac{\langle \hat{A} | g|\vec{x}_{f1} \oplus \vec{x}_{f2} \rangle}{\langle \hat{A} \rangle} W_{f1}(\vec{x}_{f1}) W_{f2}(\vec{x}_{f2}), \quad (274)$$

and the conditional expectation value $\langle \hat{A} | g|\vec{x}_{f1} \oplus \vec{x}_{f2} \rangle$ is again given by Eq. (192). The entanglement in the resulting state thus crucially depends on the exact properties of $\langle \hat{A} | g|\vec{x}_{f1} \oplus \vec{x}_{f2} \rangle$.

First of all, note that $W_{f1}(\vec{x}_{f1})$ and $W_{f2}(\vec{x}_{f2})$ are generally not pure states and as a consequence $W_{f|A}(\vec{x}_{f1} \oplus \vec{x}_{f2})$ is not a pure state either. Even though the specific structure of the Wigner function makes it a suitable case to apply the methods of Ref. [221], we follow a different route in this Tutorial by focusing on a particular example for which we can assume that $W_{f1}(\vec{x}_{f1})$ and $W_{f2}(\vec{x}_{f2})$ are pure.

Just as in Sec. 2, we concentrate on photon subtraction. To get a conceptual idea of such a setup in this specific scenario, we present two equivalent schemes in (a) and (b) of Fig. 7. Note that the equivalence stems from the fact that the beam splitters that subtract the light from the signal beams to send it to the photodetector are of extremely low reflectivity. In this limit, we can be sure that there is at most one photon in the path and when it is detected, we herald a single-photon-subtracted state. In Fig. 7(a), the combination of this heralding process and the presence of at most one photon avoids that the unmeasured output causes any losses or impurities. Nevertheless, the unmeasured output will practically change the success probability of the heralding process, such that for practical implementations Fig. 7(b) may be the preferential setup. Recall that the Wigner function for a state with a photon subtracted in

a particular mode **b** was given by Eq. (259), which here becomes

$$W^{-}(\vec{x}_{f1} \oplus \vec{x}_{f2}) = \frac{W_{f1}(\vec{x}_{f1}) W_{f2}(\vec{x}_{f2})}{\text{tr}(V_{b} - 1) + \|\tilde{\xi}_{b}\|^{2}} \times (\|B^{T}(1 - V_{f1}^{-1} \oplus V_{f2}^{-1})(\vec{x}_{f1} \oplus \vec{x}_{f2} - \tilde{\xi}_{1} \oplus \tilde{\xi}_{2}) + \tilde{\xi}_{b}\|^{2} + \text{tr}[1 - B^{T}(V_{f1}^{-1} \oplus V_{f2}^{-1})B]). \quad (275)$$

Because we consider a limit where the state is completely transmitted by the beam splitter and only a negligible amount is sent to the photon counter to subtract the photon, we can indeed assume that the state is pure. For simplicity, we also assume that $f_1$ and $f_2$ are single modes. As we did before, we now calculate the reduced state

$$W_{1}^{-}(\vec{x}_{f1}) = \int_{\mathbb{R}^{2}} d\vec{x}_{f2} W^{-}(\vec{x}_{f1} \oplus \vec{x}_{f2}). \quad (276)$$

The integral is rather tedious to evaluate, therefore we immediately jump to the result (see Ref. [153] for an alternative method that circumvents the explicit calculation of integrals):

$$W_{1}^{-}(\vec{x}_{f1}) = \frac{W_{f1}(\vec{x}_{f1})}{\text{tr}(V_{b} - 1) + \|\tilde{\xi}_{b}\|^{2}} (\|B^{T}F_{1}(1 - V_{f1}^{-1})(\vec{x}_{f1} - \tilde{\xi}_{1}) + \tilde{\xi}_{b}\|^{2} + \text{tr}[B^{T}F_{1}(1 - V_{f1}^{-1})F_{1}^{T}B] + \text{tr}[B^{T}F_{2}(V_{f2} - 1)F_{2}^{T}B]), \quad (277)$$

where we introduce the matrices $F_k$, given by

$$F_k = \begin{pmatrix} | & | \\ \vec{f}_k & \Omega \vec{f}_k \end{pmatrix}, \quad (278)$$

such that we can use the properties of the symplectic form $\Omega$ to obtain

$$B^T F_k = \begin{pmatrix} \vec{b}^T \vec{f}_k & \vec{b}^T \Omega \vec{f}_k \\ -\vec{b}^T \Omega \vec{f}_k & \vec{b}^T \vec{f}_k \end{pmatrix}. \quad (279)$$

If mode **b** is orthogonal to mode $f_1$, we find that $B^T F_1 = 0$ such that $W_1^{-}(\vec{x}_{f1}) = W_{f1}(\vec{x}_{f1})$. On the other hand, when mode **b** is exactly the same as $f_1$ we find that $B^T F_1 = 1$ such that the photon is only subtracted there. In this case $W_1^{-}(\vec{x}_{f1})$ is a pure state and no entanglement is created. In this case, one can check that $W_2^{-}(\vec{x}_{f2}) = W_2(\vec{x}_{f2})$.

To create entanglement, we are thus interested in the case where **b** is a superposition of the two modes $f_1$ and $f_2$. To keep things simple, let us assume that $\vec{b} = \cos \theta \vec{f}_1 +$
<figure><img src="image_9.png" /><figcaption>FIG. 7. Entanglement creation through photon subtraction in a superposition of uncorrelated modes $\mathbf{f}_1$ and $\mathbf{f}_2$. (a),(b) Sketches of two equivalent setups to implement a photon subtraction in the mode b, with $\vec{b} = \cos \theta \vec{f}_1 + \sin \theta \vec{f}_2$. In (c),(d), we show the created entanglement, as measure through the Rényi entropy, Eq. (256), for varying values of $\theta$. The initial Gaussian states are pure, with covariance matrices $V_{\mathbf{f}_1} = \text{diag}[s_1, 1/s_1]$ and $V_{\mathbf{f}_2} = \text{diag}[s_2, 1/s_2]$ for modes $\mathbf{f}_1$ and $\mathbf{f}_2$, respectively. The global mean field, i.e., displacement, is varied $\tilde{\xi} = 0$ (solid curves), $\tilde{\xi} = (1, 0, 0, 0)^T$ (dotted curves), and $\tilde{\xi} = (0, 0, 1, 0)^T$ (dotted curves). (c) The particular case where the squeezing is balanced, i.e., $s_1 = s_2 = 2$. (d) An unbalanced example where $s_1 = 4$ and $s_2 = 2$. All squeezing values $s_1$ and $s_2$ are measured in units of vacuum noise.</figcaption></figure>

$\sin \theta \vec{f}_2$. Because the modes $\mathbf{f}_1$ and $\mathbf{f}_2$ are orthogonal, we can use that $\hat{f}_1^\dagger \hat{f}_2 = 0$ and thus we find that $B^\dagger F_1 = \cos \theta \mathbf{1}$ and $B^\dagger F_2 = \sin \theta \mathbf{1}$. Nevertheless, the general expression for $W_1(\vec{x}_{\mathbf{f}_1})$ does not simplify much.

To acquire additional insight, let us now assume that both modes $\mathbf{f}_1$ and $\mathbf{f}_2$ have exactly the same squeezing in the same quadrature:

$$V_{\mathbf{f}_1} = V_{\mathbf{f}_2} = \begin{pmatrix} s & 0 \\ 0 & \frac{1}{s} \end{pmatrix}. \qquad (280)$$

Furthermore, let us assume that there is no mean field, such that $\tilde{\xi} = 0$. In this particular case, we find the expression

$$W_1(x_{\mathbf{f}_1}, p_{\mathbf{f}_1}) = W_1(x_{\mathbf{f}_1}, p_{\mathbf{f}_1}) \times \left[ p_{\mathbf{f}_1}^2 s + \frac{x_{\mathbf{f}_1}^2}{s} + \cos(2\theta) \left( p_{\mathbf{f}_1}^2 s + \frac{x_{\mathbf{f}_1}^2}{s} - 2 \right) \right]. \qquad (281)$$

In particular, it turns out that the purity takes a simple form, such that we can quantify the entanglement for this state as

$$\mathcal{E}_R = \log(2) - \log\left(\frac{\cos(4\theta) + 3}{2}\right). \qquad (282)$$

This shows that the maximal entanglement is reached for $\theta = \pi/4$ and—as expected—the entanglement vanishes when $\theta = 0$ and $\theta = \pi/2$, i.e., when we subtract entirely in either mode $\mathbf{f}_1$ and $\mathbf{f}_2$.

More general settings are shown in Fig. 7, where we show the entanglement creation for unbalanced squeezing, by setting

$$V_{\mathbf{f}_1} = \begin{pmatrix} s_1 & 0 \\ 0 & \frac{1}{s_1} \end{pmatrix}, \quad \text{and} \quad V_{\mathbf{f}_2} = \begin{pmatrix} s_2 & 0 \\ 0 & \frac{1}{s_2} \end{pmatrix}. \qquad (283)$$

We compare the case with $s_1 = s_2$ to the case with $s_1 \neq s_2$ and find that in absence of a mean field one can reach the same maximal amount of entanglement. However, the maximum is attained at a different value of $\theta$ when the squeezing is unbalanced. From Eq. (282) we know that in absence of a mean field, the curve for $s_1 = s_2$ does not depend on the actual value of squeezing.

Figure 7 also shows the effect of an existing mean field, by probing a mean field in mode $\mathbf{f}_1$ with $\tilde{\xi} = (1, 0, 0, 0)^T$ and in mode $\mathbf{f}_2$ with $\tilde{\xi} = (0, 0, 1, 0)^T$. Generally speaking, we observe that the mean field reduces the created entanglement. Nevertheless, the unbalance of squeezing ($s_1 \neq s_2$) also unbalances the effect of the mean field. The higher squeezing in mode $\mathbf{f}_1$ makes the entanglement creation more resilient to displacements, but a mean field in mode $\mathbf{f}_2$ will reduce the maximal attainable amount of entanglement to the same level as in the balanced case [because in both (c) and (d) mode $\mathbf{f}_2$ is squeezed with $s_2 = 2$]. In the presence of a mean field, we also find that unbalanced squeezing shifts the value $\theta$ for which most entanglement is created. In other words, to achieve maximal entanglement upon photon subtraction in two modes with unequal squeezing, one must subtract in an unbalanced superposition of these modes.

Through this example, we showed that entanglement between previously uncorrelated Gaussian states can be created by a non-Gaussian operation. This entanglement has some additional peculiarities. For example, a quick glance at how this procedure affects Eq. (190) shows that we can split the state in a Gaussian, i.e., $W_f(\vec{x}_f)$, and a non-Gaussian part, i.e., $\langle \hat{A} \rangle_{g|\vec{x}_f}/\langle \hat{A} \rangle$. In this case of Eq. (259) the Gaussian part of the state clearly remains fully separable. This means that, in this representation, all entanglement is originating from the non-Gaussian part of the state. Nevertheless, the decomposition, Eq. (190), of the state into a Gaussian and a non-Gaussian part most probably is not unique for mixed states, making it challenging to study such non-Gaussian entanglement in its most general sense.

Yet, common tools that rely on the covariance matrix, such as Refs. [216,217], to characterize entanglement in the photon subtracted states, Eq. (259), are doomed to fail. In Ref. [152] it is explicitly shown that the covariance matrix of a photon-subtracted state is given by the covariance matrix of the initial Gaussian state with a positive matrix added to it. This means that photon subtraction just adds correlated noise to the covariance matrix and if we consider a Gaussian state that has exactly this covariance matrix we can decompose it using Eq. (103). In other
words, when there is no entanglement visible in the covariance matrix of the initial Gaussian state, we do not witness any entanglement based on the covariance matrix of the photon-subtracted state. In this case, the non-Gaussian entanglement is thus genuinely non-Gaussian in the sense that it cannot be detected through Gaussian witnesses. Hence, rather than decomposing the states in a Gaussian and non-Gaussian part, as was done in Eq. (190), it may be more fruitful to define non-Gaussian entanglement as any entanglement that cannot be witnessed based solely on the covariance matrix of the state. This approach also offers a natural connection to the framework of Gaussian passivity on quantum thermodynamics [245].

Another peculiarity that was presented in Refs. [152, 153] is the intrinsic nature of this non-Gaussian entanglement. When we transform the system into a different mode basis, there will still be entanglement in the system. The entanglement is said to be intrinsic because the state is entangled in every possible mode basis. As we saw in Eq. (249) Gaussian entanglement is never intrinsic as there always exists a basis in which a Gaussian state is separable.

Figure 7(b) gives a rather interesting approach to understanding the intrinsic nature of non-Gaussian entanglement. In this sketch, the second beam splitter is intended to undo the superposition $\theta$ and return to the initial mode basis with modes $f_1$ and $f_2$. Changing this beam splitter thus implies a basis change. If we remove this beam splitter entirely, we find ourselves in the entanglement distillation scenario of Sec. 2. In this case, the photon subtraction is fully local, but it happens on a state with Gaussian entanglement. The photon subtraction can then increase the Rényi entanglement by a maximal amount of $\log 2$. When we change to a mode basis where there is no Gaussian entanglement and the entanglement is created through a nonlocal photon subtraction, we create a maximal amount of Rényi entanglement given by $\log 2$. Changing the mode basis in a different way will combine the physics of these two extreme cases such that there will always be entanglement, regardless of the basis.

Extending these ideas to more general non-Gaussian operations on more general Gaussian mixed states is a hard and currently open problem. This reflects the general status of entanglement theory in CV systems: we lack a structured theoretical understanding of this phenomenon and as a consequence we also lack good tools to detect it.

## D. Non-Gaussianity and Bell inequalities

In this final subsection of our study of quantum correlations in non-Gaussian states, we study Bell inequalities. First of all, we argue that it is impossible to violate Bell inequalities when both of the states and all the measurements involved can be described by positive Wigner functions. Then, we show that the Wigner function of the

state can itself be used to formulate a Bell inequality when we allow for nonpositive Wigner functions.

The general setup for studying nonlocality in CV revolves around a multimode state with Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ defined on a phase space $\mathbb{R}^{2m} = \mathbb{R}^{2l} \oplus \mathbb{R}^{2l'}$. Bell nonlocality entails that some measurements on this state cannot be described by a local hidden variable model of the type, Eq. (238). In a quantum framework, the local measurements with POVM elements $\{\hat{A}_j\}$ (on the modes in $f$) and $\{\hat{B}_j\}$ (on the modes in $g$) can also be described by Wigner functions $W_{A_j}(\vec{x}_f)$ and $W_{B_j}(\vec{x}_g)$. Because we are dealing with a POVM, we find that

$$ (4\pi)^l \sum_j W_{A_j}(\vec{x}_f) = (4\pi)^l \sum_j W_{B_j}(\vec{x}_g) = 1. \quad (284) $$

Note that this equality holds for all possible coordinates $\vec{x}_f$ and $\vec{x}_g$. Here we assume that the measurement outcomes $A_j$ and $B_j$ are discrete, but by correctly defining resolutions of the identity we can also deal with more general probability distributions, e.g., homodyne measurements.

The probability to get the joint measurement result $(A_j, B_k)$ is given by

$$ P(A_j, B_k) = (4\pi)^m \int_{\mathbb{R}^{2l}} \int_{\mathbb{R}^{2l'}} d\vec{x}_f d\vec{x}_g W(\vec{x}_f \oplus \vec{x}_g) W_{A_j}(\vec{x}_f) W_{B_k}(\vec{x}_g). \quad (285) $$

Now let us assume that all these Wigner functions are positive. Because they are normalized, this implies that $W(\vec{x}_f \oplus \vec{x}_g)$ is a probability distribution on the entire phase space $\mathbb{R}^{2m}$, and $W_{A_j}(\vec{x}_f)$ and $W_{B_j}(\vec{x}_g)$ are probability distributions on the reduced phase spaces $\mathbb{R}^{2l}$ and $\mathbb{R}^{2l'}$, respectively. However, the model, Eq. (238), does not require probability distributions on phase space, but rather on the possible measurement outcomes.

This is where Eq. (284) comes into play. Because $W_{A_j}(\vec{x}_f)$ and $W_{B_j}(\vec{x}_g)$ are positive, more than just treat them as probability distributions in phase space we can also consider $P_{\vec{x}_f}(A_j) = (4\pi)^l W_{A_j}(\vec{x}_f)$ and $P_{\vec{x}_g}(B_j) = (4\pi)^l W_{B_j}(\vec{x}_g)$ as the probability of getting the measurement outcomes $A_j$ and $B_k$, respectively. Because of Eq. (284) we find that these probabilities are correctly normalized

$$ \sum_j P_{\vec{x}_f}(A_j) = \sum_j P_{\vec{x}_g}(B_j) = 1, \quad (286) $$

and because the Wigner functions are positive, we also find that $P_{\vec{x}_f}(A_j), P_{\vec{x}_g}(B_j) \ge 0$. Note that the phase-space coordinates $\vec{x}_f$ and $\vec{x}_g$ are no longer treated as the variable, but rather as a label. The set $\{P_{\vec{x}_f}(A_j) | \vec{x}_f \in \mathbb{R}^{2l}\}$ denotes a family of different probability distributions on
the space of measurement outcomes {$A_1, A_2, \dots$}. The set $\{P_{\vec{x}_g}(B_j) | \vec{x}_g \in \mathbb{R}^{2l}\}$ can be interpreted analogously.

We can thus recast Eq. (285) in the following form:

$$P(A_j, B_k) = \int d\vec{x}_f d\vec{x}_g W(\vec{x}_f \oplus \vec{x}_g) P_{\vec{x}_f}(A_j) P_{\vec{x}_g}(B_k). \quad (287)$$

Because $W(\vec{x}_f \oplus \vec{x}_g)$ is a positive and normalized Wigner function, it is a joint probability distribution on the coordinates $\vec{x}_f$ and $\vec{x}_g$. These coordinates label families of probability distributions $\{P_{\vec{x}_f}(A_j) | \vec{x}_f \in \mathbb{R}^{2l}\}$ and $\{P_{\vec{x}_g}(B_j) | \vec{x}_g \in \mathbb{R}^{2l}\}$ for the measurement outcomes. The expression (287) is thus fully consistent with Bell's local hidden variable model, Eq. (238). As a consequence, we cannot violate any Bell inequalities when the system is prepared in a state with a positive Wigner function and when we only have access to POVM that have Wigner representations with positive Wigner functions.

Let us emphasize that there is generally no reason to assume that the probabilities $P_{\vec{x}_f}(A_j)$ and $P_{\vec{x}_g}(B_j)$ are also consistent with quantum mechanics. In other words, there is not necessarily any state $\hat{\rho}$ such that $P_{\vec{x}_f}(A_j) = \text{tr}[\hat{\rho} \hat{A}_j]$. However, because we are dealing with Bell nonlocality, we do not need this to be the case, since Eq. (238) allows for arbitrary local probability distributions.

To make a long story short, we have shown that Wigner negativity is necessary for witnessing Bell nonlocality. The interested reader can consult works such as Ref. [64] that relate Wigner negativity to the more general concept of quantum contextuality. However, the topic of contextuality in CV systems is still a matter of scientific debate [246].

There has been a significant body of work about the violation of Bell inequalities in CV setups [247–249]. It is evident that this is an arduous task once one approaches a realistic experimental setting [250]. Here, we focus on one particular suggestion to test Bell nonlocality based on a state’s Wigner function [251,252].

The starting point of this approach is the CHSH inequality

$$|\langle\hat{X}\hat{Y}\rangle - \langle\hat{X}'\hat{Y}'\rangle + \langle\hat{X}'\hat{Y}\rangle + \langle\hat{X}'\hat{Y}'\rangle| \le 2.$$

As we discussed in Sec. 4, this inequality relies on some assumptions for the observables $X, X', Y$, and $Y'$. In particular, we must assume that the measurement outcomes are either $-1$ or $1$. In a CV setting, where we generally deal with a continuum of possible measurement outcomes, this seems like a serious constraint. Nevertheless, we have already encountered some natural examples during this Tutorial. For example, photon counters yield a discrete number of possible measurement outcomes. Here, we choose a related observable that takes us all the way back to Sec. B, where we encountered the observable

$$\hat{\Delta}(\vec{x}) = \hat{D}(-\vec{x})(-1)^{\hat{N}} \hat{D}(\vec{x}).$$

This displaced parity operator has a rich structure, but when it comes to actual measurement outcomes is will return either $-1$ or $1$. This means that we can choose $X, X', Y$, and $Y'$ to be parity operators. First of all, let us note that

$$\hat{\Delta}(\vec{x}_f \oplus \vec{x}_g) = \hat{\Delta}(\vec{x}_f) \otimes \hat{\Delta}(\vec{x}_g). \quad (288)$$

To see this, one can first show that $(-1)^{\hat{N}_m} = (-1)^{\hat{N}_l+\hat{N}'_l} = (-1)^{\hat{N}_l} \otimes (-1)^{\hat{N}'_l}$ and subsequently use $\hat{D}(\vec{x}_f \oplus \vec{x}_g) = \hat{D}(\vec{x}_f) \otimes \hat{D}(\vec{x}_g)$ (displacements in different modes are independent from each other).

Now we can identify the observables as follows:

$$\begin{aligned} X &= \hat{\Delta}(\vec{x}_f), & X' &= \hat{\Delta}(\vec{x}'_f), \\ Y &= \hat{\Delta}(\vec{x}_g), & Y' &= \hat{\Delta}(\vec{x}'_g), \end{aligned} \quad (289)$$

and therefore the CHSH inequality is transformed into

$$|\langle\hat{\Delta}(\vec{x}_f \oplus \vec{x}_g)\rangle - \langle\hat{\Delta}(\vec{x}_f \oplus \vec{x}'_g)\rangle + \langle\hat{\Delta}(\vec{x}'_f \oplus \vec{x}_g)\rangle + \langle\hat{\Delta}(\vec{x}'_f \oplus \vec{x}'_g)\rangle| \le 2. \quad (290)$$

As a next step, we use Eq. (78) to write

$$\langle\hat{\Delta}(\vec{x}_f \oplus \vec{x}_g)\rangle = (2\pi)^m W(\vec{x}_f \oplus \vec{x}_g), \quad (291)$$

such that the inequality (290) can be recast as

$$|W(\vec{x}_f \oplus \vec{x}_g) - W(\vec{x}_f \oplus \vec{x}'_g) + W(\vec{x}'_f \oplus \vec{x}_g) + W(\vec{x}'_f \oplus \vec{x}'_g)| \le \frac{2}{(2\pi)^m}. \quad (292)$$

Any state with a Wigner function that violates this inequality for some choice of coordinates $\vec{x}_f, \vec{x}'_f, \vec{x}_g$, and $\vec{x}'_g$ possesses some form of Bell nonlocality. In Refs. [251,252] it is argued that the inequality (292) can be violated by sending a single photon through a beam splitter, but also by an EPR state. The fact that a Gaussian Wigner function suffices to violate Eq. (292) sometimes comes as a surprise, because we previously argued that one needs Wigner negativity to violate Bell inequalities. The reason why one can detect Bell nonlocality with this inequality even when the Wigner function is positive stems from our choice of observable $\hat{\Delta}(\vec{x}_f \oplus \vec{x}_g)$. The POVM elements that correspond to the measurement outcomes $1$ and $-1$ have Wigner functions that are strongly Wigner negative. As a consequence the necessary Wigner negativity is baked into Eq. (292) by construction.

In practice, the inequality (292) is highly sensitive to impurities and can often be hard to violate with experimentally reconstructed Wigner functions. The hunt for good new techniques to show Bell nonlocality in CV systems
is therefore still open. However, this subsection clearly showed us that Wigner negativity is necessary to observe one of the most exotic features in quantum physics. This negativity might be baked into the state, but it could just as well be induced by measurements. The conditional methods of Sec. B also highlight this duality, where Wigner negativity in the measurement is used to induce Wigner negativity in the state. It should come as no surprise that Wigner negativity is also a necessary ingredient for the most exotic quantum protocols. However, it should also be highlighted that even with a little extra trust, it is possible to design protocols that do not require Wigner negativity to witness quantum correlations [201]. In the next section, we discuss its importance for reaching a quantum advantage with CV systems.

## VI. NON-GAUSSIAN QUANTUM ADVANTAGES

It has been long known that systems that are entirely built with Gaussian building blocks are easy to simulate [25]. It should perhaps not come as a surprise that efficient numerical tools exist to sample numbers from a multivariate Gaussian distribution. The discrete variable analog of this result comes across as less intuitive and goes by the name “Gottesman-Knill theorem” [253]. Yet, it turns out that something stronger than mere non-Gaussian elements is required to render a system hard to simulate.

In Sec. D, we encountered the power of Wigner negativity by realizing that it is a necessary requirement for Bell nonlocality. This connection between Wigner negativity and the most exotic types of quantum correlations shows us that Wigner negativity is key to giving CV systems their most prominent quantum features. It is then perhaps not a surprise that such Wigner negativity is also a necessary requirement for implementing any type of protocol that cannot be efficiently simulated by a classical computer [33,254,255]. We thus start our discussion of quantum advantages by explaining the result of Ref. [33]. To show the necessity of Wigner negativity, we show an explicit simulation algorithm for general quantum protocols without Wigner negativity.

Any quantum protocol ultimately relies on the measurement of a certain set of measurement operators {$\hat{E}_j$} (typically a POVM) of a system prepared on a state $\hat{\rho}$. In a Wigner function formalism, we then find

$$p_j = (4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x} \, W_{E_j}(\vec{x}) W(\vec{x}). \quad (293)$$

Furthermore, the fact that the set {$\hat{E}_j$} forms a POVM implies that

$$(4\pi)^m \sum_j W_{E_j}(\vec{x}) = 1. \quad (294)$$

As we already discussed in Sec. D, surrounding Eq. (284), it is crucial that the normalization condition (294) holds for any phase-space coordinate $\vec{x}$. In the present context, we want to show that there is an efficient method for a classical device to sample values from the probability distribution {$p_j$} when all involved Wigner functions are positive.

Let us start by assuming that the Wigner functions that describe the POVM elements are all positive. When combined with the POVM condition (294), this implies that we can identify a set of probabilities $P_{\vec{x}}(e_j) = (4\pi)^m W_{E_j}(\vec{x})$ as the probability to obtain the measurement outcome $e_j$, associated with the POVM element $\hat{E}_j$. These $P_{\vec{x}}(e_j)$ depend on a parameter $\vec{x}$, we can thus form a family of probability distributions {$P_{\vec{x}}(e_j) | \vec{x} \in \mathbb{R}^{2m}$} that describe the probability of obtaining the different results $e_j$, depending on a chosen phase-space point. The normalization condition (294) now states that $\sum_j P_{\vec{x}}(e_j) = 1$ for all $\vec{x}$. Let us emphasize that this family of probabilities would not be well defined if $W_{E_j}(\vec{x})$ were not positive Wigner functions, as some of the probabilities would be negative.

Going back to the initial Eq. (293), we now find that

$$p_j = \int_{\mathbb{R}^{2m}} d\vec{x} \, P_{\vec{x}}(e_j) W(\vec{x}). \quad (295)$$

To find the actual probability of getting the jth outcome is thus given by “averaging” the probabilities $P_{\vec{x}}(e_j)$ over the different phase-space coordinates. Generally speaking, this is not a real average, unless the Wigner function $W(\vec{x})$ of the state is an actual probability distribution on phase space. The latter is exactly the case when $W(\vec{x})$ is positive. Then, we can simply think of the probability $p_j$ for obtaining event $e_j$ as $p_j = \mathbb{E}_W[P_{\vec{x}}(e_j)]$, where $\mathbb{E}_W$ is the expectation value over the probability distribution $W(\vec{x})$.

Hence, when all Wigner functions are positive, the algorithm to simulate our relevant quantum process can simply be expressed by the following steps:

1. Sample a phase-space coordinate $\vec{x}$ from the probability distribution $W(\vec{x})$.

2. Construct the probability distribution $P_{\vec{x}}(e_j)$ for the sampled value $\vec{x}$.

3. Sample an outcome $e_j$ from the probability distribution $P_{\vec{x}}(e_j)$.

Even though this is the general idea behind our sampling protocol, there are some major hidden assumptions. First, we assume here that the Wigner function for the state and the measurement are known. Furthermore, we also assume that we can simply sample points from any distribution on phase space and from any distribution of measurement outcomes $P_{\vec{x}}(e_j)$. In particular, for the sampling aspects it is not at all clear that these are reasonable assumptions to make. Standard sampling protocols for multivariate probability distributions tend to get highly inefficient once the
probability distributions become too exotic such that it is dangerous to assume that we can “just sample.”

To address this point Ref. [33] makes more assumptions on the exact setup we are trying to simulate. First, we assume that the detection is done by a series of single-mode detectors, such that our label j now become a tuple $j = (j_1, j_2, \dots, j_m)$ where $j_k$ denotes the outcome $e_{j_k}^{(k)}$ for the detector on the kth mode. We can thus write the POVM element as $\hat{E}_j = \hat{E}_{j_1}^{(1)} \otimes \dots \otimes \hat{E}_{j_m}^{(m)}$, such that

$$W_{E_j}(\vec{x}) = W_{E_{j_1}^{(1)}}(x_1, p_1) W_{E_{j_2}^{(2)}}(x_2, p_2) \dots W_{E_{j_m}^{(m)}}(x_m, p_m). \quad (296)$$

We assume that each detector has been accurately calibrated, such that all the individual Wigner functions are known. This implies that for a given point in phase space $\vec{x} = (x_1, p_1, \dots, x_m, p_m)^T$, we can simply evaluate the probabilities for each detector to produce a certain outcome. Thus, we calculate the probability distributions $P_{(x_k, p_k)}(e_j^{(k)}) = 4\pi W_{E_j^{(k)}}(x_k, p_k)$ for all the possible measurement outcomes for that specific mode. We assume that sampling outcomes $e_j^{(k)}$ from these probability distributions $P_{(x_k, p_k)}(e_j^{(k)})$ is a feasible task. For typical detectors in quantum optics experiments this is a very reasonable assumption.

The Wigner function $W(\vec{x})$ that describes the state is more subtle as it also includes all correlations between modes. If the state is Gaussian, a measurement of the covariance matrix would be sufficient to know the full Wigner function. Because of its Gaussian features, there are efficient tools to directly sample phase-space points. Yet, for more general non-Gaussian positive Wigner functions this sampling may be much harder. Therefore, Ref. [33] makes an essential assumption: it assumes that we know a protocol that combines local operations to design the state $W(\vec{x})$ from a known initial state with no correlations between the modes. The notion of “locality” should here be understood in the sense of acting on a small set of modes while leaving the others fully untouched. These local operations are also supposed to be represented by positive Wigner functions, which depend only on the phase-space coordinates of the subset of modes on which they act.

Generally speaking, such Wigner positive operations $\Xi: \mathcal{H}^{\text{in}} \to \mathcal{H}^{\text{out}}$ map a state $\hat{\rho}$ to a new state $\Xi[\hat{\rho}]$. In Ref. [33], the Choi representation [256,257] is used to represent $\Xi$ as a state on a larger Hilbert space $\mathcal{H}^{\text{in}} \otimes \mathcal{H}^{\text{out}}$. This becomes particularly appealing when we go to a Wigner representation, where the Choi representation of $\Xi$ is given by a Wigner function $W_{\Xi}(\vec{x}^{\text{in}} \oplus \vec{x}^{\text{out}})$. The action of $\Xi$ on

a state with Wigner function $W(\vec{x}^{\text{in}})$ is then given by

$$W_{\text{out}}(\vec{x}^{\text{out}}) = (4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x}^{\text{in}} W_{\Xi}(\vec{x}^{\text{in}} \oplus T\vec{x}^{\text{out}}) W(\vec{x}^{\text{in}}), \quad (297)$$

where m is the number of modes of the input state. For technical reasons, we must include the transposition operator T, Eq. (257), in the action of the channel. Because this operation must be trace preserving, we on top get the property that

$$(4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x}^{\text{out}} W_{\Xi}(\vec{x}^{\text{in}} \oplus T\vec{x}^{\text{out}}) = 1. \quad (298)$$

When we now assume that the operation $\Xi$ has a Wigner-Choi representation $W_{\Xi}(\vec{x}^{\text{in}} \oplus \vec{x}^{\text{out}})$, which is a positive function, it immediately follows that the operation $\Xi$ turns a Wigner positive initial state $W(\vec{x}^{\text{in}})$ into a Wigner positive output state $W_{\text{out}}(\vec{x}^{\text{out}})$.

It is useful to note that such operations, Eq. (297), can be trivially embedded in a larger space. Let us assume that we consider a state $W(\vec{x}_f \oplus \vec{x}_g)$, we can simply let the operation act on the modes g by taking

$$W_{\text{out}}(\vec{x}_f \oplus \vec{x}_g^{\text{out}}) = (4\pi)^{l'} \int_{\mathbb{R}^{2l}} d\vec{x}_g^{\text{in}} W_{\Xi}(\vec{x}_g^{\text{in}} \oplus T\vec{x}_g^{\text{out}}) \times W(\vec{x}_f \oplus \vec{x}_g^{\text{in}}). \quad (299)$$

Notationally, this may seem a little complicated, but, in essence, we just carry out the integration over a subset of the full phase space. We call these operations local Wigner positive operations.

In our simulation protocol, we thus assume that $W(\vec{x})$ is created by a series of such local Wigner positive operations of $\Xi_1, \dots, \Xi_t$ on a noncorrelated input state $W_{\text{in}}(\vec{x}_{\text{in}}) = W_{\text{in}}^{(1)}(x_1, p_1) W_{\text{in}}^{(2)}(x_2, p_2) \dots W_{\text{in}}^{(m)}(x_m, p_m)$.

$$W(\vec{x}) = (4\pi)^{ml} \int_{\mathbb{R}^{2m}} d\vec{x}_t \dots \int_{\mathbb{R}^{2l}} d\vec{x}_l W_{\Xi_t}(\vec{x}_t \oplus T\vec{x}) \dots \\ \times W_{\Xi_2}(\vec{x}_1 \oplus T\vec{x}_2) W_{\Xi_1}(\vec{x}_{\text{in}} \oplus T\vec{x}_1) \\ \times W_{\text{in}}(\vec{x}_{\text{in}}). \quad (300)$$

We assume on top that each operation is local over a small number of modes $l \ll m$. To model this with Eq. (299), it suffices to split $\vec{x}_{t_{k-1}} = \vec{x}_{t_{k-1}}^l \oplus \vec{x}_{t_{k-1}}^{l'}$ and $\vec{x}_{t_k} = \vec{x}_{t_k}^l \oplus \vec{x}_{t_k}^{l'}$, such that

$$(4\pi)^m W_{\Xi_{t_k}}(\vec{x}_{t_{k-1}} \oplus T\vec{x}_{t_k}) = (4\pi)^l W_{\Xi_{t_k}}(\vec{x}_{t_{k-1}}^l \oplus T\vec{x}_{t_k}^l) \delta(\vec{x}_{t_{k-1}}^l - \vec{x}_{t_k}^l). \quad (301)$$

Even though the notation is complicated, it simply describes that we act on an l-mode subspace with the operation $\Xi_{t_k}$ and leave the other l' modes untouched.
The normalization condition (298) now has an important consequence, since it allows us to identify a probability distribution on phase space $P_{\vec{x}_{t_{k-1}}}(\vec{x}_{t_k}) = (4\pi)^m W_{\Xi_{t_k}}(\vec{x}_{t_{k-1}} \oplus T\vec{x}_{t_k})$. It gives us the probability of choosing a phase space value $\vec{x}_{t_k}$, given that we know $\vec{x}_{t_{k-1}}$. Because the operations are local, Eq. (301) allows us to keep most of the phase-space coordinates constant from step to step. Furthermore, the first step is simple. Every pair $(x_k, p_k)$ of the initial coordinate $\vec{x}_{in}$ can be sampled independently because $W_{in}(\vec{x}_{in})$ factorizes. This now gives us the following new algorithm:

1. Take the initial Wigner function $W_{in}(\vec{x}^{in}) = W_{in}^{(1)}(x_1, p_1)W_{in}^{(2)}(x_2, p_2) \dots W_{in}^{(m)}(x_m, p_m)$ and sample a pair $(x_k, p_k)$ from every single-mode probability distribution $W_{in}^{(k)}(x_k, p_k)$. Put all these pairs together to obtain $\vec{x}^{in} = (x_1, p_1, \dots, x_m, p_m)$.

2. Update the coordinate by sampling new coordinates based on $P_{\vec{x}_{t_{k-1}}}(\vec{x}_{t_k}) = (4\pi)^m W_{\Xi_{t_k}}(\vec{x}_{t_{k-1}} \oplus T\vec{x}_{t_k})$. Because $\Xi_{t_k}$ are local operations, it suffices to only locally update coordinates. Let us make this clear through an example. Say we have $\vec{x}_{t_{k-1}} = (x_1^{(k-1)}, p_1^{(k-1)}, \dots, x_m^{(k-1)}, p_m^{(k-1)})^T$ and operation $\Xi_{t_k}$ acts locally on modes with labels 2, 5, and 7. Take $\vec{x}_{t_{k-1}}^l = (x_2^{(k-1)}, p_2^{(k-1)}, x_5^{(k-1)}, p_5^{(k-1)}, x_7^{(k-1)}, p_7^{(k-1)})$ and use it to evaluate $P_{\vec{x}_{t_{k-1}}}(\vec{x}_{t_k}^l) = (4\pi)^l W_{\Xi_{t_k}}(\vec{x}_{t_{k-1}}^l \oplus T\vec{x}_{t_k}^l)$. Now sample a new vector $\vec{x}_{t_k}^l = (x_2^{(k)}, p_2^{(k)}, x_5^{(k)}, p_5^{(k)}, x_7^{(k)}, p_7^{(k)})$ from this probability distribution. Then construct the new vector $\vec{x}_{t_k}$ by taking $\vec{x}_{t_{k-1}}$ and updating the coordinates associated to modes 2, 5, and 7 to the newly sampled coordinates

3. After the operations $\Xi_1, \dots, \Xi_t$ have been implemented by updating the phase-space coordinate, take the final phase-space coordinate $\vec{x} = (x_1, p_1, \dots, x_m, p_m)^T$ and the Wigner function describing the detectors $W_{E_j}(\vec{x}) = W_{E_j^{(1)}}(x_1, p_1)W_{E_j^{(2)}}(x_2, p_2) \dots W_{E_j^{(m)}}(x_m, p_m)$. For each detector $k$, use the phase-space coordinate $\vec{x}$ to generate the probability distribution $P_{(x_k, p_k)}(e_j^{(k)}) = W_{E_j^{(k)}}(x_k, p_k)$.

4. Sample an outcome $e_j^{(k)}$ from the distribution $P_{(x_k, p_k)}(e_j^{(k)})$ for every detector.

Sampling the final phase-space coordinate $\vec{x}$ by using a Monte-Carlo-style update rule is time consuming, but if the operations are local it can be done efficiently. This procedure implicitly assumes that we do not just know the state we are sampling from, but that we know the circuit of local operations that is used to create the state from local resources. Ultimately, when one considers the circuit representation of quantum algorithms, this is also how a quantum algorithm works. For example, Sec. A

exactly shows that any unitary CV circuit can be built with single- and two-mode gates. The algorithm outlined in this section shows that we can efficiently simulate any protocol where the local input state, the circuit's operations, and the measurements are described by positive Wigner functions.

One may wonder whether any positive Wigner function $W(\vec{x})$ can be constructed through such a circuit and, if so, whether there is an efficient way to design such a circuit when we know the Wigner function. If we assume that not only the state $W(\vec{x})$ but also all its marginals are known, it is possible to construct a stepwise sampling procedure through the chain rule of probability theory:

$$ W(\vec{x}) = W(x_m, p_m | x_{m-1}, p_{m-1}, \dots, x_1, p_1) \times \dots \\ W(x_3, p_3 | x_2, p_2, x_1, p_1)W(x_2, p_2 | x_1, p_1)W(x_1, p_1). \tag{302} $$

This process effectively executes a type of random walk with memory. In each step of this walk, we then sample the phase-space coordinates for one mode. Nevertheless, this process only works when we have access to all these conditional probabilities, which practically implies having access to all the marginals of the distribution. In practical setups, this will often not be the case. Nevertheless, it is quickly seen that this setup can be efficiently used to sample from Gaussian Wigner functions where these conditional distributions have a particularly simple form.

Thus, we have shown that it is impossible to obtain a quantum computational advantage by using only local states, measurements, and operations with positive Wigner functions. This means that Wigner negativity is necessary to reach a quantum advantage in such setups. However, Wigner negativity is certainly not sufficient since there are many setups of quantum systems that involve negative Wigner functions that can be efficiently simulated [258]. It is thus interesting to take the opposite approach and explore a setup that is known to lead to a quantum advantage. In the spirit of CV setups, the most logical choice for such a discussion is Gaussian boson sampling [259]. In the literature, this setup has been studied mainly from the point of view of complexity theory [260,261], but here we rather focus on its physical building blocks.

Boson sampling [262] is a problem in which one injects a set of $N$ bosons (generally photons) into an $m$-mode interferometer. On the output ports of this interferometer, photodetectors are mounted to count the particles at the output. Simulating this type of quantum Galton board is a computationally hard task, implying that a quantum advantage could be reached by implementing the setup in a quantum optics experiment. On the other hand, it turns out that the required number of photons to implement such an experiment is also hard to come by. This was the motivation for developing a new approach, where the input photons are replaced by squeezed states that are injected
in each of the interferometer inputs. Because an interferometer, built out of phase shifters and beam splitters, is a Gaussian transformation the output state will remain Gaussian. We can thus effectively say that we are sampling photons from a state with Wigner function $W_G(\vec{x})$. In addition, there is no mean field in the setup such that the entire state is characterized by its covariance matrix $V$.

When we assume that the detectors resolve photon numbers, the probability to detect a string of counts $n = (n_1, \dots, n_m)$ can we written as

$$P(\mathbf{n}) = (4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x} \, W_{\mathbf{n}}(\vec{x}) W_G(\vec{x}). \quad (303)$$

We can then use Eq. (118) to write

$$W_{\mathbf{n}}(\vec{x}) = W_{n_1}(x_1, p_1) \dots W_{n_m}(x_m, p_m). \quad (304)$$

Even though the integral (303) is hard to compute, it is insightful in the light of Eq. (295) and our discussion regarding the necessity of Wigner negativity. Indeed, we see immediately that the detectors form a crucial element in rendering the setup hard to simulate. The same holds when we replace the number-resolving detectors with their on-off counterparts [263] such that $n_k = \{0, 1\}$ and the Wigner functions are given by $W_{n_k}(x_k, p_k) = \{1 - 2 \exp[-(x_k^2 + p_k^2)/2]\}/(4\pi)$.

When we stick with number-resolving detectors that project on Fock states, it is practical to reformulate the problem in terms of $P$ functions and $Q$ functions, such that

$$P(\mathbf{n}) = \int_{\mathbb{R}^{2m}} d\vec{x} \, P_{\mathbf{n}}(\vec{x}) Q_G(\vec{x}). \quad (305)$$

For the detailed calculation, we refer to Ref. [259]. It turns out that the probabilities $P(\mathbf{n})$ can be expressed in terms of the Hafnian of a matrix [264], which establishes a connection to the problem of finding perfect matchings in graph theory. This connection has led to several suggested applications for Gaussian Boson sampling [265–267].

In the light of this Tutorial, the most interesting application of Gaussian Boson sampling is its potential role in quantum state engineering [205]. When only a subset of modes are measured, we can see Gaussian Boson sampling as a generalization of photon subtraction (and even as a generalization of “generalized photon subtraction” [157]). The idea is reasonably simply explained in the light of Sec. B: when we split the system in two parts $\mathbb{R}^{2m} = \mathbb{R}^{2l} \oplus \mathbb{R}^{2l'}$, such that the Gaussian state that comes out of the interferometer now takes the form $W_G(\vec{x}_f \oplus \vec{x}_g)$, we can postselect on a measurement outcome $\mathbf{n} = (n_1, \dots, n_l)$ for the second subsystem. We thus project on a state $W_{\mathbf{n}}(\vec{x}_g)$, which is a product of $l'$ Fock states, and from Eq. (190), we obtain that the conditional state on the remaining modes is given

by

$$W_{f|\mathbf{n}}(\vec{x}_f) = \frac{\langle [ |n_1\rangle \langle n_1| \otimes \dots \otimes |n_{l'}\rangle \langle n_{l'}| ] \rangle_{g|\vec{x}_f} W_f(\vec{x}_f)}{\langle [ |n_1\rangle \langle n_1| \otimes \dots \otimes |n_{l'}\rangle \langle n_{l'}| ] \rangle}, \quad (306)$$

with $W_f(\vec{x}_f)$ defined by Eq. (182). From Eq. (192) we recall the expression

$$\langle [ |n_1\rangle \langle n_1| \otimes \dots \otimes |n_{l'}\rangle \langle n_{l'}| ] \rangle_{g|\vec{x}_f} = (4\pi)^{l'} \int_{\mathbb{R}^{2l'}} d\vec{x}_g W_{\mathbf{n}}(\vec{x}_g) W_G(\vec{x}_g | \vec{x}_f), \quad (307)$$

and because the initial state $W_G(\vec{x}_f \oplus \vec{x}_g)$ is Gaussian, we find that the conditional probability distribution $W_G(\vec{x}_g | \vec{x}_f)$ is given by Eq. (193). Ironically, to evaluate $\langle [ |n_1\rangle \langle n_1| \otimes \dots \otimes |n_{l'}\rangle \langle n_{l'}| ] \rangle_{g|\vec{x}_f}$ and $\langle [ |n_1\rangle \langle n_1| \otimes \dots \otimes |n_{l'}\rangle \langle n_{l'}| ] \rangle$ we must essentially solve the same hard problem as for the implementation of Gaussian Boson sampling itself. Therefore, the exact description of the resulting states is generally complicated.

Nevertheless, in idealized scenarios, even small Gaussian Boson sampling circuits can be used to prepare interesting non-Gaussian states [205]. In particular, the capacity of Gaussian Boson sampling to produce GKP states has taken up a prominent place in a recent blueprint for photonic quantum computation [24]. Furthermore, the results in Sec. C suggest that states created by performing Gaussian Boson sampling on a subset of modes can have additional non-Gaussian entanglement. Yet, to be able to use this procedure to produce highly resourceful Wigner negative states, Sec. 2 highlights that the initial Gaussian state needs to be such that the modes in $f$ can steer the modes in $g$. This condition can be seen as a basic quality requirement for the Gaussian Boson samplers that are used in Ref. [24].

Finally, the experimental imperfections are also detrimental for the quantum advantage that is produced in Gaussian Boson sampling. Clearly, when the Gaussian state $W_G(\vec{x})$ can be written as a Gaussian mixture of coherent states (meaning that no mode basis exists in which the quadrature noise is below vacuum noise), the sampling can be simulated efficiently. Because multimode coherent states are always just a tensor product of single-mode coherent states, it suffices to sample a coherent state from the mixture, calculate all the individual probabilities for the output detectors, and sample independent detector outputs according to these probabilities. The presence of entanglement in the Gaussian state from which we sample is thus crucial. In addition, detector efficiencies must be sufficiently high such that their Wigner functions remain nonpositive, otherwise the protocol of Ref. [33] renders the setup easy to simulate (as explained in the first part
of this section). A more thorough analysis of how dif-
ferent experimental imperfections render Gaussian Boson
sampling easier to simulate can be found in Ref. [268].

There are clearly still many aspects of the relation
between non-Gaussian features of quantum states on the
one hand, and the ability to achieve a quantum computa-
tional advantage on the other hand, that are not yet fully
understood. The Gaussian Boson sampling setup clearly
emphasizes the importance of entanglement in combina-
tion with Wigner negativity. Furthermore, there is the
implicit fact that a simulation scheme such as Ref. [33]
requires knowledge of the actual circuit of local operations
that was used to create the state. It does make sense to
assume that we actually have some ideas of the quantum
protocol that we are attempting to simulate, but yet one
may wonder whether there could be a reasonable setting
(in the sense that we are actually implementing a well-
controlled protocol) in which the assumptions of Ref. [33]
do not hold. This clearly shows that many fundamental the-
oretical aspects of CV quantum computation remain to be
uncovered.

VII. EXPERIMENTAL REALIZATIONS

Now that we have provided an overview of some theo-
retical aspects of non-Gaussian quantum states, we inter-
pret the “where to find them” part of the title in a very lit-
eral sense. Non-Gaussian states are generally rather fragile,
as one should expect from quantum central limit theorem
and the fact that thermal states in free bosonic theories are
Gaussian. Producing and analyzing non-Gaussian states in
a laboratory setting is indeed challenging, but nevertheless
it has been done numerous times. Our main focus in Sec.
A is quantum optics, which is the historical testbed for CV
quantum physics. However, in recent years there has been
increased attention for CV approaches in other setting such
as optomechanics, superconducting circuits, and trapped
ions.

A. Quantum optics experiments

*This section provides an overview of some of the most
important milestones in the generation of non-Gaussian
states in optics. For more details, we refer the reader to*
*a specialized review [81].*

Historically, one might argue that the first experimental
realizations of non-Gaussian states in optical setups relied
on sufficiently sensitive photon detectors. Initial demon-
strations primarily used photoemission of atoms [269,270],
which are prepared in excited states (e.g., by electron bom-
bardment) or via resonance fluorescence in ions [271]. The
development of spontaneous parametric down-conversion
(SPDC) made it possible to create a single-photon state
using only bulk optical elements [202]. However, all
these early non-Gaussian states were characterized through

counting statistics, which means that we generally classify
them as DV experiments.

It is perhaps intriguing to note that SPDC is also the
process that lies at the basis of the creation of squeezed
states of light [272], which are Gaussian. These states play
a key role in the generation of single-photon states, simply
because a weakly squeezed vacuum is mainly a super-
position of vacuum and a photon pair. By detecting one
photon of the pair, the presence of the second photon is
heralded. Hence, the approach of Ref. [202] is a basic
implementation of a conditional scheme for the generation
of non-Gaussian states as presented in Sec. B.

A genuine CV treatment of such non-Gaussian states
would only be achieved much later in a work that presents
the first tomographic reconstruction of a state with Wigner
negativity in optics [273]. Due to the developments of an
easily implementable maximum-likelihood algorithm for
state reconstruction, homodyne tomography became one
of the main tools to study non-Gaussian states in CV quan-
tum optics [274]. It did not take long before this also led to
the reconstruction of a displaced single-photon Fock state
[275] and a two-photon Fock state [276]. The combination
of increased squeezing with type-II SPDC and an array of
photon detectors to increase the number of heralded pho-
tons more recently made it possible to resolve the Wigner
function of a three-photon Fock state [277]. Similar ideas
of multiplexed photon detection have also been used to
generate superpositions of Fock states [278].

For non-Gaussian states beyond Fock states, photon
subtraction, as described in Sec. 2, is a common experi-
mental tool. Its first experimental implementation success-
fully showed the capability of generating non-Gaussian
statistics in the homodyne measurements, but it failed to
demonstrate Wigner negativity [148]. Later experiments
improved the quality of the generated states, demonstrat-
ing Wigner negativity and creating so-called “Schrödinger
kittens” [129,279,280]. The terminology is chosen because
these states resemble cat states proportional to |α⟩ − |−α⟩
for small values of the mean field α. Even though such
Schrödinger kittens are ultimately not very different from
squeezed single-photon states, the nomenclature makes
more sense in the context of experiments that “breed” cat
states [281]. Here, one mixes two Schrödinger kittens on
a beam splitter and performs homodyne detection on one
output port. By conditioning on instances where this homo-
dyne detector registers values close to zero, one effectively
heralds a larger cat state (the value of α has increased). A
variation of photon subtraction has also been used to create
a type of CV qubit [282].

As an alternative to photon subtraction, one can also
add a photon [159]. Even though this operation theoreti-
cally equates to applying a creation operator on the state,
it is experimentally much harder to implement than pho-
ton subtraction as it requires nonlinear optics. However,
photon subtraction can only produce Wigner negativity
when the initial state is squeezed. Photon addition, on the other hand, provides the advantage of always creating a Wigner negative state. A simple way to see this is by applying a creation operator to the state and evaluating the $Q$ function (70). When a photon is added to the mode $g$, the $Q$ function after photon addition has the property $Q^+(\vec{\alpha}) \sim (\vec{\alpha}^T \vec{g})^2 Q_G(\vec{\alpha})$, where $Q_G(\vec{\alpha})$ is the $Q$ function of the initial Gaussian state. This relation implies automatically that the $Q$ function will be exactly zero for $\vec{\alpha} = \vec{0}$, and a zero of the $Q$ function implies Wigner negativity. This means that one can apply photon addition to highly classical states, such as a coherent state or a thermal state, and still end up creating Wigner negativity. Such photon-added coherent states were also used to experimentally measure [283] non-Gaussianity $\delta(\hat{\rho})$ as defined in Eq. (110). Remarkably, combining photon addition and photon subtraction operations in both possible orders provides a way to experimentally verify the canonical commutation relations $[\hat{a}, \hat{a}^\dagger] = 1$, as was shown in Ref. [160].

The above methods are all based on Gaussian states as initial resources to generate non-Gaussian states. The non-Gaussian states that are created as such can in turn serve as useful resources to create more intricate non-Gaussian states. Fock states are a commonly used type of input state, for example, in the first demonstration of a large Schrödinger's cat state [284]. Intriguingly, by using non-Gaussian initial states, it suffices to use homodyne detection as the conditional operation. This setup can then be extended to a cat breeding scheme [285]. Another method to create large cat states in optics relies on making the light field interact with an atom [286]. The presence of entanglement between the “macroscopic” coherent state and the “microscopic” atomic degrees of freedom make for an experiment that resembles Schrödinger’s original though experiment [106]. Once the atom and the coherent light are entangled, a spin rotation of the atom is followed by a measurement to project the state of the light field in either an even or an odd cat state. This reflects the general idea that atoms still induce much larger nonlinearities than non-linear crystals. These nonlinearities are the direct source of non-Gaussian effect, but they are also much harder to control. At present, experiments that rely on such higher-order nonlinearities to create non-Gaussian states remain rare in the optical regime.

The above methods all focus on the creation of single-mode non-Gaussian states. For multimode systems, much of the experimental progress has concentrated on two-mode systems. As we extensively discussed throughout this tutorial, an important feature in such multimode systems are quantum correlations. Some of the first experimental demonstrations of non-Gaussian quantum correlations were based on the Bell inequality (292). Homodyne tomography and a single photon, delocalized over

two modes by a beam splitter, suffices to violate the inequality [287,288]. However, these works also teach us that extreme high purities are required to do so.

Motivated by photon-subtraction experiments and challenged by the no-go theorem of [225–227], entanglement distillation soon became a new focus for non-Gaussian quantum optics experiments. Some of these experiments have focused on adding some form of non-Gaussian noise on the initial state to circumvent the no-go theorem [228, 229]. Entanglement distillation through local photon subtraction from the entangled modes of a Gaussian input state would later be demonstrated in Ref. [289]. Earlier, it had already been shown that Gaussian entanglement can be increased by photon subtraction in a superposition of the entangled mode [290]. Interestingly, in the latter case, the photon is effectively subtracted in a nonentangled mode such that the setup is essentially equivalent to mixing a squeezed vacuum and a photon-subtracted squeezed vacuum on a beam splitter. A similar photon subtraction in a coherent superposition of modes was later carried out to entangle two Schrödinger kittens [291]. This can probably be seen as the first realizations of purely non-Gaussian entanglement in CV.

Photon addition has also been considered as a tool for creating entanglement between pairs of previously uncorrelated modes [292]. The resulting state can be seen as a hybrid entangled state proportional to $|0\rangle|\alpha\rangle + |1\rangle|-\alpha\rangle$, such states have also been produced using techniques similar to photon subtraction [293]. For two modes, photon addition can be implemented in a mode-selective way [294]. This setup is particularly useful to create entanglement between coherent states by adding a photon in a superposition of displaced modes.

Going beyond two modes has always remained a challenging task. For mode-selective photon subtraction from a multimode field, one must abandon the typical implementation based on a beam splitter. For two modes, such an alternative photon subtraction scheme was, for example, realized in the time-frequency domain, by subtracting a photon from a sideband [295]. Yet, going to a genuine multimode scenario required the design of a whole new photon subtractor based on sum-frequency generation [296,297]. This finally permitted the first demonstration of multimode non-Gaussian state in a CV setting, demonstrating non-Gaussian features in up to four entangled modes [298].

Such highly multimode states of more than two modes are confronted with a considerable problem: the exponential scaling of the required number of measurements for a full state tomography. This makes it highly challenging to demonstrate non-Gaussian features such as Wigner negativity in multimode non-Gaussian states. For single-photon-subtracted states, it has been pointed out that good analytical models can be used to train machine-learning
algorithms to recognize Wigner negativity based on single-mode measurements [299]. Furthermore, the techniques of Ref. [128] combined with Ref. [300] should also make it possible to use multiplexed double homodyne detection to witness Wigner negativity in certain classes of multimode states.

In multimode systems, we are confronted with the limitations of homodyne tomography. Recently, it has been shown that machine-learning techniques can be used to implement an improved form of CV tomography based on homodyne measurements [301]. Even though this setup is computationally heavy in single-mode setups, it uses a smaller set of states as a basis for state reconstruction, which might make multimode versions of the protocol more scalable. Alternatively, one can also bypass homodyne measurements all together. Photon-number-resolving detectors such as transition edge sensors [302] make it possible to use the identity, Eq. (78), to directly measure the Wigner function [76]. Intriguingly, this implies that a photon-number-resolving detector and a setup to generate displacements of the state in arbitrary modes makes it possible to directly measure the full multimode Wigner function. Nevertheless, such a multimode protocol has so far not been realized in any experiment.

## B. Other experimental setups

Given all the experimental work in CV quantum optics, it is perhaps surprising that the first experimental demonstrations of quantum states with Wigner negativity happened in different fields. The very first realization of such a state was achieved with trapped ions. Even though one often uses the atomic transitions in these systems to isolate qubits for potential quantum computers, trapped ions also have interesting motional degrees of freedom. By exploiting a Jaynes-Cummings type interaction between the atom and the trapping field, it is possible to use the ions' internal atomic degrees of freedom to create well-controlled non-Gaussian states such as a Fock state [78] and a Schrödinger's cat state [303] in the motional degrees of freedom.

Mathematically, this setup is equivalent to cavity QED, where it was shown that photons in a cavity can be manipulated through interactions with atoms [304] and the Rabi oscillations of the injected Rydberg atoms can in turn be used to probe the field within the cavity [305]. These methods would then be combined to experimentally generate a single-photon Fock state of the microwave field in a cavity [306], confirm its Wigner negativity [79], and probe its full Wigner function [80]. A few years later, similar techniques were used to finally generate Schrödinger cat states and higher-order Fock states [307].

A third setup with very similar physics is found in circuit QED. In this field, the macroscopic microwave cavities are replaced by superconducting circuits, and

nonlinearities are induced by Josephson junctions rather than atoms [140]. Even though these setups are often used in a DV approach, the microwave fields involved can equally be treated in a CV approach. The large nonlinearities rather naturally create non-Gaussian states, but getting a good sense of control over them can be challenging. Nevertheless, a wide range of non-Gaussian states such as Fock state [308] and large Schrödinger cat states [309] have been experimentally realized. The latter have furthermore been stabilized by engineering the decoherence processes in the system [310]. Very recently these systems have also been used to demonstrate the deterministic generation of photon triplets [142].

In recent years, both, trapped ions [311] and superconducting circuits [312] were used to achieve another important milestone in CV quantum computing: the experimental generation of a GKP state. These highly non-Gaussian states are useful for encoding a fault-tolerant qubit in a CV degree of freedom. By exploiting the redundancy that is offered by the infinite dimension Hilbert space of a CV system, one can create a qubit with a certain degree of robustness. This effectively makes it possible to implement error-correction routines, as shown in Ref. [312]. In other words, these systems have managed to generate CV states that are so non-Gaussian that they can be effectively used as fault-tolerant DV states.

A final field that has shown much potential over the last decades is cavity optomechanics. Here, an optical field is injected into a cavity with one moving mirror (more generally also other types of “dynamic cavities” can be used). The goal is to cool this mirror to its ground state to observe its quantum features. This way, one hopes to create nonclassical states of motion in reasonably large objects. A wide variety of such optomechanical devices exist [313]. Several theoretical schemes have been proposed to generate non-Gaussian states in such an optomechanical setup [314,315]. Even though quantum features such as photon-phonon entanglement have been demonstrated in such systems [316,317], it remains highly challenging to obtain good experimental control over the motional quantum state. Nevertheless, some CV non-Gaussian states can be realized in the form of superpositions between vacuum and a single-phonon Fock state have been experimentally realized [318].

A common problem in these setups is the creation of entanglement between the CV degrees of freedom in different modes. Some degree of such CV entanglement has been experimentally achieved in trapped ion [319] and circuit QED setups [320]. However, the number of entangled modes is much lower than what has been achieved in optics [18–21,220], where even non-Gaussian entangled states of more than two modes have been created [298]. This shows clearly how different experimental setups have different strengths and weaknesses. Optics comes with the advantage of spatial, temporal, and spectral mode
manipulations, which allows the creation of large entan-
gled states. However, the resilience of optical setups to
decoherence is due to limited interaction with the envi-
ronment. The latter implies that it is also difficult to
find controlled ways to make these systems strongly non-
Gaussian. On the other hand, the other setups, which we
discussed, require much more significant shielding from
environmental degrees of freedom. When this coupling to
other degrees of freedom can be controlled, it provides
the means to create non-Gaussian quantum states. In this
context, it is appealing to combine the advantages of differ-
ent regimes. Optomechanics offers a potential pathway to
achieve this by converting between microwave and optical
degrees of freedom [321,322].

As a last remark, it is interesting to mention that phase-
space descriptions and non-Gaussian states also appear
in atomic ensembles. This framework relies on the fact
that an ensemble of a large number of atoms can be
described by collective observables that behave very simi-
lar to bosonic systems. The associated phase space behaves
differently from the optical phase space, in the sense that it
is compact. More specifically, the phase space will cover
a sphere and the radius of this sphere will depend on
the number of atoms. Effectively, we would recover a
bosonic system in the limit of an infinite number of atoms.
However, the compactness of phase space for a finite
ensemble comes with interesting side effects: a sufficiently
high amount of spin squeezing can create non-Gaussian
states. We do not go into details for these systems, but it
should nevertheless be highlighted that non-Gaussian spin
states have received considerable attention in the literature
[323] and have been produced in a range of experiments
[324–327].

VIII. CONCLUSIONS AND OUTLOOK

In this Tutorial, we have presented a framework
based on phase-space representations to study continuous-
variable quantum systems. We then focused on the various
aspects of non-Gaussian states, where we first represented
different ways to structure the space of continuous-variable
states in a single mode in Fig. 1. Whenever possible, we
generalized results from the literature to a multimode set-
ting. However, for certain properties such as the stellar
rank, these generalizations become insufficient to classify
all possible quantum states.

We introduced two paradigms to create non-Gaussian
states, where one is a deterministic approach based on uni-
tary transformations, reminiscent of the circuit approach
for quantum-information processing. The second approach
is conditional, in the sense that it relies on condition-
ing one part of a state on measurement outcomes for
another part of a state, which is more narrowly related
to a measurement-based approach to quantum protocols.
Throughout the remainder of the Tutorial, we have largely

focused on conditional operations, since it is the most com-
monly used approach in experiments. It also provides a
natural avenue to start studying the relation between quan-
tum correlations and non-Gaussian features. We show how
the conditional approach requires certain correlations in
the initial Gaussian state to be able to induce certain type
of non-Gaussianity in the conditional state, as summarized
in Fig. 5.

On the other hand, non-Gaussian operations can also
create a type of non-Gaussian entanglement as introduced
in Sec. C. This kind of entanglement is particular as it
can not be identified with typical techniques that rely
on the state’s covariance matrix. Nevertheless, we use
Rényi-2 entanglement as a measure to illustrate the exis-
tence of such purely non-Gaussian quantum correlations
in photon-subtracted states. Even though its existence is
known from pure-state examples, it has only received lim-
ited attention in both theoretical and experimental work.
One possible reason is the difficulty of studying this type of
entanglement for mixed states, since convex roof construc-
tions tend to become highly intractable for non-Gaussian
states.

As a final theoretical aspect of the Tutorial, we high-
light the need of Wigner negativity to achieve some of the
most striking features in quantum technologies. On the one
hand, we show that Wigner negativity in either the state or
the measurement is necessary to violate a Bell inequality.
This observation can be understood in the broader con-
text of nonlocality and contextuality: Wigner negativity is
often seen as a manifestation of the contextual behavior of
quantum systems, and nonlocality can be understood as a
type of contextuality of measurements on different subsys-
tems. On the other hand, we also present results that show
how Wigner negativity is a requirement to achieve a quan-
tum computational advantage. Intuitively, it is perhaps not
surprising that states, operations, and measurements that
can all be described by probability distributions on phase
space can be efficiently simulated on a classical computer.
However, as we showed in Sec. VI, the actual simula-
tions protocol contains many subtle points. Here, too, we
conclude that there are still many open questions surround-
ing the physics of quantum computational advantages in
continuous-variable setups.

As a last step of this Tutorial, we provided an
overview of the experimental realizations of non-Gaussian
states with continuous variables. Quantum states of light
are indeed the usual suspects for continuous-variable
quantum-information processing, but it turns out to be
remarkably challenging to engineer highly non-Gaussian
states in such setups. We highlighted how trapped ions,
cavity QED, and circuit QED have proven to be bet-
ter equipped for this task, but in return they are con-
fronted with other problems. Optomechanics presents itself
as an ideal translator between these two regimes, which
may soon make it possible to combine the scalability of
optical setups with the high nonlinearities of the microwave domain.

In a more general sense, there are definitely many open question to be resolved in the domain of continuous-variable quantum physics. In this Tutorial, we have focused extensively on questions related to non-Gaussian features, notably in multimode systems. In the greater scheme of things, this is only one of the many challenges in the field. The recent demonstration of a quantum computational advantage with Gaussian Boson sampling has set an important milestone for continuous-variable quantum technologies [23], but we are still far away from useful computational protocols as set out in the roadmap of Ref. [24]. Even though the quest for a Gottesman-Kitaev-Preskill state [31] is one of the main experimental priorities, there are still many open challenges in designing the Gaussian operations that form the basis of such a setup [328,329].

Beyond universal fault-tolerant quantum computers, there are many other potential applications for continuous-variable systems. They are widely used in quantum communications for quantum key distribution [330] and secret sharing [331]. These protocols are largely based on Gaussian states and measurements, such that also the best possible attacks to these systems are Gaussian [332]. Nevertheless, non-Gaussian protocols for quantum key distribution, based on photon subtraction, have been proposed [333]. Such non-Gaussian quantum computation protocols and their security still involve many open questions.

Continuous-variable systems also provide a natural link to other bosonic systems, which is why they have been suggested as a platform to simulate molecular vibronic spectra [334]. The continuous-variable approach also plays an important role in quantum algorithms for other chemistry-related problems such as drug discovery through molecular docking [335] and the simulation of electron transport [336].

Furthermore, the continuous-variable setting is also suitable to implement certain elements for quantum machine learning such as quantum neural networks [337]. Even though this is a promising platform for tackling a wide range of problems, the proposal is highly ambitious on several points. In the context of this Tutorial, we emphasize the need of non-Gaussian unitary transformations. In principle, neural networks require linear couplings between different “neurons,” which each implement some form of nonlinear operation. Non-Gaussian operations play the role of this nonlinear element, making them a crucial step in the scheme. To implement such continuous-variable neural networks we thus require either new developments on the implementation of non-Gaussian operations, or theoretical modifications in the protocol to make it fit for implementable conditional non-Gaussian operations. It should be highlighted that other machine-learning approaches

exist, such as reservoir computing, which can be entirely based on Gaussian states [338].

A final quantum technology that may benefit from the use of non-Gaussian states is quantum metrology, as was recently demonstrated with motional Fock states of trapped ions [339]. Even though early work has shown that there is no clear benefit in using non-Gaussian operations such as photon subtraction for parameter estimation [154], there may still be other settings where such states are beneficial. Non-Gaussian entanglement could, for example, have a formal metrological advantage that is reflected in the quantum Fisher information [340]. On the other hand, ideas from quantum metrology also provide a possible approach for measuring non-Gaussian quantum steering [341]. The effects of non-Gaussian features on the sensitivity of the state can in principle be captured by higher moments of the quadrature operators [342]. It was recently shown that postselected measurements could, indeed, offer a quantum advantage for metrology [343]. This result is narrowly connected to the field of weak measurements and makes a connection to yet another phase-space representation: the Kirkwood-Dirac distribution [70,344]. Hence, we circle back to the fundamental physics of continuous-variable systems and conclude that there are still many connections to be made.

Beyond the technological applications that continuous-variable systems may have to offer, there is an important down-to-earth perspective that must be emphasized. With the improvement of detectors throughout the years, we have reached a point where theory and experiment can be considered mature to tackle single-mode problems. In multimode systems, the same cannot be said. With the exponential scaling of standard homodyne tomography, experimental tools for studying large multimode states beyond the Gaussian regime are limited. We may have to accept that the full quantum state is out of reach for experimental measurements. Even theoretically, highly multimode Wigner functions quickly become cumbersome to handle. Treating them with numerical integration techniques becomes a near-impossible task, once the number of modes is drastically increased. This makes even numerical simulations challenging. How then can we understand and even detect the non-Gaussian features of these systems?

One clear and important future research goal in this field is to provide an answer to this question. For quantum technologies, this may provide us with new ways to benchmark our systems, but more fundamentally it might teach us something new about the physics of these systems. One place where one might look for inspiration is the field of statistical mechanics, where statistical methods show that even highly complex systems can produce clear emergent signatures. We recently took a first step in exploring such ideas by looking at emergent network structures for continuous-variable non-Gaussian states [345]. The most exciting lesson from such preliminary work is that there
is still much to be learned about non-Gaussian quantum
states.

ACKNOWLEDGMENTS

First, I want to express sincere gratitude to the anony-
mous referees for providing very thorough reports and
many useful suggestions that significantly improved the
clarity and completeness of this Tutorial. I also thank
M. Genoni for several useful suggestions. More gener-
ally speaking, the content of this text was influenced by
stimulating discussions throughout the years with many
colleagues, notably F. Grosshans, R. Filip, Q. He, M. Gess-
ner, and G. Ferrini. Furthermore, I am very grateful to M.
Fannes for teaching me the mathematical foundations that
lie at the basis bosonic quantum systems. This point of
view was complemented by colleagues in the multimode
quantum optics group of the Laboratoire Kastler Brossel,
V. Parigi, N. Treps, and C. Fabre, who have introduced me
to the wonderful world of experimental quantum optics.
Still, the main source of inspiration for this Tutorial are the
many excellent students and postdocs that I have worked
with in the last few years. Their questions and struggles
have been essential to highlight the barriers that I try to
overcome in this Tutorial. Among these students and post-
docs, I want to express explicit gratitude to U. Chabaud, G.
Sorelli, and D. Barral for their careful and detailed read-
ing of the Tutorial and for their useful comments. I also
acknowledge the many useful discussions with K. Zhang,
who notably made me aware of the possibility of reducing
entanglement through photon subtraction (here shown in
Fig. 7). Last but not least, I want to express special thanks
to C. Lopetegui, who started reading this Tutorial as a new-
comeer to the field of CV quantum optics and thus had the
perfect point of view to help fine tune the content.

APPENDIX: MATHEMATICAL REMARKS

Here we present some important well-known mathemat-
ical concepts that are regularly used in the Tutorial to make
the text more self-contained. The comments and defini-
tions given here are not very rigorous and mainly aim at
giving the reader an intuitive understanding, for a more
formal introduction one should consult a standard textbook
[41,346,347].

1. Topological vector spaces

Throughout this Tutorial, we often deal implicitly with
topological vector spaces. Vector spaces are well known
from linear algebra and can be thought of as sets of mathe-
matical objects called vectors, which can be added together
in a commutative way and multiplied by scalars. When we
consider a vector space $\mathcal{V}$ on a field $\mathcal{F}$, this means that
for any $\vec{v}_1, \vec{v}_2 \in \mathcal{V}$ and any $\alpha_1, \alpha_2 \in \mathcal{F}$ the object $\alpha_1 \vec{v}_1 +
\alpha_2 \vec{v}_2 \in \mathcal{V}$. This means that the vector space is closed under

addition and scalar multiplication. In this Tutorial, the field
$\mathcal{F}$ is either identified as $\mathbb{R}$ (for phase space) or $\mathbb{C}$ (for
Hilbert spaces).

The spaces that are considered in the Tutorial have much
more structure than what is given by the vector space.
First of all, we generally deal with normed spaces, which
means that our vector spaces are topological vector spaces
in the sense that there is a notion of distance defined upon
them. Generally speaking, topological vector spaces can
be equipped with exotic topologies, but here we simply
deal with norms. On top, we again add an additional struc-
ture when we assume that these norms are generated by
inner products (depending on exact properties, these inner
products go by different names such as “positive-definite
sesquilinear form”, which is what we typically consider in
quantum mechanics).

As we deal with infinite-dimensional spaces to describe
bosonic quantum states and Fock space, it is important
to set some terminology straight. When we talk about a
Hilbert space, there is the assumption that the space is com-
plete. In an infinite-dimensional inner-product space, we
can define sequences of elements in $\mathcal{V}$. If we consider a
sequence $(\vec{v}_j)_{j \in \mathbb{N}}$ such that for any $\epsilon$ we can find a value
$N > 0$ such that $\|\vec{v}_j - \vec{v}_k\| < \epsilon$ for all $j, k > N$, we call
the sequence a Cauchy sequence. In other words, the dis-
tance between elements in the Cauchy sequence shrinks
as we proceed further into the sequence. The fact that a
Hilbert space is closed means that all Cauchy sequences
converge in the sense that $\lim_{j \to \infty} \vec{v}_j = \vec{v} \in \mathcal{V}$. Finite-
dimensional inner-product spaces automatically have this
property, but for infinite-dimensional spaces it must be
imposed explicitly.

Another class of structured vector space, that is often
encountered in the Tutorial, is a real symplectic space.
These spaces appear when we consider phase space, and
they are given by a real vector space with an additional
symplectic form $\sigma$ instead of the usual inner product. In
the mathematical literature, one often encounters the nota-
tion $(\mathcal{V}, \sigma)$ for a symplectic space, where the symplectic
form has the following properties: we consider $\vec{v}_1 \cdot \vec{v}_2 \in \mathcal{V}$
and find that $\sigma(\vec{v}_1, \vec{v}_2) \in \mathbb{R}$, $\sigma$ is bilinear, and $\sigma(\vec{v}_1, \vec{v}_2) =
-\sigma(\vec{v}_2, \vec{v}_1)$. In all cases in this Tutorial, we also consider
that $\sigma$ is nondegenerate, which means that $\sigma(\vec{v}_1, \vec{v}_2) = 0$
for all $\vec{v}_1 \in \mathcal{V}$ if and only if $\vec{v}_2 = \vec{0}$. When the symplectic
space is finite dimensional, it is often practical to represent
the symplectic form in terms of a matrix. In the Tutorial
this is done by associating $\sigma(\vec{v}_1, \vec{v}_2) = \vec{v}_1^T \Omega \vec{v}_2$.

In principle, a real symplectic space is all that is needed
to develop the mathematical framework of the CCR alge-
bra. However, it is often natural when dealing with bosonic
systems to include an additional structure in the form of an
inner product. In the Tutorial, this is done implicitly by also
using the standard inner product $\vec{v}_1^T \vec{v}_2$ on phase space. This
allows us to ultimately get the isomorphism (44). In the
quantum statistical mechanics literature, it is common to
see references to a “pre-Hilbert space”, rather than a phase space or a symplectic space. When we refer to a pre-Hilbert space, we consider an inner-product space, which is not necessarily complete and one must consider the closure to be guaranteed to obtain a full Hilbert space. The reason is that phase space, as a real vector space $V$ with an inner product, given by a bilinear form $s(.,.)$, and a symplectic form $\sigma(.,.)$ is equivalent to a complex pre-Hilbert space $\mathcal{H}$. For finite-dimensional spaces, the equivalence between the vector spaces is obtained via isomorphism (44):

$$ \vec{f} \in \mathcal{V} \mapsto |\psi_f\rangle = \sum_j (f_{2j-1} + if_{2j}) |\varphi_j\rangle \in \mathcal{H}, \quad (A1) $$

where $|\varphi_j\rangle$ for a basis of $\mathcal{H}$. As we are talking about an isomorphism between structured vectors spaces, we also need an identity between additional structures, which is given by

$$ \langle \psi_{f_1} | \psi_{f_2} \rangle = s(\vec{f}_1, \vec{f}_2) - i\sigma(\vec{f}_1, \vec{f}_2). \quad (A2) $$

This isomorphism holds very generally and can be extended to infinite-dimensional spaces. It provides a very formal connection between the single-particle Hilbert space for a many-boson system and its phase space associated with the modes of the bosonic field. Technically, we note that the phase space is equivalent to a pre-Hilbert space, and the closure of this space is the single-particle Hilbert space. Whenever the phase space (and thus the single-particle Hilbert space) is finite dimensional, the pre-Hilbert space is closed such that phase space and single-particle Hilbert space really are equivalent. For a very rigorous treatment on all these points, we refer to Ref. [42].

## 2. Span

Throughout the Tutorial, we often refer to the “span” of a certain set of vectors. These vectors can be members of a vector space, symplectic space, topological vectors space, pre-Hilbert space, or Hilbert space, the definition of the span is always the same. Let us here assume that $V$ denotes any type of vector space over a field $\mathcal{F}$ and consider a set $\vec{v}_1, \dots, \vec{v}_n \in V$. We can now define the span of this set of vectors as the set of all linear combinations that can be made with these vectors

$$ \mathrm{span}\{\vec{v}_1, \dots, \vec{v}_n\} := \{\alpha_1 \vec{v}_1 + \dots + \alpha_n \vec{v}_n \mid \alpha_1, \dots, \alpha_n \in \mathcal{F}\}. \quad (A3) $$

We emphasize that there is no need for the set $\vec{v}_1, \dots, \vec{v}_n$ to form a basis, nor for the vectors to be linearly independent, nor for the vectors to be normalized, nor for the vectors to be orthogonal to one another.

Throughout the Tutorial, the vector spaces we encounter are either real (in the case of phase space) such that $\mathcal{F} = \mathbb{R}$

or complex (in the case of Hilbert spaces for quantum systems) such that $\mathcal{F} = \mathbb{C}$. In the case where the vector spaces have some topological structure (which we can colloquially understand as a mathematical sense of distance that allows limits to be defined), it can make sense to consider the closure of a span, denoted by

$$ \overline{\mathrm{span}\{\vec{v}_1, \dots, \vec{v}_n\}}, \quad (A4) $$

such that any convergent sequence built with elements of the span has its limit also included in the closure.

[1] E. Schrödinger, Der stetige übergang von der mikro- zur makromechanik, Naturwissenschaften 14, 664 (1926).

[2] H. Araki and E. J. Woods, Representations of the canonical commutation relations describing a nonrelativistic infinite free bose gas, J. Math. Phys. 4, 637 (1963).

[3] D. W. Robinson, The ground state of the Bose gas, Commun. Math. Phys. 1, 159 (1965).

[4] A. Verbeure, Many-Body Boson Systems: Half a Century Later, Theoretical and mathematical physics (Springer, London; New York, 2011).

[5] R. J. Glauber, Coherent and incoherent states of the radiation field, Phys. Rev. 131, 2766 (1963).

[6] E. C. G. Sudarshan, Equivalence of Semiclassical and Quantum Mechanical Descriptions of Statistical Light Beams, Phys. Rev. Lett. 10, 277 (1963).

[7] S. L. Braunstein and P. van Loock, Quantum information with continuous variables, Rev. Mod. Phys. 77, 513 (2005).

[8] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Rev. Mod. Phys. 84, 621 (2012).

[9] G. Adesso, S. Ragy, and A. R. Lee, Continuous variable quantum information: Gaussian states and beyond, Open Syst. Inf. Dynamics 21, 1440001 (2014).

[10] J. Williamson, On the algebraic problem concerning the normal forms of linear dynamical systems, American J. Math. 58, 141 (1936).

[11] S. L. Braunstein, Squeezing as an irreducible resource, Phys. Rev. A 71, 055801 (2005).

[12] M. M. Wolf, G. Giedke, and J. I. Cirac, Extremality of Gaussian Quantum States, Phys. Rev. Lett. 96, 080502 (2006).

[13] A. S. Holevo, M. Sohma, and O. Hirota, Capacity of quantum gaussian channels, Phys. Rev. A 59, 1820 (1999).

[14] M. Yukawa, R. Ukai, P. van Loock, and A. Furusawa, Experimental generation of four-mode continuous-variable cluster states, Phys. Rev. A 78, 012301 (2008).

[15] X. Su, Y. Zhao, S. Hao, X. Jia, C. Xie, and K. Peng, Experimental preparation of eight-partite cluster state for photonic qumodes, Opt. Lett. 37, 5178 (2012).

[16] D. Barral, M. Walschaers, K. Bencheikh, V. Parigi, J. A. Levenson, N. Treps, and N. Belabas, Versatile Photonic Entanglement Synthesizer in the Spatial Domain, Phys. Rev. Appl. 14, 044025 (2020).
[17] J. Roslund, R. M. de Araújo, S. Jiang, C. Fabre, and N. Treps, Wavelength-multiplexed quantum networks with ultrafast frequency combs, Nat Photon **8**, 109 (2014).

[18] M. Chen, N. C. Menicucci, and O. Pfister, Experimental Realization of Multipartite Entanglement of 60 Modes of a Quantum Optical Frequency Comb, Phys. Rev. Lett. **112**, 120505 (2014).

[19] Y. Cai, J. Roslund, G. Ferrini, F. Arzani, X. Xu, C. Fabre, and N. Treps, Multimode entanglement in recon-figurable graph states using optical frequency combs, Nat. Commun. **8**, 15645 (2017).

[20] W. Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoen-sombutamon, H. Emura, R. N. Alexander, S. Takeda, J.-i. Yoshikawa, N. C. Menicucci, H. Yonezawa, and A. Furusawa, Generation of time-domain-multiplexed two-dimensional cluster state, Science **366**, 373 (2019).

[21] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Deterministic generation of a two-dimensional cluster state, Science **366**, 369 (2019).

[22] Z. Yang, M. Jahanbozorgi, D. Jeong, S. Sun, O. Pfister, H. Lee, and X. Yi, A squeezed quantum microcomb on a chip, arXiv:2103.03380 [physics.optics] (2021).

[23] H.-S. Zhong et al., Quantum computational advantage using photons, Science 370, 1460 (2020), https://science.sciencemag.org/content/370/6523/1460.full.pdf.

[24] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, K. K. Sabapathy, N. C. Menicucci, and I. Dhand, Blueprint for a scalable photonic fault-tolerant quantum computer, Quantum **5**, 392 (2021).

[25] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Efficient Classical Simulation of Continuous Variable Quantum Information Processes, Phys. Rev. Lett. **88**, 097904 (2002).

[26] S. Lloyd and S. L. Braunstein, Quantum Computation Over Continuous Variables, Phys. Rev. Lett. **82**, 1784 (1999).

[27] N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, Universal Quantum Compu-tation with Continuous-Variable Cluster States, Phys. Rev. Lett. **97**, 110501 (2006).

[28] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and P. van Loock, Quantum computing with continuous-variable clusters, Phys. Rev. A **79**, 062318 (2009).

[29] N. C. Menicucci, Fault-Tolerant Measurement-Based Quantum Computing with Continuous-Variable Cluster States, Phys. Rev. Lett. **112**, 120504 (2014).

[30] F. Arzani, N. Treps, and G. Ferrini, Polynomial approxi-mation of non-gaussian unitaries by counting one photon at a time, Phys. Rev. A **95**, 052352 (2017).

[31] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Phys. Rev. A **64**, 012310 (2001).

[32] B. Q. Baragiola, G. Pantaleoni, R. N. Alexander, A. Karanjai, and N. C. Menicucci, All-Gaussian Universality and Fault Tolerance with the Gottesman-Kitaev-Preskill Code, Phys. Rev. Lett. **123**, 200502 (2019).

[33] A. Mari and J. Eisert, Positive Wigner Functions Render Classical Simulation of Quantum Computation Efficient, Phys. Rev. Lett. **109**, 230503 (2012).

[34] R. Hudson, When is the Wigner quasi-probability density non-negative?, Rep. Math. Phys. **6**, 249 (1974).

[35] F. Soto and P. Claverie, When is the Wigner function of multidimensional systems nonnegative?, J. Math. Phys. **24**, 97 (1983).

[36] A. Mandilara, E. Karpov, and N. J. Cerf, Extending Hudson's theorem to mixed quantum states, Phys. Rev. A **79**, 062302 (2009).

[37] L. c. v. Lachman, I. Straka, J. Hloušek, M. Ježek, and R. Filip, Faithful Hierarchy of Genuine *n*-Photon Quantum Non-Gaussian Light, Phys. Rev. Lett. **123**, 043601 (2019).

[38] U. Chabaud, D. Markham, and F. Grosshans, Stellar Representation of Non-Gaussian Quantum States, Phys. Rev. Lett. **124**, 063605 (2020).

[39] M. Walschaers, V. Parigi, and N. Treps, Practical framework for conditional non-gaussian quantum state preparation, PRX Quantum **1**, 020305 (2020).

[40] D. Petz, *An Invitation to the Algebra of Canonical Comutation Relations*, Leuven notes in mathematical and theoretical physics Series A No. 2 (Leuven Univ. Press, Leuven, 1990).

[41] O. Bratteli and D. W. Robinson, *Operator Algebras and Quantum Statistical Mechanics 1* (Springer Berlin Heidelberg, Berlin, Heidelberg, 1987).

[42] O. Bratteli and D. W. Robinson, *Operator Algebras and Quantum Statistical Mechanics Equilibrium States. Models in Quantum Statistical Mechanics* (Springer, Berlin, 1997).

[43] M. H. Stone, Linear transformations in Hilbert space: III. Operational methods and group theory, PNAS **16**, 172 (1930).

[44] J. von Neumann, Die eindeutigkeit der Schrödingerschen operatoren, Math. Ann. **104**, 570 (1931).

[45] M. H. Stone, On one-parameter unitary groups in Hilbert space, Ann. Math. **33**, 643 (1932).

[46] J. von Neumann, Über einen satz von herrn M. H. Stone, Ann. Math. **33**, 567 (1932).

[47] M. C. Tichy, F. Mintert, and A. Buchleitner, Essential entanglement for atomic and molecular physics, J. Phys. B: At. Mol. Opt. Phys. **44**, 192001 (2011).

[48] F. Benatti, R. Floreanini, F. Franchini, and U. Marzolino, Entanglement in indistinguishable particle systems, Phys. Rep. **878**, 1 (2020), entanglement in indistinguishable particle systems.

[49] N. Killoran, M. Cramer, and M. B. Plenio, Extracting Entanglement from Identical Particles, Phys. Rev. Lett. **112**, 150501 (2014).

[50] N. Killoran, F. E. S. Steinhoff, and M. B. Plenio, Converting Nonclassicality Into Entanglement, Phys. Rev. Lett. **116**, 080402 (2016).

[51] R. Lo Franco and G. Compagno, Indistinguishability of Elementary Systems as a Resource for Quantum Information Processing, Phys. Rev. Lett. **120**, 240403 (2018).

[52] B. Morris, B. Yadin, M. Fadel, T. Zibold, P. Treutlein, and G. Adesso, Entanglement between Identical Particles is a Useful and Consistent Resource, Phys. Rev. X **10**, 041012 (2020).

[53] M. C. Tichy, Interference of identical particles from entanglement to boson-sampling, J. Phys. B: At. Mol. Opt. Phys. **47**, 103001 (2014).
[54] Note that the notation is not a coincidence, we can indeed construct a formal linear map from the single-particle Hilbert space into the operator algebra of linear operators on the Fock space. This mapping takes single-particle state vectors $\psi \in \mathcal{H}$ and maps them to a creation operators. The mapping is linear in the sense that $\hat{\alpha}^{\dagger}(x\psi + y\varphi) = x\hat{\alpha}^{\dagger}(\psi) + y\hat{\alpha}^{\dagger}(\varphi)$, for all $x, y \in \mathbb{C}$.

[55] L. Mandel and E. Wolf, *Optical Coherence and Quantum Optics* (Cambridge University Press, Cambridge, 1995).

[56] M. Walschaers, Signatures of many-particle interference, *J. Phys. B: At. Mol. Opt. Phys.* **53**, 043001 (2020).

[57] W. P. Schleich, *Quantum Optics in Phase Space* (Wiley VCH, Berlin, 2001).

[58] C. Fabre and N. Treps, Modes and states in quantum optics, *Rev. Mod. Phys.* **92**, 035005 (2020).

[59] A. M. O. de Almeida, The weyl representation in classical and quantum mechanics, *Phys. Rep.* **295**, 265 (1998).

[60] In the literature, one will encounter various different choices for the symplectic structure, which correspond to different forms of ordering the amplitude and phase quadratures. In general, any matrix $J$ that satisfies the conditions $J^T = -J$ and $J^2 = -1$ defines a symplectic structure. A popular alternative to the ordering $(E_1^{(x)}, E_1^{(p)}, \dots, E_m^{(x)}, E_m^{(p)})$ that we follow in this Tutorial is the ordering $(E_1^{(x)}, \dots, E_m^{(x)}, E_1^{(p)}, \dots, E_m^{(p)})$. The choice in this Tutorial is motivated by the study of entanglement, where it is more convenient to group the quadratures that correspond to the same modes together.

[61] E. Wigner, On the quantum correction for thermodynamic equilibrium, *Phys. Rev.* **40**, 749 (1932).

[62] M. Hillery, R. O'Connell, M. Scully, and E. Wigner, Distribution functions in physics: Fundamentals, *Phys. Rep.* **106**, 121 (1984).

[63] W. P. Schleich, in *Quantum Optics in Phase Space* (John Wiley and Sons, Ltd, 2001), Chap. 3, p. 67.

[64] R. W. Spekkens, Negativity and Contextuality are Equivalent Notions of Nonclassicality, *Phys. Rev. Lett.* **101**, 020401 (2008).

[65] J. von Neumann, *Mathematical Foundations of Quantum Mechanics* (Princeton University Press, Berlin, Heidelberg, 1955).

[66] H. Weyl, Quantenmechanik und gruppentheorie, *Zeitschrift für Physik* **46**, 1 (1927).

[67] J. B. Conway, in *A Course in Functional Analysis* (Springer New York, New York, NY, 1985), p. 310.

[68] M. Plancherel, Contribution à l'étude de la représentation d'une fonction arbitraire par des intégrales définies, *Rendiconti del Circolo Matematico di Palermo* (1884-1940) **30**, 289 (1910).

[69] K. E. Cahill and R. J. Glauber, Density operators and quasiprobability distributions, *Phys. Rev.* **177**, 1882 (1969).

[70] P. A. M. Dirac, On the analogy between classical and quantum mechanics, *Rev. Mod. Phys.* **17**, 195 (1945).

[71] J. Sperling and W. Vogel, Quasiprobability distributions for quantum-optical coherence and beyond, *Physica Scripta* **95**, 034007 (2020).

[72] A. Royer, Wigner function as the expectation value of a parity operator, *Phys. Rev. A* **15**, 449 (1977).

[73] K. Banaszek and K. Wódkiewicz, Direct Probing of Quantum Phase Space by Photon Counting, *Phys. Rev. Lett.* **76**, 4344 (1996).

[74] K. Banaszek, C. Radzewicz, K. Wódkiewicz, and J. S. Krasiński, Direct measurement of the wigner function by photon counting, *Phys. Rev. A* **60**, 674 (1999).

[75] K. Laiho, K. N. Cassemiro, D. Gross, and C. Silberhorn, Probing the Negative Wigner Function of a Pulsed Single Photon Point by Point, *Phys. Rev. Lett.* **105**, 253603 (2010).

[76] R. Nehra, A. Win, M. Eaton, R. Shahrokhshahi, N. Sridhar, T. Gerrits, A. Lita, S. W. Nam, and O. Pfister, State-independent quantum state tomography by photon-number-resolving measurements, *Optica* **6**, 1356 (2019).

[77] L. G. Lutterbach and L. Davidovich, Method for Direct Measurement of the Wigner Function in Cavity qed and ion Traps, *Phys. Rev. Lett.* **78**, 2547 (1997).

[78] D. Leibfried, D. M. Meekhof, B. E. King, C. Monroe, W. M. Itano, and D. J. Wineland, Experimental Determination of the Motional Quantum State of a Trapped Atom, *Phys. Rev. Lett.* **77**, 4281 (1996).

[79] G. Nogues, A. Rauschenbeutel, S. Osnaghi, P. Bertet, M. Brune, J. M. Raimond, S. Haroche, L. G. Lutterbach, and L. Davidovich, Measurement of a negative value for the wigner function of radiation, *Phys. Rev. A* **62**, 054101 (2000).

[80] P. Bertet, A. Auffeves, P. Maioli, S. Osnaghi, T. Meunier, M. Brune, J. M. Raimond, and S. Haroche, Direct Measurement of the Wigner Function of a One-Photon Fock State in a Cavity, *Phys. Rev. Lett.* **89**, 200402 (2002).

[81] A. I. Lvovsky, P. Grangier, A. Ourjoumtsev, V. Parigi, M. Sasaki, and R. Tualle-Brouri, Production and applications of non-gaussian quantum states of light, *arXiv:2006.16985 [quant-ph]* (2020).

[82] F. Flamini, N. Spagnolo, and F. Sciarrino, Photonic quantum information processing: A review, *Rep. Prog. Phys.* **82**, 016001 (2018).

[83] R. Simon, N. Mukunda, and B. Dutta, Quantum-noise matrix for multimode systems: U(n) invariance, squeezing, and normal forms, *Phys. Rev. A* **49**, 1567 (1994).

[84] H. P. Robertson, The uncertainty principle, *Phys. Rev.* **34**, 163 (1929).

[85] Arvind, B. Dutta, N. Mukunda, and R. Simon, The real symplectic groups in quantum mechanics and optics, *Pramana* **45**, 471 (1995).

[86] J. Eisert and M. M. Wolf, in *Quantum Information with Continuous Variables of Atoms and Light* (World Scientific, 2007), p. 23.

[87] B. Demoen, P. Vanheuverzwijn, and A. Verbeure, Completely positive maps on the CCR-algebra, *Lett Math Phys* **2**, 161 (1977).

[88] B. Demoen, P. Vanheuverzwijn, and A. Verbeure, Completely positive quasi-free maps of the CCR-algebra, *Rep. Math. Phys.* **15**, 27 (1979).

[89] R. Filip and L. Mišta, Detecting Quantum States with a Positive Wigner Function beyond Mixtures of Gaussian States, *Phys. Rev. Lett.* **106**, 200401 (2011).

[90] A. Kenfack and K. Życzkowski, Negativity of the wigner function as an indicator of non-classicality, *J. Opt. B: Quantum Semiclassical Opt.* **6**, 396 (2004).
[91] J. L. W. V. Jensen, Sur les fonctions convexes et les inégalités entre les valeurs moyennes, Acta Math. **30**, 175 (1906).

[92] G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A **65**, 032314 (2002).

[93] M. B. Hastings, Superadditivity of communication capacity using entangled inputs, Nat. Phys. **5**, 255 (2009).

[94] C. D. Cushen and R. L. Hudson, A quantum-mechanical central limit theorem, J. Appl. Probab. **8**, 454 (1971).

[95] J. Quaegebeur, A noncommutative central limit theorem for ccr-algebras, J. Funct. Anal. **57**, 1 (1984).

[96] D. Goderis, A. Verbeure, and P. Vets, Non-commutative central limits, Probab. Theory Relat. Fields **82**, 527 (1989).

[97] M. G. Genoni, M. G. A. Paris, and K. Banaszek, Measure of the non-gaussian character of a quantum state, Phys. Rev. A **76**, 042327 (2007).

[98] M. G. Genoni, M. G. A. Paris, and K. Banaszek, Quantifying the non-gaussian character of a quantum state by quantum relative entropy, Phys. Rev. A **78**, 060303 (2008).

[99] M. G. Genoni and M. G. A. Paris, Quantifying non-gaussianity for quantum information, Phys. Rev. A **82**, 052341 (2010).

[100] F. Albarelli, M. G. Genoni, M. G. A. Paris, and A. Ferraro, Resource theory of quantum non-gaussianity and wigner negativity, Phys. Rev. A **98**, 052350 (2018).

[101] P. Marian and T. A. Marian, Relative entropy is an exact measure of non-gaussianity, Phys. Rev. A **88**, 012322 (2013).

[102] J. S. Ivan, M. S. Kumar, and R. Simon, A measure of non-gaussianity for quantum states, Quantum Inf. Process. **11**, 853 (2012).

[103] N. Lütkenhaus and S. M. Barnett, Nonclassical effects in phase space, Phys. Rev. A **51**, 3340 (1995).

[104] A function $f$ is said to have an order of growth $r$ when there are constant $a, b \in \mathbb{C}$ such that $|f(z)| \le a e^{\exp(b|z|^r)}$ for all $z \in \mathbb{C}$.

[105] B. Yurke and D. Stoler, Generating Quantum Mechanical Superpositions of Macroscopically Distinguishable States via Amplitude Dispersion, Phys. Rev. Lett. **57**, 13 (1986).

[106] E. Schrödinger, Die gegenwärtige situation in der quan-tenmechanik, Naturwissenschaften **23**, 807 (1935).

[107] P. T. Cochrane, G. J. Milburn, and W. J. Munro, Macroscopically distinct quantum-superposition states as a bosonic code for amplitude damping, Phys. Rev. A **59**, 2631 (1999).

[108] Z. Leghtas, G. Kirchmair, B. Vlastakis, R. J. Schoelkopf, M. H. Devoret, and M. Mirrahimi, Hardware-Efficient Autonomous Quantum Memory Protection, Phys. Rev. Lett. **111**, 120501 (2013).

[109] B. C. Sanders, Quantum dynamics of the nonlinear rotator and the effects of continual spin measurement, Phys. Rev. A **40**, 2417 (1989).

[110] A. N. Boto, P. Kok, D. S. Abrams, S. L. Braunstein, C. P. Williams, and J. P. Dowling, Quantum Interferometric Optical Lithography: Exploiting Entanglement to Beat the Diffraction Limit, Phys. Rev. Lett. **85**, 2733 (2000).

[111] H. Lee, P. Kok, and J. P. Dowling, A quantum rosetta stone for interferometry, J. Mod. Opt. **49**, 2325 (2002).

[112] B. Yurke, Input States for Enhancement of Fermion Interferometer Sensitivity, Phys. Rev. Lett. **56**, 1515 (1986).

[113] Y. Israel, I. Afek, S. Rosen, O. Ambar, and Y. Silberberg, Experimental tomography of noon states with large photon numbers, Phys. Rev. A **85**, 022115 (2012).

[114] M. G. Genoni, M. L. Palma, T. Tufarelli, S. Olivares, M. S. Kim, and M. G. A. Paris, Detecting quantum non-Gaussianity via the Wigner function, Phys. Rev. A **87**, 062104 (2013).

[115] C. Hughes, M. G. Genoni, T. Tufarelli, M. G. A. Paris, and M. S. Kim, Quantum non-Gaussianity witnesses in phase space, Phys. Rev. A **90**, 013810 (2014).

[116] J. Park, J. Zhang, J. Lee, S.-W. Ji, M. Um, D. Lv, K. Kim, and H. Nha, Testing Nonclassicality and Non-Gaussianity in Phase Space, Phys. Rev. Lett. **114**, 190402 (2015).

[117] L. Happ, M. A. Efremov, H. Nha, and W. P. Schleich, Sufficient condition for a quantum state to be genuinely quantum non-Gaussian, New J. Phys. **20**, 023046 (2018).

[118] B. Kühn and W. Vogel, Quantum non-Gaussianity and quantification of nonclassicality, Phys. Rev. A **97**, 053823 (2018).

[119] M. Bohmann, E. Agudelo, and J. Sperling, Probing non-classicality with matrices of phase-space distributions, Quantum **4**, 343 (2020).

[120] I. Straka, A. Predojević, T. Huber, L. c. v. Lachman, L. Butschek, M. Miková, M. Mičuda, G. S. Solomon, G. Weihs, M. Ježek, and R. Filip, Quantum Non-Gaussian Depth of Single-Photon States, Phys. Rev. Lett. **113**, 223603 (2014).

[121] I. Straka, L. Lachman, J. Hloušek, M. Miková, M. Mičuda, M. Ježek, and R. Filip, Quantum non-Gaussian multiphoton light, Npj Quantum Inf. **4**, 4 (2018).

[122] C. N. Gagatsos and S. Guha, Efficient representation of Gaussian states for multimode non-Gaussian quantum state engineering via subtraction of arbitrary number of photons, Phys. Rev. A **99**, 053816 (2019).

[123] U. Chabaud, G. Roeland, M. Walschaers, F. Grosshans, V. Parigi, D. Markham, and N. Treps, Certification of non-Gaussian states with operational measurements, PRX Quantum **2**, 020333 (2021).

[124] P. Lelong and L. Gruman, *Entire Functions of Several Complex Variables* (Springer-Verlag Berlin Heidelberg, 1986).

[125] U. Chabaud, G. Ferrini, F. Grosshans, and D. Markham, Classical simulation of Gaussian quantum circuits with non-Gaussian input states, Phys. Rev. Res. **3**, 033018 (2021).

[126] R. Takagi and Q. Zhuang, Convex resource theory of non-Gaussianity, Phys. Rev. A **97**, 062337 (2018).

[127] V. Veitch, S. A. H. Mousavian, D. Gottesman, and J. Emerson, The resource theory of stabilizer quantum computation, New J. Phys. **16**, 013009 (2014).

[128] U. Chabaud, P.-E. Emeriau, and F. Grosshans, Witnessing Wigner negativity, Quantum **5**, 471 (2021).

[129] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier, Generating optical Schrödinger kittens for quantum information processing, Science **312**, 83 (2006).
[130] A. Mari, K. Kieling, B. M. Nielsen, E. S. Polzik, and J. Eisert, Directly Estimating Nonclassicality, Phys. Rev. Lett. 106, 010403 (2011).

[131] R. Loudon and P. Knight, Squeezed light, J. Mod. Opt. 34, 709 (1987).

[132] L.-A. Wu, H. J. Kimble, J. L. Hall, and H. Wu, Generation of Squeezed States by Parametric Down Conversion, Phys. Rev. Lett. 57, 2520 (1986).

[133] D. P. DiVincenzo, Quantum computation, Science 270, 255 (1995).

[134] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, Experimental Realization of any Discrete Unitary Operator, Phys. Rev. Lett. 73, 58 (1994).

[135] P. Marek, R. Filip, and A. Furusawa, Deterministic implementation of weak quantum cubic nonlinearity, Phys. Rev. A 84, 053802 (2011).

[136] K. Marshall, R. Pooser, G. Siopsis, and C. Weedbrook, Repeat-until-success cubic phase gate for universal continuous-variable quantum computation, Phys. Rev. A 91, 032321 (2015).

[137] K. Miyata, H. Ogawa, P. Marek, R. Filip, H. Yonezawa, J.-i. Yoshikawa, and A. Furusawa, Implementation of a quantum cubic gate by an adaptive non-Gaussian measurement, Phys. Rev. A 93, 022301 (2016).

[138] R. Yanagimoto, T. Onodera, E. Ng, L. G. Wright, P. L. McMahon, and H. Mabuchi, Engineering a Kerr-Based Deterministic Cubic Phase Gate via Gaussian Operations, Phys. Rev. Lett. 124, 240503 (2020).

[139] Y. Zheng, O. Hahn, P. Stadler, P. Holmvall, F. Quijandría, A. Ferraro, and G. Ferrini, Gaussian conversion protocols for cubic phase state generation, PRX Quantum 2, 010327 (2021).

[140] D. Vion, A. Aassime, A. Cottet, P. Joyez, H. Pothier, C. Urbina, D. Esteve, and M. H. Devoret, Manipulating the quantum state of an electrical circuit, Science 296, 886 (2002).

[141] M. H. Devoret and R. J. Schoelkopf, Superconducting circuits for quantum information: An outlook, Science 339, 1169 (2013).

[142] C. W. S. Chang, C. Sabín, P. Forn-Díaz, F. Quijandría, A. M. Vadiraj, I. Nsanzineza, G. Johansson, and C. M. Wilson, Observation of Three-Photon Spontaneous Parametric Down-Conversion in a Superconducting Parametric Cavity, Phys. Rev. X 10, 011011 (2020).

[143] K. Kraus, General state changes in quantum theory, Ann. Phys. (N. Y) 64, 311 (1971).

[144] Statistics of quantum measurements, in Statistical Structure of Quantum Theory (Springer Berlin Heidelberg, Berlin, Heidelberg, 2001) p. 39.

[145] R. A. Horn and C. R. Johnson, Matrix Analysis (Cambridge University Press, New York, NY, 2017), 2nd ed. corrected reprint ed.

[146] L. Lami, A. Serafini, and G. Adesso, Gaussian entanglement revisited, New J. Phys. 20, 023030 (2018).

[147] D. E. Browne, J. Eisert, S. Scheel, and M. B. Plenio, Driving non-Gaussian to Gaussian states with linear optics, Phys. Rev. A 67, 062320 (2003).

[148] J. Wenger, R. Tualle-Brouri, and P. Grangier, Non-Gaussian Statistics from Individual Pulses of Squeezed Light, Phys. Rev. Lett. 92, 153601 (2004).

[149] A. Zavatta, V. Parigi, M. S. Kim, and M. Bellini, Subtracting photons from arbitrary light fields: Experimental test of coherent state invariance by single-photon annihilation, New J. Phys. 10, 123006 (2008).

[150] V. Averchenko, C. Jacquard, V. Thiel, C. Fabre, and N. Treps, Multimode theory of single-photon subtraction, New J. Phys. 18, 083042 (2016).

[151] H. M. Wiseman and G. J. Milburn, *Quantum Measurement and Control* (Cambridge University Press, Cambridge, 2009).

[152] M. Walschaers, C. Fabre, V. Parigi, and N. Treps, Entanglement and Wigner Function Negativity of Multimode Non-Gaussian States, Phys. Rev. Lett. 119, 183601 (2017).

[153] M. Walschaers, C. Fabre, V. Parigi, and N. Treps, Statistical signatures of multimode single-photon-added and -subtracted states of light, Phys. Rev. A 96, 053835 (2017).

[154] D. Braun, P. Jian, O. Pinel, and N. Treps, Precision measurements with photon-subtracted or photon-added gaussian states, Phys. Rev. A 90, 013821 (2014).

[155] K. Mølmer, Non-Gaussian states from continuous-wave Gaussian light sources, Phys. Rev. A 73, 063804 (2006).

[156] M. Dakna, T. Anhut, T. Opatrný, L. Knöll, and D.-G. Welsch, Generating Schrödinger-cat-like states by means of conditional measurements on a beam splitter, Phys. Rev. A 55, 3184 (1997).

[157] K. Takase, J.-i. Yoshikawa, W. Asavanant, M. Endo, and A. Furusawa, Generation of optical Schrödinger cat states by generalized photon subtraction, Phys. Rev. A 103, 013710 (2021).

[158] A. O. C. Davis, M. Walschaers, V. Parigi, and N. Treps, Conditional preparation of non-Gaussian quantum optical states by mesoscopic measurement, New J. Phys. 23, 063039 (2021).

[159] A. Zavatta, S. Viciani, and M. Bellini, Quantum-to-classical transition with single-photon-added coherent states of light, Science 306, 660 (2004).

[160] V. Parigi, A. Zavatta, M. Kim, and M. Bellini, Probing quantum commutation rules by addition and subtraction of single photons to/from a light field, Science 317, 1890 (2007).

[161] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).

[162] R. Uola, A. C. S. Costa, H. C. Nguyen, and O. Gühne, Quantum steering, Rev. Mod. Phys. 92, 015001 (2020).

[163] N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Bell nonlocality, Rev. Mod. Phys. 86, 419 (2014).

[164] One could introduce more general and rigorous notation based on measure theory, but here we restrict to a simpler, though less general, formulation for pedagogical purposes.

[165] R. Jozsa, Quantum effects in algorithms, Chaos, Solitons Fractals 10, 1657 (1999).

[166] This is true for general “states” in classical probability theory, but in the CV context of this Tutorial it is appealing to
phrase the argument in terms of probability distributions of phase space.

[167] R. F. Werner, Quantum states with einstein-podolsky-rosen correlations admitting a hidden-variable model, Phys. Rev. A **40**, 4277 (1989).

[168] J. B. Conway, in *A Course in Functional Analysis* (Springer New York, New York, NY, 2007), p. 63.

[169] Sohail and U. Sen, Witnessing nonseparability of bipartite quantum operations, Phys. Lett. A **404**, 127411 (2021).

[170] H. F. Hofmann and S. Takeuchi, Violation of local uncertainty relations as a signature of entanglement, Phys. Rev. A **68**, 032103 (2003).

[171] O. Gittsovich, O. Gühne, P. Hyllus, and J. Eisert, Unifying several separability conditions using the covariance matrix criterion, Phys. Rev. A **78**, 052319 (2008).

[172] A. Einstein, B. Podolsky, and N. Rosen, Can quantum-mechanical description of physical reality be considered complete?, Phys. Rev. **47**, 777 (1935).

[173] E. Schrödinger, Discussion of probability relations between separated systems, Math. Proc. Cambridge Philosophical Soc. **31**, 555 (1935).

[174] E. Schrödinger, Probability relations between separated systems, Math. Proc. Cambridge Philosophical Soc. **32**, 446 (1936).

[175] M. D. Reid, Demonstration of the einstein-podolsky-rosen paradox using nondegenerate parametric amplification, Phys. Rev. A **40**, 913 (1989).

[176] H. M. Wiseman, S. J. Jones, and A. C. Doherty, Steering, Entanglement, Nonlocality, and the Einstein-Podolsky-Rosen Paradox, Phys. Rev. Lett. **98**, 140402 (2007).

[177] D. Cavalcanti and P. Skrzypczyk, Quantum steering: A review with focus on semidefinite programming, Rep. Prog. Phys. **80**, 024001 (2017).

[178] E. G. Cavalcanti, S. J. Jones, H. M. Wiseman, and M. D. Reid, Experimental criteria for steering and the einstein-podolsky-rosen paradox, Phys. Rev. A **80**, 032112 (2009).

[179] I. Kogias, A. R. Lee, S. Ragy, and G. Adesso, Quantification of Gaussian Quantum Steering, Phys. Rev. Lett. **114**, 060403 (2015).

[180] Q. Y. He, Q. H. Gong, and M. D. Reid, Classifying Directional Gaussian Entanglement, Einstein-Podolsky-Rosen Steering, and Discord, Phys. Rev. Lett. **114**, 060402 (2015).

[181] X. Deng, Y. Xiang, C. Tian, G. Adesso, Q. He, Q. Gong, X. Su, C. Xie, and K. Peng, Demonstration of Monogamy Relations for Einstein-Podolsky-Rosen Steering in Gaussian Cluster States, Phys. Rev. Lett. **118**, 230501 (2017).

[182] C. Branciard, E. G. Cavalcanti, S. P. Walborn, V. Scarani, and H. M. Wiseman, One-sided device-independent quantum key distribution: Security, feasibility, and the connection with steering, Phys. Rev. A **85**, 010301 (2012).

[183] J. S. Bell, On the einstein podolsky rosen paradox, Physics Physique Fizika **1**, 195 (1964).

[184] J. S. BELL, On the problem of hidden variables in quantum mechanics, Rev. Mod. Phys. **38**, 447 (1966).

[185] N. D. Mermin, Hidden variables and the two theorems of john bell, Rev. Mod. Phys. **65**, 803 (1993).

[186] A. Peres, All the bell inequalities, Found. Phys. **29**, 589 (1999).

[187] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Proposed Experiment to Test Local Hidden-Variable Theories, Phys. Rev. Lett. **23**, 880 (1969).

[188] J. F. Clauser and M. A. Horne, Experimental consequences of objective local theories, Phys. Rev. D **10**, 526 (1974).

[189] S. J. Freedman and J. F. Clauser, Experimental Test of Local Hidden-Variable Theories, Phys. Rev. Lett. **28**, 938 (1972).

[190] A. Aspect, P. Grangier, and G. Roger, Experimental Tests of Realistic Local Theories via Bell's Theorem, Phys. Rev. Lett. **47**, 460 (1981).

[191] A. Aspect, J. Dalibard, and G. Roger, Experimental Test of Bell's Inequalities Using Time-Varying Analyzers, Phys. Rev. Lett. **49**, 1804 (1982).

[192] G. Weihs, T. Jennewein, C. Simon, H. Weinfurter, and A. Zeilinger, Violation of Bell's Inequality under Strict Einstein Locality Conditions, Phys. Rev. Lett. **81**, 5039 (1998).

[193] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberga, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, Loophole-free bell inequality violation using electron spins separated by 1.3 kilometres, Nature **526**, 682 (2015).

[194] M. Giustina et al., Significant-Loophole-Free Test of Bell's Theorem with Entangled Photons, Phys. Rev. Lett. **115**, 250401 (2015).

[195] L. K. Shalm et al., Strong Loophole-Free Test of Local Realism, Phys. Rev. Lett. **115**, 250402 (2015).

[196] C. Abellán et al., T. B. B. T. Collaboration, Challenging local realism with human choices, Nature **557**, 212 (2018).

[197] D. Rauch, J. Handsteiner, A. Hochrainer, J. Gallicchio, A. S. Friedman, C. Leung, B. Liu, L. Bulla, S. Ecker, F. Steinlechner, R. Ursin, B. Hu, D. Leon, C. Benn, A. Gheylene, M. Cecconi, A. H. Guth, D. I. Kaiser, T. Scheidl, and A. Zeilinger, Cosmic Bell Test Using Random Measurement Settings from High-Redshift Quasars, Phys. Rev. Lett. **121**, 080403 (2018).

[198] A. Acín, N. Brunner, N. Gisin, S. Massar, S. Pironio, and V. Scarani, Device-Independent Security of Quantum Cryptography against Collective Attacks, Phys. Rev. Lett. **98**, 230501 (2007).

[199] F. Buscemi, All Entangled Quantum States are Nonlocal, Phys. Rev. Lett. **108**, 200401 (2012).

[200] C. Branciard, D. Rosset, Y.-C. Liang, and N. Gisin, Measurement-Device-Independent Entanglement Witnesses for all Entangled Quantum States, Phys. Rev. Lett. **110**, 060405 (2013).

[201] P. Abiuso, S. Bäuml, D. Cavalcanti, and A. Acín, Measurement-Device-Independent Entanglement Detection for Continuous-Variable Systems, Phys. Rev. Lett. **126**, 190502 (2021).

[202] C. K. Hong and L. Mandel, Experimental Realization of a Localized One-Photon State, Phys. Rev. Lett. **56**, 58 (1986).

[203] G. S. Thekkadath, B. A. Bell, I. A. Walmsley, and A. I., Engineering Schrödinger cat states with a photonic even-parity detector, Quantum **4**, 239 (2020).
[204] M. Eaton, R. Nehra, and O. Pfister, Non-gaussian and gottesman-kitaev-preskill state preparation by photon catalysis, New J. Phys. 21, 113034 (2019).

[205] D. Su, C. R. Myers, and K. K. Sabapathy, Conversion of gaussian states to non-gaussian states using photon-number-resolving detectors, Phys. Rev. A 100, 052301 (2019).

[206] S. Yu, Q. Chen, C. Zhang, C. H. Lai, and C. H. Oh, All Entangled Pure States Violate a Single Bell's Inequality, Phys. Rev. Lett. 109, 120402 (2012).

[207] R. F. Werner and M. M. Wolf, Bound Entangled Gaussian States, Phys. Rev. Lett. 86, 3658 (2001).

[208] Y. Xiang, S. Liu, J. Guo, Q. Gong, N. Treps, Q. He, and M. Walschaers, Quantification of wigner negativity remotely generated via einstein-podolsky-rosen steering, arXiv:2104.00451 [quant-ph] (2021).

[209] M. D. Reid, Monogamy inequalities for the einstein-podolsky-rosen paradox and quantum steering, Phys. Rev. A 88, 062108 (2013).

[210] S.-W. Ji, M. S. Kim, and H. Nha, Quantum steering of multimode gaussian states by gaussian measurements: monogamy relations and the peres conjecture, J. Phys. A: Math. Theor. 48, 135301 (2015).

[211] X. Deng, Y. Xiang, C. Tian, G. Adesso, Q. He, Q. Gong, X. Su, C. Xie, and K. Peng, Demonstration of Monogamy Relations for Einstein-Podolsky-Rosen Steering in Gaussian Cluster States, Phys. Rev. Lett. 118, 230501 (2017).

[212] E. Shchukin and W. Vogel, Inseparability Criteria for Continuous Bipartite Quantum States, Phys. Rev. Lett. 95, 230502 (2005).

[213] J. Sperling and W. Vogel, Verifying continuous-variable entanglement in finite spaces, Phys. Rev. A 79, 052313 (2009).

[214] J. S. Kim and B. C. Sanders, Monogamy of multi-qubit entanglement using rényi entropy, J. Phys. A: Math. Theor. 43, 445305 (2010).

[215] G. Adesso, D. Girolami, and A. Serafini, Measuring Gaussian Quantum Information and Correlations Using the rényi Entropy of Order 2, Phys. Rev. Lett. 109, 190502 (2012).

[216] R. Simon, Peres-Horodecki Separability Criterion for Continuous Variable Systems, Phys. Rev. Lett. 84, 2726 (2000).

[217] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Inseparability Criterion for Continuous Variable Systems, Phys. Rev. Lett. 84, 2722 (2000).

[218] P. van Loock and A. Furusawa, Detecting genuine multipartite continuous-variable entanglement, Phys. Rev. A 67, 052315 (2003).

[219] P. Hyllus and J. Eisert, Optimal entanglement witnesses for continuous-variable systems, New J. Phys. 8, 51 (2006).

[220] S. Gerke, J. Sperling, W. Vogel, Y. Cai, J. Roslund, N. Treps, and C. Fabre, Full Multipartite Entanglement of Frequency-Comb Gaussian States, Phys. Rev. Lett. 114, 050501 (2015).

[221] A. A. Valido, F. Levi, and F. Mintert, Hierarchies of multipartite entanglement for continuous-variable states, Phys. Rev. A 90, 052321 (2014).

[222] Note that there is also a more profound connection between the transpose and time reversal through Wigner's theorem. The latter is here effectively implemented by changing the sign of the momentum variables. However, this connection is beyond the scope of this Tutorial and interested readers are invited to indulge in the literature on quantum chaos instead.

[223] C. H. Bennett, G. Brassard, S. Popescu, B. Schumacher, J. A. Smolin, and W. K. Wootters, Purification of Noisy Entanglement and Faithful Teleportation via Noisy Channels, Phys. Rev. Lett. 76, 722 (1996).

[224] C. H. Bennett, H. J. Bernstein, S. Popescu, and B. Schumacher, Concentrating partial entanglement by local operations, Phys. Rev. A 53, 2046 (1996).

[225] J. Eisert, S. Scheel, and M. B. Plenio, Distilling Gaussian States with Gaussian Operations is Impossible, Phys. Rev. Lett. 89, 137903 (2002).

[226] J. Fiurášek, Gaussian Transformations and Distillation of Entangled Gaussian States, Phys. Rev. Lett. 89, 137904 (2002).

[227] G. Giedke and J. Ignacio Cirac, Characterization of gaussian operations and distillation of gaussian states, Phys. Rev. A 66, 032316 (2002).

[228] R. Dong, M. Lassen, J. Heersink, C. Marquardt, R. Filip, G. Leuchs, and U. L. Andersen, Experimental entanglement distillation of mesoscopic quantum states, Nat. Phys. 4, 919 (2008).

[229] B. Hage, A. Samblowski, J. DiGuglielmo, A. Franzen, J. Fiurášek, and R. Schnabel, Preparation of distilled and purified continuous-variable entangled states, Nat. Phys. 4, 915 (2008).

[230] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Entanglement Purification of Gaussian Continuous Variable Quantum States, Phys. Rev. Lett. 84, 4002 (2000).

[231] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Physical implementation for entanglement purification of gaussian continuous-variable quantum states, Phys. Rev. A 62, 032304 (2000).

[232] J. Fiurášek, Distillation and purification of symmetric entangled gaussian states, Phys. Rev. A 82, 042331 (2010).

[233] E. T. Campbell, M. G. Genoni, and J. Eisert, Continuous-variable entanglement distillation and noncommutative central limit theorems, Phys. Rev. A 87, 042330 (2013).

[234] J. Dias and T. C. Ralph, Quantum repeaters using continuous-variable teleportation, Phys. Rev. A 95, 022312 (2017).

[235] G. Y. Xiang, T. C. Ralph, A. P. Lund, N. Walk, and G. J. Pryde, Heralded noiseless linear amplification and distillation of entanglement, Nat. Photonics 4, 316 (2010).

[236] T. Opatrný, G. Kurizki, and D.-G. Welsch, Improvement on teleportation of continuous variables by photon subtraction via conditional measurement, Phys. Rev. A 61, 032302 (2000).

[237] P. T. Cochrane, T. C. Ralph, and G. J. Milburn, Teleportation improvement by conditional measurements on the two-mode squeezed vacuum, Phys. Rev. A 65, 062306 (2002).
[238] S. Olivares, M. G. A. Paris, and R. Bonifacio, Teleportation improvement by inconclusive photon subtraction, Phys. Rev. A **67**, 032314 (2003).

[239] A. Kitagawa, M. Takeoka, M. Sasaki, and A. Chefles, Entanglement evaluation of non-gaussian states generated by photon subtraction from squeezed states, Phys. Rev. A **73**, 042310 (2006).

[240] Y. Yang and F.-L. Li, Entanglement properties of non-gaussian resources generated via photon subtraction and addition and continuous-variable quantum-teleportation improvement, Phys. Rev. A **80**, 022315 (2009).

[241] C. Navarrete-Benlloch, R. García-Patrón, J. H. Shapiro, and N. J. Cerf, Enhancing quantum entanglement by photon addition and subtraction, Phys. Rev. A **86**, 012328 (2012).

[242] T. Das, R. Prabhu, A. Sen(De), and U. Sen, Superiority of photon subtraction to addition for entanglement in a multimode squeezed vacuum, Phys. Rev. A **93**, 052313 (2016).

[243] M. Walschaers, S. Sarkar, V. Parigi, and N. Treps, Tailing Non-Gaussian Continuous-Variable Graph States, Phys. Rev. Lett. **121**, 220501 (2018).

[244] K. Zhang, J. Jing, N. Treps, and M. Walschaers, Maximal entanglement distillation with single-photon subtraction, arXiv:2103.09197 [quant-ph] (2021).

[245] E. G. Brown, N. Friis, and M. Huber, Passivity and practical work extraction using gaussian operations, New J. Phys. **18**, 113028 (2016).

[246] R. S. Barbosa, T. Douce, P.-E. Emeriau, E. Kashefi, and S. Mansfield, Continuous-variable nonlocality and contextuality, arXiv:1905.08267 [quant-ph] (2019).

[247] A. Acín, N. J. Cerf, A. Ferraro, and J. Niset, Tests of multimode quantum nonlocality with homodyne measurements, Phys. Rev. A **79**, 012112 (2009).

[248] R. García-Patrón, J. Fiúrášek, N. J. Cerf, J. Wenger, R. Tualle-Brouri, and P. Grangier, Proposal for a Loophole-Free Bell Test Using Homodyne Detection, Phys. Rev. Lett. **93**, 130409 (2004).

[249] D. Klyshko, The bell and ghz theorems: A possible three-photon interference experiment and the question of nonlocality, Phys. Lett. A **172**, 399 (1993).

[250] W. N. Plick, F. Arzani, N. Treps, E. Diamanti, and D. Markham, Violating bell inequalities with entangled optical frequency combs and multipixel homodyne detection, Phys. Rev. A **98**, 062101 (2018).

[251] K. Banaszek and K. Wódkiewicz, Nonlocality of the einstein-podolsky-rosen state in the wigner representation, Phys. Rev. A **58**, 4345 (1998).

[252] K. Banaszek and K. Wódkiewicz, Testing Quantum Nonlocality in Phase Space, Phys. Rev. Lett. **82**, 2009 (1999).

[253] S. Aaronson and D. Gottesman, Improved simulation of stabilizer circuits, Phys. Rev. A **70**, 052328 (2004).

[254] V. Veitch, N. Wiebe, C. Ferrie, and J. Emerson, Efficient simulation scheme for a class of quantum optics experiments with non-negative wigner representation, New J. Phys. **15**, 013037 (2013).

[255] S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Sufficient Conditions for Efficient Classical Simulation of Quantum Optics, Phys. Rev. X **6**, 021039 (2016).

[256] M.-D. Choi, Completely positive linear maps on complex matrices, Linear Algebra Appl. **10**, 285 (1975).

[257] A. S. Holevo, The choi-jamiolkowski forms of quantum gaussian channels, J. Math. Phys. **52**, 042202 (2011).

[258] L. García-Álvarez, C. Calcluth, A. Ferraro, and G. Ferrini, Efficient simulatability of continuous-variable circuits with large wigner negativity, Phys. Rev. Res. **2**, 043322 (2020).

[259] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian Boson Sampling, Phys. Rev. Lett. **119**, 170501 (2017).

[260] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Detailed study of gaussian boson sampling, Phys. Rev. A **100**, 032326 (2019).

[261] A. Deshpande, A. Mehta, T. Vincent, N. Quesada, M. Hinsche, M. Ioannou, L. Madsen, J. Lavoie, H. Qi, J. Ewert, D. Hangleiter, B. Fefferman, and I. Dhand, Quantum computational supremacy via high-dimensional gaussian boson sampling, arXiv:2102.12474 [quant-ph] (2021).

[262] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, Theory of Computing **9**, 143 (2013).

[263] N. Quesada, J. M. Arrazola, and N. Killoran, Gaussian boson sampling using threshold detectors, Phys. Rev. A **98**, 062322 (2018).

[264] E. R. Caianiello, On quantum field theory —i: Explicit solution of dyson’s equation in electrodynamics without use of feynman graphs, Il Nuovo Cimento (1943-1954) **10**, 1634 (1953).

[265] K. Brádler, P.-L. Dallaire-Demers, P. Rebentrost, D. Su, and C. Weedbrook, Gaussian boson sampling for perfect matchings of arbitrary graphs, Phys. Rev. A **98**, 032310 (2018).

[266] M. Schuld, K. Brádler, R. Israel, D. Su, and B. Gupt, Measuring the similarity of graphs with a gaussian boson sampler, Phys. Rev. A **101**, 032314 (2020).

[267] T. R. Bromley, J. M. Arrazola, S. Jahangiri, J. Izaac, N. Quesada, A. D. Gran, M. Schuld, J. Swinarton, Z. Zabaneh, and N. Killoran, Applications of near-term photonic quantum computers: Software and algorithms, Quantum Sci. Technol. **5**, 034010 (2020).

[268] H. Qi, D. J. Brod, N. Quesada, and R. García-Patrón, Regimes of Classical Simulability for Noisy Gaussian Boson Sampling, Phys. Rev. Lett. **124**, 100502 (2020).

[269] J. F. Clauser, Experimental distinction between the quantum and classical field-theoretic predictions for the photoelectric effect, Phys. Rev. D **9**, 853 (1974).

[270] H. J. Kimble, M. Dagenais, and L. Mandel, Photon Anti-bunching in Resonance Fluorescence, Phys. Rev. Lett. **39**, 691 (1977).

[271] F. Diedrich and H. Walther, Nonclassical Radiation of a Single Stored ion, Phys. Rev. Lett. **58**, 203 (1987).

[272] R. Schnabel, Squeezed states of light and their applications in laser interferometers, Phys. Rep. **684**, 1 (2017).

[273] A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller, Quantum State Reconstruction of the Single-Photon Fock State, Phys. Rev. Lett. **87**, 050402 (2001).
[274] A. I. Lvovsky and M. G. Raymer, Continuous-variable optical quantum-state tomography, Rev. Mod. Phys. 81, 299 (2009).

[275] A. I. Lvovsky and S. A. Babichev, Synthesis and tomographic characterization of the displaced fock state of light, Phys. Rev. A 66, 011801 (2002).

[276] A. Ourjoumtsev, R. Tualle-Brouri, and P. Grangier, Quantum Homodyne Tomography of a Two-Photon Fock State, Phys. Rev. Lett. 96, 213601 (2006).

[277] M. Cooper, L. J. Wright, C. Söller, and B. J. Smith, Experimental generation of multi-photon fock states, Opt. Express 21, 5309 (2013).

[278] M. Yukawa, K. Miyata, T. Mizuta, H. Yonezawa, P. Marek, R. Filip, and A. Furusawa, Generating superposition of up-to three photons for continuous variable quantum information processing, Opt. Express 21, 5529 (2013).

[279] J. S. Neergaard-Nielsen, B. M. Nielsen, C. Hettich, K. Mølmer, and E. S. Polzik, Generation of a Superposition of odd Photon Number States for Quantum Information Networks, Phys. Rev. Lett. 97, 083604 (2006).

[280] K. Wakui, H. Takahashi, A. Furusawa, and M. Sasaki, Photon subtracted squeezed states generated with periodically poled ktiopo4, Opt. Express 15, 3568 (2007).

[281] D. V. Sychev, A. E. Ulanov, A. A. Pushkina, M. W. Richards, I. A. Fedorov, and A. I. Lvovsky, Enlargement of optical schrödinger's cat states, Nat. Photonics 11, 379 (2017).

[282] J. S. Neergaard-Nielsen, M. Takeuchi, K. Wakui, H. Takahashi, K. Hayasaka, M. Takeoka, and M. Sasaki, Optical Continuous-Variable Qubit, Phys. Rev. Lett. 105, 053602 (2010).

[283] M. Barbieri, N. Spagnolo, M. G. Genoni, F. Ferreyrol, R. Blandino, M. G. A. Paris, P. Grangier, and R. Tualle-Brouri, Non-gaussianity of quantum states: An experimental test on single-photon-added coherent states, Phys. Rev. A 82, 063833 (2010).

[284] A. Ourjoumtsev, H. Jeong, R. Tualle-Brouri, and P. Grangier, Generation of optical 'schrödinger cats' from photon number states, Nature 448, 784 (2007).

[285] J. Etesse, M. Bouillard, B. Kanseri, and R. Tualle-Brouri, Experimental Generation of Squeezed cat States with an Operation Allowing Iterative Growth, Phys. Rev. Lett. 114, 193602 (2015).

[286] B. Hacker, S. Welte, S. Daiss, A. Shaukat, S. Ritter, L. Li, and G. Rempe, Deterministic creation of entangled atom-light schrödinger-cat states, Nat. Photonics 13, 110 (2019).

[287] S. A. Babichev, J. Appel, and A. I. Lvovsky, Homodyne Tomography Characterization and Nonlocality of a Dual-Mode Optical Qubit, Phys. Rev. Lett. 92, 193601 (2004).

[288] M. D'Angelo, A. Zavatta, V. Parigi, and M. Bellini, Tomographic test of bell's inequality for a time-delocalized single photon, Phys. Rev. A 74, 052114 (2006).

[289] H. Takahashi, J. S. Neergaard-Nielsen, M. Takeuchi, M. Takeoka, K. Hayasaka, A. Furusawa, and M. Sasaki, Entanglement distillation from Gaussian input states, Nat. Photon 4, 178 (2010).

[290] A. Ourjoumtsev, A. Dantan, R. Tualle-Brouri, and P. Grangier, Increasing Entanglement between Gaussian

States by Coherent Photon Subtraction, Phys. Rev. Lett. 98, 030502 (2007).

[291] A. Ourjoumtsev, F. Ferreyrol, R. Tualle-Brouri, and P. Grangier, Preparation of non-local superpositions of quasi-classical light states, Nat. Phys. 5, 189 (2009).

[292] H. Jeong, A. Zavatta, M. Kang, S.-W. Lee, L. S. Costanzo, S. Grandi, T. C. Ralph, and M. Bellini, Generation of hybrid entanglement of light, Nat. Photonics 8, 564 (2014).

[293] O. Morin, K. Huang, J. Liu, H. Le Jeannic, C. Fabre, and J. Laurat, Remote creation of hybrid entanglement between particle-like and wave-like optical qubits, Nat. Photonics 8, 570 (2014).

[294] N. Biagi, L. S. Costanzo, M. Bellini, and A. Zavatta, Entangling Macroscopic Light States by Delocalized Photon Addition, Phys. Rev. Lett. 124, 033604 (2020).

[295] T. Serikawa, J.-i. Yoshikawa, S. Takeda, H. Yonezawa, T. C. Ralph, E. H. Huntington, and A. Furusawa, Generation of a cat State in an Optical Sideband, Phys. Rev. Lett. 121, 143602 (2018).

[296] V. A. Averchenko, V. Thiel, and N. Treps, Nonlinear photon subtraction from a multimode quantum field, Phys. Rev. A 89, 063808 (2014).

[297] Y.-S. Ra, C. Jacquard, A. Dufour, C. Fabre, and N. Treps, Tomography of a Mode-Tunable Coherent Single-Photon Subtractor, Phys. Rev. X 7, 031012 (2017).

[298] Y.-S. Ra, A. Dufour, M. Walschaers, C. Jacquard, T. Michel, C. Fabre, and N. Treps, Non-gaussian quantum states of a multimode light field, Nat. Phys. 16, 144 (2020).

[299] V. Cimini, M. Barbieri, N. Treps, M. Walschaers, and V. Parigi, Neural Networks for Detecting Multimode Wigner Negativity, Phys. Rev. Lett. 125, 160504 (2020).

[300] U. Chabaud, F. Grosshans, E. Kashefi, and D. Markham, Efficient verification of boson sampling, arXiv:2006.03520 [quant-ph] (2021).

[301] E. S. Tiunov, V. V. T. (Vyborova), A. E. Ulanov, A. I. Lvovsky, and A. K. Fedorov, Experimental quantum homodyne tomography via machine learning, Optica 7, 448 (2020).

[302] A. E. Lita, A. J. Miller, and S. W. Nam, Counting near-infrared single-photons with 95% efficiency, Opt. Express 16, 3032 (2008).

[303] C. Monroe, D. M. Meekhof, B. E. King, and D. J. Wineland, A "schrödinger cat" superposition state of an atom, Science 272, 1131 (1996).

[304] M. Brune, S. Haroche, J. M. Raimond, L. Davidovich, and N. Zagury, Manipulation of photons in a cavity by dispersive atom-field coupling: Quantum-nondemolition measurements and generation of "Schrödinger cat" states, Phys. Rev. A 45, 5193 (1992).

[305] M. Brune, F. Schmidt-Kaler, A. Maali, J. Dreyer, E. Hagley, J. M. Raimond, and S. Haroche, Quantum Rabi Oscillation: A Direct Test of Field Quantization in a Cavity, Phys. Rev. Lett. 76, 1800 (1996).

[306] G. Nogues, A. Rauschenbeutel, S. Osnaghi, M. Brune, J. M. Raimond, and S. Haroche, Seeing a single photon without destroying it, Nature 400, 239 (1999).
[307] S. Deléglise, I. Dotsenko, C. Sayrin, J. Bernu, M. Brune, J.-M. Raimond, and S. Haroche, Reconstruction of non-classical cavity field states with snapshots of their decoherence, *Nature* **455**, 510 (2008).

[308] M. Hofheinz, E. M. Weig, M. Ansmann, R. C. Bialczak, E. Lucero, M. Neeley, A. D. O'Connell, H. Wang, J. M. Martinis, and A. N. Cleland, Generation of Fock states in a superconducting quantum circuit, *Nature* **454**, 310 (2008).

[309] B. Vlastakis, G. Kirchmair, Z. Leghtas, S. E. Nigg, L. Frunzio, S. M. Girvin, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Deterministically encoding quantum information using 100-photon Schrödinger cat states, *Science* **342**, 607 (2013).

[310] Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge, M. Reagor, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, Confining the state of light to a quantum manifold by engineered two-photon loss, *Science* **347**, 853 (2015).

[311] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Encoding a qubit in a trapped-ion mechanical oscillator, *Nature* **566**, 513 (2019).

[312] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, L. Frunzio, M. Mirrahimi, and M. H. Devoret, Quantum error correction of a qubit encoded in grid states of an oscillator, *Nature* **584**, 368 (2020).

[313] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, *Rev. Mod. Phys.* **86**, 1391 (2014).

[314] U. B. Hoff, J. Kollath-Bönig, J. S. Neergaard-Nielsen, and U. L. Andersen, Measurement-Induced Macroscopic Superposition States in Cavity Optomechanics, *Phys. Rev. Lett.* **117**, 143601 (2016).

[315] M. Brunelli, O. Houhou, D. W. Moore, A. Nunnenkamp, M. Paternostro, and A. Ferraro, Unconditional preparation of nonclassical states via linear-and-quadratic optomechanics, *Phys. Rev. A* **98**, 063801 (2018).

[316] T. A. Palomaki, J. D. Teufel, R. W. Simmonds, and K. W. Lehnert, Entangling mechanical motion with microwave fields, *Science* **342**, 710 (2013).

[317] R. Riedinger, S. Hong, R. A. Norte, J. A. Slater, J. Shang, A. G. Krause, V. Anant, M. Aspelmeyer, and S. Gröblacher, Non-classical correlations between single photons and phonons from a mechanical oscillator, *Nature* **530**, 313 (2016).

[318] A. P. Reed, K. H. Mayer, J. D. Teufel, L. D. Burkhart, W. Pfaff, M. Reagor, L. Sletten, X. Ma, R. J. Schoelkopf, E. Knill, and K. W. Lehnert, Faithful conversion of propagating quantum information to mechanical motion, *Nat. Phys.* **13**, 1163 (2017).

[319] J. D. Jost, J. P. Home, J. M. Amini, D. Hanneke, R. Ozeri, C. Langer, J. J. Bollinger, D. Leibfried, and D. J. Wineland, Entangled mechanical oscillators, *Nature* **459**, 683 (2009).

[320] E. Flurin, N. Roch, F. Mallet, M. H. Devoret, and B. Huard, Generating Entangled Microwave Radiation Over two Transmission Lines, *Phys. Rev. Lett.* **109**, 183901 (2012).

[321] M. Forsch, R. Stockill, A. Wallucks, I. Marinković, C. Gärtner, R. A. Norte, F. van Otten, A. Fiore, K. Srinivasan, and S. Gröblacher, Microwave-to-optics conversion using a mechanical oscillator in its quantum ground state, *Nat. Phys.* **16**, 69 (2020).

[322] R. Lescanne, S. Deléglise, E. Albertinale, U. Réglade, T. Capelle, E. Ivanov, T. Jacqmin, Z. Leghtas, and E. Flurin, Irreversible Qubit-Photon Coupling for the Detection of Itinerant Microwave Photons, *Phys. Rev. X* **10**, 021038 (2020).

[323] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, *Rev. Mod. Phys.* **90**, 035005 (2018).

[324] B. Dubost, M. Koschorreck, M. Napolitano, N. Behbood, R. J. Sewell, and M. W. Mitchell, Efficient Quantification of Non-Gaussian Spin Distributions, *Phys. Rev. Lett.* **108**, 183602 (2012).

[325] F. Haas, J. Volz, R. Gehr, J. Reichel, and J. Estève, Entangled states of more than 40 atoms in an optical fiber cavity, *Science* **344**, 180 (2014).

[326] H. Strobel, W. Muessel, D. Linnemann, T. Zibold, D. B. Hume, L. Pezzè, A. Smerzi, and M. K. Oberthaler, Fisher information and entanglement of non-Gaussian spin states, *Science* **345**, 424 (2014).

[327] R. McConnell, H. Zhang, J. Hu, S. Ćuk, and V. Vuletić, Entanglement with negative Wigner function of almost 3,000 atoms heralded by one photon, *Nature* **519**, 439 (2015).

[328] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Deterministic multi-mode gates on a scalable photonic quantum computing platform, *arXiv:2010.14422 [quant-ph]* (2020).

[329] J. M. Arrazola et al., Quantum circuits with many photons on a programmable nanophotonic chip, *Nature* **591**, 54 (2021).

[330] P. Jouguet, S. Kunz-Jacques, A. Leverrier, P. Grangier, and E. Diamanti, Experimental demonstration of long-distance continuous-variable quantum key distribution, *Nat. Photonics* **7**, 378 (2013).

[331] B. A. Bell, D. Markham, D. A. Herrera-Martí, A. Marin, W. J. Wadsworth, J. G. Rarity, and M. S. Tame, Experimental demonstration of graph-state quantum secret sharing, *Nat. Commun.* **5**, 5480 (2014).

[332] M. Navascués, F. Grosshans, and A. Acín, Optimality of Gaussian Attacks in Continuous-Variable Quantum Cryptography, *Phys. Rev. Lett.* **97**, 190502 (2006).

[333] Y. Guo, Q. Liao, Y. Wang, D. Huang, P. Huang, and G. Zeng, Performance improvement of continuous-variable quantum key distribution with an entangled source in the middle via photon subtraction, *Phys. Rev. A* **95**, 032304 (2017).

[334] J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean, and A. Aspuru-Guzik, Boson sampling for molecular vibronic spectra, *Nat. Photonics* **9**, 615 (2015).

[335] L. Banchi, M. Fingerhuth, T. Babej, C. Ing, and J. M. Arrazola, Molecular docking with Gaussian Boson sampling, *Sci. Adv.* **6** (2020).

[336] S. Jahangiri, J. M. Arrazola, and A. Delgado, Quantum algorithm for simulating single-molecule electron transport, *The Journal of Physical Chemistry Letters*, *J. Phys. Chem. Lett.* **12**, 1256 (2021).

[337] N. Killoran, T. R. Bromley, J. M. Arrazola, M. Schuld, N. Quesada, and S. Lloyd, Continuous-variable and S. Gröblacher, Microwave-to-optics conversion using a mechanical oscillator in its quantum ground state, *Nat. Phys.* **16**, 69 (2020).
quantum neural networks, Phys. Rev. Res. 1, 033063 (2019).

[338] J. Nokkala, R. Martínez-Peña, G. L. Giorgi, V. Parigi, M. C. Soriano, and R. Zambrini, Gaussian states of continuous-variable quantum systems provide universal and versatile reservoir computing, Commun. Physics 4, 53 (2021).

[339] F. Wolf, C. Shi, J. C. Heip, M. Gessner, L. Pezzè, A. Smerzi, M. Schulte, K. Hammerer, and P. O. Schmidt, Motional Fock states for quantum-enhanced amplitude and phase measurements with trapped ions, Nat. Commun. 10, 2929 (2019).

[340] M. Gessner, L. Pezzè, and A. Smerzi, Entanglement and squeezing in continuous-variable systems, Quantum 1, 17 (2017).

[341] B. Yadin, M. Fadel, and M. Gessner, Metrological complementarity reveals the Einstein-Podolsky-Rosen paradox, Nat. Commun. 12, 2410 (2021).

[342] M. Gessner, A. Smerzi, and L. Pezzè, Metrological Non-linear Squeezing Parameter, Phys. Rev. Lett. 122, 090503 (2019).

[343] D. R. M. Arvidsson-Shukur, N. Yunger Halpern, H. V. Lepage, A. A. Lasek, C. H. W. Barnes, and S. Lloyd, Quantum advantage in postselected metrology, Nat. Commun. 11, 3775 (2020).

[344] J. G. Kirkwood, Quantum statistics of almost classical assemblies, Phys. Rev. 44, 31 (1933).

[345] M. Walschaers, N. Treps, B. Sundar, L. D. Carr, and V. Parigi, Emergent complex quantum networks in continuous-variables non-gaussian states, arXiv:2012.15608 [quant-ph] (2021).

[346] H. H. Schaefer, *Topological Vector Spaces* (Springer New York, New York, NY, 1971).

[347] J. B. Conway, *A Course in Functional Analysis* (Springer New York, New York, NY, 2007).