Can a quantum PDE algorithm avoid the polynomial cost of resolving a fine spatial mesh? For standard fixed-order discretizations, direct classical methods require work polynomial in $h^{-1}$, or equivalently in the number of spatial degrees of freedom $N_h=Θ(h^{-d})$. Direct quantum implementations of a parabolic semigroup still have coherent complexity $\widetilde{\mathcal O}(\sqrt{T}/h)$, and gradient-dependent observables such as heat flux and dissipation introduce additional mesh dependence. Decay of the solution norm will further suppress the postselection probability for preparing a normalized final state.
We develop a multilevel quantum algorithm that estimates linear and quadratic observables $directly$ and places the fine--coarse cancellation inside the circuit before measurement. A contour-based LCU reconstructs each target-time correction from a coherent family of shifted resolvent differences. Rather than block encoding the fine and coarse inverses separately, we encode their difference through a shifted Ritz--Schur factorization, exposing its $\mathcal O(h_\ell^2)$ two-grid normalization. For Fourier hierarchies, the corresponding SELECT oracle consists of a quantum Fourier or sine transform, a spectral-band selector, and reversible diagonal arithmetic. We also give a non-Fourier realization based on energy-orthogonal dyadic midpoint details in one dimension, together with structured tensor-product extensions under fixed-rank coefficient and access assumptions.
For readouts with derivative order $0\leχ\le2$, optimized amplitude estimation removes $all$ polynomial dependence on the finest mesh size. Under the stated access assumptions, both linear and quadratic observables can be estimated with complexity $\widetilde{\mathcal O}(1+(Tε)^{-1})$, with only polylogarithmic dependence on $h^{-1}$.
Can a quantum PDE algorithm avoid the polynomial cost of resolving a fine spatial mesh? For standard fixed-order discretizations, direct classical methods require work polynomial in h−1, or equivalently in the number of spatial degrees of freedom Nh=Θ(h−d). Direct quantum implementations of a parabolic semigroup still have coherent complexity O(T/h), and gradient-dependent observables such as heat flux and dissipation introduce additional mesh dependence. Decay of the solution norm will further suppress the postselection probability for preparing a normalized final state.
We develop a multilevel quantum algorithm that estimates linear and quadratic observables directly and places the fine--coarse cancellation inside the circuit before measurement. A contour-based LCU reconstructs each target-time correction from a coherent family of shifted resolvent differences. Rather than block encoding the fine and coarse inverses separately, we encode their difference through a shifted Ritz--Schur factorization, exposing its O(hℓ2) two-grid normalization. For Fourier hierarchies, the corresponding SELECT oracle consists of a quantum Fourier or sine transform, a spectral-band selector, and reversible diagonal arithmetic. We also give a non-Fourier realization based on energy-orthogonal dyadic midpoint details in one dimension, together with structured tensor-product extensions under fixed-rank coefficient and access assumptions.
For readouts with derivative order 0≤χ≤2, optimized amplitude estimation removes all polynomial dependence on the finest mesh size. Under the stated access assumptions, both linear and quadratic observables can be estimated with complexity O(1+(Tε)−1), with only polylogarithmic dependence on h−1.
Spatial derivatives pose a basic challenge for quantum algorithms for partial differential equations: their standard discretizations have operator norms that grow as the mesh is refined, leading to large block-encoding normalization factors. Even when preconditioning improves access to the solution, estimating derivative-dependent physical observables can reintroduce polynomial dependence on the inverse mesh width. We develop a multilevel quantum algorithm for estimating linear and quadratic observables of elliptic equations that addresses this readout cost. The algorithm decomposes the output into corrections across nested discretizations and realizes fine--coarse cancellation coherently, before measurement. A Ritz-complement factorization connects the O(h2) size of the corrected response with a block encoding at the same scale, allowing this decay to offset growing readout normalizations. Under efficient quantum access to the corrected coordinates, input data, and readouts, observables with readout scale O(h−χ), 0≤χ≤2, can be estimated to additive accuracy ε at cost O(ε−1) using amplitude estimation or O(ε−2) using direct sampling, with only polylogarithmic dependence on the inverse finest mesh width hL−1. We give explicit constructions of the required operator oracles for one-dimensional piecewise linear finite elements with diffusion piecewise constant on a fixed coarsest partition, and for Fourier spectral hierarchies. Output-specific bias estimates connect the discrete results to continuum accuracy.
Quantum linear systems algorithms (QLSAs) aim to solve linear systems $A\bb{x}=\bb{b}$ exponentially faster than classical methods under certain conditions. In this work, we develop quantum algorithms for solving linear algebraic equations from an ODE-based perspective. Inspired by the linear combination of Hamiltonian simulation (LCHS) representation in the Fourier approach \cite{Childs2017QLSA}, we express the solution $\bb{x}$ as a linear combination of solutions to a system of linear convection equations, which become Schrödinger-type equations with unitary evolutions in the Fourier domain. We refer to this representation as LC-Schrödingerization. Based on this result, we construct an LCHS-based quantum algorithm with two LCHS instances: one for time-marching and one for numerical integration. The key construction uses the derivative of a Gaussian-smoothed hat function and recovers the solution over a fixed auxiliary interval. This permits a truncation time independent of the target accuracy and avoids the loss in success probability from selecting a single grid point. Periodization and explicit Fourier coefficients provide the corresponding projection error bounds. Under the stated oracle assumptions and given a constant-factor estimate of the solution norm, direct simulation of the select operators and block preconditioning achieve the optimal matrix-query complexity O(κAlogε1) without using variable-time amplitude amplification (VTAA). The same upper bound holds for queries to the right-hand-side preparation oracle.
Quantum PDE solvers are difficult to evaluate in practice because published studies use different discretizations, output models, reconstruction rules, and hardware assumptions. We present a reproducible, application-driven benchmark for the 1-D Dirichlet heat equation that compares eleven kernels under the same problem instances and readout contract. The benchmark covers coherent linear solvers (HHL, QSVT, and QLS-Fourier), VQLS, imaginary-time methods (QITE, var-QITE, and AVQDS), real-time Hamiltonian simulation and unitary dilations (Hamiltonian simulation, Schade-Hamiltonian, and Schr"odingerisation), and the spectral quantum simulation method (QSM). We use three initial conditions, four grid sizes from n=4 to 7 qubits (N=16 to 128), a CFL-like ratio r≈0.4, and final time T=1. Statevector, ideal-shot (105 shots per step), and noisy Aer backends separate algorithmic, sampling, and device-noise errors. On statevector, QSM and Schade-Hamiltonian reproduce the semi-discrete reference to floating-point precision, Schr"odingerisation reaches approximately 10−4 error, and QITE is the strongest non-transform method for smooth data. Under the fixed-shot setting, HHL degrades to approximately 0.79 relative ℓ2 error, while several low-depth or postselected methods become readout-limited. A norm-mismatch ablation attributes 23--29% of the n=7 smooth-initial-condition error of Hamiltonian simulation, AVQDS, and QLS-Fourier to reconstruction normalization. Compact observables, including total thermal energy and individual Fourier-mode weights, require 1--3 orders of magnitude fewer shots than full-field reconstruction. The resulting public benchmark provides a practical guide for selecting quantum PDE solvers.
This paper presents an efficient quantum spectral-Galerkin method for constant-coefficient second-order linear elliptic equations. The method employs a boundary-adapted Jacobi tensor-product basis that incorporates the boundary condition into the approximation space and produces Galerkin systems with polynomially bounded condition numbers and sparsity. For analytic solutions satisfying a suitable uniform regularity condition, spectral convergence allows the normalized Galerkin coefficient state to be prepared with complexity poly(d,log(1/ε)), where d is the spatial dimension and ε is the error tolerance, under standard sparse-access and state-preparation assumptions. At the state-preparation level, this improves the precision dependence over fixed-order finite-difference and finite-element discretizations. We further show how to estimate prescribed linear functionals of the solution, with an additional O(1/ε) factor arising from norm and inner-product estimation.
Multi-asset option pricing under local- and stochastic-volatility models leads naturally to high-dimensional parabolic PDEs. We develop an end-to-end quantum PDE framework for European option pricing under local-volatility Black--Scholes and Heston models. The framework takes classical contract and model data as input and returns classical estimates of selected option values. We solve the pricing PDEs after finite-difference discretization on spatial grids. For N=2n grid points per spatial direction and d assets, the end-to-end gate complexity for single-point recovery, counted in elementary CNOT gates and one-qubit Pauli-axis rotations, has leading grid-size dependence O(d2N2+d/2) for local-volatility Black--Scholes and O(d2Nd+2) for Heston. Relative to grid-based finite-difference baselines, these scalings correspond to polynomial improvement factors Nd/2 and Nd, respectively. These estimates translate to Clifford+T resources via standard compilation. We complement the complexity analysis with numerical benchmarks against standard classical methods. In the Heston setting, the framework recovers option prices across strikes together with the associated implied-volatility smile/skew. Overall, this work provides a complete end-to-end quantum pricing pipeline with explicit resource accounting and theoretical performance guarantees.
We present a structure-aware construction of quantum BPX preconditioners for Q1 finite element discretizations of constant-coefficient reaction--diffusion equations on the unit cube, including the Poisson equation. Building on existing quantum finite element and finite difference methods, we express the preconditioned operator through explicit one-dimensional maps, tensor products, and norm-preserving interlevel transfers. Symmetric preconditioning gives multilevel coefficients, and its favorable conditioning alone does not ensure efficient recovery of a state in the original finite element nodal basis. We instead use the local L2-orthonormal basis underlying the construction as the output representation of the same finite element function. Its coefficient norm equals the continuous L2 norm. Assuming direct access to the preconditioned right-hand side, we prepare this orthonormal coefficient state with a represented function error at most ε. The expected query complexity is O(dL3∥u∥L2∥f∥L2log(2+εL∥f∥L2)), where u is the exact solution, f is the right-hand side, d is the spatial dimension, and L=log2(1/h) for mesh width h. The bound holds for 0<ε≤∥u∥L2/2 and depends only polylogarithmically on mesh refinement for fixed continuous data. We also analyze linear-functional estimation and illustrate the construction and continuous error bounds with a one-dimensional UnitaryLab experiment.