Gaussian Boson Sampling

• Gaussian Boson Sampling is a quantum computational model that employs multimode squeezed vacuum states and linear optics to produce output distributions governed by matrix hafnians.
• Its computational hardness arises from the #P-hardness of computing hafnians, positioning GBS as a strong candidate for demonstrating quantum advantage in photonic devices.
• Experimental implementations, such as the Jiuzhang 3.0 platform, validate GBS through scalable architectures, rigorous statistical certification, and applications in graph theory and combinatorial optimization.

Gaussian Boson Sampling (GBS) is a non-universal model of quantum computation in which multimode squeezed vacuum states are propagated through a passive linear-optical network and measured by photon-number-resolving or threshold detectors. The output click patterns are distributed according to probabilities that are proportional to matrix hafnians, rendering classical simulation #P-hard in general. GBS is a leading candidate for demonstrating quantum advantage in near-term photonic devices, with recent experiments reporting photon detection rates and circuit sizes that surpass classical simulability. The architecture’s statistical foundations and complexity are now well-understood, including exact formulas for output probabilities, rigorous hardness proofs, links to graph-theoretic problems, and recent analyses of scalability under photon loss.

1. Physical Principles and Mathematical Formalism

GBS prepares $M$ optical modes in single-mode squeezed vacuum states with squeezing parameters $\{r_j\}$, which are then interfered in an $M$-port linear network represented by a unitary $U$. After transmission, the output quantum state is characterized by a zero-displacement Gaussian covariance matrix $\sigma$ (Kruse et al., 2018). The output detection is performed either by photon-number-resolving (PNR) or threshold (on–off) detectors.

The probability of observing a photon number pattern $\vec{n}$ can be written in closed form as

$$\operatorname{Pr}(\vec{n}) = \frac{1}{\prod_j n_j ! \sqrt{\det \sigma_Q}} \, |\operatorname{Haf}(A_{\vec{n}})|^2,$$

where $\sigma_Q = \sigma + \frac{1}{2} I_{2M}$ and $A_{\vec{n}}$ is a $2N \times 2N$ submatrix of state-dependent $A$ obtained by selecting rows/columns according to the detected photon pattern (Hamilton et al., 2016, Kruse et al., 2018, Zhong et al., 2019). For threshold detectors, the probability distribution is governed by the Torontonian of principal submatrices of $A$ (Quesada et al., 2018), which plays an analogous role to the hafnian in the collision-free regime.

In the case of pure SMSV and Haar-random $U$, the matrix $A$ block-diagonalizes, simplifying probability computations:

$$\operatorname{Pr}(\vec{n}) = \frac{|\operatorname{Haf}(B_S)|^2}{\prod_j n_j ! \sqrt{\det \sigma_Q}},$$

with $B = U \operatorname{diag}(\tanh r_j) U^T$ (Kruse et al., 2018, Hamilton et al., 2016).

2. Computational Complexity and Scalability

The computational hardness of GBS is attributed to the #P-complete complexity of the hafnian function. Calculating hafnians of general $2N\times2N$ matrices is exponentially hard, matching the permanent in scaling (Hamilton et al., 2016, Kruse et al., 2018). Absolute hardness is retained in the dilute regime where the number of modes scales quadratically with the photon number ($M = O(N^2)$) (Grier et al., 2021). Furthermore, approximate GBS remains #P-hard under plausible anti-concentration hypotheses for hafnians (Kruse et al., 2018).

Classical sampling difficulty has been quantitatively benchmarked: on current supercomputers, simulation is tractable up to $\sim 18$–$20$ photons, with $\sim 30$ photons marking the lower bound for quantum supremacy demonstration in GBS (1908.10070). The Titan supercomputer can, for example, simulate 800 modes postselected on 20 detector clicks, requiring about two hours per sample and massive parallel resources (Gupt et al., 2018). The best exact algorithms scale as $O(2^{N/2})$ in the number of photons detected (Kruse et al., 2018, 1908.10070).

Loss has pivotal impact: under realistic photon loss, if the number of surviving photons scales sublinearly ($N_{out} \propto \sqrt{N}$), the operator entanglement entropy (OEE) saturates at $O(\log N)$, rendering tensor-network simulation polynomial-time at fixed error (Liu et al., 2023). Conversely, quantum advantage is only possible if loss is sufficiently bounded so that $N_{out}$ grows at least as $O(N^{>1/2})$; otherwise, scalable classical simulation is feasible.

3. Experimental Implementations and Validation

GBS experiments utilize parametric down-conversion (PDC) sources for SMSV states, passive interferometric networks (chip or fiber-based), and high-efficiency detectors. Recent experimental advances include the Jiuzhang 3.0 platform, which employs pseudo-photon-number-resolving detection across $144$ modes and achieves photon click counts up to $255$ (Deng et al., 2023). Detector responses are calibrated to reconstruct effective POVMs for observed click patterns, with partial photon distinguishability and experimental loss rigorously modeled.

Experimental validation employs statistical hypothesis testing (Bayesian log-likelihood ratios, correlation function analysis) to distinguish genuine GBS output from classical mockups (thermal, squashed, distinguishable, etc). High-order click cumulants and coarse-grained feature vectors across sampling orbits serve as empirical signatures of quantum interference and computational hardness (Deng et al., 2023, Giordani et al., 2022). For instance, recent experiments achieve sampling rates and complexity levels where each sample, if classically simulated, would require hundreds to billions of years on the fastest available hardware (Deng et al., 2023).

A table summarizing experimental GBS performance benchmarks:

Platform	Modes	Max clicks	Sampling time (GBS)	Classical sim time/sample	Validation Protocol
Jiuzhang 3.0	144	255	1.27 μs	600 years (avg), $3.1\cdot10^{10}$ years (max)	Bayesian, cumulant analysis
Titan (sim)	800	20	—	2 hours	Memory, walltime scaling
Line-chip (exp)	12	5	23 kHz (5-ph event)	—	Distribution similarity, TVD

Recent experiments have also demonstrated efficient in situ Gaussian-state tomography via displacement, enabling full reconstruction of the covariance matrix and displacement vectors using only single and two-photon statistics (Thekkadath et al., 2022).

4. Connections to Graph Theory, Optimization, and Applications

GBS naturally encodes graph combinatorics: for a graph with adjacency matrix $A$, the output pattern probabilities are proportional to $|\operatorname{Haf}(A_S)|^2$, which counts perfect matchings in the corresponding subgraphs (Giordani et al., 2022, Deng et al., 2023). GBS enables proportional sampling for graph problems such as dense $k$-subgraph identification, where samples preferentially select highly interconnected subgraphs.

Hybrid quantum-classical algorithms leveraging GBS samples—e.g., random search, simulated annealing—demonstrate enhanced performance for NP-hard combinatorial optimization, outperforming uniform or thermal sampling in the number of steps required to identify optimal graph structures (Deng et al., 2023, Yin et al., 2022, Sempere-Llagostera et al., 2022). Dense subgraph certification has been achieved both in spatial-mode and time-bin-encoded interferometers, with time-bin architectures showing practical scaling to larger networks (Sempere-Llagostera et al., 2022).

GBS is also applicable to molecular vibronic spectra simulation, with output statistics mapping directly onto molecular transitions under suitable encoding.

5. Simulation Methodologies: Tensor Networks, Phase Space, Threshold Detectors

Classical simulation of GBS leverages different approaches based on the output regime and hardware parameters. Matrix Product Operator (MPO) representations of the GBS density matrix enable simulation under high loss using tensor networks, exploiting $U(1)$ symmetry to reduce memory cost and computational scaling (Liu et al., 2023). Bond dimension growth and entanglement entropy scaling dictate the tractability: polynomial scaling is observed for $N_{out}\propto\sqrt{N}$, while exponential scaling returns for higher $N_{out}$.

The positive-P and Wigner phase-space methods enable efficient computation of joint and marginal click probabilities, correlation functions, and multipartite entanglement in GBS networks with up to tens of thousands of modes (Drummond et al., 2021). These phase-space techniques bypass explicit hafnian calculation and compute grouped statistics relevant for benchmarking and certification.

Sampling with threshold detectors is governed by the Torontonian, a matrix function related to an infinite sum over hafnians (Quesada et al., 2018, Gupt et al., 2018). Exact classical algorithms for threshold GBS scale exponentially in the number of clicks but only polynomially in the number of modes, with maximum simulated clicks limited by available memory.

6. Practical Limitations: Loss, Detector Resolution, Noise

Practical scalability is predominantly limited by photon loss, detector inefficiency, partial indistinguishability, and finite photon-number resolution. The central limitation for quantum advantage is the effective transmission: unless per-mode loss is drastically reduced, the surviving photon number cannot outpace the classical simulation bottleneck (Liu et al., 2023).

Realistic detector models (finite resolution, saturation effects, dead-time) require generalized matrix functionals of the covariance matrix for accurate probability calculation, interpolating between the hafnian and Torontonian (Yeremenko et al., 5 Mar 2024). Validation protocols using orbits, Bayesian likelihoods, and $\chi^2$ statistics remain robust provided detector POVMs are faithfully incorporated.

Noise (thermal admixture, mode mismatch, phase drift) degrades quantum interference and may cause the device to emulate classical statistical mixtures if not rigorously suppressed. Experiments have shown high sampling fidelity ($>0.97$ similarity for small $N$, TVD $<0.1$) for up to $N=5$ photons and moderate loss (Zhong et al., 2019, Thekkadath et al., 2022).

7. Extensions, Alternative Models, and Outlook

GBS encompasses generalizations such as BipartiteGBS, which enables programming photonic circuits to sample output probabilities proportional to permanents of arbitrary matrices (Grier et al., 2021). This closes major complexity-theoretic loopholes and establishes hardness in the quadratic mode-number regime. Alternative measurement protocols using continuous-variable Gaussian measurements have been proven #P-hard under time-reversal symmetries (Chakhmakhchyan et al., 2017).

Current open directions include extending complexity proofs to sub-quadratic mode regimes, developing robust classical spoofing resistance under experimentally realistic noise, integrating error correction, and finding new application areas in optimization, chemistry, and machine learning. Certification via graph-theoretic feature statistics and kernel-based methods is seen as critical for cryptographically secure benchmarking of large-scale GBS devices (Giordani et al., 2022).

In summary, GBS represents a central paradigm in photonic quantum sampling, with rigorous mathematical foundations, scalable experimental architectures, and diverse applications spanning quantum advantage, graph theory, and molecular simulation. The model’s complexity, loss tolerance, simulation methodologies, and validation protocols continue to be refined as next-generation quantum photonic hardware advances (Liu et al., 2023, Deng et al., 2023, 1908.10070, Giordani et al., 2022, Yeremenko et al., 5 Mar 2024).