A Sublinear-Time Quantum Algorithm for High-Dimensional Reaction Rates

TL;DR

The paper tackles the challenge of computing reaction rates in high-dimensional Fokker-Planck dynamics by reframing the problem as estimating propagator overlaps of a self-adjoint Hamiltonian $\mathcal{H}_\beta$ obtained via a similarity transform. It introduces a Gaussian linear combination of Hamiltonian simulations (Gaussian-LCHS) to express the non-unitary propagator $e^{Ht}$ as a sublinear-time quantum sum of unitaries, and pairs this with a novel non-unitary overlap circuit that directly estimates overlaps $\langle P|e^{T\mathcal{H}_\beta}|R\rangle$ without the exponential postselection penalties typical for dissipative dynamics. The authors provide explicit block-encoding constructions for the Fokker-Planck operator, a detailed Gaussian-LCHS implementation with complexity $O\left(\alpha_A\sqrt{t\log(1/\epsilon)}\right)$ queries to the underlying block-encoding, and a multiplexed quantum signal processing framework to realize the propagator efficiently. They establish additive-error estimates for the reactive flux $\nu_{RP}$ and the time-dependent rate $k_{RP}(T)$ with gate counts scaling as $\widetilde{O}\left(\frac{\sqrt{t}\alpha_A}{\epsilon}(\eta d n + 2k d n^2 + (nd)^2)\right)$, plus state-preparation costs; under worst-case classical bounds for non-convex potentials, the quantum approach yields an exponential advantage in particle number $\eta$, a quartic improvement in $\epsilon$, and a quadratic improvement in $t$. While practical classical heuristics may close some gaps, the work demonstrates a rigorous quantum route to advantage in high-dimensional dissipative dynamics relevant to chemistry and materials science.

Abstract

The Fokker-Planck equation models rare events across sciences, but its high-dimensional nature challenges classical computers. Quantum algorithms for such non-unitary dynamics often suffer from exponential {decay in} success probability. We introduce a quantum algorithm that overcomes this for computing reaction rates. Using a sum-of-squares representation, we develop a Gaussian linear combination of Hamiltonian simulations (Gaussian-LCHS) to represent the non-unitary propagator with $O\left(\sqrt{t\|H\|\log(1/ε)}\right)$ queries to its block encoding. Crucially, we pair this with {a} novel technique to directly estimate matrix elements without exponential decay. For $η$ pairwise interacting particles discretized with $N$ plane waves per degree of freedom, we estimate reactive flux to error $ε$ using $\widetilde{O}\left((η^{5/2}\sqrt{tβ}α_V + η^{3/2}\sqrt{t/β}N)/ε\right)$ quantum gates, where $α_V = \max_{r}|V'(r)/r|$. For non-convex potentials, the {sharpest classical} worst-case analytical bounds to simulate the related overdamped Langevin {equation} scale as $O(te^{Ω(η)}/ε^4)$. This {implies} an exponential separation in particle number $η$, a quartic speedup in $ε$, and quadratic speedup in $t$. While specialized classical heuristics may outperform these bounds in practice, this demonstrates a rigorous route toward quantum advantage for high-dimensional dissipative dynamics.

Figures (9)

• Figure 1: Stochastic trajectories from an overdamped Langevin equation starting from region $R$ terminate in the configuration space with probabilities given by the Fokker-Planck equation. The stochastic trajectories terminate at some point $\mathbf{X}(t)$ in the configuration space with probability $P(\mathbf{X},t | \mathbf{x} \in R, 0)$ obtained by solving the Fokker-Planck equation with an initial condition corresponding to the equilibrium distribution over $R$.
• Figure 2: Scaling of the condition number $\kappa$ of the finite-difference discretized BKE generator for the one-dimensional double-well $(V(x) = x^4 - x^2)$ test problem as a function of the number of gridpoints per particle $N$ and inverse temperature $\beta$. Condition numbers are estimated using MATLAB's condest function MATLAB. With a barrier height of $\frac{1}{4}$ separating the local minima in the wells, our numerical estimates are consistent with the asymptotic scaling $\kappa(N,\beta) \sim C\, N^{2}e^{\beta/4}$. However, it is computationally intractable to determine the scaling with particle number $\eta$. Vertical lines in panel (a) indicate rough limits for feasible direct sparse solves on a personal workstation and on current exascale-class systems (e.g., the Frontier supercomputer at Oak Ridge National Laboratory Frontier2025). These limits are based on simple extrapolations of available memory, according to the $\sim N^{\eta}$ growth of the full configuration-space dimension. Even though exascale nodes such as those on Frontier provide on the order of a terabyte of combined CPU and high-bandwidth GPU memory, the exponential growth in state-space size quickly outstrips practical storage requirements for fully discretized many-body operators. Consequently, the "curse of dimensionality" associated with grid-based discretization of the BKE operator is typically avoided in classical computing by employing alternative high-dimensional techniques---such as tensor-network or sampling-based methods---as discussed in Sec. \ref{['sec:intro']}.
• Figure 3: A diagrammatic description of the spectral characterization of metastability. There may be many low-lying eigenstates corresponding to local minima of the potential energy surface. The effective gap $\Delta_{\text{eff}}$ corresponds to the gap between the largest of these metastable states and the next largest eigenvalue of the generator.
• Figure 4: A graphical description of how the Hamiltonian model is obtained, and how the algorithm is constructed. In the column on the left, a stochastic Langevin equation can be related to backward and forward Kolmogorov equations via the Feynman-Kac formula and Kramers-Moyal expansion respectively. By performing a similarity transformation of the generator $\mathcal{F}$ with the square root of the equilibrium density, we obtain an imaginary time Schrodinger equation with Hamiltonian $\mathcal{H}_\beta$. This Hamiltonian can then be expressed as a sum-of-squares using the decomposition given by operators $A_j$ that we define below. Then, dilating into the Hermitian matrix $\mathcal{A}_H$, and performing Hamiltonian simulation with $\mathcal{A}_H$, we can approximate time evolution under $\mathcal{H}_\beta$ via the Gaussian-LCHS formula. Finally, we can apply the non-unitary overlap estimation algorithm to estimate the reactive flux, $\bra{P}e^{\mathcal{H}_\beta t}\ket{R}$.
• Figure 5: Quantum circuit for non-unitary overlap estimation. This circuit uses the standard LCU prep construction, a controlled version of the standard sel routine, and an optional controlled phase oracle for complex and non-positive real coefficients. In conjunction with controlled state preparation, this circuit allows us to approximate $\frac{1}{\alpha_g}\operatorname{Re}\left( \bra{P}e^{t \mathcal{H}_\beta}\ket{R} \right)$.

Theorems & Definitions (47)

• Theorem 1
• Theorem 2
• Lemma 2
• Theorem 3
• Theorem 4
• Lemma 4
• Theorem 5: Complexity to perform Gaussian-LCHS
• proof
• Corollary 1: Gate complexity of Gaussian-LCHS applied to Fokker-Planck equation
• proof