Boson Sampling
The KLM protocol proved photonic quantum computing possible and priced it somewhere past the moon.
In 2011, Scott Aaronson and Alex Arkhipov asked the opposite question: strip away everything
expensive — the feedforward, the adaptive measurements, the entangling gates, even the qubits —
and what is the cheapest photonic experiment that still does something no classical
computer can? Their answer is boson sampling: shine n
single photons into an m-mode interferometer mesh, fixed in advance, and
record which detectors click. That's it. No logic, no answer to any question you asked — just
samples from a distribution. And yet, under complexity assumptions most theorists would bet their
houses on, simulating those samples classically is intractable, because hiding inside every click
pattern's probability is one of the meanest objects in computer science: the
permanent of a complex matrix.
The machine, and where the permanent comes from
The setup is Module 4's programmable mesh, verbatim: m waveguides
(typically m \gtrsim n^2), a Haar-random unitary
U compiled into MZIs, single photons injected into the first
n inputs, photon counters on every output. Each photon evolves by the
linear rule \hat a_j^\dagger \to \sum_k U_{kj}\hat a_k^\dagger — so the
n-photon input state becomes a product of
n such sums. Expand that product and collect the terms landing photons
on a chosen set of output modes S = \{s_1, \dots, s_n\}: one term for
every way of pairing input photons with output slots — every permutation
\sigma — each contributing the product of its matrix elements. Because
the photons are indistinguishable bosons, these n! histories add as
amplitudes with no signs:
A(S) \;=\; \sum_{\sigma \in S_n} \prod_{k=1}^{n} U_{s_k,\,\sigma(k)}
\;=\; \operatorname{per}\bigl(U_S\bigr),
\qquad
P(S) \;=\; \bigl|\operatorname{per}(U_S)\bigr|^2
(for collision-free outcomes; a doubly occupied mode divides by s_j!).
Here U_S is the n \times n submatrix of
U picking the occupied output rows and input columns. This is
Hong–Ou–Mandel
interference, industrialised: for n = 2 the permanent has
two terms and their cancellation is the HOM dip; for n = 30 it has
30! \approx 2.7 \times 10^{32} interfering histories.
- the permanent of an n \times n matrix is
\operatorname{per}(M) = \sum_{\sigma} \prod_k M_{k,\sigma(k)} —
the determinant's formula with every sign set to +1;
- dropping the signs is computationally catastrophic: the determinant is polynomial
(O(n^3) by elimination), but computing the permanent is
#P-hard even for 0/1 matrices (Valiant, 1979) — harder than anything in NP;
- the best general algorithm known, Ryser's formula, needs
O(2^n n) operations — about 10^{11}
for n = 30, and hopeless by n \approx 60;
- boson-sampling outcome probabilities are squared permanents of complex
n \times n submatrices of the mesh unitary.
Compute one yourself
The permanent's menace is easiest to respect after implementing it. The program below evaluates
permanents by first-row expansion (fine for toy sizes; Ryser's inclusion–exclusion is the serious
tool), contrasts the determinant, and then brute-forces an actual boson-sampling
distribution: two photons into a small three-mode mesh, every outcome probability a squared
permanent. Note the checks: the probabilities sum to 1, and a 50:50 beamsplitter's coincidence
permanent vanishes — the HOM dip, recomputed as linear algebra.
// The permanent: determinant's formula with all signs +1. Expansion along row 0.
function permanent(M: number[][]): number {
const n = M.length;
if (n === 1) return M[0][0];
let total = 0;
for (let j = 0; j < n; j++) {
const minor = M.slice(1).map((row) => row.filter((_, c) => c !== j));
total += M[0][j] * permanent(minor);
}
return total;
}
// The determinant by the same expansion — only the alternating sign differs.
function det(M: number[][]): number {
const n = M.length;
if (n === 1) return M[0][0];
let total = 0;
for (let j = 0; j < n; j++) {
const minor = M.slice(1).map((row) => row.filter((_, c) => c !== j));
total += (j % 2 === 0 ? 1 : -1) * M[0][j] * det(minor);
}
return total;
}
const A = [[1, 2, 3], [4, 5, 6], [7, 8, 10]];
console.log("per(A) = " + permanent(A) + " det(A) = " + det(A));
// --- A 3-mode mesh: two planar rotations, a baby Reck scheme (real, so unitary = orthogonal).
function rotation(a: number, b: number, t: number): number[][] {
const R = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];
R[a][a] = Math.cos(t); R[b][b] = Math.cos(t);
R[a][b] = -Math.sin(t); R[b][a] = Math.sin(t);
return R;
}
function matmul(X: number[][], Y: number[][]): number[][] {
return X.map((row, i) =>
row.map((_, j) => row.reduce((acc, _v, k) => acc + X[i][k] * Y[k][j], 0)));
}
const U = matmul(rotation(0, 1, 0.6), rotation(1, 2, 1.1));
// --- Boson sampling, brute force: photons enter modes 0 and 1.
// Outcome {j,k}: probability = per(2x2 submatrix)^2, halved if both photons share a mode.
let sum = 0;
for (let j = 0; j < 3; j++) {
for (let k = j; k < 3; k++) {
const sub = [[U[j][0], U[j][1]], [U[k][0], U[k][1]]];
let p = permanent(sub) ** 2;
if (j === k) p = p / 2; // bosonic normalisation for |2⟩ outcomes
sum += p;
console.log("photons out at modes " + j + "," + k + ": P = " + p.toFixed(4));
}
}
console.log("total probability: " + sum.toFixed(4));
// --- HOM as a permanent: coincidence amplitude of a 50:50 beamsplitter.
const s = Math.SQRT1_2;
console.log("HOM coincidence amplitude: " + permanent([[s, s], [s, -s]]));
Change the rotation angles and rerun: the distribution shifts, the total stays 1, and every
probability remains a squared permanent. Now imagine n = 50 photons in
m = 2500 modes: the same three-line formula, with each evaluation costing
more than the world's computers can supply.
What Aaronson–Arkhipov actually proved — and did not
Hardness of the permanent does not by itself make the experiment hard to simulate — the
machine only samples; nobody, machine included, learns any permanent's value. The
Aaronson–Arkhipov theorem forges the missing link: if a classical algorithm could sample the exact
boson-sampling distribution efficiently, then (via a classic counting argument) the polynomial
hierarchy — complexity theory's tower of generalised NP — would collapse to its third level, an
event theorists consider on a par with P = NP. For the realistic case of sampling with small errors,
they proved the same conclusion conditional on two plausible conjectures about permanents of
Gaussian random matrices. That "hardness of
quantum
advantage" argument style — sampling, not deciding; evidence, not proof — became the
template for every advantage claim since, including Google's circuit-sampling experiment in the
superconducting world and the photonic Gaussian-boson-sampling machines of the next lesson.
The honest ledger of what a boson sampler is not: it is not universal —
no feedforward, no adaptivity, no error correction, and no way to run Shor's algorithm; it
solves no useful problem anyone has found — the samples answer no question except
"what would this mesh do?"; and at scale it is not even verifiable — checking the
output distribution requires the very permanents no one can compute. It is a purpose-built
NISQ-era
artefact: a machine for winning an argument with classical computing, and for proving that
large-scale multi-photon interference survives contact with engineering reality.
Set the determinant and permanent side by side: identical sums over permutations, differing only in
that the determinant weights each term by the permutation's sign. That lone alternating sign is
what Gaussian elimination secretly exploits — row operations preserve the determinant because
cancellations are structured — and it is why \det costs
O(n^3) while \operatorname{per}, stripped of
the sign's algebraic structure, resists every such trick. Leslie Valiant made this precise in 1979,
inventing the class #P to say it: the permanent counts (perfect matchings, among other
things), and exact counting is brutally harder than deciding. Physics then adds the punchline:
bosons, whose exchange statistics carry no sign, compute permanents in their amplitudes; fermions,
whose antisymmetry supplies the alternating sign, compute determinants — which is why
non-interacting fermion linear optics is classically easy to simulate, and non-interacting
boson linear optics is a quantum-advantage machine. The entire gap between trivial and
intractable rides on the exchange phase from the
HOM
lesson.
The tempting syllogism — "outcome probabilities are permanents; the machine produces outcomes;
therefore the machine computes permanents" — is false, and seeing why inoculates you against a
whole family of quantum-hype errors. Each specific outcome has probability
|\!\operatorname{per}(U_S)|^2 \sim n!/m^n — exponentially small. To
estimate even one permanent from samples you would need to see the same outcome many
times, requiring exponentially many runs; the machine gives you a stream of mostly-unrepeated click
patterns, each individually uninformative. The hardness claim lives one level up: no classical
machine can produce a stream with the same statistics. Second trap, same vintage: "boson
sampling proved quantum computers beat classical ones." It provided strong evidence,
conditional on complexity conjectures, for a contrived task — and classical spoofing algorithms
have repeatedly clawed back specific experiments (especially lossy ones, as the next lesson's
Gaussian machines discovered). The argument is genuinely powerful; it is not a proof, and it is
not Shor.
Where this goes next
Boson sampling's weak point in practice is its fuel: it needs n
simultaneous, identical single photons, and heralded SPDC sources deliver each with probability
well under one — the n-fold coincidence rate dies exponentially. The
next
lesson swaps the fuel: squeezed light, the deterministic output of the same
nonlinear crystals, fed directly into the mesh. The permanent gives way to its pair-counting cousin
the hafnian, the photon budget explodes from tens toward hundreds, and the result — Gaussian boson
sampling — is the design behind the first photonic quantum-advantage claims, Jiuzhang and Borealis.