Paper deep dive
A Single Atom in Front of a Mirror is a Universal Reservoir Computer
Peter J. Ehlers, Phi Hung Nguyen, Kanu Sinha, Noelle Daigle, Travis W. Sawyer, Hendra I. Nurdin, Daniel Soh
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/13/2026, 5:20:59 AM
Summary
This paper demonstrates that a single atom in front of a mirror constitutes a universal reservoir computer. The authors prove that this minimal quantum setup can approximate any fading-memory map with arbitrary accuracy by adjusting measurement settings and mode numbers, without requiring hardware retraining. The system leverages the atom's saturable nonlinearity and the mirror-induced delay to create a rich state space, outperforming or matching classical baselines on tasks like Mackey-Glass and NARMA10.
Entities (9)
Relation Signals (8)
Single Atom in Front of a Mirror → implements → Reservoir Computing
confidence 95% · We show that universality can be associated with a single reservoir, considering a minimal setup of a single atom in front of a mirror.
Single Atom in Front of a Mirror → benchmarkson → Mackey-Glass
confidence 92% · on Mackey–Glass and NARMA10, as node count grows
Single Atom in Front of a Mirror → benchmarkson → NARMA10
confidence 92% · on Mackey–Glass and NARMA10, as node count grows
Single Atom in Front of a Mirror → approximates → Fading-Memory Maps
confidence 90% · our reservoir is a universal approximator of fading-memory maps under an operating class of checkable conditions
Standing-Wave Modes → constitute → Network State Space
confidence 90% · the standing-wave modes are the network
Saturable Response → provides → Activation Function
confidence 88% · the atom’s saturable response the activation function
Atom-Mirror Distance (L) → determines → Round-Trip Delay
confidence 85% · the round-trip delay τ = 2L/v the recurrence
Single Atom in Front of a Mirror → operatesin → Linear-Transducer Limit
confidence 85% · In its linear-transducer limit, our reservoir is a universal approximator
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Universal approximation in reservoir computing is typically associated with a class of reservoirs. We show that universality can be associated with a single reservoir, considering a minimal setup of a single atom in front of a mirror. In its linear-transducer limit, our reservoir is a universal approximator of fading-memory maps under an operating class of checkable conditions, with a rate constant measured at the operating point. A given reservoir can reach arbitrary accuracy by changing measurement settings. The proof gives an explicit recipe: for a target accuracy, it specifies the required physical resources and resonator modes. Enlarging the number of accessible modes increases the matchable kernel span without reducing capability. Beyond the linear limit, the atom's saturation replaces high-order polynomial readouts, and the device operates on real-world tasks alongside classical baselines. Our results highlight an example of universality with a minimal quantum setup.
Tags
Links
- Source: https://arxiv.org/abs/2608.10382v1
- Canonical: https://arxiv.org/abs/2608.10382v1
Trouble viewing inline? Open PDF directly →
Full Text
282,249 characters extracted from source content.
Expand or collapse full text
A Single Atom in Front of a Mirror is a Universal Reservoir Computer Peter J. Ehlers 1,*,† , Phi Hung Nguyen 1,*,† , Kanu Sinha 1,2 , Noelle Daigle 1 , Travis W. Sawyer 1 , Hendra I. Nurdin 3 , and Daniel Soh 1,* 1 Wyant College of Optical Sciences, University of Arizona, Tucson, AZ, USA 2 Department of Physics, University of Arizona, Tucson, AZ, USA 3 School of Electrical Engineering and Telecommunications, University of New South Wales, Sydney, Australia * Corresponding authors: ehlersp@arizona.edu, hpnguye7@arizona.edu, danielsoh@optics.arizona.edu † These authors contributed equally. Abstract Universal approximation in reservoir computing is typically associated with a class of reservoirs. We show that universality can be associated with a single reservoir, considering a minimal setup of a single atom in front of a mirror. In its linear-transducer limit, our reservoir is a universal approximator of fading-memory maps under an operating class of checkable conditions, with a rate constant measured at the operating point. A given reservoir can reach arbitrary accuracy by changing measurement settings. The proof gives an explicit recipe: for a target accuracy, it specifies the required physical resources and resonator modes. Enlarging the number of accessible modes increases the matchable kernel span without reducing capability. Beyond the linear limit, the atom’s saturation replaces high-order polynomial readouts, and the device operates on real-world tasks alongside classical baselines. Our results highlight an example of universality with a minimal quantum setup. 1 Introduction How much hardware does a universal learning machine require? In reservoir computing—a fixed dynamical sys- tem replacing trained recurrent connections, only a linear readout trained [1, 2]—the answer has always been: more than one machine. Universality has been proven for echo-state networks [3, 4], state-affine systems [5], stochastic echo-state net- works [6], and Ising, spin-ensemble, and Gaussian quan- tum reservoirs [7–10]—always associated with a class, by ranging over weights, couplings, or realizations, with no single member carrying the property. The two known single-member constructions [11, 12] obtain it by making the reservoir digital, so it is no longer physical hardware at all. Physically realized reservoirs [13] scale expressivity by adding components; the counter-trend, the single-node delay-line reservoir [14, 15], has been hardware reservoir computing’s workhorse for over a decade with no univer- sality guarantee ever established for it (Supplementary Ta- ble S.5). Supplementary Table S.5 states the basis of that claim from a literature search across the delay-reservoir theory line [16–19]: capacity analyses, class-level equiva- lences and folded-in-time emulations, never a universality theorem attached to one fixed physical member. In this work we show the minimal architecture can carry the guarantee: a single atom before a mirror is a provably universal reservoir computer. “Atom” means any two-level emitter coupled to the guided field: trapped atom or ion, quantum dot, colour centre, or superconduct- ing artificial atom. The device (Fig. 1a) has one geomet- ric parameter entering the feedback structure, the atom– mirror distance L, alongside the fixed band envelope and mode spacing that define the device, and one scaling re- source, the number of field modes it accesses, set by geom- etry and bandwidth. We prove universality of this single machine in a precisely stated limit with explicit rates; ca- pability monotone, and the matchable kernel span strictly growing, with mode number; and genericity of the condi- tions, holding at the exact parameters of every simulation here. It is Supplementary Table S.5’s only entry com- bining physical hardware, a single-member guarantee and strict scaling. The theorem is proven for the architecture with its node operated as a linear element and the non- linearity carried by the readout, and is available precisely because the atom–mirror loop then decomposes into inde- pendent linear modes—the structure the classical device’s nonlinear node destroys. The result is possible because quantum optics sup- plies, in one passive object, every ingredient reservoir computing otherwise assembles from parts (Box 1): the standing-wave modes are the network, the atom’s sat- urable response the activation function, the round-trip delay τ = 2L/v the recurrence [20–22]—a geometry re- alized with trapped atoms [23] and superconducting ar- tificial atoms [24]. Nothing inside is trained, wired, or manufactured as a network. We claim no quantum computational advantage: quan- tum mechanics enables minimality, collapsing the recur- rent, high-dimensional, nonlinear core of a learning ma- chine into the passive dynamics of one atom, with only a linear readout trained classically. Correspondingly the empirical standard for the demonstrations is parity with mature classical methods, which is exactly what “minimal hardware suffices” predicts. “Non-Markovian” refers to the reduced atomic dynam- ics: the atom interacts with its own past through return- ing photons, while the joint system is Markovian to the axioms’ standard (Methods; Supplement Sec. S.1). 1 arXiv:2608.10382v1 [quant-ph] 11 Aug 2026 Box 1 | The device in plain terms. A reservoir computer needs three things: a large state space to hold information, a nonlinearity to mix it, and recurrence to remember it. This device gets all three from geometry. Memory: light emitted by the atom travels to the mirror and comes back after a delay τ = 2L/v; what the atom does now therefore depends on what it did one round trip ago, and, through repeated reflections, on progressively fainter echoes of its earlier history—a memory that fades at a rate set by how fast light leaks out rather than by any clock or register. State space: between atom and mirror the field forms standing waves, one for each frequency the waveguide supports; each behaves as an independent oscillator, and together they play the role of the hidden neurons of a recurrent network—added by widening the usable bandwidth rather than by fabricating components. Nonlinearity: a two-level atom saturates (it cannot absorb a second photon while holding one) which is the activation function, built into the atom by quantum mechanics. Training touches none of this: one measures the outgoing light and fits a linear map from measurement records to targets, the only learned object in the machine. 2 Results 2.1 Mode-space representation of the non-Markovian reservoir The description on which every claim rests: the atom– mirror reservoir, with all its delayed feedback, is captured by a single Lindblad generator on the atom plus its field modes, the feedback carried exactly—relative to axioms A1–A6 of Supplement Sec. S.1 (rotating wave; flat cou- pling g(ω) ≈ p γ/2π over bandwidth B; unidirectional outcoupling; classical drive; idealized measurement). The mirror boundary selects standing-wave modes; passing to this basis is a unitary transformation of the free field—no auxiliary or lossy degrees of freedom, no pseudomodes—and the coupling acquires the form factor ̃g(ω)∝ sin(ωτ/2) sampling the standing wave at the atom. All dependence on L, hence all delay and feedback, resides in this coherent structure, and the retarded self-interaction is recovered identically from its frequency integral (Meth- ods). The reservoir is described by ̇ρ(t) = Γ g · aρ(t),a † + a,ρ(t)a † · Γ † g (1) − iηu(t)[(a † · α) + (α ∗ · a),ρ(t)], in which a is the vector of mode annihilation opera- tors, u(t) the scalar input drive with coupling ampli- tude η, and Γ g = iΩ + γ 2 (v geo ⊗ v ∗ geo ) + γ g 2 (α ⊗ α ∗ ), with Ω the diagonal matrix of standing-wave frequen- cies, α (α k ∝ sin(ω k τ/2)) the normalized emission pro- file, v geo the geometry-fixed, flat profile of the readout leakage channel (v geo,k = 1/ √ K), and γ, γ g the read- out and drive-channel rates. The anti-Hermitian part car- ries the coherent feedback; the Hermitian rank-one terms carry the Markovian measurement- and drive-induced evo- lution; the homodyne local oscillator selects only which functional of the record is computed (full reading in Sup- plement Sec. S.1). A discrete Ω is a comb of spacing ∆ 0 = πv/ℓ with quantization length ℓ an operating pa- rameter; no continuum limit enters the proof. 2.2 Universality in the linear-transducer limit The output is obtained from powers of the measured quadrature, ˆy k = N X n=0 Z T off 0 W n (t) Tr[q n geo e L 0 t ρ ζ (t k )] dt,(2) with q geo = Re[v ∗ geo ·a] the measured quadrature of the ge- ometric leakage channel, W n (t) classical readout weights, and L 0 the generator at u = 0; the drive alternates input windows T on and measurement windows T off , with read- out period T = T on + T off . The powers come from the same homodyne record but cost (2n− 1)!! (v 0 +q 2 max ) n /σ 2 shots at precision σ: the factorial price of Gaussian-limit operation, and exactly what the nonlinearity transfer of Sec. 2.5 eliminates (Methods; Supplement Sec. S.7). The proof rests on the adiabatic elimination of the atom—a reduction of the atom rather than of the environment—which locks the saturable dipole into fol- lowing the field and replaces it by a linear trans- ducer, with controlled small parameter the saturation s ∼ (ε INPUT /γ g ) 2 (Methods). We call this the linear- transducer (Gaussian) limit. Three kinds of resource. The theorem quantifies over one fixed machine: fabricated hardware never changes with the target; geometric and measurement settings are se- lected after the accuracy is announced, fabricating noth- ing. Measurement settings alone reach every accuracy down to an explicit envelope of the current geometry; the geometric dial, provably never wasted, moves the envelope itself (full taxonomy: Supplement Sec. S.2.5). A network whose width is a runtime setting still stores per-width cou- pling data; every mode-number description of this device follows by one rule from one set of physical data (Supple- ment Sec. S.1.1). Theorem 1 (Constructive universality with explicit rates). Consider the class C of reservoir computers with dynamics Eq. (1) and outputs Eq. (2) in the linear- transducer limit, realized by a single atom–mirror device operated at the explicitly constructed point of Supplement Sec. S.2.5: • Overlap by construction: τ avoids the explicit node set, so α k ̸= 0, and the geometry-fixed coupler is flat, v geo,k = 1/ √ K; the overlap conditions v ∗ geo ·v R,k ̸= 0, v ∗ L,k ·α̸= 0 (eigendata of Γ g ) then hold, with explicit margins (1 − ε v )/(2 √ K) on the coupler side and 1/(2q √ K) on the drive side certified by the comb- spacing condition of the next bullet, where ε v is the coupler-flatness error of Supplement Sec. S.1 and enters every lower bound below as a subtraction; • Weak-dressing comb: an equidistant comb with spac- ing ∆ 0 ≥ 8q √ K ̄e, where q is the prime of the delay lock of the third bullet, ̄e≥∥E∥ is the dressing con- stant of Supplement Sec. S.2.5 and E the Hermitian dressing of Γ g , large enough that second-order dress- ing shifts fall below the node tolerance, and locked to the readout period (a timing condition on T alone, 2 condition D2 of Supplement Sec. S.2.5) so the nodes e −λ k T sit at the M-th roots of unity of common ra- dius ̄r = e −λ ⋆ T , λ ⋆ = (γ + γ g )/2K; • Graded defects and generic drive spread: one exponentially graded defect hierarchy (Supplement Sec. S.2.5, D5) separates all conjugation-resolved eigenvalue-sum points of order n ≤ N with distinct signed mode content, and the drive channel’s decay- rate spread, floored by a delay lock placing the atom at a rational point of the standing wave, separates the remainder with an explicit margin, so every sum point is distinct or exactly conjugate and the higher- order selection is well posed; • Weak drive, flat coupler: conditions D4 and D6 of Supplement Sec. S.2.5 cap the drive-channel and coupler-flatness spreads within the node-placement budget; D6 alone is checked rather than set, the cou- pler flatness being fabricated data. Conditions D1–D4 are geometric or measurement settings of the one fixed device, while condition D5 specifies the resonator spectrum, likewise a geometric setting of the same fixed device, as Supplementary Table S.2 records, and each is proven satisfiable and sufficient with explicit con- stants. Then for every continuous fading-memory target y k = y(u k ,u k−1 ,... ) with a convergent Volterra expansion on K u max = [−u max ,u max ] Z − and every ε > 0, there exist finite N, memory depth M, mode number K = M, read- out period T = O(M lnM + lnε −1 )/(γ + γ g ), and weights W n (t) such that the single device achieves |y k − ˆy k |≤ ε on K u max , with the explicit bound |y k − ˆy k | ≤ δ MN + A(M,N ) 2M ̄r M + eM ̄ε , ̄r M = e −(γ+γ g )T/2 . Here δ MN is the target’s Boyd–Chua truncation error, A(M,N ) = P N n=1 n 2 2n u n max M n−1 ∥h (M,N ) n ∥ ∞ is set by the target’s truncated kernels and input bound alone, ̄ε is the node-placement error driven below any tolerance by the comb spacing, and ̄r M recovers the full trace budget over the memory span. No weight norm, revival time, or loop-stability constant enters; no infinite-mode limit is in- voked; parameters are chosen once and in order, (M,N ) then T then the comb, with no circular dependence (Sup- plement Sec. S.2.5). The rate is claimed for Volterra- convergent targets; merely continuous fading-memory tar- gets are reached through the density route of Supplement Sec. S.3.2, without a rate. The certificate’s price belongs beside the theorem. The designed operating point demands a comb spacing growing exponentially with memory depth, a relative timing pre- cision falling exponentially with it, a drive weak enough to inflate shot budgets by orders of magnitude, and a rate ordering γ g ≪ λ ⋆ (condition D4) outside the adia- batic derivation of the working generator from the physi- cal atom, so that at this point the theorem is a statement about the Gaussian generator of Eq. (1) (Methods; Supple- ment Sec. S.2.5). Its reach, bounded in Methods, does not extend beyond M = 2–3 on any platform. Every simula- tion in this paper instead operates at generic points whose distinctness conditions are verified rather than designed, a relation the Supplement makes precise: the certificate proves existence with explicit constants, and generic op- eration is how the device is used. The full proof is in Supplement Sec. S.2, its logic in Methods. Two features shape it. A uniform spectral gap is provably impossible—the Hermitian part of Γ g has rank two, so the slowest mode closes as 1/K—and the proof turns this into the rate’s engine, the fixed budget (γ + γ g )/2 over K = M modes being recovered in full over the memory span. And the bound never touches the matched weights: once kernels are matched at lags ≤ M, deeper response is fixed by the spectrum alone (extrapo- lation identity, Sec. S.2.5). Delay stability is the condi- tion under which the device has fading memory (Secs. S.4, S.10), never an input to the bound. Each condition is physical: v ∗ L,k ·α̸= 0 says no standing wave has a node at the atom, and v ∗ geo ·v R,k ̸= 0 is imme- diate for the flat coupler. Delay stability D(c 0 )—roots of s + γ 2 (1 + re iφ e −sτ ) = 0 with Re[s] ≤ −c 0 , the quantum echo-state condition, failing only at the dark state—holds with wide margins (c 0 = 0.62γ, 0.36γ at γτ = 0.5, 1) and is not consumed by the bound. A structural consequence is monotonic scaling. Corollary 1 (Strict scaling). For M ≤ M ′ and N ≤ N ′ , every accuracy achievable by the (M,N )-member of C on every target is achievable by the (M ′ ,N ′ )-member (Supplement Sec. S.3.1, capability containment): enlarg- ing the accessible mode number of the single device—the nested family of Supplement Sec. S.1.1, no spectral rela- tion across mode numbers assumed—never decreases capa- bility. The exactly matchable kernel span strictly grows: at (M ′ ,N ′ ) every block of the enlarged index set is realizable, while blocks supported outside [1,M ] n are not indexable at (M,N ) (Supplement Sec. S.3.1, Proposition 2). The statement is about spans; a capability-level separation is not claimed. Universality itself is also reachable by the standard Stone–Weierstrass route [3, 7]; we prove Theorem 1 con- structively because only the kernel route certifies scaling is never wasted, with explicit rates no density argument produces (Supplementary Discussion Sec. S.12; numerics, Sec. S.11). Gaussian-reservoir class universality [9], the nearest prior art, ranges over trainable couplings; here the swept resource is one passive geometric parameter of one fixed device, and the guarantee arrives with explicit rates and a strict-scaling corollary (Supplement Sec. S.8). What does the work is linear multimode structure with generically non-resonant frequencies, supplied here by quantization but not exclusive to it. The dichotomy is nonlinear single node with virtual nodes—the classical delay line, and why it resisted a theorem—versus linear field with independent modes, this device; one atom pack- ages encoding, mode structure and saturable nonlinearity in one passive component. The resulting objection and our answer are in Supplementary Discussion Sec. S.12. Finally, the theorem–device gap is bounded: the sat- urable device’s kernels deviate from the Gaussian-limit kernels by O(s), a computable number at each simulated operating point (Supplement Sec. S.5). 3 2.3 Genericity of the universality condi- tions Theorem 1 is only useful if realizable devices meet its con- ditions, and they hold for almost every mirror distance: each resonance function R n (τ ) = P k n k λ k (τ ) is an- alytic and provably not identically zero, with no arith- metic condition on the spectrum (the pointwise shortcut is simply false; Supplement Sec. S.6), so the bad delays form a countable union of discrete sets; overlaps fail only at standing-wave nodes, delay stability only at the dark state—likewise measure zero. We verified this on the exact device spectra: over a dense τ sweep the eigenvalues are distinct with positive real parts, overlaps stay bounded away from zero except at predicted node distances, finite-order non-resonance sur- rogates stay positive away from isolated τ, and the delay- stability margin is c 0 ≥ 0.36γ throughout the simulated regime. The Gaussian-limit devices of Fig. 2 and of the diagonalization study are verified to satisfy the overlap and non-resonance conditions on which Theorem 1 rests, while the saturable devices of the remaining figures are re- lated to the theorem through the kernel-continuity bound of Supplement Sec. S.5 rather than through membership, with measured margins reported in Supplement Sec. S.6; the designed conditions D1–D5 are a separate construc- tive certificate, and the simulated devices are generically fabricated and do not meet them. 2.4 One device, one axis: convergence with mode number The theorem and its corollary make a testable predic- tion: for a target of known Volterra structure, the sin- gle device’s error should fall as accessible mode number grows, dropping sharply once K reaches M and tracking the proof’s envelope. Figure 2 tests this in the linear- transducer regime, where the test is exact: a fixed tar- get (M = 6, N = 2, random bounded kernels), one de- vice at fixed mirror distance, mode number swept, weights trained by linear regression. Under a slope metric regis- tered before the sweep—the change in log 10 NRMSE per mode across the sampled grid—the measured error falls by a factor 11. The registered claim, that the steepest per- mode descent falls in the region K ≈ M, does not hold on the full grid: the maximum sits on the last interval, and it does so because of a single comb refinement. At K = 12 one mode lies near a standing-wave node (the per-mode overlap values are given in the caption of Fig. 2); the error rises there and the recovery beyond it produces the largest slope. Excluding refinements flagged by that overlap diagnostic alone—a criterion independent of the measured error, though one we did not register in advance—the maximum falls on the interval K = 5 → 6, exactly at K = M. We report both. The episode is it- self the genericity statement of Sec. 2.3 observed directly: overlaps fail only at isolated node distances, and one such distance is sampled here. What the device can repre- sent and what it achieves separate there: the population- optimal residual at the same readout falls by roughly six orders of magnitude across the sweep (Fig. 2b), so beyond the knee the device represents the target far more accu- rately than it achieves it and the residual error is set by the shot budget rather than by mode number. Overlaid is the same protocol on the saturable device (Bloch response, saturation depth S max = 0.2): it con- verges along the same axis, tracking the Gaussian curve closely and separating from it only beyond the knee, with an eight-fold error reduction. Its floor here is higher, as it must be: the target’s kernels are drawn from the family the Gaussian limit represents exactly, so the residual gap measures kernel distortion rather than any failure of mode scaling; on tasks outside the Gaussian family the ordering reverses (Fig. 3), which is the nonlinearity-transfer result of Sec. 2.5. Two insights surfaced (Methods): a trained readout needs only nonzero splittings, converging well below the sufficient bound, and dispersion lifts the fatal exact de- generacies. 2.5 Nonlinearity transfer: feedback re- places polynomial readout The proof lives in the Gaussian limit, where all nonlinear- ity must be supplied by the readout’s powers q n geo . Restor- ing the atom’s saturable response moves the nonlinearity from the measurement into the hardware: delayed self- interference composed with a saturable dipole generates nonlinear dependence on input history that a purely lin- ear readout harvests. Theorem and mechanism are two endpoints of one dial—where the nonlinearity resides. Figure 3 makes this quantitative: a non-Markovian reservoir with a linear readout against a Markovian atom– cavity reservoir trained to higher polynomial order, on Mackey–Glass and NARMA10, as node count grows. (One feedback node exposes two quadrature features against the comparator’s one; the resource-matched comparison, an axis relabeling anchored at the plateau, is in Meth- ods.) Despite purely linear measurements, the feedback reservoir matches or surpasses the polynomial-readout comparator, reaching by moderate node counts what the Markovian device attains only at high order or not at all; feedback strength is non-monotonic, a moderate de- lay optimal, since excessive memory retains stale informa- tion. Non-Markovian feedback thus substantially reduces, and can eliminate, the demanding requirement of high- order polynomial readouts. No tuned, noise-matched com- parison against a classical echo-state network is claimed here; the reasons, and the noiseless classical reference re- ported with the financial study, are set out in Supplement Sec. S.9. This figure also answers the sharpest objection to any minimal-hardware claim—that the regression rather than the physics is the computer. If so, the dynamics would be immaterial; instead the same linear readout on the same feature count performs very differently with feed- back (γτ = 0.01 versus 0.5, 1). 2.6 Operation on real-world tasks We evaluated the reservoir on spoken-digit recognition and tumor classification (Fig. 4b,c); financial forecasting (panel a) is analyzed in Supplement Sec. S.9, its one- step-ahead structure being closely tracked by a persis- tence predictor. Under the minimality thesis the eviden- tiary standard is parity: matching mature classical meth- 4 ods is exactly what “minimal hardware suffices” predicts, and we benchmark against standard methods and triv- ial baselines, reporting ties and shortfalls plainly. The standard is falsifiable: parity failure at achievable node counts, or performance insensitive to feedback strength, would count against the thesis. The falsifier is tested and rejected twice, by Fig. 3 and by the feedback-loss sweep of Supplement Sec. S.10, where severing the loop collapses performance to the feedback-free plateau. For speech recognition (pipeline in Methods) the reser- voir reached a word error rate of 0.142± 0.019, matching a random forest (0.143±0.009) and trailing a support-vector machine (0.096± 0.043). We report this plainly: parity with one standard method, short of the strongest. No for- mal significance claim is made, and none is attainable—at five paired folds the exact two-sided Wilcoxon test’s small- est achievable p is 0.0625, so the fold count itself precludes significance (Supplement Sec. S.9). Varying γ over 0.05– 0.25 gave word error rates of 0.174, 0.142, 0.178, 0.174, 0.192 (±0.02–0.05), with no clear monotonic trend in γ, which sets the response time. For tumor classification (multiphoton-microscopy im- ages, 48 FFPE samples near-evenly split between tumor and normal, five channels) features were prepared with the classically optimized pipeline of Ref. [25] (Haralick tex- tures, correlation pruning, exhaustive LDA subset search; Methods)—a protocol that inherently disadvantages the reservoir, which receives features selected to favor the clas- sical baseline. Performance is delay-dependent: at the robust operating point τ = 10 the reservoir classified 70.83–81.25% across subsets, below the LDA best, whereas τ = 5 ranged 58.33–87.5%, its low end no better than the majority-class baseline. On one configuration (τ = 5, one subset) it matched the best classical accuracy, 87.5% (42 of 48)—a descriptive, best-configuration tie rather than a tested equivalence (Wilson 95% interval [75.3%, 94.1%], n = 48). The paired comparison on that configuration is nonetheless exact: both classifiers scoring 42 of 48 forces the discordant pairs into balance, so McNemar’s exact test gives p = 1.0 identically, confirming no detectable differ- ence at this sample size. Because the subset search runs on the full sample, the absolute accuracies of both classi- fiers carry selection bias; the shared selection step makes the reservoir–LDA comparison, rather than either abso- lute number, the valid object (Supplement Sec. S.9). We note the limits plainly: n = 48 under one filter setting, with no stability sweep claimed. Together the tasks show a single non-Markovian atom operating alongside standard baselines: reaching the classical best on one imaging con- figuration (with the robust setting below it), and at par- ity with a random forest on acoustic data (trailing the strongest kernel method), at minimal hardware cost. 3 Discussion The results assemble into one statement about the hard- ware floor of machine intelligence. A single atom before a mirror—one active component, one geometric knob—is a provably universal approximator in its linear-transducer limit, capability never decreasing along one physical axis and the matchable span strictly growing; the saturable device converges along the same axis and outperforms the Gaussian regime on non-Gaussian benchmarks. The modes are the network, the atom the nonlinearity, the mirror delay the recurrence. The delay-line architecture, without a universality theorem in fifteen years of use, ac- quires one here for a linear-node instance with the non- linearity restored at the readout; Supplementary Discus- sion Sec. S.12 identifies the obstruction that has kept the nonlinear-node original from one. Because the scaling resource is central: modes are added by geometry, so the machine stays one piece of hardware as it scales. The theorem’s limit refines the mode spacing; a second budget governs bandwidth via the flat-coupling error ε flat (B). A superconducting atom (γ/2π ∼ 10 MHz, flat to a few percent over∼ 1 GHz) sup- ports ∼ 10 2 comb modes at percent-level ε flat , which fixes the mode spacing at ∆ 0 /2π ∼ 10 MHz and hence a quan- tization length ℓ = πv/∆ 0 ≈ 6 m at v ≈ 0.4c—an order of magnitude beyond our simulations, though not unbound- edly so (budget arithmetic, optimum, and precision costs: Methods). This arithmetic describes generic operating points, whose distinctness margins are measured rather than prescribed; the certificate’s spacing requirements are excluded from it and are unreachable at these platform parameters beyond trivial depth (Methods; Supplement Sec. S.2.5). “Arbitrarily many modes” is thus an idealiza- tion with a stated budget, as “arbitrarily many neurons” is classically. Table 1 summarizes the claim; the struc- tural point is arithmetic-free: the recurrent core has zero trained parameters and one fabricated atom. Two questions remain open: extending Theorem 1 to the saturable atom (perturbatively in s, or via analytic- ity in g), motivated by the numerical finding that this regime outperforms the Gaussian one; and whether the constructive resources—readout time T = O(M lnM + lnε −1 )/(γ + γ g ) and the exponentially graded design sep- arations whose realization costs shots—can be brought to their information-theoretic floor. The architecture maps onto existing platforms. Demonstrated jointly: the atom-before-mirror geometry with feedback-modified emission [23, 24], in the simulated regime γτ ∼ 0.5–1; separately: encoded coherent driving and homodyne detection at the required bandwidths. Tol- erances are comfortable (Supplement Sec. S.10): under 8% error at γ φ /γ = 10 −2 , jitter to ∼0.5 rad, 20 dB feedback loss recovering most of the benefit. Platform arithmetic, including the ≈ 1μs per-symbol clock, is in Methods. One structural limit belongs here. The contraction the approximation bound needs and the sub-revival horizon the uniform fading-memory bound needs cannot both hold at any mode number, the trace budget and the revival time scaling identically in it; the device is operated beyond its own revival, and the bound survives because no step of it uses fading memory (Supplement Sec. S.4). Stated conser- vatively, the contribution is threefold: an existence-and- structure theorem for one fixed machine; a generic operat- ing regime in which every simulated device meets the theo- rem’s conditions with verified margins; and a proposed ex- periment assembled from demonstrated components. The certified constants belong to the theorem rather than to the laboratory, and the difference is priced above. To close with the opening question in its defensible form: no quantum computational advantage is claimed and training remains classical, but the recurrent, high- 5 trained params (core)trained params (readout)fabricated units This device0∼ 10 2 –10 3 (linear)1 atom, 1 mirror, 1 detector, 1 fixed band envelope ESN (N=50–100) 0 (random, stored ∼ 10 3 –10 4 ) ∼ 10 2 –10 3 N digital neurons Table 1: Resource accounting for the minimality claim. The accuracy comparisons are reported in Results and are not repeated here. Against a Markovian atom–cavity comparator with polynomial readout, the device with a purely linear readout matches or surpasses ninth-order readout on Mackey–Glass and NARMA10 (Fig. 3); no tuned, noise-matched comparison against a classical echo-state network is claimed (Supplement Sec. S.9). We do not tabulate energy per inference: energy accounting for analog physical reservoirs (cryogenics amortization, detection electronics) is unsettled in the literature, and the structural comparison this table supports (what is trained, what is fabricated) does not depend on it. dimensional, nonlinear core of a learning machine can in principle be replaced by the passive dynamics of one atom and one mirror, capability purchased by bandwidth rather than component count. The purchase is metered: reach- ing accuracy through higher readout order in the Gaus- sian limit costs shots factorially in that order (Supple- ment Sec. S.7), the price the saturable device’s nonlin- earity transfer removes. The hardware floor for universal temporal learning is very low, and quantum physics puts it there. 4 Methods 4.1 Hamiltonian with non-Markovian feedback The laser-driven atom coupled to the waveguide was previ- ously modeled in [26]; we adopt the same numerical frame- work with a time-dependent input drive. The joint state evolves under iℏd t |Ψ(t)⟩ = ˆ H TOT |Ψ(t)⟩ with ˆ H TOT = ˆ H A + ˆ H F + ˆ H INT . The input is encoded in the drive on the atom, ˆ H A = ω 0 ˆσ + ˆσ − − ε INPUT 2 ε(t)(e iω L t ˆσ − + h.c.), (3) with normalized input ε(t) and strength ε INPUT . The field Hamiltonian is ˆ H F = R B dω ω ˆa † (ω)ˆa(ω) with [ˆa(ω), ˆa † (ω ′ )] = δ(ω− ω ′ ), and the interaction is ˆ H INT =i Z B dω [(g R (ω)e −iωτ/2 (4) − g L (ω)e iωτ/2 )ˆa † (ω)ˆσ − − h.c.],(5) with τ = 2L/v the round-trip delay. The idealizations re- quire the scale hierarchy ω 0 ,ω L ≫ B ≫ γ,γ g , with 1/τ in- side B so the feedback is resolved by the retained band; at the simulated parameters (γ = 0.1, τ = 10, B ∼ 2, carrier several orders above B) each inequality holds with at least an order of magnitude of margin. Regarding the input– output geometry: the mirror terminates the waveguide on one side and supports no outgoing channel; the detected field exits through the open end only, so the returning feedback (standing-wave structure) and the outgoing sig- nal (unidirectional continuum past the atom) are distinct channels, which is what licenses the exact Lindblad form of the measured channel in Sec. 2.1; residual loss, where present, is a separate unidirectional channel. Under the flat spectral density approximation g(ω) ≈ p γ/2π, mov- ing to the interaction picture and the rotating frame yields ˆ H INT,I (t) = i[( √ γ L ˆa † (t) + √ γ R ˆa † (t− τ )e iφ )ˆσ − − h.c.], (6) with ∆ = ω L − ω 0 and phase φ = π − ω L τ between the time-bin operators (φ is the round-trip phase in the frame rotating at the drive; the undriven loop’s π− ω 0 τ of Sup- plement Sec. S.4 is the same quantity, the two coinciding on resonance). The discretized propagator is ˆ U I (t k , ∆t) = exp − i ˆ H A,I (t k )∆t(7) + [( √ γ L ∆ ˆ A † (t k ) + √ γ R ∆ ˆ A † (t k−ℓ )e iφ )ˆσ − − h.c.] , with ∆ ˆ A(t k ) = R t k+1 t k dt ˆa(t) the discretized time-bin oper- ator. 4.2 Readout construction and shot bud- get In Eq. (2), t k is the start of the k-th measurement window and ρ ζ (t k ) the coherently driven reservoir state there, the subscript denoting the displaced trajectory below. The output is built from powers of the measured quadrature. The higher powers q n geo are obtained from the same ho- modyne record used for the linear quadrature, with one care required by continuous measurement: the raw pho- tocurrent contains white shot noise, whose pointwise pow- ers are not well defined, so the current is first integrated against the readout window for each sample time, yielding a well-defined, noise-corrupted outcome per window per repetition; powers of these binned outcomes are then av- eraged (Supplement Sec. S.7). Readout begins after an ini- tial washout interval—standard in reservoir computing— so that the transient of the stable dynamics has decayed and the binned outcomes carry no correlations from the initial state (Supplement Sec. S.7). No additional appara- tus is required: the shot-noise contribution to these mo- ment estimates is Gaussian and is absorbed exactly into the readout weights, and the number of shots required for precision σ on the order-n moment is bounded by (2n− 1)!! (v 0 + q 2 max ) n /σ 2 (Supplement Sec. S.7). 4.3 Adiabatic elimination and the linear- transducer limit Passing from Eq. (1) to the working generator of the proof is not an approximation of the environment (which the 6 model’s axioms already render Markovian) but the adi- abatic elimination of the atom. The induced rate is a physical quantity, fixed by the coupling and the flat band rather than by any numerical window: the band of width B has correlation time ∆t≡ 1/B, the interval over which the atom absorbs and re-emits before the field can re- spond, and the golden-rule rate of the atom-mediated drive channel over that band is γ g = 2g 2 ∆t = 2g 2 /B (correspondingly η = γ g /2g = g/B). In the simulations the time-bin width is of order this correlation time at the retained band, so the pipeline’s operational 2g 2 ∆t realizes the physical rate there, and no proof constant depends on the discretization. The limit’s hierarchy is two-sided, and we name each scale. On the fast side the band is broad, B ≫ γ,γ g : the field’s correlation time 1/B is the short- est scale, which is what licenses the golden-rule coarse- graining above, and at the simulated operating points (γ = 0.1, B ∼ 2) the ratios are B/γ = 20 and B/γ g = 3.8 at the drive-channel rate of Supplement Sec. S.2.5, so the fast-side hierarchy is the weakest of the stated inequali- ties at that operating point. On the slow side the dipole must follow the field’s envelope: its induced relaxation γ g must exceed the rates at which the driven field ampli- tude it follows evolves (the per-mode decay λ ⋆ and the drive-envelope rate 1/T on ), the same separation, in the same role, as the quasi-static condition of the Bloch route (Sec. 4.8). When both sides hold, the saturable dipole adiabatically follows the field and is replaced by a linear transducer: a c-number drive ηu(t)α together with the induced decay channel γ g 2 (α ⊗ α ∗ ). Two scopes should be kept apart. The Gaussian-limit simulations instantiate the mode-space generator directly, so for them the hier- archy is a statement about which hardware realizes that generator rather than an approximation internal to the nu- merics. The designed certificate of Supplement Sec. S.2.5 is another matter: its condition D4 drives γ g far below λ ⋆ (γ g /λ ⋆ ≈ 10 −2 at the deposited operating points), placing the certificate outside the slaving ordering; this is stated with the certificate’s other physical costs (Supple- ment Sec. S.2.5, physical resources), and generic operating points are not so constrained. The controlled small pa- rameter is the saturation s∼ (ε INPUT /γ g ) 2 , and the prop- erty the proof actually uses is the resulting Gaussianity: |ζ(t)⟩ with ζ(t k ) = P n≥1 u k−n ζ n , the moments ⟨q n geo ⟩ are polynomials in (ζ,ζ ∗ ), and the Heisenberg evolution maps a 7→ e −Γ g t · a. The retarded self-interaction responsible for the delayed feedback is recovered identically from the frequency integral of the form factor ̃g(ω) ∝ sin(ωτ/2), whose stationary-phase contributions occur at the present time and one round trip in the past. 4.4 Logic of the universality proof In the Gaussian limit the output is a finite sum of Volterra terms whose kernels factor into products of single-mode exponentials e −λ k (m−1)(T on +T off ) . The freedom in the time-dependent weights W n (t) is exactly the freedom to fix these kernels: integrating a weight against e −Λt over T off realizes any function with a well-defined Taylor ex- pansion at the eigenvalue sums, provided those sums are distinct—which the non-resonance condition guarantees; since only finitely many Taylor coefficients are ever pre- scribed, polynomial weight functions W n (t) on [0,T off ] suf- fice, and no exotic function space is invoked. (Throughout, Γ g is taken at each finite K in its biorthogonal spectral res- olution Γ g = P k λ k v R,k v ∗ L,k —guaranteed by eigenvalue distinctness under non-resonance—and no infinite-mode decomposition is invoked: the many-mode limit enters only through the uniformly decaying kernel h(t) below, so convergence of the eigenvector expansion as K →∞ is never needed; the finite-K family itself is defined from the single physical device, by a single rule and with no spectral nesting assumed, in Supplement Sec. S.1.1.) The depen- dence on memory index reduces to a Vandermonde matrix in x k = e −λ k (T on +T off ) , invertible because the eigenvalues are distinct, so all kernels up to order N and memory depth M are matched exactly whenever the device sup- ports at least M modes. Controlling the error is where the structure of the proof is decisive, and two structural facts shape it. No mode-space spectral gap uniform in K exists: the Hermitian part of Γ g has rank two, so the summed decay rates are fixed at (γ +γ g )/2 and the slowest mode closes as 1/K (Supplement, trace-identity lemma, Sec. S.2). And any bound routed through the norm of the matched weights inherits exponential growth from the node geometry (Supplement, Sec. S.2.5, Fig. S.2). The proof instead rests on a structural identity: once the kernels are matched at lags ≤ M, the realized ker- nels at all deeper lags are the spectral extrapolation of the matched block, determined by the eigenvalues alone and independent of the weights (Supplement, Extrapola- tion identity). The designed operating point—comb spac- ing locked to the readout period so the spectral nodes sit at the M-th roots of unity of a common radius ̄r = e −λ ⋆ T —then makes the extrapolation tail explicitly small, 2M ̄r M + eM ̄ε, with ̄r M = e −(γ+γ g )T/2 recovering the full trace budget over the memory span and ̄ε the node- placement error controlled by the comb spacing (Supple- ment Sec. S.2.5). The parameters are chosen once and in order—(M,N ) from the target, then the readout period, then the comb—and no infinite-mode limit is invoked: the constructed device has exactly K = M modes. 4.5 Readout resolution and the role of dispersion The constructive proof selects eigenvalue-sum points indi- vidually, which is sufficient when T off ≳ π/δ min with δ min the minimal sum splitting; a trained least-squares read- out is less demanding—it requires only that the splittings be nonzero, so that the feature span is not collapsed— and indeed Fig. 2 converges at T off = 60 although the worst-case order-two splitting of the comb corresponds to π/δ min ≈ 1.7×10 2 already at K = 6. What is fatal is exact degeneracy: for an equidistant comb the order-two sums are degenerate at leading order (split only by the non- Hermitian dressing, O((γ+γ g )/K)), the achievable feature span collapses, and performance plateaus far above the envelope. Waveguide dispersion lifts these degeneracies at leading order: a chirped comb ω k = ω 0 + ∆ 0 k + ∆ 2 k|k| splits the sums at O(∆ 2 ). In the constructive theorem this role is played instead by a single explicitly graded hi- erarchy of comb-frequency defects (Supplement Sec. S.2.5, D5), which makes the required eigenvalue-sum points with distinct signed mode content separate by construction, the generic drive-channel spread separating the remain- 7 der with measured margins; natural dispersion remains the practical mechanism on undesigned devices, and the trained-readout operation of Fig. 2 exploits it. 4.6 Detection efficiency versus feedback loss Two efficiencies must be kept apart. The feedback trans- missivity η of Supplement Sec. S.10 attenuates the return- ing field and therefore changes the dynamics—it enters r = γ R /γ in the loop characteristic equation and hence the stability margin c 0 , which is why severing it collapses the device to the feedback-free plateau. Detector inefficiency η det acts on the outgoing channel only, after the dynam- ics: it admixes vacuum into the measured field, leaving the generator, the input–output kernel and c 0 untouched, and costs only shots—the order-n moment budget of Sup- plement Sec. S.7 inflates by at most η −n det , benign at the n = 1 readout used for every real-world benchmark here (Supplement Sec. S.7). Imperfect detection therefore de- grades this device’s statistics rather than its physics; the loss that degrades its physics is loss in the feedback path, which is the quantity Supplement Sec. S.10 sweeps. 4.7 Time-bin matrix-product-state simu- lation We simulate the device with a matrix-product-state rep- resentation [27] using a time-bin (rather than spatial) decomposition, Fig. 1f. The Hilbert space is H = H S N k∈Z H k , with H S the two-level atom and each H k a Fock space per time bin, |i k ⟩ = (∆ ˆ A † ) i k |0 k ⟩/ √ ∆t i k i k !. With τ = ℓ ∆t, evolution to step k entangles time-bin pairs (t p ,t p−ℓ ) for p < k; later bins factor out as ground states, and the remaining state is decomposed into site tensors Γ and singular-value matrices Λ. The long-range propagator is implemented with SWAP gates that make each interaction local [26]. At each step the quadratures ˆ Q k = ˆ b k e −iθ + ˆ b † k e iθ and ˆ P k = i( ˆ b † k e iθ − ˆ b k e −iθ ) with ˆ b = ∆ ˆ A/ √ ∆t are extracted (physically, by homodyne de- tection) and stacked into a feature matrix X; the linear readout ˆy = WX is trained by the Moore–Penrose pseu- doinverse. The MPS convergence settings used for these simulations (maximum bond dimension χ, singular-value truncation cutoff, and time step ∆t, together with the con- vergence checks under bond-dimension and time-step re- finement) are recorded with the deposited pipeline config- uration and are not reproduced here, and this manuscript makes no claim of convergence in any of them. The com- parisons within each figure are between arms of the same pipeline at the same settings, so the orderings reported are the claims this evidence supports; the absolute error val- ues are conditional on the truncation settings, and where they stand beside classical baselines that carry no trunca- tion error the resulting parity statements should be read as conditional on those settings rather than as converged results. 4.8 Quasi-static Bloch computation of the saturable overlay (Fig. 2) The red curve of Fig. 2 (the saturable device on the con- vergence target) is computed without the MPS pipeline, by a displaced-frame quasi-static Bloch route that retains the atom’s nonlinearity while remaining closed-form in the field. We record it here so the result is reproducible from the manuscript. Displaced frame. We work in the frame displaced by the coherent field amplitude β(t) =− i g αu(t) that removes the atomic drive in favor of a field drive (as in the elim- ination of Supplement Sec. S.1), but (unlike the linear- transducer limit) we do not slave the atom to the field. Instead the atom dipole is carried as a dynamical two-level Bloch vector (⟨ˆσ x ⟩,⟨ˆσ y ⟩,⟨ˆσ z ⟩) driven by the local field, and the field modes evolve under the Gaussian generatorL 0 of Eq. (1) sourced by the atom dipole rather than by a c- number. Quasi-static approximation. When the field response time is long compared with the atom’s internal relaxation over a readout window (the same separation of scales that licenses the elimination, but retained to next order) the atom dipole follows its instantaneous steady state on the field’s slow timescale. We solve the steady-state Bloch equations for⟨ˆσ − ⟩ as a function of the instantaneous drive, including saturation through ⟨ˆσ z ⟩ = −1/(1 + S) with S the saturation set by the drive and detuning, and feed the resulting nonlinear dipole into the mode dynamics. The controlled small parameter is the ratio of field-response to atomic-relaxation time; at the operating drive the satu- ration depth is S max = 0.2 and ⟨ˆσ z ⟩ is driven from −1 to −0.833 (an excited-state population of 0.083; the atom remains predominantly in the ground state, and no in- version is implied), so the atom is driven appreciably out of the linear-transducer regime—which is the point of the overlay. Validation. The quasi-static dipole was checked against the full time-dependent solution of the coupled atom–field dynamics on a set of (K, drive) points of the sweep whose size and composition are recorded with the deposited script, so the deviation quoted is a maximum over that set rather than a bound over the sweep; the maximum relative deviation of the readout feature over the validation set is 4× 10 −3 , small compared with the eight-fold error reduction the curve reports, so the approx- imation does not affect the convergence conclusion. The script implementing this route and reproducing Fig. 2 is included in the repository script-to-figure map. 4.9 Feedback-free comparison model For the polynomial-readout comparison of Fig. 3 we use a Jaynes–Cummings reservoir in which the atom is coupled to a single cavity mode ˆc, H 0 = ω c ˆc † ˆc + ∆ˆσ + ˆσ − + g(ˆc † ˆσ − + ˆcˆσ + ),(8) H 1 (t) = i ε INPUT 2 ε(t)(ˆσ − − ˆσ + ).(9) Consistent with the mirror device in its Markovian limit (γτ ≪ 1), the returning field is re-injected essentially in- stantaneously and the cavity mode decays at the mirror- modified rate κ = γ(1 + cosφ) = γ(1− cosω L τ ),(10) so the master equation is ̇ρ =−i[H 0 +H 1 (t),ρ] +D[ √ κ ˆc]ρ with the dissipator acting on the measured cavity mode ˆc, and the readout observables ˆ Q = ˆc + ˆc † , ˆ P = i(ˆc− ˆc † ), ˆ P 2 , 8 ˆ P ˆ Q, ˆ Q 2 ,... are those of the same mode. A Lamb shift δω = γ 2 sinφ accompanies the decay, with φ = π − ω L τ the round-trip phase defined above. This makes the dis- sipative channel and the measured observable consistent, which is required for a fair comparison against the feed- back device. The comparator set tests the feedback abla- tion (this Markovian device) and a classical network ref- erence (the echo-state network of Supplement Sec. S.9); a classical single-node delay-line reservoir at matched virtual-node counts is a natural comparison not performed here, and the claim this work makes about that architec- ture concerns the availability of a single-member theorem rather than relative performance. 4.10 Scope of the model and the sense of exactness A central technical point organizes the paper. “Non- Markovian” here refers to the reduced dynamics of the atom: because photons return after a delay, the atom in- teracts with its own past. The joint atom–field system evolves under a time-independent Hamiltonian, and once the field is written in the standing-wave modes selected by the mirror, delay and feedback appear entirely as coher- ent structure among those modes. The only irreversible processes (the field leaving through the open end, and the homodyne measurement) act on unidirectional con- tinua that never return and are therefore Markovian to the same standard as the axioms themselves. Under the standard input–output (SLH) idealizations—stated as ax- ioms A1–A6 in Supplement Sec. S.1, and including the classical treatment of the input drive and the idealized continuous measurement—the device is consequently de- scribed by a single Lindblad generator that truncates none of the delayed feedback, and our theorem is a statement about that generator. We are careful about the sense of “exact”: writing a Lindbladian presupposes those axioms rather than deriving them from an unspecifiable environ- ment. What is not approximated, relative to the usual treatment of delayed feedback, is the memory itself—no Born–Markov truncation of the retarded self-interaction, no pseudomodes, no auxiliary lossy modes. The one fur- ther approximation is of the atom rather than of the en- vironment, a linear-transducer (Gaussian) limit, and we show numerically that operating outside it, where the atom’s own nonlinearity is active, reduces the hardware’s measurement requirements further on the tasks tested. Extending the proof to the saturable atom is an open prob- lem we state explicitly. 4.11 Resource matching for the Marko- vian comparator (Fig. 3) To be explicit about the resource counted, since the two devices expose different observables: for the feedback reservoir one “node” is one measured quadrature pair per time bin; for the comparator it is one polynomial read- out feature. The map is: feedback reservoir, n nodes = n time bins = 2n features (Q and P per bin); compara- tor, n nodes = n features. Resource-matching is then an axis relabeling—equal feature count places n feedback nodes against 2n comparator nodes, shifting the feedback curves by a factor of two along the node axis (rightward on the logarithmic axis of Fig. 3), which any reader can perform from this map. Since the separation at large node count is an asymptotic vertical gap (both families having plateaued), a horizontal rescaling by two leaves the order- ing, and the conclusion, unchanged. We anchor the claim at the plateau for exactly this reason: before both families saturate, a horizontal factor of two need not preserve the ordering. 4.12 Repetitions A repetition is one independent experimental shot of the input–readout cycle. Because the readout estimates mo- ments of a measured quadrature, each reported output sample is an average over R repetitions, and R is a re- source on the same footing as the node count: the shot budget of Supplement Sec. S.7 fixes R from the target rel- ative precision. Where a benchmark figure reports more than one value of R, the value used for the quoted result is stated in the figure legend. 4.13 Error metric Throughout, NRMSE denotes root-mean-square error nor- malized by the target range (maximum minus minimum) over the test set, the convention of Supplement Secs. S.9– S.10 as well, so values are comparable across figures. Range normalization yields systematically smaller val- ues than the standard-deviation normalization common in the benchmark literature, so comparison with pub- lished NARMA10 figures requires conversion. No pub- lished NARMA10 or Mackey–Glass value is invoked any- where in this work: every comparator reported here is computed within this study under the identical normal- ization. 4.14 Speech pipeline For speech recognition we used 500 utterances of the Free Spoken Digit Dataset, converting each to a mel- spectrogram reduced to per-bin mean and standard- deviation streams injected into two reservoirs, with a 144-dimensional concatenated state, ten one-vs-rest read- outs, and a winner-takes-all decision under five-fold cross- validation. Word error rates, the γ sweep, and the signifi- cance arithmetic are reported in Results, with the paired- analysis protocol of record in Supplement Sec. S.9. 4.15 Tumor feature pipeline For tumor classification, features were prepared following the feature pipeline of Ref. [25]: 13 Haralick texture fea- tures per channel, so 65 features per sample; the unit of analysis throughout is the FFPE sample, and whether sev- eral samples derive from one individual is a property of the archived dataset that this secondary analysis does not re- solve. The 65 features were first reduced to 22 by removing highly correlated features at a pairwise-correlation thresh- old of 0.45; all 22 4 = 7,315 four-feature subsets were then evaluated with linear discriminant analysis, and the four best-performing subsets were retained. Only these clas- sically optimized subsets were subsequently encoded and processed by the reservoir. 9 4.16 Physical budgets of the scaling re- source Because the scaling resource is central, we state its bud- gets in full. The theorem’s limit refines the mode spac- ing: increasing the quantization length ℓ grows the mode number K and the revival time T P = 2π/∆ 0 together, which is what places the memory horizon inside the sub- revival window where the kernel bound holds; that bud- get is set by ℓ and by coherence over the round trip. A second, independent budget governs bandwidth: the flat- coupling axiom degrades as the coupling density varies across the occupied band B = K∆ 0 , with relative error ε flat (B) ∼ sup |ω−ω 0 |≤B/2 |g(ω)/g(ω 0 ) − 1|. Refining the spacing at fixed B does not pay this cost; widening B at fixed spacing does, and it caps how fine a comb a given bandwidth supports. We are precise about where this bud- get bites: the simulations reported here impose the flat- coupling axiom exactly (g(ω) = p γ/2π over the retained band), so their operating ε flat is zero by construction, and the numbers below are a statement about a physical device rather than a correction to our numerics. The device’s in- trinsic fading-memory rate c 0 , set by loop stability rather than the spectrum, degrades with neither budget; the the- orem’s convergence rate is a third quantity—e −(γ+γ g )T/2 over the memory span—purchased by readout time, and it too is independent of both budgets. A superconducting artificial atom with γ/2π ∼ 10 MHz on a line flat to a few percent over ∼ 1 GHz supports ∼ 10 2 comb modes before ε flat reaches the percent level—an order of magnitude be- yond the K ≤ 14 of our convergence experiment and the ∼ 30 nodes of our benchmarks, though not unboundedly so. Two companion statements complete the budget. First, the trade-off has an optimum: the constructive the- orem’s node-placement error ̄ε falls as the comb spacing grows, but a wider comb occupies more band and pays ε flat (B); the optimal operating point is where the falling ̄ε meets the growing ε flat (B), and it is the honest physical boundary of the idealization: the one step of the construc- tion where the proof meets the hardware. Second, for the constructive certificate specifically: the designed operating point’s defect-protection condition demands a comb spac- ing that grows exponentially with the memory depth, and its timing condition demands a relative precision of order ̄ε/(∆ 0 T ) (Supplement Sec. S.2.5, scaling ledger); a gener- ically fabricated device pays neither price, separating the same eigenvalue sums by its own verified conditions, which is how every simulation in this paper is operated. Growing precision requirements on settings are the standard cost of universal-approximation constructions, exactly as classical weight precision grows with accuracy, and we state them rather than leave them implicit. The certificate’s reach can be stated as one number. Even at its cheapest admissible selections (N = 1, minimal readout period, maximal node budget ̄ε target = 1/4), the binding entry of D1—the radial-protection entry, which exceeds the defect-protection entry by a factor growing like (16/5) M —together with the band constraint M ∆ 0 ≤ B requires B/γ ≥ 8192M 2 q 2 ln(2M ) 4 q−2 with q the least prime ≥ 2K phys + 1: about 7.3 × 10 7 at M = 2. The superconducting example above (B/γ ≈ 10 2 ) therefore supports the certified construction at no depth M ≥ 2; realizing even M = 2 requires kHz-class linewidths under GHz-flat coupling, a combination not among the demon- strated platforms. The certificate is a mathematical exis- tence statement with explicit constants; the physically op- erable regime is the generic one, which satisfies no designed spacing entry and whose conditions every simulation veri- fies with measured margins (Supplement Sec. S.2.5, phys- ical resources). 4.17 Experimentalfeasibilityand throughput To fix the scale: the simulated feedback regime γτ ∼ 0.5–1 at the superconducting example above (γ/2π ∼ 10 MHz) corresponds to a round-trip delay τ ≈ 8–16 ns, i.e. an atom–mirror path L = vτ/2 ≈ 0.5–1 m at a typ- ical coplanar-waveguide phase velocity v ≈ 0.4c—a me- andered delay line of ordinary size rather than an ex- otic requirement. Demonstrated jointly, in single exper- iments: the atom-before-mirror geometry with distance- dependent, feedback-modified emission, in superconduct- ing circuits [24] and with a trapped atom before a mirror [23], in the same feedback regime we simulate (γτ ∼ 0.5– 1) [20–22]. Demonstrated separately on the same plat- form classes, but not yet combined with mirror feedback in one apparatus: coherent driving with encoded input se- quences, and homodyne detection at the bandwidths our readout requires. To be concrete about that requirement: the homo- dyne chain must pass the occupied comb band B = K∆ 0 about the carrier—up to ∼ 1 GHz at the 10 2 -mode ceil- ing above, correspondingly less at the simulated mode counts—together with windowed integration at the ∼ μs symbol clock, both standard for microwave homodyne de- tection. Integrating these into one operated reservoir com- puter is the experimental step this proposal calls for; no individual capability is new. The scheme does not require the coherence budget of a qubit register, and its toler- ance to imperfections is quantified directly in Supplement Sec. S.10: dephasing degrades performance gracefully (un- der 8% relative error at γ φ /γ = 10 −2 , within reach of the platforms above), phase jitter is tolerated to at least ∼0.5 rad (a round-trip path stability of ∼ λ/13, equiva- lently∼ λ/25 of the one-way distance L), and the feedback path is loss-tolerant—restoring 1% of the returning power (20 dB loss) recovers most of the benefit relative to a sev- ered loop. That last sweep uses single deterministic-input realizations and carries no seed spread, so its loss toler- ance is quoted qualitatively here and in the Supplement, with no interval attached. Two of these are more than robustness checks: the per- formance floor over static phases falls exactly at the dark- state value φ = π predicted by the delay-stability analysis, a task-level confirmation of the stability landscape; and severing the loop executes the paper’s own falsifier. Throughput completes the practical picture. The per- symbol clock is T = T on +T off , dominated by the readout window; at the operating point of Fig. 2, T off = 60/γ, which at the superconducting example’s γ/2π ∼ 10 MHz is≈ 1μs per symbol, i.e. a repetition rate of order 10 6 per second. At matched statistical precision this becomes an output rate of order 10 2 samples per second for the linear readout at one-percent relative precision, since each out- 10 put sample consumes N shots repetitions; the single-shot photonic delay-line reservoir of Ref. [15] is quoted at the former rate and the comparison should be made at the latter. Two scalings qualify the number in opposite direc- tions. Toward harder targets, the constructive theorem spends readout time T = O(M lnM + lnε −1 )/(γ +γ g ), set by the trace budget rather than the per-mode rate, so the clock grows as Θ(M lnM ) in the memory depth, as the scaling ledger records; and the readout-resolution analysis of Sec. 2.4 shows that with dispersion the required T off does not grow with K, so adding modes costs bandwidth rather than time. Toward statistics, the quoted rate is per experimental repetition: moment estimation at pre- cision σ multiplies wall-clock time by N shots (Supplement Sec. S.7), which is the same ensemble cost every analog reservoir pays and is kept at its n = 1 floor here by the nonlinearity transfer. Data Availability The datasets are publicly available.The tumor dataset is at https://doi.org/10.25422/azu.data. 29983060.The spoken-digit dataset is the Free Spoken Digit Dataset (FSDD), https://github.com/ Jakobovski/free-spoken-digit-dataset, distributed under the Creative Commons Attribution–ShareAlike 4.0 International licence; FSDD is an open, growing corpus versioned by Zenodo DOI and git tag, and the version tag used here is recorded with the analysis code so the 500- utterance subset is reproducible. Financial data are daily closing prices for the tickers GSPC (S&P 500), AAPL (Ap- ple), and IXIC (NASDAQ) over 1 January 2014 to 1 Jan- uary 2024, obtained from Yahoo Finance. Because finan- cial data providers revise historical series, the fetch script, the ticker list, the adjustment settings, the download date, and SHA-256 checksums of the exact CSV extracts used here are deposited with the code, so a reader who obtains the deposited extracts can confirm they are the ones used here by comparing those checksums. Code Availability All simulation and analysis code, including the scripts generating every figure (with fixed seeds and a script-to-figure map in the repository README), is available at https://github.com/nhula01/ nonMarkovianReservoirComputer; this manuscript does not carry a release identifier, and a reader wishing to con- firm which revision produced these results must match the deposited configuration against the repository history. 11 Figure 1: Overview of the minimal quantum reservoir computer performing real-world tasks. a, A single atom with transition energyℏω 0 sits a distance L in front of a mirror, driven by the discrete input ε(k) and decaying into the waveguide to the left and right at rates γ L ,γ R . The round-trip delay τ = 2L/v produces coherent feedback. b–c, Homodyne detection collects quadrature nodes in each time bin. d, A classical linear readout W predicts the target y k . e, The three real-world tasks: financial prediction, speech recognition, and tumor classification. f, Time-bin matrix-product-state scheme (Methods). With delay ℓ = 4, at k = 0 a sequence of SWAP gates brings time bin −4 adjacent to the atom, the local evolution ˆ U I (t 0 ) is applied, and the bins are returned to order; a finalSWAP advances k = 0→ 1. The process repeats to the final time. 12 Figure 2: Convergence of a single device with mode number. a, Test-set approximation error (median over 12 input realizations; band, interquartile range) on a target with known Volterra kernels (M = 6, N = 2, random bounded kernels) versus the mode number K of one device with a nested dispersive comb at fixed retained band and fixed mirror distance, so that ∆ 0 = B/K and the readout node count is held fixed across the sweep. Blue: linear-transducer (Gaussian) regime, exact closed-form dynamics; red: the saturable device on the identical protocol (displaced-frame quasi-static Bloch route of Methods, quasi-static Bloch computation; S max = 0.2, ⟨σ z ⟩ driven from −1 to −0.833). Both carry shot noise at N shots = 3× 10 4 , the one-percent relative precision of Methods. The non-monotonicity at K = 12 is a device effect rather than scatter: at that comb refinement one mode sits near a standing-wave node, √ K min k |α k | = 0.0067 against 0.077 at the next lowest refinement and 0.40 at K = 10, the isolated overlap failure the genericity analysis predicts. In c the two intervals adjoining that refinement are greyed; excluding it on the overlap diagnostic alone leaves the maximum on K = 5 → 6. b, The population-optimal residual over the same feature set, i.e. what the device can represent rather than what it achieves: it falls by roughly six orders of magnitude across the sweep, dropping sharply at K = 5 and running from two to five orders below the shot-limited error of a thereafter, the exact-realizability statement of the theorem. Beyond that point the measured error in a is limited by the shot budget rather than by mode number. c, The registered slope metric, −∆ log 10 NRMSE per mode across the sampled grid; the maximum falls in the region K ≈ M (orange). The metric, the grid, the seed count and the claim under test were fixed before the sweep was run; the regulariser is selected on a validation split and never on test. 13 Figure 3: Nonlinearity transfer. A non-Markovian reservoir with a purely linear readout (feedback γτ = 0.5, 1) is compared against a Markovian atom–cavity reservoir trained to increasing polynomial order (vertical guides mark the polynomial degree at which the Markovian device gains each corresponding node) and against the Markovian- limit feedback γτ = 0.01, on Mackey–Glass (top) and NARMA10 (bottom). Parameters: ∆t = 1, θ = φ = π/3, g = (γ/2(1 + cosφ),γ/4(1 + cosφ)) with φ = π−ω L τ the round-trip phase entering the comparator’s mirror-modified decay (Methods, “Feedback-free comparison model”), averaged over ε INPUT = (0.8, 1.0, 1.2, 1.4, 1.6)γ. The linear- readout feedback reservoir matches or surpasses ninth-order polynomial readout. Solid lines: mean NRMSE over the five input strengths ε INPUT ; shaded bands: ±1 s.e.m. over those strengths. This spread is over a swept device parameter rather than over repeated runs; no seed-resolved variability is reported for this figure, and the comparison is stated as an ordering rather than as a tested difference. 14 Figure 4: The single-atom reservoir across three real-world tasks. a, Financial prediction with a rolling two-year-train / one-year-test protocol, reported by NRMSE; parameters ∆t = 1, γ = 0.1, ε = 0.15, θ = φ = π/3, τ = 15, with a delay sweep τ = 5, 10, 15, 20. The quoted τ = 15 is the setting shared with panels (b,c) rather than this task’s optimum—shorter delays score better (0.042± 0.008 at τ = 5; Supplement Sec. S.9)—and is retained so the three panels report one device at one operating point. A persistence baseline ˆy t = y t−1 scores ≈ 0.037 on these identical windows, below both the reservoir and the echo-state comparator, which is why this task is a consistency check be- tween methods rather than a pillar (Supplement Sec. S.9). b, Speech recognition on the Free Spoken Digit Dataset, reported by word error rate (WER); pipeline in Methods. γ = (0.05, 0.10, 0.15, 0.20, 0.25), other parameters as in (a) with τ = 15. c, Pancreatic neuroendocrine-tumor classification from multiphoton-microscopy Haralick features (13 per channel, 65 per sample; correlation-pruned to 22, then the four LDA-optimal four-feature subsets retained and encoded), leave-one-out evaluation; γ = 0.1, ε = 0.1, θ = φ = π/3. Across the four LDA-selected subsets τ = 10 is the robust setting (70.83–81.25%) relative to τ = 5 (58.33–87.5%). Because the subsets are chosen to optimize the classical baseline, the comparison disadvantages the reservoir. 15 Declarations Ethics. The pancreatic-tumor study uses the previously published, publicly archived multiphoton-microscopy dataset cited in Methods; the original specimens were collected under the institutional approvals reported by that study: University of Michigan Endocrine Oncology Repository IRB #HUM00115310, Icahn School of Medicine at Mount Sinai Biorepository and Pathology Core IRB STUDY-12-00145, and University of Arizona Tissue Acquisi- tion and Cellular/Molecular Analysis Shared Resource IRB #0600000609; and the present work performs secondary computational analysis of the de-identified public archive only, involving no new human-subject data collection. Author contributions. D.S. conceived and supervised the project and formulated the universality claim and the design conditions. P.J.E. and P.H.N. developed the theoretical framework and the exact mode-space representation, contributed equally, and share first authorship. P.H.N. implemented and ran the matrix-product-state simulation pipeline and the benchmark experiments. P.J.E. developed the proof of universality and the associated analysis. K.S. contributed to the waveguide-QED modeling. N.D. contributed to the numerical experiments. T.W.S. provided the multiphoton-microscopy tumor dataset and its clinical and provenance context. H.I.N. provided technical analysis of the reservoir-computing theory and the moment-readout treatment. All authors discussed the results and contributed to the manuscript. Competing interests. The authors declare no competing interests. Acknowledgment. This work was supported by the U.S. Department of Energy, Office of Science, Award No. DE-SC0025910 (D.S.), and by the National Science Foundation, Award No. 2529700 (D.S.). Use of AI tools. We used a large language model (Claude, Anthropic) to verify author-developed mathematical derivations and to improve the readability and grammar of the text. All the results presented were reviewed and approved by the authors, and the authors take full responsibility for all content. References [1] Herbert Jaeger and Harald Haas. Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. Science, 304(5667):78–80, 2004. doi: 10.1126/science.1091277. URL https://w. science.org/doi/10.1126/science.1091277. [2] Wolfgang Maass, Thomas Natschläger, and Henry Markram. Real-time computing without stable states: A new framework for neural computation based on perturbations. Neural Computation, 14(11):2531–2560, 2002. doi: 10.1162/089976602760407955. [3] Lyudmila Grigoryeva and Juan-Pablo Ortega. Echo state networks are universal. Neural Networks, 108:495–508, 2018. ISSN 0893-6080. doi: https://doi.org/10.1016/j.neunet.2018.08.025. URL https://w.sciencedirect. com/science/article/pii/S089360801830251X. [4] Lukas Gonon and Juan-Pablo Ortega. Fading memory echo state networks are universal. Neural Net- works, 138:10–13, 2021. ISSN 0893-6080. doi: https://doi.org/10.1016/j.neunet.2021.01.025. URL https: //w.sciencedirect.com/science/article/pii/S0893608021000332. [5] Lyudmila Grigoryeva and Juan-Pablo Ortega. Universal discrete-time reservoir computers with stochastic inputs and linear readouts using non-homogeneous state-affine systems. Journal of Machine Learning Research, 19(24): 1–40, 2018. URL http://jmlr.org/papers/v19/18-020.html. [6] Peter J. Ehlers, Hendra I. Nurdin, and Daniel Soh. Stochastic reservoir computers. Nature Communications, 16: 3070, 2025. doi: 10.1038/s41467-025-58349-6. [7] Jiayin Chen and Hendra I Nurdin. Learning nonlinear input–output maps with dissipative quantum systems. Quantum Information Processing, 18(7):198, 2019. [8] Jiayin Chen, Hendra I. Nurdin, and Naoki Yamamoto. Temporal information processing on noisy quantum computers. Physical Review Applied, 14(2):024065, 2020. doi: 10.1103/PhysRevApplied.14.024065. [9] Johannes Nokkala, Rodrigo Martínez-Peña, Gian Luca Giorgi, Valentina Parigi, Miguel C Soriano, and Roberta Zambrini. Gaussian states of continuous-variable quantum systems provide universal and versatile reservoir computing. Communications Physics, 4(1):53, 2021. [10] Keisuke Fujii and Kohei Nakajima. Harnessing disordered-ensemble quantum dynamics for machine learning. Phys. Rev. Applied, 8:024030, Aug 2017. doi: 10.1103/PhysRevApplied.8.024030. URL https://link.aps.org/ doi/10.1103/PhysRevApplied.8.024030. 16 [11] Christa Cuchiero, Lukas Gonon, Lyudmila Grigoryeva, Juan-Pablo Ortega, and Josef Teichmann. Discrete-time signatures and randomness in reservoir computing. IEEE Transactions on Neural Networks and Learning Systems, 33(11):6321–6330, 2022. doi: 10.1109/TNNLS.2021.3076777. [12] Daniel J Gauthier, Erik Bollt, Aaron Griffith, and Wendson AS Barbosa. Next generation reservoir computing. Nature communications, 12(1):5564, 2021. [13] Gouhei Tanaka, Toshiyuki Yamane, Jean Benoit Héroux, Ryosho Nakane, Naoki Kanazawa, Seiji Takeda, Hidetoshi Numata, Daiju Nakano, and Akira Hirose. Recent advances in physical reservoir computing: A review. Neural Networks, 115:100–123, 2019. ISSN 0893-6080. doi: 10.1016/j.neunet.2019.03.005. URL https://w.sciencedirect.com/science/article/pii/S0893608019300784. [14] Lennert Appeltant, Miguel C. Soriano, Guy Van der Sande, Jan Danckaert, Serge Massar, Joni Dambre, Benjamin Schrauwen, Claudio R. Mirasso, and Ingo Fischer. Information processing using a single dynamical node as complex system. Nature Communications, 2:468, 2011. doi: 10.1038/ncomms1476. [15] Laurent Larger, Antonio Baylón-Fuentes, Romain Martinenghi, Vladimir S. Udaltsov, Yanne K. Chembo, and Maxime Jacquot. High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification. Physical Review X, 7(1):011015, 2017. doi: 10.1103/PhysRevX.7.011015. [16] Lyudmila Grigoryeva, Julie Henriques, Laurent Larger, and Juan-Pablo Ortega. Stochastic nonlinear time series forecasting using time-delay reservoir computers: performance and universality. Neural Networks, 55:59–71, 2014. doi: 10.1016/j.neunet.2014.03.004. [17] Felix Köster, Serhiy Yanchuk, and Kathy Lüdge. Master memory function for delay-based reservoir computers with single-variable dynamics. IEEE Transactions on Neural Networks and Learning Systems, 35(6):7712–7725, 2024. doi: 10.1109/TNNLS.2022.3220532. Preprint: arXiv:2108.12643. [18] Silvia Ortín, Miguel C. Soriano, Luis Pesquera, Daniel Brunner, Daniel San-Martín, Ingo Fischer, Claudio R. Mirasso, and José M. Gutiérrez. A unified framework for reservoir computing and extreme learning machines based on a single time-delayed neuron. Scientific Reports, 5:14945, 2015. doi: 10.1038/srep14945. [19] Florian Stelzer, André Röhm, Raúl Vicente, Ingo Fischer, and Serhiy Yanchuk. Deep neural networks using a single neuron: folded-in-time architecture using feedback-modulated delay loops. Nature Communications, 12: 5164, 2021. doi: 10.1038/s41467-021-25427-4. [20] U. Dorner and P. Zoller. Laser-driven atoms in half-cavities. Phys. Rev. A, 66:023816, Aug 2002. doi: 10.1103/ PhysRevA.66.023816. URL https://link.aps.org/doi/10.1103/PhysRevA.66.023816. [21] Kanupriya Sinha, Alejandro González-Tudela, Yong Lu, and Pablo Solano. Collective radiation from distant emitters. Phys. Rev. A, 102:043718, Oct 2020. doi: 10.1103/PhysRevA.102.043718. URL https://link.aps. org/doi/10.1103/PhysRevA.102.043718. [22] Tommaso Tufarelli, M. S. Kim, and Francesco Ciccarello. Non-markovianity of a quantum emitter in front of a mirror. Phys. Rev. A, 90:012113, Jul 2014. doi: 10.1103/PhysRevA.90.012113. URL https://link.aps.org/ doi/10.1103/PhysRevA.90.012113. [23] Jürgen Eschner, Christoph Raab, Ferdinand Schmidt-Kaler, and Rainer Blatt. Light interference from single atoms and their mirror images. Nature, 413(6855):495–498, 2001. doi: 10.1038/35097017. URL https://w. nature.com/articles/35097017. [24] Io-Chun Hoi, A. F. Kockum, Lars Tornberg, Arsalan Pourkabirian, Göran Johansson, Per Delsing, and C. M. Wilson. Probing the quantum vacuum with an artificial atom in front of a mirror. Nature Physics, 11(12): 1045–1049, 2015. doi: 10.1038/nphys3484. URL https://w.nature.com/articles/nphys3484. [25] Noelle Daigle, Shuyuan Guan, Suzann Duan, Thomas G. Knapp, Eungjoo Lee, Ali Azhdarinia, Sukhen C. Ghosh, Solmaz AghaAmiri, Servando Hernandez Vargas, Naruhiko Ikoma, Jeannelyn Estrella, Martin J. Schnermann, Tobias Else, Michelle Kang Kim, Juanita L. Merchant, and Travis William Sawyer. Investigating machine learning algorithms to classify label-free images of pancreatic neuroendocrine neoplasms. Biophotonics Discovery, 2(4): 045001, 2025. doi: 10.1117/1.BIOS.2.4.045001. URL https://doi.org/10.1117/1.BIOS.2.4.045001. [26] Hannes Pichler and Peter Zoller. Photonic circuits with time delays and quantum feedback. Phys. Rev. Lett., 116:093601, Mar 2016. doi: 10.1103/PhysRevLett.116.093601. URL https://link.aps.org/doi/10.1103/ PhysRevLett.116.093601. [27] Román Orús. A practical introduction to tensor networks: Matrix product states and projected entangled pair states. Annals of Physics, 349:117–158, October 2014. ISSN 0003-4916. doi: 10.1016/j.aop.2014.06.013. URL http://dx.doi.org/10.1016/j.aop.2014.06.013. 17 Supplementary Information for “A Single Atom in Front of a Mirror is a Universal Reservoir Computer” Peter J. Ehlers, Phi Hung Nguyen, Kanu Sinha, Noelle Daigle, Travis W. Sawyer, Hendra I. Nurdin, Daniel Soh This Supplement provides the full derivations underlying the main text. Section S.1 derives the mode-space Lindbladian and the linear-transducer (Gaussian) limit, states the input–output axioms A1–A6 relative to which the delayed feedback is carried without truncation, identifies the adiabatic elimination of the atom as the one further approximation, and defines the nested finite-K family from the single physical device by a single rule (Sec. S.1.1), stating explicitly that no spectral nesting across K is assumed and locating the two tiers (readout-side in-band selection, envelope-controlled band truncation) at which the operating mode number is set. Section S.2 gives the complete proof of universality. The error bound is shaped by a structural constraint on the device: the trace of the generator’s Hermitian part is pinned, so the mode-space spectrum has no uniform gap and its slowest rate closes as 1/K (Lemma 3). The proof accordingly rests on a weight-independent extrapola- tion identity (Proposition 1): once the kernels are matched at the working depth, the device’s deeper response is fixed by its spectrum alone, and an explicitly constructed operating point (Sec. S.2.5.1) makes that spectral tail exponentially small, with every constant explicit and no infinite-mode limit invoked. Delay stability is treated sep- arately in Sec. S.4, as the fading-memory property of the device itself (Lemma 11, Corollary 2); Sec. S.3.2 verifies the hypotheses of the density theorem we invoke. Section S.5 bounds the deviation between the Gaussian-limit kernels and those of the fully saturable device by the saturation parameter, and defines the invariant bounded-excitation sector on which that bound is controlled. Section S.6 proves the genericity of the universality conditions (by an argument based on functional rather than pointwise independence) and records the numerical verification, including of the delay-stability margin. Section S.7 proves that shot-noise-corrupted moment estimates suffice and bounds the required shot number. Section S.8 locates the result against prior constructions in tabular form. Section S.9 presents the benchmark numerics; Section S.10 quantifies robustness to atomic dephasing, round-trip phase jitter, and 1 feedback loss, including a direct dynamical confirmation of the delay-stability land- scape and a direct execution of the feedback-insensitivity falsifier; Section S.11 presents the property-verification numerics. Section S.12 is a Supplementary Discussion car- rying, in full, the Stone–Weierstrass route and the classical-linear-optics objection summarized in the main text. S.1 The Mode-Space Lindbladian and the Linear-Transducer Limit The full system is a single atom in a multimode field described, in the standing-wave basis selected by the mirror, by the Lindbladian below, where ρ denotes the joint atom–field density operator: ̇ρ = Γ· aρ,a † + a,ρa † · Γ † (S1) − i[gσ + (α ∗ · a) + gσ − (α· a † ) + εσ + σ − ,ρ] + u(t)[σ + − σ − ,ρ], where a is the vector of mode annihilation operators, σ ± the atomic operators, g the atom–field coupling, α the emitted-mode profile with α k ∝ sin(ω k τ/2), and ε the excited-state energy. The matrix Γ = γ 2 (v geo ⊗ v ∗ geo ) + iΩ contains the mode frequencies Ω (Hermitian, diagonal) and the non-unitary evolution from continuous weak measurement of the channel v ∗ geo · a. Epistemic status of Eq. (S1). Because the strength of the claim depends entirely on what is being held fixed, we state the model’s axioms explicitly and then say precisely what is and is not exact relative to them. The axioms are the standard input–output idealizations of waveguide QED: A1. Rotating-wave approximation for the atom–field coupling. A2. Flat coupling density over the relevant bandwidth B about the atomic transition, g(ω) ≈ p γ/2π, with the fixed Gaussian edge rolloff of Lemma 11; the residual variation across the band is the budget ε flat (B) of the main text. A3. Unidirectional outcoupling : the detected and lost channels are outgoing continua past the atom that never return to it. The quantization length ℓ is a physical termination of the guided line rather than a bookkeeping device: it is what makes the retained spectrum a discrete comb of spacing ∆ 0 = πv/ℓ, and the conditions D1, D3 and D5 are conditions on those physical resonances. The far termination is part of the coherent structure the axiom does not exclude; what the axiom excludes is return through the measured port. A4. Structureless vacuum for those outgoing channels: they are initially in vacuum, are uncorrelated with the atom–field system at the initial time, and carry no internal structure on the timescales resolved by B. A5. Classical input drive: the input field is treated as a c-number drive u(t) on the atom, i.e. a coherent input in the standard input–output sense, with its quantum fluctuations subsumed into the vacuum of the drive channel. A6. Idealized continuous measurement : the homodyne detection of the outcoupled chan- nel is described in the standard continuous-measurement framework—the detector’s 2 internal dynamics is eliminated, the local oscillator is classical, and the measured channel propagates unidirectionally to the detector. A7. Two-level validity : the emitter’s nearest non-resonant transition is detuned from the retained band by much more than the band width, B ≪ |∆ anh |, with ∆ anh the anharmonicity of an artificial atom or the nearest-level detuning of a natural emitter. The mode ceiling of a given platform is then the smaller of the flat-coupling ceiling and c|∆ anh |/∆ 0 with c < 1; for a weakly anharmonic artificial atom the second branch binds, while for a natural emitter, ion or colour centre it is inactive. These are exactly the assumptions of the SLH / input–output formalism of open quantum networks [1, 2], and we adopt them as such: the drive and measurement treatments of A5–A6 are idealizations of the same standing as A1–A4 rather than consequences of the model. Given A1–A6, two facts fix the status of Eq. (S1). First, the standing-wave decomposition is a unitary change of basis of the free field, so the delay and feedback are carried entirely by the coherent structure (Ω,α) with no auxiliary or lossy modes. Second, the measured and outcoupled channels act on unidirectional continua that never return, so their delta-correlation—and hence the Lindblad form of the Hermitian rank-one terms—follows from A3, A4, and A6 jointly, and it is a statement inside the SLH framework rather than one that survives outside it. The claim we make is therefore a relative one, and we are careful not to overstate it. Writing a Lindblad generator at all presupposes A1–A6; it is not derived from an unspecified environment, and no such derivation is available in principle, since the en- vironment of a laboratory device is never fully specifiable. In particular, the input drive as written is not an exact description of the physical input field (A5 replaces the quan- tized input by its coherent amplitude), and the weak-measurement scheme removes the dynamics of the measurement apparatus and assumes unidirectional propagation to the detector (A6); both are of SLH type and are owned as axioms here. What is not approximated in Eq. (S1), relative to the usual treatment of a delayed-feedback atom within the same SLH framework, is the memory itself: the mirror-mediated feedback is carried exactly by the coherent mode structure, with no Born–Markov truncation of the delayed self-interaction, no pseudomode fitting, and no auxiliary lossy degrees of freedom. This is the sense—and the only sense—in which we call Eq. (S1) exact, and the theorem below is a statement about that generator. How well the generator describes a given physical device is an empirical question about A1–A6, and the bud- gets that govern it (ε flat , coherence over the round trip) are stated in the Discussion of the main text. The input is the continuous-time version of the discrete drive, u(t) = ∞ X k=−∞ u k 1 (k+1)T − t 1 t− T off − kT , T = T on + T off ,(S2) so that u(t) = u k during the on-interval of step k and u(t) = 0 during the off-interval, ensuring measurements are not taken while input is injected. 3 The one approximation internal to the model: adiabatic elimination of the atom. Beyond the axioms A1–A6 themselves, the passage to the working generator involves exactly one further approximation, and it is not a Markov approximation of the envi- ronment: it is the elimination of the atom. Displacing Eq. (S1) by β d (t) = − i g αu(t) cancels the atomic drive in favor of a field drive. When g is large enough that the atom absorbs and re-emits within a window ∆t shorter than the field response time— quantified by the saturation s ∼ (ε INPUT /γ g ) 2 ≪ 1 with γ g = 2g 2 ∆t—the atom adiabatically follows the field and is reduced to a linear transducer. The window is set by the band rather than by the numerics: the flat coupling density of width B has correlation time of order 1/B, and we adopt the convention ∆t ≡ 1/B, so that γ g = 2g 2 /B is the golden-rule rate of the drive channel into the retained band up to the O(1) factor this convention fixes; any other O(1) choice of the window rescales γ g and η by the same factor, and no conclusion depends on it. In the simulations the time-bin width is of order this correlation time at the retained band, so the pipeline’s operational 2g 2 ∆t realizes the physical rate there; no proof constant depends on the discretization. Reversing the displacement gives ̇ρ(t) =L(t)ρ(t) = Γ g · aρ(t),a † + a,ρ(t)a † · Γ † g − iηu(t)[(a † · α) + (α ∗ · a),ρ(t)], (S3) Γ g = Γ + γ g 2 (α⊗ α ∗ ), γ g = 2g 2 ∆t = 2g 2 B , η = γ g 2g = g∆t = g B .(S4) The property used throughout the proof is that Eq. (S3) is Gaussian: it is quadratic in the mode operators with a c-number drive, so its driven steady state is a coher- ent state (Sec. S.2). The corrections to this limit are O(s) and reintroduce the atom’s non-Gaussian (Mollow) statistics; extending the proof to finite s is the open ana- lytic problem stated in the main text. We refer to Eq. (S3) as the linear-transducer (Gaussian) limit and reserve “non-Markovian” for the reduced atomic dynamics of Eq. (S1). S.1.1 One device, one rule: definition of the nested finite-K family Throughout this Supplement, statements are made about a family of finite-dimensional generators Γ (K) g K , one for each mode number K, while the physical claim of the main text concerns a single device. This subsection makes the relation exact: the family is not postulated member by member, but derived, by one rule, from a single set of physical data, and (this is the point of Remark 1 below) no relation between the spectral data of different members is ever assumed. The physical data. The single device is specified by: (i) a coupling density g(ω) carrying the fixed band envelope of Theorem 1, written explicitly as g(ω) 2 = γ 2π Θ B env (ω) with Θ B env (ω) = 1 for |ω − ω 0 | ≤ B env /2 and Θ B env (ω) = exp[−(|ω − ω 0 |− B env /2) 2 /2σ 2 B ] beyond it, normalized so that R Θ = B env + √ 2π σ B : flat over the occupied band, with the Gaussian edge of fixed width σ B used in Lemma 11; (i) a quantization length ℓ, fixing 4 the standing-wave comb ω k with spacing ∆ 0 = πv/ℓ; and (i) the atom–mirror delay τ , fixing the unnormalized emission profile through the mode functions at the atom, ̃α k ∝ g(ω k ) √ ∆ 0 sin(ω k τ/2)—the per-mode coupling ∼ √ ∆ 0 being the continuum- consistent scaling already used after Lemma 3. The outcoupling profile of the readout leakage channel is likewise fixed by the density and the outcoupling geometry, with the same √ ∆ 0 scaling; over the retained band it is frequency-flat, v geo,k = 1/ √ K. Nothing about it is designed, tuned, or shaped; it belongs to the physical data of the device rather than to the design, and in particular v geo,k ̸= 0 for every k, so the readout-side overlap condition holds identically. All spectral freedom of the homodyne detection resides in the local oscillator and the window functions W n (t), which act on the measurement record and never enter the generator. The rule. The member Γ (K) g of the family is obtained from this single object by retaining the K comb modes of the occupied band and repackaging at fixed total rates: writing P K for the coordinate projection onto the retained modes, α (K) = P K ̃α ∥P K ̃α∥ , v (K) geo = P K ̃v ∥P K ̃v∥ , Γ (K) g = iΩ (K) + γ 2 v (K) geo ⊗ v (K)∗ geo + γ g 2 α (K) ⊗ α (K)∗ , (S5) with Ω (K) the retained diagonal block of the comb. Two consequences should be stated plainly. First, Γ (K) g is not the literal compression P K Γ g P K of an infinite-mode opera- tor: the fixed-rate convention (which the trace identity of Lemma 3 requires) rescales the two rank-one dressings by the retained weights ∥P K ̃α∥ −2 and ∥P K ̃v∥ −2 . The dif- ference is controlled by the envelope: as the retained comb fills the fixed envelope’s support, the retained weights converge to one at the Gaussian-tail rate (bounded sym- bolically in Lemma 2 below, and consistent with the simulated comb, on which the deficit falls from O(1) to 4× 10 −4 to below 10 −8 as the retained band passes the en- velope support), so the fixed-rate members converge to the compression of the single physical generator, and no unbounded-operator limit is invoked at any point. Second, each finite K is therefore a different retained-mode description of the one device (a dif- ferent differential equation, exactly as one should expect) generated by one stated rule from one physical object; the operating mode number is set by the band the drive and outcoupling occupy, an operating parameter of the device rather than a fabrication choice. Remark 1 (No spectral nesting). No step of the proof assumes that the eigenvalues or the left/right eigenvectors of Γ (K) g and Γ (K ′ ) g are related for K ̸= K ′ , nested or otherwise, and indeed they are not: the normalization in Eq. (S5) and the rank-two dressing shift all spectral data with K (all spectral data are recomputed independently at each K). The finite-K spectral resolution (Assumption 1) is invoked member by member; the Vandermonde construction of Secs. S.2.3–S.2.4 is self-contained at each K; the strict-scaling corollary of the main text rests on kernel containment (Remark 2), which compares realized input–output kernels and never the spectra of two members; and the many-mode limit enters solely through the scalar kernel convergence h K → h ∞ 5 of Lemma 11, where the family is identified precisely by the fixed coupling-density envelope of the physical data above. The plain reading. Stated without formalism: the mode number is a setting of the machine rather than a piece of it. Nothing is fabricated when K grows; the count is fixed by how many comb lines the occupied band contains, which the quantization length ℓ controls, so changing K is turning a dial on hardware that already exists, in the same sense that a classical universality theorem selects the width of one network family after the accuracy target is announced. The family Γ (K) g is then nothing more than the list of descriptions of that one machine at its different settings, each generated by the single rule above from the single set of physical data. Remark 1 records that these descriptions owe each other nothing spectrally, and that the proof is built so this costs nothing: the only statement that ever compares two settings is the strict-scaling corollary, and it compares what the settings can do (their realized kernels, where containment holds) rather than what they are (their spectra, which reshuffle). The practical consequence is the monotonicity the main text advertises: turning the dial up can never reduce capability, and the constructive selection shows some setting of the dial reaches any accuracy, so a single device becomes more precise by being operated at more modes, with the price paid in geometry and precision budgets rather than in components. The normalisation of the emission profile is needed before the envelope tail can be estimated, and it needs no design hypothesis: it holds for the equidistant comb as such. It is also what the radial floor of Sec. S.2.5.1 compares against. Lemma 1 (Profile normalisation without the delay lock). Let ω k = ω 0 + ∆ 0 k, k = 0,...,K− 1, and θ ≡ ∆ 0 τ . Then D = K−1 X k=0 sin 2 ω k τ 2 = K 2 − 1 2 cos ω 0 τ + (K− 1)θ 2 sin(Kθ/2) sin(θ/2) ,(S6) and consequently D ≥ K/2− 1/(2| sin(∆ 0 τ/2)|). In particular sin ∆ 0 τ 2 ≥ 2 K =⇒ D ≥ K 4 .(S7) Proof. Write sin 2 x = 1 2 (1 − cos 2x), so D = K/2 − 1 2 P k cos(ω k τ ), and sum the geometric series P K−1 k=0 e i(ω 0 τ +kθ) = e iω 0 τ (e iKθ − 1)/(e iθ − 1), whose real part is the displayed Dirichlet kernel. The bound follows from | sin(Kθ/2)| ≤ 1 and the stated implication from 1/(2| sin(θ/2)|)≤ K/4. Lemma 2 (Envelope tail of the retained weights). Let the coupling density carry the envelope of the physical data above, g(ω) 2 = γ 2π Θ B env (ω) with Gaussian edge width σ B , let the retained band be B K = K∆ 0 ≥ B env , and let the comb be equidistant with | sin(∆ 0 τ/2)|≥ 2/K (Lemma 1). Then 1−∥P K ̃α∥ 2 ≤ 4 √ 2π σ B K∆ 0 erfc K∆ 0 − B env 2 √ 2σ B ,(S8) and the same bound with the factor 4 replaced by 1 holds for 1−∥P K ̃v∥ 2 . Under the delay lock D3 in place of Lemma 1, the factor 4 is replaced by q 2 . 6 Proof. For the numerator, ̃α 2 k ≤ g(ω k ) 2 ∆ 0 by sin 2 ≤ 1, and the Gaussian edge is monotone in |ω−ω 0 | beyond B env /2, so each out-of-band term is at most the integral of g 2 over its own ∆ 0 -cell; summing the two sides of the band, X k/∈K ̃α 2 k ≤ γ 2π · 2 Z ∞ B K /2 e −(u−B env /2) 2 /2σ 2 B du = γ 2π σ B √ 2π erfc K∆ 0 − B env 2 √ 2σ B . (S9) For the denominator, P k∈K ̃α 2 k ≥ γ 2π ∆ 0 D with D ≥ K/4 by Eq. (S7). Dividing gives Eq. (S8). For ̃v the mode functions are flat over the retained band, so D is replaced by K exactly and the factor 4 by 1; under D3, D ≥ K/q 2 by Lemma 5(i) and the factor is q 2 . The difference between the fixed-rate member and the compression of the single physical generator is controlled in operator norm by Eq. (S8) through the two rank- one dressings, so that difference vanishes at the Gaussian-tail rate as the retained band passes the envelope support, and no unbounded-operator limit is invoked at any point. Physical selection of the operating modes. The modes entering a given operating point are selected at two logically distinct tiers, and only the second is an approximation. Within the occupied band, the retained modes are never physically decoupled and need not be: when the device supports K in-band modes and the construction uses M ≤ K of them, the selection is performed entirely on the readout side, by the time-integration basis functions of Sec. S.2.3, which are chosen to vanish on the eigenvalue sums of every unused mode (Remark 2). No hardware acts on individual modes; the weight functions W n (t) do. At the band edge, the restriction to the occupied band is physical and approximate: the drive spectrum is confined to the band, and the coupling of out-of-band modes is suppressed by the fixed envelope (they cannot be switched off, only made envelope-small) and the error of exactly this truncation is what the envelope-convolution step of Lemma 11 absorbs, at the kernel level, uniformly in K. The finite-K statements of this Supplement are exact for the retained-band model Eq. (S5); the retained-band model approximates the physical device with the envelope-controlled error just located, and nowhere else. S.2 Proof of Universality We prove that the classC of Theorem 1 of the main text is universal. Figure S.1 maps the argument before it begins: which result feeds which, and which design condition each one consumes. The output is ˆy k = N X n=0 Z T off 0 W n (t) Tr[q n geo e L 0 t ρ ζ (t k )] dt, L 0 ρ = Γ g · aρ,a † + a,ρa † · Γ † g , (S10) with q geo = Re[v ∗ geo ·a] the measured quadrature of the geometric leakage channel and t k = kT . The profile v geo follows from the mode functions themselves. For the standing waves u n (x) = p 2/ℓ sin(k n x) of the terminated line, the atom sits at an interior point 7 Trace identity Lemma 3 no condition Spectral localization Lemma 7 D1 Profile bounds Lemma 5 D3, D1 Dissipativity Lemma 4 no condition Node geometry Lemma 8 D1, D2, D4, D5, D6 Radial floor Lemma 6 D3, D1 Distinctness Lemma 10 D1, D3, D4, D5 Extrapolation identity Proposition 1 no condition Extrapolation tail Lemma 9 D2 via Lemma 8 Error split and Step 3 parameter selection Sec. S.2.5 Theorem 1 main text Supporting results, Secs. S.4–S.7, off the S.2 spine Fading memory Corollary 2 Saturation Lemma 12 Genericity Lemma 14 Shot budget Lemma 15 Supplementary Figure S.1 Map of the universality argument. An arrow runs from a result to the result that consumes it, and the graph is acyclic; the third line of each box names the design conditions the result consumes, or records that it consumes none. Theorem 1 of the main text is closed by the parameter selection of Sec. S.2.5, which consumes the distinctness margins, the extrapolation-tail bound, and the mean decay rate the trace identity pins. Two results carry no design condition and are load bearing for that reason: Lemma 3 is a rank argument on the generator, and Proposition 1 is an algebraic identity about Vandermonde nodes. Four numbered lemmas are absent because no path from them reaches Theorem 1: Lemmas 1 and 2 support the finite-K family of Sec. S.1.1, Lemma 11 supports Corollary 2, and Lemma 13 supports Lemma 12. Sec. S.2.5 states the same dependency order in prose; this figure is that paragraph drawn, and adds no claim to it. x = ℓ a and samples u n (ℓ a )∝ sin(ω n τ/2)—the profiled emission channel α—while the outcoupling port sits at the boundary, where u n vanishes and the leakage is governed by |∂ x u n | 2 = (2/ℓ)k 2 n , identical across modes up to the slowly varying k 2 n . Hence v geo is frequency-flat over the retained band to relative error ε v ≈ B/ω 0 , the same budget as the flat-coupling axiom A2; the two profiles differ because one samples the mode function and the other its boundary derivative. The homodyne local oscillator and the window functions W n (t) carry all spectral freedom of the readout, acting on the measurement record only. Expanded in a Volterra series, ˆy k = N X n=0 X m 1 ,...,m n ≥1 ˆ h n (m 1 ,...,m n )u k−m 1 ·u k−m n . The goal is to show the weights W n (t) can realize any kernel ˆ h n with n≤ N , m i ≤ M ; if M,N grow with the mode number, arbitrarily accurate approximation follows by Theorem 1 of Cuchiero et al. [3], whose hypotheses Sec. S.3.2 verifies directly, so the citation supplies precedent rather than a load-bearing step. The Volterra route also 8 certifies that a reservoir matching kernels up to (N ′ ,M ′ )≥ (N,M ) reproduces every output of a smaller one and can match strictly more. Roadmap and notation This subsection maps the proof. Secs. S.2.2–S.2.4 derive the exact Volterra expansion and reduce kernel prescription to reduction modulo the node polynomial. Sec. S.2.5 is the quantitative heart: it opens with the trace identity, which pins the mean decay rate λ ⋆ = (γ + γ g )/2K the design spends, and continues with the error bound, the designed operating point, the parameter selection, and the scaling ledger. Sec. S.3.2 verifies the density hypotheses. Table S.1 collects the recurring symbols. Each is also defined at first use; the table exists so that no reader must reconstruct a definition from context. S.2.1 Normal ordering and Heisenberg evolution For a single mode with q = Re[a], the normal-ordered and ordinary exponentials are related through Baker–Campbell–Hausdorff by : e qt := e 1 2 a † t e 1 2 at , e qt =: e qt : e 1 8 t 2 ⇒ q n = ⌊n/2⌋ X m=0 n! 8 m (n− 2m)!m! : q n−2m : . Hence any finite P N n=0 W n q n equals P N n=0 f W n : q n : with f W n = P ⌊(N−n)/2⌋ m=0 (n+2m)! 8 m n!m! W n+2m , a triangular and hence invertible relabeling; the highest weights are unchanged and each lower one is introduced once, so the map W n ↔ f W n is a bijection for finite N . Under L 0 in the Heisenberg picture, a 7→ −Γ g · a and a † 7→−a † · Γ † g act independently on each factor of a normal-ordered product, so (a † ) m a n e L 0 t = a † · e −Γ † g t m e −Γ g t · a n , and therefore N X n=0 Z T off 0 W n (t)q n geo e L 0 t dt = N X n=0 Z T off 0 f W n (t) : Re[v ∗ geo · e −Γ g t · a] n : dt.(S11) S.2.2 Coherent driven trajectory and Volterra kernels A note on notation: we write ρ ζ (t) for the driven coherent-trajectory state (the unique attractor of the time-dependent generator under the input Eq. (S2)) and reserve ρ st for the u = 0 stationary state of L 0 . The Gaussian generator Eq. (S3) admits a coherent driven trajectory ρ ζ (t) =|ζ(t)⟩⟨ζ(t)|, where ζ(t) denotes the coherent mode-amplitude vector of the trajectory: substituting this ansatz reduces ̇ρ ζ =L(t)ρ ζ to ̇ ζ(t) =−Γ g · ζ(t)− iηu(t)α, ζ(t) =−iη Z t −∞ e −Γ g (t−τ ′ ) · αu(τ ′ ) dτ ′ .(S12) With the input Eq. (S2), and with t k denoting the end of the k-th on-window (equivalently the start of the k-th measurement window), ζ(t k ) = P n≥1 u k−n ζ n where ζ n =−iη 1− e −Γ g T on Γ g · e −Γ g (n−1)T · α.(S13) 9 SymbolFirst useMeaning KS.1.1number of retained comb modes of the finite-K description M,NS.2 intromemory depth and Volterra order of the matched kernel block TEq. (S10) readout period (one on/off drive-and-measure cycle) γ,γ g S.1outcoupling (readout) and drive-channel rates; ρ = γ g /γ v geo Eq. (S10) readout leakage profile, fixed by the outcoupling geometry, frequency-flat over the retained band (v geo,k = 1/ √ K) αS.1emission profile of the drive channel, α k ∝ sin(ω k τ/2), nor- malized E, e 0 S.2.5Hermitian dressing γ 2 v geo v ∗ geo + γ g 2 α ∗ and its operator norm e 1 , ̄eS.2.5maximal off-diagonal row sum of E (Eq. (S33), K- independent) and ̄e = maxe 0 ,e 1 λ k , v R,k , v L,k S.2.3eigenvalues and right/left eigenvectors of Γ g ω k , δ k , ∆ 0 S.2.5comb frequencies, their detunings from ω 0 , and the comb spacing g 0 , βS.2.5graded-defect amplitude and base, β = 4N + 1 (condition D5, Eq. (D5)) λ ⋆ , ̄rS.2.5mean decay rate (γ + γ g )/2K and node radius e −λ ⋆ T x k S.2.4spectral nodes e −λ k T of the Vandermonde construction ̄ε, ̄ε target S.2.5node-placementerroranditsassignedbudget minε/(4eAM ), 1/4 s α S.2.5peak-to-mean ratio K max k α 2 k of the emission profile A(M,N )Eq. (S26) target-dependent tail prefactor (kernel sup-norms and input bound) S out S.2.5extrapolation tail P m>M ∥c m ∥ 1 h K , h ∞ S.3.1finite-comb and continuum input–output kernels (distinct ob- jects; only h ∞ obeys the delay equation) c 0 , D(c 0 )S.3.1loop-stability margin and the delay-stability condition T P S.3.1revival time 2π/∆ 0 of the finite comb C h , C ψ S.3.1kernel envelope constant and band-envelope exponential mo- ment u max S.2.2uniform input bound defining K u max ηS.1input coupling amplitude of the drive channel, η = γ g /2g c on S.3.2on-window injection factor R T on 0 e −Γ g s ds acting on α, enter- ing ζ 1 Supplementary Table S.1 Recurring notation of the universality proof, with the subsection of first use. Using Eq. (S11) and ⟨ζ|(a † · )(a· )|ζ⟩ = (ζ ∗ · )(ζ· ), the output collapses to ˆy k = N X n=0 Z T off 0 f W n (t) Re[v ∗ geo · e −Γ g t · ζ(t k )] n dt,(S14) 10 so f W n (t) tunes the order-n kernels, with ˆ h n (m 1 ,...,m n ) = Z T off 0 f W n (t) n Y i=1 Re[v ∗ geo · e −Γ g t · ζ m i ] dt.(S15) S.2.3 Selecting eigenmodes by time integration We first record explicitly the spectral assumption used throughout this section. Assumption 1 (Finite-K spectral resolution). At each finite mode number K, the matrix Γ g is diagonalizable, with the biorthogonal spectral resolution Γ g = P K k=1 v R,k λ k v ∗ L,k , v ∗ L,k · v R,k ′ = δ k ′ . We fix the normalisation once, since two conventions are in play. The resolvent construction of Lemma 7 returns ̃v R,k = e k + a k and ̃v L,k = e k + b k with k-th com- ponent unity, so (a k ) k = (b k ) k = 0 and ̃v ∗ L,k ̃v R,k = 1 + b ∗ k a k with |b ∗ k a k | ≤ ∥a k ∥b k ∥: the duality defect is second order in the localisation radius, because both corrections are orthogonal to e k . We therefore set v R,k ≡ ̃v R,k and v L,k ≡ ̃v L,k /(1 + b ∗ k a k ), which is biorthogonal exactly, and every overlap bound stated for the tilde vectors trans- fers with that single bounded factor, which the kernel construction absorbs into the prescribed moment. This holds automatically under the conditions of Theorem 1: the non-resonance condition implies in particular that the eigenvalues are pairwise distinct at each finite K (the two-index case of Lemma 10(a,b) at designed operating points; at finite K it also follows from the disjoint localization disks of Lemma 7), and a finite matrix with distinct eigenvalues is diagonalizable, its left/right eigenvector systems admitting the stated biorthogonal normalization. We emphasize that the resolution is invoked only at finite K, where the sum is finite and no convergence question arises. No infinite- mode spectral decomposition of the (non-normal) generator is used anywhere in the proof: the kernel-matching construction below operates at finite K, and the many- mode limit enters solely through the scalar input–output kernel h K (t) and its uniform decay (Lemma 11). Questions about the convergence of P k v R,k λ k v ∗ L,k as K →∞— equivalently, whether the eigenvector system of the non-normal limit operator forms a Riesz basis—therefore never arise; the proof rests on the kernel precisely because operator-level arguments in the many-mode limit fail (Lemma 3). The family of finite- K generators to which this assumption applies is defined, from a single physical device by a single rule, in Sec. S.1.1; in particular, no relation between the spectral data at different K is assumed anywhere (Remark 1). Expanding Eq. (S15) in this basis, each term’s time integral has the form F (s) = R T off 0 f (t)e −st dt, whose Taylor coefficients about s = 0 are (−1) m R T off 0 f (t)t m dt. Real f enforces F (s) ∗ = F (s ∗ ). What is needed is not the prescription of Taylor coefficients at s = 0 but interpolation at the finitely many prescribed points Λ j : taking f W n in the span of t p p<J with J the number of points, the realized values are Mc with M jp = R T off 0 t p e −Λ j t dt, and prescribability is the nonsingularity of M . As T off → ∞, M jp → p!/Λ p+1 j and detM tends to a Vandermonde in the recip- rocals 1/Λ j , nonzero exactly when the Λ j are distinct and nonzero—supplied by the distinctness of Lemma 10 and by λ k ̸= 0 of Lemma 4. At finite T off the correction 11 is R ∞ T off , exponentially small once Re Λ j T off ≫ J ; nonsingularity is certified outright above that threshold. Below it nonsingularity is not lost but merely non-uniform, and a genericity argument of the same type used for the delay in Sec. S.6 settles it: each en- try R T off 0 t p e −Λ j t dt is entire in the window length, hence so is the determinant, which tends as T off → ∞ to a Vandermonde in the reciprocals 1/Λ j , nonzero whenever the Λ j are distinct and nonzero. An analytic function on (0,∞) that is not identically zero has a discrete zero set, and the union over the finitely many conjugation-resolved point sets at fixed (M,N,K) is again discrete, so prescribability holds at every win- dow length outside a discrete exceptional set. In plain terms: over a short window two nearby exponentials can look alike and the matrix can degenerate, but degeneracy is a coincidence rather than a regime, so the bad windows are isolated points on the clock exactly as the bad mirror distances are isolated points on the ruler. The deposited evaluation at the operated T off confirms that this operating point is not exceptional, and discharges no step of the argument. We state plainly why the threshold is not quantified and the genericity argument is carried instead. Certifying nonsingularity outright would require Re Λ j T off ≳ J with J the number of conjugation-resolved sum points, which grows combinatorially in (M,N ) as J ≤ P n≤N 2 n M +n−1 n , while the floor available from Lemma 4 is Re Λ j ≥ γ/18K; the resulting requirement T off ≳ 18KJ/γ is incompatible with the readout period T ≥ T off that Step 3 selects, which grows only as O(M lnM + lnε −1 ). The threshold route is therefore closed to us, and T off is deliberately not a selected parameter of the theorem: the discrete exceptional set is what carries this step, exactly as it carries the delay in Sec. S.6. Defining F n (Λ) = R T off 0 f W n (t)e −Λt dt, we only need to fix F n at the discrete points P j i=1 λ ∗ k i + P n i=j+1 λ k i . The spectral conditions P k Re[λ k ]n k ̸= 0, P k Im[λ k ]n k ̸= 0 (for P k n k = 0, P k n 2 k > 0) guarantee these points are distinct, so real basis functions F ± k 1 ...k n can be built that equal ±1,±i on the target point and 0 elsewhere while respecting F (s) ∗ = F (s ∗ ). This yields ˆ h n (m 1 ,...,m n ) = X k 1 ,...,k n Re " w Ord(k 1 ...k n ) n Y i=1 v ∗ geo · v R,k i v ∗ L,k i · ζ m i # .(S16) The overlap conditions v ∗ geo · v R,k ̸= 0 and v ∗ L,k · α̸= 0 ensure no factor vanishes. S.2.4 From kernels to weights: Vandermonde inversion By Eq. (S13), v ∗ L,k · ζ m = e −λ k (m−1)T (v ∗ L,k · ζ 1 ); every Reλ k > 0 by Lemma 4, so |e −λ k T on | < 1 and 1− e −λ k T on ̸= 0. Absorbing the k-dependent factors into W, Eq. (S16) becomes ˆ h n (m 1 ,...,m n ) = X k 1 ,...,k n Re " W Ord(k 1 ...k n ) n Y i=1 V k i m i # , V km = x m−1 k , x k = e −λ k T . (S17) V is a Vandermonde matrix in the nodes x k = e −λ k T . Injectivity of x k requires that distinct eigenvalues not be aliased by the exponential, i.e. Im[(λ k − λ k ′ )]T /∈ 2πZ for k ̸= k ′ ; this is exactly the non-resonance condition applied to the two-index tuple n k = 12 +1, n k ′ =−1 (which has P k n k = 0, P k n 2 k = 2 > 0), so distinctness of the complex eigenvalue sums guarantees it. Hence all x k ′ −x k ̸= 0, so detV = Q k<k ′ (x k ′ −x k )̸= 0 and V is invertible when exactly M modes participate. Inverting, W Ord(k 1 ...k n ) = X m ′ 1 ≤·≤m ′ n W ′ m ′ 1 ...m ′ n X P n n Y i=1 V −1 P n (m ′ ) i k i ⇒ ˆ h n (m 1 ,...,m n ) =W ′ m 1 ...m n , so all kernels of order n ≤ N and memory m i ≤ M are fixed by the weights. This requires the device to support at least M modes; no infinite-mode limit is taken, here or anywhere in the proof. Remark 2 (Exact matching persists for K > M ). For K > M we do not resort to a pseudoinverse (which would yield only least-squares matching): the eigenmode- selection step of Sec. S.2.3 allows the weights to be chosen so that the basis functions F vanish on the eigenvalue sums of any K−M unused modes, reducing the matching to a square Vandermonde system on the M selected modes. This step consumes distinctness over the full K-mode sum set—the sums involving unused modes must be resolved from those of the selected block—at the designed point this resolution is supplied by extending the defect hierarchy over every in-band mode (Remark 3) together with Lemma 6, whose prime q is sized against K phys by condition D3. The matching therefore remains exact for every K ≥ M , and the strict-scaling corollary of the main text rests on this remark. S.2.5 The error bound and the designed operating point Dependency order. The conditions and lemmas of this section form an acyclic chain, and it is worth dis- playing once. The physical data of Sec. S.1.1 fix the delay lock D3, which feeds the profile bounds and the radial floor. The rate ratio D4 and the graded defects D5 are then selected, and the comb spacing D1 last; its localization link D1a supports the spectral-localization lemma, which in turn supports both the explicit strict- dissipativity floor and the node-geometry lemma, and the latter with D2 supports the extrapolation-tail lemma. Distinctness draws on D4, D5 and the radial floor, and feeds the selection and Vandermonde steps. Four citations inside this section run forward in the text and not in the logic, and all four are recorded here rather than left for the reader to notice: the explicit floor of Lemma 4 uses the localization lemma stated below it, and Lemmas 4, 5 and 6 each cite an entry of the comb-spacing condition Eq. (D1), which is displayed in full at the localization lemma because that is where its first entry is consumed. The count is enforced by the gate check forwardrefs.py; outside this section the only such citation is Lemma 12’s use of Lemma 13. In every case the cited statement is proved without reference to the citing one, so no cycle ex- ists. Step 3 selects in the order (M,N ) → q → ̄ε target → T → ρ → g 0 → ∆ 0 . The mode-space generator obeys a conservation law that fixes the arithmetic of every rate in the proof, and we establish it before anything is built on top of it. The Hermitian part of Γ g has rank at most two, so its trace is pinned by the two physical rates and cannot grow with the mode number. Two consequences follow immediately: the mean 13 decay rate is exactly λ ⋆ = (γ + γ g )/2K, which is the quantity the designed operating point of this section spends; and no gap Re[λ k ]≥ ε D > 0 can hold uniformly in K, so spectral-gap hypotheses are unavailable to us as a matter of structure rather than of technique. Lemma 3 (Trace identity). Let Γ g = iΩ + γ 2 (v geo ⊗ v ∗ geo ) + γ g 2 (α⊗ α ∗ ) act on C K with ∥v geo ∥ =∥α∥ = 1. Then K X k=1 Re[λ k ] = Tr 1 2 (Γ g + Γ † g ) = γ + γ g 2 ,(S18) independently of K. Consequently the mean of Re[λ k ] is exactly λ ⋆ = (γ +γ g )/2K and min k Re[λ k ] ≤ λ ⋆ , so no uniform gap can hold along the continuum limit: with the Hermitian part of rank at most two, the closing of the mode-space gap is structural. Lemma 4 (Strict dissipativity at finite K). Let Γ g = iΩ + E with Ω = diag(δ k ) real and E = γ 2 (v geo ⊗ v ∗ geo ) + γ g 2 (α⊗ α ∗ ), with γ > 0, γ g ≥ 0, and the geometry- fixed flat coupler v geo,k = 1/ √ K of Sec. S.1.1. If the δ k are pairwise distinct then every eigenvalue of Γ g satisfies Reλ k > 0. If in addition the detunings are separated, |δ k − δ j | ≥ 7 8 ∆ 0 for j ̸= k — which conditions D1 and D5 supply on the designed comb, and which an undefected equidistant comb supplies with separation ∆ 0 — and ∆ 0 ≥ 8 √ K ̄e, then Reλ k ≥ γ 18K for every k.(S19) Consequently λ k ̸= 0 for every k; under the spacing hypothesis every conjugation- resolved eigenvalue sum of order n ≥ 1 obeys Re Λ (S) ≥ γ/18K, while the unconditional half Reλ k > 0—which is all that the resolvent existence of Sec. S.7 requires—holds without it. Proof. The strategy is to read the real part of any eigenvalue off the Rayleigh quotient: the frequency comb contributes nothing to it, so the whole real part comes from the Hermitian dressing, which the trace identity bounds below uniformly in K. Let Γ g x = λx with x̸= 0. Then λ = x ∗ Γ g x/(x ∗ x), and both x ∗ Ωx and x ∗ Ex are real, Ω and E being Hermitian, so Reλ = x ∗ Ex x ∗ x = γ 2 |v ∗ geo x| 2 + γ g 2 |α ∗ x| 2 x ∗ x ≥ 0.(S20) Suppose Reλ = 0. Both terms are nonnegative and γ > 0, so v ∗ geo x = 0 and Ex = 0, whence iΩx = λx and x is an eigenvector of the real diagonal Ω. Its entries being pairwise distinct, every eigenspace is one-dimensional, so x = e k for some k up to scale; but then v ∗ geo x = 1/ √ K ̸= 0, a contradiction. Hence Reλ > 0. The δ k are pairwise distinct on the designed comb:|δ k −δ j |≥ ∆ 0 −2g 0 for k ̸= j, and the second entry of Eq. (D1) with g 0 = ̄ε target /8T max , where T max is the readout-period ceiling fixed in the design step below gives ∆ 0 /g 0 ≥ 512(e 0 T max / ̄ε target ) 2 > 2048, since λ ⋆ T max = Λ max ≥ 2 ln 2M with λ ⋆ ≤ e 0 and ̄ε target ≤ 1. Here λ ⋆ ≤ e 0 because E is positive semidefinite of rank at most two with tr E = (γ + γ g )/2 (both factors being unit vectors), so e 0 =∥E∥≥ 1 2 tr E = (γ + γ g )/4≥ (γ + γ g )/2K = λ ⋆ for K ≥ 2. 14 For the explicit floor, Lemma 7 gives ∥v R,k − e k ∥ ≤ 4e 0 /∆ 0 , which under ∆ 0 ≥ 8 √ K ̄e≥ 8 √ Ke 0 is at most 1/(2 √ K)≤ 1/2. Hence ∥v R,k ∥≤ 3/2 and |v ∗ geo v R,k | ≥ |v geo,k |−∥v geo ∥v R,k − e k ∥ ≥ 1 √ K − 1 2 √ K = 1 2 √ K ,(S21) so Eq. (S20) at x = v R,k , discarding the nonnegative drive term, gives Reλ k ≥ γ 2 · 1 4K · 4 9 = γ/18K. Finally Re Λ (S) = P i Reλ k i ≥ nγ/18K, conjugation leaving real parts unchanged. The identity is immediate from the displayed trace: iΩ is anti-Hermitian and the two rank-one Hermitian terms contribute γ/2 and γ g /2; every Reλ k > 0 at finite K, with the explicit floor γ/18K, by Lemma 4, which uses only the flat coupler and the distinctness of the δ k and is independent of the emission profile. This subsection contains the quantitative heart of the proof, and we begin with a plain statement of its plan. The argument has three steps. 1. From the error to the node geometry of an ideal comb. The distance between target and device output is split into a truncation error of the target itself and the response of the matched device beyond the matched depth—the tail. An algebraic identity shows that once the kernels are matched, the tail is fixed by the spectral nodes x k = e −λ k T alone (no weight norm ever enters), and for an ideal equidistant comb with a common decay rate the tail is computable in closed form and exponentially small: ≤ 2M ̄r M . 2. The obstruction and the designed repair. The same perfect regularity that makes the ideal tail computable defeats the order-n kernel selection (degenerate eigenvalue sums) and zeroes the radial separations the selection also needs. This step repairs both, introducing each design condition at the point where the proof consumes it—graded frequency defects for the angular coincidences, the drive-channel spread for the radial ones, spacing and timing conditions as the localization and node- geometry lemmas require them—and assembles the conditions into a summary table at the end. 3. Parameter selection. With the tail bound and the distinctness margins in hand, the accuracy target ε is met by choosing, once and in order, the depth and order (M,N ), the readout period T , and the comb; every constant is explicit. A ledger of how every design parameter scales with (M,N,ε), the accuracy envelope of a device at fixed geometry, and the distinction between the designed certificate and a generically fabricated device follow at the end. Step 1: from the error split to the ideal comb Splitting the error. We bound the error against the target’s own depth-(M,N ) Boyd–Chua truncation, so that the two integers (M,N ) are chosen jointly and once, and no constant is required to stay fixed while a memory depth is enlarged. Write y (M,N ) for the truncation of the target functional y at Volterra order N and memory depth M , and let ˆy denote the realized device output whose weights match, exactly, every kernel of y (M,N ) (this 15 is the kernel-matching construction of Secs. S.2.3–S.2.4, valid for K ≥ M ). Then |y k − ˆy(u k )| ≤ |y k − y (M,N ) k | | z (i) target truncation + |y (M,N ) k − ˆy(u k )| | z (i) realized tail beyond depth M .(S22) Quantitatively, for a target with a convergent Volterra expansion, term (i) is bounded by an element δ MN of a sequence decreasing in M,N , by the Boyd–Chua approxima- tion theorem for fading-memory functionals [4]; the bound is uniform in the reservoir, and K ≥ M modes are required only so that term (i) can be realized. Term (i) is the total weight of the realized kernels beyond the matched depth—the device’s response does not stop at lag M merely because the weights were fitted there—and the next move shows that this excess is determined by the spectrum alone. The extrapolation identity. The key structural fact: once the kernels are matched at lags ≤ M , the realized kernels at all lags > M are determined by the node setx k alone, independently of the weights that realize the matching. In particular no norm of the matched weights enters the error bound. (Supplementary Fig. S.3 documents why any route through that norm fails: for nodes confined to an arc, ln∥V −1 ∥ ∞ grows linearly in M .) This is the fact that makes the bound finite-dimensional and explicit, and we now prove it. Proposition 1 (Extrapolation identity). Let x 1 ,...,x K ∈ C be distinct, let V be the Vandermonde matrix V kμ = x μ−1 k (μ = 1..K), and for each m≥ 1 let c m ∈ C K be the unique solution of (V c m ) k = x m−1 k for every k = 1,...,K(S23) (so that c m = e m , the m-th standard basis vector, for m ≤ K). Then the feature functions f m (t) = P k g k (t)x m−1 k satisfy, identically in t and for every coefficient family g k (t), f m (t) = K X μ=1 c m,μ f μ (t), m≥ 1.(S24) Consequently any output functional built multilinearly from f m and matched at lags ≤ K has its order-n kernels at arbitrary lags given by the tensor extrapolation of the matched block, ˆ h n (m 1 ,...,m n ) = K X μ 1 ,...,μ n =1 n Y i=1 c m i ,μ i H n (μ 1 ,...,μ n ),(S25) where H n (⃗μ) are the (weight-dependent, complex) matched-block moments; the extrap- olation coefficients c depend on the spectrum only. Proof. The defining equation (V c m ) k = x m−1 k reads, written out, P μ c m,μ x μ−1 k = x m−1 k for every k: the monomial of degree m− 1, evaluated on the node set, coincides with a fixed linear combination of the first K monomials evaluated on the same nodes. (For m ≤ K this combination is trivially the monomial itself, hence c m = e m ; for m > K it is the reduction of a high power modulo the degree-K polynomial vanishing on the nodes.) Multiply this scalar identity by g k (t) and sum over k: every term of f m 16 is reproduced, giving the first display. For the tensor statement, expand the product Q i f m i (t) using the first display factor by factor; multilinearity distributes the sums, and integrating against the (fixed) weights preserves the linear combination, yielding Eq. (S25). In our setting f m (t)= v ∗ geo e −Γ g t ζ m = P k g k (t)x m−1 k with g k (t)= −iη(v ∗ geo v R,k )(v ∗ L,k α)λ −1 k (1− e −λ k T on )e −λ k t and x k = e −λ k T (Eq. (S13)), so Propo- sition 1 applies verbatim. Because the output uses real weights and real parts, the matched block enters through the conjugation-resolved moments H (S) n (⃗μ), one for each conjugation pattern S ⊆ 1..n of the factors f/ ̄ f . The constructive weights of Secs. S.2.3–S.2.4 prescribe these moments directly. Two statements must be kept apart here. Prescribability —that the moments can be set to chosen real values at all— requires the eigenvalue sums to be resolved, which the design of Sec. S.2.5.1 supplies with proven margins and a generic operating point supplies with verified ones. The resulting bound |H (S) n (⃗μ)| ≤ 2 n max ⃗μ |h (M,N ) n (⃗μ)| is then combinatorial, counting the 2 n conjugation patterns across which the real target kernel is distributed (Lemma 10); it carries no design margin, and the size of the weights realising it, which does, never enters the error bound (Sec. S.2.5, resource (iv)). With Eq. (S25) in hand, term (i) reduces to a property of the node set. Write S out ≡ P m>M ∥c m ∥ 1 for the total extrapolation weight beyond the matched depth (note ∥c μ ∥ 1 = 1 for μ ≤ M ). Since y (M,N ) has no kernels beyond depth M while the realized functional carries exactly the extrapolated kernels, |y (M,N ) k − ˆy(u k )| ≤ N X n=1 u n max X ⃗m/∈[1,M ] n | ˆ h n (⃗m)| ≤ h N X n=1 n 2 2n u n max M n−1 ∥h (M,N ) n ∥ ∞ i | z =:A(M,N ) · S out ,(S26) where the second line uses the counting bound (M + S out ) n − M n ≤ nS out (M + S out ) n−1 ≤ nS out (2M ) n−1 for S out ≤ M , together with the 2 n moment bound above. The hypothesis S out ≤ M is not yet available at this point; it is a condition on the node geometry that the designed operating point delivers, and the parameter selection of Step 3 verifies it explicitly (S out ≤ 2/M at every selected operating point). The prefactor A(M,N ) is a property of the target’s own truncation, through its kernel sup- norms and the input bound; it involves no weight norms and no spectral constants. Everything now rides on making S out small, which is a pure question of where the nodes x k sit in the complex plane. For orientation on the size of A(M,N ): when the target’s kernels decay geometrically in order,∥h (M,N ) n ∥ ∞ ≤ C h ρ n h , the sum is bounded by explicit geometric terms, A(M,N ) ≤ (C h /M ) P n≤N n (4u max ρ h M ) n : summable uniformly in N when 4u max ρ h M < 1 and dominated by its last term otherwise. In either case A is explicit in the target data alone; it enters the readout period only logarithmically and the comb spacing linearly (Table S.3). 17 The ideal equidistant comb: the tail is not the problem. Suppose first that the device could be operated so that its spectral nodes sit exactly at x 0 k = ̄rω k , ω = e −2πi/M , k = 0,...,M − 1, K = M :(S27) the M -th roots of unity, shrunk onto a common radius ̄r < 1. (The design below achieves this up to a small placement error; here we take it exactly. Note that the common radius means all decay rates Reλ k coincide, so the radial separations that Lemma 10 will later require are exactly zero here: the ideal comb is an illustrative limit for the tail computation rather than an admissible operating point, and Step 2 repairs precisely this.) The Vandermonde matrix then factors into a discrete Fourier transform and a radial scaling, V 0 kμ = (x 0 k ) μ−1 = ̄r μ−1 ω k(μ−1) = (FD r ) kμ , F kμ = ω k(μ−1) , D r = diag( ̄r μ−1 ), (S28) and since F ∗ = MI, the inverse is explicit: (V 0 ) −1 = D −1 r F ∗ /M . The extrapolation coefficients can now be computed in closed form. The right-hand side of the defining equation has entries (x 0 k ) m−1 = ̄r m−1 ω k(m−1) , and applying F ∗ gives (F ∗ ξ 0 m ) μ = X k ω −k(μ−1) ̄r m−1 ω k(m−1) = ̄r m−1 X k ω k(m−μ) = ̄r m−1 M 1m≡ μ (mod M ),(S29) by the geometric-sum orthogonality of the roots of unity. Dividing by M and undoing the radial scaling, c 0 m,μ = ̄r m−μ 1m≡ μ (mod M ).(S30) Physically, Eq. (S30) says the device’s response at any lag m > M is the response at the lag μ ≡ m (mod M ) inside the matched window, damped by one factor of ̄r for every additional readout period. The angular positions repeat with period M because the nodes are roots of unity; the radius supplies pure geometric decay. Each column of the tail therefore has exactly one nonzero entry, and for m = μ + jM (j ≥ 1, 1≤ μ≤ M ), ∥c 0 m ∥ 1 = ̄r jM =⇒ X m>M ∥c 0 m ∥ 1 = X j≥1 M ̄r jM = M ̄r M 1− ̄r M ≤ 2M ̄r M (S31) for ̄r M ≤ 1/2, a condition the parameter selection of Step 3 guarantees with room to spare: the readout period there enforces ̄r ≤ 1/(2M ), whence ̄r M ≤ (2M ) −M ≤ 1/2. The tail of the ideal comb is exponentially small in M , with the exponent set by the common decay rate through ̄r = e −λ ⋆ T : the tail is not the problem. S.2.5.1 Step 2: the obstruction and the designed repair The ideal comb cannot finish the proof. This step first identifies what must change relative to the ideal configuration of Step 1, then introduces the modifications and locates the error each contributes, and only then derives the error bounds, stating each design condition at the point where a proof consumes it. The conditions are collected in Table S.2 at the end of the step. 18 Two features of the ideal comb must change. First, the kernel matching of Secs. S.2.3–S.2.4 requires more than good node positions. To prescribe an order-n ker- nel, the time-integration step of Sec. S.2.3 must resolve the individual eigenvalue sums Λ (S) ( ⃗ k) = P i∈S λ k i + P i/∈S ̄ λ k i : two kernel entries can be prescribed independently only if their sums are distinct points in the complex plane. For the perfectly equidis- tant comb, δ k = ∆ 0 k, these sums collide massively. The simplest collision already occurs at order n = 2: δ 0 + δ 2 = 0 + 2∆ 0 = ∆ 0 + ∆ 0 = δ 1 + δ 1 ,(S32) so the selection integral cannot distinguish the kernel entry at mode pair (0, 2) from the entry at (1, 1); and every arithmetic-progression relation among the δ k produces another such degeneracy. The perfect regularity that made the tail computable in Step 1 is exactly the regularity that defeats the order-n selection. Second, the common radius of the ideal nodes means all decay rates coincide, so eigenvalue sums that differ only in their mode multisets are radially degenerate as well—the failure already flagged parenthetically in Step 1. The design must break both kinds of coincidence without moving the nodes appreciably, since the node positions carry the tail bound. The design uses two mechanisms, angular and radial, both acting far below the node-placement budget; we fix the recurring symbols first, since every bound below prices its contribution against them. All symbols are fixed once: (M,T ) come from the parameter selection of Step 3; K = M ; λ ⋆ ≡ (γ + γ g )/2K is the mean decay rate that the trace identity (Lemma 3) forces; ̄r ≡ e −λ ⋆ T ; ρ ≡ γ g /γ is the drive-to-readout rate ratio; E = γ 2 (v geo ⊗v ∗ geo )+ γ g 2 (α⊗α ∗ ) is the Hermitian dressing with norm e 0 ≡∥E∥≤ (γ +γ g )/2; s α ≡ K max k α 2 k measures the peak-to-mean ratio of the emission profile (s α ≤ q 2 un- der condition D3, Lemma 5(i)); ̄ε target is the node-placement budget handed down from Step 3; and β ≡ 4N + 1. Two norms of the dressing must be kept apart, because the localization argument is a row-wise statement while the perturbation bounds are operator-norm statements. Alongside the operator norm e 0 = ∥E∥ we therefore record the maximal off-diagonal row sum e 1 ≡ max k X j̸=k |E kj | ≤ γ 2 |v geo,k |∥v geo ∥ 1 + γ g 2 |α k |∥α∥ 1 ≤ γ + γ g √ s α 2 ,(S33) where the middle inequality is general while the final bound is specific to the flat coupler,|v geo,k | = 1/ √ K with∥v geo ∥ 1 = √ K (the geometry-fixed profile of Sec. S.1.1), together with |α k | ≤ p s α /K and ∥α∥ 1 ≤ √ K∥α∥ 2 = √ K by Cauchy–Schwarz. The point of Eq. (S33) is that e 1 is bounded independently of K: this is a consequence of the rank-two structure of E with normalized factors, and it is not implied by the operator norm (a general matrix of norm e 0 has row sums as large as √ K− 1e 0 ). We write ̄e ≡ maxe 0 ,e 1 ≤ (γ + γ g max1, √ s α )/2 for the constant that the spacing condition D1, stated below at its point of use, controls. We caution against reading Eq. (S33) as an ordering of the norms themselves: it orders the two bounds. Although s α ≥ 1 always (normalization forces max k α 2 k ≥ 1/K), the flat profile α 2 k = 1/K has e 0 = (γ +γ g )/2 exactly while e 1 = (γ +γ g )(1−1/K)/2 < e 0 , and the verified operating 19 points likewise have e 0 > e 1 (numerical paragraph below); random dressings realize either ordering, so the maximum is removable in neither direction. Angular mechanism (designed): condition D5. Add to each comb frequency an exponentially graded defect, δ k = ∆ 0 k + g 0 β −k , β ≡ 4N + 1, g 0 = ̄ε target 8T max ,(D5) the amplitude g 0 being fixed by the requirement, priced below in Lemma 8, that the largest defect displace any node phase by at most g 0 T ≤ ̄ε target /8, one quarter of the placement budget’s half. Note that g 0 is an absolute frequency scale, independent of ∆ 0 : enlarging the spacing (condition D1 below) leaves the defect hierarchy, and hence the separations it guarantees, untouched. The defects act like digits of a number written in base β: an eigenvalue sum at order n≤ N picks up the defect combination g 0 P k n k β −k , where the integer n k (the signed multiplicity of mode k in the sum) obeys |n k | ≤ n ≤ N . Two different signed multiplicity vectors differ, at their first differing digit j, by at least β −j minus the largest possible carry from all deeper digits; because β = 4N + 1 exceeds four times the largest digit, the carries can never add up to a full digit, and the two combinations stay separated by at least g 0 β −j /2. This is the standard uniqueness of signed-digit representations, and it converts one continuous knob g 0 into a separation guarantee for every one of the finitely many coincidences at once. Radial mechanism (generic): conditions D3 and D4. Sums that share the same signed multiplicity vector but involve different mode multisets have identical comb and defect contributions; they are separated instead by the diagonal decay rates E k = γ/2K + (γ g /2)α 2 k , whose mode dependence enters through the emission profile α k . This spread is supplied by the drive channel, and two conditions make it usable. The profile is not free: the standing-wave geometry of Sec. S.1.1 fixes ̃α k ∝ sin(ω k τ/2) with τ = 2ℓ a /v the atom–mirror delay, ω k = (n 0 + k)∆ 0 + g 0 β −k and ∆ 0 = πv/ℓ, so that ∆ 0 τ/2π = ℓ a /ℓ. Two integers control it: the atom’s position as a fraction of the quantization length, and the band placement n 0 . Condition D3 (delay lock) fixes both. Let q be the smallest prime with q ≥ 2K phys + 1, where K phys ≥ K is the number of in-band modes; by Bertrand’s postulate q ≤ 4K phys + 2. Then ℓ a ℓ = p q ,gcd(p,q) = 1,(D3) n 0 mod q = ν with (2ν mod q)∈1,...,q− 2K phys + 1. The residue set is nonempty since q ≥ 2K phys + 1, and an admissible ν exists since q is an odd prime, so ν 7→ 2ν is a bijection modulo q; p = 1 is always admissible. Condition D3 places the atom at a rational point of the standing wave and shifts the retained band by at most q modes. It replaces a requirement that τ merely avoid a discrete set, which left the radial margin without a symbolic lower bound; Eq. (D3) supplies one (Lemmas 5 and 6). Lemma 5 (Profile bounds under the delay lock). Assume Eq. (D3) and ψ 0 ≡ g 0 τ/2≤ 1/q, which is entry (iv) of Eq. (D1). Write A k = π(n 0 +k)p/q and ψ k = g 0 β −k τ/2, so ̃α k = sin(A k +ψ k ) and α k = ̃α k / √ D with D = P j ̃α 2 j . Then for every 0≤ k ≤ K− 1: 20 (i) (n 0 + k)p ̸≡ 0 (mod q), hence | sinA k | ≥ sin(π/q) ≥ 2/q; (i) 1/q ≤ | ̃α k | ≤ 1 and K/q 2 ≤ D ≤ K; (i) |α k |≥ 1/(q √ K) and s α ≤ q 2 . Proof. (i) q is prime and gcd(p,q) = 1, so (n 0 + k)p ≡ 0 would force ν ≡ −k, hence 2ν ≡ −2k (mod q) with 0 ≤ 2k ≤ 2K − 2. The residues −2k mod q lie in 0∪q− 2K phys + 2,...,q− 1, disjoint from the admissible set of Eq. (D3). Hence | sinA k | = | sin(πm k /q)| with m k ∈ 1,...,q− 1, so | sinA k | ≥ sin(π/q) ≥ 2/q by concavity of sin on [0,π/2]. (i)|ψ k |≤ ψ 0 ≤ 1/q and| sin(A k +ψ k )− sinA k |≤|ψ k | give| ̃α k |≥ 1/q; summing over K modes bounds D. (i) Immediate from (i). Lemma 6 (Radial separation floor). Assume Eq. (D3) and ψ 0 ≤ (4N ) −(q−2) /(32N ), entry (iv) of Eq. (D1). Let d∈ Z K satisfy d̸= 0, P k d k = 0, P k |d k |≤ 2N . Then X k d k α 2 k ≥ 1 8K (4N ) −(q−2) .(S34) Proof. The strategy is to turn the radial separation into a statement about a trigono- metric sum over the comb, and then to bound that sum below by an algebraic argument at the q-th roots of unity rather than by estimating term by term. Since P k d k = 0, P k d k ̃α 2 k =− 1 2 P k d k cos(2A k + 2ψ k ), and expanding the cosine, X k d k ̃α 2 k =− 1 2 X k d k cos 2A k + 1 2 X k d k cos 2A k (1− cos 2ψ k )(S35) + 1 2 X k d k sin 2A k sin 2ψ k , the last two terms bounded by 2Nψ 2 0 and 2Nψ 0 using P k |d k |≤ 2N . Let ζ = e 2πip/q , a primitive q-th root of unity, ν = n 0 mod q, and P d (x) = P K−1 k=0 d k x k . Then P k d k cos 2A k = Re[ζ ν P d (ζ)] = 1 2 R(ζ) with R(x) ≡ x ν P d (x) + x q−ν P d (x q−1 ) reduced modulo x q − 1, an integer polynomial of degree < q. Its two summands contribute d k at exponents (ν + k) mod q and (−ν − k) mod q respec- tively; these sets are disjoint, since ν + k ≡ −ν − k ′ would give 2ν ≡ −(k + k ′ ) with 0≤ k +k ′ ≤ 2K− 2, excluded by Eq. (D3) exactly as in Lemma 5(i). No cancellation occurs, so R ̸= 0, R has at most 2K nonzero coefficients, and its coefficient sum is 2 P k |d k |≤ 4N . q being prime, Φ q = 1 + x + · + x q−1 is irreducible over Q (Eisenstein after x 7→ x + 1) and is the minimal polynomial of ζ. If R(ζ) = 0 then Φ q | R, and degR ≤ q − 1 forces R = c Φ q ; but R has at most 2K ≤ q − 1 nonzero coefficients while every coefficient of Φ q is 1, so c = 0 and R = 0, a contradiction. Hence R(ζ) is a nonzero algebraic integer of Q(ζ), of degree q− 1 over Q, with conjugates R(ζ j ) each bounded in modulus by the coefficient sum 4N . Its norm being a nonzero rational integer, 1≤ q−1 Y j=1 R(ζ j ) ≤|R(ζ)|(4N ) q−2 ,(S36) so |R(ζ)| ≥ (4N ) −(q−2) and the leading term of Eq. (S35) has modulus at least 1 4 (4N ) −(q−2) . With ψ 0 ≤ (4N ) −(q−2) /(32N ) ≤ 1 the remainders total at most 4Nψ 0 ≤ 1 8 (4N ) −(q−2) , so | P k d k ̃α 2 k | ≥ 1 8 (4N ) −(q−2) . Dividing by D ≤ K gives 21 Eq. (S34). This proves the lemma: the radial separation of the nodes is bounded be- low by a quantity fixed by the comb arithmetic alone, so the floor survives at every K and does not degrade as modes are added. Condition D4 (weak drive channel): the rate ratio obeys ρ = γ g γ ≤ min n ̄ε target 8s α Λ max , 1 8Ns α o ,Λ≡ λ ⋆ T = max n ln 2M, 1 M ln 8MA ε o , (D4) so that the very spread being used for radial separation displaces the node moduli by less than one quarter of the budget’s half (priced in Lemma 8): the drive channel is kept strong enough to separate and weak enough not to misplace. The spread is required only to exceed a perturbative remainder that the comb spacing (D1’s third entry, below) makes arbitrarily small. With the two mechanisms in place, the remaining conditions are on the comb spacing (D1) and the readout timing (D2), each stated immediately before the lemma that consumes it. The conditions divide three ways. D2 and the choices of T , N , M , the local-oscillator spectrum and the window functions are measurement settings. D1, D3, D5 and the band placement are geometric settings: dials of one fixed machine. The output coupler is neither: its profile and its flatness ε v are fabricated data. Remark 3 (Embedding the design in a device with K phys > M in-band modes). If the physical operating point retains K phys > M comb modes in the occupied band (Sec. S.1.1), the kernel matching itself is unaffected: the readout-side selection of Re- mark 2 confines every realized kernel, at all orders and all lags, to the M selected modes, so the extrapolation identity and Lemma 9(b) apply verbatim to the selected block. What changes is the rate bookkeeping. The trace budget of Lemma 3 is now shared among K phys modes, so λ ⋆ = (γ + γ g )/2K phys , and the advertised decay ̄r M = e −Mλ ⋆ T carries the factor M/K phys in the exponent relative to the K = M statement e −(γ+γ g )T/2 . The designed operating point therefore chooses the quantiza- tion length ℓ and the occupied band so that approximately M comb modes lie in-band, consistent with the main text’s statement that modes are added by geometry; when K phys > M is operated instead, the rate constant is the one just displayed, reported as such. One further bookkeeping item completes the embedding: the readout-side selec- tion of Remark 2 must resolve the eigenvalue sums involving the K phys −M unselected modes, and at the designed point this resolution is supplied by extending the graded hierarchy of D5 over every in-band mode, δ k = ∆ 0 k +g 0 β −k for k = 0,...,K phys − 1. The signed-digit separation argument of Lemma 10 is agnostic to the digit count, so it applies verbatim: the guaranteed angular floor becomes g 0 β −(K phys −1) /4, and the defect-protection entry of D1 carries the exponent K phys − 1 in place of M − 1, with the scaling ledger of Table S.3 read accordingly. We now derive the error bounds, in the order the construction forces: that the spectrum sits where the design intends (spectral localization), that the nodes land within budget of the roots-of-unity configuration (node geometry), that the tail bound of Step 1 survives the small displacements (extrapolation tail), and that the eigenvalue sums are pairwise resolved (design distinctness). 22 The localization lemma is the first consumer of a condition on the comb spacing, and we state that condition here, in full, since its five entries are consumed at five identified points. Condition D1 (equidistant comb, weak dressing, defect protection): ω k = ω 0 + ∆ 0 k, k = 0..M−1, with ∆ 0 ≥ max n 8 ̄e (D1a), 8 √ K ̄e (D1b), 8q √ K ̄e (D1c), 64e 2 0 T max ̄ε target ,(D1) 384N e 2 0 T β M−1 ̄ε target , 32πNg 0 (4N ) q−2 , 512N K e 2 0 γ g (4N ) q−2 o . The first entry keeps the dressing perturbative and is consumed twice. Its factor 8 is what makes the localization radius satisfy 8e 2 0 /∆ 0 ≤ ∆ 0 /8, an estimate the resolvent chain of Lemma 7 uses twice; its factor q √ K is what makes both overlap conditions hold with an explicit margin. The two sides are not symmetric, and the asymmetry sets the entry. Writing∥v R,k −e k ∥,∥v L,k −e k ∥≤ 4e 0 /∆ 0 (Lemma 7) and∥v geo ∥ =∥α∥ = 1, |v ∗ geo · v R,k | ≥ |v geo,k |− 4e 0 ∆ 0 = 1 √ K − 4e 0 ∆ 0 ,(S37) |v ∗ L,k · α| ≥ |α k |− 4e 0 ∆ 0 ≥ 1 q √ K − 4e 0 ∆ 0 , the drive-side floor|α k |≥ 1/(q √ K) coming from Lemma 5(i). The readout side, with its floor 1/ √ K, would close already at ∆ 0 ≥ 8 √ K ̄e, which leaves 1/(2 √ K); the drive side, whose floor is smaller by the factor q, requires ∆ 0 ≥ 8q √ K e 0 for the remainder to consume at most half of it, and then leaves 1/(2q √ K). The entry is stated at ̄e≥ e 0 so that one expression serves both. The second entry keeps the dressing’s node displacement inside the placement budget and is consumed by Lemma 8. The third entry protects the deepest graded separation against the dressing shift and is consumed by Lemma 10; which entry binds is determined once, in the reach paragraph below, and it is the fifth. The third entry is the price of the constructive certificate and grows exponentially with M (see the resource paragraph below). Because the exponential entries exceed the first by factors growing exponentially in M , the precise constant in the first entry—including its factor q ≤ 4K phys + 2—is immaterial at every realized operating point: strengthening it to secure the drive-side margin cannot change which entry binds. Lemma 7 (Spectral localization). Assume the detunings are separated,|δ k −δ j |≥ 7 8 ∆ 0 for j ̸= k, and that the dressing is weak, ∆ 0 ≥ 8 ̄e. Conditions D1 and D5 together imply the first and D1 the second, and at an undefected equidistant comb the first holds with separation ∆ 0 exactly, so the lemma applies at designed and generic operating points alike. Then the eigenvalues of Γ g = iΩ + E satisfy, for each k, λ k − iδ k − E k ≤ 8e 2 0 ∆ 0 ,(S38) each disk containing exactly one eigenvalue; moreover Reλ k ≥ 0 always, both eigen- vector systems satisfy ∥v R,k − e k ∥ ≤ 4e 0 /∆ 0 and ∥v L,k − e k ∥ ≤ 4e 0 /∆ 0 , and the diagonal decay rates are E k = γ/2K + (γ g /2)α 2 k , with mean value exactly λ ⋆ . 23 The proof introduces seven quantities local to it: the intended disk centre c k , the disk D k and its Schur function s k , the dressing blocks E k and E jk , the k-th column F k of E, and the comb part Ω k . None appears outside this proof. Proof. Throughout, write c k ≡ iδ k +E k for the intended centre of the k-th disk, and note first that the comb frequencies are separated even after the graded defects are switched on: for j ̸= k, |δ k − δ j | ≥ ∆ 0 − 2g 0 ≥ 15 16 ∆ 0 , since g 0 = ̄ε target /8T max (D5) while the second entry of Eq. (D1) forces ∆ 0 ≥ 64e 2 0 T max / ̄ε target , whence ∆ 0 /g 0 ≥ 512(e 0 T max / ̄ε target ) 2 ≥ 2048 ≥ 32; the last step uses e 0 T max ≥ λ ⋆ T max = Λ max ≥ 2 ln 2M ≥ 2 ln 4 and ̄ε target ≤ 1, valid at every operating point of Step 3. We suppress the 15 16 below and use the weaker |δ k − δ j |≥ 7 8 ∆ 0 , which is all the chain needs. First-order picture. That each eigenvalue lies near its comb frequency is, at first order, Gershgorin’s theorem, and the relevant radius is the off-diagonal row sum rather than the largest entry. By Eq. (S33) that row sum is at most e 1 , uniformly in K; since the centres c k are separated by at least 7 8 ∆ 0 ≥ 7 ̄e≥ 7e 1 under D1, the K Gershgorin disks of radius e 1 are pairwise disjoint, and by the standard homotopy argument (deform the off-diagonal part to zero; eigenvalues move continuously and cannot leave a disjoint disk) each contains exactly one eigenvalue. The K disks therefore account for the whole spectrum. We stress that the K-independence of e 1 is supplied by the rank-two structure of E with normalized factors and not by∥E∥ alone: bounding each entry by e 0 would give a row sum as large as √ K− 1e 0 , which is useless here. This picture gives radius e 1 ; the design needs the sharper radius 8e 2 0 /∆ 0 , and the rest of the proof supplies it. Reλ ≥ 0. The numerical range of Γ g has real part W (E) ⊆ [0,e 0 ] and contains the spectrum. Sharper localization. Fix k, set κ≡ 8e 2 0 /∆ 0 , and let ̄ D k =z :|z−c k |≤ κ be the closed disk. By the first entry of Eq. (D1), e 0 ≤ ̄e≤ ∆ 0 /8, hence κ = 8e 2 0 ∆ 0 ≤ 8e 0 ∆ 0 · ∆ 0 8 = e 0 ≤ ∆ 0 8 ,(S39) which is the estimate the chain below uses twice. (This is where the constant 8 rather than 4 in D1 is needed: with ∆ 0 ≥ 4e 0 one obtains only κ≤ ∆ 0 /2, and the separation estimate that follows fails.) Let ˆ A k denote the (K−1)-dimensional block of Γ g with the k-th row and column removed, i ˆ Ω k its diagonal comb part, ˆ E k its dressing part (∥ ˆ E k ∥ ≤ e 0 ), and F k the k-th column of E with its diagonal entry removed (∥F k ∥ ≤ e 0 ). We first establish invertibility of ˆ A k − z on all of ̄ D k , since this is what licenses the Schur complement, and only then form it. For z ∈ ̄ D k and j ̸= k, |iδ j − z| ≥ |δ k − δ j |−|E k |−|z− c k | ≥ 7 8 ∆ 0 − e 0 − κ ≥ 5 8 ∆ 0 ,(S40) since each subtrahend is at most ∆ 0 /8, by |E k | ≤ e 0 ≤ ∆ 0 /8 and Eq. (S39). Hence i ˆ Ω k − z is invertible with ∥(i ˆ Ω k − z) −1 ∥≤ 8/(5∆ 0 )≤ 2/∆ 0 , and ∥(i ˆ Ω k − z) −1 ˆ E k ∥ ≤ 2e 0 ∆ 0 ≤ 1 4 < 1.(S41) 24 Writing ˆ A k − z = (i ˆ Ω k − z) I + (i ˆ Ω k − z) −1 ˆ E k , the second factor is invertible by Neumann series for every z ∈ ̄ D k , so ˆ A k − z is invertible on the closed disk with ∥( ˆ A k − z) −1 ∥ ≤ 2/∆ 0 1− 1/4 = 8 3∆ 0 ≤ 4 ∆ 0 .(S42) This is the point at issue: the resolvent bound is established on the whole closed disk and not merely on its boundary, which is what both the Schur complement and the holomorphy hypothesis of Rouch ́e’s theorem require. With invertibility in hand, the Schur complement on the k-th coordinate says that z ∈ ̄ D k is an eigenvalue of Γ g if and only if s k (z) ≡ c k − z + F ∗ k ( ˆ A k − z) −1 F k = 0,(S43) and s k is holomorphic on ̄ D k because the resolvent is. The coupling correction obeys F ∗ k ( ˆ A k − z) −1 F k ≤ 4e 2 0 ∆ 0 = κ 2 < κ,(S44) so on the boundary circle |z−c k | = κ we have |s k (z)− (c k −z)|≤ κ/2 < κ =|c k −z|, with a factor-of-two margin. By Rouch ́e’s theorem s k has exactly as many zeros in the open disk as the linear function c k −z, namely one. Since the K disks ̄ D k are pairwise disjoint (their radii κ≤ ∆ 0 /8 against centre separation ≥ 7 8 ∆ 0 ) and each carries one eigenvalue, they account for the full spectrum, which is the “exactly one” claim of the statement. The right-eigenvector bound follows from v R,k ∝ e k − ( ˆ A k − λ k ) −1 F k with the same resolvent bound, λ k ∈ ̄ D k being covered by it. The left eigenvectors are the right eigenvectors of Γ † g = −iΩ + E, since E is Hermitian. That matrix has the same diagonal separation |−δ k + δ j | =|δ k − δ j |, the same dressing operator norm e 0 , and the same off-diagonal row-sum bound e 1 —Eq. (S33) is invariant under conjugate transpose, P j̸=k |E jk | = P j̸=k |E kj | for Hermitian E—so every step of this proof applies to it verbatim, giving ∥v L,k − e k ∥≤ 4e 0 /∆ 0 . Finally E k = γ 2 v 2 geo,k + γ g 2 α 2 k = γ/2K + (γ g /2)α 2 k by the geometry-fixed coupler of Sec. S.1.1, and averaging over k using P k α 2 k = 1 gives λ ⋆ . The node-geometry lemma is where the readout timing enters, and it consumes two conditions that we state here in displayed form. Condition D2 (readout-timing condition): ∆ 0 T = 2π(nM + 1) M , n∈ N,(D2) so that δ k T ≡ 2πk/M (mod 2π) for the undressed comb: the node angles are exactly the M -th roots of unity, with the comb wrapping the circle n times between readouts. Because consecutive grid values of T differ by 2π/∆ 0 , the grid always contains a point within one revival time above any required minimum, so this is a timing condition on the readout, satisfiable at fixed comb spacing, and it restricts no device parameter. The fourth and last contribution to the node budget is the readout coupler’s own flatness. The geometry-fixed profile v geo,k = 1/ √ K of Sec. S.1.1 is flat only to relative error ε v (Sec. S.2), so the diagonal decay rate carries a mode-dependent part γ 2K ν k with|ν k |≤ ε v , displacing each node modulus by λ ⋆ Tε v /(1 +ρ). Condition D6 (coupler 25 flatness) caps that displacement at the same quarter of the budget’s half that D1, D4 and D5 cap the other three: ε v ≤ ̄ε target 8 Λ max ,Λ max ≡ λ ⋆ T max .(D6) Eq. (D6) is Eq. (D4)’s first entry with s α replaced by unity and ρ by ε v : the two channels enter the node budget through the same product of a spread with the di- mensionless readout period, and are capped identically. Unlike D1–D5, D6 is not a setting. Its quantity ε v ≈ B/ω 0 is fixed once the band and the transition frequency are, so Eq. (D6) is a feasibility test evaluated at the operating point rather than a dial turned to meet it. Lemma 8 (Node geometry). Let M ≥ 2. Under D1–D5 the nodes x k = e −λ k T satisfy x k = ̄rω k (1 + ε k ) with ω = e −2πi/M , ̄r = e −λ ⋆ T (up to one global phase absorbed into the mode labeling), and |ε k | ≤ ̄ε ≡ 2 ρs α λ ⋆ T/(1 + ρ) |z drive-channel spread + 8e 2 0 T ∆ 0 |z dressing shift (S45) + g 0 T |z graded defects + λ ⋆ T ε v /(1 + ρ) |z coupler flatness ≤ ̄ε target , provided the bracket is ≤ 1/2; the final inequality holds because D4, D1, D5 and D6 cap the four contributions at ̄ε target /8 each. The fourth term is the readout coupler’s own spread: the outcoupling profile is flat only to relative error ε v (Sec. S.2), so E k inherits a mode-dependent part γ 2K ν k with|ν k |≤ ε v , and D6 caps its node displacement alongside the other three. The nodes are pairwise distinct for ̄ε < sin(π/M ). Proof. The strategy is to factor the node into the four independent contributions the design controls, bound each against the budget entry that governs it, and collect them. Write each node as a product of the four factors the design controls: x k = e −λ k T = e −i∆ 0 kT |z comb angle e −ig 0 β −k T |z defect phase e −E k T |z diagonal decay e −R k T |z dressing remainder ,(S46) with |R k | ≤ 8e 2 0 /∆ 0 by Lemma 7. Take the factors in turn. The comb angle is ex- act: by D2, ∆ 0 kT ≡ 2πk/M (mod 2π), so e −i∆ 0 kT = ω k . The defect phase deviates from unity by at most g 0 β −k T ≤ g 0 T . The diagonal decay is e −E k T = ̄re −(E k −λ ⋆ )T , and the exponent obeys |E k −λ ⋆ |T ≤ (γ g /2) max k |α 2 k − 1/K|T ≤ (γ g /2)(s α /K)T = ρs α λ ⋆ T/(1 +ρ), which condition D4 caps at ̄ε target /8 through the readout-period ceil- ing T max . The dressing remainder deviates by at most |R k |T ≤ 8e 2 0 T/∆ 0 . Combining the three small factors with the elementary bound |e w − 1|≤ 2|w| for |w|≤ 1/2 gives Eq. (S45). Distinctness: adjacent ideal nodes are separated by 2 ̄r sin(π/M ), which exceeds the total displacement 2 ̄r ̄ε. The tail beyond the matched depth is controlled at any operating point, and con- trolled far more sharply at the designed one. Both bounds are the same computation — reduction modulo the node polynomial — and differ only in the estimate available for it, so we state them together. 26 Lemma 9 (Extrapolation tail, generic and designed). Let x 1 ,...,x M ∈ C be the nodes and let c m be the extrapolation coefficients of Proposition 1. As there, the nodes are pairwise distinct; this is what makes c m well defined, since Proposition 1 fixes c m by a Vandermonde solve. Write S out for the total realized kernel weight beyond the matched depth, defined in Eq. (S49) below as a sum of moduli, where ∥·∥ 1 denotes the sum of moduli of the entries of a coefficient vector; and let r ⋆ ≡ max k |x k |(S47) denote the spectral radius of the node set. Then: (a) (Generic operating point.) If r ⋆ < 1/(M + 1), then S out ≤ M (1 + r ⋆ ) M X j≥1 r j ⋆ (M + j) M +j j j M M ≤ eM (M + 1)r ⋆ 1 + O(Mr ⋆ ) . (S48) In particular, with x k = e −λ k T and T the readout period, S out decays as e − min k Reλ k T . (b) (Designed operating point.) If in addition the nodes take the designed form x k = ̄rω k (1 + ε k ) with ω = e −2πi/M the primitive M -th root of unity, ̄r the common radius, and |ε k |≤ ̄ε≤ 1/4, and if ̄r ≤ 1/(2M ) and M ≥ 2, then S out = X m>M ∥c m ∥ 1 ≤ 2M ̄r M + eM ̄ε.(S49) The hypotheses nest: those of (b) imply that of (a), since r ⋆ ≤ ̄r(1 + ̄ε) ≤ 5 8M < 1/(M + 1) for M ≥ 2, the last step being 5(M + 1) < 8M ⇐⇒ M > 5/3. The conclusions do not nest, and (b) is not a corollary of (a): Eq. (S49) decays exponentially in M where Eq. (S48) decays only linearly in r ⋆ , and the two are reached by different arguments. The proof introduces seven quantities local to it: the reduction operatorR and the smear operator E , the polynomial coefficients a i and b i , the standard basis vectors e i and e j , and the reduced polynomial p m . None appears outside this proof. Proof. Both parts compute the same object. By Proposition 1, c m is the coefficient vector of x m−1 reduced modulo the node polynomial Π(x) = Q k (x−x k ), so S out asks how fast repeated reduction shrinks a monomial. Part (a) estimates this by a contour integral, which needs nothing of the node configuration beyond its radius. Part (b) exploits the fact that at the designed point Π is a perturbed binomial, for which reduction is a single explicit instruction. Throughout, ∥·∥ 1 on a polynomial means the sum of moduli of its coefficients. (a) Generic. Write P (x)≡ Π(x) = P M i=0 a i x i and let p m (x) = x m−1 mod P (x), so that c m is the coefficient vector of p m ; equivalently p m interpolates the entire function z 7→ z m−1 at the nodes. No inverse of the Vandermonde matrix is formed. The idea is that interpolation at nodes confined to a small disk cannot produce large coefficients, 27 and a contour integral is the cleanest way to say so. By the Hermite formula, for any circle Γ =|z| = R enclosing every node, p m (x) = 1 2πi I Γ z m−1 P (z)− P (x) (z− x)P (z) dz.(S50) Expanding P (z)−P (x) z−x = P M−1 μ=0 x μ P i>μ a i z i−1−μ and taking R < 1, so that |z| i−1−μ ≤ 1 for every i > μ, its coefficient ℓ 1 norm in x is at most M P i |a i | ≤ M (1 + r ⋆ ) M , using |a i | ≤ M i r M−i ⋆ . On Γ the denominator is bounded below, |P (z)|≥ Q k (R−|x k |)≥ (R− r ⋆ ) M . Hence for every R∈ (r ⋆ , 1), ∥c m ∥ 1 ≤ M (1 + r ⋆ ) M R m (R− r ⋆ ) M .(S51) The contour radius is a free parameter here, its freedom being exactly the statement that the estimate holds for every R ∈ (r ⋆ , 1), so we may optimise it per term rather than once. Minimising R M +j /(R− r ⋆ ) M over R gives R ⋆ = (M + j)r ⋆ /j, admissible (R ⋆ < 1) for every j ≥ 1 when (M + 1)r ⋆ < 1, and substituting yields the summand of Eq. (S48). The j = 1 term is M (1+r ⋆ ) M r ⋆ (M +1) M +1 /M M ≤ eM (M +1)r ⋆ (1+r ⋆ ) M , and the remaining terms are smaller by factors O(Mr ⋆ ). (b) Designed. Everything follows from what the node polynomial is at the designed point. The M numbers ̄rω k are precisely the roots of x M − ̄r M , which is monic of degree M , so at ̄ε = 0 the node polynomial is that binomial and reduction is a single instruction: replace x M by ̄r M . Away from ̄ε = 0 write Π(x) = x M − ̄r M + δ(x),degδ ≤ M − 1,(S52) so that δ measures, and is the only thing that measures, how far the design pushes the node polynomial off that binomial. Step 1: the size of δ. The coefficients of δ are the elementary symmetric functions e j (x k ), and these vanish identically at ̄ε = 0 for 1 ≤ j ≤ M − 1: that vanishing is the whole content of “roots of unity at a common radius”. Using| Q k∈S (1 +ε k )− 1|≤ (1 + ̄ε) j − 1≤ j ̄ε(1 + ̄ε) j−1 on each of the M j subsets of size j, and M j j = M M−1 j−1 , D ≡ ∥δ∥ 1 ≤ M ̄ε ̄r 1 + ̄r(1 + ̄ε) M−1 ≤ e 2 ̄ε,(S53) the last step using ̄r(1 + ̄ε)≤ 1/M and M ̄r ≤ 1/2, so that (1 + 1/M ) M−1 < e. The sum is dominated by its first term M ̄r ̄ε: the radius discounts the j-th symmetric function by ̄r j , and M ̄r ≤ 1/2 by the readout period. Step 2: the reduction operator does not amplify. Let R denote the operator “mul- tiply by x, then reduce modulo Π”, written with a script letter to keep it distinct from the readout period T . In the monomial basis R is the companion matrix of Π: it sends e i 7→ e i+1 for i < M − 1, columns of ℓ 1 norm exactly one, and sends e M−1 to the coefficient vector of x M mod Π = ̄r M − δ(x), of ℓ 1 norm at most ̄r M + D. Hence ∥R∥ 1→1 = max1, ̄r M + D = 1, the maximum being over columns. Multiplying by x only relabels coefficients except at the top, where it costs one reduction, and one reduction costs less than it saves. Step 3: M reductions are a scalar plus a controlled remainder. For degp≤ M − 1 the product ̄r M p is already reduced, so R M p = (x M p) mod Π = ( ̄r M − δ)p mod Π = ̄r M p− (δp mod Π),(S54) 28 that is,R M = ̄r M I +E withEp =−(δp mod Π), the scriptE again chosen to avoid the comb spacing’s ∆. Writing δp = P i b i x i gives ∥b∥ 1 ≤ D∥p∥ 1 by submultiplicativity of ℓ 1 under convolution, and (δp) mod Π = P i b i R i (e 0 ) with∥R i (e 0 )∥ 1 ≤ 1 by Step 2, so ∥E∥ 1→1 ≤ D. At the ideal point M reductions return the coefficient vector unchanged up to the factor ̄r M ; off it, each full turn also smears the answer across residue classes by δ, and the whole error analysis is the size of that smear. Step 4: summation. Write m− 1 = aM +s with 0≤ s < M . Then c m =R m−1 e 0 = (R M ) a R s e 0 , so ∥c m ∥ 1 ≤∥R M ∥ a ∥R s e 0 ∥ 1 ≤ ( ̄r M +D) a . Every m > M has a≥ 1, and each a admits exactly M values of s, whence S out ≤ M X a≥1 ( ̄r M + D) a = M ( ̄r M + D) 1− ( ̄r M + D) ≤ 2M ̄r M + D (S55) whenever ̄r M +D ≤ 1/2, which ̄r M ≤ (2M ) −M ≤ 1/16 and D ≤ e 2 ̄ε≤ 7/16 guarantee at ̄ε ≤ 1/4. Substituting Eq. (S53) gives Eq. (S49). Both sums are sums of moduli throughout, over μ inside ∥c m ∥ 1 and over m outside, as Eq. (S49) requires; no ℓ 2 quantity is used, and none is available to be traded for one. Global phase. If the ideal nodes carry a common phase e iθ , the binomial becomes x M − ( ̄re iθ ) M and every step is unchanged, | ̄re iθ | = ̄r. This proves both parts of the lemma: the extrapolation tail is bounded by the node radius alone in the generic case, and by the designed geometry in the second, in each case without any hypothesis on the weights. What the design buys is now a comparison inside one statement. At the designed point the tail decays as ̄r M = e −Mλ ⋆ T = e −(γ+γ g )T/2 , the full trace budget of Lemma 3 spent over the memory span; at a generic point it decays as e − min k Reλ k T , the slowest mode’s share alone. The design equalises the decay rates, converting the minimum into the mean. The ratio of the two exponents is bounded on both sides by results that use no design condition: min k Reλ k ≤ λ ⋆ by Lemma 3, and min k Reλ k ≥ γ/18K by Lemma 4, whence M ≤ Mλ ⋆ min k Reλ k ≤ 9M (1 + ρ γ ), ρ γ = γ g /γ.(S56) The certificate’s comb spacing therefore purchases a readout period shorter by a factor in this band, and nothing else: universality itself, with an explicit rate, is available at any admissible operating point through part (a). Lemma 10 (Design distinctness for the order-n selection). Under D1, D4, and D5, consider the conjugation-resolved eigenvalue-sum points Λ (S) ( ⃗ k) = P i∈S λ k i + P i/∈S ̄ λ k i over all orders n≤ N , mode multisets ⃗ k, and conjugation patterns S. Clas- sify each point by its signed multiplicity vector n∈ Z K (entry n k counts appearances of mode k in S minus appearances outside S; |n k |≤ N ) and its total multiplicity vector m∈ N K . Then: (a) Points with distinct signed vectors n ̸= n ′ have imaginary parts separated by at least g 0 β −(M−1) /4. (b) Points with equal signed vectors but distinct total-multiplicity vectors m ̸= m ′ are separated in their real parts: by at least γ/4K if the total orders differ, and 29 otherwise by the designed drive-channel margin ε rad ≡ min| γ g 2 P k (m k −m ′ k )α 2 k |≥ γ g 16K (4N ) −(q−2) > 0, with q the prime of condition D3. (c) The only remaining coincidences are the exact complement pairs S ↔ S c at equal multisets, which satisfy Λ (S) = Λ (S c ) identically; the reality constraint F (s) ∗ = F (s ∗ ) of Sec. S.2.3 prescribes these jointly. Consequently the constructive weights can prescribe every conjugation-resolved mo- ment H (S) n (⃗μ) to a chosen real value, and choosing these values as the symmetrized target combinations yields |H (S) n (⃗μ)|≤ 2 n max|h (M,N ) n |. Proof. The strategy is to classify every eigenvalue-sum point by its multiplicity vector, so that two points coincide only when their vectors do, and then to show the design separates the finitely many vectors that survive. Every eigenvalue decomposes, by Lemma 7, as λ k = iδ k +E k +r k with |r k |≤ 8e 2 0 /∆ 0 . A sum of n≤ N eigenvalues or conjugates therefore has Im Λ (S) = X k n k δ k + ImR,Re Λ (S) = X k m k E k + ReR, |R|≤ N 8e 2 0 ∆ 0 . (S57) (a) For two points with n̸= n ′ , write the digit-difference vector d = n− n ′ ,|d k |≤ 2N , d̸= 0. The unperturbed separation is X k d k δ k = ∆ 0 X k d k k + g 0 X k d k β −k .(S58) If the comb part is nonzero, |∆ 0 P d k k| ≥ ∆ 0 , which dwarfs everything else. If the comb part vanishes, let j be the smallest index with d j ̸= 0; then g 0 X k≥j d k β −k ≥ g 0 β −j 1− 2N β− 1 = g 0 β −j 1− 2N 4N = g 0 β −j 2 ≥ g 0 β −(M−1) 2 , (S59) where the parenthesis bounds the worst-case carry from all deeper digits by a geometric series: this is the uniqueness of signed-digit representations in base β = 4N + 1 with digits bounded by 2N . The dressing remainders shift each of the two points by at most N · 8e 2 0 /∆ 0 , and the third entry of Eq. (D1) guarantees 2N · 8e 2 0 /∆ 0 ≤ g 0 β −(M−1) /4, so half the unperturbed separation survives. (b) With n = n ′ the imaginary parts agree to within the remainders, and the real parts differ by P k (m k − m ′ k )E k = (γ/2K)(n− n ′ ) + (γ g /2) P k (m k − m ′ k )α 2 k . If the total orders differ, the first term contributes at least γ/2K while the second is at most 2Nγ g s α /2K = ρNs α γ/K, which D4 makes smaller than γ/4K; the remainders are smaller still by Eq. (D1). If the total orders agree, the separation is the pure drive- channel expression γ g 2 P k (m k −m ′ k )α 2 k with P k (m k −m ′ k ) = 0. Writing d = m− m ′ , we have d ̸= 0, P k d k = 0 and P k |d k | ≤ 2N , so Lemma 6 bounds the expression below by ε rad = γ g 16K (4N ) −(q−2) . The dressing remainders shift the two points by at most 2N · 8e 2 0 /∆ 0 in total, which the fifth entry of Eq. (D1) makes at most ε rad /2; half the separation survives. 30 (c) Complement pairs at equal multisets are exact conjugates by inspection, and the joint prescription in the F ± basis of Sec. S.2.3 is compatible with the reality con- straint; distinct-or-conjugate is exactly what that construction requires. This proves the lemma: at the designed operating point the eigenvalue-sum points of every or- der n ≤ N are distinct or exact conjugates, which is precisely the input the order-n selection of Sec. S.2.3 requires. The two separations are both design outputs: the angular floor g 0 β −(M−1) /4 from the graded defects, and the radial floor ε rad from the delay lock. Neither is a property measured at an operating point. Cond. Statement and roleConsumed by D1Comb spacing ∆ 0 bounded below by the five-entry maximum of Eq. (D1): keeps the dressing perturba- tive; caps its node displacement; protects the deepest graded separation; locks the delay; protects the radial separation (the last is the entry that binds) Lemmas 7, 8, 10 D2Readout period T on the grid ∆ 0 T = 2π(nM +1)/M : undressed node angles exactly the M -th roots of unity; a measurement setting only Lemma 8 D3Atom at a rational point ℓ a /ℓ = p/q of the standing wave, q the smallest prime ≥ 2K phys + 1, with band residue n 0 mod q admissible: fixes the emission profile, bounds s α , and floors the radial separation Lemmas 5, 6, 10 D4Rate ratio ρ bounded by Eq. (D4): the drive-channel spread separates radially without misplacing the nodes Lemmas 8, 10 D5Graded defects of Eq. (D5), base β = 4N + 1: supplies Q-independence of the detunings, which no comb poly- nomial in k provides beyond order deg +1, breaking every angular coincidence with margin g 0 β −(M−1) /4. Unlike D1–D4 this is a property of the fabricated res- onator spectrum rather than a runtime setting Lemmas 8, 10 D6Coupler flatness ε v ≤ ̄ε target /8Λ max , Eq. (D6): caps the readout coupler’s own node displacement. Physical data of the fabricated coupler rather than a setting; verified at the operating point Lemma 8 The output coupler’s profile is geometry-fixed physical data (Sec. S.1.1) and no condition sets it. Its flatness ε v is likewise physical data, but it is not unconstrained: D6 is the feasibility test that data must pass rather than a dial. D1–D5 are settings; D6 is a check. Supplementary Table S.2 The designed operating point, assembled. Each condition is stated in Step 2 at its point of use; the table exists so the complete set can be audited at a glance. Step 3: parameter selection, and the theorem’s quantifiers The proof closes by explicit parameter selection. Given a target with a convergent Volterra expansion and an accuracy ε > 0: 31 1. Choose (M,N ) so that the Boyd–Chua truncation error obeys δ MN ≤ ε/2 (term (i)); this fixes A = A(M,N ) from the target’s own kernels. Without loss of generality take M ≥ 2: enlarging M only decreases δ MN . 2. Set K = M , let q be the smallest prime ≥ 2K phys + 1, and fix the delay and band placement by condition D3. Set the node-placement budget ̄ε target = min n ε 4eAM , 1 4 o (S60) and the drive ratio ρ by D4, whose right-hand side depends only on ̄ε target , N , M and the bound s α ≤ q 2 of Lemma 5. Then set λ ⋆ = (γ + γ g )/2K. 3. Choose the readout period T on the D2 grid so that ̄r = e −λ ⋆ T ≤ min n 1 2M , ε 8MA 1/M o , ̄r M = e −(γ+γ g )T/2 ,(S61) i.e. T ≥ λ −1 ⋆ max n ln 2M, 1 M ln 8MA ε o = 2 γ + γ g max n M ln 2M, ln 8MA ε o = O M lnM + lnε −1 γ + γ g ,(S62) which the D2 grid permits at any comb spacing. The middle equality substitutes λ ⋆ = (γ + γ g )/2M at K = M , and it is worth pausing on because the factor is easy to count twice: the 1/λ ⋆ is already carried inside the displayed numerator, so the clock grows as Θ(M lnM ) in the memory depth and not as Θ(M 2 lnM ). The rate itself is the check. At the binding branch (γ + γ g )T/2 = M ln 2M , so ̄r M = e −(γ+γ g )T/2 = (2M ) −M , which is exactly the bound Step 1 uses at Eq. (S30); the discarded reading T = O(M lnM )/λ ⋆ would give (γ + γ g )T/2 = M 2 ln 2M and hence ̄r M = (2M ) −M 2 , contradicting it. The trace budget rather than the per-mode rate, sets the clock. 4. Choose the defect amplitude g 0 by D5 and then the comb spacing ∆ 0 by Eq. (D1). Then Lemma 8 gives ̄ε ≤ ̄ε target , whose two entries discharge the two standing hy- potheses: ̄ε ≤ 1/4 is the condition of Lemma 9(b), and with ̄r ≤ 1/(2M ) at M ≥ 2 it gives S out ≤ 2M ̄r M + eM ̄ε ≤ M , the counting hypothesis of Eq. (S26); while ̄ε ≤ ε/(4eAM ) gives S out ≤ ε/(2A), so term (i) ≤ ε/2 and the total error is at most ε. Every constant is explicit and no limit other than the finite construction is taken. Write T min =Λ min /λ ⋆ for the floor of item 3, with Λ min = maxln 2M, 1 M ln(8MA/ε) dimensionless, and set T max ≡ 2T min , Λ max ≡ λ ⋆ T max . Conditions D1’s second and third entries, D4, D5 and D6 are asserted at T max ; the remaining entries are T -free once g 0 is fixed, so every condition asserted at the ceiling holds at every admissible T ≤ T max . That the grid reaches into the window is then a one-line consequence rather than a hypothesis: the D2 spacing is 2π/∆ 0 , and the sec- ond entry of Eq. (D1) at T max gives ∆ 0 T min ≥ 128(e 0 T min ) 2 / ̄ε target ≥ 128 > 2π, using e 0 T min ≥ λ ⋆ T min = Λ min ≥ ln 4 > 1 and ̄ε target ≤ 1, so a grid point lies in [T min ,T max ]. 32 The selection order is acyclic: (M,N ) → q → ̄ε target → Λ min → Λ max → [check D6 against ε v ] → ρ → λ ⋆ → T min → T max → g 0 → ∆ 0 → grid → T . The bracketed step is the only one that is not a selection: Eq. (D6) is evaluated against fabricated data and either passes or fails, and if it fails the accuracy is unreachable on that hardware rather than at a different setting. ■ Scaling of the design parameters QuantityFormulaGrowing MShrinking ε Readout period TEq. (S62)Θ(M lnM )+ Θ(lnε −1 ) Nodebudget ̄ε target Eq. (S60)Θ 1/(AM ) Θ(ε) (small ε) Drive ratio ρEq. (D4), first en- try Θ 1/(AM 3 lnM ) Θ(ε/ lnε −1 ) Defect amplitude g 0 Eq. (D5)Θ 1/(AM 2 lnM ) Θ(ε/ lnε −1 ) Comb spacing ∆ 0 Eq. (D1), fifth en- try Θ AM 4 lnM (4N ) q−2 Θ(ε −1 lnε −1 ) Node radius ̄r M Eq. (S61)exponentially smallΘ(ε) Supplementary Table S.3 Scaling ledger of the designed operating point. The two regimes are growing M at fixed ε and shrinking ε at fixed (M,N ). A = A(M,N ) is the target-dependent prefactor of Eq. (S26); β = 4N + 1; q is the least prime ≥ 2M +1, and the drive-ratio row uses s α ≤ q 2 . Every formula column cites the equation it is read off, and no entry is asymptotic guesswork. In the readout-period row, λ ⋆ ∝ 1/M at the designed point K = M is already absorbed in Eq. (S62), whence Θ(M lnM ); the defect-amplitude row inherits that exponent through g 0 = ̄ε target /8T max . The exponential factor β M in the comb spacing is the defect-protection entry of Eq. (D1); it is a property of the designed certificate, and the paragraph on designed versus generic operating points below explains what replaces it on a fabricated device. At K phys > M in-band modes the exponent reads K phys − 1 (Remark 3). Reading the ledger. Each growth column is a product of factors already fixed above, and we write the chain out so that no exponent lives only inside a table cell. Take the accuracy fixed and the depth growing, on the branch ̄ε target = ε/4eAM of Eq. (S60), which binds once A > ε/(eM ); the Volterra prefactor carries M n−1 at order n, so this is every tar- get of order N ≥ 2. The budget is then Θ(1/(AM )). The clock follows from Eq. (S62) as Θ(M lnM ), so the defect amplitude g 0 = ̄ε target /8T max is Θ 1/(AM 2 lnM ) . The drive ratio divides the budget by the profile’s peak-to-mean ratio and by the dimen- sionless clock, ρ = ̄ε target /8s α Λ max with s α ≤ q 2 = Θ(M 2 ) under the delay lock D3 and Λ max = Θ(lnM ), giving Θ 1/(AM 3 lnM ) . The comb spacing is then forced, and it is worth seeing why it cannot be read independently of the row above it. Its binding entry is ∆ 0 ≥ 512NKe 2 0 (4N ) q−2 /γ g with γ g = ργ, so ∆ 0 ∝ K/ρ at fixed N and fixed e 0 = Θ(γ): the spacing is the mode count divided by the drive ratio, and its M -exponent is therefore exactly one more than the exponent of 1/ρ. With 1/ρ = Θ(AM 3 lnM ) this is Θ AM 4 lnM (4N ) q−2 , the entry tabulated. The physical reading is that the design pays twice for the same 33 radial margin: the drive channel must be weak enough not to misplace the nodes (D4), which shrinks the very separation D5 and D1 must then protect, and the comb must be opened far enough to protect it. Every polynomial factor here is dominated by (4N ) q−2 with q = Θ(M ), so the exponential entry is what a fabricated device actually feels; the polynomial prefactor matters only as an audit trail. With the depth fixed and the accuracy tightening, only two factors move: the budget is Θ(ε) and the clock acquires + Θ(lnε −1 ). Both g 0 and ρ divide a budget of order ε by a clock of order lnε −1 , so both are Θ(ε/ lnε −1 ); the comb spacing, proportional to 1/ρ, is Θ(ε −1 lnε −1 ); and ̄r M = Θ(ε) by Eq. (S61). Buying a decimal place costs a factor of ten in the comb spacing and a logarithm in the clock, and nothing in the mode count. S.3 Reach, Scaling, and the Operating Class This section collects the results that follow from the construction of Sec. S.2 without being part of the proof of Theorem 1: what a fixed geometry can reach, how the design parameters scale, the weaker operating class on which the guarantee still holds, the physical resources the certificate spends, the scaling corollary, and the density step. The accuracy envelope of a fixed geometry The selection of Step 3 can be inverted, which answers the operational question: with the hardware fabricated and the geometry already set, which accuracies are reachable by measurement settings alone? Corollary 1 (Fixed-geometry envelope). Fix the fabricated hardware and the geo- metric settings: a comb of K phys in-band modes at spacing ∆ 0 with defect amplitude g 0 and drive ratio ρ (satisfying D3–D5 at level ̄ε 0 ≡ 6 maxρs α ln 2K phys , g 0 T for the periods used). An accuracy ε is then achievable by measurement settings alone (a choice of M ≤ K phys , N , T on the D2 grid, local-oscillator spectrum, and window functions) whenever there exist (M,N ) with δ MN ≤ ε 2 , ̄ε ∆ 0 ,ρ,g 0 ,ε v ;M,T ≤ min n ε 4eA(M,N )M , 1 4 o , 2M ̄r M ≤ ε 8A(M,N ) ,(S63) with ̄ε the explicit four-term expression of Eq. (S45) evaluated at the fixed geometry, and with the distinctness margins of Lemma 10 (the graded floor g 0 β −(M−1) /4 and the radial floor ε rad of Eq. (S34), the prime q of condition D3 belonging to the fixed geometry) exceeding the dressing remainders at that geometry. The infimum of such ε is the accuracy envelope ε min of the fixed geometry. Below the envelope, two resources must move together: the geometric dial (larger quantization length ℓ, refining ∆ 0 and raising K phys ) and the retained bandwidth. They are coupled: refining ∆ 0 at fixed B = K∆ 0 lowers the spacing, while the first entry of Eq. (D1) demands ∆ 0 ≥ 8q √ K ̄e with q ≈ 2K, so that B ≳ 16K 5/2 ̄e: the bandwidth must grow at least as K 5/2 for the certificate to survive rising mode number, on top of the exponential growth of the fifth entry. Subject to that, the strict-scaling corollary of the main text guarantees the move is never wasted. 34 Proof. Immediate from Step 3 read backwards: at fixed geometry the only selected quantities are (M,N,T ) and the record-side freedoms, T ranges over the D2 grid whose spacing 2π/∆ 0 is available at fixed ∆ 0 , and the three displayed inequalities are precisely the conditions under which steps 1, 3, and 4 close with the geometric quantities held fixed. The third inequality is satisfiable at fixed λ ⋆ by taking T large on the grid, so the envelope is set by the second inequality and by M ≤ K phys . Two consequences of the corollary deserve emphasis. First, within the envelope, the retained modes are selected by the readout: the record-side machinery of Remark 2 confines the realized kernels to the M modes the measurement addresses, so no physical mode-selection mechanism is needed, and the question of how K modes are singled out from the one device has a one-word answer: measurement. Second, the envelope is a statement about a setting of the device, never about the device: every accuracy is reachable on the same fabricated hardware, and the envelope marks where the purchase switches from measurement settings to the geometric dial. The operating class, and the theorem stated on it The conditions the construction of this section supplies are sufficient, and the proof consumes less than all of them. We therefore isolate what is consumed, so that the guarantee can be stated on operating points a fabricated device can occupy. Definition S.1 (Operating class). A finite-K member Γ (K) g of the family of Sec. S.1.1, together with a readout period T , a measurement window T off and a readout order N , lies in A 0 (M,N ) when (O1) the detunings δ k are pairwise distinct; (O2) the eigenvalues λ k of Γ g are pairwise distinct; (O3) the overlaps are nonzero, with margins ω min ≡ min k |v ∗ geo v R,k | > 0 and a min ≡ min k |v ∗ L,k α| > 0; (O4) the conjugation-resolved eigenvalue-sum points of order n ≤ N , taken over all in- band modes, are pairwise distinct or exactly conjugate, with minimal separation δ Λ > 0, and the drive channel is active, γ g > 0; (O5) T off lies outside the discrete exceptional set of Sec. S.2.3; (O6) r ∗ ≡ max k |e −λ k T | < 1/(M + 1); (O7) Im(λ k − λ k ′ )T /∈ 2πZ for k ̸= k ′ . The subclass A 1 (M,N )⊂A 0 (M,N ) adds (O8) the weak-dressing condition ∆ 0 ≥ 8 √ K ̄e. Conditions D1–D5 imply (O1)–(O4) and (O7) through Lemmas 7 and 10, and D1’s coupler entry is (O8); (O6) is a condition on the readout period, and (O5) holds at every window outside a discrete set. Every entry is an open condition, so A 0 has nonempty interior. Write λ min ≡ min k Reλ k , which is strictly positive on A 0 by Lemma 4 and is bounded below by γ/18K on A 1 . Theorem S.1 (Universality with an explicit rate on the operating class). Let y be a continuous fading-memory target with a convergent Volterra expansion on K u max , let ε > 0, and choose (M,N ) so that the Boyd–Chua truncation error obeys δ MN ≤ ε/2, 35 fixing A = A(M,N ). Let the device be operated at any point of the operating class A 0 (M,N ) of Definition S.1 (whose conditions O1–O8 the proof consumes) with K ≥ M . Then for every readout period satisfying (O6)–(O7) and T ≥ 1 λ min ln 2eAM (M + 1) ε (S64) there exist weights W n (t) with |y k − ˆy k | ≤ ε on K u max . On A 1 the same holds with λ min replaced by its symbolic floor γ/18K. Proof. Kernel matching to depth M and order N is the construction of Secs. S.2.3– S.2.4, which consumes (O2), (O3), (O4) and (O7) for the selection and the Vander- monde inversion, and (O5) for prescribability at the operated window; for K > M the readout-side selection of Remark 2 confines the realized kernels to the M selected modes and consumes (O4) over all in-band modes. The extrapolation identity (Propo- sition 1) is algebraic and consumes only (O2). Lemma 9(a) then bounds the tail by S out ≤ eM (M + 1)r ∗ (1 + O(Mr ∗ )) under (O6), where r ∗ = e −λ min T is a maximum over node moduli whose exponent carries the minimum over decay rates. Requiring AS out ≤ ε/2 gives Eq. (S64), and adding the truncation error gives ε. No design condition is consumed at any step. In plain terms, the designed operating point of this section engineers the spectrum so that every constant can be written down in advance; the class above asks only that the comb lines be distinct, that the readout couple to every mode, that the eigenvalue sums not collide, and that the device be watched long enough. The price of the weaker hypothesis is that the rate constant λ min is a measured property of the operating point rather than a formula in the physical rates: one must look at the machine to know how long to watch it. Comparing the two at K = M , the designed period is Θ(M lnM/γ) and the class period is Θ(M ln(AM 2 /ε)/γ), the same order in the memory depth, so the designed point buys a bounded factor in the readout period rather than the availability of the guarantee. Designed and generic operating points The design D1–D5 is a constructive certificate: it exhibits, with every constant explicit, one operating point at which all hypotheses hold with proven margins. Its conditions are sufficient; we do not establish that any is necessary. They carry a stated price, the exponential defect-protection entry of Eq. (D1), and the guarantee does not depend on paying it (Lemma 9(a)). A generic fabricated device pays no such price: its eigenvalue sums are separated by its own spectral irregularity, with margins measured rather than prescribed, and Sec. S.6 verifies this by direct diagonalization at every operating point simulated here. Physical resources of the design The construction spends seven resources, and we state them plainly, for the designed certificate. (i) Time: T = O(M lnM + lnε −1 )/(γ +γ g ) per readout step, by Eq. (S62). (i) Bandwidth: B = K∆ 0 with ∆ 0 carrying the exponential defect-protection and 36 radial-protection factors of Eq. (D1); the flat-coupling budget ε flat (B) of the Discus- sion is charged accordingly. A generic operating point replaces this entry by its verified margins and pays only the polynomial second entry of Eq. (D1). (i) Timing preci- sion: the D2 condition must hold to a phase accuracy of order ̄ε target , which at spacing ∆ 0 means a relative timing precision of order ̄ε target /(∆ 0 T ), which falls exponentially in M through ∆ 0 , so the certificate’s clock is as far beyond laboratory reach as its comb spacing; growing precision requirements on settings are the standard price of universal-approximation constructions (classical weight precision grows likewise), and we state the requirement rather than leave it implicit. (iv) Shots: at the designed operating point the graded separations are exponentially small in M , so the construc- tive weights there, while never entering the error bound, are exponentially large, and realizing them at fixed output precision costs shots accordingly; the theorem is an ap- proximation statement, and its sample-complexity accounting is Sec. S.7’s. (v) Rate ordering : condition D4 drives the drive-channel rate far below the mean mode decay, by the two entries of D4, so the designed certificate sits outside the adiabatic-slaving ordering γ g ≫ λ ⋆ under which the linear-transducer limit is derived from the physical atom (Methods). The certificate’s requirement s≤ 10 −2 is not in conflict with this: it is imposed at the designed point, where D4 caps ρ near 10 −4 , so the admissible drive is far smaller than the simulated one and the saturation bound is slack by orders of magnitude. Which entry supplies that floor must be said, because the entries do not scale alike. The defect-protection entry gives 768M 2 ln(2M ) 5 M−1 , about 2.1× 10 4 at M = 2. Both expressions are read off Eq. (D1) by the same substitutions, and we give them once so the reader can check either: e 0 = γ/2 because the Hermitian part has rank two with unit vectors (Lemma 3); T = Λ/λ ⋆ = 2M Λ/γ at K = M ; γ g = ργ with ρ at its D4 ceiling ̄ε target /(8s α Λ max ) and s α ≤ q 2 from D3; Λ max = 2Λ min ; and B = K∆ 0 = M ∆ 0 . The third entry 384Ne 2 0 Tβ M−1 / ̄ε target then gives B/γ ≥ 192M 2 Λβ M−1 / ̄ε target , and the fifth 512NKe 2 0 (4N ) q−2 /γ g gives B/γ ≥ 2 2q+6 M 2 q 2 Λ max / ̄ε target ; at N = 1, β = 5 and ̄ε target = 1/4 these are the two displayed forms. Their ratio grows as (16/5) M when q is taken as 2M +1, and jumps with q at the primes. The radial -protection entry, the fifth, gives more: with ρ at its D4 cap, at s α ≤ q 2 and with the band constraint M ∆ 0 ≤ B, it reads B/γ ≥ 2 2q+7 ̄ε target M 2 q 2 Λ min ,Λ = λ ⋆ T, q the least prime ≥ 2K phys + 1, (S65) including the factor-two margin taken in the pair selection of Step 3. The prefactor is a single power of two because every factor entering it is one: 512 from the fifth entry of Eq. (D1), e 2 0 = γ 2 /4, the 8 of the D4 ceiling on ρ, and ̄ε target = 1/4; at that budget the prefactor reads 2 2q+9 . At the minimal readout period Λ = ln 2M this is about 7.3× 10 7 at M = 2, 6.6× 10 9 at M = 3, and 8.6× 10 12 at M = 4 (deposited script resources v36.py), and it degrades by Λ/ ln 2M when the accuracy branch of Step 3 binds. Since q ≥ 2M + 1, the fifth entry exceeds the third by a factor growing like (16/5) M , so it is the fifth entry that binds and Eq. (S65) that states the reach. The platform reading of the main text follows: the certified construction fits no depth 37 M ≥ 2 at the superconducting example (B/γ ≈ 10 2 ), and M = 2–3 only for kHz- class linewidths under GHz-flat coupling. What this floor prices is the comparison following Lemma 9: a readout period shorter by the factor of Eq. (S56) rather than the availability of the guarantee, which Lemma 9(a) supplies at any admissible operating point. Condition D6 carries a cost of the same kind, and we state it because it is the entry that decides what is buildable. Since ε v ≤ ̄ε target /8Λ max by Eq. (D6) and the geometry gives ε v ≈ B/ω 0 with ω 0 the transition frequency, the two combine into a carrier-to-linewidth requirement ω 0 γ ≥ B γ · 8Λ max ̄ε target = 2 2q+11 ̄ε 2 target M 2 q 2 Λ 2 (S66) which at ̄ε target = 1/4 reads 2 2q+15 M 2 q 2 Λ 2 , about 6.4× 10 9 at M = 2. The relax- ation of the node budget from 1/(4M 3/2 ) to 1/4, which the remainder estimate of Lemma 9(b) permits, improves this by a factor 8 at M = 2 and by 4M 3/2 in general; at the retired budget the same expression gives 5.2× 10 10 . Even relaxed, the require- ment is severe: it is the coupler’s flatness, rather than the comb spacing, that sets the hardest demand on the fabricated device, and no platform named in this work meets it at M = 2. Check of the designed operating point The design checks of Sec. S.2.5.1 are inequalities between explicit expressions and are evaluated by the deposited scripts at operating points constructed from Table S.2. They check the construction; no lemma of this Supplement rests on them. S.3.1 Capability containment and strict scaling: proof of the corollary The strict-scaling corollary of the main text asserts that enlarging the accessible mode number never decreases approximation capability, and strictly enlarges the span of exactly matchable kernel blocks. We prove the two halves separately. Proposition 2 (Capability containment and strict growth of the span). LetR(M,N ) denote the set of target-accuracy pairs (y,ε) the (M,N )-device achieves, and let S(M,N ) denote the set of kernel blocks h n (⃗m) : n ≤ N, ⃗m ∈ [1,M ] n it realizes exactly. Let M ≤ M ′ and N ≤ N ′ . Then (i) Containment. R(M,N )⊆R(M ′ ,N ′ ). (i) Strict growth of the span. S(M,N ) is the full block space of symmetric kernels, of dimension P n≤N M +n−1 n , the kernels being identified only through their symmet- ric part, as the inversion of Sec. S.2.4 over ordered indices makes explicit. Hence S(M,N ) ⊊ S(M ′ ,N ′ ) strictly whenever (M ′ ,N ′ ) ≥ (M,N ) componentwise with at least one strict inequality. Proof. (i). Let (y,ε)∈R(M,N ): the parameter selection of Step 3 chose (M,N ) from the target so that the truncation error satisfies δ MN ≤ ε/2, and met the remaining budget through Eq. (S26) by driving the extrapolation tail below ε/(2A(M,N )). We 38 Supplementary Figure S.2 Verification of the designed operating point (geometry-fixed coupler). a, Measured extrapolation tail P m>M ∥c m ∥ 1 , computed directly from the node set at the extreme admissible modulus ̄ε = 1/4 and ̄r = 1/(2M ), against the bound of Lemma 9(b) (dashed) and its two constituents. The node term dominates the bound over the whole range while the radial term 2M ̄r M is already negligible, and the measurement sits below the bound by the factors annotated, rising monotonically from 22 at M = 4 to 45 at M = 14. The plotted value is the worst case over the exhaustive deterministic phase family θ k = 2πjk/M , j = 0,...,M − 1, rather than a random sample: random draws under-sample the M -torus non-uniformly in M and produce a spurious non- monotonicity. b, Realized node-placement error against the assigned budget ̄ε target at M = 6, N = 2, with the three contributions of Eq. (S45) resolved: the placement is realized at 1 6 of budget with unit slope over three decades, consumed almost entirely by the graded-defect phase g 0 T , with the radial spread two orders below it and the dressing displacement 8e 2 0 T/∆ 0 some fifteen further orders down. Two remarks on scope. Both panels run well past the depth at which the certified operating point is physically reachable—the bandwidth requirement of Eq. (S65) already exceeds any demonstrated platform at M = 4—and every design inequality of Eq. (D1) is asserted and holds at each point plotted; the construction is therefore verified over a range where the certificate itself is not realizable, which is what makes the margins in a and b informative rather than merely reassuring. And both margins are one-sided: the bounds are conservative, the conservatism grows with M , and none of it is fitted. exhibit settings of the (M ′ ,N ′ )-device achieving the same ε in three steps. The intu- ition is that a larger device can be told to imitate a smaller one by being handed the smaller one’s kernels and zeros everywhere else, and that the zeros cost nothing. Step A (zero-padded matching). Operate the (M ′ ,N ′ )-device at any point at which its nodes are pairwise distinct and its conjugation-resolved eigenvalue sums of order n ≤ N ′ are resolved; the designed point of this section is one such, and a generically fabricated device is another, by Sec. S.6. By Lemma 8 its M ′ nodes are pairwise distinct, and by Lemma 10 at base β = 4N ′ + 1 every conjugation-resolved eigenvalue- sum point of order n ≤ N ′ needed by the matching is resolved, so the inversion of Sec. S.2 prescribes any values on the (N ′ ,M ′ ) block. Prescribe the target’s truncated kernels h (M,N ) n (⃗m) on the sub-block n≤ N , ⃗m∈ [1,M ] n , and zero on the remainder. This is an admissible choice of the weights W n (t); nothing else about the device changes. Step B (zeros extrapolate to zero). By Proposition 1 the realized kernels at all lags are the tensor extrapolation Eq. (S25) of the matched-block moments, and Eq. (S25) is linear in H n separately at each order. At orders n ∈ (N,N ′ ] the matched block is identically zero, hence H n ≡ 0, hence the realized kernels vanish at every lag: the padding introduces no spurious response at the padded orders. At orders n ≤ N the realized functional agrees with the truncated target on [1,M ′ ] n , the values on 39 Supplementary Figure S.3 Why weight-norm routes fail: motivation for the weight-independent proof. a, For nodes confined to an arc (the geometry forced when a sub-revival condition T P > MT caps the angular span), ln∥V −1 ∥ ∞ grows linearly in M (R 2 > 0.999): any bound routed through the matched-weight norm inherits exponential growth. b, The growth rate is set by the angular span rather than the radius. c, Forcing |x k | = 1 does not repair it. The proof of this subsection avoids the weight norm entirely (Proposition 1), and the design removes the sub-revival constraint that created the arc. [1,M ′ ] n \ [1,M ] n being the target’s own zeros, and deviates from it only at lags ⃗m /∈ [1,M ′ ] n . Step C (tail budget). The deviation of Step B is controlled by Eq. (S26) applied to the (M ′ ,N ′ ) matching. The zero-padded block has the target’s own kernel sup-norms, zeros not increasing a supremum, so the prefactor is A ′ ≡ P N n=1 n 2 2n u n max M ′n−1 ∥h (M,N ) n ∥ ∞ < ∞, the sum terminating at N because the padded orders carry no realized content by Step B; the tail factor is the M ′ -node S ′ out ≤ 2M ′ ̄r ′M ′ + eM ′ ̄ε ′ of Lemma 9(b). By the trace identity ̄r ′M ′ = e −(γ+γ g )T ′ /2 with T ′ free, and ̄ε ′ is driven below any tolerance by the comb spacing (Lemma 8); the parameter selection of Step 3, run verbatim at (M ′ ,N ′ ) with budget ε/(2A ′ ), there- fore delivers S ′ out ≤ ε/(2A ′ ), with the same quantifier order and no circularity. Since δ MN ≤ ε/2 is a property of the target alone and is unchanged, the total error is at most ε. (i). Surjectivity of the matching onto the full block space is the inversion itself: node distinctness (Lemma 8) makes the order-one prescription invertible, and design distinctness (Lemma 10) makes the order-n prescription well posed for all n ≤ N , so every block value is independently prescribable. The dimension count and the strict inclusion follow directly, every block supported on the enlarged index set and not on the smaller one being realizable at (M ′ ,N ′ ) and not indexable at (M,N ). No auxiliary hypothesis is required. This establishes both halves of the proposition: the capability sets are nested as the memory depth and Volterra order grow, and the nesting is strict, since each enlargement of the index set carries a kernel block the smaller device cannot index at all. Remark 4 (Sharper prefactor). Tracking the support of the nonzero matched entries (contained in [1,M ] n ) through the counting bound of Eq. (S26) replaces M ′n−1 by M n−1 in A ′ , recovering the original A(M,N ); the loose constant above is used in the proof because it requires no support bookkeeping, and either constant is finite and target-determined, which is all Step C consumes. 40 Two cautions delimit what has been proved, and they are the reason (i) and (i) are separate statements rather than one. Containment is a statement about achievable accuracies rather than about exact reproduction of output functionals: distinct mem- bers of the family have distinct spectra, no spectral nesting being assumed (Sec. S.1.1), so their extrapolation tails differ, and the (M ′ ,N ′ )-device reproduces the (M,N )- device’s outputs to every tolerance rather than identically. This is the operative sense in which scaling is never wasted, and the only sense the theorem’s own error metric defines. Strictness, meanwhile, is proved at the level of the exactly matchable span: a capability-level separation, a target and an ε achievable at (M ′ ,N ′ ) and provably not achievable at (M,N ) over all of the smaller device’s admissible settings, would require a lower bound we do not claim. Since deeper truncation strictly reduces δ MN for any tar- get with content at the added lags, the span statement is what the strict-containment claim of the abstract asserts. S.3.2 Verification of the hypotheses of the density theorem The passage from finitely many matched kernels to approximation of the full functional invokes Theorem 1 of Cuchiero et al. [3]. We verify its hypotheses in our setting. The input space is K u max = [−u max ,u max ] Z − —inputs uniformly bounded by u max , the bound already used in Lemma 13; the symbol B is reserved for the bandwidth— compact in the product topology by Tychonoff’s theorem. The target is a continuous fading-memory functional on K u max , i.e. continuous with respect to a weighted norm ∥u∥ w = sup m w m |u −m | with w m ↓ 0; onK u max this topology coincides with the product topology, so the target lies in C(K u max ). The finite Volterra polynomials u7→ Q i u −m i form a point-separating subalgebra of C(K u max ) containing the constants: two distinct sequences differ at some lag m, separated by the degree-one monomial u7→ u −m . By Stone–Weierstrass this algebra is dense in C(K u max ), which is the density hypothesis of the cited theorem; the kernel-matching construction of Secs. S.2.3–S.2.4 realizes every element of the algebra exactly (for K ≥ M ), which is its realizability hypothesis. The quantification is over the family: exact realization of the degree-(N,M ) monomial block needs K ≥ M modes, so density is attained along the nested family of operating points as M grows—the family-level, rate-free sense in which merely continuous targets are reached, exactly as Theorem 1’s statement records. We emphasize that Stone– Weierstrass enters here only to characterize the closure of the target class; the strict- scaling corollary of the main text continues to follow from kernel containment, which density arguments do not provide. Theorem 1 adds three physical statements to the abstract density results it in- vokes. Boyd–Chua [4] and Cuchiero et al. [3] establish that some functional with rich enough kernels approximates any fading-memory target; they say nothing about which hardware realizes those kernels, at what rate, or with what scaling. Here: (i) the atom– mirror filter’s own kernels are provably rich enough to match any target up to (N,M ) under generic conditions on one mirror distance; (i) the convergence rate is the trace budget of the fixed hardware; and (i) enlarging the accessible mode number never reduces capability and strictly enlarges the matchable span. 41 S.4 Fading Memory of the Device The device has fading memory, and the fact has two statements of different strength. We give the elementary one first, because it follows from the kernel bound alone and because it is all that is needed at a fixed operating point, and then prove the stronger statement that survives the many-mode limit. The elementary statement, at fixed mode number. In the linear-transducer limit the reduced generator L 0 is quadratic with a c-number drive, so a constant drive relaxes to a coherent steady state with amplitude Γ −1 g α, and the input–output kernel is the finite sum of exponentials h K (t) = v ∗ geo e −Γ g t α = P k (v ∗ geo v R,k )(v ∗ L,k α)e −λ k t . Every Reλ k > 0 at finite K (Lemma 4, which needs only the flat coupler and the distinctness of the δ k ), so |h K (t)| ≤ C(K)e −ε D (K)t with ε D (K) = min k Reλ k > 0: the kernel decays, the influence of any input recedes, and Corollary 2 follows at that K. Nothing more is required for any single operating point, and no delay equation is involved. Why the elementary statement is not enough. Its rate is not uniform in the mode number. By the trace identity of Lemma 3, ε D (K) ≤ (γ + γ g )/2K, so the elementary bound degrades as 1/K and says nothing about the family as a whole—even though the physical loop that produces the mem- ory is unchanged as the comb is refined. The discrepancy is not an artefact of a loose bound: the slow eigenmodes are the weakly coupled ones, and they carry almost none of the input–output response. The physically meaningful statement is therefore about the kernel rather than the spectrum, and it is uniform: on any sub-revival horizon the kernel decays at the loop rate c 0 set by delay stability, with constants independent of K. That is the content of the remainder of this section. S.4.1 Uniform kernel decay under delay stability Lemma 11 (Uniform kernel decay under delay stability). Two kernels appear in this lemma, and we define them separately because they obey different equations. The finite-comb kernel is the matrix element h K (t) = v ∗ geo · e −Γ g t · α,(S67) the input–output impulse response of the K-mode Gaussian-limit generator, drawn from the nested family whose coupling density carries the fixed band envelope of Theorem 1 (flat with smooth Gaussian rolloff ). The continuum kernel h ∞ (t) is its un- bandlimited limit. It is h ∞ , and only h ∞ , that obeys the delay-differential equation of the reduced linear transducer: at any finite K the comb kernel h K satisfies no delay equation (it is a finite sum of exponentials, and its revival echoes are the signature of the difference). The physical difference between the two kernels is the mirror’s apparent distance. In h ∞ the atom faces a genuine continuum: the emitted pulse reflects once, returns after the delay τ , and then escapes forever out the open end, and this one-echo causal structure is exactly what a delay-differential equation encodes. In h K the field lives on a 42 discrete comb, and a discrete comb is a periodic system: K equally spaced tones cannot dephase forever and rephase after the revival time T P = 2π/∆ 0 . Physically, finite mode spacing means finite quantization length, so the pulse has a second boundary to return from: h K shows the fast loop transient, decays as if the field were gone, and partially reconstructs at t≈ T P , the revival echo, which is the arrival of radiation a true continuum would have carried to infinity. The two kernels therefore agree on the sub- revival window, where the finite comb does not yet know it is finite, and differ precisely at and beyond T P ; this is why h K obeys no delay equation and why the uniform bound below is honestly restricted to horizons T H < T P . The characteristic equation of the continuum delay dynamics is s + γ 2 1 + re iφ e −sτ = 0, r = γ R /γ ≤ 1, φ = π− ω 0 τ.(S68) Here γ R is the emission rate into the returning mirror channel as seen at the atom, round-trip losses included, so that r = γ R /γ ≤ 1 is the returned fraction of the total decay, r = 1 the lossless ideal mirror, and r the quantity the feedback transmissivity η of Sec. S.10 attenuates (see the paragraph on detection versus feedback loss in Sec. S.7); and φ = π − ω 0 τ is the round-trip phase, the propagation phase ω 0 τ at the atomic transition frequency plus the π of the mirror reflection. (The characteristic equation describes the undriven loop, so its phase is evaluated at the atomic frequency; the main text’s interaction-picture phase π− ω L τ is the same quantity in the frame rotating at the drive, the two coinciding on resonance and differing by the detuning contribution ∆τ off it.) Consistency is internal: φ = π corresponds to sin(ω 0 τ/2) = 0, the atom at a field node, which is exactly the dark-state set excluded by the overlap conditions, while φ = 0 is the antinode of maximal collective decay, matching the static anchors of Sec. S.10. Define the delay-stability condition D(c 0 ): all roots of Eq. (S68) satisfy Re[s] ≤ −c 0 < 0. Let T P = 2π/∆ 0 be the revival time of a comb of spacing ∆ 0 . If D(c 0 ) holds, then for every horizon T H < T P there exist C h (c 0 ,T P ) <∞ and K 0 , both independent of K, such that for all K ≥ K 0 , |h K (t)| ≤ C h e −c 0 t/2 for all 0≤ t≤ T H ,(S69) i.e. on any sub-revival horizon the kernel decay is uniform in the mode number. Refin- ing the comb (∆ 0 ↓ 0) raises T P and hence extends the horizon on which the uniform bound holds; in particular, whenever the comb is refined so that T P exceeds the mem- ory horizon MT used by the readout, the bound holds with K-independent constants over the entire range of lags that enter the output. Proof. We first derive the delay equation for h ∞ , since the rest of the proof rests on it. In the continuum the linear transducer couples to the semi-infinite waveguide through the standing-wave density ̃g(ω) 2 = (γ/π) sin 2 (ωτ/2) of Sec. S.1. Writing sin 2 (ωτ/2) = 1 2 (1− cosωτ ) and passing to the frame rotating at ω 0 under the flat-band axiom A2, Wigner–Weisskopf elimination of the field gives the memory kernel as the Fourier transform of the density, Z dν γ 2π 1− cos((ν + ω 0 )τ ) e −iνt = γ δ(t)− γ 2 e −iω 0 τ δ(t− τ )− γ 2 e iω 0 τ δ(t + τ ). The rate appearing here is the atom’s mirror-modified emission into the guided field; in the mode-space generator it is the profiled channel α k ∝ sin(ω k τ/2) that carries this 43 structure, the frequency-flat readout channel being delta-correlated and contributing no delay. The causal contributions (the local term counting one half, the advanced term discarded) yield ̇ h ∞ (t) = − γ 2 h ∞ (t) + γ 2 e −iω 0 τ h ∞ (t − τ ) for t > τ , with the memoryless segment on [0,τ ]; since e −iω 0 τ =−e iφ at the round-trip phase φ = π−ω 0 τ (the π being the mirror reflection), and a lossy return scales the delayed term by r = γ R /γ, this is ̇ h ∞ (t) =− γ 2 h ∞ (t)− γ 2 re iφ h ∞ (t−τ ), whose characteristic equation is Eq. (S68). Solutions of a linear autonomous delay equation admit the spectral bound |h ∞ (t)|≤ Ce (σ 0 +ε)t for any ε > 0, with σ 0 the supremum of real parts of characteristic roots [5]; D(c 0 ) gives σ 0 ≤−c 0 , hence the bound for h ∞ at rate c 0 (take ε = c 0 /2). Uniformity in K requires relating the finite-comb kernel to h ∞ , and the relation is not a plain Riemann approximation. Two steps are involved, and we keep them separate. First, the fixed band envelope of the theorem’s nested-family setup acts on the continuum kernel by convolution: writing ψ for the (inverse) Fourier transform of the fixed envelope shape, h (env) ∞ = ψ∗ h ∞ ,i.e. h (env) ∞ (t) = Z ψ(s)h ∞ (t− s) ds,(S70) with ψ fixed once and independent of K because the envelope shape is. Second, for a comb of spacing ∆ 0 the Poisson summation formula relates the finite-comb memory kernel C K to periodic repetitions of the envelope-shaped one, C K (t) = P p∈Z C (env) ∞ (t− pT P ); the repetitions are disjoint on [0,T P ), so C K ≡ C (env) ∞ there, and since the response solves a causal Volterra equation with a unique solution, the responses agree exactly on the same window. The response is not itself a periodization of h ∞ —the map from memory kernel to response is a resolvent rather than a linear one—and the display below is written for the memory kernel, with the echo picture used only as an illustration of what happens beyond T P : h K (t) = X p≥0 h (env) ∞ t− pT P 1(t− pT P ), T P = 2π ∆ 0 ,(S71) where the p≥ 1 terms are revival echoes of the finite comb. Each echo is a copy of the decaying kernel restarted at t = pT P . Equation (S71) is exact for the equidistant comb; which is the case the lemma claims, together with the designed comb of Sec. S.2.5.1, whose defects satisfy|ω k −ω 0 −∆ 0 k|≤ g 0 . The quadratically chirped comb ∆ 0 +2∆ 2 |k| has varying local spacing, the repetitions are not exact, and it lies outside the scope of this statement. Eq. (S71) is written here for the idealized comb, in which the repetitions are exactly disjoint; with the smeared envelope the echoes acquire tails that leak marginally across t = pT P , and we account for that leakage explicitly below. The band envelope is what makes the p = 0 term decay uniformly in K, and it does so without any appeal to an ideal (brick-wall) band cutoff—whose filter would have divergent L 1 norm (the sinc kernel), and whose one-sided transform decays only as 1/|ω|, so that a naive band-limitation constant would grow like logB K rather than staying K-uniform. Instead, band-limitation with the fixed envelope is the convolution Eq. (S70), whose kernel ψ possesses a finite exponential moment, Z |ψ(s)|e c 0 |s| ds ≡ C ψ < ∞, K-independently.(S72) 44 This is where the band-edge shape enters, and the sufficient condition is exactly Eq. (S72): the envelope must possess a finite exponential moment at rate c 0 , equiva- lently its frequency profile must extend analytically to a horizontal strip of half-width c 0 . A merely smooth rolloff need not satisfy this, but the condition is far from Gaussian- specific: an exponential edge of width σ B satisfies it whenever σ B > c 0 , while a Lorentzian edge fails. The Gaussian edge of fixed width σ B used here gives a Gaus- sian ψ of width 1/σ B in time, for which every exponential moment is finite and C ψ = e c 2 0 /(2σ 2 B ) up to normalization. The constant therefore depends only on the ratio c 0 /σ B of the loop-stability rate to the fixed band-edge width—a single physical num- ber, K-independent because σ B is fixed once with the envelope shape. We note that σ B is the band-edge scale and is logically distinct from the bandwidth B over which the flat-coupling budget ε flat (B) of the Discussion is assessed: the envelope constrains how the density rolls off at the edges, ε flat how flat it is across the interior. The band-limitation step then collapses to a three-line convolution bound. Using |h ∞ (t− s)|≤ C e −c 0 (t−s) 1(t− s) from delay stability D(c 0 ) in Eq. (S70), |h (env) ∞ (t)| ≤ Z |ψ(s)||h ∞ (t− s)| ds ≤ C e −c 0 t Z s<t |ψ(s)|e c 0 s ds + (leakage from s > t) ≤ C h e −c 0 t/2 ,(S73) where the first integral is bounded by C ψ via the exponential moment Eq. (S72) (it absorbs the e c 0 s ), and the leakage from the causal tail s > t—the region where 1(t−s) has cut off h ∞ —is Gaussian-small in t because ψ inherits the Gaussian rolloff, hence absorbed into C h . On a horizon t ≤ T H < T P only the p = 0 term of Eq. (S71) is present (θ(t− pT P ) = 0 for p≥ 1), so |h K (t)| =|h (env) ∞ (t)|≤ C h e −c 0 t/2 with C h built from C ψ and C alone. The constants C h and K 0 therefore depend on (c 0 ,T P ) and the fixed envelope but not on K, which is the uniformity claimed; the pointwise exponential envelope is now derived from Eq. (S72) rather than asserted, so no separate ripple- crossing estimate is needed. Beyond one revival (t ≳ T P ) the echoes reconstruct the kernel and no fading-memory bound of this form can hold at fixed spacing—the honest content of the trace identity (Lemma 3). The uniform bound therefore holds only on horizons below the revival time. The echo leakage noted after Eq. (S71) is Gaussian- small in B K (T P −T H ) and is absorbed into C h , so the sub-revival bound is unaffected. Whether a given operating point lies below its own revival is a joint statement about the readout period and the comb spacing; we do not establish it for the simulations reported here, and the uniform-in-K constants of this lemma are correspondingly not claimed at those operating points. Remark 5 (The contraction and sub-revival requirements cannot both be met). Two requirements bear on the readout period from opposite sides, and they do not overlap. The extrapolation-tail estimate of Lemma 9(a) requires the contraction max k |e −λ k T | < 1/(M + 1), hence T > ln(M + 1)/λ min with λ min ≡ min k Reλ k the slowest decay rate of the retained block; the uniform bound of Lemma 11 requires the memory horizon to lie below the revival, MT < T P . The trace identity (Lemma 3) caps λ min ≤ λ ⋆ = (γ + γ g )/2K, and a comb of K modes filling a band B has T P = 2πK/B. Chaining 45 the three, 2KM ln(M + 1) γ + γ g < MT < 2πK B ⇐⇒ B M ln(M + 1) < π(γ + γ g ), (S74) in which the mode number cancels identically: refining the comb raises the revival time and lowers the slowest decay rate in the same proportion, so the two effects annihilate. Axiom A2 places the band far above both rates, so the right-hand side of Eq. (S74) is smaller than the left for every M ≥ 1, and no operating point of the finite comb meets both requirements. In plain terms, the device must be watched long enough for its slowest mode to decay, and that time grows with the mode number in exactly the proportion that the comb’s rephasing time does; adding modes buys both sides equally and never closes the gap. No step of the approximation bound of Sec. S.2.5 is affected, since none of it uses fading memory; what Eq. (S74) withdraws is the possibility of extending the uniform-in-K constants of Lemma 11 over the lags the readout uses. Relation to the approximation theorem. The uniform decay proven here characterizes the device: under delay stability the physical input–output kernel decays at the loop rate on sub-revival horizons, uni- formly in the mode number, which is the fading-memory property of Corollary 2 and the measurable signature of the feedback loop. It is a separate statement from the ap- proximation bound of Sec. S.2.5, which is finite-dimensional and spectral and in which neither C h nor c 0 appears. The two are complementary: the theorem certifies what a chosen operating point can approximate, and this section certifies that the machine itself forgets. The constant C h at the operating comb. The lemma asserts existence of a finite, K-independent C h on the sub-revival horizon; the constant itself should be evaluated at the operating comb, and we report its measured behaviour. On a sub-revival window T H = 0.5T P , with the comb refined as the theorem prescribes (R = T P /T growing with K), the equidistant comb gives C h = 0.9251, 0.9231, 0.9226, 0.9223, 0.9222 at K = 6, 10, 14, 20, 28: stable to four significant figures, a direct confirmation of the claimed K-uniformity. The equidistant comb is not the worst case: at K = 14 and R = 21 on T H = 0.8T P , moderate chirps raise C h to between 1.8 and 3.7 times its undispersed value, while chirped combs place their echo reconstruction at or beyond T P rather than before it. The margin to the revival also matters: at T H = 0.8T P the equidistant C h itself begins to grow with K (0.92 → 1.73 → 5.30 at K = 14, 20, 28) as the horizon approaches the revival, so the sub-revival window should be taken with room to spare rather than pressed against T P ; we accordingly operate, and recommend, T H ≤ 0.5T P . Nothing in the error bound of Sec. S.2.5 depends on C h . Exponent bookkeeping. Three rates appear, and we track the factors of two once. Delay stability D(c 0 ) places all characteristic roots at Re[s]≤−c 0 , so the continuum kernel decays at rate c 0 (up 46 to any ε > 0 in the spectral bound). Lemma 11 states its conclusion at the halved rate c 0 /2: the slack absorbs the ε of that spectral bound (we take ε = c 0 /2) and the enve- lope’s exponential-moment constant C ψ , leaving a clean pointwise bound C h e −c 0 t/2 . The fading-memory corollary (Cor. 2) halves once more, to c 0 /4: the weighted-norm statement splits the per-lag decay factor into one half used to sum the geomet- ric tail into the Lipschitz constant L and one half retained as the weight sequence w m = e −c 0 (m−1)T/4 . No step is tight in these constants, and no qualitative claim depends on the halvings. Measured envelope constants. The uniform bound is an analytic statement about the idealized many-mode nested family, proven above with no appeal to numerics. Independently, the finite-K devices actually simulated are verified to be members ofC by direct diagonalization (Sec. S.6), and on those devices the kernel envelope can be read off directly: the fitted decay rates are c eff = 0.15γ (K = 6) and 0.09γ (K = 14), and the measured envelope constants sup t≤T P |h K (t)|e +c eff t/2 are 0.56 and 0.63 respectively: finite, K-stable over the sub- revival horizon, and comfortably below the bound the lemma proves at the operating scale. These constants are a consistency check on the simulated devices rather than an ingredient of the proof. The fixed envelope disturbs nothing else in the construction: it is positive on the occupied band, so the overlap conditions v ∗ geo ·v R,k ̸= 0, v ∗ L,k ·α̸= 0 still hold; it is τ -independent, so the almost-every-τ genericity of Sec. S.6 is untouched; it preserves the normalization ∥α∥ = ∥v geo ∥ = 1, so the trace identity (Lemma 3) is unchanged; and it is a property of the many-mode family only, so the finite-K Vandermonde construction of Secs. S.2.3–S.2.4 is unaffected. Where D(c 0 ) holds. Solving Eq. (S68) via the Lambert W function, D(c 0 ) fails only at the exact anti- node resonance φ = π (mod 2π) with r = 1, where s = 0 is a root: the dark-state configuration at which the atom decouples, the same measure-zero set excluded by the overlap conditions. Numerically (Supplementary Fig. S.4), at the simulated φ = π/3 the margin is c 0 = 0.62γ at γτ = 0.5 and c 0 = 0.36γ at γτ = 1.0, and the worst case over all φ remains strictly stable for γτ ≤ 3. D(c 0 ) is a stability requirement on the feedback loop, the natural quantum analogue of the echo-state condition of classical delay-line reservoir computing, and it replaces the mode-space gap in Theorem 1. S.4.2 Fading-memory property of the composed map Corollary 2 (FMP of the realized reservoir functional). Under condition D(c 0 ), the composed input–output map ˆy :K u max → R realized with matched weights satisfies, for all u,v ∈K u max , |ˆy(u)− ˆy(v)| ≤ Lsup 1≤m≤⌊T P /T⌋ e −c 0 (m−1)T/4 |u −m − v −m |,(S75) with L = L(u max ,N,C W ,C h ,c 0 ,T ) <∞ uniform in K ≥ K 0 for combs refined so that T P > MT , where C W ≡ max n≤N sup t∈[0,T off ] | f W n (t)| is a bound on the fixed matched weight functions. Hence ˆy has the fading-memory property in the weighted-norm sense 47 Supplementary Figure S.4 The delay-stability condition D(c 0 ). Rightmost characteristic root of Eq. (S68) versus γτ at the simulated phase φ = π/3 (solid) and in the worst case over φ at r = 1 (dashed). The simulated regime γτ ∈ [0.5, 1] (shaded) has margin c 0 ≥ 0.36γ. of the echo-state literature, with weight sequence w m = e −c 0 (m−1)T/4 , uniformly along the mode-number family, over every lag the matched weights use: those are supported on the last M windows and MT < T P by hypothesis. Beyond the first revival the finite-K kernel still decays, at the rate min k Reλ k , so the all-lag statement holds with the K-dependent weight e −(min k Reλ k )(m−1)T ; by Lemma 3 no K-uniform weight is available over all lags, and none is needed, since no step of the approximation bound uses fading memory. Remark 6 (Fading memory is the device’s; truncation is the readout’s). Two distinct statements are easily conflated here, and we separate them. Corollary 2 is an all-lag statement: the weighted bound holds over every m ≥ 1, which is what the fading- memory property means in the echo-state sense, and it is a property of the device— it follows from Lemma 11, which bounds the physical kernel and knows nothing of any readout. The memory depth M , by contrast, belongs to the readout: the matched weights are supported on the last M on/off windows, and the resulting truncation is the approximation step quantified by term (i) of Sec. S.2.5. The device does not “have memory M ”; it has fading memory at rate c 0 , and the readout chooses to read M lags of it. The sub-revival caveat attaches to the certification rather than to the property: it is over the horizon MT , kept below T P by comb refinement, that the uniform-in-K constants of Lemma 11 are available. Proof. The output is a polynomial of degree ≤ N in the linear functionals ℓ t (u) = Re[v ∗ geo e −Γ g t ζ(u)] = P m≥1 Re[h K (t + (m − 1)T )c on ]u −m (with c on the fixed on-window factor of Eq. (S13)), whose coefficient at lag m is bounded by |c on |C h e −c 0 (m−1)T/2 by Lemma 11, where c on = R T on 0 e −Γ g s ds acting on α defines the fixed on-window injection factor entering ζ 1 (Eq. (S13));|c on | is a K-uniform constant by the same contractivity bound. On the bounded domain the polynomial is Lipschitz in each ℓ t with constant set by u max , N , C W ; summing the geometric tail and splitting the decay factor gives the stated weighted-sup bound. 48 S.5 Continuity of the Volterra Kernels in the Saturation Parameter The main text states that the gap between the linear-transducer limit (in which Theorem 1 holds) and the fully saturable device (which the simulations run) is a quan- tified discrepancy of order the saturation parameter s = (ε INPUT /γ g ) 2 . We record the perturbative statement here. Lemma12(Kernelcontinuityin s).Fixthereservoirparameters (Ω,α,v geo ,γ,τ,T on ,T off ) and bounded inputs |u k | ≤ u max . Let the full generator Eq. (S1) and the Gaussian-limit generator Eq. (S3) be run from a common initial state on the sector E (n max ) of Lemma 13, driven by the same input word, over a comparison horizon of M + 1 readout periods. Then the output ˆy k of the full generator has Volterra kernels at lags m i ≤ M satisfying ˆ h full n (m 1 ,...,m n ) = ˆ h G n (m 1 ,...,m n ) + O(s),(S76) where ˆ h G n are the kernels of the Gaussian-limit generator Eq. (S3), with the explicit bound ˆ h full n (⃗m)− ˆ h G n (⃗m) ≤ |s| C ext (n,u max )∥O∥ ∞ C V g √ n max + 1 (M +1)T,(S77) where ∥O∥ ∞ bounds the readout observable on the sector and C ext (n,u max ) is the coefficient-extraction constant of the input monomials. The constant depends on u max , n, m i , the fixed reservoir parameters, the comparison horizon (M +1)T , and the sector size n max of Lemma 13; it is therefore a per-operating-point statement rather than a K-uniform one. Proof. Work in the displaced frame of Lemma 13 and on the bounded-excitation sector E (n max ) there constructed, on which the physical trajectory is supported up to tails bounded by e −n max /(2β 2 max ) , absorbed into the constant below. The sector is finite- dimensional (K modes carrying at most n max total excitations), so every operator on it is bounded; write the full generator there as L full = L G + sV with ∥V∥ ≤ C V g √ n max + 1 (Lemma 13). Analyticity of the one-period propagator. The on/off drive is T -periodic, so the driven trajectory is a fixed point of the one-period propagator P s . Expanding P s in a Dyson series in sV about the Gaussian propagation, the order-n term is an n-fold time- ordered integral bounded in norm by (|s|∥V∥T ) n /n! times the uniformly bounded Gaussian segment norms, the piecewise-constant switching entering only through the segmentation of the integrals; the series converges absolutely for every s, and P s is entire in s on the sector. Comparison of the two trajectories. Write ρ full (t) and ρ G (t) for the states gener- ated by L full and L G from the common initial state under the same input word, and U full (t,τ ′ ) for the propagator of the full generator. Duhamel’s formula is exact, ρ full (t)− ρ G (t) = Z t 0 U full (t,τ ′ ) [sV]ρ G (τ ′ ) dτ ′ .(S78) Take trace norms. The propagator sits on the left of each product and is completely positive and trace preserving, so∥U full (t,τ ′ )∥ 1→1 ≤ 1 for every pair of times—a supre- mum over τ ′ rather than an average—while ∥ρ G (τ ′ )∥ 1 = 1 and ∥V∥≤ C V g √ n max + 1 49 by Lemma 13. Hence ∥ρ full (t)− ρ G (t)∥ 1 ≤ |s|C V g √ n max + 1 t,(S79) an inequality running in one direction only and containing no spectral quantity: con- tractivity of a quantum channel is available without a gap, and it is all that is used. The bound grows linearly in t, which is why the horizon is a hypothesis of the lemma rather than a remark upon it. Passage to the kernels. The output Eq. (S10) is, for each fixed weight set, a bounded linear functional of the trajectory on the sector: writing O for the readout observable and ∥O∥ ∞ for its operator norm there, duality of the trace and operator norms gives |ˆy k (s)− ˆy k (0)| ≤ ∥O∥ ∞ ∥ρ full − ρ G ∥ 1 , which Eq. (S79) bounds at t = (M +1)T . Both sides are polynomials in the input word, and the input monomials are linearly indepen- dent on K u max , a set with nonempty interior; the coefficient functionals are therefore bounded, with a constant C ext (n,u max ) depending on the monomial degree and the input bound alone. Matching coefficients yields Eq. (S77). No spectral quantity en- ters at any step, and the constant depends on K only through n max . This proves the lemma: the Volterra kernels move at most linearly in the saturation parameter, at a rate fixed by the input bound and the sector size alone. The Gaussian-limit kernels are therefore the leading term of a controlled expansion rather than a separate object. The invariant bounded-excitation sector. The proof above requires the residual superoperator V to be bounded, and it is bounded only on a sector of finitely many excitations. We supply that sector here. Lemma 13 (Invariant sector). Let ˆ N = σ + σ − + P k a † k a k and let E (n) denote the sec- tor spanned by density operators supported on eigenspaces of ˆ N with eigenvalue ≤ n. Work in the displaced frame ρ7→ D † (β d (t))ρD(β d (t)) with β d (t) =− i g αu(t), in which the input enters only through the displacement. The residual generator conserves ˆ N up to the dissipative channels, which are excitation-non-increasing; the frame displace- ment is uniformly bounded, sup t ∥β d (t)∥ ≤ u max /g, while the coherent excitation of the state itself is carried by the driven trajectory, sup t ∥ζ(t)∥≤ ζ max = ηu max ∥Γ −1 g α∥; it is the latter that the sector must contain. Hence for inputs |u k |≤ u max the physical state remains, up to a tail bounded by the Chernoff estimate P(N ≥ n)≤ e −λ (eλ/n) n at λ = ζ 2 max , in E (n max ) with n max the least integer making that tail at most the tar- get accuracy δ (equivalently n max = max⌈4ζ 2 max ⌉ + 1, ⌈log 2 δ −1 ⌉ suffices), and on E (n max ) the residual atom–field coupling satisfies ∥V∥ E(n max ) ≤ C V g √ n max + 1. The earlier prescription n max = ⌈4ζ 2 max ⌉ + 1 is recovered whenever ζ 2 max ≳ 1, where the exact Poisson tail at that n max is already below 10 −3 ; only in the weak-drive corner ζ 2 max ≲ 1/2 must n max be raised, and then by one, costing a factor p 4/3 in the bound on ∥V∥. Proof. Under the rotating-wave interaction of Eq. (S1), [ ˆ H int , ˆ N ] = 0; the drive term is removed by the displacement; the Lindblad channels v ∗ geo · a and α ∗ · a annihilate excitations. The displaced dynamics therefore mapsE (n) into itself for every n. Undo- ing the displacement adds a coherent amplitude of at most β max per mode direction, whose Poissonian number tails give the stated bound. On E (n max ), ∥a k ∥ ≤ √ n max and ∥σ ± ∥ = 1, bounding V. 50 Remark 7 (Status of the constant β max ). The constant β max = ηu max ∥Γ −1 g α∥/g is an operating-point quantity with a physical name: it is the device’s closed-loop gain, the steady displacement a bounded input can sustain against the loop’s damping. Three properties fix its status. It is finite exactly when every dressed mode is damped (Reλ k > 0 for all k), which the delay-stability and overlap conditions guarantee; at the dark- state configuration φ = π with unit feedback an undamped root appears, Γ −1 g ceases to exist, and a bounded resonant input does grow the excitation without bound, which is the same excluded set as everywhere else in the paper and the configuration probed by the φ = π task-level floor of Sec. S.10. A slowly damped, strongly driven mode makes β max large rather than infinite: growth is transient over ∼ 1/Reλ k and saturates at gain times input, which is what the sector bound then prices. Finally, β max grows along the mode-number family as the slow modes’ response grows (the trace identity closes the slowest rate as 1/K, and a near-resonant component of the gain grows accordingly), so the sector bound is a per-operating-point statement, and no uniformity in K is claimed or needed: the kernel-continuity lemma of this section is invoked at fixed operating points only. This is the sector on which Lemma 12 applies, and the constant of Eq. (S77) depends on n max , i.e. only on u max , the fixed reservoir parameters and the comparison horizon, as stated there. Remark 8 (Scope of the constant). The constant of Eq. (S77) is K-dependent: it en- ters through n max , which grows along the mode-number family with β max (Remark 7). The estimate is therefore a per-operating-point statement and no uniformity in K is claimed or needed; the uniform-in-K content of the paper lives entirely in Lemma 11, which is proven on the scalar kernel. The steady-trajectory version of the comparison is not established here, and the obstruction is structural rather than technical: the full generator Eq. (S1) carries a single dissipative channel, the γ g term of Eq. (S3) being produced by the adiabatic elimination rather than inherited from it, so the saturable device is less dissipative than its Gaussian limit and no fixed-K contraction transfers to it. That extension belongs to the finite-s problem stated in the Discussion. Remark 9. Lemma 12 does not extend Theorem 1 to finite s (that remains the open problem stated in the Discussion) but it changes the category of the theorem–device relation: the physical device’s kernels lie within an explicitly controlled, continuously tunable distance of a family proven universal, rather than in a separate regime. For the simulations of the main text the operating values of s are computable from the stated parameters. A direct numerical closure of this bound would extract the first- and second-order Volterra kernels of the simulated saturable device and overlay them on the closed-form Gaussian-limit kernels of Eq. (S15), checking agreement to O(s) at the operating drive; that is the natural verification, and the mode-number convergence study of the main text (Fig. 2, saturable overlay) already exhibits the physical device converging along the Gaussian-limit axis, consistent with the continuity asserted here. S.6 Genericity of the Universality Conditions The proof requires that the spectrum of Γ g satisfy the linear-independence (non- resonance) and overlap conditions. We show these hold for almost every mirror distance 51 and verify them numerically for the simulated devices. We first record the equivalence used implicitly above. The two non-resonance conditions are generic. The argument has two halves: the conditions are first recast as a Q-independence statement, which is then shown to hold for almost every delay. Lemma 14 (Genericity of the non-resonance conditions).(a) For K modes, the con- ditions P k Re[λ k ]n k ̸= 0 and P k Im[λ k ]n k ̸= 0 for all integern k with P k n k = 0, P k n 2 k > 0, are equivalent to linear independence over Q of Re[λ k − λ K ] k<K and of Im[λ k − λ K ] k<K respectively. Either one alone implies the complex non- degeneracy P k λ k n k ̸= 0 used in Theorem 1 (if a real or imaginary part is nonzero the complex sum cannot vanish), so the genericity result below, proven for the real part and hence for the complex sums, delivers exactly what the theorem’s hypothe- sis requires; the imaginary-part condition, which the theorem does not consume, is scoped separately in Remark 10. (b) Regard the eigenvalues λ k (τ ) of Γ g as functions of the delay τ . Away from the isolated τ at which eigenvalues cross, each λ k (τ ) is real-analytic. For every inte- ger tuple n k with P k n k = 0, P k n 2 k > 0, the resonance function R n (τ ) = P k n k λ k (τ ) is not identically zero, and neither is its real part P k n k Reλ k (τ ). This holds for every comb with pairwise-distinct frequencies—equidistant, chirped, or de- fected alike—with no disjointness condition on any frequency set. Consequently the set of τ violating the complex and real-part non-resonance conditions is a countable union of discrete sets and therefore has Lebesgue measure zero, and by Lemma 14 this alone discharges every non-resonance hypothesis of Theorem 1. Proof. (a) The constraint P k n k = 0 lets one eliminate n K = − P k<K n k , giv- ing P k<K n k (λ k − λ K ); vanishing for some nonzero integer tuple is exactly Q-linear dependence (integer and rational dependence coincide after clearing denominators). (b) Consider the one-parameter family Γ g (τ,t) = iΩ + t E(τ ), E(τ ) = γ 2 v geo ⊗ v ∗ geo + γ g 2 α(τ )⊗ α(τ ) ∗ , t∈ [0, 1], (S80) so that t = 1 is the device and t = 0 the bare comb; the ray carries the full physical dressing, both channels at their physical rates. The entries are jointly analytic in (τ,t) wherever the normalization D(τ ) = P j sin 2 (ω j τ/2) is positive, which excludes only a discrete set of delays. At t = 0 the operator is diagonal with eigenvalues iω k : pairwise distinct for any comb with distinct frequencies, independent of τ , and simple, so each λ k (τ,t) is jointly analytic near t = 0 with no crossing analysis required at the expansion point, and first-order perturbation theory at a diagonal unperturbed operator [6] returns exactly the diagonal entries of the perturbation: λ k (τ,t) = iω k + tE k (τ ) + O(t 2 ), E k (τ ) = γ 2K + γ g 2 |α k (τ )| 2 ,(S81) using the geometry-fixed flat coupler for the first term. For any admissible tuple the flat-coupler contribution cancels, P k n k γ/2K = 0, leaving ∂R n ∂t t=0 = γ g 2 X k n k |α k (τ )| 2 = γ g 2 N (τ ) D(τ ) , 52 N (τ ) = X k n k sin 2 ω k τ 2 = − 1 2 X k n k cos(ω k τ ),(S82) where P k n k = 0 removed the constant term. The frequencies ω k being pairwise distinct, the functions cos(ω k τ ) are linearly independent as functions of τ , so N ̸≡ 0: the first-order coefficient is a not-identically-zero analytic function of τ , for every comb with distinct frequencies and with no disjointness of any frequency set invoked anywhere. Joint analyticity and the identity theorem then give R n ̸≡ 0 on the connected domain, so R n (·,t) ̸≡ 0 for every ray parameter t outside a discrete exceptional set, and away from that set the zero set of R n (·,t) in τ is discrete. The union over the countably many integer tuples, together with the countable crossing set at t = 1, is measure zero, and the countably many finite K contribute countably many such sets. Finally, the first-order coefficient is real, so the identical argument applies verbatim to P k n k Reλ k (τ,t): real-part genericity is delivered outright, and with it, by Lemma 14, the complex non-degeneracy Theorem 1 consumes. This proves the lemma: the non-resonance conditions fail only on a set of mirror distances of measure zero, so a device drawn at a generic geometry satisfies them, and Theorem 1 may consume them as a hypothesis without excluding any but an exceptional set of devices. Remark 10 (Scope of the imaginary-part condition). We first locate where the con- dition is consumed. Writing a conjugation-resolved sum point over a multiset m with pattern p as Λ(p) = P k m k Reλ k +i P k (2p k −m k ) Imλ k , the real part depends on m alone; two conjugation patterns on a common multiset therefore have identically equal real parts and are separated only by P k b k Imλ k with b = p−p ′ ,|b| 1 ≤ n. This is exactly the separation Lemma 10(a) supplies at the designed point, through the graded defects and with no sum-zero restriction on b; it is what the defect hierarchy is for. Away from the design, at first order along the ray the imaginary part of R n is P k n k ω k +O(t 2 ), which on the equidistant comb vanishes identically for every tuple with P k n k k = 0: the bare equidistant comb does not satisfy imaginary-part genericity at leading order, a further instance of the regularity-versus-selection theme of Sec. S.2.5. The imaginary- part condition therefore requires either higher-order terms or a comb whose frequencies avoid the cross-frequency coincidences, and a generic chirp restores it at leading order: for ω k = ω 0 +∆ 0 k +∆ 2 k 2 , a coincidence 1 2 (ω j +ω j ′ ) = ω m requires both (j +j ′ )/2 = m and (j 2 +j ′2 )/2 = m 2 , impossible for j ̸= j ′ by strict convexity of the square, so the col- lisions occupy finitely many ∆ 2 values at each K and almost every chirp avoids them all. None of this is consumed by Theorem 1, whose hypotheses need only the complex distinctness that the real part alone supplies (Lemma 14); the numerical verification below checks the operating combs directly. Remark 11 (Why functional rather than pointwise, independence). A tempting short- cut would infer pointwise Q-linear independence of the perturbative coefficients from their distinctness; this is false in general (distinctness does not imply Q-independence, and for an equidistant comb the relevant consecutive differences are in fact rationally dependent). The proof above therefore requires no arithmetic condition at any point: non-vanishing as a function of τ , supplied by the linear independence of cosines at dis- tinct frequencies, is all the identity theorem needs. The bare equidistant comb is itself the cautionary case twice over: its arithmetic regularity defeats not only the pointwise 53 shortcut but the leading-order imaginary-part expansion (Remark 10), which is why the argument is arranged, through the diagonal expansion point, to need neither. Remark 12. Lemma 14 is what a naive cardinality argument misses: the device spectrum is not a free point in R K but an analytic curve τ 7→ λ k (τ ), and generic- ity along that curve is what must be, and is, established. The bare equidistant comb Ω = diag(ω 0 + k∆) is maximally resonant on the imaginary axis; the dressing by γ 2 (v geo ⊗v ∗ geo ) and γ g 2 (α⊗α ∗ ) lifts the degeneracy of the complex sums for almost ev- ery τ through the real parts, per the proof above, while the imaginary-axis coincidences themselves persist at leading order on the equidistant comb (Remark 10). Numerical verification. This section concerns generically fabricated devices, which do not satisfy condition D3 and are the devices the simulations of this paper use. It carries no hypothesis of Theorem 1: every condition of the theorem is discharged at the designed operating point by Lemmas 5, 6, 10 and 4. We diagonalized the exact Γ g over a dense sweep of τ for the mode numbers used in the main text. Across the swept range the eigenvalues are distinct with Re[λ k ]≥ ε D > 0; the overlaps |v ∗ geo · v R,k | and |v ∗ L,k · α| are bounded away from zero except at the isolated node distances where sin(ω k τ/2) = 0; and the finite-order non-resonance surrogate min n: P n k =0, P |n k |≤6 | P k n k λ k (τ )| stays strictly positive away from a discrete set of τ . The parameter values used for every figure in the main text and this Supplement satisfy all conditions, so those simulated reservoirs satisfy the overlap and non-resonance conditions in the linear-transducer limit; they do not satisfy the weak-dressing entry of D1, and the explicit constants that entry supplies are correspondingly unavailable at these operating points. In addition, the delay-stability condition D(c 0 ) of Lemma 11 is verified directly on the character- istic roots of Eq. (S68): at the simulated phase φ = π/3 the margin is c 0 = 0.62γ (γτ = 0.5) and c 0 = 0.36γ (γτ = 1.0), and the worst case over all phases remains strictly stable for γτ ≤ 3 (Supplementary Fig. S.4). S.7 Shot-Noise-Corrected Moment Readout From the intracavity quadrature to the detected field. The proof is written in terms of q geo = Re[v ∗ geo ·a], a quadrature of the standing-wave mode combination selected by the measured channel, whereas the laboratory observ- able is the traveling field at the detector. The two are related by the input–output boundary condition on the open end of the waveguide. The mirror terminates one side and supports no outgoing channel, so the detected continuum is the unidirectional field past the atom, and for the measured channel v geo the standard relation reads ˆ b out (t) = ˆ b in (t) + √ γ v ∗ geo · a(t) ,(S83) with ˆ b in the vacuum input to that channel, [ ˆ b in (t), ˆ b † in (t ′ )] = δ(t− t ′ ). Balanced ho- modyne detection of ˆ b out at local-oscillator phase θ measures ˆ b out e −iθ + ˆ b † out e iθ , i.e. 2 √ γ q geo plus the vacuum contribution of ˆ b in . Equation (S83) is what licenses reading the binned outcomes I j below as estimates of the windowed q geo : the signal term is the intracavity quadrature the proof uses, scaled by √ γ, and the ˆ b in term integrates 54 to the noise contribution ξ j . We stress that Eq. (S83) is an operator identity and that its two terms are correlated: the same vacuum input ˆ b in that appears explicitly also drives a(t) through the Langevin equation of the measured channel, so the fluctuation part of the windowed quadrature and the integrated shot noise of the same bin are not independent. The statistics of the binned outcomes are therefore not asserted here; they are derived in Lemma 15, where the measurement back-action is accounted for exactly and the measured output field is shown to be a displaced vacuum in the linear- transducer limit. Because ˆ b in acts on an incoming continuum that never returns to the atom (the mirror side carries no outgoing channel) Eq. (S83) introduces no additional feedback and no non-Markovian structure beyond that already carried by (Ω,α). Detection loss is not feedback loss. Two distinct efficiencies appear in this device and should not be conflated. The feed- back transmissivity η swept in Sec. S.10 attenuates the returning field and therefore alters the dynamics—it changes r = γ R /γ in the characteristic equation Eq. (S68) and with it the stability margin c 0 , which is why severing it (η = 0) collapses the device to the feedback-free plateau. Detector inefficiency η det < 1, by contrast, acts only on the outgoing channel Eq. (S83) after the physics has happened: it admixes vacuum, ˆ b det = √ η det ˆ b out + √ 1− η det ˆv, leaving the generator Eq. (S3), the kernel h K , and the rate c 0 untouched. Because a beamsplitter with a vacuum ancilla is itself a passive map, it carries the displaced vacuum of Lemma 15 to a displaced vacuum—the mean rescaled by √ η det , the fluctuations still exactly vacuum—so the binned outcome be- comes I j = √ η det γ ̄q j + ξ j with Var[ξ j ] = v 0 unchanged, and the triangular relation of Lemma 15 holds verbatim with γ 7→ η det γ: the ideal moments are still recovered exactly in expectation. The cost is entirely in the shot budget: the effective signal-to- noise scale ̄v = v 0 +q 2 max becomes v 0 +η det q 2 max , so by Lemma 15 the shot number for a target precision at order n grows by at most η −n det —a benign factor at the n = 1 read- out used for every real-world benchmark, and the reason imperfect detection degrades the statistics of this device rather than its dynamics. Shot-noise-corrupted moments. The proof uses ideal moments ⟨q n geo ⟩; experimentally one has the balanced-homodyne photocurrent, which in continuous measurement contains white shot noise. Point- wise powers of the raw current are therefore not well-defined objects (they would involve powers of white noise) so all moment estimation below is built on time-binned outcomes: for each readout sample time t j the current is first integrated against a square-integrable window w j matched to the discretized weight grid, I j = Z w j (t)dI(t),(S84) with the windows w j square-integrable and pairwise orthogonal (in practice, non- overlapping bins). Two quantities must be kept distinct in what follows. The windowed quadrature is the operator R w j (t)q geo (t) dt, whose fluctuation part is correlated with the same-bin shot noise, as noted above. The windowed amplitude ̄q j is the c-number window average of Re[v ∗ geo · ζ(t)] on the coherent driven trajectory of Sec. S.2—a 55 deterministic functional of the input history, and the quantity in terms of which the statistics of I j are stated below. Powers and empirical moments are taken of the binned outcomes I j , one per window per experimental repetition; with this definition in place, the statements below are well-posed (a loose phrasing in terms of “powers of each shot” would not be, precisely because of the white-noise issue), and v 0 ∝∥w j ∥ 2 denotes the vacuum variance of a binned homodyne outcome. The readout is a moment estimator, and its two costs—what back-action does to the estimate, and how many shots the estimate needs—are one statement. Lemma 15 (Moment readout and shot budget).(a) In the linear-transducer (Gaus- sian) limit, for pairwise orthogonal windows, and provided readout begins after a washout interval long compared with 1/c 0 —so that the state has relaxed to the co- herent driven trajectory of Sec. S.2 and any correlations carried by the initial state have decayed—the binned outcomes I j are jointly Gaussian with mean √ γ ̄q j and covariance exactly v 0 δ j ′ , independently of the drive. Consequently the empirical moments of I are linear combinations of the ideal moments⟨q j geo ⟩, j ≤ n, with fixed, known coefficients; any output linear in ⟨q j geo ⟩ is realized by a linear reprocessing of ⟨I j ⟩, and the required weights are obtained by an invertible relabeling. (b) Let ˆμ n denote the empirical n-th moment of the homodyne outcome over N shots repetitions, with vacuum-noise variance v 0 and signal bounded by q max . Then Var[ˆμ n ] ≤ (2n− 1)!! (v 0 + q 2 max ) n N shots ,(S85) so precision σ on the order-n moment requires N shots ≥ (2n− 1)!! (v 0 +q 2 max ) n /σ 2 . Writing ̄v ≡ v 0 +q 2 max for this second-moment scale, the requirement reads N shots ≥ (2n− 1)!! ̄v n /σ 2 . Composing this with the approximation bound is one line, and we record it because the two are otherwise stated separately. Writing ˆy (R) k for the R- shot estimate of the output Eq. (S10) and ∥W n ∥ 1 = R T off 0 |W n (t)| dt, linearity gives |ˆy (R) k − ˆy k |≤ P n ∥W n ∥ 1 σ n , so an end-to-end accuracy ε requires R ≥ max 0≤n≤N 4(N +1) 2 (2n− 1)!! ̄v n ∥W n ∥ 2 1 ε 2 .(S86) The weight norm does not enter the approximation bound, which is the content of Proposition 1; it does enter the sample complexity, quadratically. The reading is that the theorem’s guarantee is on expectation values, and that the cost of realis- ing it in finite time is set by the size of the weights that realise it, which at the designed operating point is governed by separations that are exponentially small in M . Bounding ∥W n ∥ 1 at a generic operating point is not attempted here. For the Gaussian-limit readout at order N the shot cost thus grows as (2N − 1)!!v N 0 ; for the saturable device with a linear readout (N =1) it is ∼ v 0 /σ 2 . The nonlinearity transfer of the main text is therefore not merely an apparatus simplification but a super-exponential reduction of the measurement budget. 56 Proof. (a) Statistics of the binned outcomes. The Gaussian-limit generator Eq. (S3) is passive—quadratic with no anomalous (a † a † ) terms—and stable (Reλ k ≥ γ/18K > 0, Lemma 4), and its only dissipation channels are the readout channel √ γ v ∗ geo · a and the drive channel √ γ g α ∗ · a, both fed by vacuum; the input enters solely as the c-number displacement. The input–output map is therefore linear scattering: in frequency, ˆ b out (ω) = S(ω) ˆ b in (ω) + d(ω) with the 2× 2 scattering matrix S(ω) = I − C (Γ g −iω) −1 C † , C = ( √ γ v ∗ geo ; √ γ g α ∗ ), and d(ω) the c-number displacement driven by u. For a passive stable system with all dissipation channels listed, the scattering matrix is Σ(ω) = I − CG(ω)C † with G(ω) = (Γ g − iω) −1 and C the 2× K coupling matrix with rows √ γv ∗ geo and √ γ g α ∗ , so that 1 2 C † C = E and Γ g + Γ † g = C † C. Then G † C † CG = G † [(Γ g − iω) + (Γ † g + iω)]G = G † + G, whence Σ † Σ = I − CG † C † − CGC † + C(G † + G)C † = I at every real ω, the resolvent existing there by Lemma 4. The identity is exact. Unitarity of the measured row means that the measured output channel, fed by vacuum on both inputs, is again an exactly δ-correlated vacuum field, displaced by d: the output state of the measured channel is a displaced vacuum. One hypothesis of the lemma enters exactly here: the frequency-domain scattering relation is a stationary statement, valid when the field is initially in vacuum or, equivalently for the readout, once the transient of the stable dynamics has decayed and the state has reached the coherent driven trajectory. At transient times the correlations carried by the initial state survive, and binned outcomes in different windows would in general remain correlated; this is why readout is taken to begin after the washout interval of the lemma’s statement—a period standard in reservoir-computing practice, and included in every simulation of this paper (Sec. S.10). With the transient washed out, integrating against orthogonal windows gives the first statement: I j jointly Gaussian, mean √ γ ̄q j , covariance v 0 δ j ′ , with v 0 the vacuum value regardless of the drive. Where the back-action went. The correlation flagged before Eq. (S83) is real: writing the windowed quadrature as ̄q j + δ ̄q j with δ ̄q j the operator fluctuation, Var[I j ] = γ Var[δ ̄q j ] + v 0 + 2 √ γ Cov[δ ̄q j ,ξ j ], and the first step fixes the total at v 0 : the back- action cross term is negative and cancels the intracavity-fluctuation excess exactly. This is the content of vacuum-in/vacuum-out for a passive network rather than an accident, and it is why no independence assumption between signal and noise appears anywhere in this proof. Moment recovery. The measured moments ⟨I n ⟩ are the moments of a Gaussian with mean √ γ ̄q and variance v 0 : ⟨I n ⟩ = P n j=0 n j γ j/2 ̄q j μ n−j with μ r the central Gaussian moments of variance v 0 (μ 0 = 1, μ 1 = 0, μ 2 = v 0 , higher μ r by Isserlis). This is lower-triangular with unit diagonal in the pair (⟨I n ⟩,γ n/2 ̄q n ), hence invertible, recovering ̄q and its powers exactly in expectation. The ideal moments⟨q n geo ⟩ are, in the coherent state of the driven trajectory, the moments of a Gaussian with mean ̄q and the state-independent vacuum quadrature variance fixed by the commutator; they are therefore obtained from the powers of ̄q by a second triangular, unit-diagonal relation with fixed coefficients. Composing the two, the ideal moments are linear combinations of the measured ones with fixed, known coefficients, and these absorb into the readout weights precisely as the normal-ordering coefficients did in Eq. (S11). No apparatus beyond the homodyne detector is required. This proves the lemma: every ideal moment 57 the construction needs is recovered exactly in expectation from the homodyne record, and the price of doing so is paid entirely in shot count rather than in bias. Remark 13 (Scope of the moment inversion; readout order for the saturable device). The displaced-vacuum argument above is a statement about the linear-transducer limit and is used only there. For the fully saturable device the output field is not Gaussian, and the signal–noise correlation of the second proof step does not cancel against a vac- uum total; no moment-inversion claim is made for it at order n≥ 2. None is needed: every real-world benchmark in this paper operates the saturable device with a purely lin- ear readout, n = 1 (Table S.4), where linearity of expectation gives⟨I j ⟩ = √ γ⟨ ̄q j +δ ̄q j ⟩ exactly (an identity that holds under any signal–noise correlation whatsoever) and the sole n = 2 usage, the Gaussian-limit convergence experiment, is run where the lemma applies. Every use of measured moments in this paper is therefore inside the lemma’s stated scope. Unbiasedness alone does not establish that finite-shot estimates suffice; the fol- lowing bound does, and quantifies the measurement cost of the moment readout. Applied to the experiments of this paper, the lemma yields the measurement bud- gets of Table S.4: every real-world benchmark operates the saturable device with a purely linear readout (n = 1), where the budget is the bare homodyne cost, while only the Gaussian-limit convergence experiment uses n = 2; a hypothetical Gaussian-limit implementation of the benchmark nonlinearity at n = 9 would pay a (2 · 9 − 1)!! ≈ 3.4 × 10 7 -fold factorial penalty. We stress that this number is comparator-relative: n = 9 is the readout order the Markovian atom–cavity compara- tor of main-text Fig. 3 requires on these tasks rather than an intrinsic scale of the problem, and the penalty should be read as the cost of matching that comparator in the Gaussian limit rather than as a property of the tasks themselves. In laboratory terms, at relative precision ε rel on the second moment (i.e. σ = ε rel ̄v) the requirement is N shots ≥ 3/ε 2 rel , independent of ̄v: 3× 10 4 repetitions at one-percent precision, or about 30 ms of wall clock per output sample at the ∼ μs symbol clock of the main text. Experimentorder n N shots penalty vs. n=1 Real-world benchmarks (saturable, linear)1 ̄v/σ 2 1 Convergence experiment (Gaussian limit)23 ̄v 2 /σ 2 3 ̄v Ninth-order comparator (Fig. 3), if Gaussian917!! ̄v 9 /σ 2 ∼ 3.4× 10 7 ̄v 8 Supplementary Table S.4 Shot budgets by experiment (Lemma 15), from N shots ≥ (2n− 1)!! ̄v n /σ 2 with ̄v ≡ v 0 + q 2 max the second-moment scale and σ the target precision. The saturable device’s nonlinearity transfer keeps all real-world tasks at the n = 1 cost. Proof. Var[ˆμ n ] = Var[I n ]/N shots ≤ E[I 2n ]/N shots ; conditioning on the bounded signal and applying Isserlis’ theorem to the centered Gaussian part, E[I 2n ]≤ (2n− 1)!! (v 0 + q 2 max ) n using E[ξ 2m ] = (2m− 1)!!v m 0 and the binomial expansion. 58 S.8 Prior Constructions Compared This table, referenced from the main text, locates the present result against prior reservoir-computing constructions. S.9 Benchmark Tasks We report no tuned, noise-matched comparison against an echo-state network; the echo-state figures that appear in the financial study below are a noiseless classical reference rather than a matched baseline, and are labelled as such. A matched compar- ison requires two things this manuscript does not supply: injection into the classical reservoir state of the shot noise implied by the device’s finite measurement bud- get, and independent hyperparameter selection for both families at each noise level. Without both, a noiseless and separately tuned classical reservoir is not a matched baseline, and a comparison against one is uninformative in either direction. The com- parators we do report are matched by construction: the atom–cavity comparator and the Markovian-limit ablation run through the same generator, the same readout and the same measurement statistics as the device, and they are what the nonlinearity- transfer claim rests on. The question a matched classical comparison would answer is also not the one this paper asks: the claim is minimality of the physical architecture rather than computational advantage, and the evidentiary standard for the real-world tasks is parity with standard classical methods. Financial forecasting study. The financial study is reported here rather than in the main text: its one-step-ahead structure on normalized prices is closely tracked by a persistence predictor, which makes it a weak discriminator between methods; we therefore treat it as a consistency check. We used daily closing prices for Apple, the S&P 500, and NASDAQ from Yahoo Finance over 1 January 2014 to 1 January 2024, normalized before injection and evaluated with a rolling two-year-train / one-year-test protocol. With ten input masks the reservoir achieved NRMSE values of 0.075± 0.058, 0.059± 0.025, and 0.079± 0.070 (mean± SD over rolling windows) for the S&P 500, Apple, and NASDAQ, against an echo-state network at 0.076± 0.066, 0.077± 0.048, and 0.072± 0.043; the differences lie well within one standard deviation, and we make no formal equivalence claim. Under the range normalization used throughout, the persistence baseline ˆy t = y t−1 scores ≈ 0.037 on these identical windows (all three series; computed on the stated protocol by the deposited script b6 persistencecheck.py), below both models: one-step-ahead prediction of slowly varying normalized prices is dominated by persistence, which is precisely why this task is retained as a consistency check between methods rather than as a benchmark, and no claim of beating persistence is made. Varying the delay over τ = 5, 10, 15, 20 gave 0.042± 0.008, 0.047± 0.011, 0.075± 0.058, and 0.101± 0.102, indicating that a moderate memory length is preferable and that excessively long feedback retains stale information. 59 Trivial baselines. For the tumor task of the main text the 48-sample FFPE set is near-evenly split between tumor and normal, so a majority-class classifier scores near chance—well below the reported 87.5% (42/48). The feature subsets were selected on the full sample by the classical pipeline of the main text’s Methods (a protocol optimized for, and shared with, the LDA baseline), so the selection step is common to both classifiers and the comparison between them is on equal footing; the Wilson 95% interval on 42/48 is [75.3%, 94.1%]. Because that subset search is performed on the full sample, the absolute accuracies of both classifiers, and the interval, carry selection bias and should not be read as generalization estimates; the reservoir–LDA comparison, whose selection step is shared, is the valid object. The paired comparison is moreover exact without any per-sample data: both classifiers scoring 42 of 48 forces the discordant pairs into balance (42 = n both + b and 42 = n both + c give b = c), and the exact two-sided McNemar test at b = c yields p = 1.0 identically—no detectable difference at this sample size, whatever the discordant count. For the speech task, no paired significance is attainable at the reported design as a matter of arithmetic: with five paired folds, the exact two-sided Wilcoxon signed-rank test’s smallest achievable p- value is 2/2 5 = 0.0625, so the fold count itself precludes significance at the conventional level regardless of the per-fold values; the deposited standby script bpairedstats.py is the protocol of record for any expanded-fold analysis. S.10 Robustness to Dephasing, Phase Jitter, and Feedback Loss This section quantifies the sensitivity of task performance to the three imperfec- tions most relevant to hardware: pure dephasing of the atom, jitter of the round-trip phase, and loss in the delayed-feedback path. Protocol common to all three sweeps: NARMA10 at the main-text parameters (γ = 0.1, τ = 10, ε INPUT = 0.1, φ = π/3 unless stated); noise simulated by stochastic unraveling of pure-state trajectories; one “node” denotes one delayed quadrature feature of the single measured output (succes- sive columns of the delay-embedded (P,Q) record; the readout dimension, as defined for Fig. 3 of the main text); features are averaged over 40 independent noise real- izations before the readout is trained, matching the ensemble averaging inherent to experimental moment estimation. Throughout this section NRMSE is normalized by the target range (max–min), the convention of the benchmark literature; uncertainty bands are 95% bootstrap confidence intervals obtained by resampling the 40 noise re- alizations with replacement before feature averaging and retraining. The three sweeps report slightly different noiseless reference errors—0.078 (dephasing), 0.0860 (jitter), and 0.084 at the loss-sweep optimum—because they were run with different train- ing and fading-memory washout points and independently configured node grids and noise-ensemble handling (the dephasing and jitter sweeps average over 40 stochas- tic trajectories, whereas the loss sweep uses single deterministic-input realizations, as noted below); the values are therefore internally consistent within each sweep but are not a single shared baseline, and comparisons are made within a sweep rather than 60 across the three reference numbers. Table S.6 collects the three reference points and their ensemble conventions in one place. Pure dephasing. Figure S.5 sweeps the dephasing rate over a logarithmically spaced grid γ φ ∈ 0, 10 −4 , 3× 10 −4 , 10 −3 , 3× 10 −3 , 10 −2 ,..., 10 −1 . Degradation is graceful, with no cliff: relative to the noiseless best error 0.078, dephasing at γ φ /γ = 10 −2 raises the error from 0.078 to 0.083, γ φ /γ = 10 −1 costs ∼25% (0.097), and at γ φ = γ the error reaches 0.158, roughly double this sweep’s own noiseless reference: full-rate dephasing erases the non-Markovian benefit. (The severed-feedback plateau is a quantity of the loss sweep and is not compared against here, per the within-sweep rule stated above.) Two structural features are visible in the node-resolved curves. First, node-scaling sat- urates progressively earlier as γ φ grows, consistent with the mode-resolution picture: dephasing broadens the dressed lines, and once γ φ exceeds the relevant splittings, addi- tional measured nodes carry no new information, giving a maximum useful node count that shrinks with γ φ . Second, the step structure of the noiseless curve (in particular the drop near 13 nodes) survives at small γ φ and washes out at large γ φ , as expected for features that live in inter-mode coherences; whether these steps admit the same polynomial-degree-boundary interpretation as the staircase of main-text Fig. 3 is sug- gestive but not exact. For the platforms discussed in the main text, atomic dephasing ratios of γ φ /γ ∼ 10 −3 –10 −2 are routine, placing them in the ≲10% penalty regime. Round-trip phase jitter. Figure S.6 adds Gaussian jitter of standard deviation δ (radians) to the round-trip phase, drawn per trajectory (modeling slow interferometric drift of the atom–mirror distance, i.e. drift on timescales long compared with a round trip; fast intra-round- trip phase noise is outside the simulated regime), together with two static anchors: the fully constructive phase φ = 0 and the fully destructive phase φ = π. Performance is statistically flat through δ = 0.5 rad—the 95% bootstrap CIs overlap the baseline across that entire range, and the shallow minimum at δ = 0.2 (0.0855 against baseline 0.0860) is well within CI, so we do not interpret it—and loses ∼10% at δ = 1 (0.095). The static anchors carry the theoretical content: the performance floor is the destruc- tive phase φ = π—precisely the dark-state configuration at which the delay-stability condition D(c 0 ) of Theorem 1 becomes marginal (Sec. S.4.1)—while the constructive phase φ = 0, which maximizes the collective decay and thus shortens the memory, is also inferior to the intermediate operating point φ = π/3. The task-level simulation thus reproduces, dynamically, the stability landscape derived analytically from the loop characteristic equation. As an engineering statement: interferometric stability of the atom–mirror path to ∼0.5 rad suffices (∼ λ/13 of the round-trip path 2L, equiv- alently ∼ λ/25 of the one-way atom–mirror distance L), a comfortable requirement on the platforms considered. The per-trajectory (slow-drift) model is the appropriate regime for those platforms: the round trip is τ = 2L/v—nanoseconds or less for the superconducting transmission lines and trapped-atom geometries cited in the main text—whereas the dominant phase noise in both, thermal and mechanical drift of the atom–mirror path length, occurs on millisecond-and-slower timescales. The round-trip 61 Supplementary Figure S.5 Pure dephasing. Top: NARMA10 test NRMSE versus number of measured nodes for logarithmically spaced dephasing rates (color bar); shaded bands, 95% bootstrap CI over the 40 noise realizations. The γ φ = 0 baseline (dashed) is deterministic and carries no band. Bottom: best (minimum-over-nodes) NRMSE versus γ φ on a logarithmic axis, with bootstrap CIs; top axis, the ratio γ φ /γ. Degradation is graceful, and the useful node count shrinks as dephasing broadens the dressed modes. phase is therefore effectively static within any single trajectory and varies between them, which is exactly the sampling this sweep implements; phase noise with appre- ciable spectral weight at the round-trip frequency would require a separate treatment and is not claimed here. Feedback loss. Figure S.7 attenuates the delayed-feedback path, with η the retained power fraction of the returning field (η = 0 severs the loop entirely). Severing the feedback collapses performance to a high plateau (0.157, the level of the feedback-free device), directly executing the falsifier stated in the main text: if performance were insensitive to the feedback, the trained readout rather than the physics would be doing the computation. It is not—and the recovery is strikingly nonlinear in η: restoring just η = 0.01 (a 20 dB round-trip loss) already recovers most of the severed-to-optimum gap (best error 0.099, against 0.084 at the optimum and 0.157 severed; recovering 0.157− 0.099 of 62 Supplementary Figure S.6 Round-trip phase jitter. Top: test NRMSE versus measured nodes for jitter widths δ (per-trajectory Gaussian draws; bands, 95% bootstrap CI over 40 realizations), with the static constructive (φ = 0) and destructive (φ = π) anchors. Bottom: best NRMSE versus δ (log axis) with bootstrap CIs against the baseline and static-anchor reference lines. The performance floor is the dark-state phase φ = π predicted by the delay-stability analysis; jitter is tolerated to at least δ ≈ 0.5 rad, the CIs overlapping the baseline throughout that range. the 0.157− 0.084 severed-to-optimum gap, a comparison we quote only as an order- of-magnitude indication because this sweep, unlike the dephasing and jitter sweeps, uses single deterministic-input realizations and therefore carries no seed spread or confidence interval), returns saturate beyond η ≈ 0.1, and the shallow minimum near η ≈ 0.2–0.4 visible in the nine-point grid is left uninterpreted, these being single deterministic-input realizations. The non-Markovian resource is therefore not only necessary for the performance but remarkably loss-tolerant: even one percent of the returning power carries most of the computational benefit, which substantially relaxes the insertion-loss budget of any experimental implementation. 63 Supplementary Figure S.7 Feedback loss. Left: test NRMSE versus measured nodes as the re- tained feedback power fraction η is varied (deterministic-input runs, one realization per η; no noise ensemble, hence no bands). Right: best NRMSE versus η over the full nine-point grid. η = 0 (severed loop) collapses to the feedback-free plateau; η = 0.01 (a 20 dB round-trip loss) already recovers most of the severed-to-optimum gap, returns saturate beyond η ≈ 0.1, and the shallow minimum near η ≈ 0.2–0.4 is left uninterpreted (single realizations). S.11 Numerical Verification of the Universality Properties Figure S.8 reports the three classical sufficient conditions for reservoir universality (separation, fading memory, and polynomial enrichment) for the non-Markovian reser- voir. They are an independent sanity check rather than ingredients of the proof of Theorem 1, which rests instead on the Volterra expansion, the resolution of the eigen- value sums, the Vandermonde inversion, the extrapolation identity and the tail bound; fading memory of the device is established separately in Sec. S.4. Separation is shown by applying small perturbations to the entire input: even minute perturbations yield distinguishable output trajectories, so nearby input histories map to distinct reser- voir states. Fading memory is shown by perturbing a single input point near t = 100: the response rises then decays, so recent inputs dominate over distant ones. Polyno- mial enrichment is shown by forming nonlinear combinations of collected states, which enlarge the readout space and improve approximation with saturating gains as the node count grows. These are numerical evidence for the ingredients that the proof establishes analytically in the linear-transducer limit. S.12 Supplementary Discussion: the Stone–Weierstrass route, and what the atom contributes This section carries, at full length, two discussions summarized in the main text: why the theorem is proven by the constructive Volterra-kernel route although the abstract density route is also available, and the strongest form of the classical-linear-optics objection together with our answer. 64 Supplementary Figure S.8 Numerical demonstration of separation, fading memory, and polyno- mial enrichment. Parameters: ∆t = 1, γ = 0.1, ε = 0.15, θ = φ = π/3, τ = 10. We emphasize that universality itself is also within reach of the standard, shorter route: the reservoir outputs contain the constants, form an algebra—sums are realized by adding weights, and products of outputs are outputs of higher readout order, under the same spectral non-degeneracy the theorem’s third condition supplies—and, ranged over the nested family of operating points, they separate distinct input histories, so the Stone–Weierstrass theorem delivers density in the continuous fading-memory func- tionals on K u max , exactly as in prior universality results for reservoir classes [11, 16]. We prove Theorem 1 of the main text constructively not because the abstract route fails here, but because of what it cannot say. The Volterra-kernel route provides this where a Stone–Weierstrass argument cannot: it certifies that scaling the one device up is never wasted, rather than merely asserting that some adequate member of a class exists—and it does so with explicit rates and constants, which no density argument produces. The containment is what is proven outright; the strictness is generic in the same sense as the theorem’s conditions. Separation, fading memory, and polynomial enrichment are additionally verified numerically for the simulated device in Sec. S.11. We are careful about what does the work in the proof, since it is easy to overstate. The essential ingredient is a linear multimode structure with generically non-resonant frequencies and a tunable readout—and linear multimode structure is not the ex- clusive property of a quantized field: classical wave optics supplies it too, in the independent spatial or spectral modes of a linear optical network. The relevant di- chotomy is therefore not quantum-versus-classical but nonlinear single node with virtual (time-multiplexed) nodes versus linear field with genuinely independent modes. The single-node delay-line reservoir [20, 21], the atom–mirror system’s closest archi- tectural relative, sits on the first side: its virtual nodes are time-multiplexed samples of one nonlinear trajectory rather than independent degrees of freedom, which is pre- cisely why it has resisted a universality theorem. The atom–mirror device sits on the second: the theorem becomes available for this instance because its delay loop decom- poses into independent linear modes—the structure the classical device’s nonlinear node destroys (the mode-space section of the main text). What quantization con- tributes specifically, in one atom, is the packaging of three ingredients into a single passive component: input encoding on the atomic drive, the mirror-selected indepen- dent mode structure that carries the recurrence, and (beyond the Gaussian limit) the atom’s saturable nonlinearity, which moves the nonlinearity into the hardware and collapses the Gaussian readout’s n = 9 shot budget to n = 1 (Sec. S.7). We prove the 65 theorem for the quantized instance and identify why the classical delay-line original cannot inherit it; we do not claim quantization is the only route to a linear multimode reservoir. That concession invites the sharpest form of the objection, and we state it in full rather than leave it implicit. If linear multimode structure is what the proof needs, and classical linear optics supplies it, then in the regime where our theorem holds (the Gaussian limit, where the atom is a linear transducer) a classical linear optical network with a comparable mode structure would appear to be an equally good instance; while in the regime where the atom’s own nonlinearity matters, we have no theorem. Our answer is not that the atom is doing something a classical field cannot, because in the Gaussian limit it demonstrably is not. The answer is that the two regimes are the same device at two settings of one dial, and that the claim being made is about packaging rather than about quantum supremacy of any kind. A classical linear multimode network can host the recurrence, but it must import its nonlinearity and its input encoding from separate components—a modulator, a nonlinear element, or a polynomial readout paid for in the factorial shot budget of Sec. S.7. The atom– mirror device carries all three in one passive object: the drive encodes, the mirror- selected modes recur, and the same atom that acts as a linear transducer at small drive becomes the saturable nonlinearity at larger drive, continuously and with a bounded, computable gap between the two (Sec. S.5). What the theorem certifies is that this packaging is not paid for in expressivity: at the setting where the device is simplest to analyze, it is already universal. The minimality claim is therefore about component count for a fixed capability, and it survives the concession intact—but it is a claim about one atom replacing an assembly rather than about quantum mechanics enabling a computation classical optics could not perform. Supplementary References [1] Gardiner, C.W., Collett, M.J.: Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation. Physical Review A 31(6), 3761–3774 (1985) [2] Combes, J., Kerckhoff, J., Sarovar, M.: The SLH framework for modeling quantum input–output networks. Advances in Physics: X 2(3), 784–888 (2017) [3] Cuchiero, C., Gonon, L., Grigoryeva, L., Ortega, J.-P., Teichmann, J.: Discrete- time signatures and randomness in reservoir computing. IEEE Transactions on Neural Networks and Learning Systems 33(11), 6321–6330 (2022) https://doi. org/10.1109/TNNLS.2021.3076777 [4] Boyd, S., Chua, L.: Fading memory and the problem of approximating nonlin- ear operators with Volterra series. IEEE Transactions on Circuits and Systems 32(11), 1150–1161 (1985) https://doi.org/10.1109/TCS.1985.1085649 [5] Hale, J.K., Verduyn Lunel, S.M.: Introduction to Functional Differential Equations. Applied Mathematical Sciences, vol. 99. Springer, New York (1993) 66 [6] Kato, T.: Perturbation Theory for Linear Operators. Classics in Mathematics. Springer, Berlin, Heidelberg (1995). Reprint of the 1980 second edition; see Ch. I §1 and Ch. VII§1 for the analyticity of Riesz projections under perturbation [7] Grigoryeva, L., Henriques, J., Larger, L., Ortega, J.-P.: Stochastic nonlinear time series forecasting using time-delay reservoir computers: performance and univer- sality. Neural Networks 55, 59–71 (2014) https://doi.org/10.1016/j.neunet.2014. 03.004 [8] Grigoryeva, L., Henriques, J., Larger, L., Ortega, J.-P.: Optimal nonlinear infor- mation processing capacity in delay-based reservoir computers. Scientific Reports 5, 12858 (2015) https://doi.org/10.1038/srep12858 [9] K ̈oster, F., Yanchuk, S., L ̈udge, K.: Master memory function for delay-based reservoir computers with single-variable dynamics. IEEE Transactions on Neu- ral Networks and Learning Systems 35(6), 7712–7725 (2024) https://doi.org/10. 1109/TNNLS.2022.3220532 . Preprint: arXiv:2108.12643 [10] Ort ́ın, S., Soriano, M.C., Pesquera, L., Brunner, D., San-Mart ́ın, D., Fischer, I., Mirasso, C.R., Guti ́errez, J.M.: A unified framework for reservoir computing and extreme learning machines based on a single time-delayed neuron. Scientific Reports 5, 14945 (2015) https://doi.org/10.1038/srep14945 [11] Grigoryeva, L., Ortega, J.-P.: Echo state networks are universal. Neural Networks 108, 495–508 (2018) https://doi.org/10.1016/j.neunet.2018.08.025 [12] Gonon, L., Ortega, J.-P.: Reservoir computing universality with stochastic inputs. IEEE Transactions on Neural Networks and Learning Systems 31(1), 100–112 (2020) https://doi.org/10.1109/TNNLS.2019.2899649 [13] Stelzer, F., R ̈ohm, A., Vicente, R., Fischer, I., Yanchuk, S.: Deep neural net- works using a single neuron: folded-in-time architecture using feedback-modulated delay loops. Nature Communications 12, 5164 (2021) https://doi.org/10.1038/ s41467-021-25427-4 [14] Gonon, L., Ortega, J.-P.: Fading memory echo state networks are universal. Neural Networks 138, 10–13 (2021) https://doi.org/10.1016/j.neunet.2021.01. 025 [15] Grigoryeva, L., Ortega, J.-P.: Universal discrete-time reservoir computers with stochastic inputs and linear readouts using non-homogeneous state-affine systems. Journal of Machine Learning Research 19(24), 1–40 (2018) [16] Chen, J., Nurdin, H.I.: Learning nonlinear input–output maps with dissipative quantum systems. Quantum Information Processing 18(7), 198 (2019) [17] Chen, J., Nurdin, H.I., Yamamoto, N.: Temporal information processing on noisy 67 quantum computers. Physical Review Applied 14(2), 024065 (2020) https://doi. org/10.1103/PhysRevApplied.14.024065 [18] Nokkala, J., Mart ́ınez-Pe ̃na, R., Giorgi, G.L., Parigi, V., Soriano, M.C., Zambrini, R.: Gaussian states of continuous-variable quantum systems provide universal and versatile reservoir computing. Communications Physics 4(1), 53 (2021) [19] Gauthier, D.J., Bollt, E., Griffith, A., Barbosa, W.A.: Next generation reservoir computing. Nature communications 12(1), 5564 (2021) [20] Appeltant, L., Soriano, M.C., Sande, G., Danckaert, J., Massar, S., Dambre, J., Schrauwen, B., Mirasso, C.R., Fischer, I.: Information processing using a single dynamical node as complex system. Nature Communications 2, 468 (2011) https: //doi.org/10.1038/ncomms1476 [21] Larger, L., Bayl ́on-Fuentes, A., Martinenghi, R., Udaltsov, V.S., Chembo, Y.K., Jacquot, M.: High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification. Physical Review X 7(1), 011015 (2017) https://doi.org/10.1103/PhysRevX.7.011015 68 Supplementary Table S.5 Reservoir-computing constructions compared on the properties this work establishes. “Physical reservoir” asks whether the reservoir is realized as hardware dynamics (with its active component count) or computed digitally. “Single-member universal?” asks whether one fixed system carries a universal-approximation guarantee: “No (class)” means the guarantee attaches only to a class ranged over, never to an individual member. “Strict scaling” asks whether increasing a single physical resource of one device is proven never to reduce its approximation capability, with the span of exactly matchable kernels strictly enlarged; a “No” records that no such statement has been made for that construction rather than that it fails—for classical network classes containment across sizes is trivially available by zero-padding weights, and what is proven here is containment for a non-nested spectral family, which zero-padding does not supply. The universality entry for this work holds under the constructive conditions of the main text’s universality theorem (Supplement Sec. S.2.5), all proven satisfiable with explicit constants; delay stability of the feedback loop, verified with margin c 0 ≥ 0.36γ throughout the simulated regime, is separately the condition for the device’s fading memory. The “No” entry for the single-node delay-line reservoir reflects a literature search across the delay-reservoir theory line: capacity analyses (including the master memory function) quantify memory rather than universal approximation [7–9]; the unified single-neuron framework of Ort ́ın, Pesquera and co-workers shows the delay node implements echo-state and extreme-learning schemes [10], whose universality is class-level [11, 12], so the equivalence transfers a class guarantee rather than a single-member one; and folded-in-time architectures inherit the guarantees of the network they emulate [13]. For the classical network classes, capability containment across sizes is trivially available by zero-padding weights; the “No” entries record that no such statement is proven in the cited works, and the nontrivial content of the present entry is containment for a family whose spectra are not nested, where zero-padding is unavailable (Sec. S.2.5). The search underlying the delay-line entries is current at the time of submission and includes the folded-in-time and deep delay-architecture line. The single-member and strict-scaling entries for this work hold on the operating class of Definition S.1, whose conditions a fabricated device can be checked against; the constructive certificate of Sec. S.2.5.1 supplies them by design and additionally specifies the resonator spectrum. ConstructionPhysical reservoir (components) Single-member universal? Strict scaling? Echo-state network classes [11, 14] Realization (O(nodes))No (class)No State-affine system classes [15] Abstract (O(dimension))No (class)No Ising / spin-ensemble quantum reservoir classes [16, 17] Spin ensemble (O(spins))No (class)No Gaussian quantum-optics classes [18] Optical network (O(modes))No (class)No Boyd–Chua / Volterra density [4] Abstract (functional class)N/A (existence only)No Signature / next-generation RC [3, 19] No (digital feature map)YesNo Single-node delay-line RC [20, 21] Yes (1 node + 1 delay loop)NoNo Dissipative quantum systems [16] YesNo (class of systems)No This workYes (1 atom, 1 mirror, 1 detector, 1 fixed band envelope) YesYes 69 SweepNoiseless ref.EnsembleUncertainty Dephasing (Fig. S.5)0.07840 noise realizations95% bootstrap CI Phase jitter (Fig. S.6)0.086040 per-trajectory draws95% bootstrap CI Feedback loss (Fig. S.7) 0.084 (optimum) single deterministic run none (no seed spread) Supplementary Table S.6 Reference baselines of the three robustness sweeps. The three noiseless references differ because the sweeps were run with different training and fading-memory washout points and independently configured node grids and ensemble handling; comparisons are made within a sweep only. The absence of a confidence interval on the loss sweep is a stated limitation of that sweep rather than of the other two. 70