A Central limit theorems reviewed

# Gaussification and entanglement distillation of continuous variable systems: a unifying picture

## Abstract

Distillation of entanglement using only Gaussian operations is an important primitive in quantum communication, quantum repeater architectures, and distributed quantum computing. Existing distillation protocols for continuous degrees of freedom are only known to converge to a Gaussian state when measurements yield precisely the vacuum outcome. In sharp contrast, non-Gaussian states can be deterministically converted into Gaussian states while preserving their second moments, albeit by usually reducing their degree of entanglement. In this work – based on a novel instance of a non-commutative central limit theorem – we introduce a picture general enough to encompass the known protocols leading to Gaussian states, and new classes of protocols including multipartite distillation. This gives the experimental option of balancing the merits of success probability against entanglement produced.

###### pacs:
03.67.Bg,42.50.Ex

Entangled quantum states are the fundamental resources that enable quantum key distribution, quantum communication, and instances of distributed quantum computing. Real physical systems are affected by decoherence and non-ideal apparatus that degrades the quality of experimentally preparable quantum states. However, entanglement distillation protocols provide a means of converting many copies of partially entangled states into a smaller number of more entangled states (1). When entanglement is required over very long distances, distillation can be implemented at regular intervals called repeater nodes (2); (3); (4). Photonic systems that carry entanglement in continuous degrees of freedom are difficult to manipulate arbitrarily. However, so-called Gaussian operations are more easily implemented by a combination of beam-splitters, phase shifters and squeezers. Furthermore, preparation of Gaussian states is routine in many laboratories, and such states are especially useful for numerous quantum information tasks.

Unfortunately, Gaussian operations are quite limited in their capacity to distill entanglement. In particular, a series of no-go theorems have shown that with only Gaussian resources, it is impossible to increase entanglement (5). These results can be circumvented when non-Gaussian resources are available. Given an appropriate resource, Gaussian operations can simultaneously increase the entanglement and make the state more Gaussian, a process we refer to as the “Gaussification protocol” (GP) (6) (for steps towards experimental realization of this and related protocols, see Refs. (7)). Alternatively, the theorem can also be circumvented by using a non-Gaussian operation implemented by photon detectors. Entanglement distillation can be achieved by a combination of photon subtraction, a de-Gaussifying operation, and Gaussification (6); (8).

While these techniques give hope for simple realizations of quantum information protocols, they do not exploit the richness of Gaussian operations. Specifically, GP utilizes only projections onto the vacuum. These projections are feasible if reliable detectors are available that distinguish zero from one or more photons. However, strictly speaking these detectors do not fall within the realm of Gaussian devices. Performing eight-port homodyne detection (9) and postselecting on the vacuum outcome achieves the same projection, but this will have zero success probability when postselecting exactly on the vacuum measurement outcomes. In this work, we prove Gaussification for a wide class of truly Gaussian protocols with nonzero, and tunable, success probabilities for multimode states. All known feasible distillation protocols, including our protocols, are iterative and consume a number of copies exponential in the iterations required. As such, our scheme’s capability to increase success probability significantly improves the prospects of distillation and repeater implementations with only modest resources. Our techniques also open up the perspective of directly distilling into multi-partite Gaussian states.

In addition to the practical applications of our results, our analysis provides a more intuitive explanation of the phenomena of Gaussification. In the original distillation protocols (6), the process of Gaussification is quite mysterious, but it is very apparent in the protocol of Ref. (10) which provides an alternative method that uses no measurements, referred to as the “extremality protocol” (EP), as it is used to show the extremality of Gaussian states with respect to several properties. Yet EP can convert many non-Gaussian states into a more Gaussian state while conserving the expectation value of observables quadratic in position and momentum. Although the EP does Gaussify, its capacity for increasing entanglement is restrained by its deterministic nature. The proof of EP elegantly employs the central limit theorem (11) that explains the ubiquity of Gaussian distributions in classical statistics and nature itself. Our approach unifies GP and EP within a comprehensive theory of “general Gaussification” (GG) founded on a non-commutative central limit theorem and so provides an intuitive mechanism for Gaussification. In addition to bipartite entanglement distillation, our approach reveals whole new classes of protocols that we discuss.

We consider GG protocols that can be implemented iteratively, with the iteration as follows:

1. Take two copies of an -mode state shared between parties;

2. Each party applies a 50:50 beam-splitter transformation between their pair of modes;

3. Every party makes a Gaussian measurement on the output of one beam-splitter port;

4. The parties compare measurement results and postselect such that the operation implemented is Gaussian;

5. The output state is used for the next iteration.

Formal description of GG protocols. There are modes involved in the protocol, and we label annihilation operators, , with two indices; The index labels the respective party and the copy at a particular node. The beam-splitters, in the Heisenberg picture, perform

 U^aj,kU†=(^aj,1+(−1)k^aj,2)/√2. (1)

The measurements at step 3 can be homodyne, eight-port homodyne, or any other Gaussian measurements projecting onto a state for measurement outcome . We are interested in Gaussian protocols that postselect on a set of measurement outcomes, and mix over all accepting outcomes,

 Missing dimension or its units for \mskip (2)

where denotes a partial trace over the second copy, with . The integral over measurement outcomes is weighted by . The weights and correspond to a rejection and an acceptance, respectively, but we also allow for probabilistic strategies where gives the probability of acceptance. The protocol is described by a single operator we call the filter

 Π=∫dmP(m)Πm. (3)

We allow for arbitrary Gaussian filters, , that are invertible and proportional to a fully separable Gaussian state with vanishing first moments and finite energy. Certain interesting cases, such as GP, EP and protocols using precise homodyne detection are included as limits within this family of filters, so keeping full generality. With this notation

 ρn+1∝tr2[U(ρn⊗ρn)U†(1l⊗Π)]. (4)

Phase space. Before we give our results, we review phase space representations for the position and momentum observables of an -mode system, labeled as

 ^R=(^R1,^R2,…,^R2m−1,^R2m)=(^X1,^P1,…,^Xm,^Pm),

where and for . The canonical commutation relations between the coordinates are embodied in the symplectic matrix . The covariance matrix of an operator records the second moments of these observables.

 (ΓA)j,k=tr({^Rj−(dA)j,^Rk−(dA)k}+A), (5)

where denotes the anticommutator and first moments are . Furthermore, we make use of characteristic functions, , that encode all the information of as , where is the displacement operator, . Such a function is said to be Gaussian when

 χA(r)=exp(ir⋅dA−rTΓAr/4). (6)

If the operator is a physical state, its covariance matrix will be real. However, an instrumental tool in our analysis is that we work with . Indeed, we will not employ characteristic functions of states satisfying the conditions of Bochner’s theorem (11); (12), but of more general objects, hence leading to more general complex-valued functions. Since is invertible, uniquely defines a quantum state .

A new non-commutative quantum central limit theorem. With these definitions at hand we can state our first result, with a stronger form of convergence demonstrated later.

###### Theorem 1 (Convergence of general Gaussifier protocols)

Consider an initial state , with associated operator such that the following conditions are satisfied: (i) ; (ii) for all ; (iii) the covariance matrix

 Γρ∞=(ΓΠ−iΣ)(ΓΠ−Γσ)−1(ΓΠ+iΣ)−ΓΠ, (7)

exists and is positive definite. Let denote the Gaussian state with covariance matrix . GG with filter causes to weakly converge to in the following sense: If and are eigenvectors of , then for all there exists a such that for all

 |⟨x|ρn|y⟩/tr(ρnΠ)−⟨x|ρ∞|y⟩/tr(ρ∞Π)|<ϵ. (8)

Proof of the statement. Much of the basic structure of the proof follows an argument of a non-commutative quantum central limit theorem for quantum states (11); (12); (10); (13). The problem once one allows for Gaussian measurements is that the characteristic function of the output is unwieldy indeed. We circumvent this problem by a bold step: We put in an additional filter at the output “by hand”, in order to exploit symmetry, at the expense of having to consider of different objects, . This will then lead to the desired result. From Eq. (4) we have,

 χσn+1(r)∝tr[(D(r)⊗1l)U(ρn⊗ρn)U†(Π⊗Π)]. (9)

Using the cyclicity of the trace and Eq. (1)

 χσn+1(r)∝tr[D(r/√2)⊗2ρ⊗2nU†Π⊗2U]. (10)

Next we recall that Gaussian states with zero first moments commute with beam-splitters, , such that

 χσn+1(r) ∝ tr{[D(r/√2)ρnΠ]⊗2}∝tr[D(r/√2)σn]2, (11) = χσn(r/√2)2=χσ(r/√N)N,

where in the last equality , iterating the formula. The key to the simplicity of this formula is to consider convergence of rather than directly . By introducing an additional projector within the trace, symmetry allows us to commute through the beam-splitter unitaries, which is the essential simplifying step. To find the limiting characteristic function, we consider a given phase space point and the function defined as . In the spirit of a classical central limit theorem (11); (10); (13) but for non-Hermitian operators we write

 fσ,r0(t√N)N ∝ (1−t2rT0Γσr04N+o(t2N))N, (12)

which converges assuming that second moments are finite, and for all . This last condition, which is always satisfied for classical characteristic functions, may be violated for non-Hermitian . However, provided for all , we find (see App. A) in the limit of large that

 limn→∞fσn+1(t) ∝ exp(−t2rT0Γσr0/4), (13)

pointwise in . Setting shows pointwise convergence of for each phase space point . Furthermore, following the reasoning of Refs. (12); (13), this entails that any operator that is absolutely integrable such that , has a convergent expectation value (see App. B). Setting where , then for large

 ⟨x|ρn|y⟩tr(ρnΠ)=⟨x|σn|y⟩λy→⟨x|σ∞|y⟩λy=⟨x|ρ∞|y⟩tr(ρ∞Π). (14)

All that remains is to show that Gaussian entails a Gaussian . First we observe (see App. C) that the product of two operators has a characteristic equation

 χσ∞(q)∝∫χρ∞(r)χΠ(q−r)exp(−irTσq/2)dr. (15)

For Gaussian and , the integral is a multi-variate Gaussian integral that evaluates (see App. C) to another Gaussian with the covariance matrix taking the form of a Schur complement (15),

 Γσ=ΓΠ−(ΓΠ+iΣ)(Γρ∞+ΓΠ)−1(ΓΠ−iΣ), (16)

Rearranging this formula for gives us the covariance matrix for the convergent state in terms of and as in the theorem.

Examples. We have introduced a broad class of protocols for which our theorem indicates Gaussification. However, for concreteness it is helpful to keep in mind a simple class of protocols. Consider when each party performs eight-port homodyne measurements that project onto a coherent state. When the party obtains outcome projecting onto coherent state vector , we declare the iteration a success with probability . The degree of postselection is quantified by a real variance . It follows that the filter is indeed proportional to a Gaussian state with covariance matrix , with . This class of protocols is important as it contains GP and EP as limits and , respectively.

We now consider the degree of entanglement that is achieved by applying our protocol with eight-port homodyne measurements. First we consider the well-studied bipartite state vector , and present the log-negativity (16) of the convergent Gaussian in Fig. 2. Varying the parameter interpolates between the entanglement achieved by GP and EP, with increased yield compensating for reductions in entanglement. We also analyzed a tripartite entangled state vector , having up to two photons in three modes, and for the log-negativity summed over all bipartitions. In fact, in this way, one can straightforwardly engineer multi-partite hybrid distillation protocols for quantum networks, giving rise to primitives in repeater architectures where entanglement is shared across many repeater nodes. These would overcome the known limitations providing a road-block against entirely Gaussian continuous-variable repeater networks (17).

Strong convergence. Our previous theorem proves a convergence result identical to that of Ref. (6) and GP. However, it would often be preferable to have convergence of to in a stronger sense (without the addition factors of and and as a convergence in trace-norm). For most physically relevant (see App. D) instances of initial quantum states, we now show that this is indeed the case. We will make use of expectation values of normally order operators,

 αx,yn=tr[V(⊗mk=1^axkk)†(⊗mj=1^ayjj)V†σn], (17)

where and is the Gaussian unitary such that is a thermal state.

###### Theorem 2 (Strong convergence)

In addition to Theorem 1, if for all the expectations values of normally ordered operators satisfy, for . It follows that for all there exists an such that for all , we have , where is the trace norm.

The conditions of the theorem are stated technically, but physically prevent overpopulation of higher Fock numbers. E.g., it is easy to check these conditions are meet in all the low photon examples analyzed in Fig 2. We begin by first showing converges to , which in terms of is

 tr(Π−1σ)=tr(Π−1ρΠ)/tr(ρΠ)=tr(ρΠ)−1. (18)

The inverse filter, , has an exponential form that can be Taylor expanded and, using the bosonic commutation relation, normally ordered such that , where . Hence, when the conditions of our theorem hold, we conclude , and we proceed by showing this holds for all . Iteratively we have (details in App. E)

 αx,yn+1=∑u≤x,v≤yCx,yu,vαu,vnαx−u,y−vn, (19)

where is a combinatorial quantity that is non-negative and real. Under our assumptions the absolute values obey

 |αx,yn+1|≤∑u≤x,v≤yCx,yu,vαu,v∞αx−u,y−v∞=αx,y∞, (20)

so initially satisfying our conditions entails satisfaction for all . Hence, for all we deduce and equivalently . Next we bound from above. Consider a finite rank projector that commutes with , then has an absolutely integrable characteristic function. Hence, for arbitrarily small there exists an such that for

 tr(Bσn)−tr(Bσ∞)=tr(Pρn)tr(Πρn)−tr(Pρ∞)tr(Πρ∞)≥−δ. (21)

Since , and can be chosen so is arbitrarily close to unity, we conclude where is again small. This inequality rearranges to , giving an arbitrarily tight bound from above. Combined with our lower bound we conclude that , under the stated assumptions, converges to . As such converges to , and as is well known this entails trace norm convergence (See Ref. (12) and App F).

Summary. We have introduced a framework for constructing a range of new protocols for entanglement distillation and manipulation. At the same time, this work provides a unified framework for existing protocols leading to Gaussian states: Notably, the mysterious emergence of Gaussian states in distillation schemes is once again related to an instance of a quantum central limit theorem, albeit for a much more broader class of protocols than previously considered. This framework also allows us to look at Gaussification in experimentally realistic acceptance windows, and to trade-off different figures of merit against each other. Such trade-off control is essential as previous proposals are so heavily postselective that over several iterations the success probability would reduce dramatically, whereas our protocols offer an arbitrarily good chance success. Potential for future research is broad as a unique protocol is defined by every separable Gaussian state. For example, our techniques can be applied to homodyne detection protocols for either multimode entanglement distillation, or single mode squeezing enhancement.

Acknowledgements. We would like to thank the EU (QESSENCE, MINOS, COMPAS), the BMBF (QuOReP), and the EURYI for support, and thank A. Mari, C. Gogolin, V. Nesme, D. E. Browne and M. Ohliger for valuable discussions.

## Appendix A Central limit theorems reviewed

Here we present additional details on the convergence of characteristic functions to a Gaussian function. Our proof follows the classical proofs, such as in Ref. (14), but since is generally neither positive nor Hermitian there are some notable subtleties and technicalities involved. Ultimately, we aim to prove the following.

###### Theorem 3 (General quantum central theorem)

If satisfies for all then in the limit of large

 fσ,r0(t/√N)N→e−νt2/4, (22)

pointwise in , with as the appropriate second moment. Specifically, pointwise convergence means that for all and all there exists a such that for all we have

 |fσ,r0(t/√N)N−e−νt2/4|≤ϵ. (23)

Setting directly entails pointwise convergence of the characteristic function , and furthermore this is uniform on compact regions, such as a ball in phase space.

###### Corollary 1 (Convergence in compact regions)

Within any compact region , converges uniformly to . That is, for any there exists an such that for all and all , we have .

Pointwise convergence is uniform on compact regions whenever the gradient is bounded within for all . To prove a gradient bound for all , we first observe that for a bounded continuos function always has a gradient bound on compact sets. That is, an can always be found such that for all we have , where denotes the 2-norm. This property is iteratively preserved for all , since

 Gn(r) = ∥∇χσ0(r/√N)N∥2, = N∥[∇χσ0(r/√N)].[χσ0(r/√N)N−1]∥2, ≤ √NG0(r/√N)≤η|r|.

Using the point in the compact region that gives the largest value giving the desired gradient bound for all Ê and all . This corollary will prove useful in App. B.

We prove several useful lemmas before directly addressing this proof. First we observe,

 fσ,r0(t) = tr[D(tr0)σ], (24) = ∫∞−∞dxeixtFσ,r0(x),

where is the representation of along a line in phase space defined by the unit vector (e.g., the position or momentum representation in case of a unit vector pointing along the axis of phase space). This is similar to a probability density for measurement of a particular quadrature; however, because is neither positive nor Hermitian, the function is not generally positive or real, which is the main cause for caution. Yet it is sufficiently well behaved, with the two following lemmas proving useful.

###### Lemma 1 (Convergence of absolute integrals)

The absolute integral satisfies .

Now clearly , such that

 Π2+ρ2≥iρΠ−iΠρ (25)

and , so that

 Π2+ρ2≥Πρ+ρΠ. (26)

Similarly, and . One finds after a few steps that, for any vector , we have

 |⟨ψ|ρΠ|ψ⟩|≤|⟨ψ|(ρ2+Π2)|ψ⟩|/√2. (27)

Hence, using the triangle inequality,

 |Fσ,r0(x)|≤|FΠ2,r0(x)|+|Fρ2,r0(x)|√2|tr(ρΠ)|. (28)

That is to say, whenever the second moments of and are finite, the absolute integral converges to a finite value.

The second lemma we require is as follows.

###### Lemma 2 (Expansion of functions)

For a given unit vector and an operator that is unit trace and zero displacement , the above function can be expressed as

 fσ,r0(t)=1−14νt2K(t) (29)

where as .

This lemma asserts that higher order moments drop off sufficiently quickly for the purposes of establishing a central limit theorem. We make use of the following formulation of the exponential function

 eitx=1+itx−12t2x2L(x,t) (30)

where

 L(x,t)=2∫10(1−u)eiutxdu, (31)

satisfying

 |L(x,t)|≤1 (32)

for . Furthermore, for any finite interval, as uniformly in . Consequently, we also have as uniformly on any finite -interval. Multiplying Eq. (30) by and integrating w.r.t.  and using the unit trace of and vanishing first moments we have

 fσ,r0(t)=1−12t2∫∞−∞x2L(x,t)Fσ,r0(x)dx, (33)

which we can rewrite as

 fσ,r0(t)=1−12t2νK(t)dx, (34)

where is the second moment defined by

 ν=2∫∞−∞x2Fσ,r0(x), (35)

and is the quantity satisfying

 νK(t)=2∫∞−∞x2L(x,t)Fσ,r0(x). (36)

If we subtract Eq. (35) from the above, we have

 ν(K(t)−1)=2∫∞−∞x2(L(x,t)−1)Fσ,r0(x). (37)

The absolute value satisfies

 ν|K(t)−1|≤2∫∞−∞x2|(L(x,t)−1)Fσ,r0(x)|. (38)

Noting that since , and dividing the integral into two parts over some finite interval and the complement we deduce that

 ν|K(t)−1| ≤ 2∫Rx2|(L(x,t)−1)Fσ,r0(x)| (39) + 4∫R⊥x2|Fσ,r0(x)|.

We now make use of Lem. (1) which tells us that the second integral always converges to a finite value. Therefore, for any the region can be choose sufficiently large that the second integral is bounded such that

 ν|K(t)−1| ≤ 2∫Rx2|(L(x,t)−1)Fσ,r0(x)|+ϵ. (40)

Recall that for finite interval, such as , uniformly goes to zero with vanishing . Hence, for any and there exists a such that for all we have

 ν|K(t)−1| ≤ 2ϵ. (41)

We conclude that converges to 1 with vanishing , proving our lemma.

We now employ the above lemmas to demonstrate pointwise convergence. We wish to bound the quantity

 δ(t,N)=|fσ,r0(t/√N)N−e−νt2/4|, (42)

which can be expanded as

 δ(t,N) = |(1−νt24NK(t/√N))N−e−νt2/4|. (43)

In the limit of large we have that

 δ(t,N) = ∣∣ ∣∣(1−νt24NK(t/√N))N−(1−νt24N)N∣∣ ∣∣.

For any and satisfying it is well known that

 |aN−bN|≤N|a−b|. (45)

Applying this inequality to our problem gives

 δ(t,N) ≤ N|(1−νt24NK(t/√N))−(1−νt24N)|, (46) ≤ |νt2(K(t/√N)−1)|/4.

For fixed , the argument of vanishes as increases, entailing pointwise convergence. Consider a point , for any we can always find a such that for all . The variable defines a , such that whenever and so

 δ(t,N>Nϵ) ≤ ϵ. (47)

Hence we have shown Thm. 3. Furthermore, this convergence is uniform on finite intervals of .

## Appendix B Convergence of expectation values

In the letter we state that convergence of to a Gaussian function with covariance matrix entails convergence for the expectation value of operators, , with absolutely integrable characteristic functions. Here we prove the relevant lemma, and also show that its conditions are meet for operators of the form where and are, up-to a Gaussian unitary, multimode Fock states. Note that, when are positive Hermitian operators, stronger conclusions can be reached by a similar means as shown in Ref. (11).

###### Lemma 3 (Convergence of expectation values)

Consider a sequence of trace class m-mode operators and a limiting operator , such that (i) for all and (ii) For any compact region, , converges uniformly to . Consider also a trace class operator satisfying and any . There exists an such that for all

 |tr(Bσn)−tr(Bσ∞)|≤ϵ. (48)

For the sequence of characteristic functions considered in the letter, property (i) is a prerequisite and property (ii) was shown in Cor. (1). We state these properties again here for clarity. First we note that for any two m-mode trace class operators, and , the characteristic functions, and , satisfy

 ∫drχB(r)χA(r)=(2π)mtr(BA). (49)

This entails that

 Dn = (2π)m|tr(Bσn)−tr(Bσ∞)| (50) = |∫χB(r)[χσn(r)−χσ∞(r)]dr|. (51)

Next we observe that for any satisfying , and for any we can find a function with compact support such that

 ∫|χB(r)−h(r)|dr≤ϵ. (52)

That is, the set of functions with compact support are dense in the set of absolutely integrable functions. Defining and , we have

 Dn = |∫[η(r)−h(r)]Δn(r)dr|, (53) ≤ |∫η(r)Δn(r)dr|+|∫h(r)Δn(r)dr|, ≤ sup(|Δn(r)|)∫|η(r)|dr+|∫h(r)Δndr|, ≤ 2ϵ+|∫h(r)Δndr|.

The third line uses a Hölder’s inequality. The final line uses Eq. (52) and . For sufficiently large , the second term can also be made arbitrarily small because has compact support and vanishes uniformly on any compact set, as required by property (ii). This completes the proof of the lemma.

Furthermore, any operator will satisfy the conditions of the lemma when and are eigenvectors of a Gaussian density operator. That is, there exists a Gaussian unitary such that and are Fock states. Clearly is trace class, and so we only need show absolute convergence of . For a 1-mode element in the Fock basis the characteristic function will be a product of a real polynomial in of finite degree and a decaying Gaussian in . Such functions are easily seen to absolutely converge, and so too shall a product of such functions. A Gaussian displacement only adds a phase factor to the characteristic function, and non-displacing Gaussian unitaries are equivalent to a unitary change of variables in the integral. Hence, Gaussian unitaries will not alter the value of the absolute integral. Indeed, most trace class operators are easily seen to have absolutely integrable characteristic functions. However, we do not currently know if this is generally true.

## Appendix C Gaussian integrals and the Schur complement

The product of two operators, e.g., has a characteristic function that equals an integral involving the characteristic functions for and as expressed by Eq. (15) of the letter. Furthermore, for Gaussian states this integral is solved by taking the Schur complement of an appropriate matrix (see Eq. (16) of the letter). These techniques are standard, but for the purpose of completeness we shall review them here.

As a side remark, note that an alternative strategy towards formulating generalized quantum central limit theorems would have been to consider objects of the form instead of . While this constitutes an alternative in principle, the resulting Schur complement expressions become much more involved, motivating the path taken above.

We begin by noting that for any trace class operator , the formula for the characteristic function can be inverted such that

 A∝∫drχA(r)D(−r). (54)

Hence, the product of two operators satisfies

 AB∝∫∫drdr′χA(r)χB(r′)D(−r)D(−r′). (55)

The Campbell-Baker-Haussdorf formula entails that

 D(−r)D(−r′)=D(−(r+r′))exp(−irTΣr′/2), (56)

and so

 AB∝∫drdr′χA(r)χB(r′)D(−(r+r′))exp(−irTΣr′/2). (57)

Changing variables to and yields,

 AB∝∫drdqχA(r)χB(q−r)D(−q)exp(−i2rTΣ(q−r)). (58)

Noting that for all eliminates part of the exponent. Next we observe that this has the form of Eq. (54) and so we can deduce that the characteristic function is

 χAB(q)∝∫drχA(r)χB(q−r)exp(−i2rTΣq). (59)

Setting and provides Eq. 15 of the letter.

Until now we have not used any properties of or , but herein assume they are Gaussian with covariance matrices and . Hence

 χAB(q) ∝ ∫drexp(−rTΓAr4) (60) × exp(−(q−r)TΓB(q−r)4−i2rTΣq).

This can be expressed more compactly using a single matrix

 χAB(q) ∝ ∫exp(−(q,r)TM(q,r)/4)dr, (61)

where is the Block matrix

 M=(ΓB−ΓB2iΣ−ΓBΓA+ΓB). (62)

Observing that allows explicit symmetrization,

 M=(ΓB−ΓB−iΣiΣ−ΓBΓA+ΓB), (63)

such that now . The Schur complements arises by decomposing the matrix as , with

 D = (Sc(M)00ΓA+ΓB), (64) N = (1l−(ΓA+iΣ)(ΓA+ΓB)−101l), (65)

where denotes the Schur complement of ,

 Sc(M)=ΓB−(ΓB+iΣ)(ΓA+ΓB)−1(ΓB−iΣ). (66)

This decomposition requires only that is symmetric and that is invertible. Returning to our Gaussian integral, a change of variables gives

 χAB(q) ∝ ∫exp(−(q,r)TD(q,r)/4)dr, (67)

and so the and variables decouple, and the integral evaluates to some number and

 χAB(q) ∝ exp(−qTSc(M)q/4)dr. (68)

Hence, we have that . By substituting in and we obtain Eq. (16) of the letter.

## Appendix D Notions of weak convergence

Here we briefly discuss technicalities related to notions of weak convergence versus convergence in trace-norm. Interestingly, these different notions can have a physical implication. While in many experimentally relevant examples the conditions of Thm. (2) of the letter are satisfied, one can concoct exotic instances for which strong convergence fails. For example, if we consider the -mode state vector

 |ψ0⟩∝|0,0⟩+0.1|1,1⟩+6|2,2⟩ (69)

and apply the GP protocol, then Thm. (2) of the letter holds for a target Gaussian

 |ψ∞⟩∝∑n(1/10)n|n,n⟩. (70)

However, once correctly normalized, straightforward numerics show the state diverges with fidelity vanishing with . This highlights the importance of proving more robust convergence results. Stronger convergence has only been established for EP, where trivially since . On a theoretical level, it is compelling to ask what the necessary and sufficient conditions are for strong convergence of GG.

## Appendix E Iterative formula

Here we go through the details of the iterative formula of Eq. (19). For brevity, we assume is diagonal in the Fock basis, and so , though the more general case can be proven by the same method but more cluttered notation. For the round we have that the state satisfies

 σx,yn+1 = tr2[U(ρn⊗ρn)U†(Π⊗Π)]tr[U(ρn⊗ρn)U†(Π⊗Π)] = tr2[U(ρn⊗ρn)U†(Π⊗Π)]tr[ρnΠ]2

and so the expectation values satisfy

 αx,yn+1 =