Efficient Non-Uniform Quantum Hermite Transform through Adaptive Sampling Nitay Mayo1 and Aryeh Lev Zabokritskiy (Yohananov)∗2
arXiv:2609.20739v1 [quant-ph] 17 Sep 2026
1
2
Department of Computer Science, Technion – Israel Institute of Technology, Haifa, Israel Department of Computer Science, Tel-Hai University of Kiryat Shmona in the Galilee; MIGAL – Galilee Research Institute
Abstract On the span of the first N oscillator modes, Gauss–Hermite quadrature gives an exact change of basis between mode coefficients and N weighted position space samples. We implement this transform with O(N polylog(N, 1/ε)) logical gates and polylogarithmic quantum width. The operator-error bound ε holds on arbitrary superpositions and includes all auxiliary registers. The construction uses signed averages on adaptive windows to convert uniform-grid samples into weighted Hermite-root samples. Their varying widths control the amplification cost, giving the near-linear bound.
1
Introduction
Gaussian quadrature allows certain continuous integrals to be evaluated exactly from finitely many samples. For arbitrary combinations of the first N oscillator modes, N Hermite roots and their quadrature weights preserve inner products and recover every continuous mode coefficient exactly. Thus the finite sample representation introduces no discretization error on this space. The question here is the quantum cost of this change of basis. It also provides a finite representation of harmonicoscillator evolution and fractional Fourier rotations, both diagonal in the oscillator basis. Golub and Welsch give the classical eigenvalue construction of Gaussian quadrature [5]. Pliś and Zak implement the finite Hermite transform on a quantum computer [13]; Webb and Maierhofer study stable classical construction of the normalized Hermite matrix [16]. Our contribution is a near-linear gate bound for the same finite transform. We convert an existing uniform-grid efficient Hermite sampler [6] by using signed averages on disjoint adaptive windows. Adapting their widths to the local root spacing controls normalization even at the outermost roots. Our conventions are as follows. Let N = 2n , where n ≥ 1 is the number of data qubits. We denote the discrete Bargmann transform from weighted root samples to oscillator coefficients by B and construct its inverse B † . Circuit errors are operator-norm errors that include every auxiliary register, with fixed phase. The input is an arbitrary quantum superposition; loading classical data into its amplitudes and reading out all coefficients are separate tasks. ∗
ORCID: 0000-0003-3151-6192
1
1.1
The finite transform and its input
We define the oscillator number basis in L2 (R) using the normalized Hermite functions: Hk (x) 2 φk (x) = e−x /2 q √ , 2k k! π
k ≥ 0,
(1)
where Hk are the physicists’ Hermite polynomials (with H0 (x) = 1). The corresponding finite subspace is HN = span{φ0 , . . . , φN −1 }. While continuous states can be formally mapped to a phase-space representation via the standard Bargmann–Fock continuous kernel [2], B outputs number-basis coefficients rather than values on a complex grid; the circuit and its reverse implement the two directions of this finite change of basis. Let x1 < · · · < xN be the roots of HN , with corresponding Gauss–Hermite quadrature weights [12]: √ 2N −1 N ! π wi = 2 > 0. (2) N HN −1 (xi )2 For a polynomial p of degree at most 2N − 1, this quadrature exactly evaluates the Gaussian integral: Z
2
e−x p(x) dx =
R
N X
wi p(xi ).
(3)
i=1
Define the elements of the transform matrix B as [13, 16]: √ wi Bk,i = q 0 ≤ k < N, √ Hk (xi ), 2k k! π
1 ≤ i ≤ N.
(4)
Rows are number-basis indices; column i is stored under computational label i − 1. Applying the exact quadrature to Hk Hm , of degree at most 2N − 2, yields orthonormality [12, Sec. 3.5(v)]: BB † = B † B = I.
(5)
The sampling convention is essential. For a continuous function with defined point values, we encode the amplitudes as weighted samples: (SN ψ)i = For any state ψ = and norms:
√
2
wi exi /2 ψ(xi ).
(6)
PN −1
k=0 ck φk ∈ HN , these weights ensure the exact recovery of both coefficients
SN ψ = B †⃗c,
BSN ψ = ⃗c,
∥SN ψ∥2 = ∥ψ∥L2 .
(7)
These identities express the core quadrature advantage: the encoded amplitudes are SN ψ, not the unweighted values ψ(xi ). The corresponding states |xi ⟩ form an orthonormal finite position basis, rather than a set of continuous position eigenstates. For inputs outside HN , truncation and sampled-tail aliasing contribute additional errors.
2
2
Main result
For the main result take N = 2n ≥ 2 and 0 < ε ≤ 1/2, and let the input and output data registers have n qubits. Let Ein and Eout denote their isometric embeddings into the full circuit space with every other register initialized to zero. Fixed register permutations identify their data labels with the mode and root labels, respectively. Theorem 1 (Near-linear implementation). There is a Clifford+T quantum circuit UN,ε such that UN,ε Ein − Eout B † ≤ ε.
(8)
The circuit uses O N polylog(N, 1/ε) logical gates, hence at most this depth, and polylog(N, 1/ε) qubits including workspace. Its reverse implements B with the same error and resources.
Here polylog(N, 1/ε) denotes a fixed polynomial in 1 + log N + log(1/ε). The bound is for quantum execution after classical compilation; it excludes arbitrary input loading and full classical readout. For N = 1 the transform is the identity. Supplementary Section S7 gives the separate continuous-input error bounds. While Theorem 1 establishes the headline complexity bound for the transform, the main construction is the efficient averaging window operator R. Generic synthesis gives a quadratic gate upper bound, but does not exploit the spacing and weights of the Hermite roots. The operator R maps uniform-grid samples to the non-uniform weighted root samples. The near-linear blockencoding of this window operator enables the efficient implementation of the transform B. The uniform-grid transform used in the proof is based on [6]; its precise sampling and phase guarantees are proved in Supplementary Theorem S1.
3
Applications and Use-Cases
The transform B connects weighted position samples to the Hermite basis, in which both harmonic evolution and fractional Fourier transformation are diagonal. We describe these two applications. The exact identities below concern normalized states in the finite space HN ; implementing them requires the finite-precision circuits of Theorem 1.
3.1
Fast-Forwarding the Quantum Harmonic Oscillator
In oscillator units, the quantum harmonic oscillator has Hamiltonian Hosc = (−d2 /dx2 + x2 )/2 and eigenvalues Ek = k + 1/2. Its restriction to HN can be evolved directly for any time t, without a time-stepping or product-formula approximation. By (7), the corresponding operator on weighted samples is N −1 UX (t) = B † D(t)B, D(t) = diag e−i(k+1/2)t k=0 . (9) To fast-forward an input state by an arbitrary time t, we map the system to the energy basis via B, execute a diagonal matrix of phase rotations corresponding to the exact QHO energy eigenvalues Ek = k + 1/2, and return to the spatial basis using B † . The exact finite-space time-evolution sequence is: |SN ψ(0)⟩
n
B
D(t)
B†
|SN ψ(t)⟩
(10)
The diagonal phases factor over the binary digits of k, together with the global phase e−it/2 . Thus the number of phase rotations need not grow with the number of elapsed oscillator periods; their 3
precision and classical angle calculation remain part of the implementation. Figure 1(a) illustrates the coherent-state trajectory. A nontrivial coherent state has infinitely many Hermite coefficients, so its use in HN requires truncation. The exact projected evolution preserves that truncation error in L2 ; transform synthesis and phase synthesis contribute separate errors. Both plots below use normalized projections onto N = 32 modes and reconstruct the continuous density |ψ(x)|2 from those modes. The omitted squared norms are approximately 1.46 × 10−18 for the coherent state and 6.87 × 10−17 for the cat state. The accompanying reproducibility package supplies the coefficients, plotted data, and generating code. Coherent-state Hermite projection
π Fractional angle (radians)
0.6
0.4 0.3 0.2
π
0.1 0.05
Density bands
Oscillator time t
0.5
0.025 0.01
0 -6
1 0.8 0.55 0.35
π/2
0.2 0.1 0.05 0.025
0 -6
0
-3 0 3 Position coordinate x
1.2
6
(a) Coherent state with real amplitude 2, evolved for 0 ≤ t√ ≤ 2π in oscillator units. The dashed curve is x = 2 2 cos t.
Density bands
2π
Even-cat Hermite projection
0
-3 0 3 Position coordinate x
6
(b) Even cat state with displacement x0 = 3, rotated through 0 ≤ α ≤ π. The momentum-domain fringes occur at α = π/2; the even position profile returns at α = π.
Figure 1: Visual demonstrations of the Quantum Hermite Transform applied to physical timeevolution (left) and continuous-variable signal processing (right). Filled bands show probabilitydensity intervals. The finite-transform identities apply to the retained Hermite modes; these illustrations are not gate-level or hardware error benchmarks.
3.2
Continuous Fractional Fourier Analysis
In signal processing, the standard quantum Fourier transform acts on a cyclic group, whereas the Hermite representation describes functions on the real line. We use the continuous Fourier R convention (F ψ)(p) = (2π)−1/2 R e−ipx ψ(x) dx. Its fractional extension is defined spectrally by Fα φk = e−ikα φk , with F0 = I and Fπ/2 = F [8, Secs. 1–3]. The fractional Fourier transform (FrFT) rotates a signal in phase space by an angle α. On −1 HN , its exact representation on weighted samples is B † diag(e−ikα )N k=0 B. Writing ψα = Fα ψ0 , this gives the circuit |SN ψ0 ⟩
n
B
diag(e−ikα )
B†
|SN ψα ⟩
(11)
The QHO evolution differs from this transform by the global phase e−it/2 when α = t [8, Eq. (16)]. Both operators preserve HN and the norm. The example in Figure 1(b) uses a bimodal Schrödinger cat state. In the position domain (α = 0), the state is described by a superposition of displaced coherent states, with real displacement 4
x0 :
1
1
2
2
ψ0 (x) ∝ e− 2 (x−x0 ) + e− 2 (x+x0 ) .
(12)
As the signal is rotated into the momentum domain (α = π/2), the Fourier transform of these displaced spatial peaks yields the interference pattern 1 2
ψ̃(p) = ψπ/2 (p) ∝ e− 2 p cos(px0 ).
(13)
These are continuous formulas. Their normalized Hermite projections yield finite-register examples, with the discarded tail and any sampled-tail aliasing treated separately. The finite operator preserves the projected norm exactly and imposes no periodic spatial boundary. Supplementary Section S8 gives a small, reproducible illustration of qubit noise viewed through Hermite-mode populations. The transform specifies the preparation and measurement bases for that observable. We pose this mode-population diagnostic as an open direction for future research to further formalize.
4
Averaging Operator R on a Uniform Grid
To implement the weighted-root sampling map B † in (6), we average uniform-grid samples near each Hermite root. Polynomial reproduction controls the error on the retained Hermite space; radii at the local root-spacing scale bound the averaging norm uniformly in N for fixed kernel degree and hence control coherent amplification. We prove Theorem 1 using the Introduction’s notation and error conventions. While coherent averaging and reversible changes between uniform sampling grids were pioneered by Kitaev and Webb [7], their construction does not map to non-uniform, weighted quadratures. Here, we extend and realize the standard polynomial kernel and amplitude amplification specifically to solve the non-uniform weighted Gauss–Hermite sampling problem.
4.1
Local weights, separation, and analytic control
The averaging row norm depends on the effective weight divided by its radius, so we need a local weight estimate. Define 2
vi = wi exi ,
qi = 2N − x2i + N 1/3 ,
c=
1 , 64
c ri = √ , qi
1 ≤ i ≤ N.
(14)
Here, vi is the effective weight, qi defines the local density scale tracking spatial root compression, c is a geometric safety constant guaranteeing disjoint windows, and ri sets the dynamically scaled safe √ radius for the averaging window. In particular (B † )i,k = vi φk (xi ), with the root-label convention above. A radius used below is always positive. Lemma 2 (Local weight and root separation). There are absolute constants cW , CW > 0 such that cW CW √ ≤ vi ≤ √ . qi qi
N 1/3 ≤ qi ≤ 3N,
(15)
Every consecutive root gap gi = xi+1 − xi satisfies −1/2
gi ≥ max{qi
−1/2
, qi+1 },
5
ri + ri+1 ≤ 2cgi .
(16)
The proof is in Supplementary Section S1.1. For parameter selection, one may use certified root and weight intervals to obtain a dyadic number √ C W ≥ max{1, max vi qi }. (17) i
Refining the intervals makes this bound a fixed-factor estimate of the maximum. Lemma 2 then bounds C W uniformly without a numerical source constant. Hermite functions are entire, but the approximation requires a bound uniform in the degree and the root, including the edges. The next lemma provides that bound on a disk scaled to ri . Cauchy’s estimate will then control all Taylor remainders with one choice of the kernel degree. Lemma 3 (Analytic bound at the local scale). Let ci ∈ R. Suppose |ci − xi | ≤ h ≤ ri and 0 < ρi ≤ ri . There is an absolute constant CM > 0 such that, for every root and every complex z with |z| ≤ 2, N −1 X √ |φk (ci + ρi z)|2 ≤ CM qi . (18) k=0
√ Consequently, there exists an absolute constant B0 > 0 such that the row vector fi (z) = ( vi φk (ci + N −1 ρi z)) √ k=0 has norm at most B0 . The analytic matrix with rows fi (z) has operator norm at most B0 N on this disk. The proof is in Supplementary Section S1.2.
4.2
The Averaging Operator R
A reproducing kernel turns uniform-grid samples into root samples by reproducing the initial Taylor terms; midpoint quadrature controls the discretization error. The Legendre reproducing kernel supplies the required polynomial weights. Let Pj be the Legendre polynomial with Pj (1) = 1. For an integer d ≥ 0 define Kd (u) =
d X
(2j + 1)Pj (0)Pj (u),
C0 = (d + 1)2 ,
C1 = (d + 1)4 .
(19)
j=0
The orthogonal-polynomial reproducing-kernel identity specializes to 1 2
Z 1 −1
Kd (u)p(u) du = p(0),
Indeed, expand p in the Legendre basis and use and |Pj′ (u)| ≤ j(j + 1)/2 on [−1, 1] give
IL = [−L/2, L/2),
(20)
R1
−1 Pj Pk = 2δj,k /(2j + 1). The bounds |Pj (u)| ≤ 1
∥Kd′ ∥∞ ≤ C1 .
∥Kd ∥∞ ≤ C0 , Let L be a power of two, put h =
deg p ≤ d.
(21)
2π/L, and define the L × N uniform sampling matrix by √ yℓ = hℓ, ℓ ∈ IL ∩ Z, JL [ℓ, k] = h φk (yℓ ). (22) p
A signed grid address ℓ is stored as its log2 L-bit residue, interpreted in [−L/2, L/2). All address arithmetic and interval tests below use this convention. Assuming the quantum grid has been made
6
sufficiently fine, for example h ≤ mini ri /16, choose an even power of two mi and an integer start address ai such that mi h ρi = ∈ (ri /4, ri ], 2
mi − 1 ci = h ai + , 2
|ci − xi | ≤ h
(23)
where ρi is the averaging window radius, and ci is the window midpoint which must be chosen to be sufficiently close to the i-th root xi . For each root i let the mi coordinates inside the i’th window be 2j + 1 − mi , 0 ≤ j < mi . (24) ui,j = mi Then ci + ρi ui,j = h(ai + j) exactly. The physical cells [ci − ρi , ci + ρi ] are disjoint by Lemma 2. Indeed, for an adjacent gap gi , the two radii and center displacements total at most (2c+c/8)gi < gi . The grid will also be chosen large enough to contain every window. The continuous average and its midpoint discretization are, respectively, √ Z 1 vi Tcont [i, k] = Kd (u)φk (ci + ρi u) du, (25) 2 −1 √ mX vi i −1 T [i, k] = Kd (ui,j )φk (ci + ρi ui,j ). (26) mi j=0 The matrix T approximates B † ; the following proposition bounds its error before circuit synthesis. Proposition 4 (Uniform approximation by adaptive windows). Under (23), define √ √ √ 3B0 2 αd = B0 C0 N = O( N (d + 1) ), αh = (1 + C0 + C1 )N. c
(27)
Then ∥T − B † ∥ ≤ Ed (h) := αd 2−d + αh h.
(28)
Both coefficients depend on N and d, through C0 = (d + 1)2 and C1 = (d + 1)4 , but are independent of h and L. Section S1.3 bounds the three errors separately. The coefficients satisfy αd = √Supplementary 2 O( N (d + 1) ) and αh = O(N (d + 1)4 ), so logarithmic degree and inverse-polynomial grid spacing suffice. To seamlessly bridge our uniform grid samples with the continuous approximation, we define a rectangular averaging matrix R : CL → CN . This matrix acts locally: it is zero everywhere outside the defined windows, while inside it applies the reproducing kernel weights: √ R[i, ai + j] =
m −1
vi Kd (ui,j ) √ , mi h
i 1 X Di = Kd (ui,j )2 , mi j=0
(29)
where Di serves as the discrete normalization constant for the i-th root. Recall that the grid addresses map exactly to our local coordinates via yai +j = ci + ρi ui,j . This alignment allows R, when applied to the uniform grid sampler JL , to perfectly reproduce our
7
discrete approximation matrix T : (RJL )[i, k] =
mX i −1 √ j=0
vi Kd (ui,j ) √ √ hφk (yai +j ) mi h
(30)
√ =
i −1 vi mX Kd (ui,j )φk (ci + ρi ui,j ) mi j=0
= T [i, k]
(31) (32)
Because the physical windows do not overlap, the rows of R are strictly orthogonal. This allows for a direct formula for the overall operator norm: RR† = diag
vi Di 2ρi
N
,
∥R∥2 = max i
i=1
vi Di . 2ρi
(33)
The integral norm of the kernel gives a sharper normalization than its pointwise bound. Orthogonality yields Z 1 1 Kd (u)2 du = Kd (0), 2 −1 (34) 1 ≤ Kd (0) ≤ d + 1. r To see the upper bound, use Legendre parity and P2r (0) = (−1)r 2r r /4 . Telescoping gives Kd (0) =
2
2r r r (2r +1)2 2r r /16 for r = ⌊d/2⌋. The ratio of successive factors r /4 is (2r −1)/(2r), so induction
2
r gives 2r r /16 ≤ 1/(2r + 1); the lower bound is the j = 0 term. Since Kd2 has derivative bounded by 2C0 C1 , midpoint quadrature gives |Di − Kd (0)| ≤ √ C0 C1 /mi ≤ 1/4 whenever mi ≥ 4C0 C1 . Combining this estimate with ρi > ri /4 and vi qi ≤ C W gives 3 5 Kd (0) ≤ Di ≤ d + , 4 4 (35) 2C 5 W 2 d+ . ∥R∥ ≤ c 4
Choose A as the least power of two satisfying
s
2C W 5 A ≥ max 2, 2 d+ , c 4 √ A = O( d + 1).
(36)
This choice precedes the final grid size and ensures A ≥ 2∥R∥ on every admissible finer grid.
4.3
Realization of the Averaging Windows Operator R
To realize R, prepare its normalized row states with equal preliminary success amplitudes, permitting a common amplification sequence. For the certified disjoint dyadic windows above, define |χi ⟩ =
mX i −1
Kd (ui,j ) √ |j⟩, mi Di j=0
W |i − 1⟩ =
8
mX i −1
Kd (ui,j ) √ |ai + j⟩. m i Di j=0
(37)
Disjoint supports and the definition of Di make W an isometry. We construct its root-to-grid extension and use the reversed circuit to average grid samples. The construction below prepares the signed states, unloads their metadata, and erases the root label coherently. Lemma 5 (Charged preparation of the window isometry). W is an isometry with a reversible unitary extension. For 0 < ζ < 1/2, its extension and its inverse can be approximated to operator error at most ζ using O N poly(d, log N, log L, log(1/ζ)) logical gates and poly(d, log N, log L, log(1/ζ)) quantum width. The same bounds hold for a unitary whose designated input-output block is R/A. Proof. We start by showing a staged construction of the operator W , detailing how the system prepares local states, amplifies them, and reversibly uncomputes auxiliary data. To clarify the internal mechanics, each stage below tracks the exact evolution of the quantum state for a simplified running example: a single input root labeled |i − 1⟩root , accompanied by initially empty auxiliary registers for metadata, local indices, the physical grid, and a success flag. 1. Metadata Loading (a) Define mi = 2si and allocate a shared quantum register of size smax = maxi si to hold the local index j. (b) Load the hardwired starting address ai , length si , and normalization constant Di by sequentially scanning the N classical records. Each scan computes equality to a label, conditionally XORs its record, and clears the comparison flag; its O(N ) times polynomial record-length gate cost is included. (c) Output State: The root register acts as a control to pull classical data into the active metadata registers. |i − 1⟩root |ai , si , Di ⟩meta |0⟩local |0⟩grid |0⟩flag
(38)
2. Local State Preparation (a) Apply Hadamard gates to precisely the lowest si bits, leaving the higher bits zero, to generate a uniform superposition over valid local addresses 0 ≤ j < mi . (b) Reversibly compute the local weight, rotate an auxiliary flag qubit to encode the target amplitude bKd (ui,j ) fi (j) = √ (39) Di and immediately uncompute the arithmetic to clear the workspace. Set t0 = 4(d + 1) P and b = sin( 4t0π+2 ). Cauchy–Schwarz gives |Kd (u)|2 ≤ Kd (0) dj=0 (2j + 1)Pj (u)2 ≤ √ Kd (0)(d + 1)2 . Thus (35) and b < 1/[4(d + 1)] imply |fi (j)| < 1/(2 3) < 1/2, leaving fixed slack for the rotation. For invalid addresses or malformed metadata define the amplitude to be zero; clamp it to [−1/2, 1/2] on the remaining register states. These reversible bounded extensions agree with the specified amplitudes on loaded valid records and define a uniformly conditioned rotation on the full register space.
9
(c) Output State: The local register is placed in a uniform superposition, with total success probability b2 . |i − 1⟩root |ai , si , Di ⟩meta |0⟩grid m −1
q i 1 X ⊗√ |j⟩local fi (j)|good⟩flag + 1 − fi (j)2 |bad⟩flag . mi j=0
(40)
3. Reduction to Exact Grover Search (a) Let A be the row-controlled state-preparation circuit, and use the Grover-type sequence (−AR0 A† Rg )t0 A.
(41)
Here Rg = I − 2Πg flips the good flag and R0 = I − 2Π0 reflects about zero preparation registers, leaving the root label and loaded metadata untouched. (b) Because the initial amplitude was artificially shrunk by the specific parameter b, this becomes an exact Grover search. Indeed, every row has good component b|χi ⟩|good⟩. The common angle θ0 = π/(4t0 + 2) rotates to (2t0 + 1)θ0 = π/2 with the same positive phase in every row [3, Sec. 2.1]. The scalar minus in the reflection product fixes this phase. (c) It perfectly isolates the target state |χi ⟩ in the ideal circuit. Arithmetic registers are zero because the preparation explicitly uncomputes them; the metadata remains loaded. This direct sum of equal-angle rotations also works on root-label superpositions. (d) Output State: The target state |χi ⟩ is cleanly isolated in the local register, the arithmetic memory is empty, and the success flag achieves exactly 100% probability. |i − 1⟩root |ai , si , Di ⟩meta |χi ⟩local |0⟩grid |good⟩flag
(42)
4. Address Mapping and Reversible Erasure (a) Shift the local index j into an absolute grid position by computing ℓ = ai + j into a new output register. XOR the computed value ℓ − ai into the local register and reverse the subtraction arithmetic. Then reverse the metadata-loading scan while the controlling root label is still present. The local and metadata registers are now both zero. (b) To safely discard the initial root label i − 1 without collapsing the quantum superposition, exploit the fact that the physical windows are strictly disjoint. Evaluate a piecewise reverse-lookup function λ(ℓ) = λ(ai + j) = i − 1 via a sequential scan, which uniquely identifies the root associated with any valid grid position. Set λ(ℓ) = 0 outside all windows. Each interval comparison conditionally XORs its hardwired label and uncomputes its test, so this defines a reversible operation on every grid input. (c) XOR this computed reverse-lookup value against the original input register. Because λ(ℓ) exactly equals i − 1, this cleanly subtracts the label from itself, restoring the root register to |0⟩ and isolating the final grid amplitudes. (d) Output State: The metadata, local index, and initial root label are perfectly uncomputed, transferring all window-state amplitudes solely into the active grid register.
mX i −1
Kd (ui,j ) √ |0⟩root |0⟩meta |0⟩local ⊗ |ai + j⟩grid |good⟩flag m D i i j=0 10
(43)
Toggle the now fixed good flag to zero. Each address operation uses O(log L) bits; comparisons and arithmetic reuse their workspace. Reversing this complete circuit gives the inverse of the clean extension of W . Having successfully constructed the isometry W , the circuit proceeds to use it to synthesize the complete adjoint averaging matrix R† . The scaled operator R† /A is a contraction, so it can be embedded as a block of a larger unitary. The unscaled R† need not shrink amplitudes. We achieve this by algebraically decomposing the target matrix into two distinct operations: s †
R = W diag(λi ),
λi =
vi Di . 2ρi
(44)
This decomposition is proven in the supplemental Sec. S2. • Step 1: The Diagonal Scaling (diag(λi )): On root label i − 1, a separate flag rotation has success amplitude λi /A ≤ 1/2. Load the row data, compute and apply the rotation, reverse its arithmetic, and unload the data while the root label is still present. Assign success to the zero flag state. The rotation angle is 2 arcsin(λi /A) in the usual Ry (θ) = e−iθY /2 convention, with a fixed flag relabeling if needed. • Step 2: The Spatial Mapping (W ): With the local quadrature scaling successfully encoded into the flag qubit’s probability amplitude, the circuit executes the exact forward sequence of W (Stages 1 through 4) to spatially distribute the root state across the physical grid. The separate row flag is left untouched by this sequence, including its inner amplification. Through block encoding, the combined effect of scaling the flag and applying W successfully embeds the target matrix into the overall quantum state. The block from root inputs with zero workspace to grid outputs with zero workspace and successful row flag is R† /A. Reversing the complete circuit and interchanging these designated subspaces gives R/A on arbitrary grid inputs. All flag amplitudes use bounded full-space extensions, as in the local preparation, and stay away from ±1. Thus amplitude error ξ gives angle error O(ξ). Polynomial evaluation, reciprocal square roots, and bounded-angle evaluation require polynomially many bits in d, log N, log L, log(1/ζ). In particular, the Legendre recurrence bounds the absolute coefficient sum of Pj by 3j , so that of Kd is at most C0 3d ; logarithmically many guard bits in this bound suffice. Reversible arithmetic and fixed-phase Clifford+T synthesis [14], with error allocated per gate occurrence, approximate the full unitary. Every inverse call reverses the actual compiled circuit. Telescoping includes the O(d + 1) inner amplification calls and all output registers. This yields a final logical gate count and total quantum width of O N poly(d, log N, log L, log(1/ζ))
and
poly(d, log N, log L, log(1/ζ)),
respectively. A rigorous bounding of the polynomial evaluation, arithmetic overhead, and unitary synthesis complexity is provided in Supplementary Section S3.
4.4
Complete amplification and choice of precision
Clean implementation also generates amplitudes outside the designated output subspace. The following two-reflection amplification estimate controls the full output with fixed phase; for projected robust amplification, see [4]. The input and output registers can have different layouts, so we use separate clean embeddings. Padding with zero registers and fixed register permutations put both layouts in one circuit space. 11
Lemma 6 (Fixed-phase amplification with all output registers). Let Q be an N × N unitary. Let U0 be a unitary with designated rank-N input and output embeddings Ein and Eout satisfying † Eout U0 Ein = Te/A,
A ≥ 2,
∥Te − Q∥ ≤ δ ≤ 1/24.
(45)
By utilizing O(A) alternating applications of U0 and U0† , alongside reflections about the designated subspaces and one additional qubit, we can construct an amplified unitary Uamp satisfying ∥Uamp Ein − Eout Q∥ ≤ (π + 2)δ.
(46)
The embeddings in the conclusion include the extra qubit in its designated fixed state. The full-output and fixed-phase estimate, including the component outside the clean output subspace, is proved in Supplementary Section S1.4. It gives a coherent implementation without postselection. To supply the input data for our amplified transform, we use the uniform-grid sampler of Supplementary Theorem S1. Based on the grid definitions in Eq. (22), for any target preparation error 0 < τ ≤ 1/10 and every power-of-two grid dimension L ≥ C(N + 1)3 τ −8 (where C is an absolute constant), this subroutine provides a coherent sampling unitary VL satisfying (47)
∥VL Ein − Egrid JL ∥ ≤ τ.
Here, Egrid embeds the L-dimensional sample register with all sampler work registers zero. This operation is efficient—requiring a circuit cost and width polynomial in log N , log L, and log(1/τ )— and maintains fixed column phases. The norm includes every output register; it therefore allows small residual workspace amplitudes. The tolerance also accounts for JL not being exactly isometric. The theorem permits a larger dyadic grid, allowing us to satisfy the window-resolution conditions simultaneously. The remaining choices separate naturally into two tasks: make the analytic and sampling errors small, then compile the resulting amplified circuit. Lemma 7 (Compatible parameter choice). For any N ≥ 2 and 0 < ε ≤ 1/2, there exist admissible choices for the polynomial degree d, grid size L, and precision tolerances such that all construction requirements are satisfied simultaneously. Specifically, √ d = O(log(N/ε)), A = O( d + 1), L = poly(N, 1/ε), (48) Ed (h) ≤ ε/24,
0 < τ ≤ 1/10,
Aτ ≤ ε/24.
(49)
The certified windows lie inside the grid, are pairwise disjoint, satisfy (23) and mi ≥ 4C0 C1 , and obey A ≥ 2∥R∥. The normalization A is fixed independently of the final grid refinement. Lemma 8 (Finite-precision compilation). Use the parameters of Lemma 7. Compose the actual sampler VL with the exact window block of Lemma 5, and apply the amplification construction of Lemma 6 with A = A. Compiling the window operations and additional coin rotation changes this complete circuit by at most ε/2 in full-unitary operator norm. The resulting Clifford+T circuit uses O N polylog(N, 1/ε) logical gates and polylog(N, 1/ε) logical qubits. The sampler circuit is identical in the comparison and compiled sequences, and every inverse is the reverse of the actual forward circuit. Classical certification and circuit generation take polynomial time in N, 1/ε and are charged separately. 12
The proofs of these two lemmas are in Supplementary Sections S4 and S5. Proof of Theorem 1. We begin by selecting compatible parameters according to Lemma 7. We compose the uniform-grid unitary VL in (47) with the ideal block for R/A from Lemma 5. The comparison circuit uses exact window rotations but the actual sampler at accuracy τ . Its good output subspace requires zero sampler and window workspace and a successful row flag. Write its block as Te/A. On ideal grid samples this block would be RJL /A = T /A by (32). Unitary composition and projection cannot increase the sampler error, so ∥Te − B † ∥ ≤ Ed (h) + Aτ =: δ ≤ ε/12 ≤ 1/24.
(50)
Padding and fixed register permutations supply the two rank-N clean embeddings. Apply Lemma 6 with Q = B † and A = A. The complete amplified output is within (π + 2)δ < ε/2 of Eout B † . Lemma 8 contributes at most ε/2 further error, proving (8) with the stated gate count and width. Tensoring with an identity includes any external reference system. † Finally, multiplying (8) on the left by UN,ε and on the right by B gives † Eout − Ein B∥ ≤ ε. ∥UN,ε
(51)
This confirms the inverse transform assertion using identical circuit resources.
5
Comparison to Known Implementations of the Transform
We compare Theorem 1 against existing complete-transform synthesis methods, targeting the transe suppresses factors polynomial form B with a clean-input operator-norm error bounded by ε. Here O in 1 + log N + log(1/ε). For our baselines, we first consider the complete reflection construction by Pliś and Zak, which requires quadratic total gates with logarithmic precision overhead [13]. We also compare against the unitary synthesis framework of Low, Kliuchnikov, and Schaeffer [11]. While their parameterization √ e 3/2 ) T -gate count by utilizing O( e N ) auxiliary width, the online coherentcan achieve an O(N e 2 ) total gate bound [11]. lookup Clifford gates still give a quadratic O(N Here λ follows Sections 2–3 of [11]: it counts b-bit lookup blocks and uses bλ borrowed qubits, with b = O(log(N/ε)). The lookup requires O(bN ) online Clifford gates per reflection-state preparation; the N -reflection synthesis therefore has a quadratic total-gate upper bound. Classical preprocessing is a separate cost. Table 1 summarizes these upper bounds. By exploiting adaptive windows, our approach reduces e ) while maintaining strictly polylogarithmic width. For the total gate complexity to near-linear O(N positive quantities, f ≍ g means cg ≤ f ≤ Cg for constants c, C > 0 independent of the varying parameters.
Conclusion and further questions By dynamically matching the averaging width to the local effective quadrature weight, we establish a near-linear finite-transform circuit with strictly polylogarithmic width. Our continuous-input analysis isolates the underlying circuit error from truncation and sampled-tail aliasing. The resource bounds are asymptotic upper bounds for this construction. Practical performance comparisons require a concrete compiler with certified numerical constants and gate counts. The oscillator and fractional Fourier examples illustrate operations diagonal in the retained Hermite space. Using them for hardware characterization will additionally require calibrated transform, 13
Construction
Online T gates
Online total gates
Total width
Adaptive windows, Theorem 1
e ) O(N
e ) O(N
polylog
Low–Kliuchnikov–Schaeffer, λ = 1
e 2) O(N
e 2) O(N
Low–Kliuchnikov–Schaeffer, √ λ≍ N
e 3/2 ) O(N
e 2) O(N
polylog √ e N) O(
Table 1: Upper bounds for the complete finite transform. All auxiliary qubits, including borrowed dirty qubits, contribute to width. preparation and measurement errors, as well as a verified noise model. Extending adaptive averaging to other quadratures is a further question: the local weight, separation and analytic estimates must be established for each new family.
References [1] Victor V. Albert, Kyungjoo Noh, Kasper Duivenvoorden, Dylan J. Young, R. T. Brierley, Philip Reinhold, Christophe Vuillot, Linshu Li, Chao Shen, S. M. Girvin, Barbara M. Terhal, and Liang Jiang. Performance and structure of single-mode bosonic codes. Physical Review A, 97:032346, 2018. https://arxiv.org/abs/1708.05010v3; https://doi.org/10.1103/PhysRevA.97.032346. [2] V. Bargmann. On a hilbert space of analytic functions and an associated integral transform part i. Communications on Pure and Applied Mathematics, 14(3):187–214, 1961. https://doi.org/10.1002/cpa.3160140303. [3] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp. Quantum amplitude amplification and estimation. https://arxiv.org/abs/quant-ph/0005055, 2000. arXiv:quant-ph/0005055v1; final version in Contemporary Mathematics 305 (2002), pp. 53–74, https://doi.org/10.1090/conm/305/05215. [4] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. https://arxiv.org/abs/1806.01838v1, 2018. arXiv:1806.01838v1; conference version: STOC 2019, pp. 193–204, https://doi.org/10.1145/3313276.3316366. [5] Gene H. Golub and John H. Welsch. Calculation of Gauss quadrature rules. Mathematics of Computation, 23(106):221–230, 1969. https://doi.org/10.1090/S0025-5718-69-99647-1. [6] Siddhartha Jain, Vishnu Iyer, Rolando D. Somma, Ning Bao, and Stephen Jordan. Efficient quantum Hermite transform. https://arxiv.org/abs/2510.04929v1, 2025. arXiv:2510.04929v1; conference version: STOC 2026, https://doi.org/10.1145/3798129.3800772. [7] Alexei Kitaev and William A. Webb. Wavefunction preparation and resampling using a quantum computer. https://arxiv.org/abs/0801.0342v2, 2009. arXiv:0801.0342v2; first version 2008.
14
[8] M. Alper Kutay and Haldun M. Ozaktas. The fractional fourier transform and harmonic oscillation. Nonlinear Dynamics, 29:157–172, 2002. https://www.ee.bilkent.edu.tr/~haldun/publications/ozaktas198.pdf. [9] A. L. Levin and D. S. Lubinsky. Christoffel functions, orthogonal polynomials, and Nevai’s conjecture for Freud weights. Constructive Approximation, 8:463–535, 1992. https://doi.org/10.1007/BF01203463. [10] A. L. Levin and D. S. Lubinsky. Erratum: Christoffel functions, orthogonal polynomials, and Nevai’s conjecture for Freud weights. Constructive Approximation, 11:417–418, 1995. https://doi.org/10.1007/BF01208563. [11] Guang Hao Low, Vadym Kliuchnikov, and Luke Schaeffer. Trading T gates for dirty qubits in state preparation and unitary synthesis. Quantum, 8:1375, 2024. https://doi.org/10.22331/q-2024-06-17-1375. [12] NIST. Digital library of mathematical functions. https://dlmf.nist.gov/. Sections 3.5(v)–(vi), 18.2, 18.7–18.9 and 18.18; accessed September 7, 2026. [13] Szymon Pliś and Emil Zak. Quantum discrete variable representations. https://arxiv.org/abs/2504.15841v1, 2025. arXiv:2504.15841v1. [14] Neil J. Ross and Peter Selinger. Optimal ancilla-free Clifford+T approximation of z-rotations. Quantum Information and Computation, 16(11–12):901–953, 2016. https://arxiv.org/abs/1403.2975; https://doi.org/10.26421/QIC16.11-12-1. [15] Gábor Szegő. Orthogonal Polynomials, volume 23 of American Mathematical Society Colloquium Publications. American Mathematical Society, Providence, RI, 4 edition, 1975. Theorem 8.22.9(a), Eq. (8.22.12). [16] Marcus Webb and Georg Maierhofer. Stable Hermite transforms via the Golub–Welsch algorithm. https://arxiv.org/abs/2604.02041v2, 2026. arXiv:2604.02041v2.
15
Supplementary proofs and additional results These supplementary proofs establish the analytic estimates, averaging-operator factorization, and assembly lemmas used in the main text. They also give the uniform-grid sampling contract, continuous-input error bounds, and introduce a mode-population noise diagnostic that we pose as an open problem for future hardware characterization. We use the article’s Hermite normalization, ordered roots x1 < · · · < xN of HN , weights wi , √ √ 2 and transform Bk,i = wi hk (xi ), with hk = Hk /(2k k! π)1/2 and φk = e−x /2 hk . Degree labels are 0 ≤ k < N and root i is stored as i − 1. Norms of finite-dimensional maps are operator norms. Circuit error includes all workspace, restored to zero in the ideal output, with fixed global phase and arbitrary external reference systems, as in article Equation (8). The uniform sampler in Section S6 uses only the analytic Mehler identity proved in Section S1.2, not its implementation theorem. Its proof therefore does not depend on the adaptive circuit whose construction uses this sampler.
S1
Supporting proofs for adaptive sampling
This appendix supplies the local weight and averaging estimates, exact normalization arithmetic, and full-output amplification bound used in main-text Section 4. The analytic estimates do not depend on the circuit construction.
S1.1
Weights and root separation
Proof of main-text Lemma 2. The roots of HN are the eigenvalues p of the real symmetric tridiagonal N × N matrix with zero diagonal and off-diagonal entries k/2, 1 ≤ k < N , by the Hermite √ recurrence [5, Sec. 2]. Every absolute row sum is below 2N , so the eigenvalue row-sum bound √ gives |xi | < 2N and the bounds on qi . Orthogonality of B gives the effective weight as vi =
N −1 X
!−1
φk (xi )
2
.
k=0 2
−1 2 −1 be the Christoffel function. Thus v = For W (x) = e−x /2 , let λN (W 2 , x) = ( N i k=0 hk (x) ) λN (W 2 , xi )/W (xi )2 . Write W = e−Q . Its Mhaskar–Rakhmanov–Saff number aN > 0 is determined by
P
N=
2 π
Z 1 0
aN tQ′ (aN t) √ dt. 1 − t2
We use the Christoffel estimate of Levin and Lubinsky [9, displayed estimate in the abstract]: for a fixed admissible Freud weight W = e−Q and fixed A > 0, −1/2 λN (W 2 , x) aN −2/3 ≍ max{N , 1 − |x|/a } , N W (x)2 N
|x| ≤ aN (1 + AN −2/3 ).
(S1)
The comparison constants are uniform in N, x. In the Gaussian case Q(x) = x2 /2 is even, Q′ (x) > 0 for x > 0, and (xQ′ (x))′ /Q′ (x) = 2, satisfying the source hypotheses. The defining integral gives N=
a2N ; 2
aN = 16
√
2N .
All roots therefore lie in the domain of (S1). An indexing convention that includes degree N gives the same value at xi , because HN (xi ) = 0. √ To compare the source’s edge scale with qi , put ui = |xi |/ 2N and δi = 1 − ui . These measure normalized position and distance from the edge. Define Mi = max{N 1/3 , N δi }. The source estimate gives vi ≍
√
(S2)
√ 2/ Mi . Since 0 ≤ ui < 1,
qi = N 1/3 + 2N δi (1 + ui ),
Mi ≤ qi ≤ 5Mi .
−1/2
Hence vi ≍ qi , proving main-text Equation (15). The root row-sum bound above and the differential-equation argument below do not use the zero-location corollary addressed in the source’s erratum [10]. To prove separation, fix one root x of a consecutive gap of length g and write q = 2N −x2 +N 1/3 . Assume g < q −1/2 , then for every coordinate y in the gap √ √ 2 2N 2 2 2N + 1 − y ≤ 2N + 1 − x + 2|x|g ≤ q + √ (S3) ≤ (1 + 2 2)q q |x−y|≤g where q 3 ≥ N was used in the last step. The Hermite equation [12, Table 18.8.1, Hermite row] gives φ′′N + (2N + 1 − y 2 )φN = 0. Multiplication by φN and integration between its two zeros, followed by the Dirichlet Poincaré inequality, gives √ π2 (2N + 1 − y 2 )|φN |2 |φ′ |2 R ≤ (1 + 2 2)q. (S4) ≤R N 2 = 2 2 |φN | |φN | g √ The denominator is nonzero. Since g 2 q < 1, this contradicts π 2 > 1 + 2 2, so the assumption is false and g ≥ q −1/2 . Applying the result at both endpoints proves main-text Equation (16). R
S1.2
R
The scaled complex Hermite bound
Proof of main-text Lemma 3. Mehler’s formula [12, Eq. 18.18.28], in the normalization of φk , gives ∞ X k=0
2
e
−tk
2
e− tanh(t/2)X +coth(t/2)Y p , |φk (X + iY )| = π(1 − e−2t ) 2
t > 0.
(S5)
Normal convergence extends the real polynomial identity to complex arguments. The generating p P 2 2zu−u k function e = k≥0 Hk (z)u /k!, Cauchy’s estimate on |u| = k/2, and Stirling’s upper bound imply, on every compact set of complex z, √ |Hk (z)| √ ≤ Ck 1/4 e 2k|z| , 2k k!
k ≥ 1.
Thus the two-variable Mehler series converges normally for |r| < 1. Apply the identity theorem in each variable, substitute X + iY and X − iY with r = e−t , and include the Gaussian factors to obtain (S5), including the positive Y 2 term. Set t = qi−1 ≤ 1 and write ci + ρi z = xi + ∆ + iY , with X = xi + ∆. The hypotheses give 3c |∆| = |(ci − xi ) + ρi Re(z)| ≤ ri + 2ri ≤ 3ri = √ qi 2c |Y | = |ρi Im(z)| ≤ 2ri = √ qi 17
(S6) (S7)
For k < N , e−tk ≥ e−N t . Multiplying the first N terms by eN t and then including the remaining nonnegative terms gives N −1 X
2
|φk (ci + ρi z)|2 ≤ eN t
k=0
e− tanh(t/2)X +coth(t/2)Y p π(1 − e−2t )
2
(S8)
The inequalities u − u3 /3 ≤ tanh u ≤ u and coth(t/2) ≤ 1 + 2/t ≤ 3/t imply 2
2
2
3 2
2
eN t−tanh(t/2)(xi +∆) +coth(t/2)Y ≤ et(N −xi /2)+t xi /24+t|xi ∆|+3Y /t √
2
≤ e1/2+1/12+3 2c+12c
−2
Here qi3 ≥ N controls the second and third terms. Concavity gives 1 − e−2t ≥ (1 − e−2 )t = 1−e for qi −1 0 ≤ t = qi ≤ 1. Combining this with the exponential bound (S8) yields main-text Equation (18), with √ 2 e7/12+3 2c+12c . (S9) CM = p π(1 − e−2 ) √ √ Multiply by vi ≤ C W / qi to obtain the row bound. The qi terms cancel, bounding the row q
vector norm by C W CM ≤ B0 . Summing the squared row bounds across all N rows bounds the √ Frobenius norm by B0 N .
S1.3
The three averaging errors
Proof of article Proposition 4. Let F (z) have rows fi (z) from main-text Lemma 3. The same dimensionless variable z parametrizes every row, even though their physical centers and radii differ. √ P j −j Write F (z) = j≥0 Fj z . These Taylor coefficients satisfy ∥Fj ∥ ≤ B0 N 2 by Cauchy’s integral formula on the radius-two disk. Thus its degree-d Taylor polynomial has uniform error at most √ −d B0 N 2 on [−1, 1]. Applying main-text Equation (20) to that polynomial to approximate F (0) and using the bound on Kd in main-text Equation (21) gives √ ∥Tcont − F (0)∥ ≤ B0 C0 N 2−d . (S10) The row lengths mi may also differ. For a real u ∈ [−1, 1] the radius-one disk about u lies in |z| ≤ 2; Cauchy’s estimate gives ∥fi′ (u)∥2 ≤ B0 . Hence gi (u) = Kd (u)fi (u) is Lipschitz with constant ∥gi′ (u)∥2 = ∥fi′ (u)Kd (u)+fi (u)Kd′ (u)∥2 ≤ B0 (C0 +C1 ). A vector-valued Lipschitz function of constant H has normalized midpoint error at most H/(2m): the integral of the distance to the midpoint of a cell of length 2/m is m−2 , and the integral is normalized by 1/2. Apply this separately to each row and then use the Frobenius norm to obtain B0 (C0 + C1 ) ∥T − Tcont ∥ ≤ 2
sX
m−2 i .
i
Finally apply main-text Lemma 3 with center xi and radius ri . Its derivative bound along the pP segment from xi to ci gives a row error at most B0 h/ri , and therefore ∥F (0)−B † ∥ ≤ (B0 h/c) i qi . Combining these three estimates gives the more detailed bound B0 h ∥T − B ∥ ≤ c †
sX
−d
qi + α d 2
i
18
B0 (C0 + C1 ) + 2
sX i
m−2 i .
(S11)
Since mi > ri /(2h) and
i qi ≤ 3N
P
2,
sX
2h m−2 i ≤
√ 2 3N h qi ≤ . c i
sX
c
i
Thus the first and third terms in (S11) sum to at most αh h, with the coefficients defined in main-text Equation (27). This proves article Equation (28).
S1.4
Amplification with fixed phase and full output
Proof of main-text Lemma 6. Put t = ⌈A⌉, θ = π/(4t + 2), and s = sin θ. Adjoin an independent coin with positive good amplitude β = As < 1; the inequality follows from As ≤ πA/(4t + 2) < 1. Require this coin to be good in the output subspace, choosing its good value to be zero. For the † † resulting unitary U , the good block is sTe. Write Πin = Ein Ein and Πout = Eout Eout , and let Rin = I − 2Πin and Rout = I − 2Πout . The amplification is Uamp = (−U Rin U † Rout )t U.
(S12)
The displayed scalar signs and the positive coin amplitude fix its global phase. The scalar −I has the exact Clifford realization σx σz σx σz = −I on one qubit. P Take a singular value decomposition Te = j σj |uj ⟩⟨vj |. Since Q is unitary, 1 − δ ≤ σj ≤ 1 + δ. Define αj = arcsin(sσj ) and the normalized bad vector |bj ⟩ =
U Ein |vj ⟩ − sin αj Eout |uj ⟩ . cos αj
The denominator is positive since s ≤ 1/2 and σj ≤ 1.1. The vectors are orthogonal to the output subspace. Subtracting the good-component Gram matrix from that of the unitary images gives ⟨bj , bk ⟩ = δj,k . Thus the planes span{Eout uj , bj } are pairwise orthogonal. In the ordered basis (Eout uj , bj ) the reflection product in (S12) is !
cos(2αj ) sin(2αj ) . − sin(2αj ) cos(2αj ) It follows, including the bad component, that Uamp Ein vj = sin((2t + 1)αj )Eout uj + cos((2t + 1)αj )bj . At σj = 1 the angle is (2t + 1)θ = π/2. On the interval of singular values in question the derivative of arcsin(sσ) is at most 2s, whence |(2t + 1)αj − π/2| ≤ 2(2t + 1)sδ ≤ πδ. The full output’s distance from Eout uj is twice the sine of half the angular deviation, at most πδ. P Orthogonality of the planes extends this bound to operator norm. Writing Qpol = j |uj ⟩⟨vj | for the unitary polar factor, we have ∥Uamp Ein − Eout Qpol ∥ ≤ πδ,
∥Qpol − Te∥ = max |1 − σj | ≤ δ. j
Together with ∥Te − Q∥ ≤ δ, this proves main-text Equation (46) for the unprojected output. Tensoring with an identity preserves operator norm and includes an external reference. 19
S2
Factorization of the averaging operator
Proof. The factorization in main-text Equation (44) follows entrywise. Recall the normalized rows and scales from article Lemma 5: W |i − 1⟩ = All coefficients are real, so
mX i −1
Kd (ui,j ) √ |ai + j⟩. m i Di j=0
Kd (ui,j ) (W † )i,ai +j = √ . mi Di
(S13)
(S14)
Since 2ρi = mi h, the row scale is s
λi =
vi Di = 2ρi
s
vi Di . mi h
(S15)
Consequently s
vi Di Kd (ui,j ) √ λi (W )i,ai +j = mi h mi Di √ vi Kd (ui,j ) √ = = R[i, ai + j]. mi h
†
(S16)
Both sides vanish outside the row window, by article Equation (29). Thus R = diag(λi )W † and
R† = W diag(λi ).
S3
Finite-Precision Synthesis and Complexity Analysis
To approximate the ideal circuit, we must first control the arithmetic error in each flag amplitude. The separation of the amplitudes from ±1 bounds the derivative of the arcsin function, ensuring that an arithmetic amplitude error ξ translates to a bounded rotation error O(ξ). Rational polynomial evaluation, reciprocal square roots on Di ≥ 1/2, and bounded-angle functions admit arithmetic circuits polynomial in d, log L, and the required output precision. For example, let Lj be the sum of the absolute monomial coefficients of the polynomial Pj . The Legendre recurrence gives Lj+1 ≤ 2Lj + Lj−1 , and hence Lj ≤ 3j . Thus, the coefficient sum of the kernel Kd is at most C0 3d . Horner evaluation on |u| ≤ 1 with absolute error ξ uses O(d + log(d + 1) + log(1/ξ)) magnitude, guard, and fractional bits, up to fixed arithmetic overhead. Integer addresses and dyadic arguments are computed exactly. The same argument applies to the row flag, using certified approximations of λi /A with its fixed slack. We choose a working-bit bound polynomial in d, log N , log L, and log(1/ζ). By retaining intermediate values, copying the output, and reversing the computation, Boolean arithmetic is made reversible in polynomial time and space. Bit-controlled rotations synthesize the angles with polynomial overhead in this bound and log(1/ζ), using fixed-phase Clifford+T approximation. We rely only on the existence of this overhead, not typical synthesis counts. At every inverse call, we use the actual synthesized circuit in reverse. We compare the resulting full unitary sequence with its ideal sequence by telescoping, allocating the error per occurrence, including the O(d+1) calls to the amplification block. This bounds the accumulated approximation 20
error on all output registers with logarithmic extra precision in the occurrence count. The bounded extensions ensure uniform approximation on every address encountered during amplification. Metadata remains loaded throughout the inner amplification sequence; loading, unloading, and the final interval scan each cost O(N ) times a polynomial in bit length. All remaining operations have polynomial cost in that length and d, proving the total gate and width bounds for the circuit and its inverse.
S4
Proof of Compatible Parameter Choice
Proof. We prove main-text q Lemma 7. For N ≥ 2 and 0 < ε ≤ 1/2, certify C W as in article
Equation (17) and fix B0 ≥ C W CM from the analytic bound. Choose an integer d ≥ 1 and then A as in article Equation (36), with √ ε B0 (d + 1)2 N 2−d ≤ ε/48, τ= . (S17) 24A √ The exponential decay permits d = O(log(N/ε)), so A = O( d + 1) and 0 < τ ≤ 1/10. These choices precede the grid size: neither error coefficient in article Equation (27) depends on L. Choose a power-of-two L ≥ C(N + 1)3 τ −8 , as required by Theorem S1, large enough that (
h ≤ min
c c cε √ √ , , √ 16 3N 8C0 C1 3N 48 3B0 (1 + C0 + C1 )N
)
,
Lh √ ≥ 2N + 2. 2
(S18)
The first bound permits the certified disjoint windows of article Equation (23). The factor-four admissible radius interval permits a dyadic choice with a margin; a root enclosure narrower than h/4 determines a half-grid√center within h of the root. The second bound gives mi ≥ 4C0 C1 , since mi > ri /(2h) and ri ≥ c/ 3N . The last condition contains every window in the grid. The third bound and (S17) give Ed (h) ≤ ε/24. Every lower bound on L = 2π/h2 is polynomial in N, 1/ε; rounding the largest to the next power of two gives L = poly(N, 1/ε) and log L = O(log(N/ε)). Reselecting the windows on this grid preserves A ≥ 2∥R∥, because the bound in article Equation (36) holds for every admissible grid. Finally Aτ = ε/24, which is the separate sampler error allocation used in the main proof.
S5
Proof of Finite Precision Compilation
Proof. We prove main-text Lemma 8 for the circuit assembled there, not for an arbitrary ideal circuit. Put ε t = ⌈A⌉, M = 2t + 1, ζ= . 4M The fixed-phase amplification sequence (S12) has exactly M forward or inverse block calls. Each block comprises the actual sampler VL , the averaging unitary, and a coin whose positive good amplitude is A sin(π/(4t + 2)) < π/4. Approximate the averaging unitary and the coin each to operator error ζ. Article Lemma 5, with the bounded extensions and guard-bit analysis of Section S3, includes the inner amplification, arithmetic, metadata scans, and cleanup in the averaging-unitary budget. The fixed coin slack makes its angle uniformly conditioned. The input and good-output reflections test specified zero registers and fixed register locations; they have exact Clifford+T implementations with reusable clean workspace. Use the same compiled VL in the comparison and implemented sequences. Its full-output accuracy τ has already entered article Equation (50); no exact implementation of the nonisometric 21
sampling matrix JL is assumed. Every inverse call is the actual reverse of the corresponding circuit, including VL . Thus each block differs from its comparison block by at most 2ζ, and full-unitary telescoping adds at most 2M ζ = ε/2 on all output registers. Variations in row success amplitudes and nonzero residual workspace are included in this norm estimate, rather than discarded by a final projection. √ Here M = O( d + 1) and log(1/ζ) = O(log(N/ε)). Each call uses charged O(N )-record scans and arithmetic of polynomial logarithmic cost; the sampler, coin, and reflections have polynomial logarithmic cost as well. The online logical Clifford+T gate count is therefore O(N polylog(N, 1/ε)), which also bounds depth. Serial scans and inverse computations reuse their workspace, giving total width polylog(N, 1/ε). The O(N ) loaded records have polynomial logarithmic length. Their classical certification uses real-root isolation for the degree-N integer polynomial HN , interval evaluation of the weight formula, and direct rational evaluation of the normalization constants Di . The coefficients of HN have O(N log N ) bits; the root-gap bound in main-text Lemma 2 and the dyadic margins above make the required isolations finite with polynomial bit cost. For the normalizations, write Kd = P/Q with integer P of degree at most d; a common denominator for Di divides Q2 m2d+1 . The Legendre i recurrence and Di ≤ d + 5/4 bound the denominator and numerator bit lengths polynomially in d + log mi . Since mi ≤ L = poly(N, 1/ε), direct summation therefore has polynomial classical bit cost. Gaussian factors can require polynomially many classical guard bits during certification; only the shorter certified records enter the online circuit. Classical certification and circuit-description generation are polynomial in N, 1/ε and are charged separately from the online gate and width bounds.
S6
Coherent sampling of Hermite expansions on a uniform grid
Jain et al. construct a uniform-grid Hermite sampler for arbitrary coefficient superpositions, including coherent degree-label erasure [6, Theorem 19, Lemma 27 and Algorithm 2]. Their factor (−1)k is removed by an input Z on the least significant degree bit when N ≥ 2, and by the identity when N = 1. We retain a direct sampled-function proof of the clean interface used in the adaptive conversion. It follows their preparation, oscillator-phase selection, and amplification strategy, without requiring the complete spectrum of a discretized p Hamiltonian. Throughout, L ≥ 8 is dyadic, m = log2 L, h = 2π/L, and IL = [−L/2, L/2). Grid labels ℓ ∈ IL ∩ Z use a fixed reversible signed encoding, such as ℓ mod L. For 1 ≤ d < L, set √ (d) (N ) JL [ℓ, k] = h φk (hℓ), 0 ≤ k < d, JL = JL , (S19) √ 2 where φk (x) = e−x /2 Hk (x)/(2k k! π)1/2 has positive leading coefficient. Norms of finitedimensional maps are operator norms; absolute constants may increase between estimates. Theorem S1 (Clean uniform-grid sampling). There is an absolute constant C with the following property. Let N = 2n ≥ 1 and 0 < τ ≤ 1/10. For every power of two L ≥ C(N + 1)3 τ −8 ,
(S20)
there is a logical quantum circuit VL with clean input and output embeddings Ein : CN −→ Hall and Egrid : CL −→ Hall such that VL Ein − Egrid JL ≤ τ. (S21) The input embedding places the degree label in a zero-padded register and initializes all other registers to zero. The output embedding retains only the grid register and sets the degree, phase, and 22
arithmetic registers to zero. The logical gate count and total number of qubits are polynomial in log L + log(1/τ ). The estimate includes the phases of all columns and holds after tensoring with an arbitrary reference system. The constant and the circuit family are qualitative asymptotic ones. The construction uses the classical bulk Hermite asymptotic and the standard phase-preserving singular-vector amplification theorem stated below; it does not specify a numerical value of C, an optimized gate exponent, or a physical-hardware compilation.
S6.1
Sampled oscillator phases
Let Ff (ξ) = (2π)−1/2 e−ixξ f (x) dx and FL [j, ℓ] = L−1/2 e−2πijℓ/L . Both Fourier transforms have negative sign. Define R
2
CL (a) = diagℓ∈IL ∩Z (e−ia(hℓ) ),
A(u) = 12 tan(u/2),
B(u) = 12 sin u,
and the bounded half-time word VL (t) = [FL∗ CL (A(t/2))FL CL (B(t/2))FL∗ CL (A(t/2))FL ]2 ,
|t| ≤ π.
(S22)
The following estimate controls its full output, including leakage from the sampled Hermite space. (d)
Lemma S2. Let J = JL and Dd (t) = diagk<d (e−it(k+1/2) ). C1 , C2 , c > 0 such that, with E(d, L) = C1 L1/4 eC2 d−cL ,
There are absolute constants
∥J ∗ J − Id ∥ + sup ∥VL (t)J − JDd (t)∥ ≤ E(d, L).
(S23)
|t|≤π
For efficiently computable t, clean-workspace implementation of the word to full-output error ζ uses poly(log L + log(1/ζ)) gates and qubits. Proof. Complex Mehler, Equation (S5) in Section S1.2 (see also [12, Eq. 18.18.28]), gives, with P a = tanh(1/2), for f = k<d αk φk and ∥α∥2 = 1, 2
2
(S24)
|f (x + iy)| ≤ Ced/2−ax /2+y /(2a) . For real v, the Hermite generating function and a Gaussian Fourier integral give 2
F(eivx /2 φk )(ξ) = s−1/2 ei arctan(v)/2
v−i s
k
2
2
e−ivξ /(2s ) φk (ξ/s),
s=
p
1 + v2.
The branch is continuous at v = 0, fixing the scalar phase. Thus Fourier transformation preserves 2 the family σ −1/2 eiκx /2 ei(kθ+ϕ) φk (x/σ), where σ > 0 and κ, θ, ϕ ∈ R, with √ 1 + κ2 σ 4 κσ 4 ′ σ = , κ′ = − . σ 1 + κ2 σ 4 Inverse Fourier transformation adds a reflection, and a chirp adds its coefficient to κ. The finitely many operations in (S22) have chirp coefficients of magnitude at most one. Their continuous parameter maps keep the widths bounded above and away from zero and the chirps bounded. Consequently every ideal continuum prefix g and its next Fourier or inverse Fourier transform 2 satisfy |g(x)| + |Fg(x)| + |F −1 g(x)| ≤ Ced/2−cx , uniformly in t and α.
23
Write SL g =
√
h(g(hℓ))ℓ∈IL ∩Z . Poisson summation, with 2π/h = Lh, gives
(FL SL g − SL Fg)j =
q √ X X h Fg(hj + qLh) − h/L g(hℓ)e−2πijℓ/L . ℓ∈I / L
q̸=0
The Gaussian bounds and |hj + qLh| ≥ (2|q| − 1)Lh/2 show that its norm is at most CL1/4 ed/2−cL ; the inverse transform obeys the same estimate. Sampling commutes exactly with a chirp. Telescope the eight Fourier steps using these bounds on ideal prefixes and unitary grid suffixes. The correb2 )/2, pb = −i d/dx, and x bf = xf : sponding continuum word equals e−itHosc , where Hosc = (pb2 + x apply the three-chirp identity [6, Theorem 3, Eq. (6)] twice at t/2, choosing the continuous branch equal to the identity at time zero. Since Hosc φk = (k + 1/2)φk , this proves the phase bound. 2 2 For the Gram bound, the entire function Gf (z) = f (z)f (z) satisfies |Gf (x+iy)| ≤ Ced−ax +y /a . 2 Shifting its Fourier contour to y = −aω/2 gives |FGf (ω)| ≤ Ced−aω /4 ; Gaussian decay removes the vertical boundary terms. Poisson summation and removal of the spatial tails therefore give P h ℓ∈IL |f (hℓ)|2 = 1 + O(ed−cL ). The uniform Rayleigh-quotient bound proves the Gram estimate. Finally, a Fourier circuit has O(m2 ) gates before finite-set synthesis. A chirp computes the signed square, applies its phase, and uncomputes. Since (hℓ)2 = O(L), coefficient accuracy O(ζ/L) suffices. Arithmetic, bounded trigonometric evaluation and rotation synthesis [14] have polynomial cost in log L + log(1/ζ); allocate error across their gates. Compute–phase–uncompute gives the stated clean-input, full-output bound. For analysis only, when γ = ∥J ∗ J − Id ∥ ≤ 1/2, set G = J ∗ J,
Q = JG−1/2 ,
qk = Q|k⟩.
Functional calculus on [1/2, 3/2] gives Q∗ Q = Id and ∥Q − J∥ ≤ Cγ. Adding and subtracting J in (S23) shows that a compiled word VeL (t) has error δ = C(E(d, L) + γ + ζ) relative to QDd (t) on this frame, including all output workspace. No commutation of G and Dd (t) is needed, and Q is not computed by the circuit.
S6.2
Preparation and clean degree-label erasure
Proof of Theorem S1. We construct smooth seeds, select their Hermite component, and amplify before erasing the retained degree label. Seeds and sampled tails. Szegő’s bulk formula [15, Theorem 8.22.9(a), Eq. (8.22.12), p. 201], in our positive-leading normalization, is √ 21/4 sin((k/2 + 1/4)(sin 2θ − 2θ) + 3π/4) φk ( 2k + 1 cos θ) = + O(k −5/4 ), (S25) √ 1/4 √ πk sin θ √ p uniformly on fixed compact subintervals of (0, π). Put ak = 3(2k + 1)/4 and wk = (10 2k + 1)−1 . Let gk equal one on |x| ≤ ak , zero on |x| ≥ ak + wk , and 1 − 10s3 + 15s4 − 6s5 in √ between, with s = (|x| − a )/w . Its support stays in a fixed angular bulk interval, since |x|/ 2k + 1 ≤ k k √ 3/2 + 1/10 < 1. For k ≥ k0 , define fk as gk times the explicit sine expression in (S25); for the finitely many k < k0 , use fk = gk φk . Set fk = 0 outside the support and evaluate the bulk expression only inside it. The cutoff is C 2 . Direct differentiation of the explicit expressions, not of the asymptotic remainder, gives (r) ∥fk ∥∞ ≤ C(k + 1)r/2−1/4 (0 ≤ r ≤ 2),
1 2
24
Z
(|fk′ |2 + x2 |fk |2 ) dx ≤ C(k + 1).
(S26)
√ Here the derivative and inverse cutoff width √ support length is O( R k 2+ 21) and both the phase R ak 2 are O( k + 1). The moment x φk = k + 1/2 implies −ak φk ≥ 1/3. The bulk remainder gives ∥gk (fek − φk )∥2 ≤ C/k, where fek is the sine expression. Choosing k0 large enough yields ⟨φk , fk ⟩ ≥ b > 0 and b ≤ ∥fk ∥2 ≤ C for all k, including zero. For a sufficiently large absolute constant Cs , set Bk = (k + 1)1/4 fk /Cs ∈ [−1/2, 1/2]. Elementary-function evaluation with angle precision O(σ/(k + 1)) and guard bits evaluates Bk to error σ with reversible resources polynomial in log L + log(k + 1) + log(1/σ). The fixed-degree exceptional branch and the C 2 cutoff joins have the same bound; no growing table is required.√ P For N ≤ d < L and 0 ≤ k < N , put Fk = h ℓ∈IL fk (hℓ)|ℓ⟩, sampling the entire seed support. P Write cj = ⟨φj , fk ⟩ and α = J ∗ Fk . The oscillator quadratic form gives j≥d |cj |2 ≤ C(k + 1)/d. (r)
The ladder relations and ∥u∥2∞ ≤ 2∥u∥2 ∥u′ ∥2 give ∥φj ∥∞ ≤ Cr (j + 1)r/2+1/4 . Leibniz’s rule and full-support trapezoidal quadrature thus yield √ C(N + 1) C N + 1(d + 1)5/4 ∥Fk ∥2 − ∥fk ∥22 ≤ , . |αj − cj | ≤ e := L L There is no endpoint term. The exact identity ∥(I−QQ∗ )Fk ∥2 = ∥Fk ∥2 −α∗ G−1 α, ∥G−1 −I∥ ≤ Cγ, and the squared-coefficient error O(de) give Fk 2 ≤C (I − QQ ) ∥Fk ∥ ∗
N + 1 (d + 1)11/4 + +γ . d L !
(S27)
For the small errors chosen below, c ≤ ∥Fk ∥ ≤ C and bk := ⟨qk , Fk /∥Fk ∥⟩ ≥ c > 0, because ⟨qk , Fk ⟩ = ck + O(e + γ). All these values and matrices are real, fixing the overlap phase. Selection and signed preparation. Let U be phase estimation using all m controlled providers it /2 j e VeL (tj ), tj = 2π2j /L, j = 0, . . . , m − 1, and the positive-sign Fourier transform on its phase register. The scalar is a relative phase on the control-one branch. Extend the ideal shifted action qk 7→ e−itj k qk by the identity on (ran Q)⊥ , obtaining an ideal circuit U0 . Character orthogonality gives U0 (|0⟩qk ) = |k⟩qk since d < L. Its intermediate operations preserve the space M with arbitrary phase register, grid in ran Q, and zero work. Unitary hybrids, including Fourier synthesis, give ∥(U ∗ − U0∗ )|M ∥ ≤ Cmδ.
∥(U − U0 )|M ∥ ≤ Cmδ,
(S28)
For the inverse use U ∗ − U0∗ = U ∗ (U0 − U )U0∗ and U0∗ M = M. Actual frame invariance is unnecessary. Retaining k, prepare a uniform dyadic √ interval of Tk grid points containing the whole seed support plus one mesh step, with Tk h = Θ( k + 1). Conservative binary endpoints suffice. Compute 2 arcsin Bk (hℓ), apply the corresponding Ry rotation to a flag initially in |0⟩, then apply X and uncompute evaluation. With Ry (2θ)|0⟩ = cos θ|0⟩ + sin θ|1⟩, the zero-flag amplitude is exactly Bk (hℓ), including its sign. The fixed slack |Bk | ≤ 1/2 bounds the derivative of arcsin, so angle approximation and rotation synthesis give full input-isometry error O(σ) with the same polylogarithmic resources. The selected vector is κk Fk + O(σ), where κk =
(k + 1)1/4 √ , Cs Tk h
0 < c ≤ κk ≤ C.
The error is independent of Tk because the initial interval is normalized. Append U , preserving the raw failure flags, and subtract k from the phase register. Call this compiled unitary A and let
25
Πgood require the zero raw flag and zero phase and work registers. Suppressing the retained label, put zk = Πgood Ak |0⟩. If the tail norm in (S27) is at most η, then ∥zk − βk (qk ⊗ |0⟩)∥ ≤ C(σ + η + mδ) =: ηg ,
βk = κk ∥Fk ∥bk ≥ a∗ > 0.
(S29)
Indeed, ideal selection extracts exactly the qk coefficient of the frame component, and unitarity bounds the tail contribution by its norm. Raw failure branches cannot enter the good subspace. Amplification and erasure. We use phase-preserving singular-vector amplification [4, Theorems 26–27 and Corollary 10, arXiv v1]: if the nonzero singular values of Πgood AΠin lie in [s∗ , 1], a common polynomial gives real overlap at least 1 − ρ with each normalized good singular vector, ∗ using O(s−1 ∗ log(1/ρ)) calls to A, A and projector phases. This external input fixes phases, not just success probabilities. Apply it to the actual A and its exact circuit inverse, with Πin selecting k < N and zero input auxiliaries. For ηg ≤ a∗ /2, the actual success amplitude sk = ∥zk ∥ ≥ a∗ /2 and direction vk = zk /sk satisfy ∥vk − qk ⊗ |0⟩∥ ≤ 4ηg /a∗ . Retained labels make these selected columns orthogonal, with singular values exactly sk . Choose ρ = ε2amp /8: real overlap 1 − ρ gives √ vector error at most 2ρ. Allocate another εamp /2 to amplification-phase synthesis. Preparation of |k⟩qk |0⟩ then has operator error at most εamp + 4ηg /a∗ , taking the maximum over orthogonal label blocks. Seed and selection errors enter once through the actual good direction, not once per amplification call. Apply the exact compiled U ∗ using the retained degree register as its phase register. For P P arbitrary coefficients λk , its ideal inverse sends k<N λk |k⟩qk to |0⟩ k<N λk qk . Equation (S28) costs Cmδ; unitarity preserves the preceding error. Replacing Q by J costs Cγ. Thus the fulloutput error is at most εamp + C(σ + η + mδ)/a∗ + Cmδ + Cγ. (S30) Parameters and resources. Choose σ, η, εamp as sufficiently small fixed multiples of τ , then take d = ⌈C3 (N + 1)τ −2 ⌉,
L ≥ C4 (d + 1)3 τ −2 ,
ζ ≤ c0 τ /(1 + m).
Increasing C3 , C4 and decreasing c0 makes (S27) at most η 2 and (S30) at most τ . Indeed, the grid bound ensures full support containment, d < L, and quadrature error O(τ 2 ), and also E(d, L), γ ≤ c0 min{τ 2 , τ /(1 + m)}. For this last assertion keep d fixed as L increases: half of −cL absorbs C2 d, and the other half dominates L1/4 (1 + log L) and the required powers of 1/τ already at the polynomial floor. These bounds decrease beyond a fixed threshold. Hence (S20), with a sufficiently large absolute constant, works for every larger dyadic grid. Each selection uses m providers, and amplification uses O(log(1/τ )) calls. Seed evaluation, label arithmetic, projector tests, finite-gate synthesis with occurrence-wise allocation, and exact circuit inverses have polynomial cost in log L + log(1/τ ). Evaluation is uncomputed; all temporary, phase and amplification registers have the same polynomial width. There is no unit-cost state-table oracle. The estimates include every auxiliary output and are operator bounds, so they survive tensoring with an arbitrary reference system. This proves the clean interface and its logical resource claims. The fixed thresholds and classical-input constants remain qualitative. Circuit-description generation is separate from the online quantum count; no numerical constant, optimized classical compiler or hardware realization is asserted. A polynomially chosen grid gives polynomial dependence on log N + log(1/τ ), while larger grids retain the explicit log L dependence.
S7
Sampling outside the finite Hermite space
An exact finite unitary need not be an exact continuous transform. Let PN be the orthogonal projection of L2 (R) onto HN . For a function ψ with well-defined samples, write r = ψ − PN ψ and 26
let c<N be the coefficient vector of PN ψ. The finite output is BSN ψ = c<N + BSN r,
∥BSN ψ − c<N ∥2 = ∥SN r∥2 .
(S31)
Moreover, interpreting the output as a function in HN gives the orthogonal error decomposition N −1 X
2
k=0
L2
= ∥SN r∥22 + ∥r∥2L2 .
(BSN ψ)k φk − ψ
(S32)
The first term is aliasing from the sampled tail; the second is truncation. These identities follow directly from (5) and orthogonal projection. They do not assume that point evaluation is bounded on L2 . Two immediate examples show why a sampling assumption is necessary. Since every xi is a root of HN , SN φN = 0. Applying the Hermite recurrence [12, Sec. 18.9] at those roots also gives s
BSN φN +1 = −
N eN −1 , N +1
(S33)
where eN −1 is the final number-basis coordinate. Thus a component just outside the cutoff can disappear entirely or alias into a retained component. There is no constant CN that bounds ∥SN r∥2 by CN ∥r∥L2 for every smooth r ⊥ HN . To see this, place a smooth bump of width h and height h−1/2 at one quadrature node, with fixed nonzero central value after rescaling. Its L2 norm is fixed and its sampled value diverges. Subtract its projection onto the fixed finite space. Each projection coefficient is O(h1/2 ), so the subtraction neither cancels the diverging sample nor makes the residual norm diverge. This is the usual unboundedness of point evaluation on L2 , specialized to the finite sampling operator. A stronger input assumption does provide a direct bound. If r ∈ H 1 (R), then its continuous representative satisfies ∥r∥2∞ ≤ ∥r∥L2 ∥r′ ∥L2 . One proof integrates (|r|2 )′ from both sides of a point P 2 and applies Cauchy–Schwarz. Consequently, with WN = i wi exi , ∥SN r∥22 ≤ WN ∥r∥L2 ∥r′ ∥L2 .
(S34)
For the oscillator Hosc = (−d2 /dx2 + x2 )/2, its quadratic form is qosc [r] = 12 ∥r′ ∥2L2 + ∥xr∥2L2 ,
r ∈ H 1 (R),
xr ∈ L2 (R).
Thus ∥r′ ∥2L2 ≤ 2qosc [r] on the form domain. On the smaller operator domain, qosc [r] = ⟨r, Hosc r⟩ as well; the sampling bound requires only the form-domain assumption. Equations (S32) and (S34) give a fully specified, although not asserted sharp, sufficient error bound. Determining sharper N dependence under an appropriate Hermite regularity assumption is a further approximation-theory question. For an actual quantum input one must normalize SN ψ. If a = ∥c<N ∥2 > 0 and e = ∥SN r∥2 < a, then the normalized finite output and the normalized projected coefficient state satisfy c<N BSN ψ − ∥SN ψ∥2 a
≤ 2
2e . a
(S35)
This follows by adding and subtracting BSN ψ/a and using the reverse triangle inequality. A circuit approximation of operator-norm error ε adds at most ε to this state-vector error. If SN ψ = 0, no normalized sampled input exists. 27
For the oscillator application, evolution and measurement already specified in the number basis can be performed directly in that basis. The scalar e−it/2 in the position-basis evolution becomes a relative control phase in a controlled implementation. More general multiplication potentials require a separate approximation: replacing a potential by its values at the quadrature nodes is a discrete-variable-representation approximation, whose error must be controlled independently of the exact finite transform [13, Appendix A, Eqs. (42)–(45)]. A direct composition check can form d = B † c for a normalized coefficient vector and compare B † D(t)Bd with B † D(t)c. This tests the conventions and composition; compiled gate counts require the separate resource analysis.
S8
A Mode-Population Noise Diagnostic
A transform–idle–inverse experiment can show how a qubit noise channel redistributes Hermitemode populations. The following protocol describes a possible diagnostic. The role of B is to choose the preparation and measurement bases, rather than to establish a noise-characterization advantage. In an ideal-transform model, noise is confined to the idle interval by assumption; a hardware implementation would also require characterization of state preparation, transform errors, and readout. For example, consider a 4-qubit register (N = 16), initially in the number state |12⟩. With the convention (4), B † maps this state to weighted position samples; after an idle interval ∆t, B returns it to the number basis for measurement. If Φ∆t is the noise channel on the sample register, the measured state is ρout = B Φ∆t B † ρin B B † , ρin = |12⟩⟨12|. (S36) Figure S1 gives this protocol. id
|0⟩
id
|0⟩ Init |12⟩
B†
B
|0⟩
id
|0⟩
id
Figure S1: Transform–idle–inverse protocol. The number state is mapped to weighted position samples by B † , subjected to noisy idle (id) operations for time ∆t, and mapped back by B for number-basis measurement. The observable populations are pk (∆t) = ⟨k|ρout |k⟩, whose sum is one for a trace-preserving channel. To illustrate this protocol, take Φt = ϕ⊗4 t , where each qubit undergoes thermal relaxation with T1 = 50 µs, coherence-decay time T2 = 30 µs, and equilibrium excited-state population peq = 0.05. For a single-qubit density matrix σ, the channel is specified by (ϕt (σ))11 = peq + (σ11 − peq )e−t/T1 , (ϕt (σ))01 = σ01 e−t/T2 ,
(S37)
with the remaining entries fixed by trace one and Hermiticity. We use the binary ordering |q3 q2 q1 q0 ⟩, with q0 least significant. Figure S2 is a direct density-matrix calculation of (S36) at 28
integer-microsecond intervals, using ideal dense transforms. Qubit relaxation in this encoding need not relax the oscillator to its ground state: even a complete reset of the sample register gives B|0⟩ in the number basis, rather than |0⟩. Here the limiting state is B diag(1 − peq , peq )⊗4 B † . The accompanying reproducibility package contains the simulator, the full population table, and numerical consistency checks. Ideal transforms; independent qubit thermal relaxation 1.0
30
0.8
20 0.6 15 0.4 10
Mode population
Idle time (microseconds)
25
0.2
5
0
0.0 0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Hermite mode index k
Figure S2: Simulated mode populations for the transform–idle–inverse protocol with N = 16, initial mode 12, and independent identical qubit channels (S37). The parameters are T1 = 50 µs, T2 = 30 µs and peq = 0.05. Each row gives all 16 populations at one idle time. The calculation uses ideal dense transforms and density matrices, without shot sampling, gate compilation, or hardware execution.
S8.1
Population and coherence tests
Population measurements determine only part of a noise channel. The distinction is already visible for bosonic pure loss, whose single lowering operator gives the generator κD[aN ] on the retained number states [1, Sec. I.A, preceding Eq. (1.3)]. Here κ ≥ 0 is the loss rate and aN =
N −1 √ X
k|k − 1⟩⟨k|,
D[L](ρ) = LρL† − 12 {L† L, ρ}.
k=1
Its infinitesimal population transition k → k − 1 has rate κk. Independent qubit relaxation instead has local lowering operators: in binary encoding, decay of bit j changes the encoded integer by 2j , at a device-dependent rate. A basis change conjugates the channel as in (S36); a model of encoded qubit noise and a model of bosonic loss therefore require distinct physical assumptions. To see why populations alone are insufficient, put Ekℓ = |k⟩⟨ℓ| and compare L1 = κD[aN ],
L2 = κ
N −1 X
kD[Ek−1,k ].
k=1
Both preserve number-diagonal density matrices and act identically on them, so every numberstate input has the same population evolution under both generators at every time. For distinct 29
1 ≤ k, ℓ < N , however,
√ (L1 − L2 )(Ekℓ ) = κ kℓ Ek−1,ℓ−1 .
(S38)
The anticommutator terms agree, while the √ collective jump transports off-diagonal coherence. For N ≥ 3 and κ > 0, the state (|1⟩ + |2⟩)/ 2 distinguishes the two generators. They coincide for N ≤ 2 or κ = 0. Number-state inputs probe population leakage, while superposition inputs with phase-sensitive readout can probe coherence transport. A zero-idle round trip supplies a baseline for the combined preparation, transform, and readout errors. These are proposed extensions of the diagnostic; the population simulation above provides neither a coherence measurement nor a device calibration.
30