Non-linear Boson Sampling

Boson Sampling is a task that is conjectured to be computationally hard for a classical computer, but which can be efficiently solved by linear-optical interferometers with Fock state inputs. Significant advances have been reported in the last few years, with demonstrations of small- and medium-scale devices, as well as implementations of variants such as Gaussian Boson Sampling. Besides the relevance of this class of computational models in the quest for unambiguous experimental demonstrations of quantum advantage, recent results have also proposed first applications for hybrid quantum computing. Here, we introduce the adoption of non-linear photon-photon interactions in the Boson Sampling framework, and analyze the enhancement in complexity via an explicit linear-optical simulation scheme. By extending the computational expressivity of Boson Sampling, the introduction of non-linearities promises to disclose novel functionalities for this class of quantum devices. Hence, our results are expected to lead to new applications of near-term, restricted photonic quantum computers.

Introduction. -Quantum technologies promise to provide speed-up in several fields, ranging from intrinsically secure long-distance quantum communication [1] to a novel generation of high-precision sensors [2], and enhanced computational and simulation capabilities [3]. Among the currently developed experimental platforms, in the last few years photonic technologies have recently experienced a technological boost in all fundamental components, namely photon sources, manipulation and detection [4].
Recent studies have focused on identifying suitable dedicated classically-hard tasks, with the aim of leveraging the necessary technological resources and system size to reach the quantum advantage regime [5,6]. Such regime corresponds to the scenario where a quantum device solves a given task faster than any classical counterpart. Within this context, a computational problem named Boson Sampling [7] has been defined as a promising approach. This problem, that consists in sampling from the output distribution of a system of n non-interacting bosons undergoing linear evolution, is a classically-hard task (in n) while it can be naturally solved by a linear optical photonic system. Such sampling problem has also subsequently inspired other classes of sampling problems [8,9] suitable to be solved with different quantum hardware [5,10,11].
In this Letter, we introduce the adoption of non-linear interactions at the few-photon level within the Boson Sampling framework as a route to increase the complexity and reduce the threshold for the quantum advantage regime. This possibility is encouraged by recent advances showing the first experimental demonstrations of non-linear photonic processes within solid state devices [61]. We will first describe the introduction of nonlinearities within the otherwise linear evolution. Then, we will provide an upper bound on the complexity of the enhanced devices via a simulation scheme based on auxiliary, linear-optical gadgets. We will discuss both the asymptotic and finite cases, leveraging results from the well-established linear Boson Sampling framework [7]. the output distribution of n indistinguishable, noninteracting photons after evolution through a m-mode linear network [see Fig. 1(a)]. Given an interferometer described by a unitary matrix U , the transition amplitude from input |S to output state |T can be written as: where per(A) = σ∈Sn n i=1 a i,σ(i) is the matrix permanent, {s i }({t i }) are the occupation numbers of states |S (|T ), and U S,T is the n×n matrix obtained by selecting rows and columns of U according to the (s 1 , . . . , s m ) and (t 1 , . . . , t m ) respectively. Calculation of permanents of matrices with complex entries is in the #P-hard computational complexity class [62]. In Eq. (1), ϕ(U ) is the unitary transformation acting on the Hilbert space H m,n of n photons in m modes, that corresponds to the linear evolution U on the optical modes. Due to the linearity of the evolution, ϕ(U ) is an homomorphism [63]. This means that, if a given evolution is the sequence of two linear networks W and V , the overall evolution can be written in terms of permanents of submatrices of U = V W .
In Ref. [7] it was shown that sampling (even approximately) from the output distribution of such a system is classically hard if (i) the input state |S = |s 1 , . . . s m has at most one photon per mode, (ii) U is drawn randomly from the uniform Haar measure, and (iii) the number of modes m and photons n satisfy m n 6 . Non-linear Boson Sampling. -Let us now consider the scheme of Fig. 1(b). An input state |S = |s 1 , . . . s m of n indistinguishable photons undergoes a m-mode evolution divided in three steps. While steps 1 and 3 are linear evolutions V and W drawn from the Haar ensemble, the intermediate step 2 now consists of a non-linear evolution N . This N transforms a state |R = |r 1 , . . . r m as |R N → R∈Φm,n N q1...qm r1...rm |Q , where |Q = |q 1 , . . . q m and Φ m,n is the set of tuples corresponding to n photons in m modes. In this equation, function N q1...qm r1...rm represents the transition amplitude A N (|R → |Q ) determined by the non-linear evolution. We assume N q1...qm r1...rm has an efficient classical description, e.g., it is given by the composition of a small number of few-mode non-linear transformations, or by a Hamiltonian with a simple form in terms of the field operators.
Let us now write the overall transformation of input state |S = |s 1 , . . . s m according to the three-step evolution W → N → V , which includes linear transformations W, V of the form given by Eq. (1), and the non-linear N q1...qm r1...rm : This amplitude is written as a Feynman path sum over all possible basis states just before and after the nonlinear evolution step. If the permanent distribution was peaked, it might be possible to obtain a good approximation to Eq. (2) by summing over only the dominant terms. Haar random matrices, however, display an anticoncentrated, relatively flat distribution [7]. In [64] we provide numerical evidence for this, showing that to account for (90%, 95%, 99%) of the total probability mass function, we need to calculate the probabilities associated with respective fractions ∼ (0.5, 0.6, 0.8) of all possible outcomes; moreover, this behavior is nearly independent of n and m.
Of course, there may be computational shortcuts to evaluating Eq. (2), other than the explicit sum over paths. For example, if we replace the non-linear term N by a linear term, the amplitude can be evaluated as a single permanent. This motivates us to investigate different ways to assess the complexity of non-linear Boson Sampling.
Single-mode non-linear phase shift gate. -Let us proceed by studying a specific example of non-linear evolution N consisting of a single non-linear phase gate introduced in mode x. The unitary operator describing this gate can be written asÛ nlp = exp(−ın 2 x φ). Its action on a generic m-mode state |R leads to a function N of the form N q1...qm r1...rm = exp(−ır 2 x φ) m i=1 δ ri,qi . Inserting this choice of non-linear evolution into the general expression (2) we obtain: Eq. (3) can be rearranged in the following form (see [64]): (4) Here,Ū is a unitary transformation composed by the sequence W , F and V , where F replaces the non-linear phase in layer N with a linear phase shift described by the operator exp(−ın x φ). Equation (4) clearly shows that the departure from linear evolution is due only to bunching terms, corresponding to more than a single photon in mode x.
Bounding complexity via linear-optical simulation using auxiliary photons. -An upper bound on the complexity of non-linear Boson Sampling can be obtained by devising a specific linear-optical simulation algorithm, that we discuss below for the case of a singlemode non-linear phase. The simulation is based on the results by Scheel et al. [65], describing how auxiliary photons and modes can be used, together with linear optics, to induce effective nonlinear gates. In particular, given a single mode state in the photon number basis |ψ in = i c i |i , it is possible to apply a polynomial of degree k in the photon number operator P k (n) to |ψ in by injecting the state in mode 1 of a suitably chosen (k + 1)-mode linear-optical gadget described by unitary U eff , where the auxiliary modes j = 2, . . . , k + 1 are injected with a single photon state |1 j . The desired output state |ψ out = P k (n)|ψ in is obtained upon conditional detection of a single photon on each of the auxiliary modes. If the input state has a maximum number of l photons |χ in = c 0 |0 +. . .+c l |l , a polynomial of degree l inn is sufficient to obtain the general evolution from |χ in to |χ out = c 0 |0 +. . .+c l |l with arbitrary coefficients {c 0 , . . . , c k } [66]. The success probability of the operation is equal to Pr succ = |per(U 0,1,...,1 eff )| 2 [65], where U 0,1,...,1 eff is the k×k submatrix of U eff obtained by removing row 1 and column 1 from the full matrix.
Finding the effective linear-optical simulation unitary U eff has been done previously only for a few types of gates and small k [65,[67][68][69], as the computational effort seems to scale exponentially with k. Nevertheless, even limited non-linear gate simulations can be quite versatile, as it is known that almost any non-linear gate can be combined with linear optics to generate arbitrary non-linear gates [70] -for details, see [64].
In Fig. 2 we describe the linear-optical, postselectionbased gadget that can be used to simulate single-mode non-linear gates. We see that the k-mode linear optical gadget (with k ≤ n) replaces the single-mode non-linear gate. In the gadget, mode x and the k single photons undergo the effective unitary U eff . This linear-optical simulation approximates the non-linear Boson Sampling evolution upon detection of k photons at the auxiliary output modes.
Using the state-of-the-art weak classical simulation algorithm of Clifford and Clifford [42], we can simulate the enlarged (n + k, m + k) linear optical system, postselect-ing only those events where a single photon is measured in each of the auxiliary modes (m+1, . . . , m+k) (see [64] for more details). This results in a classical simulation algorithm for the non-linear Boson Sampling experiment.
N -> single-mode non-linearity FIG. 2. Scheme for linear optics simulation of non-linear Boson Sampling with a single-mode non-linearity N on mode x between two Haar-random linear unitaries W and V . Mode x after linear transformation W is injected in the first port of a (k + 1)-mode linear optical gadget described by unitary U eff , whose other input ports are injected with k single photon states (k ≤ n). Detection of one photon in each of the k auxiliary modes of the gadget heralds a successful simulation of the single-mode non-linearity in mode x of the original interferometer (see inset).
Let us now discuss some issues that arise when using this scheme to simulate non-linear Boson Sampling in either the asymptotic regime of large numbers of photons/modes, or in the finite setting.
Non-linearities in the asymptotic setting. -Assuming uniformly drawn, Haar-random interferometer unitaries, it has been shown that the appropriate scaling between the original number of modes m and number of photons n will result in asymptotic suppression of multi- , and the post-selected linear-optical simulation using k < n auxiliary photons. In these plots, the exact distribution is calculated using the linear-optical simulation scheme with n auxiliary photons. P bunch is the probability of having more than k photons at the non-linearity site. (a) Parametric plot showing a correlation between TVD and P bunch , averaged over 100 different evolutions W and V for fixed φ = π/2. Each point corresponds to a different value of m, in the range [5, 7, . . . , 45]. (b)-(e) Correlation between TVD and P bunch for fixed number of modes m. Here, each point corresponds to different evolutions W and V for fixed φ = π/2. In all plots: blue points corresponds to k = 1, orange points to k = 2, and green points to k = 3.
photon collisions. More precisely: if m = O(n j/(j−1) ), then (j + 1)-fold collisions are suppressed, when n, m go to infinity [71]. In particular, this will be true for the photon occupation numbers at the non-linear gates. So, by choosing m = O(n k/(k+1) ), at most k photons will asymptotically be present at each non-linear gate, which means the linear-optical simulation (or classical simulation based on it) can be done with only k auxiliary photons per non-linear gate. As we will soon show, such a simulation for small k, e.g. k = 2, 3, 4 can be readily obtained. These simulations using k = 2 are sufficient for an asymptotically perfect simulation for the usual Boson Sampling regime of m = O(n 2 ). In other words, in this setting there is a precise correspondence between one single-mode non-linear phase gate and two extra auxiliary photons. More generally, the scaling of m with n dictates how many auxiliary photons are needed for asymptotically perfect simulation of a non-linear phase gate N .
Non-linearities in the finite setting. -The setting with finite n, m is experimentally relevant, and in this case there will be no strict suppression of multiphoton collisions at the non-linear gates. Setting k = n results in an exact classical simulation of non-linear Boson Sampling. When k < n we will have only an approximate simulation. As an example, when m ∼ n 2 numerical results suggest that a number of auxiliary photons equal to k = 2, 3 should provide a sufficiently accurate simulation given the large effective suppression of bunching at the output of Haar-random unitaries (see [64]).
There are two main features that increase the simulation complexity. First, finding an effective unitary U eff that uses k photons for a linear-optical simulation seems to require the computation of permanents of k × k matrices [64,65], which results in a classical runtime that increases exponentially with k. The other cost incurred is the postselection overhead. From all the simulated events on the enlarged linear-optical set-up with (n+k) photons, we only use events where the k auxiliary photons were detected at the linear-optical simulation gadget. This will happen with a probability Pr succ = |per(U 0,1,...,1 There is some evidence that the maximum value of p tends to decrease as k increases [72]. We have performed numerical simulations of non-linear Boson Sampling with a single non-linear phase in the finite setting, using the classical algorithm based on linearoptical simulation. The results are shown in Fig. 3 for (n = 3, m = 5, 9,16,27) and (n = 4, m = 16). As expected, having k = n results in exact sampling from the non-linear process, and in fact (once the appropriate gadget U eff has been determined) is numerically found to be more computationally effective than directly using Eq. (2). For fixed n, k, note that the simulation error decreases with increasing m, since bunching events become rarer. These results suggest that the crucial parameter for the simulation complexity is the scaling between n and m (see also [64]). The regime when m = O(n) is particularly interesting, as there is a trade-off between a faster classical simulation algorithm [46], and the increased complexity required to find the linear-optical gadget unitary U eff for larger k.
If the number of auxiliary photons k < n, the simulation scheme based on linear-optical gadgets will be only approximate, due to non-linear dynamics of more than k photons. The key open point in this scenario is to quantify the simulation error incurred. In Fig. 3 we provide a numerical study on how the simulation error depends on n, m, k, as quantified by the total variation distance (TVD) between the exact non-linear evolution and its simulation using k < n auxiliary photons. We observe a strong correlation between the TVD and the probability of bunching at the non-linearity. An open interesting research question is to obtain a quantitative description of this dependence between TVD and bunching, for instance by using bounds on bunching in the uniformly random, Haar ensemble of unitaries.
Discussion. -We have proposed the adoption of non-linear gates within the framework of Boson Sampling as a way to increase the computational complexity of the model. We have shown how to quantify this complexity using a linear-optical simulation with postselection, which itself can be simulated classically. For large number m of modes and n of photons, suppressed bunching allows asymptotically perfect simulation at a cost of two extra photons per non-linear phase gate introduced, if we assume m = O(n 2 ). For finite m, n and single-mode nonlinear phase gates, we identify the probability of bunching, governed by the scaling of m as a function of n, as the key factor affecting the complexity of our proposed simulation scheme.
The non-linear Boson Sampling model we propose is inherently more expressive than linear Boson Sampling. In light of the recent developments regarding first application of Boson Sampling and its variants for hybrid quantum computational models, we expect that having access to increased functionalities enabled by non-linearities can be turned into useful advantage for tasks solvable with linear Boson Sampling, as well as propose altogether new tasks solvable by noisy, intermediate-scale quantum (NISQ) devices. In parallel, an important research direction regards development of more efficient simulation schemes for non-linear gates, which is directly relevant not only to the model we propose, but to general photonic quantum computation.
Acknowledgments. -This work is supported by the ERC Advanced Grant QU-BOSS (QUantum advantage via non-linear BOSon Sampling, grant agreement no. 884676), by the European Union's Hori-zon 2020 research and innovation program through the FET project PHOQUSING ("PHOtonic Quantum Sam-plING machine" -Grant Agreement No. 899544) and by FCT -Fundação para Ciência e Tecnologia, via project CEECINST/00062/2018. We acknowledge useful comments from Scott Aaronson. Here we derive the transition amplitude for non-linear Boson Sampling, corresponding to Eq. (2) reported in the main text. For the case of linear dynamics, the transformation is represented by a unitary matrix U which describes the evolution of creation operators according to a † i → j U j,i a † j . A system of n indistinguishable photons evolving via a linear transformation are characterized by input-output transition amplitudes of the form: (S1) Here |S = |s 1 , . . . , s m is the input state configuration and |T = |t 1 , . . . , t m is the output state, where {s i } and {t i } are the respective lists of mode occupation numbers. Equation (S1) depends on the permanent of U S,T , which is the n × n submatrix of U obtained by choosing rows and columns according to |S and |T as described in Ref. [1]. The permanent is defined as: where the sum is extended over all permutations σ in S n . Equation (S1) corresponds also to writing the output state as: where Φ m,n is the set of tuples describing n photons in m modes, and γ U S,T = S|ϕ(U )|T . Due to the linearity of the evolution, ϕ(U ) is an homomorphism [2]. Consider now an evolution U that can be divided as the product of two matrices U = V W , thus corresponding to a sequence of two linear interferometers. Given the properties of ϕ(U ), the final output state can be written in two equivalent ways: which corresponds to the following identity for matrix permanents when U = V W : This expression is the starting point of the output state expansion in the presence of non-linearities.
Let us now consider the scenario where a non-linear layer N is inserted between two linear evolutions V and W . The non-linear evolution N transforms a state |R = |r 1 , . . . r m as: where |Q = |q 1 , . . . q m . The function N q1...qm r1...rm represents the transition amplitude A N (|R → |Q ) due to the non-linear evolution. It must in general satisfy appropriate constraints to represent a physical evolution.

arXiv:2110.13788v1 [quant-ph] 26 Oct 2021
We now write the overall evolution of input state |S = |s 1 , . . . s m according to the full evolution W → N → V . After the first (linear) transformation, the state evolves to: Then, after the non-linear term N , the state can be written as: Finally, after the third step corresponding to the (linear) evolution V we obtain the following expression: This equation can be rearranged as: Hence, the transition amplitude from input state |S to |T reads: corresponding to Eq. (2) of main text.

II. CUMULATIVE BOSON SAMPLING DISTRIBUTIONS
As a first step to analyze the case where a non-linear evolution is introduced in the linear-optical system, we numerically investigate whether a small fraction of matrix permanents can provide a significant contribution to the full transition amplitude. If that were true, a good approximation to the exact amplitude might be achievable by including only a few terms in the sum of Eq. (S11). To that end we perform a numerical simulation by (i) sampling N unit = 1000 linear transformations U according to the Haar measure; (ii) evaluating the full probability distribution P for each U ; (iii) sorting each probability distribution in decreasing order to obtain P S as a function of the fraction N/N tot of output configurations. For each P S we evaluate its cumulative probability C by summing over a given fraction of the combinations, obtained as C(X) = x≤X P S (x). This cumulative probability provides information on the effective fraction of permanents that contribute to the overall probability mass up to a chosen threshold p. The results of a numerical simulation for n = 3, 4, 5, 6 and m ∈ [5,37] are shown in Fig. S1.
We observe that the sorted probability distributions P S present similar trends, and that the cumulative probabilities C for fixed m are almost independent on the number of photons n. An analysis of the fraction N (C = p)/N tot of configurations that contribute to the full distribution up to a threshold p is then shown in the inset of Fig. S1. We find that fractions ∼ (0.5, 0.6, 0.8) of the output combinations are necessary to evaluate (90%, 95%, 99%) of the full probability mass, respectively, and that these values are almost independent of n and m. This feature is related to the anti-concentration conjecture for complex Gaussian matrices [1], which is relevant for the complexity of computing permanents of random matrices.

III. NON-LINEAR BOSON SAMPLING WITH A SINGLE-MODE NON-LINEAR PHASE
Here we perform the explicit calculation for the scenario in which a non-linear phase shift is introduced in mode x, corresponding to Eq. (4) in main text. The evolution operator for this transformation can be written asÛ nlp = exp(−ın 2 x φ). Evaluating its action on a generic m-mode Fock state |R reads: thus corresponding to a function N of the form: Substituting into the general expression (S11) we obtain: This expression can be written in a different form. We first observe that exp (−ır 2 x φ) = exp(−ır x φ) for r x ∈ {0, 1}. Equation (S14) can then be rearranged as: where the sums are extended over R ∈ Φ m,n for those terms with r x = 0, 1 and r x > 1 respectively. By summing and subtracting the following term: the expression for the transition amplitude can be rearranged as: Let us now define F as the unitary matrix associated to a linear transformation applying a phase shift on mode x: F = exp(−ın x φ). We observe that the transition amplitude between an input state |R and an output state |P which evolve according to F is written as: This allows us to write the following chain of equalities for the first term of Eq. (S17) whereŪ = W F V , and Eq. (S5) has been used to perform the last simplification. By replacing this into Eq. (S17) we obtain Eq. (4) of the main text: (S20)

IV. EVOLUTION WITH A NON-LINEAR SINGLE-MODE PHASE TERM
Here we present numerical simulations that quantify the perturbation, in the output distribution, induced by the introduction of a single-mode non-linear phase term N . We verify numerically the difference between two output distributions: the desired one (in the presence of N ) and one obtained by replacing the non-linear phase shift with a linear one defined by the same phase φ (i.e., the overall transformation U = W F V ). This is suggested by the form of Eq. (S20), which shows that a single-mode non-linear phase introduce a departure from a linear phase for those paths where multiple photons propagate in mode x. Finally, we verify numerically that replacing the non-linear phase shift with a linear one represents the correct benchmark for the analysis presented in this Section. We find that, in most cases, the closest linear evolution to the non-linear one W, N, V is that given by U .

A. Evolution in the finite setting
We now investigate the effect of introducing a non-linear phase term in a single mode. We begin by considering a scenario with n = 3 photons in m = 9, 16 modes according to the evolution corresponding to Eq. (S14). More specifically, we insert the non-linear phase in modes x = 5 and x = 9, respectively for m = 9 and m = 16, after a linear transformation W and before the second linear evolution V .
As a first step, we study the total variation distance TVD = 1/2 i |p i − q i | between the output distribution {p i }, corresponding to U = V W , and the output distribution {q i }, which includes the non-linear phase between W and V . The results of a numerical simulation with N unit = 1000 sets of unitary transformations {W, V } are shown in Fig. S2  (a-c). We observe that, as expected, the TVD increases for larger non-linear phases φ. For m = 9 and φ = π, the TVD reaches a value close to the one relative to a uniform distribution, thus showing that the insertion of a single non-linear phase introduces a significant change in the output distribution.
Having seen that a non-linear term introduces a significant change in the output distribution, as expected, we then analyze whether the evolution with the non-linear phase can be approximated by a linear transformation that is different from U . The form of Eq. (S20) suggests an ansatz for a linear approximation for the non-linear dynamics, namely, replacing the non-linear term N with a linear unitary F = exp(−in x φ), leading to an overall interferometer U = W F V . In Fig. S2 (d-f)  W, N, V ). For phase values close to φ = 0 and φ = π, the two scenarios present a small TVD, meaning that the non-linear evolution can be approximated by a simple linear phase. Indeed, for φ = 0 there is no phase term in both cases, whereas for φ = π the terms exp(−ınφ) and exp(−ın 2 φ) provide the same phase factor for all values of n. Conversely, the maximum difference is obtained for φ = π/2.
As a following step, we verify numerically that similar trends are obtained by directly looking at single amplitudes. In Fig. S3 we observe the same trend, as displayed by the TVD, but now for the differences in the moduli and in the argument of the amplitudes. More specifically, we consider the behavior of relative differences δ r abs = |(|A nlp W,N,V | − |A W,V |)|/|A W,V | and the absolute difference ∆ arg = | arg(A nlp W,N,V ) − arg(A W,V )|. Furthermore, numerical simulations show that the average trend for these quantities can be calculated by sampling N amp = 10 4 unitaries, and by calculating a single amplitude for each sampled unitary (see Fig. S4). This allows us to obtain an approximated estimate of δ r abs and ∆ arg without having to calculate the full distributions.
Following these results, we studied how the perturbation introduced by the non-linear phase scales with the number of photons n and the number of modes m. The results are shown in Fig. S5 for n = 3, 4, 5, and different values of m. We observe that the change induced by the non-linear phase relative to the linear transformation U = W F V decreases with the number of modes m. This trend is shown by dashed lines of Fig. S5 (g-h). This can be explained by the second term of Eq. (S20), which shows that the departure from the linear phase evolution depends on the weight of the bunching terms at the non-linearity. This is discussed in the main text, and will be also addressed in Sec. V.

B. Approximating non-linear evolution with linear transformations
In the previous section, we employed the linear transformation U as a benchmark to assess the disturbance induced by the non-linear phase shift. In particular, we reported numerical evidence that the insertion of exp(−ın 2 x φ) within linear evolutions W and V introduces a non-negligible departure from the linear behaviour. This choice of benchmark is suggested by the expression of the transition amplitude reported in Eq. (4) of the main text and Eq. (S20) above. There, we observe that the overall evolution can be written as a combination of a linear term, identified by unitary U , that includes a linear phase shift, and a second term which is present only when multiple photons propagate through the non-linearity. Now we report on additional numerical analysis to support the choice of benchmark evolution U . More specifically, we verified whether a different effective linear transformation H can have a smaller distance with respect to the full evolution W, N, V . We consider the case of n = 2, 3 photons and m = 9 modes, focusing on the φ = π/2 case. We then sample two Haar random unitary matrices W and V , and search numerically for a linear evolution H whose output distribution has a smaller TVD relative to the output considering the non-linear evolution W, N, V . This search is performed by randomly sampling from the Haar measure.
In Fig. S6 (a) we show the progress of the best TVD as a function of the number of sampled unitaries k up to N unit = 10 7 sampled Haar-random matrices. The result suggests that the TVD approximately saturates after N unit = 10 6 sampled matrices. We then repeated this numerical search for 100 different pairs W, V [see Fig. S6 (b)]. We observe that, for most tested cases, the linear transformation U is closer to the non-linear evolution W, N, V than the best Haar random matrix found via the brute force numerical search.
These numerical results then give supporting evidence that U , obtained by replacing the non-linear phase shift with a linear one, can represent a suitable benchmark for the chosen non-linearity as per Eq. (S20). Here, TVD(Haar) is the total variation distance of the non-linear evolution with respect to the best Haar random matrix found after Nunit iterations, and TVD(U ) is the total variation distance with respect to evolution U = W F V .

V. SIMULATION OF NON-LINEAR BOSON SAMPLING WITH LINEAR OPTICS AND ANCILLARY PHOTONS
In this Section we give more details on the algorithm for simulating non-linear Boson Sampling using a linear-optical gadget acting on auxiliary modes and photons.
We begin by arguing that, at least in principle, there should always exist some linear-optical gadget that simulates a given non-linear transformation. There are likely to be decompositions that are much more efficient than the one we describe here in terms of required ancillary photons, success probability of the gadget, computational complexity of constructing the gadget and so on. We leave this optimization as an interesting avenue for future research, focusing on small-scale numerical investigations. We then provide some numerical insight on the role of the number of auxiliary photons, k, by looking at the contribution of bunching terms in Haar-random photon number distributions. Finally we provide the explicit construction of the simulation algorithm, and discuss the calculation of the effective linearoptical matrices required for the simulation. Finally we provide some numerical analysis to complement the discussion reported in the main text, and leverage the linear-optical simulation into a classical simulation algorithm.

A. Universality of non-linear gadgets
We begin by quoting a result due to Oszmaniec and Zimboras [3] regarding the ability of some non-linear gates to extend linear optics to full universality within the state space of n photons in m modes.
In particular, let H m,n be the restriction of Fock space to n photons in m modes, and let LO b be set of all transformations generated by linear optics acting on H m,n . Let V be some unitary gate that is not in LO b (hence, a non-linear gate), and let LO b , V be the group of dynamics generated by sequences of elements from LO b ∪ {V }. and where |D l is the two-mode state with l particles in the first mode and n − l particles in the second. In our description of the result of [3] we have also skipped one set of conditions which applies to fringe cases and do not concern us here. This result allows us to check which conditions a particular (unitary) non-linear gate must satisfy in order to be able to approximately decompose any other non-linear gate. In other words, suppose we have some gate V that satisfies the above criteria. If our non-linear Boson Sampling setup includes some non-linearity N , we know that it is possible, in principle, to approximate N by a sequence of linear-optical elements interspersed with some copies of V .
By the result of [3], any non-linear unitary gate V suffices for the above purposes if m > 2, which means that we can always replace N by a 3-mode circuit consisting of beam splitters and applications of V . However, N is a single-mode non-linearity. Therefore it is likely that applying the above result for m = 3 is less efficient (in the sense of number of applications of V required to obtain an approximation of N to desired accuracy) than applying it for m = 2 whenever that option is available. Thus, we can test some reasonable non-linear gates to check whether they satisfy condition (ii) when m = 2. The authors of [3] proved that for a two-mode cross-Kerr gate, given by exp (−in 1n2 φ), this is true for some values of φ such as π/3. It is easy to check that the following gates also satisfy condition (ii): (i) The gate which we focus on in this paper, namelyÛ nlp = exp(−ın 2 x φ), for φ = π/3. If the number of photons n is odd, then this gate is also universal for φ = π/2.
(ii) The generalized non-linear sign shift gate NS k [4], which leaves all states with up to k − 1 photons invariant and applies |k → − |k .
So far we considered the ability of some non-linear gates to simulate others, but our main goal is to map this into a post-selected linear-optical gadget. Two issues arise if we try to directly apply the result of [3] to this question. The first is that the result is not constructive. It provides a test to decide whether some non-linear gate is universal when supplemented with linear optics, but does not output a sequence that approximates a specific desired gate. This can be solved by invoking the standard Solovay-Kitaev theorem [5] (with the caveat that, besides V , we also need to be able to apply its inverse, though that is trivial for the gates listed above).
The second issue is that, given a maximum number of photons present in the input state, n, we need there to always exist some unitary non-linear gate that can be implemented as a post-selected linear-optical gadget. To the best of our knowledge, this remains an open question. In subsequent sections we show that there exists a gadget for U nlp when φ = π/2 and for up to 4 photons, while previous work [4] has given evidence that the gate NS k can be implemented for arbitrary k with success probability that decreases as 1/k 2 .
We can combine the observations above into a procedure to obtain a post-selected gadget for any non-linearity N . We begin by using the Solovay-Kitaev algorithm to provide a sequence that approximates N as a sequence of beam-splitters and some standard non-linear gate such asÛ nlp or NS k . Then, we replace every copy of that gate by a postselected gadget. We do not pursue this approach as a practical means to perform this task, only as an in-principle argument that it should always be possible. The circuit generated by this reasoning is likely to lead to a very inefficient gadget for a typical N , and subsequently a very inefficient classical algorithm. However, we conjecture that is unnecessary, and that some gadget that uses only at most n additional photons should exist for any n-photon single-mode non-linearity. As evidence, throughout this work we show that, for many different gates, a direct gadget exists that does not require a Solovay-Kitaev type decomposition.

B. Haar-random distribution with photon-number threshold
As a first step, we need to evaluate to evaluate the contribution of bunching terms in the photon number distribution after the first Haar-random interferometer. This is important because the measurement-induced scheme of Ref. [6] allows us to simulate single-mode non-linearities acting on up to n input photons by using a linear-optical gadget and n auxiliary photons. Therefore, if the rate bunching is not too large after the first interferometer, a smaller number of auxiliary photons might suffice.
In the asymptotic limit, the Bosonic Birthday Paradox holds, stating that bunching contributions are negligible for m n 2 . To verify this limit for finite sizes, we performed numerical simulations that investigate the contribution of bunching terms at the output of Haar random matrices for n = 3, 4, 5, 6 and 7 ≤ m ≤ 37. For each system size (n, m) we numerically evaluated the overall average probability to obtain an output state with contributions having at most n max photons for each output mode (i.e. excluding configurations where at least one output mode has n out > n max photons). The results are shown in Fig. S7.
We observe that, for m = n 2 , there is a significant portion (∼ 5%) of the distribution probability where bunching terms with 3 photons on the same mode are occur, while higher order terms are almost negligible. This is a relevant aspect that must be taken into account for the simulation algorithm described below.

C. Simulating non-linearities with auxiliary photons and extended linear evolution
Simulation of non-linear Boson Sampling with single-mode non-linearities can be performed according to the approach reported in the main text. In particular, the scheme of Fig. 2 employs a set of ancillary photons and modes, and an auxiliary transformation U eff . Such transformation depends on the actual non-linear transformation to be simulated (more details on the calculation of U eff can be found in Sec. V D). The non-linear evolution is successfully simulated conditioned to the detection of state |1 on each of the output auxiliary modes.
The scheme of Fig. 2 thus defines both a linear-optical simulation algorithm, which can be implemented by building the corresponding experimental apparatus and thus running the device, and a classical simulation algorithm. In the second case, sampling can be performed via the Clifford and Clifford algorithm [7], with the additional ingredient of discarding those events which do not correspond to a correct post-selected configuration.
We can now define below the classical simulation algorithm for non-linear Boson Sampling. Algorithm 1. The following classical algorithm provides an approximate simulation of non-linear Boson Sampling with single-mode non-linearities: (i) Compute the overall m + k-mode linear transformation resulting from the sequential applications of W , U eff and V .
(ii) Sample a single event from the overall transformation with an input state having single photons in modes (1, . . . , n) and (m + 1, . . . , m + k), by using the exact algorithm due to Clifford and Clifford [7].
(iii) If one photon per mode is obtained for modes (m + 1, . . . , m + k), retain the event as valid. Otherwise, discard the event and repeat point (ii).
(iv) Iterate the procedure until the number of required samples N sample is obtained.

D. Calculation of the effective unitary U eff
A central step in the linear-optics based simulation algorithm os the calculation of the effective unitary U eff . This transformation depends on the employed non-linear evolution N , and a recipe for obtaining it has been reported in Ref. [6].
Consider an input state composed by an arbitrary superposition of up to k-photon Fock terms: We focus on the specific non-linear evolution N considered here, corresponding to a non-linear phase shift acting as exp(−ın 2 φ). The desired output state after transformation N is given by: S8. Scheme of the linear-optics gadget for simulation of single-mode non-linearities with an input of up to k photons.
The main idea is to consider the gadget shown in Fig. S8. It is composed by a (k +1)×(k +1) unitary transformation U eff . The input state |ψ in is injected into mode 1, while a single photon is injected in each auxiliary mode (from 2 to k + 1). Conditioned on the detection a single photon in each output mode from 2 to k + 1, the output state on mode 1 can be written as: Finding U eff such that all conditions (S25) are satisfied is not a trivial task. Indeed, evaluation of these conditions involves matrix permanents of size up to (2k) × (2k). Here, to determine the matrices U eff to simulate the non-linear phase shift evolution up to k = 4, we employed a numerical minimization approach. More specifically, we minimized the following quantity: Parametrization of U eff is obtained by using the Reck decomposition [8], which is employed as a mathematical tool to define the set of parameters (beam-splitter transmittivities and phase shifts) over which the minimization is performed.
Performing an unconstrained minimization, we found that the algorithm can collapse into a minimum (a value of U eff ) which provides a very small success probability Pr succ . To avoid this issue, we searched for U eff by performing a constrained numerical optimization, imposing that Pr succ ≥ p th . Here, p th is a threshold which has to be manually inserted before the minimization process. By trial and error, we found that the optimization process is sensitive to the value of p th , which thus should be carefully chosen.
By using this numerical approach, we calculated the effective unitaries U eff for the linear-optics simulation gadget, for k = 2, 3, 4 and φ = π/2. For k = 2 we found: We first observe that, in all cases, all matrices present a finite success probability. Furthermore, we observe a decreasing trend in Pr succ with respect to the photon number threshold k. We also find that, due to the presence of larger matrix permanents and to the increasing number of parameters, the computation becomes progressively more expensive for increasing k.
As a further analysis, we have then focused on the k = 2 scenario, in particular with reference to the results of Ref. [9]. In that paper, upper bounds on the success probabilities for measurement-induced non-linear gates have been derived. For instance, in the case of a non-linear phase gate of the form: the success probability is bounded by: This bound is derived by assuming unlimited resources for the measurement-induced scheme, that is, unbounded number of ancillary modes and photons. This result can be applied directly to the non-linear phase shift exp(−ın 2 φ). More specifically, the output state of up to 2 photon terms after this non-linear evolution is The state of Eq. (S32) can be mapped to Eq. (S30) by means of a linear phase shift operator as: Thus, the bound of Eq. (S31) applies to the non-linear phase shift considered here up to 2 input photons, by identifying ϕ = −2φ: As a first step, we focus on the φ = π/2 scenario. In this case, the bound of Eq. (S34) gives Pr succ ≤ 0.25. The effective unitary U (2) eff of Eq. (S27), found with the previous method, corresponds to Pr succ 0.209, and thus is close to the bound of Ref. [9].
We then investigated whether it is possible to use the same approach to obtain effective unitaries with finite success probability for all values of φ ∈ [0, π/2]. The results are shown in Fig. S9, and compared with the bound of Eq. (S34). We observe that, in fact, for all tested values of φ, Pr succ 0.1. Thus, it is possible to implement the simulation algorithm with non-negligible success probability. On the other hand, for small φ we also observe larger deviations from the bound of Eq. (S34). One reason could be the fact that this bound allows for measurement-induced schemes with unbounded resources, whereas the gadget we consider here uses a specific number of photons and modes (see Fig. S8). It is also relevant that at this stage we are only considering the maximum success probability for an exact implementation of the non-linear phase. Thus, even though there is a trivial implementation with Pr succ = 1 for the specific value φ = 0 (i.e. the identity matrix), as soon as φ deviates from 0 it might become necessary to have more modes and photons achieve success probability greater than ≈ 0.1.
However, in the regime where φ is small, an exact implementation might not be the best approach. As discussed in Sec. IV B, in that regime it might be better to incur a small error in TVD and use a linear phase as approximation to the non-linear one (which would also lead to a smaller runtime in the simulation, since it would not require simulation via extra photons). We leave the investigation of the tradeoff between these figures of merit-runtime and approximation error-as direction for future work.

E. Correlations between bunching and TVD
In the main text we have shown that the bunching probability at the non-linear site is related to the TVD between the exact distribution and the one obtained with the simulation scheme based on linear optics and auxiliary photons and modes. In Fig. S10 we report a similar analysis to the one shown in Fig. 3 of the main text. However, here we replace the bunching probability at the non-linear site P bunch with the overall bunching probability P bunch calculated for the full distribution.
In this scenario, we still observe an overall trend between the average TVD and the average bunching probability P bunch for varying number of modes [see Fig. S10 (a)]. However, when looking at the distribution for fixed number of modes m, we observe no correlations between TVD and the bunching probability (over all modes). This highlights the facy that, as shown in the main text and in Eq. (S20), the TVD is related to the bunching probability restricted to the mode with the non-linearity. FIG. S10. Analysis of total variation distance (TVD) and bunching probability, for n = 4 photons and different number of modes. The TVD is calculated between the exact probability distribution, whose single transition amplitudes are obtained from Eq. (S20), and the post-selected distribution corresponding to the linear optics measurement-based simulation algorithm for k < n. In these plots, the exact distribution is calculated as the output one from the simulation scheme with n ancillary photons and modes. At variance with Fig. 3 of the main text, the bunching probability P bunch corresponds to the overall probability obtained by summing up all terms having more than k photons in the full distribution. (a) Parametric plot showing the relative trend between the TVD and the bunching probability P bunch at the non-linear site, averaged over 100 different evolutions W and V for fixed φ = π/2. Each point corresponds to a different value of m, in the range [5, 7, . . . , 45]. (b)-(e) Correlation between TVD and P bunch for fixed number of modes m. Here, each point corresponds to different evolutions W and V for fixed φ = π/2. In all plots: blue points corresponds to k = 1, orange points to k = 2, and green points to k = 3.

F. Classical simulation algorithm implementation with finite sample size
In this section we report on numerical simulations performed to test the implementation of Algorithm 1 as a classical simulation for non-linear Boson Sampling. The results of these numerical simulations are shown in Fig. S11 for different values of n and m. In all cases, we considered a single non-linear phase positioned in the central mode between linear transformations W and V according to the scheme described in the main text. We then verified the convergence, in total variation distance, of a sample composed of N sample events with respect to the exact output distribution.
More specifically, for each tested value of n and m we have drawn samples by using different approaches. These include brute-force sampling from the exact distribution obtained with the Feynman integral approach, and sampling with the Algorithm 1 for k = n and k < n. From the plots of Fig. S11, we first observe that sampling from the exact distribution and from the linear-optics simulation algorithm with k = n provides the same trend as a function of N sample . This provides an additional numerical evidence that using Algorithm 1 with k = n leads to an exact simulation algorithm for non-linear Boson Sampling. Conversely, when k < n, the simulation saturates to an asymptotic value provided by the post-selected distribution at the output of the linear-optics gadget.