# Photon catalysis for general multimode multi-photon
quantum state preparation

Andrei Aralov, Émilie Gillet, Viet Nguyen, Andrea Cosentino, Mattia

Walschaers, Massimo Frigerio

► To cite this version:

Andrei Aralov, Émilie Gillet, Viet Nguyen, Andrea Cosentino, Mattia Walschaers, et al. Photon catalysis for general multimode multi-photon quantum state preparation. 2025. [hal-05363643](https://hal.science/hal-05363643v1)

HAL Id: hal-05363643

https://hal.science/hal-05363643v1

Preprint submitted on 13 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Photon catalysis for general multimode multi-photon quantum state preparation

Andrei Aralov,¹, * Émilie Gillet,¹, * Viet Nguyen,² Andrea Cosentino,² Mattia Walschaers,¹, † and Massimo Frigerio¹, ‡

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

²Sorbonne Université, F-75006 Paris, France

(Dated: August 11, 2025)

Multimode multiphoton states are at the center of many photonic quantum technologies, from photonic quantum computing to quantum sensing. In this work, we derive a procedure to generate exactly, and with a controlled number of steps, any such state by using only multiport interferometers, photon number resolving detectors, photon additions and displacements. We achieve this goal by establishing a connection between photonic quantum state engineering and the algebraic problem of symmetric tensor decomposition. This connection allows us to solve the problem by using corresponding results from algebraic geometry and unveils a mechanism of photon catalysis, where photons are injected and subsequently retrieved in measurements, to generate entanglement that cannot be obtained through Gaussian operations. We also introduce a tensor decomposition, that generalizes our method and allows to construct optimal circuits for particular classes of states. As a benchmark, we numerically evaluate our method and compare its performance with state-of-the-art results, confirming 100% fidelity on different classes of states.

# I. INTRODUCTION

Multiphoton states of a multimode system are versatile and widely studied finite linear subspaces of the Fock space of multimode bosonic systems [1] and they find applications in quantum computation [2–10], quantum metrology [11–17] quantum simulation [18–20] as well as being at the core of various types of boson sampling tasks [21–23]. Moreover, these same states have been introduced under the name of core states [24] in the stellar representation of non-Gaussian states and their characterization is therefore essential to the understanding of non-Gaussian states with finite stellar rank [25].

When the total number of photons is constrained, the corresponding subspaces are invariant under all linear optical operations (equivalently, all passive Gaussian operations). However, the vast majority of them cannot be prepared starting from factorized multimode Fock states and passive linear optics, nor, more generally, by alternating probabilistic photon additions and interferometers, as can be easily deduced by counting free parameters. This aspect of multiphoton quantum states has been implicitly known for a long time; standard references on bosonic systems often emphasize the need to consider superpositions of such states to construct the full Fock space [26–28]. The issue is also explicitly discussed in [29]. A first explicit scheme to construct arbitrary multimode multiphoton states dates back several decades [30]. This scheme relies on intense coherent states which are truncated using quantum scissors [31], and intricate Bell-type measurements [32] that must be implemented through postselection. While this approach works in principle, it is highly resource intensive and does

not come with a mathematical framework that allows to reason about its optimality. A more structured classification of multimode states with a definite number of photons was recently proposed in [33], where the authors show that only those belonging to a specific class can be generated by alternating *N* photon additions with passive linear optics transformations acting on the *M* target modes alone, but how to generate states belonging to the other classes remained an open problem, despite the fact that they constitute the vast majority of multimode multiphoton states.

In this paper, we tackle the very general problem of exactly preparing any multimode multiphoton state by considering an ancillary mode on which we perform post-selection through photon counting at the end of the process. In this way, we show that the *M* modes can be projected into *any* desired state with *N* photons, with a quantifiable cost and a fidelity of one. By adding a second ancillary mode, moreover, we can also prepare any state with non-homogeneous photons number, i.e., any multimode multiphoton state.

Our protocol can effectively be implemented with a specific realization of a boson sampling [21] or Gaussian boson sampling [34] device. In order to produce a state with up to *N* photons, the boson sampler that implements our state-engineering protocol must generally be seeded with much more than *N* photons. Even though the additional photons are retrieved by the photon counter, they serve the key purpose of creating a type of mode-intrinsic entanglement [35, 36]. This process has been known in literature under the name of photon catalysis [37]. Hitherto, it was generally applied to produce a range of single-mode states [38–41], but its role in the increase of entanglement has also been investigated [42]. Recently, photon catalysis was studied in a multimode context, but only for two-photon states [43]. Our present work can be seen as a proof by construction that photon catalysis can create any multimode multi-photon state.

<footnote>* These two authors contributed equally.</footnote>

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

<footnote>‡ massimo.frigerio@lkb.upmc.fr</footnote>
The two-photon result of [43] is retrieved as a special case of the proposed protocol.

In the alternative implementation with a Gaussian boson sampler, we inject a series of weakly squeezed states into an interferometer and subsequently count photons. This procedure is commonly used to create single-mode non-Gaussian states for continuous variable quantum computing [44–47]. While possible application for engineering multimode states are discussed in [44], no general protocol was provided. In the context of continuous-variable measurement-based quantum computing [48], it is implicitly known that we can obtain any multimode state by performing measurements on a well-chosen Gaussian state [49]. Our work can be seen as the first explicit demonstration of such a general quantum state-engineering protocol.

The layout of our paper is as follows: after an introduction to the theory of multimode multiphoton states, also known as core states, and their representation as stellar polynomials in section II, we introduce our main result in section III and prove it using Waring decomposition and monic expansion. In section IV, we explain how our protocol could be effectively implemented with common photonic setups: a weak version of a boson sampler and a Gaussian boson sampler. In section V, we provide results bounding the number of photon additions and the rank of the Fock state on which the ancillary mode has to be projected and, in section VI, we analyze some special cases in which the generating cost of some states is shown to be smaller than what is predicted by our fully general method, which is therefore not always optimal. Finally, in section VII, a software implementation of our method and its variants allows us to compare fidelity, probability of success and resource use for a selection of state preparation tasks.

## II. STELLAR REPRESENTATION AND CORE STATES

In this section we introduce the stellar representation of quantum states that will allow us to map our set of quantum operations on core states as operations on complex multivariate polynomials [24, 25]. The stellar function of a given quantum state of $M$ modes $|\psi\rangle \in \mathcal{H}^{\otimes M}$, where $\mathcal{H} \simeq \text{L}^2(\mathbb{R})$ is the Fock space of one bosonic mode, is defined as a renormalized overlap of a multimode coherent state and $|\psi\rangle$:

$$F_{\psi}^{*}(\alpha) = e^{\frac{1}{2}|\alpha|^{2}} \langle \alpha^{*} | \psi \rangle \quad (1)$$

where $\alpha \in \mathbb{C}^M$. It can be shown [25] that any state can be written as:

$$|\psi\rangle = F_{\psi}^{*}(\hat{a}^{\dagger}) |0\rangle. \quad (2)$$

where $\hat{a}^{\dagger} = (\hat{a}_{1}^{\dagger}, \dots, \hat{a}_{M}^{\dagger})$ are the creation operators of the $M$ modes, satisfying the canonical commutation relations $[\hat{a}_{j}, \hat{a}_{k}^{\dagger}] = \delta_{jk}$ with their adjoints.

In particular, for core states, which have a finite support in the Fock basis, the stellar function becomes a polynomial with complex coefficients $P_{\psi} \in \mathbb{C}[x_1, \dots, x_M]$, where from now on we denote by $\mathbb{C}[x_1, \dots, x_M]$ the ring of polynomials in $M$ variables (of arbitrary degree) on the complex field.

$$|\psi\rangle = \sum_{\mathbf{n}} \psi_{\mathbf{n}} |\mathbf{n}\rangle = \sum_{\mathbf{n}} \frac{\psi_{\mathbf{n}}}{\sqrt{\mathbf{n}!}} (\hat{a}^{\dagger})^{\mathbf{n}} |0\rangle = P_{\psi}(\hat{a}^{\dagger}) |0\rangle \quad (3)$$

$$P_{\psi}(\mathbf{x}) = \sum_{\mathbf{n}} p_{\mathbf{n}} \mathbf{x}^{\mathbf{n}}; \quad p_{\mathbf{n}} = \frac{\psi_{\mathbf{n}}}{\sqrt{\mathbf{n}!}} \quad (4)$$

Here and throughout the paper we utilize the usual multi-index notation:

$$\mathbf{n} = (n_1, \dots, n_M) \qquad \mathbf{x} = [x_1, \dots, x_M]^T \quad (5)$$

$$\mathbf{n}! = n_1! \dots n_M! \qquad |\mathbf{n}| = \sum_{i=1}^{M} n_i \quad (6)$$

$$|\mathbf{n}\rangle = |n_1 \dots n_M\rangle \qquad \mathbf{x}^{\mathbf{n}} = x_1^{n_1} \dots x_M^{n_M} \quad (7)$$

This bijective correspondence allows us to study multivariate polynomials $P_{\psi}$ and then transfer results to quantum states (Table I). For simplicity of notation, we will omit normalization factors as they are easily applied by scaling all the coefficients of the polynomial.

<table><caption>TABLE I. Correspondence of state transformations and operations on polynomials</caption><thead><tr><th>State transformation</th><th>Polynomial transformation</th></tr></thead><tbody><tr><td>$\hat{a}_k^\dagger |\psi\rangle$</td><td>$x_k P_\psi(\mathbf{x})$</td></tr><tr><td>$\hat{a}_k |\psi\rangle$</td><td>$\frac{\partial}{\partial x_k} P_\psi(\mathbf{x})$</td></tr><tr><td>$\langle n|_k |\psi\rangle$</td><td>$\frac{\partial^n}{\partial x_k^n} P_\psi(\mathbf{x})|_{x_k=0}$</td></tr></tbody></table>

Another noteworthy reason to study stellar polynomials is that one can classify core states by the number of irreducible factors of different degrees in their factorization. In particular, it is possible to prepare any core state whose stellar polynomial is homogeneous and factorizes into a product of $N$ linear terms directly acting on the $M$ target modes, starting in the vacuum state, by alternating photon additions and multiport interferometers without ancillary modes, as shown in [33]. The class of such states is denoted $[1^N]_M$.

The question remains, however, how and whether one can prepare all other classes of irreducible polynomials with these resources, by adding one ancillary mode to the $M$ target modes and conditioning on specific photon counting outcome of that mode.

## III. PREPARATION OF ARBITRARY MULTIMODE CORE STATES

In this section we present our main contribution, a method to prepare any quantum state with finite support
in the Fock basis using one or two ancillary modes, multi-port interferometers, photon additions, photon-number-resolving (PNR) detectors and displacement operations. We start from the aforementioned observation that any M-modes multiphoton state described by a homogeneous polynomial in $[1^N]_M$ can be prepared exactly with photon additions and multiport interferometers acting solely on the M modes. However, in the general case and unlike single-variable and two-variables polynomials, not all polynomials in $M \ge 3$ variables factorize into a product of linear terms, which leads to a natural question: what is the amount of resources needed to prepare an arbitrary core state? The method that we propose guarantees the preparation of an arbitrary core state of d photons in M modes with 100% fidelity in a finite number of steps, and quantifies the required resources in terms of overhead of photon additions to be performed. The proposed method is described in the following Theorem 1 and depicted in Figure 1.

**Theorem 1.** Any core state containing at most d photons in M modes, described by a polynomial $P_\psi \in \mathbb{C}[x_1, \dots, x_M]$ of degree d in M variables, whose homogenization has Waring rank r can be prepared exactly from the multimode vacuum state, using two ancillary modes in the vacuum, $d \times r$ multiport interferometers and photon additions, one photon projection conditioning on having d ($r-1$) photons, one displacement and one projection conditioning on having zero photons.

The constructive proof of this theorem is detailed in the following subsections and it is based on three main steps. First, we show how any homogeneous polynomial in M variables can be seen as the coefficient of some power of an ancillary variable $\lambda$ inside a polynomial in M + 1 variables. This coefficient, hence the target core state, can be extracted from the corresponding superposition state of M + 1 modes by post-selecting on measuring the right number of photons in the ancillary mode associated with $\lambda$. Secondly, we show that this enlarged polynomial factors as a product of linear terms, therefore the corresponding core state can be prepared with photon additions and interferometers acting solely on the M + 1 modes. Finally, we show that adding a second ancillary mode, displacing it by a fixed amount and projecting on the vacuum leads to the preparation of core states with variable numbers of photons in the superposition.

### A. Waring decomposition

As a first step, we recall the correspondence between homogeneous polynomials and symmetric tensors [50]: for any homogeneous polynomial $P_\psi$ of degree d, there is a corresponding symmetric tensor $P_\psi \in S^d(\mathbb{C}^M)$ encoding its coefficients in a monomial basis. In particular, this correspondence allows to relate the decomposition of a polynomial into a sum of powers of linear forms to

symmetric tensor decomposition. Every symmetric tensor could be decomposed as a sum of tensorial powers of some vectors $\mathbf{w}_k \in \mathbb{C}^M$, where the number of terms in the decomposition is called the symmetric rank of $P_\psi$.

$$P_\psi = \sum_{k=1}^{r} \mathbf{w}_k^{\otimes d} \qquad (8)$$

The corresponding homogeneous polynomial decomposition is called the *Waring decomposition*, and the number of terms r is the *Waring rank* of that polynomial.

$$P_{\psi}(\mathbf{x}) = \sum_{k=1}^{r} (\mathbf{w}_{k}^{T} \mathbf{x})^{d} \qquad (9)$$

### B. Monic expansion and dehomogenization

We now notice that the expression in the Waring decomposition of $P_\psi$ can be trivially rewritten as $e_1((\mathbf{w}_1^T \mathbf{x})^d, \dots, (\mathbf{w}_r^T \mathbf{x})^d)$, where:

$$e_1(\mathbf{x}) = \sum_{j=1}^{M} x_j \qquad (10)$$

is the first symmetric elementary polynomial. Pursuing this idea, we can prove the following result:

**Lemma 1.** Let $P_\psi(\mathbf{x}) \in \mathbb{C}[x_1, \dots, x_M]$ be an homogeneous polynomial of degree d, having Waring decomposition $P_\psi(\mathbf{x}) = \sum_{k=1}^r (\mathbf{w}_k^T \mathbf{x})^d$. Define

$$W = [\mathbf{w}_1, \dots, \mathbf{w}_r]^T \qquad (11)$$

$$L_{r,d,W}(\lambda, \mathbf{x}) = \prod_{k=1}^{r} \left( \lambda^d + (\mathbf{w}_k^T \mathbf{x})^d \right) \qquad (12)$$

Then the coefficient of $\lambda^{d(r-1)}$ in the expansion of $L_{r,d,W}(\lambda, \mathbf{x})$ is exactly $P_\psi(\mathbf{x})$.

*Proof.* We start by recalling the following formula for the so-called *monic expansion*, a direct corollary of Vieta's formula [51]:

$$\prod_{k=1}^{M} (y + x_k) = \sum_{k=0}^{M} y^{M-k} e_k(x_1, \dots, x_M) \qquad (13)$$

$$e_k(\mathbf{x}) = \sum_{1 \le j_1 < \dots < j_k \le M} x_{j_1} \dots x_{j_k} \qquad (14)$$

where $e_k$ is a k-th elementary symmetric polynomial. We can therefore relate the monic expansion to the Waring decomposition by applying the first to:

$$\prod_{k=1}^{r} (\lambda^d + (\mathbf{w}_k^T \mathbf{x})^d) = \\ = \sum_{k=0}^{r} \lambda^{d(r-k)} e_k ((\mathbf{w}_1^T \mathbf{x})^d, \dots, (\mathbf{w}_r^T \mathbf{x})^d) \qquad (15)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Circuit of the proposed method. By injecting $N$ photons, $n$ of which will be later projected-out, a state of stellar rank $d = N - n$ is obtained. The general method described in Theorem 1 sets $N = dr, n = d(r - 1)$; while the special case described in Corollary 3 sets $N = M, n = M - 2$. The blue part represents the preparation of the seed state, the stellar polynomial of which factorizes into a product of linear terms, the green part represents the post-selection that allows to construct an arbitrary homogeneous state and the red part represents dehomogenization.</figcaption></figure>

and noticing that the second term ($k = 1$) in the expansion in Equation 15 contains exactly the Waring decomposition of rank $r$.

$$ \lambda^{d(r-1)} e_1 ((\mathbf{w}_1^T \mathbf{x})^d, \dots, (\mathbf{w}_r^T \mathbf{x})^d) = \\ = \lambda^{d(r-1)} \left( \sum_{i=1}^{r} (\mathbf{w}_i^T \mathbf{x})^d \right) \quad (16) \quad \square $$

To obtain our result for homogeneous polynomials, the last missing step is to show that each factor $(\lambda^d + (\mathbf{w}_k^T \mathbf{x})^d)$ in Equation 15 can be decomposed into a product of linear factors:

**Lemma 2.** Let $\omega_n = e^{\frac{2\pi i}{n}}$ be the primitive $n$-th root of unity. Then we have:

$$ L_{r,d,W} = \prod_{k=1}^{r} \prod_{i=0}^{d-1} (\lambda - \omega_{2d} \omega_d^i \mathbf{w}_k^T \mathbf{x}) \in [1^{dr}]_M \quad (17) $$

We also notice that the above construction can easily be generalized to the preparation of arbitrary, non-homogeneous polynomials from homogeneous ones by a so-called dehomogenization [52].

**Lemma 3.** For any $P \in \mathbb{C}[x_1, \dots, x_M]$ there exists $P^h \in \mathbb{C}[y, x_1, \dots, x_M]$; $P^h$ - homogeneous, such that $\forall \mathbf{x} \in \mathbb{C}^M: P(\mathbf{x}) = P^h(1, \mathbf{x})$. Physically, one can prepare any core state involving superpositions of variable numbers of photons by first preparing the homogeneous state $|\psi^h\rangle$ corresponding to $P^h$, then displacing the ancillary mode used for dehomogenization by 1 and finally conditioning on measuring 0 photons in this ancillary mode: $|\psi\rangle = \langle 0|_0 \hat{D}_0(1) |\psi^h\rangle$.

### C. Proof of Theorem 1

Finally, we can piece together these results to prove our main theorem:

*Proof.* Let $P_{\psi}(\mathbf{x}) = \sum p_n \mathbf{x}^n$ of degree $d$, not necessarily homogeneous, be the stellar polynomial of the target core state $|\psi\rangle$. Denote by $P_{\psi^h}(\mathbf{y})$ its homogenization $(\mathbf{y} = [y, x_1, \dots, x_M]^T \in \mathbb{C}^{M+1})$ and compute its Waring decomposition, yielding the Waring rank $r$ and the vectors $\mathbf{w}_k \in \mathbb{C}^{M+1}$ such that

$$ P_{\psi^h}(\mathbf{y}) = \sum_{k=1}^{r} (\mathbf{w}_k^T \mathbf{y})^d \quad (18) $$

Since, by Lemma 2, with the addition of an ancillary variable $\lambda$ we have that $L_{\psi^h} = L_{r,d,W}(\lambda, \mathbf{y}) \in [1^N]_M$, we can prepare the corresponding core state using the method described in [33] (included here in Appendix C), starting from the multimode vacuum state and alternating photon additions with multiport interferometers. Then, according to Table I, conditioning the ancillary mode $\hat{a}_{\lambda}$ on having $d(r-1)$ photons corresponds to collapsing the sum onto the coefficient next to $\lambda^{d(r-1)}$, which by Equation 16 gives us $P_{\psi^h}$.

Finally, by Lemma 3, applying displacement and vacuum-projection on mode $\hat{a}_y^\dagger$ to $P_{\psi^h}(\mathbf{y})$ produces the dehomogenized polynomial $P(\mathbf{x})$. $\square$

## IV. PHOTONIC IMPLEMENTATION OF THE GENERAL SCHEME

After having outlined the way in which the preparation of a generic multimode core state can be decomposed
<figure><img src="image_3.png" /><figcaption>FIG. 2. Alternative representations of the proposed method, where approximate photon additions are implemented through injection of single photons coupled by weakly reflecting beam-splitters (Figure 2a) as in a boson sampler, or by weak single-mode squeezing operations followed by heralding of single photons (Figure 2b) as in a Gaussian boson sampler.</figcaption></figure>

into an alternating sequence of photon additions and in-
terferometers, we shall see how this construction relates
to other photonic quantum information protocols by con-
sidering in greater detail the physical implementation of
the photon addition steps. As a first remark, notice that
the multiport interferometers ($U_1$ to $U_N$ in Figure 1) ef-
fectively allow us to perform photon addition on arbi-
trary linear combination of modes, and this is their only
purpose in our construction. Moreover, a possible way
to implement a photon addition on a state $|\psi\rangle$ is to let
this state interact with the first Fock state $|1\rangle$ of another
mode through a beam-splitter, whose unitary evolution
is described by $\hat{U}_{BS}(\theta) = e^{\theta(\hat{a}_1^\dagger\hat{a}_2 - \hat{a}_1\hat{a}_2^\dagger)}$, where $\hat{a}_1, \hat{a}_2$ are
the annihilation operators of the modes described by the
state $|\psi\rangle$ and $|1\rangle$, respectively. If $\theta$ is sufficiently small
and if the output two-mode state is conditioned upon
measuring the vacuum on the second mode, the resulting
state will be (see Appendix D):

$$
\langle 0 |_B \hat{U}_{BS}(\theta) | \psi \rangle \otimes |1\rangle = \theta [\hat{a}_1^\dagger |\psi\rangle + \mathcal{O}(\theta^2)] \quad (19)
$$

which is an approximation to a photon-added state. The
probability with which this protocol will work is given by
the squared norm of this output state, $p_0 = \theta^2 (\langle \hat{a}_1^\dagger \hat{a}_1 \rangle_\psi +
1) + \mathcal{O}(\theta^3)$ and it is thus proportional to the average pho-
ton number of the state on which we wish to perform the
photon addition. Therefore, up to an error proportional
to $\theta^2$, we can replace the photon-addition elements $\hat{a}^\dagger$
in Figure 1 by adding an equal number of modes prepared
in the Fock state $|1\rangle\langle 1|$, combine all the inputs into a

single, larger interferometer with $r \times d + M + 2$ input
modes, described by the unitary $\hat{U}$, and then condition-
ing on $r \times d$ output modes to be in the vacuum state.
This rewriting of our procedure, which identifies it as
an instance of boson sampling with post-selection and is
sketched schematically in Figure 2a, makes it clear that
we are effectively injecting $r \times d$ photons in the output
and recovering $d(r-1)$ in the PNR projection on one
ancillary mode, so that exactly $d$ photons are staying
in the final state (before dehomogenization). The ex-
tra photons that have to be injected and then measured
act as *catalyzers* for the preparation of the most generic
multiphoton state, an idea upon which we shall elaborate
further after Corollary 1.

A different scheme to implement photon addition, per-
haps more commonly employed in all-optical platforms,
where single-photon sources might not be readily avail-
able, relies on a two-mode squeezing operation, per-
formed by a nonlinear crystal, and post-selection condi-
tioning on having one photon on the ancillary mode [53].
In this case, the ancillary input mode is in the vacuum
state and no additional photons have to be injected: the
energy is provided by the squeezing operation, through
the pump. This option allows us to rewrite the circuit
as a post-selected Gaussian boson sampling as depicted
in Figure 2b. Notice, in particular, that now we have
the vacuum state of $r \times d + M + 2$ modes followed by a
general multimode Gaussian unitary operation $\hat{G}$. Using
the Bloch-Messiah decomposition [54] and taking into ac-
count that any interferometer whose input modes are in
the vacuum leaves them in the vacuum, we can rewrite $\hat{G}$ acting on $|0\rangle^{\otimes (r \times d + M + 2)}$ as $r \times d$ weak single-mode squeezers acting on the vacuum states of the ancillary modes, followed by a general interferometer and finally by the conditioning projective measurements, on $|1\rangle\langle 1|$ for the $r \times d$ ancillary modes used to implement the approximate addition and on $|n\rangle\langle n|$ for the ancillary mode needed to prepare $|\psi^h\rangle$. In this case, no non-Gaussian resource is used before the measurements, but the number of catalyzing photons is doubled, since the process underlying single-mode squeezing always produces photons in pairs.

# V. ASSESSMENT OF THE EFFICIENCY OF THE METHOD

## A. Bounds on resource use

The number of photon additions to be performed is given by $r \times d$, hence the order of the PNR projection upon which one has to condition to obtain the desired state is $d(r - 1)$: in physical terms, $r \times d$ photons are injected into the interferometer and $d(r - 1)$ of them are recovered when post-selecting on the appropriate Fock state of the ancillary mode, such that the final state contains exactly $d$ photons. Both these quantities scale linearly with the Waring rank $r$ of the associated polynomial. Thus, to assess the cost of generating a state with our method we need to understand how $r$ depends on the number of modes (resp., of variables) $M$ and on the number of photons (resp., degree of the polynomial) $d$.

While the value of $r$ depends on the particular polynomial, the concept of *generic rank* provides a way of reasoning about a general or “typical” polynomial. The value of the generic rank, as a function of the degree $d$ and number of variables $M$, $\bar{r}_{\text{gen}}(M, d)$ is given by the Alexander-Hirschowitz theorem [55]:

$$ \bar{r}_{\text{gen}}(M, d) = \left\lfloor \frac{1}{M} \binom{M+d-1}{d} \right\rfloor \quad (20) $$

The precise statement of the theorem and a geometrical interpretation of the generic rank are presented in the Appendix (see Theorem 3). From a physical point of view, the generic rank can be interpreted as follows:

* A randomly selected homogeneous core state has generic rank with probability 1, but there exist states with ranks that are both higher (super-generic) or lower (subgeneric).

* Any homogeneous core state, even of a supergeneric rank, could be approximated by a homogeneous core state of generic rank with fidelity arbitrarily close to 1.

Asymptotically, for $M \gg d$, we have the following scaling for the number of photon additions:

$$ d \times \bar{r}_{\text{gen}}(M, d) = \frac{M^{d-1}}{(d-1)!} (1 + O(d^2/M)) \quad (21) $$

which is polynomial in $M$. Given the symmetry of the binomial coefficient, for $d \gg M$ one gets a similar behaviour with $d$ and $M$ exchanged, i.e. a polynomial in $d$. Finally, for $M = c \times d \gg 1$, we get:

$$ d \times \bar{r}_{\text{gen}}(M, d) = \frac{a_c}{\sqrt{2\pi d}} (k_c)^d (1 + O(d^{-1})) \quad (22) $$

where $a_c = \frac{1}{c} \sqrt{\frac{c+1}{c}}$ and $k_c = \frac{(c+1)^{c+1}}{c^c} > 1$, which is instead an exponential scaling in $d$, but only in the case where both the number of modes and the number of photons are simultaneously large.

It is worth noting that practically computing the best low rank approximation, or even the exact rank of a given tensor is an NP-hard problem (see section VII). However, the maximal rank of a tensor is always upper bounded by twice its generic rank [56] and tighter upper bounds are known for particular values of $d$ and $M$ [57]. Table II summarizes these results.

<table><caption>TABLE II. Maximal rank and corresponding number of photon additions depending on the number of photons $d$ and number of modes $M$.</caption><thead><tr><th>d</th><th>M</th><th>Maximal Rank</th><th>Photon Additions</th></tr></thead><tbody><tr><td>3</td><td>3</td><td>5</td><td>15</td></tr><tr><td>4</td><td>3</td><td>7</td><td>28</td></tr><tr><td>5</td><td>3</td><td>10</td><td>50</td></tr><tr><td>3</td><td>4</td><td>7</td><td>21</td></tr><tr><td>any</td><td>3</td><td>≤ <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">d<sup>2</sup>+6d+1</div><div style="padding-top: 1px;">4</div></div></td><td>≤ <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">d<sup>3</sup>+6d<sup>2</sup>+d</div><div style="padding-top: 1px;">4</div></div></td></tr><tr><td>any</td><td>any</td><td>≤ 2r&#773;<sub>gen</sub>(M, d)</td><td>≤ 2d&#773;<sub>gen</sub>(M, d)</td></tr></tbody></table>

Given the somewhat unfavorable scaling of the Waring rank and taking into account the numerical simulations, which confirm that the proposed method is not always optimal in terms of the number of photon additions and the order of the PNR projection to be performed on the ancillary mode (see section VII), we are prompted to look for optimal protocols to prepare particular states with less catalyzing photons; in this spirit, we introduce a first step towards estimating the optimal preparation circuit by giving another polynomial decomposition, of which Waring decomposition could be seen as a particular case through Lemma 2; that is, a decomposition in terms of elementary symmetric polynomials of linear forms. The statement of the following theorem is surprisingly equivalent to the universality of the classical algebraic circuit model introduced in [58].

**Theorem 2 (Elementary Symmetric Decomposition).**
There exists $N$ such that for any homogeneous polynomial $P \in \mathbb{C}[x_1, \dots, x_M]$ of degree $d$ there is a matrix
U : $\mathbb{C}^M \to \mathbb{C}^N$ such that the polynomial could be written as an elementary symmetric polynomial of the set of linear forms defined by this matrix:

$$P(\mathbf{x}) = e_d(U\mathbf{x}) \qquad (23)$$

The relation of Theorem 2 to the optimal preparation circuit within the constraints established in this paper could be seen from considering the most general polynomial $L \in [1^N]_{M+1}$ with one ancillary mode $\lambda$.

$$L = \prod_{k=1}^{N} L_k = \prod_{\lambda \in \text{dom } L_k} \prod_{\lambda \notin \text{dom } L_j} L_j \\ \propto \prod_{k=1}^{f_1} (\lambda + \mathbf{u}_k^T \mathbf{x}) \prod_{k=1}^{f_2} (\mathbf{v}_k^T \mathbf{x}) \quad (24)$$

If we suppose that the target polynomial $P_\psi$ is obtained as a coefficient of some power of $\lambda$ (post-selection on a Fock state) and it is irreducible over $\mathbb{C}$, then $f_2 = 0$, $f_1 = N$, therefore from Equation 13 the most general result of a projection on a power of $\lambda$ is an elementary symmetric polynomial $e_d(\mathbf{u}_1^T \mathbf{x}, \dots, \mathbf{u}_N^T \mathbf{x})$. This idea is fully developed and shown to be optimal for the practical case of stellar-rank-two states in Corollary 3.

We also have a straightforward corollary of Theorem 2 which lower bounds the number of photon additions:

**Corollary 1.** If a homogeneous multiphoton state with $d$ photons is described by $M \ge 2$ intrinsic modes, then the minimum number $N$ of photon additions required to prepare it, within the established constraints, is lower bounded by $\max(d, M)$.

Here, the number of intrinsic modes is defined as the minimum number of modes, over all possible mode bases choices, that are not in the vacuum state (see [59]). Notice that states in the class $[1^d]_M$ have at most $\min(d, M)$ intrinsic modes, since one can attach at most one independent mode to each of the $d$ photon additions. More generally, if $M > d$ and the state has $M$ intrinsic modes, Corollary 1 implies that we have to create at least $N \ge M$ photons, $N-d$ of which will have to be retrieved by a post-selection on a PNR measurement outcome; it is thus appropriate to think about our protocol as a multiphoton catalysis mediated by these $N-d$ additional photons which are used only to generate the required number of intrinsic modes, thus to prepare states outside the class $[1^d]_M$.

## B. Probability of success

It is worth observing that the coefficients of $\lambda^{d(r-1)}$ in $L_{r,d,W}$ and in $L_{r,d,\alpha W}$ are identical, up to a scaling factor, and thus that scaling-coefficient $\alpha$ is a freely adjustable parameter, which can impact the probability of success. It can effectively steer the distribution of photon number measurements on the ancillary mode to peak

<figure><img src="image_4.png" /><figcaption>FIG. 3. Probability of successfully detecting the right number of photons to obtain $|\Psi_4\rangle$, $|\Psi_8\rangle$, $|\Psi_9\rangle$ from Table III. The solid lines show the final PNR detection probability for the method described in Theorem 1 (maximized over the obtained decompositions), while the dashed lines show the PNR detection probability of the method described in Corollary 3 where it is applicable.</figcaption></figure>

at $d(r-1)$. Physically, this scaling factor corresponds to the relative weight between the ancilla and the rest of the modes in the state at the output of the addition part $L_{r,d,\alpha W} = \prod (\lambda^d + \alpha^d (\mathbf{w}_k^T \mathbf{x})^d)$, thus changing the interferometers $\hat{U}_k$ in Figure 1. This highlights that the choice of the interferometers in our state preparation protocol is not unique and can, therefore, be optimized. Practically calculating the probability of successfully projecting onto the required amount of photons in the ancillary mode can be done symbolically or numerically. For the simple but nontrivial case where $W \propto I$, the explicit formula can be provided to highlight the nonlinear behavior of such probability.

**Corollary 2.** If $W = wI_{M \times M}$ – rescaled identity matrix, $w \in \mathbb{C}$, then the probability of successfully obtaining $P_\psi$ from $L_\psi$ is given by Equation 25 and plotted in Figure 3.

$$\mathrm{Prob} = \frac{(d(M-1))! d! M w^{2d}}{\sum_{k=0}^{M} (d(M-k))! (d!)^k \binom{M}{k} w^{2dk}} \quad (25)$$

*Proof.* Denote $\|P_\psi\|_1$ – norm of the unnormalized quantum state corresponding to $P_\psi$ (this norm is equal, up to degree normalization, to the Bombieri [60] norm).

$$\|L_\psi\|_1^2 = \left\| \sum_{k=0}^{M} \lambda^{d(M-k)} e_k (wx_1, \dots, wx_M) \right\|_1^2 \quad (26)$$

$$= \sum_{k=0}^{M} (d(M-k))! (d!)^k \binom{M}{k} w^{2dk} \quad (27)$$
$$ \|\lambda^{d(M-1)} e_1(w\mathbf{x})\|_{\|}^2 = (d(M-1))! d! M w^{2d} \quad (28) $$

$$ \text{Prob} = \frac{\|L_{\psi}\|_{\|}^2}{\|\lambda^{d(M-1)} e_1(w\mathbf{x})\|_{\|}^2} \quad (29) $$

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

Such nonlinear behavior of the probability implies that in the most general case, after finding the Waring decomposition of the target polynomial, an additional optimization should be carried out to find the scaling factor $\alpha$ maximizing the probability of success.

## VI. SPECIAL CASES AND EXAMPLES

As will be shown in this section, the method given in Theorem 1 is not always optimal in terms of resources, despite the fact that it is fully general. Some classes of states could be produced within the same framework given by Lemma 2 but with a smaller number of catalyzing photons, using linear combinations of elementary symmetric polynomials, or elementary symmetric polynomials of higher degrees in place of $e_1$. We analyze in the following subsections some particular cases in which this is possible, and show how the bottleneck of the Waring decomposition can be avoided.

### A. Linear combinations of {$e_k$}

Here we present a possible specialization of Theorem 1 for the states that are described by linear combinations of elementary symmetric polynomials. We show that those states could be obtained from Equation 15 by performing a projection onto a superposition of terms on the ancillary mode, which is equivalent to a sequence of displacements and photon subtractions followed by vacuum projection.

**Example 1.** Consider the GHZ state and its homogenization.

$$ |\text{GHZ}\rangle = \frac{|0\rangle^{\otimes M} + |1\rangle^{\otimes M}}{\sqrt{2}} \quad (30) $$

$$ |\text{GHZ}^h\rangle = \frac{|M\rangle|0\rangle^{\otimes M} + |0\rangle|1\rangle^{\otimes M}}{\sqrt{2}} \quad (31) $$

$$ P_{\text{GHZ}^h} = \frac{1}{\sqrt{2}} \left( \frac{1}{\sqrt{M!}} x_0^M + x_1 \dots x_M \right) \quad (32) $$

Note that the homogenization corresponds exactly to the first and last terms in the expansion of Equation 15 for $W = I$, $r = M$, $d = 1$. Thus, by projecting this element of $[1^M]_M$ onto $|\pi\rangle \propto \left(\frac{1}{\sqrt{M!}}|M\rangle + |0\rangle\right)$, we will obtain

$|\text{GHZ}^h\rangle$. To see this explicitly, we first express:

$$ \begin{aligned} L_{M,1,I} &= \prod_{k=1}^{M} (\lambda - x_k) = \sum_{k=0}^{M} \lambda^{M-k} e_k(\mathbf{x}) \\ &= \lambda^M + x_1 \dots x_M + \sum_{k=1}^{M-1} \lambda^{M-k} e_k(\mathbf{x}). \end{aligned} \quad (33) $$

Now, we can easily evaluate the projection on $|\pi\rangle$:

$$ \begin{aligned} \langle \pi | L_{M,1,I} \rangle &\propto \left( \frac{1}{\sqrt{M!}} \langle M | + \langle 0 | \right) \\ &\left( \sqrt{M!} |M\rangle |0\rangle + |0\rangle |1\rangle + \sum_{k=1}^{M-1} \sqrt{k!} |k\rangle (\dots) \right) \\ &\propto |\text{GHZ}^h\rangle \end{aligned} \quad (34) $$

The projection $\langle \pi|$ can be realized as follows. As any single variable polynomial with complex coefficients is fully factorizable, we could write $P_\pi$ as a product of factors:

$$ P_{\pi} \propto x_{0}^{M} + M! = \prod_{j=0}^{M-1} (x_{0} - r_{j}) \quad (35) $$

where $r_j$ are the roots of $P_\pi$. Denoting by $\hat{D}(r) = \exp\{r\hat{a}_0^\dagger - r^*\hat{a}_0\}$ the displacement operator, we rewrite:

$$ |\pi\rangle \propto \prod_{j=0}^{M-1} (\hat{a}_0^\dagger - r_j) |0\rangle = \prod_{j=0}^{M-1} \hat{D}(r_j^*) \hat{a}_0^\dagger \hat{D}^\dagger(r_j^*) |0\rangle \quad (36) $$

Therefore, $|\pi\rangle$ could be prepared using $M$ displacements and $M$ photon additions. By formally reversing the circuit, we find the corresponding $\langle \pi|$, which consists of displacements and photon subtractions followed by a vacuum projection:

$$ \langle \pi | \propto \langle 0 | \hat{D}^{\dagger} (r_{M-1}^{*}) \hat{a}_{0} \hat{D} (r_{M-1}^{*}) \dots \hat{D}^{\dagger} (r_{0}^{*}) \hat{a}_{0} \hat{D} (r_{0}^{*}) \quad (37) $$

This example shows how using two terms of the monic expansion we are able to prepare the $M$-mode homogenized GHZ state using $M$ photon additions, $M$ photon subtractions, $M$ displacements and one vacuum projection. Directly using Theorem 1 would require us to use rank $r = 1 + 2^{M-1}$ (which is the rank of two non-overlapping monomials [61]) and degree $d = M$, and thus inflict a much larger PNR count $M2^{M-1}$. Both methods should be followed by the dehomogenization procedure to obtain a true GHZ state.

### B. Quadratic Polynomials

Our protocol can be simplified significantly in case of two-photon states, i.e., core states of stellar rank two.
This specific setting was already studied in literature [43] using methods based on permanents. In Corollary 3, we show that Theorem II of [43] follows very naturally from our algebraic framework, as a particular case of the decomposition described in our Theorem 2 for two-photon states. We thus find that any such state in M intrinsic modes can be produced by injecting exactly M photons. In the following corollary we explicitly give the coefficients of the corresponding decomposition by analyzing the linear transformations of quadratic forms.

**Corollary 3.** Suppose that the target state is described by $P_\psi(\mathbf{x}) \in \mathbb{C}[x_1, \dots, x_M]$ of degree 2. Then the order of the Fock state projection in the preparation of $|\psi\rangle$ can be reduced to $M' - 2$, where $M'$ is the rank of the quadratic form of the homogeneization of $P_\psi$. Furthermore, for homogeneous polynomials this is the shortest scheme involving the given resources (photon additions, multiport interferometers and PNR detection on one ancillary mode).

**Proof.** The dehomogenization step could be applied similarly to Theorem 1. Thus, suppose without loss of generality that $P_\psi$ is homogeneous with one extra mode, then it follows that it could be written as $P_\psi = \mathbf{x}^T B \mathbf{x}$ for a complex symmetric matrix $B^T = B$. Also, without loss of generality, we assume that $B$ is full rank, otherwise it is possible to apply a linear transformation that maps the variables to the essential variables – a subspace of dimension $M' < M$ such that $B'$ is full rank. We first note that, for any symmetric (possibly complex) matrix $B$ there exists a unitary matrix $U$ and a diagonal matrix $D$ such that $B = U^T DU$ (given by Autonne-Takagi decomposition [62]). By choosing a (generally non unique) complex square root of the diagonal matrix $D$, we can rewrite this factorization as $B = S_B^T S_B$ for some $S_B$. Let $S_A^T S_A$ be such a decomposition for the matrix $A$ that defines $e_2$:

$$e_2(\mathbf{x}) = \mathbf{x}^T A \mathbf{x}, \quad A = \frac{1}{2} \begin{bmatrix} 0 & 1 & 1 & \dots & 1 \\ 1 & 0 & 1 & \dots & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \dots & 1 & 1 \\ 1 & 1 & \dots & 1 & 0 \end{bmatrix} \quad (38)$$

From Equation 15 in the proof of Lemma 1, it follows for $r = M$ and $d = 1$:

$$\prod_{k=1}^{M} (\lambda + \mathbf{w}_k^T \mathbf{x}) = \sum_{k=0}^{M} \lambda^{M-k} e_k (W \mathbf{x}) = \lambda^{M-2} e_2 (W \mathbf{x}) + \cdots \quad (39)$$
other powers of $\lambda$

The coefficient of $\lambda^{M-2}$ can then be rewritten as a quadratic form related by congruence (not necessarily unitary, as the matrix $W$ is arbitrary) to the initial quadratic form of $e_2(\mathbf{x})$:

$$e_2(W\mathbf{x}) = (W\mathbf{x})^T A(W\mathbf{x}) = \mathbf{x}^T W^T A W \mathbf{x} \quad (40)$$

Therefore to prepare the desired state $|\psi\rangle$ it suffices to find $W$ such that

$$\mathbf{x}^T W^T A W \mathbf{x} = \mathbf{x}^T B \mathbf{x} \quad (41)$$

$$\Leftrightarrow W^T A W = (S_A W)^T (S_A W) = B = S_B^T S_B \quad (42)$$

$$\Leftrightarrow S_B = S_A W \quad (43)$$

$$\Leftrightarrow W = S_A^{-1} S_B \quad (44)$$

Note that $A$ is real by Equation 38 and its symmetric part can be diagonalized as $A = T^T DT$ with $D = \text{diag}[1, \frac{-1}{M-1}, \dots, \frac{-1}{M-1}]$, such that $S_A = \text{diag}[1, \frac{i}{\sqrt{M-1}}, \dots, \frac{i}{\sqrt{M-1}}] T$.

Since $S_A$ is invertible, such $W$ always exists and we conclude that any 2-photon core state in $M$ modes can be prepared using one extra mode, $M$ photon additions and multiport interferometers, by conditioning on counting $M-2$ photons in the ancillary mode. The dehomogenization procedure, when needed, just adds a second ancillary mode which is then displaced and post-selected on the vacuum.

The optimality for an arbitrary state follows from the bound $N \ge M$ in Corollary 1. $\square$

Corollary 3 could be extended to $e_k, k \ge 3$, only by increasing the number of linear factors in $L$, which is demonstrated by counting the number of available parameters. Indeed, the number of parameters in $W$ that we can tune is fixed to be $M^2$, but the total number of degree $d$ homogeneous polynomial in $M$ variables is

$$\binom{d+M-1}{d} \quad (45)$$

Already for $d = 3$ Equation 45 gives $\frac{M(M+1)(M+2)}{6} > M^2$ for $M > 2$. This highlights the importance of the decomposition introduced in Theorem 2.

**Example 2.** Consider the quantum state $|\Psi_4\rangle$ mentioned in [33] and the corresponding stellar polynomial.

$$|\Psi_4\rangle = \frac{1}{\sqrt{4}} (|2, 0, 0, 0\rangle + |0, 2, 0, 0\rangle + |0, 0, 2, 0\rangle + |0, 0, 0, 2\rangle) \quad (46)$$

$$P_{\Psi_4}(x_1, x_2, x_3, x_4) = \frac{1}{2\sqrt{2}} (x_1^2 + x_2^2 + x_3^2 + x_4^2) \quad (47)$$

We present two ways of preparing this state, first based on Theorem 1 and second based on Corollary 3. Note that the state is homogeneous so we omit the dehomogenization step and introduce only one ancillary mode.

a. Using $e_1$. Since the polynomial $P_{\Psi_4}$ is already presented as a sum of powers of its variables, we do not need to compute its Waring decomposition, we directly find $W = wI_{4 \times 4}$, where $w \approx 1.45$ according to Figure 3.
Then the seed polynomial $L_{\Psi_4}$ is:

$$L_{\Psi_4} = \prod_{k=1}^{4} (\lambda^2 + (wx_k)^2) \qquad (48)$$

$$= (\lambda + iwx_1)(\lambda - iwx_1)(\lambda + iwx_2)\dots \qquad (49)$$

$$= \sum_{k=0}^{4} \lambda^{2(4-k)} e_k \left( (wx_1)^2, \dots, (wx_4)^2 \right) \qquad (50)$$

$$= \lambda^8 + \lambda^6 w^2 (x_1^2 + x_2^2 + x_3^2 + x_4^2) + \dots \qquad (51)$$

The seed polynomial is factorizable and therefore the seed state can be prepared using photon additions and interferometers. Then, by conditioning on detecting $2(4-1) = 6$ photons on the ancillary mode $\lambda$, we obtain (up to normalization) the desired state $|\Psi_4\rangle$.

b. Using $e_2$. Denote $B$ the matrix corresponding to the quadratic form $P_{\Psi_4}(\mathbf{x}) = \mathbf{x}^T B \mathbf{x}$; $B = I_{4 \times 4}$. Then we can find $W$ such that $W^TAW = B$:

$$W = w \begin{bmatrix} i & \frac{\sqrt{3}i}{3} & \frac{\sqrt{6}i}{6} & \frac{\sqrt{6}}{6} \\ -i & \frac{\sqrt{3}i}{3} & \frac{\sqrt{6}i}{6} & \frac{\sqrt{6}}{6} \\ 0 & -\frac{2\sqrt{3}i}{3} & \frac{\sqrt{6}i}{6} & \frac{\sqrt{6}}{6} \\ 0 & 0 & -\frac{\sqrt{6}i}{2} & \frac{\sqrt{6}}{6} \end{bmatrix} \qquad (52)$$

The normalization factor $w \approx 1.25$ maximizes the success probability according to Figure 3. Finally, we can prepare the state corresponding (after normalization) to the stellar polynomial $L'_{\Psi_4}$ and condition on detecting $4-2=2$ photons on the ancillary mode.

$$L'_{\Psi_4} = \prod_{k=1}^{4} (\lambda - \mathbf{w}_k^T \mathbf{x}) \qquad (53)$$

$$= \sum_{k=0}^{4} \lambda^{4-k} e_k (W \mathbf{x}) \qquad (54)$$

The second method requires less catalyzing photons and achieves better probability of success as it appears to exploit the symmetry of the monic expansion.

## VII. NUMERICAL METHODS

While the general method presented in Theorem 1 allows the preparation of the desired state given its Waring decomposition, practically computing this decomposition can be challenging. In particular, it is known that both determining the rank of a symmetric tensor and computing its best low-rank approximation are NP-hard problems [63]. It is also worth noting that the decomposition is not unique, except for some special cases [64].

The numerous applications of tensor decomposition across different fields have given rise to various decomposition algorithms. Purely algebraic methods, based on properties of the catalectic matrix [57, 65] can provide exact solutions when the decomposition rank is low, and are thus not generally applicable. A variety of numerical

$$\begin{align*} |\Psi\rangle_1 &\propto |200\rangle + |020\rangle + |002\rangle \\ |\Psi\rangle_2 &\propto |300\rangle + |030\rangle + |003\rangle \\ |\Psi\rangle_3 &\propto |400\rangle + |040\rangle + |004\rangle \\ |\Psi\rangle_4 &\propto |2000\rangle + |0200\rangle + |0020\rangle + |0002\rangle \\ |\Psi\rangle_5 &\propto |012\rangle + |120\rangle + |201\rangle + |021\rangle + |102\rangle + |210\rangle \\ |\Psi\rangle_6 &\propto |110\rangle + |101\rangle + |011\rangle \\ |\Psi\rangle_7 &\propto |220\rangle + |202\rangle + |022\rangle \\ |\Psi\rangle_8 &\propto |2000\rangle + |0110\rangle + |0002\rangle \\ |\Psi\rangle_9 &\propto |3000\rangle + |0210\rangle + |0120\rangle + |0003\rangle \\ |\Psi\rangle_{10} &\propto |040\rangle + |121\rangle + |202\rangle \\ |R_2\rangle &\propto |300\rangle + \sqrt{3}|120\rangle + \sqrt{6}|111\rangle + \sqrt{3}|102\rangle \\ |R_4\rangle &\propto |300\rangle + |030\rangle + |003\rangle + |111\rangle \\ |R_5\rangle &\propto |210\rangle + |021\rangle \\ |K_3\rangle &\propto |3000\rangle + |2100\rangle + |2010\rangle + |2001\rangle \\ &\quad - |1110\rangle - |1101\rangle - |1011\rangle - |0111\rangle \end{align*}$$

<caption>TABLE III. States used for evaluation</caption>

algorithms have also been proposed [66, 67], covering different performance trade-off, use cases or implementation platforms.

In this study, the decompositions were performed using the Optax optimization library [68] and its implementation of the Adam gradient descent algorithm. Starting from a target state, the optimization process is as follows. First, the state is converted into the corresponding polynomial $P_{\Psi}$ in creation operators. This polynomial is represented as an order-$d$ symmetric tensor $T$ with $M^d$ complex entries. Assuming the rank $r$ is given, one can find by gradient descent the $r \times M$ matrix $W$ achieving the best approximation of $T$. Since the rank $r$ is not known, we proceed iteratively, increasing the value of $r$ until a near-perfect reconstruction of $T$ is obtained. Then, once the rank and decompositions are found, the seed state $L_{r,d,\alpha W}(\hat{\mathbf{a}}^\dagger)|0\rangle$ is considered, normalized, and its projection on $|d(r-1)\rangle\langle d(r-1)|_0$ determines the probability of success of the preparation scheme, i.e. the probability of having measured $d(r-1)$ photons on the ancillary mode. This probability depends on both $\alpha$ and the decomposition $W$. Given a decomposition $W$, the value of $\alpha$ maximizing the success probability is found by numerical optimization. Finally, in order to estimate how the probability of success can vary with the decomposition $W$, the optimization process is performed $K=25$ times with randomly seeded initial conditions. One should keep in mind that the success probability we report is based on the assumption that photon addition is a deterministic process, and that a PNR detector with ideal characteristics is available.

For the purpose of comparison, we also implemented the scheme described in appendix C of [33], which requires $M+1$ modes, $d+1$ photon additions, and conditioning on detecting exactly one photon on the ancillary
mode. We verified that the fidelity metric found by our
re-implementation matches the values reported by the
original authors, and we additionally computed its prob-
ability of success.

The states used for evaluation are presented in Table III. The states $|\Psi\rangle_1$ to $|\Psi\rangle_7$ are the examples studied in [33]. $|\Psi\rangle_8$, $|\Psi\rangle_9$ and $|\Psi\rangle_{10}$ involve a higher number of modes or photons and pose an additional challenge to preparation. States $|R2\rangle$, $|R4\rangle$ and $|R5\rangle$ are such that the corresponding tensor is respectively of subgeneric ($r=2$), generic ($r=4$), and maximal ($r=5$) rank for the $M=3, d=3$ case, and illustrate how resource use scales with Waring rank. Finally, $|K_3\rangle$ corresponds to a fully connected graph state on 3 vertices.

The results are detailed in Table IV. Since the state preparation method we propose is exact, the fidelity is always 1, with minor deviations attributable to numerical noise. While guaranteeing the preparation of an arbitrary state with 100% fidelity, our method requires more catalyzing photons. The scheme in [33] can attain perfect fidelity for a low number of modes and photons, and we conjecture that indeed it allows the preparation of any three-photon three-mode states but its fidelity degrades when a fourth mode is introduced. With both methods, the optimization process can converge to multiple solutions, which have identically high fidelities (or low reconstruction errors), but which are not equivalent in terms of success probabilities. A practical implementation should thus select the best of several candidate solutions, or include the success probability as an additional term in the function being optimized.

We also compared the scheme in [33] with the method
tailored for quadratic polynomials presented in sec-
tion VIB, for an increasing number of modes *M*. For
each value of *M*, we considered 100 two-photon states
with coefficients sampled randomly from a normal distri-
bution. The results are given in Table V and confirm the
optimality of the preparation scheme presented in Corol-
lary 3. The method we propose requires a single photon
projection for the case of three modes, but then general-
izes to two photon states in any number of modes with
linear scaling of the required target of PNR.

The software implementation of the described methods
is available on GitHub [69].

VIII. CONCLUSIONS AND OUTLOOK

In this work we have proposed a method to prepare
an arbitrary multimode multiphoton state using a fixed
set of operations that are feasible in a quantum optical
experimental setup: multi-port interferometers, photon
additions, photon subtractions and vacuum projections,
or equivalently post-selecting on the outcome of a PNR
detection. This is achieved by utilizing the symmetric
tensor decomposition. While the proposed method is not
optimal for specific classes of states, such as states having
two photons in *M* modes or states that are equivalent to

elementary symmetric polynomials under passive trans-
formations, and while practically computing the decom-
position for very large states might not be feasible, the
generality of the proposed method provides a bound on
the resources required to prepare a given state. In par-
ticular, by importing results from algebraic geometry, we
can provide bounds on the generic and maximum ranks
of symmetric tensor decompositions. Furthermore, con-
sidering that states within experimental reach still have
only a limited number of photons, we are far away from
the practical limitation of our method. It can, therefore,
be directly applied to state-of-the-art experiments.

We conjecture that further refinements of our method,
achieving higher probability of success with less resources
akin to the one provided in Corollary 3 are possible for
higher-order elementary symmetric polynomials; as well
as a generalization of Theorem 1 allowing to put uni-
tary transformations between photon subtractions as in
Example 1. The relation between the rank of the de-
composition and the physical properties of entanglement,
particularly whether it is Gaussian or not, is left to future
research. Further studying the decomposition proposed
in Theorem 2 is a necessary step towards estimating the
minimal number of catalyzing photons required to pre-
pare a generic (in the sense of Theorem 3) state with *d*
photons in *M* modes. Performances might also be evalu-
ated on mixed states under the presence of different types
of noise.

Finally, we emphasize that, in a continuous-variable
setting, any finite-stellar-rank state can be obtained by
performing a well-chosen Gaussian operation on a core
state, which is exactly a multimode multiphoton state.
This implies that our methods is straightforwardly ex-
tended to the generation of any state of finite stellar
rank by adding a Gaussian unitary. While, a priori, this
might be impractical from an experimental point of view,
adding such a Gaussian unitary transformation to the im-
plementations of section IV will give rise to a more feasi-
ble experimental design for continuous-variable quantum
information processing [5, 29].

ACKNOWLEDGMENTS

We acknowledge fruitful discussions with Valentina
Parigi, Carlos E. Lopetegui-González, and Polina Ko-
valenko. A.A. and E.G. are supported by the Euro-
pean Research Council under the Consolidator Grant
COQCOoN (Grant No. 820079). M.W. and M.F. ac-
knowledge financial support from the ANR JCJC project
NoRdiC (ANR-21-CE47-0005). M.W. received financial
support from Plan France 2030 through the project OQu-
Lus (ANR-22-PETQ-0013).
<table><caption>TABLE IV. Resources (number of photon additions, and photon number to be resolved), minimum, median and maximum probability of success, and fidelity of the scheme presented in [33] and the two proposed state preparation methods, for a selection of states. All methods require one ancillary mode.</caption><thead><tr><th rowspan="2">State</th><th rowspan="2">d</th><th rowspan="2">M</th><th rowspan="2">Rank</th><th colspan="4">Baseline</th><th colspan="4">Waring decomposition</th><th colspan="3">e<sub>2</sub> extension</th></tr><tr><th>#Add.</th><th>PNR</th><th>p<sub>min</sub></th><th>p<sub>med</sub></th><th>p<sub>max</sub></th><th>F</th><th>#Add.</th><th>PNR</th><th>p<sub>min</sub></th><th>p<sub>med</sub></th><th>p<sub>max</sub></th><th>F</th><th>#Add.</th><th>PNR</th><th>p</th><th>F</th></tr></thead><tbody><tr><td>|Ψ⟩<sub>1</sub></td><td>2</td><td>3</td><td>3</td><td>3</td><td>1</td><td>0.01</td><td>0.18</td><td>0.28</td><td>1</td><td>6</td><td>4</td><td>0.03</td><td>0.14</td><td>0.26</td><td>1</td><td>3</td><td>1</td><td>0.38</td><td>1</td></tr><tr><td>|Ψ⟩<sub>2</sub></td><td>3</td><td>3</td><td>3</td><td>4</td><td>1</td><td>0.04</td><td>0.09</td><td>0.25</td><td>1</td><td>9</td><td>6</td><td>0.17</td><td>0.17</td><td>0.17</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|Ψ⟩<sub>3</sub></td><td>4</td><td>3</td><td>3</td><td>5</td><td>1</td><td>0.03</td><td>0.14</td><td>0.33</td><td>0.95</td><td>12</td><td>8</td><td>0.09</td><td>0.09</td><td>0.09</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|Ψ⟩<sub>4</sub></td><td>2</td><td>4</td><td>4</td><td>3</td><td>1</td><td>0.06</td><td>0.19</td><td>0.26</td><td>0.75</td><td>8</td><td>6</td><td>0.09</td><td>0.19</td><td>0.28</td><td>1</td><td>4</td><td>2</td><td>0.27</td><td>1</td></tr><tr><td>|Ψ⟩<sub>5</sub></td><td>3</td><td>3</td><td>3</td><td>4</td><td>1</td><td>0.08</td><td>0.14</td><td>0.29</td><td>1</td><td>9</td><td>6</td><td>0.20</td><td>0.20</td><td>0.20</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|Ψ⟩<sub>6</sub></td><td>2</td><td>3</td><td>3</td><td>3</td><td>1</td><td>0.01</td><td>0.14</td><td>0.27</td><td>1</td><td>6</td><td>4</td><td>0.05</td><td>0.17</td><td>0.27</td><td>1</td><td>3</td><td>1</td><td>0.32</td><td>1</td></tr><tr><td>|Ψ⟩<sub>7</sub></td><td>4</td><td>3</td><td>6</td><td>5</td><td>1</td><td>0.01</td><td>0.09</td><td>0.30</td><td>1</td><td>24</td><td>20</td><td>&lt;0.01</td><td>0.01</td><td>0.02</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|Ψ⟩<sub>8</sub></td><td>2</td><td>4</td><td>4</td><td>3</td><td>1</td><td>0.02</td><td>0.15</td><td>0.25</td><td>0.83</td><td>8</td><td>6</td><td>0.07</td><td>0.20</td><td>0.28</td><td>1</td><td>4</td><td>2</td><td>0.14</td><td>1</td></tr><tr><td>|Ψ⟩<sub>9</sub></td><td>3</td><td>4</td><td>4</td><td>4</td><td>1</td><td>0.15</td><td>0.15</td><td>0.15</td><td>0.87</td><td>12</td><td>9</td><td>0.15</td><td>0.15</td><td>0.15</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|Ψ⟩<sub>10</sub></td><td>4</td><td>3</td><td>6</td><td>5</td><td>1</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td><td>0.97</td><td>24</td><td>20</td><td>&lt;0.01</td><td>0.01</td><td>0.02</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|R<sub>2</sub>⟩</td><td>3</td><td>3</td><td>2</td><td>4</td><td>1</td><td>0.05</td><td>0.31</td><td>0.66</td><td>1</td><td>6</td><td>3</td><td>0.13</td><td>0.13</td><td>0.13</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|R<sub>4</sub>⟩</td><td>3</td><td>3</td><td>4</td><td>4</td><td>1</td><td>0.02</td><td>0.18</td><td>0.29</td><td>1</td><td>12</td><td>9</td><td>0.01</td><td>0.06</td><td>0.11</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|R<sub>5</sub>⟩</td><td>3</td><td>3</td><td>5</td><td>4</td><td>1</td><td>&lt;0.01</td><td>0.10</td><td>0.23</td><td>1</td><td>15</td><td>12</td><td>0.01</td><td>0.04</td><td>0.09</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>|K<sub>3</sub>⟩</td><td>3</td><td>4</td><td>5</td><td>4</td><td>1</td><td>0.01</td><td>0.20</td><td>0.26</td><td>0.95</td><td>15</td><td>12</td><td>0.13</td><td>0.13</td><td>0.13</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>TABLE V. Resources (number of photon additions, and photon number to be resolved) and fidelity of the proposed state preparation method using e₂, compared to that of [33] and its extension, for 2-photon states in M modes with coefficients sampled randomly from a normal distribution.</caption><thead><tr><th rowspan="2">d</th><th rowspan="2">M</th><th colspan="4">Baseline</th><th colspan="4">Baseline (extra addition)</th><th colspan="3">e<sub>2</sub> extension</th></tr><tr><th>#Add.</th><th>PNR</th><th>F min</th><th>F avg</th><th>F max</th><th>#Add.</th><th>PNR</th><th>F min</th><th>F avg</th><th>F max</th><th>#Add.</th><th>PNR</th><th>F</th></tr></thead><tbody><tr><td>2</td><td>3</td><td>3</td><td>1</td><td>1</td><td>1</td><td>1</td><td>4</td><td>2</td><td>1</td><td>1</td><td>1</td><td>3</td><td>1</td><td>1</td></tr><tr><td>2</td><td>4</td><td>3</td><td>1</td><td>0.88</td><td>0.97</td><td>1</td><td>4</td><td>2</td><td>1</td><td>1</td><td>1</td><td>4</td><td>2</td><td>1</td></tr><tr><td>2</td><td>5</td><td>3</td><td>1</td><td>0.84</td><td>0.93</td><td>1</td><td>4</td><td>2</td><td>0.94</td><td>0.99</td><td>1</td><td>5</td><td>3</td><td>1</td></tr><tr><td>2</td><td>6</td><td>3</td><td>1</td><td>0.77</td><td>0.87</td><td>0.96</td><td>4</td><td>2</td><td>0.89</td><td>0.96</td><td>1</td><td>6</td><td>4</td><td>1</td></tr></tbody></table>

## Appendix A: Alexander-Hirschowitz theorem

**Theorem 3** (Alexander-Hirschowitz [55]).

Let $\mathcal{Z}_r(M, d) = \{P \in S^d(\mathbb{C}^M) | \text{rank}_S P = r\}$. r is called generic symmetric rank if it is the minimal rank such that $\mathcal{Z}_r(M, d)$ is dense in $S^d(\mathbb{C}^M)$ with respect to Zariski topology on $S^d(\mathbb{C}^M)$ viewed as a vector space over $\mathbb{C}$, i.e.

$$ \bar{Z}_r(M, d) = S^d(\mathbb{C}^M) \quad (A1) $$

For $d > 2$, the generic symmetric rank of $S^d(\mathbb{C}^M)$ is

$$ \bar{r}_{\text{gen}}(M, d) = \left[ \frac{1}{M} \binom{M+d-1}{d} \right] \quad (A2) $$

except for the cases $(M, d) \in \{(3, 5), (4, 3), (4, 4), (4, 5)\}$, where it should be increased by 1.

The geometric and physical interpretations of the generic symmetric rank follow from the following statement about Zariski topology (see [70, Proposition 4.9.5.1])

**Lemma 4.** Any Zariski closed proper subset of a projective space $\mathbb{P}V$ has measure zero with respect to any measure on $\mathbb{P}V$ compatible with its linear structure.

In particular, a tensor not of typical rank has probability zero of being selected at random.

Finally note that a proper set $S \subset V$ of measure zero, according to the measure induced by the Lebsegue one on projective spaces, has no interior points so each point $p \in S$ from such set lies on its border. It follows that any neighborhood $U \ni p$ contains points outside of the set: $U \cap (V \setminus S) \neq \emptyset$, thus it is possible to construct a sequence of points $q_i \in V \setminus S$ converging to $p$.

## Appendix B: Proofs

*Proof of Lemma 2.* The equality follows from the fact that any homogeneous polynomial in two variables is factorizable into a product of linear terms, and particularly

$$ a^d + b^d = \prod_{j=0}^{d-1} (a - \omega_{2d}\omega_d^j b) \quad (B1) $$
since $a = \omega_{2d}\omega_d^j b$ for $j = 0, \dots, d-1$ are the $d$ roots of $a^d + b^d = 0$, with $\omega_d = e^{2\pi i/d}$ is the fundamental $d$-root of unity. It follows that every factor in $L_{r,d,W}$ can be factorized in such way, thus $L_{r,d,W} \in [1^N]_M$. $\square$

**Proof of Lemma 3.** Let $P_\psi(\mathbf{x}) = \sum_{\mathbf{n}} p_{\mathbf{n}} \mathbf{x}^{\mathbf{n}}$ be the stellar polynomial of the desired non-homogeneous state. The homogenization is constructed by adding an ancillary variable to each monomial to complement its degree to the degree of the polynomial, i.e

$$P_{\psi^h}(\lambda, \mathbf{x}) = \sum_{\mathbf{n}} p_{\mathbf{n}} \lambda^{d-|\mathbf{n}|} \mathbf{x}^{\mathbf{n}} \quad (B2)$$

It is then clear that $P_{\psi^h}(1, \mathbf{x}) = P_\psi(\mathbf{x})$.

Now let's consider the action of displacement operator $\hat{D}_0(\alpha)$ on the ancillary mode of $|\psi^h\rangle$.

$$\hat{D}_0(\alpha) |\psi^h\rangle = \hat{D}_0(\alpha) P_{\psi^h}(\hat{a}_0^\dagger, \hat{a}^\dagger) U^\dagger U |00\rangle \quad (B3)$$

$$= \sum_{\mathbf{n}} p_{\mathbf{n}} \hat{D}_0(\alpha) (\hat{a}_0^\dagger)^{d-|\mathbf{n}|} (\hat{a}^\dagger)^{\mathbf{n}} U^\dagger |\alpha\rangle |0\rangle \quad (B4)$$

$$= \sum_{\mathbf{n}} p_{\mathbf{n}} (\hat{a}_0^\dagger + \alpha)^{d-|\mathbf{n}|} (\hat{a}^\dagger)^{\mathbf{n}} |\alpha\rangle |0\rangle \quad (B5)$$

$$= \sum_{\mathbf{n}} p_{\mathbf{n}} \alpha^{d-|\mathbf{n}|} (\hat{a}^\dagger)^{\mathbf{n}} |\alpha\rangle |0\rangle + \hat{a}_0^\dagger(\dots) |\alpha\rangle |0\rangle \quad (B6)$$

After the projection onto the vacuum, the terms that are multiplied by the creation operator on the ancillary mode vanish: $\langle 0|\hat{a}_0^\dagger|\phi\rangle = 0$, thus after the substitution of $\alpha = 1$ we obtain $|\psi\rangle$

$$\langle 0 | \hat{D}_0 (\alpha) | \psi^h \rangle = \sum_{\mathbf{n}} p_{\mathbf{n}} \alpha^{d-|\mathbf{n}|} (\hat{a}^\dagger)^{\mathbf{n}} |0\rangle \langle 0 | \alpha \rangle \quad (B7)$$

$$\propto \sum_{\mathbf{n}} p_{\mathbf{n}} \alpha^{d-|\mathbf{n}|} (\hat{a}^\dagger)^{\mathbf{n}} |0\rangle \quad (B8)$$

$$\stackrel{\alpha=1}{=} \sum_{\mathbf{n}} p_{\mathbf{n}} (\hat{a}^\dagger)^{\mathbf{n}} |0\rangle = |\psi\rangle \quad (B9)$$

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

**Proof of Theorem 2.** The proof proceeds by demonstrating that the Waring decomposition is a special case of the introduced elementary symmetric decomposition. Consider $L_{r,d,W}(\mathbf{x})$ and rewrite Equation 17

$$L_{r,d,W} = \sum_{k=0}^{r} \lambda^{d(r-k)} e_k \left( (\mathbf{w}_1^T \mathbf{x})^d, \dots, (\mathbf{w}_r^T \mathbf{x})^d \right) \quad (B10)$$

$$= \prod_{k=1}^{r} \prod_{i=0}^{d-1} (\lambda - \omega_{2d} \omega_d^i \mathbf{w}_k^T \mathbf{x}) \quad (B11)$$

$$= \prod_{s=1}^{rd} (\lambda + \mathbf{u}_s^T \mathbf{x}) \quad (B12)$$

$$= \sum_{s=0}^{rd} \lambda^{rd-s} e_s (\mathbf{u}_1^T \mathbf{x}, \dots, \mathbf{u}_{rd}^T \mathbf{x}) \quad (B13)$$

As the coefficients in front of different powers of $\lambda$ are unchanged, we obtain the following equality for $k=1$ in Equation B10 and $s=d$ in Equation B13

$$\begin{aligned} e_1 \left( (\mathbf{w}_1^T \mathbf{x})^d, \dots, (\mathbf{w}_r^T \mathbf{x})^d \right) &= \sum_{k=1}^{r} (\mathbf{w}_k^T \mathbf{x})^d \\ &= e_d (\mathbf{u}_1^T \mathbf{x}, \dots, \mathbf{u}_{rd}^T \mathbf{x}) = e_d (U\mathbf{x}) \quad (B14) \end{aligned}$$

Therefore, by choosing $W$ to define the Waring decomposition of the desired homogeneous polynomial $P$, we have constructed the desired $U$ with $N = rd$ rows, which proves the existence of the decomposition.

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

**Proof of Corollary 1.** Suppose that $P(\mathbf{x})$ is an homogeneous polynomial having $M$ intrinsic modes. Any polynomial in $L \in [1^N]_{M+n}$ with $n \ge 0$ ancillary modes can have at most $N$ intrinsic modes, since we can at most attach one intrinsic mode to each linear form in $L$. By conditioning on specific powers of the ancillary variables, we always obtain a linear combination of products of the linear forms appearing in $L$, so each polynomial in $M$ variables that we can obtain from the coefficients of the ancillary variables in $L$ cannot contain more than $N$ intrinsic modes. It thus follows that if the target polynomial has $M$ intrinsic modes, we must have $N \ge M$. Moreover, it is also clear that the degree $d$ of $P$ is a lower bound to the number of linear factors $N$ in $L$, if $L$ has to generate $P$. Thus $N \ge d$ and, altogether, we have $N \ge \max(d, M)$.

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

## Appendix C: Construction of the sequence of interferometers

Here we give the algorithm, first presented in [33], that could be used to obtain a sequence of unitaries $\{U_i\}$ given a sequence of linear forms in creation operators $\{\mathbf{w}_i^T \hat{a}^\dagger\}_{i=1..N}$, such that alternating photon additions on the ancillary mode $\hat{a}_0^\dagger$ with these unitaries produces the same state as $\prod_{i=1}^N \mathbf{w}_i^T \mathbf{a}^\dagger$. In other words, this is the algorithm that allows to construct a circuit that prepares the seed state corresponding to the fully factorizable stellar polynomial $L_{r,d,W}$ (depicted as the blue part in Figure 1).

<table><caption>**Algorithm 1:** Product of linear forms to unitaries</caption><tr><td colspan="2">for i from 1 to N do</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;U<sub>N-i-1</sub> &larr; UnitaryCompletion(w<sub>i</sub>)</td></tr><tr><td colspan="2">for j from i + 1 to N do</td></tr><tr><td colspan="2">&nbsp;&nbsp;&nbsp;&nbsp;| w<sub>j</sub> &larr; w<sub>j</sub>U<sub>N-i-1</sub><sup>T</sup></td></tr><tr><td colspan="2">end</td></tr><tr><td colspan="2">end</td></tr></table>
Algorithm 2: Unitary completion

Input: A vector $w \in \mathbb{C}^n$

Output: A unitary matrix $U \in \mathbb{C}^{n \times n}$ such that
$U^\dagger e_0 = w_i / ||w_i||$

$U \leftarrow 0_{n \times n} \text{ // Initialize zero matrix}$

$S \leftarrow \{i | w_i \neq 0\} \text{ // Support of } w$

$n_s \leftarrow |S|$

$V \leftarrow \mathbb{I}_{n_s} \text{ // Identity matrix of size } n_s$

for $i$ from 0 to $n_s - 1$ do
| $v_{0,i} \leftarrow w_{S[i]}$

end

$V \leftarrow \text{Gram-Schmidt}(V)$

$C \leftarrow S \text{ // Affected columns}$

$R \leftarrow [0] \cup S[1:] \text{ // Affected rows}$

for $i$ from 0 to $n_s - 1$ do

| for $j$ from 0 to $n_s - 1$ do
| | $u_{R[i],C[j]} \leftarrow v_{i,j}$

end

end

for $j$ from 1 to $n - 1$ do

| if $j \notin S$ then
| | $u_{j,j} \leftarrow 1$

end

end

if $0 \notin C$ then
| $u_{S[0,0]} \leftarrow 1$

end

return $U$

Appendix D: Asymptotic error of photon addition

a. Beam splitter. Consider $\hat{U}_{BS} = e^{\theta(\hat{a}_1^\dagger \hat{a}_2 - \hat{a}_1 \hat{a}_2^\dagger)}$.

$$ \hat{U}_{\text{BS}} = 1 + \theta (\hat{a}_1^\dagger \hat{a}_2 - \hat{a}_1 \hat{a}_2^\dagger) + \frac{\theta^2}{2} (\hat{a}_1^\dagger \hat{a}_2 - \hat{a}_1 \hat{a}_2^\dagger)^2 + \mathcal{O}(\theta^3) \quad (\text{D1}) $$

$$ \langle 0 |_2 \hat{U}_{\text{BS}} | \psi \rangle_1 | 1 \rangle_2 = \sum_{n=0} \psi_n \langle 0 |_2 \hat{U}_{\text{BS}} | n \rangle_1 | 1 \rangle_2 \quad (\text{D2}) $$

Notice that $\langle 0 |_2 (\hat{a}_1^\dagger \hat{a}_2 - \hat{a}_1 \hat{a}_2^\dagger)^2 | n \rangle_1 | 1 \rangle_2 = 0$ so

$$ \begin{aligned} \langle 0 |_2 \hat{U}_{\text{BS}} | \psi \rangle_1 | 1 \rangle_2 &= \sum_{n=0} \psi_n \langle 0 |_2 (\theta (\hat{a}_1^\dagger \hat{a}_2 - \hat{a}_1 \hat{a}_2^\dagger) + \mathcal{O}(\theta^3)) | n \rangle_1 | 1 \rangle_2 \\ &= \theta \hat{a}_1^\dagger | \psi \rangle + \mathcal{O}(\theta^3) = \theta [\hat{a}_1^\dagger | \psi \rangle + \mathcal{O}(\theta^2)] \end{aligned} \quad (\text{D3}) $$

b. Squeezer. Consider $\hat{U}_{PDC} = e^{\xi(\hat{a}_1^\dagger \hat{a}_2 - \hat{a}_1 \hat{a}_2)}$.

$$ \hat{U}_{\text{PDC}} = 1 + \xi (\hat{a}_1^\dagger \hat{a}_2^\dagger - \hat{a}_1 \hat{a}_2^\dagger) + \frac{\xi^2}{2} (\hat{a}_1^\dagger \hat{a}_2^\dagger - \hat{a}_1 \hat{a}_2^\dagger)^2 + \mathcal{O}(\xi^3) \quad (\text{D4}) $$

$$ \begin{aligned} \langle 1 |_2 \hat{U}_{\text{PDC}} | \psi \rangle_1 | 0 \rangle_2 &= \sum_{n=0} \psi_n \langle 1 |_2 \hat{U}_{\text{PDC}} | n \rangle_1 | 0 \rangle_2 \quad (\text{D5}) \\ \text{Since } \langle 1 |_2 (\hat{a}_1^\dagger \hat{a}_2^\dagger - \hat{a}_1 \hat{a}_2^\dagger)^2 | n \rangle_1 | 0 \rangle_2 = 0 \text{ as well} \end{aligned} $$

$$ \begin{aligned} \langle 1 |_2 \hat{U}_{\text{PDC}} | \psi \rangle_1 | 0 \rangle_2 &= \sum_{n=0} \psi_n \langle 1 |_2 (\xi (\hat{a}_1^\dagger \hat{a}_2^\dagger - \hat{a}_1 \hat{a}_2^\dagger) + \mathcal{O}(\xi^3)) | n \rangle_1 | 0 \rangle_2 \\ &\propto \xi \hat{a}_1^\dagger | \psi \rangle + \mathcal{O}(\xi^3) = \xi [\hat{a}_1^\dagger | \psi \rangle + \mathcal{O}(\xi^2)] \end{aligned} \quad (\text{D6}) $$

[1] J.-W. Pan, Z.-B. Chen, C.-Y. Lu, H. Weinfurter, A. Zeilinger, and M. Żukowski, Multiphoton entanglement and interferometry, Reviews of Modern Physics **84**, 777-838 (2012).

[2] E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, Nature **409**, 46 (2001).

[3] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Linear optical quantum computing with photonic qubits, Rev. Mod. Phys. **79**, 135 (2007).

[4] F. Flamini, N. Spagnolo, and F. Sciarrino, Photonic quantum information processing: a review, Reports on Progress in Physics **82**, 016001 (2018).

[5] U. Chabaud and M. Walschaers, Resources for bosonic quantum computational advantage, Phys. Rev. Lett. **130**, 090602 (2023).

[6] S. Bartolucci, P. Birchall, H. Bombín, H. Cable, C. Dawson, M. Gimeno-Segovia, E. Johnston, K. Kieling, N. Nickerson, M. Pant, F. Pastawski, T. Rudolph, and C. Sparrow, Fusion-based quantum computation, Nature Communications **14**, 912 (2023).

[7] N. Maring, A. Fyrillas, M. Pont, E. Ivanov, P. Stepanov, N. Margaria, W. Hease, A. Pishchagin, A. Lemaître, I. Sagnes, T. H. Au, S. Boissier, E. Bertasi, A. Baert, M. Valdivia, M. Billard, O. Acar, A. Brieussel, R. Mezher, S. C. Wein, A. Salavrakos, P. Sinnott, D. A. Fioretto, P.-E. Emeriau, N. Belabas, S. Mansfield, P. Senellart, J. Senellart, and N. Somaschi, A versatile single-photon-based quantum computing platform, Nature Photonics **18**, 603 (2024).

[8] L. Monbroussou, E. Z. Mamon, H. Thomas, V. Yacoub, U. Chabaud, and E. Kashefi, Toward quantum advantage with photonic state injection, Phys. Rev. Res. **7**, 033051 (2025).

[9] K. Alexander, A. Benyamini, D. Black, D. Bonneau, S. Burgos, B. Burridge, H. Cable, G. Campbell, G. Catalano, A. Ceballos, C.-M. Chang, S. S. Choudhury, C. J. Chung, F. Danesh, T. Dauer, M. Davis, E. Dudley, P. Er-Xuan, J. Fargas, A. Farsi, C. Fenrich, J. Frazer,
M. Fukami, Y. Ganesan, G. Gibson, M. Gimeno-Segovia, S. Goeldi, P. Goley, R. Haislmaier, S. Halimi, P. Hansen, S. Hardy, J. Horng, M. House, H. Hu, M. Jadidi, V. Jain, H. Johansson, T. Jones, V. Kamineni, N. Kelez, R. Koustuban, G. Kovall, P. Krogen, N. Kumar, Y. Liang, N. Li-Causi, D. Llewellyn, K. Lokovic, M. Lovelady, V. R. Manfrinato, A. Melnichuk, G. Mendoza, B. Moores, S. Mukherjee, J. Munns, F.-X. Musalem, F. Najafi, J. L. O'Brien, J. E. Ortmann, S. Pai, B. Park, H.-T. Peng, N. Penthorn, B. Peterson, G. Peterson, M. Poush, G. J. Pryde, T. Ramprasad, G. Ray, A. V. Rodriguez, B. Roxworthy, T. Rudolph, D. J. Saunders, P. Shadbolt, D. Shah, A. Bahgat Shehata, H. Shin, J. Sinsky, J. Smith, B. Sohn, Y.-I. Sohn, G. Son, M. C. M. M. Souza, C. Sparrow, M. Staffaroni, C. Stavrakas, V. Sukumaran, D. Tamborini, M. G. Thompson, K. Tran, M. Triplett, M. Tung, A. Veitia, A. Vert, M. D. Vidrighin, I. Vorobeichik, P. Weigel, M. Wingert, J. Wooding, X. Zhou, and Psi-Quantum team, A manufacturable platform for photonic quantum computing, Nature **641**, 876 (2025).

[10] A. Salavrakos, N. Maring, P.-E. Emeriau, and S. Mansfield, Photon-native quantum algorithms, Materials for Quantum Technology **5**, 023001 (2025).

[11] J. P. Dowling, Quantum optical metrology – the lowdown on high-n00n states, Contemporary Physics **49**, 125–143 (2008).

[12] U. Dorner, R. Demkowicz-Dobrzanski, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Optimal quantum phase estimation, Phys. Rev. Lett. **102**, 040403 (2009).

[13] G. Y. Xiang, H. F. Hofmann, and G. J. Pryde, Optimal multi-photon phase sensing with a single interference fringe, Scientific Reports **3**, 10.1038/srep02684 (2013).

[14] P. C. Humphreys, M. Barbieri, A. Datta, and I. A. Walmsley, Quantum enhanced multiple phase estimation, Phys. Rev. Lett. **111**, 070403 (2013).

[15] L. Zhang and K. W. C. Chan, Quantum multiparameter estimation with generalized balanced multimode non-like states, Phys. Rev. A **95**, 032321 (2017).

[16] S. Hong, J. ur Rehman, Y.-S. Kim, Y.-W. Cho, S.-W. Lee, H. Jung, S. Moon, S.-W. Han, and H.-T. Lim, Quantum enhanced multiple-phase estimation with multi-mode n00n states, Nature Communications **12**, 5211 (2021).

[17] M. Namkung, D.-H. Kim, S. Hong, Y.-S. Kim, C. Lee, and H.-T. Lim, Optimal multiple-phase estimation with multi-mode noon states against photon loss, New Journal of Physics **26**, 073028 (2024).

[18] I. Dhand, M. Engelkemeier, L. Sansoni, S. Barkhofen, C. Silberhorn, and M. B. Plenio, Proposal for quantum simulation via all-optically-generated tensor network states, Phys. Rev. Lett. **120**, 130501 (2018).

[19] B. Bartlett, O. Y. Long, A. Dutt, and S. Fan, Programmable photonic system for quantum simulation in arbitrary topologies, APL Quantum **1**, 016102 (2024).

[20] F. H. B. Somhorst, R. van der Meer, M. Correa Anguita, R. Schadow, H. J. Snijders, M. de Goede, B. Kassenberg, P. Venderbosch, C. Taballione, J. P. Epping, H. H. van den Vlekkert, J. Timmerhuis, J. F. F. Bulmer, J. Lugani, I. A. Walmsley, P. W. H. Pinkse, J. Eisert, N. Walk, and J. J. Renema, Quantum simulation of thermodynamics in an integrated quantum photonic processor, Nature Communications **14**, 3895 (2023).

[21] S. Aaronson and A. Arkhipov, The computational com-

plexity of linear optics, in *Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing*, STOC '11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 333–342.

[22] L. Chakhmakhchyan and N. J. Cerf, Boson sampling with gaussian measurements, Phys. Rev. A **96**, 032326 (2017).

[23] N. Spagnolo, D. J. Brod, E. F. Galvão, and F. Sciarrino, Non-linear boson sampling, npj Quantum Information **9**, 10.1038/s41534-023-00676-x (2023).

[24] U. Chabaud and S. Mehraban, Holomorphic representation of quantum computations, Quantum **6**, 831 (2022).

[25] U. Chabaud, D. Markham, and F. Grosshans, Stellar representation of non-gaussian quantum states, Phys. Rev. Lett. **124**, 063605 (2020).

[26] O. Bratteli and D. W. Robinson, *Operator algebras and quantum statistical mechanics*, 2nd ed., Theoretical and Mathematical Physics (Springer, Berlin, Germany, 1997).

[27] A. F. Verbeure, *Many-body boson systems*, Theoretical and mathematical physics (Springer, London, England, 2010).

[28] M. Walschaers, Non-gaussian quantum states and where to find them, PRX Quantum **2**, 030204 (2021).

[29] U. Chabaud, G. Ferrini, F. Grosshans, and D. Markham, Classical simulation of gaussian quantum circuits with non-gaussian input states, Phys. Rev. Res. **3**, 033018 (2021).

[30] J. Fiurášek, S. Massar, and N. J. Cerf, Conditional generation of arbitrary multimode entangled states of light with linear optics, Phys. Rev. A **68**, 042325 (2003).

[31] D. T. Pegg, L. S. Phillips, and S. M. Barnett, Optical state truncation by projection synthesis, Phys. Rev. Lett. **81**, 1604 (1998).

[32] M. Dušek, Discrimination of the bell states of qudits by means of linear optics, Optics Communications **199**, 161 (2001).

[33] D. A. Kopylov, C. Offen, L. Ares, B. Wembe, S. Ober-Blöbaum, T. Meier, P. R. Sharapova, and J. Sperling, Multiphoton, multimode state classification for nonlinear optical circuits, Phys. Rev. Res. **7**, 033062 (2025).

[34] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian boson sampling, Phys. Rev. Lett. **119**, 170501 (2017).

[35] J. Sperling, A. Perez-Leija, K. Busch, and C. Silberhorn, Mode-independent quantum entanglement for light, Phys. Rev. A **100**, 062129 (2019).

[36] C. E. Lopetegui, M. Isoard, N. Treps, and M. Walschaers, Detection of mode-intrinsic quantum entanglement, Optica Quantum **3**, 312 (2025).

[37] A. I. Lvovsky and J. Mlynek, Quantum-optical catalysis: Generating nonclassical states of light by means of linear optics, Phys. Rev. Lett. **88**, 250401 (2002).

[38] T. J. Bartley, G. Donati, J. B. Spring, X.-M. Jin, M. Barbieri, A. Datta, B. J. Smith, and I. A. Walmsley, Multiphoton state engineering by heralded interference between single photons and coherent states, Phys. Rev. A **86**, 043820 (2012).

[39] L.-Y. Hu, J.-N. Wu, Z. Liao, and M. S. Zubairy, Multiphoton catalysis with coherent state input: nonclassicality and decoherence, Journal of Physics B: Atomic, Molecular and Optical Physics **49**, 175504 (2016).

[40] R. J. Birrittella, M. E. Baz, and C. C. Gerry, Photon catalysis and quantum state engineering, J. Opt. Soc. Am. B **35**, 1514 (2018).

[41] M. Eaton, R. Nehra, and O. Pfister, Non-gaussian and
gottesman-kitaev-preskill state preparation by photon catalysis, New Journal of Physics **21**, 113034 (2019).

[42] T. J. Bartley and I. A. Walmsley, Directly comparing entanglement-enhancing non-gaussian operations, New Journal of Physics **17**, 023038 (2015).

[43] G. De Gliniasty, P. Bagourd, S. Draux, and B. Bourdon- cle, Simple rules for two-photon state preparation with linear optics, in *2024 IEEE International Conference on Quantum Computing and Engineering (QCE)*, Vol. 01 (2024) pp. 706–711.

[44] D. Su, C. R. Myers, and K. K. Sabapathy, Conversion of gaussian states to non-gaussian states using photon-number-resolving detectors, Phys. Rev. A **100**, 052301 (2019).

[45] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, K. K. Sabapathy, N. C. Menicucci, and I. Dhand, Blueprint for a Scalable Photonic Fault- Tolerant Quantum Computer, Quantum **5**, 392 (2021).

[46] K. Fukui, S. Takeda, M. Endo, W. Asavanant, J.-i. Yoshikawa, P. van Loock, and A. Furusawa, Efficient backcasting search for optical quantum state synthesis, Phys. Rev. Lett. **128**, 240503 (2022).

[47] M. V. Larsen, J. E. Bourassa, S. Kocsis, J. F. Tasker, R. S. Chadwick, C. González-Arciniegas, J. Hastrup, C. E. Lopetegui-González, F. M. Miatto, A. Motamedi, R. Noro, G. Roeland, R. Baby, H. Chen, P. Contu, I. Di Luch, C. Drago, M. Giesbrecht, T. Grainge, I. Krasnokutska, M. Menotti, B. Morrison, C. Puviraj, K. Rezaei Shad, B. Hussain, J. McMahon, J. E. Ort- mann, M. J. Collins, C. Ma, D. S. Phillips, M. Sey- mour, Q. Y. Tang, B. Yang, Z. Vernon, R. N. Alexander, and D. H. Mahler, Integrated photonic source of Gottes- man-Kitaev-Preskill qubits, Nature **642**, 587 (2025).

[48] N. C. Menicucci, P. van Loock, M. Gu, C. Weedbrook, T. C. Ralph, and M. A. Nielsen, Universal quantum com- putation with continuous-variable cluster states, Phys. Rev. Lett. **97**, 110501 (2006).

[49] F. Arzani, R. I. Booth, and U. Chabaud, Can effective descriptions of bosonic systems be considered complete? (2025), arXiv:2501.13857 [quant-ph].

[50] J. M. Landsberg, Multilinear algebra, in [71], Chap. 2, pp. 42–43.

[51] D. A. Cox, J. Little, and D. O’Shea, Invariant theory of finite groups, in *Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra* (Springer International Publishing, Cham, 2015) pp. 346–347.

[52] D. A. Cox, J. Little, and D. O’Shea, Projective alge- braic geometry, in *Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra* (Springer International Publishing, Cham, 2015) pp. 400–401.

[53] G. Roeland, S. Kaali, V. Roman Rodriguez, N. Treps, and V. Parigi, Mode-selective single-photon addition to a multimode quantum field, New Journal of Physics **24**,

043031 (2022).

[54] S. L. Braunstein, Squeezing as an irreducible resource, Phys. Rev. A **71**, 055801 (2005).

[55] J. M. Landsberg, Secant varieties, in [71], Chap. 5, pp. 125–126.

[56] G. Blekherman and Z. Teitler, On maximum, typical and generic ranks, Math. Ann. **362**, 1021 (2015).

[57] A. Bernardi, E. Carlini, M. V. Catalisano, A. Gimigliano, and A. Oneto, The hitchhiker guide to: Secant varieties and tensor decomposition, Mathematics **6**, 314 (2018).

[58] A. Shpilka, Affine projections of symmetric polynomials, Journal of Computer and System Sciences **65**, 639 (2002), special Issue on Complexity 2001.

[59] C. Fabre and N. Treps, Modes and states in quantum optics, Rev. Mod. Phys. **92**, 035005 (2020).

[60] B. Beauzamy, E. Bombieri, P. Enflo, and H. L. Mont- gomery, Products of polynomials in many variables, Jour- nal of Number Theory **36**, 219 (1990).

[61] K. Ranestad and F.-O. Schreyer, On the rank of a sym- metric form, Journal of Algebra **346**, 340 (2011).

[62] R. A. Horn and C. R. Johnson, Unitary congruence and complex symmetric matrices, in *Matrix Analysis* (Cam- bridge University Press, 2012) pp. 263–264, 2nd ed.

[63] C. J. Hillar and L.-H. Lim, Most tensor problems are np-hard, J. ACM **60**, 10.1145/2512329 (2013).

[64] F. Galuppi and M. Mella, Identifiability of homogeneous polynomials and cremona transformations, Journal für die reine und angewandte Mathematik (Crelles Journal) **2019**, 279 (2019).

[65] L. Oeding and G. Ottaviani, Eigenvectors of tensors and algorithms for Waring decomposition, Journal of Sym- bolic Computation **54**, 9–35 (2013).

[66] J. Kileel and J. M. Pereira, Subspace power method for symmetric tensor decomposition (2025), arXiv:1912.04007 [math.NA].

[67] J. Nie, Low rank symmetric tensor approximations, SIAM Journal on Matrix Analysis and Applications **38**, 1517 (2017), https://doi.org/10.1137/16M1107528.

[68] DeepMind, I. Babuschkin, K. Baumli, A. Bell, S. Bhu- patiraju, J. Bruce, P. Buchlovsky, D. Budden, T. Cai, A. Clark, I. Danihelka, A. Dedieu, C. Fantacci, J. God- win, C. Jones, R. Hemsley, T. Hennigan, M. Hes- sel, S. Hou, S. Kapturowski, T. Keck, I. Kemaev, M. King, M. Kunesch, L. Martens, H. Merzic, V. Miku- lik, T. Norman, G. Papamakarios, J. Quan, R. Ring, F. Ruiz, A. Sanchez, L. Sartran, R. Schneider, E. Sezener, S. Spencer, S. Srinivasan, M. Stanojević, W. Stokowiec, L. Wang, G. Zhou, and F. Viola, The DeepMind JAX Ecosystem (2020).

[69] https://github.com/EQ15T/photon-catalysis.

[70] J. M. Landsberg, Algebraic geometry for spaces of ten- sors, in [71], Chap. 4, pp. 114–115.

[71] J. M. Landsberg, *Tensors: geometry and applications*, Graduate studies in mathematics, Vol. 128 (American Mathematical Society, Providence, R.I., 2012).