Classical simulation of bosonic linear-optical random circuits beyond linear light cone
CClassical simulation of bosonic linear-optical random circuits beyond linear light cone
Changhun Oh, ∗ Youngrong Lim, † Bill Fefferman, ‡ and Liang Jiang § Pritzker School of Molecular Engineering, University of Chicago, Chicago, Illinois 60637, USA School of Computational Sciences, Korea Institute for Advanced Study, Seoul 02455, Korea Department of Computer Science, University of Chicago, Chicago, Illinois 60637, USA (Dated: February 22, 2021)Sampling from probability distributions of quantum circuits is a fundamentally and practicallyimportant task which can be used to demonstrate quantum supremacy using noisy intermediate-scale quantum devices. In the present work, we examine classical simulability of sampling fromthe output photon-number distribution of linear-optical circuits composed of random beam splitterswith equally distributed squeezed vacuum states and single-photon states input. We provide efficientclassical algorithms to simulate linear-optical random circuits and show that the algorithms’ erroris exponentially small up to a depth less than quadratic in the distance between sources using aclassical random walk behavior of random linear-optical circuits. Notably, the average-case depthallowing an efficient classical simulation is larger than the worst-case depth limit, which is linear inthe distance. Besides, our results together with the hardness of boson sampling give a lower-boundon the depth for constituting global Haar-random unitary circuits.
I. INTRODUCTION
Quantum computers are believed to provide a compu-tational power that classical computers cannot achieve.The ultimate goal in the field of quantum computationis implementing a fault-tolerant and universal quantumcomputer to solve classically intractable and practicallyimportant problems, such as integer factorization [1] andsimulation of the real-time dynamics of large quantumsystems [2]. Since a current technology cannot build ascalable fault-tolerant quantum computer, much atten-tion has been paid to demonstrate quantum supremacyusing noisy intermediate-scale quantum (NISQ) devices[3]. While numerous interesting tasks have been theoret-ically proposed and proven to be hard to simulate classi-cally, there has been a rapid development in experimentsto manipulate an intermediate size of quantum systems[4–8].Sampling from the probability distributions of ran-domly chosen quantum circuits is one of the promis-ing problems to demonstrate quantum supremacy usingNISQ devices. Especially, various sampling problemshave been proven to be hard using classical means un-less the polynomial hierarchy collapses; boson sampling[9, 10], IQP sampling [11], Fourier sampling [12], andrandom circuit sampling [13] are the representative exam-ples. Recently, random circuit sampling was experimen-tally implemented to demonstrate quantum supremacyusing the state-of-the-art superconducting qubits for thetask that they claim a classical computer would cost alarge amount of computational time [4]. On the otherhand, bosonic circuits have also been widely studied the-oretically and experimentally since Aaronson’s seminal ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] work [9], the so-called boson sampling. Boson samplingis a sampling problem based on linear optics with single-photon input states and photo-detectors. Because of itsrelatively simple structure, there have been numerousproof-of-principle experiments of boson sampling [5–8].As the experimental system size of bosonic circuits growsrapidly, establishing adequate conditions for discriminat-ing easiness and hardness of bosonic linear-optical cir-cuits becomes more crucial for the practical demonstra-tion of quantum supremacy.Circuit depth is one of the essential parameters thatdetermine a quantum circuit’s classical simulability. Inparticular, a transition of sampling complexity from eas-iness to hardness appears as a circuit depth increases.On the one hand, when a quantum circuit is shallow, theoutput state is not entangled enough and easy to classi-cally simulate [14–16]. On the other hand, when a quan-tum circuit is deep enough to implement a global Haar-random unitary circuit that generates a large amountof entanglement, sampling from the probability distribu-tions of the output state becomes classically intractableunder plausible assumptions [9, 10, 17, 18]. Indeed, underbosonic Hamiltonian dynamics, phase transition behav-ior of sampling arising from the evolution time has beendiagnosed [19, 20]. Whereas time-evolution dynamics un-der a bosonic Hamitonian using the Lieb-Robinson bound[21] has been studied in Refs. [19, 20], we focus on linear-optical circuits composed of Haar-random beam splitterarrays and employ a mapping from random linear-opticalcircuits to classical random walk [22]. Investigating ran-dom circuits enables us to analyze a generic propertyof linear-optical circuits. Remarkably, we show that forrandom linear-optical circuits, one can efficiently sampletheir photon-number outcomes for a larger depth thanthe one found in Ref. [19]. Furthermore, currently usedlinear-optical circuits for boson sampling experimentsare not sufficiently deep in implementing a global Haar-random circuit [7, 8]. Besides, reducing the circuit depthwhen implementing a classically intractable circuit is im- a r X i v : . [ qu a n t - ph ] F e b portant to minimize the adversarial effects of photon lossin experiment [23–28]. Hence, understanding how largecircuit depth is required to attain a classically intractablequantum circuit (e.g. a global Haar-random circuit [9]) iscrucial from both practical and theoretical perspectives.In this paper, we consider bosonic circuits consist-ing of Haar-random beam splitters between adjacentmodes with squeezed vacuum states and single-photonstates input. The setup’s importance stems from thecomputational-complexity hardness result of (Gaussian)boson sampling [9, 10]. Besides, this setup is phys-ically relevant in various practical situations such ascontinuous-variable random quantum networks and es-sential to study typical behaviors of generic bosonic cir-cuits [22, 29]. We provide efficient classical approximatesamplers for such bosonic random circuits and computea bound on the circuit-depth for which the probabil-ity distribution of samplers are close to the ideal dis-tribution. Notably, the proposed classical samplers’ ap-proximation error is exponentially small in the numberof sources. Therefore, this result can provide a lowerbound of linear-optical circuits’ depth to implement clas-sically intractable outcomes. Also, the result togetherwith hardness of boson sampling shows a lower bound oflinear-optical circuits’ depth to constitute a global Haar-random unitary circuit. Finally, we discuss how our re-sults relate to current boson sampling experiments.In Sec. II, we describe the problem setup. In Sec.III A, we introduce a mapping between a bosonic randomcircuit and classical random walk. Using the mapping,we present an efficient classical algorithm to sample out-comes of the problem with squeezed vacuum state inputin Sec. III B and with single-photon state input in Sec.III C. In Sec. IV, we compare our result to the previ-ously known results. Finally in Sec. V, we discuss theimplication of our results and summarize. II. PROBLEM SETUP
Let us consider a bosonic system where N identicalsources are equally distributed in M bosonic modes ofa d -dimensional lattice and evolve under random beamsplitter arrays. The problem setting is illustrated inFig. 1. Each beam splitter is characterized by an inde-pendent Haar-random unitary matrix. More specifically,applying a beam splitter on mode ˆ a and ˆ b is described as (cid:18) ˆ a ˆ b (cid:19) → (cid:18) cos θ e iφ sin θ − e − iφ sin θ cos θ (cid:19) (cid:18) ˆ a ˆ b (cid:19) (1)where θ and φ are independently sampled from a uniformdistribution on [0 , π ). Examples of the problem are bo-son sampling [9] and Gaussian boson sampling [10], whereidentical sources correspond to a single-photon state anda squeezed vacuum state, respectively. We denote sub-lattices of lattices, each of which contains a single source,as {L α } Nα =1 satisfying ∪ Nα =1 L α = M , where M containsall M bosonic modes in the problem, and sublattices are disjoint L α ∩ L β = ∅ for α (cid:54) = β . We also denotethe set of modes having sources as S = { s , . . . , s N } ,where s j represents a source mode in the j th sublattice.Thus, |M| = M , |S| = N , and |L α | = M/N for all α .Also, each sublattice is a d -cube with an edge length of L = ( M/N ) /d . For simplicity, we assume L to be apositive integer throughout the present paper.Intuitively, when the depth of beam splitter arrays isnot sufficiently large, correlations between different sub-lattices are limited. In this case, we can treat sublat-tices to be independent each other without losing muchinformation about the system. For example, the Lieb-Robinson light cone limits the correlation between sub-lattices [21]. Using this, it has been shown that up totime t ∼ ( M/N ) /d , there exists an efficient classical al-gorithm to sample from the output distribution of gen-eral quadratic bosonic Hamiltonian dynamics [19, 20]. Inthis work, we focus on a circuit model of random 2-localpassive transformation and find a depth in which classi-cal simulation of sampling is efficient. We show that forgeneric random circuits, the depth limit in which an effi-cient simulation is possible is deeper than the one foundin Ref. [19, 20].On the other hand, when the depth of beam split-ter arrays is large enough to implement a global Haar-random unitary, any approximate classical simulation ofthe relevant output probabilities is believed to be ineffi-cient for single-mode squeezed vacuum states or single-photon states input for M ∼ N because existence of anefficient simulator leads to a collapse of polynomial hi-erarchy under some plausible assumptions [9, 10]. Thus,two different aspects in terms of depths indicate that thedepth of the random circuit plays a parameter that givesa phase transition of the relevant problem. Consideringthe problem’s relation to (Gaussian) boson sampling, let M = kN γ and assume that the measurement basis isthe photon-number basis. Let us write the transforma-tion of mode operators by a given beam splitter circuitcharacterized by a unitary operator ˆ U asˆ a j → ˆ U † ˆ a j ˆ U = M (cid:88) k =1 U jk ˆ a k . (2)We note that throughout the paper, both a Haar-randomunitary and a Haar-random unitary circuit mean a Haar-random matrix of U .We assess the accuracy of our simulation by total vari-ance distance between an ideal probability distributionand the classical algorithm’s output probability distribu-tion (cid:107)D − D a (cid:107) . Explicitly, the total variance distancebetween two probability distribution P ( n ) , P ( n ) is de-fined as (cid:107)D − D (cid:107) ≡ (cid:88) n | P ( n ) − P ( n ) | . (3)Here, the probability distribution follows Born’s rule, P ( n ) = Tr[ˆ ρ ˆΠ( n )] with a positive operator-valued mea-sure (POVM) { ˆΠ( n ) } with ˆΠ( n ) = ⊗ Mj =1 | n j (cid:105)(cid:104) n j | satis- FIG. 1. (a) Initial state in 1d architecture. Black dots are occupied by sources and empty dots are vacuum. (b) Random beamsplitter arrays in 1d. (c) Initial state in 2d architecture. (d) Random beam splitter arrays in 2d architecture. A single roundconsists of four steps (1)-(4) with 2 steps for each direction. The structure can be generalized for d -dimensional architecture,where a single round consists of 2 d steps. fying ˆΠ( n ) ≥ (cid:80) n ˆΠ( n ) = ˆ . In our case, ˆ ρ rep-resents the output quantum state after a given circuit,and | n j (cid:105)(cid:104) n j | represents a projector corresponding to n j photon outcome at the j th mode.In this paper, we define an efficient simulator as fol-lows: Definition 1.
A sampling problem is easy if there existsa classical randomized algorithm that takes as input theunitary of a given circuit and outputs a sample x froma distribution D a such that the total variation distancebetween the ideal distribution D and the distribution D a satisfies (cid:107)D − D a (cid:107) ≤ O (1/poly( N )), in runtime poly( N )(See Ref. [30]). We call such a classical algorithm as efficient sampler or efficient simulator .The main result of the present work is to find an ef-ficient simulator for random passive bosonic circuits asdescribed above. III. RESULTSA. Classical random walk behavior of randomlinear-optical circuits
We consider two different types of input sources: (i)single-mode squeezed states and (ii) single-photon states in the problem. The two sources are particularly impor-tant in that they constitute Gaussian boson sampling [9]and boson sampling [10] together with a linear-opticalcircuit, respectively. Before we investigate the two casesseparately, we present a key feature of random beamsplitter arrays, which is the average propagation speed ofinput photons. As recently studied, random beam split-ter arrays can be characterized by a classical random walk[22]. More specifically, let us consider a single source at s th mode with depth D beam splitter arrays. Since theinput states are vacuum except for the source, propertiesof such a system can be captured by the dynamics of thesource, ˆ a ( D ) k = U ( D ) k,s ˆ a (0) s + vac . (4)Thus, the crucial factors are U ( D ) k,s , where the superscript( D ) is used to emphasize depth D , and vac representscontribution from vacuum states. Using a property ofa Haar-random beam splitter, one can show that afterapplying a random beam splitter on modes k and k + 1,the average of the absolute square of U ( D ) k,s transforms as(See Appendix A) E [ | U ( D +1) k,s | ] = E [ | U ( D +1) k +1 ,s | ] = E [ | U ( D ) k,s | ] + E [ | U ( D ) k +1 ,s | ]2 , (5)which exhibits the simple random walk’s property. Here, E [ · ] represents an average over random beam splitter ar-rays. The configuration of beam splitter arrays in 1darchitecture is straightforward, which is shown in Fig. 1(a). It can be generalized to 2d architecture as shownin Fig. 1 (d) and further to d -dimensional system, therandom walk behavior of which can be described by d in-dependent random walk [22]. Now, let us consider a sin-gle source in α th sublattice and vacuum states for othermodes and define a leakage rate from the source to theoutside of the sublattice as η ≡ (cid:88) j ∈M\L α | U j,s α | . (6)Using the random walk property, the following lemmashows a property of the dynamics of the source: Lemma 1. In d -dimensional architecture, when thenumber of mode is M = kN γ and L = ( M/N ) /d , theaverage leakage rate η at depth D from a source to theoutside of the sublattice of the source is bounded by E [ η ] ≤ d exp (cid:18) − L D/d (cid:19) = 2 d exp (cid:18) − dk D N d ( γ − (cid:19) . (7)Here, E [ · ] is an average over a depth D circuit com-posed of Haar-random beam splitters. Furthermore, for0 < c < dk /
8, an arbitrary c >
0, and a depth D ≤ c N d ( γ − − c , η ≤ exp( − N c ) , (8)with a probability 1 − δ over the random beam splitters,where δ is exponentially small in N . Proof.
See Appendix B.The lemma implies that up to a depth D ≤ c L − c = c N d ( γ − − c for 0 < c < dk / c >
0, only an exponentially small portion of a source canpropagate to the outside of the sublattice. The lemmaplays a crucial role in both cases in the following sections.
B. Squeezed vacuum state input
In this section, we consider a squeezed vacuum state asan input source. Since the initial state is a Gaussian statewith a zero displacement in the phase space, its covari-ance matrix fully characterizes the input state (See Refs.[31–34] for more details about Gaussian states.). The co-variance matrix of a given quantum state ˆ ρ is definedas V jk = Tr[ˆ ρ { ˆ Q j , ˆ Q k } ] / Q ≡ (ˆ x , ˆ p , . . . , ˆ x M , ˆ p M ), satisfying the canoni-cal commutation relation [ ˆ Q j , ˆ Q k ] = i Ω jk , whereΩ ≡ M ⊗ (cid:18) − (cid:19) . (9) Note that Gaussian unitary transformations applied toa Gaussian state of a zero displacement can be equiva-lently characterized by applying symplectic transforma-tions to the Gaussian state’s covariance matrix, namely,ˆ U ˆ ρ ˆ U † ⇐⇒ SV S T , where symplectic matrices conservethe canonical commutation relation, S T Ω S = Ω.We first begin with a vacuum state, the covariancematrix of which is given by M /
2. We then apply asqueezing symplectic transformation on source modes in S , which is written as S sq = ⊕ Mi =1 (cid:18) e r i e − r i (cid:19) , (10)where r i = r for i ∈ S and r i = 0 otherwise. We as-sume a momentum-squeezing operation with a real posi-tive squeezing parameter r > V in = S sq M S Tsq = 12 ⊕ Mi =1 diag( e r i , e − r i ) . (11)In addition, beam splitters are also Gaussian unitary op-erations, so that beam splitter operations (1) betweentwo modes can be characterized by symplectic matrices S BS , which is formally written as S BS = (cid:18) cos θ e iφ sin θ − e − iφ sin θ cos θ (cid:19) ⊗ . (12)The symplectic matrix corresponding to given beamsplitter arrays of depths D can be efficiently computedby matrix multiplications of 2 M × M beam splitter sym-plectic matrices.Now, based on the intuition that for low-depth cir-cuits the correlation between sublattices is bounded, wegive a classical algorithm which efficiently simulates agiven bosonic quantum circuit. The input and outputof the algorithm are random beam splitters and a sam-ple ( n , . . . , n M ) after the measurement on the photon-number basis, respectively. The algorithm runs as fol-lows: (i) Compute the symplectic transformation of beamsplitters. (ii) For j from 1 to N , iterate (iii) and (iv). (iii)Initialize the j th source to be a squeezed vacuum stateand other modes to be vacuum. Compute the final covari-ance matrix for the j th sublattice using the symplecticmatrix obtained in (i). We keep only an L × L subma-trix of the final M × M covariance on the j th sublattice.(iv) Sample ( n j +1 , . . . , n j + L ) by computing the Hafnianrelevant to the covariance matrix and using Monte-Carlosimulation [35] (See Appendix C).The Monte-Carlo simulation in step (iv) is essentiallya classical algorithm of Gaussian boson sampling, whichgenerally requires an exponential amount of time in N .An important property which guarantees that the simu-lator is efficient is that for each j , the relevant matricesfor Hafnian computation are at most rank-2 because thesystem has a single source. Remarkably, the computa-tion cost of computing Hafnian of a K × K low rankmatrix is poly( K ) [36–38]. Thus, the above classical al-gorithm costs only polynomial time in N . As a remark,a Gaussian state is a superposition of infinitely manyphoton states, we need to set a threshold for the totalphotons number for simulation, which causes a trunca-tion error. The probability of obtaining a large numberof photons is exponentially small so that the truncationerror can be properly suppressed using an appropriatethreshold (See Appendix D for details.). The presentedalgorithm can be used to simulate Gaussian boson sam-pling with threshold detectors, whose probability distri-butions are characterized by Torontonian [8, 39]. Basi-cally, since threshold detectors can only discriminate be-tween vacuum and more than or equal to a single-photon,the probability distribution for the Monte-Carlo simula-tion is coarse-grained. Specifically, the probability of adetection event from threshold detectors can be obtainedby computing Torontonian.The output distribution of the algorithm is equivalentto that obtained by replacing the true covariance matrix V out with an approximated covariance matrix as follows: V out = v x · · · x N x v · · · x N ... ... . . . ... x N x N · · · v N → V a = v (cid:48) · · · v (cid:48) · · · · · · v (cid:48) N , (13)where v j and x jk are L × L submatrices, correspond-ing to the j th sublattice’s covariance matrix and cor-relations between j th and k th sublattices, and v (cid:48) j is a L × L covariance matrix obtained in step (iii). Con-sequently, the approximated covariance matrix V a nowrepresents a product of Gaussian states on different sub-lattices, ˆ ρ a = ⊗ Nα =1 ˆ ρ α . One can expect that such anapproximation works well when the correlation betweensublattices are small enough. Now, we find a depth wherethe approximation works.Let us recall that total variance distance can bebounded by quantum infidelity 1 − F [40],12 (cid:88) n | P ( n ) − P a ( n ) | ≤ (cid:107) ˆ ρ − ˆ ρ a (cid:107) ≤ (cid:112) − F (ˆ ρ, ˆ ρ a ) . (14)The upper bound of total variance distance indicates thatsince the outcome of a quantum state of a covariancematrix V a after a photon-number measurement can beefficiently simulated, it suffices to find how close quantumstates of the covariance matrix V out and V a are in termsof quantum fidelity with respect to a quantum circuit’sdepth. Quantum fidelity between two M -mode Gaussianstates characterized by covariance matrices V , V , one ofwhich is pure, can be written as [41, 42] F ( V , V ) = 1 (cid:112) det( V + V ) , (15) where we used a covariance matrix instead of quantumstate ˆ ρ in the argument because a covariance matrix com-pletely characterizes a quantum state in our case. Weassume that the first-moment of Gaussian states is zerothroughout the paper, which is used in Eq. (15).Let us analyze the error of the approximation (13).The following lemma shows that the quantum infidelitybetween two Gaussian states characterized by V and V can be bounded by the Frobenius norm (cid:107) · (cid:107) of their dif-ference matrix X = V − V . Lemma 2.
Let V be a covariance matrix of a Gaus-sian state in bosonic modes M obtained by applyingbeam splitter arrays on single-mode squeezed states ofsqueezing parameter r and V be a covariance matrix ofa Gaussian state. If X = V − V is small, (cid:107) X (cid:107) (cid:28) − F ( V , V ) ≤ (cid:107) X (cid:107)√ N cosh 4 r + O ( (cid:107) X (cid:107) ). Proof.
See Appendix E.The lemma guarantees that if (cid:107) X (cid:107) is small enough, ourapproximation is accurate in terms of quantum fidelity.Intuitively, one can expect that such an error is smallwhen a circuit is so shallow that a source from each sub-lattice has not propagated much to other sublattices toconstitute correlations between sublattices. The follow-ing lemma shows that the correlations are bounded bythe leakage rate from a source in a sublattice to outsideof the sublattice. Lemma 3.
Let V be a covariance matrix of a Gaussianstate in bosonic mode M obtained by applying beamsplitter arrays on single-mode squeezed states of squeez-ing parameter r and V a is a covariance matrix obtainedby the proposed classical algorithm. Then, the Frobeniusnorm of X = V − V a is bounded as (cid:107) X (cid:107) ≤ e r N ( η + 2 √ η ) , (16)where η is the leakage rate from a source s in L α to othersublattices M \ L α . Proof.
See Appendix F.Therefore, the interaction between a sublattice and theoutside of the sublattice is characterized by the leakagerate of the source in the sublattice to the outside. Sincewe already have the bound of the leakage rate with re-spect to the number of sources N , modes M and depth D from lemma 1, we are ready to prove the followingtheorem. Theorem 1. (Easiness for squeezed-state input) Con-sider an M -mode bosonic system of d -dimensional archi-tecture with N number of equally distributed squeezedvacuum sources and M = kN γ with γ > D ≤ c N γ − /d − c with aprobability 1 − δ over the random beam splitters. Here, δ is exponentially small in N and for 0 < c < dk / c > Proof.
Since the proposed algorithm is efficient, it sufficesto find a depth in which the algorithm’s error is small in N . Using lemma 2 and lemma 3, we obtain the upperbound of total variance distance for an arbitrary POVMas (cid:107)D − D a (cid:107) ≤ (cid:18) N r (cid:107) X (cid:107) + O ( (cid:107) X (cid:107) ) (cid:19) / ≤ cη / N / + O ( η / ) , (17)where η is a leakage rate, and c > N , (cid:107)D − D a (cid:107) ≤ c exp (cid:18)
54 log N − N c (cid:19) , (18)for D ≤ c N d ( γ − − c with 0 < c < dk / c > − δ over the random beam splitters.Thus, the total variance distance is exponentially smallin N with probability 1 − δ .The theorem shows that our classical algorithm cansimulate a quantum circuit with an exponentially smallerror for depth D ≤ c L − c = c N d ( γ − − c in polyno-mial time in N . C. Single-photon states
In this section, we consider single-photon sources as aninitial state with the same configuration. Single-photonsources are also important because the classical simula-tion of the output distribution after a Haar-random uni-tary circuit is proven to be hard under some plausibleassumptions [9]. After a passive transformation (2), thestate evolves as | ψ (cid:105) = N (cid:89) α =1 ˆ a † s j | (cid:105) → | ψ (cid:105) = N (cid:89) α =1 (cid:32) M (cid:88) k =1 U k,s α ˆ a † k (cid:33) | (cid:105) . (19)We follow the method presented in Ref. [19] with a differ-ent correlation bound. Particularly, while Ref. [19] em-ploys the Lieb-Robinson bound [21], we use a correlationbound implied by Lemma 1. Let an M -dimensional tu-ple r represent the initial configuration of single-photonsand an M -dimensional tuple t the final configuration ofthe photons from photon-number detection. For exam-ple, r = (0 , , , ,
0) presents that we begin with singlephoton sources at the second and fourth modes when M = 5, and t = (0 , , , ,
0) presents that we detectphotons at the second and third modes. Also, we definelists in and out to represent positions of input and out-put photons, which gives in = (2 ,
4) and out = (2 ,
3) forthe above example. In our case, the input list is givenby in = ( s , . . . , s N ). The probability amplitude for theconfiguration is then written as φ ( r, t ) = 1 √ t ! (cid:88) σ N (cid:89) j =1 U out σ ( j ) , in j , (20) where t ! ≡ (cid:81) Mi =1 t i !, the summation is taken over allpossible permutations σ . The probability distributionis then given by P ( r, t ) = | φ ( r, t ) | = 1 t ! (cid:88) σ N (cid:89) j =1 | U out σ ( j ) , in j | +1 t ! (cid:88) σ (cid:54) = τ N (cid:89) j =1 U out σ ( j ) , in j (cid:32) N (cid:89) k =1 U out τ ( k ) , in k (cid:33) ∗ . (21)Our strategy is to keep the first line and ignore the secondline, which is the same algorithm used in [19]. Such anapproximation can be interpreted as treating bosons asdistinguishable particles.The algorithm runs as follows: (i) Iterate (ii) and (iii)for j from 1 to N , (ii) Compute the probability distribu-tion | U k,s j | over 1 ≤ k ≤ M , (iii) Sample k according tothe distribution and save it. After iterating (ii) and (iii)for j from 1 to N , we get N -photon outcome. It is appar-ent that step (ii) and (iii) can be efficiently performed forgiven beam splitter arrays. Thus, the algorithm’s compu-tational time cost is polynomial in N . Hence, it sufficesto show that the algorithm’s accuracy is polynomiallysmall in N .Back to error analysis from Eq. (21), by rearrang-ing permutations and using | a + b | ≤ | a | + | b | , we findan upper-bound of the total variance distance (See Ap-pendix G for details.) (cid:15) = (cid:88) t t ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) σ (cid:54) = τ N (cid:89) k =1 U out σ ( k ) , in k (cid:32) N (cid:89) l =1 U out τ ( l ) , in l (cid:33) ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (22) ≤ (cid:88) j (cid:88) ρ (cid:54) =Id (cid:89) i | U j i , in i || U j i , in ρ ( i ) | , (23)where the sum j is taken over ordered tuples ( j , . . . , j N )with 1 ≤ j l ≤ M for all 1 ≤ l ≤ N , and the sum ρ is takenover all possible permutations except for the identity Id.Equivalently, we can express the upper-bound as12 (cid:88) j (cid:88) ρ (cid:54) =Id (cid:89) i | U j i , in i || U j i , in ρ ( i ) | = 12 (cid:88) I C (cid:89) i ∈I C C i (cid:89) k ∈I D D k , (24)where the sum I C is over nonempty subsets I C of theindices of in and I D is the complement of I C , and C i ≡ (cid:88) j i ,ρ : ρ ( i ) (cid:54) = i | U j i , in i || U j i , in ρ ( i ) | , (25) D i ≡ (cid:88) j i ,ρ : ρ ( i )= i | U j i , in i || U j i , in ρ ( i ) | = 1 . (26)Equivalently, I D represents the set of fixed points. Now,we show that C i can be bounded by the leakage rate η : Lemma 4. C i ≤ (cid:112) ηkN γ +1 . (27) Proof.
Appendix G.Finally, we show that the approximation algorithm isan efficient classical simulator for D ≤ c N γ − /d − c for 0 < c < dk / c > Theorem 2. (Easiness for single-photon-state input)Consider an M -mode bosonic system of d -dimensionalarchitecture with N number of equally distributed single-photon sources and M = kN γ with γ > D ≤ c N γ − /d − c with aprobability 1 − δ over the random beam splitters, withan exponentially small δ in N and for 0 < c < dk / c > Proof.
Let c be the upper-bound of C i presented inlemma 4. From Eq. (24), (cid:15) ≤ N (cid:88) l = |I C |≥ (cid:18) Nl (cid:19) c l = 12 [( c + 1) N − N c −
1] (28) ≤ exp[2 log N + ( N −
2) log(1 + c ) + log c ] (29) ≤ exp[2 log N + 2 (cid:112) ηkN γ +3 + log(4 ηkN γ +1 )] . (30)Here, we have used Taylor’s theorem that[( c + 1) N − N c − / (cid:0) N (cid:1) (1 + h ) N − c for some h ∈ [0 , c ] for the second inequality. Using lemma 1, (cid:15) ≤ exp [ − N c + o ( N c )] , (31)for depth D ≤ c N γ − /d − c with a probability 1 − δ over random beam splitter circuits, with an exponentiallysmall δ in N and for 0 < c < dk / c >
0. Hence, the total variance distance (cid:15) isbounded by an exponentially small number in N . IV. RELATION TO PREVIOUS RESULTS
In this section, we compare our result with the re-lated previous results in the literature. First, we com-pare our result with the one obtained in Refs. [19, 20]where time evolution of single photons under spatiallylocal quadratic bosonic Hamiltonians has been stud-ied. While both study a passive nearest-neighbor-interacting quadratic Hamiltonian, a crucial differenceis that we have considered random beam splitter cir-cuits rather than time evolution under a general spa-tially local quadratic bosonic Hamiltonian. Specifically,our result shows that generic linear-optical circuits up todepth D = O ( L − (cid:15) ) = O ( N d ( γ − − (cid:15) ) are easy to clas-sically simulate except for an exponentially small frac-tion of random circuits, whereas the sampling complex-ity phase-transition of the Hamiltonian dynamics wasfound at time t = Θ( L ) = Θ( N γ − d ) in Refs. [19, 20].Combining the two results, there exists an exponentiallysmall portion of linear-optical circuits between depths D = O ( N γ − d + (cid:15) ) and D = O ( N d ( γ − − (cid:15) ) which is hard to classically simulate, where (cid:15) , > D = O (log M ).In addition, our result renders an important implica-tion on boson sampling experiments. We first note that ifone uses structured beam splitter arrays [45–47], depth D = M is sufficient to implement a Haar-random uni-tary matrix U by appropriately manipulating each beamsplitter’s transmissivity. Nevertheless, there is a prac-tical reason why we need to focus on generic circuits inexperiments. As boson sampling experiments’ circuit sizegrows to demonstrate quantum supremacy, manipulatingeach beam splitter becomes very demanding. Indeed, theunitary matrix of beam splitter circuits in current bosonsampling experiments presents a slight deviation from aHaar-random unitary matrix [6, 7]. In the case of bosonsampling experiments using 1 d architecture [5, 6], oureasiness bound indicates that the depth of circuits hasto increase as D > c N γ − − c ∼ M γ − γ − c (cid:48) unlessbeam splitters are highly well controlled. Especially for γ = 2, the depth less than a linear-scaling of M is easyto simulate. In fact, the most recent boson sampling ex-periment has employed 2 d architecture to implement aglobal Haar-random unitary [7]. If one scales up the sizeof the experiment using the same structure, one can eas-ily check that a depth increases as D ∝ √ M . Althoughtheir configuration of beam splitters and input states isdifferent from ours, it is still interesting to see that oureasiness bound, D ∝ M γ − γ − c (cid:48) , is smaller than the ex-periment for γ ≤
2. Thus, such a configuration mightbe proper to increase the system size prohibiting an ef-ficient simulation based on our results for γ ≤
2, whileone needs to invent a better configuration if the numberof modes is much larger than the number of sources suchthat M ∝ N γ with γ > V. DISCUSSIONS AND CONCLUSIONS
In the present work, we have analyzed bosonic linear-optical random circuits with squeezed vacuum states andsingle-photon states and shown an efficient classical sim-ulator for depth D = O ( N d ( γ − − (cid:15) ) with an arbitrar-ily small (cid:15) >
0. Such an easiness-depth limit showsthat if γ ≥
2, for an one(two)-dimensional configuration,bosonic circuits of a depth less than quadratic (linear)can be classically simulated. The condition γ ≥ U in Eq. (2),the depth of a beam splitter M -mode circuit needs to sat-isfy D > c N d ( γ − − c ∼ M d γ − γ − c (cid:48) , where γ ≥ c (cid:48) = c /γ > γ → ∞ , a required depth for a Haar-random unitarymatrix approaches to D = O ( M d − (cid:15) ) where (cid:15) is an arbi-trarily small number. It suggests that, for example, forone(two)-dimensional architectures, one needs to consti-tute random beam splitters of a depth at least quadratic(linear) in the number of modes.Our Lemma 1 has an implication on how fast quan-tum information propagates, which is an interesting andimportant topic from fundamental and practical perspec-tives [48–51]. While a Lieb-Robinson light cone [21] givesa general bound of information propagation speed, a so-called Frobenius light cone [52, 53] was recently shown tobe tighter than the Lieb-Robinson light cone for quantumstate-transfer tasks [54] as well as many-body quantumchaos, revealing a gap between the two light cones [53].Basically, a Frobenius light cone is obtained by averag-ing a correlation function over input states while a Lieb-Robinson light cone takes into account the worst-caseinput state. As the gap of the worst-case and averagecase was found in the literature, our average-case studyfor random circuits together with the previous worst-casestudy for Hamiltonian dynamics in Ref. [19] also clearlyreveals a gap of information propagation speeds from adifferent aspect.It is worth emphasizing that diffusive dynamics lead-ing to the sub-linearity propagation of quantum infor-mation in our system is fundamentally related to theU(1) symmetry, which conserves total excitation num-ber [22, 55, 56], whereas quantum information of typicalsystems without symmetry spreads ballistically in space-time [50, 57, 58]. Hence, our observation shows that extrasymmetry will affect the quantum information propaga-tion, and on the other hand also enable more efficientsimulation algorithms. It would be an interesting futurework to investigate different kinds of symmetry that af-fect the system’s sampling complexity.Throughout the work, we have assumed an equally dis-tributed sources and this assumption plays an importantrole in the proposed simulators. Nonetheless, once aHaar-random unitary circuit is implemented, the hard-ness of (Gaussian) boson sampling is independent of aninitial configuration, i.e., sources need not be equally dis-tributed. Thus, the implication of our result on the hard-ness of random bosonic circuits does not rely on an initialconfiguration. Still, it is an interesting open question tofind a depth for easiness assuming different configura-tions. Also, we have focused on squeezed vacuum states and single-photon states because of the relation to (Gaus-sian) boson sampling. We expect that a similar result canbe derived for a different type of states because the keyproperty of the problem was a mapping to a simple ran-dom walk presented by Lemma 1, although we leave sucha generalization as an open question.In conclusion, we have studied classical simulability ofbosonic random circuits consisting of beam splitter ar-rays. We have investigated two kinds of bosonic sources:squeezed vacuum states and single-photon states. Wehave provided efficient classical samplers that approx-imately simulate generic bosonic quantum circuits ofbeam splitters and shown a depth limit that the algo-rithms work. ACKNOWLEDGMENTS
We thank Owen Howell, Alireza Seif, RoozbehBassirian, Abhinav Deshpande for interesting and fruitfuldiscussions. C.O. and L.J. acknowledge support from theARL-CDQI (W911NF-15-2-0067), ARO (W911NF-18-1-0020, W911NF-18-1-0212), ARO MURI (W911NF-16-1-0349), AFOSR MURI (FA9550-15-1-0015, FA9550-19-1-0399), DOE (DE-SC0019406), NSF (EFMA-1640959,OMA-1936118), and the Packard Foundation (2013-39273). Y. L. acknowledges National Research Founda-tion of Korea a grant funded by the Ministry of Scienceand ICT (NRF-2020M3E4A1077861) and KIAS Individ-ual Grant (CG073301) at Korea Institute for AdvancedStudy. B.F. acknowledges support from AFOSR (YIPnumber FA9550-18-1-0148 and FA9550-21-1-0008). Thismaterial is based upon work partially supported by theNational Science Foundation under Grant CCF-2044923(CAREER).
Appendix A: Random walk behavior of randombeam splitter arrays
Let us consider a single source and random beam split-ter arrays acting on the source. Let us consider twomodes ˆ a ( D ) k and ˆ a ( D ) k +1 at depth D , which are written asˆ a ( D ) k = U ( D ) k,s ˆ a (0) s + vac , and ˆ a ( D ) k +1 = U ( D ) k +1 ,s ˆ a (0) s + vac , (A1)where vac represents the contributions from initial modesoccupied by vacuum states. After applying a beam split-ter between them, the modes are transformed asˆ a ( D +1) k = e iφ ( U ( D ) k,s cos θ + e iφ U ( D ) k +1 ,s sin θ )ˆ a s, + vac= U ( D +1) k,s ˆ a s, + vac , (A2)ˆ a ( D +1) k +1 = e iφ ( U ( D ) k +1 ,s cos θ − e − iφ U ( D ) k,s sin θ )ˆ a s, + vac= U ( D +1) k +1 ,s ˆ a s, + vac , (A3)where cos θ and sin θ represent the beam splitter’s trans-missivity and reflectivity, respectively, and U ( D +1) k,s = e iφ ( U ( D ) k,s cos θ + e iφ U ( D ) k +1 ,s sin θ ) , and U ( D +1) k +1 ,s = e iφ ( U ( D ) k +1 ,s cos θ − e − iφ U ( D ) k,s sin θ ) . (A4)After averaging the transmissivity cos θ over a uniformdistribution of θ ∈ [0 , π ), and the phases φ , φ and, φ over a uniform distribution of [0 , π ), we obtain E [ | U ( D +1) k,s | ] = E [ | U ( D +1) k +1 ,s | ] = E [ | U ( D ) k,s | ] + E [ | U ( D ) k +1 ,s | ]2 , (A5)which shows that the transmissivity of a random beamsplitter array follows a random walk behavior. Appendix B: Proof of Lemma 1
Proof.
We first note that on average the random circuitscan be characterized by a symmetric random walk andthat the goal is to find an upper-bound of the leakagerate. We observe that the leakage rate assuming an in-finite number of modes is always larger than one withboundaries because boundaries makes the walker returnto the initial lattice. Thus, it is sufficient to find anupper-bound assuming an infinite number of modes.For one-dimensional random walk, the probability ofpropagating farther than l in step t is given by P ≤ (cid:18) − l t (cid:19) . (B1)Since we are interested in the leakage rate out of a sub-lattice, we set l = L/ L = ( M/N ) /d for a d -dimensional case. Taking into account the dimension ofthe circuit, we set t = D/d . Thus, the leakage rate canbe upper-bounded as E [ η ] ≤ − (cid:20) − (cid:18) − ( L/ D/d (cid:19)(cid:21) d ≤ d exp (cid:18) − L D/d (cid:19) = 2 d exp (cid:18) − dk D N d ( γ − (cid:19) . (B2)Using Markov’s inequality, we obtain P ( η ≥ a ) ≤ E [ η ] a = 2 da exp (cid:18) − dk D N d ( γ − (cid:19) . (B3)Especially for depth D ≤ c N d ( γ − − c and a =exp( − N c ) with an arbitrary c >
0, we find P ( η ≥ exp( − N c )) ≤ d exp (cid:18) − dk c N c + N c (cid:19) , (B4)which converges to zero for large N when c < dk / Appendix C: Gaussian boson sampling classicalalgorithm
In this Appendix, we present more details about Gaus-sian boson sampling and a classical algorithm for Gaus-sian boson sampling. We consider N number of sourcesin M modes with random beam splitter arrays. Let usconsider a covariance matrix Σ ij = Tr[ˆ ρ { ˆ ξ i , ˆ ξ j } ] / ρ with ˆ ξ = (ˆ a , . . . , ˆ a M , ˆ a † , . . . , ˆ a † M ) to fol-low a notational convention used in Ref. [10, 35]. Notethat a covariance Σ can be easily obtained by a covari-ance matrix V , which we have used in the main text.In the case of general Gaussian input states of a covari-ance matrix Σ and a zero displacement, the probabilityof each outcome ( n , n , . . . , n M ) obtained by the mea-surement in photon-number basis is given by [10] P ( n , n , . . . , n M ) = 1 (cid:112) det(Σ + M /
2) Haf( A n ) n ! n ! · · · n M ! , (C1)where A = Y M [ M − (Σ + M / − ] , Y m = (cid:18) M M (cid:19) . (C2)Here, A n is a matrix obtained by repeating the j th and( j + M )th row and column of A for n j times for 1 ≤ j ≤ M .Now, let Σ ( k ) be the reduced covariance matrix ofthe first k modes. From the reduced covariance matrix,one can constitute matrices O ( k ) ≡ k − ( Q ( k ) ) − , and A ( k ) ≡ Y ( k ) O ( k ) , and the latter gives a marginal proba-bility distribution on the first k modes as P ( n , . . . , n k ) = 1 (cid:112) det(Σ ( k ) + k /
2) Haf( A ( k ) n ) n ! · · · n k ! , (C3)where A ( k ) n is the matrix obtained by repeating rows andcolumns i and i + M of the matrix A n i times for 1 ≤ i ≤ k .By directly computing Hafnian of A ( k ) n , one may usea Monte-Carlo method to sample outcomes as follows[35]: First, we compute a marginal probability distribu-tion P ( n ) for 0 ≤ n ≤ n max and sample n ∗ based onthe distribution. Since n j can be infinitely large in prin-ciple, we choose an upper-threshold of n j carefully (SeeAppendix D). After sampling n ∗ , we compute P ( n ∗ , n )for 0 ≤ n ≤ n max , sample n ∗ according to the conditionprobability distribution P ( n | n ∗ ) = P ( n ∗ , n ) P ( n ∗ ) . (C4)For k th step of the above procedure, we have( n ∗ , . . . , n ∗ k − ), so that we compute P ( n ∗ , . . . , n ∗ k − , n k )0and sample n k according to the conditional probabilitydistribution P ( n k | n ∗ , . . . , n ∗ k − ) = P ( n ∗ , . . . , n ∗ k − , n k ) P ( n ∗ , . . . , n ∗ k − ) . (C5)We iterate this procedure until we get ( n , . . . , n M ).Since the sampling procedure can be efficiently executed,the dominant computational cost is to compute the Haf-nian of the relevant matrices. when we compute theprobability distribution P ( n ), one needs to compute theHafnian of a matrix. In general, the computational costto compute the Hafnian of a complex 2 k × k matrix is O ( k k )[38], and thus the computational cost to simulateGaussian boson sampling is O ( M n n max ). However,if the rank R of a matrix is low, the computational costto compute its Hafnian can be significantly reduced as (cid:0) M + R − R − (cid:1) poly(2 M ) [37, 38]. Using this, we now showthat classical simulation of beam splitter circuits begin-ning with a single squeezed vacuum state can be effi-cient. Note that if we assume threshold detectors insteadof photon-counting detectors [39], the above probabilitydistributions are coarse-grained such that the outcomesare vacuum or more than or equal to a single-photon.The probability for the latter is obtained by computingTorontonian.Especially when the input state is comprised ofsqueezed vacuum states of real squeezing parameters r j ,we can simplify the matrix as, where A = B ⊕ B ∗ , where B = K ( ⊕ Mj =1 tanh r j ) K T . (C6)Here, K is an M × M matrix describing beam splitterarrays, and r j is a squeezing parameter, which is non-zeroonly for the sources. In this case, the covariance matrixΣ of the output state is obtained byΣ = 12 (cid:18) K K ∗ (cid:19) SS † (cid:18) K † K T (cid:19) , (C7)where S = (cid:18) ⊕ Mj =1 cosh r j ⊕ Mj =1 sinh r j ⊕ Mj =1 sinh r j ⊕ Mj =1 cosh r j (cid:19) (C8)represents the squeezing symplectic matrix. Here, weassume a single squeezed state input, so that r > r j = 0 for j > A n is a rank-2matrix for any ( n , . . . , n M ). Now, it suffices to show that A ( k ) n is also rank-2. It is known that a reduced covariancematrix from a single-source covariance matrix can alwaysbe written as [59]Σ ( k ) = 12 (cid:18) L L ∗ (cid:19) S ( k ) D ( k ) S ( k ) † (cid:18) L † L T (cid:19) , (C9) D ( k ) = diag( ν , . . . , ν k ) ⊕ diag( ν , . . . , ν k ) , (C10) S ( k ) = (cid:32) ⊕ kj =1 cosh r ( k ) j ⊕ kj =1 sinh r ( k ) j ⊕ kj =1 sinh r ( k ) j ⊕ kj =1 cosh r ( k ) j (cid:33) , (C11) where ν >
1, and ν j = 1 for 2 ≤ j ≤ k , and r ( k )1 > r ( k ) j = 0 for 2 ≤ j ≤ k . In other words, the reducedstate is a product of a thermal state and vacuum statesfollowed by a squeezing operation and beam splitter ar-rays. One can easily check that A ( k ) n is also rank-2. Thus,in this case, the computation of Hafnian is polynomial in M , so that the classical simulation using the Monte-Carlosimulation is efficient. Appendix D: Error of photon-number truncation
In this section, we analyze the influence of photon-number truncation in squeezed states. When we have N squeezed vacuum states and measure them in photon-number basis, the probability to generate a total of k photon pair events (2 k photons) is given by the negativebinomial distribution [10, 60] P N ( k ) = (cid:18) N + k − k (cid:19) sech N r tanh k r. (D1)Note that beam splitter arrays do not change the prob-ability distribution of total photon numbers. The tailprobability of the negative binomial distribution is givenby [61]Pr( k > αN sech r ) ≤ exp (cid:20) − αN (1 − /α ) (cid:21) . (D2)While increasing a constant α decreases the trunca-tion error exponentially, the order of the complexity O ( M n n max ) does not change. Any accuracy can beachieved by increasing α with a constant factor whichreduces error exponentially. In order to make the trun-cation error to be smaller than (cid:15) for sufficiently large α ,we may choose αN sech r = 2sech r log(1 /(cid:15) ) ≡ n max . Appendix E: Proof of Lemma 2
Proof.
Note that the covariance matrix V can be de-composed as V = S ( M / S T by a symplectic matrix S satisfying S Ω S T = Ω and that the symplectic matrixcan be decomposed as S = OS sq , where O representsa symplectic matrix corresponding to beam splitter ar-rays and S sq represents squeezing operators to generatesqueezed vacuum sources. If V = V + X and (cid:107) X (cid:107) (cid:28) V + V ) = det(2 V − X ) = det( M − S − XS − T ) ≤ (cid:18) M | Tr[ S − XS − T | ] (cid:19) M = M (cid:88) k =0 (cid:18) Mk (cid:19) (cid:18) M | Tr[ S − XS − T ] | (cid:19) k ≤ M (cid:88) k =0 (cid:18) Mk (cid:19) (cid:18) M (cid:107) X (cid:107) (cid:113) Tr( S − T S − ) (cid:19) k = 1 + (cid:107) X (cid:107)√ N cosh 4 r + O ( (cid:107) X (cid:107) ) , (E1) where for the first inequality, we have used the upper-bound of determinant and thatTr[ S − XS − T ] = Tr[ S − ( V − V ) S − T ]= Tr[ M / − S − V S − T ]= Tr[ M / − ˜ V ] ≤ . (E2)The Cauchy-Schwarz inequality has been used for thesecond inequality. Then, the quantum fidelity betweenGaussian states with covariance matrices V and V isapproximated as1 − F ( V , V ) = 1 − (cid:112) det( V + V ) ≤ (cid:107) X (cid:107)√ N cosh 4 r + O ( (cid:107) X (cid:107) ) . (E3) Appendix F: Proof of lemma 3
In this section, we compute (cid:107) X (cid:107) and prove Lemma 3. We first note that X is the difference matrix between twocovariance matrices V out and V a . (cid:107) X (cid:107) = N (cid:88) α =1 (cid:88) i,j ∈L α X ij + (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β X ij = 14 N (cid:88) α =1 (cid:88) i ∈L α ( δ (cid:104){ ˆ x i , ˆ x j }(cid:105) + δ (cid:104){ ˆ p i , ˆ p j }(cid:105) + 2 δ (cid:104){ ˆ x i , ˆ p j }(cid:105) )+ 14 (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β (cid:0) (cid:104){ ˆ x i , ˆ x j }(cid:105) + (cid:104){ ˆ p i , ˆ p j }(cid:105) + 2 (cid:104){ ˆ x i , ˆ p j }(cid:105) (cid:1) , (F1)where δ (cid:104){ ˆ O , ˆ O }(cid:105) ≡ (cid:104){ ˆ O , ˆ O }(cid:105) − (cid:104){ ˆ O , ˆ O }(cid:105) a with (cid:104)·(cid:105) and (cid:104)·(cid:105) a representing the expectation value for an exact andan approximate covariance matrix, respectively. We compute the elements one by one. First, we focus on the firstterms, where i, j are in the same sublattices. Usingˆ x i = 1 √ U ik ˆ a k + U ∗ ik ˆ a † k ) , and ˆ p i = i √ U ∗ ik ˆ a † k − U ik ˆ a k ) , (F2)with ˆ a k being an annihilation operator before beam splitter arrays, we obtain (cid:104){ ˆ x i , ˆ x j }(cid:105) = M (cid:88) k =1 (cid:104) U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k + U ik U ∗ jk ˆ a k ˆ a † k + U ∗ ik U jk ˆ a † k ˆ a k (cid:105) , (F3) (cid:104){ ˆ p i , ˆ p j }(cid:105) = − M (cid:88) k =1 (cid:104) U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k − U ik U ∗ jk ˆ a k ˆ a † k − U ∗ ik U jk ˆ a † k ˆ a k (cid:105) , (F4) (cid:104){ ˆ x i , ˆ p j }(cid:105) = i M (cid:88) k =1 (cid:104)− U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k (cid:105) . (F5)Here, note that the approximate covariance matrix is obtained by assuming other sublattices to be vacuum at thebeginning. Thus, for a fixed α , the difference between the exact and approximate covariance matrices is that for theapproximate covariance matrix, (cid:104) ˆ a k (cid:105) a = (cid:104) ˆ a † k (cid:105) a = (cid:104) ˆ a † k ˆ a k (cid:105) a = 0 and (cid:104) ˆ a k ˆ a † k (cid:105) a = 1 for k (cid:54)∈ L α but for the exact covariance2matrix, (cid:104) ˆ a k (cid:105) = (cid:104) ˆ a † k (cid:105) = cosh r sinh r , (cid:104) ˆ a † k ˆ a k (cid:105) = cosh r and (cid:104) ˆ a k ˆ a † k (cid:105) = sinh r for k ∈ S . Also, (cid:104) ˆ a k (cid:105) = (cid:104) ˆ a k (cid:105) a , (cid:104) ˆ a † k (cid:105) = (cid:104) ˆ a † k (cid:105) a , (cid:104) ˆ a † k ˆ a k (cid:105) = (cid:104) ˆ a † k ˆ a k (cid:105) a = 0 and (cid:104) ˆ a k ˆ a † k (cid:105) = (cid:104) ˆ a k ˆ a † k (cid:105) a for k ∈ L α . Therefore, the approximation error for a fixed α iswritten as | δ (cid:104){ ˆ x i , ˆ x j }(cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) k ∈M\L α (cid:104) U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k + U ik U ∗ jk ˆ a k ˆ a † k + U ∗ ik U jk ˆ a † k ˆ a k (cid:105) − (cid:88) k ∈M\L α U ik U ∗ jk (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F6)= (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) k ∈S\{ s α } (cid:0) U ik U jk cosh r sinh r + U ∗ ik U ∗ jk cosh r sinh r + U ik U ∗ jk sinh r + U ∗ ik U jk sinh r (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F7) ≤ (cid:88) k ∈S\{ s α } | U ik || U jk | e r sinh r ≤ e r (cid:88) k ∈S\{ s α } | U ik || U jk | . (F8)Now, we compute the sum of the absolute value for different i and j , N (cid:88) α =1 (cid:88) i,j ∈L α | δ (cid:104){ ˆ x i , ˆ x j }(cid:105)| ≤ e r N (cid:88) α =1 (cid:88) i,j ∈L α (cid:88) k,l ∈S\{ s α } | U ik || U jk || U il || U jl | (F9) ≤ e r N (cid:88) α =1 (cid:88) k,l ∈S\{ s α } (cid:115) (cid:88) i ∈L α | U ik | (cid:88) i ∈L α | U il | (cid:88) j ∈L α | U jk | (cid:88) j ∈L α | U jl | (F10) ≤ e r N (cid:88) α =1 (cid:88) k,l ∈S\{ s α } η = e r N ( N − η . (F11)Here, we have used the Cauchy-Schwarz inequality and defined η to be a leakage rate from a source to other sublattices.Following the same procedure, one can find the same upper-bound for momentum quadratures, N (cid:88) α =1 (cid:88) i,j ∈L α | δ (cid:104){ ˆ p i , ˆ p j }(cid:105)| ≤ e r N ( N − η . (F12)Similarly, we get | δ (cid:104){ ˆ x i , ˆ p j }(cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) k ∈M\L α (cid:104)− U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) k ∈M\L α ( − U ik U jk + U ∗ ik U ∗ jk ) cosh r sinh r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F13) ≤ (cid:88) k ∈S\{ s α } | U ik || U jk | cosh r sinh r ≤ e r (cid:88) k ∈S\{ s α } | U ik || U jk | . (F14)Again, using the Cauchy-Schwarz inequality, we obtain N (cid:88) α =1 (cid:88) i,j ∈L α | δ (cid:104){ ˆ x i , ˆ p j }(cid:105)| ≤ e r N ( N − η . (F15)Now, we obtain the upper-bound of the second terms in Eq. (F1), where i and j are in different sublattices Here,we approximate correlations as zero, which leads to | δ (cid:104){ ˆ x i , ˆ x j }(cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M (cid:88) k =1 (cid:104) U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k + U ik U ∗ jk ˆ a k ˆ a † k + U ∗ ik U jk ˆ a † k ˆ a k (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F16)= (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) k ∈S (cid:0) U ik U jk cosh r sinh r + U ∗ ik U ∗ jk cosh r sinh r + U ik U ∗ jk cosh r + U ∗ ik U jk sinh r (cid:1) + (cid:88) k (cid:54)∈S U ik U ∗ jk (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F17)= (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) k ∈S (cid:0) U ik U jk cosh r sinh r + U ∗ ik U ∗ jk cosh r sinh r + U ik U ∗ jk sinh r + U ∗ ik U jk sinh r (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F18) ≤ e r sinh r (cid:88) k ∈S | U ik || U jk | ≤ e r (cid:88) k ∈ S | U ik || U jk | . (F19)3After taking summation over different pairs of α and β , we obtain (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β | δ (cid:104){ ˆ x i , ˆ x j }(cid:105)| ≤ e r (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β (cid:88) k,l ∈S | U ik || U jk || U il || U jl | (F20) ≤ e r (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) k,l ∈S (cid:115) (cid:88) i ∈L α | U ik | (cid:88) i ∈L α | U il | (cid:88) j ∈L β | U jk | (cid:88) j ∈L β | U jl | (F21)= e r (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) k ∈S (cid:115) (cid:88) i ∈L α | U ik | (cid:88) j ∈L β | U jl | (cid:88) l ∈S (cid:115) (cid:88) i ∈L α | U il | (cid:88) j ∈L β | U jk | . (F22)Here, (cid:88) k ∈S (cid:115) (cid:88) i ∈L α | U ik | (cid:88) j ∈L β | U jk | = (cid:88) k ∈S\{ s α ,s β } (cid:115) (cid:88) i ∈L α | U ik | (cid:88) j ∈L β | U jk | + (cid:115) (cid:88) i ∈L α | U is α | (cid:88) j ∈L β | U js α | (F23)+ (cid:115) (cid:88) i ∈L α | U is β | (cid:88) j ∈L β | U js β | (F24) ≤ (cid:88) k ∈S\{ s α ,s β } ( η + 2 √ η ) = ( N − η + 2 √ η ) . (F25)Thus, the sum of the approximation errors for different pairs of α and β is bounded as (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β | δ (cid:104){ ˆ x i , ˆ x j }(cid:105)| ≤ e r N ( N − N − ( η + 2 √ η ) . (F26)Again, the same upper-bound can be found for momentum quadratures, (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β | δ (cid:104){ ˆ p i , ˆ p j }(cid:105)| ≤ e r N ( N − N − ( η + 2 √ η ) . (F27)Finally, the correlation of position and momentum quadratures is bounded as | δ (cid:104){ ˆ x i , ˆ p j }(cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M (cid:88) k =1 (cid:104)− U ik U jk ˆ a k + U ∗ ik U ∗ jk ˆ a † k (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = cosh r sinh r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) k ∈ S ( − U ik U jk + U ∗ ik U ∗ jk ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (F28) ≤ r sinh r (cid:88) k ∈S | U ik || U jk | ≤ e r (cid:88) k ∈S | U ik || U jk | , (F29)which leads to the sum of the errors as (cid:88) ≤ α (cid:54) = β ≤ N (cid:88) i ∈L α (cid:88) j ∈L β | δ (cid:104){ ˆ x i , ˆ p j }(cid:105)| ≤ e r N ( N − N − ( η + 2 √ η ) . (F30)Hence, we prove lemma 3 (cid:107) X (cid:107) ≤ e r (cid:2) N ( N − η + N ( N − N − ( η + 2 √ η ) (cid:3) ≤ e r N ( η + 2 √ η ) . (F31)4 Appendix G: Single-photon state
In this section, we analyze the error of an approxima-tion by treating single-photons as distinguishable parti-cles. From Eq. (22), (cid:15) = (cid:88) t t ! (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) σ (cid:54) = τ N (cid:89) k =1 U out σ ( k ) , in k (cid:32) N (cid:89) l =1 U out τ ( l ) , in l (cid:33) ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (G1) ≤ (cid:88) t (cid:88) σ (cid:54) = τ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:89) k =1 U out σ ( k ) , in k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:89) l =1 U out τ ( l ) , in l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (G2)= 12 (cid:88) t (cid:88) σ,ρ N (cid:89) k =1 | U out σ ( k ) , in k U out σ ( k ) , in ρ ( k ) | (G3)= 12 (cid:88) j (cid:88) ρ (cid:54) =Id (cid:89) i | U j i , in i || U j i , in ρ ( i ) | , (G4)where we have used | a + b | ≤ | a | + | b | for the first inequalityand rearranged the summation in the last equality. Here,the sum j is taken over ordered tuples ( j , . . . , j N ) with1 ≤ j l ≤ M for all 1 ≤ l ≤ N , and the sum ρ is takenover all possible permutations except for the identity Id.Now, we prove lemma 4. Let us find an upper-boundof C i : C i = (cid:88) j i ,ρ : ρ ( i ) (cid:54) = i | U j i ,i || U j i ,ρ ( i ) | = (cid:88) j : | i − j |≤ L/ (cid:88) k ∈ in ,k (cid:54) = i | U j,i || U j,k | + (cid:88) j : | i − j | >L/ (cid:88) k ∈ in ,k (cid:54) = i | U j,i || U j,k | . (G5)Here, a distance | i − j | is l ∞ distance between i th positionand j th position. Consider the first term of C i : (cid:88) j : | i − j |≤ L/ | U j,i | (cid:88) k ∈ in ,k (cid:54) = i | U j,k | (G6) ≤ (cid:88) j : | i − j |≤ L/ | U j,i | (cid:115) ( N − (cid:88) k ∈ in ,k (cid:54) = i | U j,k | (G7) ≤ (cid:88) j : | i − j |≤ L/ | U j,i | (cid:112) ( N − η (G8) ≤ (cid:115) (cid:88) j : | i − j |≤ L/ | U j,i | (cid:112) ηM ( N − /N (G9) ≤ (cid:112) ηM ( N − /N , (G10)where we have used the Cauchy-Schwarz inequality forthe first and second inequalities. Now, the second term: (cid:88) j : | i − j | >L/ | U j,i | (cid:88) k ∈ in ,k (cid:54) = i | U j,k | (G11) ≤ (cid:88) j : | i − j | >L/ | U j,i |√ N − ≤ (cid:115) ( M − N )( N − (cid:88) j : | i − j | >L/ | U j,i | (G13) ≤ (cid:112) η ( M − N )( N − , (G14)where we have used the Cauchy-Schwarz inequality forthe first and the second inequalities. Thus, C i ≤ (cid:112) ηM ( N − /N + (cid:112) η ( M − N )( N − ≤ (cid:112) ηM N . (G15)5 [1] P. W. Shor, Algorithms for quantum computation: dis-crete logarithms and factoring, in Proceedings 35th an-nual symposium on foundations of computer science (Ieee, 1994) pp. 124–134.[2] S. Lloyd, Universal quantum simulators, Science ,1073 (1996).[3] J. Preskill, Quantum computing in the NISQ era andbeyond, Quantum , 78 (2018).[4] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin,R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao,D. A. Buell, et al. , Quantum supremacy using a pro-grammable superconducting processor, Nature , 505(2019).[5] H. Wang, Y. He, Y.-H. Li, Z.-E. Su, B. Li, H.-L. Huang,X. Ding, M.-C. Chen, C. Liu, J. Qin, et al. , High-efficiency multiphoton boson sampling, Nat. Photonics , 361 (2017).[6] H. Wang, W. Li, X. Jiang, Y.-M. He, Y.-H. Li, X. Ding,M.-C. Chen, J. Qin, C.-Z. Peng, C. Schneider, et al. ,Toward scalable boson sampling with photon loss, Phys.Rev. Lett. , 230502 (2018).[7] H. Wang, J. Qin, X. Ding, M.-C. Chen, S. Chen, X. You,Y.-M. He, X. Jiang, L. You, Z. Wang, et al. , Boson sam-pling with 20 input photons and a 60-mode interferome-ter in a 10 -dimensional Hilbert space, Phys. Rev. Lett. , 250503 (2019).[8] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C.Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, et al. ,Quantum computational advantage using photons, Sci-ence , 1460 (2020).[9] S. Aaronson and A. Arkhipov, The computational com-plexity of linear optics, in Proceedings of the forty-thirdannual ACM symposium on Theory of computing (2011)pp. 333–342.[10] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen,C. Silberhorn, and I. Jex, Gaussian boson sampling,Phys. Rev. Lett. , 170501 (2017).[11] M. J. Bremner, R. Jozsa, and D. J. Shepherd, Classicalsimulation of commuting quantum computations impliescollapse of the polynomial hierarchy, Proceedings of theRoyal Society A: Mathematical, Physical and Engineer-ing Sciences , 459 (2011).[12] B. Fefferman and C. Umans, The power of quan-tum Fourier sampling, arXiv preprint arXiv:1507.05592(2015).[13] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush,N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, andH. Neven, Characterizing quantum supremacy in near-term devices, Nat. Phys. , 595 (2018).[14] G. Vidal, Efficient classical simulation of slightly entan-gled quantum computations, Phys. Rev. Lett. , 147902(2003).[15] J. Napp, R. L. La Placa, A. M. Dalzell, F. G. Bran-dao, and A. W. Harrow, Efficient classical simulationof random shallow 2d quantum circuits, arXiv preprintarXiv:2001.00021 (2019).[16] H. Qi, D. Cifuentes, K. Br´adler, R. Israel, T. Kala-jdzievski, and N. Quesada, Efficient sampling from shal-low Gaussian quantum-optical circuits with local inter-actions, arXiv preprint arXiv:2009.11824 (2020).[17] D. J. Brod, The complexity of simulating constant-depth bosonsampling, Phys. Rev. A , 042316 (2015).[18] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani,On the complexity and verification of quantum randomcircuit sampling, Nat. Phys. , 159 (2019).[19] A. Deshpande, B. Fefferman, M. C. Tran, M. Foss-Feig,and A. V. Gorshkov, Dynamical phase transitions in sam-pling complexity, Phys. Rev. Lett. , 030501 (2018).[20] N. Maskara, A. Deshpande, A. Ehrenberg, M. C. Tran,B. Fefferman, and A. V. Gorshkov, Complexity phasediagram for interacting and long-range bosonic hamilto-nians, arXiv preprint arXiv:1906.04178 (2019).[21] E. H. Lieb and D. W. Robinson, The finite group ve-locity of quantum spin systems, in Statistical mechanics (Springer, 1972) pp. 425–431.[22] B. Zhang and Q. Zhuang, Entanglement formation incontinuous-variable random quantum networks, arXivpreprint arXiv:2005.12934 (2020).[23] S. Aaronson and D. J. Brod, Bosonsampling with lostphotons, Phys. Rev. A , 012335 (2016).[24] J. Renema, V. Shchesnovich, and R. Garcia-Patron, Clas-sical simulability of noisy boson sampling, arXiv preprintarXiv:1809.01953 (2018).[25] M. Oszmaniec and D. J. Brod, Classical simulation ofphotonic linear optics with lost particles, New J. Phys. , 092002 (2018).[26] R. Garc´ıa-Patr´on, J. J. Renema, and V. Shchesnovich,Simulating boson sampling in lossy architectures, Quan-tum , 169 (2019).[27] H. Qi, D. J. Brod, N. Quesada, and R. Garc´ıa-Patr´on,Regimes of classical simulability for noisy Gaussian bosonsampling, Phys. Rev. Lett. , 100502 (2020).[28] C. Oh, K. Noh, B. Fefferman, and L. Jiang, Classicalsimulation of lossy boson sampling using matrix productoperators, arXiv preprint arXiv:2101.11234 (2021).[29] Q. Zhuang, T. Schuster, B. Yoshida, and N. Y. Yao,Scrambling and complexity in phase space, Phys. Rev.A , 062334 (2019).[30] The definition of the efficient sampler is weaker than theone that requires for total variance distribution to be (cid:15) -close in runtime poly( n, /(cid:15) ). The depth bound efficientlysimulated depends on the choice of the definition.[31] X.-B. Wang, T. Hiroshima, A. Tomita, and M. Hayashi,Quantum information with Gaussian states, Phys. Rep , 1 (2007).[32] C. Weedbrook, S. Pirandola, R. Garc´ıa-Patr´on, N. J.Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussianquantum information, Rev. Mod. Phys. , 621 (2012).[33] G. Adesso, S. Ragy, and A. R. Lee, Continuous vari-able quantum information: Gaussian states and beyond,Open Syst. Inf. Dyn. , 1440001 (2014).[34] A. Serafini, Quantum continuous variables: a primer oftheoretical methods (CRC press, 2017).[35] N. Quesada and J. M. Arrazola, Exact simulation ofGaussian boson sampling in polynomial space and ex-ponential time, Phys. Rev. Res. , 023005 (2020).[36] A. I. Barvinok, Two algorithmic results for the travelingsalesman problem, Mathematics of Operations Research , 65 (1996).[37] R. Kan, From moments of sum to moments of product,J. Multivar. Anal. , 542 (2008).[38] A. Bj¨orklund, B. Gupt, and N. Quesada, A faster hafnian formula for complex matrices and its benchmarking ona supercomputer, Journal of Experimental Algorithmics(JEA) , 1 (2019).[39] N. Quesada, J. M. Arrazola, and N. Killoran, Gaussianboson sampling using threshold detectors, Phys. Rev. A , 062322 (2018).[40] C. A. Fuchs and J. Van De Graaf, Cryptographic dis-tinguishability measures for quantum-mechanical states,IEEE Transactions on Information Theory , 1216(1999).[41] G. Spedalieri, C. Weedbrook, and S. Pirandola, A limitformula for the quantum fidelity, J. Phys. A Math. Theor. , 025304 (2012).[42] L. Banchi, S. L. Braunstein, and S. Pirandola, Quantumfidelity for arbitrary gaussian states, Phys. Rev. Lett. , 260501 (2015).[43] R. Movassagh, Quantum supremacy and random circuits,arXiv preprint arXiv:1909.06210 (2019).[44] A. Bouland, B. Fefferman, Z. Landau, and Y. Liu, Noiseand the frontier of quantum supremacy, arXiv preprintarXiv:2102.01738 (2021).[45] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani,Experimental realization of any discrete unitary opera-tor, Phys. Rev. Lett. , 58 (1994).[46] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S.Kolthammer, and I. A. Walmsley, Optimal design for uni-versal multiport interferometers, Optica , 1460 (2016).[47] N. J. Russell, L. Chakhmakhchyan, J. L. O’Brien, andA. Laing, Direct dialling of haar random unitary matri-ces, New J. Phys. , 033007 (2017).[48] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-robinson bounds and the generation of correlations andtopological quantum order, Phys. Rev. Lett. , 050401(2006).[49] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller,R. Blatt, and C. F. Roos, Quasiparticle engineering andentanglement propagation in a quantum many-body sys- tem, Nature , 202 (2014).[50] A. Nahum, S. Vijay, and J. Haah, Operator spreading inrandom unitary circuits, Phys. Rev. X , 021014 (2018).[51] X. Mi et al. , Information scrambling in computa-tionally complex quantum circuits, arXiv preprintarXiv:2101.08870 (2021).[52] D. A. Roberts and B. Swingle, Lieb-robinson bound andthe butterfly effect in quantum field theories, Phys. Rev.Lett. , 091602 (2016).[53] M. C. Tran, C.-F. Chen, A. Ehrenberg, A. Y. Guo,A. Deshpande, Y. Hong, Z.-X. Gong, A. V. Gorshkov,and A. Lucas, Hierarchy of linear light cones with long-range interactions, Phys. Rev. X , 031009 (2020).[54] J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi,Quantum state transfer and entanglement distributionamong distant nodes in a quantum network, Phys. Rev.Lett. , 3221 (1997).[55] V. Khemani, A. Vishwanath, and D. A. Huse, Opera-tor spreading and the emergence of dissipative hydrody-namics under unitary evolution with conservation laws,Physical Review X , 031057 (2018).[56] T. Rakovszky, F. Pollmann, and C. Von Keyserlingk, Dif-fusive hydrodynamics of out-of-time-ordered correlatorswith charge conservation, Physical Review X , 031058(2018).[57] H. Kim and D. A. Huse, Ballistic spreading of entan-glement in a diffusive nonintegrable system, Phys. Rev.Lett. , 127205 (2013).[58] S. Xu and B. Swingle, Accessing scrambling using matrixproduct operators, Nat. Phys. , 199 (2020).[59] A. Botero and B. Reznik, Modewise entanglement ofGaussian states, Phys. Rev. A , 052311 (2003).[60] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen,C. Silberhorn, and I. Jex, Detailed study of Gaussianboson sampling, Phys. Rev. A , 032326 (2019).[61] D. G. Brown, How I wasted too long finding a con-centration inequality for sums of geometric variables, https://cs.uwaterloo.ca/~browndg/negbin.pdfhttps://cs.uwaterloo.ca/~browndg/negbin.pdf