LCHS

LCHS (linear combination of Hamiltonian simulation) computes the solution at time \(T\) of the linear ODE \(du/dt=-Au+b\) for a constant matrix \(A\), an initial vector \(u(0)\) and an optional constant source \(b\). It splits \(A=L+iH\) into the Hermitian parts \(L=(A+A^\dagger)/2\) and \(H=(A-A^\dagger)/(2i)\) and, for positive semidefinite (PSD) \(L\) and the default near-optimal kernel \(g_\beta\), writes the propagator as a weighted integral of unitary evolutions,

\[ e^{-AT}=\int_{\mathbb R}g_\beta(k)\,e^{-iT(kL+H)}\,dk, \]

which NWQLib truncates to \(|k|\le K\) and replaces by a finite quadrature sum (An, Childs and Lin, arXiv:2312.03916v2, here ACL, Eqs. (3)–(4), (6)–(7) and Theorem 6, stated as Result 25). A constant source adds the Duhamel term \(\left(\int_0^Te^{-As}\,ds\right)b\) of their Eq. (2).

Use LCHS for the state of a linear dynamical system at a given time. For a steady state \(Ax=b\), use QLS. For the default kernel–quadrature pair, approximation_tolerance=0.01 limits two discretization error components, the cutoff tail and the k quadrature. It does not bound the total error of the output (Selection and accuracy).

Solve a linear ODE

from nwqlib import LinearDynamics, solve
from nwqlib.algorithms import LCHS

problem = LinearDynamics(A=[[.4, .15], [.05, .25]],
                         initial_state=[1, 0], time=.1)
result = solve(problem, method=LCHS())
print(result.solution)

[ 0.96006038-1.54102084e-12j -0.00513885+3.61167323e-13j]

The matrix exponential gives \(e^{-0.1A}(1,0)^T\approx(0.96082565,-0.00484022)\), so this default finite approximation has absolute L2 discrepancy about .000821472 (measured values). Replace A, initial_state and time to solve your own system. The default output is the physical solution in the original coordinate order.

Each term \(c_je^{-iT(k_jL+H)}\) of the finite sum is a branch. The quantum circuit applies the sum as a linear combination of unitaries (LCU) block (ACL Appendix A.3, Lemma 24, Eq. (178)), in which a SELECT operation applies the branch whose index a binary address register holds. Each computational basis state of the address register is an address slot, and slots without a branch carry zero coefficient. In a flat layout, which lists the branches consecutively, the number of address slots is the branch count rounded up to a power of two. Problems without a source and the branch-controlled SELECT, the only one that dense_exact uses, have a flat layout. Other SELECT implementations of a constant-source problem use a field-padded layout instead, which pads the k-node and Duhamel-node fields separately and, when both inputs are present, adds one bit to distinguish the initial state from the source. The default plan of this example builds one nine-qubit dense SELECT with 204 branches and 256 address slots. NWQLib computes each branch matrix classically from its node's Hermitian eigensystem and inserts it as a controlled circuit, a finite construction for small systems. execution="classical" evaluates the planned finite sum classically and reports classical outputs. Larger problems need a supported structured construction with its own limit check.

A constant source fits the default limits in classical execution:

source_problem = LinearDynamics(A=[[.6, .1], [.1, .3]],
    initial_state=[1., 0.], time=.1, source=[.2, -.1])
source_result = solve(source_problem, method=LCHS(), execution="classical")
print(source_result.solution)

[ 0.9601867 -3.83384053e-17j -0.01971285-2.65777763e-18j]

It evaluates a finite LCHS and Duhamel sum. The closed-form solution \(e^{-AT}u_0+\left(\int_0^Te^{-As}\,ds\right)b\) is approximately \((0.96127277,-0.0195092)\), and the observed absolute L2 discrepancy is about .00110499 (measured values). Its quantum plan needs 1836 branches, 2048 address slots and 12 qubits. Building its SELECT takes 42,558,888 work units (NWQLib's count of planned classical operations) and 266,359,472 known peak bytes (the declared peak allowance for the numerical arrays and circuit data counted during SELECT construction), both within the default construction limits. The classical route evaluates the finite sum directly.

The linear dynamics notebook (examples/lchs_linear_dynamics_intro.ipynb) applies LCHS to advection-diffusion. It separates product-formula error from quadrature error and compares kernels, quadratures, and exact and matrix-product-state (MPS) state preparation by their errors, success probabilities and gate counts. The LCHS scientific notebook (examples/lchs_scientific.ipynb) applies a constant source and imaginary boundary penalty to four-coordinate heat flow. Running its cells performs one ten-qubit quantum solve, within the 12 qubits that the example notebooks simulate at most, and separates LCHS approximation error from finite-penalty error. The notebook is generated from examples/generators/lchs_scientific.py.

Supported problems and limits

• \(A\) and the source are time-independent. General time-dependent \(A\) or source is not implemented (open work).
• With \(A\) given as a periodic stencil, the problem must have no source.
• \(L\) must be PSD. By default (make_l_psd=True) NWQLib shifts \(L\) and restores the physical growth factor (Inputs). With make_l_psd=False, a negative eigenvalue beyond the numerical window psd_tolerance*||L||_2 is rejected.
• The default dense_exact SELECT is refused above max_dense_select_slots=4096 address slots, independently of its work and byte limits. Planning never changes the grid, duhamel_nodes or the tolerance to fit, because fewer Duhamel or k nodes would change the approximation. Each dense branch is a controlled, classically computed matrix exponential. The refusal names the time, tolerance, limit or backend choices that you can change.
• For dense \(A\) with nonzero \(L\), the spectral check accepts physical dimension at most 464 under the default max_spectral_work=100_000_000. A periodic stencil has an analytically PSD \(L\) with a known norm bound, so its compact construction does no dense spectral solve and this limit does not apply.
• approximation_tolerance and the reported kernel and quadrature bounds are component bounds. They do not bound the total physical-output error, which also depends on input scaling, PSD recovery, preparation, evolution and numerical errors.

The technical roadmap lists the current limitations of LCHS and the planned scalable SELECT.

What fits the default limits

Case	System qubits	Branches (address slots) or dimension	Default limit	Outcome
First example, default settings	1	204 (256), nine qubits in total	max_dense_select_slots=4096	fits
First example, approximation_tolerance=.001	1	396 (512), ten qubits in total	max_dense_select_slots=4096	SELECT phase fits at 7,729,128 work units and 48,453,952 known peak bytes
Constant-source example, quantum	1	1836 (2048), 12 qubits in total	max_dense_select_slots=4096	SELECT phase fits, classical route also available
dense_exact SELECT on a signed-binary k grid, dense \(A\) with generic branches	3, 4, 5, 6	at most 1024, 256, 64, 16 power-of-two branch counts	max_select_work=1_000_000_000, max_bytes=10_000_000_000, 4096 address slots	SELECT-phase upper limits, with the next branch count exceeding both work and bytes at q=5 and q=6, where bytes is reported first
Spectral check of dense \(A\) with nonzero \(L\) (none for a periodic stencil)	8, or 9 after padding	physical dimension at most 464	max_spectral_work=100_000_000	a power-of-two input reaches 256 coordinates, dimensions 257 to 464 pad to 512
Pauli decomposition of dense \(L\) and \(H\) for product formulas	at most 10	385,888,256 units at \(d=1024\), 1,686,179,840 at \(d=2048\)	max_select_work=1_000_000_000	eleven qubits are refused at this stage

Work units do not predict runtime. These stage limits set no memory aside by themselves, and a problem within them is not guaranteed to run to completion. The generic dense SELECT rows have respective work counts 756,523,008, 702,578,688, 514,477,056 and 492,054,784, and known peak bytes 3,884,752,224, 3,308,166,112, 8,597,770,464 and 1,294,063,840. They exclude the separate spectral and coefficient-preparation (PREP) checks. The figures depend on system dimension, node and branch counts and address width, and use dense_control_route="auto". Planning limits gives each stage's work.

Inputs

LinearDynamics field	Meaning
A	Finite square dense matrix, or a periodic stencil. LCHS has no sparse or Pauli route. A SciPy sparse \(A\) can be passed as A.toarray() and a Qiskit SparsePauliOp as A.to_matrix(), and the dense limits then apply
initial_state	Physical \(u\) at initial_time, including magnitude and phase
time	Final time, at least initial_time
source	Optional constant physical source \(b\). None means zero
initial_time	Initial time, default 0
time_unit	Optional label. The numerical product of A and time must be dimensionless

Known negative spectral endpoints of \(L\) follow make_l_psd. The default shifts \(L\) by \(s\) and multiplies the result by the physical growth factor, \(e^{-AT}=e^{sT}e^{-(A+sI)T}\) (Result 25). make_l_psd=False rejects a violation beyond the numerical PSD window psd_tolerance*||L||_2, the scale of the eigensolver's rounding (LAPACK Users' Guide, 3rd ed., Sec. 4.7), so the decision does not depend on the time unit. No Method option bypasses a known false PSD assumption. The lower-level cartesian_decomposition(..., eigcheck=False) only decomposes \(A\). It does not establish PSD or replace the Method's check.

For a dimension that is not a power of two, the encoded problem is \(A\oplus0\), with zero-padded initial and source vectors. PSD shifting includes the dummy coordinates coherently. The original spectral endpoints and the analytic zero endpoint determine the padded problem, so no second eigensolve is needed. Results are projected back to the original coordinates. Constant sources keep their Duhamel times, weights, physical input norms and elapsed-time PSD compensation.

An exactly zero \(L\) has known spectral endpoints (0,0), so its eigensolve and the cubic-work limit check of that eigensolve are skipped. A near-zero nonzero \(L\) follows the ordinary spectral and PSD rules. Small nonzero dissipation remains a quadrature problem, while exactly zero \(L\) uses the unitary identity.

When time equals initial_time, there is no evolution. Known input vectors and norms are used directly. A requested quantum observable or sample still uses its state preparation and readout. With classical execution, an observable of an initial state given as an explicit vector is computed once, and analysis and loading a saved result reuse that value. A zero initial state given in a compact representation rather than as a vector has no normalized state to prepare. Normalized outputs are then undefined, and a requested physical solution is formed once as a zero array, within its limits, and saved.

Outputs

Output	Value	Notes
Solution (default), StateVector(normalization="physical")	\(u(T)\) in the original coordinates	Keeps the physical phase and norm
StateVector(normalization="unit"), NormSquared, QuadraticForm, NormalizedExpectation	The state direction or a scalar of \(u(T)\)	Each comes from its own readout. Normalized outputs are unavailable when the physical mass is zero
Samples	Measured original-coordinate indices and counts	Excludes dummy coordinates. Counts do not recover a complex vector

With shots=None (exact readout), a scalar observable is evaluated from one simulation of the circuit without shot sampling, so it has no finite-shot error. Circuit, observable-action and normalization rounding remain separate error contributions. With finite shots, Pauli terms that commute qubit-wise share a measurement basis, and shots is the number of shots per group. Readout gives the details, and readout tolerance branches explains how the masses of an exact readout are checked.

plan(...) returns a Plan, the chosen construction and its costs, computed before any circuit exists. The fields coefficient_l1_norm, physical_branches and padded_branches of plan.reconstruction describe its coefficient table. These fields differ between the classical and quantum Plan objects of a constant-source problem. The classical Plan of the example above evaluates the initial state and each of the eight Duhamel nodes with the same 204 kernel coefficients, and it reports their dimensionless one-norm 1.40368 and 204 branches. A quantum Plan prepares all 1836 branches in one SELECT, each kernel coefficient multiplied by the initial-state norm or by a Duhamel weight times the source norm, so its one-norm, approximately 1.43507 at the defaults, is in the solution's unit.

Each discretization bound is reported in two forms:

Name	Where	Meaning
kernel_approximation_bound, quadrature_bound	Fields of plan.reconstruction	Operator norm of the error for a unit input vector, before the PSD growth factor. The Result summary labels them "per unit input before PSD growth"
kernel_approximation, k_quadrature	Entries of plan.facts (values that planning attaches to the Plan, here error bounds), and of result.facts for a classical Result	Physical L2 bounds. Each multiplies the field by the input norm, the Duhamel weight and exp(psd_shift*t) of each operator application (the propagator applied to the initial state, or to the source at one Duhamel node) and sums over the applications

For initial_state=[20, 0] with the first matrix, the facts are twenty times the fields. The kernel term is the error of replacing the exact propagator by the finite kernel integral. For the ACL Eq. (7) kernel (arXiv:2312.03916v2) this is the cutoff tail. For the Low–Somma kernel it also includes the approximate-identity term, independent of the cutoff \(R\), of Theorem 1, Eq. (8), so increasing the cutoff alone does not remove that second contribution.

plan(...) keeps the nodes, coefficients, PSD conversion, input directions and numerical products it chose. It does not compute a solution or run a reference comparison. estimate(plan) adds up the resource formulas of the planned circuit without building it. Numerical errors whose bounds are unavailable remain unknown.

Selection and accuracy

The default kernel is near_optimal_eq7 with beta=.75, composite Gauss quadrature, truncation_multiplier=1 and approximation_tolerance=.01. Beta lies in ACL's practical range in Sec. 2.2. The construction gives half the tolerance to the cutoff tail and half to the k quadrature. It chooses \(K\) by inverting the finite tail bound from ACL Eq. (186), independently of the panel width. An explicit multiplier of at least one can enlarge this cutoff.

The closed tail bound is \(2e^{-x}/(\beta C_\beta x)\) with \(x=\cos(\beta\pi/2)K^\beta\) and \(C_\beta=2\pi e^{-2^\beta}\), the normalization constant of the kernel \(g_\beta\), from ACL Eq. (186) and \(E_1(x)=\int_x^\infty e^{-v}\,dv/v\le e^{-x}/x\) (Proposition 26 gives the derivation). For \(K\ge1\), ACL Eq. (62) may tighten it.

For the default kernel–quadrature pair, Q-point Gauss panels are selected with the Bernstein-ellipse bound of Trefethen, ATAP Theorem 19.3 (ISBN 978-1-61197-239-9), Eq. (19.8), whose index is \(n=Q-1\). ACL's own panel prescription (arXiv:2312.03916v2, Lemma 11, Eqs. (64)–(65)) assumes \(T\max\|L\|\ge32/e\). The first example has \(T\|L\|=.045\), so that assumption fails there, while the ellipse bound holds for every \(T\|L\|\ge0\). Let \(m\) be the number of panels on each half of \([-K,K]\), so each panel has width \(h=K/m\). With margin \(d_0=.1\) from the boundaries of the analytic strip \(|\operatorname{Im}k|<1\), the ellipse lies in \(|\operatorname{Im}k|\le b_0\) with \(b_0=1-d_0=.9\). Its parameter \(\rho\) satisfies \(\log\rho=\operatorname{asinh}(2b_0/h)\), and the bound is

\[ \frac{64K\,e^{b_0T\|L\|}}{15C_\beta d_0(\rho^2-1)\rho^{2(Q-1)}}. \]

It holds for noncommuting \(L\) and \(H\) and has no dimension factor (Proposition 27). The finite scalar search chooses the fewest nodes that satisfy this bound and builds only that grid. For each successive number \(m\) of panels per half-axis, it chooses the least number \(Q\ge2\) of Gauss nodes per panel that satisfies the bound, giving \(2mQ\ge4m\) nodes in total. Once \(4m\) exceeds the smallest total node count already found, neither this nor any larger \(m\) can improve that count, so the search stops. The field quadrature.h1 of the coefficient table that resolve_lchs_coefficient_plan returns is that grid's panel width. Rules other than Gauss have no panel width. Bounds are binary64 evaluations, not interval arithmetic. A positive bound that underflows is rounded upward and never reported as an exact zero.

For the first example, \(K\) is about 27.51574, the tail bound is .005 and the quadrature bound is about .00431140. With this unit input and no PSD shift, the facts kernel_approximation and k_quadrature of the Plan equal these unit-input bounds. The observed absolute L2 discrepancy is about .000821472. Changing the tolerance to .001 selects 396 branches, with tail and quadrature bounds about .0005 and .000450744. Its classical finite sum has discrepancy about .0000532605. That finer dense SELECT needs 512 address slots and ten qubits and fits the construction limits. The two component bounds do not bound the total physical-output error, which also depends on input scaling, PSD recovery, preparation, evolution and numerical errors.

The quadrature bound can be conservative even when its numerical reference resolves the component error. For \(L=.5I+.2Z\), \(H=.3X\), \(T=1\), \(\beta=.75\) and tolerance .01, the 216-node rule uses \(K=27.515743\), nine panels per half-axis and 12 points per panel. Its finite Gauss sum differs from an adaptive integral over \([-K,K]\) by about \(1.41\times10^{-10}\) in operator norm, against the quadrature bound .00416129. The adaptive reference's error estimate is \(1.65\times10^{-13}\), which is numerical resolution evidence rather than a deterministic enclosure. The same finite sum differs from \(e^{-T(L+iH)}\) by .000731460, against the combined tail and quadrature bound .00916129. The finite-interval comparison isolates quadrature error, while the propagator comparison also includes the cutoff approximation.

The component bounds describe the finite sum in exact arithmetic. The norms of its terms add up to \(\alpha\) times the input norm, where \(\alpha\) is the coefficient one-norm, while the exact solution of the shifted problem is no larger than the input norm. A large \(\alpha\) therefore means cancellation, and binary64 evaluation carries rounding of order \(u_{\rm fp}\alpha\) times the input norm, with unit roundoff \(u_{\rm fp}=2^{-53}\). Because LCHS reports floating-point error as unknown, planning refuses a coefficient table whose \(u_{\rm fp}\alpha\) exceeds 0.1 times approximation_tolerance, the error budget of each construction stage. This affects Low–Somma with a large shift \(c\), whose \(\alpha\) is close to \(e^c\operatorname{erfc}(1/(2\gamma))\), with \(\gamma\) the kernel width that Theorem 2 selects, about 1.2e16 at \(c=50\).

Low–Somma's node savings can increase coefficient mass. At the same \(T=\|L\|=1\), tolerance \(10^{-4}\) and symmetric trapezoidal rule, \(c=1\) uses 203 nodes with \(\alpha\approx2.26846\), while \(c=5\) uses 79 with \(\alpha\approx54.6904\). Their \(1/\alpha^2\) scales are .194329 and .000334332. For an ideal LCU with unit input and contractive target \(x\), a finite vector error \(\|\widehat x-x\|\le\eta\) implies \(p_{\rm succ}\le\min(1,(1+\eta)^2/\alpha^2)\). The scale is not an equality, and sampled scalar discrepancies do not establish the vector-error assumption uniformly over matrices.

Error components fall into two classes. Discretization error (the kernel-integral error, k quadrature and, with a source, Duhamel quadrature) separates the exact propagator from the finite sum. Construction error (product-formula or quantum signal processing (QSP) synthesis and approximate state preparation) separates that finite sum from the block the circuit applies. Each stage becomes a physical-vector bound by its own rule. The kernel-integral, k-quadrature and product-formula stages are weighted per application by the input norm, the Duhamel weight and exp(shift*elapsed). QSP synthesis and coefficient preparation are scaled by the recovery scale of the prepared coefficient state. The Duhamel and input-preparation stages are already physical bounds. The available stages are summed once and then converted to the requested output, and each component is reported with its own value or the reason it is unavailable.

For NormalizedExpectation, a complete physical L2 bound \(\delta\), physical norm \(r>\delta\ge0\) from exact readout and finite \(w\ge\|O\|\) give the sharp bound \(2w\delta/r\) for a Hermitian observable (Proposition 30). Missing component, norm or observable evidence leaves that conversion unknown. A sampled norm estimate does not supply the radius required by this deterministic conversion. This exact-arithmetic relation is evaluated with ordinary scalar rounding, and the floating-point errors of the computed vector are treated separately.

Planning fixes the inputs of the kernel-integral, k-quadrature and budgeted product-formula bounds, so plan(...).facts already holds their values. A classical Plan lists these stages, and its Result records the same values as sums over the executed applications. A quantum Plan also lists the preparation stages and, for QSP, the QSP synthesis stage. For a dense \(A\), the synthesis stage of fixed-step trotter and the Duhamel stage of a nonzero constant source stay unknown on both execution routes until LCHSRefinement evaluates them. Request them after the solve with result.verify(checks=LCHSRefinement(components=("duhamel",))) for the Duhamel stage or components=("fixed_pf",) for the product-formula stage, with LCHSRefinement imported from nwqlib.algorithms.lchs. Refinement evaluates them at the saved nodes and step counts and chooses no new grid or step count. The one exception is quantum qsp_block_encoding with \(A=0\), whose constant Duhamel integrand gives an exactly zero Duhamel stage.

The function compare plans several configured LCHS Methods for one Problem and estimates their resources before execution. Kernels compared on identical signed-binary nodes isolate the kernel change. Low–Somma selects its own symmetric trapezoidal grid from its tolerance, so a comparison that includes it changes the kernel–quadrature pair. Comparisons at fixed construction settings do not measure cost at matched accuracy. The MPS analysis studies compression and circuit fitting for one kernel–quadrature pair.

Error budget for heat flow on a ring

For heat flow on a ring, plan(...).facts gives every error bound before anything runs, and only the cutoff tail, the k quadrature and the Strang product have nonzero bounds. The problem is \(du/dt=-Lu\) with \(L=\tfrac14(2I-S-S^\dagger)\) on \(2^q\) sites, where \(S\) shifts the site index by one, from a unit point source at site 0 until \(T=1\), with \(L\) given as a periodic stencil. approximation_tolerance=0.005 is split equally between the cutoff tail and the k quadrature, and trotter_synthesis_tolerance=0.005 sets the Strang step count. With the default kernel–quadrature pair, the three bounds together bound the ideal L2 error of the physical vector by 0.01 for this unit input, before rounding in the circuit and sampling.

import numpy as np
from scipy.linalg import expm
from math import fsum
from nwqlib import LinearDynamics, plan
from nwqlib.algorithms import LCHS
from nwqlib.algorithms.lchs import resolve_lchs_coefficient_plan
from nwqlib.operators import PeriodicStencil
from nwqlib.problems import ingest_occupation

q = 3  # eight coordinates for the ideal finite-branch comparison
heat = LinearDynamics(A=PeriodicStencil(q, mass=0.0, diffusion=0.25),
                      initial_state=ingest_occupation("0" * q, num_qubits=q),
                      time=1.0)
method = LCHS(hamiltonian_evolution_backend="trotter",
              approximation_tolerance=0.005, trotter_steps=None,
              trotter_synthesis_tolerance=0.005)
selected = plan(heat, method=method, seed=7)
for item in selected.facts:
    print(item.fact.quantity, item.fact.value.value)

coefficients = resolve_lchs_coefficient_plan(selected)
r = selected.reconstruction.step_counts[0]
n = 2**q
I = np.eye(n)
S = np.roll(I, 1, axis=0)
E = I[np.arange(n) ^ 1]
O = S @ E @ S.T
terms = []
for k, c in zip(coefficients.nodes, coefficients.coefficients, strict=True):
    angle = k * .25 / r
    half = np.cos(angle / 2) * I + 1j * np.sin(angle / 2) * E
    full = np.cos(angle) * I + 1j * np.sin(angle) * O
    step = half @ full @ half
    terms.append(c * np.exp(-.5j * k) * np.linalg.matrix_power(step, r)[:, 0])
ideal = np.array([complex(fsum(v[j].real for v in terms),
                         fsum(v[j].imag for v in terms)) for j in range(n)])
exact = expm(-(2 * I - S - S.T) / 4)[:, 0]
print(np.linalg.norm(ideal - exact))

kernel_approximation 0.002499999999999994
k_quadrature 0.002304654032281964
trotter_synthesis 0.004970884480827643
lcu_coefficient_preparation 0.0
initial_state_preparation 0.0
0.00026603733009987544

The coefficient and initial-state preparation bounds are zero because both preparations are direct. The last line is the numerical L2 discrepancy of the ideal finite-branch construction from the heat-equation solution. It sums the 270 planned branches with 8×8 matrices and compares the physical vector with scipy.linalg.expm(-L)[:, 0]. It includes finite kernel, quadrature and Strang approximations. Circuit formation, execution and sampling errors are outside this comparison. The value is about .000266037 at 23 steps, evaluated with Python 3.12.14, NumPy 2.5.2 and SciPy 1.18.1 on macOS arm64.

Quantity	Value	How it is obtained
Kernel	ACL Eq. (7) kernel near_optimal_eq7, \(\beta=0.75\)	Default
Cutoff \(K\)	32.1994	Inverts the tail bound of ACL Eq. (186) at the tail's share, 0.0025
Gauss rule	18 panels of 15 points	Fewest nodes whose quadrature bound meets the quadrature's share, 0.0025
Branches (address slots)	270 (512)	One branch per node. Address slots are the branch count rounded up to a power of two
Coefficient 1-norm	1.40522	Sum of the branch-weight magnitudes, not a success probability
Cutoff tail bound, kernel_approximation	0.0025	Closed-form bound from ACL Eq. (186), evaluated in binary64
Quadrature bound, k_quadrature	0.00230465	Bernstein-ellipse bound of Trefethen, ATAP Theorem 19.3, evaluated in binary64
Strang coefficient \(B\)	2.62960	Weighted matching-commutator coefficient for \(q\ge3\)
Strang steps \(r\)	23	Smallest \(r\) with \(B/r^2\le0.005\)
Strang bound, trotter_synthesis, at \(r-1=22\) and at \(r\)	0.00543305 and 0.00497088	\(B/r^2\) in exact rational arithmetic, recorded rounded upward. 22 steps exceed 0.005
Coefficient and initial-state preparation bounds	0 and 0	Direct preparation
Sum of the bounds	0.00977554	Ideal L2 error bound of the physical vector for this unit input, against the target 0.01
Numerical L2 discrepancy at \(q=3\)	0.000266037	Ideal finite-branch action evaluated with 8×8 matrices, compared with \(e^{-L}e_0\)

The diamond marks the numerical discrepancy of the ideal action above. docs/scripts/error_budget_figures.py draws the figure from this calculation. Its ordinary invocation runs only the LCHS calculation, and --qls separately requests the QLS examples.

For \(q\ge3\), the branches, step count and bounds depend on \(T\lVert L\rVert=1\) and the diffusion coefficient and are independent of ring size. The diffusion-only matchings commute at \(q=1,2\), where the ideal synthesis component is zero and one step suffices. Planning and resource estimation build no circuit and extend the example to rings of up to \(2^{100}\) sites:

System qubits \(q\)	Total qubits	Branches (address slots)	Strang steps \(r\)	Sum of the bounds	CX upper bound per attempt
3	12	270 (512)	23	0.00977554	37,686
10	19	270 (512)	23	0.00977554	46,518
20	29	270 (512)	23	0.00977554	74,578
40	49	270 (512)	23	0.00977554	185,898
100	109	270 (512)	23	0.00977554	961,458

Each row plans the problem above at that q and counts its CX gates with estimate(selected, context=ResourceContext(basis="cx")) (Resource inspection). The CX upper bound counts the CX gates of one coherent attempt before gate cancellation or routing, with no repeat-until-success factor. ResourceContext is imported from nwqlib.resources.

estimate reports the T count of these plans as unavailable in the default, cx and clifford_t bases, because the coefficient preparation, its inverse and the periodic phase and step blocks have no T-count rule (Choose the gate basis). Estimate fault-tolerant resources gives T counts of small compiled circuits.

Measured values

Value	Quantity	Problem and settings	Environment
.000821472	Absolute L2 discrepancy against expm	First example, default LCHS(), quantum execution on the default local Aer statevector	Python 3.12.14, NumPy 2.5.2, SciPy 1.18.1, Qiskit 2.6.0, Aer 0.18.0, macOS arm64
.0000532605	Absolute L2 discrepancy against expm	First example, LCHS(approximation_tolerance=.001), execution="classical"	Same
.00110499	Absolute L2 discrepancy against closed_form	Constant-source example, default LCHS(), execution="classical"	Same
up to 4.2e-5	Entry error of Qiskit's two-qubit synthesis of a unitary near the identity	60 random generators (dependency issues)	Qiskit 2.5.2, mechanism reproduced with 2.6.0

Evolution backends

Select the branch evolution with LCHS(hamiltonian_evolution_backend=...):

Backend	Branch evolution	Step count
dense_exact (default)	Controlled dense branch matrix from the node's Hermitian eigensystem	None, no product formula. Accepts only trotter_order=2
trotter	Fixed-step sorted-Pauli Lie-1 or Suzuki-2 product formula	trotter_steps, or on the periodic route trotter_synthesis_tolerance
trotter_error_budgeted	Product formula with the smallest step count per node that fits an error budget	Chosen per node. Requires trotter_steps=1
qsp_block_encoding	QSP evolution of a joint generator built from L and H block encodings	None, no product formula. Accepts only trotter_order=2

trotter_order selects the Lie-1 or Suzuki-2 formula of either product-formula backend. psd_tolerance is the numerical PSD window. Kernel and k-grid choices use lchs_kernel and k_quadrature. The LCHS API explains the MPS state-preparation controls and the work and byte limits. Those limits do not bound the total physical error.

dense_exact and qsp_block_encoding do not use trotter_steps. The Method sets that field to one before computing its content hash, so Methods that differ only in an ignored step count are the same Method. A nondefault trotter_steps is ignored with a UserWarning at planning. A copy made with model_copy() and a loaded Method hold trotter_steps=1 and no record of the supplied value, so planning them gives no warning (saved run mechanics).

When you change the backend of an existing Method to "trotter", set trotter_steps in the same change. Changing only the backend keeps trotter_steps=1. trotter_error_budgeted chooses per-node counts from its separate error budget, based on approximation_tolerance.

• dense_exact computes each small branch matrix from its node's Hermitian eigensystem and controls the corresponding circuit on the address bits. dense_control_route selects how the address bits control a branch on two or more system qubits. The default "auto" synthesizes the whole controlled branch for one address bit and controls each synthesized gate for more (dense control route). On the gate-wise route, before adding the address controls, NWQLib replaces the dense branch unitary by a circuit of CX and one-qubit gates that equals it to rounding (subroutines/_dense_synthesis.py). Generic control would otherwise use Qiskit 2.6.0's own synthesis, which replaces a two-qubit unitary close to the identity by a local gate. A branch \(e^{-it(kL+H)}\) is close to the identity when \(t(kL+H)\) is small, for example after a short elapsed time, and for such a branch that replacement gave entry errors measured up to 4.2e-5 (measured values). Controlled dense unitaries describes the construction, and Planning limits its cost.
• trotter uses fixed-step sorted-Pauli Lie-1 or Suzuki-2 evolution. Its order and per-node repetitions are saved.
• trotter_error_budgeted gives each node the smallest step count \(r\) for which \(r\) one-step bounds at step length \(t/r\), plus the reported pruning bound \(tm\) for an upper bound \(m\) on the pruned Pauli mass, fit the per-node error budget of 0.1 times approximation_tolerance, with \(t\) the node's elapsed time (CSTWZ doi:10.1103/PhysRevX.11.011020, Sec. V B). It does not read trotter_steps, so a Method that pairs it with a value other than the default 1 is refused. The one-step bound is CSTWZ Proposition 9, Eq. (120), for trotter_order=1 and Proposition 10, Eq. (121), for order 2, with each commutator norm replaced by its triangle-inequality sum over Pauli terms. Each one-step bound is a constant times the square (order 1) or cube (order 2) of the step length, so order-1 step counts grow as the inverse of the per-node error budget and order-2 counts as its inverse square root. The budget bounds only the product-formula synthesis error of each node.
• qsp_block_encoding chooses the L and H block encodings, joint-generator diagonal weights, polynomial and phases once. Execution and loading a saved result reuse them. The success projector also covers its synthesis ancillas. An exactly zero joint generator gives the identity evolution directly, without a zero block-encoding scale or a phase fit.

lcu_select_implementation="auto" uses only the SELECT implementations that apply to the problem, with their cost rules. An explicit structured, multiplexor or branch_controlled choice is refused where it does not apply, and the lcu_select_implementation entry of the LCHS API lists where each applies. structured requires proof of an eligible affine structure, computed from the k-node addresses and the generator's real Pauli coefficients.

Periodic stencils

For \(A\) given as a PeriodicStencil, planning keeps the compact input. This route has a homogeneous Strang construction only, so a problem with a source is refused (periodic.plan_periodic). The Strang path keeps the one-qubit double-wrap edge convention. For \(a_j=|k_j|\,\mathrm{diffusion}\) and \(p=|\mathrm{potential}|\), its weighted unit-input synthesis bound is \(B/r^2\), where \(B=T^3\sum_j\overline c_j W_q(a_j,p)\), \(\overline c_j\ge|c_j|\) and \(W_q(a,p)=g_q a^3/2+4a^2p/3+ap^2/3\). Here \(g_q=0\) for \(q=1,2\) and one for \(q\ge3\). The proof follows the order in which the circuit applies the potential, even-matching and odd-matching factors (Proposition 36). The weighted bound concerns the finite physical action on a unit input. trotter_steps fixes \(r\). LCHS(hamiltonian_evolution_backend="trotter", trotter_steps=None, trotter_synthesis_tolerance=eps) instead selects the smallest \(r\) with \(B/r^2\le\) eps, then applies max_trotter_steps, max_bytes and max_select_work. A Method that sets both count and tolerance, or neither, is refused. Other operators refuse this tolerance because their fixed-step routes have no planning bound to invert.

For heat flow with PeriodicStencil(q, mass=0.0, diffusion=0.25), \(T=1\) and both tolerances .005, \(q\ge3\) gives \(B\approx2.62960\) and \(r=23\). At \(q=1,2\), the matchings commute and \(B=0\), so automatic selection chooses one step. Nonzero potential keeps both cross terms at every width, and zero diffusion gives \(W_q=0\) for every potential, while the step selection still requires mass+4*diffusion>0. A zero ideal product error leaves numerical formation, execution and readout errors separate. The two tolerances together bound the ideal vector error of a unit input, before circuit rounding and sampling. Error budget for heat flow on a ring shows the three components against the .01 target.

A periodic construction refusal reports its work and byte requirements alongside max_select_work and max_bytes. Raise the limiting budget or reduce the construction. A smaller fixed trotter_steps, or a larger trotter_synthesis_tolerance in automatic mode, can reduce the step count and construction cost, while increasing or preserving the Strang error bound. If even one step cannot fit, change the construction or its limits. In automatic mode, a step count above max_trotter_steps is refused before the construction limits are checked.

Explicit checks and saved results

from nwqlib.algorithms.lchs import LCHSVerification
from nwqlib.evidence import Certificate

checks = LCHSVerification(reference="expm", metric="absolute_l2", threshold=.01)
receipt, facts = result.verify(checks=checks)  # The only reference computation.
print(facts[0].fact.value.value)  # This comparison emits a Float64 discrepancy.
certificate = Certificate(plan_id=result.plan_id, result_id=result.content_id,
    assessment=result.assess(absolute_tolerance=.01), checks=())
certificate = certificate.with_verification(
    result, options=checks, evidence=facts)
print(certificate.checks[0].status)  # PASS: .000821... is below .01.
result.save("lchs-result")

0.0008214720329548587
PASS

Checks are explicit extra computations. Choose the reference with reference:

reference	Compares the result with	Applies to
expm	\(e^{-AT}u_0\)	Problems without a source
ivp	An RK45 solve of \(du/dt=-Au+b\), with explicit rtol, atol and max_rhs_calls	Problems with a constant source
closed_form	\(e^{-AT}u_0+\left(\int_0^Te^{-As}\,ds\right)b\), both terms of ACL arXiv:2312.03916v2, Eq. (2)	Problems with a constant source
ivp_closed_form	Both of the above, with their signed consistency	Problems with a constant source
selected_grid	The exact saved finite sum, including the saved product-formula steps, so its discrepancy excludes the kernel-integral, k- and Duhamel-quadrature and product-formula errors	Problems with or without a source

Every reference needs a dense physical A and never densifies one implicitly. closed_form evaluates one exponential of the augmented matrix [[-A*T, b*T], [0, 0]], whose top-right block is the Duhamel term. It equals \(A^{-1}(I-e^{-AT})b\) for invertible \(A\), is defined for singular \(A\) as well and avoids the cancellation of forming \(I-e^{-AT}\) when \(\|A\|T\) is small. The expm and closed_form references use scipy.linalg.expm, or scipy.sparse.linalg.expm for a triangular matrix, because the dense routine scipy.linalg.expm loses accuracy on triangular input with close diagonal entries (dependency issues). A reference whose matrix has a one-norm above 2**37, or whose exponential is not finite in binary64, is reported as unknown. The quantum selected_grid reference evolves the input vector of each operator application through the decomposed branch circuit. Reference discrepancies are empirical numerical evidence, not a complete proven error bound. Reference work gives the work each reference counts.

LCHSRefinement requests component-bound analysis with its own explicit limits. Its default max_structural_work=1_000_000_000 covers the shared Pauli commutation table, which lists the anticommuting Pauli label pairs and the triples with a nonzero nested commutator over which the triangle-inequality commutator bound sums (Product-formula planning), and its contractions at each node, while max_dense_work=100_000_000 separately limits requested dense validation. Its verification record contains all observations, while the returned facts contain only output-error facts suitable for assessment, so a request for "spectral_norms" alone leaves facts empty.

result.verify returns the verification record receipt and the error facts facts. The check assessment above applies its threshold to the same facts and Result, without calling verify again. These facts apply only to the options in checks. Another threshold, metric or reference needs result.verify(checks=...) with the new options, which computes the reference again (Check accuracy and verify a result). The aggregate accuracy assessment can independently remain INCONCLUSIVE. See Check against a tolerance.

load_result(...) restores the original Method, Plan, computed output, saved phases, MPS cores and numerical data. The saved plan names a separate JSON file with NPY arrays for the data that grows with the grid and schedule: node and coefficient tables, LCHSData.host_actions (classical product-formula evaluation) and SELECT angle tables. Loading gives each circuit block the same data that planning gave it. SELECT uses those tables and the Problem's elapsed time, the periodic Strang blocks use the saved periodic Strang data, the input preparation uses the Problem's initial state, and a state-preparation block uses its saved tensor. Loading does not fit phases, solve an eigensystem, recompute a reference or reselect a grid.

Coefficient MPS analysis

from nwqlib.algorithms.lchs import resolve_lchs_coefficient_plan

coefficients = resolve_lchs_coefficient_plan(problem)
decomposition = coefficients.decompose_mps(max_bond_dim=4)

resolve_lchs_coefficient_plan(problem) returns the coefficient table of the finite sum for a problem without a source, and resolve_lchs_coefficient_plan(plan) the table that a Plan of such a problem executes. Neither call builds a circuit. decompose_mps returns an MPSDecomposition with the TT (tensor-train) cores, the bond dimensions and the discarded weight, the sum of the squared singular values discarded across the TT-SVD sweep of the normalized coefficient preparation tensor. To prepare the coefficient state with an MPS circuit in a solve, pass LCHS(lcu_state_preparation="mps_circuit", lcu_mps_max_bond_dim=...).

The coefficient state-preparation target is the ordered normalized magnitude tensor sqrt(abs(c) / sum(abs(c))). Complex coefficient phases stay in SELECT. A kernel, quadrature or tolerance choice passed to resolve_lchs_coefficient_plan together with a Plan is refused, so another coefficient table requires a new Plan. Positive scaling permits reuse of identical normalized tensor data, while a decomposition passed as reuse for a reordered or changed tensor raises ValueError. TT compression error and layered MPS circuit fidelity are different quantities. Another kernel–quadrature pair has its own coefficient tensor and needs its own decomposition. Constant-source quantum construction supports direct state preparation only. MPS state preparation there would need its own decomposition of the source-branch coefficient tensor and of both input directions.

Finite-layer MPS state preparation may change the LCU block even when the TT-SVD discarded weight is zero. Its circuit error remains unevaluated until an explicit validation of that circuit, and the summary that print(result) shows states this scope.

That summary reports the saved MPS layer count, largest bond dimension, core bytes and construction work. Circuit fidelity validation still performs an explicit statevector simulation, and explicit MPS reference analysis may reconstruct the tensor. Printing the summary calls neither. max_svd_work bounds the prospective TT-SVD and the layered circuit construction, whose work formula mps.layered_construction_size follows the scikit_tt sweeps that mps_circuit.mps_to_circuit repeats after every extracted gate. The direct MPS builder checks both before it starts. LCHS planning checks the layered construction of each mps_circuit state preparation against the Method's max_svd_work and max_bytes, the limits that building the circuit applies again. The internal workspace of scikit_tt remains unknown.

Resource inspection

estimate(plan) adds up the resource formulas of the plan. For a separate count on representative built circuits:

from nwqlib.algorithms.lchs.quantum_resources import sample_resources

samples = sample_resources(result.plan, max_qubits=13)

samples["representatives"] lists each representative circuit with its name, its multiplicity in the construction and its operation counts. samples["weighted_totals"] gives, for each readout setting (a sampled group, a mass setting or the exact reduction), the operation counts of every core representative (every representative other than a readout representative) multiplied by its multiplicity, plus the counts of that setting's readout representative. A representative's operation names stand for every block it represents, so a weighted count under a node-specific name can cover several nodes. Weighted totals are estimates. They do not model angle-specific cancellation, cross-block optimization or the working memory of Qiskit, and they do not combine depths.

A dense exact SELECT builds no representative. Its samples["structural_bounds"] record uses the value-independent upper counts of the library synthesis and control route for every branch. Each record names its basis, route and assumptions, and has count_kind="structural_upper_bound". It is not a count of an observed circuit, and the weighted totals do not include it. Counts for any explicitly requested compilation require inspection of that compiled circuit. For a problem with a source or a compact periodic stencil, inspect the full circuit with prepare(plan).inspect_resources(index=...). The library function inspect_circuit_resources runs no simulation and transpiles only when asked to. Import the inspection function with from nwqlib.backends.inspection import inspect_circuit_resources. estimate can instead count CX gates with ResourceContext(basis="cx"), as shown in resources, without preparing a circuit. No resource calculation establishes solution accuracy.

The compiled SELECT of qsp_block_encoding is one QSP Hamiltonian evolution of the joint generator built by build_control_diagonal_generator_encoding. One iteration of oblivious amplitude amplification uses the combined cosine and sine block twice and its adjoint once, separated by two reflections. With cos_degree and sin_degree denoting the cosine and sine polynomial degrees, the evolution therefore makes 3*(cos_degree + sin_degree) queries to the joint block encoding, and each query costs one call per child encoding plus one multiplexed RY of 2**control_qubits rotation angles per child, where control_qubits is the width of the address register. These counts are structural, rotation angles and child calls before gate synthesis, so they do not depend on the Qiskit version. In Pocrnic et al., arXiv:2506.20760v2, Lemma 7, p. 15, each effective-Hamiltonian query uses one U_A plus one U_A^dagger call and 2M multi-controlled rotations, where M is the number of quadrature nodes. Here each query uses child_count child calls and child_count * P multiplexed rotation angles. child_count is 2 for the \(L\) and \(H\) branches and 1 when \(H=0\), and P = 2**control_qubits is the number of address slots after power-of-two padding, so P >= M. With two children, the 2P multiplexed rotation angle entries per joint-generator query are at least the paper's 2M multi-controlled rotation calls per effective-Hamiltonian query. With one child, NWQLib uses one child call and P angle entries per query. These are counts of different rotation constructions and child encodings, so the comparison does not order their synthesized gate costs. The paper's generalized QSP (GQSP) and qubitization query count and amplification constants are not part of this count.

Planning limits and numerical scope

The Method's work and byte limits have the defaults max_select_work=1_000_000_000, max_spectral_work=100_000_000, max_quadrature_work=100_000_000, max_svd_work=100_000_000, max_readout_work=2_000_000_000 and max_bytes=10_000_000_000. Work units count planned classical operations and do not predict runtime.

max_admission_steps (default 1,000,000) limits the work of checking each planned circuit description and does not change the quantum operations. If planning refuses with this field, raise it at least to the reported count. A refused count of the fields of the circuit description is complete, while a refused count of checking work is a lower bound that may need raising further, and later preparation can need more in either case (planning work limit).

Quadrature selection

For the composite Gauss rule, planning first completes the scalar panel search. It then checks max_quadrature_work against the search work, 32 units per chosen node and the cubic Gauss-rule construction work, before allocating the node arrays. A refusal names the complete work of that rule. The limit therefore covers the search plus the construction of the chosen rule, not every scalar probe.

Dense SELECT construction

On the gate-wise route, the CX formula of the dense_exact backend bounds the gate count of the exact synthesis and counts each gate at Qiskit's CX count for that gate with the address controls added. On the whole-matrix route it is the CX bound of the whole-matrix synthesis on the system qubits and the address bits. Planning counts the classical work and bytes of each node eigensystem and branch matrix (Engineering constants), of each branch synthesis and of Qiskit's control of the synthesized gates with the address bits against max_select_work and max_bytes (exact synthesis check). This work depends on the system dimension, the node and branch counts and the address width, not on the norm of the branch generator. What fits the default limits gives the resulting branch counts for generic dense branches. A dense-dilation child of qsp_block_encoding is counted the same way, and dense_control_route also selects how the joint generator's combine qubit controls it. With one address bit the whole-matrix route counts the whole-matrix synthesis on the system qubits plus one instead.

The dense_exact route computes one Hermitian eigensystem for each distinct \(kL+H\) and reuses it for every elapsed time, including the source nodes. Classical evaluation applies the eigenphases to vectors without forming each dense propagator. A quantum SELECT forms one elapsed-time branch matrix at a time and synthesizes it within its work and byte limits. Its memory count includes the eigensystem factors in memory, matrix-formation workspace, synthesis and control workspace and completed branch circuits. Classical vector actions and explicit branch matrices have different operation counts, and the resource record counts the operations each performs. PSD recovery uses each application's elapsed time. "Dense exact" names the ideal matrix function. Its binary64 evaluation and circuit synthesis carry eigensolver, phase-evaluation, contraction and synthesis errors, which are separate from the real-arithmetic LCHS approximation budget.

Product-formula planning

Dense product-formula planning decomposes \(L\) and \(H\) once with the recursive Pauli block transform. Its work and output conversion grow as \(qd^2\) for \(d=2^q\). At max_select_work=1_000_000_000, decomposition sets aside 385,888,256 units at \(d=1024\) and 1,686,179,840 at \(d=2048\), accepting at most ten system qubits at this stage. At ten qubits, decomposition plus conservative second-order relaxed SELECT coefficient selection allows a union of 1231 nonidentity Pauli labels for 200 distinct nodes, with either 256 or 4096 address slots. These figures assume no identity coefficient, \(p\) entries in each L/H support and the full pair upper bound. They include support alignment and the work set aside for dropped mass. Angle-table work, bytes and circuit construction remain separate limits on a complete Plan.

Product-formula selection builds one Pauli commutation table over the ordered Pauli labels, shares it across k-nodes and evaluates each node's bound coefficient once. Initial and source applications reuse that coefficient at their own elapsed times, pruning bounds and step counts. The full Pauli-triangle expression is used when its remaining work and storage fit. Otherwise selection uses the second-order degree-capped suffix relaxation bound_variant="relaxed_prefix", or stops before the refused stage. For anticommuting i<j it uses min(S_i,D_i+D_j-a_i), with later mass S_i and weighted anticommuting degree D_v. Degrees are rebuilt from the rows of the Pauli commutation table restricted to each node's kept labels. The packed pair tests of a shared support are counted once, while each of the \(K_n\) distinct nodes adds 2p+4E contraction visits, with E the number of anticommuting label pairs. The byte check includes the combined-node records and the arrays in memory. Padding address slots add no nodes.

For order two, consider 200 distinct active k-nodes, nine active applications and one shared sorted union of \(p\) nonidentity labels, with no identity coefficient. Use \(p\) entries in each L/H dictionary and at every node, with full pair and triple upper bounds. The table below uses one billion work units for each stage. Fixed-step refinement uses LCHSRefinement.max_structural_work, while the classical and SELECT stages use LCHS.max_select_work. For \(P=p(p-1)/2\), \(J=p(p-1)(2p-1)/6\) and these widths, the full Pauli commutation table and its contractions cost \(201(P+J)\), and the relaxed variant, whose table holds pairs only, costs \(P+200(2p+4P)\) with its contractions. Each stage adds its preparation and application work.

System qubits	Stage	Largest \(p\), full / relaxed
20	Proven fixed-step bounds	244 / 1498
20	Classical stage	246 / 1573
20	SELECT coefficient selection (4096 address slots)	246 / 1570
9	Proven fixed-step bounds	245 / 1540
9	Classical stage	246 / 1573
9	SELECT coefficient selection (4096 address slots)	246 / 1573

These are work-only stage upper limits, before decomposition, the classical actions, SELECT angle tables and circuit construction. The full limit check also includes the data in memory and the storage of the Pauli commutation table, so these figures do not imply that a complete dense Plan fits. In particular, the current dense arrays at 20 qubits exceed the default byte limit. A homogeneous 200-node SELECT uses 256 address slots, while the field-padded layout with an initial state and eight source applications uses 4096. Padding adds address-slot records and angle storage but no k nodes. The work and byte formulas of the Pauli commutation table are applied to the shared structure and distinct node contractions, adding 9 bytes per row of the larger of the pair and reserved triple counts for one node's row selection (Engineering constants). Each caller adds its own preparation and application work.

The dropped coefficient mass is an upper bound for the terms omitted. Its product with the evolution time is reported, rounded upward, as the pruning bound. The product formula is selected against the requested error budget minus that pruning bound, rounded downward. Its formula bound and the combined pruning-plus-formula bound are reported rounded upward. These directions ensure that the exact sum of the reported contributions and the reported total do not exceed the requested error budget. Selection can require additional steps or fail when no positive representable formula error budget remains.

Classical product-formula evaluation

For a classical Pauli product formula, planning records the coefficient table once per k node and the elapsed time, step count and route per operator application. Each rotation angle is computed when the action is applied, as coefficient times elapsed time divided by steps, in that order. The limit check takes the largest of the decomposition, node-table construction and action phases, and every saved table still in memory stays counted. Vector actions use cos(a)*v - i*sin(a)*P@v. A large repetition count can instead select one dense step and binary powering. With dimension \(D\), \(r\) effective repetitions and \(R\) rotations per step, the work estimates of the vector route and the dense route are r*R*(W_action + 2*D), where W_action is the work count of one single-label Pauli action on a D-entry vector, and R*D**2 + (floor(log2(r))+popcount(r)-1)*D**3 + D**2. The vector update adds one D-entry scaling and one D-entry scaled addition per rotation, with the scalar trigonometric preparation included in W_action and in those two passes. Effective repetitions are max(1, steps) for the saved step count. Only routes whose known arrays and work fit the Method limits are eligible, and equal work chooses the vector route. These are operation-size estimates, not timing predictions.

The node table, per-node steps, routes and route work are saved in LCHSData.host_actions of the Plan. Construction work and aggregate action work are in plan.reconstruction.classical_work. result.applications holds one record per operator application with its schedule and component error facts. The classical evaluation gives an upper bound on the JSON size of these records, so the run sets aside their bytes together with the outputs. Reading a report or a saved result neither chooses routes again nor evaluates an action. Constant-source nodes use their own elapsed time and PSD recovery. The dense_exact backend instead forms one Hermitian eigensystem per node and uses it for every operator application, the nonzero initial state and each constant-source node, each at its own elapsed time.

Readout

An exact scalar readout (shots=None) fixes the method's success and physical-condition bits and excludes dummy coordinates. A Pauli observable with \(n_{\rm P}\) terms then needs \(O(n_{\rm P}N)\) classical work on the \(N\)-coordinate encoded system. The normalized expectation divides the projected quadratic moment by the physical-slice mass. A physical quadratic form also applies the method's recovery scale.

Exact readout reduces the saved state to scaled scalar statistics. The reduction scans the complete mass (the squared norm of the whole simulated state) and the success mass of that state in bounded tiles, shares the physical mass when the two sets of outcomes coincide, and evaluates the observable on the success-and-physical slice. The limit check includes the saved state, the observable's serialized, parsed and packed representations, and the workspace for gathering, scaling and grouped action. Scaled statistics carry the recovery information needed for reanalysis. Readout tolerance branches explains which tolerance checks the masses. max_select_work covers SELECT construction. max_readout_work, whose default is 2,000,000,000, covers sampled-readout grouping comparisons and exact projected-reduction input visits. The reduction's known data and workspace must also fit max_bytes before the circuit runs.

Sampled observable readout groups qubit-wise-commuting Pauli terms into shared local bases. shots is the number of shots per group. Terms in a group share outcomes, so uncertainty calculations use the weighted group outcome and its covariance. Grouping reduces the number of settings. Its variance at fixed total shots depends on the state, coefficients and allocation. A padded normalized output also measures its physical-coordinate mass in an unrotated basis. Every counts setting measures all success and system qubits into one shared classical layout, and a group decodes each label's parity from the coordinates of that label's support. The grouping comparisons count against max_readout_work. The Result records each group's labels, basis, returned shots, shots postselected on success and weighted moments. Samples holds its original-coordinate indices and counts as two integer arrays. Sampled group decoding and Samples read outcomes of at most 64 measured bits.

Arrays and numerical range

Classical observable outputs use one action of the observable in its own dense, sparse or Pauli storage, with binary scaling combined with the vector's physical recovery before producing the scalar. A representable quadratic form can therefore remain available when its normalized moment or physical norm squared is not. If a common observable scale would erase a nonzero entry, the original unscaled action is used, chosen before execution, with no second action and no implicit dense conversion. Binary64 cancellation and heterogeneous-range limitations remain, and an unrepresentable result is unavailable rather than a clipped value. Work and workspace include the observable copy and scans.

The PSD recovery exp(shift*T) can exceed binary64 while the physical solution does not, for example when strong growth acts on a tiny input. The coefficients and arrays then hold each recovery divided by 2**e, with e = floor(shift*T/ln 2), and 2**e joins the binary output scale, so a representable solution is still returned. The weighted error components need the recovery as a number and are unknown with the reason unrepresentable_psd_recovery:applications.

NWQLib forms the Cartesian parts \(L=(A+A^\dagger)/2\) and \(H=(A-A^\dagger)/(2i)\) from one conjugate-transpose workspace. The parts are private arrays and become read-only after the PSD decision and any diagonal shift. Input whose dimension is not a power of two is padded as \(A\oplus0\), with the known zero eigenvalue added to the original spectral endpoints. Sequential padding and in-place scaling bound the declared array storage by four encoded complex matrices plus the vector workspace explicitly set aside. When only the coefficient table is requested, as by resolve_lchs_coefficient_plan(problem), the spectral endpoints of the unpadded \(L\) come from the same spectral function as in planning, and the padding zeros are added analytically. Its limit check counts the known NumPy arrays and excludes the linked eigensolver's queried workspace. PSD compensation remains part of the coefficients.

Reference work

The classical selected_grid reference reads the coefficient table at most once per distinct k-node with a nonzero coefficient. Fixed-step refinement reads every distinct k-node of the table once. Reading a row of the shared table without copying it costs one logical visit, and reading the identity coefficient and the upper bound on the dropped mass costs two visits per node. The table is counted in the memory of the phase that reads it. Preparation of the Pauli commutation table, application phases and bounds, vector actions and requested dense validation have separate work counts. Fixed-step refinement evaluates the saved step counts using the saved coefficients and pruning mass. The work of the quantum selected_grid reference follows each executed gate's contraction arity, and its workspace includes state rearrangement copies.

Source map

ACL is An, Childs and Lin, arXiv:2312.03916v2. Low–Somma is arXiv:2508.19238v2. CSTWZ is Childs, Su, Tran, Wiebe and Zhu, Phys. Rev. X 11, 011020 (2021), doi:10.1103/PhysRevX.11.011020. Code paths are in src/nwqlib/algorithms/lchs/ unless they name another directory. "NWQLib" marks a derivation stated in the docstring of the function in the Code column.

Scientific step	Source	Location	Code
A = L + iH and the PSD assumption	ACL	Eqs. (3)-(4), Lemma 21 Eq. (162)	time_independent_terms.cartesian_decomposition
Numerical PSD window and shift, recovery exp(shift*T)	NWQLib, from ACL Eq. (4), window from LAPACK Users' Guide, 3rd ed., Sec. 4.7	exp(-AT) = exp(sT) exp(-(A+sI)T), window psd_tolerance*norm(L, 2)	solution_error_budget._numerical_psd_decision, providers._homogeneous_coefficient_plan
Recovery beyond binary64 carried as a power of two	NWQLib	exp(sT) = (exp(sT)*2**-e)*2**e, e = floor(sT/ln 2)	solution_error_budget.psd_recovery_exponent, solution_error_budget.psd_recovery_part
Rounding refusal for a large coefficient 1-norm	NWQLib, with Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., doi:10.1137/1.9780898718027	Eq. (4.4)	providers._admit_coefficient_rounding
Encoding A ⊕ 0 for a non-power-of-two dimension	NWQLib, from ACL Eq. (4)	L ⊕ 0 and H ⊕ 0, zero endpoint added without a second eigensolve	time_independent_terms._prepare_decomposition
LCHS identity with kernel f(k)/(1-ik)	ACL	Eq. (6), Theorem 6 Eq. (12), Eq. (60)	providers.eq7_coefficient
Near-optimal kernel and C_beta	ACL	Eq. (7), same as Eq. (32), C_beta from Eq. (33)	providers.eq7_kernel_function, providers.eq7_cbeta
Default beta=.75	NWQLib, the midpoint of ACL's practical range	Sec. 2.2, p. 14, range [0.7, 0.8]	providers._resolve_eq7
Cauchy kernel	An, Liu and Lin, PRL 131, 150603 (2023), doi:10.1103/PhysRevLett.131.150603, restated by ACL	ACL Eqs. (5) and (181)	providers._cauchy_provider_coefficient
Equal tail and quadrature budgets, epsilon/2 each	NWQLib default	none	providers._build_composite_gauss
Finite-budget cutoff K	ACL bound inverted by NWQLib	Eq. (186), with E1(x)<=exp(-x)/x	providers.eq7_truncation_range
Tail bound	ACL, closed form by NWQLib (Proposition 26)	Eqs. (185)-(186), and Lemma 10 Eq. (62) for K >= 1	providers.eq7_tail_bound
Composite Gauss nodes and weights	ACL	Eq. (61), panel map of Appendix A.4	providers.composite_gauss_grid
Gauss panel selection	Trefethen, ATAP ISBN 978-1-61197-239-9, Chapter 19, applied by NWQLib (strip bound and panel sum, Proposition 27)	Theorem 19.3, Eq. (19.8) for the (n+1)-point Gauss rule, n=Q-1	providers._ellipse_rule
Allowed kernel-quadrature pairs and their bound scope	NWQLib	none	providers.LCHS_KERNEL_QUADRATURE_COMPATIBILITY
Low–Somma kernel j=2, y=1	Low–Somma	Eq. (6)	providers.low_somma_fhat_2
Low–Somma default c=1	Low–Somma	Theorem 2 (any c > 0), Figure 1 as discussed on p. 6	providers._resolve_low_somma_f2
gamma, R, h, J and the epsilon/3 split	Low–Somma	Theorems 2-4, pp. 5-6	providers._resolve_low_somma_pair_profile
Symmetric uniform trapezoid grid	Low–Somma	Theorem 3, Eq. (11)	providers._build_symmetric_uniform_trapezoid
Signed-binary grid	NWQLib	none	providers._build_signed_binary
LCU block sum_j c_j U_j / alpha	ACL, with the complex phases moved from the PREP pair into SELECT by NWQLib	Appendix A.3, Lemma 24, Eq. (178)	quantum.plan_quantum, native.construct_select
Duhamel integral and its quadrature	ACL	Eq. (2), Eq. (72) with one panel	source_selection.select_source_branches, inhomogeneous_theory.duhamel_quadrature
Duhamel remainder	DLMF §3.5 and ACL	DLMF Eqs. 3.5.19 and 3.5.21, ACL Eq. (162)	inhomogeneous_theory._duhamel_quadrature_error_bound
Product-formula step bounds and budgeted step counts	CSTWZ	Propositions 9-10, Eqs. (120)-(121), and the smallest-r rule of Sec. V B	time_independent_terms._fixed_trotter_certificate_records, time_independent_terms._budgeted_trotter_node_record
Dense branch evolution from one Hermitian eigensystem per node	ACL, evaluated by NWQLib	Eq. (70), U(T,s,k)=exp(-i(T-s)(kL+H)), with kL+H=V diag(lambda) V†	time_independent_terms._node_eigensystem, time_independent_terms._branch_unitary, time_independent_terms.spectral_lchs_sum
Shared Pauli commutation table across k-nodes	NWQLib, on the CSTWZ triangle expansions	pruned labels have zero magnitude, so each node's pair and triple rows are the shared rows restricted to its kept labels	time_independent_terms._node_bound_coefficients, time_independent_terms._lchs_census_choice
Pruning bound T*m and its upward-rounded composition	NWQLib	operator norm of exp(-iTH)-exp(-iTH') at most T times that of H-H', the triangle inequality, and directed rounding of the reported pruning, remainder and total	time_independent_terms._budgeted_trotter_node_record, error_budget.upper_dropped_mass, error_budget.upper_pruning_error, error_budget.upper_combined_error
Per-branch identity phase exp(-i t c_I)	NWQLib	c_I depends on k through k*tr(L)/2**q	time_independent_terms._combined_trotter_terms
Shared product-formula schedule and repetition blocks	NWQLib, on Qiskit's Lie-1 and Suzuki-2 templates	angle m*t_j*c_P(k_j)/r_j of address slot j	time_independent_terms._build_product_formula_select_plan, time_independent_terms._repetition_blocks
Structured SELECT on the signed-binary grid	NWQLib	k_j affine in the address bits, so each angle table is affine	time_independent_terms._attach_affine_pauli_structure
Classical Pauli actions	NWQLib	exp(-iaP) = cos(a) I - i sin(a) P	host_pf.apply_node
Periodic Strang bound	CSTWZ, specialized by NWQLib to matching-operator commutators	Proposition 10, Eq. (121), with the factors in circuit order (potential, even matching, odd matching), Proposition 36	periodic.select_periodic_parameters
Periodic Strang step count from trotter_synthesis_tolerance	CSTWZ smallest-r rule, applied by NWQLib to the weighted periodic bound	Sec. V B, r = ceil(sqrt(ceil(B/eps))) in exact rational arithmetic	periodic.select_periodic_parameters, subroutines/trotterization/error_budget._smallest_step_count
Joint-generator SELECT for QSP	Pocrnic et al., arXiv:2506.20760v2, adapted to separate L and H encodings	Section IV, Eqs. (61)-(65)	parameters.select_parameters, compiled_selection._compiled_select_qsp_plan
Source elapsed times in the QSP diagonals	NWQLib	D_L = (t_j/T) k_j, D_H = t_j/T	parameters.select_parameters
CX formula of the compiled QSP SELECT	NWQLib, from the child gate counts as Qiskit 2.6.0 unrolls and controls them, a dense child through its exact synthesis	Qiskit 2.6.0 qiskit/circuit/_add_control.py::apply_basic_controlled_gate, derivation in the docstring	compiled_selection.compiled_select_cx_projection, compiled_selection._qsp_child_gate_census
Exact synthesis of a dense branch before its address controls	Two-qubit KAK and block-ZXZ Shannon decompositions, with the sources in the exact dense synthesis source map	See that map	subroutines/qiskit_compat.controlled, subroutines/_dense_synthesis.dense_unitary_circuit
CX formula of the controlled dense branches	NWQLib, from that synthesis's gate count and Qiskit 2.6.0's add_control	_dense_synthesis.py recursion and Qiskit 2.6.0 qiskit/circuit/_add_control.py::apply_basic_controlled_gate, derivation in the docstring	native.dense_branch_select_cx
Constant-source closed-form reference	ACL, integrated by NWQLib for constant A	Eq. (2): exp(-AT)u0 + (integral_0^T exp(-As) ds) b, which is A^{-1}(I - exp(-AT))b for invertible A, evaluated as expm([[-AT, bT], [0, 0]])	references.closed_form_reference
Matrix exponential of the expm and closed_form references	Al-Mohy and Higham, SIAM J. Matrix Anal. Appl. 31, 970 (2009), doi:10.1137/09074721X	Algorithm 5.1, and Code Fragment 2.1 for triangular input	references._exponential
Propagation of the component errors to the requested output	NWQLib	Norm inequalities and sharp projector trace-norm bridge in Proposition 30	solution_error_budget.propagate_physical_error