Can lossy Gaussian boson sampling be simulated beyond the square-root regime?

The square-root regime is a limit on one method

The setup is standard lossy Gaussian boson sampling. N squeezed vacuum states with a common fixed squeezing magnitude r, plus M−N vacuum modes, pass through uniform pure loss with photon survival probability η and a passive interferometer V, then reach photon-number-resolving detectors. The target P is the full photon-count distribution over all M output modes, and accuracy is measured by total variation distance TV(P,Q), the standard half-sum of absolute count-probability differences. The mean surviving photon number is N_out = Nη·sinh²r.

A widely used classical approach replaces each lossy input with the closest classical Gaussian state and treats the rest exactly. That method stays accurate while N_out grows no faster than the square root of N at fixed squeezing and accuracy. The paper's first result shows this is not an artifact of choosing a Gaussian reference. Theorem 1 constructs explicit interferometers — pairs of equally squeezed inputs sent through a beam splitter, with the remaining outputs left in vacuum — for which every state with a nonnegative Glauber-Sudarshan P representation, Gaussian or not and however correlated across modes, has TV error at least 1 − exp(−c_r Nη²). On that family the optimal error within the P-classical set is Θ_r(min{1, Nη²}).

The mechanism is a sign constraint on classical count probabilities. For a product coherent state with intensities I and J, the pair probabilities satisfy q₁₁ − q₂₀ − q₀₂ = −½e^(−I−J)(I−J)² ≤ 0, and every P-classical mixture inherits this inequality. The lossy pair violates it, with p₁₁ − p₂₀ − p₀₂ = n₀η²/[1 + n₀(2η − η²)]² > 0 at order η². Summing the constraint across pairs and applying Cauchy-Schwarz turns this local violation into the global lower bound.

The scope matters as much as the bound. Theorem 1 is a statement about a method, not a complexity lower bound: it does not cover references that are Gaussian but nonclassical, nor the nonclassical core retained by tensor-network approaches, and it says nothing about the hardness of the distribution itself. What it rules out is any algorithm that replaces the input by classical light and then detects.

Approximating detection escapes the barrier

Because the barrier is a property of the approximation class rather than of classical computation, a method that keeps the input and approximates the detection is not bound by it. The paper rewrites photon-number-resolving detection as a fictitious absorption process: each output mode is coupled to a perfect absorber that removes photons one at a time and records which mode each removal came from. Running until every photon is absorbed, the number of removals recorded in mode i is exactly its photon number, so the record reproduces the PNR distribution. The fictitious time only orders the clicks and is integrated out; nothing is applied in the experiment.

Along this process the conditional state follows two rules. While no click occurs it is damped by e^(−t n̂/2), where n̂ is the total photon-number operator — the state conditioned on no photon having been removed yet. At a click in mode i it becomes b_i ρ b_i†/λ_i, a photon subtraction, with λ_i = ⟨b_i† b_i⟩ the rate of that mode. Neither rule is unitary. For a Gaussian state both have closed forms: with Σ = I/2 + K, the no-click evolution gives K̇ = −K − K², and a click updates K → K + Δ_i with Δ_i = KΠ_iK/λ_i and λ_i = ½Tr(Π_iK). The waiting time to the next click also has a closed form, S_K(t) = det[I + (1 − e^(−t))K]^(−1/2), whose limit det(I + K)^(−1/2) is the vacuum probability of the current state and therefore the probability that no photon is left.

The sampler stores only K. Starting from K₀ = ηO diag(p_r I_N, −d_r I_N)Oᵀ, where O is the realification of V and p_r, d_r are fixed by r, it draws whether the record has ended and otherwise the time of the next click, chooses the clicked output in proportion to the rates λ_i, and applies the update. The clicks per output form the sample. The no-click step is exact; the single approximation is that after a click the exact state, which is no longer Gaussian, is replaced by the centered Gaussian with the same covariance I/2 + K + Δ_i. No hafnian is needed, and neither photon number nor covariance is capped.

The error scales as Nη³, not Nη²

Theorem 2 states that for every fixed r > 0, every M ≥ N, and every passive isometry V, there is a randomized classical algorithm with polynomial expected arithmetic cost whose output distribution Q satisfies TV(P, Q) ≤ min{1, C₃(r)Nη³} for 0 ≤ η ≤ 1, with C₃(r) independent of M, N, and V. The theorem applies to the squeezed and vacuum inputs specified above, with no restriction on interferometer depth or connectivity, and output collisions are included.

The cubic scale has a simple reading. At low transmission each output holds O(η) photons, so a subtraction changes the state by Δ_i = O(η), and along typical records the Gaussian replacement costs O(η²) per click. Clicks occur at a total rate of order Nη, and the product is Nη³. The proof also covers rare records with large covariance changes, using a shift Y_t = K_t + b e^(−t)P with b = η/(1−η) and P = OOᵀ, and a nonnegative functional Φ_t = Tr(Y_t K_t²). The inequalities Y_t ⪰ 0, 0 ⪯ Δ_i ⪯ 2Y_t, and Σ_i λ_iΔ_i = K_t² rely on the paired quadrature structure of squeezed inputs rather than on a lower bound on any output's occupation.

The cost is polynomial but not small. A direct implementation has expected cost O(MN + (1 + Nηn₀)(N³ + MN²)) and memory O(MN + N²), with scalar inversion of S_K counted as a primitive, and the expected number of clicks is at most Nηn₀. The factor N rather than M in the error bound comes from the occupied covariance subspace.

What Nη³ buys in surviving photons

Setting Nη³ = O(1) at fixed squeezing and accuracy permits η = c₀N^(−1/3) for a sufficiently small fixed c₀, and therefore a mean surviving photon number N_out proportional to N^(2/3). That is the headline improvement over the square-root regime, where Nη² = O(1) forces N_out = O(√N).

The two theorems separate on the same target family. Choose η = N^(−a) with 1/3 < a < 1/2. Then Nη² → ∞ while Nη³ → 0, so every classical-light count distribution becomes asymptotically maximally distant from the target while the detection-approximation sampler stays accurate. This is a direct separation on one family of interferometers, not a claim that all lossy GBS instances are classically simulable at this loss.

A separate pair construction reaches further. Using pair generation and amplification with an exact treatment of a few pair sectors, the total variation distance is bounded by (3/4)√(N x⁴) with x = (n₀ + 1)η, so for fixed r and accuracy η proportional to N^(−1/4) gives a mean surviving count proportional to N^(3/4). This construction applies when every input mode is equally squeezed. The paper is explicit that the finitely-many-exact-sectors bound does not by itself establish the N^(3/4) exponent, because the core pair number grows with N in that regime.

Where the detection approximation stops working

The Gaussian reset is not uniformly small. If the displacement intensity |α_i|² is comparable to or smaller than the core occupation n_c,i, the reset can have order-one error: at zero displacement it replaces S(u)|1⟩ by a squeezed thermal state at trace distance 3/4. In the initial displacement ensemble this region has probability of order √η and the detection rate there is of order η², suggesting a contribution of order η^(5/2), but the paper states that these initial rates do not bound the full trajectory error.

A displaced variant Q′_G combines the exact displacement decomposition with Gaussian resets, retaining the sampled displacement and matching both mean and covariance after each click. Theorem 2 applies to the covariance-only sampler Q_G, not to this variant. The negative term −v_i v_iᵀ/λ_i² in the conditional-moment update prevents direct use of the relative jump bound and cubic budget for the displaced case, and no error analysis is given for a two-quadrature-noise variant that might remove the restriction to one noise quadrature.

The numerical tests are correspondingly narrow. They use one or two squeezed inputs in up to three modes, with r in {0.3, 0.5, 1} and η in {0.1, 0.2, 0.3}. At η = 0.5, Q_G becomes less accurate than the classical reference above r ≈ 1.1. For random three-mode interferometers at η = 0.3, the empirical errors of the reference, Q_G, Q′_G, and Q_pair are (0.058, 0.058, 0.014, 0.015) at r = 1 and (0.062, 0.136, 0.04, 0.080) at r = 1.5. For Q′_G both empirical distances are comparable to sampling noise and do not resolve the true error. These few-mode tests do not establish a uniform error bound and do not compare runtime with tensor-network methods.

The sampler trades accuracy for a lower mean

One numerical check shows the approximation failing in a specific, measurable way. For one mode at r = 1 and η = 0.3, a run with 10⁶ records gives mean 0.3795 ± 0.0010, compared with the exact 0.4143, and probability 1.3 × 10⁻⁴ for more than 16 clicks, compared with 2 × 10⁻⁶ for the target. The smaller sampled mean is consistent with the theoretical result, which gives only an upper bound on the expected count rather than equality. The paper notes that these numerical values do not establish failure of uniform integrability.

The formulas themselves were checked independently of the numerical tests. The reduced covariance representation, summed jump identity, cubic generator, and local cost inequality were verified in exact rational arithmetic for 18 Gaussian prefixes and 54 detector jumps, with 42 noncommuting jumps and complex rational isometries at (N, M) = (1, 2), (2, 3), and (2, 4). In 18 single-mode Fock calculations, increasing the dimension from 40 to 56 changes the local trace distance by at most 1.4 × 10⁻¹⁶. The analytical bounds are proved independently of these checks, which do not rigorously bound truncation or rounding errors.

What decides whether this method applies

The practical question for anyone assessing a lossy GBS advantage claim is which approximation is being made, because the two choices have different ceilings. If the method replaces the lossy input with classical light, the relevant threshold is Nη² = O(1), and the mean surviving photon number cannot exceed O(√N) in the worst case on the explicit interferometer family, regardless of whether the classical reference is Gaussian or non-Gaussian and correlated. If the method keeps the input exact and approximates the detection, the threshold becomes Nη³ = O(1), permitting N_out ∝ N^(2/3), and N^(3/4) for the equally-squeezed pair construction.

The separation is sharpest in the window 1/3 < a < 1/2 with η = N^(−a), where the classical-light error approaches its maximum while the detection-approximation error vanishes. That window, not the asymptotic exponents alone, is where the two approaches can be told apart on the same target family.

Two boundaries should travel with any use of these results. The lower bound is a statement about the classical-light replacement method, not a complexity lower bound, so it does not cover Gaussian-but-nonclassical references or retained nonclassical cores. And the upper bound is proved for the covariance-only sampler Q_G; the displaced variant Q′_G, which showed lower empirical error in the small tests, has no corresponding theorem. The most decision-relevant comparison in the paper is therefore not the exponent alone but the pairing of threshold and method: Nη² for input replacement, Nη³ for detection approximation, with the N^(3/4) pair construction restricted to equally squeezed inputs.

Sources used in this answer