Introduction

Imagine shooting individual photons — particles of light — one by one into a network of beam splitters and mirrors. Each beam splitter partly reflects and partly transmits its incoming light. When two photons meet at a beam splitter, something purely quantum happens: they refuse to behave as independent particles. Instead they bunch together, anti-bunch, or interfere in patterns governed by the laws of quantum mechanics.

At the output ports you observe which photons land where. The probability of any particular outcome is linked to a quantity called the permanent of a matrix — a sum over all permutations that looks like a determinant but without the minus signs. Computing permanents is one of the hardest tasks in all of combinatorics.

In 2011, Scott Aaronson and Alex Arkhipov proposed boson sampling as a deliberately minimal quantum experiment: it does not try to be a universal quantum computer. It just routes bosons (here, photons) through a linear-optical network and samples from the resulting distribution. The striking conjecture is that no classical computer can efficiently simulate this sampling, making even small experiments potential demonstrations of quantum advantage beyond anything classical hardware can match.

Boson sampling thus sits at a fascinating crossroads: it is experimentally accessible with today's photonic hardware, yet its output is believed to be classically intractable — a window into the power of quantum mechanics without needing full fault-tolerant qubits.

Try It: 3-Photon Network

Below is a tiny boson sampling experiment: 3 photons entering 3 input ports of a linear-optical network. The network is described by a 3×3 unitary matrix (you can randomize it). The simulator computes the exact output probability for every way 3 photons can be distributed across 3 output ports using the permanent of the relevant submatrix.

<p class="hint">A <b>3-photon boson sampling</b> experiment: photons enter ports 0, 1, 2 and exit across 3 output ports. Probabilities are computed from the <b>permanent</b> of the relevant submatrix. Randomize the network to see how drastically the distribution changes.</p>
<div class="matrix-area">
  <div class="mat-label">Network matrix U (real part shown):</div>
  <div id="matrix-display" class="matrix-grid"></div>
</div>
<div class="btns">
  <button id="randomize" type="button">Randomize Network</button>
  <button id="hadamard" type="button" class="ghost">Reset to Balanced</button>
</div>
<div class="chart-label">Output distribution (3 photons in 3 modes):</div>
<div id="chart" class="chart"></div>
<div class="status" id="status"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; font-size: 14px; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .75rem; line-height: 1.45; }
.matrix-area { margin-bottom: .6rem; }
.mat-label { font-size: .8rem; color: #666; margin-bottom: .3rem; }
.matrix-grid { display: grid; grid-template-columns: repeat(3, 54px); gap: 3px; width: max-content; }
.mat-cell { background: #e8eef3; border: 1px solid #cdd9e3; border-radius: 5px;
            text-align: center; padding: 4px 2px; font: 600 12px ui-monospace, monospace;
            color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .75rem; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .85rem;
         border: 1px solid #1d3557; background: #1d3557; color: #fff;
         border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.chart-label { font-size: .8rem; color: #666; margin-bottom: .3rem; }
.chart { display: flex; flex-direction: column; gap: 4px; }
.bar-row { display: flex; align-items: center; gap: 6px; }
.bar-key { font: 600 11px ui-monospace, monospace; color: #1d3557; width: 80px; flex-shrink: 0; text-align: right; }
.bar-track { flex: 1; background: #e8eef3; border-radius: 4px; height: 20px; overflow: hidden; }
.bar-fill { background: #1d3557; height: 100%; border-radius: 4px; transition: width .3s; }
.bar-pct { font: 600 11px ui-monospace, monospace; color: #555; width: 44px; flex-shrink: 0; }
.status { font-size: .85rem; color: #444; margin-top: .5rem; min-height: 1.2em; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Boson sampling simulation: 3 photons in 3 modes.
// Input state: one photon in each of modes 0, 1, 2.
// Output basis: all ways to place 3 photons across 3 modes (multisets of size 3).
// For output pattern s = [s0, s1, s2] (sum = 3), the probability is
//   P(s) = |Perm(U_S)|^2 / (s0! * s1! * s2!)
// where U_S is the 3x3 submatrix with row i repeated s_i times.
// For this demo all inputs have one photon each, so the input factor is 1.

// Compute permanent of a square matrix (Ryser's formula for small n).
// Perm(A) = (-1)^n * sum_{S subset [n], S nonempty} (-1)^|S| * prod_i (sum_{j in S} A[i][j])
function permanent(M) {
  const n = M.length;
  let total = 0;
  // iterate over all non-empty subsets of columns
  for (let mask = 1; mask < (1 << n); mask++) {
    const bits = [];
    for (let i = 0; i < n; i++) if (mask & (1 << i)) bits.push(i);
    const k = bits.length;
    // row products: for each row, sum of M[row][col] for col in bits
    let prod = 1;
    for (let r = 0; r < n; r++) {
      let rowSum = 0;
      for (const c of bits) rowSum += M[r][c];
      prod *= rowSum;
    }
    // inclusion-exclusion sign: (-1)^|subset|
    const sign = (k % 2 === 0) ? 1 : -1;
    total += sign * prod;
  }
  // outer factor: (-1)^n
  return ((n % 2 === 0) ? 1 : -1) * total;
}

function factorial(n) { let f = 1; for (let i = 2; i <= n; i++) f *= i; return f; }

// Generate all multisets of size n from {0..modes-1}
function outputPatterns(modes, n) {
  const results = [];
  function gen(start, remaining, current) {
    if (remaining === 0) { results.push([...current]); return; }
    for (let i = start; i < modes; i++) {
      current.push(i); gen(i, remaining - 1, current); current.pop();
    }
  }
  gen(0, n, []);
  return results;
}

// Build output occupation vector from pattern list
function patternToOccupation(pattern, modes) {
  const occ = new Array(modes).fill(0);
  for (const p of pattern) occ[p]++;
  return occ;
}

// Build the submatrix for input modes inputOcc and output occupation outputOcc
// For input all-single-photon: select row i repeated inputOcc[i]=1 times -> just U itself
// For output occupation outputOcc: repeat column j outputOcc[j] times
function buildSubmatrix(U, inputOcc, outputOcc) {
  const rows = [];
  for (let i = 0; i < U.length; i++) {
    for (let k = 0; k < inputOcc[i]; k++) rows.push(i);
  }
  const cols = [];
  for (let j = 0; j < U.length; j++) {
    for (let k = 0; k < outputOcc[j]; k++) cols.push(j);
  }
  return rows.map(r => cols.map(c => U[r][c]));
}

// Current network matrix (3x3 real unitary from Hadamard-like construction)
let U = [];

function balancedMatrix() {
  // 3x3 DFT-like real unitary (from Hadamard construction)
  const s = 1 / Math.sqrt(3);
  return [
    [s,  s,  s],
    [s, -s * Math.cos(2 * Math.PI / 3) - s * Math.sin(2 * Math.PI / 3),
         s * Math.sin(2 * Math.PI / 3) - s * Math.cos(2 * Math.PI / 3)],
    [s,  s * Math.cos(4 * Math.PI / 3) - s * Math.sin(4 * Math.PI / 3),
        -s * Math.sin(4 * Math.PI / 3) - s * Math.cos(4 * Math.PI / 3)]
  ];
}

function randomUnitary3() {
  // Random 3x3 real orthogonal matrix via Gram-Schmidt
  function randVec() { return [Math.random()-0.5, Math.random()-0.5, Math.random()-0.5]; }
  function dot(a, b) { return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; }
  function sub(a, b, s) { return [a[0]-s*b[0], a[1]-s*b[1], a[2]-s*b[2]]; }
  function norm(a) { return Math.sqrt(dot(a, a)); }
  function normalize(a) { const n = norm(a); return [a[0]/n, a[1]/n, a[2]/n]; }
  let v0, v1, v2;
  // Generate with retry to avoid near-zero norms
  do { v0 = normalize(randVec()); } while (norm(randVec()) < 0.01);
  let w1 = randVec(); w1 = sub(w1, v0, dot(w1, v0)); v1 = normalize(w1);
  let w2 = randVec(); w2 = sub(w2, v0, dot(w2, v0)); w2 = sub(w2, v1, dot(w2, v1)); v2 = normalize(w2);
  return [v0, v1, v2]; // rows
}

function computeDistribution(matrix) {
  const modes = 3;
  const n = 3; // photons
  const inputOcc = [1, 1, 1]; // one photon per input mode
  const patterns = outputPatterns(modes, n);
  const probs = [];
  let total = 0;
  for (const pattern of patterns) {
    const occ = patternToOccupation(pattern, modes);
    const sub = buildSubmatrix(matrix, inputOcc, occ);
    const perm = permanent(sub);
    const denom = occ.reduce((acc, k) => acc * factorial(k), 1);
    const prob = (perm * perm) / denom;
    probs.push({ occ, prob });
    total += prob;
  }
  // Normalize (should already be ~1 for unitary, but floating point)
  for (const p of probs) p.prob /= total;
  return probs;
}

function renderMatrix(matrix) {
  const grid = document.getElementById('matrix-display');
  grid.innerHTML = '';
  for (let r = 0; r < 3; r++) {
    for (let c = 0; c < 3; c++) {
      const cell = document.createElement('div');
      cell.className = 'mat-cell';
      cell.textContent = matrix[r][c].toFixed(3);
      grid.appendChild(cell);
    }
  }
}

function renderChart(dist) {
  const chart = document.getElementById('chart');
  chart.innerHTML = '';
  // Sort by probability descending
  const sorted = [...dist].sort((a, b) => b.prob - a.prob);
  const maxProb = sorted[0].prob;
  for (const { occ, prob } of sorted) {
    const label = '(' + occ.join(',') + ')';
    const pct = (prob * 100).toFixed(1);
    const barWidth = maxProb > 0 ? (prob / maxProb * 100) : 0;
    const row = document.createElement('div');
    row.className = 'bar-row';
    row.innerHTML =
      '<div class="bar-key">' + label + '</div>' +
      '<div class="bar-track"><div class="bar-fill" style="width:' + barWidth + '%"></div></div>' +
      '<div class="bar-pct">' + pct + '%</div>';
    chart.appendChild(row);
  }
}

function update() {
  renderMatrix(U);
  const dist = computeDistribution(U);
  renderChart(dist);
  const permSum = dist.reduce((s, d) => s + d.prob, 0);
  document.getElementById('status').textContent =
    fmt('Total probability: %v%% (10 output patterns, computed: %length%)', { v: (permSum * 100).toFixed(1), length: dist.length });
  // Log top outputs for node verification
  const top = [...dist].sort((a, b) => b.prob - a.prob).slice(0, 3);
  console.log('Top outputs:', top.map(d => '(' + d.occ.join(',') + ')=' + (d.prob*100).toFixed(2) + '%').join(' '));
}

document.getElementById('randomize').addEventListener('click', () => {
  U = randomUnitary3();
  update();
});
document.getElementById('hadamard').addEventListener('click', () => {
  U = balancedMatrix();
  update();
});

U = balancedMatrix();
update();
Result
Console

Notice how changing the matrix changes the output distribution dramatically. The permanent calculation here is trivial for 3 photons — but for 50 photons it would require summing over 50! ≈ 3 × 106410^{64} permutations. That exponential blowup is exactly why classical computers cannot efficiently simulate large boson sampling experiments.

The Real Complexity

Why is boson sampling believed to be hard for classical computers? The key is the matrix permanent.

  • What is a permanent? For an n×n matrix M, its permanent is Σ (over all permutations σ) of the product M[i][σ(i)]. Unlike the determinant, no minus signs appear, so no Gaussian-elimination trick works. Computing the permanent exactly is #P-complete (Valiant, 1979) — as hard as counting the solutions to any NP problem.
  • Why does it appear? When n indistinguishable photons enter n input ports and you ask for the probability that each output port receives exactly one photon, the amplitude is |Perm(U_S)|² / n! where U_S is the submatrix of the network's unitary matrix U selected by the input and output modes.
  • The hardness argument (Aaronson–Arkhipov, 2011). If there were a classical algorithm that could efficiently sample from the boson sampling distribution (even approximately), then — under two widely-believed conjectures about the permanent — the polynomial hierarchy would collapse to the third level. That is considered extremely unlikely by complexity theorists, so the existence of such an algorithm is conjectured to be false.
  • Status: open conjecture. The classical hardness of boson sampling is not a proven theorem — it rests on the conjectures that (1) estimating permanents of random Gaussian matrices is #P-hard even approximately, and (2) the output distribution of a random boson sampler cannot be spoofed efficiently. Disproving either conjecture would remove the complexity argument, but no one has succeeded in doing so.

This places boson sampling in a unique spot: unlike P vs NP, the hardness claim is about sampling rather than decision problems, and it is directly testable with quantum optics hardware. Related problems like counting #P underlie the same mathematical challenge.

Where It Matters

Boson sampling started as a theoretical curiosity but has grown into a rich research programme:

  • Quantum supremacy benchmarks: experiments by groups at USTC (China, 2020–2022) and elsewhere have run boson sampling variants with 50–100 photons, claiming outputs that classical supercomputers cannot reproduce in practical time.
  • Molecular vibronic spectra: the same permanent formula governs the Franck–Condon factors that describe how molecules absorb light and vibrate. A boson sampler can in principle simulate vibronic spectra of molecules — a direct application to quantum chemistry.
  • Graph problems via the hafnian: a variant called Gaussian boson sampling computes the hafnian (a relative of the permanent for symmetric matrices), which counts perfect matchings in graphs. This links photonic hardware to combinatorial optimization.
  • Randomness certification: the unpredictability of boson sampling outputs can serve as a source of certified quantum randomness, useful in cryptography.
  • Foundational tests: boson sampling experiments probe Hong–Ou–Mandel interference (the photon bunching effect) at scale, testing whether quantum mechanics continues to hold in complex multi-photon regimes.

The field is moving fast: each new experiment pushes the photon count higher and sharpens the question of whether classical simulators can ever catch up.

Conclusion

Boson sampling is a beautiful reminder that quantum mechanics is not just strange — it may be computationally powerful in ways classical physics cannot match. A handful of photons routed through a passive optical network can, in principle, produce a distribution whose simulation requires resources that scale exponentially on any classical machine.

The experiment does not need qubits, error correction, or a universal gate set. It needs only indistinguishable photons and well-calibrated beam splitters. And yet the output is connected to one of the hardest computational objects we know — the matrix permanent, a #P-complete quantity tied to counting, permutations, and the deep structure of counting problems.

Whether classical computers will ever find a shortcut to simulate large boson samplers remains an open question, but the burden of proof now lies with the classical side. As photonic hardware matures, boson sampling experiments will keep raising the bar — and keep asking whether quantum computers truly offer an exponential edge.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/boson-sampling/Content licensed under CC BY-NC 4.0.