# Two quantum algorithms for solving the one-dimensional advection-diffusion
equation

Julia Ingelmann  Sachin S. Bharadwaj  Philipp Pfeffer  Katepalli R.
Sreenivasan  Jörg Schumacher

###### Abstract

Two quantum algorithms are presented for the numerical solution of a linear
one-dimensional advection-diffusion equation with periodic boundary
conditions. Their accuracy and performance with increasing qubit number are
compared point-by-point with each other. Specifically, we solve the linear
partial differential equation with a Quantum Linear Systems Algorithms (QLSA)
based on the Harrow–Hassidim–Lloyd method and a Variational Quantum Algorithm
(VQA), for resolutions that can be encoded using up to 6 qubits, which
corresponds to $N=64$ grid points on the unit interval. Both algorithms are of
hybrid nature, i.e., they involve a combination of classical and quantum
computing building blocks. The QLSA and VQA are solved as ideal statevector
simulations using the in-house solver QFlowS and open-access Qiskit software,
respectively. We discuss several aspects of both algorithms which are crucial
for a successful performance in both cases. These are the sizes of an
additional quantum register for the quantum phase estimation for the QLSA and
the choice of the algorithm of the minimization of the cost function for the
VQA. The latter algorithm is also implemented in the noisy Qiskit framework
including measurement and decoherence circuit noise. We reflect the current
limitations and suggest some possible routes of future research for the
numerical simulation of classical fluid flows on a quantum computer.

††journal: Computers & Fluids

[label1]organization=Institute of Thermodynamics and Fluid Mechanics,
addressline=Technische Universität Ilmenau, P.O.Box 100565, city=Ilmenau,
postcode=D-98684, country=Germany

[label2]organization=Tandon School of Engineering, addressline=New York
University, city=New York City, postcode=11201, state=NY, country=USA

[label3]organization=Courant Institute of Mathematical Sciences,
addressline=New York University, city=New York City, postcode=10012, state=NY,
country=USA

[label4]organization=Department of Physics, addressline=New York University,
city=New York City, postcode=10012, state=NY, country=USA

[label5]organization=Center for Space Science, addressline=New York University
Abu Dhabi, city=Abu Dhabi, postcode=129188, country=United Arab Emirates

##  1 Introduction

Quantum computing has the potential to open new ways to classify, generate,
and process data [1, 2] thus changing paradigms in many application fields,
such as material science, renewable energy technology, and finance. The reason
for the expected advantage over classical algorithms is the physical
foundation of quantum computing. Quantum algorithms are capable of encoding
information in superposition states and of combining several such quantum
states into tensorial product states which span high-dimensional spaces. They
can perform unitary transformations (quantum gates) on these product states in
parallel rather than on individual bits, as done in classical computers. In
this way, $n$ qubits—the smallest units of quantum information—span a
$2^{n}$-dimensional Hilbert space. This parallelism is tightly connected to
the possibility of entangling qubits, representing inseparable correlations
between qubits, which is absent in classical bit registers [3]. Already these
two properties suggest faster solutions of problems with high computational
complexity, as has been demonstrated for operations such as prime number
factorization [4], data search [5], and data sampling [6]; see ref. [7] for a
discussion. Still open is the question of whether similar advantages survive
the application of quantum algorithms to solutions of nonlinear ordinary and
partial differential equations.

Fluid dynamics comprises many applications with high computational effort, for
instance the modeling of flows over complex objects such as airplanes and the
Direct Numerical Simulation (DNS) of turbulent flows [8] that resolves all
physically relevant flow scales from the system size down to those dominated
by viscous and diffusive effects. The nonlinear partial differential equations
(PDEs) relevant to us are the Navier-Stokes equations for the flow, and
(simultaneously) the advection-diffusion equation for the transport of the
scalar field such as a substance concentration and temperature. The numerical
effort to resolve all these spatial scales increases as $N^{3}$ for the three-
dimensional case, which varies at least as fast as $Re^{9/4}$. Here, $N$ the
number of mesh points along one spatial direction and $Re$ is the flow
Reynolds number that quantifies the vigor of the fluid turbulence. In many
technological applications for which the geometry of the flow domain is
complex, one requires in addition adaptive refinements of the computational
meshes. Consequently, resource limits are reached quickly, even on the largest
state-of-the-art supercomputers. The present solution to this problem is to
model the small-scale part, e.g., in the form of Reynolds-averaged Navier-
Stokes equation models or large eddy simulations.

A further possible solution might be the transformation of classical fluid
flow problems on a quantum computer to make use of the parallelism that
originates from the quantum mechanical foundations. As one example, a single
velocity component of a DNS of homogeneous isotropic turbulence in a periodic
box with $N^{3}=8192^{3}\approx 5.5\times 10^{11}$ grid points [9, 10] could
be encoded theoretically in less than 40 qubits, which should be eminently
doable since the biggest quantum chip contains 433 qubits. This motivates our
present work.

Several approaches have been suggested in the past years to study fluid flows
on quantum computers. They include a transformation into a quantum computing-
inspired tensor product framework with an effective mapping of the excited
degrees of freedom of a three-dimensional turbulent flow [11] or the mapping
of specific classical flow problems to a Schrödinger-type quantum dynamics
[12, 13, 14]. They include also a surrogate modeling of thermally driven flows
within quantum machine learning frameworks, such as hybrid quantum-classical
reservoir computing [15, 16]. Implementations of mostly one-dimensional flow
problems on a quantum computer in the form of pure or hybrid algorithms
comprise quantum linear systems algorithms for steady pipe Poiseuille, plane
Couette flows and Burgers equation [17, 18, 19], quantum amplitude estimation
for one-dimensional gas dynamics [20], Variational Quantum Algorithms for the
one-dimensional nonlinear Burgers equation [21, 22], advection-diffusion
problems [23, 24, 25], and quantum lattice Boltzmann methods [26, 27]. See
also ref. [28] for a recent perspective.

The potential of quantum computing algorithms for solving advection-diffusion
problems has been investigated in different ways recently. One approach is the
decomposition of the PDE into finite differences such that the resulting
system of linear equations can be solved. For sparse linear equation systems,
the Harrow-Hassidim-Lloyd (HHL) algorithm can provide a exponential speed up
in comparison to classical computation [29] under certain caveats [30, 31]. A
further approach are variational methods. Different versions, such as the
variational quantum imaginary time evolution [25], the Variational Quantum
Linear Solver (VQLS) [23], or the Variational Quantum Algorithm (VQA) [21, 32,
33], have been used, even for two-dimensional problems, such as the heat
equation [33].

The present work compares these two popular hybrid quantum-classical
algorithms for a standard benchmark problem in fluid mechanics, which is the
one-dimensional advection-diffusion equation with a constant advection
velocity $U$, described by a linear partial differential equation. To this
end, we will compare one-to-one a hybrid quantum-classical Variational Quantum
Algorithm (VQA) with a Quantum Linear Systems Algorithm (QLSA). The purpose of
the present study is to explore the scalability of both algorithms up to mesh
grids which will be encoded in registers consisting of $n\leq 6$ qubits giving
resolutions of $N\leq 64$ grid points. Furthermore, we identify the
bottlenecks that exist in both cases for some of their main building blocks:
for the VQA scheme, this turns out to be the classical optimization algorithm
for the minimization of the cost function; for the QLSA it is the quantum
phase estimation routine—an approximate method to find eigenvalues of a
unitary matrix. Several classical optimization algorithms are therefore
compared in the VQA case. Here, we also investigate the role of the depth of
the parametric quantum circuit on the performance of the VQA algorithm and
report the impact of measurements for data readout on the overall performance.
In case of QLSA, the underlying hybrid algorithm presented here (which in
itself preserves the speed-up [18]) is customized carefully for the advection-
diffusion problem. We analyse the algorithm’s performance after prescribing
specific strategies for accurate eigenvalue estimation. We also evaluate its
dependence on the number of qubits, preconditioning and measurement. To keep
our manuscript self-contained and accessible to the fluid dynamics community,
we provide compact introductions to quantum computing as well as the two
algorithms. Finally, we critically assess both algorithms for this simple
fluid mechanical problem and thereby discuss possible limitations of quantum
algorithms for (nonlinear) fluid flow problems in one red(and higher-
dimensional) cases.

The article is structured as follows. First, the analytical solution for the
one-dimensional advection-diffusion equation is obtained as the basis for the
comparison of the quantum algorithms (Sec. 2). Second, the numerical scheme of
the finite differences approach is given for forward and backward Euler
stepping, which is the groundwork for the quantum algorithms considered here
(Sec. 3). Then, the quantum algorithms are introduced in detail (Sec. 4). The
comparison of both quantum algorithms is shown in Sec. 5 on aspects such as
the time evolution of the concentration profiles, dependence on the number of
qubits, dependence on parameter $T_{0}$ and the realisation on Noisy
Intermediate Scale Quantum (NISQ) devices. The results are summarized and
discussed in Sec. 6.

##  2 One-dimensional advection-diffusion equation

We demonstrate and compare the performance of the two quantum computing
algorithms considered here for the advection-diffusion equation given by

$\displaystyle\partial_{t}c=D\nabla^{2}c-{\bm{u}}\cdot\nabla c,$ (1)

where $c({\bm{x}},t)$ is the concentration field of the solvent, $D$ is the
diffusion constant and ${\bm{u}}({\bm{x}},t)$ is the velocity vector field
that advects the solvent. This equation describes the transport of the
solvent, such as a dye or a cloud of tracer particles subject to diffusion and
advection with the velocity field. In this paper, we consider the simplest
case of a one-dimensional linear equation, which is given by

$\displaystyle\partial_{t}c(x,t)=D\partial^{2}_{x}c(x,t)-U\partial_{x}c(x,t).$
(2)

Here, the advection velocity $U$ in the $x$-direction is taken as a constant.
The problem is discretized in space and time. For the spatial discretization,
the interval $x\in[-L,L]$ is divided into $N$ segments of width $\Delta
x=2L/N$. The time evolution is also discretized uniformly, such that
$t=m\tau$, where $\tau$ is the time step. For the analytical solution, the
wave-like ansatz $c(x,t)=\text{exp}\left(\omega t+i\lambda x\right)$ is chosen
such that $\omega=-D\lambda^{2}-iU\lambda$ follows from eq. (2). Hence, the
concentration profile takes the form of

$c(x,t)=\left[a\cos{\lambda(x-Ut)}+b\sin{\lambda(x-Ut)}\right]\text{exp}(-D\lambda^{2}t)\,.$
(3)

Periodic boundary conditions are imposed, such that $c(x=0,t)=c(x=N,t)$.
Consequently, the wavenumber $\lambda=k\pi/L$ with $k\in\mathbb{N}$. Thus
there follows the general solution to the problem in the form of a series
expansion

$\displaystyle c(x,t)$
$\displaystyle=\sum_{k=0}^{\infty}\left[a_{k}\cos{\left(\frac{k\pi}{L}(x-Ut)\right)}+b_{k}\sin{\left(\frac{k\pi}{L}(x-Ut)\right)}\right]$
$\displaystyle\times\text{exp}\left(-D\left(\frac{k\pi}{L}\right)^{2}t\right)$
(4)

As initial condition the delta function is applied such that
$c(x,0)=\delta(x)$. The delta function is standard and defined to be
$\delta(x)=0$ for $x\neq 0$ and $\int_{-\infty}^{\infty}\delta(x)\text{d}x=1$.
The initial condition specifies the expansion coefficients in the general
solution as

$\displaystyle
c(x,0)=\frac{a_{0}}{2}+\sum_{k=1}^{\infty}a_{k}\cos\left(\frac{k\pi}{L}x\right)+b_{k}\sin\left(\frac{k\pi}{L}x\right).$
(5)

The Fourier coefficients are given by

$\displaystyle a_{0}$
$\displaystyle=\frac{1}{L}\int_{-L}^{L}\delta(x)\text{d}x=\frac{1}{L},$ (6)
$\displaystyle a_{k}$
$\displaystyle=\frac{1}{L}\int_{-L}^{L}\delta(x)\cos\left(\frac{k\pi}{L}x\right)\text{d}x=\frac{1}{L}\text{
and}$ (7) $\displaystyle b_{k}$
$\displaystyle=\frac{1}{L}\int_{-L}^{L}\delta(x)\sin\left(\frac{k\pi}{L}x\right)\text{d}x=0,$
(8)

such that the initial condition can written as

$\displaystyle
c(x,0)=\frac{1}{2L}+\sum_{k=1}^{\infty}\frac{1}{L}\cos\left(\frac{k\pi}{L}x\right).$
(9)

Considering these coefficients, the analytical solution can be found to be

$\displaystyle c(x,t)$
$\displaystyle=\frac{1}{2L}+\sum_{k=1}^{\infty}\left[\frac{1}{L}\cos{\left(\frac{k\pi}{L}(x-Ut)\right)}\right]\text{exp}\left(-D\left(\frac{k\pi}{L}\right)^{2}t\right)\,.$
(10)

Equation (10) describes a Gauss-shaped pulse that diffuses while moving to the
right given that $U>0$. For the following, lengths are measured in units of
the interval length $2L$. Times can be expressed in units of either the
advection time $\tau_{a}=2L/U$ or the diffusive time $\tau_{d}=(2L)^{2}/D$. If
not stated otherwise, we will use $\tau_{a}$.

##  3 Finite difference methods with Euler time stepping

The numerical solution of the advection-diffusion equation (2) can be obtained
by a finite difference method (FDM). In the simplest case, these are Euler
methods, either an explicit forward or an implicit backward Euler time step
method. For this method, the partial differential equation is approximated by
a system of algebraic discretization equations. Furthermore, the problem is
discretized in space and time uniformly, such that $x_{i}=x_{0}+i\Delta x$ and
$t_{m}=t_{0}+m\tau$ with $x_{0}=0$ and $t_{0}=0$. Indices $i=0,\dots,N-1$ and
$m=0,\dots,M$. When the forward difference in time and the centered difference
in space is taken, one gets the forward in time and centered in space (FCTS)
method. It is of 1st order accuracy in time, of 2nd order accuracy in space,
and given by

$\displaystyle\frac{c_{i}^{m+1}-c_{i}^{m}}{\tau}=D\frac{c_{i+1}^{m}-2c_{i}^{m}+c_{i-1}^{m}}{(\Delta
x)^{2}}-U\frac{c_{i+1}^{m}-c_{i-1}^{m}}{2\Delta x}\,,$ (11)

and thus

$c_{i}^{m+1}=\left(\frac{D\tau}{(\Delta x)^{2}}-\frac{U\tau}{2\Delta
x}\right)c_{i+1}^{m}+\left(1-2\frac{D\tau}{(\Delta
x)^{2}}\right)c_{i}^{m}+\left(\frac{D\tau}{(\Delta
x)^{2}}+\frac{U\tau}{2\Delta x}\right)c_{i-1}^{m}\,.$

We define the following abbreviations

$r=\frac{D\tau}{(\Delta x)^{2}}\quad\mbox{and}\quad s=\frac{U\tau}{2\Delta
x}\,,$ (12)

where $s$ is the parameter of the convective part and $r$ is the stability
parameter. For this explicit scheme $r\leq 1/2$ should hold. The scheme can be
expressed as a system of linear equations via

$\displaystyle A{\bm{c}}^{m}$ $\displaystyle={\bm{c}}^{m+1}$ (13)
$\displaystyle\text{with}\ \ A$
$\displaystyle=\left[\begin{array}[]{cccccc}1-2r&r-s&0&\ldots&0&s+r\\\
s+r&1-2r&r-s&&&0\\\ 0&s+r&1-2r&r-s&&0\\\
\vdots&&\ddots&\ddots&\ddots&\vdots\\\ 0&&&s+r&1-2r&r-s\\\
r-s&0&\ldots&0&s+r&1-2r\\\ \end{array}\right].$ (20)

In case of the implicit backward Euler scheme (BTCS), the system of linear
equation follows such that the matrix $A$ has to be inverted to find the
desired solution, since this method imposes the expression

$\displaystyle\frac{c_{i}^{m+1}-c_{i}^{m}}{\tau}=D\frac{c_{i+1}^{m+1}-2c_{i}^{m+1}+c_{i-1}^{m+1}}{(\Delta
x)^{2}}-U\frac{c_{i+1}^{m+1}-c_{i-1}^{m+1}}{2\Delta x}\,,$ (21)

which can be reformulated to

$c_{i}^{m}=\left(-\frac{D\tau}{(\Delta x)^{2}}+\frac{U\tau}{2\Delta
x}\right)c_{i+1}^{m+1}+\left(1+2\frac{D\tau}{(\Delta
x)^{2}}\right)c_{i}^{m+1}+\left(-\frac{D\tau}{(\Delta
x)^{2}}-\frac{U\tau}{2\Delta x}\right)c_{i-1}^{m+1}\,.$

In other words, the scheme can be expressed as a system of linear equations
via

$\displaystyle A{\bm{c}}^{m+1}$ $\displaystyle={\bm{c}}^{m}$ (22)
$\displaystyle\text{with}\ \ A$
$\displaystyle=\left[\begin{array}[]{cccccc}1+2r&-r+s&0&\cdots&0&-r-s\\\
-r-s&1+2r&-r+s&&&0\\\ 0&-r-s&1+2r&-r+s&&0\\\
\vdots&&\ddots&\ddots&\ddots&\vdots\\\ 0&&&-r-s&1+2r&-r+s\\\
-r+s&0&\cdots&0&-r-s&1+2r\\\ \end{array}\right]$ (29)

The comparison of the analytical solution (ANA) with those of FCTS and BCTS is
shown in Fig. 1. The comparison is made via the mean squared error (MSE)
defined as

$\displaystyle{\rm
MSE}(t_{m})=\frac{1}{N}\sum_{i=0}^{N-1}\left[c_{i}^{\text{ANA}}(t_{m})-c_{i}^{\text{FDM}}(t_{m})\right]^{2}\,,$
(30)

where $c_{i}=c(x_{i},t)$. In Figs. 1(a) and 1(b), it can be seen that the
numerical methods approximate the analytical solution sufficiently well. This
could be different when nonlinear equations have to be solved with VQA.

Figure 1: Time evolution of the concentration profiles of the analytical
solution (ANA) and the results of the classical numerical methods, namely the
explicit method (FTCS), the implicit method (BTCS) and the midpoint method
(MP) for $N=32$, $D=1$, $U=10$. Panels (a) and (b) compare the concentration
profiles of the methods at two time instants. The corresponding mean squared
error (MSE) of the results of the classical numerical methods to the
analytical solution is shown in panel (c).

##  4 Quantum algorithms

This section describes both quantum algorithms, namely the VQA and the QLSA.
The quantum part of the VQA is implemented in the quantum simulation
environment Qiskit [34]. The QLSA is done with QFlowS, a C++ based simulation
package [18]. For the direct comparison of both algorithms, an ideal
statevector simulation will be used. In the following, we will briefly
introduce the basics of both quantum algorithms. The building block of both
algorithms are the qubits, the smallest information units in a quantum
algorithm. While a single classical bit can take two discrete values only,
namely $\\{0,1\\}$, a qubit is a superposition of the two basis states of the
Hilbert space $\mathbb{C}^{2}$

$|q_{1}\rangle=\alpha_{1}|0\rangle+\alpha_{2}|1\rangle=\alpha_{1}\left(\begin{array}[]{c}1\\\
0\\\ \end{array}\right)+\alpha_{2}\left(\begin{array}[]{c}0\\\ 1\\\
\end{array}\right)\,,$ (31)

with $\alpha_{1},\alpha_{2}\in\mathbb{C}$ and
$\|q_{1}\|_{2}=\sqrt{|\alpha_{1}|^{2}+|\alpha_{2}|^{2}}=1$ and basis vectors
$|0\rangle$ and $|1\rangle$ in Dirac’s notation [3]. It can be combined into
an $n$-qubit system, also denoted as an $n$-qubit quantum register, by
successive tensor products of qubits. An unentangled two-qubit state vector is
the tensor product of two single-qubit vectors,

$|q_{1}\rangle\otimes|q_{1}^{\prime}\rangle\in\mathbb{C}^{2}\otimes\mathbb{C}^{2}\,.$
(32)

The basis of this tensor product space is given by 4 vectors, usually
formulated in integer or binary bit-string notation:
$|j_{1}\rangle=|0\rangle\otimes|0\rangle=|00\rangle$,
$|j_{2}\rangle=|01\rangle$, $|j_{3}\rangle=|10\rangle$, and
$|j_{4}\rangle=|11\rangle$. The $n$-qubit quantum state $|c(t)\rangle$ at a
time $t$ is consequently defined in a $2^{n}$-dimensional tensor product
Hilbert space ${\cal H}=(\mathbb{C}^{2})^{\otimes n}$ and given by

$\lvert c(t)\rangle=\sum\limits_{k=1}^{2^{n}}c_{k}(t)\lvert j_{k}\rangle\hskip
10.00002pt\text{with}\hskip
10.00002pt\sum\limits_{k=1}^{2^{n}}|c_{k}(t)|^{2}=1\,.$ (33)

In other words, when connecting this formalism to the present flow problem,
the discretization of the concentration profile $c(x,t)$ on $N=2^{n}$ grid
points at time $t$ is obtained by an $n$-qubit quantum state vector. In eq.
(33), the quantum state vector is normalized to 1. Technically, the square
magnitude of each coefficient represents the probability of measuring the
respective basis state. Thereby, they naturally have to sum up to 1. For a
classical concentration profile this does not have to be the case; see
subsection 4.2. The time evolution of the state vector in quantum algorithms
is established by unitary transformations or operators $\hat{U}$ with
$\hat{U}^{-1}=\hat{U}^{\dagger}$, realized on a quantum computer by a sequence
of quantum gates. These gates can be viewed as rotation operators on quantum
state vectors that can also generate entanglement between qubit states [3].

Note that the accuracy of the quantum algorithms depends on the accuracy of
the numerical input data described in Sec. 3. For the comparison of the
quantum algorithms, the analytical solution, which we derived in Sec. 2, is
also considered.

###  4.1 Quantum Linear Systems Algorithm (QLSA)

(a)

(b)

Figure 2: (a) A sketch of the quantum circuit depicting the QLSA-1
architecture. The horizontal lines in the diagram stand either for individual
qubits or collections of several qubits, so-called quantum registers. Boxes
across qubits stand for parametric circuits, collections of multiple unitary
transformations on a single or multiple qubits. These are from left to right
Quantum State Preparation (QSP), Quantum Phase Estimation (QPE), Conditional
Rotation and inverse Quantum Phase Estimation (IQPE). (b) A sketch depicting
the hybrid workflow of the QLSA algorithm. The classical preparation of the
matrix elements and pre-conditioning is done classically. Then, on the quantum
simulator the quantum states are prepared and processed with quantum state
preparation (QSP) and QLSA, respectively. The solutions are measured or post
processed at the quantum level and read into the classical machine. The
residue and convergence checks are made classically and the solutions are
post-processed when the iterations have finished or converged.

Quantum algorithms which solve a linear system of equations of the form,
$A\mathbf{x}=\mathbf{b}$, belong to the category of Quantum Linear Systems
Algorithm (QLSA). All such algorithms [29, 35, 36] (excluding variational
methods, which will be described subsequently in Sec. 4.2) can be broadly
categorized into two approaches which compute a quantum-numerical
approximation to $A^{-1}\mathbf{b}$ (BTCS) or $A\mathbf{x}$ (FTCS). The
approach presented here, which we call QLSA-1, is a modified version of the
original HHL algorithm [29]. Here, we compute the eigenvalues ($\sigma_{j}$)
of the matrix $A$ and thereby approximate the solution $A^{-1}\mathbf{b}$. The
central computational issue here is to identify the eigenvalues of $A$ and the
following evaluation of their inverse. An alternative algorithm we call QLSA-2
proceeds by approximating the action of the matrix $A$ (or $A^{-1}$) as purely
a matrix-vector multiplication operation, implemented by decomposing the
matrix into a Linear Combination of Unitary (LCU) quantum gates [35], acting
on a suitably prepared quantum state. The central goal in that case is to find
the best unitary basis to produce a probabilistic implementation of the
matrix. Both methods have been implemented on QFlowS in [18] to solve laminar
Poiseuille and Couette flows. The solution $c(t)$ can be obtained either
iteratively at every time-step or by one-shot QLSA algorithms that would offer
higher quantum advantage [37, 18]. However, the latter strategy can be
computationally expensive to simulate large system sizes over long integration
times. For our present purposes, we present results for QLSA-1 using the
former approach.

The QLSA-1 algorithm is implemented as a full gate-level circuit simulation
with at most single qubit or (double) controlled NOT gates [3] to integrate
eq. (2) using the BTCS method. The outline of the algorithm’s work flow and
its circuit is shown in Fig. 2. It comprises the steps or quantum sub-routines
briefly outlined below (whose details can be found in [18]). The herein flow
problem has the matrix A of eq. (29), which is not Hermitian if advection is
present, namely $U\neq 0\ (\Rightarrow s\neq 0)$. Since the algorithm admits
only Hermitian matrices, the matrix $A$ is first extended to an Hermitian
classically as

$\tilde{A}=\begin{pmatrix}0&A\\\ A^{{\dagger}}&0\end{pmatrix}.$ (34)

The implementation involves the following steps.

Step 1 - Quantum State Preparation (QSP): The concentration field at every
time step $m$ is loaded onto an $(n+1)$-qubit ($=n_{b}$ from here on) state
proportional
$\tilde{\mathbf{c}}^{m}=[\mathbf{c}^{m};0]\,\|\mathbf{c}^{m}\|_{2}^{-1}$ to
make it compatible with eq. (34) (and therefore one expects the solution state
in the form $\mathbf{x}=[0;\mathbf{c}^{m+1}]$). As will be described shortly,
the algorithm also requires an additional $n_{q}+1$ ancillary qubits. The
latter are helper qubits or in short ancillas. All these qubits that are
initially set to basis state $|0\rangle$, are then initialized using either
the functional form type state preparation or the sparse-state preparation
oracle $\hat{U}_{\rm QSP}$ (see Sec. 3 of SI Appendix in [18]),

$|\psi_{\text{STEP1}}\rangle=\ \hat{U}_{\rm QSP}|0\rangle_{n_{b}+n_{q}+1}=\
|\tilde{c}(t)\rangle_{n_{b}}\otimes|0\rangle_{n_{q}}\otimes|0\rangle.$ (35)

Step 2 - Quantum Phase Estimation (QPE): Given a linear operator $U$, if
$e^{i\pi\sigma}$ is an eigenvalue, the QPE essentially estimates the phase
angle $\sigma$ as a binary representation $n_{q}$-bit
$|\varphi_{1}\varphi_{2}\cdots\varphi_{n_{q}}\rangle$,
$\forall\varphi_{k}\in\\{0,1\\}$. Using this algorithm, an $n_{q}$-bit binary
approximation to the eigenvalues $\tilde{\sigma}_{j}$ of $\tilde{A}$ is
computed. For this purpose, we first rescale the matrix by a suitable value so
that its eigenvalues lie in a range that is optimal for the algorithm’s
performance [18, 36], and, in addition, is a subset of $[-0.5,0.5]$, to obtain
the matrix $\bar{A}$. To now invoke QPE, this matrix is exponentiated as
$e^{i\bar{A}T_{0}}$ to form a linear unitary operator, where $T_{0}$ is the
Hamiltonian simulation time [29, 3]. This parameter can be regarded as a
scaling parameter that rescales the eigenvalues of $\bar{A}$ such that the
eigenvalues $\bar{\sigma}_{j}$ can be represented nearly exactly using an
$n_{q}$-bit binary state with minimal truncation error. The matrix $\bar{A}$
can be expanded in the eigenbasis $|v_{j}\rangle\langle v_{j}|$ such that

$e^{i\bar{A}T_{0}}:=\sum\limits_{j=1}^{2^{n_{b}}}e^{i\bar{\sigma}_{j}T_{0}}|{v_{j}}\rangle\langle{v_{j}}|.$
(36)

Following this, the QPE then produces the state proportional to

$|\psi_{\text{STEP2}}\rangle=\sum\limits_{j=1}^{2^{n_{b}}}\hat{\mathbf{c}}_{j}^{m}|v_{j}\rangle_{n_{b}}\otimes|\bar{\sigma}_{j0}\rangle_{n_{q}}\otimes|0\rangle,$
(37)

where $\bar{\sigma}_{j0}=\bar{\sigma}_{j}T_{0}$ are the binary represented
eigenvalues of A rescaled by $T_{0}$ while $\hat{\mathbf{c}}_{j}^{m}$ are the
coefficients of the normalized $\tilde{\mathbf{c}}^{m}$ generated by rotating
into the basis of A’s eigenvectors $|v_{j}\rangle$.

Step 3 - Conditional Rotation: Here we apply a relative rotation operator on
the last ancilla qubit, conditioned on $\bar{\sigma}_{j}$ to compute the
inverse $1/\bar{\sigma}$,

$|\psi_{\text{STEP3}}\rangle=\sum\limits_{j=1}^{2^{n_{b}}}\hat{\mathbf{c}}_{j}^{m}|v_{j}\rangle_{n_{b}}\otimes|\bar{\sigma}_{j0}\rangle_{n_{q}}\otimes\Bigg{(}\sqrt{1-\frac{K^{2}}{\bar{\sigma}_{j0}^{2}}}|0\rangle+\frac{K}{\bar{\sigma}_{j0}}|1\rangle\Bigg{)}$
(38)

where $K$ is a suitably chosen normalization constant.
Step 4: Finally, we perform the inverse QPE (IQPE) operation to reset $n_{q}$
to $|0\rangle$, and follow it up by a measurement of the last ancilla qubit in
the computational basis, producing a state proportional to

${|\mathbf{c}^{m+1}\rangle}\sim
R\times\sqrt{\dfrac{1}{\sum\limits_{j=0}^{2^{n_{b}-1}}|b_{j}|^{2}/|\bar{\sigma}_{j0}|^{2}}}\,\sum\limits_{j=1}^{2^{n_{b}}}\frac{\hat{\textbf{c}}^{m}_{j}}{\bar{\sigma}_{j0}}|v_{j}\rangle_{n_{b}}\otimes|0\rangle_{n_{q}}\otimes|1\rangle\,,$
(39)

where $R$ is the corresponding rescaling constant to extract the solution
appropriately. The solution can now either be read into classical formats by
sampling every component of the wavefunction from performing multiple runs of
the circuit (so-called shots), or the state can also be post-processed within
the quantum simulator to estimate linear and nonlinear functions of the
concentration field as shown in [18]. This would help conserve a certain
degree of quantum advantage. In any case, the results are then finally
assimilated in the classical device for post processing and output.

###  4.2 Variational Quantum Algorithm (VQA)

The variational quantum algorithm (VQA) is a hybrid quantum-classical
algorithm where a parameterized cost function $C$ is minimized by an optimizer
[38]. While the cost function is evaluated by a quantum circuit composed of
single- and two-qubit gates, the optimization is performed classically. This
defines the hybrid character of the algorithm. The general principle of the
VQA is shown in Fig. 3. The initial parameter set
$(\lambda_{0},{\bm{\lambda}})_{\text{init}}$, which consists of the
normalization factor $\lambda_{0}$ and the angles of the single-qubit unitary
rotation gates ${\bm{\lambda}}=(\lambda_{1},\lambda_{2},\dots)$, is the input
to the algorithm. Then, a cost function $C(\lambda_{0},{\bm{\lambda}})$, which
is parameterized with $(\lambda_{0},{\bm{\lambda}})$, is evaluated on a
quantum device. For our approach, a classical device adds results of multiple
quantum circuits together to generate the final cost. The minimum of the cost
function corresponds to the solution of the considered problem. This solution
is modeled by a quantum ansatz function $\hat{U}({\bm{\lambda}})$ which is
initialized with the parameter set ${\bm{\lambda}}$. Measurements of the
quantum circuits evaluate the costs. These costs are minimized with an
classical optimizer. The optimal parameter set
$(\lambda_{0}^{*},{\bm{\lambda}}^{*})$ initializes the ansatz function such
that the solution of the given problem can be observed [38].

Figure 3: Principal sketch of the Variational Quantum Algorithm (VQA) which
interacts between a quantum computational cost evaluation and a classical
supervised parameter update. Here, $(\lambda_{0},{\bm{\lambda}})$ is the
parameter vector which is optimized, $C(\lambda_{0},{\bm{\lambda}})$ is the
cost function, $\hat{U}({\bm{\lambda}})$ is the quantum ansatz function and
$(\lambda_{0}^{*},{\bm{\lambda}}^{*})$ is the optimal parameter vector, which
corresponds with the minimum of the cost function, hence with the solution of
the considered problem. The output of the quantum simulator results from
measurements of the quantum Hadamard test circuit which is indicated in the
quantum simulator block.

Note that VQA can also straightforwardly be applied to nonlinear equations
[21, 22]. In order to derive the cost function for the present advection-
diffusion problem, the discrete concentration profile is transformed in vector
notation such that eq. (2) can be written explicitly as in Euler type FTCS
methods by

$\displaystyle|c(t+\tau)\rangle=(\mathds{1}+\tau\hat{O})|c(t)\rangle,$ (40)

with the linear operator $\hat{O}=D\partial^{2}_{x}-U\partial_{x}$ and the
identity operator $\mathds{1}$. Then, the corresponding cost function $C$ can
be found as

$\displaystyle
C(|c(t+\tau)\rangle)=\||c(t+\tau)\rangle-(\mathds{1}+\tau\hat{O})|c(t)\rangle\|_{2}^{2},$
(41)

where the minimum of the cost function $C$ corresponds to the solution
$|c(t+\tau)\rangle$. Following Lubasch et al. [21], we define

$\displaystyle|c(t+\tau)\rangle$
$\displaystyle=\lambda_{0}|\psi({\bm{\lambda}})\rangle=\lambda_{0}\hat{U}({\bm{\lambda}})|0\rangle,$
(42) $\displaystyle|c(t)\rangle$
$\displaystyle=\tilde{\lambda}_{0}|\tilde{\psi}\rangle=\tilde{\lambda}_{0}\hat{\tilde{U}}|0\rangle,$
(43)

where ${\bm{\lambda}}$ is the parameter vector which initializes the quantum
ansatz for the solution $|c(t+\tau)\rangle$. The quantum states $|\psi\rangle$
are normalized, such that $\|\psi\|_{2}^{2}=1$, while for the concentration
profile holds

$\int_{-L}^{L}c(x,t)dx=\mbox{const}\quad\Rightarrow\sum_{k=1}^{2^{n}}c_{k}=1\,.$
(44)

The constant is 1 for the present case due to $c(x,0)=\delta(x)$. In order to
fulfill both constraints, normalization parameters, $\lambda_{0}$ and
$\tilde{\lambda}_{0}$ are introduced. In the present work we will rescale our
solution to an $L_{2}$-norm of 1 to be directly comparable to the QLSA case.
Thus eq. (41) results in

$\displaystyle
C(\lambda_{0},{\bm{\lambda}})=\|\lambda_{0}|\psi({\bm{\lambda}})\rangle-\tilde{\lambda}_{0}(\mathds{1}+\tau\hat{O})|\tilde{\psi}\rangle\|_{2}^{2}.$
(45)

The norm is evaluated by the scalar product and gives

$\displaystyle C(\lambda_{0},{\bm{\lambda}})$
$\displaystyle=\lambda_{0}^{2}\underbrace{\langle\psi({\bm{\lambda}})|\psi({\bm{\lambda}})\rangle}_{=1}-2\lambda_{0}\tilde{\lambda}_{0}\langle\psi({\bm{\lambda}})|(\mathds{1}+\tau\hat{O})\tilde{\psi}\rangle$
(46)
$\displaystyle+\underbrace{\tilde{\lambda}_{0}^{2}\langle\tilde{\psi}|(\mathds{1}+\tau\hat{O})^{\dagger}(\mathds{1}+\tau\hat{O})|\tilde{\psi}\rangle}_{=\text{const.}}$

The last term is constant for each time step because it depends only on
$|\tilde{\psi}\rangle$ and $\tilde{\lambda}_{0}^{2}$, which are fixed from the
previous time step. A further decomposition of the scalar product leads to

$\displaystyle C(\lambda_{0},{\bm{\lambda}})$
$\displaystyle=\lambda_{0}^{2}-2\lambda_{0}\tilde{\lambda}_{0}\left[\langle\psi({\bm{\lambda}})|\tilde{\psi}\rangle+\tau\langle\psi({\bm{\lambda}})|\hat{O}\tilde{\psi}\rangle\right]$
$\displaystyle+\tilde{\lambda}_{0}^{2}\left[1+2\tau\langle\tilde{\psi}|\hat{O}\tilde{\psi}\rangle+\tau^{2}\langle\tilde{\psi}|\hat{O}^{\dagger}\hat{O}\tilde{\psi}\rangle\right]\,.$
(47)

The linear operator $\hat{O}$ consists of 2 terms, the diffusion term and the
advection term. Rather than implementing these terms directly, a 2nd-order
finite difference discretization of both operators is used, which is in line
with the discussion in Sec. 3. For this, the unitary shift operators
$\psi_{i+1}=\hat{S}_{+}\psi_{i}$ and $\psi_{i-1}=\hat{S}_{-}\psi_{i}$ are
defined; see A for details. Note that thereby, the periodic boundary
conditions are imposed. With the definition of these shift operators and
$2l/\Delta x=N$ one gets

$\displaystyle\hat{O}$
$\displaystyle=DN^{2}(\hat{S}_{+}-2\mathds{1}+\hat{S}_{-})-U\frac{N}{2}(\hat{S}_{+}-\hat{S}_{-})$
$\displaystyle=(\underbrace{DN^{2}-U\frac{N}{2}}_{=\alpha})\hat{S}_{+}-\underbrace{2DN^{2}}_{=\beta}\mathds{1}+(\underbrace{DN^{2}+U\frac{N}{2}}_{=\gamma})\hat{S}_{-}\,.$
(48)

Consequently, the cost function can be written as a further decomposition of
the scalar products, leading to the following summary of the cost function:

$\displaystyle C(\lambda_{0},{\bm{\lambda}})=\lambda_{0}^{2}$
$\displaystyle-2\lambda_{0}\tilde{\lambda}_{0}\bigg{[}(1-\tau\beta)\langle\psi({\bm{\lambda}})|\tilde{\psi}\rangle+\tau\alpha\langle\psi({\bm{\lambda}})|\hat{S}_{+}\tilde{\psi}\rangle+\tau\gamma\langle\psi({\bm{\lambda}})|\hat{S}_{-}\tilde{\psi}\rangle\bigg{]}$
$\displaystyle+\tilde{\lambda}_{0}^{2}\bigg{[}1+2\tau\alpha\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle-2\tau\beta+2\tau\gamma\underbrace{\langle\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle}_{=\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}\bigg{]}$
$\displaystyle+\tilde{\lambda}_{0}^{2}\tau^{2}\bigg{[}\alpha^{2}\underbrace{\langle\hat{S}_{+}\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}_{=1}-\alpha\beta\underbrace{\langle\hat{S}_{+}\tilde{\psi}|\tilde{\psi}\rangle}_{=\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}+\alpha\gamma\underbrace{\langle\hat{S}_{+}\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle}_{=\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle}\bigg{]}$
$\displaystyle+\tilde{\lambda}_{0}^{2}\tau^{2}\bigg{[}-\beta\alpha\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle+\beta^{2}\underbrace{\langle\tilde{\psi}|\tilde{\psi}\rangle}_{=\mathds{1}}-\beta\gamma\underbrace{\langle\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle}_{=\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}\bigg{]}$
$\displaystyle+\tilde{\lambda}_{0}^{2}\tau^{2}\bigg{[}\gamma\alpha\underbrace{\langle\hat{S}_{-}\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}_{=\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle}-\gamma\beta\underbrace{\langle\hat{S}_{-}\tilde{\psi}|\tilde{\psi}\rangle}_{=\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}+\gamma^{2}\underbrace{\langle\hat{S}_{-}\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle}_{=1}\bigg{]}\,.$
(49)

We use that $\hat{S}_{+}=\hat{S}_{-}^{-1}=\hat{S}_{-}^{\dagger}$ and
$\hat{S}_{-}=\hat{S}_{+}^{\dagger}$. This leads to the final cost function

$\displaystyle C(\lambda_{0},{\bm{\lambda}})$
$\displaystyle=\lambda_{0}^{2}-2\lambda_{0}\tilde{\lambda}_{0}\bigg{[}(1-\tau\beta)\underbrace{\langle\psi({\bm{\lambda}})|\tilde{\psi}\rangle}_{=C_{\mathds{1}}}+\tau\alpha\underbrace{\langle\psi({\bm{\lambda}})|\hat{S}_{+}\tilde{\psi}\rangle}_{=C_{S_{+}}}+\tau\gamma\underbrace{\langle\psi({\bm{\lambda}})|\hat{S}_{-}\tilde{\psi}\rangle}_{=C_{S_{-}}}\bigg{]}$
$\displaystyle+\tilde{\lambda}_{0}^{2}\bigg{[}1-2\tau\beta+2\tau(\alpha+\gamma)\underbrace{\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}_{=\tilde{C}_{S_{+}}}\bigg{]}$
$\displaystyle+\tilde{\lambda}_{0}^{2}\tau^{2}\bigg{[}\alpha^{2}+\beta^{2}+\gamma^{2}-2\beta(\alpha+\gamma)\underbrace{\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle}_{=\tilde{C}_{S_{+}}}+2\alpha\gamma\underbrace{\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle}_{=\tilde{C}_{S_{++}}}\bigg{]},$
(50)

where $C_{\mathds{1}}$ expresses the contribution of the identity part and
$C_{S_{+/-}}$ the contribution of the shift parts. The contributions
$\tilde{C}_{S_{+}}$ and $\tilde{C}_{S_{++}}$ to the cost function depend on
the solution of the previous time step only, and are hence constants. Note
that these different terms are evaluated separately and summed classically to
give the cost function. This means that one re-prepares the parametrized state
a few times every integration time step.

The cost function is evaluated by an adaption of a fundamental quantum
circuit, the so-called Hadamard test. In general, the Hadamard test provides
an expectation value $\Re\langle\psi|\hat{U}|\psi\rangle$ for any variable
$\hat{U}|\psi\rangle$ (see B). The measurement on the ancilla qubit $q_{0}$
delivers a measure for the manipulation on the lower qubits $q_{1}$ to
$q_{n-1}$. This measurement is performed such that $p_{0}-p_{1}$ is evaluated,
where $p_{0}$ and $p_{1}$ is the probability to measure the standard basis
eigenstates $|0\rangle$ and $|1\rangle$ at the ancilla qubit $q_{0}$,
respectively. In order to evaluate $C_{\mathds{1}}$ which is
$\langle\psi({\bm{\lambda}})|\tilde{\psi}\rangle$, the parameterized quantum
ansatz for the solution $\hat{U}({\bm{\lambda}})$ and the inverse of previous
time step $\tilde{U}^{\dagger}$ are implemented as controlled gates. If
$\hat{U}({\bm{\lambda}})$ initializes a state which is completely removed by
$\tilde{U}^{\dagger}$, the probability to measure the $|0\rangle$ state would
be $p_{0}=1$ because just the Hadamard gates by themselves, cancel their
effects causing no net rotation in total. For the evaluation of the $C_{S+/-}$
which is $\langle\psi({\bm{\lambda}})|\hat{S}_{+/-}\tilde{\psi}\rangle$, the
shift operation is added by implementing controlled NOT gates (CNOT) and
Toffoli gates which are organized in a particular way. For the $C_{S+}$ case,
the CNOT and Toffoli gates are organized in reverse order compared to
$C_{S-}$. The structure is shown in Fig. 4. The evaluation of
$\tilde{C}_{S_{+}}$ and $\tilde{C}_{S_{++}}$ requires the implementation of
$\hat{\tilde{U}}$ instead of $\hat{U}({\bm{\lambda}})$. In order to realize
the double shift operation for $\tilde{C}_{S_{++}}$, the $\hat{S}_{+}$ block
can be either implemented twice, or more efficiently, the processing structure
starts one qubit lower such that $q_{1}$ is not affected by the shift
operator.

$q_{0}$$q_{1}$$\vdots$$q_{n}$$H$$H$$\hat{U}({\bm{\lambda}})$$\hat{S}_{-}$$\ddots$$\hat{S}_{+}$$\iddots$$\hat{\tilde{U}}^{\dagger}$
Figure 4: Quantum circuits for the evaluation of the main cost contributions
$C_{\mathds{1}}$, $C_{S_{+}}$ and $C_{S_{-}}$. For the evaluation of
$C_{\mathds{1}}$, the both shift blocks are neglected, while the evaluation of
$C_{S_{-/+}}$ requires the $\hat{S}_{-/+}$ block. Qubit $q_{0}$ is an
ancillary (in short ancilla) qubit which collects the information for the
measurement to the right.

The quantum ansatz function can be designed problem-specific or generic
without any knowledge about the form of the solution. In this work, a quantum
ansatz with an universal function is used which is shown in Fig. 5. The ansatz
is defined by a special structure of $R_{y}$ rotation gates and CNOT gates.
The $R_{y}$ gates are parameterized with the parameter set and perform
rotations by $\lambda_{i}/2$ about the y axis of the qubit. CNOT gates negate
the state of the target qubit whenever the control qubit is in state
$|1\rangle$. This ansatz allows to implement all possible quantum states. The
trade-off is, that the solution can always be found, but the optimization has
as many parameters to tune as reasonably possible. The usage of other ansatz
functions showed no improvement in performance, as discussed in Sec. 5.4.
However, the inherent disadvantage of the considered universal ansatz function
is the circuit depth which would diminish a possible quantum advantage. This
ansatz requires $2^{n}-1$ parameterized gates and thus, $2^{n}$ parameters
(one additional parameter $\lambda_{0}$ for normalization purpose) need to be
optimized which leads to a high computational effort in circuit execution and
optimization.

$q_{0}$$q_{1}$$q_{2}$$\lambda_{1}$$R_{y}$$\lambda_{2}$$R_{y}$$\lambda_{3}$$R_{y}$$\lambda_{4}$$R_{y}$$\lambda_{5}$$R_{y}$$\lambda_{6}$$R_{y}$$\lambda_{7}$$R_{y}$
Figure 5: Example of a universal quantum ansatz function
$\hat{U}({\bm{\lambda}})$ for a qubit amount of $n=3$ with parameterized
$R_{y}(\lambda_{i})$ and CNOT gates.

The Nelder-Mead algorithm [39] is chosen as the classical optimization
algorithm. This algorithm is designed to solve unconstrained optimization
problems by a geometric method. For this, the function values of the cost
function are evaluated at some points. These points define the so-called
simplex. For a two-dimensional data space, a simplex would correspond with a
triangle. In the optimization process, the simplex is transformed by
reflections and expansions or contractions of the sides of the triangle.

We also tested other classical optimization algorithms and found that the the
Nelder-Mead algorithm is most suitable for the present problem in low-
dimensional data spaces. Thus, it is used mainly. The comparison of our
results with the other popular classical optimization algorithms is reported
in detail in C.

Figure 6: Comparison of the analytical solution to the results of the quantum
algorithms. The results are shown as contour plots in $x$ direction over the
time $t$ for (a) the analytical solution (ANA), (b) the QLSA results and (c)
the VQA results. The mean squared error (MSE) is evaluated to determine the
deviation to the analytical solution. For this, the MSE is evaluated with
respect (d) to the classical backward (red) and forward (blue) Euler method,
(e) to the QLSA results and (f) to the VQA results. Times are given in units
of $\tau_{a}=2L/U$ and length in units of $2L$. The $x$-axis was resolved with
4 qubits in both cases.

##  5 Comparison of quantum algorithms

###  5.1 Time evolution of concentration profile

In this section, the time evolution of the concentration profiles of the QLSA
and the VQA is shown and compared to the analytical solution, cf. eq. (10).
For this comparison, the 4-qubit case with $N=16$ is chosen. The parameters
are time step width $\tau=0.001$ s, diffusion constant $D=1$ m2/s, unit length
$2L=1$ m, and constant advection velocity $U=10$ m/s. The characteristic time
scales of the problem are the advection and the diffusion times. They are
given by $\tau_{a}=0.1$ s and $\tau_{d}=1$ s. From now on, we proceed with the
dimensionless form. Characteristic scales are combined in the dimensionless
Péclet number which is given by

${\rm Pe}=\frac{2UL}{D}=\frac{\tau_{d}}{\tau_{a}}=10\,.$ (51)

Figure 6 provides a first impression of the dynamical evolution of the
concentration profile in the form of a contour plot. We provide the analytical
solution together with those obtained from QLSA and VQA. The bottom row shows
the time evolution of the corresponding mean squared errors which will be
detailed below.

Figure 7: Comparison of the concentration profiles of the analytical solution
(ANA) and the results of the VQA and the QLSA for $N=16$ and a time step width
$\tau=0.01\tau_{a}$. (a)–(e) illustrate the time evolution of the
concentration profiles for increasing time. In panel (f), the corresponding
MSE is shown where the results of the quantum methods are compared to the
analytical solution. The VQA algorithm used 4 qubits for the spatial
discretization plus 1 ancilla qubit for the Hadamard tests which gives $n_{\rm
tot}=5$. The QLSA algorithm used 4 qubits for the spatial discretization plus
1 ancilla qubit for converting the matrix $A$ into a hermitian matrix. A
further ancilla qubit is used together with 8 additional qubits for the
register of the QPE. Thus, a total of $n_{\rm tot}=14$ qubits are used in
QLSA.

The corresponding concentration profiles are plotted in Fig. 7 where the time
interval is approximately $1/30\tau_{d}$ or $1/3\tau_{a}$. The concentration
profiles of the VQA reproduce the advection and diffusion process as expected
but the accuracy is limited by the Euler method used in the cost function
considered; see subsection 4.2. Especially for the early time steps, the
performance of the VQA is excellent (see Figs. 7a and 7b). In the course of
the time evolution, the advection-diffusion process starts to depart slightly
in comparison to the analytical solution: see Figs. 7(c)-(e).

This behaviour can also be seen in the time evolution of the mean squared
error (MSE) in Fig. 7(f) where the VQA result is compared to the analytical
solution. The MSE is given by, cf. eq. (30),

${\rm MSE}(t_{m})=\frac{1}{N}\sum_{i=1}^{N}\left[c_{i}^{m}-c_{i}^{\rm
ANA}(t_{m})\right]^{2}\,.$ (52)

The MSE curve decreases first, but starts to increase for $t\gtrsim 0.1$
$L/U$. The reason for this behaviour can be assigned to the phase when the
largest fraction of the concentration profile crosses the periodic boundary
for the first time. This is connected with stronger changes in the
parameterized vector components $(\lambda_{0},{\bm{\lambda}})$. In more
detail, comparing the iterative update of the parameter set, for a time step
where the maximum concentration is away from the periodic boundaries, the
components $(\lambda_{0},{\bm{\lambda}})$ instead change rapidly when the
boundary is crossed by the bulk of distribution. This is a specific property
of the geometric Nelder-Mead algorithm. We continued with this algorithm since
it still gave the lowest MSE amplitudes for the present problem, boundary
conditions, and qubit number range. See also Appendix C for more discussion.
With progressing time evolution, the problem fades out and the MSE curve
decreases again. This non-monotoneous behaviour was observed for system sizes
$N\geq 16$; here the number of optimization parameters was always similar to
the number of grid points $N$.

The concentration profiles computed with the QLSA using the implicit Euler
method (BTCS) also capture the physics of the advection-diffusion process very
well, both qualitatively and quantitatively. It should be reiterated here that
the error in the QLSA solutions (or from any algorithm for that matter) with
respect to the analytical solution is bounded from below by the error of the
classical solutions from the same underlying numerical scheme, in this case
the implicit Euler method. With this in mind, we can now observe from Figs.
7(a) and (b) that the QLSA, in contrast to the VQA, deviates from the
analytical solution only during the initial few time steps (for small $t$).
However, this is natural since the classical implicit Euler solution also
deviates almost exactly in the same manner the QLSA solution does. In fact,
the QLSA performs excellently when compared to the classical BTCS solution
alone, which we shall discuss more closely in the next subsection. This
behaviour is anticipated from Fig. 1(c) where, for the problem under
discussion, the MSE of the BTCS is in general higher than the FTCS scheme
which forms the basis to the VQA solutions. Proceeding further, the QLSA
performs progressively better for increasing $t$, as can be seen in Figs.
7(c)-(e), when it begins to closely follow the analytical solution. This is
quantified by observing the monotonic decay in the MSE of QLSA with evolving
time, as shown in Fig. 7(f).

In contrast to the VQA, the performance of QLSA improves with increasing
system size as one would naturally expect from higher degrees of resolution.
On the other hand, similar to the blockades posed by the parametric
optimization in the VQA, the accuracy of QLSA critically depends on (a) large
enough registers $n_{q}$ for the Quantum Phase Estimation and (b) the right
choice of $T_{0}$, see eq. (36) in Step 2. For instance, though the MSE of
QLSA asymptotes closely with the MSE of VQA for large $t$ in Fig. 7 (f), these
MSE values of the QLSA can, in general, be further lowered by providing a
higher $n_{q}$, without any increase in the finite difference resolution. We
investigate all such dependencies more closely in the following sections.

The maximum number of grid points is $N=64$ which corresponds to 6 qubits. In
this case, $64$ parameters have to be optimized, as described in C for a
detailed explanation of the classical optimization. The corresponding
concentration profiles are shown in Fig. 8. Within the short time interval
considered, the advection-diffusion dynamics can be reproduced very well by
the VQA.

Figure 8: Concentration profiles for $N=64$ of the analytical solution (ANA),
the VQA and QLSA results for (a) $t=6.1\times 10^{-3}\tau_{a}$ and (b)
$t=1.22\times 10^{-2}\tau_{a}$ where the time step width is $\tau=6.1\times
10^{-4}\tau_{a}$. The VQA algorithm used 6 qubits for the spatial
discretization plus 1 ancilla qubit which gives $n_{\rm tot}=7$. The QLSA
algorithm used 6 qubits for the spatial discretization plus 1 ancilla qubit
for the converting the matrix $A$ into a hermitian matrix. A further ancilla
qubit is used together with 6 additional qubits for the register of the QPE.
This results to $n_{\rm tot}=14$ qubits for QLSA.

###  5.2 Dependence on the number of qubits

With the number of qubits the resolution of the spatial discretization $N$
increases exponentially as $N=2^{n}$. Apart from the spatial resolution $N$,
the total number of qubits $n_{\rm tot}$ in both algorithms is as follows:

$\displaystyle n_{\rm tot}^{\rm VQA}$ $\displaystyle=\log_{2}(N)+1\,,$ (53)
$\displaystyle n_{\rm tot}^{\rm QLSA}$ $\displaystyle=\log_{2}(N)+2+n_{q}\,.$
(54)

In case of QLSA, an additional register with $n_{q}$ qubits is required which
corresponds to the QPE qubits.111In case of QLSA-2, the need for $n_{q}$ can
be eliminated given the absence of QPE. They determine the accuracy of the
eigenvalue estimation. The specific dependence on $n_{q}$ will be dealt with
in subsection 5.3. The discussion in this section focuses on the dependence on
the number of qubits associated with resolution alone.

For this investigation, the diffusion constant and the velocity are fixed to
$D=1$ and $U=10$ and the cases for $N=8,16,32$ grid points are analyzed. In
case of the VQA, the time constant $\tau$ is adapted for each discretization.
The cost function (see eq. (4.2)), which is derived in Sec. 4.2, includes the
prefactor $(1-2DN^{2}\tau)$. In order to include the term with this prefactor,
the time step width $\tau<1/2DN^{2}$ besides the Courant-Friedrichs-Lewy (CFL)
condition, $\tau<1/(NU)$. Thus, the time steps are $\tau=4\times 10^{-3}$ for
$N=8$, $\tau=10^{-3}$ for $N=16$, and $\tau=2.4\times 10^{-4}$ for $N=32$. If
we take a CFL number of 0.5, the time steps are smaller by about a factor of
1.5, 3, and 7 than would be classically possible for the first-order scheme.

For fair comparison, the same set of system parameters is prescribed for the
QLSA simulations as well. In the classical finite difference method, which is
the basis for the cost function of the VQA, a finer resolution results in a
decreased error. For the VQA, this means a finer resolution increases the
number of states, hence an increase of the qubit amount. For a time $t\leq
0.12\tau_{a}$, it can be shown that the error decreases for cases with a
higher number of qubits. The evaluation of the MSE over the time
$t=0.04-0.24\tau_{a}$ is shown in Fig. 9(a). It can be seen that a larger
number of qubits lead to smaller errors. However, the error for $N=8$
decreases while the curves for $N=16$ and $32$ show a rapid increase. The
reason for these inaccuracies in the concentration profiles, which lead to the
increased MSE, is the crossing of the periodic boundary of the bulk of the
concentration profile as discussed above. This is also shown in Fig. 9(b).
Here, it can be seen that the case for $N=8$ reaches lower cost than the cases
for $N=16$, $32$. Thus, it can be assumed that the global minimum in the
higher-dimensional parameter space of the optimization is harder to find and
hence, the concentration profiles differ slightly from the analytical
solution. This can be seen in Figs. 9(c)–(h). It can also be seen that
problems in reproducing the analytical solution appear especially after
crossing the boundary.

Figure 9: Comparison of VQA results for $N=8,16,32$ to the analytical solution
(ANA) with (a) the Mean Squared Error (MSE) over time and (b) the cost
function over time. Panels (c) and (d) show the concentration profiles for
$N=8$, panels (e) and (f) for $N=16$ and panels (g) and (h) for $N=32$ where
$t=0.12\tau_{a}$ and $t=0.24\tau_{a}$ are considered for each discretization
case.

In the case of QLSA, we compare its performance with respect to both the
analytical solution and the classical BTCS solution. It should be noted here
that, in order to study the effect of resolution alone, we assign a fixed,
sufficient number of $n_{q}=\log_{2}(N)+2$ qubits for each case for QPE module
of the algorithm. This corresponds to the minimum number of qubits required
for any given $N$, such that the solutions are of comparable performance and
$n_{q}$ affects every case somewhat similarly, although for increasing
resolutions one would need far larger $n_{q}$ registers and therefore lack of
which would still bear a weak effect. For $N=8$, $16$, and $32$, one would
therefore need to assign a total of $n_{\rm tot}=2\log_{2}(N)+3$ qubits, that
is, $n_{\rm tot}=10$, $12$ and $14$ qubits, respectively.

More generally, one needs to assign at least
$n_{q}=\max\\{\log_{2}(N)+1,\log_{2}(\kappa)\\}$, where $\kappa$ is the
condition number of $\tilde{A}$ (see eq. (34)). That said, the effect of
increasing resolution has a clear consequence of lowering the MSE with respect
to the analytical solution for all $t$, which is shown in Fig. 10(b). This can
also be qualitatively observed from Figs. 10(c)–(h). It can also be seen that
in all those figures that QLSA follows the BTCS very closely; however, when
quantified by computing the MSE with respect to classical BTCS solution as
shown in Fig. 10(a), we observe a rather non-trivial trend with increasing
resolution. Firstly, this MSE in this case is overall lower than the MSE in
Fig. 10 (b), suggesting that QLSA solutions are performing extremely well in
closely reproducing the classical BTCS, which forms the basis of QLSA
discretization. However, one can see that, overall, the $N=8$ case performs
the best followed by $N=32$ and $N=16$ which, more or less, have close time
evolution trends. This behavior is purely an artifact of the $n_{q}$ assigned
in each case. The $n_{q}$ provided for increasing $N$ is progressively
inadequate to foster an accurate eigenvalue estimation.

This can be seen more pronounced in Figs. 8(a)-(b) which shows the $N=64$ case
for small $t$. The QLSA solution though effectively reproduces the analytical
solution barring a modest quantitative error which is seen as spurious
oscillations around zero. This quantitative deviation can again be attributed
to two factors – (1) Inadequate $n_{q}$ (which also causes sign flips around
zero for small values of the solution field) which in turn causes improper
sign handling and evaluation of negative eigenvalues. The seemingly small and
negative concentration values, are not essentially non-physical, but are just
values with the wrong sign, which once measured can readily be flipped to
positive values classically. However, we still show this to highlight that the
sign handling quantum subroutines also suffer with insufficient $n_{q}$. (3)
The expected errors of solutions from BTCS based schemes in the initial few
time steps. The consequence of such insufficient resource allocation and its
remedy is further detailed in subsection 5.3. In essence, we can summarize
from the above that the performance and accuracy of QLSA when compared to the
analytical solution clearly increases with increasing resolution, when
provided with adequate algorithmic resources performs very well as already
seen in Figs. (6) and (10).

Figure 10: Comparison of QLSA results for $N=8,16,32$ and their MSE over time
with respect to (a) implicit time stepping (BTCS) classical solution and (b)
the analytical solution (ANA). Panels (c) and (d) show the concentration
profiles for $N=8$ at times $t=0.08\tau_{a}$ and $t=0.16\tau_{a}$, panels (e)
and (f) for $N=16$ at times $t=0.08\tau_{a}$ and $t=0.15\tau_{a}$ and panels
(g) and (h) for $N=32$ at times $t=0.0768\tau_{a}$ and $t=0.1488\tau_{a}$,
respectively.

###  5.3 Accurate eigenvalue estimation with QLSA

QLSA, specifically the QLSA-1 discussed in this work, relies heavily on
accurate estimation by QPE of the eigenvalues of the matrices under
discussion. The errors in this module occur from two primary sources:

(1) Numerical truncation: We recall from expression (37), that the process of
estimating eigenvalues in the QPE module requires an intermediate encoding of
those values into a binary format using $n_{q}$ qubits. For a given value, an
insufficient $n_{q}$ will naturally cause truncation errors of the order
$O(2^{-(n_{q}+1)})$. However, given an $n_{q}$, the eigenvalues can always be
scaled by choosing an appropriate $T_{0}$ such that $\sigma_{j}T_{0}$ can be
represented with the required accuracy. Unfortunately, since the eigenvalues
are unknown a priori, the choice of both $n_{q}$ and $T_{0}$ becomes elusive.
Figure 11(e) depicts the intricate connection between the two quantities and
their effect on the MSE. A similar contour can be made for the fidelity of the
solution as well. The fidelity $F$ would be given by

$F=\dfrac{\sum_{i=1}^{N}|c_{i}^{\rm QLSA}c_{i}^{\rm BTCS}|}{\|c^{\rm
QLSA}\|_{2}\,\|c^{\rm BTCS}\|_{2}}\,.$ (55)

$F$ is a measure of overlap instead of the difference. However, as noted in
[18], it would provide only a rough indication of the QLSA performance. So,
for brevity of discussion, we limit ourselves to the computing of MSE only.
Though some recommendations for the choice of $T_{0}$ can be made by bounding
the minimum and maximum eigenvalues, with functions of either the condition
number $\kappa$ or the trace of the matrix, they would still be rough
estimates.

The optimal choice of $T_{0}$ would be such that (a) all eigenvalues are
almost exactly represented, and (b) the MSE (with respect to analytical
solution) of the concentration field, should neither diverge nor oscillate
with time, and decrease with increasing $n_{q}$. This estimation of $T_{0}$ is
described in [18]. In summary, it requires one to first compute the behavior
of condition number $\kappa$ with increasing system sizes. If accurate
estimation of $\kappa$ turns out to be expensive, they can also be estimated
by tight, theoretical upper bounds (which, of course, would give less accurate
results). With this relation in hand, the system of equations is then solved
with QLSA for a smaller range of system sizes ($N,t$), $n_{q}$ and
$T_{0}\in[0,2\kappa(1-2^{-n_{q,max}})]$. From these results an MSE is computed
with respect to classical or analytical solution (available in this case) for
every combination of $n_{q}$ and $T_{0}$ as shown in Fig. 11(e), for the $N=8$
case integrated up to $t=0.1\tau_{a}$. The MSE could be of either the entire
concentration field (as is the case here), or of a function of the
concentration field, such as the scalar dissipation computed using the by
Quantum Post Processing (QPP) protocol [18], as denoted in Fig. 2(b).
Computing the latter is more efficient and speed-up preserving since it avoids
measuring the entire field – which is a $\mathcal{O}(N)$ operation, and also
minimizes the measurement errors associated with it.

Proceeding further, the trajectory of the minimum MSE is traced for every
$n_{q}$ and $T_{0}$ as shown in Fig. 11(e) (cyan dotted line) to find a
$T^{*}_{0}$ for which most eigenvalues are accurately represented with $n_{q}$
qubits. Finally, using the previously computed $\kappa-N$ relation, a new
relation between $N$ and $T^{*}_{0}$ is determined (power-law like behavior),
with which one can predict with nominal confidence a $T^{*}_{0}$ for all large
system sizes. Note that, for a given problem, this exercise needs to be
performed only once and larger system sizes can thereafter be solved with
minimal classical precomputing. With the right choice of $T_{0}$ we now solve
the system for $N=16$, $\tau=0.001$ with increasing $n_{q}$—and therefore $n$
is given by eq. (54). In this case, $n_{q}\in[4,8]$ and thus $n\in[10,14]$.
The MSE is computed with respect to analytical and BTCS solution as shown in
Figs. 11(a)–(c). Three observations are possible:

(i) The overall magnitudes of MSE between QLSA and BTCS is lower than with
QLSA and the analytical solution. This is expected since QLSA is based on the
BTCS scheme and thus follows the classical solution closely.

(ii) The MSE seems to exhibit a non-monotonic trend with time. The error is
initially high as expected when estimating a delta peak, consistent with Fig.
1(c). It decreases as the field tends to become more uniform due to diffusion.
On the other hand, the initial few time steps pose a beneficial setting to
QLSA, since many values of $c(t)$ are close to 0, and therefore errors in
eigenvalue estimation are somewhat diminished (except sign issues, which can
be corrected easily).

As the concentration field becomes more uniform and mixed, inaccuracies in
$T_{0}$ estimation manifest as a slight increase in MSE. It has to be
emphasized here that these errors and their trends also have contributions
from errors due to finite differences that are of the order $O(\tau,(\Delta
x)^{2})$. Another reason is the following, from eq. (39) the solution is of
the form $b_{j}/\sigma_{j}$. So the maximum error stems from the smallest
eigenvalue. If the $b_{j}$ associated with the smallest eigenvalue is
negligible, then the error from that is also relatively smaller. However, as
the concentration peak is advected in space and more $b_{j}$ become finite,
the error from the smallest eigenvalue becomes magnified. Further as the field
tends to diffuse, again the inaccuracy in eigenvalue estimation is smaller.
However, the final large $t$ value of MSE is bounded by $2^{-n_{q}}$. To
accurately estimate very small values of $c(t)$ would require a larger
$n_{q}$.

(iii) Finally, the effect of $n_{q}$ is clear from Figs. 11 (a) and (b).
Increasing $n_{q}$ tends to lower MSE in general, but it is in step-like
fashion as shown in Fig. 11(c) plotted for final time $t=0.4\tau_{a}$. This is
because increasing $n_{q}$ in small steps (of $\mathcal{O}(1)$) does not lower
the least count significantly (in $\log_{10}$ or $\log_{e}$). We also plot the
residue, given by $RES=|c(t)-c(t-\tau)|$, as a function of $t$ as shown in
Fig. 11 (d). The monotonic decay in residue symbolizes two aspects: (a) The
numerical method and choice of $\tau$ and $\Delta x$ produces stable non-
diverging solutions. (b) When the residue falls below a threshold (which can
be set arbitrarily small), a steady-state of the solution is reached. The
overall residue also decreases with increasing refinement as expected.

(2) Sign flips: Finally, a finer observation is that of the somewhat erratic
behavior in MSE in the initial few time steps as shown Fig. 11 (b) for the
$n_{\rm tot}=14$ case. This is because of improper handling of very small
negative eigenvalues. The eigenvalues are generally mirrored about 0, to lie
in the range [-0.5,0.5]. If the eigenvalues get too close to 0 or when
improperly scaled with $T_{0}$, the signs might flip causing a rough MSE
profile. This also manifests physically in the concentration field, as
depicted by the tiny dark peaks for $t$ close to 0 in Fig. 6(b). This can be
minimized by adding another sign control qubit or by increasing $n_{q}$ and
estimating $T_{0}$ more accurately.

(a)

(b)

(c)

(d)

(e)

Figure 11: (a) and (b) depict the evolution in time of the MSE from QLSA
computed with respect to the analytical and BTCS solutions respectively,
plotted for increasing number of qubits. (c) compares the MSE at $t=0.4$ (in
$\tau_{a}$) with respect to analytical and BTCS solutions as function of $n$.
(d) shows the time-decay of residue for $N=8,16$ and $32$. (e) Shows the
contour of MSE of QLSA with respect to analytical solution as function of both
$n_{q}$ and $T_{0}$. The dotted line (cyan) plots the trajectory of minimum
MSE for every combination of $(T_{0},n_{q}$). The color bar shows MSE in
logarithmic scale (of base 10).

###  5.4 Comparison of different ansatz functions for VQA

The quantum ansatz $\hat{U}({\bm{\lambda}})$ has to meet several requirements.
The ansatz should be able to construct the unknown solution for the next time
step in the given problem. Secondly, the amount of rotation and entanglement
gates should be reduced to a minimum in order to get an efficient and shallow
parametric circuit $\hat{U}({\bm{\lambda}})$. The efficiency of the ansatz is
a key factor in determining whether a quantum advantage can be achieved at all
or not. The currently applied universal ansatz (see Fig. 5) contains $2^{n}-1$
parameterized gates to generate the next $n$-qubit state. However with ${\cal
O}(N)$ parameterized gates, one cannot obtain a quantum advantage. In order to
improve the efficiency of the algorithm, further ansatz functions are now
tested. To this end, we analyse the performance of shallow tensor networks
(TNs) where $R_{y}$ and CNOT gates are structured in a staggered way, see e.g.
ref. [40]. These TNs are visualized in Fig. 12. The TN1 ansatz shows a generic
structure where the marked code block can be repeated as often as desired in
order to build $\hat{U}({\bm{\lambda}})$ with a different number of gates. In
TN2, a row of $R_{y}$ gates is added which is inspired by the universal ansatz
and should enable an easier generation.

The quality of an ansatz is quantified by a Hadamard test-like quantum circuit
which is shown in Fig. 4. In the $\hat{U}({\bm{\lambda}})$ code block, the
parameterized ansatz is initialized. Then, the inverse of the wanted solution
which is given as the classical FTCS solution of a certain time step $t$ is
initialized. The ancilla qubit is measured to obtain the probability of state
$|1\rangle$. In this way we determine the match (or overlap) of the generated
with the desired state. In detail, if the probability for the state
$|1\rangle$ is zero, the ansatz $\hat{U}({\bm{\lambda}})$ generates a state
which is perfectly uncomputed by the inverse of the wanted FTCS solution. In
other words, the ansatz allows to reconstruct the considered quantum state
exactly. If the probability for the state $|1\rangle$ is greater than zero, we
obtain a measure for the deviation between both states. This measure is called
identity costs $C_{\text{id}}$.

In this work, the considered ansatz structures are evaluated for quantum
registers of size $n=4$ and 6. For this, three concentration profiles are
chosen which capture the significant shapes of the advection-diffusion
problem. The corresponding identity costs $C_{\text{id}}$ are evaluated for
these chosen concentration profiles and for a varying number of parameterized
$R_{y}$ gates. In Fig. 13, the results for the $n=4$ qubits are shown. We can
observe that the universal ansatz leads to low identity costs of $\approx
10^{-11}$–$10^{-5}$, such that this ansatz is suitable to construct the wanted
concentration profile, but the number of used parameterized gates is ${\cal
O}(2^{n})$. The considered tensor networks (TN1, TN2) lead to increased
identity costs of $\approx 10^{-4}$–$10^{-2}$. The identity costs decrease
slowly for a higher number of parameterized gates, but the costs cannot reach
the level of the universal ansatz. Thus, one observes that the investigated
TNs are less suitable as an ansatz function for the considered advection-
diffusion problem.

Interestingly, even if the number of gates is similar to the universal ansatz,
the tensor networks cannot reproduce the given concentration profiles well.
The evaluation of both investigated TNs results in similar identity costs,
whereby TN1 tends to be more suitable for sharp Gaussian shaped concentration
profiles and TN2 seems to be more appropriate for concentration profiles which
are further decayed. Similar results are obtained for the $n=6$ qubit case
(see Fig. 14). The identity costs for TN1 and TN2 differ marginally. The
evaluation of the universal ansatz results in significantly lower costs for
No. of $\lambda_{i}\to N$. Furthermore, a reduction of the parameter space is
not possible, particularly for the reconstruction of the concentration
profiles in early times a high amount of parameterized gates is necessary (see
Fig. 14(d)).

To conclude, the investigated TN structures for $n=4$ and 6 qubits could not
achieve the required accuracy for the state vector generation. Thus, we
proceed with the universal ansatz for this investigation. It is also worth
noting here that, in contrast to the ansatz used here, the general quantum
state preparation (QSP) on the other hand (as used in QLSA), prepares a state
exactly, however preparing arbitrary states requires $O(N)$ depth as well.
Nevertheless efficient state preparation algorithms exist that prepare states
which either have functional forms (such as Gaussian-like or wave-like forms
as encountered in the current problem under discussion) or when they are
sparse states [18].

(a)
TN1$q_{0}$$q_{1}$$q_{2}$$q_{3}$$\lambda_{1}$$R_{y}$$\lambda_{2}$$\lambda_{3}$$\circlearrowright$(b)
TN2$q_{0}$$q_{1}$$q_{2}$$q_{3}$$\lambda_{1}$$R_{y}$$\lambda_{2}$$R_{y}$$\lambda_{3}$$R_{y}$$\lambda_{4}$$R_{y}$$\lambda_{5}$$R_{y}$$\lambda_{6}$$\lambda_{7}$$\circlearrowright$
Figure 12: Example of tensor networks (TN) as generic quantum ansatz functions
$\hat{U}({\bm{\lambda}})$ for $n=4$ qubits with parameterized
$R_{y}(\lambda_{i})$ and CNOT gates. The parameters are continuously indexed
for repeated blocks ($\circlearrowright$). (a) Full generic ansatz TN1 and (b)
Ansatz TN2 which includes an additional row of $R_{y}$ gates. Figure 13:
Comparison of the performance of the different investigated ansatz structures
for $n=4$ qubits. The upper panel shows the concentration profiles of the FTCS
solution for the time steps (a) $t=0.01\tau_{a}$, (b) $t=0.1\tau_{a}$ and (c)
$t=0.4\tau_{a}$ with $\tau_{a}=2L/U$. In the lower panel, the corresponding
identity costs $C_{\text{id}}$ of the investigated ansatz structures are shown
for a varying number of parameters $\\#\lambda$.  Figure 14: Comparison of the
performance of the different investigated ansatz structures for $n=6$ qubits.
The upper panel shows again the concentration profiles of the FTCS solution
for the time steps (a) $t=6.1\times 10^{-4}\tau_{a}$, (b) $t=5.54\times
10^{-2}\tau_{a}$ and (c) $t=1.7\times 10^{-1}\tau_{a}$. In the lower panel,
the corresponding identity costs $C_{\text{id}}$ of the investigated ansatz
structures are shown for a varying number of parameters $\\#\lambda$. We do
not show the individual grid points in (a)–(c) which sum to $N=64$.

###  5.5 Qiskit implementation of VQA with circuit noise

In this section, the VQA implementation is performed with a noisy quantum back
end, which is close to a real noisy intermediate scale quantum device. The
difference from the ideal state vector (SV) simulation is that it includes a
measurement process with sampling noise and error rates of the gates, such as
bit flips or phase flips. For this, the QASM simulator environment in Qiskit
is used which is a noisy quantum circuit back end. A single measurement of an
$n$-qubit quantum state on a quantum computer is a random projection on one of
the $2^{n}$ eigenstates with respect to an observable. This observable is the
$Z$ matrix on each qubit. In order to obtain the full quantum state vector,
such measurements have to be repeated many times to sample all eigenstates
sufficiently well. These repetitions of identically prepared quantum
simulations of each integration time step of the advection-diffusion equation
are termed shots. In this investigation, the number of shots is fixed to
$N_{S}=2^{20}$. This sampling error of the shots decreases with
$1/\sqrt{N_{S}}$.

Real quantum computers are never perfectly isolated from the environment; thus
many different types of decoherence errors appear at each of the individual
gates. They are smaller for single qubit gates than for entanglement gates. In
the simulation software Qiskit, the decoherence noise model implemented is
such that customized quantum errors can be set. Thereby, the probabilities for
the appearance of quantum gate errors ($p_{\text{gate}}$), errors in
measurement ($p_{\text{meas}}$) and in resetting ($p_{\text{reset}}$) of
qubits are defined. We have done a study for the case with $N=8$ and compared
the results to the corresponding ideal SV simulations reported earlier. To
this end, the noise model is implemented as follows. We choose the
probabilities $p_{\text{gate}}=0.008$, $p_{\text{meas}}=0.03$, and
$p_{\text{reset}}=0.0003$ that a gate error, a measurement error, and a qubit
reset error appear in the course of the quantum simulation.

Furthermore, the evaluation of the cost function is simplified to reduce the
appearance of decoherence noise in the quantum circuits. As discussed in
subsection 4.2, the costs are calculated on the basis of eq. (46). When the
last term is dropped, which is always a constant term, the number of quantum
circuits for the evaluation of the overlap terms can be reduced from 5 to 3.
Thus, the minimum of $C(\lambda_{0},{\bm{\lambda}})$ is technically no longer
at zero, but at a negative constant value.

The direct comparison with the ideal state vector simulation is shown in Figs.
15(a)–(b). It can be seen that the concentration profile with the QASM
simulator can reproduce advection and diffusion, but the profile differs
slightly from the those of the ideal simulation and the analytical solution.
The MSE is evaluated again as a measure of the deviation from the analytical
solution and shown in Fig. 15(c). As expected, the MSE for the QASM simulation
case is higher than that of the ideal simulation. Furthermore, it increases
with respect to time. This can be explained by the error propagation from the
previous step, which is included in this iterative framework. The cost
function of the QASM simulation case is increased in comparison with the state
vector simulation case, see Fig. 15(d), because the additional noise
complicates the optimization procedure.

Figure 15: Comparison between the shot based simulator (QASM) with statistical
noise where shot number is $2^{20}$ and the ideal state vector (SV) simulator
results and the analytical solution (ANA) with the concentration profiles for
(a) $t=0.04\tau_{a}$ and (b) $t=0.2\tau_{a}$ for the the advection-diffusion
case. In (c), the mean squared errors (MSE) of the different simulation
methods are compared from which the deviation from the analytical solution is
evaluated; in (d) the corresponding cost functions are shown.

##  6 Final discussion and outlook

The goal of the present work has been to present a one-to-one comparison of
two quantum algorithms to simulate a simple linear flow problem numerically.
In the present work we considered a time-dependent linear and one-dimensional
advection-diffusion problem on the unit interval with periodic boundary
conditions. The time evolution of this fluid mechanical problem is not unitary
and hence requires specific steps to be taken in both algorithms. A passive
scalar or concentration profile $c(x,t)$ is advected by a constant velocity
$U>0$ and subject to a constant molecular diffusion $D$. The Péclet number is
$Pe=10$.

The two algorithms chosen were a Quantum Linear Systems Algorithm (QLSA) and a
Variational Quantum Algorithm (VQA), both of which are hybrid quantum-
classical in nature. We have investigated their performance on computational
grids varying between $N=8$ and 64, which correspond to 3 and 6 qubits,
respectively. We were able to show that both algorithms perform well for the
numerical solution of the fluid mechanical problem with a first order time
integration scheme, using either a backward or a forward Euler integration
scheme. The accuracy of the time evolution in both quantum algorithms, i.e.,
the forward Euler (FTCS) for VQA and backward Euler (BTCS) for QLSA is bounded
from below by the round-off errors of the corresponding classical integration
schemes. Accuracy was quantified by a mean squared error (MSE).

We have shown that both algorithms involved detailed pre-conditioning with
respect to specific aspects; this was the major part of this work and could,
in our view, be of interest to other users of these specific quantum
algorithms. In QLSA, the central point that required comprehensive
investigations is related to the approximate determination of eigenvalues of
the unitary matrix $\hat{U}(t)=\exp(i\tilde{A}t)$ in the Quantum Phase
Estimation (QPE) stage. It is demonstrated that the number of additional
qubits $n_{q}$ needed for this task and appropriate pre-conditioning is key to
the accuracy of the QLSA method. In case of the VQA, the classical
optimization algorithm to determine the minimum of the cost function
$C(\lambda_{0},{\bm{\lambda}})$ turns out to be the bottleneck. In the present
work, we found that the geometric Nelder-Mead algorithm gave the best results,
despite a non-monotonic time evolution of the MSE; see also Appendix C. This
result holds for the present benchmark task and the chosen boundary
conditions. For example, Dirichlet boundary conditions of the concentration
profile at the wall might eliminate this behaviour.

Our results suggest immediate directions for future research for both
algorithms. For QLSA, the ongoing and upcoming work focuses on developing
algorithms, which are mainly based on the concept of Linear Combination of
Unitaries (LCU) [35, 18] (QLSA-2) and eliminate the need for QPE. This
ameliorates the higher circuit depths and gate count encountered in QLSA-1,
making it more suitable for implementation on NISQ devices. Extending these
tools to solve nonlinear flow problems by new embedding techniques such as
Homotopy Analysis forms a major part of the future work [19]. In the case of
VQA, surrogate algorithms for the global minimum search of the cost function
have been suggested recently [41]. Finding minima in high-dimensional
parameter spaces is a general problem of quantum algorithms. This includes
quantum machine learning where barren plateaus limit the efficiency of
implementation [38]. The strength of the VQA might become better visible for
nonlinear problems already attempted, e.g. in Lubasch et al. [21]. These
problems will, however, require higher-order time integration scheme to avoid
numerical instabilities, e.g. when time-dependent nonlinear Schrödinger
equations have to be solved. These equations describe also nonlinear wave
phenomena in fluid mechanics [42].

In the end, we wish to provide a few more general comments on the subject as a
whole. The numerical implementations of the classical fluid flow problems as a
quantum algorithm have so far not gone beyond the proof-of-concept level. We
discuss mostly one-dimensional linear and nonlinear problems while the
realistic flows are two- and three-dimensional. Many studies, including most
of the existing ones, are implemented in ideal quantum simulation frameworks,
thus avoiding the decoherence problems of real and noisy quantum devices that
are state-of-the art today. These considerations appear to hinder
demonstration of true quantum advantage. In case of the variational methods,
most algorithms do not come with theoretical guarantees of quantum advantage
(or complexity). The advantage is contingent on problem-specific
implementation of the ansatz as well as the parametrization and the
optimization methods. On the other hand, QLSA algorithms come with theoretical
guarantees of quantum complexity and advantage. However, these algorithms tend
to be very sensitive to parameters such as sparsity $S$ of the linear systems
matrix $\tilde{A}$, its condition number $\kappa$, or the choice of unitary
bases in case of methods of LCU, making it hard to project their performance
on real quantum devices. In both approaches, one also needs to account for the
number of shots needed to sample the final quantum state. Therefore careful
implementation of Quantum Amplitude Amplification [43] is necessary such that
one obtains the solution while maintaining quantum advantage.

A desired quantum advantage will most possibly require us to rethink the
solution of classical flow problems even more as a quantum mechanical problem.
This might be obtained by transforming a nonlinear problem, which is
numerically formulated in a finite-dimensional space (for example by a
Galerkin method), to a corresponding linear problem in a much higher-
dimensional (theoretically infinite-dimensional) Hilbert space. In the latter,
the encoding capacity of quantum algorithms would fully unfold. One possible
pathway in this respect can be provided by the quantum mechanical
implementation of Carleman embeddings [37, 44], the Koopman operator formalism
[45, 46, 47] or the Homotopy analysis method. Apart from these, the quantum
volume222Representing the combined measure of the size of quantum circuits
(qubits and depth) that can be reliably used [48]. of the current and near-
term quantum devices is an important consideration while designing algorithms
in the hope of harnessing any quantum advantage. Future investigations will
probably show us if these routes are indeed successful, and leave us with new
scenarios in this research field.

Acknowledgements
We wish to thank Georgy Zinchenko for insightful discussions. The work of J.I.
is funded by the European Union (ERC, MesoComp, 101052786). Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union or the European Research Council. Neither
the European Union nor the granting authority can be held responsible for
them. P.P. is supported by the project no. P2018-02-001 ”DeepTurb – Deep
Learning in and of Turbulence” of the Carl Zeiss Foundation, Germany.

##  Appendix A Shift operations for cost function
$C(\lambda_{0},{\bm{\lambda}})$

In the following, the shift operations are explained by a two-qubit example.
These operations have been used in the evaluation of the cost function in the
VQA. For this, a quantum register with the qubits $q_{1}$ and $q_{2}$ is
considered to be in the initial state $|01\rangle$. For a shift to the left
which is defined as $\hat{S}_{-}$ operation, an $X$ gate is implemented on the
first qubit and afterwards, a controlled NOT gate (CNOT) acts on the register
as it is shown in Fig. 16(a). Consequently, the register is in the following
state:

$\displaystyle|10\rangle\xlongrightarrow[]{B^{\prime}}|11\rangle\xlongrightarrow[]{B}|01\rangle$
(56)

For the $\hat{S}_{+}$ operation, the gates are organized reversely as it is
shown in Fig. 16(b). Then, the following states can be found:

$\displaystyle|10\rangle\xlongrightarrow[]{C^{\prime}}|10\rangle\xlongrightarrow[]{C}|11\rangle$
(57)

In the case of the considered cost function (see 4.2), the shift operations
are applied to fixed quantum states $|\tilde{\psi}\rangle$ within a Hadamard
test which is analogous to the evaluation of a scalar product in classical
computation. Now we want to show by an example that an application of a
$\hat{S}_{+}$ operations
($\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle$) equals the application
of any single shift operation for this special case. For this, the two-qubit
quantum state is $|\tilde{\psi}\rangle=(a,b,c,d)^{T}$. Then follows,

$\displaystyle\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle$
$\displaystyle=(a,b,c,d)\cdot(d,a,b,c)^{T}=ad+ba+cb+dc$ (58)
$\displaystyle\langle\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle$
$\displaystyle=(a,b,c,d)\cdot(b,c,d,a)^{T}=ab+bc+cd+da$ (59)
$\displaystyle\langle\hat{S}_{+}\tilde{\psi}|\tilde{\psi}\rangle$
$\displaystyle=(d,a,b,c)\cdot(a,b,c,d)^{T}=da+ab+bc+cd$ (60)
$\displaystyle\langle\hat{S}_{-}\tilde{\psi}|\tilde{\psi}\rangle$
$\displaystyle=(b,c,d,a)\cdot(a,b,c,d)^{T}=ba+cb+dc+ad$ (61)

It can be observed that
$\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle=\langle\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle=\langle\hat{S}_{+}\tilde{\psi}|\tilde{\psi}\rangle=\langle\hat{S}_{-}\tilde{\psi}|\tilde{\psi}\rangle$.
Furthermore, these shift operations can be applied to both factors of the
scalar product. If the same shift operation is applied to both scalar product
entries, e.g.,
$\langle\hat{S}_{+}\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle$, the identity
will be computed because

$\displaystyle\langle\hat{S}_{(+/-)}\tilde{\psi}|\hat{S}_{(+/-)}\tilde{\psi}\rangle=\langle\tilde{\psi}|\tilde{\psi}\rangle=\sum_{i}\tilde{\psi}_{i}^{2}=1.$
(62)

Furthermore, $\langle\hat{S}_{-}\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle$
can be rewritten to
$\langle\tilde{\psi}|\hat{S}_{+}\hat{S}_{+}\tilde{\psi}\rangle=\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle$,

$\displaystyle\langle\hat{S}_{-}\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle$
$\displaystyle=(b,c,d,a)\cdot(d,a,b,c)^{T}=bd+ca+db+ac$ (63)
$\displaystyle\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle$
$\displaystyle=(a,b,c,d)\cdot(c,d,a,b)^{T}=ac+bd+ca+db$ (64)

Analogously, it holds
$\langle\tilde{\psi}|\hat{S}_{--}\tilde{\psi}\rangle=\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle$.

$q_{1}=|0\rangle$$q_{2}=|1\rangle$$A$$X$$B^{\prime}$$B$ (a) $\hat{S}_{-}$
operation.

$q_{1}=|0\rangle$$q_{2}=|1\rangle$$A$$X$$C^{\prime}$$C$ (b) $\hat{S}_{+}$
operation.

Figure 16: Definition of the shift operations for two qubits with $X$ and CNOT
gates.

##  Appendix B The Hadamard test

In general, the Hadamard test is a method which allows to find the expectation
value $\Re\langle\varphi|\hat{U}|\varphi\rangle$. For this, a unitary gate
$\hat{U}$ acts on the qubit $q_{1}$ which is in the state $|\varphi\rangle$.
The corresponding quantum circuit is shown in Fig. 17. First, the quantum
circuit is in the state

$\displaystyle|q_{1}q_{0}\rangle_{A}=|\varphi\rangle\otimes|0\rangle.$ (65)

The first Hadamard gate acts on the zeroth qubit such that

$\displaystyle|q_{1}q_{0}\rangle_{B}$
$\displaystyle=|\varphi\rangle\otimes\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)=\frac{1}{\sqrt{2}}\left(|\varphi\rangle\otimes|0\rangle+|\varphi\rangle\otimes|1\rangle\right).$
(66)

The unitary gate is implemented on the qubit $q_{1}$ and is controlled to the
ancilla qubit $q_{0}$ such that follows:

$\displaystyle|q_{1}q_{0}\rangle_{C}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|\varphi\rangle\otimes|0\rangle+\hat{U}|\varphi\rangle\otimes|1\rangle\right).$
(67)

Considering the second Hadamard gate, the state of the quantum circuit changes
to the following one:

$\displaystyle|q_{1}q_{0}\rangle_{D}$
$\displaystyle=\frac{1}{2}\left(|\varphi\rangle\otimes(|0\rangle+|1\rangle)+\hat{U}|\varphi\rangle\otimes(|0\rangle-|1\rangle\right)$
(68)
$\displaystyle=\frac{1}{2}\left((\mathds{1}+\hat{U})|\varphi\rangle\otimes|0\rangle+(\mathds{1}-\hat{U})|\varphi\rangle\otimes|1\rangle\right)$

The measurement is performed in the standard $Z$ basis such that

$\displaystyle p_{0}-p_{1}$
$\displaystyle=\frac{1}{4}\langle\varphi|(\mathds{1}+\hat{U})^{\dagger}(\mathds{1}+\hat{U})|\varphi\rangle-\frac{1}{4}\langle\varphi|(\mathds{1}-\hat{U})^{\dagger}(\mathds{1}-\hat{U})|\varphi\rangle$
(69)
$\displaystyle=\frac{1}{2}\langle\varphi|\hat{U}^{\dagger}+\hat{U}|\varphi\rangle=\Re\langle\varphi|\hat{U}|\varphi\rangle\,.$

$q_{0}=|0\rangle$$q_{1}=|\varphi\rangle$$A$$H$$\hat{U}$$H$$B$$C$$D$ Figure 17:
Two qubit quantum circuit which defines the general Hadamard test. Two
Hadamard gates are applied on the ancilla qubit $q_{0}$ and a controlled
unitary transformation $\hat{U}$ is applied between these Hadamard gates. The
measurement of the ancilla qubit returns the expectation value
$\Re\langle\varphi|\hat{U}|\varphi\rangle$ for the variable
$\hat{U}|\varphi\rangle$. This Hadamard test principle is used in the VQA to
evaluate the cost function.

##  Appendix C Classical optimization methods for the VQA

In this appendix, the classical optimization methods for the VQA are
introduced and compared. The optimization in the considered advection-
diffusion problem is challenging due to different aspects. First, the number
of parameters which are optimized and hence, the complexity of the
optimization problem, scales with the qubit number. Consequently, a fine
discretization with a large number of qubits and the chosen ansatz function
lead to a high-dimensional parameter space and a complex-shaped cost function
for the classical optimization. Secondly, vanishing gradients which drive the
search in a local minimum, also known as barren plateaus [49], complicate the
search for the global minimum.

Furthermore, the imposed periodic boundary conditions can induce rapidly
changing parameter sets at the boundaries. This aspect is visualized by the
simple example of a triangle function which is moving by one cell per time
step in Fig. 18. The movement far away from the boundaries results in small
changes of the parameter set $(\lambda_{0},{\bm{\lambda}})$. If the periodic
boundary is crossed and entries at the other boundary appear, the state vector
which models the concentration profile changes strongly and hence, the
corresponding parameter set ${\bm{\lambda}}$ shows major modifications. This
aspect is specific to geometric optimization algorithms. Lastly, the existence
of noise in the evaluation of quantum circuits contributes to the challenges
of the classical optimization. In this work, ideal simulations were considered
and hence, the impact of noise is neglected in the selection of the
optimization algorithm. For the comparison of the classical optimization
algorithms, the VQA is applied to the one-dimensional advection-diffusion
equation for the case with $N=16$. The parameters are $D=1$, $u=10$, and
$\tau=0.001$ and the computation is performed for a total time $T=30\tau$.

Figure 18: Comparison of the changing parameter set ${\bm{\lambda}}$ for an
advection crossing the periodic boundary by the example of a simple advection
of a hat signal for $N=8$. For advection far from the boundaries, the
concentration profile for a time (a) $t_{1}$ and (b) $t_{2}$ is shown with (c)
the corresponding distribution of the parameter set
${\bm{\lambda}}=(\lambda_{1},\dots,\lambda_{7})$. The concentration profiles
for a time (d) $t_{3}$ and (e) $t_{4}$ show the case where the periodic
boundaries are crossed and the corresponding distribution of the parameters
${\bm{\lambda}}$ can be found in panel (f). The comparison between the panels
(c) and (f) show that crossing the boundary layer results in larger changes of
the parameter vector which can be problematic for the optimization routine.

The Nelder-Mead algorithm (NM) [39] or downhill-simplex algorithm is designed
to solve classical unconstrained optimization problems without any gradient
approximation. The algorithm only uses the function values at some points
which construct the simplex in the hyperplane. This simplex is transformed by
geometric operations such as reflection, expansion, contraction and shrinking.
The Nelder-Mead algorithm uses a geometric method in order to find the minimum
of the given cost function. The application of the Nelder-Mead algorithm in
the considered advection-diffusion problem showed that this algorithm is
robust for small search spaces ($N\leq 16$). Moreover, the mean computation
time per time step is small in comparison to other methods. However, problems
appear if the qubit amount is increased or there are entries at the elements
at the boundary. Reasonable for this might be the fact that the parameter set
changes rapidly at the periodic boundary in comparison to the changes far away
from the boundary, see again Fig. 18. Consequently, other optimization
algorithms need to be considered for an increased qubit amount.

The Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS) applies a quasi-Newton
method for solving unconstrained, nonlinear optimization problems. Thereby,
the Hessian matrix of the cost function is approximated by the evaluation of
the gradients (or the approximated gradients) in order to find the descent
direction in the hyperparameter landscape. In this work, an adaption of the
BFGS algorithm is used. The Limited memory-BFGS algorithm for bound
constraints (L-BFGS-B) [50] uses a limited amount of computer memory which
makes the algorithm suitable for large search spaces. Furthermore, it can
handle bound constraints. In this investigation, the L-BFGS-B algorithm cannot
find the global minimum such that the mean MSE is approximately $10^{-4}$. A
possible reason for this is the disadvantageous initial parameter set. In
order to improve the performance, the L-BFGS-B method is combined with other
preceding optimization methods which aim at finding an appropriate region for
further optimization.

The first one is the combination of the Bayesian optimization and L-BFGS-B
algorithm (BO+L-BFGS). Bayesian optimization [51] is suitable for optimization
problems where the costs are given as black box functions, are expensive to
evaluate or include noise. This method approximates the cost function by a
Gaussian process regression based on previous observations. An acquisition
function determines the next samples for the observations whereby random
exploration steps can be added in order to include a wide range of
observations for the fitting of the cost landscape. This aims at finding a
promising region for further optimization with the L-BFGS-B algorithm. With
this preceding Bayesian optimization, the results of L-BFGS-B algorithm could
be improved such that the mean MSE is approximate $10^{-5}$, but the
computation time is increased. However, the application of this combination of
methods is reasonable if the system size is increased. For this, the test case
was expanded to $N=64$. Thereby, the Nelder-Mead optimization cannot find the
global minimum of the cost function ($C\approx 10^{-2}$) and hence, the
optimization fails. In contrast, the combination of Bayesian optimization and
L-BFGS-B algorithm shows small costs ($C\approx 10^{-10},\dots,10^{-5}$) and
good results in accuracy. This is shown qualitatively with the concentration
profiles in Fig. 19a and 19b and with the comparison of the cost functions
(Fig. 19c).

Secondly, the Adaptive moments algorithm (Adam) [52] is combined with the
L-BFGS-B algorithm. The Adam algorithm uses a gradient-based method to
determine the descent direction in the hyperparameter landscape. It includes
an adaptive learning rate and momentum for each update step of the parameter
which improves the performance in cases of sparse gradients and non-stationary
problems. Furthermore, it is suitable for the optimization of large parameter
sets. In this investigation, the combination of Adam and L-BFGS-B algorithm
can process the test case with an accuracy $\approx 10^{-5}$, but the
computational effort is too high to use this method for increased system
sizes. Reasonable for this high computation time is the large amount of
required iteration steps which all include the calculation of the gradients.

The Simultaneous Perturbation Stochastic Approximation (SPSA) [53] is an
optimization algorithm which uses a stochastic method to approximate the
gradient of the cost function. Thereby, the cost function is evaluated twice
with completely perturbed parameter sets. The parameters are chosen randomly
using a zero-mean distribution. This algorithm is robust to noise. In this
work, the SPSA optimization could not find the minimum such that the costs
were found to be $\approx 10^{-2}$ which results in high mean squared errors
($\approx 10^{-3}$).

In conclusion, the Nelder-Mead algorithm can be recommended in cases of low
parameter spaces ($N\leq 16$) due to its accuracy and the computation time. If
the system is increased it is advisable to chose the combination of Bayesian
optimization and L-BFGS-B algorithm. The comparison of the mean squared errors
and the cost functions of the VQA with the presented optimization algorithms
are shown in Fig. 20. Furthermore, an overview of the used optimization
algorithms, the corresponding methods, accuracy and computation efforts is
presented in table 1.

Figure 19: Comparison of the concentration profiles of the analytical solution (ANA) and the VQA results with the Nelder-Mead optimization (NM) and the combined Bayesian and Broyden-Fletcher-Goldfarb-Shanno algorithm optimization (BO+BFGS-L) for $N=64$. In (a) the concentration profile of the for $t=2\tau$ is shown and in (b) for $t=4\tau$. Here, the parameters are chosen to be $D=1$, $U=10$, $\tau=6.1\times 10^{-5}$. In (c) the corresponding cost function is shown. Figure 20: Comparison of investigated optimization methods with (a) the mean squared error (MSE) and (b) the cost function over time. Overall the combined optimization with Adam and L-BFGS algorithm shows the best performance for the considered test case with $N=16$, $D=1$, $U=10$ and $\tau=0.001$. Optimizer | Method | Accuracy | Computation time
---|---|---|---
NM | geometric | $\sim 10^{-5}$ | $0.75h$
L-BFGS-B | gradient-based | $\sim 10^{-4}$ | $0.75h$
BO | approximation | $\sim 10^{-5}$ | $2.25h$
& L-BFGS-B | gradient-based |  |
Adam | gradient-based | $\sim 10^{-5}$ | $7.5h$
& L-BFGS-B | gradient-based |  |
SPSA | stochastic | $\sim 10^{-3}$ | $6.25h$
| approximation |  |
Table 1: Comparison of classical optimization methods for VQA. The accuracy is
given as the magnitude of the MSE, see eq. (52) and the computation time is
the mean value for the computation of one time step.

## References

  * Preskill [2018] J. Preskill,  Quantum computing in the NISQ era and beyond,  Quantum 2 (2018) 79\.
  * Deutsch [2020] I. H. Deutsch,  Harnessing the power of the second quantum revolution,  PRX Quantum 1 (2020) 020101.
  * Nielsen and Chuang [2011] M. Nielsen, I. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition, Cambridge University Press, 2011\.
  * Shor [1997] P. W. Shor,  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,  SIAM J. Comput. 26 (1997) 1484–1509.
  * Grover [1997] L. K. Grover,  Quantum mechanics helps in searching for a needle in a haystack,  Phys. Rev. Lett. 79 (1997) 325–328.
  * Deng et al. [2023] Y.-H. Deng, Y.-C. Gu, H.-L. Liu, S.-Q. Gong, H. Su, Z.-J. Zhang, H.-Y. Tang, M.-H. Jia, J.-M. Xu, M.-C. Chen, J. Qin, L.-C. Peng, J. Yan, Y. Hu, J. Huang, H. Li, Y. Li, Y. Chen, X. Jiang, L. Gan, G. Yang, L. You, L. Li, H.-S. Zhong, H. Wang, N.-L. Liu, J. J. Renema, C.-Y. Lu, J.-W. Pan,  Gaussian boson sampling with pseudo-photon-number-resolving detectors and quantum computational advantage,  Phys. Rev. Lett. 131 (2023) 150601.
  * Choi et al. [2023] S. Choi, W. S. Moses, N. Thompson, The Quantum Tortoise and the Classical Hare: A simple framework for understanding which problems quantum computing will accelerate (and which it won’t), 2023. arXiv:2310.15505.
  * Moin and Mahesh [1998] P. Moin, K. Mahesh,  Direct Numerical Simulation: A tool in turbulence research,  Annu. Rev. Fluid Mech. 30 (1998) 539–578.
  * Iyer et al. [2019] K. P. Iyer, J. Schumacher, K. R. Sreenivasan, P. K. Yeung,  Scaling of locally averaged energy dissipation and enstrophy density in isotropic turbulence,  New J. Phys. 21 (2019) 033016.
  * Buaria et al. [2019] D. Buaria, A. Pumir, E. Bodenschatz, P. K. Yeung,  Extreme velocity gradients in turbulent flows,  New J. Phys. 21 (2019) 043004.
  * Gourianov et al. [2022] N. Gourianov, M. Lubasch, S. Dolgov, Q. Y. van den Berg, H. Babaee, P. Givi, M. Kiffner, D. Jaksch,  A quantum-inspired approach to exploit turbulence structures,  Nat. Comput. Sci. 2 (2022) 30–37.
  * Meng and Yang [2023] Z. Meng, Y. Yang,  Quantum computing of fluid dynamics using the hydrodynamic Schrödinger equation,  Phys. Rev. Research 5 (2023) 033182.
  * Jin et al. [2023] S. Jin, N. Liu, Y. Yu,  Quantum simulation of partial differential equations: Applications and detailed analysis,  Phys. Rev. A 108 (2023) 032603.
  * Succi and Tiribocchi [2023] S. Succi, A. Tiribocchi, Navier-Stokes-Schrödinger equation for dissipative fluids, 2023. arXiv:2308.05879.
  * Pfeffer et al. [2022] P. Pfeffer, F. Heyder, J. Schumacher,  Hybrid quantum-classical reservoir computing of thermal convection flow,  Phys. Rev. Research 4 (2022) 033176.
  * Pfeffer et al. [2023] P. Pfeffer, F. Heyder, J. Schumacher,  Reduced-order modeling of two-dimensional turbulent Rayleigh-Bénard flow by hybrid quantum-classical reservoir computing,  Phys. Rev. Research 5 (2023) 043242.
  * Bharadwaj and Sreenivasan [2020] S. S. Bharadwaj, K. R. Sreenivasan,  Quantum computation of fluid dynamics,  Indian Academy of Sciences Conference Series 3 (2020) 77–96.
  * Bharadwaj and Sreenivasan [2023] S. S. Bharadwaj, K. R. Sreenivasan,  Hybrid quantum algorithms for flow problems,  Proc. Natl. Acad. Sci. USA 120 (2023) e2311014120.
  * Bharadwaj et al. [2023] S. S. Bharadwaj, B. Nadiga, S. Eidenbenz, K. R. Sreenivasan,  Quantum computing of nonlinear flow problems with a homotopy analysis algorithm,  Bull. Am. Phys. Soc. ZC17 (2023) 002.
  * Gaitan [2021] F. Gaitan,  Finding flows in a Navier-Stokes fluid through quantum computing,  npj Quantum Inf. 6 (2021) 61.
  * Lubasch et al. [2020] M. Lubasch, J. Joo, P. Moinier, M. Kiffner, D. Jaksch,  Variational quantum algorithms for nonlinear problems,  Phys. Rev. A 101 (2020) 010301. doi:10.1103/PhysRevA.101.010301.
  * Pool et al. [2022] A. J. Pool, A. D. Somoza, M. Lubasch, B. Horstmann,  Solving partial differential equations using a quantum computer,  2022 IEEE International Conference on Quantum Computing and Engineering (2022) 864–866.
  * Demirdjian et al. [2020] R. Demirdjian, D. Gunlycke, C. A. Reynolds, J. D. Doyle, S. Tafur,  Variational quantum solutions to the advection–diffusion equation for applications in fluid dynamics,  Quantum Inf. Process. 21 (2020) 322.
  * Leong et al. [2022] F. Y. Leong, W.-B. Ewe, D. E. Koh,  Variational quantum evolution equation solver,  Sci. Rep. 12 (2022) 10817.
  * Leong et al. [2023] F. Y. Leong, D. E. Koh, W.-B. Ewe, J. F. Kong,  Variational quantum simulation of partial differential equations: Applications in colloidal transport,  Int. J. Numer. Method H. 33 (2023) 3669–3690.
  * Todorova and de Steijl [2020] B. N. Todorova, R. de Steijl,  Quantum algorithm for the collisionless Boltzmann equation,  J. Comp. Phys. 409 (2020) 109347.
  * Budinski [2021] L. Budinski,  Quantum algorithm for the advection–diffusion equation simulated with the Lattice Boltzmann method,  Quantum Inf. Process. 20 (2021) 57.
  * Succi et al. [2023] S. Succi, W. Itani, K. R. Sreenivasan, R. Steijl,  Quantum computing for fluids: Where do we stand?,  Europhys. Lett. 144 (2023) 10001.
  * Harrow et al. [2009] A. H. Harrow, A. Hassidim, S. Lloyd,  Quantum algorithm for linear systems of equations,  Phys. Rev. Lett. 103 (2009) 150502. doi:10.1103/PhysRevLett.103.150502.
  * Aaronson [2015] S. Aaronson,  Read the fine print,  Nature Physics 11 (2015) 291–293.
  * Montanaro and Pallister [2016] A. Montanaro, S. Pallister,  Quantum algorithms and the finite element method,  Phys. Rev. A 93 (2016) 032324.
  * Guseynov et al. [2023] N. M. Guseynov, A. A. Zhukov, W. V. Pogosov, A. V. Lebedev,  Depth analysis of variational quantum algorithms for the heat equation,  Phys. Rev. A 107 (2023) 052422. doi:10.1103/PhysRevA.107.052422.
  * Liu et al. [2023] Y. Y. Liu, Z. Chen, C. Shu, S. C. Chew, B. C. Khoo, X. Zhao, Y. D. Cui,  Application of a variational hybrid quantum-classical algorithm to heat conduction equation and analysis of time complexity,  Phys. Fluids 34 (2023) 117121.
  * Qis [2023] Qiskit version 0.23.2 (2023).
  * Childs et al. [2017] A. M. Childs, R. Kothari, R. D. Somma,  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision,  SIAM J. Comput. 46 (2017) 1920–1950.
  * Childs et al. [2021] A. M. Childs, J.-P. Liu, A. Ostrander,  High-precision quantum algorithms for partial differential equations,  Quantum 5 (2021) 574.
  * Liu et al. [2021] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, A. M. Childs,  Efficient quantum algorithm for dissipative nonlinear differential equations,  Proc. Natl. Acad. Sci. USA 118 (2021) e2026805118.
  * Cerezo et al. [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, P. J. Coles,  Variational quantum algorithms,  Nat. Rev. Phys. 3 (2021) 625–644.
  * Nelder and Mead [1965] J. Nelder, R. Mead,  A simplex method for function minimization,  The Computer Journal 7 (1965) 308–313. doi:10.1093/comjnl/7.4.308.
  * Barratt et al. [2021] F. Barratt, J. Dborin, M. Bal, V. Stojevic, F. Pollmann, A. G. Green,  Parallel quantum simulation of large systems on small NISQ computers,  npj Quantum Inf. 7 (2021) 79.
  * Shaffer et al. [2023] R. Shaffer, L. Kocia, M. Sarovar,  Surrogate-based optimization for variational quantum algorithms,  Phys. Rev. A 107 (2023) 032415.
  * Dudley et al. [2019] J. M. Dudley, G. Genty, A. Mussot, A. Chabchoub, F. Dias,  Rogue waves and analogies in optics and oceanography,  Nat. Rev. Phys. 1 (2019) 675–689.
  * Brassard et al. [2002] G. Brassard, P. Hoyer, M. Mosca, A. Tapp,  Quantum amplitude amplification and estimation,  Contemp. Math. 305 (2002) 53–74.
  * Engel et al. [2021] A. Engel, G. Smith, S. E. Parker,  Linear embedding of nonlinear dynamical systems and prospects for efficient quantum algorithms,  Phys. Plasmas 28 (2021) 062305.
  * Joseph [2020] I. Joseph,  Koopman–von Neumann approach to quantum simulation of nonlinear classical dynamics,  Phys. Rev. Research 2 (2020) 043102.
  * Giannakis et al. [2022] D. Giannakis, A. Ourmazd, P. Pfeffer, J. Schumacher, J. Slawinska,  Embedding classical dynamics in a quantum computer,  Phys. Rev. A 105 (2022) 052404.
  * Lin et al. [2022] Y. T. Lin, R. B. Lowrie, D. Aslangil, Y. Subaşi, A. T. Sornborger, Koopman von Neumann mechanics and the Koopman representation: A perspective on solving nonlinear dynamical systems with quantum computers, 2022. arXiv:2202.02188.
  * Cross et al. [2019] A. W. Cross, L. S. Bishop, S. Sheldon, P. D. Nation, J. M. Gambetta,  Validating quantum computers using randomized model circuits,  Phys. Rev. A 100 (2019) 032328.
  * Uvarov and Biamonte [2021] A. V. Uvarov, J. D. Biamonte,  On barren plateaus and cost function locality in variational quantum algorithms,  Journal of Physics A: Mathematical and Theoretical 54 (2021). doi:10.1088/1751-8121/abfac7.
  * Liu and Nocedal [1989] D. Liu, J. Nocedal,  On the limited memory BFGS method for large scale optimization,  Mathematical Programming 45 (1989) 503–528. doi:10.1007/BF01589116.
  * Mockus et al. [1978] J. Mockus, V. Tiesis, A. Zilinskas,  The application of Bayesian methods for seeking the extremum,  Towards Global Optimization 2 (1978) 117–129.
  * Kingma and Ba [2015] D. K. Kingma, J. L. Ba,  Adam: A method for stochastic optimization,  3rd International Conference for Learning Representations (2015). doi:10.48550/arXiv.1412.6980.
  * Spall [1998] J. C. Spall,  An overview of the simultaneous perturbation method for efficient optimization,  Johns Hopkins Apl Technical Digest 19 (1998).