Conditional non-Gaussian quantum state preparation

We develop a general formalism, based on the Wigner function representation of continuous-variable quantum states, to describe the action of an arbitrary conditional operation on a multimode Gaussian state. We apply this formalism to several examples, thus showing its potential as an elegant analytical tool for simulating quantum optics experiments. Furthermore, we also use it to prove that EPR steering is a necessary requirement to remotely prepare a Wigner-negative state.


I. INTRODUCTION
In continuous-variable (CV) quantum physics, Gaussian states have long been a fruitful topic of research [1][2][3][4][5][6][7][8][9][10]. They appear naturally as the ground states of systems of many noninteracting particles in the form of thermal states [11], or as the coherent states that describe the light emitted by a laser [3]. Through nonlinear processes, it is possible to reduce the noise beyond the shot noise limit (at the price of increased noise in a complementary observable), and create squeezed states [12][13][14][15][16][17]. For the purpose of metrology, such squeezed states are often enough to obtain a significant boost in performance [18][19][20][21].
On theoretical grounds, Gaussian states are relatively easy to handle [8,9]. The quantum statistics of the continuousvariable observables (e.g., the quadratures in quantum optics) are described by Gaussian Wigner functions. All interesting quantum features can be deduced from the covariance matrix that characterises this Gaussian distribution on phase space. Hence, whenever the number of modes remains finite, the techniques of symplectic matrix analysis are sufficient to study Gaussian quantum states. This has generated an extensive understanding of the entanglement properties of Gaussian states [22][23][24][25][26][27], and recently it has also lead to the development of a measure for quantum steering (see [28]) of Gaussian states with Gaussian measurements [29][30][31][32], which we refer to as Einstein-Podolsky-Rosen (EPR) steering.
Even though they have many advantages, Gaussian states are of limited use to quantum technologies beyond sensing. They have been shown to be easily simulated on classical devices [33], and in particular Wigner negativity is known to be a necessary resource for reaching a quantum computation advantage [34]. However, it should be stressed that recent work has found large classes of Wigner negative states that can also be simulated easily [35]. In other words, Wigner negativity is necessary but not sufficient to reach a quantum computation advantage [36].
In the particular case of CV quantum computation, Gaussian states play an essential role in the measurement-based approach [37]. In this paradigm, one establishes large Gaussian entangled states, known as cluster states, which form the backbone of the desired quantum routine [38]. Several recent * mattia.walschaers@lkb.upmc.fr breakthroughs have led to the experimental realisation of such states [39][40][41][42][43]. Nevertheless, to execute quantum algorithms that cannot be simulated efficiently, one must induce Wigner negativity. In the spirit of measurement-based quantum computation, this feature is induced by measuring non-Gaussian observables, e.g., the number of photons, on a subset of modes [44][45][46][47]. Such a measurement then projects the remainder of the system into a non-Gaussian state. The exact properties of the resulting state depend strongly on the result of the measurement.
Such a conditional preparation of non-Gaussian quantum states is common procedure in quantum optics experiments [48]. Basic examples include the heralding of single-photon Fock states after parametric down-conversion [49][50][51], photon addition and subtraction [52][53][54][55][56][57], and known schemes to prepare more exotic states such as Schrödinger-cat [58,59] or Gottesman-Kitaev-Preskill states [60]. Remarkably, though, a practical framework to describe the effect of such conditional operations on arbitrary Gaussian states is still lacking. Notable exceptions where one does study arbitrary initial states usually rely on specific choices for the conditional measurement.
Here, in Section III, we introduce a practical framework to describe the resulting Wigner function for a quantum state that is conditionally prepared by measuring a subset of modes of a Gaussian multimode state. The techniques used in this work are largely based on classical multivariate probability theory and provide a conceptually new understanding of these conditioned states. In Section IV, we unveil the most striking consequence of this new framework: we can formally prove that EPR steering in the initial Gaussian state is a necessary requirement for the conditional preparation of Wigner-negative states, regardless of the measurement upon which we condition. This solidifies a previously conjectured general connection between EPR steering and Wigner-negativity. As shown in Section V, our framework reproduce a range of known state-preparation schemes and can be used to treat more advanced scenarios, which could thus far not be addressed by other analytical methods. First, however, we review the phase space description of multimode CV systems in Section II.

II. PHASE SPACE DESCRIPTION OF MULTIMODE CONTINUOUS-VARIABLE SYSTEMS
The CV approach studies quantum systems with an infinitedimensional Hilbert space H based on observables,x andp that have a continuous spectrum and obey the canonical commutation relation [x,p] = 2i (the factors two is chosen to normalise the vacuum noise to one). Common examples include the position and momentum operators in mechanical systems, or the amplitude and phase quadratures in quantum optics. In this work, we will use quantum optics terminology, but the results equally apply to any other system that is described by the algebra of canonical commutation relations (i.e., any bosonic system).
In a single-mode system, the quadrature observablesx and p determine the optical phase space. The latter is a twodimensional real space, where the axes denote the possible measurement outcomes forx andp. It is common practice to represent a given stateρ by means of its measurement statistics forx andp on this optical phase space, as in statistical physics. However, becausex andp are complementary observables, they cannot be measured simultaneously, and thus, a priori, we cannot construct a joint probability distribution of phase space that reproduces the correct marginals to describe the measurement statistics the quadratures. Therefore, the phase space representation of quantum states are quasiprobability distributions. The quasi-probability distribution that reproduces the measurement statistics of the quadrature observables as its marginals, is known as the Wigner function [61][62][63] W(x, p) = 1 (2π) 2 For some quantum states, this function has the peculiar property of reaching negative values. This Wigner-negativity is a genuine hallmark of quantum physics, and it is understood to be crucial in reaching a quantum computational advantage.
Here, we will consider a multimode system comprising m modes. Every mode comes with its own infinite-dimensional Hilbert space, associated to a two-dimensional phase space, and observablesx j andp j . The total optical phase space is, thus, a real space R 2m with a symplectic structure Ω = m ω, where the two-dimensional matrix ω is given by Therefore, Ω has the properties Ω 2 = −1 and Ω T = −Ω. Any normalised vector f ∈ R 2m defines a single optical mode with an associated phase space span{ f , Ω f } (i.e., when f generates the phase space axis associated with the amplitude quadrature of this mode, Ω f generates the axis for the associated phase quadrature). Henceforth, we will refer to the subsystem associated with the phase space span{ f , Ω f } as "the mode f ". Every point α ∈ R 2m can also associated with a generalised quadrature observablê These observables satisfy the general canonical commutation relation [q( α),q( β)] = −i α T Ω β. Physically, such observablê q( α) can be measured with a homodyne detector by selecting the mode that is determined by the direction of α, and multiplying the detector outcome by α . In our theoretical treatment, such generalised quadratures are useful to define the quantum characteristic function of any multimode stateρ for an arbitrary point α in phase space. The multimode Wigner function of the state is then obtained as the Fourier transform of the characteristic function where x ∈ R 2m can, again, be any point in the multimode phase space, and the coordinates of x represent possible measurement outcomes forx j andp j . The Wigner function can be used to represent and characterise an arbitrary quantum state of the multimode system. In the same spirit, we can also define the phase space representation of an arbitrary observableÂ as such that we can fully describe the measurement statistics of an arbitrary quantum observable on phase space, by invoking the identity to evaluate expectation values. In practice, it is often challenging to obtain Wigner functions for arbitrary states or observables, but in some cases they can take convenient forms. A particular class of convenient states are Gaussian states, where the Wigner function W( x) is a Gaussian. As a consequence, the Wigner function is positive, and can thus be interpreted as a probability distribution. This Gaussian distribution is completely determined by a covariance matrix V, and mean-field ξ, such that the Wigner function takes the form This forms the basis of our preparation procedure for non-Gaussian states as we assume that our initial multimode system is prepared in such a Gaussian state.
To perform the conditional state-preparation, we divide the m-mode system in two subsets of orthogonal modes, f = { f 1 , . . . , f l } and g = {g 1 , . . . , g l }, with l + l = m, and perform a measurement on the modes in g. We can then describe EPR Steering Conditioning Figure 1. Sketch of the conditional state preparation scenario: a multimode quantum state with density matrixρ is separated over two subsets of modes, f and g. A measurement is performed on the modes in g, yielding a result associated with a POVM elementÂ. Conditioning on this measurement outcome "projects" the the subset of mode f into a stateρ f|Â . The directional EPR steering, discussed in Section IV, is highlighted.
the subsystems of modes f and g by phase spaces R 2l and R 2l , respectively. As such, the joint phase space can be mathematically decomposed as R 2m = R 2l ⊕ R 2l . A general point x in the multimode phase space R 2m can thus be decomposed as x = x f ⊕ x g , where x f and x g describe the phase space coordinates associated with the sets of modes f and g, respectively.
In particular, x f can be expanded in a particular modes basis f 1 , . . . , f l as x f = (x f 1 , p f 1 , . . . , x f l , p f l ), where the coordinates x f 1 and p f 1 are obtained as A completely analogous treatment is possible for the coordinates associated with the set of modes g.

III. CONDITIONAL OPERATIONS IN PHASE SPACE
In quantum optics, we associate a Hilbert space (more precisely a Fock space) to each of these modes. The Hilbert space H of the entire system can then be structured as H = H f ⊗H g , where H f ( H g ) describes that quantum states of the set of orthogonal modes f (g). Formally, the state of our full m-mode system is then described by a density matrixρ that acts on H.
Within this manuscript, we perform a conditional operation in the set of modes g, which we describe through a (not necessarily normalised) set of Kraus operators [64]X j that act on Such a conditional operation naturally arises as a postmeasurement state, whenX j is a projector, or whenÂ = jX † jX j is a more general POVM element. The positive semidefinite operatorÂ is useful to express the reduced state of the set of modes f:ρ where tr g denotes the partial trace of the Hilbert space H g associated with the set of mode g. Our general goal is to understand the properties of the stateρ f|Â .
As we are interested in the Wigner function for the state of the subset of modes f, we translate (12) to its phase-space representation. We initialize the total system in a Gaussian state with Wigner function W( x). Subsequently we also define the Wigner function WÂ( x g ) of the positive operatorÂ, which is a function that is defined according to (6) on the phase space that describes the subset of modes g. As such, we find that BecauseÂ is a positive semi-definite operator, the denominator is a positive constant.
As presented in (13), the Wigner function W f|Â ( x f ) is impractical to use and its properties are not apparent. Hence, we now introduce some mathematical tools to obtain a more insightful expression for W f|Â ( x f ). First, we use that, for Gaussian states, W( x) is a probability distribution on phase space, such that we can define the conditional probability distribution through where W f ( x f ) is the reduced Gaussian state for the set of modes f, Because W( x) is a Gaussian probability distribution, the conditional probability distribution W( x g | x f ) is also a Gaussian distribution [66] with covariance matrix where V g and V f are the covariance matrices describing the subsets of modes g and f in the initial state, whereas V gf describes all the initial Gaussian correlations between those subsets. Note that this covariance matrix is the same far all points x f ∈ R 2l , which is a particular property of Gaussian conditional probability distributions. Furthermore, the distribution W( x g | x f ) also contains a displacement where ξ g and ξ f describe the displacements of the initial state in the sets of modes g and f, respectively. Generally, the phase space probability distribution W( x g | x f ) is not a valid Wigner function of a well-defined quantum state, in the sense that it would violate the Heisenberg inequality. However, it does remain a well-defined probability distribution, i.e., it is normalised and positive. Thus, it still has interesting properties that we can exploit to formulate a general expression for W f|Â ( x f ). Let us therefore define which is the expectation value of the phase-space representation ofÂ with respect to the probability distribution W( x g | x f ). We can then use (14) and (18) to recast (13) in the following form: where we introduce the notation The major advantage of this formulation is that Â g| x f represents the average with respect to a Gaussian probability distribution, such that one can use several computational techniques that are well-know for Gaussian integrals. A notable property is the factorisation of higher moments in multivariate Gaussian distributions, such that Â g| x f can generally be expressed algebraically in terms of the components of V g| x f and ξ g| x f .
Finally, we remark that Â g| x f = Â in absence of correlations between the set of modes g that are conditioned upon and the set modes f for which we construct the reduced state. This result is directly responsible for the previously obtained results related to the spread of non-Gaussian features in cluster states [67].

IV. EINSTEIN-PODOLSKY-ROSEN STEERING AND WIGNER-NEGATIVITY
When two systems are connected through a quantum correlation, on can, in some cases, perform quantum steering [28]. Colloquially, we say that a subsystem X can steer a subsystem Y when measurements of certain observables in X can influence the conditional measurement statistics of observables in Y beyond what is possible with classical correlations. Ultimately, in quantum steering one studies properties of conditional quantum states as compared to a local hidden variable model for any observables X and Y, acting on X and Y, respectively. Contrary to the case of Bell non-locality, quantum steering considers an asymmetric local hidden variable model: (21) where one assumes that the probability distributions P Q (Y = y | λ) of steered party Y follow the laws of quantum mechanics. For the party X, which performs the steering, no such assumption is made and any probability distribution is allowed. Such local hidden variable model can typically be falsified, either by brute force computational methods [68] or via witnesses [69]. These methods have been applied in a variety of contexts to experimentally observe quantum steering [31,32,[70][71][72][73][74][75][76].
A paradigmatic example is found when performing homodyne measurements on EPR state [77]: when the entanglement in the system is sufficiently strong, one can condition thê x andp quadrature measurements in Y on the outcome of the same quadrature measurement in X. The obtained conditional probability distributions for the quadrature measurements in Y can violate the Heisenberg inequality, even when averaged over all measurement outcomes in X. The violations of such a conditional inequality is impossible with classical correlations, but is a hallmark of quantum steering.
Quantum steering can occur in all types of quantum states, with all kinds of measurements. In CV, one often refers to the particular case of Gaussian states that can be steered through Gaussian measurements as EPR steering. Recently, other forms of steering for Gaussian states have been developed under the name of non-classical steering [78]. In this approach, one checks whether Gaussian measurements in X can induce a nonclassical conditional state in Y. Throughout our current work, the focus lies on EPR steering, where the systems X and Y are the sets of modes f and g, respectively.
In previous work, we showed EPR steering is a necessary prerequisite to remotely generate Wigner negativity through photon subtraction [79]. More precisely, when a photon is subtracted in a mode g, the reduced state Wigner function of a correlated mode f can only be non-positive if mode f is able to steer mode g. When one allows for an additional Gaussian transformation on mode g prior to photon subtraction, we found that EPR steering from f to g is also a sufficient condition to reach Wigner negativity in mode f .
The formalism that was developed in the previous section allows us to generalize this previous result to arbitrary conditional operations on an arbitrary number of modes: Theorem 1. For any initial Gaussian stateρ and any conditional operationÂ in (12), EPR steering between the set of modes f and the set of modes g is necessary to induce Wigner negativity in W f|Â ( x f ).
Proof. Gaussian EPR steering is generally quantified through the properties of V g| x f . In particular, one can show that the set of modes in f can jointly steer the set of modes g if and only if V g| x f violates the Heisenberg inequality [29,30]. The crucial consequence is that W( x g | x f ), as defined in (14), is itself a well-defined Gaussian quantum state when the modes in f cannot steer the modes g. For all possible x f , we can thus associate this Gaussian quantum state with a density matrix The crucial observation is that, for any x f , is the expectation value ofÂ in a well-defined quantum state ρ g| x f . BecauseÂ is a positive semi-definite operator, we directly find that Therefore, the overall conditional Wigner function W f|Â ( x f ) in (19) is non-negative. We can only achieve Â g| x f < 0 for certain points x f ∈ R 2l when V g| x f violates the Heisenberg inequality. This concludes that in absence of EPR steering Note that the steps in this proof rely heavily on the fact that the initial state is Gaussian. For other types of quantum states, we cannot directly relate quantum steering to the properties of W( x g | x f ).

A. Heralding
In the first example, we consider a scenario where a photonnumber revolving measurement is performed on one of the output modes, which can be considered a special case of the situation considered in [47]. Heralding is ubiquitous in quantum optics, as it is one of the most common tools to generate single photon Fock states [49][50][51].
To study heralding, we use (19) where a measurement of the number of photon n in a single mode g is performed. We assume that this measurement is optimal, and, thus, that we project on a Fock state |n . In this case, we setÂ = |n n|, and therefore we obtain that where we used the closed form of the Laguerre polynomial. Hence, we can now use this expression to calculate |n n| g| x f . It is convenient to explicitly write and we can recast After a substitution in the integral, we then find that where we defined σ = V g| x f (1 + V g| x f ), which is now the covariance matrix of a new Gaussian probability distribution. The final expression is then determined by the moments of the Gaussian distribution with covariance matrix σ and displacement ξ g| x f . Even though this expression is relatively elegant, it can be remarkably tedious to compute for larger values of n.
First, let us focus on the experimentally relevant case where n = 1 as an illustration. The evaluation of (26) is than conducted by calculating the second moments of a Gaussian distribution, such that we ultimately find where we set ξ g = 0, thus assuming that there is no mean field in mode g. We note that this function reaches negative values if and only of tr[(1 + V g| x f ) −1 V g| x f ] < 1. Using Williamson's decomposition as we did in [79], it can be shown that this condition can only be fulfilled when the set of modes f can perform EPR steering in mode g, or, in other words, when V g| x f violates the Heisenberg inequality. This is exactly what we can expect from our general result in Section IV.
In general we know that the Wigner function (19) can only be negative when V g| x f is not a covariance matrix of a welldefined quantum state. However, determining the existence of zeroes of this Wigner function is a cumbersome task. For heralding with n > 1 we, therefore, restrict to numerical simulations using a specific initial state.
This specific initial state is generated by mixing two squeezed thermal states on a balance beamsplitter, where one of the output modes will serve as f , and the other as g. In the limiting case where the initial thermal noise vanishes, we recover the well-known EPR state which manifests perfect photon-number correlations between modes f and g. In this case, it is clear a detection of n photons in mode g will herald the state |n in mode f . However, by introducing thermal noise the photon-number correlations fade and the properties of the heralded state in mode f are less clear. Thermal noise will also gradually reduce the EPR steering in the system, such that the Wigner negativity in mode f will vanish when the thermal noise becomes too strong. Hence, with this example we can study the interplay between Wigner negativity and EPR steering in a controlled setting.
The squeezed thermal state is characterised by a covariance matrix V = diag[δ/s, δs], where δ denotes the amount of initial thermal noise, and s is the squeezing parameter. We initially start with two copies of such a state, and rotate the phase of one of them by π/2 (see Fig. 2). When both modes are mixed on a beamsplitter, the resulting state manifests EPR steering depending on parameters δ and s, which can be quantified through [30] where we explicitly use the fact that V g| x f is a two-dimensional matrix. When we then post-select on the number of photons, n, measured in one output mode, we herald a conditional non-Gaussian state in the other mode. In Fig. 2, we show the resulting Wigner functions for the case where the detected number of photons is n = 5. When the amount of EPR steering is varied (note that µ = 0.55 corresponds to the pure state), we see that the resulting Wigner function rapidly loses Wigner negativity. In full agreement with our general result of the previous section, we also find that the Wigner negativity vanishes when there is no EPR steering. A more quantitative study of the Wigner negativity can be found in Panel (c) of Fig. 2, where we vary, both, the amount of steering µ and the number of detected photons n. The Wigner negativity is measured by the quantity [80][81][82] When the state is pure (here for µ = 0.55), a detection of n photons in one mode herald a Fock state |n in the other mode and the Wigner negativity, thus, increases with n. However, once the state is no longer pure and the steering decreases, we observe the existence of an optimal value n for which the maximal amount of Wigner negativity is obtained. For very weak EPR steering (e.g. µ = 0.08 in this calculation), this optimal value is obtained for n = 1. This numerical study shows the fruitfulness of our presented framework to study a very concrete heralding scheme. Furthermore, the example confirms the relationship between Gaussian EPR steering and Wigner negativity.

B. Photon-added and -subtracted states
Ideal photon addition and subtraction are defined by acting with a creation operatorâ † or annihilation operatorâ, respectively, on the quantum state. In practice, these operations are often realised by using some form of heralding [52], which we treated in the previous example. However, it tends to be more convenient to use the idealised model, based on creation and annihilation operators, and it has been shown experimentally that this model is highly accurate. This model also fits the conditional state framework of (11), where we setX j to be a creation or annihilation operator.
In multimode systems, photon addition and subtraction have been considered for their entanglement properties, which sprouted a range of theoretical [83][84][85][86][87][88][89][90] and experimental [91][92][93] results. Many of the obtained theoretical results rely on the purity of the initial Gaussian state, and are hard to generalise to arbitrary Gaussian states. In recent years, there has been some progress in developing analytical tools to describe general photon subtracted states [90,94], but it remain challenging to use these techniques to evaluate entanglement measures. Therefore, one has investigated related questions, such as, for example, the spread of non-Gaussian features in multimode systems [67,79,95,96].
The framework presented in this manuscript is particularly fruitful to investigate the spread of non-Gaussian features through photon addition or subtraction. We will first show how the results of [79] can be recovered via (19). Then, we use the present framework to provide analytical results for the states that can be obtained by subtracting multiple photons in a multimode system.

Adding or subtracting a single photon
We will start by studying the addition and subtract of a single photon. The scenario for photon-subtracted states was studied in detail in [79] and out goal in this example is to show how these previous results can be obtain in the context of our present framework. Furthermore, we also study photon addition, which not yet been considered in the context of the remote generation of Wigner negativity.
Creating and annihilation operators are by construction operators that act on a single mode g. In the single photon scenario, we find the photon-subtracted statê and the photon-added statê These states are clearly fit the framework of (11). In the context of (12), the reduced state of the set of modes f is obtained by choosingÂ =n g andÂ =n g + 1 for photon subtraction and addition, respectively. We can then use (19) to obtain the Wigner function in the subset of modes f, for which we must evaluate n g g| x f . To this goal, we evaluate the Wigner function of the number operator, which is given by such that we directly find that where the dependence on x f comes from ξ g| x f . Thus, we find for the photon-subtracted state that

< l a t e x i t s h a 1 _ b a s e 6 4 = " Q X Y C H 6 k 7 b v a p e h I W z f Z i X X I 1 T 8 8 = " > A A A B 8 X i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B 0 7 I r F v U g F L 1 4 r G A / s F 1 K N s 2 2 o U l 2 S b J C W f o v v H h Q x K v / x p v / x r T d g 7 Y + G H i 8 N 8 P M v D D h T B v P + 3 Y K K 6 t r 6 x v F z d L W 9 s 7 u X n n / o K n j V B H a I D G P V T v E m n I m a c M w w 2 k 7 U R S L k N N W O L q d + q 0 n q j S L 5 Y M Z J z Q Q e C B Z x A g 2 V n r s i h R d I 8 + t V n v l i u d 6 M 6 B l 4 u e k A j n q v f J X t x + T V F B p C M d a d 3 w v M U G G l W G E 0 0 m p m 2 q a Y D L C A 9 q x V G J B d Z D N L p 6 g E 6 v 0 U R Q r W 9 K g m f p 7 I s N C 6 7 E I b a f A Z q g X v a n 4 n 9 d J T X Q Z Z E w m q a G S z B d F K U c m R t P 3 U Z 8 p S g w f W 4 K J Y v Z W R I Z Y Y W J s S C U b g r / 4 8 j J p n r n + u X t 1 f 1 6 p 3 e R x F O E I j u E U f L i A G t x B H R p A Q M I z v M K b o 5 0 X 5 9 3 5 m L c W n H z m E P 7 A + f w B X M W P b Q = = < / l a t e x i t >
µ = 0.37

< l a t e x i t s h a 1 _ b a s e 6 4 = " P 8 r 5 8 / g m a T S B U G a h F W G / e X U P D W s = " > A A A B 8 X i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B 0 7 K r h e p B K H r x W M F + Y L u U b J p t Q 5 P s k m S F s v R f e P G g i F f / j T f / j W m 7 B 6 0 + G H i 8 N 8 P M v D D h T B v P + 3 I K K 6 t r 6 x v F z d L W 9 s 7 u X n n / o K X j V B H a J D G P V S f E m n I m a d M w w 2 k n U R S L k N N 2 O L 6 Z + e 1 H q j S L 5 b 2 Z J D Q Q e C h Z x A g 2 V n r o i R R d I c 8 9 r / X L F c / 1 5 k B / i Z + T C u R o 9 M u f v U F M U k G l I R x r 3 f W 9 x A Q Z V o Y R T q e l X q p p g s k Y D 2 n X U o k F 1 U E 2 v 3 i K T q w y Q F G s b E m D 5 u r P i Q w L r S c i t J 0 C m 5 F e 9 m b i f 1 4 3 N d F F k D G Z p I Z K s l g U p R y Z G M 3 e R w O m K D F 8 Y g k m i t l b E R l h h Y m x I Z V s C P 7 y y 3 9 J 6 8 z 1 q + 7 l X b V S v 8 7 j K M I R H M M p + F C D O t x C A 5 p A Q M I T v M C r o 5 1 n 5 8 1 5 X 7 Q W n H z m E H 7 B + f g G X M O P b Q = = < / l a t e x i t >
µ = 0.22

< l a t e x i t s h a 1 _ b a s e 6 4 = " R o x q 2 2 i P 5 b o W d 0 H M D c 9 d l 2 N t A U U = " > A A A B 8 X i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B 0 7 J b C u p B K H r x W M F + Y L u U b J p t Q 5 P s k m S F s v R f e P G g i F f / j T f / j W m 7 B 2 1 9 M P B 4 b 4 a Z e W H C m T a e 9 + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U U v H q S K 0 S W I e q 0 6 I N e V M 0 q Z h h t N O o i g W I a f t c H w 7 8 9 t P V G k W y w c z S W g g 8 F C y i B F s r P T Y E y m 6 R p 5 b r f b L F c / 1 5 k C r x M 9 J B X I 0 + u W v 3 i A m q a D S E I 6 1 7 v p e Y o I M K 8 M I p 9 N S L 9 U 0 w W S M h 7 R r q c S C 6 i C b X z x F Z 1 Y Z o C h W t q R B c / X 3 R I a F 1 h M R 2 k 6 B z U g v e z P x P 6 + b m u g y y J h M U k M l W S y K U o 5 M j G b v o w F T l B g + s Q Q T x e y t i I y w w s T Y k E o 2 B H / 5 5 V X S q r p + z b 2 6 r 1 X q N 3 k c R T i B U z g H H y 6 g D n f Q g C Y Q k P A M r / D m a O
f F e X c + F q 0 F J 5 8 5 h j 9 w P n 8 A U 6 q P Z w = = < / l a t e x i t > µ = 0.08  From this result, we immediately observe that the potential Wigner-negativity of these states depends on whether or not trV g| x f < 2. In [79] it was shown through the Williamson decomposition that trV g| x f 2 det V g| x f . This directly implies that EPR steering (28) is a necessary condition to reach Wigner-negativity. It is instructive to emphasise that

< l a t e x i t s h a 1 _ b a s e 6 4 = " R F N A S X U P u K c C C c / 8 x Q + w x 0 f e e g c = " > A A A B 8 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s x I w b o Q i m 5 c V r A P b I e S S T N t a J I Z k o x Q h v 6 F G x e K u P V v 3 P k 3 p t N Z a O u B k M M 5 9 3 L v P U H M m T a u + + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U V t H i S K 0 R S I e q W 6 A N e V M 0 p Z h h t N u r C g W A a e d Y H I 7 9 z t P V G k W y Q c z j a k v 8 E i y k B F s r P T Y F w m 6 R m 7 V r Q / K F f t l Q K v E y 0 k F c j Q H 5 a / + M C K J o N I Q j r X u e W 5 s / B Q r w w i n s 1 I / 0 T T G Z I J H t G e p x I J q P 8 0 2 n q E z q w x R G C n 7 p E G Z + r s j x U L r q Q h s p c B m r J e 9 u f i f 1 0 t M W P d T J u P E U E k W g 8 K E I x O h + f l o y B Q l h k 8 t w U Q x u y s i Y 6 w w M T a k k g 3 B W z 5 5 l b Q v q l 6 t e n V f q z R u 8 j i K c A K n c A 4 e X E I D 7 q A J L S A g 4 R l e 4 c 3 R z o v z 7 n w s S g t O 3 n M M f + B 8 / g B Z u I 9 r < / l a t e x i t > µ = 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 y f s F u U X Q 1 Q G D g c b w 7 e P m C 8 3 U R s = " > A A A B 7 n i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 9 m V g n o Q i l 4 8 V r A f 0 C 4 l m 2 b b 0 C Q b k q x Q l v 4 I L x 4 U 8 e r v 8 e a / M W 3 3 o K 0 P B h 7 v z T A z L 1 K c G e v 7 3 1 5 h b X 1 j c 6 u 4 X d r Z 3 d s / K B 8 e t U y S a k K b J O G J 7 k T Y U M 4 k b V p m O e 0 o T b G I O G 1 H 4 7 u Z 3 3 6 i 2 r B E P t q J o q H A Q 8 l i R r B 1 U r s n U n S D / H 6 5 4 l f 9 O d A q C X J S g R y N f v m r N 0 h I K q i 0 h G N j u o G v b J h h b R n h d F r q p Y Y q T M Z 4 S L u O S i y o C b P 5 u V N 0 5 p Q B i h P t S l o 0 V 3 9 P Z F g Y M x G R 6 x T Y j s y y N x P / 8 7 q p j a / C j E m V W i r J Y l G c c m Q T N P s d D Z i m x P K J I 5 h o 5 m 5 F Z I Q 1 J t Y l V H I h B M s v r 5 L W R T W o V a 8 f a p X 6 b R 5 H E U 7 g F M 4 h g E u o w z 0 0 o A k E x v
from which one ultimately retrieves the expression which is the result that was derived in [79].
For the photon-added state, we can perform a completely analogous computation with from which we find that This result immediately shows that this Wigner function is always positive, which implies that it is impossible to remotely create Wigner negativity through photon addition.
In previous work, we highlighted that photon addition always creates Wigner negativity in the mode where the photon is added [90]. What we observe in (37) can be understood as the complementary picture for the other modes. This result also highlights an operational difference between photon subtraction and addition: photon additional is a more powerful tool to locally create Wigner negativity, whereas photon subtraction has the potential to create Wigner negativity nonlocally (i.e. in modes that can steer the mode in which the photon is subtracted).

Subtracting multiple photons
When multiple photons are added or subtracted, or when we chain combinations of addition and subtraction operations, the evaluation of Â g| x f ( x f ) will rapidly become more complicated. A general strategy to approach this problem avoids the explicit evaluation of WÂ( x g ), but rather uses standard techniques for the evaluation of moments of multivariate Gaussian distributions. This ultimately boils down to applying Wick's theorem [97] and summing over all matchings (see Appendix A for details). Even though this task can be implemented numerically, the corresponding analytical expressions quickly become intractable.
To illustrate this method, we consider the multimode scenario where two photons are subtracted in different orthogonal modes, g 1 and g 2 , which implies that the conditioning implements the following map ρ →â This implies that we must apply our formalism withÂ = n g 1n g 2 . To treat this problem with the technique of matchings, we use the Gaussian identity (not that we do not explicitly write the dependence on x f to simplify notation) where . . . g| x f denotes the non-displaced version of the distribution. We can immediately identify subsequently, we obtain from (33) that and finally we find new types of terms, that are given by and Computation required to obtain the final result is tedious but straightforward. We find that where we have defined the submartix C as the off-diagonal block of V g| x f via Non-zero entries in the block C can occur due to various caused. First of all, it can be due to a correlation between the modes g 1 and g 2 in the initial Gaussian state (as seen from the term V g in (16)). However, non-trivial entries in C also arise when modes g 1 and g 2 are both correlated to the same modes in f, which is induced by the term V gf V −1 f V T gf in (16). The result (44) directly show the appearance of a trivial term, (trV , which multiplies the effect of photon subtraction in g 1 with that of photon subtraction in g 2 . However, when both modes are sufficiently "close" to each other, we find the additional terms 2tr(C T C) + 4 ξ T g 1 | x f C ξ g 2 | x f , which can be interpreted as some form of interference between the two photon subtractions. Fig. 3 provides an illustration, where we inject three pure squeezed vacuum states into a series of beamsplitters to generate an entangled three-mode state from which we subtract two photons. The first two squeezed vacuum states have 5dB squeezing in opposite quadratures and are mixed on beamsplitter with 75% transmittance. One of the output ports will serve as mode g 1 , whereas the other is injected into a section beamsplitter of 25% transmittance. In the other input port of this beamplitter, we inject the third squeezed vacuum state, which is also squeezed by 5dB. One of the output ports of the 25% transmittance beamsplitter serves as mode g 2 , and in the other output port we find mode f , which is the mode for which we reconstruct the output Wigner function using (44). Photon subtraction is represented by a highly transmitting beamsplitter which sends a small amount of light to a photon detector. Two-photon subtraction then happens when both detectors click at the same time, and we can condition the state in mode f upon this detection outcome. This postselection scheme effectively implements the operatorsâ g 1 and a g 2 on modes g 1 and g 2 , respectively.
We observe that the conditional state W f |Â ( x f ), witĥ A =n g 1n g 2 reaches negative values in two distinct regions of phase space. Indeed, with the Williamson decomposition of V g| x f we can quantify [30] the strength of EPR steering from mode f to the set of modes g to be µ = 0.548. Furthermore, the fact that there are two negativity regions is a hallmark of the subtraction of two photons. This example shows that our framework is a highly versatile tool for CV quantum state engineering.
Finally, we consider the complementary scenario where two photons are subtracted from one mode. In this case, we can still use the perfect matching technique (39), when creation and annihilation operators are in normal ordering. In this case, we obtainÂ =â † gâ † gâgâg , and analogously to (39), we find that This result can then directly be inserted in (19) to obtain the final conditional state for the set of modes f when two photons T = 0.75 T = 0.25 Figure 3. Conditional state Wigner function W f |Â ( x f ), obtained by subtracting a photon in two of the three modes in a three-mode entangled state. This entangled state is generated by mixing three squeezed vacuum states in a sequence of beamsplitters with transmittance of 75% (left) and 25% (right). Two of the squeezed vacuum states are squeezed by 5dB in the x-quadrature (left, right) and one is squeezed by 5dB along the p quadrature (middle). The photon subtraction is represented by highly transmitting beamsplitters which send a small fraction of light to a photon detector, which effectively implements the operatorsâ g 1 andâ g 2 on the mode g 1 and g 2 , respectively.
are subtracted in mode g. As expected, the subtraction of two photons can induce Wigner negativity only when there is EPR steering from the modes f to mode g.
As such, we have shown that our framework allows us to analytically describe conditional non-Gaussian states in a regime which is highly challenging for many other methods. For example, it is highly challenging to approach the problem with the correlation function methods of [90], even though this method was highly successful for single-photon subtraction in multimode states.
These methods can in principle be extended to deal with higher numbers of added and/or subtracted photons in various modes. However, it must be emphasised that one will quickly encounter practical boundaries as finding all possible matchings is a computationally hard problem [98]. Finding an exact description of the Wigner function that is obtained by subtracting a large number of photons from a subset of an entangled Gaussian state seems to be a computationally hard problem that has its roots in graph theory. The problem of finding all matchings also lies at the basis of Gaussian boson sampling [99,100], and it is not expected to be easy to overcome. The problem of Gaussian boson sampling can in turn also be related to CV sampling from photon-added orsubtracted states [101].

VI. CONCLUSIONS
We presented a general framework that describes the Wigner function that is obtained by applying an arbitrary operation on a subset of modes of a multimode Gaussian state, and conditioning the remaining modes on this operation. The most natural way of interpreting this scenario is by considering this operation to be a measurement, such that the state of the remaining modes is obtained by post-selecting on a specific measurement outcome, as is the case for heralding. However, this framework can also be used to study the non-local effects of photon addition and subtraction.
Our framework relies heavily on classical probability theory, and in particular on properties of conditional probability distributions (14). We use the fact that Gaussian states have positive Wigner functions, such that associated conditional probability distributions on phase space are well-defined as probability distribution (but not necessarily as quantum states, because they can violate the Heisenberg inequality). In this regard, our general results (18 -20) are valid for all initial states with a positive Wigner function.
Gaussian states are not only the most relevant initial states from an experimental point of view, they also have the theoretical advantage of leading to a Gaussian conditional probability distribution. The latter is an enormous advantage for evaluating the crucial quantity Â g| x f , as defined in (18). On a more fundamental level, we note that the covariance matrix (16) of this Gaussian conditional probability distribution is essential in the theory of Gaussian EPR steering. This observation allows us to directly prove that Gaussian EPR steering is a necessary prerequisite for the conditional preparation of Wigner-negativity, regardless of the conditional operation that is performed.
In previous work, we already showed that Gaussian EPR steering is also a sufficient ingredient for the remote preparation of Wigner-negativity, in the sense that there always exists a combination of a Gaussian operation and photon subtraction in the modes g that induces Wigner-negativity in the modes f. We thus establish a fundamental relation between Gaussian EPR steering and the ability to prepare a Wigner-negative state in correlated modes. This result is particularly important in the light of measurement-based quantum computation, where large Gaussian cluster states form the backbone for implementing a quantum algorithm. The actual computation is then executed by performing measurements (or more general operations) on some modes of the cluster, in order to project the remainder of the system in a desired quantum state. To claim that such a computation is universal, one must be able to induce Wigner-negativity. Our results therefore show that EPR steering is an essential figure of merit in these cluster states in order to claim that a cluster state is suitable for universal quantum computation.
From a practical point of view, the examples in Section V show that our framework is highly versatile. However, it also highlights the boundaries of analytical treatments. Even though the obtained expression (18) for Â g| x f is easy to interpret conceptually, the actual evaluation can still be challenging. Regardless, we must emphasise that the elegance and simplicity of our framework does allow us to obtain results with far greater ease than previously possible. Many of the methods known in literature are either hard to generalise to arbitrary Gaussian initial states [44,59], focused on one par-ticular measurement or operation [47,90,94], or are just generally hard to interpret or use analytically. Our framework can be applied to any initial Gaussian state, and any conditional operation, provided the Wigner function ofÂ is known.
As such, our results provide a starting point for investigating a wide range of new questions related to multimode conditional preparation of non-Gaussian states. By establishing a fundamental relation between EPR steering and Wignernegativity, we specifically highlight that this framework is also suited to obtain general analytical results, which is often challenging in the study of states that are, both, highly non-Gaussian and highly multimode.