Polynomial-time classical simulation of boson and fermion sampling under thermal loss

The paper in detail

There are \(m\) bosonic or fermionic modes, a passive interferometer \(U\), and occupation-number measurements. Each mode is coupled to its own thermal environment with mean occupation \(\beta\) by \(\hat a_i\mapsto\sqrt\eta\,\hat a_i+\sqrt{1-\eta}\,\hat b_i\). After averaging, number operators transform as \(\hat n_i\mapsto\eta\hat n_i+(1-\eta)\beta\). Uniform loss commutes with passive optics and with dephasing in the occupation basis, so loss after each layer equals a single loss with transmissivity \(\eta_0^L\) at the end.

The output distribution is expanded in occupation polynomials orthogonal for a reference distribution \(\tau\): the uniform distribution over \(n\)-particle configurations, or a thermal product state for Gaussian boson sampling. The bases are Young polynomials for fermions, Hahn polynomials for Fock-state bosons and Meixner polynomials for Gaussian boson sampling. A polynomial's degree matches the irreducible components of the adjoint action of passive transformations, and degree \(k\) carries \(k\)-particle correlations. Loss multiplies the degree-\(k\) part by a factor at most \(q^k\) with \(q<1\), so truncating at degree \(D=O(\log m)\) costs little when the output is flat.

Each retained degree-\(k\) contribution reduces to an overlap on the \(k\)-particle space, computable exactly in time \(m^{O(k)}\). That gives the quasipolynomial theorem. For Haar-random interferometers, unbiased randomised estimators adapted from colour coding cost \(\operatorname{poly}(m)\,2^{O(k)}\), with variances controlled on average over \(U\), which gives polynomial time. Rejection sampling from the noisy reference, which is easy to sample, turns the signed approximation into samples.

With only a bounded number of empty input modes, the reduced state of any fixed number of particles approaches the maximally mixed state. The noisy output then concentrates around the loss applied to the uniform occupation state, and the degree-zero approximation alone suffices for large systems. The authors leave open peaked output distributions, the optimal classical cost at finite size, and matching lower bounds.

Transferable methods

1. Loss damps high-degree correlations

   Loss thins each mode's occupation number, and a \(k\)-particle correlation survives only if all \(k\) particles do. In a basis graded by degree, noise therefore acts diagonally with factors shrinking geometrically in \(k\), and truncation error is controlled by how flat the ideal output is.

   The polynomials satisfy \(\operatorname{Tr}(\tau X_iX_j)=\delta_{ij}h_j\), and noise preserves this orthogonality with factors \(\Gamma_j\le q^{\deg X_j}\). By Parseval, \(\chi^2(\rho_U\|\tau)=\sum_{j\ge1}|\operatorname{Tr}(\rho_UX_j)|^2/h_j\). Keeping degrees up to \(D\) before the noise channel gives \[\bigl\|\mathcal N(\rho_U)-\widetilde\rho_{D,U}\bigr\|_1^2\le\sum_{\deg X_j>D}\Gamma_j\frac{|\operatorname{Tr}(\rho_UX_j)|^2}{h_j}\le q^{D+1}\chi^2(\rho_U\|\tau),\] so \(D=O\bigl(1+\log[(1+C)/\varepsilon^2]/\log(1/q)\bigr)\) suffices when \(\chi^2\le C\).

   The paper's Table 1 constants give \(\kappa_{\rm F}\le\eta\) for fermions, \(q_{\rm B}\le2\eta\) for Fock-state bosons and \(q\le(1+a)\eta\) for Gaussian boson sampling, with \(q_{\rm B}=s(m+n)/(m+sn)<1\) in general. Editorial check of the mechanism: for pure loss, binomial thinning multiplies the factorial moments \(\mathbb E[(n)_k]\) by exactly \(\eta^k\), which a direct computation confirmed for \(k=1,2,3\) at \(\eta=0.3\) and \(0.7\).

   Editorial toy check, not a test of the theorem: for a Haar-random interferometer with \(m=5\) modes and \(n=3\) photons under pure loss, the total-variation distance between the noisy output and the noisy uniform reference, which is the degree-zero approximation, was 0.326 at \(\eta=1\), 0.232 at 0.8, 0.135 at 0.5, 0.083 at 0.3 and 0.030 at 0.1. Interference information fades steadily as loss grows.

   How to reuse it

   For a noisy sampling task, look for a basis graded by correlation order in which the noise acts diagonally and contracts. Low-degree truncation plus a flatness bound then gives classical simulation, exactly as Pauli weight does for noisy random circuits.

2. From a signed approximation to samples

   The truncated distribution can be negative, so it cannot be sampled directly. The paper instead proposes configurations from the noisy reference distribution and accepts them according to the approximation's relative density, which only requires evaluating the approximation at single outcomes.

   Let \(\mu=\mathcal N^{\otimes m}(\tau)\) and \(f_{D,U}(x)=\langle x|\widetilde\rho_{D,U}|x\rangle/\mu(x)\). Draw \(x\sim\mu\) and accept with probability \(\min\{1,\max\{0,f_{D,U}(x)\}/M\}\), where \(M=(1+C)/\varepsilon\). Lemma D.1 bounds the total-variation error by the approximation error plus clipping and stopping errors, each \(O(\varepsilon)\) once the number of proposals is capped.

   Evaluating \(f_{D,U}\) needs the degree-\(k\) overlaps \[p_{k,U}(x)\propto\operatorname{Tr}\bigl[\rho_x^{(k)}U_kT_kU_k^\dagger\bigr],\qquad T_k=\Pi_k\bigl(\rho_{\rm in}^{(k)}\bigr),\] which live on the \(k\)-particle space. Exact tables cost \(m^{O(k)}\). Random operators with \(\mathbb E\widehat\rho_{\rm in}^{(k)}=\rho_{\rm in}^{(k)}\) and \(\mathbb E\widehat\rho_x^{(k)}=\rho_x^{(k)}\) give unbiased estimates at cost \(\operatorname{poly}(m)\,2^{O(k)}\), computed with determinant formulas.

   Editorial remark: this is the whole difference between the two main theorems. With \(D=O(\log m)\), exact evaluation costs \(m^{O(\log m)}\), which is quasipolynomial, while the randomised cost \(2^{O(D)}\) is polynomial. Haar averaging is needed only to control the estimators' variance.

   How to reuse it

   Rejection sampling with clipping turns any pointwise-computable signed approximation that is close in \(L^1\) into a sampler. The same wrapper works for Pauli-path, Majorana-propagation or tensor-network approximations of noisy circuits.

Source: Section V.1-V.5 and Appendix H.3 (Table 1). Based on the linked paper version. Derivations marked as editorial, and suggested uses in other problems, are editorial guidance rather than results of the paper.

Context & limitations

Why it matters. The results suggest that superpolynomial advantage from Haar-random sampling needs a loss rate that vanishes as the system grows, even when a constant fraction of photons is detected and collisions dominate. One framework covers bosons and fermions, with representation theory identifying the "degree" that loss suppresses. For fermions, combined with Mathur and Kerenidis's anticoncentration bound, it makes logarithmic-depth random fermionic circuits with per-layer loss classically simulable.

Limitations. The guarantees are asymptotic. The runtime exponents grow as the loss rate falls, so substantial polynomial advantage at finite size is not excluded. Haar randomness or flatness is required, and peaked distributions, including recent proposals for verifiable advantage, fall outside the bounds. Only homogeneous thermal loss is treated, not partial distinguishability or detector noise. Concurrent work by Go and Jeong and by Lim overlaps. The authors state that ChatGPT (GPT-5.2, 5.4, 5.5, 5.6 Sol and 6 Astra) was used from January to October 2026 to help identify the polynomial basis and literature, prove technical results and draft parts of the manuscript, and that they devised the proof strategy and verified all derivations.

Description based on the full arXiv HTML of v1, read 2026-11-06: Sections I to V (framework, theorems, outlook and methods), Table 1 and the constants in Appendix H.3, the note added and the AI statement. The model-specific proofs in Appendices B to G were not checked.