Paired Quantum Trajectories: Common Random Numbers for Protocol Selection in Continuously Monitored Quantum Systems

1. Introduction

An open quantum system under continuous monitoring evolves coherently between detection events and jumps at random times. The quantum jump unraveling of the master equation makes this structure exact: the waiting time until the next jump is drawn from the decay of the no-jump norm, and the jump then selects a dissipation channel [1,2,3,4]. Radaelli, Landi and Binder made the correspondence explicit by recasting jump trajectories as a Gillespie process [5]: the same algorithm that underlies chemical kinetics [7] and, in its scheduling structure, discrete-event simulation generally. Under counting-type unravelings the detection record is not merely simulable by discrete events. It is a marked point process exactly, and not as an approximation of some underlying continuous description: it is what the detector produces. (Diffusive unravelings, homodyne and heterodyne, are outside this correspondence and outside this paper.)

This paper takes the correspondence seriously in both directions, and each direction produces results the other cannot.

In one direction, we build the laboratory inside a discrete-event engine designed for verification (Section 2). The dominant modelling errors in this domain are structural, not numerical: a decay path omitted from the operator set, or a detector inefficiency modelled as a loss channel and then also counted in the total rate. Neither produces a crash or an implausible trajectory. Both produce a smooth, converged, tightly error-barred, wrong ensemble average, and more replications never catch them. What catches them is exhaustive bookkeeping: machine invariants checked on every transition of every trajectory (the channel rates sum to the declared decay; the conditional state has unit norm at every anchor; photons reconcile exactly against clicks, inefficiency misses, and dead-time absorptions; every click is consumed by exactly one controller rule), and load-time gates on the declared wiring (every jump channel has exactly one detection destination; every emitted event kind is covered). The engine machine-checks the event structure and the bookkeeping, not the Hamiltonian: the coherent evolution between jumps is closed-form numerics that no model checker touches, and we validate it instead against analytic known answers and against QuTiP [9] (Section 3).

In the other direction, we import the discrete-event community’s standard comparison discipline into a field that does not use it. The decision problem in monitored quantum systems is protocol selection: which control pulse, which feedback policy, which measurement strategy performs best under realistic detector imperfections. Current practice simulates independent trajectory ensembles per protocol and compares their means [9], which discards the fact that both protocols face the same stochastic environment. Common random numbers (CRN) make the comparison paired [13,14,15]: the shared stochastic backbone cancels out of the difference. We are not aware of any use of this technique in the quantum trajectory literature (Section 6), and the nearest quantum relative, correlated sampling in ground-state quantum Monte Carlo [11,12], lives in a different domain. Our implementation addresses every random draw by event ordinal (the k-th jump, the k-th detector arrival), so two protocol arms consume identical uniforms up to the ordinal where their event histories diverge, and a mechanical audit verifies after the fact that the shared draws were byte-identical while reporting the arm-only tails (Section 5).

Our contributions:

1.

A verified event-structure laboratory for monitored quantum systems (Section 2, Section 3). Gillespie waiting times computed in closed form, detector dead time and inefficiency as journaled bookkeeping, click-driven feedback with stale-jump cancellation, all invariant-checked per transition; validated against exact finite-horizon renewal expectations and QuTiP (twelve assertions on three levels, all passing), with the laboratory’s own estimator defect, caught by external review, as the demonstration that the target defect class is real (Section 3.1).

2.

An exact dead-time correction for antibunched resonance fluorescence, for an effect reaching 21% (Section 4). The textbook non-paralyzable correction $r/(1+r{\tau }_d)$ assumes Poisson arrivals. Sub-Poissonian arrivals are regularized, fewer land in the dead window, and the measured rate sits up to 21% above the Poisson prediction. That dead time modifies non-Poisson counting statistics is established; what we add is the exact delayed-renewal expectation for this source, reproduced by the laboratory across ${\tau }_d=0.25$ to 8 decay times.

3.

Paired trajectories, with the audit as the instrument (Section 5). The variance of the estimated protocol gap falls by $40\times$ or more for near-identical protocols (40 to 430 across five seeds: the ratio itself is heavy-tailed at the smallest gap, and we report it as a floor), decaying to $1.4\times$ as the protocols diverge, and the CRN audit’s shared-stream ledger measures the alignment mechanism directly, closely enough that, calibrated on one protocol family, it predicts another family’s variance ratios to a median factor of $1.3$ (Section 5.5). At every seed the paired design resolves protocol differences the unpaired analysis of the same data cannot.

4.

The rate-optimal feedback pulse is not $\pi$, and not because of latency (Section 5.4). The mean wait after the pulse is $a_0+a_{1}cos\theta +a_{2}sin\theta$ exactly, for every drive, detuning and latency, so the optimum is $atan2(-a_2,-a_1)$ and equals $\pi$ only when $a_2=0$. At our drive, which is the exceptional point of the effective Hamiltonian, the zero-latency optimum is exactly $3\pi /4$ and the $\pi$ pulse gives up exactly $\sqrt{2}-1=41.4\%$ of the rate; the measurement at the shortest representable latency returns $0.75\pi$ and 40% ($\pi$ is recovered only as $\Omega \to 0$, which is the statement $a_2/a_1\to 0$); latency lowers the optimum along a measured curve ($0.68\pi$ at $L=0.5/\gamma$, $0.63\pi$ at $L=1/\gamma$; measured baselines are the prediction rounded to the grid of the sweep that produced them); and the optimum does not maximize post-pulse excitation. Locating one optimum took $n=400$ paired trajectories per arm; the unpaired equivalent would need roughly an order of magnitude more.

2. Model and Simulation Framework

2.1. Physics and the Gillespie Step

The system is a driven two-level emitter under photon counting, in the rotating frame:

$$H=\Delta {\sigma }^{+}{\sigma }^{-}+\frac{\Omega }{2}({\sigma }^{+}+{\sigma }^{-}),L_k=\sqrt{{\gamma }_k}{\sigma }^{-}(k=1\dots K),$$

(1)

with total decay $\gamma ={\sum }_k{\gamma }_k$ split across K detection ports. Every result in this paper is on resonance, $\Delta =0$. Between events the conditional (unnormalized) state evolves under $H_{\mathrm{eff}}=H-\frac{i\gamma }{2}{\sigma }^{+}{\sigma }^{-}$ as $|\tilde{\psi }\left(\tau \right)\rangle =e^{-i{H}_{\mathrm{eff}}\tau }|{\psi }_0\rangle$, computed in closed form by $2\times 2$ complex eigendecomposition. No integrator and no tolerance knob appear anywhere. The no-jump survival is $S\left(\tau \right)=\langle \tilde{\psi }|\tilde{\psi }\rangle$, monotone from 1 toward 0 under driving.

The Gillespie step [5,7]: at every anchor (the start of time, each jump, each control pulse) the emitter receives two uniforms $(u_t,u_c)$ and decides deterministically. The waiting time solves $S\left({\tau }^{*}\right)=u_t$ by bracket and bisection on the closed form. The channel follows from $u_c$ against the weights ${\gamma }_k/\gamma$ ($u_c$ is drawn even when $K=1$: the fixed-draw rule below). The jump is then scheduled as a future event. After a ${\sigma }^{-}$ jump the state is $|g\rangle$ exactly. A control pulse arriving before the scheduled jump propagates the conditional state to the pulse time, applies its unitary, cancels the pending jump, and re-draws: the old waiting time was conditioned on the pre-pulse state, and reusing it would be a conditioning error that no ensemble average would flag. Time is discretized at ${10}^{-3}$ decay times per tick and waiting times round to at least one tick. The suite bounds the resulting bias by halving the tick and re-running the known-answer gates; it sits below the statistical floor of every measurement in this paper.

2.2. Machine Architecture, Invariants, and Randomness Allocation

The laboratory is composed from three machine types. The emitter holds the conditional state and the declared channel partition. The photodetector (one per port) applies efficiency $\eta$ and non-paralyzable dead time ${\tau }_d$: a photon arriving during the dead window is a journaled refusal with a named reason, counted in absorbed, never silently dropped; a live arrival is resolved against one uniform, becoming a click (which extends the dead window and notifies the controller) or a counted miss. The controller consumes clicks under a declared policy: none, or pulse-on-click with a declared angle, latency, and refractory window. A click in the refractory window is a journaled refusal. The topology declares which dissipation channel is monitored by which detector and which click drives which controller action, and the engine’s load gates check the wiring is total.

The division of responsibilities defines the architecture: responders draw, machines decide, the kernel schedules. Machines own no randomness. A responder hands the emitter its two uniforms at the current jump ordinal, from a named stream (jump:{node}:n{k}); a second responder hands each detector one uniform per arrival (eff:{node}:n{k}). Two rules make this a comparison instrument rather than a convenience. The fixed-draw rule: every ordinal consumes the same number of uniforms in every scenario, unconditionally. The ordinal-addressing rule: streams are named by event index, not by simulation time, so two protocols under one seed share their stochastic backbone draw for draw until their event histories diverge. Section 5 builds on both.

The bookkeeping invariants run on every transition of every trajectory: the channel rates sum to the declared $\gamma$ exactly; the anchor state has unit norm; photons = clicks + misses + absorptions + in-flight resolutions, exactly; clicks seen = consumed + refused; pulses fired never exceed clicks consumed. Across every run reported in this paper these invariants evaluated with zero violations, and each run reproduces its journal digest byte for byte from its content-addressed identity (world hash, scenario hash, seed, replication index).

3. Verification Against Known Answers and QuTiP

The two-level emitter is chosen because its ground truth is independently known. The steady state gives the photon rate $R=\gamma {\langle {\sigma }^{+}{\sigma }^{-}\rangle }_{\mathrm{ss}}$ with ${\langle {\sigma }^{+}{\sigma }^{-}\rangle }_{\mathrm{ss}}=({\Omega }^2/4)/({\Delta }^2+{\gamma }^2/4+{\Omega }^2/2)$, and renewal theory makes the mean inter-jump waiting time exactly $1/R$. At a finite horizon the correct comparison target is not R itself: for an ordinary renewal process started at a renewal point, $\mathbb{E}\left[N\left(T\right)\right]/T=R+({\sigma }^2-{\mu }^2)/\left(2{\mu }^{2}T\right)+o(1/T)$, and for sub-Poissonian intervals ($\sigma <\mu$) that constant is negative. We therefore gate against the exact finite-horizon renewal expectation $\mathbb{E}\left[N\right(T\left)\right]/T$, computed deterministically from the no-jump survival by a grid-and-FFT renewal sum (make_figs.py). Table 1 shows the laboratory against both targets at $n=40$ trajectories of $T=200$ decay times per cell (± values here and throughout are $2$SE, ≈95% CIs).

3.1. Case Study: An Estimator Bias Identified by External Review

The first draft of this paper carried a rate estimator that normalized clicks by the time of the last recorded event rather than by the declared horizon. For a renewal process that is the classic inspection-paradox estimator: $\mathbb{E}[N/t_N]=(1/\mu )(1+\mathbb{E}\left[A\right]/T)+o(1/T)$, with $\mathbb{E}\left[A\right]=\mathbb{E}\left[{\tau }^2\right]/2\mu$ the mean backward recurrence time, a deterministic $+0.7$ to $+1\%$ bias in every absolute rate. Every gate passed anyway: at $n=40$ the bias sat inside $2$SE, and three same-sign one-sigma excesses were absorbed by three separate error bars. This is exactly the defect class the introduction names (a smooth, converged, tightly error-barred, slightly wrong number), and it was caught not by more replications but by external review computing the backward-recurrence correction from the machine’s own survival function and reproducing the excess to four decimal places. The fix is one line (divide by the horizon); the residual finite-horizon physics is the renewal constant above, now part of the gates; and the paired-comparison results of Section 5 moved by less than a tenth of their interval widths, because the bias was common mode across arms.

We retain the episode because of what it delimits. The defect did not live in a machine: it lived in a probe, and probes are pure folds over the journal that sit outside the invariant perimeter entirely. The machine invariants cannot see an estimator, so no amount of per-transition checking was ever going to flag this, and the cross-validation did not either, because QuTiP was refereeing the same biased quantity. The constructive lesson is that estimator identities are themselves checkable, and cheaply: the horizon-normalized and last-event-normalized rates differ by precisely the bias in question, so their discrepancy is a standing gate, and the laboratory now carries it. The episode supports a narrow claim: because the event structure, the draw sites, and the estimator were all declared artifacts rather than implementation details of a bespoke simulator, an external reader could localize the defect analytically from the machine’s own survival function, without rerunning anything. Reproducibility of the event structure is what made the error findable by someone other than its author.

Figure 1 shows the full waiting-time distribution at $\Omega =2$: the antibunching hole at the origin ($W\left(\tau \right)\to {\tau }^2$: a jump resets the emitter to $|g\rangle$ and re-emission requires re-excitation [3,18]), the Rabi lobe structure under the decay envelope, and the same-mean exponential reference cutting through both. None of this structure is put in by hand. It emerges from the event process, which is the point of validating against it.

QuTiP is the external referee [9] (worlds/monitored/qutip_cross.py), on three levels mirroring the cross-validation discipline of our companion study [24]. The twelve assertions counted here are the QuTiP comparisons only; the laboratory’s own test suite is separate and is described under Reproducibility. Analytic: QuTiP’s steadystate matches the machine’s closed form to ${10}^{-9}$ at all three drives. Distributional: collapse times recorded by mcsolve give waiting-time means and short-wait fractions consistent with the laboratory’s (for example 2.240 versus 2.252 at $\Omega =2$). Statistical: the master-equation click rates sit inside the laboratory’s CIs at every drive. Twelve assertions across the three levels, all passing. Detector dead time and click-triggered feedback are outside QuTiP’s machinery, which is precisely the gap the laboratory exists to fill; their referees are the exact renewal expectations of Section 4 and Section 5.4.

One boundary should be stated where the reader will look for it. The renewal-expectation curves in Figure 5 and Figure 2 are computed in an independent numpy implementation of the closed forms (a separate codebase from the machine’s JavaScript), so agreement checks both implementations of the model. It does not check the model against nature; that is what the analytic steady state, QuTiP, and the antibunching structure of Figure 1 are for. And all parameter sweeps in this paper run a single detection port ($K=1$); the multi-channel partition machinery is exercised by the load gates, the unit suite, and one reported two-port measurement (k2-check.mjs: ports at $0.7/0.3$ of $\gamma$, per-port click rates $0.3167\pm 0.0121$ and $0.1271\pm 0.0069$ against expected shares $0.3106$ and $0.1331$ of the exact finite-horizon expectation, photons reconciling exactly per port in all 40 runs, zero invariant findings; and the paired latency comparison of Section 5.3 rerun on the two-port world resolves at a $7.9\times$ variance ratio with the audit verifying both detectors’ streams).

4. Antibunched Light and the Dead-Time Formula

Every photon-counting experiment corrects for detector dead time. The standard non-paralyzable correction [19] is

$$r_{\mathrm{meas}}={\displaystyle \frac{r}{1+r{\tau }_d}}.$$

(2)

That formula assumes Poisson arrivals. That dead time modifies the counting statistics of non-Poisson light is an established literature: dead-time-corrected photocounting distributions go back to Cantor and Teich [20], dead-time-modified count moments and their observation to Vannucci and Teich [21], and corrections for correlation measurements through modern single-photon detectors [22,23]. Clark et al. model count-rate saturation for bunched, Poissonian, and antibunched sources through detectors with a measured efficiency-recovery profile [23]; our step-function dead window is the non-paralyzable idealization of that recovery, and the machine’s journaled design would absorb a graded profile as a time-dependent $\eta$. The light in this laboratory is antibunched: arrivals are regularized and short gaps are suppressed (Figure 1), so fewer photons arrive inside the dead window, and the loss is smaller than Eq. (2) predicts. The contribution here is not the existence of the effect but its exact, gate-grade computation for resonance fluorescence, verified in a laboratory where dead time is a first-class journaled event.

The exact prediction is renewal-theoretic. At $\eta =1$ the jumps form a renewal process with waiting density $W\left(\tau \right)=-S^{\prime }\left(\tau \right)$ from $|g\rangle$. An accepted click at a renewal point opens a dead window ${\tau }_d$; the next accepted click is the first jump after it, so the inter-accepted intervals are i.i.d. with $I={\tau }_d+Y\left({\tau }_d\right)$, Y the forward recurrence time at lag ${\tau }_d$, while the first interval of a trajectory is an ordinary wait (the first click pays no dead time). The accepted-click process is therefore a delayed renewal process, and we compute its exact finite-horizon expectation $\mathbb{E}\left[N_{\mathrm{acc}}\left(T\right)\right]/T$ by the same grid-and-FFT renewal sum as the gates of Section 3. Figure 2 shows the laboratory against this prediction and against Eq. (2) across ${\tau }_d=0.25$ to 8 decay times: every measured point sits within $0.71$ interval half-widths of the exact curve, and above the Poisson formula everywhere, by up to 21% at ${\tau }_d\approx 1/\gamma$, converging only in the long-dead-time limit where both scale as $1/{\tau }_d$.

The practical implication is direct. Dead-time corrections calibrated on Poisson sources over-correct antibunched sources. The size of the effect at ${\tau }_d\sim 1/\gamma$ is tens of percent, and the correction requires the waiting-time distribution, not just the rate. The laboratory computes it because dead time is a first-class, journaled event in the detector machine, not a post-hoc analysis step.

5. Paired Trajectories

5.1. Paired-Comparison Method and Audit Procedure

Comparing feedback protocols with independent ensembles estimates ${\overline{Y}}_A-{\overline{Y}}_B$ with variance $\mathrm{Var}\left(Y_A\right)/n+\mathrm{Var}\left(Y_B\right)/n$. Pairing replication i of protocol A with replication i of protocol B on the same draws subtracts $2\mathrm{Cov}(Y_A,Y_B)/n$ [13,15]. The technique is decades old in operations research. It is not guaranteed to be a gain: the classical theory gives monotonicity conditions under which the covariance is positive [13], and a recent planning result exhibits decision processes where full pairing increases the variance [17]. Nothing forces $\mathrm{Cov}(Y_A,Y_B)>0$ for quantum protocol arms either. The variance ratio is therefore a measured output of every comparison in this paper. Every ratio we measure is above one. The quantum trajectory literature evaluates protocols on independent realizations (Section 6), and the reason is partly mechanical: pairing requires that “the same randomness” be well defined across two systems whose event counts differ, and a global RNG sequence loses alignment at the first divergent event.

Ordinal stream addressing solves the alignment problem. The k-th jump decision of the emitter draws from a stream named by k in both arms; the k-th arrival at a detector likewise. Two protocols under one seed therefore consume identical uniforms at every ordinal both of them reach. When the histories diverge (one protocol jumps more often), the arms consume different total draw counts, and that is not an error: it is a measured property of the comparison. After every paired run the engine’s CRN audit compares the two arms’ per-stream draw ledgers and refuses the paired analysis on any mismatch, while reporting the arm-only tails.

Byte-identity of shared streams in a correct engine is close to tautological, since streams derive deterministically from seed and name. The audit’s content is threefold: it catches fixed-draw-rule violations, an arm whose draw count per ordinal became state- or scenario-dependent, which is exactly the silent CRN-killer in hand-rolled trajectory pairings; it catches ordinal misaddressing and stream-name collisions; and its shared/arm-only ledger is the alignment instrument of Section 5.3. A mutation test in the laboratory’s suite injects a responder that consumes one extra uniform past a scenario-controlled jump ordinal (both arms still run; both produce plausible physics), and the audit refuses the comparison with the mismatched streams named (CRN_AUDIT_FAILED, monitored.test.mjs).

Protocol arms are scenario documents (data, not code), so a protocol comparison is a declared experiment: paired scenarios, a declared endpoint (the accepted-click rate), paired-t intervals, and the audit, in one report.

5.2. Distant Protocols: A Negative Result

Pairing is not automatically beneficial, and the first measurement establishes its scope. Comparing feedback off against pulse-on-click at $\Omega =0.5\gamma$ ($n=60$), the paired and unpaired interval widths nearly coincide: variance ratios 1.1 to 1.4 across latencies. The reason is physical. The protocols differ so strongly (the pulse arm emits three times faster) that their event histories decorrelate after the first click, the outcome correlation collapses, and the covariance term buys almost nothing. The audit ledger reflects this directly: 81 shared streams against up to 380 arm-only streams. For protocol screening across qualitatively different policies, independent ensembles are essentially optimal. This is the scope boundary of the method.

5.3. Close Protocols: The Variance Curve

The selection decisions that matter are between close variants: this latency against that one, this pulse angle against that one. There the trajectories stay aligned and pairing pays exactly where the decision is hard. Figure 3 sweeps the latency gap $\Delta L$ between two pulse-on-click protocols at $n=400$: the variance ratio falls monotonically from $40\times$ at $\Delta L=0.01$ decay times to $1.4\times$ at $\Delta L=1.5$, and the audit’s shared-stream count falls along the same curve (319 shared streams per pair at the smallest gap, 171 at the largest). The arm-only tails are instructively asymmetric (Figure 3b): the baseline-only count grows (150 to 254) while the variant-only count collapses toward zero (90 to 3), because the slower high-latency arm reaches fewer event ordinals and its consumed set becomes a near-subset of the baseline’s. Protocol distance shows up not as symmetric divergence but as one arm’s ordinal range outrunning the other’s, and the ledger distinguishes the two. For close protocols, correlation, alignment, and the audit ledger are three views of the same phenomenon, and Section 5.5 turns the third into a quantitative predictor of the first.

The operational cell is $\Delta L=0.05$: the paired estimate of the protocol gap is $-0.0099\pm 0.0017$, resolved decisively, while the unpaired analysis of the same 800 trajectories reports $\pm 0.0077$ and resolves it only marginally (the point estimate exceeds the interval by a factor of 1.3); at $\Delta L=0.01$ the paired design resolves $-0.0024\pm 0.0012$ (a 0.7% relative difference) where the unpaired interval $\pm 0.0078$ contains zero with room to spare. In terms of computational cost, an unpaired study would need at least $40\times$ the trajectories to make the $\Delta L=0.01$ decision the paired study makes at $n=400$. One caution on breadth: the sweep runs eight gaps at $\alpha =0.05$ with no multiplicity correction; the resolved/unresolved contrasts quoted here are the two pre-declared headline cells, not the product of multiple-comparison selection.

The $40\times$ should be read as a floor, not a point value, and the distinction is itself a measured property of the method (seed-rep.mjs). Replicating both headline cells across five seeds leaves the paired point estimate essentially unchanged ($-0.0020$ to $-0.0024$ at $\Delta L=0.01$) and the decision pattern stable within each cell: at $\Delta L=0.01$ the paired design resolves at every seed and the unpaired design at none, while at $\Delta L=0.05$ both resolve at every seed. The marginal unpaired resolution reported above is a property of that gap, not of this paper’s seed. The variance ratio at the smallest gap ranges from 40 to $433\times$, with this paper’s seed at the minimum. The mechanism is the same alignment physics read at one more level: at $\Delta L=0.01$ most replication pairs never invert their event order and contribute nearly zero difference, so the paired variance is carried by the few pairs that do, and how many such pairs a seed contains fluctuates. The difference distribution is heavy-tailed exactly where pairing pays most, which is also a caution about paired-t intervals in that regime; at the operational gap $\Delta L=0.05$ the ratio is stable, 20 to $44\times$ across the same five seeds.

5.4. The Pulse-Angle Optimum Is Not $\pi$, Even at Zero Latency

The applied payoff is protocol optimization. The intuitive click-triggered protocol is a $\pi$ pulse: the emitter just jumped to $|g\rangle$, so re-excite it fully. The intuition fails in two separable ways, and neither is a latency artifact.

The first failure needs no latency at all, and it is available in closed form. Write the post-pulse state as $|\psi \rangle =R_x\left(\theta \right)|{\psi }_L\rangle$, with $R_x\left(\theta \right)=cos(\theta /2)\mathbb{1}-isin(\theta /2){\sigma }_x$ and $|{\psi }_L\rangle$ the normalized conditional state at the end of the latency window. The no-jump survival $S\left(\tau \right)=\langle \psi \left|M\right(\tau \left)\right|\psi \rangle$, with $M\left(\tau \right)=e^{i{H}_{\mathrm{eff}}^{†}\tau }{e}^{-i{H}_{\mathrm{eff}}\tau }$, is a quadratic form in $|\psi \rangle$, so the half-angle products collapse onto $\{1,cos\theta ,sin\theta \}$ and the mean wait is

$$\mathbb{E}\left[\tau \right]\left(\theta \right)=a_0+a_{1}cos\theta +a_{2}sin\theta ,$$

(3)

with $a_0=\frac{1}{2}\int (\langle M\rangle +\langle {\sigma }_{x}M{\sigma }_x\rangle )$, $a_1=\frac{1}{2}\int (\langle M\rangle -\langle {\sigma }_{x}M{\sigma }_x\rangle )$ and $a_2=\int \mathrm{Im}\langle M{\sigma }_x\rangle$, all evaluated in $|{\psi }_L\rangle$ and independent of $\theta$. Equation (3) is exact for every drive, detuning and latency. The mean cycle is $L+\mathbb{E}\left[\tau \right]$, so the rate optimum and its value are

$${\theta }^{*}=atan2(-a_2,-a_1),\mathbb{E}{\left[\tau \right]}^{*}=a_0-\sqrt{a_1^2+a_2^2},$$

(4)

and ${\theta }^{*}=\pi$ if and only if $a_2=0$. The coefficient $a_2$ carries the off-diagonal response to ${\sigma }_x$, and it vanishes only when the pre-pulse state leaves the drive no dipole to act on. That is the first failure in one line. The $\pi$ pulse is rate-optimal exactly when the drive does no work during the wait, and the weak-driving recovery of $\pi$ is the statement $a_2/a_1\to 0$: ${\theta }^{*}(L=0)$ runs from $0.94\pi$ at $\Omega =0.1\gamma$ to $0.65\pi$ at $\Omega =\gamma$ (Figure 4b). Coherence with the drive, not population, sets the re-emission clock.

At the operating drive the coefficients are elementary, because $\Omega =\gamma /2$ is the exceptional point of $H_{\mathrm{eff}}$: the eigenvalues ${\lambda }_{\pm }=\left[-i\gamma /2\pm {({\Omega }^2-{\gamma }^2/4)}^{1/2}\right]/2$ coalesce, the propagator acquires a nilpotent term, and the survival is $(1+b\tau +c{\tau }^2)e^{-\gamma \tau /2}$ rather than a sum of two exponentials. At $L=0$ this gives $(a_0,a_1,a_2)=(4,2,-2)/\gamma$, so

$$\mathbb{E}\left[\tau \right]\left(\theta \right)={\displaystyle \frac{4+2cos\theta -2sin\theta }{\gamma }},{\theta }^{*}={\displaystyle \frac{3\pi }{4}},\mathbb{E}{\left[\tau \right]}^{*}={\displaystyle \frac{4-2\sqrt{2}}{\gamma }},$$

(5)

and the $\pi$ pulse gives up exactly $\sqrt{2}-1=41.4\%$ of the rate. The zero-latency optimum is $3\pi /4$, not a value near it.

Two statements bound what this covers. Equation (3) is the asymptotic cycle; the finite-horizon optimum from the renewal sum of Section 3 at $T=200$ sits $+0.002\pi$ above it uniformly in L, inside the $0.01\pi$ grid of the sweeps. And because Figure 4b carries $\Omega$ from $0.1\gamma$ to $\gamma$ across the exceptional point with no feature in ${\theta }^{*}(\Omega )$, the structure of the optimum is not an artifact of operating there.

The engine represents a latency of one tick (${10}^{-3}/\gamma$, the quantization floor), so the zero-latency limit is measurable and not merely derived. At one-tick latency (theta-zero.mjs, $n=400$/arm) the $0.75\pi$ protocol is resolved above both declared neighbors ($0.70\pi$ at $-0.0240\pm 0.0030$, $0.80\pi$ at $-0.0223\pm 0.0031$) and above the $\pi$ pulse by $-0.3349\pm 0.0069$, a 40% rate deficit for the textbook choice against the derived 41.4%.

The second failure is the latency pre-rotation. During L the conditional state evolves under the drive and acquires part of the rotation, $|{\psi }_L\rangle$ acquires a dipole, $a_2/a_1$ grows, and ${\theta }^{*}$ falls along a monotone curve, from $0.75\pi$ at $L\to 0$ to $0.58\pi$ at $L=1.5/\gamma$ (Figure 4a and Table 2). The measurement tracks the curve (theta-L.mjs, $n=400$/arm): at $L=0.1$, $0.25$, and $1.0/\gamma$ the predicted-optimum baseline beats both its $\pm 0.05\pi$ neighbors by the paired design (all six comparisons resolved below, the tightest at $-0.0071\pm 0.0010$), just as it does at the operating point of Figure 5.

At the operating point ($\Omega =0.5\gamma$, $L=0.5/\gamma$) the click process under the protocol is a delayed renewal process (a long un-pulsed first wait from $|g\rangle$, then cycles of length $L+{\tau }_{\mathrm{pulsed}}$), and its exact finite-horizon expectation, computed as in Section 3, has its maximum at ${\theta }^{*}\approx 0.68\pi$. The post-pulse excitation peaks elsewhere, near $0.93\pi$ (Figure 5b): the rate optimum deliberately sacrifices population for coherence with the drive.

Figure 5a shows the measurement: with the predicted optimum as the paired baseline, every measured angle is resolved below it, including the two pre-declared nearest neighbors, $0.65\pi$ at $-0.0030\pm 0.0017$ (a 0.5% relative difference) and $0.75\pi$ at $-0.0201\pm 0.0020$ (the other five angles are also resolved, but the localization claim rests on the two declared neighbors, so the seven tests are not the product of multiple-comparison selection). The whole optimization cost $n=400$ trajectories per arm; the variance ratios against the baseline run 3 to 17 across the sweep, with the two nearest neighbors at 16 and 17: an unpaired-equivalent cost of roughly 7,000 trajectories per arm for the decisions that localize the optimum.

The absolute agreement with the theory curve is weaker than the localization result, and it separates into two parts. Every measured point at $L=0.5/\gamma$ sits within one interval half-width of the exact expectation ($+0.21$ to $+0.95$ half-widths). All eight residuals are positive, but the eight arms are paired against a common baseline, so they carry one common-mode realization between them rather than eight independent ones. Reading the paired differences separates the two contributions.

A common-mode offset of $+0.0013$ to $+0.0018$ is present at every angle, and the during-latency emission branch accounts for it: a jump inside the latency window (probability $1-S\left(L\right)=0.22\%$ per cycle) yields an extra counted click and delivers the pulse to a freshly reset state, contributing $+0.0011$ to $+0.0013$ to the rate across the sweep. Net of that offset the paired differences match the theory to within $1\times {10}^{-4}$ from $0.55\pi$ to $0.75\pi$, which is the range the localization claim occupies.

A separate effect appears at $\theta \ge 0.8\pi$, reaching $+0.0046$ in the paired difference at $0.9\pi$. It does not follow the angle. It follows post-pulse excitation, peaking where the excitation peaks ($0.998$ at $0.9\pi$) and falling back at $\pi$, where the excitation falls to $0.988$. Three candidate sources are excluded. Quadrature error: Richardson extrapolation at $h=0.005$ and at $h=0.0025$ agrees to five decimals at every angle, and the residual grid drift falls with $\theta$ while the discrepancy grows. Solver tolerance: the survival inversion runs bracket doubling then bisection to machine precision. Tick quantization: the probability of a wait shorter than half a tick is $5\times {10}^{-4}$ even at full inversion, an order of magnitude too small and of the wrong sign. The effect is not identified. The one-tick-latency arms show an offset of the opposite sign (measured 2.5% below the prediction).

These bound what the theory comparison establishes: the two implementations of the renewal model agree with each other to quadrature accuracy, and the machine reproduces the model to $1\times {10}^{-4}$ in paired difference across the localization range and to about 1 to 2% in absolute rate at the largest angles. The localization claims are unaffected, because they rest on paired differences within a single arm-pair, where a common-mode model offset cancels; this is the same cancellation the paper’s method exploits throughout, here visible as a property of the comparison rather than of the model.

5.5. How Much of the Achievable Precision Does Pairing Recover?

Section 5.3 and Section 5.4 report variance ratios against the unpaired design, which answers what pairing buys relative to current practice but not what it leaves on the table. The click record is a renewal process with a known interval density, so the question has an absolute answer: the Fisher information of that record bounds any estimator of the protocol gap, and the achieved variance can be placed against it.

For the pulse-on-click protocol the inter-click interval is $I=L+\tau$, with $\tau$ drawn from the waiting density of the pulsed state $\psi (\theta ,L)$. The interval density $f_I(x;L)=W(x-L;\psi (\theta ,L))$ depends on the latency twice over, by translation and through the pre-pulse evolution, so the information is not the trivial shift information. At the operating point ($\Omega =0.5\gamma$, $L=0.5/\gamma$, $\theta =\pi$) numerical quadrature on the closed forms of Section 3 gives ${\mathcal{I}}_1=401$ per interval and, at $rT\approx 66$ intervals per trajectory, $\mathcal{I}=2.65\times {10}^4$ per trajectory. By the delta method the Cramér-Rao floor on the standard error of the rate gap at $n=400$ is $|{\partial }_{L}r|/\sqrt{n\mathcal{I}}=6.5\times {10}^{-5}$.

Table 3 places four estimators of the same decision against that floor. The rows between the two extremes are computed in an idealized renewal replica in which both arms draw the same uniform at every ordinal, which is the strongest coupling ordinal addressing can express. The replica is not a prediction of the laboratory: it reproduces the laboratory’s baseline rate to $0.2\%$, its paired gap to $2\%$ and its unpaired standard error to $3\%$, but its arms never lose alignment, so its paired row is a ceiling rather than a measurement.

For the two coupled rows the third column is a ratio and not an efficiency. A coupled design is a stronger experiment than the independent-arms floor describes, so the paired score at 68% of that floor does not license the reading that the estimator is close to optimal: the bound for a coupled design sits lower, and where it sits is the open question at the end of this section.

Two separate losses account for the distance between current practice and the bound, and the table separates them.

The first is the endpoint. The count is a lossy summary of a renewal record: it discards the interval structure that carries most of the information about the latency. Replacing the count by the score ${\sum }_k{\partial }_{L}log{f}_I(I_k;L)$, at the same coupling and the same n, moves the estimator from 36% to 68% of the bound, a further factor of $3.5$ in variance. Pairing and the endpoint are therefore independent axes: pairing removes the common noise, the endpoint decides how much of the remaining signal is read.

The second is alignment, and the laboratory already measures it. Under ideal coupling the paired count reaches a variance ratio of 451; the laboratory achieves 40 at the same gap, and the shortfall is trajectory divergence. That invites reading the audit ledger of Section 5.3 quantitatively rather than as illustration. Across all seventeen paired comparisons in this paper, taking the overlap relative to the larger arm, $\omega =s/max(s+b,s+v)$, the variance ratio is exponential in the overlap (Figure 6):

$$\mathrm{VR}\approx exp\left[13.96(\omega -0.407)\right].$$

(6)

The fit is taken over the fourteen close-protocol comparisons, but the substantive test is out of sample. The pulse-angle family occupies $\omega \in [0.504,0.571]$, inside the latency family’s $[0.402,0.680]$, so the relation can be calibrated on the latency family alone and used to predict the angle family cold. Fitted on the eight latency comparisons, Eq. (6) predicts the six angle comparisons to a median factor of $1.29$ and a worst case of $1.58$, across a perturbation type the fit never saw: a latency change translates the interval density, a pulse-angle change reshapes it. The ledger is therefore not merely a correctness check on the pairing. It is a predictor of what the pairing will buy, and it carries across the kind of perturbation within a calibrated range of overlap.

That range qualification is load-bearing. Fitted separately the two families have slopes $13.5\pm 1.2$ (latency, $n=8$) and $21.8\pm 3.4$ (angle, $n=6$), differing by $2.3$ standard errors, so they do not lie on a common line; predicting in the other direction, out of the angle family’s narrow $\omega$ span, degrades to a median factor of $1.84$. Within each family the comparisons also share a baseline arm, so the fourteen points carry fewer independent degrees of freedom than their count suggests. What the measurement supports is an interpolating predictor over a calibrated range of overlap, not a universal exponent.

The exception is the distant-protocol family of Section 5.2, which sits off the line and whose measured ratios fall below the prediction by about half at comparable overlap. That is the regime in which the arms differ threefold in rate, and it is the regime the same section already recommends running unpaired, so the predictor fails exactly where its answer would not be used.

Two consequences follow. Because the dependence is exponential, overlap is a sensitive design lever: the difference between $\omega =0.5$ and $\omega =0.6$ is a factor of four in the achievable precision, which argues for choosing protocol parameterizations that keep arms on shared ordinals. And because the relation is calibrated rather than derived, a single pilot pair fixes $\omega$ and predicts the benefit of a planned comparison before the full study is committed. Deriving Eq. (6) from the coupling structure, rather than fitting it, remains open.

Three scope statements. The information is the classical Fisher information of the jump record under this unraveling, not the quantum Fisher information over all measurement strategies; the latter is the subject of Radaelli et al. [6], whose renewal and Markov-chain machinery computes it properly, and whose bound would sit below ours. The quadrature uses the ordinary-renewal approximation, so the delayed first interval and horizon censoring enter at order $1/66$. And the bound is stated for independent arms; establishing the corresponding bound for a coupled design, where the shared randomness is part of the experiment rather than a nuisance, is the open theoretical question this section raises.

6. Related Work

Quantum jump unravelings go back to Dalibard, Castin and Mølmer [1,2] and Carmichael [3], with Plenio and Knight’s review [4] and Wiseman and Milburn’s treatment of measurement and feedback [8] as standard references; QuTiP’s mcsolve is the field’s reference implementation [9]. Radaelli, Landi and Binder recast the jump unraveling as a Gillespie process [5]; we build directly on that observation, executing the Gillespie step inside a general discrete-event engine rather than a bespoke simulator, which is what makes the event structure checkable and the comparison discipline available. The same group has since carried the formalism from simulation to inference: the Fisher information of jump records, computed through their renewal and Markov-chain structure, with a Gillespie-based sampler [6]. That work asks how much one protocol’s click record reveals about a parameter; ours asks which of two protocols is better. The questions are complementary and share the renewal machinery, and their Fisher information bounds the precision any estimator of a rate difference, paired or not, can reach.

Variance reduction exists in trajectory simulation on an orthogonal axis: choosing the unraveling to reduce single-ensemble estimator variance [10]. That technique and ours compose (one picks the representation, the other correlates the arms) and answer different questions. The nearest quantum use of correlated randomness is correlated sampling in ground-state quantum Monte Carlo, where synchronized walkers estimate energy differences between nearby trial functions or geometries [11,12]. The mechanism is the same; the domain (imaginary-time projection for electronic structure) and the failure mode of alignment are not. In classical simulation, CRN theory and practice are mature [13,14,15], including recent empirical analyses of when seed pairing pays [16] and, in simulation-based planning, a depth-dependent rollout estimator that shares randomness only where the compared policies’ histories align [17]: an independent arrival at the same structural insight as our ordinal addressing, that the value of pairing lives in the shared prefix before the compared systems diverge. Our contribution on this axis is the transplant into quantum trajectories, the ordinal-stream construction that keeps arms aligned, and the mechanical audit that verifies the pairing instead of assuming it. We are not aware of a prior CRN-paired protocol comparison for quantum trajectory ensembles; QuTiP documents seeds for reproducibility of a single calculation only.

7. Discussion and Limitations

The methodological claim travels beyond this laboratory. Any trajectory simulator that addresses its randomness by event ordinal (which the Gillespie formulation makes natural [5]) can run paired protocol comparisons and, with a draw ledger, audit them. The measured scoping rule is simple: pairing pays in proportion to trajectory alignment, so it is the tool for the fine end of protocol selection (latencies, pulse angles, thresholds), and independent ensembles remain appropriate for coarse screening. The audit matters as much as the pairing: a comparison whose arms silently consumed different shared draws is wrong in a way no error bar reveals, and the refusal semantics (shared streams must be identical; arm-only tails are reported, not refused) are, to our knowledge, the correct and previously unstated contract for state-dependent event processes.

Limitations. The physics here is a two-level emitter with radiative channels only: no cavity mode, no dephasing channel, no thermal excitation, a single detection port in every parameter sweep (one two-port check is reported in Section 3), and jump unraveling only (homodyne and heterodyne unravelings are diffusive and need a different event discipline). The endpoint is the click rate; state-dependent endpoints (fidelity, concurrence) need the same treatment but were not run. The dead-time result is for non-paralyzable detectors at unit efficiency, gated against the exact delayed-renewal expectation. No warm-up deletion is used, and none is owed: at $\eta =1$ without feedback every inter-click interval starts at a jump anchor in $|g\rangle$, so the click process is exactly renewal from the start of time, and under feedback every cycle likewise starts at a renewal point; the finite-horizon effects that remain are computed, not deleted. The variance curve is measured for one protocol family on one observable; the decay of alignment with protocol distance should be generic, but its scale is not. Equation (6) is calibrated rather than derived, its two close-protocol families have slopes that differ by $2.3$ standard errors, and it is tested only by interpolation within the overlap range the latency family spans. The information bound of Section 5.5 is the classical Fisher information of the jump record under this unraveling rather than the quantum Fisher information over all measurement strategies, is computed in the ordinary-renewal approximation, and is stated for independent arms; the ideal-coupling rows of Table 3 come from a renewal replica whose arms never lose alignment, and are a ceiling rather than a measurement. The coherent evolution between events is trusted numerics validated against known answers and QuTiP, not verified logic: the engine machine-checks the event structure and the bookkeeping, and that boundary is stated wherever it binds.