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

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.

[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 Preskill (2018); Deutsch (2020) 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 Nielsen and Chuang (2011). Already these two properties suggest faster solutions of problems with high computational complexity, as has been demonstrated for operations such as prime number factorization Shor (1997), data search Grover (1997), and data sampling Deng et al. (2023); see ref. Choi et al. (2023) 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 Moin and Mahesh (1998) 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 Iyer et al. (2019); Buaria et al. (2019) 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 Gourianov et al. (2022) or the mapping of specific classical flow problems to a Schrödinger-type quantum dynamics Meng and Yang (2023); Jin et al. (2023); Succi and Tiribocchi (2023). They include also a surrogate modeling of thermally driven flows within quantum machine learning frameworks, such as hybrid quantum-classical reservoir computing Pfeffer et al. (2022, 2023). 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 Bharadwaj and Sreenivasan (2020, 2023); Bharadwaj et al. (2023), quantum amplitude estimation for one-dimensional gas dynamics Gaitan (2021), Variational Quantum Algorithms for the one-dimensional nonlinear Burgers equation Lubasch et al. (2020); Pool et al. (2022), advection-diffusion problems Demirdjian et al. (2020); Leong et al. (2022, 2023), and quantum lattice Boltzmann methods Todorova and de Steijl (2020); Budinski (2021). See also ref. Succi et al. (2023) 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 Harrow et al. (2009) under certain caveats Aaronson (2015); Montanaro and Pallister (2016). A further approach are variational methods. Different versions, such as the variational quantum imaginary time evolution Leong et al. (2023), the Variational Quantum Linear Solver (VQLS) Demirdjian et al. (2020), or the Variational Quantum Algorithm (VQA) Lubasch et al. (2020); Guseynov et al. (2023); Liu et al. (2023), have been used, even for two-dimensional problems, such as the heat equation Liu et al. (2023).

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 Bharadwaj and Sreenivasan (2023)) 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

$$\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

$$\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

$$c(x,t) =\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]$$

$$\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

$$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

$$a_{0} =\frac{1}{L}\int_{-L}^{L}\delta(x)\text{d}x=\frac{1}{L},$$ (6)

$$a_{k} =\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)

$$b_{k} =\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

$$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

$$c(x,t) =\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

$$\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

$$A{\bm{c}}^{m} ={\bm{c}}^{m+1}$$ (13)

$$\text{with}\ \ A =\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

$$\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

$$A{\bm{c}}^{m+1} ={\bm{c}}^{m}$$ (22)

$$\text{with}\ \ A =\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

$${\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 Qis (2023). The QLSA is done with QFlowS, a C++ based simulation package Bharadwaj and Sreenivasan (2023). 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 Nielsen and Chuang (2011). 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 Nielsen and Chuang (2011).

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 Harrow et al. (2009); Childs et al. (2017, 2021) (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 Harrow et al. (2009). 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 Childs et al. (2017), 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 Bharadwaj and Sreenivasan (2023) 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 Liu et al. (2021); Bharadwaj and Sreenivasan (2023). 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 Nielsen and Chuang (2011) 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 Bharadwaj and Sreenivasan (2023)). 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 Bharadwaj and Sreenivasan (2023)),

$$|\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 Bharadwaj and Sreenivasan (2023); Childs et al. (2021), 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 Harrow et al. (2009); Nielsen and Chuang (2011). 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 Bharadwaj and Sreenivasan (2023). 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 Cerezo et al. (2021). 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 Cerezo et al. (2021).

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 Lubasch et al. (2020); Pool et al. (2022). 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

$$|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

$$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. Lubasch et al. (2020), we define

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

$$|c(t)\rangle =\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

$$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

$$C(\lambda_{0},{\bm{\lambda}}) =\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)

$$+\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

$$C(\lambda_{0},{\bm{\lambda}}) =\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]$$

$$+\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

$$\hat{O} =DN^{2}(\hat{S}_{+}-2\mathds{1}+\hat{S}_{-})-U\frac{N}{2}(\hat{S}_{+}-\hat{S}_{-})$$

$$=(\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:

$$C(\lambda_{0},{\bm{\lambda}})=\lambda_{0}^{2} -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{]}$$

$$+\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{]}$$

$$+\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{]}$$

$$+\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{]}$$

$$+\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

$$C(\lambda_{0},{\bm{\lambda}}) =\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{]}$$

$$+\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{]}$$

$$+\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.

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.

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 Nelder and Mead (1965) 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$ m${}^{2}$/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:

$$n_{\rm tot}^{\rm VQA} =\log_{2}(N)+1\,,$$ (53)

$$n_{\rm tot}^{\rm QLSA} =\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. 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 Bharadwaj and Sreenivasan (2023), 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 Bharadwaj and Sreenivasan (2023). 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 Bharadwaj and Sreenivasan (2023), 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. Barratt et al. (2021). 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 Bharadwaj and Sreenivasan (2023).

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) Childs et al. (2017); Bharadwaj and Sreenivasan (2023) (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 Bharadwaj et al. (2023). In the case of VQA, surrogate algorithms for the global minimum search of the cost function have been suggested recently Shaffer et al. (2023). 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 Cerezo et al. (2021). The strength of the VQA might become better visible for nonlinear problems already attempted, e.g. in Lubasch et al. Lubasch et al. (2020). 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 Dudley et al. (2019).

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 Brassard et al. (2002) 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 Liu et al. (2021); Engel et al. (2021), the Koopman operator formalism Joseph (2020); Giannakis et al. (2022); Lin et al. (2022) or the Homotopy analysis method. Apart from these, the quantum volume 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. 15(a). Consequently, the register is in the following state:

$$|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. 15(b). Then, the following states can be found:

$$|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,

$$\langle\tilde{\psi}|\hat{S}_{+}\tilde{\psi}\rangle =(a,b,c,d)\cdot(d,a,b,c)^{T}=ad+ba+cb+dc$$ (58)

$$\langle\tilde{\psi}|\hat{S}_{-}\tilde{\psi}\rangle =(a,b,c,d)\cdot(b,c,d,a)^{T}=ab+bc+cd+da$$ (59)

$$\langle\hat{S}_{+}\tilde{\psi}|\tilde{\psi}\rangle =(d,a,b,c)\cdot(a,b,c,d)^{T}=da+ab+bc+cd$$ (60)

$$\langle\hat{S}_{-}\tilde{\psi}|\tilde{\psi}\rangle =(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

$$\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$,

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

$$\langle\tilde{\psi}|\hat{S}_{++}\tilde{\psi}\rangle =(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$.

(a) $\hat{S}_{-}$ operation.

(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

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

The first Hadamard gate acts on the zeroth qubit such that

$$|q_{1}q_{0}\rangle_{B} =|\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:

$$|q_{1}q_{0}\rangle_{C} =\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:

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

$$=\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

$$p_{0}-p_{1} =\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)

$$=\frac{1}{2}\langle\varphi|\hat{U}^{\dagger}+\hat{U}|\varphi\rangle=\Re\langle\varphi|\hat{U}|\varphi\rangle\,.$$

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 Uvarov and Biamonte (2021), 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) Nelder and Mead (1965) 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) Liu and Nocedal (1989) 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 Mockus et al. (1978) 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) Kingma and Ba (2015) 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) Spall (1998) 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.