# Non-Gaussian Quantum States and Where to Find Them

Mattia Walschaers

► To cite this version:

Mattia Walschaers. Non-Gaussian Quantum States and Where to Find Them. PRX Quantum, 2021, 2, pp.030204. 10.1103/PRXQuantum.2.030204. hal-03216896

HAL Id: hal-03216896

https://hal.science/hal-03216896

Submitted on 29 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-Gaussian Quantum States and Where to Find Them

Mattia Walschaers*

Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, 4 place Jussieu, Paris F-75252, France

<figure><img src="image_2.png" /></figure>

(Received 2 May 2021; published 28 September 2021)

Gaussian states have played an important role in the physics of continuous-variable quantum systems. They are appealing for the experimental ease with which they can be produced, and for their compact and elegant mathematical description. Nevertheless, many proposed quantum technologies require us to go beyond the realm of Gaussian states and introduce non-Gaussian elements. In this Tutorial, we provide a roadmap for the physics of non-Gaussian quantum states. We introduce the phase-space representations as a framework to describe the different properties of quantum states in continuous-variable systems. We then use this framework in various ways to explore the structure of the state space. We explain how non-Gaussian states can be characterized not only through the negative values of their Wigner function, but also via other properties such as quantum non-Gaussianity and the related stellar rank. For multimode systems, we are naturally confronted with the question of how non-Gaussian properties behave with respect to quantum correlations. To answer this question, we first show how non-Gaussian states can be created by performing measurements on a subset of modes in a Gaussian state. Then, we highlight that these measured modes must be correlated via specific quantum correlations to the remainder of the system to create quantum non-Gaussian or Wigner-negative states. On the other hand, non-Gaussian operations are also shown to enhance or even create quantum correlations. Finally, we demonstrate that Wigner negativity is a requirement to violate Bell inequalities and to achieve a quantum computational advantage. At the end of the Tutorial, we also provide an overview of several experimental realizations of non-Gaussian quantum states in quantum optics and beyond.

DOI: 10.1103/PRXQuantum.2.030204

## CONTENTS

I. INTRODUCTION 2

II. CONTINUOUS-VARIABLE QUANTUM STATES 3

A. Fock space 3

B. Phase space 5

C. Discrete and continuous variables 11

D. Gaussian states 12

III. NON-GAUSSIAN QUANTUM STATES 15

A. Gaussian versus Non-Gaussian 15

B. Examples of non-Gaussian states 18

C. Quantum non-Gaussianity 20

D. Stellar rank 22

E. Wigner negativity 24

IV. CREATING NON-GAUSSIAN STATES 26

A. Deterministic methods 26

B. Conditional methods 27

1. General framework 28

2. An example: photon subtraction 29

V. NON-GAUSSIAN STATES AND QUANTUM CORRELATIONS 32

A. Quantum correlations: a crash course 32

1. Correlations 32

2. Quantum entanglement 33

3. Quantum steering 34

4. Bell nonlocality 36

B. Non-Gaussianity through quantum correlations 36

1. Quantum non-Gaussianity and entanglement 37

2. Wigner negativity and Einstein-Podolsky-Rosen steering 39

C. Quantum correlations through non-Gaussianity 39

1. Entanglement measures on phase space 40

2. Entanglement increase 41

3. Purely non-Gaussian quantum entanglement 43

D. Non-Gaussianity and Bell inequalities 46

<footnote>*mattia.walschaers@lkb.upmc.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<table><tbody><tr><td>VI. NON-GAUSSIAN QUANTUM ADVANTAGES</td><td>48</td></tr><tr><td>VII. EXPERIMENTAL REALIZATIONS</td><td>52</td></tr><tr><td>A. Quantum optics experiments</td><td>52</td></tr><tr><td>B. Other experimental setups</td><td>54</td></tr><tr><td>VIII. CONCLUSIONS AND OUTLOOK</td><td>55</td></tr><tr><td>ACKNOWLEDGMENTS</td><td>57</td></tr><tr><td>APPENDIX: MATHEMATICAL REMARKS</td><td>57</td></tr><tr><td>1. Topological vector spaces</td><td>57</td></tr><tr><td>2. Span</td><td>58</td></tr><tr><td>REFERENCES</td><td>58</td></tr></tbody></table>

# I. INTRODUCTION

Gaussian states have a long history in quantum physics, which dates back to Schrödinger's introduction of the coherent state as a means to study the harmonic oscillator [1]. In later times, Gaussian states rose to prominence due to their importance in the description of Bose gases [2–4] and in the theory of optical coherence [5,6]. With the advent of quantum-information theory, the elegant mathematical structure of Gaussian states made them important objects in the study of continuous-variable (CV) quantum-information theory [7–9]. In this Tutorial, we focus on bosonic systems, which means that the continuous variables of interest are field quadratures. Gaussian quantum states are then defined as the states for which measurement statistics of these field quadratures is Gaussian.

Gaussian states can be fully described by their mean field and covariance matrix, and, due to Williamson's decomposition [10], the latter can be studied using a range of tools from symplectic vector spaces. As such, one can directly relate quadrature squeezing to Gaussian entanglement via the Bloch-Messiah decomposition [11]. In the full state space of CV systems, Gaussian states are furthermore known to play a specific role: of all possible states with the same covariance matrix, the Gaussian state will always have the weakest entanglement [12] and the highest entropy [13]. From a theoretical point of view, Gaussian quantum states provide, thus, an elegant and highly relevant framework for quantum-information theory. On an experimental level, CV quantum information has long been motivated by advances in quantum optics, due to the capability of on-demand generation of ever larger entangled states using either spatial modes [14–16] or time-frequency modes [17–22]. Furthermore, Gaussian states also play a key role in the recent demonstration of a quantum advantage with Gaussian Boson sampling [23]. These developments have made the CV quantum optics an important platform for quantum computation [24].

Regardless of all the experimental and theoretical successes of Gaussian states, they have a major shortcoming in the context of quantum technologies: all Gaussian measurements of such states can be efficiently simulated [25]. In pioneering work on CV quantum computation,

it is already argued that a non-Gaussian operation is necessary to implement a universal quantum computer in CV [26]. Later works that laid the groundwork for CV measurement-based quantum computing have left the question of this non-Gaussian operation somewhat in the open [27–29]. Common schemes, based on the cubic phase gate, turn out to be particularly hard to implement in realistic setups [30]. Furthermore, these protocols require highly non-Gaussian states, such as Gottesman-Kitaev-Preskill (GKP) states [31], to encode information. Even though such states could also serve as a non-Gaussian resource for implementing non-Gaussian gates [32], these states remain notoriously challenging to produce. In spite of the practical problems involved with non-Gaussian states, one is obliged to venture into non-Gaussian territory to reach a quantum computational advantage in the CV regime [33]. This emphasizes the importance of a general understanding of non-Gaussian states and their properties. In this Tutorial, we attempt to provide a roadmap to navigate within this quickly developing field.

In Sec. II, we take an unusual start to introduce CV systems. We first present some elements of many-boson physics, by treating Fock space. This mathematical environment is probably familiar to most readers to describe photons. We then explain how such a Fock space can also be described in phase space, which is the more natural framework from CV quantum optics. We introduce phase-space representations of states and observables in CV systems such as the Wigner function, and to familiarize the reader with the language of multimode systems. By first reviewing the basics of Fock space, we can make interesting connections between what is known as the discrete-variable (DV) approach and the CV approach to quantum optics. We see that there is often a shady region between these two frameworks, where techniques that are typically associated with one framework can be applied in the other. We finally argue that the main distinction lies in whether one measures photons (DV) or field quadratures (CV).

In Sec. III, we provide the reader with an introduction to some of the different structures that can be identified in the space of CV quantum states. When pure states are considered, all non-Gaussian states are known to have a nonpositive Wigner function [34,35], but this no longer holds when mixed states enter the game [36]. In the entirety of the state space, non-Gaussian states occupy such a vast territory that it is impossible to describe all of them within one single formalism. Nevertheless, there has recently been considerable progress in the classification of non-Gaussian states [37,38]. We introduce some key ideas behind quantum non-Gaussianity, the stellar rank, and Wigner negativity as tools to characterize non-Gaussian states.

Section IV introduces two main families of techniques to create non-Gaussian states starting from Gaussian inputs. The first approach concentrates on deterministic
methods, which rely on the implementation of non-Gaussian unitary transformations. We show how such transformations can be built by using a specific non-Gaussian gate. We then introduce the second class of techniques, which are probabilisitic and rely on performing non-Gaussian measurements on a Gaussian state and conditioning on a certain measurement result. We introduce our recently developed approach to describe these systems [39] and present mode-selective photon subtraction as a case study.

Then all the pieces are set to discuss the interplay between non-Gaussian effects and quantum correlations in Sec. V. First, we consider the resources that are required to conditionally prepare certain non-Gaussian states. The conditional scheme relies on performing a non-Gaussian measurement on one part of a bipartite Gaussian state, and will show that the nature of the quantum correlations in this bipartite state is essential. We show that we can only generate quantum non-Gaussian states if the initial bipartite state is entangled. Furthermore, to conditionally generate Wigner negativity we even require quantum steering. In the second part of Sec. V, we show how non-Gaussian operations can in return enhance or create quantum correlations. Finally, we show that Wigner negativity (in either the state or the measurement) is necessary to violate Bell inequalities in CV systems.

In a similar fashion, we spend most of Sec. VI explaining the result of Ref. [33], which shows that Wigner negativity is also necessary to reach a quantum advantage. To show this, we explicitly construct a protocol to efficiently simulate the measurement outcomes of a setup with states, operations, and detectors that are described by positive Wigner functions. In the remainder of the section, we provide comments on the quantum computational advantage reached with Gaussian Boson sampling.

Finally, in Sec. VII, we provide a quick overview of non-Gaussian states in CV experiments. Due to the author's background, the first half of this overview focuses on quantum optics. In the second part, we also discuss some key developments in other branches of experimental quantum physics. Readers should be warned that this is by no means an extensive review of all the relevant experimental progress. A more general conclusion and outlook on what the future may have in store is presented in Sec. VIII.

## II. CONTINUOUS-VARIABLE QUANTUM STATES

Before we can start our endeavor to classify non-Gaussian states of CV systems and study their properties, we must develop some basic formalism for dealing with multimode bosonic systems. At the root of bosonic systems lies the canonical commutation relation, $[\hat{x}, \hat{p}] \sim i\mathbf{1}$, which can be traced back to the early foundations of quantum mechanics. The study of the algebra of such noncommuting observables has given birth to rich branches of

mathematics and mathematical physics that ponder on the subtleties of these observables and their associated states. In this Tutorial, we keep a safe distance from the representation theory of the associated $C^*$ algebras that describe bosonic field theories in their most general sense. We do refer interested readers to a rich but technical literature [4,40–42].

In this Tutorial, we do exclusively work within the Fock representation, which implies that we consider systems with a finite expectation value for the number of particles. In quantum optics, this assumption translates to the logical requirement that energies remain finite. There are many approaches to mathematically construct such systems (luckily for us they are all equivalent [43–46]). Here, we briefly present two such approaches that nicely capture one of the key dualities on quantum physics. First we take the particle approach by introducing the Fock space that describes identical bosonic particles in Sec. A. Subsequently, in Sec. B, we take the approach that starts out from a wave picture, by concentrating on the phase-space representation of the electromagnetic field. Here we also introduce the phase-space representations of CV quantum states that proves to be crucial tools in the remainder of this Tutorial. We show how these approaches are quite naturally two sides of the same coin. In Sec. C, we briefly discuss the concept of modes and the role they play in CV quantum systems. This subsection is both intended to provide some clarification about common jargon and to eliminate common misconceptions. We finish this section by presenting a brief case study of Gaussian states in Sec. D, reviewing some key results. After all, it is difficult to appreciate the subtleties of non-Gaussian states without having a flavor from their Gaussian counterparts.

### A. Fock space

In typical quantum mechanics textbooks, the story of identical particles usually starts by considering a set of $n$ particles, which are each described by a quantum state vector in a single-particle Hilbert space $\mathcal{H}$, thus for the $i$th particle we ascribe a state vector $|\psi_i\rangle \in \mathcal{H}$. The joint state of these $n$ particles is then given by the tensor product of the state vectors $|\psi_1\rangle, \dots, |\psi_n\rangle$. However, if the particles are identical in all their internal degrees of freedom, we should be free to permute them without changing the observed physics. Formally, such permutation is implemented by a unitary operator $U_\sigma$, for the permutation $\sigma \in S_n$, which acts as

$$U_{\sigma} |\psi_{1}\rangle \otimes \dots \otimes |\psi_{n}\rangle = |\psi_{\sigma(1)}\rangle \otimes \dots \otimes |\psi_{\sigma(n)}\rangle. \quad (1)$$

Invariance of physical observables under such permutations can be achieved by either imposing the $n$-particle state vector to be fully symmetric (bosons) or fully antisymmetric (fermions) under these permutations of particles. In this Tutorial, we focus exclusively on bosons, and
thus the condition that must be imposed to obtain a bosonic n-particle state is

$$U_{\sigma} |\Psi^{(n)}\rangle = |\Psi^{(n)}\rangle. \qquad (2)$$

Because these are the only states that are permitted to describe the bosonic system, we commonly use the Hilbert space $\mathcal{H}_s^{(n)}$, which is a subspace of $\mathcal{H}^{\otimes n}$ that contains only those states that fulfil Eq. (2). It is usually convenient to generate these spaces with a set of elementary tensors, known as Fock states, which we define as

$$|\psi_1\rangle \vee \dots \vee |\psi_n\rangle := \sum_{\sigma \in S_n} |\psi_{\sigma(1)}\rangle \otimes \dots |\psi_{\sigma(n)}\rangle, \quad (3)$$

such that

$$\mathcal{H}_s^{(n)} = \overline{\mathrm{span}\{ |\psi_1\rangle \vee \dots \vee |\psi_n\rangle \mid |\psi_i\rangle \in \mathcal{H} \}}, \quad (4)$$

where we refer to the Appendix for some further details on the span. This fully describes a system of n bosonic particles in what is often referred to as first quantization. It is interesting to note that these identical particles appear to be entangled with respect to the tensor product structure of $\mathcal{H}^{\otimes n}$. There is still debate on whether this is a mathematical artefact of our description or rather a genuine physical feature of identical particles. Even though there is still debate about how to exactly define entanglement between indistinguishable particles [47,48], several authors have shown how these symmetrizations can [49–52] induce useful entanglement. Furthermore, it is undeniable that this structure leads to physical interference phenomena that do not exist for distinguishable particles [53].

The name “first quantization” suggests the existence of a second quantization, which turns out to be more appropriate for this Tutorial. Second quantization finds its origins in models where particle numbers are not fixed or conserved. This formalism is largely based on creation and annihilation operators, denoted $\hat{a}^\dagger$ and $\hat{a}$, respectively, that add or remove particles. To accommodate these operators in our mathematical framework, we must equip our Hilbert space to describe a varying number of particles. Therefore, we introduce the Fock space

$$\Gamma(\mathcal{H}) := \mathcal{H}_s^{(0)} \oplus \mathcal{H}_s^{(1)} \oplus \mathcal{H}_s^{(2)} \oplus \dots, \quad (5)$$

where the single-particle Hilbert space is given by $\mathcal{H}_s^{(1)} = \mathcal{H}$. Furthermore, we retrieve a peculiar component $\mathcal{H}_s^{(0)}$, which describes the fraction of the system that contains no particles at all. On its own, $\mathcal{H}_s^{(0)}$ is thus populated by a single state $|0\rangle$ that we refer to as the vacuum. This implies that technically $\mathcal{H}_s^{(0)} \cong \mathbb{C}$ the zero-particle Hilbert space is just described by a complex number that corresponds to the overlap of the state with the vacuum. A general pure state

in Fock space $|\Psi\rangle \in \Gamma(\mathcal{H})$ can then be described using the structure, Eq. (5), as

$$|\Psi\rangle = \Psi^{(0)} \oplus \Psi^{(1)} \oplus \Psi^{(2)} \oplus \dots, \quad (6)$$

where $\Psi^{(i)} \in \mathcal{H}_s^{(i)}$ are non-normalized vectors (and therefore we omit the $|\cdot\rangle$) in the i-particle Hilbert space. Because $|\Psi\rangle$ is a state, we must impose the normalization condition $\|\Psi\|^2 = \sum_{i=0}^{\infty} \|\Psi^{(i)}\|^2 = 1$

We can now define a creation operator $\hat{a}^\dagger(\varphi)$ for every $\varphi \in \mathcal{H}$ [54], which acts as

$$\hat{a}^{\dagger}(\varphi) |\Psi\rangle = 0 \oplus (\Psi^{(0)} |\varphi\rangle) \oplus (|\varphi\rangle \vee \Psi^{(1)}) \oplus (|\varphi\rangle \vee \Psi^{(2)}) \oplus \dots \quad (7)$$

In the same spirit, it is possible to provide an explicit construction of the annihilation operators $\hat{a}(\varphi)$, but here we content ourselves by just introducing the annihilation operator as the hermitian conjugate of the creation operator. Just as the creation operator that literally adds a particle to the system, the annihilation operator literally removes one. One additional property of the annihilation operators is that they destroy the vacuum state:

$$\hat{a}(\varphi) |0\rangle = 0. \quad (8)$$

We can now use creation and annihilation operators to build an arbitrary Fock state by creating particles on the vacuum state

$$|\psi_1\rangle \vee \dots \vee |\psi_n\rangle = \hat{a}^{\dagger}(\psi_1) \hat{a}^{\dagger}(\psi_2) \dots \hat{a}^{\dagger}(\psi_n) |0\rangle \quad (9)$$

and by considering superpositions of such Fock states, we can ultimately generate the entire Fock space. By considering any basis of the single-particle Hilbert space $\mathcal{H}$ and constructing all possible Fock states of all possible lengths that can be formed by generating particles in these basis vectors we construct a basis of the Fock space $\Gamma(\mathcal{H})$. We refer to this basis as the Fock basis.

The beauty of second quantization lies in the natural appearance of states, which have no fixed particle number. The most important example is the coherent state

$$|\alpha\rangle := e^{-\frac{\|\alpha\|^2}{8}} \sum_{j=0}^{\infty} \frac{[\hat{a}^{\dagger}(\alpha)]^j}{2^j j!} |0\rangle, \quad (10)$$

where $\alpha \in \mathcal{H}$ is a non-normalized vector in the single particle Hilbert space. One can, indeed, simply generalise (7) to non-normalized vectors in $\mathcal{H}$ which we use explicitly in Eq. (10). Second, we note that an unusual factor 2 is included to make the definition consistent with Eq. (62).

In quantum optics, these coherent states are crucial objects as they describe perfectly coherent light [5]. It is
important to remark that a coherent state is always generated by a single vector in the single-particle Hilbert space. Coherent states often provide a good approximation for the state that is produced by a single-mode laser far above threshold [55]. More generally, the study of laser light is a whole field in its own right and often the light deviates from the fully coherent approximation.

The creation and annihilation operators are not only important objects because they populate the Fock space; they are also of key importance for describing observables in a many-boson system. These operators are the generators of the algebra of observables that represents the canonical commutation relations on Fock space. This implies that any observable can ultimately be approximated by a polynomial of creation and annihilation operators. At the heart of this mathematical formalism lies the canonical commutation relation (CCR):

$$[\hat{a}(\varphi), \hat{a}^{\dagger}(\psi)] = \langle \varphi | \psi \rangle, \quad (11)$$

which describes the algebra of observables. Note that this relation holds for any vectors $|\varphi\rangle$ and $|\psi\rangle$ in the single-particle Hilbert space $\mathcal{H}$. These vectors should not form a basis, nor should they be orthogonal. When $|\varphi\rangle = |\psi\rangle$, we find that $[\hat{a}(\psi), \hat{a}^{\dagger}(\psi)] = 1$. On the other hand, when the single-particle states $|\varphi\rangle$ and $|\psi\rangle$ are fully orthogonal, we find that $[\hat{a}(\varphi), \hat{a}^{\dagger}(\psi)] = 0$. In these cases we recover the typical creation and annihilation operators for harmonic oscillators. However, by introducing the creations and annihilation operators through Eq. (7), we can also deal with more general cases. Furthermore, all definitions and the form of the CCR are still valid when $\varphi$ and $\psi$ are unnormalized vectors in $\mathcal{H}$. A more detailed discussion can be found in Ref. [56].

When we leave the realm of pure states, the description of quantum states becomes tedious. Commonly, one uses a density operator $\hat{\rho}$ with tr $\hat{\rho} = 1$ to formally describe a state. However, we can generally think of these density operators as infinite-dimensional matrices with an infinite number of components in the Fock basis. In other words, this is not necessarily a convenient description. In an operational sense, any state is considered to be characterized when we know all the moments of all the possible observables. Because the creation and annihilation operators generate the algebra, one knows all the moments of all the observables if one knows all the correlation functions tr[$\hat{\rho}\hat{a}^{\dagger}(\psi_1)\dots\hat{a}^{\dagger}(\psi_n)\hat{a}(\varphi_1)\dots\hat{a}(\varphi_m)$], for all possible lengths $n$ and $m$. Even though this might seem like an equally challenging endeavor, much of quantum statistical mechanics boils down to finding expressions of the correlation functions for relevant classes of states.

## B. Phase space

In the previous subsection, we started our analysis by extending a system of one quantum particle to a system of

many quantum particles. Here we follow a different route, where we start by considering the classical electric field. With some effort, we can apply such an analysis to any bosonic field, but in this Tutorial we focus on quantum optics as our main field of application. For a more extensive introduction from a quantum optics perspective we recommend Refs. [55,57,58], whereas a general introduction to quantum physics in phase space can be found in Ref. [59].

A traveling electromagnetic wave is described by a solution of Maxwell's equations. As is commonly the case in optics, we focus on the complex representation of the electric field, which is generally given by $E^{(+)}(\mathbf{r}, t)$. It is related to the real-valued electric field $E(\mathbf{r}, t)$ that is encountered in standard electrodynamics textbooks by $E(\mathbf{r}, t) = E^{(+)}(\mathbf{r}, t) + [E^{(+)}(\mathbf{r}, t)]^*$. To express the electric field, it is useful to introduce an orthonormal mode basis $\{\mathbf{u}_i(\mathbf{r}, t)\}$. These modes are solutions to Maxwell's equations

$$\nabla \cdot \mathbf{u}_i(\mathbf{r}, t) = 0, \quad (12)$$

$$\left( \Delta - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right) \mathbf{u}_i(\mathbf{r}, t) = 0. \quad (13)$$

The orthogonalization property is implemented by the following condition:

$$\frac{1}{V} \int_V d^3\mathbf{r} [\mathbf{u}_i(\mathbf{r}, t)]^* \mathbf{u}_j(\mathbf{r}, t) = \delta_{ij}, \quad (14)$$

where $V$ is some large volume that contains the entire physical system. This assumption serves the practical purpose of allowing us to consider a discrete mode basis and on top it makes physical sense. Note that we do not integrate over $t$, which implies that at every instant of time $t$ we consider a mode basis that is normalized with respect to the spatial degrees of freedom. It is practical to assume that all relevant physics can be described by a (possibly large) finite number of modes $m$. These modes now form a basis in which we can expand any solution to Maxwell's equations and thus we may write

$$\mathbf{E}^{(+)}(\mathbf{r}, t) = \sum_{j=1}^{m} \mathcal{E}_j \mathbf{u}_j(\mathbf{r}, t), \quad (15)$$

where $\mathcal{E}_j$ are a set of complex numbers, which can be written in terms of the real and imaginary parts

$$\mathcal{E}_j = E_j^{(x)} + iE_j^{(p)}. \quad (16)$$

These real and imaginary parts of the field are known as the amplitude and phase quadrature, respectively. We can interpret these quantities $\vec{E} := (E_1^{(x)}, E_1^{(p)}, \dots, E_m^{(x)}, E_m^{(p)}) \in \mathbb{R}^{2m}$ as the coordinate in optical phase space that describes the light field.
The space of solutions of Maxwell's equations forms a Hilbert space, which we call the mode space $M$ and the mode basis chosen to describe this space is far from unique. As with all Hilbert spaces, we can define unitary transformations and use them to change from one basis to another. As such, let us introduce the unitary operator $U$ to change between bases

$$ \mathbf{u}_i(\mathbf{r}, t) = \sum_{j=1}^{m} U_{ji} \mathbf{v}_j(\mathbf{r}, t), \quad (17) $$

$$ \mathbf{v}_i(\mathbf{r}, t) = \sum_{j=1}^{m} U_{ji}^{\dagger} \mathbf{u}_j(\mathbf{r}, t), \quad (18) $$

where we can in principle obtain $U$ as an infinite-dimensional matrix with

$$ U_{ji} = \frac{1}{V} \int_V d^3\mathbf{r} [\mathbf{v}_j(\mathbf{r}, t)]^* \mathbf{u}_i(\mathbf{r}, t), \quad (19) $$

which remarkably does not depend on time due to the normalization properties of the mode bases. We can analogously expand the electric field in the new mode basis

$$ \mathbf{E}^{(+)}(\mathbf{r}, t) = \sum_{i} \mathcal{E}'_{i} \mathbf{v}_{i}(\mathbf{r}, t), \quad (20) $$

where $\mathcal{E}'_i = \sum_j U_{ij} \mathcal{E}_j$. This observation is of great importance when we quantize the electric field. The change of mode basis also imposes a change of coordinates in the optical phase space. Like in Eq. (16) the new components can also be divided in real and imaginary parts, which leads to a new coordinate $\vec{E}'$. Because the coordinate vectors in optical phase space are real $2m$-dimensional vectors, we obtain

$$ \vec{E}' = O\vec{E}, \quad (21) $$

where $O$ is an orthonormal transformation. However, the orthogonal transformation $O$ on the phase space must correspond to the unitary transformation $U$ on the modes, which imposes the constraint

$$ O_{2i-1,2j-1} = \frac{1}{2}(U_{ij} + U_{ij}^*), \quad (22) $$

$$ O_{2i-1,2j} = -\frac{1}{2i}(U_{ij} - U_{ij}^*), \quad (23) $$

$$ O_{2i,2j-1} = \frac{1}{2i}(U_{ij} - U_{ij}^*), \quad (24) $$

$$ O_{2i,2j} = \frac{1}{2}(U_{ij} + U_{ij}^*). \quad (25) $$

This imposes a symplectic structure to the transformation $O$ such that the optical phase space, just like the phase

space of analytical mechanics, can be treated as a symplectic space. The conserved symplectic structure associated with this space is given by

$$ \Omega = \bigoplus_{j=1}^{m} \omega, \quad \text{with } \omega = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}, \quad (26) $$

such that $O^T \Omega O = \Omega$. Note that $\Omega$ can be interpreted as a matrix representation of the imaginary $i$, in the sense that it has the properties $\Omega^T = -\Omega$ and $\Omega^2 = -1$ [60].

In quantum optics, the electric field of light is treated as a quantum observable $\hat{\mathbf{E}}^{(+)}(\mathbf{r}, t)$. In this quantization, the modes, i.e., the normalized solutions to Maxwell's equations, remain classical objects and all the quantum features are absorbed in the coefficients. We can thus write

$$ \hat{\mathbf{E}}^{(+)}(\mathbf{r}, t) = \sum_{i} \mathcal{E}_{i}^{(1)} \frac{\hat{x}_{i} + i\hat{p}_{i}}{2} \mathbf{u}_{i}(\mathbf{r}, t), \quad (27) $$

where $\mathcal{E}_i^{(1)}$ is a constant that carries the dimensions of the field, which can be interpreted as the electric field of a single photon. Glossing over many subtleties of the quantization of the electromagnetic field, we remind the reader that any system that is described on phase space can be quantized through canonical quantization. The quadrature operators $\hat{x}_j$ and $\hat{p}_k$ therefore follow the canonical commutation relations $[\hat{x}_j, \hat{p}_k] = 2i\delta_{j,k}$, such that they satisfy the Heisenberg relation $\Delta\hat{x}\Delta\hat{p} \ge 1$. As they are introduced above, the quadrature operators are specifically related to the specific mode basis. Indeed, $\hat{x}_j$ and $\hat{p}_j$ are the quadrature operators that describe the field in mode $\mathbf{u}_j(\mathbf{r}, t)$. Thus, when we change the basis of modes, we should change the quadrature operators accordingly in line with Eq. (21). To overcome these difficulties, it is often convenient to introduce a basis-independent expression for the quadrature operators, which can be done by mapping any point in the optical phase space $\vec{f} \in \mathbb{R}^{2m}$ to an observable $\hat{q}(\vec{f})$, given by

$$ \hat{q}(\vec{f}) := \sum_{j=1}^{m} f_{2j-1} \hat{x}_j + f_{2j} \hat{p}_j. \quad (28) $$

These quadrature operators follow a generalized version of the CCR, given by

$$ [\hat{q}(\vec{f}_1), \hat{q}(\vec{f}_2)] = -2i \hat{f}_1^T \Omega \hat{f}_2, \quad \text{for all } \vec{f}_1, \vec{f}_2 \in \mathbb{R}^{2m}. \quad (29) $$

We highlight the particular case where $[\hat{q}(\vec{f}), \hat{q}(\Omega\vec{f})] = 2i\|\vec{f}\|^2$, such that we recover the typical form of the CCR for $\|\vec{f}\| = 1$. This highlights that $\Omega$ maps an amplitude quadrature to its associated phase quadrature. From a mathematical point of view, everything is perfectly well
defined for arbitrary $\vec{f} \in \mathbb{R}^{2m}$ and no normalization conditions have to be imposed. From Eq. (28) we can see that the norm of $\vec{f}$ can be factored out, such that it serves as a general rescaling factor of the quadrature operator. In a physical context, when a quadrature is measured, it is common to renormalize measurements to units of vacuum noise, which practically means that we set $\|\vec{f}\| = 1$. Unless explicitly stated otherwise, we assume that $\|\vec{f}\| = 1$ throughout this Tutorial. We can use these general quadrature operators to express electric field operator as

$$ \hat{\mathbf{E}}^{(+)}(\mathbf{r}, t) = \sum_{j=1}^{m} \mathcal{E}_{j}^{(1)} \frac{\hat{q}(\vec{e}_{j}) + i\hat{q}(\Omega\vec{e}_{j})}{2} \mathbf{u}_{j}(\mathbf{r}, t), \quad (30) $$

and we can use the basis transformation (21) to equivalently express the electric field operator in a different mode basis as

$$ \hat{\mathbf{E}}^{(+)}(\mathbf{r}, t) = \sum_{j=1}^{m} \mathcal{E}_{j}^{(1)} \frac{\hat{q}(\Omega\vec{e}_{j}) + i\hat{q}(\Omega\Omega\vec{e}_{j})}{2} \mathbf{v}_{j}(\mathbf{r}, t). \quad (31) $$

This procedure shows us that optical elements, that change the mode basis, change the associated quadrature operators accordingly.

Equation (30) shows us explicitly that quadrature operators $\hat{q}(\vec{f})$ and $\hat{q}(\Omega\vec{f})$ correspond to the same mode, regardless of the mode basis. This reflects the fact that $\vec{f}$ generates one axis in the optical phase space and $\Omega\vec{f}$ generates the second axis that corresponds to the same mode. As such, any arbitrary mode comes with an associated two-dimensional phase space that mathematically can be denoted as $\text{span}(\vec{f}, \Omega\vec{f})$. Because this phase space is uniquely associated with a specific mode, we introduce the notation

$$ \mathbf{f} = \text{span}(\vec{f}, \Omega\vec{f}), \quad (32) $$

and we refer to this as “mode f.” This allows us to concentrate on the multimode quantum states within this Tutorial, while the specifications of the modes can be left ambiguous. The modes can be seen as the physical implementations of the quantum system and are of major importance in the experimental setting as multimode quantum optics experiments rely on the manipulation of these modes.

Multimode quantum states define expectation values of the field, and when we consider CV quantum optics, we primarily focus on the expectation values of the quadrature operators $\hat{q}(\vec{f})$. These operators are unbounded and have a continuous spectrum. The measurement of a field quadrature thus leads to a continuum of possible outcomes and the continuous-variable approach to quantum optics implies that this characterizes quantum properties of light through the measurement of such quadrature operators.

Formally, we can again describe a quantum state on such a system by a density operator $\hat{\rho}$ but this description is rather inconvenient. It turns out that the quadrature operators $\hat{q}(\vec{f})$ generate the algebra of observables for the quantum system that is comprised within our multimode light. In other words, any observable can be approximated by a polynomial of quadrature operators. This generally implies that we can fully characterize the quantum state $\hat{\rho}$ by correlation functions of the type $\text{tr}[\hat{\rho}\hat{q}(\vec{f}_1)\dots\hat{q}(\vec{f}_n)]$. When we know these correlation functions for all lengths $n$ and normalized vectors in phase space, we have fully characterized the state.

To go beyond the information that is contained in correlation functions, it is often convenient to consider probability distributions as a whole. For a single quadrature $\hat{q}(\vec{f})$ we can introduce the characteristic function for any $\lambda \in \mathbb{R}$ as

$$ \chi(\lambda) = \text{tr}[\hat{\rho} e^{i\lambda\hat{q}(\vec{f})}] = \sum_{n=0}^{\infty} \frac{(i\lambda)^n}{n!} \text{tr}[\hat{\rho}\hat{q}(\vec{f})^n], \quad (33) $$

which is clearly related to the moments $\text{tr}[\hat{\rho}\hat{q}(\vec{f})^n]$. The characteristic function is the Fourier transform of the probability distribution of the outcomes of observable $\hat{q}(\vec{f})$. We can thus obtain the probability distribution as

$$ p(x) = \frac{1}{2\pi} \int_{\mathbb{R}} d\lambda \, \chi(\lambda) e^{-ix\lambda}. \quad (34) $$

This approach can be readily generalized to the joint probability distribution for a set of commuting quadrature operators. We thus consider $\vec{f}_1, \dots, \vec{f}_n$ with $[\hat{q}(\vec{f}_j), \hat{q}(\vec{f}_k)] = 0$ for all $j, k$, and we define for all $\vec{\lambda} = \lambda_1\vec{f}_1 + \lambda_2\vec{f}_2 + \dots + \lambda_n\vec{f}_n$ (note that $\vec{\lambda}$ is not normalized). We can then use the properties of the quadrature operators to construct $\hat{q}(\vec{\lambda}) = \sum_{k=1}^n \lambda_k \hat{q}(\vec{f}_k)$ and define the function

$$ \chi(\vec{\lambda}) = \text{tr}[\hat{\rho} e^{i\hat{q}(\vec{\lambda})}]. \quad (35) $$

This function generates all the correlations between observables $\hat{q}(\vec{f}_1), \dots, \hat{q}(\vec{f}_n)$ and it can be used to obtain the multivariate probability distribution

$$ p(\vec{x}) = \frac{1}{(2\pi)^n} \int_{\mathbb{R}^n} d\vec{\lambda} \, \chi(\vec{\lambda}) e^{-i\vec{\lambda}\cdot\vec{x}}, \quad (36) $$

where $d\vec{\lambda} = d\lambda_1 \dots d\lambda_n$. The function $p(\vec{x})$ describes the probability density to jointly obtain $x_1, \dots, x_n$ as measurement outcomes for the measurements of $\hat{q}(\vec{f}_1), \dots, \hat{q}(\vec{f}_n)$, respectively. This approach relies on the fact that commuting observables can be jointly measured and what we presented to derive Eqs. (34) and (36) is ultimately just classical probability theory. However, not all quadrature operators commute such that joint measurements are not
always possible. This implies that a quantum state cannot be straightforwardly defined by a probability distribution of the optical phase space.

Intriguingly, we can carry out the same procedure for a full multimode system over a set of m modes. To this goal, let us define the quantum characteristic function

$$ \chi : \mathbb{R}^{2m} \to \mathbb{C} : \vec{\lambda} \mapsto \chi(\vec{\lambda}) := \mathrm{tr}[\hat{\rho} e^{i\hat{q}(\vec{\lambda})}]. \quad (37) $$

This function, defined on the full optical phase space, can be used to generate all correlation functions between all quadrature operators. As such, it does characterize the full quantum state, but it is common practice to rather study its inverse Fourier transform, which is known as the Wigner function [61–63]

$$ W(\vec{x}) := \frac{1}{(2\pi)^{2m}} \int_{\mathbb{R}^{2m}} d\vec{\lambda} \chi(\vec{\lambda}) e^{-i\vec{\lambda}^T \vec{x}}. \quad (38) $$

This function has many appealing properties even though it is not a probability distribution but rather a quasiprobability distribution. First of all, the Wigner function is normalized, i.e., $\int_{\mathbb{R}^{2m}} d\vec{x} W(\vec{x}) = 1$. Furthermore, its marginals consistently describe all the joint probability distributions for sets of commuting quadratures in the system. Formally, this implies that $p(\vec{x})$ of Eq. (36) can be obtained by integrating over all the phase-space axes that are not contained within $\mathrm{span}(\vec{f}_1, \dots, \vec{f}_n)$. To do so, let us introduce the n-dimensional vector $\vec{x}_M$ that is associated with the measured quadratures, and the $2m-n$ dimensional vectors $\vec{x}_c$, which describe all other axes in phase space. An arbitrary point in phase space can thus be written as $\vec{x} = \vec{x}_M \oplus \vec{x}_c$. Then we find that

$$ P(\vec{x}_M) = \int_{\mathbb{R}^{2m-n}} d\vec{x}_c W(\vec{x}_M \oplus \vec{x}_c). \quad (39) $$

Finally, the Wigner function also produces the correct expectation values

$$ \int_{\mathbb{R}^{2m}} d\vec{x} \vec{f}_1^T \vec{x} \dots \vec{f}_n^T \vec{x} W(\vec{x}) = \mathrm{Re}\{\mathrm{tr}[\hat{\rho}\hat{q}(\vec{f}_1) \dots \hat{q}(\vec{f}_n)]\}. \quad (40) $$

Note that considering the real part of $\mathrm{tr}[\hat{\rho}\hat{q}(\vec{f}_1) \dots \hat{q}(\vec{f}_n)]$ is essentially equivalent to considering symmetric ordering of the operators. Regardless of these nice properties the Wigner function is by itself not a well-defined probability distribution. Due to complementarity, the function can reach negative values for some states. This Wigner negativity is consistent with the impossibility to jointly describe the measurement statistics of all quadratures while also complying with the laws of quantum physics (notably the Heisenberg relation). The profound relation between negativity of the Wigner functions and joint measurability is perhaps most strikingly illustrated by its connection to contextuality [64]. The formalism of Wigner functions can be

used to construct phase-space representations of arbitrary observables by introducing

$$ \chi_A(\vec{\lambda}) = \mathrm{tr}[\hat{A} e^{i\hat{q}(\vec{\lambda})}], \quad (41) $$

such that the Wigner representation is given by

$$ W_A(\vec{x}) = \frac{1}{(2\pi)^{2m}} \int_{\mathbb{R}^{2m}} d\vec{\lambda} \chi_A(\vec{\lambda}) e^{-i\vec{\lambda}^T \vec{x}}. \quad (42) $$

These Wigner representations have the appealing property that

$$ \mathrm{tr}[\hat{A}\hat{\rho}] = (4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x} W_A^*(\vec{x}) W(\vec{x}). \quad (43) $$

When $\hat{A}$ is an observable and thus has $\hat{A} = \hat{A}^\dagger$, its Wigner function will be real such that $W_A^*(\vec{x}) = W_A(\vec{x})$. However, it may sometimes be useful to extend the formalism to more general operators. As such the entire theory of continuous-variable quantum systems can be developed using Wigner functions.

Several aspects of the phase-space representations in this section are reminiscent of earlier results in Sec. A, which was fully developed in a language of particles (also known as a discrete-variable approach). Indeed, the algebra of operators that is generated by the creation and annihilation operators is actually the same as the algebra generated by the quadrature operators. To formalize this, we must first stress that the optical phase space is isomorphic to an m-dimensional complex Hilbert space, which can equally be interpreted as the single-particle Hilbert space of a photon. Formally, this equivalence is constructed through the bijection (see also the Appendix)

$$ \vec{f} \in \mathbb{R}^{2m} \mapsto \sum_j (f_{2j-1} + if_{2j}) |\varphi_j\rangle \in \mathcal{H}, \quad (44) $$

where $\{|\varphi_j\rangle\}$ is an arbitrary basis of $\mathcal{H}$. We can introduce the operators

$$ \begin{aligned} \hat{a}(\vec{f}) &= \frac{1}{2}[\hat{q}(\vec{f}) + i\hat{q}(\Omega\vec{f})], \\ \hat{a}^\dagger(\vec{f}) &= \frac{1}{2}[\hat{q}(\vec{f}) - i\hat{q}(\Omega\vec{f})]. \end{aligned} \quad (45) $$

By using Eq. (44), we can naturally associate these operators to creation and annihilation operators on the single-particle Hilbert space. We retrieve the canonical commutation relation

$$ [\hat{a}(\vec{f}_1), \hat{a}^\dagger(\vec{f}_2)] = \vec{f}_1^T \vec{f}_2 - i \vec{f}_1^T \Omega \vec{f}_2, \quad (46) $$

which can be connected to the inner product on the Hilbert space $\mathcal{H}$ via Eq. (44).
The definition of creation and annihilation operators allows us to make sense of the vacuum state in our phase-space picture. The vacuum state is completely characterized by the property

$$ \hat{a}(\vec{f}) |0\rangle = 0, \text{ for all } \vec{f} \in \mathbb{R}^{2m}. \quad (47) $$

This simple fact can be used to evaluate the quantum characteristic function

$$ \chi_0(\vec{\lambda}) = \text{tr}(|0\rangle \langle 0| e^{i[\hat{a}^\dagger(\vec{\lambda}) + \hat{a}(\vec{\lambda})]}) \quad (48) $$

$$ = \sum_{n=0}^{\infty} \frac{i^n}{n!} \langle 0 | [\hat{a}^{\dagger}(\vec{\lambda}) + \hat{a}(\vec{\lambda})]^n | 0 \rangle \quad (49) $$

$$ = \sum_{n=0}^{\infty} - \frac{\|\vec{\lambda}\|^2 n}{2^n n!} \quad (50) $$

$$ = \exp \left[ - \frac{\|\vec{\lambda}\|^2}{2} \right]. \quad (51) $$

To obtain Eq. (50) we need a considerable amount of combinatorics to evaluate $\langle 0| [\hat{a}^\dagger(\vec{\lambda}) + \hat{a}(\vec{\lambda})]^n |0\rangle$. In general, it can be shown that $\langle 0| \hat{a}^\dagger(\vec{\lambda}_1) \dots \hat{a}^\dagger(\vec{\lambda}_k) \hat{a}(\vec{\lambda}_{k+1}) \dots \hat{a}(\vec{\lambda}_{k+l}) |0\rangle = 0$. Thus it suffices to cast $[\hat{a}^\dagger(\vec{\lambda}) + \hat{a}(\vec{\lambda})]^n$ in normal ordering and extract the term proportional to identity. Even though straightforward, this calculation is quite cumbersome and thus we do not present the details.

From Eq. (51) the Wigner function can be obtained via an inverse Fourier transformation that leads to

$$ W_0(\vec{x}) = \frac{e^{-\frac{1}{2} \|\vec{x}\|^2}}{2\pi}. \quad (52) $$

This Wigner function describes a Gaussian distribution on the phase space with unit variance along every axis. We can thus use Eq. (40) to see that the vacuum state saturates Heisenberg's inequality, i.e., $\Delta\hat{q}(\vec{f})\Delta\hat{q}(\Omega\vec{f}) = 1$.

The quadrature operators thus generate the same algebra of observables as the creation and annihilation operators. However, both sets of observables tend to cause mathematical problems because they are unbounded operators [42,65–67]. The unboundedness means that, when $|\Psi\rangle$ is contained in the Fock space, there is no guarantee that $\hat{q}(\vec{f}) |\Psi\rangle$ will also be contained in the Fock space. One way of solving this problem explicitly is by only considering states for which $\langle\Psi|\hat{q}(\vec{f})^2|\Psi\rangle < \infty$, such that $\hat{q}(\vec{f}) |\Psi\rangle$ is a well-defined state. Physically this assumption makes sense, as it ultimately implies that we consider only states with finite energies. However, the unboundedness of quadrature operators also disqualifies them as well-defined generators of the $C^*$ algebra of observables (since elements of such algebras must be bounded). $C^*$ algebras are essential tools as they allow reconstruction

of the whole framework of Hilbert spaces based on representation theory of abstract algebras (which is essentially the idea of canonical quantization). A highly formal and detailed treatment that considers all these subtleties for bosonic systems is found in Ref. [41]. The key idea is to rather consider a set of bounded operators that describes the same algebra of observables [65,66] and are known as the displacement operators:

$$ \hat{D}(\vec{\alpha}) = e^{-i\hat{q}(\Omega\vec{\alpha})/2}, \quad (53) $$

where $\vec{\alpha} \in \mathbb{R}^{2m}$ need not be normalized. Again, we can use the isomorphism (44) to identify the displacement operator on the quantized phase space to a displacement operator on the Fock space. These operators can be seen as generators of the quadrature operators and they act in a very natural way on them:

$$ \hat{D}^{\dagger}(\vec{\alpha}) \hat{q}(\vec{f}) \hat{D}(\vec{\alpha}) = \hat{q}(\vec{f}) + \vec{\alpha}^T \vec{f}, \quad (54) $$

which means that the value $\vec{\alpha}^T \vec{f}$ is added to the measurement outcomes of $\hat{q}(\vec{f})$. The displacement operator can be combined according to the rule

$$ \hat{D}(\vec{\alpha}_1) \hat{D}(\vec{\alpha}_2) = \hat{D}(\vec{\alpha}_1 + \vec{\alpha}_2) e^{i \vec{\alpha}_1^T \Omega \vec{\alpha}_2}. \quad (55) $$

This rule is yet another representation of the canonical commutation relation and it generates the same algebra of observables. This implies that any observable $\hat{A}$ can be written as a linear combination of displacement operators. We use the Hilbert-Schmidt inner product $\langle \hat{A}, \hat{B} \rangle_{\text{HS}} = \text{tr}[\hat{A}^\dagger \hat{B}]$ to make this explicit

$$ \begin{aligned} \hat{A} &= \int_{\mathbb{R}^{2m}} d\vec{\lambda} \langle \hat{D}(2\Omega\vec{\lambda}), \hat{A} \rangle_{\text{HS}} \hat{D}(2\Omega\vec{\lambda}), \\ &= \int_{\mathbb{R}^{2m}} d\vec{\lambda} \text{tr}[\hat{A}\hat{D}(-2\Omega\vec{\lambda})]\hat{D}(2\Omega\vec{\lambda}), \end{aligned} \quad (56) $$

and we can readily identify that

$$ \text{tr}[\hat{A}\hat{D}(-2\Omega\vec{\lambda})] = \chi_A^*(\vec{\lambda}). \quad (57) $$

It can then directly be seen that

$$ \text{tr}[\hat{A}\hat{\rho}] = \int_{\mathbb{R}^{2m}} d\vec{\lambda} \chi_A^*(\vec{\lambda}) \text{tr}[\hat{D}(2\Omega\vec{\lambda})\hat{\rho}], \quad (58) $$

$$ = \int_{\mathbb{R}^{2m}} d\vec{\lambda} \chi_A^*(\vec{\lambda}) \chi(\vec{\lambda}). \quad (59) $$

And we immediately obtain Eq. (43) via Plancherel's theorem [67,68].

The displacement operators also implement a unitary operation on a quantum state. This unitary operation has a remarkably simple effect when it is expressed on the level
of the Wigner function. Via the property, Eq. (55), we can calculate that

$$ \hat{\rho} \mapsto \hat{D}(\vec{\alpha}) \hat{\rho} \hat{D}^{\dagger}(\vec{\alpha}) \implies \chi(\vec{\lambda}) \mapsto \chi(\vec{\lambda}) e^{-i\vec{\alpha}^{\mathrm{T}}\vec{\lambda}}. \quad (60) $$

Performing the inverse Fourier transform of these quantum characteristic functions leads to

$$ W(\vec{x}) \underset{\vec{\alpha}}{\stackrel{D(\vec{\alpha})}{\mapsto}} W(\vec{x} - \vec{\alpha}). \quad (61) $$

The displacement operator thus literally implements a displacement of the Wigner function by a vector $\vec{\alpha} \in \mathbb{R}^{2m}$ in phase space.

Displacement operators are also well known as the generators of the coherent states that were introduced in Eq. (10). We can combine the bijection between phase space and Hilbert space, Eq. (44), the expression of creation and annihilation operators in terms of quadratures, Eq. (45), and the definition of the displacement operator, Eq. (53), to derive that

$$ |\vec{\alpha}\rangle = \hat{D}(\vec{\alpha}) |0\rangle. \quad (62) $$

By combining Eqs. (52) and (61), we immediately see that the Wigner function for such a coherent state is given by

$$ W_{\alpha}(\vec{x}) = W_{0}(\vec{x} - \vec{\alpha}) = \frac{e^{-\frac{1}{2}\|\vec{x} - \vec{\alpha}\|^2}}{2\pi}. \quad (63) $$

We emphasize that there is a slight difference between the coherent states as defined here, and coherent states as sometimes introduced in the literature. The difference is a factor of 2, which appears because we normalized the shot noise to 1 rather than to 1/2. As such, our coherent states have an energy in mode $f$, which is given by

$$ \langle\vec{\alpha}|\hat{a}^{\dagger}(\vec{f})a(\vec{f})|\vec{\alpha}\rangle = \frac{1}{4}[(\vec{\alpha}^{\mathrm{T}}\vec{f})^{2} + (\vec{\alpha}^{\mathrm{T}}\Omega\vec{f})^{2}]. \quad (64) $$

The coherent states lead us to two other representations of quantum states and observables: the $Q$ function and $P$ function. The definition of the $P$ function is related to the idea that coherent states form an overcomplete basis of Fock space. This implies, notably, that for a $m$-dimensional single-particle Hilbert space

$$ \frac{1}{(4\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\alpha} |\vec{\alpha}\rangle \langle\vec{\alpha}| = 1, \quad (65) $$

with $1$ the identity operator. We can then show that any observable can be written as [5,6]

$$ \hat{A} = \frac{1}{(4\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\alpha} P_A(\vec{\alpha}) |\vec{\alpha}\rangle \langle\vec{\alpha}|, \quad (66) $$

where we refer to $P_A(\vec{\alpha})$ as the $P$ function of the observable $\hat{A}$. Similarly, we can represent a density of operator $\hat{\rho}$

by its $P$ function $P(\vec{\alpha})$. The reader should be warned that $P$ functions often have rather unpleasant mathematical properties. In particular, they often are not actual functions and can be highly singular.

The $P$ function naturally comes with a dual representation that is known as the $Q$ function. As often in quantum physics, what actually counts is the expectation value of an observable in a specific state. We can use the $P$ function to write

$$ \mathrm{tr}[\hat{A}\hat{\rho}] = \frac{1}{(4\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\alpha} P_A(\vec{\alpha}) \langle\vec{\alpha}|\hat{\rho}|\vec{\alpha}\rangle \quad (67) $$

$$ = \frac{1}{(4\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\alpha} P(\vec{\alpha}) \langle\vec{\alpha}|\hat{A}|\vec{\alpha}\rangle. \quad (68) $$

This naturally introduces the $Q$ function, given by

$$ Q_A(\vec{\alpha}) = \frac{1}{(4\pi)^m} \langle\vec{\alpha}|\hat{A}|\vec{\alpha}\rangle, \quad (69) $$

and, in particular, for the quantum state $\hat{\rho}$ we find that

$$ Q(\vec{\alpha}) = \frac{1}{(4\pi)^m} \langle\vec{\alpha}|\hat{\rho}|\vec{\alpha}\rangle. \quad (70) $$

The latter is of particular interest because it represents the quantum state $\hat{\rho}$ as an actual probability distribution. This leads us to the general identity that

$$ \mathrm{tr}[\hat{A}\hat{B}] = \int_{\mathbb{R}^{2m}} d\vec{\alpha} P_A(\vec{\alpha}) Q_B(\vec{\alpha}) = \int_{\mathbb{R}^{2m}} d\vec{\alpha} Q_A(\vec{\alpha}) P_B(\vec{\alpha}). \quad (71) $$

Thus finishing our introduction to the various descriptions of the quantum states and observables of bosonic many-particle systems.

The $Q$ function has a clear physical interpretation. It is directly proportional to the fidelity of the state $\hat{\rho}$ with respect to a target coherent state $|\vec{\alpha}\rangle$. Furthermore, it is always positive, which implies that it is a well-defined probability distribution. Because we can write $\langle\vec{\alpha}|\hat{\rho}|\vec{\alpha}\rangle = \mathrm{tr}[\hat{\rho}|\vec{\alpha}\rangle\langle\vec{\alpha}|]$, we can use Eqs. (43) and (52) to express the $Q$ function in terms of the Wigner function as

$$ Q(\vec{\alpha}) = \int_{\mathbb{R}^{2m}} d\vec{x} W(\vec{x}) W_0(\vec{x} - \vec{\alpha}). \quad (72) $$

To satisfy both Eqs. (43) and (71), we find that

$$ W(\vec{x}) = \frac{1}{(4\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\alpha} P(\vec{\alpha}) W_0(\vec{x} - \vec{\alpha}), \quad (73) $$

$$ = \frac{1}{(4\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\alpha} P(\vec{\alpha}) \frac{e^{-\frac{1}{2}\|\vec{x} - \vec{\alpha}\|^2}}{(2\pi)^m}. \quad (74) $$
In turn, this implies that

$$Q(\vec{\alpha}) = \frac{1}{(8\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\beta} P(\vec{\beta}) \frac{e^{-\frac{1}{4}\|\vec{\beta}-\vec{\alpha}\|^2}}{(2\pi)^m}. \quad (75)$$

These results thus show that all these phase-space representations are ultimately related to one another through convolution or deconvolution with a Gaussian [recall that the Wigner function of the vacuum Eq. (52) is a Gaussian distribution on phase space]. One can now follow Ref. [69] to define a continuous family of phase-space representations $W_\sigma$ for $\sigma \in [-1, 1]$

$$W_\sigma(\vec{\alpha}) = \left(\frac{1}{4\pi[1-\sigma]}\right)^m \int_{\mathbb{R}^{2m}} d\vec{\beta} P(\vec{\beta}) \frac{e^{-\frac{1}{2[1-\sigma]}\|\vec{\beta}-\vec{\alpha}\|^2}}{(2\pi)^m}, \quad (76)$$

where we convolute the $P$ function with an ever-increasing Gaussian, smoothening its features. We can then see that

$$\mathrm{tr}[\hat{A}\hat{\rho}] = (4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{\alpha} W_{A,-\sigma}(\vec{\alpha}) W_\sigma(\vec{\alpha}). \quad (77)$$

We find, notably, that $W(\vec{x}) = W_{\sigma=0}(\vec{x})$, $Q(\vec{\alpha}) = W_{\sigma=-1}(\vec{\alpha})$, and $P(\vec{\alpha}) = (4\pi)^m W_{\sigma=1}(\vec{\alpha})$. This shows that the phase-space representation becomes more regular when decreasing $\sigma$. Other generalized probability distributions have been considered in the literature [70,71], often to circumvent the unappealing properties of the $P$ function. In this Tutorial, we mainly use the Wigner function and (to a lesser extent) the $Q$ function, as they are suitable tools to classify non-Gaussian quantum states. The $P$ function is often used in the literature to characterize the nonclassicality of a state, where the intuition is that classical light is a mixture of coherent states (and thus its $P$ function is a probability distribution) [5,6].

Before we close this introductory section on the phase-space description of CV quantum systems, we introduce one final tool that often comes in handy. The Wigner function can itself be obtained as the expectation value of an operator [72]. Formally, we write

$$W(\vec{x}) = \frac{1}{(2\pi)^m} \mathrm{tr}[\hat{\rho}\hat{\Delta}(\vec{x})]. \quad (78)$$

Using linearity and Eq. (38), we obtain the special case

$$\hat{\Delta}(\vec{0}) = \frac{1}{(2\pi)^m} \int_{\mathbb{R}^{2m}} d\vec{\lambda} e^{i\hat{q}(\vec{\lambda})}. \quad (79)$$

By using techniques based on Eqs. (54) and (55), we can show that

$$\hat{\Delta}(\vec{0})\hat{q}(\vec{f})\hat{\Delta}(\vec{0}) = -\hat{q}(\vec{f}). \quad (80)$$

This means that $\hat{\Delta}(\vec{0})$ is the parity operator. Its eigenstates are the Fock states, since

$$\hat{\Delta}(\vec{0})a^{\dagger}(\vec{f}_1) \dots a^{\dagger}(\vec{f}_n) |0\rangle = (-1)^n a^{\dagger}(\vec{f}_1) \dots a^{\dagger}(\vec{f}_n) |0\rangle. \quad (81)$$

Thus we can formally identify

$$\hat{\Delta}(\vec{0}) = (-1)^N, \quad (82)$$

where $\hat{N}$ is the number operator. We define this operator by introducing a mode basis $\{\vec{e}_1, \Omega\vec{e}_1, \dots, \vec{e}_m, \Omega\vec{e}_m\}$ of the optical phase space, such that $\hat{N} := \sum_{j=1}^m a^\dagger(\vec{e}_j)a(\vec{e}_j)$. This definition can be combined with the properties of the displacement operator to obtain that

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

Note that what we just obtained is the operator equivalent of a $\delta$ function, which becomes even more explicit when we explicitly write down its Wigner representation

$$W_{\Delta(\vec{x}')}(\vec{x}) = \frac{1}{(4\pi)^m} \delta(\vec{x} - \vec{x}'), \quad (84)$$

which follows directly from Eqs. (43) and (78).

This result may seem somewhat artificial, but it turns out to be extremely useful. The observable $\hat{\Delta}(\vec{x})$ can be measured experimentally by counting photons, which means that the combination of photon counting and displacements directly allows us to reconstruct the Wigner function of the quantum state [73]. Until recently, the lack of good photon-number-resolving detectors in the optical frequency range has long made this method unfeasible for most states. Even though there was an early demonstration of the method for coherent states [74], it is only due to recent developments in detector technologies that the method can be applied to more general states [75,76]. The idea was also applied in other settings [77], and was used in pioneering CV experiments with trapped ions, such as Ref. [78], and in cavity QED [79,80].

### C. Discrete and continuous variables

In Sec. A, we have introduced a many-boson system, regardless of the physical realization of these bosons. Such a many-boson system and its Fock space are built upon the structure that is determined by the single-particle Hilbert space $\mathcal{H}$. The Fock space that is constructed accordingly has a rich structure that is further explored in the Tutorial [56]. In optics, the bosons that we consider are photons, and quantum optics can thus be seen as the theory of a many-boson system in the context of Sec. A. This approach to quantum optics is referred to as the DV approach.
In Sec. B, we contrast this with the CV approach to quantum optics. This approach relies on the measurement of the field quadratures, and can thus be seen as a bosonic quantum field theory. Therefore we started this approach by introducing the classical electric field and its modes, which we subsequently quantized through canonical quantization. We introduced the notion of optical phase space as a general way of describing CV quantum systems. The phase space is directly related to the modes of the field and manipulations of the modes also cause changes in the optical phase space. Nevertheless, any system with a phase space can be described by these techniques.

Both of these approaches are ultimately equivalent. Bosonic creation and annihilation operators describe the same algebra of observables as bosonic quadrature operators, which means that on the level of mathematical structure, both approaches can be interchanged and even mixed. This is strikingly clear when the Wigner function, i.e., the phase-space representation of quantum states and observables that is most naturally associated with field quadratures, turns out to be directly measurable by counting photons. Notably, this implies that when it comes to mathematical structures, bosonic particles such as atoms can also be described on phase space.

The real difference between CV and DV approaches is of an experimental nature. What is important is not the observables that are technically present in the quantum system, but the observables that are practically measured in the lab. For the CV approach we typically use homodyne detection to measure quadratures [7,81], whereas in DV approaches we count photons [82].

A common source of misunderstanding between the DV and CV community stems from the role they attribute to the single-particle Hilbert space and optical phase space, respectively. As we argued, both spaces are (at least for a finite-dimensional number of modes) isomorphic, see the Appendix for some additional mathematical intuition. However, the Hilbert space of a photon, which is inherently a quantum particle, is often interpreted as a quantum object. At the same time, the optical phase space represents the field quadratures of optical modes and is thus rather considered to be a classical object. The origin of this confusion lies in the fact that the optical modes, i.e., normalized solutions of Maxwell's equations, also form a Hilbert space that has its origins entirely in classical physics.

The optical modes are the vessels that contain photons much in the same way as a set of electrons contains spins. The crucial difference is that optical modes are not uniquely defined, we can manipulate them, transform them from one mode basis to another with an interferometer and thus consider new superpositions of modes. In typical experimental settings, one would not consider a superposition of two electrons a new well-defined electron.

Because creation and displacement operators always act in one specific mode (i.e., they are generated by a single vector on the single-photon Hilbert space), single-photon states and coherent states are always single-mode states. We may expand this single mode in a different mode basis, which can even be done physically by sending the state through a beam splitter, to create some form of entanglement in the quantum states. However, this entanglement is just a manifestation of the fact that we are not considering the optimal mode basis. In the CV approach, this has led to the notion of “intrinsic” properties [58], which are those properties of quantum states that are independent of the chosen mode basis. The purity and entropy of a state are notable examples, but one can also introduce a notion of “intrinsic entanglement” to refer to a state that is entangled in any possible mode basis. In the next subsection, we introduce Gaussian states, which will later be shown to never be intrinsically entangled.

### D. Gaussian states

Now that we have introduced phase-space representations for states and observables of CV quantum systems, we still need one building block before we can tackle multimode non-Gaussian states: a good understanding of Gaussian states. It is not the goal of this subsection to delve deep into decades worth of research on Gaussian states. We rather highlight a few key results that set apart Gaussian quantum states from the rest of the vast states' space. For more extended reviews, we refer the reader to Refs. [8,9]. These states are also extensively studied in the mathematical physics literature under the name “quasifree states of the CCR algebra.”

Gaussian states are by definition states that have a Wigner function, which is a Gaussian:

$$W_G(\vec{x}) = \frac{e^{-\frac{1}{2}(\vec{x}-\vec{\xi})^T V^{-1} (\vec{x}-\vec{\xi})}}{(2\pi)^m \sqrt{\det V}}, \quad (85)$$

where $\vec{\xi}$ is referred to as the mean field (or displacement) and $V$ is known as the covariance matrix. With Eq. (40), we can verify that the mean field indeed corresponds to the expectation value of the field quadrature

$$\mathrm{tr}[\hat{\rho}\hat{q}(\vec{f})] = \vec{\xi}^T\vec{f}, \quad (86)$$

similarly, we find for the covariance matrix

$$\mathrm{tr}[\hat{\rho}\hat{q}(\vec{f}_1)\hat{q}(\vec{f}_2)] - \mathrm{tr}[\hat{\rho}\hat{q}(\vec{f}_1)]\mathrm{tr}[\hat{\rho}\hat{q}(\vec{f}_2)] \\ = \vec{f}_1^T V \vec{f}_2 - i \vec{f}_1^T \Omega \vec{f}_2. \quad (87)$$

These quantities can thus be obtained for arbitrary quantum states, but for Gaussian states the covariance matrix and the mean field also determine all higher-order expectation values. The most elegant way to see this is via the
multivariate cumulants (also known as truncated correlation functions), which vanish beyond order two [4]. This fact implies that all properties of Gaussian states can ultimately be deduced from their mean field and—often more importantly—from their covariance matrix.

Hitherto, we have encountered the vacuum state $|0\rangle$ and the coherent states $|\vec{\alpha}\rangle$ as examples of Gaussian states. Both of the examples have a covariance matrix $V = 1$. However, there is a much larger range of possible covariance matrices available and they have to satisfy certain constraints [83]. At first instance, we note that a covariance matrix must be positive. An additional constraint is obtained by imposing that the variance $\Delta^2\hat{q}(\vec{f}) \ge 0$ for all $\vec{f}$ in phase space. Equation (87) then directly yields that $\vec{f}^T(V - i\Omega)\vec{f} \ge 0$, which implies that $V \ge 0$ and suggests that $(V - i\Omega) \ge 0$. However, the latter is not obvious, since $\vec{f}$ are real vectors, whereas $(V - i\Omega)$ is a complex matrix. We thus need an additional ingredient: the Heisenberg inequality. Formally, this inequality can be obtained through Robertson's more general inequality [84], such that we find

$$ \Delta^2 \hat{q}(\vec{f}_1) \Delta^2 \hat{q}(\vec{f}_2) \ge \frac{1}{4} \left| \text{tr}\{\hat{\rho}[\hat{q}(\vec{f}_1), \hat{q}(\vec{f}_2)]\} \right|^2. \quad (88) $$

We can now apply the CCR (29) to obtain the general form

$$ \Delta^2 \hat{q}(\vec{f}_1) \Delta^2 \hat{q}(\vec{f}_2) \ge |\vec{f}_1^T \Omega \vec{f}_2|^2. \quad (89) $$

On the other hand, the definition (87) of the covariance matrix can be used to translate this result to

$$ \vec{f}_1^T V \vec{f}_1 \vec{f}_2^T V \vec{f}_2 \ge |\vec{f}_1^T \Omega \vec{f}_2|^2. \quad (90) $$

This identity can then be used to prove $(\vec{f}_1^T - i\vec{f}_2^T)(V - i\Omega)(\vec{f}_1 + i\vec{f}_2) \ge 0$ for all $\vec{f}_1, \vec{f}_2 \in \mathbb{R}^{2m}$. As a consequence, we find that

$$ V - i\Omega \ge 0, \quad (91) $$

an important constraint on the covariance matrix $V$, which can be understood as combining the positivity conditions and the Heisenberg inequality.

To further understand the structure of covariance matrices and the Gaussian states that they describe, we highlight some important results on symplectic matrices. The first of these results is Williamson's decomposition [10], which states that any positive-definite real matrix $V$ can be diagonalized by a symplectic matrix $S$ (i.e., a matrix with $S^T \Omega S = \Omega$):

$$ V = S^T N S, \text{ with } N = \text{diag}[\nu_1, \nu_1, \nu_2, \nu_2, \dots, \nu_m, \nu_m]. \quad (92) $$

The values $\nu_1, \dots, \nu_m$ are also known as the symplectic spectrum of $V$. From Heisenberg's relation, we then

find the additional constraint that $\nu_1, \dots, \nu_m \ge 1$, in other words, the values in the symplectic spectrum are larger than shot noise. It now becomes straightforward to see the Heisenberg's relation also implies that

$$ \det V \ge 1. \quad (93) $$

It thus becomes apparent that the Heisenberg inequality is saturated when $\det V = 1$. The states for which this is the case must have a covariance matrix $V = S^T S$.

The Gaussian states for which the Heisenberg inequality is saturated turn out to be the pure Gaussian states. Recall that the purity of a quantum state is given by $\mu = \text{tr}[\hat{\rho}^2]$. This quantity can be directly calculated from the Wigner function via Eq. (43). We then find for an arbitrary Gaussian state

$$ \mu_G = (4\pi)^m \int_{\mathbb{R}^{2m}} d\vec{x} \, W_G(\vec{x})^2 = \frac{1}{\sqrt{\det V}}. \quad (94) $$

Alternatively, we may use the symplectic spectrum to express $\mu_G = \prod_{k=1}^m \nu_k^{-1}$. This shows us that a Gaussian state is pure if and only if its covariance matrix is a positive symplectic matrix, i.e., it can be written as $V = S^T S$.

The class of states with a covariance matrix given by $V = S^T S$ is much larger than just the vacuum and coherent states with $V = 1$. The additional states turn out to have asymmetric noise in their quadratures, and because the Heisenberg inequality is saturated this implies that some quadratures have less noise than the vacuum state. The states with such covariance matrices are therefore known as squeezed states. To formalize this intuition, we consider the Bloch-Messiah decomposition (which is known in mathematics and classical mechanics as Euler's decomposition) [11,85]. Any symplectic matrix $S$ can be decomposed as follows:

$$ S = O_1 K O_2, \text{ with } K = \text{diag}[s_1^{1/2}, s_1^{-1/2}, \dots, s_m^{1/2}, s_m^{-1/2}], \quad (95) $$

where $O_1$ and $O_2$ are orthogonal symplectic matrices, i.e., $O_j^T O_j = 1$ and $O_j^T \Omega O_j = \Omega$. We can then see that for any pure Gaussian state, we find

$$ V = S^T S = O^T K^2 O. \quad (96) $$

We have already encountered orthogonal symplectic transformations in Eq. (21), where we associated them with transformations of mode bases. Thus, if we find a set of optical modes that are prepared in a pure Gaussian state, we can always find a different mode basis in which the state
is given by

$$V' = OVO^T = \begin{pmatrix} s_1 & & & & \\ & 1/s_1 & & & \\ & & \ddots & & \\ & & & s_m & \\ & & & & 1/s_m \end{pmatrix}. \quad (97)$$

This means that we can always find a set of symplectic eigenvectors {$\vec{e}_1, \Omega\vec{e}_1, \dots, \vec{e}_m, \Omega\vec{e}_m$} of a pure Gaussian state's covariance matrix, which have the properties that $\Delta^2\hat{q}(\vec{e}_j) = s_j$ and $\Delta^2\hat{q}(\Omega\vec{e}_j) = 1/s_j$, such that the Heisenberg relation is saturated: $\Delta^2\hat{q}(\vec{e}_j)\Delta^2\hat{q}(\Omega\vec{e}_j) = 1$. At the same time, we find that clearly either $\Delta^2\hat{q}(\Omega\vec{e}_j)$ of $\Delta^2\hat{q}(\Omega\vec{e}_j)$ is smaller than one (and thus below shot noise).

Gaussian states naturally come with the notion of Gaussian channels [86], they are the completely positive trace-preserving transformations that map Gaussian states into other Gaussian states. We have already seen that the displacement operators are unitary transformations that fulfil this condition. Because any Gaussian transformation $\Gamma$ preserves the general shape of the Wigner function (85), we can simply describe the Gaussian channel $\Gamma$ in terms of its actions on the mean field and the covariance matrix:

$$V \stackrel{\Gamma}{\mapsto} XVX^T + V_c, \quad (98)$$

$$\vec{\xi} \stackrel{\Gamma}{\mapsto} X\vec{\xi} + \vec{\alpha}. \quad (99)$$

The vector $\vec{\alpha}$ simply serves to displace the entire Gaussian to a different location in phase space. On the level of the covariance matrix, $X$ transforms and reshapes the initial covariance matrix, whereas $V_c$ describes the addition of Gaussian classical noise. Both can a priori be any real matrices, as long as they satisfy the constraint

$$V_c - i\Omega + iX\Omega X^T \ge 0. \quad (100)$$

This constraint derives from the demand that $XVX^T + V_c$ is a well-defined covariance matrix, and therefore $XVX^T + V_c - i\Omega \ge 0$. Because $V$ is a well-defined covariance matrix, $X(V - i\Omega)X^T \ge 0$ and thus it can be seen that $XVX^T + V_c$ is also a well-defined covariance matrix whenever Eq. (100) holds. This simple argument proves Eq. (100) is a sufficient condition for $\Gamma$ to transform the covariance matrix of the initial state into a new bona fide covariance matrix.

An important case is obtained when we impose that $\Gamma$ conserves the purity of the state and is thus a unitary transformation. It then immediately follows that $V_c = 0$, since there cannot be any classical noise. The displacement $\vec{\alpha}$ is

simply implemented by a displacement operator, and the constraint (100), combined with the demand that purity is conserved implies that $X$ is a symplectic matrix. In other words, a Gaussian unitary transformation $\hat{U}_G$ satisfies $V \mapsto S^T VS$. Another relevant example is the case of uniform Gaussian losses, where we set $X = \sqrt{1 - \eta}1$, $V_c = \eta1$ and $\vec{\alpha} = 0$, with the positive value $\eta \le 1$ denoting the amount of loss.

More generally, the action of a Gaussian channel on an arbitrary state can be understood from its action on $\exp[i\hat{q}(\vec{\lambda})]$, which can be proven to take the form

$$\exp[i\hat{q}(\vec{\lambda})] \stackrel{\Gamma}{\mapsto} \exp \left[ i\hat{q}(X^T\vec{\lambda}) + i\vec{\alpha}^T\vec{\lambda} - \frac{1}{2}\vec{\lambda}^T V_c \vec{\lambda} \right]. \quad (101)$$

We can then calculate the quantum characteristic function and use some properties of Fourier transforms to find that the Wigner function transforms as

$$W(\vec{x}) \stackrel{\Gamma}{\mapsto} \int_{\mathbb{R}^{2m}} d\vec{y} \, W(X^{-1}\vec{x} - \vec{y}) \frac{e^{-\frac{1}{2}(\vec{y}-\vec{\alpha})^T V_c^{-1} (\vec{y}-\vec{\alpha})}}{(2\pi)^m \sqrt{\det V_c}}. \quad (102)$$

For Gaussian unitary transformations we find the appealing result that $W(\vec{x}) \mapsto W[S^{-1}(\vec{x} - \vec{\alpha})]$. This means that a Gaussian unitary transformation is simply a coordinate transformation on phase space.

Proving that any completely positive Gaussian channel $\Gamma$ is of the form Eq. (102) with condition (100), is a challenging task. The result was first obtained in Refs. [87,88], using the language of $C^*$ algebras. The proof is rather technical and we do not go into details here.

The paradigm of Gaussian channels is also useful to structure general Gaussian states. One may, for example, wonder which Gaussian channel would transform the vacuum state into the Gaussian state with covariance matrix $V$. In general, there is no unique solution to this question, but there is a straightforward route to find an answer. First, take any symplectic matrix $S$ that satisfies $V - S^T S \ge 0$ (the Williamson decomposition guarantees that this is always possible). This implies that there is a positive-definite matrix $V_c$ such that $V = S^T S + V_c$. As such, an arbitrary Gaussian state can always be decomposed as

$$W_G(\vec{x}) = \int_{\mathbb{R}^{2m}} d\vec{y} \, W_0(S^{-1}\vec{x} - \vec{y}) e^{-\frac{1}{2}(\vec{y}-\vec{\alpha})^T V_c^{-1} (\vec{y}-\vec{\alpha})}. \quad (103)$$

The symplectic operation that is applied to the vacuum is known as multimode squeezing in optics. These transformations are fully equivalent to Bogoliubov transformations that are regularly used in condensed-matter physics [4,8]. Combined with a displacement, this operation provides the most general operation that maps quadrature operators into well-defined quadrature operators.

Now that we have introduced the basic concepts of Gaussian states, we are equipped to start exploring their
non-Gaussian counterparts. Several other important properties of Gaussian states will be introduced along the way to stress just how peculiar these Gaussian states are compared to the rest of state space.

### III. NON-GAUSSIAN QUANTUM STATES

Contrary to Gaussian states with their elegant Wigner function and properties that can be nicely deduced from the covariance matrix, the set of non-Gaussian states is vast and wild. Literally all states with Wigner functions that are not Gaussian are contained within this class. To give an idea of the enormous variety, one can consider that highly exotic states such as Gottesman-Kitaev-Preskill states [31] and Schrödinger cat states inhabit the set of non-Gaussian states together with the states that describe single photons and even certain convex mixtures of Gaussian states. Throughout the years, there have been considerable efforts to structure the set of non-Gaussian states. We introduce the notion of quantum non-Gaussian states [89] and then extend it to a hierarchy based on stellar rank [38]. A different approach is provided by considering that the negativity of the Wigner function can be used as a genuine signature of nonclassicality [90]. However, before we attack these different measures to structure non-Gaussian quantum states, we contrast some properties of Gaussian and non-Gaussian states.

Figure 1 provides an overview that can be used as a brief guide to understand the structure of non-Gaussian states. We attempt to highlight how the different quantities used to structure the non-Gaussian part of state space are interconnected.

#### A. Gaussian versus Non-Gaussian

Gaussian states have many extraordinary properties that set them apart from non-Gaussian states. First of all, pure Gaussian states turn out to be the only quantum states that saturate the uncertainty relation. The easiest way to see this is by describing arbitrary pure states in terms of their wave functions. The wave functions associated with amplitude quadratures $\hat{q}(\vec{f})$ and those associated with phase quadratures $\hat{q}(\Omega\vec{f})$ are related by a Fourier transform. This fact can then be used to show that only Gaussian wave functions saturate the Heisenberg inequality. The extension to arbitrary mixed states can be achieved via Jensen's inequality, which emphasizes that no mixed states can saturate the uncertainty relation. Let us consider a mixed state $\hat{\rho} = \sum_k p_k |\Psi_k\rangle\langle\Psi_k|$ with variances $\Delta^2\hat{q}(\vec{f})$. We also introduce the variances $\Delta_k^2\hat{q}(\vec{f})$ for the pure states $|\Psi_k\rangle$. From Jensen's inequality [91], it follows that

$$ \Delta^2 \hat{q}(\vec{f}) \geq \sum_k p_k \Delta_k^2 \hat{q}(\vec{f}). \quad (104) $$

<figure><img src="image_3.png" /><figcaption>FIG. 1. Overview of the different types of non-Gaussian states that can be found in state space. The different aspects will all be considered throughout Sec. III. Here we attempt to show the stellar hierarchy and how it differentiates itself from the convex hull (mixtures) of Gaussian states. Furthermore, we emphasize that the stellar rank and Wigner negativity are different quantifiers of non-Gaussianity. It should be noted that all non-Gaussian pure states are Wigner-negative states, but we can find states that are not mixtures of Gaussian states without Wigner negativity. For all the classes, we provide examples of states that belong to this group, the Wigner functions for several of these examples are shown in Fig. 2.</figcaption></figure>
For Heisenberg's inequality, we calculate

$$ \begin{aligned} \Delta^2 \hat{q}(\vec{f}) \Delta^2 \hat{q}(\Omega \vec{f}) &\geq \sum_k p_k^2 \Delta_k^2 \hat{q}(\vec{f}) \Delta_k^2 \hat{q}(\Omega \vec{f}) \\ &\quad + \sum_{k \neq l} p_k p_l \Delta_k^2 \hat{q}(\vec{f}) \Delta_l^2 \hat{q}(\Omega \vec{f}) \\ &\geq 1. \end{aligned} \quad (105) $$

The presence of cross terms highlights that even when all the pure states in the mixture saturate the inequality, the mixture does not. The only possible exception is the case where the state is pure.

That only pure Gaussian states saturate the Heisenberg inequality may seem like an innocent observation, but it has an important implication for non-Gaussian states. The Heisenberg inequality can be formulated entirely in terms of the covariance matrix. We showed in Eq. (92) that the inequality is saturated if and only if the covariance matrix is symplectic, i.e., $V = S^T S$. Furthermore, we showed in Eq. (94) that a Gaussian state is pure if and only if its covariance matrix is symplectic $V = S^T S$. The fact that no non-Gaussian states can saturate the inequality thus implies that non-Gaussian states can never have a symplectic covariance matrix $V = S^T S$. This is a first hint of the special role played by Gaussian states.

A more general result along these lines states that for all states $\hat{\rho}$ with the same covariance matrix $V$, the Gaussian state always has the highest von Neumann entropy [13]. First of all, note that entropy $-\text{tr}[\hat{\rho} \log \hat{\rho}]$ is conserved under unitary transformations. Due to the Williamson decomposition (92), we can write any Gaussian state as

$$ \hat{\rho}_G = \hat{U}_G \bigotimes_{j=1}^{m} \hat{\rho}_{\bar{n}_j} \hat{U}_G^{\dagger}, \quad (106) $$

where $\hat{\rho}_{\bar{n}_j}$ is a thermal state of the Hamiltonian $\hat{a}_j^\dagger \hat{a}_j$ with average particle number $\bar{n}_j = (v_j - 1)/2$. From statistical mechanics, we know that thermal states are the quantum states that maximize the von Neumann entropy for a given temperature (here fixed by the occupations $\bar{n}_j$).

It turns out that Gaussian states are limiting cases for many quantities [12]. This result shows that for a range of functionals $f$ on the state space, we find that $f(\hat{\rho}) \ge f(\hat{\rho}_G)$, where $\hat{\rho}_G$ is the Gaussian state with the same covariance matrix as $\hat{\rho}$. Apart from some more technical aspects such as continuity, $f$ must have two important features: it must be conserved under (a certain class of) unitary operations $f(\hat{U}\hat{\rho}\hat{U}^\dagger) = f(\hat{\rho})$ and it must be strongly superadditive $f(\hat{\rho}) \ge f(\hat{\rho}_1) + f(\hat{\rho}_2)$ (note that $\hat{\rho}_1$ and $\hat{\rho}_2$ are marginals of $\hat{\rho}$). The equality must be saturated for product states, i.e., $f(\hat{\rho}_1 \otimes \hat{\rho}_2) = f(\hat{\rho}_1) + f(\hat{\rho}_2)$. For strongly subadditive functions with $f(\hat{\rho}) \le f(\hat{\rho}_1) + f(\hat{\rho}_2)$ the same result implies $f(\hat{\rho}) \le f(\hat{\rho}_G)$, (after all, in that case $-f$

is a strongly superadditive function). It is clear that the von Neumann entropy fulfils the latter conditions and is maximized for Gaussian states. For superadditive entanglement measures, this result can be used to show that for all states with the same covariance matrix, Gaussian states are the least entangled ones (entanglement is much more extensively discussed in Sec. V). However, several common entanglement measures, e.g., the logarithmic negativity [92] and the entanglement of formation [93], are not superadditive.

At the heart of these extremal properties lies the central limit theorem [4,94–96]. There are many versions of the central limit theorem in quantum physics, but we stick to what is probably the simplest one. As always, we consider our optical phase space $\mathbb{R}^{2m}$, but this time, we take $N$ copies of it, which implies that we are dealing with a phase space $\mathbb{R}^{2Nm} = \mathbb{R}^{2m} \oplus \cdots \oplus \mathbb{R}^{2m}$ for the full system. We can then embed a vector $\vec{\lambda} \in \mathbb{R}^{2m}$ in the $j$th of these $N$ copies via $\vec{\lambda}_j := \vec{0} \oplus \cdots \oplus \vec{0} \oplus \vec{\lambda} \oplus \vec{0} \cdots \oplus \vec{0}$ and introduce the new averaged operator

$$ \bar{q}_N(\vec{\lambda}) := \frac{1}{\sqrt{N}} \sum_{j=1}^{N} \hat{q}(\vec{\lambda}_j). \quad (107) $$

It is rather straightforward to see that these observables follow the canonical commutation relation. We can now restrict ourselves to studying the algebra that is generated entirely by such averaged quadrature operators. When we then assume that the different copies of the system are “independently and identically distributed” we must set the overall state to be $\hat{\rho}^{(N)} = \hat{\rho}^{\otimes N}$. We then find the characteristic function of the algebra of averaged observables by

$$ \chi_N(\vec{\lambda}) = \text{tr}[\hat{\rho}^{\otimes N} e^{i\bar{q}_N(\vec{\lambda})}]. \quad (108) $$

The following pointwise convergence can be shown:

$$ \chi_N(\vec{\lambda}) \xrightarrow{N \to \infty} \chi_G(\vec{\lambda}), \quad (109) $$

where $\chi_G(\vec{\lambda})$ is the characteristic function of the Gaussian state $\hat{\rho}_G$ that has the same covariance matrix as $\hat{\rho}$. This means that the non-Gaussian features in any state $\hat{\rho}$ can be coarse grained away by averaging sufficiently many copies of the state. Note that this result considers $N$ copies of an arbitrary $m$-mode state. The single-mode version of this result was proven in Ref. [94], whereas a much more general versions are derived in Refs. [95,96]. In Ref. [12] the central limit theorem is combined with invariance under local unitary transformations to prove the final extremality result, we do not review these points in detail.

The extremality of Gaussian states and the associated central limit theorem highlight why Gaussian states are important in quantum-information theory and quantum
statistical mechanics. It also shows that Gaussian states have some particular properties compared to non-Gaussian states. It thus should not come as a surprise that some of these properties can be used to measure the degree of non-Gaussianity of the state [97–99]. As we mentioned before, for a fixed covariance matrix V the von Neumann entropy is maximized by the Gaussian state. This suggest that we can use the difference in von Neumann entropy as a measure for non-Gaussianity. To formalize things, let us consider an arbitrary state $\hat{\rho}$ with covariance matrix V and mean field $\bar{\xi}$ (this quantities can be derived, respectively, for the second and first moments of the quadrature operators). We then construct a Gaussian state $\hat{\sigma}_V$, which has the same covariance matrix and the mean field. In the spirit of extremality, we then define

$$ \delta(\hat{\rho}) = S(\hat{\sigma}_V) - S(\hat{\rho}), \qquad (110) $$

where $S(\hat{\rho}) := -\text{tr}[\hat{\rho} \log \hat{\rho}]$. Because von Neumann entropy is constant under unitary transformations, for a Gaussian state it depends only on the symplectic spectrum $v_1, \dots, v_m$. In other words, we can calculate $S(\hat{\sigma}_V)$ directly by using the Williamson decomposition (92) on V. We find from Ref. [13] that

$$ S(\hat{\sigma}_V) = \sum_{j=1}^{m} \left[ \frac{v_j + 1}{2} \log \frac{v_j + 1}{2} - \frac{v_j - 1}{2} \log \frac{v_j - 1}{2} \right]. \qquad (111) $$

However, it should be noted that the entropy of the non-Gaussian states $S(\hat{\rho})$ is generally harder to calculate unless we can accurately approximate the state by a finite density matrix in the Fock basis. Furthermore, if the state $\hat{\rho}$ is pure, we simply find that $\delta(\hat{\rho}) = S(\hat{\sigma}_V)$.

Due to extremality of Gaussian states it directly follows that $\delta(\hat{\rho}) \ge 0$, but this does not necessarily mean that $\delta(\hat{\rho})$ is a good measure for non-Gaussianity. References [13,98] establish that

$$ \delta(\hat{\rho}) = S(\hat{\rho} | |\hat{\sigma}_V), \qquad (112) $$

where $S(\hat{\rho} || \sigma_V) := \text{tr}[\hat{\rho}(\log \hat{\rho} - \log \hat{\sigma}_V)]$ is the quantum relative entropy between $\hat{\rho}$ and reference state $\hat{\sigma}_V$. The quantum relative entropy allows us to connect $\delta(\hat{\rho})$ to a range of interesting properties, as shown in Ref. [98]. For example, it directly follows that $\delta(\hat{\rho}) = 0$ if and only if $\hat{\rho} = \hat{\sigma}_V$. Furthermore, the measure $\delta(\hat{\rho})$ inherits convexity and monotonicity properties from the relative entropy. These are exactly the properties that made this measure a useful ingredient in the resource theory for quantum non-Gaussianity presented in Ref. [100].

Thus, the connection between Eq. (110) and relative entropy shows that $\delta(\hat{\rho})$ can indeed be used as a measure for non-Gaussianity in the sense that it measures “entropic

distance” between $\hat{\rho}$ and $\hat{\sigma}_V$. Yet, there is one important question that remains to be answered: is $\sigma_V$ indeed the closest Gaussian state to $\hat{\rho}$? An affirmative answer to this question was provided in Ref. [101], where it was shown that

$$ \delta(\hat{\rho}) = \min_{\hat{\rho}_G} S(\hat{\rho} || |\hat{\rho}_G|), \qquad (113) $$

where we minimise over all possible Gaussian states $\hat{\rho}_G$. The main idea of the proof is to show that $S(\hat{\rho} || \hat{\rho}_G) - S(\hat{\rho} || \hat{\sigma}_V) = S(\hat{\sigma}_V || \hat{\rho}_G) \ge 0$ such that the smallest relative entropy in indeed achieved for $\hat{\sigma}_V$. For the technical details, we refer the interested reader to Ref. [101]. Furthermore, we note that a similar non-Gaussianity measure was introduced by using the Wehrl entropy (based on the Q function) rather than the von Neumann entropy [102].

We have thus shown that Gaussian states are special in the sense that they minimize entanglement and maximize entropy as compared to other states with the same covariance matrix. Another profound distinction can be found when comparing pure Gaussian states to pure non-Gaussian states. In this case, there is a seminal result by Hudson [34] that was extended by Soto and Claverie to multimode systems [35], which states that a pure state can have a non-negative Wigner function if and only if the state is Gaussian. In other words, all non-Gaussian pure states exhibit Wigner negativity.

Here, we follow the approach of Ref. [103] to prove this result. First of all, we introduce the function

$$ F_{\Psi}^{*}(\vec{\alpha}) := \langle \vec{\alpha} | \Psi \rangle e^{\frac{1}{8} \|\vec{\alpha}\|^2}, \qquad (114) $$

such that the Q function (70) of the state $|\Psi\rangle$ is given by

$$ Q(\vec{\alpha}) = \frac{1}{(4\pi)^m} |F_{\Psi}^{*}(\vec{\alpha})|^2 e^{-\frac{1}{4} \|\vec{\alpha}\|^2}. \qquad (115) $$

From Eq. (114), we directly find that

$$ |F_{\Psi}^{*}(\vec{\alpha})|^2 \le e^{\frac{1}{4} \|\vec{\alpha}\|^2}. \qquad (116) $$

Next, we observe that a Q function that reaches zero implies a negative Wigner function, which can be seen from Eq. (72). Thus, demanding that the state has a positive Wigner function implies demanding that $Q(\vec{\alpha}) > 0$, and thus that $F_{\Psi}^{*}(\vec{\alpha})$ has no zeros. Using the equivalence between 2m-dimensional phase space and a complex m-dimensional Hilbert space, we can use the multidimensional but restricted version of the Hadamard theorem [35], which states that any entire function $f: \mathbb{C}^m \to \mathbb{C}$ without any zeros and with order of growth [104] r is an exponential $f(z) = \exp g(z)$, where g(z) is a polynomial of degree $s \le r$. We then note that $F_{\Psi}^{*}$ is an entire function of maximal growth $r = 2$ as given by Eq. (116).
Hadamard's theorem then tells us that $F_{\psi}^{*}(\vec{\alpha})$ is Gaussian. In other words, if a pure quantum state $|\Psi\rangle$ has a positive Wigner function $F_{\psi}^{*}(\vec{\alpha})$ must be Gaussian. The only states for which this is the case are Gaussian states.

In summary, we have seen that Gaussian states inherit a particular extremal behavior from the central limit theorem. This should not come as a surprise given that they are the Gibbs states of a free bosonic field at finite temperature. Thus, a non-Gaussian state can be expected to have more “exotic” features than the Gaussian state with the same covariance matrix. This also formalizes the intuition that Gaussian states are more classical states. This idea is further established by the fact that all pure Gaussian states are the only possible pure states that have a positive Wigner function.

The fact that non-Gaussian states automatically have nonpositive Wigner functions no longer holds when mixed states are considered. A simple example is that state $\hat{\rho} = [0|0\rangle\langle0| + a^{\dagger}(\vec{f})|0\rangle\langle0|a(\vec{f})]/2$, which is clearly non-Gaussian but also has a positive Wigner function. There have been considerable efforts to extend Hudson's theorem in some form to mixed states [36]. However, in what follows, we see that there are many ways for a state to be non-Gaussian. This makes it particularly hard to connect a measure such as Eq. (110) to more operational interpretations. In the next section, we start by showing some examples of different non-Gaussian states to make the reader appreciate their variety.

## B. Examples of non-Gaussian states

An overview of the different examples discussed in this section is shown in Fig. 2.

A first important class of non-Gaussian states are Fock states, generated by acting with creation operators $a^{\dagger}(\vec{f})$ on the vacuum state

$$|n_f\rangle := \frac{1}{\sqrt{n!}} [\hat{a}^{\dagger}(\vec{f})]^n |0\rangle, \qquad (117)$$

which is a state of $n$ photons in mode $f$. These states are inherently single mode, even though they can be embedded in a much larger multimode state space. The Wigner function for such states is commonly found in quantum optics textbooks, but deriving it using Eq. (38) is a good exercise. Here we simply state the result:

$$W_{n_f}(\vec{x}) = \sum_{k=0}^{n} \binom{n}{k} \frac{(-1)^{n+k} \|\vec{x}_f\|^2 k}{k!} \frac{e^{-\frac{1}{2} \|\vec{x}\|^2}}{2\pi}. \qquad (118)$$

In the most general sense, we write $\vec{x}_f = (\vec{f}^T \vec{x})\vec{f} + (\vec{f}^T \Omega \vec{x})\Omega \vec{f}$ as the projector of $\vec{x}$ on the phase space of mode $f$, but it is most practical to use the coordinate representation $\vec{x}_f = (x_f, p_f)^T$ where $x_f = \vec{f}^T \vec{x}$ and $p_f = \vec{f}^T \Omega \vec{x}$ to describe the two-dimensional phase space associated with mode $f$.

Multimode Fock states can be obtained by acting on the vacuum with different creation operators in different modes. Even though these different modes do not necessarily have to be orthogonal [56], we here focus on the case where they are. For example, in the $m$-mode system, we can choose a basis $\{\vec{e}_1, \Omega \vec{e}_1, \dots, \vec{e}_m, \Omega \vec{e}_m\}$ of the phase

<figure><img src="image_4.png" /><figcaption>FIG. 2. Several examples of Wigner functions for single-mode non-Gaussian states. The Wigner functions for the Fock states are obtained from Eq. (118) and the mixture of a Fock state and a vacuum is given by Eq. (124). The mixture of coherent states is expressed in Eq. (122) where we set $\|\vec{\alpha}\| = 4$. The Wigner function for the cat state is given by Eq. (126) where we choose $\|\vec{\alpha}\| = 6$. Finally, for the GKP state, Eq. (130), we numerically integrated a wave-function expression of the state with $s = 2$ and $\delta = 0.3$ to obtain the Wigner function.</figcaption></figure>
space $\mathbb{R}^{2m}$ and define multimode Fock states as

$$|n_{e_1}\rangle \otimes \dots \otimes |n_{e_m}\rangle \\ := \frac{1}{\sqrt{n_1! \dots n_m!}} [\hat{a}^\dagger(\vec{e}_1)]^{n_1} \dots [\hat{a}^\dagger(\vec{e}_m)]^{n_m} |0\rangle, \quad (119)$$

where the kth mode in the basis contains $n_k$ photons. Note that for $n_k = 0$ we have a vacuum mode. For the Wigner function, this implies

$$W_{n_{e_1}, \dots, n_{e_m}}(\vec{x}) = W_{n_{e_1}}(\vec{x}_{e_1}) \dots W_{n_{e_m}}(\vec{x}_{e_m}). \quad (120)$$

Here we use that the phase-space point $\vec{x}$ can be expressed as $\vec{x} = \vec{x}_{e_1} \oplus \dots \oplus \vec{x}_{e_m}$, where $\vec{x}_{e_k}$ is the phase-space coordinate within the subspace spanned by $\vec{e}_k$ and $\Omega \vec{e}_k$. We can note $\vec{x}_{e_k} = (x_k, p_k)^T$, such that we find the coordinate representation $\vec{x} = (x_1, p_1, \dots, x_m, p_m)^T$ in the chosen basis of phase space.

Non-Gaussian states do not necessarily have to be pure, they can also come in the form of statistical mixtures. The most basic example of such a state is a non-Gaussian mixture of Gaussian states. As a simple example, let is consider a mixture of two coherent states

$$\hat{\rho} = \frac{1}{2} (|\vec{\alpha}\rangle \langle\vec{\alpha}| + |-\vec{\alpha}\rangle \langle-\vec{\alpha}|). \quad (121)$$

Even though this is a highly classical state, it is still non-Gaussian as clearly seen from its Wigner function

$$W(\vec{x}) = \frac{1}{4\pi} \left( e^{-\frac{1}{2} \|\vec{x}-\vec{\alpha}\|^2} + e^{-\frac{1}{2} \|\vec{x}+\vec{\alpha}\|^2} \right). \quad (122)$$

Another important example of a non-Gaussian mixed state is

$$\hat{\rho}_{\lambda} = \lambda |0\rangle \langle 0| + (1-\lambda) |1_f\rangle \langle 1_f|. \quad (123)$$

In the mode $\vec{f}$, the Wigner function of the state behaves as

$$W_{\lambda}(\vec{x}) = [(1 - \lambda)\|\vec{x}_{\mathrm{f}}\|^2 + 2\lambda - 1] \frac{e^{-\frac{1}{2}\|\vec{x}\|^2}}{2\pi}. \quad (124)$$

This Wigner function reaches negative values as long as $\lambda < 1/2$ and subsequently becomes positive. Nevertheless, it will remain non-Gaussian until $\lambda = 1$. As we see in Sec. C, even when the Wigner function is positive, it is not always possible to describe this state as a mixture of Gaussian states.

Generally speaking, non-Gaussian states can come in a wide variety of shapes, which can be much more exotic than the examples discussed above. A popular class of Gaussian states is obtained by taking coherent superpositions of Gaussian states. As we see in Sec. VII, there has been a strong experimental focus on two specific types of

such states: Schrödinger's cat states [105] and GKP states [31]. The former are obtained by coherently superposing two coherent states, and often are split in even $|\mathrm{cat}_+\rangle$ and odd $|\mathrm{cat}_-\rangle$ cat states:

$$|\mathrm{cat}_{\pm}\rangle := \frac{1}{\mathcal{N}} (|\vec{\alpha}\rangle \pm |-\vec{\alpha}\rangle), \quad (125)$$

where $\mathcal{N} = \sqrt{2(1 \pm \exp[-\|\vec{\alpha}\|^2])}$ is the normalization coefficient, which depends on the displacement $\vec{\alpha}$. The latter is often referred to as “the size of the cat.” The Wigner function of these states resembles that of Eq. (122) but has an additional interference term

$$W_{\mathrm{cat}_{\pm}}(\vec{x}) = \frac{e^{-\frac{1}{2}\|\vec{x}-\vec{\alpha}\|^2} + e^{-\frac{1}{2}\|\vec{x}+\vec{\alpha}\|^2} \pm \cos(\sqrt{2}\vec{\alpha}^T\vec{x})e^{-\frac{1}{2}\|\vec{x}\|^2}}{4\pi(1 \pm e^{-\|\vec{\alpha}\|^2})}. \quad (126)$$

The appearance of these interference terms creates several regions in phase space where the Wigner function attains negative values. The term “Schrödinger’s cat state” has historically grown from the idea that coherent states describe classical electromagnetic fields and can thus be considered “macroscopic,” in particular, for large values of $\|\vec{\alpha}\|$. However, one should honestly admit that they fail to capture an important point of Schrödinger’s thought experiment [106]: the entanglement with a microscopic quantum system (i.e., the decay event that triggers the smashing of the vial of poison). Nevertheless, the term “cat state” has established itself firmly in the CV jargon, and now also lies at the basis of derived concepts such “cat codes” for error correction [107,108].

Finally, there are the GKP states. In their idealized form, they rely on eigenvectors of the quadrature operators (sometimes also known as infinitely squeezed states). To keep notation simple, we restrict to the single mode with quadrature operators $\hat{x}$ and $\hat{p}$. The eigenvectors of these operators are then formally written as

$$\hat{x} |x\rangle = x |x\rangle, \quad \text{and} \quad \hat{p} |p\rangle = p |p\rangle. \quad (127)$$

Furthermore, we have the relations $\langle x' | x \rangle = \delta(x' - x)$, $\langle p' | p \rangle = \delta(p' - p)$, and $\langle p | x \rangle = e^{-ipx}/\sqrt{2\pi}$. GKP states are constructed by considering a grid of such states to create a qubit, by identifying the two following GKP vectors:

$$|\mathrm{GKP}_0\rangle := \sum_{k \in \mathbb{Z}} |x = 2k\sqrt{\pi}\rangle, \quad (128)$$

$$|\mathrm{GKP}_1\rangle := \sum_{k \in \mathbb{Z}} |x = (2k + 1)\sqrt{\pi}\rangle. \quad (129)$$

Clearly, these states are not normalizable and not physical as they would require infinite energy to be created. Thus, it is common to construct approximate GKP states,
by replacing the states $|x\rangle$ with displaced squeezed states, and by truncating the summation by adding a Gaussian envelope:

$$|GKP_0\rangle := \mathcal{N}_0 \sum_{k \in \mathbb{Z}} e^{-2\pi[k\delta]^2} \hat{D}[2k\sqrt{\pi}] |s\rangle, \quad (130)$$

$$|GKP_1\rangle := \mathcal{N}_1 \sum_{k \in \mathbb{Z}} e^{-2\pi[(k+1/2)\delta]^2} \hat{D}[(2k+1)\sqrt{\pi}] |s\rangle. \quad (131)$$

Here $\mathcal{N}_{0,1}$ are normalization constants and $|s\rangle$ is a single-mode squeezed vacuum, which implies that its Wigner function is given by Eq. (85) with $\tilde{\xi} = (0,0)^T$ and $V = \text{diag}[1/s, s]$ for $s > 1$. To get a good GKP state for quantum error correction, we generally need that $s \gg \sqrt{\pi}$. The Wigner function of an ideal GKP state is a grid of delta functions, which again highlights that it is a nonphysical state. The more realistic states $|GKP_0\rangle$ and $|GKP_1\rangle$ have well-defined Wigner functions, even though they are not very insightful to write down explicitly. In Fig. 2, we plot an example that was calculated numerically by taking into account only the first few terms around $k=0$ in the sum.

GKP states may seem a little artificial at first glance, but they have been developed with a very clear purpose: to encode a qubit in a harmonic oscillator [31]. This encoding implies a notion of fault tolerance as these states are designed to be very efficient at correcting displacement errors. The more realistic incarnations of these states, Eq. (130), are therefore often proposed as candidates for encoding the information in CV quantum computation protocols [29]. Furthermore, it was shown that these states can also be used as the sole non-Gaussian resource to implement a CV quantum computer [32].

Once we progress into the realm of multimode states, the class of non-Gaussian states becomes even more vast. In Sec. 2, we present a dedicated introduction to multimode photon-subtracted states, which is a useful state to illustrate several of the concepts treated in this Tutorial. Furthermore, these states have a particular importance in CV quantum optics experiments. As a final example, we introduce another class of multimode non-Gaussian states, which have been highly relevant for quantum metrology: $N00N$ states [109–111]. Even though these states are very promising for quantum sensing with optical setups, the general idea that underlies these states was first introduced for fermions [112] in an attempt to mimic the advantage that is provided by squeezing in optics.

$N00N$ states are two-mode entangled states defined in a pair of orthogonal modes $g_1$ and $g_2$. The state contains exactly $N$ photons, and is a superposition of a state with all photons being mode $g_1$ and a state with all photons in

mode $g_2$:

$$|N00N\rangle := \frac{1}{\sqrt{2}} (|N_{g_1}\rangle + |N_{g_2}\rangle). \quad (132)$$

Here we recall that the state $|N_{g_1}\rangle$ can be trivially embedded to the full multimode space by adding vacuum in all other modes. We study these states in more detail for $N=2$ in our discussion of the Hong-Ou-Mandel effect surrounding Eq. (165). However, here we highlight already that the Wigner function of $|N00N\rangle$ is not simply the sum of Wigner functions of the form Eq. (118). The entanglement will create additional interference terms, just like we saw in Eq. (126). In the present case, these interferences are genuinely multimode, and thus related to quantum correlations.

Experimentally, these states have been created and analyzed using a DV approach [113]. As we highlighted in Sec. C, the distinction between DV and CV is somewhat subtle and mainly depends on what is measured. Because $N00N$ states are built from Fock states and have a well-defined total photon number, they are most natural to analyze using photon-number-resolving detectors.

### C. Quantum non-Gaussianity

Non-Gaussian states come in a wide variety, which means also that some of them are more exotic than others. Non-Gaussian states that are of limited interest, are those which are convex combinations of Gaussian states. Gaussian states do not form a convex set, after all, we can immediately see that, e.g., $[W_0(\vec{x} - \vec{\alpha}_1) + W_0(\vec{x} - \vec{\alpha}_2)]/2$ is not a Gaussian function even though it is a convex combination of Gaussian states.

The fact that the set of non-Gaussian states contains mixtures of Gaussian states may lead one to suspect that any mixed state with a positive Wigner function can be written as a well-chosen mixture of Gaussian states. After all, Gaussian states are the only pure states with positive Wigner functions. This intuition turns out to be false [89], which means that the set of states with a nonpositive Wigner function is not the same as the set of states that lie outside of the convex hull of Gaussian states $\mathcal{G}$. More formally, let us define

$$\mathcal{G} := \{ \hat{\rho} | \hat{\rho} = \int d\gamma p(\gamma) \hat{\rho}_G(\gamma) \}, \quad (133)$$

where $\gamma$ is some arbitrary way of labelling Gaussian states $\hat{\rho}_G(\gamma)$ and $p(\gamma)$ is a probability distribution on these labels. Note that Eq. (103) tells us that we can generate all Gaussian states by taking convex combinations of displaced squeezed states, and thus we can limit ourselves to $\hat{\rho}_G(\gamma) = |\Psi_G(\gamma)\rangle\langle\Psi_G(\gamma)|$ in the definition of $\mathcal{G}$.

Any quantum state $\hat{\rho}$ that is not contained in the convex hull of Gaussian states, i.e., $\hat{\rho} \notin \mathcal{G}$, is referred to
as a “quantum non-Gaussian” state. The intuition behind this terminology is that Gaussian pure states are less quantum than non-Gaussian pure states that boast a non-positive Wigner function. A mixed state that is quantum non-Gaussian may have a positive Wigner function, but it cannot be created without adding states with nonpositive Wigner functions into the pure-state decomposition. Hence, these states are more quantum than the states that are in the convex hull of Gaussian states $\mathcal{G}$.

Next, one may wonder how to differentiate between states that are quantum non-Gaussian and states which are in the convex hull $\mathcal{G}$. Throughout the last decade, many methods have been developed to answer this question. We start by introducing the main idea of Ref. [114] because it is based on the Wigner function. The key idea is that Gaussian distributions have tails, which means that we can take an arbitrary pure Gaussian state $W_0[S^{-1}(\vec{x} - \vec{\alpha})]$, and evaluate the Wigner function at to origin of phase space:

$$W_0(S^{-1}\vec{\alpha}) = \frac{e^{-\frac{1}{2}\|S^{-1}\vec{\alpha}\|^2}}{(2\pi)^2}. \quad (134)$$

Clearly, when $\|S^{-1}\vec{\alpha}\|^2 \to \infty$, we do find that $W_0(S^1\vec{\alpha}) \to 0$. This limit essentially corresponds to a system with infinite energy. It is thus natural to try to bound the value of the Wigner function in the origin by a function that depends on the energy $\bar{N}$ of the state. For an arbitrary pure Gaussian state, we find that

$$\bar{N} = \sum_{j=1}^{m} \mathrm{tr}[\hat{\rho} \hat{a}^{\dagger}(\vec{e}_j) \hat{a}(\vec{e}_j)] = \frac{1}{4} (\mathrm{tr}[S^T S - 1] + \|\vec{\alpha}\|^2). \quad (135)$$

Using the properties of the operator norm, we write $\|S^{-1}\vec{\alpha}\|^2 = \|S^{-1}\mathcal{O}\vec{\alpha}\|^2$, where $\mathcal{O}$ is a symplectic orthogonal transformation. Furthermore, we note that $\|\vec{\alpha}\|^2 = \|\mathcal{O}\vec{\alpha}\|^2$. It is then useful to explicitly write the coordinate representation of the vector $\mathcal{O}\vec{\alpha} = (\alpha_1^{(x)}, \alpha_1^{(p)}, \dots, \alpha_m^{(x)}, \alpha_m^{(p)})^T$, such that

$$\bar{N} = \sum_{j=1}^{m} \frac{1}{4} \left( s_j + \frac{1}{s_j} + (\alpha_j^{(x)})^2 + (\alpha_j^{(p)})^2 - 2 \right) \quad (136)$$

$$= \sum_{j=1}^{m} \bar{n}_j. \quad (137)$$

At the same time, we expand

$$\|S^{-1}\vec{\alpha}\|^2 = \sum_{j=1}^{m} s_j (\alpha_j^{(x)})^2 + \frac{(\alpha_j^{(p)})^2}{s_j}, \quad (138)$$

and with a little algebra we can show that

$$\frac{1}{2} \|S^{-1}\vec{\alpha}\|^2 \le \sum_{j=1}^{m} 4\bar{n}_j (2\bar{n}_j + 1) \le 4\bar{N}(2\bar{N} + 1), \quad (139)$$

such that

$$W_0(S^{-1}\vec{\alpha}) \ge \frac{1}{(2\pi)^m} e^{-4\bar{N}(2\bar{N}+1)}. \quad (140)$$

This means that the value of the Wigner function of a pure Gaussian state in the origin of phase space is bounded from below by a function of the average number of particles $\bar{N}$. However, it is not obvious that we can extend this bound to arbitrary mixtures of Gaussian states. Let us assume that $\hat{\rho} \in \mathcal{G}$, then the Wigner function in the origin is given by

$$W(\vec{0}) = \int d\gamma p(\gamma) W_0(S_{\gamma}^{-1}\vec{\alpha}_{\gamma}), \quad (141)$$

where we saw that $W_0(S_{\gamma}^{-1}\vec{\alpha}_{\gamma})$ is the value of a pure Gaussian state's Wigner function in the origin and $\gamma$ is some arbitrary label for the Gaussian states in the mixture. Therefore we can bound the states in the convex combination

$$W(\vec{0}) \ge \frac{1}{(2\pi)^m} \int_0^\infty d\bar{N}_\gamma \tilde{p}(\bar{N}_\gamma) e^{-4\bar{N}_\gamma(2\bar{N}_\gamma + 1)}, \quad (142)$$

where we introduce a probability distribution $\tilde{p}$ on the average particle numbers of the pure Gaussian states in the mixture. The overall average number of particles in the state $\hat{\rho}$ is then given by $\bar{N} = \int_0^\infty d\bar{N}_\gamma \tilde{p}(\bar{N}_\gamma)\bar{N}_\gamma$. The final element that we require is the fact that $\exp[-4\bar{N}_\gamma(2\bar{N}_\gamma + 1)]$ is a convex function, such that we can apply Jensen's inequality to find that

$$\hat{\rho} \in \mathcal{G} \implies W(\vec{0}) \ge \frac{1}{(2\pi)^m} e^{-4\bar{N}(2\bar{N}+1)}. \quad (143)$$

This means that we can simply use the total energy of the state to construct a witness for quantum non-Gaussianity. This clearly shows that there are quantum non-Gaussian states with positive Wigner functions. An explicit example can be constructed by tuning the $\gamma$ in the state $[(1 - \gamma)|0\rangle\langle 0| + \gamma|1_\bar{f}\rangle\langle 1_\bar{f}|]$ to $1/2 > \gamma > 1/2 - e^{-4\gamma(2\gamma+1)}$ (where we use that $\gamma$ is also the average particle number in this particular state).

Nevertheless, there are many quantum non-Gaussian states that do not violate inequality (143). After all, why would the origin of phase space be the most interesting point? A first solution is provided in Ref. [114], where it is
argued that

$$ \hat{\rho} \in \mathcal{G} \implies W_{\Gamma}(\vec{0}) \ge \frac{1}{(2\pi)^m} e^{-4\overline{N}_{\Gamma}(2\overline{N}_{\Gamma}+1)}, \quad (144) $$

for all Gaussian channels $\Gamma$ that act on the Wigner function of $\hat{\rho}$ as $W(\vec{x}) \stackrel{\Gamma}{\mapsto} W_{\Gamma}(\vec{x})$. Recall that the action of a Gaussian channel was defined in Eq. (102). The quantity $\overline{N}_{\Gamma}$ then denotes the average number of particles in the state $\Gamma(\hat{\rho})$. Further generalizations of this scheme have been worked out in Ref. [115]. Moreover, Ref. [116] has considered combinations of the value of the Wigner function in several points to reach better witnesses for quantum non-Gaussianity. Further progress has been made by identifying observables that are more easily measurable with typical CV techniques [117]. Others have considered other phase-space representations of the state to identify witnesses of quantum non-Gaussianity [118,119].

Quantum non-Gaussianity has been investigated with a wide range of tools. In this Tutorial we have limited ourselves to phase-space methods in the spirit of the CV approach. However, there is also a significant body of work on quantum non-Gaussianity using techniques that are more typical in DV quantum optics. The earliest works on the subject used photon statistics to distinguish quantum non-Gaussian states from convex mixtures of Gaussian states [89]. This research line has been continued in recent years to uncover new aspects of quantum non-Gaussian states, such as the “non-Gaussian depth” [120] and techniques to differentiate different types of multiphoton states [121]. Ultimately, these photon-counting techniques were extended to develop a whole hierarchy of quantum non-Gaussian states [37]. These ideas have been further formalized and generalized through the notion of the “stellar rank” of a quantum state.

## D. Stellar rank

An interesting starting point to introduce the stellar representation is the method [103] to prove Hudson’s theorem. We recall the definition

$$ F_{\Psi}^{*}(\vec{\alpha}) := \langle \vec{\alpha} | \Psi \rangle e^{\frac{1}{8} \|\vec{\alpha}\|^2}, $$

of what we henceforth refer to as the stellar function. To avoid technical complications, let us now restrict ourselves to single-mode systems such that the optical phase space is $\mathbb{R}^2$. Note that, in a single-mode system, there is only one creation operator $\hat{a}^\dagger$ with associated annihilation operators $\hat{a}$. We can then follow Ref. [38] to introduce the stellar representation of single-mode quantum states. Note that some similar ideas are also present in other works [122].

First, we develop the stellar representation for pure states, which will then be used to generalize the framework to mixed states in Eq. (148). We use the definition of the

displacement operator to show that

$$ |\Psi\rangle = F_{\Psi}^{*}(\hat{a}^{\dagger}) |0\rangle, \quad (145) $$

which immediately implies that the stellar representation is unique. In other words, if $F_{\Psi}^* = F_{\Phi}^*$ it follows that $|\Psi\rangle = |\Phi\rangle$ up to a phase. In our proof of Hudson’s theorem, we have already highlighted that $F_{\Psi}^*$ satisfies the property (116), which means that it is an entire function with growth order $r = 2$. Because in this single-mode setting $F_{\Psi}^*$ can be interpreted as a function of a single complex variable, the Hadamard-Weierstrass theorem implies that $F_{\Psi}^*$ can be fully represented by its zeros (one can consider this as a generalization of the fundamental theorem of calculus). This thus implies that a single-mode state is completely determined by the zeros of the $F_{\Psi}^*$, and thus by the zeros of the $Q$ function in Eq. (115).

It is thus natural to use these zeros in order to classify pure single-mode quantum states and thus the stellar rank is introduced. Ultimately, the stellar rank is simply given by the number of zeros of $F_{\Psi}^*$ or alternatively the number of zeros of the $Q$ function. Because in practice zeros may coincide, one should also consider the multiplicity of the zeros. We thus define the stellar rank $r^*(\Psi)$ of $|\Psi\rangle$ as the number of zeros of $F_{\Psi}^*$: $\mathbb{C} \mapsto \mathbb{C}$ counted with multiplicity. Alternatively one may count the zeros of the $Q$ function with multiplicity and divide by two.

The fact that a state is fully characterized by its stellar representation $F_{\Psi}^*$ can be made more explicit by considering the roots $\{\vec{\alpha}_1, \dots, \vec{\alpha}_{r^*(\Psi)}\}$ of the $Q$ function, which represents $|\Psi\rangle$ [note that we use (44) to interchange between phase-space representation and complex Hilbert-space representation]. The single-mode state $|\Psi\rangle$ can then be used to express

$$ |\Psi\rangle = \frac{1}{N} \prod_{j=1}^{r^*(\Psi)} \hat{D}^{\dagger}(\vec{\alpha}_j) \hat{a}^{\dagger} \hat{D}(\vec{\alpha}_j) |\Psi_G\rangle, \quad (146) $$

where $|\Psi_G\rangle$ is a pure Gaussian state and $N$ a normalization constant. We can then use the stellar rank to induce some further structure in the set of states by defining

$$ \mathcal{R}_N := \{|\Psi\rangle | r^*(\Psi) = N\}. \quad (147) $$

that groups all states of stellar rank $N$. Note that the Hadamard-Weierstrass theorem also considers functions with an infinite amount of zeros and the case $N = \infty$ is thus mathematically well defined. It turns out that this case is not just a pathological limit. An evaluation of the $Q$ function shows that Gottesman-Kitaev-Preskill states and Schrödinger cat states inhabit the set $\mathcal{R}_\infty$.

Clearly, all that was introduced so far only works for pure states. We can naturally extend this result via a convex
roof construction, by defining

$$r^*(\hat{\rho}) := \inf_{\{p(\gamma), |\Psi(\gamma)|\}} \sup_{\gamma} \{r^*[\Psi(\gamma)]\}, \quad (148)$$

where the infimum is considered of all probability distributions on the set of pure states that lead to $\hat{\rho} = \int d\gamma p(\gamma) |\Psi(\gamma)\rangle\langle\Psi(\gamma)|$. In words, there are many ways to decompose the state $\hat{\rho}$ in pure states and we consider all of them. For each decomposition, we define the stellar rank as the highest rank of the states in the decomposition. Then we minimize these values over all possible decompositions to arrive at the stellar rank of $\hat{\rho}$.

Convex roof constructions are commonly used to treat mixed states as they are easy and natural to formally define. However, they are often much harder to calculate in practice. This is where the stellar representation unveils its most remarkable property: stellar robustness. To formalize this idea, Ref. [38] introduced the robustness as the trace distance between the state and the nearest possible state of lower stellar rank.

$$R^*(\Psi) := \inf_{r^*(\hat{\rho})<r^*(\Psi)} \frac{1}{2} \text{tr} \sqrt{|\Psi\rangle\langle\Psi|-\hat{\rho}|^2}. \quad (149)$$

And it can be shown that

$$R^*(\Psi) = \sqrt{1 - \sup_{r^*(\hat{\rho})<r^*(\Psi)} \langle\Psi|\hat{\rho}|\Psi\rangle}, \quad (150)$$

where $\langle\Psi|\hat{\rho}|\Psi\rangle$ is the fidelity of $\hat{\rho}$ with target state $|\Psi\rangle$. Remarkably, it can be shown that $R^*(\Psi) > 0$ when $r^*(\Psi) < \infty$ [38]. This means that any state that is sufficiently close to a pure state of $r^*(\Psi)$ is also of rank $r^*(\Psi)$ or higher.

The idea of stellar robustness can be generalized [123] by introducing $k$ robustness. For any $k < r^*(\Psi)$, we define

$$R_k^*(\Psi) := \inf_{r^*(\Phi) \le k} \sqrt{1 - |\langle\Phi|\Psi\rangle|^2}, \quad (151)$$

and show subsequently that

$$R_k^*(\Psi) = \sqrt{1 - \sup_{r^*(\hat{\rho}) \le k} \langle\Psi|\hat{\rho}|\Psi\rangle}. \quad (152)$$

The $k$ robustness can thus be interpreted as the nearest distance from a state $|\Psi\rangle$ at which we can find any state of stellar rank $k$, provided $k < r^*(\Psi)$. Beyond showing that $R_k^*(\Psi)$ is nonzero when $|\Psi\rangle$ is of finite stellar rank, Ref. [123] also provides an explicit method to calculate $R_k^*(\Psi)$. Thus, for whichever state $\hat{\rho}$ is available in an experiment one can attempt to find a pure target state $|\Psi\rangle$ for which $\langle\Psi|\hat{\rho}|\Psi\rangle > 1 - R_k^*(\Psi)^2$ to prove that $\hat{\rho}$ is at least of stellar rank $k$.

We note that for pure states $r^*(\Psi) = 0$ implies that the state is Gaussian (this is essentially what is proven in Hudson's theorem). From the definition (148) we can then deduce that

$$r^*(\hat{\rho}) = 0 \iff \hat{\rho} \in \mathcal{G}, \quad (153)$$

where $\mathcal{G}$ denotes, again, the convex hull of Gaussian states. On the other hand, states for which $r^*(\hat{\rho}) > 0$ cannot be written as a mixture of Gaussian states and are thus quantum non-Gaussian.

This idea can be extended by using the stellar $k$ robustness, Eq. (152), as a witness of quantum non-Gaussianity. When we want to check whether $\hat{\rho}$ is quantum non-Gaussian, it suffices to find a pure target state $|\Psi\rangle$ such that $\hat{\rho}$ is closer to $|\Psi\rangle$ than the 1 robustness $R_k^*(\Psi)$. More formally written, whenever a pure state $|\Psi\rangle$ exists with the following property:

$$\langle\Psi|\hat{\rho}|\Psi\rangle > 1 - R_1^*(\Psi)^2 \implies \hat{\rho} \notin \mathcal{G}, \quad (154)$$

and $\hat{\rho}$ is quantum non-Gaussian. This may seem like a complicated challenge, but for a single-photon state $|\Psi\rangle = |1\rangle$ we find that $1 - R_1^*(\Psi)^2 \approx 0.478$ [123]. This means that any state that has a fidelity of more than 0.478 with respect to a Fock state is quantum non-Gaussian. This idea can be extended to higher stellar ranks: whenever we find a target state $|\Psi\rangle$ such that $\langle\Psi|\hat{\rho}|\Psi\rangle > 1 - R_k^*(\Psi)^2$, the state $\hat{\rho}$ is at least of stellar rank $k$. Note that the fidelity of an experimentally generated state $\hat{\rho}$ with any pure target state $|\Psi\rangle$ can be calculated from double homodyne measurements on $\hat{\rho}$ [123]. There is no need to experimentally create the pure state $|\Psi\rangle$, the latter is just theoretical input needed to analyze the data.

Obviously, the stellar rank imposes a lot of additional structure on the state space. It rigorously orders all states that can be achieved by combining a finite number of creation operators and Gaussian transformations. The creation operator serves as a tool to increase the stellar rank by one and the stellar rank actually corresponds to the minimal number of times the creation operator must be applied to obtain the state, together with Gaussian operations. The stellar rank remains unchanged under Gaussian unitary transformations, which makes sense for a measure of the non-Gaussian character of the state, and thus it falls within the set of intrinsic properties of a state as discussed in Sec. C. Furthermore, the class of states with infinite stellar rank can be understood as the set that contains the most exotic states. However, it is lonely at the top as it can be shown that $R_\infty^*(\Psi) = 0$ for states of infinite rank. This means that we can find states of finite stellar rank arbitrarily close to a state of infinite stellar rank. As stressed in Ref. [38], this implies that finite-rank states are dense in the full Fock space and any state of infinite rank can be arbitrarily well approximated by finite-rank states. Whereas a finite
stellar rank k of any experimental state can be certified by achieving a sufficiently high fidelity to a target state $|\Psi\rangle$ to fall within the range given by its k robustness $R_k^*(\Psi)$, a similar procedure is impossible for infinite stellar ranks. This means, in practice, that genuinely infinite-rank states are impossible to certify in experiments because one never achieves perfect fidelity. Nevertheless, different states of infinite rank may differ significantly in the values $R_k^*(\Psi)$ for $k < \infty$.

Many of the results on stellar rank rely on the Hadamard-Weierstrass theorem that allows one to uniquely factorize $F_\Psi^*(\vec{\alpha})$ as a Gaussian and a polynomial, where the roots of the polynomial are the roots of $F_\Psi^*$ and thus also the roots of the Q function. Sadly, this theorem cannot be straightforwardly generalized to a multimode setting, which is known in mathematics as Cousin's second problem [124]. Notable progress was made in Ref. [125] where one studies multimode stellar functions, which are polynomials and it was shown that there is no straightforward generalization of Eq. (146).

## E. Wigner negativity

Hudson's theorem shows us that all pure non-Gaussian states have nonpositive Wigner functions, which sets them apart from normal probability distributions of phase space. For mixed states, this no longer holds and thus we spent the previous two sections developing methods to characterize the non-Gaussian features of these states. Whether it is through quantum non-Gaussianity or the more refined stellar rank, these methods focus on characterizing the non-Gaussian resources that are required to generate a certain state. In this subsection, we change the perspective and focus rather on negative values of the Wigner function ("Wigner negativity" in short) as a resource of interest.

Wigner negativity has the advantage of being a clear quantum feature, it reflects that different quadratures in the same mode cannot be jointly measured and thus goes hand in hand with the principle of complementarity. More formally, it has even been connected to the principle of quantum contextuality [64]. Indeed, in Sec. VI we elaborate on the fact that Wigner negativity is a necessary resource for reaching a quantum advantage, i.e., performing a task that cannot be efficiently simulated by a classical computer. However, the idea of using Wigner negativity as a signature of nonclassicality was already around before it was connected to a quantum computational advantage. An important step to formalize this idea was the introduction of a measure for Wigner negativity [90], which lies at the basis of recent resource theories of Wigner negativity [100,126].

A priori, there are several natural measures than can be used for Wigner negativity. It is therefore useful to consider some desirable properties that are required for a measure of Wigner negativity. First of all, we want the

measure to be zero if and only if the Wigner function is positive. It seems natural to demand that, furthermore, Wigner negativity remains unchanged under Gaussian unitary transformations. It is then tempting to simply consider the absolute value of the lowest possible value of the Wigner function, but this would have some unnatural outcomes. It would mean that a single-photon state would have more Wigner negativity than a two-photon state. We thus need to look for a different measure.

The starting point of Ref. [90] is that the normalization of the Wigner function implies that

$$ \int_{\mathbb{R}^{2m}} d\vec{x} |W(\vec{x})| \ge 1, \qquad (155) $$

and that the inequality is strict whenever there is Wigner negativity. Furthermore, Liouville's theorem implies that integrals over phase space are unchanged by Gaussian transformations. A first possible way of measuring Wigner negativity is through the negativity volume

$$ \mathcal{N}(\hat{\rho}) := \int_{\mathbb{R}^{2m}} d\vec{x} |W(\vec{x})| - 1. \qquad (156) $$

This measure has the major advantage of being convex due to the triangle inequality, which means that for $\hat{\rho} = \int d\gamma p(\gamma) \hat{\rho}(\gamma)$ we find that

$$ \mathcal{N}(\hat{\rho}) \le \int d\gamma p(\gamma) \mathcal{N}[\hat{\rho}(\gamma)]. \qquad (157) $$

However, this measure is not additive, i.e., $\mathcal{N}(\hat{\rho}_1 \otimes \hat{\rho}_2) \ne \mathcal{N}(\hat{\rho}_1) + \mathcal{N}(\hat{\rho}_2)$. To circumvent this shortcoming, another measure for Wigner negativity has been introduced [100, 126,127]:

$$ \mathfrak{N}(\hat{\rho}) := \log \int_{\mathbb{R}^{2m}} d\vec{x} |W(\vec{x})|. \qquad (158) $$

Clearly, $\mathfrak{N}(\hat{\rho}_1 \otimes \hat{\rho}_2) = \mathfrak{N}(\hat{\rho}_1) + \mathfrak{N}(\hat{\rho}_2)$ making this measure additive. However, the introduction of the logarithm destroys the convexity of the measure. Note that the two measures are closely related by $\mathfrak{N}(\hat{\rho}) = \log[\mathcal{N}(\hat{\rho}) + 1]$. Thus when $\mathcal{N}(\hat{\rho}_1) > \mathcal{N}(\hat{\rho}_2)$, we also find that $\mathfrak{N}(\hat{\rho}_1) > \mathfrak{N}(\hat{\rho}_2)$.

The single-mode examples that are considered in Ref. [90] lead to some interesting observations. First of all, they show that for Fock states Wigner negativity increases with the photon number. Furthermore, they show that for Schrödinger cat states the integral is bounded from above by a value smaller than the Wigner negativity of a two-photon state. Even though Fock states of increasing stellar rank have increasing Wigner negativity, there is no clear relation between stellar rank and Wigner negativity for more general classes of states. For example, Schrödinger
cat states are of infinite stellar rank, suggesting that they are in this regard the most exotic states, but they manifest only a limited amount of Wigner negativity.

As a case study, let us briefly concentrate on the Wigner negativity of Fock states, Eq. (117). One can now evaluate the Wigner negativity of such states to find that

$$ \mathcal{N}(|1\rangle) \approx 0.42612 \quad \text{and} \quad \mathfrak{N}(|1\rangle) \approx 0.354959, \quad (159) $$

$$ \mathcal{N}(|2\rangle) \approx 0.72899 \quad \text{and} \quad \mathfrak{N}(|2\rangle) \approx 0.547537, \quad (160) $$

$$ \mathcal{N}(|3\rangle) \approx 0.97667 \quad \text{and} \quad \mathfrak{N}(|3\rangle) \approx 0.681415, \quad (161) $$

which shows that the negativity does not simply increase linearly with the number of photons even for the additive measure $\mathfrak{N}$. However, let us now look at a multimode n-photon state where each photon occupies a different mode, i.e., a Fock state generated by creation operators in $\vec{f}_1, \dots, \vec{f}_n$ with $\text{span}\{\vec{f}_j, \Omega\vec{f}_j\} \neq \text{span}\{\vec{f}_k, \Omega\vec{f}_k\}$ for all $j \neq k$,

$$ \hat{a}^{\dagger}(\vec{f}_1) \dots \hat{a}^{\dagger}(\vec{f}_n) |0\rangle = |1_{f_1}\rangle \otimes \dots \otimes |1_{f_n}\rangle. \quad (162) $$

The Wigner function for this state can be shown to be [showing this based on (38) is again a good exercise]

$$ W_{1_{f_1}, \dots, 1_{f_n}}(\vec{x}_{f_1} \oplus \dots \oplus \vec{x}_{f_n}) = \prod_{k=1}^{n} W_{1_{f_k}}(\vec{x}_{f_k}). \quad (163) $$

Either by explicitly using the expression of the Wigner function, or by using the additivity property, we find that

$$ \mathfrak{N}(|1_{f_1}\rangle \otimes \dots \otimes |1_{f_n}\rangle) = n\mathfrak{N}(|1\rangle). \quad (164) $$

Numerically, we can show that $n\mathfrak{N}(|1\rangle) > \mathfrak{N}(|n\rangle)$ and thus we can generally conclude that $n$ photons in different modes hold more Wigner negativity than $n$ photons in the same mode.

Let us now concentrate on the case where $n=2$. We showed that two photons in different modes are more Wigner negative than two photons in the same mode, and now we combine this finding with the idea that Wigner negativity remains unchanged under Gaussian transformations. A particularly simple Gaussian transformation is a balanced beam splitter, which ultimately just implements a change in mode basis that we describe by an orthonormal transformation $O_{BS}$. When we mix two photons, prepared in orthogonal modes $f_1$ and $f_2$ by such a balanced beam splitter, we see the Hong-Ou-Mandel effect in action (more details can be found in Ref. [56] where a similar notation is used):

$$ |1_{f_1}\rangle \otimes |1_{f_2}\rangle \xrightarrow{O_{BS}} \frac{1}{\sqrt{2}} (|2_{g_1}\rangle - |2_{g_2}\rangle) := |\text{HOM}\rangle, \quad (165) $$

where $f_1, f_2$ and $g_1, g_2$ are the input and output modes of the beam splitter, respectively. The Hong-Ou-Mandel output state $|\text{HOM}\rangle$ is thus a superposition of two photons in

mode $g_1$ and two photons in mode $g_2$. We can now use the simple fact that Wigner negativity is unchanged under Gaussian unitary transformations to show that

$$ \mathfrak{N}(|\text{HOM}\rangle) = \mathfrak{N}(|1_{f_1}\rangle \otimes |1_{f_2}\rangle) = 2\mathfrak{N}(|1\rangle) > \mathfrak{N}(|2\rangle), \quad (166) $$

this then also implies that $\mathcal{N}(|\text{HOM}\rangle) > \mathcal{N}(|2\rangle)$. At first sight, this is somewhat of a peculiar finding: by taking a superposition of two states with the same Wigner negativity one finds a state with a higher Wigner negativity.

An explicit look at the Wigner function of the Hong-Ou-Mandel state $|\text{HOM}\rangle$ provides some insight. We find that this Wigner function can be written as (yet again a good exercise to show this explicitly)

$$ W_{\text{HOM}}(\vec{x}_{g_1} \oplus \vec{x}_{g_2}) = \frac{1}{2}[W_{2g_1}(\vec{x}_{g_1}) + W_{2g_2}(\vec{x}_{g_2})] + W_{\text{int}}(\vec{x}_{g_1} \oplus \vec{x}_{g_2}), \quad (167) $$

where $W_{\text{int}}$ is the contribution to the Wigner function that contains all the interference terms that are induced by the superposition. We can calculate that

$$ \mathcal{N}\left(\frac{1}{2}[W_{2g_1}(\vec{x}_{g_1}) + W_{2g_2}(\vec{x}_{g_2})]\right) \le \mathcal{N}(|2\rangle), \quad (168) $$

and thus, by additionally applying the triangle inequality, we can understand that the additional negativity in the Hong-Ou-Mandel state is due to the term $W_{\text{int}}$.

In the Hong-Ou-Mandel effect, it is common to talk about interference between particles, but in a more general CV language this interference will be equivalent to some form of entanglement, which is exactly described by the Wigner-function contribution $W_{\text{int}}$. In other words, the superposition between $|2_{g_1}\rangle$ and $|2_{g_2}\rangle$ has more Wigner negativity than each of its two constituents because it creates entanglement between the modes $g_1$ and $g_2$. This is a first indication that there is a connection between quantum correlations and non-Gaussian features of the Wigner function. We explore this connection in further detail in Sec. V.

Even though Wigner negativity is an important non-Gaussian feature, it is often hard to witness [128]. The most common experimental technique is homodyne tomography [129] to fully reconstruct the quantum state. These methods come with the inconvenience that it is hard to set good error bars. Techniques to circumvent the need for a full tomography have been developed based on homodyne [130] and double-homodyne (or heterodyne) measurements [123,128]. These methods come with the advantage of permitting to put a degree of confidence on the proclaimed Wigner negativity.
#### IV. CREATING NON-GAUSSIAN STATES

In Sec. III we have discussed the many ways of characterizing non-Gaussian quantum states and their properties. In this section, we explore the different theoretical frameworks for creating these states. An overview of some important experimental advances to put these theoretical techniques into practice is left for Sec. VII.

Gaussian quantum states can in some sense be understood as naturally occurring states. The foundational work of Planck that lies at the basis of all of quantum mechanics provides a first description of the thermal states of light that describe black-body radiation. In a more modern language, we refer to this as the thermal states of an ensemble of quantum harmonic oscillators or a free bosonic field. It has long been understood that these states are Gaussian [2–4]. Creating this kind of Gaussian states of light is thus literally as simple as switching on a light bulb.

When we turn towards more sophisticated light sources such as lasers, we can encounter coherent light that is described by coherent states [5,6]. Generating squeezed light becomes much harder and typically requires nonlinear optics [131]. Nevertheless, pumping a nonlinear crystal with a coherent pump generally suffices to deterministically create a squeezed state [132]. Recall from the end of Sec. D that from a theoretical point of view all these pure Gaussian states can be created by applying Gaussian unitary transformation to the vacuum state.

From an experimental point of view, the creation of non-Gaussian states is much harder than the creation of their Gaussian counterparts. Nevertheless, we start by introducing an ideal theoretical approach that is not too different from Gaussian states. In essence, it suffices to apply a non-Gaussian unitary operation to the state to create a non-Gaussian state. In Sec. A we dig deeper into the desired structure of such non-Gaussian unitary transformations that would in principle allow for the deterministic generation of non-Gaussian quantum states. In experiments (in particular, those in optics) such non-Gaussian unitary transformations are hard to come by, which is why one very often uses different preparation schemes. In Sec. B, we provide a general introduction into the conditional preparation of non-Gaussian quantum states, where one measures part of the system and conditions on a certain measurement outcome. This process projects the remainder of the system into a new non-Gaussian state.

##### A. Deterministic methods

To introduce some further structure in the sets of Gaussian and non-Gaussian unitary transformations, it is useful to take a quantum computation approach that is inspired by Ref. [26,133]. The central idea of this work is that Gaussian unitary transformations are always generated by “Hamiltonians” that are at most quadratic in the quadrature operators (or equivalently in the creation and annihilation

operators). Let us denote that as

$$ \hat{U}_G = \exp\{i\mathcal{P}_2(\hat{q})\} \qquad (169) $$

where the polynomials $\mathcal{P}_2(\hat{q})$ are generated by combining terms of the types $1$, $\hat{q}(\vec{f})$, and $\hat{q}(\vec{f}_1)\hat{q}(\vec{f}_2)$. A remarkable property of these three types of observables is that they are closed under the action of a commutator. Indeed, using the canonical commutation relation (29) and the general properties of commutators, we can show that

$$ [\hat{q}(\vec{f}_1), \hat{q}(\vec{f}_2)] \sim 1, \qquad (170) $$

$$ [\hat{q}(\vec{f}_1), \hat{q}(\vec{f}_2)\hat{q}(\vec{f}_3)] \sim \hat{q}(\vec{f}'), \qquad (171) $$

$$ [\hat{q}(\vec{f}_1)\hat{q}(\vec{f}_2), \hat{q}(\vec{f}_3)\hat{q}(\vec{f}_4)] \sim \sum \hat{q}(\vec{f}_1')\hat{q}(\vec{f}_2'). \qquad (172) $$

Thus, we can use the Baker-Campbell-Hausdorff formula to show that the combination of two Gaussian unitaries $\hat{U}_G \hat{U}_G'$ is again a Gaussian unitary.

This notion lies at the basis of universal gate sets in the CV approach. Using typical techniques from Lie groups, we can look for a minimal set of Gaussian unitaries than can be combined to generate all possible Gaussian unitary transformations. Generally, such a set is clearly not unique, but there are some natural choices. For example, we previously saw that a Gaussian unitary transformation is a combination of displacement operations and symplectic transformations. Furthermore, the Bloch-Messiah decomposition (95) shows us that any symplectic transformation can be decomposed into a combination of multimode interferometers and single-mode squeezing. In turn, interferometers can be decomposed as a combination of beam splitters and phase shifters [134]. Indeed, we can choose the set of Gaussian gates to be

$$ \hat{U}_D(\vec{\lambda}) := \hat{D}(\vec{\lambda}), \qquad (173) $$

$$ \hat{U}_S(\vec{\lambda}) := \exp i[\hat{q}(\vec{\lambda})\hat{q}(\Omega\vec{\lambda}) + \hat{q}(\Omega\vec{\lambda})\hat{q}(\vec{\lambda})], \qquad (174) $$

$$ \hat{U}_P(\vec{\lambda}) := \exp i[\hat{q}(\vec{\lambda})^2 + \hat{q}(\Omega\vec{\lambda})^2], \qquad (175) $$

$$ \hat{U}_{BS}(\vec{\lambda}_1, \vec{\lambda}_2) := \exp i[\hat{q}(\vec{\lambda}_1)\hat{q}(\vec{\lambda}_2) + \hat{q}(\Omega\vec{\lambda}_1)\hat{q}(\Omega\vec{\lambda}_2)], \qquad (176) $$

where we note that $\vec{\lambda}, \vec{\lambda}_1, \vec{\lambda}_2 \in \mathbb{R}^{2m}$ are not normalized and $\vec{\lambda}_1 \perp \vec{\lambda}_2$. These unitary operators describe a displacement, a squeezer, a phase shifter, and a beam splitter, respectively. We note that all these operations act on a single mode, except for the beam splitter, which connects a pair of modes. These transformations are referred to as a Gaussian gate set; when we can implement all these gates in all the modes of some mode basis, we can generate any multimode Gaussian transformation, and thus any Gaussian state.

To generate non-Gaussian unitary transformations and thus non-Gaussian states, we need to add more unitary
gates to the Gaussian gate set. The relevant question is thus how many gates one should add and which gates are the best choices. The answer to the first question is surprising: one needs to add just one single gate [27]. The argument is simple, when we consider the operators $\hat{q}(\vec{\lambda})^3$, we find that

$$[\hat{q}(\vec{\lambda}_1)^3, \hat{q}(\vec{\lambda}_2)^n] = \hat{q}(\vec{\lambda}_1)^2[\hat{q}(\vec{\lambda}_1), \hat{q}(\vec{\lambda}_2)^n] \\ + [\hat{q}(\vec{\lambda}_1), \hat{q}(\vec{\lambda}_2)^n]\hat{q}(\vec{\lambda}_1)^2 \\ + \hat{q}(\vec{\lambda}_1)[\hat{q}(\vec{\lambda}_1), \hat{q}(\vec{\lambda}_2)^n]\hat{q}(\vec{\lambda}_1). \quad (177)$$

With the canonical commutation relations we can show that $[\hat{q}(\vec{\lambda}_1), \hat{q}(\vec{\lambda}_2)^n] \sim \hat{q}(\vec{\lambda}_2)^{n-1}$, which can be inserted into Eq. (177) to obtain

$$[\hat{q}(\vec{\lambda}_1)^3, \hat{q}(\vec{\lambda}_2)^n] \sim \hat{q}(\vec{\lambda}_1)^2 \hat{q}(\vec{\lambda}_2)^{n-1} + \hat{q}(\vec{\lambda}_2)^{n-1} \hat{q}(\vec{\lambda}_1)^2 \\ + \hat{q}(\vec{\lambda}_1) \hat{q}(\vec{\lambda}_2)^{n-1} \hat{q}(\vec{\lambda}_1). \quad (178)$$

This calculation thus shows that commutation with the operators $\hat{q}(\vec{\lambda})^3$ increases the order of the quadrature operators. Thus, with the quadratic Hamiltonians to generate all operations that conserve the order of polynomials of quadrature operators and $\hat{q}(\vec{\lambda})^3$ to increase the order of the polynomial by one, we can ultimately generate the full algebra of observables. On the level of unitary gates, this implies that a full universal gate set is given by

$$U = \{\hat{U}_D(\vec{\lambda}), \hat{U}_S(\vec{\lambda}), \hat{U}_P(\vec{\lambda}), \hat{U}_{BS}(\vec{\lambda}_1, \vec{\lambda}_2), \hat{U}_C(\vec{\lambda})\}, \quad (179)$$

$$\text{with } \hat{U}_C := \exp i\hat{q}(\vec{\lambda})^3. \quad (180)$$

In other words, combining sufficiently many of these gates allows us to build any arbitrary unitary transformation generated by a Hamiltonian, which is polynomial in the quadrature operators.

The non-Gaussian gate $\hat{U}_C$ is known as the cubic phase gate. The argument above shows that any experiment that can implement Gaussian transformations and a cubic phase gate can in principle generate any arbitrary non-Gaussian state. Even though many protocols have been proposed to experimentally realize a cubic phase gate [28,135–138], any convincing implementations have yet to be demonstrated. One of the key problems is that experimental imperfections and finite squeezing are detrimental for the most commonly proposed methods [30].

In principle, there is no particular reason to limit our attention to cubic phase gates. Already in the very first work on the subject it is argued that essentially any Hamiltonian of a higher than quadratic order can be used as a generator [26]. Thus, optical processes that perform photon triplet generation can also be used as a non-Gaussian gate, which can even be converted into the cubic phase gate [139]. This requires well-controlled high $\chi^{(3)}$ nonlinearities, which are generally only achieved by using exotic

nonlinear crystals or well-controlled individual atoms. Handling such setups with a sufficient degree of control to actually implement a quantum gate is extremely challenging.

Other CV systems are more appropriate for the implementation of non-Gaussian unitary transformations. In particular, the systems used in circuit QED have such non-Gaussian contributions in their Hamiltonians [140,141], which suggests that they may be more capable of deterministically generating non-Gaussian states than their optical counterparts. Still, the characterization, detection, and control of such states is expected to be challenging. Recently, some important progress was made by demonstrating triplet generation in these systems [142].

## B. Conditional methods

The experimental difficulties that are encountered when trying to implement non-Gaussian unitary transformations can be circumvented by abandoning the demand of unitarity. This implies that we no longer consider operations that can be implemented deterministically, but rather resort to what can be broadly referred to as conditional operations. This idea was formalized by Kraus when characterizing the most general ways of manipulating quantum states [143].

In the most general sense, we can implement a conditional operation by taking a set of linear operators on Fock space $\hat{X}_1, \hat{X}_2, \dots$ and acting on the state in the following way:

$$\hat{\rho} \mapsto \frac{\sum_j \hat{X}_j \hat{\rho} \hat{X}_j^\dagger}{\text{tr}[\hat{\rho} \sum_j \hat{X}_j^\dagger \hat{X}_j]}. \quad (181)$$

This formalism is typically implemented by performing some form of generalized measurement on the state $\hat{\rho}$ [144]. When $\hat{\rho}$ is a deterministically generated Gaussian state, the action of a well-chosen set of operators $\hat{X}_1, \hat{X}_2, \dots$ can turn it into a non-Gaussian state. In optics, two of the most well-known examples of this technique are single-photon addition and subtraction. In both cases, there is only a single operator $\hat{X}_1$. For photon addition, we implement $\hat{X}_1 = \hat{a}^\dagger(\vec{f})$, whereas photon subtraction requires the realization of a case where $\hat{X}_1 = \hat{a}(\vec{f})$.

In many physical setups, and, in particular, in optics, the problem is that measurements are destructive and a measurement effectively removes the measured mode from the system. Therefore, it is common to prepare large multimode Gaussian states of which a subset of modes is measured in order to conditionally prepare a non-Gaussian state in the remaining modes. We now introduce a general framework to describe the non-Gaussian Wigner functions that are created accordingly [39].
<figure><img src="image_5.png" /><figcaption>FIG. 3. Sketch representation of the conditional preparation scheme for creating the non-Gaussian states described by Eq. (190). Note that both **f** and **g** can be highly multimode. The Wigner function shown on the right is obtained by a conditional protocol shown in Ref. [39].</figcaption></figure>

1. General framework

First of all, let us consider a general multimode phase space and separate it into two subsystems, i.e., $\mathbb{R}^{2m} = \mathbb{R}^{2l} \oplus \mathbb{R}^{2l'}$, where we perform some generalized measurement on the $l'$ modes and leave the remaining $l$ modes untouched. This introduces a general structure in the points of phase space $\vec{x} \in \mathbb{R}^{2m}$, which can now be written as $\vec{x} = \vec{x}_f \oplus \vec{x}_g$ with $\vec{x}_f \in \mathbb{R}^{2l}$ and $\vec{x}_g \in \mathbb{R}^{2l'}$. The general procedure is schematically outlined in Fig. 3 and we present the details step by step.

Any state $\hat{\rho}$ on this system then comes with a Wigner function $W(\vec{x}) = W(\vec{x}_f \oplus \vec{x}_g)$ that is defined on the global phase space. This state can be reduced to one of the two subsystems by tracing out the other subsystem, which can be described on the level of the Wigner function by the following integrals:

$$W_f(\vec{x}_f) := \int_{\mathbb{R}^{2l'}} d\vec{x}_g W(\vec{x}_f \oplus \vec{x}_g), \quad (182)$$

$$W_g(\vec{x}_g) := \int_{\mathbb{R}^{2l'}} d\vec{x}_f W(\vec{x}_f \oplus \vec{x}_g). \quad (183)$$

When the state is Gaussian and the Wigner function is given by Eq. (85), the structure of the phase space is reflected in the mean field vector $\vec{\xi}$ and in the covariance matrix $V$:

$$\vec{\xi} = \vec{\xi}_f \oplus \vec{\xi}_g, \quad (184)$$

$$V = \begin{pmatrix} V_f & V_{fg} \\ V_{gf} & V_g \end{pmatrix}, \quad (185)$$

with $V_{fg} = V_{gf}^T$. The matrices $V_f$ and $V_g$ describe all the variances and correlations of the modes within $\mathbb{R}^{2l}$ and $\mathbb{R}^{2l'}$, respectively. In addition, the submatrix $V_{gf}$ contains all the correlations between the modes in the different subspaces, which will be important for conditional state preparation. One can show that for such Gaussian states, the reduced states are also Gaussian, for the modes in $\mathbb{R}^{2l}$ given by

$$W_f(\vec{x}_f) = \frac{e^{-\frac{1}{2}(\vec{x}_f - \vec{\xi}_f)^T V_f^{-1} (\vec{x}_f - \vec{\xi}_f)}}{(2\pi)^m \sqrt{\det V_f}}, \quad (186)$$

and analogously for the modes in $\mathbb{R}^{2l'}$. As Gaussian states are the states that are least challenging to produce, they form the starting point of the conditional state preparation scheme.

As a next step, we must implement some form of operation on the modes that correspond to the phase space $\mathbb{R}^{2l'}$. To do so, we consider the action of a general positive operator-valued measure (POVM) element $\hat{A} \ge 0$ that corresponds to a specific measurement outcome. We can then obtain a conditional state via

$$\hat{\rho}_{f|A} := \frac{\mathrm{tr}_g[\hat{A}\hat{\rho}]}{\mathrm{tr}[\hat{A}\hat{\rho}]}. \quad (187)$$

The partial trace $\mathrm{tr}_g[\hat{A}\hat{\rho}]$ runs only over the modes in $\mathbb{R}^{2l'}$ because the other modes are left untouched. The denominator $\mathrm{tr}[\hat{A}\hat{\rho}]$ renormalizes the state and gives the probability of actually obtaining the measurement result that corresponds to $\hat{A}$. In an actual experiment, this operation is implemented by many repeated measurements of the modes in $\mathbb{R}^{2l'}$ and $\hat{A}$ corresponds to a specific detector output of these measurements. The nonmeasured part of the state is only used when the detector indicates this specific output, otherwise it is simply discarded. This conditional selection of the state significantly changes the properties of the state in a way that is strongly influenced by $\hat{A}$.

As we described in Eq. (42), the operator $\hat{A}$ comes with an associated phase-space representation $W_A(\vec{x}_g)$, which can be used to formally describe the phase-space representation of $\hat{\rho}_{f|A}$:

$$W_{f|A}(\vec{x}_f) = \frac{\int_{\mathbb{R}^{2l'}} d\vec{x}_g W_A(\vec{x}_g) W(\vec{x}_f \oplus \vec{x}_g)}{\int_{\mathbb{R}^{2l''}} d\vec{x}_g W_A(\vec{x}_g) W_g(\vec{x}_g)}. \quad (188)$$

There is a more practical way of expressing this Wigner function by exploiting the fact that the initial multimode Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ is positive and therefore describes a well-defined probability distribution on phase space. This implies that the conditional probability distribution

$$W(\vec{x}_g | \vec{x}_f) := \frac{W(\vec{x}_f \oplus \vec{x}_g)}{W_f(\vec{x}_f)}, \quad (189)$$

is also a well-defined probability distribution, which is obtained when we fix one point in phase space $\vec{x}_f \in \mathbb{R}^{2l}$ and look at the probability distribution for the remaining modes in $\mathbb{R}^{2l'}$. We can then use this conditional probability distribution to write $W(\vec{x}_f \oplus \vec{x}_g) = W(\vec{x}_g | \vec{x}_f) W_f(\vec{x}_f)$, which can be inserted in Eq. (188) to find

$$W_{f|A}(\vec{x}_f) = \frac{\langle \hat{A} \rangle_{g|\vec{x}_f} W_f(\vec{x}_f)}{\langle \hat{A} \rangle}, \quad (190)$$
where we define

$$ \langle \hat{A} \rangle := (4\pi)^{\prime} \int_{\mathbb{R}^{2\prime}} d\vec{x}_g W_A(\vec{x}_g) W_g(\vec{x}_g), \quad (191) $$

$$ \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). \quad (192) $$

The quantity $\langle \hat{A} \rangle$ is simply the expectation value of the observable $\hat{A}$ in the state $\hat{\rho}$. $\langle \hat{A} \rangle_{g|\vec{x}_f}$, on the other hand, is the expectation value of the function $W_A(\vec{x}_g)$ where $\vec{x}_g$ is distributed according to the distribution $W(\vec{x}_g | \vec{x}_f)$, which makes $\langle \hat{A} \rangle_{g|\vec{x}_f}$ a function of the selected phase-space point $\vec{x}_f$. However, even though $W(\vec{x}_g | \vec{x}_f)$ is a well-defined probability distribution on phase space, it does not necessarily correspond to a quantum state. Indeed, $W(\vec{x}_g | \vec{x}_f)$ may violate the Heisenberg inequality, which will be of vital importance in Sec. V as it is narrowly connected to quantum steering.

In the specific case where $W(\vec{x}_f \oplus \vec{x}_g)$ is Gaussian, we find that $W(\vec{x}_g | \vec{x}_f)$ is also a Gaussian probability distribution, given by

$$ W(\vec{x}_g | \vec{x}_f) = \frac{\exp\left[-\frac{1}{2}(\vec{x}_g - \vec{\xi}_{g|\vec{x}_f})^T V_{g|\vec{x}_f}^{-1}(\vec{x}_g - \vec{\xi}_{g|\vec{x}_f})\right]}{(2\pi)^{\prime} \sqrt{\det V_{g|\vec{x}_f}}}. \quad (193) $$

Using the notation of Eq. (185), we express its covariance matrix

$$ V_{g|\vec{x}_f} = V_g - V_{gf}V_f^{-1}V_{gf}^T, \quad (194) $$

and mean field vector

$$ \vec{\xi}_{g|\vec{x}_f} = \vec{\xi}_g + V_{gf}V_f^{-1}(\vec{x}_f - \vec{\xi}_f). \quad (195) $$

The covariance matrix $V_{g|\vec{x}_f}$ is known in the mathematics literature [145] as the Schur complement of $V$. The Schur complement has interesting properties, for example, $V$ is a positive matrix if and only if the same holds for the Schur complement $V_{g|\vec{x}_f}$. This immediately implies that the Gaussian probability distribution in Eq. (193) is well defined. Furthermore, the Schur complement also plays an important role in the theory of Gaussian quantum correlations [146]. It should be noted that $V_{g|\vec{x}_f}$ does not actually depend on the chosen value for $\vec{x}_f$. Thus, the conditional expectation value $\langle \hat{A} \rangle_{g|\vec{x}_f}$ depends only on the phase-space point $\vec{x}_f$ through the displacement $\vec{\xi}_{g|\vec{x}_f}$. This is a particular feature of Gaussian states.

Finally, remark that the derivation of Eq. (190) holds true for all initial states with a positive Wigner function. Whenever the initial multimode Wigner function $W(\vec{x}_f \oplus \vec{x}_g)$ is positive, it follows that $W_f(\vec{x}_f)$ is also positive. Furthermore, given that $\langle \hat{A} \rangle$ is the quantum expectation

value of a positive semidefinite operator it clearly also is a positive quantity. Hence, Wigner negativity is entirely contained with $\langle \hat{A} \rangle_{g|\vec{x}_f}$. The fact that $\langle \hat{A} \rangle_{g|\vec{x}_f}$ can take negative values is exactly due to $W(\vec{x}_g | \vec{x}_f)$ not being the Wigner function of a quantum state. Furthermore, Eq. (192) teaches us that the conditionally generated Wigner function $W_{f|A}(\vec{x}_f)$ can only achieve negative value when $W_A(\vec{x}_g)$ is nonpositive.

Thus, in order to conditionally prepare a state with Wigner negativity, one faces strict requirements, on both the POVM element $\hat{A}$ that is conditioned upon and on the conditional probability distribution $W(\vec{x}_g | \vec{x}_f)$ that is obtained from the initial multimode state. We discuss this point in greater detail in Sec. B. For a more experimentally inclined perspective on the production of non-Gaussian states, we refer to Ref. [81].

Before we move on to consider photon subtraction as an example of conditional creation of non-Gaussian states, let us take for a moment the opposite process: Gaussianification. The authors of Ref. [147] consider several copies of an initial non-Gaussian state, which are mixed through linear optics and subsequently some output modes are measured with on-off detectors. The conditioning is done of the events where no photons are detected, and such that we can interpret $\hat{A}$ as a projector in vacuum. By repeating several iterations of this scheme (assuming many successful conditioning events), the initial non-Gaussian state is converted into a Gaussian state. The Gaussianification process thus relies on starting from a non-Gaussian state and conditioning by projecting on a Gaussian state: the vacuum. This point of view nicely complements our approach to create non-Gaussian states.

## 2. An example: photon subtraction

Single-photon subtracted states are theoretically obtained by acting with an annihilation operator on the state. Their density matrices are given by

$$ \hat{\rho}^{-} = \frac{\hat{a}(\vec{b})\hat{\rho}\hat{a}^{\dagger}(\vec{b})}{\text{tr}[\hat{a}^{\dagger}(\vec{b})\hat{a}(\vec{b})\hat{\rho}]}, \quad (196) $$

if the photon is subtracted in one specific mode b. In practice [129,148,149], we can implement this operation on the state $\hat{\rho}$ through a mode-selective beam splitter $\hat{U} = \exp\{\theta[\hat{a}^{\dagger}(\vec{g})\hat{a}(\vec{b}) - \hat{a}^{\dagger}(\vec{b})\hat{a}(\vec{g})]\}$, that couples the mode b to an auxiliary mode g, which is prepared in a vacuum state. We thus describe the action of the beam splitter on the system of interest and the auxiliary mode as $\hat{U}(\hat{\rho} \otimes |0\rangle\langle 0|)\hat{U}^{\dagger}$. As a next step, we mount a photon detector on one of the output modes of the beam splitter. This detector is crucial to make sure that no information is lost, without it we would effectively trace out the mode and the beam splitter would simply induce losses. In contrast, we condition on the specific events where the detector counts a single
photon, we generate the state

$$ \hat{\rho}_{\theta}^{-} = \frac{\mathrm{tr}_{\mathrm{g}}[\hat{U}(\hat{\rho} \otimes |0\rangle\langle0|) \hat{U}^{\dagger}(1 \otimes |1\rangle\langle1|)]}{\mathrm{tr}[\hat{U}(\hat{\rho} \otimes |0\rangle\langle0|) \hat{U}^{\dagger}(1 \otimes |1\rangle\langle1|)]} \quad (197) $$

The reader can now recognize Eq. (187). As a next step, we assume that the beam splitter is transmitting nearly all the incoming light, such that $\theta \to 0$. We can then approximate $\hat{U} \approx 1 + \theta[\hat{a}^{\dagger}(\vec{g})\hat{a}(\vec{b}) - \hat{a}^{\dagger}(\vec{b})\hat{a}(\vec{g})]$. Then, when we insert this approximation in the expression for $\hat{\rho}_{\theta}^{-}$, we find that only the terms proportional to $\theta^2$ survive such that

$$ \hat{\rho}^{-} = \lim_{\theta \to 0} \hat{\rho}_{\theta}^{-} = \frac{\hat{a}(\vec{b})\hat{\rho}\hat{a}^{\dagger}(\vec{b})}{\mathrm{tr}[\hat{a}^{\dagger}(\vec{b})\hat{a}(\vec{b})\hat{\rho}]} \quad (198) $$

A much more detailed analysis of multimode photon sub-traction with imperfect mode selectivity can be found in Ref. [150]. We note that through this approach, photon subtraction can be understood as a weak measurement of the number of photons [151].

We can now derive the Wigner function of a single-photon-subtracted state through Eq. (193) by following the idea of Eq. (197). We initially start from a Gaussian state with covariance matrix $V_f$ and one auxiliary mode that is prepared in the vacuum

$$ V_{\text{ini}} = \begin{pmatrix} V_f & 0 \\ 0 & 1 \end{pmatrix}. \quad (199) $$

We then implement a mode-selective beam splitter that mixes one specific mode **b** with the auxiliary vacuum mode, following the scheme outlined in Fig. 4. An effective way to describe such a transformation is by designing a new mode basis $\mathcal{B}$, which has **b** as one of the modes in the mode basis. We complete the basis with complementary modes $b_1^c, \dots, b_{m-1}^c$, such that the modes basis of phase space is given by $\mathcal{B} = \{\vec{b}_1^c, \Omega\vec{b}_1^c, \dots, \vec{b}_{m-1}^c, \Omega\vec{b}_{m-1}^c, \vec{b}, \Omega\vec{b}\}$. Thus, we can perform such a basis change as

$$ \begin{pmatrix} V_f & 0 \\ 0 & 1 \end{pmatrix} \mapsto \begin{pmatrix} O_{\mathcal{B}}^T & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} V_f & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} O_{\mathcal{B}} & 0 \\ 0 & 1 \end{pmatrix}, \quad (200) $$

where the matrix of basis change is given by

$$ O_{\mathcal{B}} = \begin{pmatrix} | & | & | & | & | & | \\ \vec{b}_1^c & \Omega\vec{b}_1^c & \dots & \vec{b}_{m-1}^c & \Omega\vec{b}_{m-1}^c & \vec{b} & \Omega\vec{b} \end{pmatrix}. \quad (201) $$

It is now instructive to explicitly write the rows and columns corresponding to mode **b**:

$$ O_{\mathcal{B}}^{T} V_{\mathrm{f}} O_{\mathcal{B}} = \begin{pmatrix} V_{\mathrm{f}}^{c} & V_{\mathrm{f}}^{cb} \\ V_{\mathrm{f}}^{bc} & V_{\mathrm{f}}^{b} \end{pmatrix}. \quad (202) $$

Note that $V_f^b$ is the $2 \times 2$ matrix that describes the initial state covariances of mode **b**, while $V_f^c$ is the $(m-1) \times (m-1)$ that describes all the covariances in the complementary modes. The rectangular matrices $V_f^{cb}$ and $V_f^{bc}$ contain all the correlations between the mode **b** and the complementary modes in the basis. Now, we mix the mode **b** and the auxiliary vacuum mode on a beam splitter. This beam splitter is implemented by the transformation

$$ V_{\mathrm{BS}}^{(\mathcal{B})} = O_{\mathrm{BS}}^{(\mathcal{B})} \begin{pmatrix} V_{\mathrm{f}}^{c} & V_{\mathrm{f}}^{cb} & 0 \\ V_{\mathrm{f}}^{bc} & V_{\mathrm{f}}^{b} & 0 \\ 0 & 0 & 1 \end{pmatrix} O_{\mathrm{BS}}^{(\mathcal{B})}^{T}, \quad (203) $$

where $O_{\mathrm{BS}}^{(\mathcal{B})}$ is given by

$$ O_{\mathrm{BS}}^{(\mathcal{B})} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \theta 1 & -\sin \theta 1 \\ 0 & \sin \theta 1 & \cos \theta 1 \end{pmatrix}. \quad (204) $$

As a final step, we change the basis back to the original basis, such that the final state's covariance matrix becomes

$$ V = \begin{pmatrix} O_{\mathcal{B}} & 0 \\ 0 & 1 \end{pmatrix} V_{\mathrm{BS}}^{(\mathcal{B})} \begin{pmatrix} O_{\mathcal{B}}^{T} & 0 \\ 0 & 1 \end{pmatrix}. \quad (205) $$

We can now rewrite this entire transformation such that the matrix $V$ in Eq. (185) is given by

$$ V = O_{\mathrm{BS}} V_{\mathrm{ini}} O_{\mathrm{BS}}^{T}, \quad (206) $$

with

$$ \begin{aligned} O_{\mathrm{BS}} &= \begin{pmatrix} O_{\mathcal{B}} & 0 \\ 0 & 1 \end{pmatrix} O_{\mathrm{BS}}^{(\mathcal{B})} \begin{pmatrix} O_{\mathcal{B}}^{T} & 0 \\ 0 & 1 \end{pmatrix} \\ &= \begin{pmatrix} (\cos \theta - 1) BB^{T} + 1 & \sin \theta B \\ -\sin \theta B^{T} & \cos \theta 1 \end{pmatrix}. \end{aligned} \quad (207) $$

We introduce the $2m \times 2$ matrix $B$, which implements the mode selectivity of the beam splitter in mode **b** and is defined as

$$ B = \begin{pmatrix} | & | \\ \vec{b} & \Omega\vec{b} \end{pmatrix}. \quad (208) $$

Hence, we can simply use $O_{\mathrm{BS}}$ as a mode-selective beam splitter that mixes one specific mode of a multimode state with the auxiliary mode. We should highlight that $O_{\mathrm{BS}}$ ultimately turns out to be independent of the complementary modes $b_1^c, \dots, b_{m-1}^c$. This means that the finer details of the interferometer $O_{\mathcal{B}}$ are not important for the final $O_{\mathrm{BS}}$, the key point is that $O_{\mathcal{B}}$ changes towards a mode basis in which $\vec{b}$ and $\Omega\vec{b}$ are basis vectors of the phase space.

For the particular case of photon subtraction, we consider a very weak beam splitter, such that we consider the limit $\theta \to 0$. In this case, we can express the conditional
<figure><img src="image_6.png" /><figcaption>FIG. 4. Schematic representation of an implementation of mode-selective photon subtraction. See main text for details. For illustration, the initial state on the left is a product of single-mode squeezed vacuum states, but the protocol can in principle be applied to any Gaussian state.</figcaption></figure>

mean field, Eq. (195), and covariance matrix, Eq. (194), of the auxiliary mode g by

$$ \begin{aligned} \vec{\xi}_{g|\vec{x}_f} & \approx \theta(V_{bf} - B^T)V_f^{-1}(\vec{x}_f - \vec{\xi}_f) + \theta\vec{\xi}_b \\ & = \theta B^T(1 - V_f^{-1})(\vec{x}_f - \vec{\xi}_f) + \theta\vec{\xi}_b \end{aligned} \quad (209) $$

$$ \begin{aligned} V_{g|\vec{x}_f} & \approx 1 + \theta^2 [V_b - 1 - (V_{bf} - B^T)V_f^{-1}(V_{fb} - B)], \\ & = 1 + \theta^2 (1 - B^T V_f^{-1} B), \end{aligned} \quad (210) $$

where we introduce the matrices $V_b = B^T V_f B$, $V_b f = B^T V_f f$, and $V_{fb} = V_f B$ as well as the vector $\vec{\xi}_b = B^T \vec{\xi}_f$. We can then use these quantities to evaluate that

$$ W(\vec{x}_g | \vec{x}_f) \approx \frac{e^{-\frac{1}{2} \|\vec{x}_g\|^2}}{(2\pi)^m} \left( 1 + \theta \vec{x}_g^T B^T (1 - V_f^{-1}) (\vec{x}_f - \vec{\xi}_f) \right. \\ \left. + \theta \vec{x}_g^T \cdot \vec{\xi}_b + \frac{\theta^2}{2} \left[ (\vec{x}_g^T \cdot B^T (1 - V_f^{-1}) (\vec{x}_f - \vec{\xi}_f) \right. \right. \\ \left. \left. + \vec{x}_g^T \cdot \vec{\xi}_b \right)^2 - \|B^T(1 - V_f^{-1})(\vec{x}_f - \vec{\xi}_f) + \vec{\xi}_b\|^2 \right. \\ \left. + \vec{x}_g^T (1 - B^T V_f^{-1} B) \vec{x}_g \right] + \mathcal{O}(\theta^3) \right). \quad (211) $$

As a next step, we must choose a POVM element $\hat{A}$ to measure. In the case of photon subtraction, we mount a photon counter on the auxiliary mode and for single-photon subtraction we condition on the event where this detector detects exactly one photon. Because we use a very weakly reflective beam splitter, the probability of obtaining such an event is small but when it occurs, we have created a photon subtracted state on the remaining modes.

On a theoretical level, mounting a photon counter and conditioning on a single photon is translated to choosing

$\hat{A} = |1_g\rangle\langle1_g|$. We already encountered the corresponding Wigner function in Eq. (118), and thus we can combine this with Eq. (211) to obtain

$$ \begin{aligned} \langle \hat{A} \rangle_{g|\vec{x}_f} &= 4\pi \int_{\mathbb{R}^2} d\vec{x}_g W_{1\vec{g}}(\vec{x}_g) W(\vec{x}_g | \vec{x}_f) \\ & \approx \frac{\theta^2}{2} (\|B^T(1 - V_f^{-1})(\vec{x}_f - \vec{\xi}_f) + \vec{\xi}_b\|^2 \\ & \qquad + \mathrm{tr}[1 - B^T V_f^{-1} B]) + \mathcal{O}(\theta^3), \end{aligned} \quad (212) $$

and in a similar fashion we find that

$$ \begin{aligned} 0 \langle \hat{A} \rangle &= 4\pi \int_{\mathbb{R}^2} d\vec{x}_g W_{1\vec{g}}(\vec{x}_g) W_g(\vec{x}_g) \\ & \approx \frac{\theta^2}{2} \left( \mathrm{tr}[V_b - 1] + \|\vec{\xi}_b\|^2 \right) + \mathcal{O}(\theta^3). \end{aligned} \quad (213) $$

The actual evaluation of these integrals is not completely straightforward. As a key idea, we use that the integral takes the form of a polynomial multiplied by a Gaussian. We can thus evaluate the expectation value of the polynomial with respect to this Gaussian distribution. In essence, this boils down to calculating a set of moments of a Gaussian probability distribution. We see rather quickly that the lowest orders in $\theta$ vanish, such that the leading order is $\theta^2$. Putting everything together, we find that

$$ \begin{aligned} \lim_{\theta \to 0} W_{f|A}(\vec{x}_f) &= \lim_{\theta \to 0} \frac{\langle \hat{A} \rangle_{g|\vec{x}_f}}{\langle \hat{A} \rangle} W_f(\vec{x}_f) \\ &= \frac{\|B^T(1 - V_f^{-1})(\vec{x}_f - \vec{\xi}_f) + \vec{\xi}_b\|^2 + \mathrm{tr}[1 - B^T V_f^{-1} B]}{\mathrm{tr}(V_b - 1) + \|\vec{\xi}_b\|^2} W_f(\vec{x}_f). \end{aligned} \quad (214) $$

As such, we obtain the Wigner function for a multi-mode photon-subtracted state. This Wigner function can be obtained using several different methods, ranging from algebraic [152,153] to analytical [154]. The difference between those approaches and our method here is that we do not directly use the properties of the annihilation operator, but rather model the exact experimental setup, while relying entirely on phase-space representations.

The methods presented here for treating photon-subtracted states can straightforwardly be extended to the subtraction of multiple photons in different modes and we can easily replace photon-number-resolving detection with an on-off detector by setting $\hat{A} = 1 - |0\rangle\langle 0|$. The techniques used in the calculations remain essentially the same and it yields the same result in the $\theta \to 0$ limit (doing this calculation may prove to be a good exercise for the motivated reader). However, any real implementation of a photon-subtraction experiment will use a beam splitter with finite reflectivity, such that there will be a
difference between on-off detectors and photon-number-resolving detectors due to the small contributions of higher order terms in $\theta$. In practice, one chooses the reflectivity of the beam splitter with respect to the energy content of the initial state to effectively suppress all higher-order terms in Eq. (211). In the single-mode case, an early thorough analysis of the implementation of photon subtraction can be found in Ref. [155]. There are also proposals in the literature to use a photon-subtraction setup with larger values of $\theta$ to gain an additional advantage in quantum state preparation [156–158].

As a final note, we point out that a similar treatment can be used to describe photon-added states, which are also relevant in experiments [159,160]. It is perhaps surprising that such a state can be obtained by performing a measurement on a part of a Gaussian state, but it suffices to replace the beam splitter in Eq. (197) with a two-mode squeezer. In other words, we set $\hat{U} = \exp\{\theta[\hat{a}^\dagger(\vec{g})\hat{a}^\dagger(\vec{b}) - \hat{a}(\vec{b})\hat{a}(\vec{g})]\}$, and consider again the limit where the parameter $\theta$ is small, i.e., weak squeezing. Even though this is a simple step in theory, it is much harder in an actual experimental setting. Photon subtraction can be implemented with a passive linear optics element, while photon addition always requires squeezing and thus a nonlinear optics implementation.

# V. NON-GAUSSIAN STATES AND QUANTUM CORRELATIONS

In this section, we explore the interplay between non-Gaussian effects and quantum correlations. First, in Sec. A, we provide a crash course to introduce the unfamiliar reader to the most important types of quantum correlations: entanglement, steering, and Bell nonlocality. In Sec. B we subsequently highlight how certain types of quantum correlations can be used to create certain types of non-Gaussian states via the methods of Sec. B. In Sec. C, we then explore how non-Gaussian operations can create or enhance quantum correlations by focusing on photon-subtracted states. Finally, we explore the role that is played by non-Gaussian states in Bell inequalities in Sec. D.

## A. Quantum correlations: a crash course

We start by giving a quick introduction to the different kinds of common quantum correlations. Readers who want to get a more thorough overview on these subjects are referred to Refs. [161–163] as natural starting points.

In this Tutorial, we solely consider bipartite quantum correlations. This implies that we structure the system in a similar way as in Sec. B and divide the m-mode system in two parts, each with their own phase space, i.e., $\mathbb{R}^{2m} = \mathbb{R}^{2l} \oplus \mathbb{R}^{2l'}$. It is noteworthy that the corresponding Fock space takes the structure $\Gamma(\mathcal{H}_m) = \Gamma(\mathcal{H}_l) \otimes \Gamma(\mathcal{H}_{l'})$, where we again use the mapping, Eq. (44), between the phase space $\mathbb{R}^{2k}$ and the k-dimension Hilbert space $\mathcal{H}_k$.

These structures are crucial to understand quantum correlations.

## 1. Correlations

To better understand quantum correlations, it is useful to start by generally defining what a correlations is. In a statistical sense, two stochastic variables X and Y are correlated when the expectation values have the following property:

$$ \mathbb{E}(XY) \neq \mathbb{E}(X)\mathbb{E}(Y). \qquad (215) $$

This can be translated to the level of probability distributions by stating that the joint probability distribution for outcomes X = x and Y = y is not the product of the marginals

$$ P(x,y) \neq P(x)P(y), \qquad (216) $$

where

$$ P(x) = \int_{\mathcal{Y}} dy P(x,y), \quad \text{and} \quad P(y) = \int_{\mathcal{X}} dx P(x,y). \qquad (217) $$

Here, $\mathcal{X}$ and $\mathcal{Y}$ denote the possible outcomes of the stochastic variables X and Y, respectively [164].

When we talk about quantum systems, there are generally many observables that can be considered. When we consider a global multimode system with phase space $\mathbb{R}^{2m}$ and two subsystems with phase spaces $\mathbb{R}^{2l}$ and $\mathbb{R}^{2l'}$, there is a whole algebra of observables involved. The role of the stochastic observables X and Y will be taken up by local observables $\hat{X}$ and $\hat{Y}$ that are contained in the observable algebra generated by, respectively, $\hat{q}(\vec{f})$ and $\hat{q}(\vec{g})$, with $\vec{f} \in \mathbb{R}^{2l}$ and $\vec{g} \in \mathbb{R}^{2l'}$. These local observables are correlated when

$$ \mathrm{tr}(\hat{X} \otimes \hat{Y}\hat{\rho}) \neq \mathrm{tr}(\hat{X}\hat{\rho}_f)\mathrm{tr}(\hat{Y}\hat{\rho}_g), \qquad (218) $$

where $\hat{\rho}_f$ and $\hat{\rho}_g$ are the marginals (or reduced states) of $\hat{\rho}$ for the subsystems $\mathbb{R}^{2l}$ and $\mathbb{R}^{2l'}$.

When we talk about correlated systems rather than correlated observables, we consider that there exists a pair of local observables such that Eq. (218) holds. Thus, if two systems are not correlated, it follows that for all possible observables $\mathrm{tr}(\hat{X} \otimes \hat{Y}\hat{\rho}) = \mathrm{tr}(\hat{X}\hat{\rho}_f)\mathrm{tr}(\hat{Y}\hat{\rho}_g)$. This lack of correlations can be expressed on the level of the quantum state by the identity $\hat{\rho} = \hat{\rho}_f \otimes \hat{\rho}_g$. On the level of Wigner functions, we can therefore say that a state contains correlations if the Wigner function satisfies

$$ W(\vec{x}_f \oplus \vec{x}_g) \neq W_f(\vec{x}_f)W_g(\vec{x}_g), \qquad (219) $$

where the marginal Wigner functions are defined as in Eqs. (182), (183).
It is clear that correlations between systems can occur, both, in the context of classical probability theory and in quantum theory. However, we already established that quantum physics imposes additional constraints on the statistics of observables, which ultimately make it impossible to describe CV quantum systems in terms of probability distributions on phase space. Similarly, quantum physics leads to new features for the correlations of subsystems. Thus, in our study of quantum correlations we explore correlated systems, in the sense of Eq. (219), and we seek to differentiate between correlations that are of classical origin and those that can be attributed to a quantum origin.

## 2. Quantum entanglement

Quantum entanglement is probably the most well-known type of quantum correlation. The notion of entanglement derives directly from the structure of the quantum state space and is related to the contrast between pure states in classical and quantum physics.

To understand this contrast, we loosely follow the idea of Ref. [165]. Let us be a bit more precise as to what is meant with pure states in classical physics in the context of CV systems. Classically, in a context of statistical mechanics, any CV system can be described by a probability distribution on phase space. From a mathematical point of view, this means that the space containing all the possible classical states is a convex set because any convex combination of two probability distributions is again a probability distribution. Pure states are formally defined as the extreme points of the convex set, i.e., the states that cannot be decomposed as being a convex combination of two other states. In a classical theory, where states can unambiguously be represented by probability distributions on phase space, the pure states are delta functions centered on the different points of phase space. From a physical point of view, this corresponds to the intuition that pure states are “the least noisy” states, which simply corresponds to a single point in phase space.

For our phase space $\mathbb{R}^{2m} = \mathbb{R}^{2l} \oplus \mathbb{R}^{2l'}$ these delta functions factorize with respect to the subsystems, i.e., $\delta(\vec{x} - \vec{x}') = \delta(\vec{x}_f - \vec{x}'_f)\delta(\vec{x}_g - \vec{x}'_g)$, with $\vec{x}, \vec{x}' \in \mathbb{R}^{2m}$, $\vec{x}_f, \vec{x}'_f \in \mathbb{R}^{2l}$, and $\vec{x}_g, \vec{x}'_g \in \mathbb{R}^{2l'}$. In the light of Eq. (216) we thus conclude that pure states of classical systems are always uncorrelated [166]. Any correlations that are present in classical states are thus obtained by taking a convex combination of uncorrelated pure states.

In quantum systems, pure states are represented by state vectors in a Hilbert space (in our case Fock space). They also can be seen as the extreme points of a convex set of states that contains all density matrices $\hat{\rho}$. As we saw in the example where we discussed the Hong-Ou-Mandel state $|\text{HOM}\rangle$ in Eq. (165), pure quantum states can actually be correlated in the sense of Eq. (219). This crucial difference

between classical and quantum pure states lies at the basis of quantum entanglement.

The notion of entanglement derives directly from the structure of the quantum state and is defined as the opposite of a separable state. For pure states, separable states $|\Psi\rangle \in \Gamma(\mathcal{H}_m)$ are the pure states that are uncorrelated and can thus be written as $|\Psi\rangle = |\Psi_l\rangle \otimes |\Psi_{l'}\rangle$ with $|\Psi_l\rangle \in \Gamma(\mathcal{H}_l)$ and $|\Psi_{l'}\rangle \in \Gamma(\mathcal{H}_{l'})$. All other pure states are said to be entangled. They possess correlations that are not due to some type of convex combination of uncorrelated states, something which is impossible for classical pure states.

The situation is more subtle when considering mixed states, i.e., convex combinations of pure states. Convex mixtures of classical pure states can also show correlations, and it is therefore crucial to make a distinction between this type of classical correlations and quantum correlations. Due to the structure of classical pure states, we find that any classical joint probability distribution on phase space can be written as a convex combination of local probability distributions

$$P(\vec{x}_f \oplus \vec{x}_g) = \int d\gamma p(\gamma) W_f^{(\gamma)}(\vec{x}_f) W_g^{(\gamma)}(\vec{x}_g), \quad (220)$$

where $\gamma$ is some arbitrary way of labeling states, distributed according to distribution $p(\gamma)$. This notion of classical correlations can directly be generalized to quantum states [167], and thus a mixed state is said to be separable when all of its correlations are classical, i.e., when it is a convex mixture of product states

$$\hat{\rho} = \int d\gamma p(\gamma) \hat{\rho}_f^{(\gamma)} \otimes \hat{\rho}_g^{(\gamma)}. \quad (221)$$

In the language of Wigner functions, the separability condition translates to

$$W(\vec{x}_f \oplus \vec{x}_g) = \int d\gamma p(\gamma) W_f^{(\gamma)}(\vec{x}_f) W_g^{(\gamma)}(\vec{x}_g), \quad (222)$$

where we again use the definitions of Eqs. (182), (183). Quantum states that cannot be described by a Wigner function of the form Eq. (222) are not separable and are said to be entangled.

Hence, quantum entanglement describes the origin of the quantum correlations rather than their properties. Nevertheless, the set of separable states is a convex set and thus the Hahn-Banach separation theorem [168,169] teaches us that it is in principle possible to use observables to distinguish between separable and entangled states. In this sense the difference between entangled and separable states is measurable. For the sake of uniformity, we highlight that separable states lead to the following measurement
statistics of local observables $\hat{X}$ and $\hat{Y}$:

$$P(x,y) = \int d\gamma p(\gamma) P_{\hat{\rho}}^{(\gamma)}(x) P_{\hat{\rho}}^{(\gamma)}(y). \quad (223)$$

It is crucial to emphasize that the distributions of measurement outcomes $P_{\hat{\rho}}^{(\gamma)}(x)$ and $P_{\hat{\rho}}^{(\gamma)}(y)$ are governed by the laws of quantum physics. Formally, we can use the spectral theorem to write

$$\hat{X} = \int_{\mathcal{X}} dx \, x \hat{E}_x, \quad \text{and} \quad \hat{Y} = \int_{\mathcal{Y}} dy \, y \hat{E}_y, \quad (224)$$

such that $\hat{E}_x$ and $\hat{E}_y$ are the POVM elements that correspond to the measurement outcomes $x$ and $y$ for the measurement of the (generalized) observables $\hat{X}$ and $\hat{Y}$, respectively. The probability distribution $P_{\hat{\rho}}^{(\gamma)}(x)$ is then given by

$$P_{\hat{\rho}}^{(\gamma)}(x) = \text{tr}[\hat{E}_x \hat{\rho}_{\hat{f}}^{(\gamma)}] = (4\pi)^l \int_{\mathbb{R}^{2l}} d\vec{x}_f W_{E_x}(\vec{x}_f) W_{\hat{f}}(\vec{x}_f), \quad (225)$$

and analogously for $P_{\hat{\rho}}^{(\gamma)}(y)$.

For separable states, Eq. (223) with local probability distribution given by Eq. (225) holds for any arbitrary pair of local observables. The model that is described by these equations is known as a local hidden variable model for quantum entanglement, where $\gamma$ is the hidden variable. We may not necessarily know the origins and behavior of $\gamma$, but the model generally captures two important elements. First, all correlations are induced by the common variable $\gamma$ that governs the convex mixture. Second, the local probability distributions $P_{\hat{\rho}}^{(\gamma)}(x)$ and $P_{\hat{\rho}}^{(\gamma)}(y)$ have a quantum origin. For CV systems the latter point, for example, implies that these local probability distributions must satisfy the Heisenberg inequality. These quantum constraints on the local probability distributions $P_{\hat{\rho}}^{(\gamma)}(x)$ and $P_{\hat{\rho}}^{(\gamma)}(y)$ are typically useful for the falsification of the local hidden variable model, Eq. (223), and thus prove the presence of quantum entanglement [170,171].

### 3. Quantum steering

In a formal sense, quantum steering is a rather recent addition to the family of quantum correlations. Nevertheless, it is exactly this phenomenon that lies at the basis of the Einstein-Podolsky-Rosen (EPR) paradox [172]. Schrödinger's response to the EPR paper [173,174] lies at the basis of what we now call quantum steering, but the broader implications of these results were only sporadically discovered and formalized [175,176].

Just like for quantum entanglement, a system is said to be steerable if the measurement statistics cannot be explained in terms of a local hidden variable model. A

peculiarity of quantum steering is that it involves a certain directionality, where one of the subsystems is said to “steer” the other subsystem. This asymmetry is represented in the local hidden variable model, which takes the following form:

$$P(x,y) = \int d\gamma p(\gamma) P^{(\gamma)}(x) P_{\hat{\rho}}^{(\gamma)}(y), \quad (226)$$

where we emphasize the striking resemblance to Eq. (223). Note that, contrary to the case of quantum entanglement, we now allow the probability distribution $P^{(\gamma)}(x)$ for the first subsystem to be arbitrary and thus do not impose any constraints of quantum theory on it. If there exist observables $\hat{X}$ and $\hat{Y}$ for which the probability distribution is not consistent with the model, Eq. (226), the subsystems with phase space $\mathbb{R}^{2l}$ is able to steer the subsystem with phase space $\mathbb{R}^{2l'}$.

Quantum steering is perhaps most logically explained in terms of conditional states and probability distributions. For nonsteerable states, the local hidden variable model, Eq. (226), must hold for all observables, which in turn imposes conditions on the level of states. Here these conditions manifest on the level of conditional states of the type Eq. (187). To see this, we consider the conditional probability distribution associated with Eq. (226):

$$P(y|x) = \frac{\int d\gamma p(\gamma) P^{(\gamma)}(x) P_{\hat{\rho}}^{(\gamma)}(y)}{P(x)}, \quad (227)$$

where the probability to obtain a certain outcome $\hat{X} = x$ is given by

$$P(x) = \int d\gamma p(\gamma) P^{(\gamma)}(x). \quad (228)$$

Note that for any $x$ the function

$$\tilde{P}(\gamma | x) := \frac{p(\gamma) P^{(\gamma)}(x)}{P(x)} \quad (229)$$

is a well-defined probability distribution. Furthermore, if we demand that Eq. (227) holds for all observables $\hat{Y}$, we find the following condition for the conditional state:

$$\hat{\rho}_{g|\hat{X}=x} = \int d\gamma \, \tilde{P}(\gamma | x) \hat{\rho}_{g}^{(\gamma)}. \quad (230)$$

Because quantum steering is a property of the state, we again require Eq. (230) to hold for all observables $\hat{X}$ for a state to not be steerable.

The local hidden variable model, Eq. (226), and the consequence for the conditional state, Eq. (230), may seem stringent, but it is often intricate to formally prove that such a model cannot explain observed data. It turns out
that computational methods based on semidefinite programming [177] are well suited to prove that the set of all possible conditional states is inconsistent with Eq. (230). A more physical point of view is based on developing steering inequalities [178]. As a notable example, one can derive a type of conditional Heisenberg inequality for states of the form Eq. (230). The local hidden variable model, Eq. (226), assumes that the laws of quantum physics constrain the measurement statistics in the second subsystem. We can then define the conditional variance of an arbitrary observable $\hat{Y}$

$$ \Delta^2(\hat{Y} | \hat{X} = x) := \mathrm{tr}[\hat{Y}^2 \hat{\rho}_{g|\hat{X}=x}] - \mathrm{tr}[\hat{Y} \hat{\rho}_{g|\hat{X}=x}]^2, \quad (231) $$

which leads to the “average inference variance”

$$ \Delta_{\mathrm{inf}}^2(\hat{Y}) := \int_X dxP(x) \Delta^2(\hat{Y} | \hat{X} = x), \quad (232) $$

that characterizes the precision with which we can infer the measurement outcome of $\hat{Y}$, given a measurement outcome of $\hat{X}$. Under the assumption that Eq. (230) holds, we can then prove the inference Heisenberg inequality [178]

$$ \Delta_{\mathrm{inf}}^2(\hat{Y}_1) \Delta_{\mathrm{inf}}^2(\hat{Y}_2) \geq \frac{1}{2} \int_{\mathcal{X}_3} dxP(x) \left| \mathrm{tr} \left( [\hat{Y}_1, \hat{Y}_2] \hat{\rho}_{g|\hat{X}_3=x} \right) \right|^2, \quad (233) $$

where $\Delta_{\mathrm{inf}}^2(\hat{Y}_1)$ and $\Delta_{\mathrm{inf}}^2(\hat{Y}_2)$ can be conditioned on any observables $\hat{X}_1$ and $\hat{X}_2$, respectively.

Thus, whenever one performs a series of conditional measurements that violate the inference Heisenberg inequality (233), the assumption (230) cannot hold and thus the measurements in the subsystem with phase space $\mathbb{R}^{2l}$ have steered those in the subsystem with phase space $\mathbb{R}^{2l'}$. In more colloquial terms, the inequality (233) sets a limit on the precision with which classical correlations between observables can be used to infer measurement outcomes of one quantum subsystem, based on measurement outcome of the other subsystem (regardless of whether it is quantum or not). Quantum correlations allow us to outperform these bounds and provide better inference than classically possible, and this phenomenon is the essence of quantum steering.

Now let us now express Eq. (233) for quadrature operators:

$$ \Delta_{\mathrm{inf}}^2[\hat{q}(\vec{g}_1)] \Delta_{\mathrm{inf}}^2[\hat{q}(\vec{g}_2)] \geq |\vec{g}_1^T \Omega \vec{g}_2|^2, \quad (234) $$

where $\vec{g}_1, \vec{g}_2 \in \mathbb{R}^{2l'}$. As a next step, we must understand the properties of the average inference variance $\Delta_{\mathrm{inf}}^2[\hat{q}(\vec{g}_1)]$, which we obtain by conditioning on a quadrature observable in the other subsystem's phase space $\mathbb{R}^{2l'}$. More specifically let us assume that we condition on measurements of

$\hat{q}(\vec{f}_1)$, such that we must evaluate the conditional variance $\Delta^2[\hat{q}(\vec{g}_1) | \hat{q}(\vec{f}_1) = x]$. The conditional variance $\Delta^2[\hat{q}(\vec{g}_1) | \hat{q}(\vec{f}_1) = x]$ is then given by the matrix element of the covariance matrix that describes $W(\vec{x}_g | x\vec{f}_1)$ as defined in Eq. (193):

$$ \Delta^2[\hat{q}(\vec{g}_1) | \hat{q}(\vec{f}_1) = x] = \vec{g}_1^T \left[ V_g - \frac{V_g \vec{f}_1 \vec{f}_1^T V_{fg}}{\vec{f}_1^T V_{fg} \vec{f}_1} \right] \vec{g}_1, \quad (235) $$

because the quantity does not depend on the actual outcome that is postselected upon, we find that

$$ \Delta_{\mathrm{inf}}^2[\hat{q}(\vec{g}_1)] = \vec{g}_1^T \left[ V_g - \frac{V_g \vec{f}_1 \vec{f}_1^T V_{fg}}{\vec{f}_1^T V_{fg} \vec{f}_1} \right] \vec{g}_1, \quad (236) $$

$$ \Delta_{\mathrm{inf}}^2[\hat{q}(\vec{g}_2)] = \vec{g}_2^T \left[ V_g - \frac{V_g \vec{f}_2 \vec{f}_2^T V_{fg}}{\vec{f}_2^T V_{fg} \vec{f}_2} \right] \vec{g}_2. \quad (237) $$

From Eqs. (236) and (237) we can deduce that $\Delta_{\mathrm{inf}}^2[\hat{q}(\vec{g})] \geq \vec{g}^T V_{g|\vec{x}_g} \vec{g}$ for all $\vec{g} \in \mathbb{R}^{2l'}$ regardless of the $\hat{q}(\vec{f})$ that is conditioned upon. Thus, if $V_{g|\vec{x}_g}$ satisfies the Heisenberg inequality the inference Heisenberg inequality (234) is also satisfied.

The setting with homodyne measurements, or more general Gaussian measurements, is close to the system that is discussed in Ref. [172]. For this reason, we refer to quantum steering with Gaussian measurements as EPR steering in contrast to more general quantum steering. This type of steering has been studied extensively in the literature, e.g., [176,179–181] and will be a key element in Sec. 2.

Note that both subsystems clearly play a very different role in this setting. The first subsystem simply produces measurement results of different observables. The information of these measurements in the first subsystem is then used to infer measurement results in the second subsystem, which is assumed to be a quantum system. In a quantum communication context, this asymmetry corresponds to a level of trust: we position ourselves in the steered system and trust that our system is a well-behaved quantum system, but we do not trust the party that controls the other subsystem (up to a point where we do not even want to assume that the data that are communicated to us come from an actual quantum system). The violation of a steering inequality practically allows verification in such a setting that there is indeed a quantum correlation between the two subsystems [182].

The inference Heisenberg inequality (233) shows that quantum steering describes certain properties of the quantum correlations. States that can perform quantum steering thus possess correlations that can be used to infer measurement outcomes better than any classical correlations
could. These correlations cannot be described by a hidden variable model of the form Eq. (226), which is more general than the model, Eq. (223). Thus, all states that produce statistics consistent with Eq. (223) are also consistent with Eq. (226) such that states that can perform steering must be entangled. However, there are states that produce statistics that is consistent with Eq. (226), but inconsistent with Eq. (223). In other words, not all entangled states can be used to perform quantum steering. In this sense, quantum steering can be said to be “stronger” than quantum entanglement.

### 4. Bell nonlocality

To date, the seminal work of John S. Bell on the Einstein-Podolsky-Rosen paradox [183] is probably one of the most remarkable findings on the foundations of quantum physics. What most had long taken for granted, the existence of local hidden variables to explain the probabilistic nature of quantum physics, turned out to be inconsistent with the theoretical quantum formalism. It is here that we find the real historical origin of the concept of quantum correlations as something fundamentally different from classical ones.

As for quantum entanglement and steering, the story of Bell nonlocality starts from a local hidden variable model that bears a strong resemblance to Eqs. (223) and (226). In this case, the model attempts to describe the joint measurement statistics of $\hat{X}$ and $\hat{Y}$ as

$$P(x,y) = \int d\gamma p(\gamma) P^{(\gamma)}(x) P^{(\gamma)}(y). \quad (238)$$

The key observation is that now all the quantum constraints on the probability distributions have been dropped and both the $P^{(\gamma)}(x)$ and $P^{(\gamma)}(y)$ can be any mathematically well-defined probability distributions. Even though the difference between Eqs. (223) and (226) on the one hand, and Eq. (238) on the other hand, may appear small, the impact of dropping the constraints on the local distributions is enormous. Think, for example, of the Hahn-Banach separation theorem that is invoked to define entanglement witnesses, this crucially relies on the Hilbert-space structure of the state space. Think for example of (233) which crucially depends on the fact that quantum probabilities are constrained by the Heisenberg inequality. Abandoning all connections that tie probabilities to operator algebras on Hilbert spaces deprives quantum mechanics of their toolbox. Nevertheless, it turns out that some quantum states induce statistics that is inconsistent with Eq. (238).

Again, we note that states that can be described by the models, Eqs. (223) or (226), can also be described by the model, Eq. (238). Bell's local hidden variable model, Eq. (238), is thus the most general one and the class of states that lead to measurement statistics that cannot be described by it is the smallest. Therefore, we say that the correlations

that lead to a violation of the mode, Eq. (238), also known as Bell nonlocality, are the strongest types of quantum correlations.

The key insight of Bell's work [183,184] is that Eq. (238) puts constraints on the correlations of different combinations of observables in the subsystems. These constraints, cast in the form of Bell inequalities can be violated by certain quantum states. The inconsistency of quantum physics with the model, Eq. (238), can in itself be seen as a special case of contextuality [185]. Over the decades, many different kinds of Bell inequalities have been derived (see, for example, Ref. [186]). Here we restrict to presenting one of the most commonly used incarnations: the Clauser-Horne-Shimony-Holt (CHSH) inequality [187]. This inequality relies on the measurement of four observables: $\hat{X}$ and $\hat{X}'$ on the first subsystem and $\hat{Y}$ and $\hat{Y}'$ on the second subsystem. Furthermore, we consider that the observables can take two possible values: -1 or 1. Assuming the model in Eq. (238) it is then possible to derive

$$|\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, \quad (239)$$

where $\langle\cdot\rangle$ denotes the expectation value. In this Tutorial we skip the derivation of this result, but the interested reader is referred to Ref. [188] for a detailed discussion. Remarkably, certain highly entangled states can violate this inequality.

The experimental violation of Bell's inequalities formally shows that quantum correlations are profoundly different than classical correlations [189–192]. However, one needs clever combinations of several observables in both subsystems to actually observe the difference. With most experimental loopholes now closed [193–197], Bell inequalities can now in principle be used to impose a device-independent level of security on various quantum protocols [198].

As a concluding remark, it is interesting to highlight the existence of a semidevice-independent framework for testing quantum correlations [199,200]. The key idea is that nothing is assumed about the measurement devices nor about the states, much like in the scenario of Bell inequalities. Yet, in the framework of Refs. [199,200] one does add an additional level of trust in the sense that one assumes that the inputs of the measurement device can be controlled and trusted. In a way, this additional intermediate level of trust is somewhat reminiscent of quantum steering. This framework was very recently extended to the CV setting [201].

### B. Non-Gaussianity through quantum correlations

In Sec. B, we explained how conditional operations can be used to create non-Gaussian quantum states. The presence of correlations plays an essential role in this