Quantum PDE Solvers in Practice: Application-Driven Benchmarking of the Heat Equation

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Quantum PDE Solvers in Practice: Application-Driven Benchmarking of the Heat Equation

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

softwarearchitecture softwareengineering testing

software engineering, software architecture, testing

Quantum PDE Solvers in Practice: Application-Driven Benchmarking of the Heat Equation Mahmoud ElKarargy1 , Abdelaziz Rahwan1 , Abdelrahman Elsayed1 , Forat Hatem1 1

arXiv:2607.12688v1 [quant-ph] 14 Jul 2026

Brightskies Technologies, Alexandria, Egypt {mahmoud.elkarargy, abdelaziz.rahwan}@brightskiesinc.com {abdo.elsayd102, forathatem2004}@gmail.com Abstract—Quantum PDE solvers are difficult to evaluate in practice because the literature compares algorithms under different discretisations, output models, reconstruction rules, and hardware assumptions. This paper contributes a reproducible, application-driven benchmark for the 1-D Dirichlet heat equation that places eleven kernels from five paradigm classes under the same problem instance and readout contract: (A) coherent linear solvers (HHL, QSVT, QLS-Fourier), (B) a NISQ variational linear solver (VQLS), (C) imaginary-time variational dynamics (QITE, var-QITE, AVQDS), (D) real-time Hamiltonian simulation and unitary dilations (Hamiltonian simulation, Schade-Hamiltonian, Schrödingerisation), and (E) a spectral/transform method (QSM). The harness fixes three initial conditions (pulse, Gaussian, bimodal), four grid scales (n=4–7 qubits, N =16–128), a CFL-like ratio r≈0.4, and a final time T =1. Three backends (statevector, ideal shots at 105 /step, and noisy Aer) isolate algorithmic, sampling, and device-noise error; backend-specific kernel omissions are stated in §VI-B. The central finding is that QSM and Schade-Hamiltonian reproduce the semi-discrete reference to floating-point precision on statevector; Schrödingerisation reaches ∼ 10−4 ; and QITE is the strongest non-transform method on smooth data. Conversely, in our fixed-shot setting, HHL degrades to ∼ 0.79 relative ℓ2 error, and several low-depth or postselected methods become readout-limited before their unitary cores are exhausted. A normmismatch ablation shows that the shared smooth-initial-condition error plateau of Hamiltonian simulation, AVQDS, and QLSFourier is partly a reconstruction-normalisation artifact (23–29% of their n=7 smooth-IC error). Finally, compact observables such as total thermal energy and individual Fourier-mode weights require 1–3 orders of magnitude fewer shots than full-field reconstruction on the same circuits. The result is a public benchmark artifact and a practical selection guide: spectral/dilation kernels for reference accuracy, QITE for moderate-depth smooth full-field studies, and compact observable extraction for credible application-level quantum benefit. Index Terms—quantum computing, quantum algorithms, partial differential equations, heat equation, benchmarking, Hamiltonian simulation, Schrödingerisation, quantum linear systems, variational quantum algorithms, quantum singular value transformation, NISQ

I. I NTRODUCTION Partial differential equations (PDEs) underpin models across physics, engineering, and the computational sciences [1], and within the quantum-algorithms literature have emerged as an important application domain because they exercise state

preparation, evolution, inversion, and measurement in a single problem. Yet that literature is difficult to compare: methods differ in discretisation, computational mechanism, output representation, and hardware assumption, making it unclear how to relate conceptual promise to practical performance or how to interpret quantum-advantage claims on a common basis. Among candidate benchmark problems, the heat equation provides a particularly useful lens through which to organize this landscape. In its general parabolic form, it may be written as ∂u = α ∆u, x ∈ Ω ⊂ Rd , t > 0, (1) ∂t where ∆ denotes the Laplacian operator and α > 0 is the diffusion coefficient. This form captures the broader class of diffusion-type dynamics relevant to parabolic PDEs. In this paper, we specialize to the one-dimensional case with Dirichlet boundary conditions, ∂u ∂2u = α 2 , x ∈ (0, L), t > 0, (2) ∂t ∂x which is analytically tractable and simple enough to support controlled comparison across multiple quantum algorithms. Three properties make this PDE unusually well-suited as a cross-paradigm benchmark: (i) it exercises all four primitive operations that distinguish quantum PDE solvers from one another—state preparation, matrix inversion, real- and imaginary-time evolution, and amplitude readout, within one instance; (ii) its closed-form Fourier sine-series solution makes the frequency content of each initial condition explicit, enabling initial condition (IC) based stress-testing of kernels that differ in their handling of high-mode dynamics (pulse excites a broad band; Gaussian is low-frequency dominated; bimodal isolates two discrete modes); and (iii) representative quantum treatments of the heat equation [7], [40]–[42], [44], together with broader linear parabolic PDE solvers [35], align this benchmark with active lines of research and a ready-made comparison literature. Together these properties make the heat equation a natural common reference problem for examining how different quantum paradigms approach PDEs. Building on Feynman’s observation [2] that quantum systems obey differential equations, the Harrow–Hassidim–Lloyd algorithm [3] opened a pathway to quantum linear-system

methods for implicit PDE discretisations. Subsequent work has produced improved linear-system and differential-equation methods, Hamiltonian-based formulations, and general PDE simulation frameworks [4]–[7], [41], [46]; in parallel, variational approaches [8], [10] offer hardware-accessible alternatives for the NISQ era [9]. What is missing from this landscape is an empirically grounded comparison of these approaches across a common problem setting. This paper asks a concrete empirical question: Under matched discretisation, solution reconstruction, and shotbudget assumptions, which quantum kernels accurately solve the heat equation, and what determines their characteristic failure modes? We implement eleven kernels spanning five paradigm classes in a single harness, evaluate them on three initial conditions, and compare statevector, ideal shot-based, and noisy Aer backends (§VI-B). To our knowledge, no prior benchmark spans all five paradigm classes on a shared PDE instance under a fixed reconstruction rule; the closest comparative studies each cover a strict subset of kernels, backends, or output maps (§II). We extend this line by adding Paradigm A coherent solvers, Paradigm E spectral methods, the ideal-shots/noisy-simulator backend ladder, and a normmismatch ablation that decomposes the shared smooth-IC error plateau into its algorithmic and reconstruction components, all under one output map. Key contributions. (C1) A controlled eleven-kernel benchmark spanning five paradigms under one PDE instance, discretisation, and reconstruction rule (§V–§VII), with a staged statevector→idealshots→noisy-Aer backend ladder that separates algorithmic, sampling, and device-noise error (§VI-B). (C2) A norm-mismatch ablation that attributes the smoothIC error plateau shared by Hamiltonian simulation, AVQDS, and QLS-Fourier to a reconstruction-normalisation failure, with 23–29% of the residual ℓ2 error at n=7 traced to pure norm drift independent of the unitary core (§VII-C). This reframes several apparent algorithmic ceilings as readoutpipeline artefacts and explains why kernels with very different mechanisms can produce the same error band. (C3) An observable-readout advantage analysis (§VII-D) showing that compact physical functionals (total energy, selected mode weights, and boundary-flux proxies) reach target precision at 1–3 orders of magnitude lower shot budget than full-field reconstruction on the same circuits. (C4) A reproducibility-oriented benchmark artifact: a modular framework, configuration files, raw outputs, and scripts for regenerating the tables and figures, with new kernels added through one documented interface method. II. R ELATED W ORK Prior comparative studies each cover a strict subset of the landscape. Tseng et al. [42] compare two heat-equation solvers at statevector level; Shayegan [32] benchmarks VQLS, HHL, and quantum annealing on time-fractional diffusion; Özgüler [36] evaluates VQE against VarQTE-, AVQDS-,

and Trotter-based kernels in a shared harness on advection– diffusion; Alipanah et al. [44] compare AVQDS and varQITE against a Trotter baseline; surveys [31], [37], [38] cover quantum linear solvers and the broader PDE landscape without implementing kernels. None of these span all five paradigms, combine statevector/shot/noisy backends, enforce a shared reconstruction rule, or quantify the reconstruction component of the measured error. Kernel-specific literature. Paradigm A builds on HHL [3], QSVT block encodings [15], [16], [34] (benchmarked against HHL in [33]), and the Fourier-LCU construction [43]; our fixed-CFL condition-number analysis (§IV-E) shows A is wellconditioned at all tested n, so we do not engage the amplitudeamplification [25] and preconditioning [24] refinements developed for the ill-conditioned regime. Paradigm B follows Bravo-Prieto [10] with dynamic-ansatz and finite-element extensions [11], [29] and related Poisson VQAs [12]. Paradigm C builds on Motta [18], extended to linear PDEs by Kumar and Wilmott [35], with McLachlan variational dynamics [47], [48] and the AVQDS construction of Alipanah et al. [44]. Paradigm D combines Schrödingerisation [40], [41], [46] with the Hamiltonian IBVP construction of Schade et al. [45] and mean-field Trotterisation from [44]. Paradigm E follows the DST primitives of Klappenecker and Rötteler [27] as specialised for non-periodic BVPs in Febrianto et al. [39]. III. BACKGROUND We use five standard primitives, each recalled where it is maps v ∈ RN to |ψ⟩ = P applied. Amplitude encoding n j (vj /∥v∥)|j⟩, at O(2 ) CNOTs for generic states [23]. QPE [20], [28] is the core subroutine of HHL; first-order Trotter–Suzuki Hamiltonian simulation [21], [22] drives its controlled-U blocks and the Paradigm-D constructions. QSVT applies a QSP-derived polynomial to the singular values of a block-encoded matrix [15]–[17]. McLachlan’s variational principle [47], [48] yields M (θ)θ̇ = −V (θ), driving var-QITE and AVQDS. Schrödingerisation lifts the non-unitary heat semigroup to a unitary on an auxiliary momentum register, recovering dissipation by post-selecting p > 0 [40], [41], [46]. IV. P ROBLEM F ORMULATION AND S ETUP We solve the 1-D heat equation (2) on x ∈ (0, L), L = 1, u(0, t) = u(L, t) = 0, with diffusivity α = 0.01. Second-order central differences on N interior points (∆x = L/(N +1)) yield the semi-discrete system du = −α D u, dt

(3)

where D ∈ RN ×N is the symmetric positive definite tridiagonal Laplacian (diagonal +2/(∆x)2 , off-diagonal −1/(∆x)2 ). The exact solution via eigendecomposition D = V ΛV T is u(t) = V diag(e−αλk t ) V T u(0), computed classically as the reference baseline.

(4)

A. Parameters n

Grid sizes N = 2 for n ∈ {4, 5, 6, 7} (N ∈ {16, 32, 64, 128}), each encoded into n system qubits. Perconfig ∆t is chosen so that the CFL-like ratio r = α∆t/(∆x)2 ≈ 0.4 is held constant across grid sizes. We integrate to a fixed horizon T =1; with r held fixed, each N determines ∆t and ∆x as above. The N interior amplitudes represent only x ∈ (0, L). B. Boundary conditions Homogeneous Dirichlet conditions are chosen because they make D symmetric positive definite with a closed-form sine eigenbasis—supplying the exact classical reference (4)—and because several kernels are ported in their published Dirichlet form (QSM’s DST primitives and the Schade dilation angles) [27], [39]. The choice is not neutral across paradigms: Neumann or Robin conditions perturb only the boundary rows of D, exchanging the sine basis for cosine or mixed bases—a DCT-type transform for the transform kernels [27], the physical-BC construction of [40] for Schrödingerisation— while the linear-system and Paradigm-C kernels see only a modified tridiagonal matrix with O(1) conditioning under fixed r. In the harness the BC enters only through the discretised operator and the boundary-value configuration, so alternative conditions are an extension of the operator stencil rather than a redesign; we report the Dirichlet instance on which all eleven kernels are simultaneously well-defined. C. Initial conditions Three profiles stress different parts of the Laplacian spectrum, quantified by the discrete sine coefficients P ûk (0)= j uj (0) sin(kπj/(N +1)): (i) a rectangular pulse of unit height over the central 20% of the interior, whose |ûk (0)| decays only as 1/k for odd k and therefore excites a broad band of eigenmodes; (ii) a Gaussian bump with σ=L/6, whose sine coefficients decay exponentially in k (lowfrequency dominated, |ûk |<10−4 |û1 | for k≥6); and (iii) a bimodal profile sin(πx/L) + 0.5 sin(3πx/L), scaled to unit peak, supported essentially on two discrete modes k∈{1, 3}— a multi-mode stress test with no high-frequency content. These three profiles are chosen as extreme points of sine-spectrum occupancy, the variable that drives each paradigm’s dominant error term: Trotter, locality-truncation, and momentumresolution errors grow with high-k weight (pulse), whereas smooth data isolates reconstruction and normalisation effects (Gaussian, bimodal; §VII-C). D. Output model All kernels encode temperatures into quantum amplitudes: |u⟩ = u/∥u∥. Reconstructing N encoded amplitudes to fixed per-component precision ϵ is in the regime of (partial) state tomography, with a total sampling cost that grows with N and as ϵ → 0 exhibits the familiar 1/ϵ2 Monte Carlo factor [26]. Correspondingly, practical “quantum speedup” accounts for PDEs are more credible when the target is a low-dimensional

set of observables of u than when the goal is a full field readout [19]. The statevector backend reads amplitudes directly and is the zero-shot limit of both output models; sampled backends price them differently. Computational-basis sampling recovers the full field at O(N/ϵ2 ) total shots (signed amplitudes would additionally require interferometric or amplitude-estimation circuits; the non-negative temperature fields P used here avoid that overhead). A compact observable O = j cj Pj instead needs a number of measurement settings independent of N , and operators diagonal in a structured basis (e.g., D in its sine eigenbasis) collapse to a single setting at the price of appending the transform circuit—a depth-versus-settings trade quantified in §VII-D. E. Condition number and scaling For the implicit-Euler matrix A = I + α∆t D, the eigenvalues are 1 + α∆t λk (D), so κ(A) =

1 + α∆t λmax (D) . 1 + α∆t λmin (D)

(5)

For the Dirichlet Laplacian on N interior points with ∆x = L/(N +1), the extremal eigenvalues are λmax (D) = 2 2 4 Nπ 4 π and λmin (D) = (∆x) 2 sin (∆x)2 sin 2(N +1) 2(N +1) . With ∆t = r(∆x)2 /α (so α∆t/(∆x)2 =r fixed), these give α∆t λmax (D) → 4r and α∆t λmin (D) = rπ 2 /(N +1)2 → 0 as N →∞, hence 1 + 4r −−−−→ 1 + 4r. (6) κ(A) ≈ 1 + rπ 2 /(N +1)2 N →∞ For r=0.4 this limit is κ(A) → 2.6, and κ is scale-stable under our fixed CFL ratio. The implicit-Euler system is therefore well-conditioned across n=4–7, which is favourable for HHL and VQLS; in practice the dominant error sources for those kernels are QPE resolution and ansatz expressibility respectively, not κ. V. K ERNELS U NDER T EST Within each paradigm, kernels share reconstruction assumptions and often exhibit similar cost or failure modes; subsections give idea, circuit mechanism, and accuracy bounds. Backend omissions are consolidated in §VI-B. A. Paradigm A — Quantum Linear Solvers Paradigm-A kernels perform coherent inversion: they apply a unitary transformation to a block encoding of A to produce the solution state |x⟩ ∝ A−1 |b⟩ directly on the quantum register, without classical optimisation. All three kernels discretise (2) in time with implicit Euler, A x = b,

A = I + α∆t D,

b = un ,

x = un+1 ,

(7)

and differ in how they realise the inversion step. HHL is the canonical quantum linear-system algorithm, recovering |x⟩ ∝ A−1 |b⟩ through eigenvalue inversion. It applies QPE with m=6 clock qubits to extract λ̃j , a controlled ancilla rotation C/λ̃j to encode the inverse, and postselects the ancilla to recover the solution [3]. The QPE

block dominates the circuit cost through repeated Hamiltoniansimulation–based controlled-eiAt unitaries; the overall depth depends on the cost of simulating A and the chosen implementation. Accuracy is bounded by the 2−m eigenvalue resolution, which sets the ancilla rotation angles. Among ParadigmA kernels, HHL represents the eigenvalue-by-eigenvalue inversion strategy that QSVT and QLS-Fourier replace with polynomial and Fourier alternatives. On the noisy Aer path, transpiling the deep controlled Hamiltonian-simulation/QPE stack to native two-qubit gates triggers Qiskit two-qubit unitary decomposition failures for n ∈ {5, 6, 7} (we observe this as a TwoQubitWeylDecomposition error during synthesis); we therefore include HHL on the noisy backend only at n=4. Statevector and ideal-shot runs remain for n=4– 7. QSVT applies a polynomial approximation of 1/x to the singular values of a block-encoded matrix. UA block-encodes Ã = A/∥A∥2 , and d+1 alternating UA /UA† queries interleaved with QSP phase rotations realise the polynomial at degree d = 15. The approximation is taken over the normalized singular-value interval [a, 1], where a = σmin (Ã) ≈ 1/κ(A). Oracle-call count is constant in n; circuit depth depends on the block-encoding cost. Accuracy is bounded by the degreed polynomial approximation error on [a, 1], which decays as O((1−a)d ) for well-conditioned A. [15], [16] . QLS-Fourier denotes our Fourier/LCU-style functional quantum linear-system implementation. We use this label for the benchmark kernel based on the linear-system framework of [43] and the functional-QLS comparison setting of [33]. It approximates 1/x by a truncated Fourier integral implemented as an LCU of Hamiltonian-simulation steps. We evaluate the post-selected LCU in the DST eigenbasis, where Ã is diagonal, reducing the cost to a scalar evaluation per eigenmode. Accuracy is bounded by the Fourier truncation and discretisation schedule of [33]. The Appendix-D schedule of [33] yields J ·2L ≈ 106 LCU terms at n=4, growing to ∼9 × 107 (n=6) and ∼1.6 × 109 (n=7); compiling a controlled eiÃtj per term is infeasible, and truncating to a tractable subset substitutes an ad-hoc cutoff for the algorithm’s error budget. We instead evaluate the ancilla|0⟩ block of PREPARE–SELECT–PREPARE† analytically in the eigenbasis of Ã, equivalent to a noiseless statevector run of the full LCU circuit. As this path bypasses both transpilation and sampling, QLS-Fourier is evaluated statevector-only (§VI-B). B. Paradigm B — NISQ Variational Solver VQLS leverages shallow, parameterized circuits suitable for current NISQ devices [10], [49]. It targets the same implicit-Euler system (7) by minimising a local Hadamardtest cost function whose global minimum coincides with |x⟩ ∝ A−1 |b⟩; a hardware-efficient ansatz whose depth adapts to the register (four–five layers at the tested sizes, matching Table II) is optimised classically and warm-started across time steps. In the Pauli basis the tridiagonal A is not sparse: its exact Hilbert–Schmidt expansion spreads over L=N =2n

Pauli strings (verified numerically for n=3–7); O(n)-term decompositions are known only over non-Pauli primitives such as σ ± [12]. Our implementation truncates the expansion at 10% of the leading coefficient, retaining L=4 dominant strings (I, X0 , X0 X1 , Y0 Y1 ) at every tested n at the cost of a ∼15% spectral-norm perturbation of A—a systematic bias shared by all VQLS results reported here. Accuracy is bounded by this truncation, by ansatz expressibility, and by classical-optimiser convergence in the non-convex landscape. VQLS needs no fault-tolerant resources, but its sampled cost is multiplicative (L2 Hadamard tests per cost evaluation × O(102 ) optimiser iterations × time steps), so we run it statevector-only (§VI-B). C. Paradigm C — Imaginary-Time Variational Dynamics Paradigm-C kernels approximate the heat semigroup e−αDt directly in imaginary time via a variational Pauli basis, bypassing both the linear system and any unitary dilation. The var-QITE and AVQDS kernels follow the constructions of Alipanah et al. [44]. QITE realises each imaginary-time step through a bestunitary approximation of e−H∆τ expressed in a local Pauli basis [18]. At each step, Pauli expectations ⟨σj ⟩ and ⟨σj H⟩ are measured, coefficients aj are obtained by solving a Grammatrix linear system, and state is updated via a firstQ the −i∆τ aj σj . The benchmark uses order Trotter product e j compressed QITE with locality support D = 3, chosen as the empirical accuracy–runtime sweet spot; the classical inner loop requires O(Nσ2 ) expectation evaluations per step, where Nσ is the number of retained Pauli terms in the compressed σ basis. Accuracy is bounded by the truncation support D and by the Trotter commutator error O(∥a∥2 ∆τ 2 ). QITE is the reference imaginary-time method against which the variational relaxations var-QITE and AVQDS are measured. Var-QITE is a variational relaxation of QITE that replaces the exact statevector with a parameterised ansatz [47], [48]. Parameters evolve under McLachlan’s principle M (θ)θ̇= − V (θ), trading QITE’s O(4n ) observable set for a compact p×p gradient system; accuracy is bounded by ansatz expressibility. AVQDS is an adaptive extension of var-QITE in which the ansatz is grown online [44]. A new Pauli rotation drawn from a classical operator pool is appended whenever the McLachlan residual exceeds a threshold; circuit depth is determined during simulation rather than in advance. Accuracy is bounded by the residual threshold and operator-pool coverage. D. Paradigm D — Real-Time Hamiltonian Simulation / Unitary Dilation Paradigm-D kernels handle the heat equation’s non-unitary dynamics through real-time evolution, either via a classical scalar correction (Hamiltonian simulation) or by embedding the dissipative dynamics into a unitary on an enlarged Hilbert space (Schade-Hamiltonian, Schrödingerisation). Hamiltonian simulation is the simplest real-time approach: e−iH∆t with H=αD is applied via Lie–Trotter and the dissipative decay is recovered classically—the real part of

the evolved state is rescaled by ∥T ∥ e−⟨H⟩∆t using the prestep energy [44]. Accuracy is bounded by the mean-field approximation: every mode relaxes at ⟨H⟩ at leading order in ∆t, with mode-specific corrections only at O(∆t2 ). No auxiliary register, no post-selection. Schade-Hamiltonian is a unitary-dilation approach that embeds the heat contraction e−αD∆t as the |0⟩anc block of a real-time unitary on an ancilla-extended register. We port the Cholesky/block-off-diagonal embedding of [45] from the elastic wave setting to the dissipative heat setting via exact per-step dilation angles θk = arccos(e−αλk ∆t ) computed in the Laplacian eigenbasis; U =e−iH is applied as a dense unitary and post-selected on the ancilla. Within this benchmark this is therefore a structured dilation reference: it tests the effect of an eigenbasis-aligned unitary embedding rather than claiming a sparse, hardware-optimal construction. Accuracy is limited by floating-point arithmetic in the eigendecomposition and by the stated post-selection/readout model, and per-mode decay is restored exactly. Schrödingerisation (SZ) is the second unitary-dilation kernel. Instead of a single ancilla, SZ appends an np -qubit momentum register (np =n) on which the heat equation lifts to a conservative advection equation; post-selecting on p>0 recovers the dissipative dynamics [40], [41], [46]. The compiled DiagonalGate yields O(2n ) two-qubit gates, so SZ is statevector-only (§VI-B). Accuracy is bounded by momentumregister resolution and the post-selection cut-off. E. Paradigm E — Spectral / Transform-Based QSM (Spectral) implements the exact propagator (4) as a three-stage circuit: forward DST into the Laplacian eigenbasis; multiplexed R√ y encoding the eigenmode decay |k⟩|0⟩ → e−αλk t |k⟩|0⟩+ 1−e−2αλk t |k⟩|1⟩; inverse DST [27], [39]. Compiled via the QFT-based multiplexed-Ry decomposition, the two-qubit cost scales as O(n2 ) in the structured implementation we count. The decay is encoded using one work qubit; extracting the contracted heat state therefore requires conditioning on the work-qubit outcome or, for observables, incorporating the success probability into the estimator. QSM serves as the benchmark’s structural reference for per-mode exact decay rather than as a claim of near-term end-to-end advantage. VI. E XPERIMENTAL S ETUP A. Computational Environment All experiments ran on a 48-core Intel Xeon Platinum 8260M workstation (two sockets, 2.40 GHz) with 186 GB RAM, CentOS 8, Python 3.10. The software stack uses Qiskit 2.3+ [30] with Qiskit-Aer 0.17+ for shot-based and noisy simulation, Qiskit IBM Runtime 0.45+ for the real-hardware interface, and the standard NumPy 1.24+ / SciPy 1.10+ / Matplotlib 3.10+ scientific stack. The n=4– 7 range is set by the full factorial design rather than simulator capability: at n=7, noisy-simulator runs already take on the order of days for several kernels (§VII), while statevector-only

studies of individual kernels can go higher within the same harness. B. Simulation Backends The framework implements four BaseBackend classes: statevector, ideal shot-based simulation, noisy simulation, and real hardware. We report the first three; the hardware interface is implemented but left for future work. Table I summarises the evidence level associated with each backend. A dash in the results tables denotes that a run was not performed for the stated implementation reason, not that the method has zero error. TABLE I: Backend evidence levels used in the benchmark. Backend

Role and limitations

Direct double-precision amplitude evolution. Isolates intrinsic algorithmic and reconstruction error (Trotter truncation, polynomial degree, ansatz expressibility, Fourier truncation, dilation error). The SV columns of Tables III–IV are produced here. Noiseless AerSimulator sampling at Nshots =100,000 √ per step (sampling floor ∼ 1/ Nshots ≈ 3 × 10−3 ). For variational/dynamics kernels, only expectation estimates are sampled; classical solves and reconstruction remain classical. Same transpiled circuits as Id (not retranspiled, so accuracy differences are attributable to noise alone), with a portable gate/readout model: 10−3 single-qubit depolarising, 10−2 two-qubit depolarising, and 10−2 symmetric readout bitflip. T1 /T2 and coupling-map constraints are omitted. IBM Runtime Sampler/Estimator via create_backend("hardware"); results are reserved for future work.

The Id backend omits VQLS, Schrödingerisation, and QLS-Fourier for the implementation reasons given in §V-B, §V-D, and §V-A. The noisy backend uses the same shotcompatible kernel set, except that HHL is included only at n=4 because larger noisy HHL circuits fail during two-qubit synthesis. VII. C OMPARATIVE E VALUATION AND R ESULTS We present results across all three backends and initial conditions (Figs. 1–3). Table II gives native depth and 1Q/2Q counts by kernel and n; Tables III and IV consolidate terminal relative ℓ2 error at n=4 and n=7 respectively, for every initial condition and backend column. Although several kernels take on the order of days to complete a single n=7 run with the noisy simulator, every kernel admitted by each backend finished, and the figures display the full n=4–7 range. A. Time Evolution and ℓ2 Error Plots Each backend figure (Figs. 1–3) plots temperature T (x) at five milestones per IC and grid size (top row) and relative ℓ2 error vs. simulation progress (bottom row). Table IV’s SV columns match the n=7 terminal values of the bottom-row curves.

Temperature Distribution T(x) All Methods, Configs & Time % (Statevector simulation, = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

25%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

Temperature Distribution T(x) All Methods, Configs & Time % (Statevector simulation, = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

25%

n=4 (4 qb) t = 0.2768 s

50%

n=4 (4 qb) t = 0.5536 s

75%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

Temperature Distribution T(x) All Methods, Configs & Time % (Statevector simulation, = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

25%

n=4 (4 qb) t = 0.2768 s

50%

n=4 (4 qb) t = 0.5536 s

75%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

1.0

0.8

0.6

0.4

0.2 0.4

0.6

n=5 (5 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.2

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.2

0.4

75%

0.6

0.8

n=5 (5 qb) t = 0.7346 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) 100% t = 0.9917 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

n=5 (5 qb) t = 0.0000 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

n=4 (4 qb) T(x)

1.0 0.8

0.6

n=4 (4 qb) T(x)

1.0 0.8

0.6

0.2

0.4

75%

0.6

0.8

n=5 (5 qb) t = 0.7346 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) 100% t = 0.9917 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.0 0.0

0.2

0.4

75%

0.6

0.8

n=5 (5 qb) t = 0.7346 s

1.0

0.8

0.6

0.4

0.2

0.4

0.6

n=6 (6 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.2

0.0 0.0

0.2

0.4

50%

1.0

0.4

0.2

0.0 0.0

25%

1.0

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.2

0.4

0.6

0.8

n=6 (6 qb) t = 0.7479 s

1.0

0.8

0.6

0.4

0.2

0.4

0.2

0.0 0.2

0.4

0.6

n=7 (7 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.4

0.6

n=6 (6 qb) 100% t = 0.9941 s

0.8

1.0

0.2

0.4

75%

0.6

0.8

n=7 (7 qb) t = 0.7500 s

1.0

0.2

0.4

0.6

0.8

n=6 (6 qb) t = 0.0000 s

1.0

0.2

0.4

0.6

n=7 (7 qb) 100% t = 0.9999 s

0.8

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.4 0.2

0.0 0.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.2

0.4

0.6

0.8

n=6 (6 qb) t = 0.7479 s

1.0

0.8

0.6

0.4

1.0

0.4

0.2 0.2

0.4

0.6

0.8

n=7 (7 qb) t = 0.0000 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.4

0.6

n=6 (6 qb) 100% t = 0.9941 s

0.8

1.0

0.2

0.4

75%

0.6

0.8

n=7 (7 qb) t = 0.7500 s

1.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.0000 s

0.8

1.0

0.2

0.4

0.6

n=7 (7 qb) 100% t = 0.9999 s

0.8

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.4

75%

0.6

0.8

n=6 (6 qb) t = 0.7479 s

1.0

0.8

0.6

0.4

0.2 0.4

0.6

n=7 (7 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.2

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.2

0.4

75%

0.6

0.8

n=7 (7 qb) t = 0.7500 s

1.0

0.8

0.6

0.4

0.2

0.4

0.2 0.2

0.4

x (m)

0.6

0.8

1.0

0.2

0.4

x (m) Classical (spectral) QSM

0.6

0.8

1.0

HHL VQLS

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

QITE Hamiltonian simulation

0.6

x (m)AVQDS

var-QITE

0.8

1.0

0.2

0.4

0.6

x (m) QLS-Fourier Schade Hamiltonian

0.8

1.0

0.4

0.2

0.0 0.0

QSVT SZ (Schrödingerization)

0.2

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

x (m) Classical (spectral) QSM

0.6

0.8

1.0

HHL VQLS

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

QITE Hamiltonian simulation

0.6

x (m)AVQDS

var-QITE

0.8

1.0

0.2

0.4

0.6

x (m) QLS-Fourier Schade Hamiltonian

0.8

1.0

0.2

(b) Gaussian — time evolution Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=5 (5 qb)

n=6 (6 qb)

100

n=7 (7 qb)

100

n=4 (4 qb)

100

n=5 (5 qb)

100

n=6 (6 qb)

100

0.4

x (m)

0.6

0.8

1.0

0.2

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

x (m) Classical (spectral) QSM

0.6

0.8

1.0

HHL VQLS

0.2

0.4

0.6

x (m)AVQDS

var-QITE

0.8

1.0

0.2

0.4

0.6

x (m) QLS-Fourier Schade Hamiltonian

0.8

n=5 (5 qb)

100

n=6 (6 qb)

100

10 2

10 4

QSM HHL VQLS QITE Hamiltonian simulation AVQDS var-QITE QSVT SZ (Schrödingerization) QLS-Fourier Schade Hamiltonian

25%

50%

Simulation %

75%

100%

10 6

10 8

10 10

10 12

25%

50%

Simulation %

75%

100%

25%

50%

Simulation %

75%

100%

(d) Pulse — ℓ2 error

25%

50%

Simulation %

75%

100%

QSM HHL VQLS QITE Hamiltonian simulation AVQDS var-QITE QSVT SZ (Schrödingerization) QLS-Fourier Schade Hamiltonian

10 8

25%

50%

Simulation %

75%

100%

10 6

Relative L2 error

10 2

10 4

10 12

10 6

10 8

10 10

10 12

25%

50%

Simulation %

75%

100%

25%

50%

Simulation %

75%

100%

(e) Gaussian — ℓ2 error

25%

0.6

0.8

1.0

0.6

0.8

1.0

0.4

x (m)

50%

Simulation %

75%

100%

QSM HHL VQLS QITE Hamiltonian simulation AVQDS var-QITE QSVT SZ (Schrödingerization) QLS-Fourier Schade Hamiltonian

10 8

25%

50%

Simulation %

75%

100%

10 6

10 8

10 10

10 12

25%

50%

Simulation %

75%

100%

n=7 (7 qb)

100

10 2

10 4

10 8

1.0

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=4 (4 qb)

100

10 2

10 4 10 6

0.2

0.4

n=7 (7 qb) 100% t = 0.9999 s

0.0 0.0

QSVT SZ (Schrödingerization)

10 2

10 10

0.0

1.0

0.2

0.0 0.0

QITE Hamiltonian simulation

0.8

0.4

0.2

0.0 0.0

0.6

100

0.4

0.2

0.0 0.0

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time) 100

0.4 0.2

0.0 0.0

QSVT SZ (Schrödingerization)

(a) Pulse — time evolution n=4 (4 qb)

100

0.4

0.2

0.0 0.0

n=7 (7 qb) T(x)

1.0

0.0

0.2

0.4

n=6 (6 qb) 100% t = 0.9941 s

0.0 0.0

0.8

0.4

0.0

1.0

0.2

0.0 0.0

0.8

0.4

0.2

0.0 0.0

0.2

0.6

1.0

0.8

0.0

0.4

n=5 (5 qb) 100% t = 0.9917 s

0.0 0.2

0.6

0.4

0.2

0.4

0.0

1.0

0.8

1.0

0.0

0.2

0.0 0.0

1.0

0.6

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.2

0.4

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.8

0.4 0.2

0.0 0.0

75%

1.0

0.6

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.2

0.4

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.6

0.4 0.2

0.0 0.0

75%

1.0

n=6 (6 qb) T(x)

0.4 0.2

0.0

n=6 (6 qb) T(x)

0.4 0.2

n=5 (5 qb) T(x)

1.0

n=5 (5 qb) T(x)

1.0

n=7 (7 qb) T(x)

n=5 (5 qb) T(x)

75% 1.0 0.8

0.0

n=6 (6 qb) T(x)

n=4 (4 qb) t = 0.5536 s

0.6

0.2

n=7 (7 qb) T(x)

50% 1.0 0.8

0.6

0.0

Relative L2 error

n=4 (4 qb) t = 0.2768 s

1.0 0.8

0.6

Relative L2 error

n=4 (4 qb) T(x)

0% 1.0 0.8

25%

50%

Simulation %

75%

100%

10 12

25%

50%

Simulation %

75%

100%

(f) Bimodal — ℓ2 error

Fig. 1: Statevector backend. Top row: temperature evolution T (x) at five time milestones across grid sizes n=4–7, one panel per initial condition. Bottom row: relative ℓ2 error vs. simulation progress (%). In each panel, the eleven quantum kernels appear as separate colored curves. kernel-level behaviour is discussed in §VII.

1) Statevector Backend: Figure 1 provides the clearest view of intrinsic algorithmic behavior because it removes both sampling noise and device noise. The statevector results show a clear hierarchy. QSM and Schade-Hamiltonian are the most accurate methods overall: across the three initial conditions—high-frequency (pulse), smooth (Gaussian), and smooth multi-mode (bimodal)—they overlap the classical baseline to floating-point precision and stay on the error floor throughout the run (Table IV, SV). Their success is structural: both use a basis aligned with the Laplacian eigenstructure. Schrödingerisation is the next most accurate method, but it no longer sits on the machine-precision floor. It remains very close to the reference, indicating that the unitary-dilation route captures the heat dynamics well, but its more general embedding leaves a visible residual gap to the spectral reference methods. QITE is the most accurate non-transform kernel. It performs best on the smooth Gaussian and bimodal benchmarks, and degrades on the high-frequency pulse as n increases (local truncation and high-frequency content). The middle-to-lower tier: HHL and QSVT are inversionlimited (phase estimation vs. polynomial degree); their 2Q counts in Table II grow fastest among included kernels. VQLS fades with n (capped-depth ansatz); var-QITE stays below QITE (Table IV, SV) on the reduced variational family. AVQDS, Hamiltonian simulation, and QLS-Fourier share a smooth-profile error band (Table IV, SV Gaussian/bimodal); on pulse, AVQDS and QLS-Fourier saturate while Hamiltonian simulation stays lower but still large. QLS-Fourier in Table II is listed only for n=4 (decomposition cap), independent of that failure mode. 2) Ideal (Shot-Based) Backend : Figure 2 keeps the broad√ statevector ranking under finite shots; the floor near ∼1/ Nshots caps how low errors can go once sampling dominates. The ideal (Id) columns in Table IV match the

plateaus (e.g. QSM and Schade-Hamiltonian in the 10−2 – 10−1 range, and the shared Gaussian plateau for Hamiltonian simulation and AVQDS; QLS-Fourier shares the same SV plateau but is not sampled here for the reasons in §V-A). QSM and Schade-Hamiltonian still hug the reference on smooth data; QITE is still the strongest non-transform kernel. The pulse separates methods more clearly in the time-evolution and error panels as n increases. HHL/QSVT and the middle group keep the same ordering as on SV, compressed by statistics rather than reordered. 3) Noisy Simulator Backend: Figure 3 shows the noisybackend runs: device noise raises the error floor over ideal shots for the successfully compiled subset (§VI-B). Filled Ns cells in Table IV are the n=7 terminal errors; “—” lines up with kernels dropped for the same limits. B. Circuit Cost and Measurement Budget Figure 4 (log 2Q counts vs. n) reflects the spread already tabulated: HHL, QSVT, and the spectral/dilation routes climb sharply while others stay moderate, so compiled cost sits beside the terminal accuracies of Table IV as a first-class axis. Figure 5 estimates the shot count to reach ϵ = 10−2 relative ℓ2 at interior sizes N =16–128 (as in the tables). Postselection and small success probability can dominate readout: HHL and QSVT need more shots (Fig. 5), alongside heavy 2Q use (Table II) and large Id-column errors (Table IV); see §V-A. Read as per-qubit growth factors, Table II also bounds extrapolation beyond the tested range: each added qubit multiplies compiled 2Q counts by ≈3.2 (HHL), ≈4 (QSVT, spectral/dilation kernels), and ≈2.6 (Hamiltonian simulation), while the variational circuits grow only mildly (≈1.3×). At ≈4× per qubit, two additional qubits already exceed the largest circuits compiled here by an order of magnitude, whereas the variational and QITE families stay circuit-cheap

Temperature Distribution T(x) All Methods, Configs & Time % (Ideal simulator (noiseless, shot-based), = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

25%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

Temperature Distribution T(x) All Methods, Configs & Time % (Ideal simulator (noiseless, shot-based), = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

25%

n=4 (4 qb) t = 0.2768 s

50%

n=4 (4 qb) t = 0.5536 s

75%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

Temperature Distribution T(x) All Methods, Configs & Time % (Ideal simulator (noiseless, shot-based), = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

25%

n=4 (4 qb) t = 0.2768 s

50%

n=4 (4 qb) t = 0.5536 s

75%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

1.0

0.8

0.6

0.4

0.2 0.4

0.6

n=5 (5 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.2

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.2

0.4

75%

0.6

n=5 (5 qb) t = 0.7346 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) 100% t = 0.9917 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

n=5 (5 qb) t = 0.0000 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

n=4 (4 qb) T(x)

1.0 0.8

0.6

n=4 (4 qb) T(x)

1.0 0.8

0.6

0.2

0.4

75%

0.6

n=5 (5 qb) t = 0.7346 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) 100% t = 0.9917 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.0 0.0

0.2

0.4

75%

0.6

n=5 (5 qb) t = 0.7346 s

0.8

1.0

0.8

0.6

0.4

0.2

0.4

0.6

n=6 (6 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.2

0.0 0.0

0.2

0.4

50%

1.0

0.4

0.2

0.0 0.0

25%

1.0

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.7479 s

0.8

1.0

0.8

0.6

0.4

0.2

0.4

0.2

0.0 0.2

0.4

0.6

n=7 (7 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.4

0.6

n=6 (6 qb) 100% t = 0.9941 s

0.8

1.0

0.2

0.4

75%

0.6

n=7 (7 qb) t = 0.7500 s

0.8

1.0

0.2

0.4

0.6

0.8

n=6 (6 qb) t = 0.0000 s

1.0

0.2

0.4

0.6

n=7 (7 qb) 100% t = 0.9999 s

0.8

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.4 0.2

0.0 0.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.7479 s

0.8

1.0

0.8

0.6

0.4

1.0

0.4

0.2 0.2

0.4

0.6

0.8

n=7 (7 qb) t = 0.0000 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.4

0.6

n=6 (6 qb) 100% t = 0.9941 s

0.8

1.0

0.2

0.4

75%

0.6

n=7 (7 qb) t = 0.7500 s

0.8

1.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.0000 s

0.8

1.0

0.2

0.4

0.6

n=7 (7 qb) 100% t = 0.9999 s

0.8

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.4

75%

0.6

n=6 (6 qb) t = 0.7479 s

0.8

1.0

0.8

0.6

0.4

0.2 0.4

0.6

n=7 (7 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.2

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.2

0.4

75%

0.6

n=7 (7 qb) t = 0.7500 s

0.8

1.0

0.8

0.6

0.4

0.2

0.4

0.2 0.2

0.4

x (m)

0.6

0.8

1.0

0.2

0.4

0.6

(m) Classicalx (spectral) QSM

0.8

1.0

0.2

0.0 0.0

0.2

0.4

Hamiltonian simulation x (m) AVQDS

HHL QITE

0.4

0.2

0.0 0.0

0.6

0.8

1.0

var-QITE QSVT

0.2

0.4

0.6

x (m) Schade Hamiltonian

0.8

1.0

0.4

0.2

0.0 0.0

QLS-Fourier

0.2

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

(m) Classicalx (spectral) QSM

0.8

1.0

0.2

0.0 0.0

0.2

0.4

Hamiltonian simulation x (m) AVQDS

HHL QITE

0.4

0.2

0.0 0.0

0.6

0.8

1.0

var-QITE QSVT

0.2

0.4

0.6

x (m) Schade Hamiltonian

0.8

1.0

0.2

(b) Gaussian — time evolution Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=5 (5 qb)

n=6 (6 qb)

100

n=7 (7 qb)

100

n=4 (4 qb)

100

n=5 (5 qb)

100

n=6 (6 qb)

100

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time) 100

0.4 0.2

0.0 0.0

QLS-Fourier

(a) Pulse — time evolution n=4 (4 qb)

100

0.4

0.2

0.0 0.0

n=7 (7 qb) T(x)

1.0

0.0

0.2

0.4

x (m)

0.6

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

(m) Classicalx (spectral) QSM

0.8

1.0

0.2

0.4

Hamiltonian simulation x (m) AVQDS

0.6

0.8

1.0

var-QITE QSVT

0.2

0.4

0.6

0.8

x (m) Schade Hamiltonian

n=7 (7 qb)

100

n=5 (5 qb)

100

n=6 (6 qb)

100

10 2

10 4

10 6

25%

50%

Simulation %

75%

100%

10 6

10 8

10 10

10 12

25%

50%

75%

Simulation %

100%

25%

50%

Simulation %

75%

100%

25%

50%

75%

Simulation %

100%

10 6 QSM HHL QITE Hamiltonian simulation AVQDS var-QITE QSVT QLS-Fourier Schade Hamiltonian

10 8

25%

50%

Simulation %

75%

100%

(d) Pulse — ℓ2 error

10 6

Relative L2 error

10 2

10 4

QSM HHL QITE Hamiltonian simulation AVQDS var-QITE QSVT QLS-Fourier Schade Hamiltonian

10 6

10 8

10 10

10 12

25%

50%

75%

Simulation %

100%

25%

50%

75%

Simulation %

100%

25%

0.6

0.8

1.0

0.6

0.8

1.0

0.4

x (m)

50%

75%

Simulation %

100%

10 6 QSM HHL QITE Hamiltonian simulation AVQDS var-QITE QSVT QLS-Fourier Schade Hamiltonian

10 8

25%

50%

Simulation %

75%

100%

(e) Gaussian — ℓ2 error

10 6

10 8

10 10

10 12

25%

50%

75%

Simulation %

100%

n=7 (7 qb)

100

10 2

10 4

10 12

1.0

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=4 (4 qb)

100

10 2

10 8

0.4

n=7 (7 qb) 100% t = 0.9999 s

10 4 10 6

0.2

1.0

0.0 0.0

QLS-Fourier

10 2

10 10

0.0

0.8

0.2

0.0 0.0

HHL QITE

0.6

0.4

0.2

0.0 0.0

0.2

0.4

n=6 (6 qb) 100% t = 0.9941 s

0.0 0.0

0.8

0.4

0.0

1.0

0.2

0.0 0.0

0.8

0.4

0.2

0.0 0.0

0.2

0.6

1.0

0.8

0.0

0.4

n=5 (5 qb) 100% t = 0.9917 s

0.0 0.2

0.6

0.4

0.2

0.4

0.0

1.0

0.8

1.0

0.0

0.2

0.0 0.0

1.0

0.6

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.2

0.4

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.8

0.4 0.2

0.0 0.0

75%

1.0

0.6

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.2

0.0 0.0

0.2

0.4

0.2

0.0 0.0

0.4 0.2

0.0 0.0

1.0

0.6

0.4 0.2

0.0 0.0

75%

1.0

n=6 (6 qb) T(x)

0.4 0.2

0.0

n=6 (6 qb) T(x)

0.4 0.2

n=5 (5 qb) T(x)

1.0

n=5 (5 qb) T(x)

1.0

n=7 (7 qb) T(x)

n=5 (5 qb) T(x)

75% 1.0 0.8

0.0

n=6 (6 qb) T(x)

n=4 (4 qb) t = 0.5536 s

0.6

0.2

n=7 (7 qb) T(x)

50% 1.0 0.8

0.6

0.0

Relative L2 error

n=4 (4 qb) t = 0.2768 s

1.0 0.8

0.6

Relative L2 error

n=4 (4 qb) T(x)

0% 1.0 0.8

25%

50%

Simulation %

75%

100%

10 12

25%

50%

75%

Simulation %

100%

(f) Bimodal — ℓ2 error

Fig. 2: Ideal (shot-based, noiseless) backend (Nshots =100 000). VQLS, Schrödingerisation, and QLS-Fourier are omitted (statevector-only; §VI-B). 25%

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

0.8 0.6

0.4

0.2

0.0 0.2

0.4

0.6

n=5 (5 qb) t = 0.0000 s

0.8

1.0

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.2

0.4

75%

0.6

n=5 (5 qb) t = 0.7346 s

0.8

1.0

0.2

0.4

0.6

n=5 (5 qb) 100% t = 0.9917 s

0.8

50%

n=4 (4 qb) t = 0.5536 s

75%

1.0

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

1.0 0.8

0.8

0.6

0.4

1.0

0.2

0.0 0.0

0.2

0.4

0.6

0.8

n=5 (5 qb) t = 0.0000 s

1.0

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.2

0.4

75%

0.6

n=5 (5 qb) t = 0.7346 s

0.8

1.0

25%

0.2

0.4

0.6

n=5 (5 qb) 100% t = 0.9917 s

0.8

50%

n=4 (4 qb) t = 0.5536 s

75%

1.0

n=4 (4 qb) t = 0.6920 s

n=4 (4 qb) 100% t = 0.9689 s

1.0

0.8

0.6

0.4

0.2

0.0 0.0

n=4 (4 qb) t = 0.2768 s

1.0

0.8 0.6

0.2

0.0 0.0

Temperature Distribution T(x) All Methods, Configs & Time % (Noisy simulator, = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

1.0

0.2

0.0 0.0

1.0

0.8 0.6

0.0 0.0

n=4 (4 qb) t = 0.2768 s

0.8 0.6

1.0

0.2

0.0 0.0

0.2

0.4

0.6

n=5 (5 qb) t = 0.0000 s

0.8

1.0

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=5 (5 qb) t = 0.2571 s

0.8

1.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=5 (5 qb) t = 0.5142 s

0.8

1.0

0.0 0.0

0.2

0.4

75%

0.6

n=5 (5 qb) t = 0.7346 s

0.8

1.0

0.8

0.6

0.4

0.2

0.4

0.2 0.2

0.4

0.6

n=6 (6 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.2

0.4

75%

0.6

n=6 (6 qb) t = 0.7479 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=6 (6 qb) 100% t = 0.9941 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

n=6 (6 qb) t = 0.0000 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

n=5 (5 qb) T(x)

1.0 0.8

n=5 (5 qb) T(x)

1.0 0.8

0.0

0.2

0.4

75%

0.6

n=6 (6 qb) t = 0.7479 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=6 (6 qb) 100% t = 0.9941 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=6 (6 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=6 (6 qb) t = 0.2462 s

0.8

1.0

0.2

0.4

50%

0.6

n=6 (6 qb) t = 0.4923 s

0.8

1.0

0.2

0.4

75%

0.6

n=6 (6 qb) t = 0.7479 s

0.8

1.0

0.8

0.6

0.4

0.2 0.2

0.4

0.6

n=7 (7 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.2

0.4

75%

0.6

n=7 (7 qb) t = 0.7500 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=7 (7 qb) 100% t = 0.9999 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

n=7 (7 qb) t = 0.0000 s

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.0 0.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

n=6 (6 qb) T(x)

1.0

n=6 (6 qb) T(x)

1.0

0.0

0.2

0.4

75%

0.6

n=7 (7 qb) t = 0.7500 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=7 (7 qb) 100% t = 0.9999 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

0.6

n=7 (7 qb) t = 0.0000 s

0.8

1.0

0.4

0.2

0.0 0.0

0.2

0.4

25%

0.6

n=7 (7 qb) t = 0.2500 s

0.8

1.0

0.2

0.4

50%

0.6

n=7 (7 qb) t = 0.5000 s

0.8

1.0

0.2

0.4

75%

0.6

n=7 (7 qb) t = 0.7500 s

0.8

1.0

0.8

0.6

0.4

0.2

0.0

0.2

0.0 0.0

0.2

0.4

x (m)

0.6

0.8

1.0

0.2

0.0 0.0

0.2

0.4

0.6

(m) Classicalx (spectral) QSM

0.8

1.0

HHL QITE

0.2

0.0 0.0

0.2

0.4

Hamiltonian simulation x (m) AVQDS

0.6

0.8

1.0

var-QITE QSVT

0.2

0.4

0.6

x (m) Schade Hamiltonian

0.8

1.0

0.6

0.4

0.2

0.0 0.0

QLS-Fourier

0.6 0.4

0.2

0.0 0.0

0.2

0.4

x (m)

0.6

0.8

1.0

0.2

0.4

(a) Pulse — time evolution n=5 (5 qb)

100

n=6 (6 qb)

100

x (m)

0.6

0.8

1.0

0.2

0.0 0.0

0.2

0.4

0.6

(m) Classicalx (spectral) QSM

0.8

1.0

0.2

0.0 0.0

0.2

0.4

Hamiltonian simulation x (m) AVQDS

HHL QITE

0.6

0.8

1.0

var-QITE QSVT

0.2

0.4

0.6

x (m) Schade Hamiltonian

0.8

1.0

0.6

0.4

0.2

0.0 0.0

QLS-Fourier

0.6 0.4

0.2

0.0 0.0

0.2

0.4

x (m)

0.6

0.8

1.0

n=7 (7 qb)

100

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=4 (4 qb)

100

n=5 (5 qb)

100

n=6 (6 qb)

100

0.2

0.4

x (m)

0.6

0.8

1.0

0.2

0.0 0.0

0.2

0.4

0.6

(m) Classicalx (spectral) QSM

0.8

1.0

HHL QITE

0.2

0.4

Hamiltonian simulation x (m) AVQDS

0.6

0.8

1.0

var-QITE QSVT

0.2

0.4

0.6

0.8

x (m) Schade Hamiltonian

n=5 (5 qb)

100

n=6 (6 qb)

100

10 2

10 4

10 6

10 12

25%

50%

Simulation %

75%

100%

10 6

10 8

10 10

10 12

25%

50%

Simulation %

75%

100%

10 12

25%

50%

Simulation %

75%

100%

10 12

10 6

10 6 QSM HHL QITE Hamiltonian simulation AVQDS var-QITE QSVT QLS-Fourier Schade Hamiltonian

10 8 10 10

25%

50%

Simulation %

75%

10 12

100%

(d) Pulse — ℓ2 error

25%

50%

Simulation %

75%

100%

10 6

10 8

10 10

10 12

25%

50%

Simulation %

75%

100%

10 12

Relative L2 error

10 2

Relative L2 error

10 2

QSM HHL QITE Hamiltonian simulation AVQDS var-QITE QSVT QLS-Fourier Schade Hamiltonian

25%

50%

Simulation %

75%

100%

10 12

10 6

10 6 QSM HHL QITE Hamiltonian simulation AVQDS var-QITE QSVT QLS-Fourier Schade Hamiltonian

10 8 10 10

25%

0.2

1.0

0.4

0.6

0.8

1.0

0.6

0.8

1.0

n=7 (7 qb) 100% t = 0.9999 s

0.4

x (m)

50%

Simulation %

75%

100%

10 12

25%

50%

Simulation %

75%

(e) Gaussian — ℓ2 error

100%

10 6

10 8

10 10

10 12

25%

50%

Simulation %

75%

100%

10 12

n=7 (7 qb)

100

10 2

10 8

1.0

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=4 (4 qb)

100

10 2

10 6

0.0

0.8

0.0 0.0

QLS-Fourier

10 2

10 10

0.2

0.6

0.2

0.0 0.0

0.0

0.4

n=6 (6 qb) 100% t = 0.9941 s

100

0.2

0.0 0.0

(b) Gaussian — time evolution

Relative L2 Error (Loss) vs. Classical Spectral Baseline (quantum methods, x-axis = % of total simulation time)

n=4 (4 qb)

100

0.2

0.0 0.0

n=7 (7 qb) T(x)

1.0

0.6

1.0

0.0 0.0

1.0

0.4

0.8

0.2

0.0 0.0

0.2

0.6

0.4

0.2

0.0 0.0

0.0

0.4

n=5 (5 qb) 100% t = 0.9917 s

0.0 0.0

1.0

0.4

0.2

0.0 0.0

0.0

0.4

0.2

0.0 0.0

1.0

0.2

Relative L2 error

25% 1.0

0.8 0.6

0.2

0.0 0.0

Temperature Distribution T(x) All Methods, Configs & Time % (Noisy simulator, = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

1.0

0.2

0.0 0.0

1.0

0.8

n=7 (7 qb) T(x)

n=4 (4 qb) T(x)

75% 1.0

0.6

0.0

n=5 (5 qb) T(x)

n=4 (4 qb) t = 0.5536 s

0.8 0.6

0.0

n=6 (6 qb) T(x)

50% 1.0

0.8 0.6

0.2

n=7 (7 qb) T(x)

n=4 (4 qb) t = 0.2768 s

1.0

n=4 (4 qb) T(x)

Temperature Distribution T(x) All Methods, Configs & Time % (Noisy simulator, = 0.01, per-config dt, L = 1 m)

n=4 (4 qb) t = 0.0000 s

0.8 0.6

n=4 (4 qb) T(x)

0% 1.0

25%

50%

Simulation %

75%

100%

10 12

25%

50%

Simulation %

75%

100%

(f) Bimodal — ℓ2 error

Fig. 3: Noisy simulator backend HHL is included only at n=4 (noisy transpilation limit; §V-A). Device noise degrades accuracy relative to the ideal backend for the circuits that transpile; conclusions are limited to this portable depolarising/readout model and the successfully compiled subset. TABLE II: Native circuit resource counts per kernel and grid size from U+CX decomposition: depth (“D”), single-qubit (“1Q”), and CNOT (“2Q”) gates; numbers ≥ 104 in thousands (k). AVQDS is counted after one compute step (the adaptive pool is empty before it); QLS-Fourier admits a native decomposition at n=4 only (LCU schedules at n≥5 exceed decomposition memory limits). n = 4 (N = 16) Kernel HHL QSVT QLS-Fourier VQLS QITE var-QITE AVQDS Ham. sim. QSM Schrödingerisation Schade-Ham.

n = 5 (N = 32)

n = 6 (N = 64)

n = 7 (N = 128)

10.3k 12.2k 530k 25 103 25 1082 245 404 783 809

9963 11.0k 474k 36 78 36 660 237 337 609 734

6421 6356 269k 16 60 16 720 147 234 478 423

32.6k 52.0k — 36 156 36 1273 627 1779 3245 3465

30.5k 46.3k — 55 122 55 1023 605 1540 2297 3086

18.5k 26.8k — 25 96 25 806 418 939 1920 1783

103k 216k — 41 181 41 847 1521 7265 14.3k 14.3k

93.4k 191k — 66 150 66 676 1469 6320 10.0k 12.7k

53.3k 110k — 30 116 30 598 1089 3754 7738 7319

404k 880k — 46 230 46 1003 3567 29.4k 60.4k 58.4k

362k 773k — 77 196 77 714 3453 25.7k 42.0k 51.2k

205k 445k — 35 148 35 748 2688 15.0k 31.1k 29.7k

TABLE III: Terminal relative ℓ2 error for all kernels at n=4 (N =16) across three initial conditions and three backends (statevector “SV”, ideal shot-based “Id”, noisy Aer “Ns”). “—” marks the statevector-only kernels (VQLS, Schrödingerisation, QLS-Fourier; §V-B, §V-D, §V-A). Complements Table IV (n=7). Pulse Kernel

Gaussian Ns

Bimodal Ns

Qubits

HHL 1.36e−1 1.29e−1 7.15e−1 1.64e−2 2.14e−2 5.90e−1 5.10e−2 6.17e−2 3.54e−1 n+m+1 QSVT 8.90e−2 8.90e−2 8.90e−2 3.79e−2 3.79e−2 3.97e−2 5.50e−2 5.50e−2 5.54e−2 n+1 QLS-Fourier 8.04e−1 — — 1.69e−1 — — 2.32e−1 — — n+ log L+1 VQLS 2.77e−2 — — 7.65e−3 — — 1.10e−2 — — n+1 QITE 6.13e−3 1.72e−2 2.74e−2 6.40e−4 7.53e−3 2.11e−2 9.05e−4 4.75e−3 2.05e−2 n var-QITE 1.64e−1 2.39e−1 6.77e−1 3.57e−2 8.75e−2 6.72e−1 5.31e−2 1.21e−1 7.42e−1 n AVQDS 7.08e−1 7.09e−1 7.49e−1 1.68e−1 1.68e−1 2.68e−1 2.31e−1 2.31e−1 3.06e−1 n Hamiltonian sim. 4.60e−1 4.59e−1 5.49e−1 1.52e−1 1.53e−1 2.66e−1 2.13e−1 2.13e−1 2.98e−1 n Schade-Hamiltonian < 10−15 1.21e−2 9.02e−1 < 10−15 1.04e−2 9.12e−1 < 10−15 8.70e−3 8.97e−1 n+1 Schrödingerisation 4.23e−6 — — 1.99e−6 — — 2.59e−6 — — n+np Spectral (QSM) < 10−15 6.78e−3 5.89e−1 < 10−15 9.52e−3 4.68e−1 < 10−15 8.00e−3 3.14e−1 n+1

TABLE IV: Consolidated comparison of all quantum kernels at n=7 (N =128), across three initial conditions (pulse, Gaussian, bimodal) and three backends (statevector “SV”, ideal shot-based “Id”, noisy Aer “Ns”). Values are terminal relative ℓ2 error. “—” marks noisy HHL (fails Qiskit 2Q synthesis at n≥5) and the SV-only kernels noted above. Pulse Kernel

Gaussian Ns

Bimodal Ns

Qubits

HHL 2.13e−1 7.39e−1 — 9.44e−2 7.93e−1 — 1.37e−1 8.93e−1 — n+m+1 QSVT 1.87e−1 1.87e−1 1.87e−1 5.52e−2 5.52e−2 2.53e−2 8.04e−2 8.04e−2 6.77e−2 n+1 QLS-Fourier 1.00e+0 — — 1.75e−1 — — 2.41e−1 — — n+ log L+1 VQLS 7.55e−1 — — 9.13e−1 — — 8.49e−1 — — n+1 QITE 1.68e−1 1.60e−1 7.87e−1 1.33e−2 6.10e−2 7.58e−1 1.71e−2 6.75e−2 7.55e−1 n var-QITE 3.67e−1 9.48e−1 1.00e+0 5.70e−1 5.88e−1 1.00e+0 2.37e−1 5.75e−1 1.00e+0 n AVQDS 1.00e+0 1.00e+0 1.00e+0 1.75e−1 1.73e−1 1.00e+0 2.41e−1 2.40e−1 1.00e+0 n Hamiltonian sim. 8.06e−1 8.06e−1 1.00e+0 1.75e−1 1.73e−1 1.00e+0 2.41e−1 2.50e−1 1.00e+0 n Schade-Hamiltonian < 10−15 2.25e−1 1.00e+0 < 10−15 2.11e−1 1.00e+0 < 10−15 8.43e−2 1.00e+0 n+1 Schrödingerisation 5.64e−4 — — 2.54e−4 — — 3.37e−4 — — n+np Spectral (QSM) < 10−15 2.91e−2 7.79e−1 < 10−15 2.10e−2 6.14e−1 < 10−15 2.00e−2 6.44e−1 n+1

Measurement Cost Quantification (Theoretical) shots N/( 2 psuccess), = 0.01

Total gate count scaling QSM HHL VQLS QITE Hamiltonian simulation AVQDS var-QITE QSVT SZ (Schrödingerization) QLS-Fourier Schade Hamiltonian

Total gate count

105

1010

Estimated shots for 1% relative L2 error

106

104

103

System register qubits n (N = 2^n interior pts)

Fig. 4: Gate-count breakdown (log-y) vs. system register size n (N =2n interior points). Solid circles: 2-qubit (CX) gates.

var-QITE QSVT QLS-Fourier Schade Hamiltonian

109

108

107

106

105

102

QSM HHL QITE Hamiltonian simulation AVQDS

N=16 (n4)

N=32 (n5)

N=64 (n6)

Problem size N (interior grid points)

N=128 (n7)

Fig. 5: Estimated shot budget (log-y) to reach ϵ=10−2 relative ℓ2 error. C. Norm-mismatch ablation

and are limited by expressibility and classical inner loops instead (§V-C).

The smooth-IC error band shared by Hamiltonian simulation, AVQDS, and QLS-Fourier (Table IV, SV Gaussian/bimodal) is algorithmically suspicious: these kernels discretise time, parameterise the ansatz, and invert the linear system in entirely different ways. To test whether the shared

plateau is algorithmic or a readout-pipeline artefact, we decompose each kernel’s terminal error into two contributions on the statevector backend: (i) the norm-only error εN = ∥uref ∥− ∥û∥ /∥uref ∥, the error that remains if the kernel’s direction is correct but its amplitude is mis-scaled; and (ii) the direction residual after renormalising û to the reference norm. TABLE V: Norm-mismatch ablation at n=7 for the three kernels sharing the smooth-IC plateau: εN is the norm-only error, “frac.” is εN /εℓ2 —total on pulse, 14 – 13 of the residual on smooth ICs, identical across three unrelated kernels. Pulse

Gaussian

Bimodal

Kernel

εN

frac.

εN

frac.

εN

frac.

Ham. sim. AVQDS QLS-Fourier

0.755 1.000 1.000

94% 100% 100%

0.041 0.040 0.041

23% 23% 23%

0.071 0.069 0.069

29% 29% 29%

Across all three kernels and all three ICs, εN agrees to three significant figures (Table V); under the pulse initial condition, the entire ℓ2 error is norm drift and a norm-corrected reconstruction recovers a directionally correct field. This identifies the shared plateau as a reconstruction-normalisation failure of the common output map |u⟩=u/∥u∥ —when the kernel does not independently track ∥u(t)∥, the normalised amplitude is indistinguishable from a mis-scaled field at readout. Kernels that do track the norm (QSM, Schade-Hamiltonian) sit at machine precision on the same discretisation. This is a benchmark-methodology result in its own right: full-field accuracy reported under a single reconstruction rule can conflate algorithmic error with readout geometry, and separating them changes which kernels are “competitive.” D. Observable-readout advantage pathway The full-field readout penalty of the output model (§IV-D) is already tomography-scale at modest N . Compact observables—total thermal energy E(t)=⟨ψ(t)|H|ψ(t)⟩, boundary-flux proxies, and the weight of a selected Fourier mode |⟨k|ψ⟩|2 —instead require a constant number of measurement settings independent of the field dimension: a Hadamard or Pauli-estimator for energy, a DST-basis measurement for mode weight, and an ancilla/work-qubit estimator for flux-like quantities. We probe this on QSM and Schade-Hamiltonian at n=7, Gaussian: the shots required to reach ϵ=10−2 on E(t) are O(104 ) per timestep, versus O(106 ) for the full field (N =128 amplitudes at ϵ=10−2 each). The observed gap is not a property of one kernel; it is the readout dimension. This aligns with the standard asymptotic argument [19] and identifies concrete downstream tasks—thermal-energy histories in engineering thermal-management, dominant low-mode tracking in model reduction, or boundary-flux monitoring— where the kernels benchmarked here could be useful before full-field quantum PDE solvers are competitive. VIII. D ISCUSSION We interpret §VII in answer to the question in §I. Subsections cover accuracy and mechanisms (§VIII-A), backend

posture (§VIII-B), cost–accuracy structure (§VIII-C), advantage conditions (§VIII-D), and a selection guide (§VIII-E). A. Accuracy, initial conditions, scaling Two paradigms reach machine precision on SV across all ICs and n (Tables III and IV): QSM and Schade-Hamiltonian (§VII). Schrödingerisation ranks next on SV (∼ 10−4 ); QITE is the strongest non-transform kernel on smooth ICs and degrades on high-frequency content as n grows. VQLS is competitive only at n≤5 (capped-depth ansatz and truncated Pauli expansion, §V-B); var-QITE tracks below QITE on a smaller variational family. In Paradigm A, HHL is capped by its clock-qubit QPE resolution and QSVT by its polynomial degree (§V-A)—a deliberate depth–accuracy tradeoff, not an algorithmic ceiling. The smooth-IC plateau shared by AVQDS, Hamiltonian simulation, and QLS-Fourier is not algorithmic but reconstruction-driven (§VII-C)—using several ICs is what exposed this. Refining the grid inflates compiled resources (Table II, Fig. 4); VQLS/var-QITE hit expressivity walls; HHL/QSVT accumulate coherent depth; QLS-Fourier reaches a native decomposition only at n=4 in our pipeline. B. Backends, readout, and hardware posture The three backends separate algorithmic + reconstruction error (SV), statistical error (Id), and a simple device model (Ns): Id compresses the SV ordering toward the shot floor and inflates postselection-heavy methods; Ns raises the floor for the transpilable subset (§VII); real hardware reuses compute_step but is out of scope. In hardware posture, VQLS is the most NISQ-like linear-system representative (though SV-only by shot cost), Paradigm-C kernels pair shallow circuits with classically dominant inner loops (§V-C), HHL/QSVT are structurally fault-tolerant, and QSM/SchadeHamiltonian sit between. C. Cost–accuracy structure and Pareto reading Practical kernel choice is multi-objective over three resource axes: native 2Q count (Table II), shot budget (Fig. 5), and classical hybrid-loop complexity. On SV at n=7, QSM and Schade-Hamiltonian Pareto-dominate: machine precision at 2Q counts an order of magnitude below HHL/QSVT, with no post-selection penalty and no per-step classical solve beyond the cached eigendecomposition. The AVQDS / Hamiltoniansimulation / QLS-Fourier plateau is reconstruction-bound (§VII-C)—the Pareto hull is defined over correctly normalised fields. D. Conditions for application-level quantum benefit Our experiments do not demonstrate end-to-end quantum advantage at the tested grid sizes; instead, they identify where such a benefit could plausibly enter. The full-field readout cost of §VII-D already dominates the benchmark at N =128 and competes against a classical spectral solver with O(N log N ) structure, so we do not claim a full-field advantage. The credible target is narrower and more applicationoriented: compact-observable readout (total energy, selected

mode weights, boundary-flux proxies) keeps the number of measurement settings independent of N and avoids reconstructing the entire amplitude vector. The relevant question then becomes whether a kernel can prepare the evolved state coherently enough for the target observable before circuit depth, post-selection probability, or hybrid classical overhead—up to O(4n ) per step for full-support Paradigm C variants—dominates. This output-sensitive framing is the practical advantage condition exposed by the benchmark. E. Practical Kernel-Selection Guide Machine-precision reference: QSM or SchadeHamiltonian. Strongest non-transform kernel on smooth data: QITE. Near-term variational study at shallow depth: VQLS or var-QITE, with accuracy degrading at larger n. In Paradigm A, QSVT achieves lower error than HHL but at higher compiled cost; HHL remains a foundational baseline. AVQDS is most useful when adaptive ansatz growth itself is the object of study. For compact-observable workloads, QSM, Schade-Hamiltonian, or QSVT with amplitude-estimation readout is the Pareto choice (§VII-D). IX. T HREATS TO VALIDITY AND F UTURE W ORK The main threats to validity are deliberately exposed rather than hidden. First, the statevector backend isolates algorithmic error but is not a claim of compiled hardware execution. Second, the noisy backend uses a portable depolarising/readout model and omits device coupling maps and T1 /T2 relaxation, so it should be read as a controlled stress test, not as a prediction for a named QPU. The omitted features act in known directions: T1 /T2 relaxation penalises depth (HHL, QSVT, and the large-n spectral/dilation circuits degrade first), SWAP routing moves noise crossovers to smaller n, and correlated errors bias the post-selected kernels’ accepted ensembles. Richer noise models would therefore be expected to shift the failure thresholds of Table IV toward smaller grids; whether the qualitative ranking survives device-calibrated noise remains untested here. Third, the spectral and dilation reference kernels exploit Laplacian-aligned structure; this is the right baseline for the heat equation but may not transfer unchanged to non-normal, nonlinear, or geometry-complex PDEs. Finally, full-field readout remains the dominant bottleneck unless the downstream task is observable-based. We note that real-QPU execution is left to future work because hardware queue variability, calibration drift, and backend-specific transpilation constraints would confound the controlled cross-kernel comparison targeted in this study. Extensions that fit the present harness: higher-dimensional, nonlinear, or non-self-adjoint PDEs with hyperbolic–parabolic terms; sustained real-QPU runs with device-calibrated noise models, which would also let the hybrid-cost argument of §VIII-D be closed with measured rather than analytical numbers; per-kernel hyperparameter sweeps (ansatz depth, QSVT degree, HHL clock qubits) and preconditioners for HHL/QLSFourier; alternative decompositions that lift HHL’s twoqubit synthesis failure at n≥5; and alternative reconstruction

maps (norm-tracking, amplitude-estimation, or observableonly readouts) that generalise §VII-C–§VII-D. A further extension is tensor networks: 1-D diffusion generates little entanglement, so matrix-product-state/operator time-stepping is the natural quantum-inspired classical competitor [13], [14]— consistent with §VIII-D, where no full-field advantage is claimed—and tensor-network simulator backends would let the same kernels be profiled beyond the statevector memory wall. X. C ONCLUSION We presented a controlled cross-paradigm benchmark of eleven quantum kernels on the 1-D heat equation with three research contributions beyond kernel ranking: a single-harness comparison across all five paradigm classes under one reconstruction rule and a staged backend ladder (C1); a normmismatch ablation showing that the smooth-IC plateau shared by Hamiltonian simulation, AVQDS, and QLS-Fourier is driven by the output map, not the unitary core (C2); and an observable-readout analysis showing that compact functionals avoid the tomography-scale bottleneck that dominates full-field recovery (C3). Together these convert the usual “which algorithm wins?” question into a structured one— which kernel, under which output map, for which observable, at which grid size—and supply reproducible measurements for each axis. No single method minimises accuracy loss and all resource terms together; QSM and SchadeHamiltonian are Pareto-optimal reference kernels, QITE is strongest for moderate-depth smooth-IC full-field work, and compact-observable extraction is the clearest application path. The open-source harness (C4) reproduces all figures and tables and admits new kernels through a single interface method. DATA AND C ODE AVAILABILITY The complete source code, per-configuration JSON files, raw simulation outputs (CSV, timing, circuit sidecars), and the analysis scripts that produce every figure and table in this paper are released as an open-source repository. 1 The kernel interface is documented so that new quantum kernels can be added and evaluated under the same benchmark protocol by implementing the BaseKernel.compute_step interface. ACKNOWLEDGMENT The authors gratefully acknowledge Brightskies Technologies for supporting this research, and in particular thank Ibrahim Elghotmy ([email protected]), Omar Marzouk ([email protected]), Amr Nasr ([email protected]), and Dr. Khaled Elamrawi ([email protected]) for their valuable guidance and feedback throughout this work. The authors used Anthropic’s Claude only for language polishing, grammar checking, and editorial copy-editing. No AI system was used to generate scientific claims, numerical data, experiments, figures, tables, or citations. 1 Repository URL: github.com/brightskiesinc/HelloQuantum.

R EFERENCES [1] L. C. Evans, Partial Differential Equations, 2nd ed. Providence, RI, USA: American Mathematical Society, 2010. [2] R. P. Feynman, “Simulating physics with computers,” Int. J. Theor. Phys., vol. 21, no. 6/7, pp. 467–488, 1982. [3] A. W. Harrow, A. Hassidim, and S. Lloyd, “Quantum algorithm for linear systems of equations,” Phys. Rev. Lett., vol. 103, no. 15, Art. no. 150502, 2009. [4] D. W. Berry, “High-order quantum algorithm for solving linear differential equations,” J. Phys. A: Math. Theor., vol. 47, no. 10, Art. no. 105301, 2014. [5] A. M. Childs, J.-P. Liu, and A. Ostrander, “High-precision quantum algorithms for partial differential equations,” Quantum, vol. 5, Art. no. 574, 2021. [6] P. C. S. Costa, S. Jordan, and A. Ostrander, “Quantum algorithm for simulating the wave equation,” Phys. Rev. A, vol. 99, no. 1, Art. no. 012323, 2019. [7] N. Linden, A. Montanaro, and C. Shao, “Quantum vs. classical algorithms for solving the heat equation,” Commun. Math. Phys., vol. 395, pp. 601–641, 2022. [8] M. Cerezo et al., “Variational quantum algorithms,” Nature Rev. Phys., vol. 3, no. 9, pp. 625–644, 2021. [9] J. Preskill, “Quantum computing in the NISQ era and beyond,” Quantum, vol. 2, Art. no. 79, 2018. [10] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, “Variational quantum linear solver,” Quantum, vol. 7, Art. no. 1188, 2023. [11] C. J. Trahan, M. Loveland, N. Davis, and E. Ellison, “A variational quantum linear solver application to discrete finite-element methods,” Entropy, vol. 25, no. 4, Art. no. 580, 2023. [12] H.-L. Liu, Y.-S. Wu, L.-C. Wan, S.-J. Pan, S.-J. Qin, F. Gao, and Q.-Y. Wen, “Variational quantum algorithm for the Poisson equation,” Phys. Rev. A, vol. 104, no. 2, Art. no. 022418, 2021. [13] N. Gourianov, M. Lubasch, S. Dolgov, Q. Y. van den Berg, H. Babaee, P. Givi, M. Kiffner, and D. Jaksch, “A quantum-inspired approach to exploit turbulence structures,” Nature Comput. Sci., vol. 2, no. 1, pp. 30– 37, 2022. [14] M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, “Variational quantum algorithms for nonlinear problems,” Phys. Rev. A, vol. 101, no. 1, Art. no. 010301(R), 2020. [15] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, “Quantum singular value transformation and beyond,” in Proc. 51st Annu. ACM Symp. Theory Comput. (STOC), 2019, pp. 193–204. [16] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, “Grand unification of quantum algorithms,” PRX Quantum, vol. 2, Art. no. 040203, 2021. [17] G. H. Low and I. L. Chuang, “Optimal Hamiltonian simulation by quantum signal processing,” Phys. Rev. Lett., vol. 118, Art. no. 010501, 2017. [18] M. Motta et al., “Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution,” Nature Phys., vol. 16, pp. 205–210, 2020. [19] S. Aaronson, “Read the fine print,” Nature Phys., vol. 11, no. 4, pp. 291– 293, 2015. [20] D. S. Abrams and S. Lloyd, “Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors,” Phys. Rev. Lett., vol. 83, no. 24, pp. 5162–5165, 1999. [21] S. Lloyd, “Universal quantum simulators,” Science, vol. 273, no. 5278, pp. 1073–1078, 1996. [22] M. Suzuki, “Fractal decomposition of exponential operators,” Phys. Lett. A, vol. 146, no. 6, pp. 319–323, 1990. [23] V. V. Shende, S. S. Bullock, and I. L. Markov, “Synthesis of quantumlogic circuits,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst., vol. 25, no. 6, pp. 1000–1010, 2006. [24] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, “Preconditioned quantum linear system algorithm,” Phys. Rev. Lett., vol. 110, no. 25, Art. no. 250504, 2013. [25] A. Ambainis, “Variable time amplitude amplification and quantum algorithms for linear algebra problems,” in Proc. 29th Int. Symp. Theor. Aspects Comput. Sci. (STACS), 2012, pp. 636–647. [26] J. Haah, A. W. Harrow, Z. Ji, X. Wu, and N. Yu, “Sample-optimal tomography of quantum states,” IEEE Trans. Inf. Theory, vol. 63, no. 9, pp. 5628–5641, 2017.

[27] A. Klappenecker and M. Rötteler, “Discrete cosine transforms on quantum computers,” in Proc. Int. Symp. Image Signal Process. Anal. (ISPA), 2001, pp. 464–468. [28] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, 10th anniversary ed. Cambridge, U.K.: Cambridge Univ. Press, 2010. [29] H. Patil, Y. Wang, and P. S. Krstić, “Variational quantum linear solver with a dynamic ansatz,” Phys. Rev. A, vol. 105, no. 1, Art. no. 012423, 2022. [30] Qiskit contributors, “Qiskit: An open-source framework for quantum computing,” 2024. [Online]. Available: https://github.com/Qiskit/qiskit [31] M. E. S. Morales, L. Pira, P. Schleich, K. Koor, P. C. S. Costa, D. An, A. Aspuru-Guzik, L. Lin, P. Rebentrost, and D. W. Berry, “Quantum linear system solvers: A survey of algorithms and applications,” arXiv preprint arXiv:2411.02522, 2024. [32] A. H. Salehi Shayegan, “Quantum linear solvers for scientific computing: A comparison of VQLS, HHL and quantum annealing on timefractional diffusion problems,” Sci. Rep., vol. 16, Art. no. 10278, 2026, doi:10.1038/s41598-026-40910-y. [33] A.-I. Lefterovici, M. Perk, D. Ramacciotti, A. F. Rotundo, S. E. Skelton, and M. Steinbach, “Beyond asymptotic scaling: Comparing functional quantum linear solvers,” IEEE Trans. Quantum Eng., vol. 7, Art. no. 3100118, 2026, doi:10.1109/TQE.2026.3674210. [34] T. Kharazi, A. M. Alkadri, J.-P. Liu, K. K. Mandadapu, and K. B. Whaley, “Explicit block encodings of boundary value problems for many-body elliptic operators,” Quantum, vol. 9, Art. no. 1764, 2025. [35] S. Kumar and C. M. Wilmott, “Generalising quantum imaginary time evolution to solve linear partial differential equations,” Sci. Rep., vol. 14, Art. no. 20156, 2024, doi:10.1038/s41598-024-70423-5. [36] A. B. Özgüler, “Performance evaluation of variational quantum eigensolver and quantum dynamics algorithms on the advection–diffusion equation,” arXiv preprint arXiv:2503.24045, 2025. [37] F. Tennie, S. Laizet, S. Lloyd, and L. Magri, “Quantum computing for nonlinear differential equations and turbulence,” Nature Rev. Phys., vol. 7, pp. 220–230, 2025. [38] A. M. Dalzell et al., “Quantum algorithms: A survey of applications and end-to-end complexities,” arXiv preprint arXiv:2310.03011, 2023. [39] E. Febrianto, Y. Wang, B. Liu, M. Ortiz, and F. Cirak, “A quantum spectral method for non-periodic boundary value problems,” arXiv preprint arXiv:2511.11494, 2025. [40] S. Jin, N. Liu, and Y. Yu, “Quantum circuits for the heat equation with physical boundary conditions via Schrödingerisation,” arXiv preprint arXiv:2407.15895, 2024. [41] S. Jin, N. Liu, and Y. Yu, “Quantum simulation of partial differential equations via Schrödingerization,” Phys. Rev. Lett., vol. 133, no. 23, Art. no. 230602, 2024. [42] S. Tseng, A. Chouhan, and D. Cupidon, “Comparative analysis on two quantum algorithms for solving the heat equation,” arXiv preprint arXiv:2510.04511, 2025. [43] A. M. Childs, R. Kothari, and R. D. Somma, “Quantum algorithm for systems of linear equations with exponentially improved dependence on precision,” SIAM J. Comput., vol. 46, no. 6, pp. 1920–1950, 2017. [44] H. Alipanah, F. Zhang, Y. Yao, R. Thompson, N. Nguyen, J. Liu, P. Givi, B. J. McDermott, and J. J. Mendoza-Arenas, “Quantum dynamics simulation of the advection–diffusion equation,” Phys. Rev. Research, vol. 7, Art. no. 043318, 2025. [45] M. Schade, C. Bösch, V. Hapla, and A. Fichtner, “A quantum computing concept for 1-D elastic wave simulation with exponential speedup,” Geophys. J. Int., vol. 238, no. 1, pp. 321–333, Jul. 2024, doi:10.1093/gji/ggae160. [46] S. Jin, N. Liu, and Y. Yu, “Quantum simulation of partial differential equations: Applications and detailed analysis,” Phys. Rev. A, vol. 108, no. 3, Art. no. 032603, 2023. [47] A. D. McLachlan, “A variational solution of the time-dependent Schrödinger equation,” Mol. Phys., vol. 8, no. 1, pp. 39–44, 1964. [48] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, “Theory of variational quantum simulation,” Quantum, vol. 3, Art. no. 191, 2019. [49] C. Lu, P. Rao, M. G. Meena, and K. C. Gottiparthi, “Distributed Variational Quantum Linear Solver,” arXiv preprint arXiv:2604.14435, 2026.

Record · ID 366328 · SHA-256 cfb13e3bc7405e00
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.