Limitations and open work¶

The priorities below address executable scientific workflows, scalable constructions and accuracy at finite measurement budgets. Current capabilities and limitations describe the implementation available now. The guides for LCHS, QPE, QLS, GCiM, QFD and QHD give the details of each method.

Open work¶

Priority	Work	When to proceed
1	One hardware scientific workflow, with resource-model calibration	Begin with one existing method, one backend and an independently checkable instance.
2	Structured inputs and scalable LCHS SELECT	Pursue as a focused algorithm and resource project for one useful problem family.
3	Near-term measurement efficiency and statistical decisions	Develop together with the addition of new near-term methods.

Hardware-backed scientific validation¶

Start with one existing method, one backend and a small scientific instance whose requested output has an independent analytical or classical reference. Define the comparison and applicable statistical acceptance criteria before collecting data. Complete preparation, submission, collection and interpretation of the result, recording the backend configuration, calibration provenance, circuits, shots and execution costs. Use the observed error and cost to choose the next instance and scale.

Qualification applies to the method, backend and workload that were checked. Larger instances need a recorded input source, a reproducible runtime budget and an independent reference appropriate to the requested output.

Resource-model calibration¶

For the selected hardware study, compare the planned resource estimate with transpiled circuits and observed execution costs. Fix the gate set, coupling map, transpiler configuration and versioned hardware assumptions, and include preparation and measurement costs. Keep measured backend costs distinct from portable algorithmic resource formulas. Address device-fit and run-time forecast gaps when they limit that comparison.

Scalable LCHS SELECT¶

The default quantum route limits LCHS to small systems. That route, dense_exact, computes each branch exponential classically and synthesizes it as a controlled dense unitary, so its cost grows with the number of branches and exponentially with the number of system qubits. The qsp_block_encoding backend already implements the whole SELECT at once. It evolves, with QSP Hamiltonian simulation, the joint generator sum_j |j><j| (x) (k_j L + H) of Pocrnic et al., arXiv:2506.20760v2, Sec. IV, which NWQLib builds from separate L and H block encodings instead of one encoding of A. It has a query-cost formula and an error allowance within the LCHS budget. It must be selected explicitly, its polynomial degree is limited by max_qsp_degree (256 by default) while the degree it needs grows with the kernel cutoff times the evolution time, and a dense input can give it block-encoding children that need dense synthesis.

Begin with one scientifically useful problem family whose L and H have compact access. Derive its normalization, required polynomial degree, ancillary workspace and circuit costs, then determine whether dense child construction or work at that degree is the limiting factor. Extend qsp_block_encoding or replace its polynomial transformation on that evidence. Generalized QSP (Motlagh and Wiebe, arXiv:2308.01501v2), which Classiq's LCHS notebook uses, is one candidate.

The work needs access to L and H that does not start from a dense matrix, a way to reach the degrees that long evolutions need, and a default route that depends on the problem size. A replacement transformation also needs its own query-cost formula and error composition with the existing kernel and quadrature budgets. Small cases must agree with dense_exact and with independent matrix-exponential references.

Alongside new near-term methods¶

Develop the following measurement, uncertainty and adaptive-decision work when adding new near-term methods. Select a requested output and an execution budget, then compare implementations by the error and total cost of that output.

For GCiM ground-state correlation recovery, the GCiM options page separates alternative acquisition circuits from changes to the basis, selection rule or Hamiltonian. It records current availability, additional input requirements and the evidence needed to compare results.

Sampling and measurement implementations for near-term methods¶

Lanczos already offers pilot-based sensitivity allocation, and SPE uses importance sampling over Fourier frequencies. The sampled Pauli expectation workflow currently uses uniform per-setting shot counts for its observable measurements. Sampled GCiM and QFD use them by default, and PencilSampling instead chooses their main-stage shots from a pilot stage (Pilot and main shot allocation).

Evaluate importance sampling and shot allocation based on variance or output sensitivity for these workflows and the new methods. Include covariance between quantities measured in the same shots. Compare the final output's bias and root-mean-square error at matched total shot budgets, including the pilot, and report circuit and classical costs. For nonlinear outputs such as projected energies, establish whether improved measurement precision actually improves the final estimate.

GCiM and QFD offer grouped joint-state measurements and prefix pair preparation. ADAPT-GCiM also offers vacuum-based pair preparations when its generators satisfy the documented conditions. Compare these routes and ADAPT-GCiM's explicit termwise Hadamard tests on concrete workloads, and assess alternatives for other near-term methods. Direct measurements, other overlap or interference circuits, and simpler controlled operations are possible choices when they compute the same requested quantity. Compare ancilla requirements, compiled circuit depth and two-qubit gate counts, measurement settings and shots, and sensitivity to hardware noise. State each alternative's input assumptions and preserve any required relative phase and real and imaginary components. The choice of implementation needs an error-versus-cost comparison for the workload.

GCiM physical-output error bounds and automatic assessment¶

Start from an energy comparison, a measurement-budget decision or an adaptive stopping decision. Derive the uncertainty statements needed for that decision before expanding the total-error assessment.

FixedGCIM currently computes the projected energy, H/S matrices, overlap spectrum, selected rank, conditioning and numerical residual/backward-error diagnostics. Its default eigen_error_model declares six unresolved error sources, native_preparation and native_simulation (preparing and simulating the backend circuit), sampling, projected_solve, ground_identification and physical_model. Neither the quantum nor classical FixedGCIM analysis path automatically supplies the corresponding physical-output error facts. A default total result.assess(absolute_tolerance=...) therefore remains INCONCLUSIVE, including for a full trial basis. This describes an unfinished automatic bound-generation capability, rather than a failure to compute the projected energy.

The existing facts= argument accepts supplied bounds whose output, error frame (the quantity, metric, unit and scope that a bound refers to) and Result context match those of the Result. It does not derive those bounds. Absolute-error assessment does not require the exact ground energy. Supplying reference= alone does not resolve missing error contributions. A known reference is useful for an independent development comparison, but it must not become a prerequisite for solving unknown scientific problems.

For that decision, connect method-specific numerical and statistical analysis to the existing error model:

Derive bounds for the measurement of the chosen matrix elements and for the numerical reduction, including their units, assumptions and the observations that contribute to them. Distinguish a measured residual or variance estimate from a bound on the reported energy.
Propagate supported H/S perturbations and overlap-cutoff effects to the selected projected energy, accounting for conditioning. Handle finite-shot covariance and confidence explicitly. Treat trial-space approximation and ground identification as separate questions whose bounds require additional problem-specific assumptions.
Add applicable bounds to the existing error evidence of the Result (facts) and its saved data. Validate analytically checkable full-basis and partial-basis cases, an excited eigenstate with zero residual, and an ill-conditioned Gram matrix. A generic status change or an unconditional zero error term would not complete this work.

The bounds would be computed in algorithms/_eigen_support.py and the FixedGCIM planning and analysis code, and read by evidence/error_model.py and the projected checks recorded by evidence/verification.py. ADAPT and Lanczos share the error-model constructor, so a change to it needs a separate check of the code that computes their bounds and of its assumptions. GCiM diagnostics remain usable while automatic physical-output bounds are developed.

ADAPT-GCiM stopping beyond the flat counter¶

Replace the consecutive-small-energy-change heuristic with a stopping rule supported by evidence about remaining improvement. A flat energy sequence can reflect an unproductive sequence of selected generators or insufficient measurement precision while the physical energy error remains large. The current min(T_auto, T_usr) rule therefore provides an operational stop, not an accuracy certificate.

Evaluate possible rules at matched measurement and classical-work budgets, including temporary plateaus, near-dependent basis vectors and finite-shot uncertainty. Distinguish a small projected eigensolver residual from convergence in the full Hilbert space. The design target is 20–25 qubits. Prefer matrix elements already measured and structural evidence, account for work repeated at every iteration, and require explicit opt-in for additional full-space residual calculations, dense expansion or simulation. The stopping criterion is open research work and is coupled to the statistical questions below.

ADAPT-GCiM finite-shot adaptive statistics¶

The current quantum counts path returns point estimates and compares noisy gradients against a fixed selection floor. It does not quantify selection error across the pool. An exact-null or near-tie screen can choose sampling noise. Earlier experiments over random seeds give no coverage guarantee for the current code.

A future result must record the complete gradient vector, per-setting counts or probabilities, independent seed streams, variance estimates, confidence intervals, the multiple-comparison method, and the separation between the winner and its closest competitor. It must report an unresolved selection outcome when the evidence does not separate candidates.

Finite-shot residual certification remains unresolved. Current sampled-pencil diagnostics keep the record of their sources and do not become deterministic spectral or convergence guarantees. A future residual certificate needs a joint statistical model for every contributing quantity. Clipping an inconsistent estimate cannot supply that model.

Within one ADAPT run, gradient and matrix-element circuits consume a shared increasing sequence of simulator seeds. Distinct seeds do not establish statistical independence. A certified study still needs justified independent substreams or a covariance model that accounts for shared randomness.

Quantum exact and counts execution paths remain available. Finite-shot adaptive selection and finite-shot residual certification remain statistically unresolved and are not certified convergence evidence.

Application-driven extensions¶

Start these extensions when an identified scientific problem needs them. Establish the input, requested output, strongest relevant classical comparison and resource budget before expanding the supported model.

Time-dependent LCHS¶

Start from a concrete dynamics application that needs a time-varying generator or source, and specify the requested observable, input access and execution budget.

The quadrature, evolution ordering, input records, and error composition must be designed together. Reusing the time-independent records without those semantics would create a misleading API.

Scalable input and execution models¶

Users evaluating asymptotic advantage need problem-loading and Hamiltonian-access models that do not begin with a dense matrix.

Each algorithm needs an explicit oracle specification, cost model, and small independently checkable construction. A generic “sparse” flag would hide algorithm-specific assumptions and is not sufficient.

The augmented-Lagrangian layer and box refinement run one ordinary QHD Plan per round or level and search the finite grid of each box. They follow Wu et al. arXiv:2605.12066v1, Sec. III.B and Sec. V. Choose a concrete optimization task and a competitive classical baseline before expanding these methods. The items below are candidates when that comparison identifies the corresponding limitation. An item that would change a default needs an independent classical oracle or a derived bound first.

Narrower eligibility for a quantum automatic default. Qualify a restricted class beginning with one affine inequality, KineticGroundState, joint-mode readout, no refinement and exact noiseless readout. Executed evidence currently covers only the three n = 3, K = 4 grid configurations under the settings that the guide's measured-evidence row for the comparison states, while the cases up to fourteen variables show that planning accepts them and give their resource counts. Qualification needs converted quantum trajectories that preserve feasibility and pass a predeclared best-feasible objective margin against the default Powell–Hestenes–Rockafellar (PHR) inequality form at matched budgets, checked against independent feasible-objective references. The resource evidence must also include backend-circuit costs, representation-selection overhead, observed simulator time and total process memory. Quantum means, quantum refinement, finite shots, other initial states and other schedules remain untested by this study, and the classical sensitivity checks do not establish quantum time-step convergence.
Scales. The augmented-Lagrangian layer takes the objective and constraint scales from the caller, with default 1. Estimating them from the problem, for example with the gradient-based factors that Algencan computes once at the starting point (Birgin and Martínez, doi:10.1137/1.9781611973365, Eq. (10.5)), or rescaling between rounds, would reduce how much the run's decisions depend on the units in which a constraint is written.
Finite grid and continuous problem. Stopping statuses describe computed feasibility and complementarity at the selected points, and constrained_grid_minimum compares evaluated finite-grid values. A verifiable error relation between the finite-grid result and the continuous constrained optimum would let a grid result support a statement about the continuous problem. A narrower open question concerns the grid problem itself. With best_observed under exact readout, every grid point supported and no box refinement, each round's point is the grid point with the least binary64 table value of L_k, the smallest index on ties. Whether the convergence results of Birgin and Martínez, doi:10.1137/1.9781611973365, Chapter 5 (Assumption 5.1, Theorems 5.1 and 5.2), then apply to the finite-grid problem, with the grid taken as the book's lower-level set Omega, is not established, and max_penalty and the default absence of a safeguard depart from Algorithm 4.1 as those results assume it. The layer claims none of these results.
Retry an unfinishable round or level. Resuming a saved run cannot give an unfinishable inner Run a fresh attempt. Add an explicit option that creates a new Run while keeping the old attempt's counted or uncertain work and counting the new attempt against the cumulative limits, so retrying cannot erase its cost.
Finite-shot decisions. With shots, every round's point and every level's intervals are read from samples. Apart from the Hoeffding screen that accepts a stall-split valley, the probability that sampling changed one of these decisions is not quantified, and the shots per round are not chosen adaptively.
Feasibility restoration and polishing. A grid point rarely satisfies an equality exactly. A local feasibility restoration or classical polishing step would be explicit post-processing with its own budget, reported apart from the QHD result so that it does not enter the performance of the quantum part.
Shared evidence and reports. An augmented-Lagrangian result keeps its own record and inner Results and has its own report and archive. It does not supply error evidence (facts) to Certificate or ErrorModel, and the CLI report command does not read its archive.
Early-trigger refinement split. stall_split="best_region" acts only when the whole box stops changing, so for the same observations the best relative objective is never worse than that of the run without a split. Splitting as soon as one axis is kept whole and its marginal has two separated peaks did better on the double well of the QHD guide. With the search model at potential gain 1, the uniform initial state, max_levels=10 and max_no_improve=10, it split at the first level and reached best-point error 0.000857 in ten solves after discarding 5.7% of that level's distribution, where the stall split reached 0.0108 in ten levels. At the default gain 8, with the other settings unchanged, the double well does not stall within ten levels, so the stall split does not act, and the run reaches 0.00325. The early trigger changes the run before any stall and so loses that guarantee. It needs its own risk statement and an explicit decision before it becomes an option.
Best-point Gaussian width and regime. level_initial_state="best_point_gaussian" has been compared only at width 1/6 of the level box, on one problem where it helped and on three two-variable problems where the kinetic ground state at every level did better (box refinement). A comparison over widths and problem classes that shows when the Gaussian helps would support a rule for choosing it, or a width set from the level's distribution.
Application instances. An instance generator for the AC optimal power flow problems of Wu et al. would need an independent implementation and validation.

QHD core models¶

Prioritize the model or numerical limitation that affects a selected application or accuracy claim. The items below describe unresolved cases, with their own mathematical and resource requirements.

Walsh reconstruction in IR comparisons. Propagate R_W into the readout accuracy of the ir_product reference, the mode tie window and the verification uncertainty. When the resulting uncertainty interval crosses the requested tolerance, report an inconclusive comparison so the decision does not exceed the available numerical evidence.
Constant-phase report. Report the objective constant's phase contribution separately from the rest of the phase-sensitive error record. This would help readers distinguish a global-phase contribution from errors that can change readout probabilities.
Classical time integration. The integrated coefficient rule removes the quadrature error of the schedule coefficients but keeps the second Magnus term, and on the measured model, with s = 1 for the cubic schedules, both coefficient rules converge at second order with errors within 1% of each other. Whether a higher-order Magnus or interaction-picture classical model is needed depends on accuracy and cost evidence on a stated model.
Kinetic terms below the normal binary64 range. Planning refuses a kinetic coefficient, energy, hopping term or one-hot kinetic phase whose exact value would fall below 2**-1022, and names the grid spacing and schedule weight as the remedy (algorithms/qhd/validation.py::_kinetic_range_error). With schedule weights of order one, only grids with spacing near 1e150 or wider reach this, for example K = 64 on a periodic binary box of width 6.4e154. The omission rule for objective data would extend to them by replacing an analytic Walsh coefficient or dense energy with zero and bounding its error by the radius of its enclosure about zero, adding each omitted one-hot hopping occurrence's exact intended exponent magnitude to the omitted-block error bound, and adding each omitted one-hot kinetic diagonal event's exact identity-phase magnitude to the identity-phase error bound. This needs those radii in the binary angle-formation bound, the counts in QHDRangeOmissions and a public test per case.
Binary encoding on Dirichlet grids. A sine-transform circuit for the Dirichlet kinetic term would let the binary encoding use the Dirichlet grids, which the current QHD limitations exclude.
Odd K on the one-hot periodic grid. The second-order one-hot product splits each variable's links into two layers of disjoint links (kinetic.KineticCompiler). On a cycle of odd length K the wrap link (0, K-1) falls in the even layer and shares point 0 with link (0, 1), so the one-hot periodic grid requires an even K (method.QHD._periodic_grid). The rule also covers first order and the classical variants, so that the periodic grids a Method accepts do not depend on its product order or its execution. A compiler that splits odd cycles into three layers would accept odd K. It needs the third layer in the compiled step, its terms in the splitting bound of evolution_bounds, its CX and rotation counts, and a comparison of the emitted blocks with the circulant stencil at odd K.
Start-vector bound on very large grids. The construction error bound of a product initial state includes (d - 1) ceil(sqrt(K**d)) 2**-1074 for gradual underflow in forming its tensor product (initial_state.restricted_state_error). This term grows with sqrt(K**d). Planning therefore refuses a Gaussian state or the Dirichlet kinetic ground state when the bound in units of u exceeds the binary64 range, including from d = 2035 variables at K = 4, and recommends UniformState(), whose bound does not depend on the grid size. This shared planning check also stops resource-only and circuit Plans, although circuit preparation acts on each register separately without forming the full tensor product. The Gaussian state_preparation error estimate uses the aggregate construction bound, while the Dirichlet kinetic-ground estimate uses per-register bounds. Supporting these Plans at larger grid sizes needs a planning check appropriate to circuit preparation and a derived Gaussian preparation estimate independent of forming the full tensor product. Classical evaluations that form that product still need its construction bound.
Split-step budget beyond first order. The Plan and observed state budgets of the split-step variant are first order in u (split_step.state_error, split_step.evolve). A budget valid at every order in u would remove the first-order approximation from its mass and tie windows. This requires the higher-order terms of every error contribution in both the Plan and observed budgets, together with an upper bound computed during planning on the tie window obtained from the new observed budget.
Dense phase-diagonal circuit error. The circuit error record (algorithms/qhd/resources.py::circuit_resources) has no bound for converting a dense phase diagonal into gates (subroutines/_multiplexors.py::append_control_diagonal_phases). This conversion wraps each phase with NumPy's complex exponential and argument and forms recursive means and Gray-code rotation angles. Its diagonal_wrap entry is therefore unavailable for a one-hot circuit with a number projector on 4 to 6 qubits and for a binary circuit with a dense diagonal. A bound on this contribution needs qualified error bounds for NumPy's complex exp and angle and a rounding bound for the multiplexor recursion. A total circuit-error bound would also need the other error sources.
Periodic Gaussian initial state. GaussianState uses the chart distance |x - c| of the coordinates in [lower, upper) on the periodic grid, so for a center near one end of the period the points across the wrap link get the amplitudes of their chart distance. A variant that moves with the center across the period would use the minimum image min_k |x - c + k L| for the period L, or a smoothly periodized Gaussian sum_k exp(-(x - c + k L)**2/(2 sigma**2)). Its application is the best-point Gaussian start of BoxRefinement(level_initial_state="best_point_gaussian") on a periodic grid, whose best point can lie near either end of the level box. It would be a separate option, so that states supplied with the current definition keep their meaning.
Tie window with an excluded preparation. When Qiskit's state preparation is built from an excluded unitary, the preparation record has no state-error bound and the most probable point is chosen by exact comparison (Readout). A qualified window needs preparation-time evidence that the errors in constructing and applying the supplied matrices fit the target's per-instruction state-error allowance, or a separate bound eta_excluded on the excluded instructions, giving delta = (C G_covered/2 + 5) u + eta_excluded for the covered instruction count G_covered and the target's per-instruction constant C. That delta is relative to the saved backend circuit, so a comparison with the finite model also needs the error of preparing its initial state and of compiling the model into the backend circuit.
Backend readout cases without a check. Two cases of the QHD readout windows have no test or check (Readout). The tie window of exact NWQ-Sim probabilities takes its state budget from the per-instruction constant NWQSIM_CPU_SV_ROUNDOFF_PER_OPERATION, and no QHD run on NWQ-Sim has compared its chosen point and window with an independent high-precision evaluation of the same circuit. The readout roundoff of exact backend probabilities uses COMPONENT_SQUARES_ROUNDOFF, derived for a bin that holds one amplitude, which holds when the backend circuit has no qubits beyond its register. The first case needs such a run with a qualified NWQ-Sim runner. The second needs a check of which QHD constructions add ancilla qubits and, for any that do, a readout window for bins that sum several amplitudes.
State-dependent Schrodinger tie window. The schrodinger variant computes a first-order state-error contribution for each expm_multiply call from bounds on its shifted generator norm and its number of nonzeros per row. A budget that follows the state, as the split-step variant's does, needs the Taylor terms of SciPy's algorithm (Al-Mohy and Higham, doi:10.1137/100788860, Algorithm 3.2). With Z the scaled shifted generator and theta >= ||Z||_2, the term errors satisfy e_j <= (theta/j) e_(j-1) + r_j with a sparse-product defect r_j per term, and, when theta < m + 2, a geometric bound covers the tail after the Taylor degree m that the algorithm uses. SciPy returns neither the terms nor the degrees and scaling counts it chooses, so the bound needs an instrumented implementation with several O(D) reductions per Taylor term, or one more sparse matrix-vector product per term for the componentwise defect.
Verification bytes from observed outcomes. For QHD verification from observed counts or probabilities without a kept state, QHDVerification.max_bytes does not bound the success-mass step at any readout width. In the measured 68-bit case, the memory used exceeded the bytes counted against the limit. Count the bytes of persistent histogram arrays, first-access construction or loading buffers, index expansion, decoding masks and coordinates, temporary arrays created during selection, the arrays of the previous block of observations and accumulated divided passing weights, including storage for invalid outcomes. Read count lengths and classical-layout widths, or probability-array entry counts and widths, from metadata before building or loading a histogram array. Add persistent payloads to the live grid data and combine the successive grid-evaluation and histogram scratch phases by a maximum, without counting shared payloads twice or treating the grid slab size as a bound on histogram entries. Qualify the object constants and reader buffers on the chosen 64-bit CPython and NumPy versions by tracing components on narrow and wide readouts, invalid outcomes and noncached local indices, using both cached histograms and their first construction or load from recorded blocks of observations. These allowances cover the named allocations, not the resident memory (RSS) of the process. The kept-state route keeps its current byte formula. The entry-dependent work count remains open together with the byte formula.

Structured Dalzell coefficient recovery¶

The Dalzell kernel-reflection polynomial has known even structure, but the current constructor samples it on a degree-dependent grid and recovers coefficients with dense Chebyshev least squares. A recurrence or stable orthogonal-transform construction could reduce classical coefficient-generation cost at large degree.

A replacement must preserve the coefficient convention, even parity, the K(0)=1 anchor, and the off-kernel bound -1 <= K(x) <= -1 + 4 eta/(1 + eta) for 1/kappa <= |x| <= 1, with independent large-degree stability comparisons. The existing fit remains the validation-scale fallback until those checks establish equivalence. This item is separate from the shortcut Program's composition of child block-encoding queries.

Maintenance and targeted accuracy work¶

Address these items when they affect an actual workflow, resource limit or requested accuracy claim. Their scope and evidence requirements are specific to the affected construction.

Controlled dense-query construction¶

Exact dense synthesis describes how NWQLib synthesizes and controls dense unitaries and how each construction checks that work against its limits. The following points remain open.

Dense input does not scale. The synthesis takes of order M³ classical work and M² gates for an M-square unitary, so a dense encoding serves small validation systems. A larger system needs a structured encoding with its own oracle, normalization, query count and ancillary workspace (scalable input and execution models).
The gate-wise control formula counts the synthesized gates of a controlled composite. Qiskit also controls the other gates of the composite, such as the projector phases of a QSP pass or the input preparation of a controlled LCHS source branch, and a backend then compiles the multi-controlled gates of the controlled circuit. These costs have no formula.
The parity control of a QSP pass controls the whole pass and is gate-wise on every route. Controlling the pass instruction by instruction would let the whole-matrix route reach a dense child at every query, and it would need its own CX and work formulas.
No limit of its own bounds the memory of the Run's synthesis cache. The Run sets aside each synthesis's work from max_synthesis_work before caching its circuit. At the default limit of 10**9 work units, a count of scalar operations derived from the matrix sizes rather than measured CPU time, the byte allowances of the cached circuits plus the matrices that serve as cache keys total at most about 2.1 GB, or about 0.79 GB when every cached circuit has at least three qubits and about 0.41 GB when every one has at least four.

End-to-end certified resource estimation¶

Choose one method, requested physical output, total error target and set of assumptions before developing a certificate. Current estimates identify their basis, precision, the circuits and shots they count, and their model scope. Several families bound some error sources, but no family composes every source of its output into one certified budget.

The selected case needs a complete error budget, applicable rotation and T-count formulas, and interval arithmetic so that rounding cannot move the certified bound. Extend to another family when a concrete application needs it. GCiM physical-output bounds belong to the work alongside new near-term methods.

Planning and resource estimation have separate work limits. QPE checks the work of computing commutator bounds for its error-budgeted product formula against max_work. With max_work=100000000 and 20 system qubits, the initial work bounds for this calculation alone permit at most 531 terms for the full expression and 6316 for the relaxed expression. An exact-trajectory Plan must also fit the work set aside for choosing its common time step, which depends on its power schedule and inputs, so these counts alone are not term limits for an exact-trajectory Plan. Resource estimation limits the work needed to total a Program's resources to max(max_steps, 24*a), where a is the work of checking the Program during planning. The current accounting is described in Resource estimate bookkeeping.

QHD objective expansion check¶

Before expanding a QHD objective, planning computes a nominal expansion-term count (objective.expansion_term_charge) and compares it with work and 16-byte-per-term allowances. The count bounds formal terms for numeric constants, symbols, Add, Mul and nonnegative integer Pow. Atoms contribute one, sums add, products multiply and an m-th power of an n-term base has at most binomial(n+m-1,m) products. It is not a bound on general SymPy rewriting time, temporary storage or coefficient size. After expansion, objective.node_count counts the work of table evaluation. A complete pre-expansion resource model remains open for the following cases.

It does not count coefficient digits, so the default max_work accepts the expansion of (x+y)**1000000, whose binomial coefficients reach about 301,000 digits.
A product of negative powers forms the product of its expanded denominators, which the nominal count does not multiply. 1/((x+y)**400*(z+w)**400) has a largest partial count of 401 but receives count one, and its expansion forms 160,801 denominator terms.
Function rewrites can create powers or sums outside the counted polynomial skeleton. With positive variables, log(x*y*z)**3 has count one and expands to ten terms, while exp((w+4)*log(x+y+z)) has count one and expands to fifteen. Disabling the log, power or basic expansion hints alone does not prevent the growth that SymPy's expression constructors can cause in general.
An unevaluated Sum or Product in the objective is lambdified as a loop over its range but counted as a few nodes, in planning and in the grid_minimum verification.

Preparing every setting on remote backends¶

prepare(plan, settings="all") prepares every setting that is known before running only on local Aer, local NWQ-Sim and classical execution (Run on a backend). Extending it to IBM Runtime, IonQ, Nexus, Slurm and other remote backends needs two decisions. The first concerns cost. A remote preparation can be a compilation that the provider performs and may charge for, so preparing every setting in advance commits that cost for settings that a later cancellation or failure would never submit, and the caller needs to see that cost before the request. The second concerns pending preparations. A remote compilation can remain pending as a PendingPreparation, and the Run never repeats a preparation attempt on its own. Preparing every setting therefore needs defined behavior while some settings are still pending, namely what Prepared.circuits and Prepared.setting_names expose and how a later submit or a reopened Run advances the pending compilations without a second attempt and without submitting a setting from prepare.

Device fit and run-time forecasts¶

Check device fit and run time answers whether a Plan fits a declared device and how long it may take only in part. Three gaps remain:

The built-in AER_STATEVECTOR_TARGET declares its readouts but not its artifacts, instructions or program_nodes (Program steps), so its capability check is unknown until the user declares them with AER_STATEVECTOR_TARGET.revise(...) (Assess a plan against a device).
A finite time model needs an exact logical operation count. A Plan whose operation count is unavailable, such as the QLS example of Estimate resources, gets an unavailable time prediction (Finite time models).
A residual between a forecast and a measured Aer time needs the forecast's runtime and compiler to equal those in the preparation record. The Aer adapter builds these sources internally and NWQLib exports no constant for them, so a user copies them from run.data.receipts[0] (Compare forecasts with measured times).

Closing these gaps needs built-in targets that declare what Aer supports, logical operation counts for the Plans that lack them, and public runtime and compiler sources for Aer.

Possible extensions¶

These extensions require a concrete use case, input data, and an independent validation method before work begins:

Additional QPE estimators or broad estimator benchmark matrices.
Additional low-level synthesis backends.
New chemistry excitation pools beyond the fixed and adaptive pools that NWQLib provides.
Provider-specific convenience layers beyond the common checked execution path.
Trial dependency upgrades across release series that compare every saved scientific reference value and ignore the version fields that only record provenance.
Visualization and report front ends beyond the structured records and JSON reports that workflows already produce.

Before implementing an extension, define its scientific input, public interface, independent validation method, resource limits, and effects on existing checks. Listing an extension here does not promise future compatibility.

Current capabilities¶

Every method follows one workflow. A Problem and a Method give a Plan, from which NWQLib estimates resources, runs the circuits as a Run and attaches the Result, and saved data is loaded only on request.

Area	What it does	What has been checked
Finite Pauli expectation	Exact readout by group, in which the simulator computes the expectation values without shot sampling (`shots=None`), subject to numerical rounding, or sampled readout without error mitigation that measures the whole register once per qubit-wise commuting group. Binary readout calibration uses parity shots. Point estimates, variances and fixed-time, time-uniform and Beta intervals stay separate quantities.	Small cases against independently computed scalars and injected test data. Shared calibration covariance stays explicit. An interval uses only a chronological prefix of complete sets of fresh shots, so data that omit an earlier set or keep a chosen subset leave it unavailable (Accumulated and imported data). The statistical assumptions and the transfer through the readout channel do not certify complete physical accuracy.
Chebyshev Lanczos	Configurable initial state and Krylov dimension, a shared signed-walk Plan, exact marginals or counts, resource totals of the planned construction, an explicit two-stage sensitivity analysis and a classical recurrence.	Small cases against independent algebra, and checks of the Chebyshev chain that was executed. Projected values keep the original target, the smallest eigenvalue, and leave the `ground_identification`, `native_preparation` and `native_simulation` error sources unknown. The finite basis, the empirical shot allocation and the cutoff do not certify physical energy accuracy.
LCHS	Time-independent homogeneous and constant-source inhomogeneous workflows. Classical and quantum Methods with explicit verification, finite scalar refinement, coefficient MPS analysis and resource operations. Dense, product-formula and compiled-QSP evolution.	Small dense references, independent product-formula references, agreement of circuits with theory, and the saved scientific examples.
QPE	QCELS, SPE, RFE and adaptive RWPE share one planning and execution workflow, and coherent QPE is separate. Controlled powers of dense unitary or Hermitian matrices and of real `PauliSum` Hamiltonians.	Tiny cases against analytic values, injected data and backend circuits, tracing of each reported scalar to the observations it came from, explicit checks of the nominal eigenvalue cluster, and reopening of saved RWPE runs. Local intervals do not establish a complete accuracy claim.
GCiM	FixedGCIM and explicit ADAPT on shared planning, execution and report code, with several operator pools, optional reoptimization and verification. Chemistry inputs, which need the optional chemistry packages, and conversion of a supplied real double-factorized (DF) Hamiltonian feed the existing Jordan–Wigner energy evaluation.	Small cases against independent matrices, gradients and backend circuits. Sampled and exact pencils stay distinct, and saved runs keep the current basis apart from its history. Chemistry and larger trajectories keep their separately recorded scope.
QFD	Quantum filter diagonalization in a basis generated by integer powers of one fixed product-formula step applied to one reference state. Pauli-term and fermionic-block evolution, exact readout or grouped shot measurements, optional pilot-based shot allocation, and saved-result reanalysis at another overlap cutoff share the public workflow.	Small cases against independent matrices, circuit phase and particle-number checks, save/load and interrupted-run checks, and compact 26-qubit planning. The QFD guide explains the projected energy and the conditions of these checks.
QLS	Three solvers, the polynomial inverse and two shortcuts, with support for the structured periodic family, shared classical and quantum execution, and explicit resource and verification operations.	Small cases against independent numerical results and backend circuits, with each physical output tied to its source. Classical computations used by the method, conditioning on measurement outcomes and unsupported capabilities are stated explicitly.
QHD	One-hot construction on a Dirichlet or periodic grid, or binary construction on a periodic grid, with compact Qiskit circuit blocks and bounded classical evolution. Quadratic, cubic and shifted-cubic schedule records with midpoint or integrated coefficients, uniform, kinetic-ground and Gaussian initial states, classical Schrodinger, circuit-product (`ir_product`) and split-step variants (`QHD(theory_flavor=...)`), readout of the candidate (the observed valid grid point with the least objective value) and of the most probable point with conditional position statistics, and explicit verification. Rotation counting formulas, budgeted T estimates and separate per-circuit error sources for one circuit, with rotation and T totals over an augmented-Lagrangian run or a refinement. The augmented-Lagrangian layer solves a `ConstrainedOptimization` by a sequence of QHD Plans and has an explicit evaluated finite-grid reference. Box refinement repeats QHD on shrinking boxes in a physical or dimensionless search model, with an optional stall split. Both layers can save their state in a directory. Resuming reopens the recorded inner Run with its recorded Plan and run log. A recoverable Run can continue, including a remote result that becomes available after an interruption. Some interrupted local preparation, data collection or classical evolution cannot be completed from the saved Run.	Small cases against independent grid, action and outcome-distribution values, grid minimization as the reference for the augmented-Lagrangian trajectories, and independent interval, box and invariance checks for refinement. With the default `trotter_order=2` and `rotation_threshold=0`, which drops no rotation, the one-hot circuit's probabilities approach those of the finite model at second order in the step length, in a small periodic case with a constant Hamiltonian (Scientific model and ordering, Proposition 40). This convergence is distinct from discretization, floating-point error and global optimality. Stopping statuses and refined points describe finite-grid searches, not continuous or global optima.
Shared subroutines	State preparation, LCU, block encodings, Hamiltonian evolution, QSP and QSVT, QPE primitives, resource records, reporting, validation, and the interfaces to local and cloud backends.	Unitary and statevector tests, independent numerical references, mutation tests of scientific formulas, and mocked providers in CI.
Saved data and public workflows	Typed verification records and certificates, size-limited saved metadata and loading of the binary arrays a Plan uses, explicit reanalysis, method cards and configuration schemas, the CLI, and conformance checks with explicit scientific cases for Method authors.	Loading a saved `Source` imports and runs nothing, a missing array stays unavailable, and reanalysis runs no new measurement. Checks of an installed package establish only the package and runtime combination they ran on and the behavior of records in the current format.

The top-level nwqlib package exposes the common workflow and search functions, the Problem and output records and its modules, and imports each module on first use. Methods and their Result classes are imported from their own packages, and the API reference documents the exports.

Limitations¶

Chebyshev Lanczos¶

The Lanczos guide documents the one method, chebyshev_lanczos. Its signed Pauli encoding does not require THC factorization. It is not a THC-BLISS implementation or a claim of minimum hardware cost.
Exact and counts-based quantum execution and classical execution each have their own numerical and circuit limits. ExecutionLimits bounds the simulator width and the size of the recorded data, and raising those limits does not make execution scalable.
Pilot allocation, overlap regularization, and projected numerical diagnostics are empirical choices. Final main-stage moments exclude pilot data. Nonlinear Ritz energies have no finite-budget unbiasedness or physical-accuracy certificate.
Fermionic Hamiltonians read in XACC format with read_xacc or parse_xacc keep the literal order of their ladder products, with at most four ladder operators per term. Reference occupations, electron sector, spin ordering, and units are explicit metadata, never inferred chemistry conventions.
Counts readout adds a label-controlled basis transform for odd-degree moment settings. Fewer measurement settings do not establish lower compiled gate cost or improved hardware accuracy. Broader encodings, certified adaptive statistics, and Hamiltonian-specific factorization improvements require a concrete workload and independent evidence.

LCHS¶

Time-dependent LCHS is not implemented. The public workflows and shared term construction implement time-independent homogeneous systems and time-independent systems with a constant source term.
Dense input and dense-exact SELECT provide explicit finite constructions, not scalable oracle models. Product-formula and compiled-QSP (qsp_block_encoding) choices expose more structured circuit paths, but end-to-end scale is still limited by simulation and synthesis. Scalable LCHS SELECT describes the planned work.
The general dense-input compiled-QSP route derives its Pauli and block-encoding decomposition from that input. Structured families chosen separately use their documented access and error formulas. Neither route establishes an arbitrary sparse-input oracle.
A comparison of the LCHS quadrature node count with an asymptotic cutoff recipe needs a stated comparator. Before the guide compares the two, the recipe, its parameters and the evidence must be named.
The contribution of the positive endpoint of the LCHS Low–Somma quadrature is checked only through an assertion that also allows for the surrounding tolerance. A sharper independent check that separates the endpoint's contribution from that tolerance is an open item.

QPE¶

By default a dense target uses an explicit dense power circuit, and a real PauliSum Hamiltonian uses compact product-formula powers directly. The classical route remains a bounded nominal-overlap model. Arbitrary unitary targets do not acquire a scalable power oracle.
Estimator confidence statements are method-specific. They should not be compared as if all methods returned the same statistical object.
Broad performance benchmarking across hardware providers is outside the current validation suite.

GCiM¶

Chemistry input requires the optional chemistry packages and is unavailable when they are not installed.
FixedGCIM reports projected matrices and numerical diagnostics, but nothing computes the required physical-output error bounds automatically for its default total-error assessment. See GCiM physical-output error bounds and automatic assessment for the open work and the guide for quantities available now.
Adaptive pool growth and optimizer behavior are research choices rather than universal convergence guarantees.
Circuit-level validation is limited to small active spaces. Larger chemistry studies require an external resource and accuracy assessment.
Exact quantum ADAPT execution reads Pauli values from the backend, execution with shots measures the planned groups, and classical execution evaluates its queries without running circuits. The three paths share one projected solver and keep sampled values distinct from exact algebraic ones.
The ADAPT iteration state keeps the completed basis and its history apart from the current analysis throughout a Run. Reopened scalar results need no rebuilt pencil for reporting or scalar comparison, while state verification needs the current pencil and the planned data access. A noisy enclosure diagnostic is distinct from a refusal by the solver.

QFD¶

The QFD result is a Ritz value in the trial space generated from the chosen reference. The basis can miss the ground state, and neither the projected solve nor the sampled enclosure diagnostic establishes full-space ground-energy accuracy. No energy-error bound is computed automatically.
The product-formula settings and time schedule define the trial space. QFD does not evaluate the distance of its evolution step from exact time evolution or propagate an evolution-error bound into the reported energy.
Fermionic-block evolution preserves particle number. The Plan states the reference's particle-number sector only for an occupation reference or a product of computational basis states. Conservation of spin projection and total spin, a symmetry of the target Hamiltonian, and molecular accuracy each need separate evidence. The fermionic route distinguishes the evolution generator from the Hamiltonian whose energy is reported.

QLS¶

Oracle-first structured inputs do not support every automatic norm or condition-number strategy. Unsupported combinations fail explicitly.
Analytical resource formulas describe their planned constructions and declared domains. They are not optimized hardware counts.
Polynomial fits and shortcut constructions remain sensitive to the supplied spectral interval and error budget.
The classical QLS model records its implementation as qls.classical_model version "2", while the arguments that a classical QLS Result records in applications[0].arguments have changed within that version. To close this, define whether an implementation version covers these arguments. If it does, increase the version when they next change, otherwise state that they are not part of the version.

QHD¶

The QHD core solves a supplied bounded objective using the one-hot encoding with a Dirichlet or periodic finite-difference kinetic term, or the binary encoding with a periodic finite-difference or spectral kinetic term applied through the QFT. The one-hot periodic grid needs an even number of points of at least 4, a limit of the one-hot link compiler rather than of the periodic Laplacian. The classical split-step variant applies the kinetic term in its eigenbasis and also offers the periodic spectral kinetic term for the one-hot model.
The binary encoding supports the periodic grid only. A Dirichlet binary kinetic term needs a sine-transform circuit. Its diagonal syntheses are chosen table by table, and parity-network synthesis that shares CX between Walsh strings is not implemented. Its structured preparation prepares the uniform state only, and other initial states use Qiskit's state preparation.
Statevector and shot validation is limited to small grids. The classical route evaluates the planned finite model. It is not evidence of scalable circuit execution.
The augmented-Lagrangian layer and box refinement search the finite grid of each box, except that the mode_or_mean point rule can read an off-grid mean position. The default stopping status feasible_complementary does not include the projected stationarity of Birgin and Martínez, doi:10.1137/1.9781611973365, Eq. (10.6), so it is not convergence and does not assess optimality. Finite-grid multipliers are dual iterates of the grid problem, not continuous KKT multipliers. A problem with equality constraints must set feasibility_tolerance, because the residual a grid can reach depends on the grid, and when no grid point satisfies an equality within that tolerance, no penalty value lets a grid point pass the feasibility test.
Constraint and objective scales default to 1, which assumes a dimensionless problem. The layer does not estimate them, and an unscaled constraint factor of 1000 multiplies the quadratic penalty by 10**6. A larger penalty can increase the range of the augmented objective and, on the Schrodinger path, the planning work counted from its generator norm. The range need not grow monotonically because the objective, linear multiplier term and quadratic penalty can cancel. A round after the first that the planning checks refuse ends the run with inner_failed, while a refused first round raises its original error.
BoxRefinement defaults to the search model with potential gain 8, chosen from a two-variable comparison of the physical model with gains 1 and 8 of the search model, so no gain is established as best in general. The step count must resolve both the kinetic and potential evolution. A larger gain can increase potential phases and product-formula error, so the default gain does not establish that a chosen step count is adequate. The gain-8 convergence evidence uses classical Schrodinger evolution. For a circuit-product comparison, also check the product approximation and the numerical-reference limitations of ir_product. The search model solves a different dimensionless Hamiltonian at each level and makes no claim about the dynamics of the original QHD Hamiltonian. stall_split="best_region" discards the probability of one region on the evidence of two objective values, so it is not the default. With a positive split budget it refuses a periodic grid, because one cut opens a periodic axis but does not divide it into two regions.
The Schrodinger and split-step variants need different step counts for the same final-state error. In the measured comparison the split-step variant reached each error target with less work and wall time than the Schrodinger variant, by a factor that depended on the grid, the target and the penalty, so no general rule selects the cheaper variant for a new problem.
Runs of the problems of Wu et al., arXiv:2605.12066v1, and Liu et al., arXiv:2607.16996v1, use NWQLib's numerical choices and are cross-checks of those papers, not reproductions of their computations. Explicit comparison settings with ShiftedCubicSchedule(s=1), the form that Liu et al.'s code runs as qhd-c, reproduce the binary CX and Rz counts of their Table II, whose caption names QHD-C, which their Eq. (92) defines as the cubic schedule (circuit synthesis). The optimization tables of Wu et al. and the dynamics figures of Liu et al. are not reproduced. Wu et al.'s results used a left-endpoint time rule, level rescaling rules and left-endpoint boxes that NWQLib does not offer, because the time rule is first order, the boxes are shifted by half a cell and the rescaling couples the potential gain to an evolution speed that depends on the coordinate units. The arrays and the executed schedule behind Liu et al.'s figures are not available.
The size and work counts that planning checks for an objective miss some expansions and evaluations (QHD objective expansion check).

Cross-cutting¶

Dense simulation is used deliberately as an independent validation layer in several algorithms. A successful dense run is not a scalability claim.
Cloud backends share one checked workflow to prepare circuits, submit jobs, refresh job status and collect observations. Checks with the public SDKs and injected responses cover the documented mappings of counts and estimates. Live credentials, queues, compilation and hardware have not been qualified. See backends, IonQ and Nexus.
Resource estimates identify their basis, precision, the set of circuits and shots they count, and their model scope. Analytical, representative and observed backend counts must not be treated as interchangeable.
The resource total of a Plan describes the workload the Plan represents, while preparation and execution keep their own records of preparations and events. Resource operations do not submit jobs or turn unknown backend costs into zero.
Saved reference outputs and example checks protect recorded behavior and complement independent numerical tests without replacing them.
Independent numerical checks that compare the operator matrices of the circuits built by the block-encoding, LCU and QSP subroutines with the documented wire order and controlled relative phase, and checks of compact planning at 10 and 40 qubits beside the 80 to 100 of the resource notebook, have not been run. Each is a small script that takes minutes per revision. The notebook workloads agree bitwise between releases, and these checks would cover inputs the notebooks do not reach.
load_run reopens a classical Run, and a Run saved on the default AerBackend(), without a backend argument. A Run saved with another configuration needs its backend passed in. Closing this needs a saved configuration that names its backend class for every provider, so that no Run needs its backend restated.
The errors that reject an unsupported saved format of a Run, a Result, a method's archive and a QHD refinement or augmented-Lagrangian directory use different messages. One message form for all of them would give the format label found, the supported label and, as the remedy, the release that wrote the data, and would name the matching resume function when such a QHD directory is passed to the resume function of the other layer.
QasmMaterializationBudget requires max_dynamic_visits and max_operations. Defaulting them to the QasmWriteReceipt fields dynamic_visits and expanded_operations, which the field descriptions name as their source, would let a user set only the byte and width limits.
A Plan built with a caller-supplied encoding or preparation gets a new content_id each time that encoding or preparation is chosen. A deterministic circuit digest would make the content_id reproducible, and QPY bytes cannot serve as that digest because they differ between processes. The encodings that a Method chooses are already deterministic.