A Derivation truncated correlation functions

# Statistical signatures of multimode single-photon added and subtracted states of light

## Abstract

The addition or subtraction of a photon from a Gaussian state of light is a versatile and experimentally feasible procedure to create non-Gaussian states. In multimode setups, these states manifest a wide range of phenomena when the photon is added or subtracted in a mode-tunable way. In this contribution, we derive the truncated correlations, which are multimode generalisations of cumulants, between quadratures in different modes as statistical signatures of these states. These correlations are then used to obtain the full multimode Wigner function, the properties of which are subsequently studied. In particular we investigate the effect of impurity in the subtraction or addition process, and evaluate its impact on the negativity of the Wigner function. Finally, we elaborate on the generation of inherent entanglement through subtraction or addition of a photon from a pure squeezed vacuum.

\LetLtxMacro\ORIGselectlanguage\ORIGselectlanguage

english

## I Introduction

Continuous variable (CV) quantum optics has ample advantages for quantum information processing. The most notable strength of general optical systems is their resilience against decoherence, which proves useful for quantum protocols. In CV quantum optics, states with arbitrary many entangled modes (normalised solutions to Maxwell’s equations) can be deterministically generated Chen et al. (2014); Yoshikawa et al. (2016). However, these experimentally generated states are Gaussian, i.e. they can be described by a multivariate Gaussian probability distribution on the optical phase space (for formal details, see Section II). Because Gaussian statistics can easily be simulated with classical computation resources Rahimi-Keshari et al. (2016), the use of these states in quantum computation is limited.

To reach full universal quantum computation, CV setups require at least one non-Gaussian ingredient. There have been several theoretical and experimental proposals to achieve this, ranging from ancillary Gottesman-Kitaev-Preskill (GKP) states Menicucci (2014) to specific non-Gaussian gates, e.g. Arzani et al. (2017). Within this paper, we focus on photon addition and subtraction as de-Gaussification techniques Parigi et al. (2007). Such techniques have, for example, proven their worth for entanglement distillation Ourjoumtsev et al. (2007); Takahashi et al. (2010); Navarrete-Benlloch et al. (2012). In particular the subtraction of a photon is in essence a simple procedure which, as originally proposed Dakna et al. (1997), only requires a beamsplitter and a photodetector. However, because beamsplitters are not mode-selective, this simple photon subtraction scheme will increase the impurity of the state. To avoid such incoherent mixing of modes, theory for coherent mode-selective photon subtraction was recently developed Averchenko et al. (2014, 2016). Considerable steps have already been undertaken to implement this coherent mode-dependent photon subtraction in a quantum frequency comb Ra et al. (2017).

In this paper, we extend the theoretical framework for such multimode photon-added and -subtracted states. Our central achievement is the derivation of the general Wigner function Braunstein and van Loock (2005); Weedbrook et al. (2012) for these states. Even though the Wigner functions for (multi)photon-added and -subtracted states are known in several specific setups, e.g. Agarwal and Tara (1991); Dakna et al. (1997); Kim et al. (2005); Biswas and Agarwal (2007); Kim (2008), general multimode results were still lacking, even in the case of single-photon addition and subtraction.

We approach this problem from a statistical mechanics perspective, by deriving truncated correlation functions Robinson (1965a); Petz (1990); Bratteli and Robinson (1997); Verbeure (2011) for these particular case of mode-selective photon addition and subtraction from multimode Gaussian states. Truncated correlations as such are useful witnesses for the Gaussianity of states, but they are also connected to phase-space representations. Specifically, we employ the truncated correlations to derive the characteristic function, which upon Fourier transformation gives us the Wigner function –this key result is shown in (50). In the remainder of the work, we investigate the negativity of this Wigner function, which is an important indicator of the non-classicality of the state Hudson (1974); Soto and Claverie (1983); Mandilara et al. (2009); Kenfack and Życzkowski (2004); Mari and Eisert (2012) from a quantum probability theory perspective. Finally, we also investigate the entanglement properties that can be deduced from the Wigner function, which are ultimately the features that we want to exploit in future application in quantum technologies. This work elaborates on the details behind Walschaers et al. (2017) and generalises the results to non-pure photon addition and subtraction.

The paper is structured in three major parts. In Section II, we introduce the mathematical formalism and concepts which mix techniques from quantum statistical mechanics Petz (1990); Bratteli and Robinson (1997); Verbeure (2011) and quantum optics Braunstein and van Loock (2005); Weedbrook et al. (2012). These techniques are applied in Section III to investigate multimode mode-selective photon addition and subtraction. To make our abstract results more concrete, we finally study two examples in Section IV: the subtraction and addition from the two-mode symmetrically squeezed vacuum, and from an experimentally obtained state. The latter is an extension of the results of Walschaers et al. (2017).

## Ii Multimode quantum optics on phase space

### ii.1 Optical phase space and quadrature operators

The study of continuous variable multimode quantum optics is in essence a study of quantum physics in a high-dimensional phase space. For optical systems, the relevant phase space is generated by the real and imaginary parts of the contributing electric fields, the amplitude and phase quadratures, respectively.

The modal structure of light is essential in the present contribution. A mode is simply a normalised 1 solution to Maxwell’s equation, which has both a spatial and a temporal structure, as indicated by the arguments and , respectively. A general complex electric field can than be represented in terms of a mode basis as

 E(r,t)=m∑j=1(xj+ipj)uj(r,t), (1)

where the and are the amplitude and phase quadratures, respectively. Therefore, every vector can be associated with a set of phase and amplitude quadratures in the specific mode basis. The vectors space which is generated in this way is the optical phase space. As such, any vector represents a classical electromagnetic field. When, in addition, is normalised, this classical electric field is associated with a new mode. Thus, it is useful to introduce

 N(R2m)={f∈R2m∣∥f∥=1}, (2)

to describe modes. However, the dimension of is larger than the number of modes . This is a consequence of the complex amplitude of the field, which associates two quadratures to every mode. To faithfully reproduce the properties of these complex amplitudes in (1), the phase and amplitude quadratures are connected through a symplectic structure represented by a matrix which acts on the optical phase space , and has the properties:

 J2=−1, (3) (f1,Jf2)=−(f2,Jf1)for all f1,f2∈R2m, (4)

where denotes the standard inner product on . This structure implies that the optical phase space is a phase space as also studied in analytical mechanics. An important consequence of (4) is that for every . The orthogonal vectors are associated with the same mode , such that the space generated by and is the two-dimensional phase space which describes all possible electromagnetic fields in mode .

One can always construct an orthonormal symplectic basis of the optical phase space, where is the basis vector which generates the phase space axis which denotes the amplitude quadrature of mode , whereas generates the associated phase quadrature. The symplectic basis of the optical phase space is directly associated with a mode basis . Hence, a change of basis in the optical phase space implies a change in mode basis.

When we combine the above optical phase space with the framework of statistical mechanics, we can describe classical optics setups. However, to treat problems in multimode quantum optics, we must go through the procedure of canonical quantisation. To do so, we associate a quadrature operator to each . These operators fulfil the crucial mathematical property

 Q(af1+bf2)=aQ(f1)+bQ(f2), (5)

for all and with . This property implies that the operator is independent of the basis chosen to express . Moreover, these operators are governed by the canonical commutation relation von Neumann (1931); Petz (1990):

 [Q(f1),Q(f2)]=−2i(f1,Jf2)for all f1,f2∈N(R2m). (6)

We have defined (6) such that the operator corresponds to a quadrature operator with the shot noise equal to one.

The linearity condition (5) can be extended to all to define operators for non-normalised This generalisation does not lead to any mathematical problems and (6) still holds. Physically, such different norms of can be associated with rescaled quadrature measurements. In this article, the generalised quadratures will be used to limit notational overhead in the definition of the displacement operator (7).

### ii.2 Representing quantum states

Because quantum physics is a statistical theory, we require a mathematical object to describe the statistics of measurements: the quantum state. We focus on systems which can accurately be represented in a Hilbert space , on which is an (unbounded) operator. The quantum state can then be represented by a density operator , which is positive and has 2.

However, the density operator is not the most practical tool to characterise a state of a continuous variable system. Quasi-probability distributions on phase space are a common and practical alternative, not only due to their importance to interpret the fundamental physics of the state, but also because they can be measured experimentally. Throughout this paper, we will particularly emphasise the Wigner function as an important tool, because –at least for small mode numbers– it can be experimentally reconstructed through tomographic methods.

From the mathematical point of view, we start by constructing the characteristic function in order to derive the Wigner function. To do so, we first define the displacement operator

 D(α)≡exp(−iQ(Jα)/2),α∈R2m. (7)

Importantly, is generally not normalised as its norm dictates the distance of the displacement. Indeed, this operator’s action on a quadrature operator is given by

 D(−α)Q(f)D(α)=Q(f)+(f,α), (8)

such that the strength and direction of the displacement are given by and respectively. Via the displacement operator, we can introduce characteristic function as

 χ(α)≡tr{ρD(2Jα)}=tr{ρexp(iQ(α))}. (9)

Subsequently, one can construct the Wigner function by a multi-dimensional Fourier transformation

 W(β)=1(2π)2m∫R2mdαχ(α)e−i(α,β),for β∈R2m, (10)

where indicates a point in phase space, hence it needs not be normalised. The Wigner function shares the normalisation properties of a probability distribution and its marginals are probability measures. However, the full Wigner function is merely a quasi-probability distribution in the sense that it can assume negative values for some regions of phase space. It is this negativity which sets quantum mechanics apart from classical probabilistic theories on phase space. As such, negativity can be seen as a genuine sign of “quantumness”, which was also associated with quantum supremacy Mari and Eisert (2012); Rahimi-Keshari et al. (2016).

Furthermore, the characteristic function can be directly linked to the cumulants of a specific quadrature measurement statistics. Indeed, is also known as the cumulant-generating function, which implies that

 ∂n∂λnlogχ(λf)∣∣λ=0≡⟨Q(f)n⟩T, (11)

with and . denotes the th cumulant for the measurement of the quadrature . A straightforward calculation Petz (1990); Verbeure (2011) now shows that one can recast the characteristic function in the form

 χ(α)=exp{∞∑n=1inn!⟨Q(α)n⟩T},α∈R2m. (12)

Thus, knowledge of all the cumulants for all the different quadratures, i.e. for all and all orders , implies full knowledge of the quantum state.

We now introduce these cumulants in a more explicit form, by treating them as a special case of truncated correlation functions.

### ii.3 Truncated correlation functions

The cumulants of (11) are related to the statistics of a single mode and do not explicitly elucidate how different modes are correlated. However, to study such questions the cumulant can be generalised to a multimode form that is commonly referred to as the truncated correlation Robinson (1965a) between different quadratures 3.

Truncated correlation functions are the multivariate extensions of cumulants and describe how quadratures for different modes are correlated. They, too, are generated using displacement operators (7). In general, we obtain

 ⟨Q(f1)…Q(fn)⟩T (13) ≡∂logtr(ρD(−2λ1Jf1)…D(−2λnJfn))∂λ1…∂λn∣∣λ1=⋯=λn=0,

which can be related to in (9) through the identity .

In an experimental setting, it is more practical to obtain the truncated correlation functions by jointly measuring distinct quadratures and following a recursive recipe:

 ⟨Q(f1)⟩T=tr{ρQ(f1)} (14) ⟨Q(f1)Q(f2)⟩T=tr{ρQ(f1)Q(f2)}−⟨Q(f1)⟩T⟨Q(f2)⟩T ⟨Q(f1)Q(f2)Q(f3)⟩T=tr{ρQ(f1)Q(f2)Q(f3)} −⟨Q(f1)Q(f2)⟩T⟨Q(f3)⟩T −⟨Q(f1)Q(f3)⟩T⟨Q(f2)⟩T −⟨Q(f2)Q(f3)⟩T⟨Q(f1)⟩T −⟨Q(f1)⟩T⟨Q(f3)⟩T⟨Q(f2)⟩T …\em et cetera.

These truncated correlation functions are experimentally measurable through, for example, multimode homodyne measurement Beck (2000); Armstrong et al. (2012); Cai et al. (2017). By expanding in a specific mode basis in (12), the role of truncated correlations becomes apparent. Thus, the set of truncated correlations forms an important tool for characterisation. More specifically, one may wonder to what order one needs to measure these correlations to extract a given property of the state. Such a property, which is of special interest throughout this text, is the state’s Gaussianity.

### ii.4 Gaussian states

Of particular importance in quantum optics are the Gaussian states. In the broad sense, a state of a CV system is said to be Gaussian if it induces Gaussian statistics in all modes. This implies that the function , and hence also the Wigner function, is a multivariate Gaussian Braunstein and van Loock (2005); Weedbrook et al. (2012):

 χG(α) =exp{−(α,Vα)2+i(ξ,α)}, (15) WG(β) =exp{−12((β−ξ),V−1(β−ξ))}(2π)m√detV (16)

where is a vector which describes the states displacement, and is referred to as the covariance matrix. Therefore is positive semi-definite matrix on , which describes the correlations between different field quadratures in a specific mode basis. A crucial demand for this to be associated with a well-defined quantum state is given by Petz (1990)

 (f1,Vf1)(f2,Vf2)⩾|(f1,Jf2)|2,for all f1,f2∈N(R2m), (17)

which is the multimode version of Heisenberg’s uncertainty relation. Alternatively, the properties of can also be expressed by the condition .

The insertion of (15) in (13) imposes important conditions upon the truncated correlations of Gaussian states. At first it is directly obtained that

 ⟨Q(f)⟩T=(ξ,f), (18) ⟨Q(f1)Q(f2)⟩T=(f1,Vf2)−i(f1,Jf2). (19)

Furthermore, we deduce the general condition that for a Gaussian state

 ⟨Q(f1)…Q(fn)⟩T=0,n>2, (20)

for all . The implication of (20) is that all non-Gaussian states must have non-zero truncated correlations. Furthermore, it was shown Robinson (1965b) that for non-Gaussian states, there is never an order from which onward the truncated correlation functions become zero. Therefore these functions are ideal tools for the operational characterisation of non-Gaussian states. Specifically in multimode systems where full tomographies tend to be completely unfeasible, they are an experimentally accessible alternative.

### ii.5 Entanglement

In the context of quantum physics, one often associates correlations to the study of entanglement which is commonly seen as an important resource for quantum computation and quantum communication. The profound advantage of CV quantum optics, is the simplicity with which Gaussian entanglement between modes can be generated. CV entanglement is strongly dependent of the mode basis in which the problem is described. To see this, it is instructive to consider an arbitrary symplectic basis of and express in (10) in this basis:

 β=m∑i=1ζ(i)xb(i)+ζ(i)pJb(i). (21)

In this mode basis, the Wigner function is a function of the -variables, i.e.  4.

We refer to a CV state as fully separable in the mode basis when its Wigner function can be written as

 W(ζ(1)x,…,ζ(m)x,ζ(1)p,…,ζ(m)p)=∫dλp(λ)m∏i=1Wλ(ζ(i)x,ζ(i)p), (22)

a statistical mixture of a product of single-mode Wigner functions. To obtain a statistical mixture, must correspond to a way of labelling states, and is a probability distribution on this set of labels. Any state for which (22) does not hold is said to be entangled in mode basis . Note that one can introduce more refined terminology for multimode entanglement Gerke et al. (2015). Analysing such different types of multipartite entanglement, however, falls beyond the scope of this work.

It is natural to ask whether there always exists a mode basis in which the state is separable. We will provide a negative answer to this question by showing that this is generally not the case for photon-added and subtracted states. If we can construct a mode basis for which (22) holds, we will refer to the state as passively separable, to highlight that any entanglement present in a specific mode basis can be undone by passive linear optics. States which are not passively separable are now referred to as inherently entangled.

We show now that the Gaussian states of Section II.4 are always passively separable, by using the properties of their covariance matrices. The Wigner function (16) of a non-displaced (i.e. ) Gaussian state is completely governed by the positive semidefinite covariance matrix , which can be decomposed as through the Williamson decomposition. Here is a symplectic matrix and is diagonal (the diagonal elements of are known as the symplectic spectrum). Because is symplectic, it can be further decomposed with the Bloch-Messiah decomposition: , where and are orthogonal and symplectic, and is diagonal and symplectic. This now allows us to rewrite , where is the covariance matrix of a thermal state. We use this structure to separate the covariance into a pure part and added classical noise,

 V=Vs+Vc. (23)

Here, is the covariance matrix of a pure squeezed vacuum state , and is the covariance matrix of the additional noise. Note that, a priori, does not fulfil (17) and is therefore not the covariance matrix of a quantum state. We can think of the state characterised by as being generated by injecting into a noisy Gaussian channel Wolf et al. (2004); Eisert and Wolf (2007). We obtain that

 ρG=∫R2md2mξ′D(ξ′)ρsD(−ξ′)exp{−(ξ,V−1cξ)2}(2π)m√detVc, (24)

which implies that the Wigner function for can be represented by

 WG(β)=∫d2mξWs(β−ξ)pc(ξ), (25)

where

 pc(ξ)=exp{−(ξ,V−1cξ)2}(2π)m√detVc, (26)

and

 Ws(β)=exp{−12(β,V−1sβ)}(2π)m√detVs. (27)

The Bloch-Messiah decomposition naturally gives a specific mode basis (obtained through the orthogonal transformation ) in which factorises for any . In this mode basis, the Wigner function (25) has the form (22), which implies that the state is passively separable. Thus, any entanglement that is present in the original mode basis can be undone by a passive linear optics circuit that is described by , or by measuring quadratures in mode basis associated with .

Thus we provided an explicit construction of a linear optics operation to render a given Gaussian state separable. For mixed states this is typically not the only linear optics operation that can undo entanglement. Indeed, the core ingredient of the decomposition (24) is (23), such that for every pure state covariance matrix , we can set in (23). Because characterises a pure, Gaussian state, we can find a mode basis of symplectic eigenvectors of . We can thus simply diagonalise , where describes an alternative linear optics circuit that can undo entanglement in the Gaussian state. The above method, using Williamson and Bloch Messiah, shows that such a always exists.

For simplicity, we assumed that was non-displaced. However, the argument is straightforwardly extended to displaced Gaussian states by letting the displacement operator act on .

## Iii Single-photon added and subtracted Gaussian states

### iii.1 Induced correlations

As we argued in the introduction of this paper, non-Gaussianity is a crucial ingredient to achieve universal quantum computation. Moreover, for technological applications, it is essential that the complexity of any quantum device can be increased, hence requiring a sense of scalability. In CV quantum optics, we first and foremost consider such scalability in the number of modes. Therefore, we must consider a multimode setup, in which we can incorporate a non-Gaussian operation. From an experimental perspective, a promising procedure to fulfil these conditions is mode-selective photon subtraction Averchenko et al. (2016, 2014); Ra et al. (2017) or addition Zavatta et al. (2007); Parigi et al. (2007).

In this contribution, we limit ourselves to the subtraction or addition of a single photon in a setup with an arbitrary mode number . To effectively model the associated subtraction and addition procedures, we must introduce the annihilation and creation operators for an arbitrary vector in phase space , and , respectively. In our framework, they are defined as

 Missing or unrecognized delimiter for \big (28)

from which we directly obtain an alternative version of the canonical commutation relation (6),

 [a(g1),a†(g2)]=(g1,g2)+i(g1,Jg2). (29)

These operators create or annihilate photons in a specific mode, represented by . However, is a vector in the dimensional phase space, whereas there are only modes. Therefore, we stress that , such that the photons created by the operators and clearly only differ by a global phase. As global phases have no physical importance in quantum physics, the creation operators and really create a photon in the same mode.

We now focus on an arbitrary non-displaced Gaussian state, which we formally describe by a density matrix , and convert it to a non-Gaussian state by means of mode-selective photon-addition or -subtraction in a mode . As was argued in Section II.4, this non-displaced Gaussian state can be completely characterised by its covariance matrix . The new, photon-added and -subtracted states’ density operators are then given by

 ρ−=a(g)ρGa†(g)⟨^n(g)⟩G, (30)

for subtraction, and

 ρ+=a†(g)ρGa(g)⟨^n(g)⟩G+1, (31)

for addition. We introduced the notation for the expectation values in the state , and for the number operator that counts the number of photons in the mode .

In order to characterise the non-Gaussian features of the system, we follow the ideas of Section II.3 and evaluate the truncated correlation functions. If the state is, indeed, non-Gaussian, we should obtain non-zero values for the truncated correlation functions of some order beyond than two. However, because we intend to use the recursive procedure of Section II.3 to evaluate the correlations, it is instructive to start by evaluating the two-point correlation , for arbitrary . Because the state is non-displaced, by definition and we obtain that for the photon-subtracted (hence the superscript “”) state

 ⟨Q(f1)Q(f2)⟩−T =tr{ρ−Q(f1)Q(f2)} (32) =⟨a†(g)Q(f1)Q(f2)a(g)⟩G⟨^n(g)⟩G, (33)

where we used (30) and the cyclic property of the trace. Analogously, for photon-addition we obtain

 ⟨Q(f1)Q(f2)⟩+T=⟨a(g)Q(f1)Q(f2)a†(g)⟩G⟨^n(g)⟩G+1, (34)

The property (20) for non-displaced Gaussian states implies that expectation values of products of quadrature operators factorises in pairs Glauber (1963); Robinson (1965a); Petz (1990); Verbeure (2011). Combining this with the definition (28) for the creation and annihilation operators in terms of quadratures, and with the linearity of the trace, we find

 ⟨Q(f1)Q(f2)⟩−T= ⟨a†(g)Q(f1)⟩G⟨Q(f2)a(g)⟩G⟨^n(g)⟩G +⟨a†(g)Q(f2)⟩G⟨Q(f1)a(g)⟩G⟨^n(g)⟩G (35) +⟨Q(f1)Q(f2)⟩G.

and

 ⟨Q(f1)Q(f2)⟩+T= ⟨a(g)Q(f1)⟩G⟨Q(f2)a†(g)⟩G⟨^n(g)⟩G+1 +⟨a(g)Q(f2)⟩G⟨Q(f1)a†(g)⟩G⟨^n(g)⟩G+1 (36) +⟨Q(f1)Q(f2)⟩G.

To proceed in the evaluation, we use (19) and (28) to obtain

 ⟨a†(g)Q(f)⟩G Missing or unrecognized delimiter for \Big ⟨Q(f)a†(g)⟩G Missing or unrecognized delimiter for \Big ⟨Q(f)a(g)⟩G Missing or unrecognized delimiter for \Big ⟨a(g)Q(f)⟩G =12((f,[V+1]g)+i[(f,[V+1]Jg)]), ⟨^n(g)⟩G =14((g,Vg)+(Jg,VJg)−2). (37)

When we insert these results in (III.1), we ultimately obtain that for the photon-subtracted non-displaced Gaussian state

 ⟨Q(f1)Q(f2)⟩±T=⟨Q(f1)Q(f2)⟩G+(f1,A±gf2), (38) with(f1,A−gf2)≡⟨a†(g)Q(f1)⟩G⟨Q(f2)a(g)⟩G⟨^n(g)⟩G +⟨a†(g)Q(f2)⟩G⟨Q(f1)a(g)⟩G⟨^n(g)⟩G, and(f1,A+gf2)≡⟨a(g)Q(f1)⟩G⟨Q(f2)a†(g)⟩G⟨^n(g)⟩G+1 +⟨a(g)Q(f2)⟩G⟨Q(f1)a†(g)⟩G⟨^n(g)⟩G+1,

is a matrix which acts on the space . Inserting (37) in (38) directly leads to

 A±g=2(V±1)(Pg+PJg)(V±1)tr{(V±1)(Pg+PJg)}. (39)

Here we introduced and as the projectors on the vectors and , respectively. This implies that is the projector on the two-dimensional phase space, associated with the mode in which the photon was subtracted. It can directly be verified that describes additional correlations between quadratures that are induced by the photon-subtraction or -addition process. Ultimately, these additional correlations are completely determined by the mode from which the photon is subtracted, and the correlation matrix of the initial non-displaced Gaussian state .

Experimentally, however, it is hard to guarantee that (30) and (31) are the exact states which we obtain. In general, the subtraction process adds some degree of impurity. There are various sources of impurities in an experimental context, ranging from photon-losses to contributions of higher photon-numbers in the subtraction Ra et al. (2017), which go beyond the scope of the present work. Nevertheless, we consider one important type of impurity in the subtraction process, related to lack of control of the mode-selectivity Averchenko et al. (2016). In the most extreme case, one may think of photon subtraction by means of a beamsplitter, where it is impossible to infer from which mode the photon originated in the case of co-propagating modes. In general terms, it is hard to control exactly in which mode the photon is added or subtracted Averchenko et al. (2016). This implies that we have to deal with a mixture of the form

 ρ−=∑kγka(gk)ρGa†(gk)∑kγk⟨^n(gk)⟩G,with∑kγk=1, and γk⩾0, (40)

for subtraction, or

 ρ+=∑kγka†(gk)ρGa(gk)1+∑kγk⟨^n(gk)⟩G,with∑kγk=1, and γk⩾0, (41)

for addition. The details of the participating modes and the depend strongly on the experimental setup and can be estimated through a detailed modelling Averchenko et al. (2016, 2014).

Through the linearity of the expectation value, we can directly verify that

 ⟨Q(f1)Q(f2)⟩±T=⟨Q(f1)Q(f2)⟩G+(f1,A±mixf2), (42)

where

 A±mix=2(V±1)∑kγk(Pgk+PJgk)tr{(V±1)∑kγk(Pgk+PJgk)}(V±1). (43)

We use (42) to evaluate the higher-order truncated correlations in Appendix A. This leads to the remarkable result that these truncated correlations, too, are governed by the matrix . As a final result, we obtain

 Unknown environment '%' (44)

for all . indicates the set of all pair-partitions, i.e. all the ways of dividing the set up in pairs. In literature, e.g. Mansour (2013), this partition is also known as a perfect matching. For even orders, the truncated correlations (44) are generally non-zero. This is a clear statistical signature of the non-Gaussian character of the state.

### iii.2 Phase space representations

#### Multimode Wigner function

To highlight that the above truncated correlation functions grant us full knowledge of the quantum state, we use them to construct the quantum characteristic function for the photon-subtracted state (30) by virtue of (12). To do so, we need to know the state’s cumulants, which are the truncated correlations (44) for . Central in this evaluation is that every pair-partition in (44) contributes the same term, , because contribution of each mode is the same. We only need to count the number of pair-partitions to know which combinatorial factor to add. We find that the cumulant is given by

 ⟨Q(f)2k⟩T=(−1)k−1(2k−1)!2k−1(f,A±mixf)k (45) +(f,Vf)δk,1, ⟨Q(f)2k−1⟩T=0. (46)

It now remains to evaluate the series in Eq. (12), for which we find that

 ∞∑n=1inλnn!⟨Q(f)n⟩T=−λ2(f,Vf)2+∞∑k=1i2kλ2k(2k)!(−1)k−1(2k−1)!2k−1(f,A±mixf)k=−λ2(f,Vf)2−∞∑k=11k(λ2(f,A±mixf)2)k, (47)

where the series was already rewritten to only sum over the even cumulants since all odd contributions are zero.The final series in (47) is subtle because it does not necessarily converge. However, maps the points of phase space to the complex plane, as such we may resort to an analytical continuation of the series to obtain that

 χ(λf) =exp{−λ22(f,Vf)−∞∑k=11k(λ2(f,A±mixf)2)k} =(1−λ2(f,A±mixf)2)exp{−λ22(f,Vf)}. (48)

Therefore we have obtained the quantum characteristic function (9) to fully characterise the states. In principle, this also allows us to derive the Wigner function by means of a multi-dimensional Fourier transformation. We may formally write the Wigner function as

 Missing or unrecognized delimiter for \Bigg (49)

This Fourier transform is explicitly computed in Appendix B, and leads to

 W±(β)=12((β,V−1A±mixV−1β)−tr(V−1A±mix)+2)WG(β), (50)

where is the Wigner function of the initial Gaussian state (16) before the addition or subtraction of the photon. This now gives us the full, multimode, Wigner functions of a non-displaced photon-added or -subtracted state. We observe that the general structure of the Wigner function is given by a multivariate polynomial of order two, multiplied by the Gaussian Wigner function of the initial state.

#### Negativity

The negativity of the Wigner function is often seen as a genuine quantum feature in CV systems. With (50) we have all the tools at hand to analyse such features in the Wigner function.

At first, we note that the Wigner function (50) is negative if and only if there are vectors for which

 (β,V−1A±mixV−1β)−tr(V−1A±mix)+2⩽0. (51)

However, it is directly verified that is a positive-semidefinite matrix, hence

 (β,V−1A±mixV−1β)⩾0, for all β∈R2m. (52)

Therefore, the necessary and sufficient condition for the existence of negative values of the Wigner function is

 tr(V−1A±mix)⩾2. (53)

By setting in (51) we clearly see that (53) is, indeed, a sufficient condition. Through (43), we can rephrase this condition as

 ∑kγk[(gk,V−1gk)+(Jgk,V−1Jgk)]>2,for subtraction,∑kγk[(gk,V−1gk)+(Jgk,V−1Jgk)]>−2,for addition, (54)

which is automatically fulfilled for photon addition. Hence, we formally show that photon addition to a non-displaced Gaussian state always induces a negative Wigner function, even when the initial state and the addition process are mixed. On the other hand, for photon subtraction the condition for negativity of the Wigner function can be violated when there is too much thermal noise compared to the amount of squeezing (see, e.g., the example in Section IV.1).

Finally, we emphasise that the equation

 (β,V−1A±mixV−1β)=tr(V−1A±mix)−2 (55)

defines the manifold of zeros of the Wigner function. Specifically equation (55) generates a multidimensional ellipsoid. The details of the manifold depend strongly on the details of the subtraction or addition process, and on the covariance matrix . However, as expected, the general condition for equation (55) to have solutions is also given by (53).

#### Entanglement

In this section we elaborate on the passive separability of the Wigner function (50). First, we prove that, whenever a photon is added or subtracted to or from a mode which is not entangled to any other modes in the initial Gaussian state, the resulting photon-added or -subtracted state will remain passively separable. We then prove for pure states, that subtraction or addition of a photon in any other mode renders the state inherently entangled.

For any possible decomposition (23) of the Gaussian state’s covariance matrix , we may use (24) to write the photon-subtracted state as

 ρ−=a(g)ρGa†(g)⟨^n(g)⟩G=1⟨^n(g)⟩G∫d2mξa(g)D(ξ)ρsD(−ξ)a†(g)pc(ξ), (56)

and the photon-added state as

 ρ+=a†(g)ρGa(g)⟨^n(g)⟩G+1=1⟨^n(g)⟩G+1∫d2mξa†(g)D(ξ)ρsD(−ξ)a(g)pc(ξ), (57)

where we initially focus on the pure subtraction of a photon from mode . The evaluation of is cumbersome and is therefore left for Appendix C, where we describe the general subtraction/addition of a photon from/to a displaced state. In general, we can use (127) to write the Wigner function of in (56) and (57) as

 W−(β)=∫d2mξW−ξ(β)×⟨^n(g)⟩s+14[(ξ,g)2+(ξ,Jg)2]⟨^n(g)⟩Gpc(ξ), (58)

and

 Unknown environment '%' (59)

with

 W±ξ(β)= Ws(β−ξ)tr((Vs+∥ξ∥2Pξ±1)(Pg+PJg)) (60) ×(∥(Pg+PJg)(1±V−1s)(β−ξ)∥2 +2(ξ,(Pg+PJg)(1±V−1s)(β−ξ)) +tr((Pg+PJg)(∥ξ∥2Pξ−V−1s∓1))).

Note that and denote the Wigner function and covariance matrix, respectively, of , as introduced in (24, 25). The passive separability of the Wigner functions (58) and (59) now depends on two aspects. Firstly, it hinges on the factorisability of the pure state Wigner function , as given by (60). Secondly, we require that

 p−c(ξ)≡⟨^n(g)⟩s+14[(ξ,g)2+(ξ,Jg)2]⟨^n(g)⟩Gpc(ξ), (61)

and

 p+c(ξ)≡⟨^n(g)⟩s+1+14[(ξ,g)2+(ξ,Jg)2]⟨^n(g)⟩G+1pc(ξ), (62)

are well-defined probability distributions.

Probability distributions It is straightforwardly verified that are well-defined probability distributions. Because we know that and are positive functions, it suffices to validate their normalisation. To do so, we evaluate

 ∫d2mξ[(ξ,g)2+(ξ,Jg)2]pc(ξ)=∫d2mξ[(ξ,g)2+(ξ,Jg)2]exp{−(ξ,V−1cξ)2}(2π)m√detVc=(g,Vcg)+(Jg,VcJg)=tr{(Pg+PJg)Vc}, (63)

where we used that

 ∫d2mξPξ∥ξ∥2exp{−(ξ,V−1cξ)2}(2π)m√detVc=Vc. (64)

This implies that

 ∫d2mξp−c(ξ)=⟨^n(g)⟩s+14tr{(Pg+PJg)Vc}⟨^n(g)