Non-producibility of arbitrary non-Gaussian states using zero-mean Gaussian states and partial photon number resolving detection

Gaussian states and measurements collectively are not powerful-enough resources for quantum computing, as any Gaussian dynamics can be simulated efficiently, classically. However, it is known that any one non-Gaussian resource -- either a state, a unitary operation, or a measurement -- together with Gaussian unitaries, makes for universal quantum resources. Photon number resolving (PNR) detection, a readily-realizable non-Gaussian measurement, has been a popular tool to try and engineer non-Gaussian states for universal quantum processing. In this paper, we consider PNR detection of a subset of the modes of a zero-mean pure multi-mode Gaussian state as a means to herald a target non-Gaussian state on the undetected modes. This is motivated from the ease of scalable preparation of Gaussian states that have zero mean, using squeezed vacuum and passive linear optics. We calculate upper bounds on the fidelity between the actual heralded state and the target state. We find that this fidelity upper bound is $1/2$ when the target state is a multi-mode coherent cat-basis cluster state, a resource sufficient for universal quantum computing. This proves that there exist non-Gaussian states that are not producible by this method. Our fidelity upper bound is a simple expression that depends only on the target state represented in the photon-number basis, which could be applied to other non-Gaussian states of interest.

Gaussian states and Gaussian unitaries, produced by the action of linear and quadratic Hamiltonians on the vacuum state, have efficient and complete mathematical representations [13][14][15]. Non-Gaussian states is a vast set-it consists of states generated via the action, on the multi-mode vacuum state, of a unitary with Hamiltonian that is a third or higherorder polynomial in the field operators. Therefore, non-Gaussian states are inherently under-explored and their general representations less-understood.
Deterministic realization of non-Gaussian unitary operations, such as the self-Kerr gate [16,17] and the cubic-phase gate [18] is near impossible at optical frequencies [19]. The extreme resource inefficiency resulting from this deficiency, combined with the fact that Gaussian states and Gaussian unitaries are a classically-simulable resource [20], have kept allphotonic quantum computing from being pursued as one of the top contenders for quantum computing, for decades since their invention, despite their obvious importance in optical communications and sensing applications, and it not requiring quantum transduction for networking far-flung quantum processors-a major benefit unique to photonic quantum encodings.
Recent advances in discrete variable (DV), i.e., single-photon-qubit based, quantum computing [16] have revealed that deterministic production of even small non-Gaussian resource states (e.g., a 3-photonentangled GHZ state) can enable resource-efficient universal photonic quantum computing, despite twoqubit gates being inherently probabilistic [1,21]. However, a systematic understanding of efficient production of even such simple non-Gaussian states as GHZ states and realization of simple two-qubit non-Gaussian measurements (e.g., Bell-state measurements) required for DV quantum computing has proven extremely difficult [22].
A major attraction of continuous variable (CV) quantum computing [18] is that large Gaussian entangled (cluster) states [23] can be produced experimentally in a one-shot deterministic fashion [24,25]. Further, CV qubit states, such as the GKP qubit is known to be the most loss-resilient encoding of the qubit in a bosonic mode [26] that admit deterministic Clifford gates using Gaussian unitaries, and there are native CV quantum codes to correct for loss errors [27]. Since Gaussian states are not uni-versal [20], one needs a non-Gaussian operation to enable universal quantum computing [28]. Experimentally, the most readily-available non-Gaussian resource is photon number resolving (PNR) detection [29]. One common modality in which a PNR detector can be used to probabilistically engineer non-Gaussian states is photon subtraction [30,31], which also is known to increase entanglement [32]. Photon subtraction from multi-mode Gaussian states has been achieved experimentally [33][34][35] and several theoretical aspects of photon subtraction have been studied [36][37][38][39][40][41][42][43]. In this paper, we consider a more general and simpler to describe set-up, namely partial PNR detection [44,45], i.e., employing PNR detection on a subset of the modes of a multi-mode Gaussian state, to herald the undetected modes in a desired state, conditioned on the PNR detectors' click pattern on the detected modes. The heralded state is non-Gaussian unless all PNR detectors register zero photon clicks. This is because the projection on vacuum is a Gaussian operation and therefore will not impart any non-Gaussianity to the heralded state). We point out that partial measurements has been used before in different contexts, e.g., to realize minimal-disturbance measurements experimentally [46].
The most general multimode Gaussian state is described by a covariance matrix and a coherent displacement vector. Zero-mean Gaussian states are those whose displacement vector, or the mean field amplitude, is zero. A general K-mode Gaussian state can be produced by passing K displaced-squeezed states through a linear-optical unitary transformation, which in turn admits a systematic design in terms of K(K − 1) 50-50 beamsplitters and an equal number of unspecified phase elements [13,14]. Experimentally, the most challenging part in the above is preparing a displaced squeezed state. Preparation of squeezed vacuum state, or two-mode squeezed vacuum state of light-both of which are zeromean Gaussian states-on the other hand, is routinely performed using spontaneous parametric downconversion (SPDC), e.g., using a χ (2) non-linear medium. Recent experiments have demonstrated on-chip squeezed-vacuum generation [47,48]. Further, since fully-programmable linear optical circuits have also been realized on-chip [49], scalable generation of arbitrary multi-mode zero-mean Gaussian states is well-within the reach of modern technology. This is the reason why we focus in this paper, on evaluating whether arbitrary non-Gaussian states can be prepared just by partial-PNR detection on zero-mean Gaussian states. This paper is organized as follows: In Sec. II we explain how photon subtraction and addition can be seen as a special case of partial PNR. Said section can be skipped by the experienced reader but it possesses some pedagogical value, and it sets notation. In Sec. III we review briefly the mathematical description of partial PNR (see also App. A and [50]) and we also provide a simple proof that any pure zero-mean Gaussian state engineered with partial PNR will necessarily give a zero-mean non-Gaussian state. In Sec. IV we present the fidelity between the heralded state and a given non-Gaussian (target) state. The idea behind the upper bound on said fidelity is based on the Cauchy-Schwartz inequality. In Sec. V, we calculate said bound for any product state of single-mode superposition of the binary phase-shift keyed (BPSK) coherent states | − γ , |γ and give a few examples. In Sec. VI we calculate fidelity upper bounds for multimode entangled states that are superpositions of multimode coherent states where each mode is from the BPSK constellation. The first special case of an entangled state of this type that we consider is the coherent GHZ state, for which our fidelity upper bound comes out to 1 (a trivial upper bound). However, for the coherent-cat basis cluster states (CCCS), our fidelity upper bound evaluates to 1/2, showing such a state cannot be prepared by partial PNR on a zero-mean Gaussian state. Finally, in Sec. VII, we summarize our results: we discuss how the fidelity upper bound relates to the absence of a non-zero mean (or coherent displacement), we put our findings in context with the literature on non-Gaussian quantum state preparation, present further intuition, and briefly discuss future directions of research.

II. PARTIAL PNR AS A GENERALIZATION OF PHOTON SUBTRACTION AND ADDITION
Let us consider photon subtraction on the most general single-mode Gaussian pure state: a squeezedcoherent state |α 1 , r 1 , where α 1 , r 1 ∈ C are the displacement and squeezing parameters respectively. The state interacts with vacuum on a beamsplitter of transmissivity τ ∼ 1. On the reflective (lowtransmissivity) output port of the beamsplitter, PNR detection is performed, which registers, say, n 1 photons. This heralds the subtraction of n 1 photons from the input state. This conditional photon-subtracted state |Φ n1 -i.e., the state heralded on the transmitted port of the beamsplitter-is non-Gaussian whenever n 1 ≥ 1. One can write down the probability of detecting n 1 photons (and hence producing |Φ n1 ) as a function of α 1 , r 1 and τ [45]. A natural generalization of this setup is to allow for further Gaussian resources by substituting the vacuum state at the other input of the beamsplitter with another single mode Gaussian pure state |α 2 , r 2 , α 2 , r 2 ∈ C, and then proceed with PNR detection on one output port, and considering the heralded state on the other output port if n 1 photons are detected. No matter what figure of merit we might choose on the quality of the heralded state (e.g., fidelity to some target state), it can only improve, or at worst remain the same as compared to that using photon subtraction. This is because vacuum is a trivial special case of a pure Gaussian state. We term this setup partial PNR detection [44]: PNR detection on one mode of a two-mode general Gaussian pure state, to seek a desired post-selected non-Gaussian state on the undetected mode.
Next, let us consider photon subtraction on all modes of a K-mode pure Gaussian state (which in general is entangled), by coupling the i-th mode with vacuum on a beamsplitter of transmissivity τ i , i = 1, . . . , K. If we count all the ancillary vacuum states as input modes, we have a 2K-mode Gaussian state K modes of which are detected with PNR, resulting in a K-mode (generally non-Gaussian) state. An obvious generalization is to consider an N -mode Gaussian pure state (where N can be even or odd) and to apply PNR detection on N −M modes, resulting in an M -mode state |Φ n M +1 ,...,n N , conditioned on the PNR pattern (n M +1 . . . , n N ).
Let us note that partial PNR detection schemes incorporates multiple photon addition as well. Photon addition is modeled utilizing a beam splitter whose upper input is the state |Ψ (or a mode of the state) which will undergo photon addition while in the beam splitter's lower input port a Fock state |n is injected. A PNR detector is applied on the lower output port which if heralds m ≤ n photons, then photons have been added to the state |Ψ , resulting to a state |Φ which in general will be non-Gaussian. To produce a Fock state |n , i.e., the input to the lower port of the beam splitter, one can consider a two-mode squeezed vacuum state (TMSV) whose one mode is detected using a PNR detector. Then, this Fock state can be used for the photon addition task. Equivalently, one can include several TMSV states as part of a general Gaussian state and all the PNR detectors to be included in the very last step of the state engineering protocol.

III. MATHEMATICAL DESCRIPTION OF PARTIAL PNR ON ZERO-MEAN GAUSSIAN STATES
Let us consider a zero-mean N -mode Gaussian (in general entangled) pure state |Ψ , prepared by mixing N single-mode squeezed vacuum states of squeezing parameters r i , 1 ≤ i ≤ N , in a linear optical Nmode unitary operation U . We consider partial PNR detection on N − M of those modes. We proved that (see [45] and Section II of [50], also included as App. A of this paper in uniform notation, for the benefit of the reader), conditioned on the PNR detection pattern (n M +1 . . . , n N ), the M -mode heralded state |Φ n M +1 ,...,n N ≡ |Φ can be written in the Fock basis as: where, The probability of obtaining the PNR pattern (n M +1 , . . . , n N ) is given by, and where the loop-Hafnian in Eq. (4) is evaluated for a matrix σ, whose elements are given in the Appendix, in Eq. (A37). Per Eq. (A35), σ is directly related to the covariance matrix of the Q function of the state |Ψ on which the partial PNR is performed. As such, the matrix elements of σ are functions of the squeezing parameters r i and the entries of U . Without loss of generality, we set the modes undergoing PNR to be the 'last' N − M modes of |Ψ . Further, we consider r i > 0 to be real-valued, or equivalently, all phases are pushed into the passive interferometer U that entangles the squeezed vacuum states to create the resource Gaussian state |Ψ .
For the main result of this paper, we will not need to invoke the explicit dependence of σ on |Ψ , through the squeezing parameters r i and the parameters of the entangling passive linear-optical unitary U . The property of importance to us will be the parity of the PNR detector's pattern in Eq. (4).

IV. FIDELITY UPPER BOUND ON THE CONDITIONAL STATE
In Sec. III we proved that a zero-mean Gaussian pure state under partial PNR will necessarily give a zero-mean conditional state on the unmeasured modes. Therefore, it is natural to anticipate that any non-Gaussian target state with an arbitrary nonzero mean-field would not have a fidelity arbitrarily close to 1 with a non-Gaussian state engineered using partial PNR on zero-mean Gaussian states. However, the question of what the highest said fidelity can be, remains open. In this section, we provide a general recipe to find an upper bound to the fidelity for any target state. In subsequent sections, we will apply this technique to evaluate our fidelity upper bound on specific non-Gaussian states of interest.
The fidelity F = | Φ t |Φ | 2 between the conditional state |Φ of Eq. (1) and a non-Gaussian target state It is apparent that if we use the Cauchy-Schwartz inequality on Eq. (5), we will get F ≤ 1. However, we will see that under the constraint (Eqs. (2) and (4)) n 1 + . . . + n N = even, the Cauchy-Schwartz inequality gives a non-trivial upper bound. Then n 1 + . . .+n M = even if the summation of the PNR pattern (n M +1 , . . . , n N ) is even, and n 1 + . . . + n M = odd if the summation of the PNR pattern (n M +1 , . . . , n N ) is odd. Therefore, we rewrite Eq. (5) as, where, Let us consider the case where n 1 + . . . + n M = even. Then F odd = 0 and we can use the Cauchy-Schwartz inequality to get, Finally, exploiting the fact that the state |Φ has nonzero coefficients under the constraint n 1 + . . . + n M = even, we write, and we get, Similarly, for the complementary case where n 1 + . . . + n M = odd, we have that F even = 0 and Four observations are necessary here. First, we observe that 0 ≤ u even ≤ 1 and 0 ≤ u odd ≤ 1 and both bounds depend only on the target state, therefore they are easy to compute. Second, the non-Gaussian target state is normalized, therefore, It is possible that for the desired non-Gaussian target state, u even and u odd to be unequal. In that case, we will use as upper bound the larger among the two, and herald on the PNR pattern whose parity corresponds to that of the higher upper bound. Third, let us note that we view fidelity as necessary criterion for successful non-Gaussian state engineering. For example, a coherent cat state N −1 0 (|γ +|−γ ) (where | ± γ is a coherent state and N 0 is normalization) can have high fidelity with vacuum for small, albeit non-zero, γ amplitude. However, vacuum and small coherent cat states are inherently different. On the other hand, if one derives a low enough upper bound for the fidelity, then the impossibility of producing the state under consideration is certain. Last, we observe that assuming a zero-mean Gaussian resource state, resulted to imposing a specific parity on the PNR pattern. The question now is how this parity constraint impacts the state engineering performance.

V. FIDELITY UPPER BOUNDS FOR COHERENT CAT PRODUCT STATES
Consider a single mode state |c , which is a superposition of two coherent states | ± γ , where b 1 , b 2 ∈ C satisfy, so that c|c = 1. Let us calculate the upper bounds of Eqs. (10) and (11) for the product state |c ⊗M . Following Eq. (10), we have, which can be rewritten as, By separating the fraction of Eq. (15) and using the fact that state |c ⊗M is normalized we get, Using Eq. (13), the Fock basis expansion of a coherent state |γ = exp(−|γ| 2 /2) ∞ n=0 γ n / √ n!|n , and Eq. (14), we find, and using Eq. (12) we get, We observe that u even and u odd for the state of Eq. (13), depend only on the absolute values of the the state's coefficients when expressed as a coherent states' superposition. We note that for our N -mode Gaussian state, the M -mode produced state, and the M target states, we allow N and M, 1 ≤ M < N , to be arbitrary. As applications, we will consider the following target states, where N k = 2[+(−1) k e −2|γ| 2 ], k = 0, 1. The states of Eqs. (19), (20), (21), (22), are the computationaland rotated-basis qubit states corresponding to the coherent cat-basis qubit-one of the leading qubit candidates for all-photonic quantum computing. The states |0 and |1 form the so-called logical qubit basis, while the states |+ and |− are derived by the action of a Hadamard gate (defined on the qubit basis) on the logical qubit basis' kets. For the |0 , |1 states we find, Consistent with the parity of the |0 , |1 states, we see that a PNR pattern whose summation is odd cannot herald the state |0 , while the state |0 is not impossible to be engineered if the summation of the PNR pattern is an even number. Also, if M is an odd (even) number, |1 cannot be heralded if the PNR pattern is summed to an even (odd) number. We note that an upper bound equal to 1 does not mean that the state can be engineered with perfect fidelity. However, high fidelity for generating approximations of the |0 , |1 states has been found in the literature [45,50,51] using the partial PNR method, even with the resource Gaussian state being zero-mean. For the |+ , |− states we find, Since |+ , |− are not zero-mean states, we expect that the upper bound should reflect that by being less than 1. In fact, the upper bounds are low enough to conclude that the |+ , |− states cannot be heralded no matter what the summation of the PNR pattern is. Let us assume that we can engineer the |0 or |1 state with perfect fidelity from a zero-mean Gaussian state using partial PNR. Then, since it is impossible to engineer the |+ , |− states utilizing a zero-mean Gaussian state, we conclude that any optical implementation of a Hadamard gate (defined in the qubit space) based on Gaussian resources and partial PNR, must necessarily include displacements, in accordance with the setups presented in [52].

VI. FIDELITY UPPER BOUNDS FOR COHERENT GHZ AND CLUSTER STATES
Consider a non-Gaussian target state that is the multi-mode superposition, where b l ∈ C are such that state |C is normalized and |γ (l) is a product of M coherent states |γ , | − γ , or any combination thereof (there exist 2 M such product states). We can rewrite Eq. (33) as, where |γ ≡ |γ ( Working out Eq. (35) and using (34) we get, where n = (n 1 , . . . , n M ). It is hard to write Eq. (36) in closed form, however if one specifies the coefficients b l , the summation is rendered computable. One could write a similar to Eq. (36) expression for a state like (33) but with different coherent amplitudes per mode, however states with equal coherent amplitudes are relevant to quantum computing. We remind the reader that the upper bound u odd is always given by u odd = 1 − u even as per Eq. (12).

A. GHZ states
Let us examine the following GHZ states, where N ± = 2 ± 2e −2M |γ| 2 are the normalization constants. Applying Eq. (36) and (12) we find, which is again consistent with the fact that the mean filed amplitude of GHZ states is zero and with the parity of the PNR patterrn imposed by the absence of displacement in the resource Gaussian state. We note again, that we do not prove that our upper bound is attainable. However, it has been shown that GHZ states can be produced with high fidelity [45] even with zero-mean resource Gaussian states.

B. Coherent Cat-basis Cluster State
Let us move to a more interesting case. Consider the CZ 2 gate whose action is defined as CZ 2 |00 = |00 , CZ 2 |01 = |01 , CZ 2 |10 = |10 , and CZ 2 |11 = −|11 , and therefore is an entangling operation when it acts on | + + . In this work, we denote as CZ any product consisting of multiple two-mode CZ 2 gates, acting on any two qubits of a multi-qubit product state. In fact, we consider that CZ acts on the state |+ ⊗M , i.e., CZ|+ ⊗M , to create entanglement between any possible couple of |+ states at the same time, therefore creating a cluster state on the coherent-cat basis, i.e, is any coherent cat-basis cluster state. Using Eq.  √ n!, we see that any cross-term will be proportional to Therefore, the only non-zero terms in Eq. (44) are of the form |s i n 1 . . . n M |q 1 . . .q M | 2 = | n 1 . . . n M |q 1 . . .q M | 2 , since |s i | 2 = 1. Therefore, we have, From Eq. (10) and following the methods of Section V, we find that, Finally, from Eqs. (45), (46), and (12) we find, Any |CCCS has inherently non-zero mean-field amplitude because the |+ has non-zero displacement. The upper bound of Eqs. (47) and (48) quantifies the damage of not considering displacement as a resource. An upper bound equal to 1/2 on the fidelity with any produced state shows that any |CCCS state is well beyond reach with a zero-mean Gaussian resource state.

VII. CONCLUSIONS AND DISCUSSION
Partial PNR is the new trend for non-Gaussian bosonic state engineering because essentially it circumvents the technical difficulties of constructing non-Gaussian optical unitary operations. However, the are two main drawbacks in said approach: (i) optimization methods are needed to reveal an optimal resource Gaussian state that maximizes the fidelity and probability of occurrence of an acceptable produced state (ii) fidelity is merely a necessary criterion. Any numerical optimization typically does not give intuition on the underlying physics of state engineering. In this paper we asked what would happen if we forbid our resource state to posses any displacements and therefore reveal the implications on coherent-cat basis clusters under any optimization of such resource. We recognized that zero displacement restricts the parity of the observed PNR pattern and therefore it restricts the Fock expansion coefficients (modulo squared) one should sum up to derive a fidelity upper bound, yielding a hard 1/2 upper bound for target states with non-zero mean field amplitude such as the |+ , |− , and |CCCS states.
As a byproduct, we argued that any optical implementation, i.e., based on Gaussian resources and partial PNR, of a qubit Hadamard gate (an operation transforming |0 → |+ ) must necessarily include displacements. This Hadamard gate could be a separate primitive consisting m displaced squeezed states as inputs to a passive n-mode interferometer (m < n). The rest of the n − m input modes could be the output state of another partial PNR based scheme which produces the |0 or |1 states. The Hadamard optical primitive and the |0 or |1 state generator could be combined into a single interferometer, with single mode dispalced squeezed inputs, and an array of PNR detectors at the output, some of which control the production of |0 or |1 and another PNR subset the realization of the Hadamard gate.
It is known that PNR detectors and Gaussian states comprise a universal resource set [28]. Therefore, by working with a general pure Gaussian state, i.e., including displacements, universality must be restored. Apparently, a displacement D(α) on the undetected output would not suffice as it can be easily seen that for example D(α)|0 = |+ . All dis-placements must be applied on the squeezed single mode states going into the passive interferometer, or just before partial PNR (i.e. equivalently partially projecting a zero-mean Gaussian state onto displaced Fock states). However, a constructive way of designing partial PNR based schemes which would attain universality is still elusive.
The holy grail of this line of research would be a systematic theory for non-Gaussian state engineering for specific classes of states that are useful in various quantum information processing tasks such as cluster states for quantum computing, all-optical quantum repeaters, metrologically-optimal states in distributed quantum sensing, etc. One specific interesting question that relates to the states considered in this paper is: Whether the GHZ states considered in this work can be transformed into the CCCS by using local unitaries (e.g., it is known that a star-topology cluster state and a GHZ state are local-Hadamard equivalent), where the local unitaries are themselves realized by post-selected non-Gaussian ancilla states which in turn were heralded using Gaussian states and PNR detectors [52].
Such questions could be answered by expanding the mathematical formalism developed in [45,50] (also given as App. A of this work) to include displacements. This could catalyze further progress toward the non-Gaussian state engineering, if not in providing constructive ways for attaining universality, but at least for constructing optical implementations for specific useful to quantum computation primitives.

ACKNOWLEDGMENTS
CNG was supported by Xanadu Quantum Technologies. CNG and SG acknowledge Xanadu Quantum Technologies for supporting multiple useful discussions on this topic. CNG also acknowledges partial funding support from the National Science Foundation, award number 2122337.

Appendix A: Derivations of probabilities and Fock coefficients of engineered non-Gaussian states
The following can be found as Section II of [50] co-authored by CNG, SK, and collaborators. Said subsection was authored by CNG.
Here we briefly review the results in Ref. [45] and then evolve those to further worked-out formulas. Among other things, in Ref. [45] it was proven that any N -mode pure Gaussian state |Ψ with covariance matrix (CM) V and displacement vector x β can be written in the coherent basis | α as where with Γ = V + I/2, where A = A T , B = B T , and C are defined as the blocks of Γ −1 as follows: Note that we have simplified the expressions compared to Ref. [45]. We note that since the CM V is symmetric, Γ and Γ −1 are also symmetric. We work with the convention = 1 (therefore the CM of vacuum is I/2) and consider the qqpp representation where vectors are defined as x T α = ( q T α , p T α ) with q T α = (q α1 , . . . , q α N ) and p T α = (p α1 , . . . , p α N ) the canonical position and momentum vectors. The volume element for integration is then defined as d 2N x α = dq α1 . . . dq α N dp α1 . . . dp α N , and α i = (q αi + ip αi )/ √ 2. The coherent basis representation is a valuable tool for working on photon-subtraction-based or, more generally, partial PNR detection schemes aimed at engineering Gaussian states into desired non-Gaussian states. Photon subtraction can be modeled either (i) as a beam-splitter whose two input ports are fed with the ith mode of |Ψ and vacuum |0 , respectively, followed by PNR detection on the lower output port; or (ii) simply by acting the annihilation operatorâ i , where the index i refers to the mode, on |Ψ . Therefore, the photon subtraction operator will act only on the basis vectors of the state, i.e., coherent states in this instance. The action of beam-splitters or annihilation operators on coherent states is straightforward, making this basis particularly efficient for analytical or numerical evaluation. The situation is similar for partial PNR detection on a Gaussian state written as a coherent state expansion; the projection of a coherent state on a Fock state is the well known expression n|α = exp(−|α| 2 /2)α n / √ n!. In Ref. [45] it was shown that the probability of a length-N PNR pattern for an N -mode Gaussian state |Ψ with zero displacements, i.e., x β = 0 in Eq. (A1), is given by where In this work, we will derive the explicit relation of the matrix σ to the matrix H −1 and consequently to matrices Γ and V . We also give the expressions for the Fock expansion coefficients of the produced non-Gaussian states and simplify further the expressions. The following subsections summarize new simplifications, observations, and new results which improve on Eqs. (A5-A8).

The determinant and inverse of Γ
The matrix Γ is defined as Γ = V + I/2, where V is the CM and I the identity matrix. Since V corresponds to a pure Gaussian state, it can be written as V = S p V 0 S T p , where S p is an orthogonal symplectic matrix for a general passive transformation (beamsplitters and phase rotations, but not squeezers) and V 0 is the CM for a product of N single mode squeezed vacuum states, i.e., the diagonal matrix where r 1 , . . . , r N are the real and positive squeezing parameters for each of the N single-mode squeezed vacuum states (note that the phase of the squeezing has been absorbed into the orthogonal symplectic transformation S p ).
We have the following relation, In the case where the input squeezing is the same among all single mode squeezed vacuum states, i.e. r 1 = . . . = r N = r, Eq. (A14) reduces to det Γ = cosh 2N r. Now let us simplify Eq. (A5). We can write Γ = S p (V 0 + I/2)S T p , and since S T −1 p = S p is a symplectic orthogonal matrix we have The symplectic orthogonal matrix S p has the following block matrix structure and properties: Moreover, since V 0 is diagonal we can write where T = diag (tanh r 1 , . . . , tanh r N ). In virtue of Eqs. (A16), (A17), and (A19), we find that in Eq. (A5) Therefore, in the most general case possible, Eq. (A5) is simplified to where A and C are given in Eqs. (A23) and (A24), respectively, as functions of the passive symplectic transformation S p and the input squeezing parameters. Consequently, matrix B of Eq. (A3) simplifies to

The determinant and inverse of H
The matrix H appearing in Eq. (A8) is defined as We find it easier if we transform asH = W † HW using the unitary matrix W defined as Utilizing Eqs. (A27), (A28), and (A29) we find from which we see that detH = det I = 1. Since | det W | 2 = 1, we have detH = det H and conclude that det H = 1.
Therefore, Eqs. (A8) and (A15) are further simplified to P n1...n N = I n1...n N condition can be verified in practice by successively increasing the limits and observing no change to P .

Fock expansion coefficients of the produced state
The non-Gaussian state |Φ on the M undetected modes (see Fig. A1), can be written as a partial projection on Fock states of the detected modes: where P is given in Eq. (A38) and |Ψ is the input N -mode Gaussian state. The Fock expansion coefficients of heralded state |Φ are c n1...n M = n 1 . . . n M |Φ . Using Eqs. (A7) and (A39) we find where the numerator is given by Eq. (A7). Therefore, for any given partial PNR pattern (n M +1 , . . . , n N ) one can compute the Fock expansion coefficients of the produced state |Φ , which can be benchmarked against a target non-Gaussian state |Φ t through direct comparison of Fock coefficients or collectively through fidelity F = | Φ t |Φ | 2 .