Quantum computing research,
built on 23 years of engineering.

The KLM gates (Knill, Laflamme, Milburn, Nature 409, 46, 2001). The nonlinear sign gate, from KLM's own 3-mode unitary and from the beamsplitter form (angles 22.5°, 65.5302°, −22.5° and a π phase): on 50 random inputs α₀|0⟩ + α₁|1⟩ + α₂|2⟩ it gives α₀|0⟩ + α₁|1⟩ − α₂|2⟩ with probability 1/4 and fidelity 1 − 10⁻¹⁵. Two of them make a CZ on dual-rail qubits (8 modes, 4 photons): process fidelity 1 − 3 × 10⁻¹⁶, success exactly 1/16.

Checks that bite, and one that did not. The middle angle off by 0.5° lowers the CZ fidelity to 0.9997. Leaving out the π phase was expected to break the gate; it does not — it gives diag(1, −1, −1) per NS gate, and the circuit becomes CZ followed by Z on both qubits, still a CZ up to single-qubit phases. So that is no power check, and we say so.

Fusion (Browne and Rudolph, PRL 95, 010501, 2005). Type-I fusion, simulated photon by photon — polarising beamsplitter, 45° rotation, one polarisation-resolving detector — acts on the two input qubits as (|0⟩⟨00| ± |1⟩⟨11|)/√2 (to 10⁻¹⁶), succeeds with the input's weight on |HH⟩ and |VV⟩ (½ for cluster qubits), and its failures are Z measurements of both. Type II succeeds with probability ½ and projects onto one Bell state, Φ⁺ or Ψ⁺: linear optics tells apart two of the four. On state vectors, fusing the ends of linear clusters of n and m qubits (2–5) gives the cluster of n + m − 1 with fidelity 1 − 2 × 10⁻¹⁶.

Growing a chain. A long linear cluster grows by fusing k-qubit pieces onto its end, each piece itself built from Bell pairs by fusion; on a failure the chain loses its end qubit and the piece its first. Bell pairs used per qubit gained (Monte Carlo, 20 × 20,000 attempts):

Bell pairs alone (k = 2) never grow the chain — success adds one qubit, failure takes one away. Three-qubit pieces cost 8 pairs per qubit, or 7 when the two-qubit leftover of a failure is reused: the review's "2 × 4 − 1 = 7". Its "6.5" we do not find in any strategy we simulated; the best are 7.0, with three- or four-qubit pieces and reuse (data, fusion checks, /quantum/optics/fusion).

Measurement-based computation (on state vectors). One-bit teleportation through a two-node cluster and a rotation along a chain (adaptive angles, Pauli byproducts) reach fidelity 1 on random inputs. The review's four-node example does not hold as drawn: its Fig. 10 chain (CZ₁₂CZ₂₃CZ₃₄) gives a CZ between the teleported input and the fourth qubit — CNOT-class by its local invariants, but needing gates before as well as after; the couplings its Eq. 56 writes (1–2, 2–3, 2–4) do give a CNOT for every outcome and for superposition inputs, with Pauli byproducts when built from CZ gates (data).

Detectors and loss. A detector of efficiency η is a beamsplitter onto an empty mode before an ideal detector: its outcome probabilities match the bucket and number-resolving formulas (the review's Eqs. 80–81) to 4 × 10⁻¹⁶. Through such detectors on the heralds, lost photons make false heralds; with bucket detectors, which cannot tell one photon from two, the gates fail even at η = 1:

And a photon lost from a dual-rail qubit leaves exactly ηρ + (1 − η)|vac⟩⟨vac| (to 2 × 10⁻¹⁶): an erasure, heralded by the missing photon — the error model of our next step, decoding erasures (data, CSV).

Erasure decoding — the part that is ours. A lost photon is a heralded erasure: the qubit is known to hold a random Pauli. Our union-find decoder, QxDecode, takes erased edges as already grown and weightless in the peeling tree — for pure erasure the peeling decoder of Delfosse and Zémor, which is maximum-likelihood. At code capacity (the rotated surface code, perfect syndrome, 20,000 shots a point, PyMatching with the erased edges at weight 0 on the same shots beside it), pure erasure crosses at 50 % for d = 3 to 9, the percolation threshold, and the two decoders agree:

With bit flips on the other qubits the crossing falls as the erased share grows — near 8–10 % of flips without erasure, about 8 % with 10 % erasure, about 6 % with 20 % — in line with the usual linear trade-off; QxDecode is a little less accurate than matching without erasure (as union-find is) and close to it once erasures dominate:

At circuit level, QxStab samples heralded erasure (Stim's HERALDED_ERASE semantics; on a d = 3 memory with 2 % erasure its 51 detectors and 27 heralds agree with Stim's sampler, exact tests, smallest corrected p 0.89). Surface-code memories, d rounds, Pauli noise 0.1 % on every operation and heralded erasure on every data qubit each round; the decoding graph from Stim's error model, each herald mapped to the edges its erasure can flip; 20,000 shots a point:

Knowing where the losses were is worth about twenty times in logical error at d = 5 and 4 % erasure, and keeps larger codes better even at 8 % per round — where, without the heralds, the larger code is already worse. QxDecode uses the heralds about as well as matching does (data, pure-erasure data, mixed data, tools/erasure/).

The 2.7 % loss figure. Romero and Milburn's review (p. 27) says fusion-based computation "can tolerate a 2.7 % probability of suffering a single-photon loss". It comes from Bartolucci et al. (Nat. Commun. 14, 912, 2023, Sec. VI and Appendix B): the 6-ring fusion network with (2,2)-Shor-encoded resource states and fusions boosted with a Bell pair (failure 1/4, four photons per fusion). We reproduce it from their formulas — 2.71 % per photon, which is 10.4 % that at least one of the four photons is lost, and their 43.2 % failure threshold at zero loss — taking their simulated erasure threshold (11.98 %) as given; we have not simulated the fusion network itself. Unencoded, the same network tolerates no loss at failure 1/4; the 4-star network, encoded, about 1.0 % (tools/optics/fbqc_loss.py).

One mode: GKP and cat qubits (Gottesman, Kitaev, Preskill, PRA 64, 012310, 2001; Ralph et al., PRA 68, 042319, 2003). A third simulator, for one optical mode in continuous variables, again written twice: cv.py in Python works in the Fock basis and on a grid in q, QxCV in Elixir on the grid alone, its p side summed analytically over the peaks instead of by FFT. Squeezed vacuum from the Fock basis matches the closed-form wavefunction to 8 × 10⁻¹⁵, and Eq. 24 holds to 9 × 10⁻¹⁴. Eq. 26 holds only if the phases of Eqs. 22 and 25 agree: with ζ = re^(−2iϕ) as printed, the rotation must be e^(+iϕ), not Eq. 25's e^(−iϕ) (flipping either sign works); as printed, the variances are off by up to 5.6 when ϕ ≠ 0.

GKP states (Eq. 38; α = √π, peaks squeezed 10 dB, Δ² = 0.01). The logical Z̄ = e^(iπq/α) reads ±0.924 on |0̃⟩ and |1̃⟩, equal to an analytic sum over the Gaussian peaks to 10⁻¹⁰; X̄ = e^(−ipα) takes |0̃⟩ to |1̃⟩ with fidelity 0.984 — the envelope does not move with the peaks — also as predicted. Read out after a Gaussian shift of σ, the logical error from the exact integral, from 200,000 sampled readouts, and from the closed form for an ideal comb with σ² plus the peak's variance agree (the closed form to 1.5 × 10⁻⁷), in q and in p. Eq. 39's bounds show as an error of ½ (to 2 × 10⁻⁵) at a shift of α/2 in q or π/2α in p. QxCV's X readout (the Fourier sum) and Python's (the FFT) agree to 10⁻¹³:

The review's Fig. 5, read from the vector curve in its PDF. The peaks sit at multiples of √π, as Eq. 38 gives with α = √π; each peak's variance fits 0.11 — near e^(−2r) = 0.09 for r = 1.2, while the caption says e^(2r) = 11; and the |1̃⟩ envelope falls off four times faster than the |0̃⟩ one (4.00), which a single Δ in Eq. 38 cannot give. With the caption's Δ̃ = 0.01, amplitudes e^(−Δ̃(sα)²) for |0̃⟩ and e^(−Δ̃[(2s+1)α]²) for |1̃⟩ fit both panels: the |0̃⟩ panel looks drawn with s where Eq. 38 has 2s.

Cat states (Eq. 75). From the Fock basis and from the grid they agree to 10⁻¹⁵; they hold only even or only odd photon numbers, and the Wigner function at the origin is ⟨Π⟩/π = ±1/π (to 2 × 10⁻¹⁶). Eq. 40's N±⁻¹ leaves the state with norm² 1/N±; Eq. 75's N±^(−1/2) is right. Photon loss is the cat code's error — one lost photon flips the parity — and after a loss 1 − η, ⟨Π⟩ = cosh(β²(1 − 2η))/cosh β² for the even cat (sinh for the odd): the generating function (QxCV) and the Fock amplitudes thinned photon by photon (Python) agree to 2 × 10⁻¹⁶.

A comb heralded by homodyne detection (Eqs. 66–73). With Eq. 69's Kraus amplitude the outcomes y sum to one for every q (to 10⁻⁹); Eq. 70 as printed sums to 1/√2 — it lacks the ½ in the exponent. A phase-squeezed input with the outcome y = 0 gives peaks every π/θ (√π here: a GKP-like comb); other outcomes give pairs, with gaps 2 arccos(y/2β)/θ and the rest of the period, as sampled:

checks · Elixir against Python · readout · Eq. 39 shifts · Fig. 5 as read · comb · tools/optics/cv_checks.py, xcheck_cv.py · /quantum/optics/gkp, /cat, /comb

From the review we started from (Romero and Milburn, arXiv:2404.03367): its Fig. 8 gives the third angle as +22.5°, which does not make an NS gate — only the text's −22.5° reproduces KLM's unitary; we also found slips in its Eqs. 13, 22, 40 and 70, a sign between Eqs. 22 and 25, its growth arithmetic, its Fig. 5 (caption and envelopes) and its Fig. 10 (above), listed in the repository's NOTES.md. Every number here comes from the primary paper or our own derivation.

HOM data · KLM checks · against Perceval · tools/optics/ · /quantum/optics/hom?gt=1, /quantum/optics/ns, /quantum/optics/cz?input=11&shots=10000

Question. Under concurrent load, does a union-find decoder on the BEAM — every syndrome its own lightweight process, many decoded at once on all cores — keep its tail latency bounded, at a throughput near a compiled decoder's?

Baselines, to separate the algorithm from the language. Union-find in a compiled language (Rust or C++), and in plain Python; and PyMatching 2.4 (minimum-weight matching, C++) — each run under the same load as ours: a worker pool on the same cores, fed from a queue, with 1, 100 and 10,000 syndromes in flight. Same machine, same syndromes (sampled by QxStab and by Stim from the threshold-study circuits, d = 3 to 9, p = 0.1 to 1 %). The decoding graph comes from Stim's detector error model for every decoder, ours included: there our work depends on Stim.

The deadline, and why the BEAM. The BEAM is soft real-time: it will not decode for superconducting hardware, whose rounds take about a microsecond. Trapped-ion and neutral-atom machines run rounds of a millisecond or more, and the software that orchestrates any machine works in milliseconds; we test against a 1 ms deadline at d = 5. The question is whether preemptive per-process scheduling keeps the tail there when the queue grows, where a single-threaded decoder's queue waits pile up.

Hypothesis, and what falsifies it. At d = 5: our 99th-percentile latency stays under 1 ms with 10,000 syndromes in flight, growing no more than 2× from 1 to 10,000; and our throughput per core is within 2× of compiled union-find. Either number missing falsifies it. We expect to lose against PyMatching, which decodes such syndromes in microseconds, and we report by how much — and the logical error rate, since union-find is a little less accurate than matching. A partial result still says something: bounded latency at a lower throughput would describe what the BEAM is for.

Protocol. One warm-up pass; 10 repeated runs per point, medians with 95 % intervals by bootstrap; garbage collection left on and reported (BEAM processes collect their own heaps; we record the time spent); scripts, raw timings and syndromes published with the results.

QxDecode (lib/qx_decode.ex): union-find after Delfosse and Nickerson, weighted (edge lengths from the error probabilities, in quarter steps of the lightest edge), peeled over a minimum-weight spanning tree of each cluster. Its first version peeled over a breadth-first tree and failed a simple pattern — two detectors next to each other and the boundary — ten times too often at low noise; comparing with PyMatching shot by shot found it. On the same 100,000 syndromes per point, both on Stim's detector error model:

What this says. Accuracy: our logical error rate is 1.06 to 1.37 times PyMatching's where both have enough failures to tell — union-find's usual distance from matching; McNemar's exact test on the shots where they disagree finds the difference real at most points. Speed: per syndrome on one core we are 14 to 50 times slower than PyMatching; on all 20 cores our throughput is about one PyMatching core. That is the throughput half of the hypothesis already answered for PyMatching (we lose, as expected); the comparison against compiled union-find and the latency under load follow below (data, tools/dem_graph.py, tools/decode_runs.exs, tools/decoder_compare.py).

First, our design had a flaw. With a fixed number of syndromes in flight, Little's law makes the mean latency that number divided by the throughput: 10,000 in flight at 1 ms needs 10 million syndromes a second from any decoder. That measures throughput, not scheduling. We kept the closed-loop runs, and added the usual test: syndromes arriving at random (Poisson) at 50, 80 and 95 % of what each decoder sustains through its own harness, and the latency from arrival to answer.

The setup. The same 8 fast cores (of a 20-core Arm machine that also runs other services) for every decoder; the same graph and syndromes at p = 0.4 %. QxDecode on 8 BEAM schedulers, one process per syndrome, released by a dispatcher that spins on the clock; union-find in C++ (tools/load/uf.cpp, the same algorithm: identical predictions on 99.94–100 % of shots), a spinning dispatcher and 7 worker threads; PyMatching and union-find in Python, a dispatcher and 7 worker processes. Each harness also ran with a decoder that does nothing, to show its own cost. 10 seeds at d = 5, 3 at d = 3 and 7; medians.

Against the hypothesis, at d = 5. Throughput within 2× of compiled union-find: no — 12 times less decoding alone, 20 times less through the harness. 99th percentile under 1 ms: at half load 1.24 ms (95 % interval 894 µs–1.73 ms), at 80 % 13.57 ms, at 95 % 8.53 ms. Growth bounded: no. The hypothesis is falsified.

What it does show. A process per syndrome costs the BEAM more than the decoding: its harness alone, with no decoding, has a 99th percentile of 2.41 ms at half load for d = 3 — spawning and messaging at a few hundred thousand a second is the bottleneck, not the scheduler's fairness. C++'s tails are bimodal across seeds — microseconds in some runs, milliseconds in others — which fits interference from the machine's other work on shared cores; a clean tail measurement needs a machine to itself. Python's harness costs more than any deadline here (its no-op p99 is tens of milliseconds), so for PyMatching and Python union-find only throughput is meaningful: about 124k and 93k syndromes a second through it at d = 5.

Next. Decode many syndromes per process (batches) rather than one each, to separate the BEAM's scheduling from its spawning cost; a dedicated machine for the tails; then the question again. Data: summary, every run; tools/load/.