Two Quantum Algorithms for Nonlinear Reaction-Diffusion Equation using Chebyshev Approximation Method

TL;DR

The paper introduces two quantum algorithms that solve nonlinear reaction–diffusion dynamics by first embedding the nonlinear ODEs from a spatial discretization (e.g., Fisher–KPP) into a linear Carleman system and then applying Chebyshev-based methods. A key theoretical advance is a set of sufficient conditions for diagonalizing the Carleman matrix, plus an iterative procedure to construct the diagonalizing transform, enabling global Chebyshev approximations rather than local Taylor-based schemes. Method I (matrix exponentiation) and Method II (quantum spectral method) achieve gate complexities that scale polylogarithmically in the Carleman matrix size and polynomially in time $T$ and precision $1/ obreak\\varepsilon$, with performance comparable to the truncated Taylor-series approach in leading terms. The work also quantifies the impact of scaling, no-resonance constraints, and the norm/dissipation structure on the complexity, and discusses practical shortcomings such as the absence of a universal bound on the diagonalizing condition number. Overall, this Chebyshev-based framework broadens quantum options for nonlinear PDEs by leveraging Carleman linearization and structured matrix techniques to enable efficient quantum simulation of reaction–diffusion dynamics.

Abstract

We present two new quantum algorithms for reaction-diffusion equations that employ the truncated Chebyshev polynomial approximation. This method is employed to numerically solve the ordinary differential equation emerging from the linearization of the associated nonlinear differential equation. In the first algorithm, we use the matrix exponentiation method (Patel et al., 2018), while in the second algorithm, we repurpose the quantum spectral method (Childs et al., 2020). Our main technical contribution is to derive the sufficient conditions for the diagonalization of the Carleman embedding matrix, which is indispensable for designing both quantum algorithms. We supplement this with an efficient iterative algorithm to diagonalize the Carleman matrix. Our first algorithm has gate complexity of O(d$\cdot$log(d)+T$\cdot$polylog(T/$\varepsilon$)). Here $d$ is the size of the Carleman matrix, $T$ is the simulation time, and $\varepsilon$ is the approximation error. The second algorithm is polynomial in $log(d)$, $T$, and $log(1/\varepsilon)$ - the gate complexity scales as O(polylog(d)$\cdot$T$\cdot$polylog(T/$\varepsilon$)). In terms of $T$ and $\varepsilon$, this is comparable to the speedup gained by the current best known quantum algorithm for this problem, the truncated Taylor series method (Costa et.al., 2025). Our approach has two shortcomings. First, we have not provided an upper bound, in terms of d, on the condition number of the Carleman matrix. Second, the success of the diagonalization is based on a conjecture that a specific trigonometric equation has no integral solution. However, we provide strategies to mitigate these shortcomings in most practical cases.

Figures (16)

• Figure 1: Major steps in producing a quantum state encoding the solution of the $\mathsf{Fisher}$-$\mathsf{KPP}$ equation. We have investigated the possibility and advantage of using the Chebyshev approximation method to solve the emerging linear ODE -(i) $\mathsf{matrix\ exponentiation}$ method, and (ii) $\mathsf{quantum\ spectral}$ method. We overcome the main hurdle to adopting these two quantum ODE solvers — the diagonalization of the Carleman matrix ($A$) — by providing a necessary and sufficient condition for it.
• Figure 2: Three point stencil method for special discretization
• Figure 3: Three key steps to solve the quadratic ODE problem. Existing algorithms have attempted to solve the linear ODE (step 2) in different ways. Liu_An_Fang_Wang_Low_Jordan_2023 and costa2025further have solved it using the Euler and Taylor series methods, respectively. We have explored two different methods based on the Chebyshev series method.
• Figure 4: (An iterative procedure for diagonalization of the Carleman matrix) It begins with the uppermost three block matrices (red) and uses \ref{['lem:diag_P']} to get its diagonalized form. Subsequently, it moves to the next block (blue) for diagonalization. This iterative procedure terminates at the last block $A_N^N$.
• Figure 5: (Structure of similarity transformation matrix for $A$) Although matrix $\textbf{A}$ is a block bi-diagonal matrix, the matrix $V$ and $V^{-1}$ are block upper triangular. The product of these three matrices will give the diagonal matrix $\Lambda$.

Theorems & Definitions (35)

• Theorem 2
• Lemma 1
• Theorem 3
• Definition 1: No-resonance condition
• Definition 2: No-resonance condition
• Theorem 4
• Lemma 2
• Lemma 3
• Theorem 5
• Definition 3: Rescaled Quadratic ODE problem