Classical simulation of linear optics subject to nonuniform losses

TL;DR

This paper demonstrates that non-uniform losses in linear optical networks enable efficient classical simulation of boson sampling, showing that certain loss patterns negate quantum advantage in these systems.

Contribution

It proves that non-uniform, path-dependent losses can be exploited to classically simulate linear optical processes, extending existing algorithms to more general input states.

Findings

Non-uniform losses can be commuted to the input layer based on network geometry.

Classical simulation algorithms can be extended to arbitrary n-photon Fock states.

Boson sampling becomes classically simulable under specific input photon distributions.

Abstract

We present a comprehensive study of the impact of non-uniform, i.e.\ path-dependent, photonic losses on the computational complexity of linear-optical processes. Our main result states that, if each beam splitter in a network induces some loss probability, non-uniform network designs cannot circumvent the efficient classical simulations based on losses. To achieve our result we obtain new intermediate results that can be of independent interest. First, we show that, for any network of lossy beam-splitters, it is possible to extract a layer of non-uniform losses that depends on the network geometry. We prove that, for every input mode of the network it is possible to commute $s_i$ layers of losses to the input, where $s_i$ is the length of the shortest path connecting the $i$th input to any output. We then extend a recent classical simulation algorithm due to P. Clifford and R. Clifford to allow for arbitrary $n$-photon input Fock states (i.e. to include collision states). Consequently, we identify two types of input states where boson sampling becomes classically simulable: (A) when $n$ input photons occupy a constant number of input modes; (B) when all but $O(\log n)$ photons are concentrated on a single input mode, while an additional $O(\log n)$ modes contain one photon each.

Figures

Equations

$$p_U\left(\vec{S}\to \vec{T}\right)=\frac{{\left|\mathrm{Per}(U_{ST})\right|}^2}{{\prod }_{i=1}^m{s}_i!{\prod }_{i=1}^m{t}_i!},$$

$$\rho \mapsto {\Lambda }_{\eta }\left(\rho \right):=\sum _{l=0}^n\left(\genfrac{}{}{0ex}{}{n}{l}\right){\eta }^l(1-\eta {)}^{n-l}{\mathrm{tr}}_{n-l}\left(\rho \right),$$

$${\Lambda }_{\mathcal{N}}={\Lambda }_{\tilde{\mathcal{N}}}\circ {\Lambda }_{\vec{\eta }}.$$

$${\Lambda }_{tot}={\Lambda }_{{\tilde{\mathcal{N}}}^{\prime }}\circ {\Lambda }_{{\vec{\eta }}^{\prime }},$$

$${\eta }_i=\underset{L_i}{\max }\prod _{j\in L_i}{\eta }_{i,j},$$

$$\mathcal{T}=nm\prod _{j=1}^m(s_j+1{)}^2{\alpha }_S,$$

$$\text{Per}(U_{ST})=\frac{1}{2^n}\sum _{v_1=0}^{s_1}\cdots \sum _{v_m=0}^{s_m}(-1{)}^{N_v}(\genfrac{}{}{0ex}{}{s_1}{v_1})\cdots (\genfrac{}{}{0ex}{}{s_m}{v_m})\prod _{i=1}^m{\left[\sum _{j=1}^m(s_j-2v_j)U_{jk}\right]}^{t_i},$$

$$\prod _{i=1}^m{\left[\sum _{j=1}^m(s_j-2v_j)U_{jk}\right]}^{t_i}.$$

$${\tau }_{S,T}=\min \left(\prod _{i=1}^m(s_i+1),\prod _{j=1}^m(t_j+1)\right){\alpha }_S{\alpha }_T.$$

$$p_U(\vec{S}\to \vec{r})=\frac{1}{n!}\frac{1}{{\prod }_{i=1}^m{s}_i!}{\left|\text{Per}{U}_{S,\vec{r}}\right|}^2,$$

$$p(\vec{r})=p(r_1)p(r_2|r_1)p(r_3|r_1,r_2)\dots p(r_n|r_1,\dots ,r_{n-1}).$$

$$\tau =n\prod _{i=1}^m(s_i+1){\alpha }_S.$$

$$m\cdot \sum _{l=1}^n{N}_l(\vec{S})\cdot \tau ,$$

$${\mathcal{T}}_A$$

$${\mathcal{T}}_B$$

$${\mathcal{T}}_A$$

$${\mathcal{T}}_B$$

$$\tau =\text{O}\left[{\left(\frac{n}{\log n}\right)}^{\log n}n\log n\right],$$

$$d_{\mathrm{TV}}\left(\{p_x\},\{q_x\}\right)=\frac{1}{2}\sum _x|p_x-q_x|.$$

$$d_{\mathrm{TV}}\left(\{p_x^{\rho }\},\{p_x^{\sigma }\}\right)\le d_{\mathrm{tr}}(\rho ,\sigma ).$$

$$|{\Psi }_0⟩=|\stackrel{n}{\overbrace{1\cdots 1}}\stackrel{m-n}{\overbrace{0\cdots 0}}⟩$$

$$\vec{\eta }=(\stackrel{k}{\overbrace{1,\dots ,1}},\stackrel{m-k}{\overbrace{\eta ,\dots ,\eta }}).$$

$$d_{\mathrm{TV}}\left(\{p_x^{\Lambda }\},\{q_x^{\Lambda }\}\right)\le \frac{{\eta }^2(n-k)}{2}+\frac{\eta (1-\eta )}{2}.$$

$$\tilde{\sigma }=\sum _{\alpha }{\tilde{P}}_{\alpha }|{\varphi }_{\alpha }⟩⟨{\varphi }_{\alpha }{|}^{\otimes n_{\alpha }},$$

$$d_{\mathrm{tr}}(\tilde{\rho },\tilde{\sigma })\le \tilde{\Delta }=\frac{{\eta }^2\tilde{n}}{2}+\frac{\eta (1-\eta )}{2}.$$

$$\sigma =\sum _{\alpha }{P}_{\alpha }|{\psi }_{\alpha }⟩⟨{\psi }_{\alpha }|,$$

$$d_{\mathrm{tr}}(\rho ,\sigma )\le \Delta =\frac{{\eta }^2(n-k)}{2}+\frac{\eta (1-\eta )}{2}.$$

$$|{\Psi }_0⟩=|\stackrel{n}{\overbrace{1\cdots 1}}\stackrel{m-n}{\overbrace{0\cdots 0}}⟩.$$

$$d_{\mathrm{TV}}\left(\{p_x^{\mathcal{N}}\},\{q_x^{\mathcal{N}}\}\right)={\Delta }_{\mathcal{N}}\le \frac{{\eta }_{\mathrm{eff}}^2(n-k)}{2}+\frac{{\eta }_{\mathrm{eff}}(1-{\eta }_{\mathrm{eff}})}{2}$$

$$\vec{\eta }=(\stackrel{k}{\overbrace{1,\dots ,1}},\stackrel{m-k}{\overbrace{{\eta }_{\mathrm{eff}},\dots ,{\eta }_{\mathrm{eff}}}}).$$

$$|{\vec{T}}\rangle\mathrel{:=}\binom{n}{t_{1},t_{2},\ldots t_{m}}^{{\color[rgb]{0,0,0}1/2}}\mathbb{P}_{\textrm{sym}}^{n,m}|{\vec{r}}\rangle,$$

Full text

To achieve our result we obtain new intermediate results that can be of independent interest. First we show that, for any network of lossy beam-splitters, it is possible to extract a layer of non-uniform losses that depends on the network geometry. We prove that, for every input mode of the network it is possible to commute $s_{i}$ layers of losses to the input, where $s_{i}$ is the length of the shortest path connecting the $i$th input to any output. We then extend a recent classical simulation algorithm due to P. Clifford and R. Clifford to allow for arbitrary $n$-photon input Fock states (i.e. to include collision states). Consequently, we identify two types of input states where boson sampling becomes classically simulable: (A) when $n$ input photons occupy a constant number of input modes; (B) when all but $O(\log n)$ photons are concentrated on a single input mode, while an additional $O(\log n)$ modes contain one photon each.

1 Introduction and relation to previous works

The recent paradigm of quantum computational advantage (or supremacy), is regarded as a promising route towards demonstrating that quantum computers are more powerful than their classical counterparts [1]. Although the systems used in this approach are not expected to be universal for quantum computation, and are often not known to perform explicitly useful tasks, they allow us to test quantum mechanics in the limit of high computational complexity. The pursuit of a near-term demonstration of quantum superiority is also aligned with the overarching goals of the field of quantum computing, as it pushes the development of technologies that will undoubtedly be necessary for scalable universal quantum computers, and allows us to better understand the effects of real-world imperfections in intermediate-size quantum systems.

One candidate for demonstrating quantum advantage is nonadaptive linear optics, or boson sampling [2]. In this model, an $n$-photon $m$-mode Fock state evolves according to a passive $m$-mode linear-optical transformation $U$ and is subsequently measured by particle-number resolving detectors [see fig. 1(a)]. This physical process leads to a probability distribution $p_{U}(\vec{T})$. The seminal contribution of [2] was to show that, subject to certain complexity-theoretic conjectures, for generic transformations $U$ it is hard for a classical computer to generate samples from a distribution suitably close to $p_{U}(\vec{T})$. Besides being an elegant and physically-motivated computational model and, arguably, the first proposal of quantum supremacy as we understand it today, boson sampling also benefits from the technology and expertise developed by the quantum optics community in the last decades [3].

Since the original proposal several small-scale boson sampling experiments have been reported [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], with state-of-the-art ranging up to five photons with near-deterministic quantum dot sources [12, 13, 14]. From the opposite end, there has also been intense activity in the development of classical simulation algorithms [16, 17, 18], and current supercomputers are expected to simulate 50-photon experiments without much difficulty [19, 17]. A more refined analysis of the complexity-theoretic arguments underpinning boson sampling has suggested 90 photons as a concrete milestone for the demonstration of quantum computational advantage [20].

Many technological and theoretical challenges remain that can undermine the scalability of boson sampling. Several sources of imperfection affect linear-optical experiments, and it is essential to understand which can be mitigated and which degrade the computational power of the model. This boundary between classical simulability and quantum advantage is an area of intense investigation, with losses in particular receiving most attention. In [21] it was shown that boson sampling retains its computational power if a constant number of photons is lost. In the other extreme, recent papers showed that lossy boson sampling can be efficiently simulated when less than $\sqrt{n}$ photons are left [22, 23] for arbitrary interferometers, or even when a suitably high constant fraction of the photons is lost for typical Haar-random interferometers [24] (see also a recent follow-up [25]). Boson sampling has also been investigated under the effect of fabrication noise in the linear-optical components [26, 27], losses combined with dark counts [28], partial photon distinguishability [29], and Gaussian noise in the experimental data [30].

In this work we build on our previous findings of [23] and solve some of the open problems posed there. The main result of [23] stated that, when less than $\sqrt{n}$ out of $n$ photons are left, it is possible to approximate a lossy boson sampling state by a state of distinguishable photons, which is known to be classically simulable. This approximation is done at the level of the input state, and holds for arbitrary linear-optical experiments. For this result to work, it was necessary to assume losses happen at the input–a reasonable assumption when the interferometer is balanced, but which might fail for unbalanced network geometries [31]. To justify this, in [23] we also showed it is possible to commute $s$ layers of mode-independent losses to the input of the network, where $s$ is the smallest number of beam splitters in any path connecting any input to any output. However, it is possible to construct a universal linear-optical transformation where $s=1$ [31], suggesting a possible route to circumvent the algorithm of [23].

Here, we extend the results of [23] in different directions, with the goal of adapting them for very unbalanced network geometries [31]. In particular, our main result is that an unbalanced network geometry cannot be used to circumvent classical simulation due to losses. Our intermediate results regarding the accumulation of losses in linear-optical networks and an extension of the algorithm for simulation of boson sampling [17] can be also be of independent interest to the quantum optics and quantum complexity theory communities.

The paper is organized as follows. First, in section 2 we give an informal overview of our results. Sections 3-6 contain formal statements and proofs of our main results. Finally, appendix A contains the technical details needed in the proof of theorem 2.

Note added. After completion of this work, the first experimental demonstration of quantum computational advantage was realised in a chip consisting of 53 superconducting qubits [32]. The experiment implemented random circuit sampling, which constitutes another paradigm for realising quantum computational superiority [33]. Around the same time, a boson sampling experiment with 20 photons occupying 60 optical modes was reported in [34].

Notation and main concepts. Throughout the paper we use the following notation. Given two positive-valued functions $f$ and $g$, we write $f=\textrm{o}(g)$ if $\lim_{x\rightarrow\infty}f(x)/g(x)=0$ and $f=\textrm{O}(g)$ if $\lim_{x\rightarrow\infty}f(x)/g(x)<\infty$. Likewise, $g=\omega(f)$ if $f=\textrm{o}(g)$. Finally, $f\approx g$ when $\lim_{x\rightarrow\infty}f(x)/g(x)=1$. By $\lfloor x\rfloor$ ($\lceil x\rceil$) we denote the largest (smallest) integer smaller (larger) or equal to $x$. We denote the total number of particles by $n$, the total number of modes by $m$, and by $[m]$ the $m$-element set $\{1,2,\ldots,m\}$. We alternatively refer to input modes as “bins”, whenever there are many photons initialized in that mode (accordingly, we refer to this situation as bunching or binning). We also refer to any state where more than one photon occupies the same mode as a collision state.

The Hilbert space of photonic states in $m$ modes is spanned by Fock states, characterized by a definite photon number $s_{i}$ in each optical mode $i$. In what follows we denote Fock states by $|{\vec{S}}\rangle$, where $\vec{S}=(s_{1},\ldots,s_{m})$ is a vector of occupation numbers. In the formalism of second quantization we write $|{\vec{S}}\rangle=\prod_{i=1}^{m}\frac{(a^{\dagger}_{i})^{s_{i}}}{\sqrt{s_{i}!}}|{0}\rangle$, where $a^{\dagger}_{i}$ is the standard creation operator associated to bosonic mode $i$ and $|{0}\rangle$ is the vacuum state. Passive linear optical transformations are associated with $m\times m$ unitary matrices $U$, which induce the transformation of modes $a^{\dagger}_{i}\mapsto a^{\prime\dagger}_{i}=\sum_{j=1}^{m}U_{ji}a^{\dagger}_{i}$. When this transformation is applied to input state $|{\vec{S}}\rangle$ the output probabilities are given by

[TABLE]

where $\mathrm{Per}(A)$ denotes the permanent of $A$ and $U_{ST}$ is an $n\times n$ matrix constructed according to the following prescription. We first construct an $m\times n$ matrix $U_{S}$ by taking $s_{i}$ copies of the $i$th column of $U$, then construct an $n\times n$ matrix $U_{S,T}$ by taking $t_{j}$ copies of the $j$th row of $U_{T}$. 111The arguments for computational hardness of Boson sampling rely on the fact that computation of the permanent is a $\#\mathrm{P}$-hard problem. Importantly, however, a device implementing boson sampling does not, in general, have the power to accurately compute or estimate permanents that appear in Eq. (1).. Boson sampling can also be framed in the language of first quantization, which we briefly use (see e.g. [35] or  the introductory section of [23] for the translation between the two descriptions).

2 Overview of main results

In this Section we give an informal summary of our main findings and discuss the relation to previous work.

Result 1: Extracting nonuniform losses from the network. We improve the procedure of [23] to commute losses to the input of a network. For simplicity, we make the following assumptions regarding losses. We assume that all losses happen within the network and are located at beam splitters, and that every beam splitter carries the same amount of loss. Although losses due to imperfect sources and detectors are experimentally important [14], they are in principle constant, whereas losses due to linear-optical elements dictate how overall transmissivity of the network scales as the experiments become more complex. Furthermore, the geometry of a network is the main source of loss non-uniformity which we wish to study. Experimentally, beam splitter losses are also relevant e.g. in integrated photonic devices, since waveguide bends required to build directional couplers can cause photons to leave in unguided modes [5]. We represent a lossy linear-optical network as in fig. 2(a).

The main figure of merit for the result of [23] is the length of input-output paths inside a network, measured in the number of beam splitters a photon has to traverse along that path. Previously it was known that, if the shortest input-output path has length $s$, we can rewrite the network in such a way that it has $s$ layers of mode-independent losses at the input. However this is not particularly useful for unbalanced networks such as that of Reck et al [31], depicted in fig. 2(a), as in that case we can only pull out one layer of uniform losses.

Our new result improves on this by showing how to extract nonuniform loss layers. More specifically, let $s_{i}$ be the shortest path between input $i$ and any output within some network. We show how to write an equivalent network that is preceded by a round of losses where input $i$ suffers the effect of $s_{i}$ consecutive loss elements. The precise formulation of this result can be found in theorem 1. Our proof is efficient and constructive, and works for networks of arbitrary geometry. The final result is illustrated in fig. 2(b) for the case of the network of Reck [31]. Overall, our result gives a more efficient description of the network where asymmetry is taken into account, and which we use to prove result 4 described later.

Result 2: Classical simulation of boson sampling with collision input states. It is well known that boson sampling is efficiently classically simulable if all $n$ photons are initialized in a single mode, i.e., input state $|{n,0,0,\ldots,0}\rangle$. This follows from the fact that the permanent of a matrix of $n$ repeated columns can be computed trivially in linear time [2], or alternatively from the fact that if all photons are initialized in the same mode they behave as distinguishable particles [23]. Here we prove a stronger version of this result, by identifying two new types of input state which allow for efficient classical simulation of boson sampling [see fig. 1(b) and fig. 1(c)]. The first (type A) corresponds to $n$ photons concentrated in a constant number of input modes. The second (type B) is when all but $O(\log n)$ photons are concentrated on a single input mode, while an additional $O(\log n)$ modes contain one photon each.

Concretely, we provide a strong simulation of boson sampling (i.e. computing outcome probabilities) by showing how an expression for the permanent function from [36] can be computed efficiently for inputs of types A and B (though a similar scaling follows from an Appendix of [37]). We then generalize a result of [17] for weak classical simulation (i.e. producing samples from the correct distribution) to allow for bunched inputs, showing that it is efficient for inputs of types A and B. Our extension of the result of [17] is also based on a physical description of the state (via first quantization) rather than combinatorial considerations, and we believe it could make the original result more transparent for researchers with a stronger physics background. The technical formulation of this result can be found in theorem 2 and in corollary 1.

Besides giving us a new understanding of the regimes where efficient classical simulability of boson sampling is possible, and ruling out attempts at demonstrating a quantum advantage with highly concentrated bosonic input states, this result also serves as a nontrivial stepping stone towards result 3 described below.

Result 3: Classical simulation of lossy boson sampling with a few lossless modes. We improve the classical simulation algorithm of [23] by allowing the lossy boson sampling instance to be coupled to some lossless modes. Concretely, suppose that $n$ photons pass through a lossy channel in such a way that less than $\sqrt{n}$ survive, and are then combined with up to $c\log n$ photons input in modes that have perfect transmission. We show that the classical simulation based on the ideas similar those of [23] works for this setting as well. The technical formulation of this result can be found in theorem 3.

The classical simulation of [23] was based on approximating a lossy bosonic state $\rho$ by another state $\sigma$ where photons effectively behave as distinguishable particles, then leveraging known results for classical simulation in that case [2]. Our new result is based on appending up to $c\log n$ single-photon states to $\sigma$, and using this state to simulate the state obtained by also appending $c\log n$ photons to $\rho$. The state $\sigma$ is a convex combination of states of type B so we can leverage result 2 above. The distance between the two states (and hence the error bound on the simulation) is the same as in [23], but the simulation is more costly due to the use of theorem 2. This result is also essential to prove result 4 that follows.

Result 4: Nonuniform losses do not avoid classical simulation. By combining results 1–3 we obtain our main result. Informally it states that the classical simulation algorithms of [23, 22] cannot be circumvented by the design of unbalanced networks such as that of Reck et al [31] (see theorem 4 for the formal statement). This was left as an open question in [23]. The idea was that, in an unbalanced network, there might be some input modes where losses are very mild, and the high chance of survival of these photons might break the assumptions behind the simulation algorithm in [23, 22]. Here we show that classical simulation remains possible as long as only $O\log(n)$ inputs have small optical path to the output. Our result can be immediately applied to lossy interferometers with the geometry proposed in [31]. This closes one avenue for scaling up boson sampling experiments, but is also relevant for current experiments since this geometry is common in integrated-photonic implementations [10, 5].

In what follows we outline our argument for the simulability of the triangular construction of [31] [cf. fig. 2(b)]. This construction has the property that the shortest path from mode $i$ to any output has length $i$ (if we label modes starting from the lowest one). Suppose the input state we wish to simulate occupies the bottom $n$ modes in an attempt to avoid losses as much as possible. To simplify the argument, suppose now that the bottom $c\log n$ modes are lossless (or include their loss channel as part of the network rather than the input state). By result 1 we can extract more than $c\log n$ loss elements for all modes above mode $c\log n+1$. However, as discussed in [23], given a per-beam splitter loss parameter $\eta$, there is always some $c$ such that, if an $n$-photon state suffers $c\log n$ rounds of losses, we expect less than $\sqrt{n}$ photons to survive on average. Therefore, the input to the network is amenable to approximation by a convex combination of states of type B, as described in result 3, and so efficient classical simulation is possible.

3 Extracting nonuniform losses from a lossy network

Linear-optical networks are the natural implementation of general multimode interferometers. They are typically built by composition of simple few-mode transformations, often easier to realize in practice, which are then arranged in some network layout (see [31, 9] for the commonly used geometries, [38, 39, 40] for mathematical perspective, and [41] for the generalization to non-unitary transformations). At the same time, each small transformation also introduces some loss. The purpose of this Section is to provide an effective description of lossy linear-optical networks that takes into account their geometric properties.

Our focus is on losses induced by linear-optical elements within the circuit, such as beam-splitters, which determine how losses scale with the size of experiments. We disregard losses in detectors and sources, as well as uniform losses associated with free-space or waveguide propagation. These do not depend on the network layout and can be easily incorporated into our description as additional losses at the input of the network (alternatively, any path-dependent transmission loss can be incorporated into the nearest beam splitter). A loss channel that acts uniformly on $m$ modes is parameterized by a transmission probability $\eta$, and its action on an $n$-photon state can be written as:

[TABLE]

where $\binom{n}{l}\eta^{l}(1-\eta)^{n-l}$ is the probability that $l$ out of $n$ particles remain. This is equivalent to adding an extra beam splitter that, with amplitude $\sqrt{1-\eta}$, diverts each photon into some inaccessible mode [see fig. 3(a)]. It can be easily shown that $\Lambda_{\eta_{2}}\circ\Lambda_{\eta_{1}}=\Lambda_{\eta_{1}\eta_{2}}$.

It is well-known that a layer of identical loss elements commutes with any linear-optical network acting only on the same modes. In fig. 3(b) this is represented for $m=2$. This is the rationale for a standard assumption that losses happen at the input of a network, as it is usually assumed that they are approximately mode-independent. In practice, however, losses can also occur in nontrivial optical elements composing a network. In what follows, we assume for simplicity that each beam splitter carries two identical loss elements at its inputs. In [23] we showed the following. Given some linear-optical network built out of the building block of fig. 3(b), let $s$ be the smallest number of beam splitters that a photon has to traverse within the network starting from any input to reach any output. It is possible to rewrite this network as a different (possibly still containing some losses) linear-optical network preceded by $s$ uniform layers of loss elements. This formalized the assumption that losses happen at the input by basing it on a geometrical property of the network, namely the smallest number of loss elements that will affect any photon as it propagates inside it.

The above result works very well for symmetric networks, such as that of Clements et al [42], depicted in fig. 4(a), but not for unbalanced network such as the triangular decomposition of Reck et al [31], depicted in fig. 4(b) (often used in integrated photonics [10, 5]). For the construction of Reck et al specifically, the previous result only allows commuting a single layer of losses to the input. This led us to conjecture it might be possible to circumvent the efficient classical simulation of [23] by using purposely unbalanced networks. Motivated by this we now prove a stronger result that allows us to treat each mode individually. Namely we show that for, each input mode $i$, it is possible to pull $s_{i}$ loss elements to the beginning of the network, where $s_{i}$ is the smallest number of beam splitters in any path between mode $i$ and any output mode.

Theorem 1**.**

Let $\mathcal{N}$ be an $m$-mode linear-optical network composed of beam splitters, each associated with a transmission probability $\eta$. Let the length of a path connecting some input and output modes be the number of beam splitters that a photon encounters along that path. Let $s_{i}$ be the shortest path between input mode $i$ and any output mode. The quantum channel associated with the entire network can be decomposed as follows:

[TABLE]

Here $\Lambda_{\vec{\eta}}$ is a nonuniform loss channel described by $\vec{\eta}=(\eta_{1},\eta_{2},\ldots,\eta_{m})$, where $\eta_{i}=\eta^{s_{i}}$ is the transmission probability for mode $i$, and $\Lambda_{\tilde{\mathcal{N}}}$ is another (possibly still lossy) linear-optical network that can be described explicitly and efficiently (in the number of modes and beam splitters in the network).

Proof.

We assume without loss of generality that the network is arranged in $D$ layers, where each layer contains at most $\left\lfloor{m/2}\right\rfloor$ beam splitters acting on disjoint pairs of modes. Under our assumption that each beam splitter is preceded by one loss element in each of its inputs, the length of some input-output path corresponds to the number of loss elements the photon encounters in that path.

We prove the Theorem by iteratively applying the identity of fig. 3(b). We begin with the last layer in the network, and build it by progressively adding each previous layer one beam splitter at a time. At each step we commute losses through the new beam splitter to make sure the new network has a specific form.

Consider initially the more general situation of fig. 5(a). We have some linear-optical network $\tilde{\mathcal{M}}$. It is composed of a network $\mathcal{M}$, that can contain any number of beam splitters and loss elements distributed in an arbitrary manner, preceded by a layer of loss elements with transmissivity $\mu_{i}$ in mode $i$. We add to the beginning of $\mathcal{M}$ a new beam splitter, $\mathcal{B}$, acting on modes $i$ and $j$ (and two corresponding loss elements with transmissivity $\eta$). We cannot necessarily commute $\mu_{i}$ and $\mu_{j}$ through $\mathcal{B}$, as they may not be equal, but we can commute losses corresponding to the smallest of them. Suppose $\mu_{i}\leq\mu_{j}$, and write $\mu_{j}=(\mu_{j}/\mu_{i})\mu_{i}$. We can commute two loss elements $\mu_{i}$ through $\mathcal{B}$, leaving behind a loss element $\mu_{j}/\mu_{i}$ in mode $j$. We are now left with the situation of fig. 5(b).

Let us now return to our network $\mathcal{N}$ and apply the above procedure. Start with the last layer of beam splitters. Clearly it can be written in the form of $\tilde{\mathcal{M}}$ of fig. 5(a). We now add each previous layer, one beam splitter at a time, always applying the argument of the previous paragraph and fig. 5. Suppose that, at step $k$ of this iteration, we have some intermediate network $\tilde{\mathcal{M}}_{k}$ as in fig. 5(a). Suppose also that $\mu_{i}=\eta^{t_{i}}$ where $t_{i}$ is the length of the shortest path between input mode $i$ of $\mathcal{M}_{k}$ and any output. In other words, we assume that at step $k$ the intermediate network satisfies the statement of our theorem.

Consider now the next beam splitter to be added to the network, $\mathcal{B}$, acting between modes $i$ and $j$. By the reasoning of fig. 5 we can now commute $\min{(t_{i},t_{j})}$ loss elements through both modes of $\mathcal{B}$. The new network consisting of $\tilde{\mathcal{M}}_{k}$ and $\mathcal{B}$ has three properties:

• (i)

the shortest path starting from either input $i$ or $j$ has length $\min{(t_{i},t_{j})}+1$,

• (ii)

there are $\min{(t_{i},t_{j})}+1$ loss elements at both inputs modes $i$ and $j$, and

• (iii)

the action of the network on modes other than $i$ and $j$ is unchanged.

Clearly the new network, $\tilde{\mathcal{M}}_{k+1}$, also satisfies the statement of our theorem. Therefore, by induction, this holds all the way until the network $\mathcal{N}$ is complete. Furthermore, this procedure is constructive by also giving explicitly the form of $\tilde{\mathcal{N}}$. We conclude the proof by noting that the above scheme is also efficient, in the sense of requiring only polynomial effort in the number of modes and beam splitters in the network.

∎

Remark 1**.**

As stated above, our result can be easily extended to include losses in the preparation and measurement stages of the experiment. Indeed, let $\Lambda_{\vec{\eta}_{in}}$ and $\Lambda_{\vec{\eta}_{out}}$ be channels describing (possibly non-uniform) losses occurring during the state-preparation and readout phases respectively. The quantum channel describing the whole physical process is then given by $\Lambda_{tot}=\Lambda_{\vec{\eta}_{out}}\circ\Lambda_{\mathcal{N}}\circ\Lambda_{\vec{\eta}_{in}}$. Combining this with (3) we get a decomposition

[TABLE]

where $\Lambda_{\tilde{\mathcal{N}}^{\prime}}=\Lambda_{\vec{\eta}_{out}}\Lambda_{\tilde{\mathcal{N}}}$ and $\vec{\eta}^{\prime}=(\eta_{in,1}\eta_{1},\eta_{in,2}\eta_{2},\ldots,\eta_{in,m}\eta_{m})$.

Remark 2**.**

Inspection of the proof of theorem 1 shows that the decomposition of Eq. (3) also holds in a network in which different optical elements induce different loss probabilities. In that case, it is still possible to write a decomposition similar to Eq. (3) where, for each input $i$, we pull out a factor

[TABLE]

where $L_{i}$ are all paths connecting input $i$ to some output, $j$ indexes all linear-optical elements connected to input $i$, each with transmissivity $\eta_{i,j}$, and the product is taken over all $j$’s in path $L_{i}$. More explicitly, we pull out sufficient losses such that the final network $\tilde{\mathcal{N}}$ has one completely lossless path from input $i$ to some output. To see how this follows from the proof of the Theorem, note two things about the procedure of fig. 5: (i) it does not require each beam splitter to have equal loss, and (ii) at every iterative step, it leaves one output port of $\mathcal{B}$ lossless. Therefore, it will create one completely lossless path from each input to some output bu construction.

Finally we point out that, if each beam splitter has a different amount of loss, this might affect conclusions about finite-sized experiments drawn from our results, but it will not change any of our asymptotic conclusions.

4 Simulation of boson sampling with binned inputs

In this Section we describe two new types of input states for which boson sampling is easy to simulate classically. One of them (type B) will be used later in a classical simulation for optical networks with nonuniform losses. We start by providing a generalization of the algorithm by P. Clifford and R. Clifford [17]. This generalization, contrary to its predecessor, works for bunched input states, and similarly avoids an exponential cost in the number of modes $m$. The proof of the following Theorem follows a similar reasoning to that of [17]. The main difference lays in the extensive usage of the first quantization formalism, which makes it easier to perform the necessary combinatorics and, in our opinion, can make the algorithm more approachable for physicists.

Theorem 2** (Generalization of the algorithm of P. Clifford and R. Clifford to collision input states).**

Consider an $n$-photon, $m$-mode boson sampling instance described by linear optical transformation $U$ and an input state $|{\vec{S}}\rangle=|{s_{1},s_{2},\ldots,s_{m}}\rangle$. We do not assume $s_{i}\leq 1$. There exists a classical algorithm that generates a sample from the output probability distribution $p_{U}\left(\vec{S}\rightarrow\vec{T}\right)$ in time upper bounded by

[TABLE]

where $\alpha_{S}$ is the number of nonzero $s_{i}$’s. Moreover, the algorithm also outputs the probability corresponding to the sampled outcome.

Proof.

Let $\vec{S}=(s_{1},s_{2},\ldots s_{m})$ encode the input state with $s_{i}$ photons in mode $i$ such that $|\vec{S}|=\sum s_{i}=n$. Define similar quantities for the output state $\vec{T}=(t_{1},t_{2}\ldots t_{m})$. In what follows the initial configuration of photons $\vec{S}$ is fixed while $\vec{T}$ denote the possible outcomes of boson sampling experiment.

We split the proof into two parts. First, we recall the results of [36] to give bounds on the computation cost on individual probability amplitudes $p_{U}\left(\vec{S}\rightarrow\vec{T}\right)$. Then we show, by generalizing the ideas from [17], how to sample from the output probability distribution by sampling the evolution of the individual particles one at the time, sidestepping the need to compute exponentially many individual probabilities.

We start by reviewing results of [36]. According to Eq. (1), the probability of some outcome $\vec{T}$ is proportional to $\lvert\textrm{Per}(U_{ST})\lvert^{2}$. This permanent can be expanded as follows:

[TABLE]

where $N_{v}=\sum_{i=1}^{m}v_{i}$. Let us compute the cost of calculating the permanent in this manner. First note there are $(s_{i}+1)$ terms for each of the $v_{i}$ sums. For each term in that multiple sum we need to compute

[TABLE]

The product requires computing $\alpha_{T}$ terms (since those with $t_{k}=0$ do not contribute), and the sum has $\alpha_{S}$ terms for similar reasons (i.e., if $s_{i}=0$ then $v_{i}=0$ as well). Since the absolute value of the permanent is invariant under exchanging $S$ and $T$, computing $\textrm{Per}(U_{ST})$ according to the expansion in eq. 6 has a runtime of

[TABLE]

For the second part of the proof, we move from the usual notation of second quantization (or mode occupation lists) where a state is represented as a list of occupation numbers [i.e. $\vec{T}=(t_{1},t_{2},\ldots t_{m})$], to a notation native to first quantization (or mode assignment lists). In this notation states are represented by multisets of $n$ elements in $[m]$, which we represent as tuples written in nondecreasing order, i.e. $\vec{z}=(z_{1}\ldots z_{n})$. Each $z_{i}$ represents in which mode photons $i$ is. Mapping back to the previous notation is easy: just identify $t_{i}$ as the multiplicity of mode $i$ in $\vec{z}$. More details on the mapping between these descriptions can be found in appendix A.

The first step in constructing the sampling algorithm is to enlarge the sample space. Instead of sampling output states as multisets $\vec{z}$, we can sample a tuple $\vec{r}\in[m]^{n}$ and rewrite it in nondecreasing order to obtain the output $\vec{z}$. It is easy to show [see Eq.(27) in appendix A] that sampling from the distribution

[TABLE]

and then reordering $\vec{r}$ is equivalent to sampling from the correct boson sampling distribution. Notice that, in the above equation, we have expanded the sample space for the outputs but not for the inputs, since those are fixed (and of type A or B). We have also defined $U_{S,\vec{r}}$ in a similar spirit to $U_{S,T}$ just above eq. 6, but now for each $r_{i}$ in $\vec{r}$ we take row $r_{i}$ of $U_{S}$. The claim that sampling from eq. 8 is equivalent to sampling from the correct distribution follows from the fact that the permanent is invariant under permutations of rows and columns, and so it is equal for all $\vec{r}$’s that, when reordered, lead to the same $\vec{z}$.

The second step is to find the marginal probability mass function (pmf) of leading subsequences of $(r_{1},r_{2},\ldots r_{n})$. The combinatorial factors that appear in the expression for this pmf make it cumbersome to write explicitly. In appendix A we show how to compute these marginal pmfs, and show that they have two properties crucial for the simulation algorithm (which they share with the corresponding pmfs for no-collision inputs, i.e. Lemma 1 in [17]):

• (i)

the probability of outcome $(r_{1},r_{2},\ldots r_{l})$, for $l\leq n$, is the sum over probabilities corresponding to all inputs consistent with choosing $l$ photons from the actual $n$-photon input state $S$, and

• (ii)

each term in that sum corresponds to an input with the same binning structure as $S$.

Therefore, if for some input $S$ the permanents $\textrm{Per}(U_{ST})$ can be efficiently computed, then the same is true for individual marginal probabilities of the distribution $p_{U}(\vec{S}\rightarrow\vec{r})$.

We now describe the sampling algorithm. This is done by sampling the marginal distributions one at a time. That is, we use the chain rule for conditional probabilities to write

[TABLE]

It is easy to see that, for fixed $r_{1}$, sampling $r_{2}$ according to $p(r_{2}|r_{1})$ is equivalent to sampling it according to (unnormalized) weights $p(r_{1},r_{2})$, and similarly for the subsequent $r_{i}$’s. Therefore, the above expression allows us to sample the elements of $(r_{1},r_{2},\ldots r_{n})$ in order, one at a time, by computing the corresponding marginal probabilities.

All that remains is to determine the running time necessary to generate a sample with the above procedure. There are $n$ steps. In step $l$, there are $m$ probabilities $p(r_{1},\ldots,r_{l})$ we need to compute according to the possible values of $r_{l}$ (recall that $r_{1}$ to $r_{l-1}$ are fixed by previous steps). Each probability requires computing some number of permanents, as per property (ii) above [cf. the discussion after eq. 29 in appendix A]. Each of these permanents can be computed in time equal to at most

[TABLE]

This follows from Eq. (7), and details can be found in appendix B. Therefore, the runtime of the algorithm can be upper bounded by

[TABLE]

where $N_{l}(\vec{S})$ is the number of $l$-photon configurations that can result by removing $n-l$ photons from the initial configuration $\vec{S}$ [see discussion preceding Eq.(30)]. Finally, note that $\sum_{l=1}^{n}N_{l}(\vec{S})=\prod_{i=1}^{m}(s_{i}+1)-1$, which concludes the proof. ∎

A straightforward inspection of Eq. (5) shows that, for the standard collision-free boson sampling input (i.e. $s_{i}=1,\ i=1,\ldots,n$), our algorithm remains with an exponential (in $n$) running time. However, using this formula we can identify two new families of input states for which boson sampling can be efficiently simulated classically.

Corollary 1** (Efficient classical simulation of boson sampling for certain classes of unbalanced inputs).**

Consider the following two types of input states:

• A

: A Fock state with photons distributed arbitrarily in $k=$O$(1)$ different input modes (or bins).

• B

: A Fock state with photons distributed in $k$=$c\log n$* bins such that $k-1$ modes contain a single photon and all remaining photons are in the *$k$th mode.

Then, from (5) it follows that there exists an efficient classical algorithm to simulate the corresponding instances of boson sampling. The runtimes of this algorithm for inputs of type A and B are, respectively,

[TABLE]

If $c$ and $k$ are constants, and in the standard boson sampling assumption that $m=$O$(n^{2})$, we have

[TABLE]

Note that $\mathcal{T}_{B}$ is not a consequence only of having $\textrm{O}(\log n)$ bins. To see this, suppose we have $\log n$ bins with $n/\log n$ photons each. In this case the running time for the computation of a single permanent according to Eq. (10) is

[TABLE]

which is super-polynomial in $n$. Therefore, the restriction that inputs of type B have most photons starting in a single mode is important.

We finally remark that types A and B of inputs states do not exhaust all the cases for which the reasoning given in the proof of Theorem 2 gives an efficient classical sampling algorithm. In particular, it was pointed out to us by Alex Moylett that the following class of input states, interpolating between types A and B, is also classically simulable:

• C

: An $n$-photon Fock state with photons distributed in $k=k_{A}+k_{B}$ input modes, where $k_{A}=\textrm{O}(1)$ and $k_{B}=c\log n$. Moreover, photons located in modes “A” are distributed arbitrarily and remaining photons follow the distribution specified in the definition of input of type B.

5 Efficient simulation of lossy boson sampling with some lossless inputs

We now combine the ideas underlying the classical simulation of lossy boson sampling from [23] with the results from the preceding section to expand the known scenarios for which boson sampling becomes classically simulable. We begin by recalling the notion of total variation distance—the figure of merit typically used to assess the accuracy of boson sampling [2]. The total variation distance between probability distributions $\{p_{x}\}$ and $\{q_{x}\}$ is given by

[TABLE]

Importantly, the total variation distance between probabilities $\{p^{\rho}_{x}\}$ and $\{p^{\sigma}_{x}\}$ generated by measuring two quantum states $\rho$ and $\sigma$, respectively, with the same measurement $\mathbf{M}=\{M_{x}\}$ is always upper bounded by the trace distance between $\rho$ and $\sigma$,

[TABLE]

We are now ready to formally state our result.

Theorem 3** (Classical simulation of lossy states when photons in some modes are not lost).**

Consider an $m$-mode $n$-photon input state

[TABLE]

and a linear-optical transformation $\Lambda=\Lambda^{\prime}\circ\Lambda_{k,\eta}$, where $\Lambda^{\prime}$ is some (perhaps lossy) linear optical transformation, and $\Lambda_{k,\eta}$ is a pure loss channel such that the first $k$ modes have zero loss, while the remaining $m-k$ modes have transmission probability $\eta$, i.e. we have

[TABLE]

Let $\{p^{\Lambda}_{x}\}$ be the probability distribution obtained by measuring the output state $\Lambda\left(|{\Psi_{0}}\rangle\langle{\Psi_{0}}|\right)$ with particle-number measurements. Then, provided $k=\textrm{O}(\log n)$, for any $\Lambda$ there exists a (weakly) classically simulable distribution $\{q^{\Lambda}_{x}\}$ which satisfies

[TABLE]

Consequently, if $\eta=o(1/\sqrt{n-k})$, the error incurred in approximating $\{p^{\Lambda}_{x}\}$ with the (classically easy) distribution $\{q^{\Lambda}_{x}\}$ goes to [math] as the number of particles $n$ tends to infinity.

Proof sketch.

The proof generalizes the main idea from [23] where the classical simulation of lossy linear optics was based on approximating the lossy bosonic state on $\tilde{n}$ particles by particle-separable states, i.e. states of the form

[TABLE]

Importantly, state $|{\phi_{\alpha}}\rangle\!\langle{\phi_{\alpha}}|^{\otimes n_{\alpha}}$ can be generated from $|{n_{\alpha}0\cdots 0}\rangle$ by application of linear-optical transformation $\Lambda_{U_{\alpha}}$.

The lossy state considered in [23] was of the form $\tilde{\rho}=\Lambda_{\eta}(|{\Phi_{0}}\rangle\!\langle{\Phi_{0}}|)$, where $\Lambda_{\eta}$ is the channel that describes uniform losses [see Eq.(2)] and $|{\Phi_{0}}\rangle\!\langle{\Phi_{0}}|$ is a standard boson sampling input state on $\tilde{n}$ particles. It was shown that there exists a particle-separable state $\tilde{\sigma}$ and a probability distribution $\{\tilde{P}_{\alpha}\}$ that can be sampled efficiently, and such that the trace distance to the target state $\tilde{\rho}$ is upper bounded by

[TABLE]

From the above and the fact that boson sampling with input states $|{\psi_{\alpha}}\rangle\!\langle{\psi_{\alpha}}|^{\otimes n_{\alpha}}$ can be simulated classically it follows that any linear-optical process $\tilde{\Lambda}$ applied to $\tilde{\rho}$ and terminating with particle-number measurements can also be (weakly) simulated classically within total variation distance $\tilde{\Delta}$.

Here we generalize the above argument to the situation at hand. The key difference is that, as the (partially) lossy input state, we take the state $\rho=\Lambda_{k,\eta}(|{\Psi_{0}}\rangle\!\langle{\Psi_{0}}|)$ for which no particles are lost in the first $k$ modes. We can approximate $\rho$ by

[TABLE]

where, in contrast to Eq. (17), states entering in the above convex decomposition are of the type $\Lambda_{U_{\alpha}}(|{\Psi_{k,\alpha}}\rangle\!\langle{\Psi_{k,\alpha}}|)$, where $|{\Psi_{k,\alpha}}\rangle=|{\overbrace{1\cdots 1}^{k},n_{\alpha},0\cdots 0}\rangle$ and $\Lambda_{U_{\alpha}}$ is a suitable linear optical unitary transformation that acts as the identity on the first $k$ optical modes. If $k=\textrm{O}(\log n)$, states of type $|{\Psi_{k,\alpha}}\rangle$ lead to classically easy boson sampling instances (they belong to type B appearing in Theorem 2). Repeating the considerations that led to Eq.(18) we conclude it is possible to choose a distribution $\{P_{\alpha}\}$ which can be sampled efficiently and for which we have

[TABLE]

By the assumptions stated above, sampling from $\Lambda^{\prime}(\sigma)$ can be done efficiently, which concludes the proof. ∎

6 Classical simulation of linear optics in unbalanced lossy networks

In this Section we combine results from Theorems 1 and 3 to show that, for arbitrary lossy linear networks, boson sampling becomes classically simulable provided there are at most O($\log n$) input modes for which the input-output paths are short.

Theorem 4** (Classical simulation in unbalanced lossy networks).**

Consider an $m$-mode lossy optical linear network $\mathcal{N}$ acting on the standard input $n$-photon state of the form

[TABLE]

Assume every beam splitter introduces losses, characterized by transmissivity $\eta$, on the modes on which it acts. Moreover, assume that at most $\textrm{O}(\log n)$ inputs are such that their shortest input-output path has length smaller than $c\log n$ for some suitable constant $c>0$. Then there exists a classical algorithm to sample from a $\{q^{\mathcal{N}}_{x}\}$ that approximates the target output distribution $\{p^{\mathcal{N}}_{x}\}$, obtained by evolving the input state through network $\mathcal{N}$, to accuracy

[TABLE]

where $\eta_{\mathrm{eff}}=\eta^{c\log(n)}=n^{-c\log(1/\eta)}$. In particular, if $c>\frac{1}{2\log(1/\eta)}$, then $\Delta_{\mathcal{N}}\rightarrow 0$ as the number of photons $n$ increase.

Proof.

The result follows directly from Theorems 1 and 3. First, note that at most $k$ input modes have their shortest input-output path of length $s_{i}\leq c\log n$. Then, by Theorem 1, we can extract from the network a layer of nonuniform losses that, up to relabeling of modes, is of the form

[TABLE]

Next, note that assumptions of Theorem 3 are met since $k=\textrm{O}(\log n)$ and the vector of losses from Eq. (23) directly corresponds to Eq. (15). By directly invoking of Theorem 3 we obtain Eq. (22), which concludes the proof. ∎

Example 1** (Triangular decomposition of linear-optical networks [31]).**

Consider the standard construction of linear-optical networks given in Fig. 2. Assume that every beam splitter in the network is lossy. The number of beam-splitters photons need to cross, as they propagate from input $i$ (counting from the bottom of the network) scales linearly with $i$. Therefore, a lossy network in this geometry satisfies the assumptions of Theorem 4 and the corresponding boson sampling instance can be simulated as described in the statement of the Theorem.

The conclusions of theorem 4 are of an asymptotic nature. We note that the error incurred in the classical simulation described here is the same as in [23], and we direct the interested reader to that reference for a discussion on how this error scales for finite-sized experiments.

7 Conclusions and outlook

We presented a comprehensive study of the effect of nonuniform losses on classical simulability of boson sampling. Our main conclusion is that linear-optical networks for which losses exhibit high degree of non-uniformity cannot evade a classical simulation strategies based on ideas from [22, 23]. Our main contributions are twofold. First, in theorem 1 we established the quantitative relation between the geometric properties of lossy networks and the amount of losses that can be pulled-out to the input of the network. Second, in theorem 4 we provided sufficient conditions for effective classical simulability of of boson sampling subject to nonuniform losses (that can result, for example, from the complicated network structure). To establish the second result we extended the classical simulation algorithm from [17] to instances of boson sampling with binned inputs. This allowed us to identify two classes of input states for which the simulation of boson sampling becomes efficient: (A) when $n$ input photons occupy a constant number of input modes; (B) when all but $O(\log n)$ photons are concentrated on a single input mode, while an additional $O(\log n)$ modes contain one photon each.

We believe that our results can help in understanding which linear-optical designs can offer scalable platforms for demonstration of quantum advantage based on boson sampling. Moreover, the classical simulation techniques presented here can be leveraged to benchmark the performance of photonic networks with complicated input states in the near future. We would like to conclude by stating a number of interesting open questions.

First, in a recent work [24] Renema et al. proved that, for typical (i.e. Haar-random) linear-optical interferometers $U$, boson sampling can be classically simulated if only a constant fraction of input photons are lost and every optical mode suffers from identical losses. The technique used in [24] was based on approximating the lossy output probability distributions by distributions originating from interference of partially distinguishable photons (with increasing degree of indistinguishability). It is an interesting question whether input states of type A or type B can be used to obtain similar results for generic multi-mode interferometers.

Another interesting question is what are the consequences of our result for Gaussian boson sampling [43], a variant of boson sampling that uses squeezed Gaussian states as input. A recent work [44] studied under which conditions lossy Gaussian boson sampling becomes easy to simulate. We conjecture that our results from section 4 can be adapted to the case of Gaussian boson sampling to prove, for example, that it is also classically simulable if the squeezed states are inputs only in $k=$O$(1)$. We also conjecture that our results of section 6 can be adapted to show that the conclusions of [44] also hold for unbalanced networks (there, the authors made the standard assumption of uniform losses).

Lastly, nonuniform losses clearly seem to undermine the quantum features of nonclassical states. It is natural to explore the consequences of our results regarding extraction of nonuniform losses from a lossy network (theorem 1) for the nonclassical features of quantum states generated by such networks [45].

We thank Michał Horodecki for interesting and fruitful discussions. M. O. acknowledges the financial support by TEAM-NET project (contract no. POIR.04.04.00-00-17C1/18-00). D. J. B. acknowledges funding from Instituto Nacional de Ciência e Tecnologia de Informação Quântica (INCT-IQ/CNPq).

Appendix A Marginal probabilities for partial photon outputs

In this Appendix, we show how to compute the marginal probabilities necessary for theorem 2. This generalizes Lemma 1 of [17] by allowing for collision inputs. The proof technique, however, is slightly different from that of [17], as we use the physical description of bosonic states via first quantization rather than combinatorial arguments. Therefore, we hope it might also make the original Lemma 1 of [17] more transparent for physicists.

Throughout this section we denote by $|{\vec{T}}\rangle$ an $n$-boson, $m$-mode Dicke state. We can write it as

[TABLE]

where $\mathbb{P}_{\textrm{sym}}^{n,m}$ is the projector onto the $n$-particle, $m$-mode symmetric subspace, and $|{\vec{r}}\rangle\mathrel{:=}|{r_{1}}\rangle\ldots|{r_{n}}\rangle$. The state $|{\vec{T}}\rangle$ is such that $t_{i}$ is the number of times state $|{i}\rangle$ appears in $|{\vec{r}}\rangle$. Conversely, $\vec{r}$ are the tuples as defined just above eq. 8.

Let now $|{\vec{S}}\rangle$ and $|{\vec{T}}\rangle$ be the (not necessarily collision-free) input and output states to a linear-optical circuit $U$. In this formalism, we can write the corresponding transition probability as

[TABLE]

where $\varphi_{U}$ is the unitary transformation induced by $U$ in the $n$-particle Hilbert space, i.e. $\varphi_{U}\mathrel{:=}U^{\otimes n}$. The above expression coincides with the standard definition of bosonic transition probabilities in terms of permanents of [2].

Expanding the output state as in eq. 24 we can write

[TABLE]

where we have used the cyclic property of the trace and the fact that $\mathbb{P}_{\textrm{sym}}^{n,m}$ acts trivially on $\varphi_{U}|{\vec{S}}\rangle\!\langle{\vec{S}}|\varphi_{U}^{\dagger}$. Equation 26 is similar to eq. 8. It represents the fact that, due to the symmetries of the input state and the dynamics (or, equivalently, the symmetries of the permanent function), in order to sample from the output distribution over symmetric states (or multisets $\vec{z}$) we can equivalently sample from the distribution over tuples $\vec{r}$ given by

[TABLE]

and reorder the outcome in nondecreasing manner.

Equation 27 gives us an expression for the full pmf over tuples $\vec{r}$. We now need to compute the marginal probabilities over leading subsequences of $\vec{r}$. To that end, let $\vec{r}_{l}=(r_{1},r_{2}\ldots r_{l})$ denote the subsequence of the first $l$ elements of $\vec{r}$, for $l\leq n$, and let $\vec{r}_{l,C}=(r_{l},r_{l+1}\ldots r_{n})$ denote the remaining elements of the list over which we wish to marginalize. We can write the marginal pmf for that subsequence as

[TABLE]

The sum inside the trace equals simply to the identity on the subspace of particles $l+1$ to $n$. We can thus write the above as

[TABLE]

where by $\mathrm{tr}_{n-l}$ we mean tracing over the last $n-l$ particles, and we have used the property that $\mathrm{tr}_{n-l}(U^{\otimes n}\rho(U^{\dagger})^{\otimes n})=U^{\otimes l}\mathrm{tr}_{n-l}(\rho)(U^{\dagger})^{\otimes l}$, as discussed in [23].

The result we need already follows from eq. 29. To see this note that $(\mathrm{tr}_{n-l}|{\vec{S}}\rangle\!\langle{\vec{S}}|)$ is a convex combination over all $l$-photon Dicke states compatible with choosing $l$ out of $n$ photons from input state $|{\vec{S}}\rangle$ i.e. states of the form

[TABLE]

where $\sum_{i=1}^{m}K_{i}=n-l$ and $K_{i}\leq S_{i}$. The probability that a particular state $|{\vec{S}-\vec{K}}\rangle$ appears in the decomposition of $(\mathrm{tr}_{n-l}|{\vec{S}}\rangle\!\langle{\vec{S}}|)$ is proportional to the multinomial coefficient $\binom{n-l}{K_{1},K_{2},\ldots K_{m}}$. Therefore, eq. 29 satisfies properties (i) and (ii) laid out above eq. 9. The fact that each probability in eq. 29 can be computed efficiently for input state of types A and B follows from the discussion in the proof of theorem 2.

Appendix B Upper bound on the cost of permanents involved in theorem 2

The maximal runtime associated to the computation of the individual permanent in the simulation algorithm given in the proof of theorem 2 can be computed via the following expression.

[TABLE]

where

[TABLE]

and the internal maximization is over all $l$-photon input configurations $\vec{S}^{\prime}$ that are compatible with with the initial state $\vec{S}$ and over $l$-photon outputs $\vec{T}$. Let us note that $\tau_{S,T}$ is a non-decreasing function of individual components of both inputs $\vec{S}$ and $\vec{T}$. Therefore, for a specific number of particles we have $\tau_{S}=\max_{T,|T|=n}\tau_{S,T}=\prod_{i=1}^{m}(s_{i}+1)\alpha_{S}n$. This concludes the proof of the upper bound from Eq.(10).