# Efficient quantum algorithm for nonlinear
reaction-diffusion equations and energy estimation

Jin-Peng Liu1,2,3,4, Dong An1, Di Fang5, Jiasu Wang3, Guang Hao Low6, Stephen
Jordan6,7
1 Joint Center for Quantum Information and Computer Science, University of
Maryland, MD
2 Simons Institute for the Theory of Computing, Berkeley, CA
3 Department of Mathematics, University of California, Berkeley, CA
4 Center for Theoretical Physics, Massachusetts Institute of Technology,
Cambridge, MA
5 Department of Mathematics and Duke Quantum Center, Duke University, Durham,
NC
6 Microsoft Quantum, Redmond, WA
7 Google Quantum AI, Santa Barbara, CA

###### Abstract

Nonlinear differential equations exhibit rich phenomena in many fields but are
notoriously challenging to solve. Recently, Liu et al. [1] demonstrated the
first efficient quantum algorithm for dissipative quadratic differential
equations under the condition $R<1$, where $R$ measures the ratio of
nonlinearity to dissipation using the $\ell_{2}$ norm. Here we develop an
efficient quantum algorithm based on [1] for reaction-diffusion equations, a
class of nonlinear partial differential equations (PDEs). To achieve this, we
improve upon the Carleman linearization approach introduced in [1] to obtain a
faster convergence rate under the condition $R_{D}<1$, where $R_{D}$ measures
the ratio of nonlinearity to dissipation using the $\ell_{\infty}$ norm. Since
$R_{D}$ is independent of the number of spatial grid points $n$ while $R$
increases with $n$, the criterion $R_{D}<1$ is significantly milder than $R<1$
for high-dimensional systems and can stay convergent under grid refinement for
approximating PDEs. As applications of our quantum algorithm we consider the
Fisher-KPP and Allen-Cahn equations, which have interpretations in classical
physics. In particular, we show how to estimate the mean square kinetic energy
in the solution by postprocessing the quantum state that encodes it to extract
derivative information.

##  1 Introduction

Nonlinear partial differential equations (PDEs) of reaction-diffusion type are
widespread and have many applications, ranging from biology and ecology to
data science. Exhibiting rich phenomena, reaction-diffusion equations have
been applied to describe biological transport networks such as leaf venations
and blood flow [2, 3, 4, 5, 6, 7, 8, 9, 10, 11], predator-prey dynamics in
interacting populations [12, 13, 14], prediction of brain functions and tumor
growth [15, 16], the formation of the Turing patterns in tissues and organs
[17, 18, 19, 20], dendritic colony growth [21, 22], complex chemical processes
such as combustion [23, 24, 25, 26] and calcium dynamics [27]. Reaction-
diffusion equations have also been applied to data classification [28, 29,
30], and image segmentation and inpainting [31, 32, 33, 34]. In many cases,
the underlying model can be viewed as an energy optimization procedure, with
the reaction-diffusion equations as the gradient flow. Such reaction-diffusion
equations inherit the property that energy decays with time. Moreover, the
maximum principle is satisfied, which states that if the initial and boundary
conditions are bounded by a certain constant, then the entire solution remains
bounded (in the $L^{\infty}$ sense) for all time. When designing numerical
approximation schemes it is often of great interest to maintain these
properties exactly.

Attempts to solve such PDEs on classical computers are hampered by the so-
called curse of dimensionality, in which computational complexity grows
exponentially with spatial dimension [35]. For example, in $d$ dimensions, if
each coordinate is discretized by $n$ grid points, the grid will have
$\Omega(n^{d})$ grid points.

Recent advances in quantum computing offer a fresh approach to the efficient
solution of such high-dimensional problems. Quantum algorithms have been
developed to prepare a quantum state encoding the solution to an
$n^{d}$-dimensional linear system, while in some cases only requiring quantum
circuits of complexity $\operatorname{poly}(d,\log n)$ [36, 37, 38, 39, 40,
41, 42, 43, 44]. Such quantum algorithms have been applied to address high-
dimensional problems governed by linear ODEs [45, 46, 47, 48] and PDEs [49,
50, 51, 52, 53, 54, 55].

It has been a longstanding open problem to understand the capability of
quantum computers to solve nonlinear differential equations. An early work
proposed a quantum algorithm for ODEs that simulates polynomial nonlinearities
by employing multiple copies of the solution [56]. For $n$-dimensional systems
of polynomial ordinary differential equations, this quantum algorithm scales
as $\operatorname{poly}(\log n,1/\epsilon^{T})$. Finding quantum algorithms
with polynomial scaling in $T$ for solving nonlinear differential equations
remained an open problem. Furthermore, complexity-theoretic arguments indicate
that this should not be achievable in the most general case [57, 58, 59], but
rather will require exploiting specific properties of restricted classes of
nonlinear differential equations.

Recently, in [1], a quantum algorithm based on Carlemann linearization [60,
61, 62] was proposed for solving a class of nonlinear differential equations

$\frac{dU}{dt}=F_{1}U+F_{2}U^{\otimes 2}+F_{0}(t).$ (1.1)

Here we assume $F_{1}$ is dissipative, i.e. all the eigenvalues of $F_{1}$ are
negative. We are given an initial condition $U(0)=U_{\mathrm{in}}$, an error
tolerance $\epsilon$, and a time-duration $T$. The efficiency of the algorithm
depends on $R$, defined as

$R=\frac{\|F_{2}\|}{|\lambda_{1}|}\|U_{\mathrm{in}}\|,$ (1.2)

where $\lambda_{1}$ is the largest eigenvalue of $F_{1}$. The quantity $R$ is
used to quantify the relative strength of the nonlinearity and forcing to the
linear dissipation according to the $\ell_{2}$ norm. A more general definition
of $R$ for polynomial differential equations is given in (3.10). Under the
condition $R<1$, the algorithm has complexity
$\mathcal{O}\left(\frac{T^{2}q}{g\epsilon}\operatorname{poly}(\log T,\log
n,\log 1/\epsilon)\right)$, where $q=\|U_{\mathrm{in}}\|/\|U(T)\|$, and
$g=\|U(T)\|$. This quadratic scaling with $T$ was an exponential improvement
over prior quantum algorithms for solving nonlinear differential equations.
The error dependence of quantum Carleman linearization was subsequently
improved from $\operatorname{poly}(1/\epsilon)$ to $\operatorname{poly}(\log
1/\epsilon)$ in [63], by assuming the log-norm of the dissipation matrix is
negative rather than exploiting a diagonalizability condition.

Various quantum algorithms for nonlinear differential equations have also been
investigated recently based on Koopman-von Neumann linearization [64, 65, 66,
67, 68] and the related level set formalism [69]. Others have been proposed
based on the non-Hermitian Hamiltonian approach [64, 66, 70] and the homotopy
perturbation approach [71]. Carleman linearization can be treated as a
particular Koopman-von Neumann linearization, while the non-Hermitian
Hamiltonian approach is inspired by quantum simulation. The quantum algorithm
of [71] combines the homotopy perturbation method with the high-precision
quantum linear ODE solver of [46], achieving complexity that scales linearly
in $T$ and polylogarithmically in $1/\epsilon$. This is shown under the
condition $K=4R<1$, which is stricter than the condition $R<1$ used in [1].
The quantum algorithm of [69] is based on the level set method, which maps a
nonlinear differential equation into a linear differential equation describing
the dynamics of the level sets of the solution to the original nonlinear
differential equation. Given a specific construction for the encoding of
initial data, the quantum algorithm of [69] encodes the level set function
from which physical observables can be estimated corresponding to multiple
initial conditions.

Many of the quantum algorithms proposed for solving differential equations
depend on solving a high dimensional linear system, and their complexities are
thus determined by the condition number of this system. Deriving bounds on
this condition number based on the properties of the original differential
equations is challenging and largely unsolved. Quantum complexity lower bounds
on simulating nonlinear quantum dynamics [59] or classical dynamics [1] show
that this condition number becomes exponential in the worst case.

In this paper, we extend upon the quantum algorithm of [1] by adapting the
Carlemann linearization approach to the context of reaction-diffusion PDEs. We
also show that reaction-diffusion equations are still tractable on quantum
computers even for larger $R$ under the Maximum Principle, ruling out the
worst-case exponential time complexity in [1]. Finally, we conduct several
numerical experiments for Fisher-KPP equations and Allen-Cahn equations to
verify the convergence rate and efficiency of the improved Carleman
linearization.

We compare our improved quantum Carleman linearization algorithm to the
original one [1] in Table 1. Both quantum Carleman linearization methods solve
an $n_{d}$-dimensional system of ordinary differential equations with initial
condition $U_{\mathrm{in}}$, for a given evolution time $T$ and normalized
$\ell_{2}$ error tolerance $\epsilon$. (In the present work, we consider this
set of ODEs as arising from the discretization of a PDE.) The quantum
algorithm in [1] solves nonlinear dissipative differential equations of the
form in (1.1). Our new quantum algorithm solves reaction-diffusion equations,
where $F_{1}$ is Laplacian and $F_{0}=0$, but the $F_{2}U^{\otimes 2}$ term is
instead allowed to be a high-degree polynomial. Thus, the class of problems we
consider here is neither a strict generalization nor a strict special case of
that considered in [1].

Our new quantum algorithm produces a Feynman-Kitaev history state encoding the
full time-evolution of the solution $U(t):t\in[0,T]$, whereas the algorithm of
[1] produces a quantum state encoding the final value $U(T)$. The history
state we produce corresponds to the gradient flow of the energy functional. We
can post-process this state to extract derivative information which can be
interpreted as kinetic energy in classical physical systems modeled by the
nonlinear PDE.

The quantum algorithm of [1] has time-complexity proportional to $q/g$ where
$q=\|U_{\mathrm{in}}\|/\|U(t)\|$ and $g=\|U(T)\|$ is the final norm of the
solution. Our new algorithm instead has complexity proportional to
$\|U_{\mathrm{in}}\|/G$, where $G=\frac{1}{T}\int_{0}^{T}\|U(t)\|$ is the
time-averaged norm of the solution. In some cases, $1/G$ can be much smaller
than $1/g$. For example, the solution $u(t)=e^{-t}$ arises in many homogeneous
dissipative differential equations. In this case $1/g=e^{T}$ and
$1/G=\Omega(T)$.

Polynomial complexity is here shown under the assumption $R_{D}<1$, where
$R_{D}$ is a ratio of nonlinearity to dissipation in $\ell_{\infty}$ norm,
whereas the algorithm of [1] requires $R<1$, where $R$ is a ratio of
nonlinearity to dissipation in $\ell_{2}$ norm. The latter is a stronger
assumption, not well suited for the solution of high-dimensional PDEs because
it grows under grid refinement, whereas $R_{D}$ converges to a constant.
Specifically, in the limit where each spatial dimension is discretized into
$n\to\infty$ steps, the number of lattice sites scales as $n^{d}$ and the
$\ell_{2}$ norm of the discretized solution vector scales as $n^{d/2}$, which
leads to divergent $R$ (see (1.12) and (1.13) for detailed discussion).

The solutions to general nonlinear differential equations can have
exponentially growing norms, which would result in an exponential complexity
for the algorithm introduced here. However, we rule this out for reaction-
diffusion equations by establishing an upper bound on the $\ell_{\infty}$ norm
of the solution independent of $T$ as shown in Theorem 5.1.

We also study the extraction of classical information of practical interest
from the history state produced by our algorithm. First, from the quantum
state, we can directly estimate the mean square amplitude over a specific sub-
domain, which can be understood as the portion of a physical observable on
this sub-domain. Our approach is a direct application of amplitude estimate
technique [72] and can achieve a quadratic speedup in precision over standard
classical Monte Carlo sampling. Secondly, we show how to estimate the portion
of the kinetic energy on a specific sub-domain by developing a quantum
algorithm that can transfer a quantum state with function values to a quantum
state encoding its partial derivatives. This algorithm is based on the
discrete Fourier transform. It only requires $\mathcal{O}(1)$ uses of quantum
Fourier transform (QFT) and input oracle of a diagonal matrix, and can
potentially be of independent interest in other problems such as quantum
optimization algorithms. Our main results are summarized in Table 2. Second,
we briefly discuss the potential advantages of the history state compared to
the final state. In particular, the history state structure allows us to
estimate the time when the system reaches equilibrium and run a pre-diagnosis
procedure to avoid possible exponential overhead brought by the fast decay of
the solution.

Algorithm | Model | Output | Condition | Complexity | Grid refinement
---|---|---|---|---|---
[1, Theorem 1] | dissipative quadratic ODEs | final state |  $R<1$ (3.10) | $T^{2}q/(g\epsilon)$ | $R=\mathcal{O}(n_{d}^{1/2})$
Theorem 4.1 |  polynomial R-D systems (2.16) | history state |  $R_{D}<1$ (3.14) | $T^{2}\|U_{\mathrm{in}}\|^{2N}/(G\epsilon)$ | $R_{D}=\mathcal{O}(1)$
Table 1: Comparison between original and improved results of quantum Carleman linearization (QCL). We consider an $n_{d}$-dimensional system of ODEs, given the initial condition $U_{\mathrm{in}}$, the evolution time $T$ and normalized $\ell_{2}$ error tolerance $\epsilon$. We denote $g$, $q$, and $G$ as $\ell_{2}$ norm of the solution, $\ell_{2}$ relation between initial and final solutions, and time-average $\ell_{2}$ norm of the solution. R-D refers to reaction-diffusion. The condition $R<1$ implies the weaker condition $R_{D}<1$. When $R_{D}<1\leqslant R$, we usually choose a suitably small value of the truncation number of Carleman linearization $N$ to reduce the cost resulting from the prefactor $\|U_{\mathrm{in}}\|^{2N}$. Grid refinement refers to the case that the system of ODEs is discretized from a PDE using an increasing number $n_{d}$ of grid points. Logarithmic factors in the complexities are omitted.  Quantum algorithm | Output | Query complexity
---|---|---
Theorem 6.2 | history state of derivatives | $\mathcal{\widetilde{O}}\left(\sqrt{d}\|\vec{f}\|/\|\vec{\nabla f}\|\right)$
Corollary 6.3 | kinetic energy ratio | $\mathcal{\widetilde{O}}\left(\sqrt{d}\|\vec{f}\|/(\|\vec{\nabla f}\|\epsilon)\right)$
Table 2: Summary of potential applications of a quantum state encoding the
information of a smooth function in its amplitudes. Query complexity is the
number of queries to the circuit preparing the quantum state encoding the
function $f$. Here $d$ is the spatial dimension of the function, $\epsilon$ is
the tolerated level of error. We use $\vec{f}$ and $\vec{\nabla f}$ to denote
the unnormalized vectors of $f$ and $\nabla f$ evaluated at discrete grid
points, respectively.

The paper is organized as follows. Section 2 introduces the background of
reaction-diffusion equations. Section 3 develops the Carleman linearization
with $\ell_{2}$ and $\ell_{\infty}$ convergence analysis. Section 4 presents
the problem model and gives the quantum algorithm with a detailed complexity
analysis. Section 5 establishes lower bound results. Section 6 describes how
our approach could be applied to kinetic energy estimation problems. Finally,
we conclude with a discussion of the results and some possible future
directions in Section 7.

###  1.1 Preliminaries

Here we denote the domain, boundaries, functions, and norms as follows.

We consider a $d$-dimensional hypercube as the spatial domain, denoting as
${\mathscr{D}}\coloneqq[0,1]^{d}$. We denote the spatial and time domain as
${\mathscr{D}}_{T}\coloneqq[0,1]^{d}\times(0,T]$. We also denote
$\partial{\mathscr{D}}$ and $\partial{\mathscr{D}}_{T}$ as boundary domains of
${\mathscr{D}}$ and ${\mathscr{D}}_{T}$, respectively.

We consider a uniform spatial discretization on ${\mathscr{D}}$ and introduce
$n$ discretization nodes for each coordinate. To represent it, we denote
$[{n}]_{0}\coloneqq\\{0,1,\ldots,n-1\\}$ and a set of multi-indices as

$\mathcal{I}\coloneqq[{n}]_{0}^{d}=\Bigl{\\{}l=(l_{1},\ldots,l_{d})~{}\Big{|}~{}l_{j}\in[{n}]_{0}\Bigr{\\}}.$
(1.3)

We then denote the set of uniform nodes as

$\chi\coloneqq\Bigl{\\{}\chi_{l}~{}\Big{|}~{}l\in\mathcal{I}\Bigr{\\}},$ (1.4)

where $\chi_{l}$ maps index $l$ to the discretization node. The exact
expression for $\chi_{l}$ depends on the boundary condition. For periodic
boundary condition, $\chi_{l}$ is defined as

$\chi_{l}\coloneqq\left(\frac{l_{1}}{n},\ldots,\frac{l_{d}}{n}\right),$ (1.5)

while for Dirichlet boundary condition, it is given by

$\chi_{l}\coloneqq\left(\frac{l_{1}+1}{n+1},\ldots,\frac{l_{d}+1}{n+1}\right).$
(1.6)

For convenience, we also introduce the set of boundary indices, which is
defined as

$\mathcal{B}\coloneqq\Bigl{\\{}l=(l_{1},\ldots,l_{d})~{}\Big{|}~{}\chi_{l}\in\partial{\mathscr{D}}\Bigr{\\}}.$
(1.7)

Let $u:{\mathscr{D}}_{T}\to\mathbb{R}$ be the solution to a PDE. We can
discretize $u(x,t)$ on the set of uniform nodes $\chi$ to obtain an
$n_{d}$-dimensional vector $U(t)$, where $n_{d}=n^{d}$. The vector’s entries
$U_{1}(t),U_{2}(t),\ldots,U_{n_{d}}(t)$ are the elements of
$\\{u(\chi_{l},t)|l\in\mathcal{I}\\}$, arranged according to the lexicographic
order on $\mathcal{I}$.

We now discuss our notations for norms. For a vector
$a=[a_{1},a_{2},\ldots,a_{n}]\in{\mathbb{R}}^{n}$, we denote the vector
$\ell_{p}$ norm as

$\|a\|_{p}:=\left({\sum_{k=0}^{n-1}|a_{k}|^{p}}\right)^{1/p}.$ (1.8)

For a matrix $A\in{\mathbb{R}}^{n\times n}$, we denote the operator norm
$\left\lVert\cdot\right\rVert_{p,q}$ induced by the vector $\ell_{p}$ and
$\ell_{q}$ norms as

$\|A\|_{p,q}:=\sup_{x\neq 0}\frac{\left\lVert Ax\right\rVert_{q}}{\left\lVert
x\right\rVert_{p}},\quad\|A\|_{p}:=\|A\|_{p,p}.$ (1.9)

For a continuous scalar function $f:[0,T]\rightarrow{\mathbb{R}}$, we denote
the $L^{\infty}$ norm as

$\|f\|_{\infty}:=\max_{t\in[0,T]}|f(t)|.$ (1.10)

For a continuous scalar function
$u:\overline{{\mathscr{D}}}_{T}\rightarrow{\mathbb{R}}$, for a fixed $t$, the
$L^{p}$ norm of $u(\cdot,t)$ is given by

$\|u(\cdot,t)\|_{L^{p}({\mathscr{D}})}:=\left(\int_{\mathscr{D}}|u(x,t)|^{p}\,dx\right)^{1/p}.$
(1.11)

In particular, when no subscript is used, we mean
$\left\lVert\cdot\right\rVert=\left\lVert\cdot\right\rVert_{2}$ for vector and
matrix norms by default, and
$\left\lVert\cdot\right\rVert=\left\lVert\cdot\right\rVert_{L^{2}}$ for
function norm by default.

For a continuous scalar function $u:[0,1]\times[0,T]\rightarrow{\mathbb{R}}$,
which is discretized in space using uniform interpolation nodes, obtain an
estimate of its $L^{p}$ norm as in a Riemann sum:

$\|u(\cdot,t)\|_{L^{p}({\mathscr{D}})}^{p}=\int_{\mathscr{D}}|u(x,t)|^{p}\,dx\approx\sum_{k=0}^{n-1}\left(\left|u\left(\frac{k}{n},t\right)\right|^{p}\frac{1}{n}\right).$
(1.12)

If we denote
$U(t)=[u(\frac{0}{n},t),u(\frac{1}{n},t),\ldots,u(\frac{n-1}{n},t)]$, the RHS
of (1.12) is $\frac{1}{n}\|U(t)\|_{p}^{p}$. This indicates that

$\|u(\cdot,t)\|_{L^{p}({\mathscr{D}})}\approx\frac{1}{n^{1/p}}\|U(t)\|_{p}.$
(1.13)

Thus, if $u(\cdot,t)$ is a given function in continuous one-dimensional space,
then the $\ell_{p}$ norm of its spatially discretized function vector $U(t)$
increases under grid refinement as $n^{1/p}$. Similarly, for a
$u:\overline{{\mathscr{D}}}_{T}\rightarrow{\mathbb{R}}$ with a general spatial
dimension $d$, the $\ell_{p}$ norm of the spatially discretized function
vector $U(t)$ increases as $n^{d/p}$.

Note that when $p=\infty$,

$\|u(\cdot,t)\|_{L^{\infty}({\mathscr{D}})}\approx\|U(t)\|_{\infty}.$ (1.14)

That means the $\ell_{\infty}$ norm of the spatially discretized function
vector $U(t)$ stays convergent in the continuum limit of $n\to\infty$.

For real functions $f,g:{\mathbb{R}}\rightarrow{\mathbb{R}}$, we write
$f=\mathcal{O}(g)$ if there exists $c>0$, such that $|f(\tau)|\leqslant
c|g(\tau)|$ for all $\tau\in{\mathbb{R}}$. We write $f=\Omega(g)$ if
$g=\mathcal{O}(f)$, and $f=\Theta(g)$ if both $f=\mathcal{O}(g)$ and
$g=\mathcal{O}(f)$. We use $\widetilde{\mathcal{O}}$ to suppress logarithmic
factors in the asymptotic expression, i.e., $f=\widetilde{\mathcal{O}}(g)$ if
$f=\mathcal{O}\Bigl{(}g\operatorname{poly}(\log g)\Bigr{)}$. We write $f=o(g)$
if $\limsup_{\tau\to\infty}\frac{|f(\tau)|}{|g(\tau)|}=0$, where $g$ is
nonzero.

##  2 Problem Settings

In this section, we introduce the class of nonlinear PDEs that we focus on and
discuss its spatial discretization, as well as a priori bounds on its
solutions. We then introduce the problem statement with input and output
settings.

###  2.1 Reaction-diffusion equation

We focus on a class of nonlinear PDEs – the reaction-diffusion equations

$\frac{\partial u}{\partial t}(x,t)=D\Delta u(x,t)+f(u(x,t)),$ (2.1)

where $u$ is a real-valued scalar function at position
$x\in{\mathscr{D}}=[0,1]^{d}$ and time $t\in{\mathbb{R}}^{+}$, $f(u)$ is the
nonlinear term and $D$ is a positive number. Without loss of generality,
denoting ${\mathscr{D}}={\mathscr{D}}_{1}\times{\mathscr{D}}_{2}$, we are
given homogeneous Dirichlet boundary conditions imposed on $d_{1}$ dimensions
($x\in{\mathscr{D}}_{1}$) and periodic boundary conditions imposed on $d_{2}$
dimensions ($x\in{\mathscr{D}}_{2}$, $d_{1}+d_{2}=d$)

$\displaystyle u(x,t)$ $\displaystyle=0,\qquad\qquad\qquad
x\in\partial{\mathscr{D}}_{1},$ (2.2) $\displaystyle u(x,t)$
$\displaystyle=u(x+v,t),\qquad v\in
0^{d_{1}}\times{\mathbb{Z}}^{d_{2}},~{}x\in\partial{\mathscr{D}}_{2}.$ (2.3)

The solution $u(x,t)$ in (2.1) is the $L^{2}$ gradient flow of the free energy
functional

$E(u)=\frac{D}{2}\int|\nabla u|^{2}\mathrm{d}x+\int F(u)\mathrm{d}x$ (2.4)

with a potential satisfying

$\frac{\partial F}{\partial u}(u)=-f(u).$ (2.5)

In other words, $f$ is a field driven by the potential $F$.

In this paper, we focus on the following specific reaction-diffusion equations

$\frac{\partial u}{\partial t}(x,t)=D\Delta u(x,t)+au(x,t)+bu^{M}(x,t),$ (2.6)

with the integer $M\geqslant 2$. Without loss of generality, we assume
$|a|,|b|=o(d)$. The motivation to consider this type of PDEs is two-fold:
first, in physical and biological applications, a nonlinearity of this form is
frequently encountered. For example, in the phase transition model (the so-
called Allen-Cahn equation), $M=3$ [73], while $M=2$ corresponds to the
Fisher-KPP equation [74, 17]. Furthermore, on quantum computers, it is a
reasonable task to construct tensor powers of a quantum state, such as
$u^{\otimes M}$, which exactly corresponds to the polynomial nonlinearity in
(2.6). Although we do not do so in this paper, one might also consider an
input model in which a more general nonlinearity $f$ is specified by an oracle
$U_{f}$.

###  2.2 Spatial Discretization

Our approach to solving reaction-diffusion PDEs (2.6) on quantum computers
starts by performing spatial discretization to reduce to a problem of solving
a system of nonlinear ODEs. Specifically, we apply the central difference
discretization on (2.1) to obtain the $n^{d}$-dimensional polynomial ODE

$\frac{\mathrm{d}U_{j}}{\mathrm{d}t}=D\sum_{k}(\Delta_{h})_{jk}U_{k}+f(U_{j}),\quad
j\in\mathcal{I},$ (2.7)

where $\Delta_{h}$ stands for the central difference of the Laplacian with
homogeneous Dirichlet boundary condition or periodic boundary condition,
defined as

$\Delta_{h}=D_{h}\otimes\underbrace{I\otimes\cdots\otimes I}_{d-1\text{
terms}}+I\otimes D_{h}\otimes\underbrace{I\otimes\cdots\otimes I}_{d-2\text{
terms}}+\cdots+\underbrace{I\otimes\cdots\otimes I}_{d-1\text{ terms}}\otimes
D_{h}.$ (2.8)

Here $D_{h}$ is the one-dimensional discrete Laplacian operator. For
homogeneous Dirichlet boundary conditions, $D_{h}$ is

$D_{h}=D_{h}^{\text{Dir}}:=(n+1)^{2}\left(\begin{array}[]{ccccc}-2&1&&&\\\
1&-2&1&&\\\ &\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ &&&1&-2\\\
\end{array}\right)_{n\times n}.$ (2.9)

We denote the eigenvalues of $D_{h}$ as
$\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}$. Specifically,
$\mu_{1}=-4(n+1)^{2}\sin^{2}\left(\frac{\pi}{2n+2}\right)$. For periodic
boundary conditions, the one-dimensional discrete Laplacian operator is

$D_{h}=D_{h}^{\text{per}}:=n^{2}\left(\begin{array}[]{ccccc}-2&1&&&1\\\
1&-2&1&&\\\ &\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ 1&&&1&-2\\\
\end{array}\right)_{n\times n}.$

In this case, the largest eigenvalue of $D_{h}^{\text{per}}$ is
$\mu_{1}^{\text{per}}=0$.

It is worth pointing out that besides serving as a numerical discretization of
the corresponding PDE, the discrete reaction-diffusion equation is of interest
unto itself. For example, the discrete Allen-Cahn equation has been applied to
unsupervised and semi-supervised graph classification, graph cut minimization,
social network segmentation, and image inpainting [28, 29, 75, 76, 77, 30].

We also introduce bounds on the solution to the discrete reaction-diffusion
equations (2.7), which are used in the proof of later theorems.

######  Lemma 2.1 ($\ell_{\infty}$ A priori Bounds on the Solution).

Assume $f\in C^{\infty}(\mathbb{R})$ and has at least two distinct real-valued
roots. Denote any two distinct roots of $f$ as $\gamma_{1},\gamma_{2}$ with
$\gamma_{1}<\gamma_{2}$. Consider the solution
$U(t)=(U_{j})_{j\in\mathcal{I}}$ to (2.7) with initial condition
$U_{j}(0)=u_{0}(x_{j})$ for all $j\in\mathcal{I}$.
(i) Comparison Principle. If the initial condition satisfies

$\gamma_{1}\leqslant U_{j}(0)\leqslant\gamma_{2},\quad\text{for all
$j\in\mathcal{I}$,}$ (2.10)

and so does the solution on the boundary indices $\mathcal{B}$, then the
solution $U_{j}(t)$ remains bounded, that is,

$\gamma_{1}\leqslant U_{j}(t)\leqslant\gamma_{2},\quad\text{for all $t\geq 0$
and all $j\in\mathcal{I}$.}$ (2.11)

(ii) Maximum Principle. In particular, we denote $\gamma$ as the largest
absolute value of roots of $f$. If the initial condition satisfies

$\left\lVert U(0)\right\rVert_{\infty}\leqslant\gamma,$ (2.12)

then

$\left\lVert U(t)\right\rVert_{\infty}\leqslant\gamma,\quad\text{for all
$t\geq 0$.}$ (2.13)

We present the proof of Lemma 2.1 in Section A. It is worth remarking that
this result asserts that the solution $U(t)$ stays in the invariant set
$[\gamma_{1},\gamma_{2}]$, which is different from the type of estimate where
$\left\lVert U(t)\right\rVert_{\infty}\leqslant\left\lVert
U(0)\right\rVert_{\infty}$. In fact, the solution $\left\lVert
U(t)\right\rVert_{\infty}$ can increase in time as depicted in Figure 1. We
also point out that the formation of this invariant region is precisely due to
the nonlinear terms in $f$, and hence the nonlinear parts of the differential
equation cannot be neglected, even if the solution remains small.

For the case of the specific reaction-diffusion equation (2.6), the roots of
$f(u)=au+bu^{M}=0$ gives the explicit expression
$\gamma=(\frac{|a|}{|b|})^{\frac{1}{M-1}}$.

(a) $f(u)=0.4u-u^{2}$ and $u_{0}(x)=0.1-0.1\cos(4\pi x)$.

(b) $f(u)=0.16u-u^{3}$ and $u_{0}(x)=0.1\sin(4\pi x)$.

Figure 1: Solutions $U=(U_{j}(t))$ to (2.7) with homogeneous Dirichlet
boundary condition for different nonlinear terms $f(u)$ and initial conditions
$U_{j}(0)=u_{0}(x_{j})$, where $D=0.005$.

###  2.3 Problem statement

We are interested in solving high-dimensional reaction-diffusion equations
with quantum computers. Given the initial condition described by a quantum
state, we aim to provide a quantum state description of the solution given the
evolution time $T$.

The main computational problem we consider is as follows.

######  Problem 1.

We consider an initial value problem of an $n_{d}$-dimensional polynomial ODE
on $[0,T]$ as in (2.14)

$\frac{\mathrm{d}{U}}{\mathrm{d}{t}}=F_{1}U+F_{M}U^{\otimes M},\qquad
U(0)=U_{\mathrm{in}},$ (2.14)

Here $U=[U_{1},\ldots,U_{n_{d}}]^{T}\in{\mathbb{R}}^{n_{d}}$, $U^{\otimes
M}=[U_{1}^{M},U_{1}^{M-1}U_{2},\ldots,U_{n_{d}}^{M-1}U_{n_{d}-1},U_{n_{d}}^{M}]^{T}\in{\mathbb{R}}^{n_{d}^{M}}$.
We assume $F_{1}$, $F_{M}$ are $s$-sparse111have at most $s$ nonzero entries
in each column and row., $F_{1}$ is symmetric diagonalizable and eigenvalues
of $F_{1}$ are negative, and $\|F_{M}\|\leqslant|\lambda_{1}|$ by
rescaling222given any nonlinear ODE, we can rescale $u\to\alpha u$ with a
proper $\alpha$ to ensure $\|F_{M}\|\leqslant|\lambda_{1}|$.. We also know an
a priori bound on the solution $\max_{t}\|U(t))\|_{\infty}\leqslant\gamma$ as
in (2.13). We have oracles $O_{F_{1}}$, $O_{F_{M}}$ that provide the locations
and values of the nonzero entries of $F_{1}$, $F_{M}$. We also know
$\|U_{\mathrm{in}}\|=\|U(0)\|$ and have an oracle $O_{x}$ that maps $|00\ldots
0\rangle\in{\mathbb{C}}^{n}$ to a quantum state proportional to
$U_{\mathrm{in}}$. Our goal is to produce a quantum state as a superposition
of the solution at different timesteps

$|\hat{Y}_{\mathrm{evo}}\rangle=\frac{1}{\|\hat{Y}_{\mathrm{evo}}\|}\hat{Y}_{\mathrm{evo}}=\frac{1}{\|\hat{Y}_{\mathrm{evo}}\|}\sum_{k=0}^{m}\hat{y}_{1}(kh)|k\rangle$
(2.15)

with a sufficiently large $m$, where $\hat{y}_{1}(t)\in{\mathbb{R}}^{n_{d}}$
is a vector function that approximates $U(t)$,
$\hat{Y}_{\mathrm{evo}}=\sum_{k=0}^{m}\hat{y}_{1}(kh)|k\rangle$ is a
superposition of $\hat{y}_{1}(t)$ at different timesteps, and
$\|\hat{Y}_{\mathrm{evo}}\|=\sqrt{\sum_{k=0}^{m}\|\hat{y}_{1}(kh)\|^{2}}$ is
the normalization factor.

For the reaction-diffusion equation, we have the specific form of (2.14) as

$\frac{\mathrm{d}{U}}{\mathrm{d}{t}}=(D\Delta_{h}+aI)U+bU^{.M},\qquad
U(0)=U_{\mathrm{in}}.$ (2.16)

The corresponding $F_{1}$ and $F_{M}$ in (2.14) satisfy
$F_{1}=D\Delta_{h}+aI\in{\mathbb{R}}^{n_{d}\times n_{d}}$, and
$F_{M}\in{\mathbb{R}}^{n_{d}\times n_{d}^{M}}$ maps $u^{\otimes M}$ to
$bu^{M}$, with
$U^{.M}=[U_{1}^{M},U_{2}^{M},\ldots,U_{n_{d}}^{M}]^{T}\in{\mathbb{R}}^{n^{M}}$,
and henceforth $F_{1}$ and $F_{M}$ are $s$-sparse with $s=O(1)$. The
representation of the Laplacian matrix $\Delta_{h}$ with mixed boundary
conditions refers to (2.3), with $d_{1}$-dimensional Dirichlet boundary
conditions and $d_{2}$-dimensional periodic boundary conditions
($d_{1}+d_{2}=d$). We require the eigenvalues $\lambda_{j}$ of $F_{1}$ are
negative, i.e., $4Dd_{1}(n+1)^{2}\sin^{2}\left(\frac{\pi}{2n+2}\right)>a$.

The quantum state (2.15) corresponding with $[U(0),U(h),\ldots,U(mh)]$, also
known as the history state in [78], describes a quantum-state evolution for
the gradient flow of the free energy (2.4).

##  3 Carleman linearization

We aim to perform Carleman linearization on discretized nonlinear PDEs (2.14)
and then use quantum linear system solvers to obtain quantum states proportion
to the solutions. In this section, we revisit the Carleman linearization
procedure and then introduce the improved convergence result.

Defining $\widetilde{y}_{j}:=U^{\otimes j}$ for $j=1,2\cdots$, one has

$\frac{d}{dt}\widetilde{y}_{j}=A_{j}^{j}\widetilde{y}_{j}+A_{j+M-1}^{j}\widetilde{y}_{j+M-1},$
(3.1)

where

$\displaystyle A_{j}^{j}=\sum_{\nu=1}^{j}\overset{\text{j
factors}}{\overbrace{\mathbb{I}_{n_{d}\times
n_{d}}\otimes\cdots\otimes\underset{\underset{\nu-\text{th
position}}{\uparrow}}{F_{1}}\otimes\cdots\otimes\mathbb{I}_{n_{d}\times
n_{d}}}},$ (3.2) $\displaystyle A_{j+M-1}^{j}=\sum_{\nu=1}^{j}\overset{\text{j
factors}}{\overbrace{\mathbb{I}_{n_{d}\times
n_{d}}\otimes\cdots\otimes\underset{\underset{\nu-\text{th
position}}{\uparrow}}{F_{M}}\otimes\cdots\otimes\mathbb{I}_{n_{d}\times
n_{d}}}}.$ (3.3)

Therefore, the Carleman linearization procedure gives rise to the following
infinite-dimensional system
$\widetilde{y}^{\prime}(t)=\widetilde{A}\widetilde{y}(t),$ where
$\widetilde{A}$ is an infinite-dimensional block upper-triangular matrix

$\widetilde{A}:=\setcounter{MaxMatrixCols}{11}\begin{pmatrix}A_{1}^{1}&0&\cdots&0&A_{M}^{1}&0&\ldots&0&0&0&\ldots\\\
0&A_{2}^{2}&0&\cdots&0&A_{M+1}^{2}&\ldots&0&0&0&\ldots\\\
0&0&A_{3}^{3}&0&\cdots&0&A_{M+2}^{3}\ldots&0&0&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\\\
\end{pmatrix}.$

It follows from the definition of $A_{k}^{j}$ that the following inequalities
are satisfied.

######  Lemma 3.1.

For all $j\geq 1$,

$\left\lVert A^{j}_{j+M-1}\right\rVert_{2},\left\lVert
A^{j}_{j+M-1}\right\rVert_{\infty}\leqslant j|b|,$ (3.4) $\left\lVert
e^{tA_{j}^{j}}\right\rVert_{2}\leqslant e^{j\lambda_{1}t},$ (3.5)

where $\lambda_{1}\coloneqq Dd_{1}\mu_{1}+a$,
$\mu_{1}=-4(n+1)^{2}\sin^{2}\left(\frac{\pi}{2n+2}\right)$.

We then truncate the above infinite-dimensional system of linear ODEs at order
$N$, thereby obtaining a finite system

$\frac{\mathrm{d}{\hat{y}}}{\mathrm{d}{t}}=A\hat{y},\qquad\hat{y}(0)=\hat{y}_{\mathrm{in}}$
(3.6)

with the upper triangular block structure

$\frac{\mathrm{d}{}}{\mathrm{d}{t}}\begin{pmatrix}\hat{y}_{1}\\\
\hat{y}_{2}\\\ \vdots\\\ \vdots\\\ \vdots\\\ \hat{y}_{N-1}\\\ \hat{y}_{N}\\\
\end{pmatrix}=\begin{pmatrix}A_{1}^{1}&0&\cdots&0&A_{M}^{1}&&\\\
&A_{2}^{2}&\ddots&&\ddots&\ddots&\\\ &&\ddots&\ddots&&\ddots&A_{N}^{N-M+1}\\\
&&&\ddots&\ddots&&0\\\ &&&&\ddots&\ddots&\vdots\\\ &&&&&A_{N-1}^{N-1}&0\\\
&&&&&&A_{N}^{N}\\\ \end{pmatrix}\begin{pmatrix}\hat{y}_{1}\\\ \hat{y}_{2}\\\
\vdots\\\ \vdots\\\ \vdots\\\ \hat{y}_{N-1}\\\ \hat{y}_{N}\\\ \end{pmatrix}.$
(3.7)

Here, $\hat{y}_{j}=U^{\otimes j}\in{\mathbb{R}}^{n_{d}^{j}}$,
$\hat{y}_{\mathrm{in}}=[U_{\mathrm{in}};U_{\mathrm{in}}^{\otimes
2};\ldots;U_{\mathrm{in}}^{\otimes N}]$, and
$A_{j}^{j}\in{\mathbb{R}}^{n_{d}^{j}\times n_{d}^{j}}$,
$A_{j+1}^{j}\in{\mathbb{R}}^{n_{d}^{j}\times n_{d}^{j+1}}$ for $j\in[{N}]$ are
defined as (3.2). Note that $A$ is an $(Ns)$-sparse matrix, where $s$ is the
largest nonzero number of each column and row in $F_{1}$ and $F_{M}$. The
dimension of (3.6) is denoted as

$\mathcal{N}_{d,N}\coloneqq
n_{d}+n_{d}^{2}+\cdots+n_{d}^{N}=\frac{n_{d}^{N+1}-n_{d}}{n_{d}-1}=\mathcal{O}(n^{Nd}).$
(3.8)

Denote the solution to the truncated system as $\hat{y}_{j}$ ($j=0,1,\cdots$)
and define the error resulting from the truncation as

$\eta_{j}(t)\coloneqq\widetilde{y}_{j}(t)-\hat{y}_{j}(t)=U^{\otimes
j}(t)-\hat{y}_{j}(t).$ (3.9)

In particular,
$\eta_{1}(t)=\widetilde{y}_{1}(t)-\hat{y}_{1}(t)=U(t)-\hat{y}_{1}(t)$ is the
error due to the Carleman linearization procedure.

######  Theorem 3.2 ($\ell_{2}$ Convergence of the Carleman Linearization).

For the discrete reaction-diffusion equation (2.16) with mixed boundary
conditions in (2.3), as originally proposed in [1], we define

$\displaystyle R$
$\displaystyle=\frac{\|F_{M}\|}{|\lambda_{1}|}\|U_{\mathrm{in}}\|^{M-1},$
(3.10) $\displaystyle\overline{R}$
$\displaystyle=\frac{\|F_{M}\|}{|\lambda_{1}|}\max_{t}\|U(t)\|^{M-1}.$

Suppose that the largest eigenvalue of $D\Delta_{h}+aI$, denoted by
$\lambda_{1}$, is negative. Then the approximation error of the Carleman
linearization satisfies

$\left\lVert\eta_{j}(t)\right\rVert\leqslant\max_{t}\|U(t)\|^{j}\overline{R}^{\lceil\frac{N-j+1}{M-1}\rceil}\left(1-e^{j\lambda_{1}t}\right)$
(3.11)

for $t\geq 0$. In particular, if $N$ is some integer multiple of $M-1$, then
the error of the solution
$\eta_{1}(t)=y_{1}(t)-\hat{y}_{1}(t)=U(t)-\hat{U}(t)$ satisfies

$\left\lVert\eta_{1}(t)\right\rVert\leqslant\max_{t}\|U(t)\|\overline{R}^{\frac{N}{M-1}}\left(1-e^{\lambda_{1}t}\right)$
(3.12)

for $t\geq 0$. Furthermore, if $R\leq 1$, then $R=\overline{R}$.

The detailed proof is presented in Section D. This theorem for polynomial
differential equations is a straightforward extension of the quadratic case in
[1, Corollary 1] and implies exponential convergence in the order of
truncation $N$ in terms of $\ell_{2}$ norm as long as $R\leq 1$, that is

$\frac{\|F_{M}\|}{|\lambda_{1}|}\|U_{\mathrm{in}}\|^{M-1}=\frac{|b|}{|\lambda_{1}|}\|U_{\mathrm{in}}\|^{M-1}<1,$
(3.13)

as shown in (3.10). However, if $R>1$, we only have $R\leqslant\overline{R}$,
and $\overline{R}$ is the exact convergence radius.

######  Theorem 3.3 ($\ell_{\infty}$ Convergence of the Carleman
Linearization).

For the discrete reaction-diffusion equation (2.16) with mixed boundary
conditions as proposed in (2.3), we define

$\displaystyle R_{D}\coloneqq\frac{|b|}{|\lambda_{1}|}\gamma^{M-1}C(\lambda).$
(3.14)

Here $\gamma=(\frac{|a|}{|b|})^{\frac{1}{M-1}}$, and $\gamma$ upper bounds
$\|U(t)\|_{\infty}$ for all $t\geqslant 0$, as stated in Lemma 2.1; the
constant $C(\lambda)$ has the form

$\displaystyle
C(\lambda)\coloneqq\begin{cases}\frac{|\lambda_{1}|}{a}(e^{\frac{\ln(3)d_{1}}{2(\lambda-\lambda_{1})}a}-1)+\frac{|\lambda_{1}|}{|\lambda|},&a\neq
0,\\\
\frac{\ln(3)d_{1}}{2(\lambda-\lambda_{1})}|\lambda_{1}|+\frac{|\lambda_{1}|}{|\lambda|},&a=0,\end{cases}$
(3.15)

where $\mu_{1}=-4(n+1)^{2}\sin^{2}\left(\frac{\pi}{2n+2}\right)$,
$\lambda_{1}=Dd_{1}\mu_{1}+a<0$, and $\lambda$ is an arbitrary value
satisfying $\lambda_{1}<\lambda<0$. Then the approximation error of the
Carleman linearization satisfies

$\left\lVert\eta_{j}(t)\right\rVert_{\infty}\leqslant\gamma^{j}R_{D}^{\lceil\frac{N-j+1}{M-1}\rceil}$
(3.16)

for $t\geq 0$, with $R_{D}$ defined as (3.14). In particular, if $N$ is some
integer multiple of $M-1$, then the error of the solution
$\eta_{1}(t)=y_{1}(t)-\hat{y}_{1}(t)=U(t)-\hat{U}(t)$ satisfies

$\left\lVert\eta_{1}(t)\right\rVert_{\infty}\leqslant\gamma
R_{D}^{\frac{N}{M-1}}$ (3.17)

for $t\geq 0$.

This theorem implies an alternative exponential convergence in the order of
truncation $N$ in terms of $\ell_{\infty}$ norm as long as $R_{D}<1$, that is,

$\frac{|b|}{|\lambda_{1}|}\gamma^{M-1}C(\lambda)<1,$ (3.18)

as shown in (3.14).

###### Proof.

The truncation error $\eta_{j}$ satisfies the equation

$\eta^{\prime}_{j}(t)=A^{j}_{j}\eta_{j}(t)+A^{j}_{j+M-1}\left(y_{j+M-1}(t)-\hat{y}_{j+M-1}(t)\delta_{j+M-1\leqslant
N}\right),\qquad 1\leqslant j\leqslant N.$ (3.19)

Applying the variation of constants formula [79] to (3.19), one has

$\displaystyle\eta_{j}(t)=\int_{0}^{t}e^{A^{j}_{j}(t-s)}A^{j}_{j+M-1}y_{j+M-1}(s)\,ds,\quad
N-M+2\leqslant j\leqslant N.$ (3.20)

Note that it follows from Lemma 2.1 that

$\left\lVert y_{j+M-1}(s)\right\rVert_{\infty}=\left\lVert
U^{\otimes(j+M-1)}(s)\right\rVert_{\infty}\leqslant\left\lVert
U(s)\right\rVert^{j+M-1}_{\infty}\leqslant\gamma^{j+M-1}.$

Therefore, we have for $N-M+2\leqslant j\leqslant N$,

$\displaystyle\left\lVert\eta_{j}(t)\right\rVert_{\infty}$
$\displaystyle\leqslant\int_{0}^{t}\left\lVert
e^{A^{j}_{j}(t-s)}\right\rVert_{\infty}\left\lVert
A^{j}_{j+M-1}\right\rVert_{\infty}\left\lVert
y_{j+M-1}(s)\right\rVert_{\infty}\,ds$ (3.21) $\displaystyle\leqslant
j|b|\gamma^{j+M-1}\int_{0}^{t}\left\lVert
e^{A^{j}_{j}(t-s)}\right\rVert_{\infty}\,ds$ $\displaystyle\leqslant
j|b|\gamma^{j+M-1}\int_{0}^{t}\|e^{j(t-s)(D\Delta_{h}+a)}\|_{\infty}\,ds.$

For simplicity, we denote

$C_{j}(t)\coloneqq
j|\lambda_{1}|\int_{0}^{t}\|e^{j(t-s)(D\Delta_{h}+a)}\|_{\infty}\,ds$ (3.22)

such that for $N-M+2\leqslant j\leqslant N$,

$\displaystyle\left\lVert\eta_{j}(t)\right\rVert_{\infty}\leqslant\frac{|b|}{|\lambda_{1}|}\gamma^{j+M-1}C_{j}(t).$
(3.23)

Next, for $N-2M+3\leqslant j\leqslant N-M+1$,

$\displaystyle\left\lVert\eta_{j}(t)\right\rVert_{\infty}$
$\displaystyle\leqslant\int_{0}^{t}\left\lVert
e^{A^{j}_{j}(t-s)}\right\rVert_{\infty}\left\lVert
A^{j}_{j+M-1}\right\rVert_{\infty}\left\lVert\eta_{j+M-1}(s)\right\rVert_{\infty}\,ds$
(3.24)
$\displaystyle\leqslant\frac{|b|}{|\lambda_{1}|}\gamma^{j+2M-2}C_{j+M-1}(t)\int_{0}^{t}\left\lVert
e^{A^{j}_{j}(t-s)}\right\rVert_{\infty}\left\lVert
A^{j}_{j+M-1}\right\rVert_{\infty}\,ds$
$\displaystyle\leqslant\left(\frac{|b|}{|\lambda_{1}|}\right)^{2}\gamma^{j+2M-2}C_{j+M-1}(t)C_{j}(t).$

One can continue by mathematical induction for every group of $M-1$ terms and
obtain

$\displaystyle\left\lVert\eta_{j}(t)\right\rVert_{\infty}$
$\displaystyle\leqslant\int_{0}^{t}\left\lVert
e^{A^{j}_{j}(t-s)}\right\rVert_{\infty}\left\lVert
A^{j}_{j+M-1}\right\rVert_{\infty}\left\lVert\eta_{j+M-1}(s)\right\rVert_{\infty}\,ds$
(3.25)
$\displaystyle\leqslant\left(\frac{|b|}{|\lambda_{1}|}\right)^{\lceil\frac{N-j+1}{M-1}\rceil}\gamma^{j+(M-1)\lceil\frac{N-j+1}{M-1}\rceil}\prod_{k=1}^{\lceil\frac{N-j+1}{M-1}\rceil}C_{j+(M-1)(k-1)}(t).$

We now consider an upper bound on $C_{j}(t)$. By computing the integration of
$\|e^{j(t-s)(D\Delta_{h}+a)}\|_{\infty}$ in Lemma B.6, we have

$\displaystyle\int_{0}^{t}\|e^{j(t-s)(D\Delta_{h}+a)}\|_{\infty}\,ds\leqslant\begin{cases}\frac{1}{ja}\left(e^{\frac{\ln(3)d_{1}}{2(\lambda-\lambda_{1})}a}-1\right)+\frac{1}{j|\lambda|},&a\neq
0,\\\
\frac{\ln(3)d_{1}}{2j(\lambda-\lambda_{1})}+\frac{1}{j|\lambda|},&a=0.\end{cases}$
(3.26)

Therefore,

$\displaystyle C_{j}(t)$
$\displaystyle=j|\lambda_{1}|\int_{0}^{t}\|e^{j(t-s)(D\Delta_{h}+a)}\|_{\infty}\,ds\leqslant\begin{cases}\frac{|\lambda_{1}|}{a}\left(e^{\frac{\ln(3)d_{1}}{2(\lambda-\lambda_{1})}a}-1\right)+\frac{|\lambda_{1}|}{|\lambda|},&a\neq
0,\\\
\frac{\ln(3)d_{1}}{2(\lambda-\lambda_{1})}|\lambda_{1}|+\frac{|\lambda_{1}|}{|\lambda|},&a=0,\end{cases}$
(3.27)

where the right-hand side is exactly $C(\lambda)$ in (3.15) and is independent
of both $j$ and $t$.

Finally, substituting $C_{j}(t)$ with its upper bound $C(\lambda)$ in (3.25)
gives

$\displaystyle\left\lVert\eta_{j}(t)\right\rVert_{\infty}\leqslant\left(\frac{|b|}{|\lambda_{1}|}\right)^{\lceil\frac{N-j+1}{M-1}\rceil}\gamma^{j+(M-1)\lceil\frac{N-j+1}{M-1}\rceil}C(\lambda)^{\lceil\frac{N-j+1}{M-1}\rceil}\leqslant\gamma^{j}R_{D}^{\lceil\frac{N-j+1}{M-1}\rceil},$
(3.28)

where we use $R_{D}=\frac{|b|}{|\lambda_{1}|}\gamma^{M-1}C(\lambda)$ as
defined in (3.14). Therefore, a sufficient condition for the convergence of
$\left\lVert\eta_{j}(t)\right\rVert_{\infty}$ in $N$ is $R_{D}<1$.

In practice, we can set $N$ as some integer multiple of $M-1$, and thus

$\left\lVert\eta_{1}\right\rVert_{\infty}\leqslant\gamma
R_{D}^{\frac{N}{M-1}}.$ (3.29)

This completes the proof of the desired result. ∎

Remark. According to Theorem 3.2 and Theorem 3.3, the Carleman linearized
solution $\hat{y}_{1}(t)$ approximates the exact solution $U(kh)$ of the
original reaction-diffusion equations with exponential convergence rate in
terms of the convergence radius $R$ or $R_{D}$.

The quantities $R$ and $R_{D}$ are used to characterize the ratios of reaction
and diffusion strengths in terms of $\ell_{\infty}$ and $\ell_{2}$ norms. Here
we briefly discuss the relationship between these. In particular, we are
interested in the case where the convergence radius $R_{D}$ has an advantage
over $R$. Note that $R_{D}\leq R$ is equivalent to
$C(\lambda)\leq\left(\frac{\|U_{\mathrm{in}}\|}{\gamma}\right)^{M-1}$ and that
$\lambda$ can be an arbitrary number between $\lambda_{1}$ and 0. Hence the
regime of our interest turns out to be

$\min_{\lambda_{1}<\lambda<0}C(\lambda)\leq\left(\frac{\|U_{\mathrm{in}}\|}{\gamma}\right)^{M-1},$
(3.30)

which holds true for a large regime of parameters in the high-dimensional or
finely discretized scenarios because the right-hand side is likely to grow
rapidly in $n$ and $d$ while the left-hand side only has a weak dependence.

Specifically, according to Section E, solving the optimization problem on the
left-hand side of (3.30) helps to obtain a sharper estimate of $R_{D}$. When
$a=0$, the minimum of $C(\lambda)$ has an explicit expression. As for $a\neq
0$, we advise tuning $\lambda$ for a sharper estimate of $R_{D}$ in real
applications, since the optimization problem is hard to solve explicitly.
Nevertheless, in both cases, we can show that there exists an upper bound of
$\min_{\lambda_{1}<\lambda<0}C(\lambda)$ which is also independent of $n$. For
sufficiently large $n$, the quantity $\min_{\lambda_{1}<\lambda<0}C(\lambda)$
is $\mathcal{O}(d)$, where $d$ is the spatial dimension.

For $n_{d}$-dimensional vectors, $\|U_{\mathrm{in}}\|$ can be significantly
larger than $\gamma$ due to the inequality
$\|U_{\mathrm{in}}\|_{\infty}\leqslant\|U_{\mathrm{in}}\|\leqslant\sqrt{n^{d}}\|U_{\mathrm{in}}\|_{\infty}$.
First, when $d$ is fixed, $\|U_{\mathrm{in}}\|$ has a polynomial growth with
$n$ in the worst case and then $R_{D}<R$ for large $n$. Second, when $n$ is
large enough and fixed, $\|U_{\mathrm{in}}\|$ increases exponentially with $d$
in the worst case while $\min_{\lambda_{1}<\lambda<0}C(\lambda)$ grows at most
linearly with $d$. Therefore, $R_{D}$ is smaller than $R$ for large $d$ as
well.

For the case of grid refinement, the above results also show that $R_{D}$
stays bounded in the continuum limit $n_{d}\to\infty$, while $R$ diverges.

###  3.1 Numerical results

In this subsection, we present some numerical results to examine the
effectiveness of Carleman linearization.

Figure 2: Applying Carleman linearization to (2.1) on a classical computer,
where $D=0.2$, $f(u)=0.2u-u^{2}$, the initial condition is
$u(x,0)=0.1\left(1-\cos(2\pi x)\right)$ and homogeneous Dirichlet boundary
conditions are imposed. Left: $l_{\infty}$ norm of the absolute error between
the Carleman solutions at various truncation levels $N$. Right: the
convergence of the corresponding time-maximum error. Figure 3: Applying
Carleman linearization to (2.1) on a classical computer, where $D=0.1$,
$f(u)=0.16u-u^{3}$, the initial condition is $u(x,0)=0.2\sin(2\pi x)$ and
homogeneous Dirichlet boundary conditions are imposed. Left: $l_{\infty}$ norm
of the absolute error between the Carleman solutions at various truncation
levels $N$. Right: the convergence of the corresponding time-maximum error.

In order to demonstrate the convergence of Carleman linearization (Theorem
3.3), we apply our algorithm to (2.1) with different types of nonlinearity
$f(u)$. In the first example, the nonlinear term is $f(u)=0.2u-u^{2}$, the
Fisher-KPP type. We assume $D=0.2$, choose $u(x,0)=0.1\left(1-\cos(2\pi
x)\right)$ as the initial condition, and impose homogeneous Dirichlet boundary
conditions. In our second example, $f(u)=0.16u-u^{3}$, $D=0.1$,
$u(x,0)=0.2\sin(2\pi x)$, and homogeneous Dirichlet boundary conditions are
again used. The numerical results for these two examples are depicted in
Figure 2 and Figure 3, respectively. We see from the error convergence plots
that the absolute error, maximized over $t\in[0,1]$, decreases exponentially
as the truncation level $N$ is increased. As a function of the time $t$, the
absolute error first increases and then decreases exponentially due to the
decay of the exact solution. In particular, in Figure 3, the absolute error
curve depicting the absolute error for $N=1$ agrees with that for $N=2$. That
is because, according to (3.7), only $\hat{y}_{1}$ takes part in the time
evolution of $\hat{y}_{1}$, no matter whether $N=1$ or $2$. A similar argument
holds for the agreement of two curves for $N=3$ and $N=4$.

(a) $D=0.01$, $f(u)=0.25u-u^{3}$ and $u(x,0)=0.08\sin(2\pi x)$.

(b) $D=0.012$, $f(u)=0.0196u-u^{3}$ and $u(x,0)=0.14\sin(2\pi x)$ .

Figure 4: Examples of the application of Carleman linearization to (2.1) with
homogeneous Dirichlet boundary conditions for different parameters and initial
conditions. Left: The eigenvalues of $F_{1}$ are no longer all positive.
Right: $R>1$ while $R_{D}<1$.

Based on our numerical tests, we also find that Carleman linearization works
for more general cases. In Figure 4(a), we relax the requirement that the
eigenvalues of $F_{1}$ are all negative. We test (2.1) with homogeneous
boundary conditions and choose $D=0.01$, $f(u)=0.25u-u^{3}$, and
$u(x,0)=0.08\sin(2\pi x)$ as initial condition. Carleman linearization still
has good numerical performance, as implied by the absolute error plot. Figure
4(b) illustrates the advantage of $R_{D}$ over $R$. In that example, we
consider (2.1) with homogeneous boundary conditions and assume $D=0.012$,
$f(u)=0.0196u-u^{3}$ and $u(x,0)=0.14\sin(2\pi x)$. We discretize the spatial
domain into $15$ sub-intervals, i.e., the value of $n$ is 16. By computation,
$R=1.4924$ and $R_{D}=0.9299$ where we choose $\lambda=\lambda_{1}/2.3$. This
illustrates that error remains well-controlled under $R_{D}<1$, which is a
milder condition than $R<1$.

##  4 Quantum algorithm

We now describe an efficient quantum algorithm for computing the numerical
solution of the linearized ODEs (3.6). Our main algorithmic result is stated
as follows.

######  Theorem 4.1 (Quantum Carleman linearization).

Consider an instance of the quantum ODE problem as defined in Problem 1, with
its $N$-th Carleman linearization as defined in (3.6). We denote a parameter

$G\coloneqq\frac{1}{\sqrt{m+1}}\|\hat{Y}_{\mathrm{evo}}\|=\sqrt{\frac{\sum_{k=0}^{m}\|\hat{y}_{1}(kh)\|^{2}}{m+1}},$
(4.1)

which parameterized the average $\ell_{2}$ norm of the evolution of the $N$-th
Carleman linearized solution $\hat{Y}_{\mathrm{evo}}$. Assume $R_{D}<1$. Then
there exists a quantum algorithm that produces a state that approximates
$\hat{Y}_{\mathrm{evo}}$ in terms of the $\ell_{2}$ normalized error
$\epsilon\leqslant 1/2$ succeeding with probability $\Omega(1)$, with a flag
indicating success. The query complexity (to the oracles
$O_{F_{1}},O_{F_{M}}$, and $O_{x}$) is

$\displaystyle\frac{1}{G^{2}\epsilon}sT^{2}D^{2}d^{2}n^{4}N^{3}\|U_{\mathrm{in}}\|^{2N}\cdot\operatorname{poly}\biggl{(}\log\frac{aDdMnNsT}{G\epsilon}\biggr{)}.$
(4.2)

The gate complexity is larger than its query complexity by logarithmic
factors.

We next describe the quantum algorithm in detail, including ingredients such
as state preparation, quantum linear system algorithm, and measurement. We
conclude the proof of Theorem 4.1 at the end of this section.

Remark. We notice that a prefactor $\|U_{\mathrm{in}}\|^{2N}$ in the
complexity of the $N$-th Carleman linearization. This cost is similar to
prefactor $5^{2k}$ in the complexity of the $k$-th product formula [80]. In
practice, we usually choose a suitably small value of $N$, such as $1,2,3$, to
reduce the cost.

###  4.1 State preparation

We first recall a lemma used in [1] for preparing a quantum state
corresponding to the initial vector
$\hat{y}_{\mathrm{in}}=[U_{\mathrm{in}};U_{\mathrm{in}}^{\otimes
2};\ldots;U_{\mathrm{in}}^{\otimes N}]$, given the value $\|U_{\mathrm{in}}\|$
and the ability to prepare a quantum state proportional to $U_{\mathrm{in}}$.

######  Lemma 4.2 (Lemma 5 of [1]).

Assume we are given the value $\|U_{\mathrm{in}}\|$, and let $O_{x}$ be an
oracle that maps $|00\ldots 0\rangle\in{\mathbb{C}}^{n}$ to a normalized
quantum state $|U_{\mathrm{in}}\rangle$ proportional to $U_{\mathrm{in}}$.
Then the quantum state $|\hat{y}_{\mathrm{in}}\rangle$ proportional to
$\hat{y}_{\mathrm{in}}$ can be prepared using $\mathcal{O}(N)$ queries to
$O_{x}$, and the gate complexity is larger by an $\mathcal{O}(\text{poly}(\log
N,\log n))$ factor.

We remark that, in fact, we embed $\hat{y}_{\mathrm{in}}$ into a slightly
larger space with a more convenient tensor product structure. Further details
refer to Section 4.3 of [1].

###  4.2 Quantum linear system algorithm

After the state preparation of the initial condition, we perform the forward
Euler method to discretize the time interval $[0,T]$ into $m=T/h$ sub-
intervals, and construct a system of linear equations as

$\frac{y^{k+1}-y^{k}}{h}=Ay^{k}$ (4.3)

where $y^{k}\in{\mathbb{R}}^{\mathcal{N}_{d,N}}$ approximates $\hat{y}(kh)$
for each $k\in[{m+1}]_{0}\coloneqq\\{0,1,\ldots,m\\}$, with
$y^{0}=y_{\mathrm{in}}\coloneqq\hat{y}(0)=\hat{y}_{\mathrm{in}}$. This gives
an $(m+1)\mathcal{N}_{d,N}\times(m+1)\mathcal{N}_{d,N}$ linear system

$L|Y\rangle=|B\rangle,$ (4.4)

where

$L=\sum_{k=0}^{m+1}|k\rangle\langle k|\otimes
I-\sum_{k=1}^{m+1}|k\rangle\langle k-1|\otimes(I+Ah).$ (4.5)

(4.4) encodes (4.3) and uses it to produce a numerical solution at time $T$.
Observe that the system (4.4) has the lower triangular structure

$\begin{pmatrix}I&&&&\\\ -(I\\!+\\!Ah)&I&&&\\\ &\ddots&\ddots&&\\\
&&-(I\\!+\\!Ah)&I&\\\ &&&-(I\\!+\\!Ah)&I\\\
\end{pmatrix}\begin{pmatrix}y^{0}\\\ y^{1}\\\ \vdots\\\ y^{m-1}\\\ y^{m}\\\
\end{pmatrix}=\begin{pmatrix}y_{\mathrm{in}}\\\ 0\\\ \vdots\\\ 0\\\ 0\\\
\end{pmatrix}.$ (4.6)

For each $\mathcal{N}_{d,N}$-dim vector $y^{k}$ with $k\in[{m+1}]_{0}$, its
first $n$ components (i.e., $y_{1}^{k}$) approximate the exact solution
$U(kh)$, up to normalization. We apply the high-precision quantum linear
system algorithm (QLSA) [38] to (4.4) and postselect to produce
$\sum_{k}|y_{1}^{k}\rangle|k\rangle$ for representing the gradient flow
evolution. We would like to note that a more advanced QLSA with block-encoding
input models was recently proposed in [81]. However, for technical simplicity,
in this work we still employ the algorithm described in [38]. This is because
the improvements introduced in [81] over [38] are relatively minor, only
affecting logarithmic factors. Additionally, the input model used in this
work, which involves sparse input oracles, is more consistent with that
employed in [38].

In Theorem 4.1, the solution error has two contributions: the error in the
time discretization of (3.6) by the forward Euler method, and the error from
the QLSA. Since the QLSA produces a solution with error at most $\epsilon$
with complexity $\operatorname{poly}(\log(1/\epsilon))$ [38], we focus on
bounding the first contribution.

We provide an upper bound for the error incurred by approximating (3.6) with
the forward Euler method. The following proof basically follows [1, Lemma 3].

######  Lemma 4.3.

Consider an instance of the quantum ODE problem as defined in Problem 1, with
$R_{D}<1$ as defined in (3.14). Choose a time step

$h\leqslant\frac{1}{N^{2}[4Dd(n+1)^{2}+a]}.$ (4.7)

Then the global error from the forward Euler method (4.3) on the interval
$[0,T]$ for (3.6) satisfies

$\|\hat{y}_{1}(T)-y^{m}_{1}\|\leqslant\frac{N^{2}Th}{2}[4Dd(n+1)^{2}+a+b]^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|.$
(4.8)

###### Proof.

First of all, we establish the following bound

$\|I+Ah\|\leqslant 1.$ (4.9)

We decompose $I+Ah$ as

$I+Ah=\sum_{j=1}^{N}K_{j}$ (4.10)

where

$\displaystyle K_{j}$ $\displaystyle=\frac{1}{N}I+|j\rangle\langle j|\otimes
A_{j}^{j}h+|j\rangle\langle j+M-1|\otimes A_{j+M-1}^{j}h,\quad j\in[{N-M+1}],$
(4.11) $\displaystyle K_{j}$ $\displaystyle=\frac{1}{N}I+|j\rangle\langle
j|\otimes A_{j}^{j}h,\quad j\in[{N-1}]\backslash[{N-M+1}].$ (4.12)

All eigenvalues of $\frac{1}{N}I+|j\rangle\langle j|\otimes A_{j}^{j}h$ range
from $\frac{1}{N}+j\lambda_{n}h$ to $\frac{1}{N}+j\lambda_{1}h$. Here we
require that these eigenvalues lie in $[0,1]$, given by
$N^{2}h\|F_{1}\|=N^{2}h[4Dd(n+1)^{2}+a]\leqslant 1$ in (4.7). The norm of
$A_{j+M-1}^{j}$ is bounded by $j\|F_{M}\|$. So we have the bound

$\|K_{j}\|\leqslant\frac{1}{N}-j|\lambda_{1}|h+j\|F_{M}\|h,\quad
j\in[{N-M+1}],$ (4.13)

Then $\|F_{M}\|\leqslant|\lambda_{1}|$ in Problem 1 gives

$\|K_{j}\|\leqslant\frac{1}{N},\quad j\in[{N-M+1}].$ (4.14)

It also holds for the case $j\in[{N-1}]\backslash[{N-M+1}]$. Henceforth,

$\|I+Ah\|\leqslant\sum_{j=1}^{N}\|K_{j}\|\leqslant 1.$ (4.15)

We then define the global error by

$g^{k}\coloneqq\|\hat{y}(kh)-y^{k}\|,$ (4.16)

where $\hat{y}(kh)$ is the exact solution of (3.6), and $y^{k}$ is the
numerical solution of (4.3). Note that $g^{m}=\|\hat{y}(T)-y^{m}\|$.

The stable condition (4.15) implies the local truncation error from the
forward Euler method is non-increasing, and the global error increase at most
linear in time. Following the standard procedure of the global error estimate
(we refer it to the proof of [1, Lemma 3]), the global error is bounded by

$g^{k}\leqslant\frac{kh^{2}}{2}\max_{t\in[0,T]}\|\hat{y}^{\prime\prime}(t)\|,\quad
k\in[{m+1}]_{0},$ (4.17)

where we have the following estimate

$\max_{t\in[0,T]}\|\hat{y}^{\prime\prime}(t)\|=\max_{t\in[0,T]}\|A^{2}\hat{y}(t)\|\leqslant\|A\|^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|,$
(4.18) $\|A\|=\biggl{\|}\sum_{j=1}^{N}|j\rangle\langle j|\otimes
A_{j}^{j}+\sum_{j=1}^{N-M+1}|j\rangle\langle j+M-1|\otimes
A_{j+M-1}^{j}\biggr{\|}\leqslant N(\|F_{1}\|+\|F_{M}\|).$ (4.19)

Sequentially, we conclude that

$\displaystyle\|\hat{y}^{1}(T)-y^{m}_{1}\|$ $\displaystyle\leqslant
g^{m}\leqslant\frac{N^{2}Th}{2}(\|F_{1}\|+\|F_{M}\|)^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|$
(4.20)
$\displaystyle\leqslant\frac{N^{2}Th}{2}[4Dd(n+1)^{2}+a+b]^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|.$

∎

Given the above linear system, we can upper bound the condition number that
affects the complexity of the quantum linear system algorithm. Under the same
construction of the matrix $L$, we can follow the same estimate proposed by
[1, Lemma 4] to claim the following result.

######  Lemma 4.4.

Consider an instance of the quantum ODE problem as defined in Problem 1. Apply
the forward Euler method (4.3) with time step size (4.7) to the Carleman
linearization (3.6). Then the condition number of the matrix $L$ defined in
(4.5) satisfies

$\kappa\leqslant 2(m+1).$ (4.21)

###### Proof.

We begin by upper bounding $\|L\|$. We write

$L=L_{1}+L_{2},$ (4.22)

where

$\displaystyle L_{1}$ $\displaystyle=\sum_{k=0}^{m}|k\rangle\langle k|\otimes
I,$ (4.23) $\displaystyle L_{2}$
$\displaystyle=-\sum_{k=1}^{m}|k\rangle\langle k-1|\otimes(I+Ah).$ (4.24)

Clearly $\|L_{1}\|=1$. Furthermore, $\|L_{2}\|\leqslant\|I+Ah\|\leqslant 1$ by
(4.15). Therefore,

$\|L\|\leqslant\|L_{1}\|+\|L_{2}\|\leqslant 2.$ (4.25)

Next we upper bound $\|L^{-1}\|$. We notice that $L^{-1}$ can be directly
written as

$L^{-1}=\begin{pmatrix}I&&&&\\\ (I\\!+\\!Ah)&I&&&\\\
(I\\!+\\!Ah)^{2}&(I\\!+\\!Ah)&I&&\\\ \ldots&\ddots&\ddots&\ddots&\\\
(I\\!+\\!Ah)^{m}&\cdots&(I\\!+\\!Ah)^{2}&(I\\!+\\!Ah)&I\\\ \end{pmatrix}.$
(4.26)

So that

$\|L^{-1}\|\leqslant\|I\|+\|I+Ah\|+\ldots+\|(I+Ah)^{m}\|.$ (4.27)

Since $\|I+Ah\|\leqslant 1$ by (4.15), we have

$\|L^{-1}\|\leqslant\|I\|+\|I+Ah\|+\ldots+\|(I+Ah)\|^{m}=m+1.$ (4.28)

Finally, combining (4.25) with (4.28), we conclude

$\kappa=\|L\|\|L^{-1}\|\leqslant 2(m+1)$ (4.29)

as claimed. ∎

###  4.3 Measurement probability

After applying the QLSA to (4.4), we perform a measurement to extract a final
state of the desired form. We now consider the probability of this measurement
succeeding. Differing from [1, Lemma 6], we are interested in providing a
history state $\sum_{k}|y_{1}^{k}\rangle|k\rangle$ instead of a final state
$|y_{1}^{m}\rangle$. Thus the measurement probability does not include the
$\ell_{2}$ norm of the final state as well as the $\ell_{2}$ scaling of the
initial and final states (i.e., the parameters $g$ and $q$ as in [1, Lemma
6]).

######  Lemma 4.5.

Consider an instance of the quantum ODE problem defined in Problem 1, with the
QLSA applied to the linear system (4.4) using the forward Euler method (4.3)
with time step (4.7). Suppose the global error from the forward Euler method
as defined in Lemma 4.3 is bounded by

$\|\hat{y}(kh)-y^{k}\|\leqslant\frac{1}{2}G.$ (4.30)

Then the probability of measuring a quantum state $|y^{k}_{1}\rangle$ for
$k\in[{m+1}]_{0}$ satisfies

$P_{\mathrm{measure}}\geqslant\frac{2G^{2}}{16\max_{t\in[0,T]}\|\hat{y}(t)\|^{2}+G^{2}}.$
(4.31)

###### Proof.

The quantum state produced by the QLSA applied to (4.4) has the form

$|Y\rangle=\frac{1}{\||Y\rangle\|}\sum_{k=0}^{m}y^{k}|k\rangle=\frac{1}{\|Y\|}\sum_{k=0}^{m}\sum_{j=1}^{N}y_{j}^{k}|j\rangle|k\rangle$
(4.32)

where the normalization factor satisfies
$\||Y\rangle\|^{2}=\sum_{k=0}^{m}\|y^{k}\|^{2}=\sum_{k=0}^{m}\sum_{j=1}^{N}\|y_{j}^{k}\|^{2}$.

We aim to obtain the target quantum state as the form

$|Y_{\mathrm{evo}}\rangle=\frac{1}{\|Y_{\mathrm{evo}}\|}\sum_{k=0}^{m}y^{k}_{1}|k\rangle.$
(4.33)

which corresponds to the gradient flow evolution state (2.15) We measure the
register $|j\rangle$, $j\in[{N}]$ and extract $|Y_{\mathrm{target}}\rangle$
from $|Y\rangle$ when $k=1$. The success probability is lower bounded as
below.

According to the Cauchy–Schwarz inequality,

$\|\hat{y}_{1}(kh)\|^{2}=\|\hat{y}_{1}(kh)-y^{k}_{1}+y^{k}_{1}\|^{2}\leqslant
2\|\hat{y}_{1}(kh)-y^{k}_{1}\|^{2}+2\|y^{k}_{1}\|^{2},$ (4.34)

so that

$\||y^{k}_{1}\rangle\|^{2}\geqslant\frac{1}{2}\|\hat{y}_{1}(kh)\|^{2}-\|\hat{y}_{1}(kh)-y^{k}_{1}\|^{2}\geqslant\frac{1}{2}\|\hat{y}_{1}(kh)\|^{2}-\frac{1}{4}G^{2}.$
(4.35)

Summing $k$ from $0$ to $m$, and using the definition of $G$ in (4.1), we have

$\displaystyle\||Y_{\mathrm{evo}}\rangle\|^{2}$
$\displaystyle=\sum_{k=0}^{m}\|y_{1}^{k}\|^{2}\geqslant\frac{1}{2}\sum_{k=0}^{m}\|\hat{y}_{1}(kh)\|^{2}-\frac{m+1}{4}G^{2}$
(4.36)
$\displaystyle\geqslant\frac{m+1}{2}G^{2}-\frac{m+1}{4}G^{2}=\frac{m+1}{4}G^{2}.$

Second, we use (4.30) and the parallel inequality again to upper bound
$\|y^{k}\|^{2}$ by

$\displaystyle\|y^{k}\|^{2}\leqslant
2\|\hat{y}(kh)\|^{2}+2\|\hat{y}(kh)-y^{k}\|^{2}\leqslant
2\|\hat{y}(kh)\|^{2}+\frac{1}{2}G^{2}.$ (4.37)

Therefore

$\displaystyle\||Y\rangle\|^{2}=\sum_{k=0}^{m}\|y^{k}\|^{2}\leqslant
2\sum_{k=0}^{m}\|\hat{y}(kh)\|^{2}+\frac{m+1}{2}G^{2}\leqslant
2(m+1)\max_{t\in[0,T]}\|\hat{y}(t)\|^{2}+\frac{m+1}{2}G^{2}.$ (4.38)

Finally, using (4.36) and (4.38), we have

$P_{\mathrm{measure}}=\frac{\||Y_{\mathrm{evo}}\rangle\|^{2}}{\||Y\rangle\|^{2}}\geqslant\frac{G^{2}}{8\max_{t\in[0,T]}\|\hat{y}(t)\|^{2}+2G^{2}}$
(4.39)

as claimed.

∎

Using amplitude amplification, $\mathcal{O}(\sqrt{1/P_{\mathrm{measure}}})$
iterations suffice to succeed with constant probability.

###  4.4 Proof of Theorem 4.1

###### Proof.

We first present the quantum Carleman linearization (QCL) algorithm and then
analyze its complexity.

#### The QCL algorithm.

We introduce the choice of parameters as follows. Given an error bound
$\epsilon\leqslant 1$ and $G$, we define

$\delta\coloneqq\frac{G\epsilon}{1+\epsilon},$ (4.40)

which satisfies $\delta\leqslant G/2$ for any $t\in[0,T]$.

Now we discuss the choice of $h$. On the one hand, $h$ must follow (4.7) to
satisfy the conditions of Lemma 4.3 and Lemma 4.4. On the other hand, we
choose

$\displaystyle
h\leqslant\min\biggl{\\{}\frac{1}{N^{2}[4Dd(n+1)^{2}+a]},\frac{G\epsilon}{N^{2}T[4Dd(n+1)^{2}+a+b]^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|}\biggr{\\}}$
(4.41)

Then according to the requirement (4.7) in Lemma 4.3, and for
$k\in[{m+1}]_{0}$,

$\|\hat{y}_{1}(kh)-y^{k}_{1}\|\leqslant\|\hat{y}(kh)-y^{k}\|\leqslant\frac{N^{2}Th}{2}[4Dd(n+1)^{2}+a+b]^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|\leqslant\delta.$
(4.42)

It also leads to $\|\hat{y}(kh)-y^{k}\|\leqslant G/2$ used as a condition in
Lemma 4.5.

We now consider the error between the exact and numerical gradient flow
evolutions

$|\hat{Y}_{\mathrm{evo}}\rangle=\frac{1}{\|\hat{Y}_{\mathrm{evo}}\|}\hat{Y}_{\mathrm{evo}}=\frac{1}{\|\hat{Y}_{\mathrm{evo}}\|}\sum_{k=0}^{m}\hat{y}_{1}(kh)|k\rangle$
(4.43)

and (as denoted in (4.33))

$|Y_{\mathrm{evo}}\rangle=\frac{1}{\|Y_{\mathrm{evo}}\|}Y_{\mathrm{evo}}=\frac{1}{\|Y_{\mathrm{evo}}\|}\sum_{k=0}^{m}y^{k}_{1}|k\rangle,$
(4.44)

where $\|\hat{Y}_{\mathrm{evo}}\|$ and $\|Y_{\mathrm{evo}}\|$ are
normalization factors. Recall the definition of $G$ in (4.1)

$G=\frac{1}{\sqrt{m+1}}\|\hat{Y}_{\mathrm{evo}}\|=\sqrt{\frac{\sum_{k=0}^{m}\|\hat{y}_{1}(kh)\|^{2}}{m+1}},$
(4.45)

The $\ell_{2}$ normalized error can be controlled by

$\biggl{\|}|\hat{Y}_{\mathrm{evo}}\rangle-|Y_{\mathrm{evo}}\rangle\biggr{\|}\leqslant\frac{\|\hat{Y}_{\mathrm{evo}}-Y_{\mathrm{evo}}\|}{\min\\{\|\hat{Y}_{\mathrm{evo}}\|,\|Y_{\mathrm{evo}}\|\\}}\leqslant\frac{\|\hat{Y}_{\mathrm{evo}}-Y_{\mathrm{evo}}\|}{\|\hat{Y}_{\mathrm{evo}}\|-\|\hat{Y}_{\mathrm{evo}}-Y_{\mathrm{evo}}\|}.$
(4.46)

Then using (4.42), since

$\|\hat{Y}_{\mathrm{evo}}-Y_{\mathrm{evo}}\|^{2}\leqslant\sum_{k=0}^{m}\|\hat{y}_{1}(kh)-y^{k}_{1}\|^{2}\leqslant(m+1)\delta^{2},$
(4.47)

we have

$\|\hat{Y}_{\mathrm{evo}}-Y_{\mathrm{evo}}\|\leqslant\sqrt{m+1}\delta,$ (4.48)

which gives

$\biggl{\|}|\hat{Y}_{\mathrm{evo}}\rangle-|Y_{\mathrm{evo}}\rangle\biggr{\|}\leqslant\frac{\sqrt{m+1}\delta}{\|\hat{Y}_{\mathrm{evo}}\|-\sqrt{m+1}\delta}=\frac{\delta}{G-\delta}=\epsilon,$
(4.49)

i.e., $\epsilon$ upper bounds the $\ell_{2}$ normalized error between
$|\hat{Y}_{\mathrm{evo}}\rangle$ and $|Y_{\mathrm{evo}}\rangle$.

We follow the procedure in Lemma 4.2 to prepare the initial state
$|\hat{y}_{\mathrm{in}}\rangle$. We apply the QLSA [38] to the linear system
(4.4) with $m=\lceil T/h\rceil$, giving a solution $|Y\rangle$. By Lemma 4.5,
the probability of obtaining a state is
$\sum_{k=0}^{m}|y_{1}^{k}\rangle|k\rangle$

$P_{\mathrm{measure}}\geqslant\frac{2G^{2}}{16\max_{t\in[0,T]}\|\hat{y}(t)\|^{2}+G^{2}}.$
(4.50)

By amplitude amplification, we can achieve success probability $\Omega(1)$
with $\mathcal{O}(\max_{t\in[0,T]}\|\hat{y}(t)\|/G)$ repetitions of the above
procedure.

#### Analysis of the complexity.

By Lemma 4.2, the initial state $|\hat{y}_{\mathrm{in}}\rangle$ can be
prepared with $\mathcal{O}(N)$ queries to $O_{x}$, with gate complexity larger
by a $\mathcal{O}(\text{poly}(\log N,\log n))$ factor. The matrix $L$ is an
$(m+1)\mathcal{N}_{d,N}\times(m+1)\mathcal{N}_{d,N}$ matrix with
$\mathcal{O}(Ns)$ nonzero entries in any row or column. By Lemma 4.4 and our
choice of parameters, the condition number of $L$ is at most

$\displaystyle 2(m+1)$ (4.51)
$\displaystyle\quad=O\biggl{(}\frac{1}{G\epsilon}N^{2}T^{2}[4Dd(n+1)^{2}+a+b]^{2}\max_{t\in[0,T]}\|\hat{y}(t)\|+N^{2}T[4Dd(n+1)^{2}+a]\biggr{)}$
$\displaystyle\quad=O\Bigl{(}\frac{1}{G\epsilon}N^{2}T^{2}D^{2}d^{2}n^{4}\max_{t\in[0,T]}\|\hat{y}(t)\|\Bigr{)}.$

Here we use $(\|F_{1}\|+\|F_{M}\|)^{2}=[4Dd(n+1)^{2}+a+b]^{2}$, and
$|a|,|b|=o(d)$. Consequently, by Theorem 5 of [38], the QLSA produces the
state $|Y\rangle$ with

$\displaystyle\frac{1}{G\epsilon}sT^{2}D^{2}d^{2}n^{4}N^{3}\max_{t\in[0,T]}\|\hat{y}(t)\|\cdot\operatorname{poly}\biggl{(}\log\frac{aDdnNsT}{G\epsilon}\biggr{)}$
(4.52)

queries to the oracles $O_{F_{1}}$ and $O_{F_{2}}$. Using
$O(\max_{t\in[0,T]}\|\hat{y}(t)\|/G)$ steps of amplitude amplification to
achieve success probability $\Omega(1)$, the overall query complexity of our
algorithm is

$\displaystyle\frac{1}{G^{2}\epsilon}sT^{2}D^{2}d^{2}n^{4}N^{3}\max_{t\in[0,T]}\|\hat{y}(t)\|^{2}\cdot\operatorname{poly}\biggl{(}\log\frac{aDdnNsT}{G\epsilon}\biggr{)}$
(4.53)

and its gate complexity is larger than its query complexity only by
logarithmic factors, based on the gate-efficient algorithm in Theorem 5 of
[38].

We now estimate the quantity $\max_{t\in[0,T]}\|\hat{y}(t)\|$. By the
definition of $\eta_{j}(t)$, Theorem 3.3, and $R_{D}<1$, we have

$\|\hat{y}_{j}(t)\|=\|U^{\otimes j}(t)-\eta_{j}(t)\|\leqslant\|U^{\otimes
j}(t)\|+\|\eta_{j}(t)\|\leqslant\max_{t\in[0,T]}\|U(t)\|^{j}+\gamma^{j}\leqslant
2\max_{t\in[0,T]}\|U(t)\|^{j},$ (4.54)

so that

$\max_{t\in[0,T]}\|\hat{y}(t)\|\leqslant
2\sum_{j=1}^{N}(\max_{t\in[0,T]}\|U(t)\|)^{j}.$ (4.55)

Based on the $\ell_{2}$ estimate of the solution in Lemma C.4, when $R_{D}<1$,
we have

$\max_{t\in[0,T]}\|U(t)\|\leqslant\|U_{\mathrm{in}}\|.$ (4.56)

Therefore, the overall query complexity of our algorithm is

$\displaystyle\frac{1}{G^{2}\epsilon}sT^{2}D^{2}d^{2}n^{4}N^{3}\|U_{\mathrm{in}}\|^{2N}\cdot\operatorname{poly}\biggl{(}\log\frac{aDdMnNsT}{G\epsilon}\biggr{)}$
(4.57)

and the gate complexity is larger than its query complexity by logarithmic
factors as claimed. ∎

##  5 Efficiency and lower bound results

The reference [1] established a limitation on the ability of quantum computers
to solve the quadratic ODE problem when the nonlinearity is sufficiently
strong. In other words, general nonlinear differential equations are
intractable on quantum computers when $R\geqslant\sqrt{2}$. However, we can
rule out such a worst-case by assuming that the initial condition of reaction-
diffusion equations fulfills the maximum principle, and thus show the problem
is still tractable on quantum computers.

In the following, we state and prove our hardness and efficiency results. Part
(i) focuses on the hardness result when $R\geq\sqrt{2}$, which has been
studies [1, Theorem 2] by leveraging the hardness result of quantum state
discrimination. However, there is a technical flaw in the original proof in
[1]. The hardness result for quantum state discrimination used in [1] only
assumes multiple copies of the input states at the beginning and does not
allow access to the state during the algorithm. But, in most quantum ODE
algorithms, including the Carleman linearization method, we indeed have a
stronger assumption that we assume the state preparation oracle for the input
state and its _inverse_ , and we frequently apply those during the
implementation of the algorithm. Therefore the existing lower bound in [1] has
not yet fully ruled out the possibility of efficient algorithms with strong
oracle assumptions. We fix this gap in part (i) by applying a recent lower
bound for amplifiers [82, Theorem 13], where the state preparation oracles are
assumed. Part (ii) shows that the worst-case scenario can be precluded by
assuming the maximum principle, implying that our maximum principle analysis
captures the underlying reason for the efficiency of Carleman linearization
method.

######  Theorem 5.1.

We consider the same assumptions in Problem 1.
(i) Assume $R\geqslant\sqrt{2}$, and the initial condition satisfies
$\|U_{\mathrm{in}}\|_{\infty}>\gamma$. Then there is an instance of the
quantum quadratic ODE problem defined in Problem 1 such that any quantum
algorithm for producing a quantum state approximating $u(T)/\|u(T)\|$ with
bounded error must have worst-case query complexity exponential in $T$ to the
input state preparation oracle.
(ii) If the initial condition satisfies the maximum principle
$\|U_{\mathrm{in}}\|_{\infty}\leqslant\gamma$ as (2.12), then such a worst-
case example can be precluded even $R\geqslant\sqrt{2}$.

###### Proof.

We consider the lower bound result when $R\geqslant\sqrt{2}$ and
$\|U_{\mathrm{in}}\|_{\infty}>\gamma$. The same as Theorem 2 of [1], we
consider a $2$-dimensional system of the form

$\displaystyle\frac{\mathrm{d}{u_{1}}}{\mathrm{d}{t}}$
$\displaystyle=-u_{1}+Ru_{1}^{2},$ (5.1)
$\displaystyle\frac{\mathrm{d}{u_{2}}}{\mathrm{d}{t}}$
$\displaystyle=-u_{2}+Ru_{2}^{2},$

with two single-qubit states as initial conditions

$|\phi(0)\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ (5.2)

and

$|\psi(0)\rangle=v_{0}|0\rangle+w_{0}|1\rangle:=\cos\Bigl{(}\theta+\frac{\pi}{4}\Bigr{)}|0\rangle+\sin\Bigl{(}\theta+\frac{\pi}{4}\Bigr{)}|1\rangle,$
(5.3)

where $\theta\in(0,\pi/4)$, $2\sin^{2}\frac{\theta}{2}=\epsilon$, with an
arbitrary small $\epsilon>0$. Then the overlap between the two initial states
is

$\langle\phi(0)|\psi(0)\rangle:=\cos\theta=1-\epsilon.$ (5.4)

We let $U(t)=[v(t);w(t)]$ denote the solution evolved from
$U(0)=[v_{0};w_{0}]$. According to Lemma 8 of [1],
$|\phi(t)\rangle=|\phi(0)\rangle$ is the fixed state; but if $1/R\geqslant
1/\sqrt{2}$, $w(t)$ increases with $t$ and goes to infinity after

$t>t^{\ast}\coloneqq\log\biggl{(}\frac{R}{R-1/w_{0}}\biggr{)}.$ (5.5)

The overlap of $|\phi(T)\rangle$ and $|\psi(T)\rangle$ is no larger than a
constant (e.g., $\frac{3}{\sqrt{10}}$ used in [1]) after a short evolution
time

$T<t^{\ast}=\log\biggl{(}\frac{R}{R-\frac{1}{w_{0}}}\biggr{)}<\log\biggl{(}\frac{\sqrt{2\epsilon-\epsilon^{2}}+1-\epsilon}{\sqrt{2\epsilon-\epsilon^{2}}-\epsilon}\biggr{)}=\mathcal{O}(\log({1}/{\epsilon})).$
(5.6)

It was shown in [82, Theorem 13] that, if a quantum algorithm with oracle
input model can amplify the infidelity of two quantum states from $\epsilon$
to a constant level, then it must use $\Omega(1/\sqrt{\epsilon})$ queries in
the worst case. By applying this result and noticing that
$1/\sqrt{\epsilon}=e^{\Omega(T)}$, we directly obtain that when
$R\geqslant\sqrt{2}$, there is an instance of the quantum quadratic ODE
problem that any quantum algorithm must have worst-case time complexity
exponential in $T$.

In our paper, the ODE system (5.1) is a reduced example of reaction-diffusion
equations (2.6) with $D=0$, $a=-1$, $b=R$, $M=2$, and $R_{D}=1$. Besides,
$\|U_{\mathrm{in}}\|_{\infty}=\sin\Bigl{(}\theta+\frac{\pi}{4}\Bigr{)}$
satisfies

$1/\sqrt{2}<\sin\Bigl{(}\theta+\frac{\pi}{4}\Bigr{)}=(\sqrt{2\epsilon-\epsilon^{2}}+1-\epsilon)/\sqrt{2}<(1+\sqrt{2\epsilon})/\sqrt{2},$
(5.7)

where $\sin\Bigl{(}\theta+\frac{\pi}{4}\Bigr{)}$ is close to $1/\sqrt{2}$ when
$\epsilon$ is close to $0$. Notice that this example disobeys (2.12), the
condition of the Maximum Principle Lemma 2.1, because

$\|U_{\mathrm{in}}\|_{\infty}>\gamma=1/\sqrt{2}.$ (5.8)

Secondly, we consider an upper bound on $\|U_{\mathrm{in}}\|$ given the
maximum principle $\|U_{\mathrm{in}}\|_{\infty}\leqslant\gamma$. Then we have

$\|U_{\mathrm{in}}\|\leqslant\sqrt{n_{d}}\gamma.$ (5.9)

Substituting this estimate into the complexity in Theorem 4.1, we can upper
bound the query complexity by

$\displaystyle\frac{1}{G^{2}\epsilon}sT^{2}D^{2}d^{2}n^{4}N^{3}n_{d}^{N}\gamma^{2N}\cdot\operatorname{poly}\biggl{(}\log\frac{aDdMnNsT}{G\epsilon}\biggr{)}.$
(5.10)

Notice that the upper bound of the complexity still depends exponentially on
$N$. However, according to Theorem 3.3, the Carleman error converges
exponentially in $N$ and can be bounded independently of $T$. So $N$ can be
chosen independently of $T$ as well, therefore our algorithm does not have the
exponential overhead in $T$ stated in part (i). ∎

The upper bounds on the query and gate complexity demonstrate that the quantum
algorithm we develop has a roughly quadratic dependence on $T$ when $R_{D}<1$
and $\|U_{\mathrm{in}}\|_{\infty}\leqslant\gamma$, regardless of the value of
$R$. Such a loose upper bound includes a polynomial dependence on $n$,
revealing that the quantum algorithm does not have a potential exponential
speedup in the dimension.

But if we are given an additional assumption

$\|U_{\mathrm{in}}\|=o(n_{d}),$ (5.11)

then we can upper bound the query complexity by

$\displaystyle\frac{1}{G^{2}\epsilon}sT^{2}D^{2}d^{2}n^{4}N^{3}\cdot\operatorname{poly}\biggl{(}\log\frac{aDdMnNsT}{G\epsilon}\biggr{)},$
(5.12)

and the gate complexity has an upper bound that is larger by logarithmic
factors as claimed. In this case, our algorithm still maintains the potential
exponential speedup in the dimension over classical algorithms.

##  6 Applications

In this section, we show how the quantum state obtained by solving Problem 1
can be used to compute quantities of practical interest. For generality, in
this section, we consider the applications of a quantum state in the form
specified in Problem 1, but not limited to the output by particular
algorithms. More specifically, let $l=(l_{1},\ldots,l_{d})$ with
$l_{j}\in[{n}]$, and let $f(t,x)$ be a function defined on
$[0,T]\otimes[0,1]^{d}$. We assume that there exists a quantum algorithm
$\mathcal{A}(\epsilon)$ which can prepare the quantum state

$|\hat{f}\rangle\sim\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}\hat{f}_{k,l}{|{k}\rangle}|l_{1}\rangle\ldots|l_{d}\rangle$
(6.1)

proportional to the vector $(f(kT/m,l_{1}/n,\cdots,l_{d}/n))$ within some
prescribed error tolerance $\epsilon>0$ in $\ell^{2}$ norm. Here
$\hat{f}_{k,l}$ represents an approximation of the function $f$ evaluated at
$(kT/m,l_{1}/n,\cdots,l_{d}/n)$. In the context of this paper, the algorithm
$\mathcal{A}$ is the quantum Carleman linearization method, but the discussion
in this section works for any quantum algorithm which can encode a function
evaluated at discrete grid points.

###  6.1 Mean square amplitude

One quantity of potential practical interest is the fraction of the squared
amplitude contained in a sub-domain ${\mathscr{D}}_{t}\times{\mathscr{D}}_{x}$
defined by ${\mathscr{D}}_{t}\subset[0,T]$ and
${\mathscr{D}}_{x}\subset[0,1]^{d}$. This can be described as the ratio

$\frac{\int_{{\mathscr{D}}_{t}}\int_{{\mathscr{D}}_{x}}|f(t,x)|^{2}dxdt}{\int_{0}^{T}\int_{[0,1]^{d}}|f(t,x)|^{2}dxdt}.$
(6.2)

In the context of quantum mechanics such quantities are motivated by Born’s
rule, whereas in the context of classical wave mechanics such quantities are
motivated by notions of energy.

In the spatial discretization context, we can approximate the integrals via
numerical quadrature with equidistant nodes, and thus we are interested in
computing the ratio

$\frac{\sum_{kT/m\in{\mathscr{D}}_{t}}\sum_{(l_{1}/n,\cdots,l_{d}/n)\in{\mathscr{D}}_{x}}|f(kT/m,l_{1}/n,\cdots,l_{d}/n)|^{2}}{\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}|f(kT/m,l_{1}/n,\cdots,l_{d}/n)|^{2}}.$
(6.3)

Note that the difference between (6.2) and (6.3) scales as
$\mathcal{O}(T(T/m+d/n))$ [83], so it can be reduced to the level of
$\mathcal{O}(\epsilon)$ with arbitrarily small $\epsilon$ by refining time and
spatial discretization as
$m=\mathcal{O}(T^{2}/\epsilon),n=\mathcal{O}(dT/\epsilon)$. For simplicity,
here we fix the grids for discretization and focus on computing the
discretized ratio as shown in (6.3). This quantity can be easily estimated
given the quantum state in the form of (6.1). Let $P$ be the projector onto
the space spanned by the computational basis corresponding the $\Omega$, i.e.,
$P=\sum_{kT/m\in{\mathscr{D}}_{t}}\sum_{(l_{1}/n,\cdots,l_{d}/n)\in{\mathscr{D}}_{x}}({|{k}\rangle}{\langle{k}|})\otimes({|{l_{1}}\rangle}{\langle{l_{1}}|})\otimes\cdots\otimes({|{l_{d}}\rangle}{\langle{l_{d}}|})$.
Then

$\begin{split}&\quad\frac{\sum_{kT/m\in{\mathscr{D}}_{t}}\sum_{(l_{1}/n,\cdots,l_{d}/n)\in{\mathscr{D}}_{x}}|f(kT/m,l_{1}/n,\cdots,l_{d}/n)|^{2}}{\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}|f(kT/m,l_{1}/n,\cdots,l_{d}/n)|^{2}}\\\
&\approx\frac{\sum_{kT/m\in{\mathscr{D}}_{t}}\sum_{(l_{1}/n,\cdots,l_{d}/n)\in{\mathscr{D}}_{x}}|\hat{f}_{k,l}|^{2}}{\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}|\hat{f}_{k,l}|^{2}}\\\
&=\langle\hat{f}|P|\hat{f}\rangle\end{split}$ (6.4)

and thus can be estimated by the amplitude estimate technique [72].

The complexity of such an algorithm is given in the following theorem, and the
proof can be found in the appendix Section G.

######  Theorem 6.1.

Assume that we are given an algorithm $\mathcal{A}(\epsilon^{\prime})$
preparing (6.1) and a unitary transform $(I-2P)$ for the projector $P$
associated with the domain ${\mathscr{D}}$. Then, for any $0<\epsilon<1$,
$0<\delta<1$, there exists a quantum algorithm which can output an
approximation of (6.3) within tolerated error $\epsilon$ and probability at
least $(1-\delta)$, using $\mathcal{A}(\epsilon/4)$ and $(I-2P)$ for
$\mathcal{O}((1/\epsilon)\log(1/\delta))$ times.

Theorem 6.1 proves an almost linear sampling scaling in terms of the precision
for the quantum approach to estimate the ratio. Compared with the most widely
used classical approach for estimating high dimensional integral, standard
Monte Carlo methods, which typically scale quadratically in the precision, we
can obtain a quadratic speedup for the sampling step.

###  6.2 Derivatives and kinetic energy

As a particular example of practical interest, we would like to study the
dynamical kinetic energy ratio of the system on the domain
${\mathscr{D}}_{t}\times{\mathscr{D}}_{x}\subset[0,T]\times[0,1]^{d}$, which
is defined to be

$\frac{\int_{{\mathscr{D}}_{t}}\int_{{\mathscr{D}}_{x}}|\nabla_{x}f(t,x)|^{2}dxdt}{\int_{0}^{T}\int_{[0,1]^{d}}|\nabla_{x}f(t,x)|^{2}dxdt}\approx\frac{\sum_{kT/m\in{\mathscr{D}}_{t}}\sum_{(l_{1}/n,\cdots,l_{d}/n)\in{\mathscr{D}}_{x}}|\nabla_{x}f(kT/m,l_{1}/n,\cdots,l_{d}/n)|^{2}}{\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}|\nabla_{x}f(kT/m,l_{1}/n,\cdots,l_{d}/n)|^{2}}.$
(6.5)

To apply our result in Theorem 6.1, we first study how to prepare a quantum
state which is proportional to the partial derivative of the function $f$.
That is, our goal is to prepare a quantum state which is a good approximation
of

${|{\nabla_{x}f}\rangle}\sim\sum_{j=1}^{d}\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}\nabla_{x_{j}}f(kT/m,l_{1}/n,\cdots,l_{d}/n){|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.6)

Compared to (6.1), here we introduce one more register to simultaneously
encode all the partial derivatives in a single quantum state. In this
subsection, we further assume that the function $f(t,x)$ satisfies periodic
boundary conditions for spatial variable $x$. Notice that the periodic
boundary condition here is only for technical simplicity and not essential
because a general function can be smoothly extended to a larger space domain
with periodic boundary conditions. Our algorithm will also work by requiring
access to another projector connecting the original and the extended space
domains.

We first study how to prepare a quantum state encoding the derivative
information within amplitudes. The idea is using the discrete Fourier
transform to transform the state to the frequency domain, multiplying the
frequency in this domain and then transforming back to the space domain by
inverse discrete Fourier transform. More specifically, let $\mathcal{F}$
denote the one-dimensional quantum Fourier transform333Note that, following
the standard convention of choosing signs in [84], the quantum Fourier
transform exactly corresponds to the inverse discrete Fourier transform in the
classical setting. with $n$ nodes. Furthermore, for any positive integer
$\theta\leq n/2$, let $l=(l_{1},\cdots,l_{d})$ and $D_{j,\theta}$ be a
diagonal matrix

$D_{j,\theta}=2\pi
i\theta\sum_{l_{1},\cdots,l_{d}=0}^{n-1}D_{j,\theta}(l){|{l}\rangle}{\langle{l}|}$
(6.7)

where

$D_{j,\theta}(l)=\begin{cases}l_{j}/\theta,&\text{ if }1\leq
l_{j}\leq\theta,\\\ (l_{j}-n)/\theta,&\text{ if }n-\theta\leq l_{j}\leq
n-1,\\\ 0,&\text{ else}.\end{cases}$ (6.8)

Notice that we define $D_{j,\theta}(l)$ such that $|D_{j,\theta}|\leq 1$.
Then, for any smooth function $g(x)$ defined on $[0,1]^{d}$, we have

$\begin{split}&\quad(\otimes_{j-1}I\otimes\mathcal{F}\otimes_{d-j}I)D_{j,\theta}(\otimes_{j-1}I\otimes\mathcal{F}^{-1}\otimes_{d-j}I)\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}g(l_{1}/n_{1},\cdots,l_{d}/n_{d}){|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}\\\
&\approx\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}\partial_{x_{j}}g(l_{1}/n_{1},\cdots,l_{d}/n_{d}){|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.\end{split}$
(6.9)

Such an approach has been widely used and is regarded as the standard way to
compute derivatives in classical scientific computing, and we briefly
illustrate the reasoning in the appendix Section F.

Now we discuss how to implement this approach on a quantum device. In the
general high-dimensional case, we introduce another ancilla register with
$\mathcal{O}(\log d)$ qubits, which we will refer to as the dimension register
later, and start with the state

${|{0}\rangle}{|{f}\rangle}=\frac{1}{\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}f(kT/m,l_{1}/n,\cdots,k_{d}/n){|{0}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.10)

Applying the Hadamard gates to the dimension register, we obtain

$\frac{1}{\sqrt{d}\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}f(k_{1}/n,\cdots,k_{d}/n){|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.11)

The discrete Fourier transform can be efficiently implemented via quantum
Fourier transform. Specifically, we apply the operation
$(\sum_{j=0}^{d-1}{|{j}\rangle}{\langle{j}|}\otimes
I_{n_{s}}\otimes\mathcal{F}_{j}^{-1})$ and denote the resulting state as

$\frac{1}{\sqrt{d}\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}\tilde{f}(j,k,l){|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.12)

For the multiplication of the matrix $D_{j,\theta}$ (we will show later that
$\theta$ can be chosen as an $\mathcal{O}(1)$ parameter for smooth functions),
we assume that we are given an oracle of the mapping

$O_{D}:{|{0}\rangle}{|{j}\rangle}{|{l}\rangle}\rightarrow{|{D_{j,\theta}(l)}\rangle}{|{j}\rangle}{|{l}\rangle}.$
(6.13)

Then the multiplication of the matrix $D_{j,\theta}$ can be implemented as
follows. We first add two ancilla registers, one as the rotation register on
which we will perform conditional rotation later and the other as the
$D$-register for encoding $D_{j,\theta}$. Applying $O_{D}$ to encode
$D_{j,\theta}(l)$ in the $D$-register gives

$\frac{1}{\sqrt{d}\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}\tilde{f}(j,k,l){|{0}\rangle}{|{D_{j,\theta}(l)}\rangle}{|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.14)

Performing a rotation on the rotation register conditioned on
${|{D_{j,\theta}(l)}\rangle}$ yields

$\frac{1}{\sqrt{d}\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}\tilde{f}(j,k,l)\left(D_{j,\theta}(l){|{0}\rangle}+\sqrt{1-D_{j,\theta}(l)^{2}}{|{1}\rangle}\right){|{D_{j,\theta}(l)}\rangle}{|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.15)

Uncomputing the $D$-register gives the state

$\frac{1}{\sqrt{d}\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}D_{j,\theta}(l)\tilde{f}(j,k,l){|{0}\rangle}{|{0}\rangle}{|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}+{|{\perp}\rangle},$
(6.16)

where the first part is the desired outcome after the diagonal transformation,
and ${|{\perp}\rangle}$ represents a quantum state with the rotation register
being ${|{1}\rangle}$. Finally, applying
$(\sum_{j=0}^{d-1}{|{j}\rangle}{\langle{j}|}\otimes\mathcal{F}_{j})$ on the
registers $({|{j}\rangle}\otimes{|{l}\rangle})$ completes the operation for
computing the partial derivatives as discussed before, which yields an
approximation of

$\displaystyle\frac{1}{2\pi
i\theta}\frac{1}{\sqrt{d}\|\vec{f}\|}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}\partial_{x_{j}}f(kT/m,l_{1}/n,\cdots,l_{d}/n){|{0}\rangle}{|{0}\rangle}{|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}+{|{\perp}\rangle}.$
(6.17)

By measuring the ancilla rotation register to get $0$ and discarding the
$D$-register, we get a quantum state approximately proportional to

$\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\sum_{j=0}^{d-1}\partial_{x_{j}}f(k_{1}/n,\cdots,k_{d}/n){|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle},$
(6.18)

which encodes the partial derivatives in the amplitude controlled by a
dimension register. The entire quantum circuit is summarized in Fig. 5, and
the overall complexity estimate is given in the following theorem, of which
the proof can be found in the appendix Section G.

@R=1em @C=1em Rotation & R ${\langle{0}|}$
$D$ 2O_D -1 2O_D^-1
Dimension ⊗H 1 O_D O_D^-1 1
Space 1A F_j^-1 O_D O_D^-1 F_j
Time A

Figure 5:  Quantum circuit for preparing a quantum state encoding partial
derivatives of a known function in amplitudes. Here $\mathcal{A}$ is the
algorithm for the state encoding the function, $\mathcal{H}$ represents the
Hadamard gate, $\mathcal{F}_{j}$ represents the one-dimensional quantum
Fourier transform acting on the $j$-th direction, $O_{D}$ is the oracle
specified in (6.13), and $R$ is the rotation operation.

######  Theorem 6.2.

Let $f(t,x)$ be a smooth function and $\vec{f}$ be a possibly unnormalized
vector

$\vec{f}=\sum_{k=0}^{m-1}\sum_{l_{1}=0}^{n-1}\cdots\sum_{l_{d}=0}^{n-1}f(kT/m,l_{1}/n_{1},\cdots,l_{d}/n_{d}){|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle}.$
(6.19)

Assume that $f$ satisfies periodic boundary condition for $x$ and
$\sup_{j,p}(\|\partial_{x_{j}}^{p}f\|_{\infty})^{1/p}<\pi n$. Then for any
$0<\epsilon<1$, $0<\delta<1$, there exists a quantum algorithm which, with
probability at least $(1-\delta)$, outputs an $\epsilon$-approximation of the
state proportional to

$\vec{\nabla
f}=\sum_{j=0}^{d-1}\sum_{k=0}^{m-1}\sum_{l_{1},\cdots,l_{d}=0}^{n-1}\partial_{x_{j}}f(kT/m,l_{1}/n,\cdots,l_{d}/n){|{j}\rangle}{|{k}\rangle}{|{l_{1}}\rangle}\cdots{|{l_{d}}\rangle},$
(6.20)

using queries to $\mathcal{A}(\epsilon/Q)$ and $O_{D}$ for
$\mathcal{O}(Q\log(1/\delta))$ times and additional $\mathcal{O}(d(\log
n)^{2})$ gates, where

$Q=\frac{4\sqrt{d}\|\vec{f}\|(\sup_{j,p}(\|\partial_{x_{j}}^{p}f\|_{\infty})^{1/p}+1)}{\|\vec{\nabla
f}\|}.$

We briefly compare our result with the standard classical approach of
computing gradient in terms of the dimension parameters, including $n$ and
$d$. On the one hand, the cost for computing gradient evaluated at all
discrete grid points typically scale $\mathcal{O}(dn^{d})$, since the sizes of
both the vector storing the information of the gradients and the matrices
related to finite difference and discrete Fourier transform scale
$\mathcal{O}(dn^{d})$. On the other hand, the corresponding scaling for our
quantum approach is subtler since it depends on the scaling of the quantity
$Q$. Notice that this quantity can be very large if the function $f$ is close
to a constant function. However, it can also be of $\mathcal{O}(1)$ if at
least one of the non-trivial Fourier components in $f$ is on the order of
$\Omega(1)$, since in this case $\|\vec{f}\|\sim\mathcal{O}(\sqrt{mn^{d}})$
and $\|\vec{\nabla f}\|\sim\Omega(\sqrt{dmn^{d}})$. In this scenario, the
overall complexity scales only polynomially in terms of $d$ and poly-
logarithmically in $n$, which achieves exponential speedup in terms of the
dimension parameters.

The cost of computing the ratio of the kinetic energy can be directly
estimated by combining Theorem 6.1 and Theorem 6.2. As discussed before, we
can get a quadratic speedup in terms of precision compared to the classical
Monte Carlo method. The result is summarized in the following Corollary.

######  Corollary 6.3.

Let $f(t,x)$ be a smooth function such that $f$ satisfies periodic boundary
condition for $x$, and
$\sup_{j,p}(\|\partial_{x_{j}}^{p}f\|_{\infty})^{1/p}<\pi n$. Assume that we
are given a unitary transform $(I-2P)$ for the projector $P$ associated with
the domain ${\mathscr{D}}$. Then for any $0<\epsilon<1$, $0<\delta<1$, there
exists a quantum algorithm which can output an approximation of (6.5) within
tolerated error $\epsilon$ and probability at least $(1-\delta)$, using
queries to $(I-2P)$ for $\mathcal{O}((1/\epsilon)\log(1/\delta))$ times,
queries to $\mathcal{A}(\epsilon/(4Q))$ and $O_{D}$ for
$\mathcal{O}(Q(\log(1/\delta))^{2}/\epsilon)$ times and additional
$\mathcal{O}(d(\log n)^{2}\log(1/\delta)/\epsilon)$ gates, where

$Q=\frac{4\sqrt{d}\|\vec{f}\|(\sup_{j,p}(\|\partial_{x_{j}}^{p}f\|_{\infty})^{1/p}+1)}{\|\vec{\nabla
f}\|}.$

###  6.3 History state and decay of kinetic energy

Unlike the existing work on quantum differential equation solvers that
typically output a final state encoding the solution at the final time, our
Carleman linearization algorithm and the input model assumed in this section
take a more general history state encoding the solutions at all time steps. In
this subsection, we briefly discuss how the general history state structure
may broaden the application of our algorithm.

One potential application is to study the kinetic energy curve. As discussed
in [85, 11], applied mathematicians are interested in the curve that describes
the decay of the kinetic energy of the gradient flow (2.4), and particularly
the time of reaching the equilibrium, _i.e._ , the time when the kinetic
energy almost stops changing. We define this time to be the equilibrium time
$t^{*}$, and it can be easily estimated by combining the history state and our
algorithm for computing gradients in Fig. 5. Specifically, we first run the
algorithm in Fig. 5 to get an approximation of the history state of gradients
of the solution as in (6.6), and then measure the time register
${|{k}\rangle}$ to obtain an integer $k$. Notice that after the equilibrium
time $t^{*}$ when the kinetic energy almost stops changing, the corresponding
partial derivatives are very close to 0. This implies that we have almost no
measurement outcomes after $t^{*}$. By repeating such a procedure and taking
the maximum of the measure outcomes $k_{\max}$, we can use $k_{\max}h$ to
estimate $t^{*}$ with high probability.

The history state structure can also allow us to overcome potential
exponential cost in the differential equation solvers caused by the decay of
the solution. To the extent of our knowledge, most of the existing quantum
differential equation solvers which output a final state [45, 46, 47, 1] scale
at least linearly in terms of the parameter
$\|u_{\text{in}}\|/\|u_{\text{out}}\|$, where $u_{\text{in}}$ and
$u_{\text{out}}$ denote the unnormalized solutions of the differential
equation at the initial and final time, respectively. Such a linear dependence
is typically caused by post-processing a quantum state obtained from solving
specific linear systems of equations to get the desired final state and may
introduce extra exponential time dependence if the solution of the
differential equations experiences rapid decay. A simple example is the
imaginary time evolution

$\frac{du}{ds}=-Hu$ (6.21)

where $H$ is a positive definite matrix. Here
$\|u(T)\|=\|\exp(-HT)u(0)\|\leq\exp(-\lambda T)\|u(0)\|$ with $\lambda$ being
the smallest eigenvalue of $H$, and thus $\|u(0)\|/\|u(T)\|\geq\exp(\lambda
T)$, leading to an extra exponentially large term in $T$. Another example is
the nonlinear ordinary differential equation with no constant term, which can
be explicitly written down as

$\frac{du}{ds}=F_{1}u+F_{2}u^{\otimes 2}.$ (6.22)

It is proved in [1] that, if $F_{1}$ only has negative eigenvalues and the
nonlinearity is relatively weak in certain sense (namely the parameter $R$
defined in this work is smaller than 1), then $\|u(t)\|$ decays exponentially
in terms of $t$, which also leads to an exponentially large
$\|u_{\text{in}}\|/\|u_{\text{out}}\|$ in $T$.

With a history state, we can run a “pre-diagnosis” to first identify whether
the final state is close to 0. The key observation is that if the state has
sufficiently decayed such that the solution is very close to 0, then the
success probability of getting the corresponding time step by measuring the
time register is exponentially small. In particular, assume that we are
interested in obtaining an approximation of the final solution $u(T)$ when
$\|u(T)\|$ is very small or the corresponding state $u(T)/\|u(T)\|$ when
$\|u(T)\|$ is reasonably away from 0. We can repeatedly prepare a history
state, measure the time register, and obtain the output $k$. If for all $k$’s,
we have $kh<T$, then the final solution is expected to be exponentially small
with high probability, so we can stop here and directly use 0 to be the
approximation of the solution. On the other hand, if there exists a $k$ such
that $kh=T$, then this is a reasonable indication that the quantity
$\|u_{\text{in}}\|/\|u_{\text{out}}\|$ is not quite large, and we can follow
the standard post-processing procedure to obtain the final state
$u(T)/\|u(T)\|$. The history state here helps us determine which scenario the
differential equation is in without exponential cost in the evolution time
$T$.

##  7 Discussion

We have presented an efficient quantum algorithm for the gradient flow
evolution of reaction-diffusion equations. We improve the Carleman
linearization under a condition $R_{D}<1$. It relaxes the previous condition
$R<1$ in [1] for high-dimensional systems of nonlinear differential equations.
Besides, we discussed estimating the mean square amplitude and ratios of the
kinetic energy of the gradient flow as potential applications.

This work raises several natural open problems. The first aspect is regarding
further improvement of our algorithm. Though our work focuses on improving the
convergence condition for the Carleman linearization, it is also interesting
to seek further improvement of the dependence on other parameters in the
complexity, such as the evolution time and the error tolerance. Another
related topic is to obtain meaningful classical outputs with super-polynomial
quantum speedups over the best-known classical algorithms. It is also an
interesting question to generalize our algorithms or design new quantum
algorithms dealing with other types of nonlinear differential equations. Our
analysis relies heavily on the Maximum Principle and some good regularity of
the solutions, which are essential properties of reaction-diffusion equations.
On a high level, the maximum principle controls the norm of the solution that
shares the same spirit of spectral norm preserving properties of Hamiltonian
simulation. It is thus an interesting direction to consider other norm-
controlled problems, such as the gradient flow structured on certain norm
preserving manifold, and see whether such confinement helps conquer the
nonlinearity at hand.

However, relaxing the regularity assumptions of the solutions seems to be a
fundamentally difficult problem. All linearization-based techniques require
the solutions to be well-posed and regular, which is not the case for
applications such as the conservation laws, fluid dynamics, and Hamilton-
Jacobi equation, where the solutions blow up in finite time. Therefore, other
approaches beyond the linearization framework are desired for such nonlinear
problems, requiring some new insights. We also point out that the quantum
Carleman linearization based approaches may suffer from an overhead sensitive
to the spatial grid refinement for partial differential equations. Though our
improved $\ell_{\infty}$ framed convergence criterion can concur the
sensitivity to grid refinements, this sensitivity can be reintroduced through
the $\ell_{2}$ norm dependence when implementing the quantum algorithm since
the $\ell_{2}$ norm of the solution discretized by finite difference still
grows as more grid points are used. In particular, it appears when getting the
solution from the huge Carleman state vector. It is our future work to make
the quantum algorithm fully insensitive to the spatial grid refinements,
probably by trying other spatial discretization which can preserve the
$\ell_{2}$ norm in grid refinements such as Fourier discretization.

In the application section, we use the amplitude estimate technique to obtain
classical information beyond the quantum state output and study the kinetic
energy distribution by preparing a state with derivative information. The
techniques we propose in this section do not rely on specific models or
differential equations. Thus it is interesting to find applications of our
technique to output classical information for other problems, such as phase
separation and transition, chemical reactions, and self-organized biological
patterns. We want to remark that all the classical outputs we study in this
work are in terms of ratios. While obtaining the absolute value rather than
the ratio seems to require accurate computation of the observable on the
entire domain, which might incur exponential overhead, it is still very
interesting to further study whether the absolute value of observables can be
approximated with only a relatively small overhead. Our algorithm for
preparing quantum states of derivatives requires the regularity of the
function, and we want to understand its performance for non-smooth functions
as well. It may also be of interest in some quantum optimization problems
which require gradient information to optimize the objective function and will
be our future work.

## Acknowledgments

We thank Andrew M. Childs and Lin Lin for valuable discussions. DA acknowledge
the hospitality of the Simons Institute for the Theory of Computing in
Berkeley. DA, DF and JW acknowledge the Challenge Institute for Quantum
Computation funded by NSF through grant number OMA-2016245, and the Department
of Energy under grant No. DE-SC0017867. DA acknowledges the support by the
Department of Defense through the Hartree Postdoctoral Fellowship at QuICS,
and the NSF under Grant No. DMS-1652330. DF is supported by the NSF grant
number DMS-2208416. JPL acknowledges support from the Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research, Quantum
Algorithms Teams and Accelerated Research in Quantum Computing programs, the
National Science Foundation award (CCF-1813814, DMS-2008568), and from the
National Science Foundation Quantum Information Science and Engineering
Network (QISE-NET) triplet award (DMR-1747426).

## References

  * [1] Jin-Peng Liu, Herman Øie Kolden, Hari K Krovi, Nuno F Loureiro, Konstantina Trivisa, and Andrew M Childs.  Efficient quantum algorithm for dissipative nonlinear differential equations.  Proceedings of the National Academy of Sciences, 118(35), 2021.  arXiv:2011.03185.
  * [2] Yuankai Lu and Dan Hu.  Optimisation of biological transport networks.  East Asian Journal on Applied Mathematics, 12(1):72–95, 2022.
  * [3] Dan Hu and David Cai.  Adaptation and optimization of biological transport networks.  Physical Review Letter, 111:138701, Sep 2013.
  * [4] Jan Haskovec, Peter Markowich, and Benoit Perthame.  Mathematical analysis of a PDE system for biological network formation.  Communications in Partial Differential Equations, 40(5):918–956, 2015.
  * [5] Jan Haskovec, Peter Markowich, Benoît Perthame, and Matthias Schlottbom.  Notes on a PDE system for biological network formation.  Nonlinear Analysis, 138:127–155, 2016.
  * [6] Giacomo Albi, Marco Artina, Massimo Foransier, and Peter A. Markowich.  Biological transportation networks: modeling and simulation.  Analysis and Applications, 14(1):185–206, 2016.
  * [7] Jan Haskovec, Lisa Maria Kreusser, and Peter Markowich.  ODE and PDE based modeling of biological transportation networks.  arXiv: 1805.08526, May 2018.
  * [8] Martin Burger, Jan Haskovec, Peter Markowich, and Helene Ranetbauer.  A mesoscopic model of biological transportation networks.  arXiv: 1806.00120, May 2018.
  * [9] Jan Haskovec, Lisa Maria Kreusser, and Peter Markowich.  Rigorous continuum limit for the discrete network formation problem.  arXiv: 1808.01526, August 2018.
  * [10] Giacomo Albi, Martin Burger, Jan Haskovec, Peter Markowich, and Matthias Schlottbom.  Continuum modeling of biological network formation, pages 1–48.  Model. Simul. Sci. Eng. Technol. Birkhäuser/Springer, Cham, 2017.
  * [11] Di Fang, Shi Jin, Peter Markowich, and Benoît Perthame.  Implicit and semi-implicit numerical schemes for the gradient flow of the formation of biological transport networks.  SMAI Journal of Computational Mathematics, 5:229–249, 2019.  https://hal.sorbonne-universite.fr/hal-01984371/file/caihu_v7_Small.pdf.
  * [12] Marcus R. Garvie.  Finite-difference schemes for reaction–diffusion equations modeling predator–prey interactions in m atlab.  Bulletin of mathematical biology, 69(3):931–956, 2007.
  * [13] Horst Malchow.  Spatiotemporal patterns in ecology and epidemiology: theory, models, and simulation.  Chapman and Hall/CRC, 2007.
  * [14] Sergei V. Petrovskii and Horst Malchow.  A minimal model of pattern formation in a prey-predator system.  Mathematical and Computer Modelling, 29(8):49–63, 1999.
  * [15] Julien Lefèvre and Jean-François Mangin.  A reaction-diffusion model of human brain development.  PLoS computational biology, 6(4):e1000749, 2010.
  * [16] Salman Habib, Carmen Molina-París, and Thomas S. Deisboeck.  Complex dynamics of tumors: modeling an emerging brain tumor system with coupled reaction–diffusion equations.  Physica A: Statistical Mechanics and its Applications, 327(3-4):501–524, 2003.
  * [17] James D. Murray.  Mathematical biology II: spatial models and biomedical applications, volume 3.  Springer New York, 2001.
  * [18] James D. Murray.  Mathematical biology I : An introduction. interdisciplinary applied mathematics.  Mathematical Biology, Springer, 2002.
  * [19] Stéphane Genieys, Vitaly Volpert, and Pierre Auger.  Pattern and waves for a model in population dynamics with nonlocal consumption of resources.  Mathematical Modelling of Natural Phenomena, 1(1):63–80, 2006.
  * [20] Hans Meinhardt.  Models of biological pattern formation.  New York, 118, 1982.
  * [21] Ido Golding, Yonathan Kozlovsky, Inon Cohen, and Eshel Ben-Jacob.  Studies of bacterial branching growth using reaction–diffusion models for colonial development.  Physica A: Statistical Mechanics and its Applications, 260(3-4):510–554, 1998.
  * [22] Masayasu Mimura, Hideo Sakaguchi, and Mitsugu Matsushita.  Reaction–diffusion modelling of bacterial colony patterns.  Physica A: Statistical Mechanics and its Applications, 282(1-2):283–303, 2000.
  * [23] Henri Berestycki, Basil Nicolaenko, and Bruno Scheurer.  Traveling wave solutions to combustion models and their singular limits.  SIAM Journal on Mathematical Analysis, 16(6):1207–1242, 1985.
  * [24] I. A. Zeldovich, Grigory Isaakovich Barenblatt, V. B. Librovich, and G. M. Makhviladze.  Mathematical theory of combustion and explosions, 1985.
  * [25] Thierry Poinsot and Denis Veynante.  Theoretical and numerical combustion.  RT Edwards, Inc., 2005.
  * [26] Benoıt Perthame.  Growth, reaction, movement and diffusion from biology.  Lecture Notes, University Paris, 6, 2012.
  * [27] Shawn Means, Alexander J. Smith, Jason Shepherd, John Shadid, John Fowler, Richard J. H. Wojcikiewicz, Tomas Mazel, Gregory D. Smith, and Bridget S. Wilson.  Reaction diffusion modeling of calcium dynamics with realistic er geometry.  Biophysical journal, 91(2):537–557, 2006.
  * [28] Andrea L. Bertozzi and Arjuna Flenner.  Diffuse interface models on graphs for classification of high dimensional data.  Multiscale Modeling & Simulation, 10(3):1090–1118, 2012.  https://www.math.ucla.edu/ bertozzi/papers/MMS2012-final.pdf.
  * [29] Andrea L. Bertozzi and Arjuna Flenner.  Diffuse interface models on graphs for classification of high dimensional data.  SIAM Review, 58(2):293–328, 2016.  https://www.math.ucla.edu/ bertozzi/papers/SIGEST2016.pdf.
  * [30] Ekaterina Merkurjev, Tijana Kostic, and Andrea L. Bertozzi.  An MBO scheme on graphs for classification and image processing.  SIAM Journal on Imaging Sciences, 6(4):1903–1930, 2013.  https://www.math.ucla.edu/ bertozzi/papers/MBOgraphsFinal.pdf.
  * [31] Andrea L. Bertozzi, Selim Esedoglu, and Alan Gillette.  Inpainting of binary images using the Cahn–Hilliard equation.  IEEE Transactions on image processing, 16(1):285–291, 2006.
  * [32] Julia A. Dobrosotskaya and Andrea L. Bertozzi.  A wavelet-Laplace variational technique for image deconvolution and inpainting.  IEEE Transactions on Image Processing, 17(5):657–663, 2008.  https://www.math.ucla.edu/ bertozzi/papers/julia-short-07.pdf.
  * [33] Selim Esedoglu and Riccardo March.  Segmentation with depth but without detecting junctions.  Journal of Mathematical Imaging and Vision, 18(1):7–15, 2003.  http://dept.math.lsa.umich.edu/ esedoglu/Papers_Preprints/nms_esedoglu_march.pdf.
  * [34] Selim Esedog, Yen-Hsi Richard Tsai, et al.  Threshold dynamics for the piecewise constant Mumford–Shah functional.  Journal of Computational Physics, 211(1):367–384, 2006.  http://dept.math.lsa.umich.edu/ esedoglu/Papers_Preprints/esedoglu_tsai.pdf.
  * [35] Richard Bellman.  Dynamic programming, 1957.  Princeton University Press, New Jersey.
  * [36] Andris Ambainis.  Variable time amplitude amplification and quantum algorithms for linear algebra problems.  In 29th Symposium on Theoretical Aspects of Computer Science, volume 14, pages 636–647. LIPIcs, 2012.  arXiv:1010.4458.
  * [37] Dong An and Lin Lin.  Quantum linear system solver based on time-optimal adiabatic quantum computing and quantum approximate optimization algorithm.  ACM Transactions on Quantum Computing, 3(2):1–28, Jun 2022.
  * [38] Andrew M. Childs, Robin Kothari, and Rolando D. Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.  arXiv:1511.02306.
  * [39] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.  arXiv:1806.01838.
  * [40] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15):150502, 2009.  arXiv:0811.3171.
  * [41] Lin Lin and Yu Tong.  Optimal quantum eigenstate filtering with application to solving quantum linear systems.  Quantum, 4:361, 2020.  arXiv:1910.14596.
  * [42] Yiğit Subaşı, Rolando D. Somma, and Davide Orsucci.  Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing.  Physical Review Letters, 122(6):060504, 2019.  arXiv:1805.10549.
  * [43] Yu Tong, Dong An, Nathan Wiebe, and Lin Lin.  Fast inversion, preconditioned quantum linear system solvers, fast Green’s-function computation, and fast evaluation of matrix functions.  Physical Review A, 104(3):032422, 2021.  arXiv:2008.13295.
  * [44] Pedro C. S. Costa, Dong An, Yuval R. Sanders, Yuan Su, Ryan Babbush, and Dominic W. Berry.  Optimal scaling quantum linear systems solver via discrete adiabatic theorem, 2021.  arXiv:2111.08152.
  * [45] Dominic W. Berry.  High-order quantum algorithm for solving linear differential equations.  Journal of Physics A: Mathematical and Theoretical, 47(10):105301, 2014.  arXiv:1010.2745.
  * [46] Dominic W. Berry, Andrew M. Childs, Aaron Ostrander, and Guoming Wang.  Quantum algorithm for linear differential equations with exponentially improved dependence on precision.  Communications in Mathematical Physics, 356(3):1057–1081, 2017\.  arXiv:1701.03684.
  * [47] Andrew M. Childs and Jin-Peng Liu.  Quantum spectral methods for differential equations.  Communications in Mathematical Physics, 375:1427–1457, 2020.  arXiv:1901.00961.
  * [48] Di Fang, Lin Lin, and Yu Tong.  Time-marching based quantum solvers for time-dependent linear differential equations.  Quantum, 7:955, March 2023.
  * [49] B. David Clader, Bryan C. Jacobs, and Chad R. Sprouse.  Preconditioned quantum linear system algorithm.  Physical Review Letters, 110(25):250504, 2013.  arXiv:1301.2340.
  * [50] Yudong Cao, Anargyros Papageorgiou, Iasonas Petras, Joseph Traub, and Sabre Kais.  Quantum algorithm and circuit design solving the Poisson equation.  New Journal of Physics, 15(1):013021, 2013.  arXiv:1207.2485.
  * [51] Ashley Montanaro and Sam Pallister.  Quantum algorithms and the finite element method.  Physical Review A, 93(3):032324, 2016.  arXiv:1512.05903.
  * [52] Pedro C. S. Costa, Stephen Jordan, and Aaron Ostrander.  Quantum algorithm for simulating the wave equation.  Physical Review A, 99(1):012323, 2019.  arXiv:1711.05394.
  * [53] Andrew M Childs, Jin-Peng Liu, and Aaron Ostrander.  High-precision quantum algorithms for partial differential equations.  Quantum, 5:574, 2021.  arXiv:2002.07868.
  * [54] Alexander Engel, Graeme Smith, and Scott E. Parker.  Quantum algorithm for the Vlasov equation.  Physical Review A, 100(6):062315, 2019.  arXiv:1907.09418.
  * [55] Noah Linden, Ashley Montanaro, and Changpeng Shao.  Quantum vs. classical algorithms for solving the heat equation.  arXiv:2004.06516.
  * [56] Sarah K. Leyton and Tobias J. Osborne.  A quantum algorithm to solve nonlinear differential equations, 2008.  arXiv:0812.4423.
  * [57] Daniel S. Abrams and Seth Lloyd.  Nonlinear quantum mechanics implies polynomial-time solution for NP-complete and #P problems.  Physical Review Letters, 81(18):3992, 1998.  arXiv:quant-ph/9801041.
  * [58] Scott Aaronson.  NP-complete problems and physical reality.  ACM SIGACT News, 36(1):30–52, 2005.  arXiv:quant-ph/0502072.
  * [59] Andrew M. Childs and Joshua Young.  Optimal state discrimination and unstructured search in nonlinear quantum mechanics.  Physical Review A, 93(2):022314, 2016.  arXiv:1507.06334.
  * [60] Torsten Carleman.  Application de la théorie des équations intégrales linéaires aux systèmes d’équations différentielles non linéaires.  Acta Mathematica, 59(1):63–87, 1932.
  * [61] Krzysztof Kowalski and Willi-Hans Steeb.  Nonlinear Dynamical Systems and Carleman Linearization.  World Scientific, 1991.
  * [62] Marcelo Forets and Amaury Pouly.  Explicit error bounds for Carleman linearization, 2017.  arXiv:1711.02552.
  * [63] Hari Krovi.  Improved quantum algorithms for linear and nonlinear differential equations.  Quantum, 7:913, February 2023.
  * [64] Ilya Y. Dodin and Edward A. Startsev.  On applications of quantum computing to plasma simulations.  Physics of Plasmas, 28(9):092101, 2021.  arXiv:2005.14369.
  * [65] Ilon Joseph.  Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics.  Physical Review Research, 2(4):043102, 2020.  arXiv:2003.09980.
  * [66] Seth Lloyd, Giacomo De Palma, Can Gokler, Bobak Kiani, Zi-Wen Liu, Milad Marvian, Felix Tennie, and Tim Palmer.  Quantum algorithm for nonlinear differential equations, 2020.  arXiv:2011.06571.
  * [67] Alexander Engel, Graeme Smith, and Scott E. Parker.  Linear embedding of nonlinear dynamical systems and prospects for efficient quantum algorithms.  Physics of Plasmas, 28(6):062305, 2021.  arXiv:2012.06681.
  * [68] Cesare Tronci and Ilon Joseph.  Koopman wavefunctions and Clebsch variables in Vlasov-Maxwell kinetic theory, 2021.  arXiv:2105.00294.
  * [69] Shi Jin and Nana Liu.  Quantum algorithms for computing observables of nonlinear partial differential equations, 2022.  arXiv:2202.07834.
  * [70] Ilya Y. Dodin and Edward A. Startsev.  Quantum computation of nonlinear maps, 2021.  arXiv:2105.07317.
  * [71] Cheng Xue, Yu-Chun Wu, and Guo-Ping Guo.  Quantum homotopy perturbation method for nonlinear dissipative ordinary differential equations.  New Journal of Physics, 23(12):123035, 2021.  arXiv:2111.07486.
  * [72] Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  Contemporary Mathematics, 305:53–74, 2002.  arXiv:quant-ph/0005055.
  * [73] Samuel M. Allen and John W. Cahn.  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening.  Acta Metallurgica, 27(6):1085–1095, 1979.
  * [74] R. A. FISHER.  The wave of advance of advantageous genes.  Annals of Eugenics, 7(4):355–369, 1937.
  * [75] Yves Van Gennip, Andrea L Bertozzi, et al.  $\Gamma$-convergence of graph Ginzburg-Landau functionals.  Advances in Differential Equations, 17(11/12):1115–1180, 2012.  arXiv:1204.5220.
  * [76] Yves Van Gennip, Nestor Guillen, Braxton Osting, and Andrea L. Bertozzi.  Mean curvature, threshold dynamics, and phase field theory on finite graphs.  Milan Journal of Mathematics, 82(1):3–65, 2014.  https://link.springer.com/article/10.1007/s00032-014-0216-8.
  * [77] Xiyang Luo and Andrea L. Bertozzi.  Convergence of the graph Allen–Cahn scheme.  Journal of Statistical Physics, 167(3-4):934–958, 2017.  https://www.math.ucla.edu/ bertozzi/papers/JSTATPHYS2017-final.pdf.
  * [78] Bobak T. Kiani, Giacomo De Palma, Dirk Englund, William Kaminsky, Milad Marvian, and Seth Lloyd.  Quantum advantage for differential equation analysis, 2020.  arXiv:2010.15776.
  * [79] F. Brauer and J.A. Nohel.  The Qualitative Theory of Ordinary Differential Equations: An Introduction.  Dover Books on Mathematics. Dover Publications, 2012.
  * [80] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  Communications in Mathematical Physics, 270:359–371, 2007.  arXiv:quant-ph/0508139.
  * [81] Pedro C.S. Costa, Dong An, Yuval R. Sanders, Yuan Su, Ryan Babbush, and Dominic W. Berry.  Optimal scaling quantum linear-systems solver via discrete adiabatic theorem.  PRX Quantum, 3:040303, Oct 2022.
  * [82] Dong An, Jin-Peng Liu, Daochen Wang, and Qi Zhao.  A theory of quantum differential equation solvers: limitations and fast-forwarding, 2023.
  * [83] R. L. Burden, J. D. Faires, and A. C. Reynolds.  Numerical analysis.  Brooks Cole, 2000.
  * [84] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2000.
  * [85] Benoît Perthame.  Parabolic equations in biology.  Lecture Notes on Mathematical Modelling in the Life Sciences. Springer, Cham, 2015.  Growth, reaction, movement and diffusion.
  * [86] Mark R. Jerrum, Leslie G. Valiant, and Vijay V. Vazirani.  Random generation of combinatorial structures from a uniform distribution.  Theoretical computer science, 43:169–188, 1986.

##  Appendix A Proof of Lemma 2.1

We now discuss the proof of the a priori estimate of the solution. Before
that, we introduce the comparison principle lemma for the discrete reaction-
diffusion equation (2.7), which implies the a priori estimate as a direct
consequence.

######  Lemma A.1 (Comparison principle).

Assume $f(u)\in C^{\infty}(\mathbb{R})$ and $U(t)=[U_{1},\cdots,U_{n}]^{T}$,
$V(t)=[V_{1},\cdots,V_{n}]^{T}$ are continuous functions that satisfy

$\frac{\mathrm{d}U_{j}}{\mathrm{d}t}-D\left(\Delta_{h}U\right)_{j}-f(U_{j})\leqslant\frac{\mathrm{d}V_{j}}{\mathrm{d}t}-D\left(\Delta_{h}V\right)_{j}-f(V_{j}),$

for $t\in(0,T]$ and all the multi-indices $j\in\mathcal{I}$. Furthermore,
$U_{j}(0)\leqslant V_{j}(0)$ for all $j\in\mathcal{I}$ and $U_{j}(t)\leqslant
V_{j}(t)$ for $j\in\mathcal{B}$ and $t\in[0,T]$. Then it holds that

$U_{j}(t)\leqslant V_{j}(t)$

for all $j$ and all time $t\in(0,T]$.

We remark that Lemma A.1 immediately implies Lemma 2.1, because both
$U_{j}(t)=\gamma_{1}$ and $U_{j}(t)=\gamma_{2}$ all $j$ and $t$ are solutions
to (2.7) and hence comparison principle can be applied to the solution of
interests and these constant-valued equilibrium solutions, which yields the
desired result.

###### Proof.

Our proof can be split into the following three steps: in the first two steps,
we consider a linear operator of $U$ given by
$\frac{\mathrm{d}U_{j}}{\mathrm{d}t}-D\left(\Delta_{h}U\right)_{j}+\tilde{C}_{j}U_{j}$
for some vector $\tilde{C}=\left(\tilde{C_{j}}\right)$, and show that a
maximum principle result for this linear operator; and in the last step, we
prove the comparison principle for the nonlinear problem as considered in
(A.1). It is worth pointing out that although a linear problem is considered
at first, there is no linearization procedure introducing any extra error
here.

First, we claim that if

$\frac{\mathrm{d}V_{j}}{\mathrm{d}t}-D\left(\Delta_{h}V\right)_{j}+\tilde{C}_{j}V_{j}<0,$
(A.1)

for some $\tilde{C}$ with all entries positive, then it holds that

$\displaystyle\max_{t\in[0,T],j\in\mathcal{I}}V_{j}(t)=\displaystyle\max_{(j,t)\in(\mathcal{B}\times(0,T])\cup(\mathcal{I}\times\\{0\\})}V^{+}_{j}(t)$

where $\mathcal{B}$ denotes the boundary indices and
$V^{+}_{j}:=\max(V_{j},0)$.

Suppose the claim does not hold, then there exists some time $t_{0}\in(0,T]$
and $j_{0}\not\in\mathcal{B}$ such that $V_{j}(t)$ attains the positive
maximum value at $(t_{0},j)$. Note that the definition (2.8) of $\Delta_{h}$
can be written in the following equivalent form

$\displaystyle\Delta_{h}v_{j_{1},\cdots,j_{n}}:=$
$\displaystyle\frac{v_{j_{1}+1,\cdots,j_{n}}-2v_{j_{1},\cdots,j_{n}}+v_{j_{1}-1,\cdots,j_{n}}}{h^{2}}+\frac{v_{j_{1},j_{2}+1,\cdots,j_{n}}-2v_{j_{1},j_{2},\cdots,j_{n}}+v_{j_{1},j_{2}-1\cdots,j_{n}}}{h^{2}}$
$\displaystyle+\cdots+\frac{v_{j_{1},\cdots,j_{n}+1}-2v_{j_{1},\cdots,j_{n}}+v_{j_{1},\cdots,j_{n}-1}}{h^{2}},$
(A.2)

and hence one has $-\left(\Delta_{h}V\right)_{j_{0}}\geq 0$. Meanwhile, if
$t_{0}\in(0,T)$, then

$\frac{\mathrm{d}V_{j_{0}}}{\mathrm{d}t}(t_{0})=0;$

otherwise, $t_{0}=T$, and one has

$\frac{\mathrm{d}V_{j_{0}}}{\mathrm{d}t}(t_{0})\geq 0.$

Therefore, one has

$\frac{\mathrm{d}V_{j}}{\mathrm{d}t}-D\left(\Delta_{h}V\right)_{j}+\tilde{C}_{j}V_{j}\geq
0,$

which is a contradiction. This completes the proof of the claim.

In the second step, we relax the condition of the claim from the following two
angles: the condition on $C$ is changed from positive to bounded from below;
and the equality is allowed. To be precise, we shall show that for $C$
satisfying that $C_{j}\geq c_{\rm min}$ for all $j\in\mathcal{I}$, if

$\frac{\mathrm{d}V_{j}}{\mathrm{d}t}-D\left(\Delta_{h}V\right)_{j}+\tilde{C}_{j}V_{j}\leqslant
0,$ (A.3)

then

$\displaystyle\max_{t\in[0,T],j\in\mathcal{I}}V_{j}(t)=\displaystyle\max_{(j,t)\in(\mathcal{B}\times(0,T])\cup(\mathcal{I}\times\\{0\\})}V^{+}_{j}(t).$
(A.4)

To prove this, we define $\tilde{V}(t)$ by the change of variable
$\tilde{V}(t)=e^{c_{\rm min}t}V(t)-\delta t$. A straightforward calculation
yields

$\frac{\mathrm{d}\tilde{V}_{j}}{\mathrm{d}t}=c_{\rm
min}\tilde{V}_{j}+e^{c_{\rm
min}t}\frac{\mathrm{d}V_{j}}{\mathrm{d}t}-\delta\leqslant
D\left(\Delta_{h}\tilde{V}\right)_{j}-\left(\tilde{C}_{j}-c_{\rm
min}\right)\tilde{V}_{j}-\delta.$

Thus $\tilde{V}_{j}$ satisfies (A.1). Letting $\delta\to 0_{+}$ yields the
desired result.

The last part of the proof is to show the comparison principle. Let $W:=U-V$.
One has $W_{j}(0)\leqslant 0$ for $j\in\mathcal{I}$ and $W_{j}(t)\leqslant 0$
for all $t\in(0,T]$ and $j\in\mathcal{B}$ so that the right-hand-side of (A.4)
is non-positive. Moreover, $W$ satisfies

$\frac{\mathrm{d}W_{j}}{\mathrm{d}t}-D\left(\Delta_{h}W\right)_{j}-f(U_{j})+f(V_{j})\leqslant
0,$

and by mean value theorem $f(U_{j})-f(V_{j})=f^{\prime}(\xi_{j})$ with
$\xi_{j}$ in between $U_{j}$ and $V_{j}$, we arrive at

$\frac{\mathrm{d}W_{j}}{\mathrm{d}t}-D\left(\Delta_{h}W\right)_{j}-f^{\prime}(\xi_{j})W_{j}\leqslant
0.$

Applying the result of the second step yields

$W_{j}(t)\leqslant 0,$

for all $j\in\mathcal{I}$ and $t\in[0,T]$, which completes the proof of this
lemma. ∎

##  Appendix B Matrix inequality

######  Lemma B.1 (Discrete maximal principle).

Let $D_{h}$ be the discrete $1$-dimensional Laplacian operator with
homogeneous Dirichlet boundary conditions, then

$\left\lVert e^{tD_{h}}\right\rVert_{\infty}\leq 1,\quad\forall~{}t>0.$ (B.1)

###### Proof.

In this proof, we will apply the results,

$e^{tD_{h}}=\lim_{k\to\infty}\left(1+\frac{t\Delta_{h}}{k}\right)^{k}.$ (B.2)

First, we restrict $t\leq t_{0}:=\frac{1}{(n+1)^{2}}$. Then when $k$ is larger
enough, i.e. $k\geq 2t_{0}(n+1)^{2}$, one has

$\left\lVert
1+\frac{tD_{h}}{k}\right\rVert_{\infty}\leq\left(1-2\frac{t}{k}(n+1)^{2}\right)+\frac{t}{k}(n+1)^{2}+\frac{t}{k}(n+1)^{2}\leq
1.$ (B.3)

Hence,

$\left\lVert\left(1+\frac{tD_{h}}{k}\right)^{k}\right\rVert_{\infty}\leq\left\lVert
1+\frac{tD_{h}}{k}\right\rVert_{\infty}^{k}\leq 1.$

By continuity, one can get that $\left\lVert
e^{tD_{h}}\right\rVert_{\infty}\leq 1$. As for $t>t_{0}$, we choose positive
integer $l$ such that $t/l\leq t_{0}$ and get

$\left\lVert
e^{tD_{h}}\right\rVert_{\infty}=\left\lVert\left(e^{tD_{h}/l}\right)^{l}\right\rVert_{\infty}\leq\left\lVert
e^{tD_{h}/l}\right\rVert_{\infty}^{l}\leq 1.$

In this way, we obtain the desired results. ∎

######  Lemma B.2.

Let $D_{h}^{\text{per}}$ be the discrete $1$-dimensional Laplacian operator
with periodic boundary conditions, then

$\left\lVert
e^{tD_{h}^{\text{per}}}\right\rVert_{\infty}=1,\quad\forall~{}t>0.$ (B.4)

###### Proof.

Using the same argument in Lemma B.1, we can prove that

$\left\lVert e^{tD_{h}^{\text{per}}}\right\rVert_{\infty}\leq
1,\quad\forall~{}t>0.$

Now we only need to prove the inequality in the opposite direction. Denote $v$
to be the eigenvector of $D_{h}^{\text{per}}$ associated with eigenvalue $0$.
Then we know that $e^{tD_{h}^{\text{per}}}v=v$, which implies $\left\lVert
e^{tD_{h}^{\text{per}}}\right\rVert_{\infty}\geq 1$. In this way, we obtain
that $\left\lVert e^{tD_{h}^{\text{per}}}\right\rVert_{\infty}=1$. ∎

######  Lemma B.3.

Let $D_{h}^{\text{Dir}}$ be the discrete $1$-dimensional Laplacian operator
with homogeneous Dirichlet boundary conditions, and denote
$\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}$ to be the eigenvalues of
$D_{h}^{\text{Dir}}$. Then

$\left\lVert
e^{tD_{h}^{\text{Dir}}}\right\rVert_{\infty}\leq\left(\frac{4}{\pi}+\frac{2e^{\mu_{2}t}}{e^{-2\mu_{1}t/\pi}-1}\right)e^{\mu_{1}t},\quad\forall~{}t>0.$
(B.5)

###### Proof.

We use the decomposition of $D_{h}^{\text{Dir}}$:

$D_{h}^{\text{Dir}}=U\text{diag}(\mu_{1},...,\mu_{n})U^{\top}\text{ with
}U_{kl}=\sqrt{\frac{2}{n+1}}\sin\left(\frac{kl\pi}{n+1}\right).$

Here $\mu_{k}=-4(n+1)^{2}\left(\sin\left(\frac{k\pi}{2n+2}\right)\right)^{2}$.
Besides, we know that for any $3\leq k\leq n$,

$\begin{split}\mu_{k}-\mu_{2}&=-4(n+1)^{2}\left(\sin\left(\frac{k\pi}{2n+2}\right)\right)^{2}+4(n+1)^{2}\left(\sin\left(\frac{\pi}{n+1}\right)\right)^{2}\\\
&=2(n+1)^{2}\left(\cos\left(\frac{k\pi}{n+1}\right)-\cos\left(\frac{2\pi}{n+1}\right)\right)\\\
&=-4(n+1)^{2}\sin\left(\frac{(k-2)\pi}{2n+2}\right)\sin\left(\frac{(k+2)\pi}{2n+2}\right)\\\
&\leq-4(n+1)^{2}\sin\left(\frac{\pi}{2n+2}\right)\sin\left(\frac{(k+1)\pi}{2n+2}\right)\\\
&\leq\mu_{1}\frac{\sin\left(\frac{(k+1)\pi}{2n+2}\right)}{\sin\left(\frac{\pi}{2n+2}\right)}\leq\mu_{1}\frac{\frac{k+1}{n+1}}{\frac{\pi}{2n+2}}=\frac{2\mu_{1}}{\pi}(k+1).\end{split}$
(B.6)

In the last step, we use the inequality $\frac{2}{x}\leq\sin(x)\leq x$ for
$x\in[0,\frac{\pi}{2}]$. Note that

$\left\lVert e^{tD_{h}^{\text{Dir}}}\right\rVert_{\infty}=\max_{1\leq k\leq
n}\sum_{l=1}^{n}\left\lvert\sum_{j=1}^{n}U_{kj}e^{\mu_{j}t}U_{lj}\right\rvert.$
(B.7)

For any $1\leq k\leq n$, we have

$\begin{split}&\sum_{l=1}^{n}\left\lvert\sum_{j=1}^{n}U_{kj}e^{\mu_{j}t}U_{lj}\right\rvert\\\
&=\frac{2}{n+1}\sum_{l=1}^{n}\left\lvert\sum_{j=1}^{n}e^{\mu_{j}t}\sin\left(\frac{kj\pi}{n+1}\right)\sin\left(\frac{lj\pi}{n+1}\right)\right\rvert\\\
&\leq\frac{2}{n+1}\sum_{l=1}^{n}\sum_{j=1}^{n}e^{\mu_{j}t}\left\lvert\sin\left(\frac{kj\pi}{n+1}\right)\sin\left(\frac{lj\pi}{n+1}\right)\right\rvert\\\
&\leq\frac{2n}{n+1}\sum_{j=2}^{n}e^{\mu_{j}t}+\frac{2}{n+1}e^{\mu_{1}t}\sin\left(\frac{k\pi}{n+1}\right)\sum_{l=1}^{n}\sin\left(\frac{l\pi}{n+1}\right).\end{split}$
(B.8)

We notice that for any $k$,

$\begin{split}\frac{2}{n+1}e^{\mu_{1}t}\sin\left(\frac{k\pi}{n+1}\right)\sum_{l=1}^{n}\sin\left(\frac{l\pi}{n+1}\right)&\leq\frac{2}{n+1}e^{\mu_{1}t}\sin\left(\frac{k\pi}{n+1}\right)\frac{1}{\tan\left(\frac{\pi}{2n+2}\right)}\\\
&\leq\frac{4}{\pi}e^{\mu_{1}t}\sin\left(\frac{k\pi}{n+1}\right)\leq\frac{4}{\pi}e^{\mu_{1}t},\end{split}$
(B.9)

where we use the equality
$\sum_{l=1}^{n}\sin\left(\frac{l\pi}{n+1}\right)=\frac{1}{\tan\left(\frac{\pi}{2n+2}\right)}$
and $\tan(x)\geq x$ for $x\in[0,\frac{\pi}{2})$.

We also notice that

$\begin{split}\frac{2n}{n+1}\sum_{j=2}^{n}e^{\mu_{j}t}&=\frac{2n}{n+1}e^{\mu_{2}t}\sum_{j=2}^{n}e^{(\mu_{j}-\mu_{2})t}\leq\frac{2n}{n+1}e^{\mu_{1}t}\sum_{j=2}^{n}e^{2(k+1)\mu_{1}t/\pi}\\\
&\leq
2e^{\mu_{2}t}\frac{e^{6\mu_{1}t/\pi}\left(1-e^{2(n-1)\mu_{1}t/\pi}\right)}{1-e^{2\mu_{1}t/\pi}}\leq\frac{2}{e^{-2\mu_{1}t/\pi}-1}e^{(\mu_{1}+\mu_{2})t}.\end{split}$
(B.10)

In this way, we obtain

$\left\lVert
e^{tD_{h}^{\text{Dir}}}\right\rVert_{\infty}\leq\left(\frac{4}{\pi}+\frac{2e^{\mu_{2}t}}{e^{-2\mu_{1}t/\pi}-1}\right)e^{\mu_{1}t}.$
(B.11)

∎

Combining Lemma B.1 and Lemma B.3, we obtain the following results.

######  Lemma B.4.

Let $\mu_{1}$ be the largest eigenvalue of $1$-dimensional
$D_{h}^{\text{Dir}}$ with homogeneous Dirichlet boundary conditions. Given
$t\in{\mathbb{R}}^{+}$, for arbitrary $\mu>\mu_{1}$, we have the decay
estimate

$\left\lVert
e^{tD_{h}^{\text{Dir}}}\right\rVert_{\infty}\leq\begin{cases}1,&0<t<\frac{\ln(3)}{2(\mu-\mu_{1})},\\\
e^{t\mu},&t\geq\frac{\ln(3)}{2(\mu-\mu_{1})}.\end{cases}$ (B.12)

###### Proof.

The case when $t<\frac{\ln(3)}{2(\mu-\mu_{1})}$ directly follows from Lemma
B.1. When $t\geq\frac{\ln(3)}{2(\mu-\mu_{1})}$, recall that

$\mu_{2}=-4(n+1)^{2}\left(\sin\left(\frac{2\pi}{2n+2}\right)\right)^{2}=4\mu_{1}\left(\cos\left(\frac{\pi}{2n+2}\right)\right)^{2}\leq
3\mu_{1}.$ (B.13)

Then we have