N single photons enter N distinct input modes of an M-mode linear-optical network built from random beamsplitters. The network implements an M×M unitary matrix U. For an output pattern with occupation numbers (t₁,…,t_M), the exact quantum probability is:
P_quantum(t) = |Per(U_ST)|² / (t1!·t2!·…·tM!)
where U_ST is the N×N submatrix formed by the N input columns and the output rows (repeated per occupation number), and Per(·) is the matrix permanent — the same formula as a determinant but with every term added instead of alternating in sign.
If the particles were classical and distinguishable instead of identical bosons, the same formula holds with the permanent taken over |U_ST|² (probabilities, no interference):
P_classical(t) = Per(|U_ST|²) / (t1!·t2!·…·tM!)
- The permanent has no efficient known classical algorithm (it is #P-hard) — this gap is exactly why sampling from P_quantum is believed to be classically intractable for large N, M, the basis of the original "quantum supremacy" boson-sampling proposal.
- Total variation distance between the quantum and classical distributions measures how much genuine multi-photon interference reshapes the output — it is the multi-photon generalization of the Hong-Ou-Mandel dip.
- Fire one shot draws one output pattern from the currently selected distribution and adds it to the tally; repeat to watch the sampled histogram converge to the theoretical curve.
Real-world relevance: this is the same computational model behind photonic quantum-advantage demonstrations (e.g. Jiuzhang, Borealis) — one of the hardware platforms national and corporate quantum strategies fund alongside superconducting and trapped-ion qubits.