Time complexity analysis of quantum algorithms via linear representations for nonlinear ordinary and partial differential equations

Shi Jin [email protected] School of Mathematical Sciences, Institute of Natural Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R. China. Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China. Ministry of Education, Key Laboratory in Scientific and Engineering Computing, Shanghai Jiao Tong University, Shanghai 200240, China Nana Liu [email protected] Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China. Ministry of Education, Key Laboratory in Scientific and Engineering Computing, Shanghai Jiao Tong University, Shanghai 200240, China University of Michigan-Shanghai Jiao Tong University Joint Institute, Shanghai 200240, China Yue Yu ¹¹1Corresponding author. [email protected] School of Mathematical Sciences, Institute of Natural Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R. China.

Abstract

We construct quantum algorithms to compute the solution and/or physical observables of nonlinear ordinary differential equations (ODEs) and nonlinear Hamilton-Jacobi equations (HJE) via linear representations or exact mappings between nonlinear ODEs/HJE and linear partial differential equations (the Liouville equation and the Koopman-von Neumann equation). The connection between the linear representations and the original nonlinear system is established through the Dirac delta function or the level set mechanism. We compare the quantum linear systems algorithms based methods and the quantum simulation methods arising from different numerical approximations, including the finite difference discretisations and the Fourier spectral discretisations for the two different linear representations, with the result showing that the quantum simulation methods usually give the best performance in time complexity. We also propose the Schrödinger framework to solve the Liouville equation for the HJE with the Hamiltonian formulation of classical mechanics, since it can be recast as the semiclassical limit of the Wigner transform of the Schrödinger equation. Comparsion between the Schrödinger and the Liouville framework will also be made.

Keywords: Linear representation methods; Liouville representation; Koopman-von Neumann representation; Semiclassical Schrödinger equation; Quantum linear systems algorithms.

1 Introduction

Some of the most important problems in physics, chemistry, engineering, biology and finance are modelled by nonlinear ordinary and partial differential equations (ODEs and PDEs). Prominent examples include climate modelling, aircraft design, molecular dynamics and drug design, deep learning neural networks and mean-field games in mathematical finance. Among the most important examples of nonlinear ODEs include Newton’s equations in molecular dynamics, while examples of nonlinear PDEs include the Euler and Navier-Stokes equations in fluid dynamics, the Boltzmann equations in rarified gas, and the Hamilton-Jacobi equations in geometric optics, front propagation, mean-field games and optimal control. In spite of tremendous progresses in developing classical algorithms for solving these equations, there remain major challenges that are difficult to be handled by classical algorithms, for examples the curse-of-dimensionality, multiple scales, strong nonlinearity, and large numbers of degree of freedoms. On the other hand, quantum algorithms, due to their potential polynomial and even exponential advantages, could be the game changer to deal with some of these difficulties. However, since quantum algorithms are based on the principle of quantum mechanics, which is fundamentally linear (as far as we know), so far the development of quantum algorithms have been mostly confined to linear problems, with the most notable in linear algebra [62, 57, 60, 59, 22, 65, 33, 3, 14, 20, 66, 30, 51, 28]. For linear ODEs and PDEs, once they are numerically discretised, they became linear algebra problems which can then be handled by quantum linear algebra solvers (e.g. [19, 15, 21, 55, 26, 11, 24]).

Since most natural phenomena are nonlinear, the ability of quantum computing to solve nonlinear problems will significantly extend the horizon of quantum computing. The most natural idea of handling nonlinear problems by a quantum computer is to represent the nonlinear problem in a linear way, so quantum algorithms for linear problems can be used. There are two approaches in recent literatures. One is to approximate the nonlinear problem through linearisation of the nonlinearity, or truncation of the equation, which is referred to as linear approximation methods. The second is coined as linear representation methods [38], which tries to find a map from the linear to the nonlinear systems, and usually yields a system in the phase space which is equivalent (without any approximation) to the original system.

In the linear approximation approach, in which one linearises the nonlinear system, modeling errors are inevitable, hence the approach may only be valid for a short time, weak or special nonlinearities, and consequently after a long time solutions may lose significant nonlinear features. A more appealing approach uses the Carlemann linearisation in [58] or the techniques in [61, 50], but they are restricted to polynomial nonlinearities, and need to truncate the system since only finite number of equations is used. Such a truncation may lead to loss of some important nonlinear features of the original system. It is similar to moment closure technique in kinetic theory, in which one attempts to close the moment system with finite number of moments but often ends up with a closed system that has mathematical stability problems [32, 8] or physical realizibility issues [49]. This is similarly true for methods in [61, 50], see [54].

For the linear representation methods, there are two current approaches. One is the Koopman-von Neumann (KvN) approach [44] (see also [23]), which was introduced for nonlinear ODEs. The other is the level set method, first introduced in [38], applicable to both nonlinear ODEs, and nonlinear PDEs – more specifically the Hamilton-Jacobi equations and scalar nonlinear hyperbolic equations. Both approaches introduce equations in the phase space but the extra dimensionality that is difficult for classical algorithms can be significantly eased by quantum algorithms, hence quantum advantages in most numerical parameters can still be achieved, even including the measurements of physical observables as analysed in [38].

The goal of this paper is to compare the two linear representation methods – the Liouville equation in the level set approach and the KvN equation, their variants and their different numerical approximations. The KvN is, so to speak, the square of the Liouville equation, and, when the solution is smooth, they are equivalent if the force is divergence free (for example in the case of Hamilton system). While both equations are of linear transport nature, which can be solved similarly by classical algorithms, the KvN equation, due to its unitary structure, can be directly solved by quantum simulations, while the Liouville equation was solved by quantum linear system solvers in [38]. Despite the absence of the unitary structure, we still proposed a “quantum simulation” algorithm for the Liouville representation in Appendix A by using the dimensional splitting Trotter based approximation. The basic idea is to transform the asymmetric evolution in each direction into a symmetric one, which requires only a simple variable substitution with the transformation matrix being diagonal. However, different from the traditional time-marching Hamiltonian simulation, non-unitary procedures for the variable substitution are involved, which may lead to exponential increase of the cost arising from multiple copies of initial quantum states at every time step as pointed out in Remark A.1.

The connection between the Liouville or KvN equation with the original nonlinear system can be made through the Dirac delta-function $\delta(x)$ , which is naturally defined in the weak sense for the Liouville equation which solves for the probability distribution of particles. However, to connect the KvN model – which computes, so to speak, the square root of the probability density distribution, one needs to use $\sqrt{\delta}$ , which is not well-defined mathematically, even in the weak sense, so one needs to be more careful in interpreting its solution (as will be discussed in Section 3.1.2) and the consequent numerical convergence in a suitable solution space. In addition, since the KvN is not in conservation form, it also requires higher smoothness for the force field than the Liouville approach.

Another alternative approach to the Liouville equation (in the classical Hamiltonian formalism) is to use the linear Schrödinger equation, which approximates the Liouville equation in the classical limit by sending the “Planck constant” to zero. Here the Planck constant is an artificial small parameters, which can be chosen to depend on $\varepsilon$ , the numerical precision. The Schrödinger equation can be solved by quantum simulation techniques [37, 62], or by quantum linear algebraic solvers after spatial and temporal discretisations. We will compare all these different models and their different approximations: spectral methods vs. finite difference methods; and quantum simulation vs. quantum linear algebra solvers, and identify the pros and cons of each of these solvers in order to identify the best possible linear representation method.

Our results on time complexities for the computation of physical observables are summarised in Tab. 1. For nonlinear ODEs, the quantum simulation method has less computational cost than the quantum linear systems algorithm (QLSA) based approaches for both Liouville and KvN representations, and the cost of computing the physical observables for the Liouville representation has a multiplicative factor squared times as larger as the one for the KvN representation (if the number of copies needed of the initial state is neglected).

For nonlinear Hamilton-Jacobi equations (in the classical Hamiltonian formalism), the Schrödinger framework/representation, using quantum simulation, has advantages over other approaches in time complexity, for both $d$ (the space dimension) and $\varepsilon$ . However, we would like to point out that the solution to the Schrödinger equation has oscillations of frequency of $\mathcal{O}(1/\sqrt{\varepsilon})$ , as shown in Fig. 2. If one wants high resolution results – for example oscillations free, then the Liouville representation with spectral approximation and QLSA has edges on $d$ , and $\varepsilon$ for smooth solution, while the Liouville replantation with finite difference approximation and QLSA has edges in $\varepsilon$ if solution is less smooth (in Sobolev space $H^{l}$ for $l\leq 4$ ).

For general scalar hyperbolic equation, one can still use the Liouville representation [38] but the Schrödinger representation is not available. It is also unclear how to devise the KvN approach, as will be discussed in Section 6.

We also note that we do not compare these quantum alrogithms with the corresponding classical algorithms, which were partly made, for example, in [38]. In particular, for nonlinear ODEs, these algorithms could be more expensive than the classical solvers [38].

It should be pointed out that the optimal QLSA with query complexity $Q=\mathcal{O}(s\kappa\log(1/\varepsilon))$ is presented in [20], where $s$ is the sparsity, meaning it has at most $s$ nonzero entries per row and column, $\kappa$ is the condition number of the coefficient matrix. On the other hand, the gate complexity may be quantified by $\mathcal{O}(Q\text{poly}(\log Q,\log N))$ , which is larger than the query complexity only by logarithmic factors [14, 52, 20]. For simplicity, we use $\widetilde{\mathcal{O}}(Q)$ to denote the case where all logarithmic factors are suppressed.

The outline of the paper is as follows. In Section 2 we give a more detailed review of the several techniques for solving the nonlinear differential equations. In Section 3 we summarize the general framework for computing the nonlinear ODEs by using the linear representation approaches, and analyse in detail the time complexity of the QLSA based methods and the quantum simulation methods. We propose a simple algorithm to compute the observables of integral form and figure out the multiplicative factor for the sampling procedure. In Section 4 we propose and analyse several quantum algorithms for solving the nonlinear Hamilton-Jacobi PDEs. By using the level set mechanism proposed in [38], we map the nonlinear PDEs of $(d+1)$ -dimension to a linear $(2d+1)$ -dimensional Liouville equation, referred to as the Liouville representation for nonlinear PDEs, and repeat the analysis of the quantum algorithms for the associated Liouville equation as in the previous section. In Section 5, we propose a Schrödinger framework to solve the nonlinear Hamilton-Jacobi PDEs (in the classical Hamiltonian formalism) since the Schrödinger equation can be transformed into the quantum Liouville equation via the Wigner transform, which in turn leads to the Liouville equation when taking the semiclassical limit. We present the quantum interpretation of the classical time-splitting Fourier spectral method proposed for the Schrödinger equation, and give a detailed discussion about the associated sampling law and the gate complexity for the computation of the expectation of observables. In Section 6 we briefly study scalar nonlinear hyperbolic equations. Discussion and summary are given in Section 7.

2 Review of quantum nonlinear differential equation solvers

In [50] a quantum algorithm, perhaps the first nonlinear differential equation solver, is described to solve sparse systems of nonlinear differential equations whose nonlinear terms are quadratic polynomials:

\frac{\mathrm{d}\bm{z}(t)}{\mathrm{d}t}=F(\bm{z}),\qquad f_{\alpha}=\sum\limits_{k,l}a_{kl}^{(\alpha)}z_{k}z_{l},

where $\bm{z}(t)=[z_{1}(t),z_{2}(t),\cdots,z_{n}(t)]^{T}$ and $F=[f_{1}(\bm{z}),f_{2}(\bm{z}),\cdots,f_{n}(\bm{z})]^{T}$ . In the algorithm at least two copies of the state vector are required at every iteration step, which means the complexity scales exponentially with the number of steps to represent the nonlinearity throughout the evolution. Note that polynomials of degree higher than two (or more general nonlinearities) can be reduced to the quadratic case by introducing additional variables. However, the approach is only effective for low-order polynomials and is unlikely to be suitable for practical ODE solvers. Ref. [58] presents a quantum Carleman linearisation algorithm for a class of quadratic nonlinear differential equations

\frac{\mathrm{d}u}{\mathrm{d}t}=F_{2}u^{\otimes 2}+F_{1}u+F_{0}(t),\qquad u^{\otimes 2}=[u_{i}u_{j}].

Compared to the approach of [50], their algorithm improves the complexity from an exponential dependence on the evolution time to a nearly quadratic dependence. This algorithm is efficient in the regime $R<1$ , where the quantity $R$ characterizes the relative strength of the nonlinear and dissipative linear term. As an alternative, Ref. [67] uses a different truncation method - the homotopy perturbation method - to transform the problem into a series of nonlinear ODEs with a special structure, which are further embedded into linear ODEs with a technique similar to the Carleman linearisation. The dependence on the error is exponentially better than the quantum Carleman linearisation algorithm, but the ratio of nonlinearity to dissipation is required to be smaller. Moreover, the algorithm is only discussed for homogeneous differential equations. In [46], the author presents substantially generalised and improved quantum algorithms over prior work for inhomogeneous linear ODEs and applies it to nonlinear ODEs using Carleman linearisation. Ref. [27] describes a strategy of embedding the nonlinear system into an infinite-dimensional linear system with several specific mappings presented, which is then truncated to finite dimension. In particular, the linear Carleman embedding is also included in this approach and connections between several embedding techniques are discussed. Ref. [61] describes a quantum algorithm, referred to as the quantum nonlinear Schrödinger linearisation technique, to solve the following nonlinear ODEs

\frac{\mathrm{d}f}{\mathrm{d}t}+f(x)x=b(t),

where $x\in\mathbb{C}^{d}$ and $f(x)$ is a $d\times d$ matrix that is an order $m$ polynomial function of the vectors $x$ and $x^{\dagger}$ . After technical treatment, $f(x)$ can be assumed in the form of $f(x)=x^{{\dagger}\otimes m}Fx^{\otimes m}$ for a suitable tensor $F$ . In the approach, the forward Euler discretisation $x(t+\Delta t)=x-\Delta tf(x)x$ is approximated by $N$ -body Schrödinger evolution:

x-\Delta tf(x)x\approx\mathrm{e}^{-iH\Delta t}x^{\otimes N},\qquad H=-i\begin{pmatrix}N\\ m\end{pmatrix}^{-1}\sum\limits_{j_{1}\cdots j_{m}}F_{j_{1}\cdots j_{m}},

with

F_{j_{1}\cdots j_{m}}=I\otimes\cdots\otimes F\otimes I\otimes\cdots\otimes F\otimes\cdots\otimes I,

where the right-hand side has $N$ matrices and $F$ s are placed in the $j_{l}$ -th positions.

The above methods can be categorised as linear approximation approaches as in the introduction, in which the nonlinear term is linearised, so the approach may only be valid for a short time, for weak nonlinearities and consequently may lose significant nonlinear features of the problem after a long time [38]. The linear approximation approach is very different from the linear representation method for the ODE, where the original nonlinear ODE is transformed into a linear counterpart by using the level set mechanism without any modelling errors, and for general nonlinearity. On the other hand, when reduced to the linear representations, the impact of the nonlinearity will be reflected in the CFL condition (or the numerical stability condition) as well as the accumulation of the local truncation error for the upwind discretizations. For example, the CFL condition in Theorem 3.1 and the truncation error in Remark 3.6.

3 Linear representation methods for nonlinear ODEs: Liouville vs. KvN representations

In this section we consider the following nonlinear ODEs

\frac{\mathrm{d}q_{j}}{\mathrm{d}t}=F_{j}(q_{1},\cdots,q_{d}),\quad q_{j}(0)=q_{0,j},\quad j=1,\cdots,d,

where $q_{0,j}$ are initial data and $F_{j}$ are real-valued, which can be written in vector form as

\frac{\mathrm{d}q(t)}{\mathrm{d}t}=F(q(t)),\quad q(0)=q_{0},\quad q=[q_{1},\cdots,q_{d}]^{T}.

(3.1)

3.1 Linear representation methods

3.1.1 The Liouville representation

For $x=(x_{1},\cdots,x_{d})$ , let $\delta(x)=\Pi_{i=1}^{d}\delta(x_{i})$ be the Dirac delta distribution. The Liouville equation corresponding to (3.1) can be derived by considering a function $\rho(t,x):\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}$ , defined by

\rho(t,x)=\delta(x-q(t)),

(3.2)

which represents the probability distribution in space $x$ that corresponds to the solution $x=q$ . By the properties of the delta function, one obtains the solution of (3.1) by taking the moment:

q(t)=\int x\delta(x-q(t))\mathrm{d}x=\int x\rho(t,x)\mathrm{d}x.

(3.3)

Other quantity of interest $G(q(t))$ can be obtained by the moment

G(q(t))=\int G(x)\delta(x-q(t))\mathrm{d}x=\int G(x)\rho(t,x)\mathrm{d}x.

(3.4)

To this end, we can characterize the dynamics of $\rho(t,x)$ and find the solution $q(t)$ via (3.3). One can check that $\rho$ satisfies, in the weak sense, the linear $(d+1)$ -dimensional PDE [23, 38]

\begin{cases}\partial_{t}\rho(t,x)+\nabla_{x}\cdot[F(x)\rho(t,x)]=0,\\ \rho_{0}(x):=\rho(0,x)=\delta(x-q_{0}).\end{cases}

(3.5)

Since the initial data involves a delta function, one can consider the following problem with the smoothed initial data

\begin{cases}\partial_{t}\rho^{\omega}(t,x)+\nabla\cdot[F(x)\rho^{\omega}(t,x)]=0,\\ \rho^{\omega}_{0}(x):=\rho^{\omega}(0,x)=\delta_{\omega}(x-q_{0})\end{cases}

(3.6)

when solving (3.5) by a classical or quantum algorithm, where $\omega$ is a smoothing parameter of the delta function [38]. For example, in one-dimensional case, one can choose

\delta_{\omega}=\begin{cases}\frac{1}{\omega}\beta(x/\omega),\quad&|x|\leq\omega,\\ 0,\quad&|x|>\omega,\end{cases}

where typical choices of $\beta$ include

\beta(x)=1-|x|\qquad\mbox{and}\qquad\beta(x)=\frac{1}{2}(1+\cos(\pi x)).

Here $\omega=mh$ and $m$ is the number of mesh points within the support of $\delta_{\omega}$ . For $d$ dimensions, one defines

\delta_{\omega}(x)=\Pi_{i=1}^{d}\delta_{\omega}(x_{i}),\quad x=(x_{1},\cdots,x_{d}).

In addition, the periodic boundary conditions can be imposed since $\rho(0,x)$ or $\rho^{\omega}(0,x)$ has compact support and the solutions to problems (3.5) and (3.6) propagate with finite speed.

To compare with the KvN approach to be introduced next, note the equation in (3.5) is usually transformed to the (some classical) analogue of the Schrödinger equation

\mathrm{i}\partial_{t}\rho=L\rho,

(3.7)

where $L$ is referred to as the Liouville operator, satisfying

L\rho:=-\mathrm{i}\sum\limits_{j}\Big{(}F_{j}\frac{\partial}{\partial x_{j}}+\frac{\partial F_{j}}{\partial x_{j}}\Big{)}\rho.

The operator $L$ is generally not a Hermitian operator, and thus cannot be directly simulated by quantum Hamiltonian unless ${\rm div}F=0$ (In this case one has

L\rho:=-\mathrm{i}\sum\limits_{j}\Big{(}F_{j}\frac{\partial}{\partial x_{j}}+\frac{\partial F_{j}}{\partial x_{j}}\Big{)}\rho=-\mathrm{i}\sum\limits_{j}\Big{(}F_{j}\frac{\partial}{\partial x_{j}}+\frac{1}{2}\frac{\partial F_{j}}{\partial x_{j}}\Big{)}\rho,

which is the symmetric KvN operator defined later).

3.1.2 The Koopman-von Neumann representation

The starting point of the Koopman-von Neumann (KvN) approach to classical mechanics is the introduction of a Hilbert space of complex and square integrable functions $\psi$ referred to as the KvN wave functions such that $\rho:=|\psi|^{2}$ can be interpreted as the probability density of finding a particle at the point of the phase space. To this end, it is desirable to find a complex-valued function $\psi$ which satisfies a dynamical behavior similar to that of the Schrödinger equation, and such that $\rho=|\psi|^{2}=\psi^{\dagger}\psi$ is a solution of (3.7).

Formally, one can verify that $\psi$ satisfies [9, 10, 44, 54]

\mathrm{i}\partial_{t}\psi=\mathcal{H}_{\text{KvN}}\psi,

where the KvN operator is

\mathcal{H}_{\text{KvN}}=-\mathrm{i}\sum\limits_{j}\Big{(}F_{j}\frac{\partial}{\partial x_{j}}+\frac{1}{2}\frac{\partial F_{j}}{\partial x_{j}}\Big{)}.

It can be written in the following symmetric form

\mathcal{H}_{\text{KvN}}\psi=\frac{1}{2}\sum\limits_{j}\Big{(}F_{j}(x)P_{j}\psi+P_{j}(F_{j}(x)\psi)\Big{)}=:\sum\limits_{j=1}^{d}H_{j}\psi,

where

H_{j}\psi=\frac{1}{2}\Big{(}F_{j}(x)P_{j}\psi+P_{j}(F_{j}(x)\psi)\Big{)},\quad P_{j}=-\mathrm{i}\frac{\partial}{\partial x_{j}}.

By introducing the position operator $\hat{x}_{j}$ and the momentum operator $\hat{P}_{j}$ , the KvN representation can be rewritten as

\mathrm{i}\partial_{t}|\psi\rangle=\frac{1}{2}\sum\limits_{j}(F_{j}(\hat{x})\hat{P}_{j}+\hat{P}_{j}F_{j}(\hat{x}))|\psi\rangle=:\sum\limits_{j=1}^{d}\hat{H}_{j}|\psi\rangle,

(3.8)

where the notation $F_{j}(\hat{x})$ denotes a nonlinear map of the position operator which resembles the nonlinear flow, potentially through series expansions. Unlike the Liouville operator, the KvN operator is Hermitian and thus allows for quantum Hamiltonian simulations.

To obtain the quantity of interest $G(q(t))$ one uses

G(q(t))=\int G(x)\delta(x-q(t))\mathrm{d}x=\int G(x)|\psi(t,x)|^{2}\mathrm{d}x.

(3.9)

The setting of Eq. (3.2) may be essential, since it allows to determine the solution of the ODEs via (3.3). To be consistent, one may need to set $\psi(t,x)=\delta^{1/2}(x-q(t))$ , the square root of the delta function. Namely, one needs to solve

\begin{cases}\mathrm{i}\partial_{t}\psi=\mathcal{H}_{\text{KvN}}\psi,\\ \psi(0,x)=\delta^{1/2}(x-q_{0}).\end{cases}

(3.10)

The solution to the above problem, in general, cannot be defined mathematically since one cannot define $\delta^{1/2}$ , let alone its derivative, even in the weak sense. However, by connecting it with (3.5) by $\rho=|\psi|^{2}$ one can make sense of the solution to (3.10) since $|\psi|^{2}$ satisfies the Liouville equation (3.5) which can be defined in the weak sense [56, 29]. This connection is quite important, especially if one discretises (3.10) and hopes the numerical solution will converge, as will be elaborated more later on.

At the discrete level, as was done in [23], one could instead consider smoothing the $\delta$ -function and obtain

\begin{cases}\mathrm{i}\partial_{t}\psi^{\omega}=\mathcal{H}_{\text{KvN}}\psi^{\omega},\\ \psi^{\omega}(0,x)=\delta_{\omega}^{1/2}(x-q_{0}),\end{cases}

(3.11)

which satisfies

\rho^{\omega}=|\psi^{\omega}|^{2}=(\psi^{\omega})^{\dagger}\psi^{\omega}.

(3.12)

Although for every fixed $\omega$ , the problem (3.11) is well defined, it will not be possible to define the $\omega\to 0$ limit which is needed to make mathematical sence of the KvN equation with initial data $\sqrt{\delta(x-q_{0})}$ . Furthermore, when discretising (3.11) numerically, one will not be able to prove the convergence of the numerical approximation (when $\omega$ , mesh size and time step all go to zero) toward the solution of the KvN since the solution of KvN is not suitably defined. In such a case, a different choice of $\omega$ or mesh size could pick up different numerical solutions, which is a well-known phenomenon when numerically approximating the (weak) solution of hyperbolic conservation laws [48]. These solutions may not be unique, unless stronger conditions, such as the entropy condition, are satisfied by the numerical approximations.

Even if one disregards the $\sqrt{\delta}$ problem, when using the Liouville equation (3.5), one only needs $F$ to be Lipschitz continuous to define its solution, like for the original ODEs (3.1). For the KvN equation, since it is not in conservation form, one needs ${\text{div}}F$ to be Lipschitz continuous. Thus this is more restrictive than the case of the Liouville representation.

Unlike in the Liouville representation, in the KvN framework one cannot get the ensemble average defined as $\frac{1}{M_{0}}\sum_{j=1}^{M_{0}}a_{j}$ for $M_{0}$ different initial data [38], where $a(x)$ is a quantity of interest (say physical observables). Rather it gives $(\frac{1}{M_{0}}\sum_{j=1}^{M_{0}}\sqrt{a_{j}})^{2}$ , Thus one loses quantum advantage in $M_{0}$ if one is interested in the ensemble average with $M_{0}\gg 1$ initial data. These two kinds of ensemble average observables are relevant to the quantum sampling setting, referring to [34, 17].

3.1.3 Computation of the quantity of interest and error analysis

After regularising the $\delta$ function by $\delta_{\omega}$ , one now needs to approximate the physical quantity of interest, which are

•

For the Liouville approximation:

$\langle G_{\rho^{\omega}}(t)\rangle=\int G(x)\rho^{\omega}(t,x)\mathrm{d}x,$ (3.13)

where $\rho^{\omega}$ is the solution of (3.6);
•

For the KvN approximation:

$\langle O_{\psi^{\omega}}(t)\rangle=\int G(x)|\psi^{\omega}(t,x)|^{2}\mathrm{d}x,$ (3.14)

where $\psi^{\omega}$ is the solution of (3.11) satisfying $\rho^{\omega}=|\psi^{\omega}|^{2}$ for every fixed $\omega>0$ ,

by some quadrature rules.

By choosing specific $G$ one recovers physical quantities of interest (usually observables). For example, if $G(x)=x$ , (3.13) and (3.14) gives $q$ , the solution of the original equation (3.1). For other choices of $G$ one could recover other physical quantities such as momentum and energy for PDE problems such as the Hamilton-Jacobi equations. See [38] and Eq. (4.8).

For this quantity, we have two ways to approximate it. For the Liouville approximation, one can compute the integral by using the numerical quadrature rule

	$\displaystyle\langle G(t_{n})\rangle$	$\displaystyle=\int_{[0,1]^{d}}G(x)\rho(t_{n},x)\mathrm{d}x\approx\int_{[0,1]^{d}}G(x)\rho^{\omega}(t_{n},x)\mathrm{d}x$
		$\displaystyle=:\langle G_{\rho^{\omega}}(t_{n})\rangle\approx\frac{1}{M^{d}}\sum_{\bm{j}}G_{\bm{j}}\rho^{\omega}_{\bm{j},n}=:\langle G_{\rho^{\omega},n}\rangle,$		(3.15)

where, $G_{\bm{j}}=\bm{\omega}_{\bm{j}}G(x_{\bm{j}})$ with $\bm{\omega}_{\bm{j}}$ being the weight, $\bm{j}=(j_{1},\cdots,j_{d})$ and $M$ is the number of points in each dimension of the $d$ -dimensional space. Throughout the paper, we only consider the trapezoidal rule: with $w=[\frac{1}{2},1,\cdots,1,\frac{1}{2}]^{T}$ , the weight vector can be arranged as $\sum\nolimits_{\bm{j}}\bm{\omega}_{\bm{j}}|\bm{j}\rangle=w\otimes\cdots\otimes w$ . Accordingly, the solution vector is denoted as

\bm{\rho}_{n}^{\omega}=\sum\nolimits_{\bm{j}}\bm{\rho}^{\omega}_{\bm{j},n}|\bm{j}\rangle=\sum\limits_{j_{1},\cdots,j_{d}}\bm{\rho}^{\omega}_{j_{1},\cdots,j_{d},n}|j_{1}\rangle\otimes\cdots\otimes|j_{d}\rangle.

That is, the $n_{\bm{j}}$ -th entry of $\bm{\rho}_{n}^{\omega}$ is $\bm{\rho}^{\omega}_{\bm{j},n}$ , with the global index given by

n_{\bm{j}}:=j_{1}2^{d-1}+\cdots+j_{d}2^{0}.

(3.16)

Note that for periodic boundary conditions, $\rho$ can be assumed to be periodic since the solution to the Liouville equation is essentially zero outside a compact support. In such a case, the trapezoidal rule is of spectral accuracy [2].

Lemma 3.1 (Error of the Liouville representation).

Let $\rho^{\omega}$ be the analytical solution of the Liouville representation (3.5) with the smoothed initial data, and $\rho^{\omega}_{h}$ the numerical solution of $\rho^{\omega}$ . Then

e_{\rho}:=|\langle G(t_{n})\rangle-\langle G_{\rho^{\omega}_{h},n}\rangle|\leq C(\omega\mathrm{e}^{t_{n}\|\text{div}F\|_{\infty}}+d\Delta x^{\ell}/\omega^{\ell+1}+e_{\rho,h}),

(3.17)

where $\|\text{div}F\|_{\infty}=\text{sup}_{q}|\text{div}F(q)|$ , $\ell$ is the Sobolev regularity of $\rho^{\omega}$ (namely $\rho^{\omega}\in C^{\ell}$ ), and $e_{\rho,h}=|\bm{\rho}^{\omega}_{n}-\bm{\rho}^{\omega}_{h,n}|$ is the (relative) discretisation error for the linear Liouville equation, given by for examples:

•

For the first-order upwind finite difference scheme, one has

e_{\rho,h}\leq C\Big{(}\frac{\Delta t}{\omega}+\frac{d\Delta x}{\omega^{2}}\Big{)}=\mathcal{O}\Big{(}\frac{d\Delta x}{\omega^{2}}\Big{)}

(3.18)

with the CFL condition $d\lambda=\mathcal{O}(1)$ (for numerical stability), where $\lambda=\Delta t/\Delta x$ .

•

For the Fourier spectral discretisation, one has

$e_{\rho,h}\leq C\Big{(}\frac{\Delta t^{\alpha}}{\omega^{\alpha}}+\frac{d\Delta x^{\ell}}{\omega^{\ell+1}}\Big{)},$ (3.19)

where $\alpha$ depends on the accuracy of the temporal discretisation.

Proof.

The error can be split as

$\displaystyle e_{\rho}$	$\displaystyle=\|\langle G(t_{n})\rangle-\langle G_{\rho^{\omega}_{h},n}\rangle\|=\|\langle G_{\rho}(t_{n})\rangle-\langle G_{\rho^{\omega}_{h},n}\rangle\|$
	$\displaystyle\leq\|\langle G_{\rho}(t_{n})\rangle-\langle G_{\rho^{\omega}}(t_{n})\rangle\|+\|\langle G_{\rho^{\omega}}(t_{n})\rangle-\langle G_{\rho^{\omega},n}\rangle\|+\|\langle G_{\rho^{\omega},n}\rangle-\langle G_{\rho_{h}^{\omega},n}\rangle\|$
	$\displaystyle=:I_{1}+I_{2}+I_{3}.$	(3.20)

For $I_{1}$ , one can apply the the method of characteristics as done in [64, 38]. To do so, we introduce the characteristics of (3.5) as

\begin{cases}\dfrac{\mathrm{d}X}{\mathrm{d}t}=F(X),\quad X\in\mathbb{R}^{d},\\ X(s)=x,\end{cases}

with the solution denoted by $X(t;x,s)$ . Let the Jacobian determination of the map from $x$ to $X$ be

J(t;x,s)=\text{det}\Big{(}\frac{\partial X_{i}}{\partial x_{j}}(t;x,s)\Big{)}.

Then one has

J(t;x,s)>0,\qquad J(t;x,s)=\exp\Big{(}\int_{s}^{t}\nabla\cdot F(X(\sigma;x,s))\mathrm{d}\sigma\Big{)}.

By the method of characteristics (see Eq. (1.11) of [64] or Appendix J of [38]), the solution to (3.5) can be given by

\rho(t,x)=\rho_{0}(X(0;x,t))J(0;x,t)=\rho_{0}(X(0;x,t))\exp\Big{(}-\int_{0}^{t}\nabla\cdot F(X(\sigma;x,t))\mathrm{d}\sigma\Big{)}.

Similarly, the solution to (3.6) is

\rho^{\omega}(t,x)=\rho_{0}^{\omega}(X(0;x,t))J(0;x,t)=\rho_{0}^{\omega}(X(0;x,t))\exp\Big{(}-\int_{0}^{t}\nabla\cdot F(X(\sigma;x,t))\mathrm{d}\sigma\Big{)}.

Therefore,

	$\displaystyle I_{1}$	$\displaystyle=\|\langle G_{\rho}(t_{n})\rangle-\langle G_{\rho^{\omega}}(t_{n})\rangle\|=\Big{\|}\int_{[0,1]^{d}}G(x)(\rho(t_{n},x)-\rho^{\omega}(t_{n},x))\mathrm{d}x\Big{\|}$
		$\displaystyle=\Big{\|}\int_{[0,1]^{d}}G(x)\Big{(}\rho_{0}(X(0;x,t_{n}))-\rho_{0}^{\omega}(X(0;x,t_{n}))\Big{)}J(0;x,t_{n})\mathrm{d}x\Big{\|}$
		$\displaystyle\leq C\omega\max_{x\in[0,1]^{d}}J(0;x,t_{n})\leq C\omega\mathrm{e}^{t_{n}\\|\text{div}F\\|_{\infty}},$

where $\|\text{div}F\|_{\infty}=\text{sup}_{q}|\text{div}F(q)|$ .

The second term is just the error of the quadrature rule, hence $I_{2}\leq Cd\Delta x^{\ell}/\omega^{\ell+1}$ , where the $1/\omega^{\ell+1}$ factor comes from the $\ell$ -th derivative of $\delta_{w}$ . Obviously, $I_{3}\leq Ce_{\rho,h}$ . The final result comes from the standard error analysis for the linear hyperbolic equation [48]. ∎

Remark 3.1.

Note that for the upwind scheme, due to the CFL condition $\Delta t=\mathcal{O}(\Delta x/d)$ , one has

\mathcal{O}(\Delta t/\omega)=\mathcal{O}(\Delta x/d\omega)={\scriptstyle\mathcal{O}}(d\Delta x/\omega^{2}),

hence the second equality in (3.18) holds. For the spectral discretisation in Subsect. A.1, we require that

d\Delta t^{\alpha}/\omega^{\alpha}\sim d\Delta x^{\ell}/\omega^{\ell+1}\sim\varepsilon,

see (A.10) for example, where $d$ for time comes from the dimension splitting in (A.2). This means

\Delta t^{\alpha}/\omega^{\alpha}\leq d\Delta t^{\alpha}/\omega^{\alpha}\lesssim d\Delta x^{\ell}/\omega^{\ell+1}.

That is, we can still combine the two terms on the right hand side of (3.19). We leave it as is in what follows.

For the KvN approximation, the quadrature rule gives

	$\displaystyle\langle O_{\psi^{\omega}}(t_{n})\rangle$	$\displaystyle=\int_{[0,1]^{d}}G(x)\|\psi^{\omega}(t_{n},x)\|^{2}\mathrm{d}x\approx\frac{1}{M^{d}}\sum_{\bm{j}}G_{\bm{j}}\|\psi^{\omega}_{\bm{j},n}\|^{2}$
		$\displaystyle=\frac{1}{M^{d/2}}(\bm{\psi}_{n}^{\omega})^{\dagger}G_{M}\bm{\psi}_{n}^{\omega}=:\langle O_{\psi^{\omega},n}\rangle,$		(3.21)

where the elements of the vector $\bm{\psi}_{n}^{\omega}$ are arranged as $\bm{\psi}_{n}^{\omega}=\sum_{\bm{j}}\psi^{\omega}_{\bm{j},n}|\bm{j}\rangle$ , and $G_{M}=\text{diag}(\bm{g})$ is a diagonal matrix with $\bm{g}=\sum_{\bm{j}}G_{\bm{j}}/M^{d/2}|\bm{j}\rangle$ satisfying $\|\bm{g}\|\sim 1$ .

Lemma 3.2 (Error of KvN representation).

Let $\psi^{\omega}$ be the analytical solution of the KvN representation (3.11), and $\psi^{\omega}_{h}$ the numerical solution of $\psi^{\omega}$ . Then

e_{\psi}:=|\langle G(t_{n})\rangle-\langle O_{\psi^{\omega}_{h},n}\rangle|\leq C(\omega\mathrm{e}^{t_{n}\|\text{div}F\|_{\infty}}+d\Delta x^{\ell}/\omega^{\ell+1}+e_{\psi,h}),

(3.22)

where $\ell$ is the regularity of $\rho^{\omega}$ , and $e_{\psi,h}=||\bm{\psi}_{n}^{\omega}|^{2}-|\bm{\psi}_{h,n}^{\omega}|^{2}|$ is the (relative) discretisation error for the KvN equation (3.8).

•

For the first-order upwind finite difference scheme, one has

$e_{\psi,h}\leq C\Big{(}\frac{d\Delta x}{\omega^{2}}\Big{)},$

with the CFL condition $d\lambda=\mathcal{O}(1)$ , where $\lambda=\Delta t/\Delta x$ .
•

For the Fourier spectral discretisation, one has

$e_{\psi,h}\leq C\Big{(}\frac{\Delta t^{\alpha}}{\omega^{\alpha}}+\frac{d\Delta x^{\ell}}{\omega^{\ell+1}}\Big{)},$

where $\alpha$ depends on the accuracy of the temporal discretisation.

Proof.

Noting that $\rho^{\omega}=|\psi^{\omega}|^{2}$ , we can split the error as

	$\displaystyle e_{\psi}$	$\displaystyle=\|\langle G(t_{n})\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|$
		$\displaystyle\leq\|\langle G_{\rho}(t_{n})\rangle-\langle O_{\psi^{\omega}}(t_{n})\rangle\|+\|\langle O_{\psi^{\omega}}(t_{n})\rangle-\langle O_{\psi^{\omega},n}\rangle\|+\|\langle O_{\psi^{\omega},n}\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|$
		$\displaystyle=\|\langle G_{\rho}(t_{n})\rangle-\langle G_{\rho^{\omega}}(t_{n})\rangle\|+\|\langle O_{\psi^{\omega}}(t_{n})\rangle-\langle O_{\psi^{\omega},n}\rangle\|+\|\langle O_{\psi^{\omega},n}\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|$
		$\displaystyle=:I_{1}+I_{2}+I_{3}.$

Here, $I_{1}$ is exactly the first term in (3.20) for the Liouville representation, so $I_{1}\leq C\omega\mathrm{e}^{t_{n}\|\text{div}F\|_{\infty}}$ . The second term is the (relative) error of the quadrature rule. One again has $I_{2}\leq Cd\Delta x^{\ell}/\omega^{\ell+1}$ since $\rho^{\omega}=|\psi^{\omega}|^{2}$ , where the $1/\omega^{\ell+1}$ factor comes from the $\ell$ -th derivative of $\delta_{w}$ . For the last term, one has

	$\displaystyle I_{3}$	$\displaystyle=\|\langle O_{\psi^{\omega},n}\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|=\|(\bm{\psi}_{n}^{\omega})^{\dagger}G_{M}\bm{\psi}_{n}^{\omega}-(\bm{\psi}_{h,n}^{\omega})^{\dagger}G_{M}\bm{\psi}_{h,n}^{\omega}\|$
		$\displaystyle=\frac{1}{M^{d}}\Big{\|}\sum_{j}G_{\bm{j}}\|\psi_{\bm{j},n}^{\omega}\|^{2}-\sum_{j}G_{\bm{j}}\|(\psi_{h}^{\omega})_{\bm{j},n}\|^{2}\Big{\|}\lesssim\|\|\bm{\psi}_{n}^{\omega}\|^{2}-\|\bm{\psi}_{h,n}^{\omega}\|^{2}\|.$

This completes the proof. ∎

Remark 3.2.

When discretising the Liouville equation by the upwind scheme, since it is in conservative form, one can show that the $l_{1}$ norm at later time is bounded by the norm at $t=0$ when the CFL condition is satisfied (namely it is $l_{1}$ contracting). In addition, the accumulation of the local truncation error which grows linearly in $t$ . However, this is not true for the KvN representation. If one uses the upwind scheme to discretize the transport (spatial derivative) term, since the $1/2$ term is a forcing term, it will contribute to the $l_{1}$ error an exponentially growing term like $\mathrm{e}^{t_{n}\|\text{div}F\|_{\infty}}$ . When $t_{n}\leq T=\mathcal{O}(1)$ this is not a big issue. But it is worth pointing out the difference with the Liouville representation here, since if one wants to compute the long time solution, the $l_{1}$ contraction of the upwind scheme for the Liouville representation gives much smaller error. The above discussion will be further elaborated later.

For simplicity, we only consider $t_{n}\leq T=\mathcal{O}(1)$ in this article.

3.2 Finite difference discretisation for the Liouville representation

3.2.1 The QLSA for the finite difference discretisation

The Liouville representation can be rewritten as

\begin{cases}\frac{\partial w}{\partial t}+\sum\limits_{i=1}^{d}\frac{\partial}{\partial x_{i}}(F_{i}(x)w(t,x))=0,\\ w_{0}(x)=\delta_{\omega}(x-x_{0}),\end{cases}

(3.23)

where we have assumed the smoothed initial data $\delta_{\omega}$ .

Denote $\bm{e}_{i}=[0,\cdots,1,\cdots,0]$ to be the unit vector with the $i$ -th entry being 1. Let $\bm{j}=(j_{1},\cdots,j_{i},\cdots,j_{d})$ and $x_{i}=j_{i}\Delta x$ are the $i$ -th components of $x_{\bm{j}}$ . The first-order upwind discretisation at $(t_{n},x_{\bm{j}})$ takes the following form [38, Appendix K]:

	$\displaystyle\partial_{t}w\longrightarrow\frac{w_{\bm{j}}^{n+1}-w_{\bm{j}}^{n}}{\Delta t},$
	$\displaystyle\frac{\partial}{\partial x_{i}}(F_{i}(x)w(t,x))\longrightarrow\frac{1}{\Delta x}\Big{(}\Big{\{}F_{i}(x_{j_{i}+1/2})\Big{\}}_{\bm{j}}^{-}w^{n}_{\bm{j}+\bm{e}_{i}}-\Big{\{}F_{i}(x_{j_{i}-1/2})\Big{\}}_{\bm{j}}^{+}w^{n}_{\bm{j}-\bm{e}_{i}}\Big{)}$
	$\displaystyle\quad\hskip 85.35826pt+\frac{1}{\Delta x}\Big{(}\Big{\{}F_{i}(x_{j_{i}+1/2})\Big{\}}_{\bm{j}}^{+}-\Big{\{}F_{i}(x_{j_{i}-1/2})\Big{\}}_{\bm{j}}^{-}\Big{)}w_{\bm{j}}^{n},$

where

\alpha^{+}=\max\{\alpha,0\}=\frac{\alpha+|\alpha|}{2}\geq 0,\quad\alpha^{-}=\min\{\alpha,0\}=\frac{\alpha-|\alpha|}{2}\leq 0,

and

\Big{\{}F_{i}(x_{j_{i}\pm 1/2})\Big{\}}_{\bm{j}}=\frac{1}{2}(F_{i}(x_{\bm{j}\pm\bm{e}_{i}})+F_{i}(x_{\bm{j}}))=\frac{1}{2}(F_{i}(x_{j_{1},\cdots,j_{i}\pm 1,\cdots,j_{d}})+F_{i}(x_{j_{1},\cdots,j_{i},\cdots,j_{d}})).

For convenience we introduce the following notation

a_{\bm{j}}^{i,\pm}=\Big{\{}F_{i}(x_{j_{i}+1/2})\Big{\}}_{\bm{j}}^{\pm},\qquad b_{\bm{j}}^{i,\pm}=\Big{\{}F_{i}(x_{j_{i}-1/2})\Big{\}}_{\bm{j}}^{\pm}.

The discrete scheme can be written as

w_{\bm{j}}^{n+1}-\Big{[}1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i,+}-b_{\bm{j}}^{i,-})\Big{]}w^{n}_{\bm{j}}+\lambda\sum\limits_{i=1}^{d}\Big{[}a_{\bm{j}}^{i,-}w^{n}_{\bm{j}+\bm{e}_{i}}-b_{\bm{j}}^{i,+}w^{n}_{\bm{j}-\bm{e}_{i}}\Big{]}=0,

(3.24)

where $\lambda=\Delta t/\Delta x$ . In matrix form one has

\bm{w}^{n+1}-B\bm{w}^{n}=\bm{f}^{n+1},\quad n=0,1,\cdots,N_{t}-1,

with $\bm{f}^{i}$ being the terms resulting from the initial and boundary conditions, where the nodal values at $t=t_{n}$ are arranged as $\bm{w}^{n}=\sum_{\bm{j}}w_{\bm{j}}^{n}|\bm{j}\rangle$ . With the help of the global index (3.16), the non-zero entries of $B$ can be given by

B_{n_{\bm{j}},n_{\bm{j}}}=1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i,+}-b_{\bm{j}}^{i,-}),\qquad B_{n_{\bm{j}},n_{\bm{j}+\bm{e}_{i}}}=-\lambda a_{\bm{j}}^{i,-},\qquad B_{n_{\bm{j}},n_{\bm{j}-\bm{e}_{i}}}=\lambda b_{\bm{j}}^{i,+}.

By introducing the notation $\bm{w}=[\bm{w}^{1};\cdots;\bm{w}^{N_{t}}]$ , where “;” indicates the straightening of $\{\bm{w}^{i}\}_{i\geq 1}$ into a column vector, one obtains the following linear system

L\bm{w}=F,

(3.25)

where

L=\begin{bmatrix}I&&&\\ -B&I&&\\ &\ddots&\ddots&\\ &&-B&I\\ \end{bmatrix},\qquad F=\begin{bmatrix}\bm{f}^{1}\\ \bm{f}^{2}\\ \vdots\\ \bm{f}^{N_{t}}\\ \end{bmatrix}.

For periodic boundary conditions, one has $\bm{f}^{1}=B\bm{w}^{0}$ and $\bm{f}^{i}=\bm{0}$ for $i\geq 2$ .

Theorem 3.1.

Suppose that $\lambda=\Delta t/\Delta x$ satisfies the following CFL condition

\lambda\sum\limits_{i=1}^{d}\sup_{x}|F_{i}(x)|\leq 1.

(1)

The condition number and the sparsity of $L$ satisfy $\kappa=\mathcal{O}(1/\Delta t)$ and $s=\mathcal{O}(d)$ .
(2)

For fixed spatial step $\Delta x$ , let $\Delta t=\mathcal{O}(\Delta x/d)$ and $\omega=(d\Delta x)^{1/3}$ . Given the error tolerance $\varepsilon$ , the gate complexity of the QLSA (for the problem in Eq. (3.25)) is

$N_{\text{Gates}}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3}}{\varepsilon^{3}}\log\frac{1}{\varepsilon}\Big{)}.$

Proof.

1) We claim that $\|B\|_{2}\leq 1+\Delta t\|{\rm div}F\|_{\infty}$ . To this end, one can show that the summation of the absolute values of row or column entries is not greater than $1+\Delta t\|{\rm div}F\|_{\infty}$ .

Let $R_{i}$ be the sum of absolute values of row entries. Since the $i$ -th equation in $Ax=b$ is actually the $i$ -th row of $A$ multiplied by $x$ , the entries of the $i$ -th row of $A$ are completely determined by the $i$ -th equation. On the other hand, the grid values corresponding to the initial or boundary conditions in the $i$ -th equation will increase the values of the summation of the absolute values of row entries since they will be moved to the right-hand side. That is, $R_{i}\leq\tilde{R}_{i}$ , where $\tilde{R}_{i}$ includes the contribution from the initial-boundary values.

With the CFL condition, one has

c:=1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i,+}-b_{\bm{j}}^{i,-})\geq 0.

The argument is as follows. Without loss of generality, we set $d=1$ and obtain

c=1-\frac{\lambda}{4}(F_{j+1}-F_{j-1}+|F_{j+1}+F_{j}|+|F_{j-1}+F_{j}|)=:1-\frac{\lambda}{4}c_{F},

where $F_{j}=F(x_{j})$ . This reduces to verify that $c_{F}$ contains only four terms of the grid values of $F$ :

•

For $F_{j-1}+F_{j}\geq 0$ , one has

$c_{F}=F_{j+1}-F_{j-1}+|F_{j+1}+F_{j}|+F_{j-1}+F_{j}=F_{j+1}+F_{j}+|F_{j+1}+F_{j}|,$

as required.

•

For $F_{j-1}+F_{j}<0$ , one has

	$\displaystyle c_{F}$	$\displaystyle=F_{j+1}-F_{j-1}+\|F_{j+1}+F_{j}\|-F_{j-1}-F_{j}$
		$\displaystyle=\begin{cases}F_{j+1}-F_{j-1}+F_{j+1}-F_{j-1},\qquad&\mbox{if $F_{j+1}+F_{j}\geq 0$},\\ -F_{j-1}-F_{j}-F_{j-1}-F_{j},\qquad&\mbox{if $F_{j+1}+F_{j}<0$},\end{cases}$

as required.

The above argument implies that

	$\displaystyle R_{i}$	$\displaystyle\leq\tilde{R}_{i}\leq 1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i,+}-b_{\bm{j}}^{i,-})+\lambda\sum\limits_{i=1}^{d}(b_{\bm{j}}^{i,+}-a_{\bm{j}}^{i,-})$
		$\displaystyle=1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i,+}+a_{\bm{j}}^{i,-}-b_{\bm{j}}^{i,+}-b_{\bm{j}}^{i,-})=1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i}-b_{\bm{j}}^{i})$
		$\displaystyle=1-\lambda\sum\limits_{i=1}^{d}\Big{(}\Big{\{}F_{i}(x_{i+1/2})\Big{\}}_{\bm{j}}-\Big{\{}F_{i}(x_{i-1/2})\Big{\}}_{\bm{j}}\Big{)}.$

Since $F$ is smooth, we obtain from the mean value theorem that

	$\displaystyle R_{i}$	$\displaystyle\leq 1-\lambda\sum\limits_{i=1}^{d}\Big{(}\Big{\{}F_{i}(x_{i+1/2})\Big{\}}_{\bm{j}}-\Big{\{}F_{i}(x_{i-1/2})\Big{\}}_{\bm{j}}\Big{)}$
		$\displaystyle=1-\lambda\sum\limits_{i=1}^{d}\partial_{x_{i}}F_{i}(\xi_{i})\Delta x\qquad(x_{i-1/2}<\xi_{i}\leq x_{i+1/2})$
		$\displaystyle=1-\Delta t\sum\limits_{i=1}^{d}\partial_{x_{i}}F_{i}(\xi_{i})\leq 1+\Delta t\\|{\rm div}F\\|_{\infty},$

where $F_{i}(\xi_{i}):=F_{i}(x_{j_{1}},\cdots,\xi_{i},\cdots,x_{j_{d}})$ .

Let $C_{j}$ be the absolute column under discussion. The $j$ -th column of $A$ is exactly the collection of the coefficients of $x_{j}$ in each equation of $Ax=b$ , hence the absolute column sum is simply the sum of the absolute values of the coefficients with respect to $x_{j}$ . Consider the variable $w_{\bm{j}}^{n}$ in (3.24). Notice that the other subscripts are only shifted left and right once in some direction, and thus the other elements of the corresponding column are only changed by the upper index $i$ . This again implies that

\displaystyle C_{j}\leq 1-\lambda\sum\limits_{i=1}^{d}(a_{\bm{j}}^{i,+}-b_{\bm{j}}^{i,-})+\lambda\sum\limits_{i=1}^{d}(b_{\bm{j}}^{i,+}-a_{\bm{j}}^{i,-})\leq 1+\Delta t\|{\rm div}F\|_{\infty}.

2) By definition, $\sigma_{\min}(L)=1/\sigma_{\max}(L^{-1})$ . After simple algebra, one has

{\scriptsize L^{-1}=\begin{bmatrix}I&&&&\\ B&I&&&\\ B^{2}&\ddots&\ddots&&\\ \vdots&\ddots&\ddots&\ddots&\\ B^{N_{t}-1}&\cdots&B^{2}&B&I\\ \end{bmatrix}=\begin{bmatrix}I&&&&\\ &I&&&\\ &&\ddots&&\\ &&&\ddots&\\ &&&&I\\ \end{bmatrix}+\begin{bmatrix}&&&&\\ B&&&&\\ &\ddots&&&\\ &&\ddots&&\\ &&&B&\\ \end{bmatrix}+\cdots,}

which gives

\displaystyle\sigma_{\max}(L^{-1})=\|L^{-1}\|_{2}\leq\|I\|_{2}+\|B\|_{2}+\|B\|_{2}^{2}+\cdots+\|B\|_{2}^{N_{t}-1}.

According to the previous analysis, one has

\|B\|_{2}\leq 1+\Delta t\|{\rm div}F\|_{\infty}\leq c=\begin{cases}1,\quad&\|{\rm div}F\|_{\infty}=0,\\ 1+\Delta t,\quad&0<\|{\rm div}F\|_{\infty}<1,\\ 1+\Delta t\|{\rm div}F\|_{\infty},\quad&\|{\rm div}F\|_{\infty}\geq 1,\end{cases}

and hence

\sigma_{\max}(L^{-1})\leq 1+c+c^{2}+\cdots c^{N_{t}-1}=\frac{c^{N_{t}}-1}{c-1}.

Noting that $(1+x/n)^{n}\leq\mathrm{e}^{x}$ holds for any real number, we then have

\sigma_{\max}(L^{-1})\leq\frac{1}{\Delta t}\times\begin{cases}1,\quad&\|{\rm div}F\|_{\infty}=0\\ \mathrm{e},\quad&0<\|{\rm div}F\|_{\infty}<1\\ \exp(\|{\rm div}F\|_{\infty}),\quad&\|{\rm div}F\|_{\infty}\geq 1\end{cases},

which can be simply written as

\sigma_{\max}(L^{-1})\leq\exp(\|{\rm div}F\|_{\infty}+1)\frac{1}{\Delta t},

hence

\sigma_{\min}(L)\geq\frac{\Delta t}{\exp(\|{\rm div}F\|_{\infty}+1)}.

By the Gershgorin-type theorem for singular values [63, 35],

\sigma_{\max}(L)\leq 1+\|B\|_{2}\leq 2+\Delta t\|{\rm div}F\|_{\infty}.

which gives

\kappa(L)\leq(2+\Delta t\|{\rm div}F\|_{\infty})\exp(\|{\rm div}F\|_{\infty}+1)\frac{1}{\Delta t}\lesssim\exp(\|{\rm div}F\|_{\infty})\frac{1}{\Delta t}.

3) In view of the CFL condition, we set $\Delta t=\mathcal{O}(\Delta x/d)$ . According to Lemma 3.1, the classical error of the numerical approximation (if the first order upwind scheme is used to discretize the Liouville equation) to the observables is $\mathcal{O}(\omega+\Delta t/\omega+d\Delta x/\omega^{2})$ . One can choose $\omega$ such that $\omega\sim d\Delta x/\omega^{2}$ or $\omega=(d\Delta x)^{1/3}$ so the numerical error becomes $O((d\Delta x)^{1/3})$ . To reach the precision $\mathcal{O}(\varepsilon)$ , we choose $\Delta x\sim\varepsilon^{3}/d$ , and hence $\Delta t/\omega\sim\varepsilon^{2}/d^{2}\leq\varepsilon$ . This naturally leads to the query complexity [20]

Q=\mathcal{O}\Big{(}s\kappa\log\frac{1}{\varepsilon}\Big{)}=\mathcal{O}\Big{(}\frac{d^{3}}{\varepsilon^{3}}\log\frac{1}{\varepsilon}\Big{)}.

The gate complexity is larger than the query complexity only by logarithmic factors. ∎

Remark 3.3.

Since the Liouville equation is in conservative form, one can get the $l_{1}$ contracting of the upwind scheme. Without loss of generality we set $d=1$ . The upwind scheme in (3.24) is then given by

w_{j}^{n+1}-(1-\lambda(a_{j}^{+}-b_{j}^{-}))w^{n}_{j}+\lambda(a_{j}^{-}w^{n}_{j+1}-b_{j}^{+}w^{n}_{j-1})=0,

where $a_{j}^{\pm}=F(x_{j+1/2})^{\pm}$ and $b_{j}^{\pm}=F(x_{j-1/2})^{\pm}$ . Since $b_{j}^{\pm}=a_{j-1}^{\pm}$ , the scheme can be rewritten as

w_{j}^{n+1}-(1-\lambda(a_{j}^{+}-a_{j-1}^{-}))w^{n}_{j}+\lambda(a_{j}^{-}w^{n}_{j+1}-a_{j-1}^{+}w^{n}_{j-1})=0.

Under the CFL condition given in Theorem 3.1, one has

	$\displaystyle\\|\bm{w}^{n+1}\\|_{1}$	$\displaystyle=\sum\limits_{j}\|w_{j}^{n+1}\|\leq\sum\limits_{j}\Big{[}(1-\lambda(a_{j}^{+}-a_{j-1}^{-}))\|w^{n}_{j}\|-\lambda a_{j}^{-}\|w^{n}_{j+1}\|+\lambda a_{j-1}^{+}\|w^{n}_{j-1}\|\Big{]}$
		$\displaystyle=\sum\limits_{j}\Big{[}(1-\lambda(a_{j}^{+}-a_{j-1}^{-}))\|w^{n}_{j}\|-\lambda a_{j-1}^{-}\|w^{n}_{j}\|+\lambda a_{j}^{+}\|w^{n}_{j}\|\Big{]}$
		$\displaystyle=\sum\limits_{j}\|w_{j}^{n}\|=\\|\bm{w}^{n}\\|_{1},$

as required.

3.2.2 The algorithm for the computation of the observable

The quantum algorithm to approximate physical observables is presented in [38] with a detailed analysis on the gate complexity, where the observable is computed by using the amplitude estimation algorithm with block-encoding techniques augmented by amplitude amplification. Amplitude amplification was used to achieve optimal scaling of the query complexity with respect to the error $\varepsilon$ while measuring an expectation value. However, since our purpose here is only to compare the strengths of the Liouville approach versus the Koopman-von Neumann approach, for simplicity we can instead compute the observable with a more straightforward means, without using amplitude amplification. In this paper, we first obtain the quantum state proportional to the solution of the problem, either with QLSA or with quantum simulation, then compute the observable afterwards.

The expectation of the observable

In order to measure observables, we need to express them as Hermitian operators. Let $\bm{w}_{\bm{j},n}:=\bm{\rho}^{\omega}_{\bm{j},n}$ be the solution of the upwind finite difference method (with the smoothed initial data $\delta_{\omega}$ ). One has

|\psi\rangle=\frac{1}{N_{\psi}}\sum\limits_{\bm{j},n}\bm{w}_{\bm{j},n}|\bm{j}\rangle|n\rangle,

where the normalization constant $N_{\psi}=\|\bm{w}\|$ . With $G_{\bm{j}}$ in (3.1.3), we define the state

|G_{n}\rangle:=\frac{1}{N_{G}}\sum_{\bm{j}}G^{\dagger}_{\bm{j}}|\bm{j}\rangle|n\rangle,

where $N_{G}=(\sum_{\bm{j}}|G^{\dagger}_{\bm{j}}|^{2})^{1/2}$ is the normalisation constant. We assume we are given the density matrix $\mathcal{G}:=|G_{n}\rangle\langle G_{n}|$ and we define $\Upsilon:=\langle\psi|\mathcal{G}|\psi\rangle$ . Simple algebra yields

\langle G(t_{n})\rangle\approx\langle G_{\rho^{\omega},n}\rangle=n_{\psi}n_{G}|\sqrt{\Upsilon}|,

(3.26)

where $n_{G}=N_{G}/M^{d/2}=\mathcal{O}(1)$ is known and $n_{\psi}=N_{\psi}/M^{d/2}$ may be unknown. We further define

\langle O\rangle=\langle G_{\rho^{\omega},n}\rangle^{2}=(n_{G}n_{\psi})^{2}\Upsilon:=\langle\psi|O|\psi\rangle,\qquad O=(n_{G}n_{\psi})^{2}\mathcal{G}.

(3.27)

Then one only needs to estimate it to precision $\varepsilon$ since

|\langle G_{\rho^{\omega},n}\rangle-\langle G_{\rho^{\omega},n}\rangle_{\text{app}}|=\frac{1}{|\langle G_{\rho^{\omega},n}\rangle+\langle G_{\rho^{\omega},n}\rangle_{\text{app}}|}|\langle O\rangle-\langle O\rangle_{\text{app}}|

and $\langle G_{\rho^{\omega},n}\rangle$ and $\langle G_{\rho^{\omega},n}\rangle_{\text{app}}$ can be considered as $\mathcal{O}(1)$ , where the subscript “app” refers to the approximations.

The problem then reduces to approximating the normalisation constant $N_{\psi}=\|\bm{w}\|$ to a desired precision, where $\bm{w}=L^{-1}F$ . This can be referred to as the amplitude estimation or linear equation norm estimation in quantum computing [55, 13]. As shown in Fig. 1, for the QLSA of the upwind discretisation, we first compute an approximation $\tilde{N}_{\psi}$ of $N_{\psi}$ by using amplitude estimation, and then construct the approximate observable with $n_{\psi}$ replaced by $\tilde{n}_{\psi}$ .

Refer to caption — Fig. 1: Construction of the observables when $\|\bm{w}\|$ is unknown.

The general sampling law

Since the measurement outcome is probabilistic in general, we have to evaluate the expectation value via sampling. Let $O$ be an observable with $\mu:=\langle O\rangle=\langle\psi|O|\psi\rangle$ being the expectation value, where $|\psi\rangle$ is a quantum state. Suppose that we conduct $n$ experiments with the outcomes labelled $\mu_{1},\cdots,\mu_{n}$ . By the law of large numbers

{\rm P}_{\text{r}}\Big{(}\Big{|}\frac{\mu_{1}+\cdots+\mu_{n}}{n}-\mu\Big{|}<\varepsilon\Big{)}\geq 1-\frac{\text{Var}(O)}{n\varepsilon^{2}},

where $\text{Var}(O)$ is the variance. For a given lower bound $p$ , the number of samples required to estimate $\langle O\rangle$ to additive precision $\varepsilon$ satisfies

1-\frac{\text{Var}(O)}{n\varepsilon^{2}}\geq p\quad\Longrightarrow\quad n\geq\frac{1}{1-p}\frac{\text{Var}(O)}{\varepsilon^{2}}.

This implies a multiplicative factor $\text{Var}(O)/\varepsilon^{2}$ in the total gate complexity [44, 51], which is referred to as the “general sampling law” in this article. We remark that in many cases of interest, the number of repetitions can be reduced to $\mathcal{O}(1/\varepsilon)$ , up to polylogarithmic factors. For example, quantum algorithms based on amplitude amplification and estimation are able to compute numerical approximations to sums and integrals with a quadratic speedup over classical probabilistic algorithms, so that the number of repetitions of the quantum simulation follows the “quantum sampling law” $\mathcal{O}(1/\varepsilon)$ , up to polylogarithmic factors [44]. When $O=U$ is a unitary gate acting on the space of $|\psi\rangle$ , one can apply the Hadamard test to create a random variable whose expected value is the expected real part $\text{Re}\langle\psi|U|\psi\rangle$ or imaginary part $\text{Im}\langle\psi|U|\psi\rangle$ [1]. The above sampling law is still valid. Such a technique is also employed in [25, 53] for ground-state preparation and energy estimation, where $U=\mathrm{e}^{-\mathrm{i}\tau H}$ is the time evolution operator. In this paper, we only assume the general sampling law $\mathcal{O}(\text{Var}(O)/\varepsilon^{2})$ . Noting that

{\rm Var}(O)=(n_{G}n_{\psi})^{4}{\rm Var}(\mathcal{G})\lesssim n_{\psi}^{4},

we may need to include this multiplicative factor $n_{\psi}^{4}$ in the query complexity. Below we will show that in fact we can replace this factor instead by $n_{\psi_{0}}$ where $n_{\psi_{0}}=N_{\psi_{0}}/M^{d/2}=\|\bm{w}^{0}\|/M^{d/2}$ as defined in [38].

For the QLSA of the upwind discretisation, referring to the linear system (3.25), and noting that $\|B\|\lesssim 1+\Delta t$ and $\|L^{-1}\|=\sigma_{\max}(L^{-1})\lesssim N_{t}$ , we have

N_{\psi}=\|\bm{w}\|=\|L^{-1}F\|\leq N_{t}\|F\|\lesssim N_{t}\|\bm{w}^{0}\|,\qquad\mbox{and}\qquad n_{\psi}\lesssim N_{t}n_{\psi_{0}},

which gives the multiplicative factor $N_{t}^{4}n_{\psi_{0}}^{4}$ . However this $N_{t}^{4}$ factor can be removed as addressed in [5, 51] by adding $N_{t}$ copies of the final state $\bm{w}^{N_{t}}$ . That is, we add the following additional equations

\bm{w}^{n+1}-\bm{w}^{n}=0,\quad n=N_{t},\cdots,2N_{t}

(3.28)

in (3.25), which is referred to as the dilation procedure. For simplicity, we assume that $\|\bm{w}^{1}\|=\cdots=\|\bm{w}^{N_{t}}\|=\|\bm{w}^{0}\|=N_{\psi_{0}}$ . Let the padded state vector be

\widehat{\bm{w}}=[\widehat{\bm{w}}^{1};\cdots;\widehat{\bm{w}}^{N_{t}};\widehat{\bm{w}}^{N_{t}},\cdots,\widehat{\bm{w}}^{N_{t}}]=|0\rangle\otimes\bm{x}+|1\rangle\otimes\bm{y},

(3.29)

where $|0\rangle=[1,0]^{T}$ , $|1\rangle=[0,1]^{T}$ , and the unnormalized vectors are

\bm{x}=[\widehat{\bm{w}}^{1};\cdots;\widehat{\bm{w}}^{N_{t}}],\qquad\bm{y}=[\widehat{\bm{w}}^{N_{t}},\cdots,\widehat{\bm{w}}^{N_{t}}],

satisfying

\|\widehat{\bm{w}}^{N_{t}}\|^{2}=\frac{1}{2N_{t}}=\frac{1}{2N_{t}N_{\psi_{0}}^{2}}\|\bm{w}^{N_{t}}\|^{2}.

Let us block the matrix $O$ as $(O_{ij})$ according to the structure of $\bm{w}$ . One easily finds that $O_{ij}=\bm{O}$ are zeros matrices when $(i,j)\neq(N_{t},N_{t})$ . Then,

$\displaystyle\langle O\rangle$	$\displaystyle=\langle\psi\|O\|\psi\rangle=\frac{1}{N_{\psi}^{2}}\langle\bm{w}\|O\|\bm{w}\rangle=\frac{1}{N_{\psi}^{2}}(\bm{w}^{N_{t}})^{\dagger}O_{N_{t},N_{t}}\bm{w}^{N_{t}}$	(3.30)
	$\displaystyle=\frac{2N_{t}N_{\psi_{0}}^{2}}{N_{\psi}^{2}}(\widehat{\bm{w}}^{N_{t}})^{\dagger}O_{N_{t},N_{t}}\widehat{\bm{w}}^{N_{t}}=\frac{2N_{\psi_{0}}^{2}}{N_{\psi}^{2}}(\bm{y})^{\dagger}O_{\bm{y}}\bm{y}^{N_{t}}$
	$\displaystyle=\frac{N_{\psi_{0}}^{2}}{M^{d/2}}\frac{M^{d/2}}{N_{\psi}^{2}}\langle\widehat{\bm{w}}\|\widehat{O}\|\widehat{\bm{w}}\rangle=n_{\psi_{0}}^{2}\frac{1}{n_{\psi}^{2}}\langle\widehat{\bm{w}}\|\widehat{O}\|\widehat{\bm{w}}\rangle,$

where $O_{\bm{y}}=\text{diag}(O_{N_{t},N_{t}},\cdots,O_{N_{t},N_{t}})$ , and $\widehat{O}=\text{diag}(\bm{O},\cdots,\bm{O},O_{\bm{y}})$ . It is evident that $\text{Var}(\widehat{O})\sim\text{Var}(O)$ , and hence $\text{Var}(\widehat{O}/n_{\psi}^{2})\lesssim 1$ , which implies that the new multiplicative factor is $n_{\psi_{0}}^{4}$ , as expected. It’s worth pointing out that the solution vector in (3.29) only requires one ancilla qubit.

Boosting the success probability of the final state projection

It can be seen from (3.30) that the observable is an expectation value taken with respect to the value of the state at the final time. When the solution decays exponentially in time, the success probability of projecting the history state onto the final state is exponentially small. One can raise the success probability via the amplitude amplification as in [7, 18]. This implies a multiplicative factor $g=\max_{t\in[0,T]}\|\bm{w}(t)\|/\|\bm{w}(T)\|$ in the time complexity, which characterises the decay of the final state relative to the initial state.

In conclusion, the additional multiplicative factor in the time complexity is $gn_{\psi_{0}}^{4}/\varepsilon^{2}$ for the computation of the observable. In the subsequent discussion, we will ignore the parameter $g$ for convenience. For the range of $n_{\psi_{0}}$ , please refer to Lemma 14 in [38] for some discussions.

3.2.3 The gate complexity of computing the observable

The corresponding result for the quantum state is given in Theorem 3.1. We now consider the computation of the observable. In this case, we have to discuss the additional cost in estimating the norm of $\bm{w}$ . Let $\langle\tilde{O}\rangle_{\text{app}}$ be the approximate value of $\langle\tilde{O}\rangle$ . Then the total error Err satisfies

\text{Err}:=|\langle O\rangle-\langle\tilde{O}\rangle_{\text{app}}|\leq|\langle O\rangle-\langle\tilde{O}\rangle|+|\langle\tilde{O}\rangle-\langle\tilde{O}\rangle_{\text{app}}|=:I_{1}+I_{2},

(3.31)

where the first term is for the estimation of $N_{\psi}=\|\bm{w}\|$ , and the second one is for the sampling.

We first need to evaluate $N_{\psi}=\|\bm{w}\|$ , with the result described as follows. One can refer to [55, Theorem 16] and [13, Corollary 32] for details.

Lemma 3.3 (Estimation of $\|A^{-1}b\|$ ).

Let $Ax=b$ for an $N\times N$ matrix with sparsity $s$ and condition number $\kappa$ . Then there exists a quantum algorithm that outputs $\tilde{\alpha}$ such that

|\tilde{\alpha}-\|x\||\leq\eta\|x\|

with probability at least 0.99, in time

\mathcal{O}\Big{(}(T_{U}+T_{b})\frac{\kappa}{\eta}\log^{3}\kappa\log\log\frac{\kappa}{\eta}\Big{)},

where $T_{b}$ is the time of constructing the state $|b\rangle=\frac{1}{\|b\|}\sum b_{i}|i\rangle$ , and

T_{U}=\log N\Big{(}\log N+\log^{2.5}\frac{s\kappa\log(\kappa/\eta)}{\eta}\Big{)}\log^{2}\frac{\kappa}{\eta}.

The cost $T_{b}$ is neglected throughout the paper. With the help of the above result, we are able to bound the gate complexity of computing the observable for the QLSA of the Liouville representation.

Theorem 3.2.

Suppose the condition of Theorem 3.1 is satisfied and $\sup_{x}|F_{i}(x)|=\mathcal{O}(1)$ for $i=1,\cdots,d$ . Given the error tolerance $\varepsilon$ , if the QLSA for the upwind finite difference discretisation is used, then the observable of the Liouville representation (3.5) can be computed with gate complexity given by

N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{L}^{4}d^{3}}{\varepsilon^{5}}\log\frac{1}{\varepsilon}\Big{)},

where $n_{L}=\|(\bm{\rho}^{\omega})^{0}\|/M^{d/2}$ .

Proof.

(1) For the error $I_{1}$ in (3.31), let $\alpha$ and $\tilde{\alpha}$ be the exact and approximate norms of $\bm{w}$ , respectively. Denote $\tilde{n}_{\psi}=\tilde{N}_{\psi}/M^{d/2}=\tilde{\alpha}/M^{d/2}$ and $n_{\psi}=N_{\psi}/M^{d/2}=\alpha/M^{d/2}$ , where $|\tilde{\alpha}-\alpha|\leq\eta\alpha$ . Then

|(\tilde{n}_{\psi})^{2}-(n_{\psi})^{2}|\leq\eta n_{\psi}(n_{\psi}+\tilde{n}_{\psi})\leq\eta(2+\eta)(n_{\psi})^{2},

and the error

\displaystyle I_{1}=|(n_{G}\tilde{n}_{\psi})^{2}\Upsilon-(n_{G}n_{\psi})^{2}\Upsilon|\leq\eta(1+\eta)(n_{G}n_{\psi})^{2}\Upsilon=\eta(2+\eta)\langle O\rangle.

This suggests to take $\eta=\mathcal{O}(\varepsilon)$ since $\langle O\rangle=\mathcal{O}(1)$ .

For the error $I_{2}$ , by the general sampling law, we can obtain an approximation $\langle\tilde{O}\rangle_{\text{app}}$ to precision $\varepsilon$ by repeating the quantum algorithm $k=\mathcal{O}(n_{\psi_{0}}^{4}/\varepsilon^{2})$ times, where $n_{\psi_{0}}=n_{L}:=\|(\bm{\rho}^{\omega})^{0}\|/M^{d/2}$ .

(2) According to Theorem 3.1, one has

M=1/\Delta x\sim d/\varepsilon^{3},\qquad\Delta t\sim\Delta x/d,\qquad\kappa\sim 1/\Delta t\sim d^{2}/\varepsilon^{3},\qquad s\sim d.

(3) Let $T_{1}$ be the gate complexity of obtaining the estimation of $N_{\psi}$ . By Lemma 3.3,

T_{1}=\widetilde{\mathcal{O}}(\kappa/\eta)=\widetilde{\mathcal{O}}(d^{2}/\varepsilon^{4}).

Let $T_{2}$ be the gate complexity of the QLSA. Then,

T_{2}=\widetilde{\mathcal{O}}\Big{(}s\kappa\log\frac{1}{\varepsilon}\Big{)}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3}}{\varepsilon^{3}}\log\frac{1}{\varepsilon}\Big{)}.

The overall gate complexity is

T=T_{1}+kT_{2}=\widetilde{\mathcal{O}}\Big{(}\frac{n_{\psi_{0}}^{4}d^{3}}{\varepsilon^{5}}\log\frac{1}{\varepsilon}\Big{)}.

The proof is completed. ∎

Remark 3.4.

As observed in the proof, the overall complexity $T$ of the algorithm is dominated by the complexity $kT_{2}$ of sampling. This implies, when computing the observables, we just need to multiply the original gate complexity under an appropriate mesh strategy by the sampling factor $k=\mathcal{O}({\rm Var}(O)/\varepsilon^{2})$ . We further remark that the classical cost contains exponential terms in dimension like $d^{d}$ and $(1/\varepsilon)^{d}$ , which is absent in applications where $n_{L}$ does not grow so quickly.

Remark 3.5.

Despite the absence of the unitary structure, we can still propose a “quantum simulation” algorithm for the Liouville representation as shown in Appendix A by using the dimensional splitting Trotter based approximation. The basic idea of the algorithm is to transform the asymmetric evolution in each direction into a symmetric one, which requires only a simple variable substitution with the transformation matrix being diagonal. However, unlike the traditional time-marching Hamiltonian simulation, non-unitary procedures for the variable substitution are involved, which leads to exponential increase of the cost arising from multiple copies of initial quantum states at every time step as pointed out in Remark A.1.

3.3 Finite difference discretisation for the KvN representation

3.3.1 The QLSA for the finite difference discretisation

We consider the upwind finite difference discretisation for the KvN representation, which can be written as

\begin{cases}\partial_{t}u+\sum\limits_{i=1}^{d}F_{i}\frac{\partial u}{\partial x_{i}}+\frac{1}{2}({\rm div}F)u=0,\qquad u=\psi^{\omega},\\ u(x,0)=\psi_{0}^{\omega}.\end{cases}

(3.32)

The scheme reads

\displaystyle u_{\bm{j}}^{n+1}-\Big{[}1-\frac{1}{2}\Delta t({\rm div}F)_{\bm{j}}-\lambda\sum\limits_{\ell=1}^{d}(b_{\bm{j}}^{\ell,+}-b_{\bm{j}}^{\ell,-})\Big{]}u_{\bm{j}}^{n}+\lambda\sum\limits_{\ell=1}^{d}\Big{[}b_{\bm{j}}^{\ell,-}u_{\bm{j}+\bm{e}_{\ell}}^{n}-b_{\bm{j}}^{\ell,+}u_{\bm{j}-\bm{e}_{\ell}}^{n}\Big{]}=0,

(3.33)

where

b_{\bm{j}}^{k,\pm}=\Big{\{}F_{k}\Big{\}}_{\bm{j}}^{\pm},\qquad\bm{j}=(j_{1},\cdots,j_{d}).

In matrix form one has

\bm{u}^{n+1}-B\bm{u}^{n}=0,\quad n=0,1,\cdots,N_{t}-1.

(3.34)

The final coefficient matrix $L$ is of the same form as in Eq. (3.25).

Theorem 3.3.

Suppose $\lambda=\Delta t/\Delta x$ satisfies the following CFL condition

\lambda\sum\limits_{i=1}^{d}\sup_{x}|F_{i}(x)|\leq 1.

(1)

The condition number and the sparsity of $L$ satisfy $\kappa\lesssim{1}/{\Delta t}$ and $s=\mathcal{O}(d)$ .
(2)

For fixed spatial step $\Delta x$ , let $\Delta t=\mathcal{O}(\Delta x/d)$ and $\omega=(d\Delta x)^{1/3}$ . Given the error tolerance $\varepsilon$ , the gate complexity of the QLSA is

$N_{\text{Gates}}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3}}{\varepsilon^{3}}\log\frac{1}{\varepsilon}\Big{)}.$

Proof.

The proof is similar to the argument in Theorem 3.1, so we omit the details. ∎

Remark 3.6.

For the upwind discretisation of the KvN representation, the forcing term $\frac{1}{2}({\rm div}F)u$ will contribute to the $l_{1}$ error an exponentially growing term like $\mathrm{e}^{t_{n}\|\text{div}F\|_{\infty}}$ . In fact, one easily obtains from (3.33) that

	$\displaystyle\\|\bm{u}^{n+1}\\|_{1}$	$\displaystyle\leq\Big{(}1+\frac{1}{2}\Delta t\\|\text{div}F\\|_{\infty}\Big{)}\\|\bm{u}^{n}\\|_{1}\leq\Big{(}1+\frac{1}{2}\Delta t\\|\text{div}F\\|_{\infty}\Big{)}^{n}\\|\bm{u}^{0}\\|_{1}$
		$\displaystyle=\Big{(}1+\frac{1}{2}\frac{t_{n}}{n}\\|\text{div}F\\|_{\infty}\Big{)}^{n}\\|\bm{u}^{0}\\|_{1}\leq\mathrm{e}^{t_{n}\\|\text{div}F\\|_{\infty}}\\|\bm{u}^{0}\\|_{1}.$

As a comparison, if one uses the Liouville equation, the $l_{1}$ norm of the error is contracting, as shown in Remark 3.3.

3.3.2 The gate complexity of computing the observable

Let $n=N_{t}$ for convenience. For the KvN approximation, from (3.1.3) we know that the observable at $t=t_{n}$ is

\langle O_{\psi^{\omega},n}\rangle=\frac{1}{M^{d/2}}(\bm{\psi}_{n}^{\omega})^{\dagger}G_{M}\bm{\psi}_{n}^{\omega}\approx\frac{1}{M^{d/2}}(\bm{u}^{n})^{\dagger}G_{M}\bm{u}^{n}=\frac{1}{M^{d/2}}\bm{u}^{\dagger}\mathcal{G}_{M}\bm{u},

where $\bm{u}$ is the solution of (3.34), and $\mathcal{G}_{M}=\text{diag}(\bm{O},\cdots,\bm{O},G_{M})$ . Let

|\psi\rangle=\frac{1}{N_{\psi}}\sum\limits_{\bm{j},n}\bm{u}_{\bm{j},n}|\bm{j}\rangle|n\rangle,

where the normalisation $N_{\psi}=\|\bm{u}\|$ . The expectation of the observable can be defined as

\langle O_{\psi^{\omega},n}\rangle=\frac{1}{M^{d/2}}\langle\psi|O|\psi\rangle=:\langle O\rangle,\qquad O=\frac{N_{\psi}^{2}}{M^{d/2}}\mathcal{G}_{M}.

One easily finds from (3.1.3) that $\text{Var}(\mathcal{G}_{M})$ is bounded since

\text{Var}(\mathcal{G}_{M})=\langle\mathcal{G}_{M}^{2}\rangle-\langle\mathcal{G}_{M}\rangle^{2}\leq\|\mathcal{G}_{M}\psi\|^{2}\lesssim 1.

Following the similar analysis in Subsect. 3.2.2 for the Liouville representation, one obtains the multiplicative factor in the time complexity can be given by $n_{K}^{4}/\varepsilon^{2}$ , where $n_{K}=\|\bm{u}^{0}\|/M^{d/4}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/4}$ , where we have omitted the parameter $g$ which characterises the decay of the final state relative to the initial state. The arguments for computing the observable of the Liouville representation also apply to the KvN representation. The corresponding result is described in the following theorem.

Theorem 3.4.

Suppose the condition of Theorem 3.3 is satisfied and $\sup_{x}|F_{i}(x)|=\mathcal{O}(1)$ for $i=1,\cdots,d$ . Given the error tolerance $\varepsilon$ , if the QLSA for the upwind finite difference discretisation is used, then the observable of the KvN representation can be computed with gate complexity given by

N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{K}^{4}d^{3}}{\varepsilon^{5}}\log\frac{1}{\varepsilon}\Big{)},

where $n_{K}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/4}$ .

Proof.

The KvN representation has the same error estimate as the Liouville representation, which leads to the same mesh strategy. From Theorem 3.3, we also observe the same condition number and sparsity for the associated coefficient matrix. We therefore obtain the same gate complexity with the multiplicative factor replaced by $n_{K}^{4}/\varepsilon^{2}$ . ∎

Remark 3.7.

In view of the relation (3.12), one easily finds that

n_{L}=\frac{\|(\bm{\rho}^{\omega})^{0}\|}{M^{d/2}}=\frac{\|(\bm{\psi}^{\omega})^{0}\|^{2}}{M^{d/2}}=n_{K}^{2}.

3.4 Spectral discretisation for the KvN representation

The KvN representation can be solved by quantum Hamiltonian simulation directly since the evolutionary operator is Hermitian, where the Hamiltonian simulation can be realised by the quantum version of the classical Fourier spectral method. On the other hand, we can also develop the QLSA based method for the spectral discretisation.

3.4.1 The notations

We consider the Fourier spectral discretisation. To this end, we first introduce some notations frequently used in this article.

For one-dimensional problems we choose a uniform spatial mesh size $\Delta x=1/M$ for $M=2N=2^{m}$ with $m$ an positive integer and the time step $\Delta t$ , and we let the grid points and the time step be

x_{j}=j\Delta x,~{}~{}t_{n}=n\Delta t,\quad j=0,1,\cdots,N,~{}~{}n=0,1,\cdots.

We consider the periodic boundary conditions. For $x\in[0,1]$ , the 1-D basis functions for the Fourier spectral method are usually chosen as

\phi_{l}(x)=\mathrm{e}^{\mathrm{i}\mu_{l}x},\quad\mu_{l}=2\pi l,\quad 1=-N,\cdots,N-1.

For convenience, we adjust the index as

\phi_{l}(x)=\mathrm{e}^{\mathrm{i}\mu_{l}x},\quad\mu_{l}=2\pi(l-N-1),\quad 1\leq l\leq M=2N.

The approximation in the 1-D space is

u(t,x)=\sum\limits_{l=1}^{M}c_{l}(t)\phi_{l}(x),\quad x=x_{j},~{}~{}j=0,1,\cdots,M-1.

(3.35)

which can be written in vector form, $\bm{u}(t)=\Phi\bm{c}(t)$ , where

\bm{u}(t)=(u(t,x_{j}))_{M\times 1},\quad\bm{c}=(c_{l})_{M\times 1},\quad\Phi=(\phi_{jl})_{M\times M}=(\phi_{l}(x_{j}))_{M\times M}.

The $d$ -dimensional grid points are then given by ${x}_{\bm{j}}=(x_{j_{1}},\cdots,x_{j_{d}})$ , where $\bm{j}=(j_{1},\cdots,j_{d})$ , and

x_{j_{i}}=j_{i}\Delta x,\quad j_{i}=0,1,\cdots,M-1,\quad i=1,\cdots,d.

We use the notation $1\leq\bm{j}\leq M$ to indicate $1\leq j_{i}\leq M$ for every component of $\bm{j}$ . The multi-dimensional basis functions are written as $\phi_{\bm{l}}({x})=\phi_{l_{1}}(x_{1})\cdots\phi_{l_{d}}(x_{d})$ , where $\bm{l}=(l_{1},\cdots,l_{d})$ and $1\leq\bm{l}\leq M$ . The corresponding approximate solution is $u(t,{x})=\sum\nolimits_{\bm{l}}c_{\bm{l}}(t)\phi_{\bm{l}}({x})$ , with the coefficients determined by the exact values at the grid or collocation points ${x}_{\bm{j}}$ . These collocation values will be arranged as a column vector:

\bm{u}(t)=\sum\limits_{\bm{j}}u(t,{x}_{\bm{j}})|j_{1}\rangle\otimes\cdots\otimes|j_{d}\rangle.

That is, the $n_{\bm{j}}$ -th entry of $\bm{u}$ is $u(t,{x}_{\bm{j}})$ , with the global index given by

n_{\bm{j}}:=j_{1}2^{d-1}+\cdots+j_{d}2^{0},\qquad\bm{j}=(j_{1},\cdots,j_{d}).

Similarly $c_{\bm{l}}$ is written in a column vector as $\bm{c}=\sum\nolimits_{\bm{l}}c_{\bm{l}}|l_{1}\rangle\otimes\cdots\otimes|l_{d}\rangle$ .

To determine the transformation matrix between $\bm{u}$ and $\bm{c}$ , let $c_{\bm{l}}=c_{l_{1}}\cdots c_{l_{d}}$ . Then

u(t,{x}_{\bm{j}})=\sum\limits_{\bm{l}}c_{l_{1}}\cdots c_{l_{d}}\phi_{l_{1}}(x_{j_{1}})\cdots\phi_{l_{d}}(x_{j_{d}}).

(3.36)

The direct calculation gives

\displaystyle\sum\limits_{\bm{j}}u(t,{x}_{\bm{j}})|j_{1}\rangle\otimes\cdots\otimes|j_{d}\rangle=\Big{(}\sum\limits_{j_{1},l_{1}}c_{l_{1}}\phi_{l_{1}}(x_{j_{1}})|j_{1}\rangle\Big{)}\otimes\cdots\otimes\Big{(}\sum\limits_{j_{d},l_{d}}c_{l_{d}}\phi_{l_{d}}(x_{j_{1}})|j_{d}\rangle\Big{)},

which implies

\displaystyle\bm{u}=(\Phi\bm{c}^{(1)})\otimes\cdots\otimes(\Phi\bm{c}^{(d)})=(\Phi\otimes\cdots\otimes\Phi)(\bm{c}^{(1)}\otimes\cdots\otimes\bm{c}^{(d)})=\Phi^{\otimes^{d}}\bm{c},

where

\Phi^{\otimes^{d}}=\underbrace{\Phi\otimes\cdots\otimes\Phi}_{\text{$d$ matrices}},\qquad\bm{c}^{(i)}=(c_{l_{i}})_{M\times 1},

\bm{c}=\bm{c}^{(1)}\otimes\cdots\otimes\bm{c}^{(d)}=\sum\limits_{\bm{l}}c_{\bm{l}}|l_{1}\rangle\otimes\cdots\otimes|l_{d}\rangle.

(3.37)

This shows that by arranging ${x}_{\bm{j}}$ in the order of $|j_{1}\rangle\otimes\cdots\otimes|j_{d}\rangle$ , and ${c}_{\bm{l}}$ in the order of $|l_{1}\rangle\otimes\cdots\otimes|l_{d}\rangle$ , the corresponding coefficient matrix is exactly the tensor product of the matrices in one dimension.

For later use, we next determine the transitions between the position operator $\hat{x}_{j}$ and the momentum operator $\hat{P}_{j}=-\mathrm{i}\frac{\partial}{\partial x_{j}}$ in discrete settings.

We first consider the one-dimensional case. Let $u(x)$ be a function in one dimension and $\bm{u}=[u(x_{0}),\cdots,u(x_{M-1})]^{T}$ be the mesh function with $M=2N$ . The discrete position operator $\hat{x}^{\rm d}$ of $\hat{x}$ can be defined as

\hat{x}^{\rm d}:\bm{u}=\Big{(}u(x_{i})\Big{)}\quad\to\quad\Big{(}x_{i}u(x_{i})\Big{)}=D_{x}\bm{u}\qquad\mbox{or}\qquad\hat{x}^{\rm d}\bm{u}=D_{x}\bm{u},

where $D_{x}=\text{diag}(x_{0},x_{1},\cdots,x_{M-1})$ is the matrix representation of the position operator in $x$ -space. By the discrete Fourier expansion in (3.35), the momentum operator can be discretised as

	$\displaystyle\hat{P}u(x)$	$\displaystyle\approx\hat{P}\sum\limits_{l=1}^{M}c_{l}\phi_{l}(x)=\sum\limits_{l=1}^{M}c_{l}\hat{P}\phi_{l}(x)=\sum\limits_{l=1}^{M}c_{l}(-\mathrm{i}\partial_{x}\phi_{l}(x))$
		$\displaystyle=\sum\limits_{l=1}^{M}c_{l}\mu_{l}\phi_{l}(x),\quad\mu_{l}=2\pi(l-N-1)$

for $x=x_{j}$ , $j=0,1,\cdots,M-1$ , which is written in matrix form as

\hat{P}^{\rm d}\bm{u}=\Phi D_{\mu}\Phi^{-1}\bm{u}=:P_{x}\bm{u},\qquad D_{\mu}=\text{diag}(\mu_{1},\cdots,\mu_{M}),

where $\hat{P}^{\rm d}$ is the discrete momentum operator. The matrices $D_{\mu}$ and $P_{x}$ can be referred to as the matrix representation of the momentum operator in $p$ -space and $x$ -space, respectively, and are related by the discrete Fourier transform.

For $d$ dimensions, we still denote $\bm{u}=\sum\nolimits_{\bm{j}}u(x_{\bm{j}})|j_{1}\rangle\cdots|j_{d}\rangle$ . Let

u(x_{\bm{j}})=u(x_{j_{1}},\cdots,x_{j_{d}})=u^{(1)}(x_{j_{1}})\cdots u^{(d)}(x_{j_{d}}),

where $\bm{u}^{(l)}=\Phi\bm{c}^{(l)}$ . One has $\bm{u}=\bm{u}^{(1)}\otimes\cdots\otimes\bm{u}^{(d)}$ . The discrete position operator $\hat{x}_{l}^{\rm d}$ is defined as

\hat{x}_{l}^{\rm d}:\bm{u}=\bm{u}^{(1)}\otimes\cdots\otimes\bm{u}^{(d)}\quad\to\quad\bm{u}^{(1)}\otimes\cdots\otimes\tilde{\bm{u}}^{(l)}\otimes\cdots\otimes\bm{u}^{(d)},

where

\tilde{\bm{u}}^{(l)}:=\Big{(}x_{j_{l_{i}}}u^{(l)}(x_{j_{l_{i}}})\Big{)}=D_{x}\bm{u}^{(l)}.

Then,

\hat{x}_{l}^{\rm d}\bm{u}=(I^{\otimes^{l-1}}\otimes D_{x}\otimes I^{\otimes^{d-l}})\bm{u}=:\bm{D}_{l}\bm{u}.

Using the expansion in (3.36), one easily finds that

\hat{P}_{l}^{\rm d}\bm{u}=(I^{\otimes^{l-1}}\otimes P_{x}\otimes I^{\otimes^{d-l}})\bm{u}=:\bm{P}_{l}\bm{u}.

Note that

(\Phi^{\otimes^{d}})^{-1}\bm{P}_{l}\Phi^{\otimes^{d}}=I^{\otimes^{l-1}}\otimes D_{\mu}\otimes I^{\otimes^{d-l}}=:\bm{D}^{\mu}_{l}.

(3.38)

3.4.2 The QLSA for the spectral discretisation

Let us consider the matrix representation of the operator $\hat{H}_{j}$ in (3.8), where

\hat{H}_{j}=\frac{1}{2}(F_{j}(\hat{x})\hat{P}_{j}+\hat{P}_{j}F_{j}(\hat{x})).

For clarity, we still use $\bm{u}$ to denote the mesh function of $\psi$ (see the notations in Subsect. 3.4.1). When performing series expansion on $F$ , one has

	$\displaystyle F_{j}(\hat{x}^{\rm d})\bm{u}$	$\displaystyle:=\sum\limits_{\bm{l}}a_{\bm{l}}(\hat{x}_{1}^{\rm d})^{l_{1}}\cdots(\hat{x}_{d}^{\rm d})^{l_{d}}\bm{u}=\sum\limits_{\bm{l}}a_{\bm{l}}(\bm{D}_{1}^{l_{1}}\cdots\bm{D}_{d}^{l_{d}})\bm{u}$
		$\displaystyle=\sum\limits_{\bm{l}}a_{\bm{l}}(D_{x}^{l_{1}}\otimes\cdots\otimes D_{x}^{l_{d}})\bm{u}=:\bm{F}_{j}\bm{u}$

in the discrete setting, where $\bm{F}_{j}$ is clearly a diagonal matrix. We assume that the series expansion is accurate enough to simplify the discussion. Then one has

\displaystyle\hat{H}_{j}^{\rm d}\bm{u}

\displaystyle=\frac{1}{2}(F_{j}(\hat{x}^{\rm d})\hat{P}_{j}^{\rm d}+\hat{P}_{j}^{\rm d}F_{j}(\hat{x}^{\rm d}))\bm{u}=\frac{1}{2}(\bm{F}_{j}\bm{P}_{j}+\bm{P}_{j}\bm{F}_{j})\bm{u}=:\bm{H}_{j}\bm{u}.

One easily finds that the sparsity of $\bm{H}_{j}$ is $\mathcal{O}(M)$ . The resulting system of ordinary differential equations is

\begin{cases}\frac{\mathrm{d}}{\mathrm{d}t}\bm{u}(t)=A\bm{u}(t),\qquad A=-\mathrm{i}\sum\limits_{j=1}^{d}\bm{H}_{j},\\ \bm{u}(0)=(\psi^{\omega}(0,x_{\bm{j}})).\end{cases}

(3.39)

The analytic solution is obviously given by

\bm{u}(t)=\mathrm{e}^{At}\bm{u}(0)=\exp\Big{(}-\mathrm{i}\sum\limits_{j=1}^{d}\bm{H}_{j}\Big{)}\bm{u}(0),

which implies $\|\bm{u}(t)\|=\|\bm{u}(0)\|$ for any $t\geq 0$ since $\bm{H}_{j}$ are real symmetric matrices. Let $n=N_{t}$ and denote

|\psi\rangle=\frac{1}{N_{\psi}}\sum\limits_{\bm{j}}\bm{u}_{\bm{j}}^{n}|\bm{j}\rangle,\qquad N_{\psi}=\|\bm{u}^{n}\|.

Since $N_{\psi}=\|\bm{u}^{n}\|=\|\bm{u}^{0}\|=:N_{\psi_{0}}$ , the observable can be reformulated as

\langle O_{\psi^{\omega},n}\rangle=\frac{1}{M^{d/2}}(\bm{\psi}_{n}^{\omega})^{\dagger}G_{M}\bm{\psi}_{n}^{\omega}\approx\frac{1}{M^{d/2}}(\bm{u}^{n})^{\dagger}G_{M}\bm{u}^{n}=\frac{N_{\psi_{0}}^{2}}{M^{d/2}}\langle\psi|G_{M}|\psi\rangle=:\langle\psi|O|\psi\rangle.

The ODEs in (3.39) can be solved by the quantum differential equations solver reported in [5, 7, 18]. Here we consider the one in [7] with the result described below. For convenience, we still refer it to as the QLSA based method since the approach in [7] applies the QLSA.

Lemma 3.4.

Suppose $A=V^{-1}DV$ is an $N\times N$ diagonalizable matrix, where $D=\text{diag}(\lambda_{1},\cdots,\lambda_{N})$ satisfies $\text{Re}(\lambda_{j})\leq 0$ for any $j\in\{1,\cdots,N\}$ . In addition, suppose $A$ has at most s nonzero entries in any row and column, and we have an oracle $O_{A}$ that computes these entries. Suppose $x_{in}$ and $b$ are $N$ -dimensional vectors with known norms and that we have two controlled oracles, $O_{x}$ and $O_{b}$ , that prepare the states proportional to $x_{in}$ and $b$ , respectively. Let $x$ evolve according to the differential equation

\frac{\mathrm{d}x}{\mathrm{d}t}=Ax+b

with the initial condition $x(0)=x_{in}$ . Let $T>0$ and $g=\max_{t\in[0,T]}\|x(t)\|/\|x(T)\|$ . Then there exists a quantum algorithm that produces a state $\epsilon$ -close to $x(T)/\|x(T)\|$ in $l^{2}$ norm, succeeding with probability $\Omega(1)$ , with a flag indicating success, using

\mathcal{O}\Big{(}s\kappa_{V}\|A\|gT\cdot\text{Poly}(\log(s\kappa_{V}\|A\|gT\beta/\varepsilon))\Big{)}

queries to $O_{A}$ , $O_{x}$ , and $O_{b}$ , where $\kappa_{V}$ is the condition number of the transformation matrix $V$ , $g$ characterises the decay of the final state relative to the initial state, and $\beta=(\|x_{in}\|+T\|b\|)/\|x(T)\|$ . The gate complexity of this algorithm is larger than its query complexity by a factor of $\text{Poly}(\log(sN\kappa_{V}\|A\|gT\beta/\varepsilon))$ .

Note that the parameter $g$ can be dropped if we only output the quantum state $|x\rangle$ , not the projection $|x(T)\rangle$ . For the approximate evolutionary operator in (A.2), we are ready to quantify the gate complexity for the QLSA.

Remark 3.8.

The authors in [7] utilised the matrix exponential to construct a linear system for the ODEs and solved the linear system by using the QLSA proposed in [20]. As claimed in [20], the gate complexity exceeds the query complexity by a multiplicative factor $\mathcal{O}(\log N+\log^{2.5}(s\kappa/\varepsilon))$ , where $N=\mathcal{O}(M^{d})$ is the order of the matrix $A$ . This implies the linear dependence of the dimension $d$ when considering the gate complexity with respect to the matrix order.

Theorem 3.5.

Assume that $\max_{1\leq j\leq d}\|\bm{F}_{j}\|=\mathcal{O}(1)$ and $T=\mathcal{O}(1)$ .

(1)

There exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{u}(T)/\|\bm{u}(T)\|$ with the gate complexity given by

$N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)}.$
(2)

The observable of the KvN representation can be computed with gate complexity given by

$N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{K}^{4}d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)},$

where $n_{K}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/4}$ .

Proof.

(1) Let $A=V^{-1}DV$ with $D=\text{diag}(\lambda_{1},\cdots,\lambda_{N})$ , where $N=\mathcal{O}(M^{d})$ is the order of $A$ . Since $\bm{H}_{j}$ are real symmetric matrices, the matrix $\sum_{j=1}^{d}\bm{H}_{j}$ has only real eigenvalues and the transformation matrix $V$ can be chosen as an orthogonal matrix. This implies $\text{Re}(\lambda_{j})=0$ for any $j\in\{1,\cdots,N\}$ and $\kappa_{V}=1$ . According to Lemma 3.4, there exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{u}(T)/\|\bm{u}(T)\|$ with the gate complexity given by

N_{Gates}=\mathcal{O}\Big{(}s\kappa_{V}\|A\|T\cdot\text{Poly}(\log(sN\kappa_{V}\|A\|T/\varepsilon))\Big{)},

where we have omitted the parameter $g$ that characterises the decay of the final state relative to the initial state.

It is evident that the sparsity of $A$ is $\mathcal{O}(dM)=\mathcal{O}(d/\Delta x)$ . The norm of $A$ satisfies

	$\displaystyle\\|A\\|$	$\displaystyle\leq\sum\limits_{j=1}^{d}\\|\bm{H}_{j}\\|\leq\sum\limits_{j=1}^{d}\\|\bm{F}_{j}\\|\\|\bm{P}_{j}\\|\leq\sum\limits_{j=1}^{d}\\|\bm{F}_{j}\\|\\|\bm{D}_{j}^{\mu}\\|$
		$\displaystyle\leq M\sum\limits_{j=1}^{d}\\|\bm{F}_{j}\\|\leq dM\max_{1\leq j\leq d}\\|\bm{F}_{j}\\|=d/\Delta x\cdot\max_{1\leq j\leq d}\\|\bm{F}_{j}\\|.$

According to Lemma 3.2, the error of the spectral discretisation is $\mathcal{O}(\omega+\Delta t^{\alpha}/\omega^{\alpha}+d\Delta x^{\ell}/\omega^{\ell+1})$ , where $\alpha$ is for the precision of the temporal discretisation which has been considered in the quantum algorithm in [7]. To reach a precision of $\varepsilon$ , one just needs to set $\omega\sim d\Delta x^{\ell}/\omega^{\ell+1}\sim\varepsilon$ , and gets $\Delta x\sim\varepsilon^{1+2/\ell}/d^{1/\ell}$ . Therefore, we have

N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2}}{\Delta x^{2}}\Big{)}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)}.

where in the last equal sign we have included the additional factor $d$ arising from the matrix order (see Remark 3.8).

(2) For the spectral discretisation, the constant $N_{\psi}$ is known. The desired estimate follows from the general sampling law. ∎

3.4.3 The quantum simulation for the spectral discretisation

The ODEs (3.39) can also be solved by quantum Hamiltonian simulations.

Theorem 3.6.

Given the error tolerance $\varepsilon$ , assume that $\max_{1\leq j\leq d}\|\bm{F}_{j}\|_{\infty}=\mathcal{O}(1)$ and the simulation time $t=\mathcal{O}(1)$ .

(1)

The semi-discrete problem (3.39) obtained from the spectral discretisation of the KvN representation can be simulated with gate complexity given by

$N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{3+4/\ell}}\Big{)}.$
(2)

The observable of the KvN representation can be computed with gate complexity given by

$N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{K}^{4}d^{3+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)},$

where $n_{K}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/4}$ .

Proof.

(1) Let $H=\sum_{j=1}^{d}\bm{H}_{j}$ . Then the evolution of (3.39) can be written as $|\psi(t)\rangle=\mathrm{e}^{-\mathrm{i}Ht}|\psi(0)\rangle$ . According to Theorem 1 in [6], $\mathrm{e}^{-\mathrm{i}Ht}$ can be simulated within error $\varepsilon$ with

\mathcal{O}\Big{(}\tau(m_{d}+\log^{2.5}(\tau/\varepsilon))\frac{\log(\tau/\varepsilon)}{\log\log(\tau/\varepsilon)}\Big{)}=\mathcal{O}(\tau m_{d}\cdot\text{polylog})

2-qubits gates, where $\tau=s\|H\|_{\max}t$ , $s$ is the sparsity of $H$ and $\|H\|_{\max}$ denotes the largest entry of $H$ in absolute value, and

\text{polylog}\equiv\log^{2.5}(\tau/\varepsilon)\frac{\log(\tau/\varepsilon)}{\log\log(\tau/\varepsilon)}.

This result is near-optimal by Theorem 2 therein.

The sparsity of $H$ is $s=\mathcal{O}(dM)$ . According to the proof of Theorem 3.5, the mesh strategy is $M=1/\Delta x=d^{1/\ell}/\varepsilon^{1+2/\ell}$ , and hence the number of qubits per dimension is

m=\mathcal{O}(\log M)=\mathcal{O}\Big{(}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\Big{)}.

The total number of qubits is $m_{d}=dm$ . With these settings, noting that

H_{\max}\leq\sum\limits_{j=1}^{d}\|\bm{H}_{j}\|_{\infty}\lesssim M\sum\limits_{j=1}^{d}\|\bm{F}_{j}\|_{\infty}\leq dM\cdot\max_{1\leq j\leq d}\|\bm{F}_{j}\|_{\infty},

one has

\tau=\mathcal{O}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)},\qquad\tau/\varepsilon=\mathcal{O}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{3+4/\ell}}\Big{)}.

The gate complexity for solving the ODEs is then given by

N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)}.

(2) The gate complexity for computing the observable is obtained from the general sampling law. ∎

One can also run the simulation along each direction by using the Trotter based approximation. The evolution of (3.8) can be written as

|\psi(t+\Delta t)\rangle=\mathrm{e}^{-\mathrm{i}(\hat{H}_{1}+\cdots+\hat{H}_{d})\Delta t}|\psi(t)\rangle.

Let

U_{\Delta t}=\mathrm{e}^{-\mathrm{i}\hat{H}_{d}\Delta t}\cdots\mathrm{e}^{-\mathrm{i}\hat{H}_{1}\Delta t}.

(3.40)

One has [62, 16]

\mathrm{e}^{-\mathrm{i}(\hat{H}_{1}+\cdots+\hat{H}_{d})\Delta t}=U_{\Delta t}+C_{H}\Delta t^{2},

(3.41)

where $C_{H}$ depends on the operator $\hat{H}=\hat{H}_{1}+\cdots+\hat{H}_{d}$ or the matrix $H$ , considered as $\mathcal{O}(1)$ in the following. Therefore, the problem is reduced to the simulation of each $\hat{H}_{j}$ .

Remark 3.9.

One can clearly make the time discretization second order by using Strang’s splitting. Since other methods use first order time discretization, in order to compare the time complexities on equal footing we also use first order time discretization, namely the simple splitting, here.

First, we determine the mesh strategy. According to the error estimate in Lemma 3.2 and noting Eq. (3.41), one has the error estimate

e_{\psi}\leq C(\omega+\Delta t/\omega+d\Delta x^{\ell}/\omega^{\ell+1}).

The above error bounds suggest the following mesh strategy:

M=1/\Delta x=\mathcal{O}(d^{1/\ell}/\varepsilon^{1+2/\ell}),\qquad\Delta t\sim\varepsilon^{2}.

(3.42)

Second, we quantify the number of gates used in the quantum simulation. According to Theorem 1 in [6], $\mathrm{e}^{-\mathrm{i}\hat{H}_{j}\Delta t}$ can be simulated within error $\eta$ with

\mathcal{O}\Big{(}\tau(m_{d}+\log^{2.5}(\tau/\eta))\frac{\log(\tau/\eta)}{\log\log(\tau/\eta)}\Big{)}

2-qubits gates, where $\tau=s\|\bm{H}_{j}\|_{\max}\Delta t$ , $s$ is the sparsity of $\bm{H}_{j}$ and $\|\bm{H}_{j}\|_{\max}$ denotes the largest entry of $\bm{H}_{j}$ in absolute value. One can check that the sparsity of $\bm{H}_{j}$ is $s=\mathcal{O}(M)$ . Therefore, $U_{\Delta t}$ defined in (3.40) can be simulated within error $\mathcal{O}(d\eta)$ [51, Proposition 1.12] with

N_{\text{Gates}}(U_{\Delta t})=\mathcal{O}(d\tau m_{d}\cdot\text{polylog}),

where

\tau=M\tilde{H}_{\max}\Delta t,\qquad\tilde{H}_{\max}=\max_{j}\|\bm{H}_{j}\|_{\max}\lesssim M\max_{j}\|\bm{F}_{j}\|_{\max},

\text{polylog}=\log^{2.5}(\tau/\eta)\frac{\log(\tau/\eta)}{\log\log(\tau/\eta)}.

We also need $d\eta=\mathcal{O}(\Delta t^{2})$ or $\eta=\mathcal{O}(\varepsilon^{4}/d)$ , and the number of qubits per dimension is $m=\mathcal{O}(\log(d^{1/\ell}/\varepsilon^{1+2/\ell}))$ . The total number of qubits is $m_{d}=dm$ . With these settings, we obtain

\tau=\mathcal{O}\Big{(}\Delta t\frac{d^{2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)},\qquad\tau/\eta=\mathcal{O}\Big{(}\frac{d^{1+2/\ell}}{\varepsilon^{6+4/\ell}}\Big{)},

and the total number of gates required to iterate to the $n$ -th step is

\displaystyle N_{\text{Gates}}=nN_{\text{Gates}}(U_{\Delta t})=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)}.

Compared with the one-step implementation in Theorem 3.6, the complexity is improved in sparsity.

Remark 3.10.

From time $t=t_{n}$ to time $t=t_{n+1}$ , we solve $\bm{u}^{n+1}=U_{\Delta t}^{\rm d}\bm{u}^{n}$ , where $U_{\Delta t}^{\rm d}$ is the discrete version of $U_{\Delta t}$ in (3.40), given by $U_{\Delta t}^{\rm d}=\mathrm{e}^{-\mathrm{i}\hat{H}_{d}^{\rm d}\Delta t}\cdots\mathrm{e}^{-\mathrm{i}\hat{H}_{1}^{\rm d}\Delta t}$ . According to the previous discussions, one has $U_{\Delta t}^{\rm d}\bm{u}^{n}=B\bm{u}^{n}$ , with $B=\mathrm{e}^{-\mathrm{i}\bm{H}_{d}\Delta t}\cdots\mathrm{e}^{-\mathrm{i}\bm{H}_{1}\Delta t}$ . One can alternatively solve a linear system $LU=F$ in the form of Eq. (3.25), which, however, is not a suitable algorithm because the sparsity of $L$ grows exponentially with the number of dimensions, i.e., $s(L)=\mathcal{O}(M^{d})$ .

4 Liouville representation for nonlinear Hamilton-Jacobi PDEs: finite difference vs. spectral approximations

In [38], the level set method was used to map the nonlinear Hamilton-Jacobi equation into linear Liouville equation in the phase space, based on which quantum algorithms were then constructed.

Hamilton-Jacobi equations take the following general form

		$\displaystyle\partial_{t}S+H(\nabla S,x)=0,$		(4.1)
		$\displaystyle S(0,x)=S_{0}(x)$

with $t\in\mathbb{R}^{+}$ , $x\in\mathbb{R}^{d}$ , $S(t,x)\in\mathbb{R}$ . Define $u=\nabla S\in\mathbb{R}^{d}$ . Then $u$ solves a hyperbolic system of conservation laws in gradient form:

		$\displaystyle\partial_{t}u+\nabla H(u,x)=0,$		(4.2)
		$\displaystyle u(0,x)=\nabla S_{0}(x).$

The level set function $\phi_{i}(t,x,p)$ can be defined by

\displaystyle\phi_{i}(t,x,p=u(t,x))=0,

where $i=1,\cdots,d$ and $\,x,p\in\mathbb{R}^{d}$ , and $u(t,x)$ is the solution of Eq. (4.2). The zero level set of $\phi$ is the set $\{(t,x,p)|\phi_{i}(t,x,p)=0\}$ . Since $u(t,x)$ solves Eq. (4.2), one can show that $\phi=(\phi_{1},\cdots,\phi_{d})\in\mathbb{R}^{d}$ solves a (linear!) Liouville equation [43]

\partial_{t}\phi+\nabla_{p}H\cdot\nabla_{x}\phi-\nabla_{x}H\cdot\nabla_{p}\phi=0.

(4.3)

The initial data can be chosen as

\phi_{i}(0,x,p)=p_{i}-u_{i}(0,x),\quad i=1,\cdots,d.

(4.4)

Then $u$ can be recovered from the intersection of the zero level sets of $\phi_{i}\,(i=1,\cdots,d)$ , namely

u(t,x)=\{p(t,x)|\,\phi_{i}(t,x,p)=0,\,i=1,\cdots,d\}.

To retrieve physical observables (and to avoid finding the zero level set of $\phi$ which is challenging) later, [38] proposed to solve for $\psi$ , defined by the following problem

		$\displaystyle\partial_{t}\psi+\nabla_{p}H\cdot\nabla_{x}\psi-\nabla_{x}H\cdot\nabla_{p}\psi=0,$		(4.5)
		$\displaystyle\psi(0,x,p)=\prod_{i=1}^{d}\delta(p_{i}-u_{i}(0,x)),$

whose analytical solution is $\psi(t,x,p)=\delta(\phi(t,x,p))$ . We have thus transformed a $(d+1)$ -dimensional nonlinear Hamilton-Jacobi PDE to a $(2d+1)$ -dimensional linear PDE – the Liouville equation, without any approximations or constraints on the nonlinearity. The mapping is exact, but at the expense of doubling the spatial dimension. In the following, we denote $w(t,x,p)$ to be the solution to (4.5) with smoothed initial data, i.e., $w=\psi^{\omega}$ .

4.1 Finite difference discretisation for the Liouville equation

4.1.1 The QLSA for the finite difference discretisation

Consider $x,p$ together as a new variable, and write $y=(x,p)=(x_{1},\cdots,x_{d},p_{1},\cdots,p_{d})=(y_{1},\cdots,y_{2d})$ . Use the same uniform mesh in each $y_{i}$ direction. Let $\bm{j}=(j_{1},\cdots,j_{d},j_{d+1},\cdots,j_{2d})$ . Then the upwind discretisations for each term of the equation in (4.5) are

	$\displaystyle\partial_{t}w$	$\displaystyle\longrightarrow\frac{w_{\bm{j}}^{n+1}-w_{\bm{j}}^{n}}{\Delta t},$
	$\displaystyle\frac{\partial H}{\partial p_{i}}\frac{\partial w}{\partial x_{i}}$	$\displaystyle\longrightarrow\frac{1}{\Delta x}\Big{\{}\frac{\partial H}{\partial p_{i}}\Big{\}}_{\bm{j}}^{-}(w^{n}_{\bm{j}+\bm{e}_{d+i}}-w^{n}_{\bm{j}})+\frac{1}{\Delta x}\Big{\{}\frac{\partial H}{\partial p_{i}}\Big{\}}_{\bm{j}}^{+}(w^{n}_{\bm{j}}-w^{n}_{\bm{j}-\bm{e}_{d+i}}),$
	$\displaystyle-\frac{\partial H}{\partial x_{k}}\frac{\partial w}{\partial p_{k}}$	$\displaystyle\longrightarrow-\frac{1}{\Delta x}\Big{\{}\frac{\partial H}{\partial x_{k}}\Big{\}}_{\bm{j}}^{+}(w^{n}_{\bm{j}+\bm{e}_{k}}-w^{n}_{\bm{j}})-\frac{1}{\Delta x}\Big{\{}\frac{\partial H}{\partial x_{k}}\Big{\}}_{\bm{j}}^{-}(w^{n}_{\bm{j}}-w^{n}_{\bm{j}-\bm{e}_{k}}),$

where

\alpha^{+}=\max\{\alpha,0\}=\frac{\alpha+|\alpha|}{2},\quad\alpha^{-}=\min\{\alpha,0\}=\frac{\alpha-|\alpha|}{2}.

For convenience we introduce the following notation

a_{\bm{j}}^{k,\pm}=\Big{\{}\frac{\partial H}{\partial x_{k}}\Big{\}}_{\bm{j}}^{\pm},\qquad b_{\bm{j}}^{i,\pm}=\Big{\{}\frac{\partial H}{\partial p_{i}}\Big{\}}_{\bm{j}}^{\pm}.

The discrete scheme can be written as

	$\displaystyle w_{\bm{j}}^{n+1}-w_{\bm{j}}^{n}$	$\displaystyle+\lambda\sum\limits_{i=1}^{d}\Big{[}b_{\bm{j}}^{i,-}(w^{n}_{\bm{j}+\bm{e}_{d+i}}-w^{n}_{\bm{j}})+b_{\bm{j}}^{i,+}(w^{n}_{\bm{j}}-w^{n}_{\bm{j}-\bm{e}_{d+i}})\Big{]}$
		$\displaystyle-\lambda\sum\limits_{k=1}^{d}\Big{[}a_{\bm{j}}^{k,+}(w^{n}_{\bm{j}+\bm{e}_{k}}-w^{n}_{\bm{j}})+a_{\bm{j}}^{k,-}(w^{n}_{\bm{j}}-w^{n}_{\bm{j}-\bm{e}_{k}})\Big{]}=0,$

	$\displaystyle w_{\bm{j}}^{n+1}$	$\displaystyle-\Big{[}1-\lambda\sum\limits_{\ell=1}^{d}(b_{\bm{j}}^{\ell,+}-b_{\bm{j}}^{\ell,-}+a_{\bm{j}}^{\ell,+}-a_{\bm{j}}^{\ell,-})\Big{]}w_{\bm{j}}^{n}$
		$\displaystyle+\lambda\sum\limits_{\ell=1}^{d}\Big{[}b_{\bm{j}}^{\ell,-}w^{n}_{\bm{j}+\bm{e}_{d+\ell}}-b_{\bm{j}}^{\ell,+}w^{n}_{\bm{j}-\bm{e}_{d+\ell}}-a_{\bm{j}}^{\ell,+}w^{n}_{\bm{j}+\bm{e}_{\ell}}+a_{\bm{j}}^{\ell,-}w^{n}_{\bm{j}-\bm{e}_{\ell}}\Big{]}=0.$		(4.6)

In matrix form one has

\bm{w}^{n+1}-B\bm{w}^{n}=\bm{f}^{n+1},\quad n=0,1,\cdots,N_{t}-1,

with $\bm{f}^{i}$ being the terms resulting from the initial and boundary conditions, where the nodal values at $t=t_{n}$ are arranged as

\bm{w}^{n}=\sum\limits_{\bm{j}}w_{\bm{j}}^{n}|j_{1}\rangle\otimes\cdots\otimes|j_{d}\rangle\otimes\cdots\otimes|j_{2d}\rangle.

That is, the $n_{\bm{j}}$ -th entry of $W^{n}$ is $w_{\bm{j}}^{n}$ , with the global index given by $n_{\bm{j}}:=j_{1}2^{2d-1}+\cdots+j_{2d}2^{0}$ . The non-zero entries of $B$ can be provided by using the global index as before. The resulting linear system is

L\bm{w}=F,

(4.7)

where

\bm{w}=[\bm{w}^{1};\cdots;\bm{w}^{N_{t}}],\qquad F=[\bm{f}^{1};\bm{f}^{2};\cdots;\bm{f}^{N_{t}}].

The coefficient matrix $L$ is of the same form as in Eq. (3.25). Note that we consider the Liouville equation with the smoothed initial data, and for periodic boundary conditions, one has $\bm{f}^{1}=B\bm{w}^{0}$ and $\bm{f}^{i}=\bm{0}$ for $i\geq 2$ .

Theorem 4.1.

Suppose $\lambda=\Delta t/\Delta x$ satisfies the following CFL condition

\lambda\sum\limits_{i=1}^{d}\Big{(}\sup_{x,p}|\partial_{x_{i}}H|+\sup_{x,p}|\partial_{p_{i}}H|\Big{)}\leq 1,

and assume

\sup_{x,p}|\partial_{x_{i}}H|+\sup_{x,p}|\partial_{p_{i}}H|=\mathcal{O}(1),\quad i=1,\cdots,d.

Then the condition number and the sparsity of $L$ satisfy $\kappa=\mathcal{O}(1/\Delta t)$ and $s=\mathcal{O}(d)$ . For fixed spatial step $\Delta x$ , let $\Delta t=\mathcal{O}(\Delta x/d)$ and $\omega=(d\Delta x)^{1/3}$ . Given the error tolerance $\varepsilon$ , the gate complexity of the quantum difference method is

N_{\text{Gates}}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3}}{\varepsilon^{3}}\log\frac{1}{\varepsilon}\Big{)}.

Proof.

The proof is similar to the argument in Theorem 3.1, so we omit the details. ∎

4.1.2 The computation of the physical observables

In the following, we consider the computation of the physical observables for the Liouville equation and assume the periodic boundary conditions. Then physical observables are defined as

\langle G(t,x)\rangle=\int_{\mathbb{R}^{d}}G(p)\psi(t,x,p)\mathrm{d}p.

(4.8)

For example, $G(p)=1,p,|p|^{2}/2$ yield density, momentum and kinetic energy respectively [38]. As in (3.1.3), one can compute the integral (4.8) by using the numerical quadrature rule

\langle G(t_{n},x_{\bm{j}})\rangle=\int_{\mathbb{R}^{d}}\psi(t_{n},x_{\bm{j}},p)\mathrm{d}p\approx\frac{1}{M^{d}}\sum_{\bm{l}}G_{\bm{l}}\psi^{\omega}_{\bm{j},\bm{l},n}=\frac{1}{M^{d}}\sum_{\bm{l}}G_{\bm{l}}w_{\bm{j},\bm{l},n}=:\langle G^{\omega}_{n,\bm{j}}\rangle,

(4.9)

where, $G_{\bm{l}}$ are the weights, $\bm{j}=(j_{1},\cdots,j_{d})$ and $\bm{l}=(l_{1},\cdots,l_{d})$ and $M$ is the number of points in each dimension of the $2d$ phase space.

Let $\bm{w}_{\bm{j},\bm{l},n}$ be the solution of the classical spectral method or the upwind finite difference method (for the smoothed initial data). Then for the QLSA one has

|\psi\rangle=\frac{1}{N_{\psi}}\sum\limits_{\bm{j},\bm{l},n}\bm{w}_{\bm{j},\bm{l},n}|\bm{j}\rangle|\bm{l}\rangle|n\rangle,

where the normalisation $N_{\psi}=\|\bm{w}\|$ . For the quantum simulation method, one can just remove the “time register” (in this case $N_{\psi}$ is for time $t=t_{n}$ ). With $G_{\bm{l}}$ in (4.9), we define the state

|G_{n,\bm{j}}\rangle:=\frac{1}{N_{G}}\sum_{\bm{l}}G^{\dagger}_{\bm{l}}|\bm{j}\rangle|\bm{l}\rangle|n\rangle

where $N_{G}=\sqrt{\sum_{\bm{l}}|G^{\dagger}_{\bm{l}}|^{2}}$ is the normalisation. Given the density matrix $\mathcal{G}:=|G_{n,\bm{j}}\rangle\langle G_{n,\bm{j}}|$ , we define $\Upsilon:=\langle\psi|\mathcal{G}|\psi\rangle$ . A simple algebra yields

\langle G(t_{n},x_{\bm{j}})\rangle\approx\langle G^{\omega}_{n,\bm{j}}\rangle=n_{\psi}n_{G}|\sqrt{\Upsilon}|,

where $n_{G}=N_{G}/M^{d/2}=\mathcal{O}(1)$ is known and $n_{\psi}=N_{\psi}/M^{d/2}$ may be unknown. We further define

\langle O\rangle=\langle G^{\omega}_{n,\bm{j}}\rangle^{2}=(n_{G}n_{\psi})^{2}\Upsilon:=\langle\psi|O|\psi\rangle,\qquad O=(n_{G}n_{\psi})^{2}\mathcal{G},

as in Subsect. 3.2.2.

Theorem 4.2.

Suppose the condition of Theorem 4.1 is satisfied. Given the error tolerance $\varepsilon$ , if the QLSA for the upwind finite difference discretisation is used, then the observable of the Liouville equation (4.5) can be computed with gate complexity given by

N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{H}^{4}d^{3}}{\varepsilon^{5}}\log\frac{1}{\varepsilon}\Big{)},

where $n_{H}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/2}$ .

Proof.

According to Remark 3.4, one just needs to multiply the gate complexity in Theorem 4.1 by the factor $n_{H}^{4}/\varepsilon^{2}$ . ∎

Remark 4.1.

According to Lemma 14 in [38], one has $n_{H}=\mathcal{O}(\beta^{d}M^{d/2})$ if we assume the initial data has support in a box of size $\beta$ .

4.2 Spectral discretisation for the Liouville equation

For the spectral discretisation, we only consider an important case, namely $H(x,p)=\frac{1}{2}|p|^{2}+V(x)$ , which is the total energy (kinetic and potential energy) of classical particle. The Liouville equation is then rewritten as

\begin{cases}\partial_{t}w+p\cdot\nabla_{x}w-\nabla_{x}V(x)\cdot\nabla_{p}w=0,\\ w(0,x,p)=\psi^{\omega}(0,x,p)=\prod_{i=1}^{d}\delta_{\omega}(p_{i}-u_{i}(0,x)),\end{cases}

(4.10)

where we have assumed the smoothed initial data.

4.2.1 The QLSA for the spectral discretisation

Now we consider solving the Liouville equation in (4.10) by using the Fourier spectral methods. For simplicity, the periodic boundary conditions are used for the spectral discretisation. To this end, we introduce some notations. We always assume that $x=(x_{1},\cdots,x_{d})\in[0,1]^{d}$ and $p=(p_{1},\cdots,p_{d})\in[0,1]^{d}$ . Introduce a new variable ${y}=({x},{p})=(x_{1},\cdots,x_{d},p_{1},\cdots,p_{d})$ , and set

w(t,{x},{p})=w(t,{y})=\sum\limits_{\bm{l}}c_{\bm{l}}(t)\phi_{\bm{l}}({y}),\quad\bm{l}=(l_{1},\cdots,l_{d},l_{d+1},\cdots,l_{2d}).

(4.11)

The collocation points are denoted by ${y}_{\bm{j}}$ with $\bm{j}=(\bm{j}_{x},\bm{j}_{p})$ . As in (3.37), we define $\bm{c}=\bm{c}_{x}\otimes\bm{c}_{p}$ , where $\bm{c}_{x}=\bm{c}^{(1)}\otimes\cdots\otimes\bm{c}^{(d)}$ and $\bm{c}_{p}=\bm{c}^{(d+1)}\otimes\cdots\otimes\bm{c}^{(2d)}$ . We also introduce the notation $\bm{w}=\bm{w}_{x}\otimes\bm{w}_{p}$ , where $\bm{w}_{x}=\bm{w}^{(1)}\otimes\cdots\otimes\bm{w}^{(d)}$ , $\bm{w}_{p}=\bm{w}^{(d+1)}\otimes\cdots\otimes\bm{w}^{(2d)}$ , and $\bm{w}^{(l)}=\Phi\bm{c}^{(l)}$ can be viewed as the approximate solution of $w$ in $y_{l}$ direction.

According to the discussion in Subsect. 3.4.1, the first term can be discretised as

	$\displaystyle p\cdot\nabla_{x}w$	$\displaystyle=\mathrm{i}\sum\limits_{l=1}^{d}y_{l+d}(-\mathrm{i}\partial_{y_{l}})w\longrightarrow\mathrm{i}\sum\limits_{l=1}^{d}\hat{y}_{l+d}\hat{P}_{l}^{\text{d}}(\bm{w}_{x}\otimes\bm{w}_{p})$
		$\displaystyle=\mathrm{i}\sum\limits_{l=1}^{d}(I^{\otimes^{d}}\otimes\bm{D}_{l})(\bm{P}_{l}\otimes I^{\otimes^{d}})(\bm{w}_{x}\otimes\bm{w}_{p})=\mathrm{i}\sum\limits_{l=1}^{d}(\bm{P}_{l}\otimes\bm{D}_{l})(\bm{w}_{x}\otimes\bm{w}_{p}).$

For the second term, one has

	$\displaystyle\nabla_{x}V({x})\cdot\nabla_{p}w$	$\displaystyle=(\partial_{1}V({x}),\cdots,\partial_{d}V({x}))\cdot\nabla_{p}w=:(v_{1}({x}),\cdots,v_{d}({x}))\cdot\nabla_{p}w$
		$\displaystyle=\sum\limits_{l=1}^{d}v_{l}(x)\partial_{y_{l+d}}w=\mathrm{i}\sum\limits_{l=1}^{d}v_{l}(x)(-\mathrm{i}\partial_{y_{l+d}})w\longrightarrow\mathrm{i}\sum\limits_{l=1}^{d}v_{l}(\hat{x}^{\text{d}})\hat{P}_{l+d}^{\text{d}}(\bm{w}_{x}\otimes\bm{w}_{p})$
		$\displaystyle=\mathrm{i}\sum\limits_{l=1}^{d}(\bm{V}_{l}\otimes I^{\otimes^{d}})(I^{\otimes^{d}}\otimes\bm{P}_{l})(\bm{w}_{x}\otimes\bm{w}_{p})=\mathrm{i}\sum\limits_{l=1}^{d}(\bm{V}_{l}\otimes\bm{P}_{l})(\bm{w}_{x}\otimes\bm{w}_{p}),$

where

\bm{V}_{k}=\text{diag}(\bm{v}_{k}),\quad\bm{v}_{k}=\sum\limits_{\bm{j}_{x}}v_{k}({x}_{\bm{j}_{x}})|j_{1}\rangle\cdots|j_{d}\rangle=\sum\limits_{\bm{j}_{x}}\partial_{k}V({x}_{\bm{j}_{x}})|j_{1}\rangle\cdots|j_{d}\rangle.

Let

A=-\mathrm{i}\sum\limits_{l=1}^{d}(\bm{P}_{l}\otimes\bm{D}_{l}-\bm{V}_{l}\otimes\bm{P}_{l})=:-\mathrm{i}\tilde{A}.

(4.12)

Note that $\tilde{A}$ is a real symmetric matrix. The resulting ODEs is

\begin{cases}\frac{\mathrm{d}}{\mathrm{d}t}\bm{u}(t)=A\bm{u}(t),\\ \bm{u}(0)=(\psi^{\omega}(0,y_{\bm{j}})).\end{cases}

(4.13)

We are ready to apply the quantum algorithm in [7] to solve the above ODEs, with the time complexity described below.

Theorem 4.3.

Assume that $\max_{1\leq l\leq d}\|\bm{V}_{l}\|=\mathcal{O}(1)$ and $T=\mathcal{O}(1)$ .

(1)

There exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{u}(T)/\|\bm{u}(T)\|$ with the gate complexity given by

$N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)}.$
(2)

The observable of the Liouville equation can be computed with gate complexity given by

$N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{H}^{4}d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)},$

where $n_{H}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/2}$ .

Proof.

The argument is similar to that of Theorem 3.5.

(1) Let $A=V^{-1}DV$ with $D=\text{diag}(\lambda_{1},\cdots,\lambda_{N})$ . Since the matrix $\tilde{A}$ in (4.12) is a real symmetric matrix, $\text{Re}(\lambda_{j})=0$ for any $j\in\{1,\cdots,N\}$ and $\kappa_{V}=1$ . According to Lemma 3.4, there exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{u}(T)/\|\bm{u}(T)\|$ with the gate complexity given by

N_{Gates}=\mathcal{O}\Big{(}s\kappa_{V}\|A\|T\cdot\text{Poly}(\log(sN\kappa_{V}\|A\|T/\varepsilon))\Big{)},

where we have omitted the parameter $g$ that characterises the decay of the final state relative to the initial state.

It is evident that the sparsity of $A$ is $\mathcal{O}(M)=\mathcal{O}(1/\Delta x)$ . The norm of $A$ satisfies

	$\displaystyle\\|A\\|$	$\displaystyle\leq\sum\limits_{l=1}^{d}\\|\bm{P}_{l}\otimes\bm{D}_{l}-\bm{V}_{l}\otimes\bm{P}_{l}\\|\leq\sum\limits_{l=1}^{d}(\\|\bm{P}_{l}\\|\\|\bm{D}_{l}\\|+\\|\bm{V}_{l}\\|\\|\bm{P}_{l}\\|)$
		$\displaystyle\lesssim M\sum\limits_{l=1}^{d}(\\|\bm{D}_{l}\\|+\\|\bm{V}_{l}\\|)\lesssim dM\max_{1\leq l\leq d}(1+\\|\bm{V}_{l}\\|).$

N_{Gates}=\mathcal{O}\Big{(}\frac{d}{\Delta x^{2}}\Big{)}=\mathcal{O}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)},

where in the last equal sign we have included the additional factor $d$ arising from the matrix order (see Remark 3.8).

(2) For the spectral discretisation, the constant $n_{\psi}=N_{\psi}/M^{d/2}$ is known since $N_{\psi}=\|\bm{u}^{N_{t}}\|=\|\bm{u}^{0}\|=N_{\psi_{0}}$ . The desired estimate follows from the general sampling law. ∎

4.2.2 The quantum simulation for the spectral discretisation

One can directly solve the ODEs (4.13) by using the quantum simulation since $\tilde{A}$ is real symmetric. We in the following consider the Fourier spectral methods based on the time-splitting approximations.

From time $t=t_{n}$ to time $t=t_{n+1}$ , the Liouville equation is solved in two steps: One solves

\partial_{t}w+p\cdot\nabla_{x}w=0

(4.14)

for one time step, followed by solving

\partial_{t}w-\nabla_{x}V({x})\cdot\nabla_{p}w=0

(4.15)

again for one time step.

Step 1. According to the previous discussion, one has

\displaystyle p\cdot\nabla_{x}w\longrightarrow\mathrm{i}\sum\limits_{l=1}^{d}(\bm{P}_{l}\otimes\bm{D}_{l})(\bm{w}_{x}\otimes\bm{w}_{p}).

Since $\bm{P}_{l}\bm{w}_{x}=\Phi^{\otimes^{d}}\bm{D}_{l}^{\mu}\bm{c}_{x}$ , the first step (4.14) gives

\frac{\mathrm{d}}{\mathrm{d}t}(\bm{c}_{x}\otimes\bm{w}_{p})+\mathrm{i}\sum\limits_{l=1}^{d}(\bm{D}_{l}^{\mu}\otimes\bm{D}_{l})(\bm{c}_{x}\otimes\bm{w}_{p})=\bm{0},

which can be written as

\frac{\mathrm{d}}{\mathrm{d}t}(\bm{c}_{x}\otimes\bm{w}_{p})+\mathrm{i}L(\bm{c}_{x}\otimes\bm{w}_{p})=\bm{0},

where

L=D_{\mu}\otimes I^{\otimes^{d-1}}\otimes D_{p}\otimes I^{\otimes^{d-1}}+\cdots+I^{\otimes^{d-1}}\otimes D_{\mu}\otimes I^{\otimes^{d-1}}\otimes D_{p}

is a diagonal matrix. Therefore the intermediate solution of the first step is

(\bm{c}_{x}\otimes\bm{w}_{p})^{*}=\mathrm{e}^{-\mathrm{i}L\Delta t}(\bm{c}_{x}\otimes\bm{w}_{p})^{n}.

Step 2. The second step is to solve (4.15), i.e.,

	$\displaystyle 0$	$\displaystyle=\partial_{t}w-\nabla_{x}V({x})\cdot\nabla_{p}w=\partial_{t}w-(\partial_{1}V({x}),\cdots,\partial_{d}V({x}))\cdot\nabla_{p}w$
		$\displaystyle=:\partial_{t}w-(v_{1}({x}),\cdots,v_{d}({x}))\cdot\nabla_{p}w.$

Similar to the first step, one has

(\bm{w}_{x}\otimes\bm{c}_{p})^{n+1}=\mathrm{e}^{\mathrm{i}U\Delta t}(\bm{w}_{x}\otimes\bm{c}_{p})^{*},

where

U=\bm{V}_{1}\otimes D_{\mu}\otimes I^{\otimes^{d-1}}+\cdots+\bm{V}_{d}\otimes I^{\otimes^{d-1}}\otimes D_{\mu}

is a diagonal matrix, and

\bm{V}_{k}=\text{diag}(\bm{v}_{k}),\quad\bm{v}_{k}=\sum\limits_{\bm{j}_{x}}v_{k}({x}_{\bm{j}_{x}})|j_{1}\rangle\cdots|j_{d}\rangle=\sum\limits_{\bm{j}_{x}}\partial_{k}V({x}_{\bm{j}_{x}})|j_{1}\rangle\cdots|j_{d}\rangle.

Given the initial state of $\bm{w}^{0}$ , applying the inverse QFT to the $x$ -register, one gets $(\bm{c}_{x}\otimes\bm{w}_{p})^{0}$ . At each time step, one needs to consider the following procedure

(\bm{c}_{x}\otimes\bm{w}_{p})^{n}\xrightarrow{\mathrm{e}^{-\mathrm{i}L\Delta t}}(\bm{c}_{x}\otimes\bm{w}_{p})^{*}\xrightarrow{F_{x}\otimes F_{p}^{-1}}(\bm{w}_{x}\otimes\bm{c}_{p})^{*}\xrightarrow{\mathrm{e}^{\mathrm{i}U\Delta t}}(\bm{w}_{x}\otimes\bm{c}_{p})^{n+1}\xrightarrow{F_{x}^{-1}\otimes F_{p}}(\bm{c}_{x}\otimes\bm{w}_{p})^{n+1},

where $F_{x}=F_{p}=\Phi^{\otimes^{d}}$ .

Theorem 4.4.

Given the error tolerance $\varepsilon$ , assume that $S_{0}(x)$ , $A_{0}(x)$ and $V(x)$ are smooth enough.

(1)

The Liouville equation can be simulated with gate complexity given by,

$N_{\text{Gates}}=\mathcal{O}\Big{(}\frac{d}{\varepsilon^{2}}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\Big{)}.$

(2)

The observable of the Liouville equation can be computed with gate complexity given by

N_{\text{Gates}}(\langle O\rangle)=\mathcal{O}\Big{(}\frac{n_{H}^{4}d}{\varepsilon^{4}}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\Big{)},

where $n_{H}=\|(\bm{\psi}^{\omega})^{0}\|/M^{d/2}$ .

Proof.

When $S_{0}(x)$ and $V(x)$ are smooth, the time-splitting spectral method has the error estimate

\|w^{n}(\cdot)-w(t_{n},\cdot)\|\leq C_{\ell}\Big{(}\omega+\frac{\Delta t}{\omega}+\frac{d\Delta x^{\ell}}{\omega^{\ell+1}}\Big{)},

(4.16)

where $\omega^{\ell+1}$ comes from the $\ell$ -th order derivative of $w:=w^{\omega}$ , and $C_{\ell}$ is an $\mathcal{O}(1)$ constant. Then one can implement the following meshing strategy

\omega\sim\varepsilon,\quad\Delta t\sim\varepsilon^{2},\quad\Delta x\sim\varepsilon^{1+2/\ell}/d^{1/\ell}

(4.17)

by forcing both error terms to be of order $\varepsilon$ . Thus,

M=L/\Delta x=2^{m}\quad\Longrightarrow\quad m=\log\frac{L}{\Delta x}=\mathcal{O}\Big{(}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\Big{)},

where $m$ is the number of qubits per dimension, and the total number of qubits is $m_{2d}=2dm$ . The diagonal unitary operators $\mathrm{e}^{-\mathrm{i}L\Delta t}$ and $\mathrm{e}^{\mathrm{i}U\Delta t}$ can be implemented using $J(m_{2d})=\mathcal{O}(m_{2d})$ gates, and the quantum Fourier transforms $F_{x}$ or $F_{p}$ can be implemented using $d\mathcal{O}(m\log m)$ gates. Therefore, the gate complexity required to iterate to the $n$ -th step is

	$\displaystyle N_{\text{Gates}}$	$\displaystyle=2n(J(m_{2d})+2d\mathcal{O}(m\log m))=2n\mathcal{O}(2dm+2dm\log m)$
		$\displaystyle=\mathcal{O}(4ndm\log m)=\mathcal{O}\Big{(}\frac{d}{\varepsilon^{2}}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\Big{)}.$

The gate complexity for the observable is obtained from the sampling law. ∎

Remark 4.2.

We remark that this section is the follow-up to [38]. Here we make precise statements about concrete implementations or running times of quantum algorithms. Compared to [38], we present a different (and simpler) algorithm for computing the physical observables and include the discussion of the spectral discretization. For the finite difference discretization, we also give details on the estimation of the condition number.

5 The Schrödinger framework

In this section we propose another framework to solve the important example of the nonlinear Hamilton-Jacbobi equations with $H(x,p)=\frac{1}{2}|p|^{2}+V(x)$ based on solving the Schrödinger equation, since the Liouville equation is the classical limit of the Schrödinger equation. The idea is to choose a semiclassical parameter, still denoted by $\hbar$ here, sufficiently small, so the solution of the Schrödinger equation is close to that of the Liouville equation.

Since the error between the expectation of the wave function and its classical counterpart (the physical observables of the Liouville equation) is of $O(\hbar^{2})$ [47], one can take $\hbar=\mathcal{O}(\sqrt{\varepsilon})$ , to maintain the computational precision of $\mathcal{O}(\varepsilon)$ for this framework.

We consider the Schrödinger equation in the semiclassical regime

\begin{cases}{\rm i}\hbar\partial_{t}u(t,{x})=-\frac{\hbar^{2}}{2}\Delta u(t,{x})+V({x})u(t,{x})\quad\mbox{in}~{}~{}\Omega=(a,b)^{d},\quad t>0,\\ u(0,{x})=u_{0}({x})\end{cases}

(5.1)

with periodic boundary conditions, where ${x}=(x_{1},x_{2},\cdots,x_{d})\in\mathbb{R}^{d}$ , $u(t,x):=u^{\hbar}(t,x)$ is the complex-valued wave function, $V({x})$ is the external potential and $\hbar=\mathcal{O}(\sqrt{\varepsilon})$ with $\varepsilon\ll 1$ being the precision. Without loss of generality, we always set $a=0$ and $b=1$ . The initial condition in (5.1) is chosen in a WKB form,

u_{0}({x})=A_{0}({x}){\rm e}^{{\rm i}\frac{S_{0}({x})}{\hbar}},

(5.2)

with $A_{0}$ and $S_{0}$ independent of $\hbar$ , real-valued and smooth. The periodic boundary conditions, for example, in one-dimensional case can be written as

u(t,a)=u(t,b),\quad u_{x}(t,a)=u_{x}(t,b),\quad t\geq 0.

The problem (5.1) will be solved by the classical time-splitting Fourier spectral method [4], which, as described below, can be interpreted as the Trotter based Hamiltonian simulation [37, 51].

5.1 The semiclassical approximation

We first recall the WKB analysis, which assumes that the solution remains the same form as the initial data at later time:

u=A(t,x){\rm e}^{{\rm i}\frac{S(t,x)}{\hbar}},

where $A(t,x)$ and $S(t,x)$ are the amplitude and phase respectively. Substituting this into (5.1), and separating the real and imaginary parts, one gets

	$\displaystyle A\partial_{t}S+\frac{1}{2}\|\nabla_{x}S\|^{2}+AV=\frac{\hbar^{2}}{2}\Delta A,$
	$\displaystyle\partial_{t}A+\nabla_{x}A\cdot\nabla_{x}S+\frac{1}{2}A\Delta S=0.$

Ignoring the $\mathcal{O}(\hbar^{2})$ terms, and multiplying the second equation by $A$ , one gets

	$\displaystyle\partial_{t}\|A\|^{2}+\nabla_{x}(\|A\|^{2}\nabla_{x}S)=0,$		(5.3)
	$\displaystyle\partial_{t}S+\frac{1}{2}\|\nabla_{x}S\|^{2}+V=0.$		(5.4)

The first equation (5.3) is a transport equation, and the second one (5.4) is the eikonal equation, which is exactly the Hamilton-Jacobi equation (4.1). Note that the eikonal equation admits solutions $S$ with discontinuous derivatives (usually referred to as the caustic) even if the initial data of $u$ is smooth. Thus the WKB analysis is only valid up to the time when the first caustic forms. Beyond caustics, the solution becomes multi-valued [36, 42].

In contrast to that, the Wigner transform technique yields the Liouville equation on phase space, in the semiclassical limit $\hbar\to 0$ , whose solution does not exhibit caustics, hence is valid globally in time. The Wigner transform of $u$ is defined as [36, 42]

w^{\hbar}(t,x,p)=w^{\hbar}[u](t,x,p):=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}u\Big{(}x+\frac{\hbar}{2}\eta\Big{)}\overline{u}\Big{(}x-\frac{\hbar}{2}\eta\Big{)}{\rm e}^{{\rm i}p\cdot\eta}{\rm d}\eta.

(5.5)

Applying this transformation on the Schrödinger equation (5.1), one obtains the Wigner equation (also called the quantum Liouville equation):

\partial_{t}w^{\hbar}+p\cdot\nabla_{x}w^{\hbar}-H_{V}w^{\hbar}=0,\quad w^{\hbar}(0,x,p)=w_{\rm in}(x,p),

where

H_{V}w^{\hbar}=\frac{{\rm i}}{(2\pi)^{d}}\iint_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\delta V(x,y)w^{\hbar}(x,p^{\prime}){\rm e}^{{\rm i}\eta(p-p^{\prime})}{\rm d}\eta{\rm d}p^{\prime},

\delta V=\frac{1}{\hbar}\Big{(}V(x-\frac{\hbar}{2}y)-V(x+\frac{\hbar}{2}y)\Big{)}.

When $\hbar\to 0$ , the Wigner equation becomes the classical Liouville equation on the phase space:

\partial_{t}w+p\cdot\nabla_{x}w-\nabla_{x}V(x)\cdot\nabla_{p}w=0,

which corresponds to the important case $H(x,p)=\frac{1}{2}|p|^{2}+V(x)$ . One easily finds that the Liouville equation can be written as $\partial_{t}w+\{w,H\}=0$ , where $\{\cdot,\cdot\}$ is the Poisson bracket, defined as

\{w,H\}=\nabla_{p}H\cdot\nabla_{x}w-\nabla_{x}H\cdot\nabla_{p}w.

When $u_{0}$ is given in WKB form (5.2), the corresponding Wigner measure is found to be

w^{\hbar}[u_{0}]\xrightarrow{\hbar\to 0}w_{0}=|A_{0}(x)|^{2}\delta(p-\nabla S_{0}(x)),

see [42, Eq. (3.9)]. It should be pointed out that Eqs. (5.3) and (5.4) can be deduced from the moment-closure of the Liouville problem

\begin{cases}\partial_{t}w+\nabla_{p}H\cdot\nabla_{x}w-\nabla_{x}H\cdot\nabla_{p}w=0,\\ w(0,x,p)=|A_{0}(x)|^{2}\delta(p-\nabla S_{0}(x))\end{cases}

(5.6)

with mono-kinetic ansatz $w(t,x,p)=|A(t,x)|^{2}\delta(p-\nabla_{x}S(t,x))$ , but are not valid beyond caustics since $\delta(p-\nabla_{x}S(t,x))$ is not well-defined when $\nabla_{x}S(t,x)$ becomes discontinuous. However, the equation in (5.6) is valid globally in time, since it unfolds the caustics in the phase space [36, 42, 38]. For this reason, we instead solve (5.6) in the semiclassical regime.

Clearly, (5.6) is the level set formulation (4.5) if $A_{0}(x)\equiv 1$ . Here we leave the general $A_{0}(x)$ . One can solve the problems in (5.6) by upwind finite difference methods or spectral methods as shown in the previous section.

5.2 Quantum simulations for the spectral discretisation

5.2.1 The time-splitting spectral approximations

From time $t=t_{n}$ to time $t=t_{n+1}$ , the Schrödinger equation is solved in two steps [4]: One solves

\hbar u_{t}-\mathrm{i}\frac{\hbar^{2}}{2}\Delta u=0

(5.7)

for one time step, followed by solving

\hbar u_{t}+\mathrm{i}V({x})u=0

(5.8)

again for one time step. Equation (5.7) will be discretised in space by the Fourier pseudo-spectral method and integrated in time exactly. The ODE (5.8) will then be solved exactly.

Remark 5.1.

We remark that usually the Trotter or Strang splitting is used for quantum simulation of the Schrödinger equation, which is second order in time rather than first order in the above simple splitting. We use the first order one in order to compare, on the equal footing, with other methods since all other methods use the first order approximation in time. For time complexity of the Trotter or Strang splitting see [37].

Let $u_{\bm{j}}^{n}$ be the numerical solution at $t=t_{n}$ and $u_{\bm{j}}^{*}$ the solution given by the first step for $0\leq\bm{j}\leq M-1$ .

•

For the first step, according to the previous discussion one easily obtains

\frac{\mathrm{d}}{\mathrm{d}t}\bm{u}(t)+\mathrm{i}\hbar/2\cdot(\bm{P}_{1}^{2}+\cdots+\bm{P}_{d}^{2})\bm{u}(t)=\bm{0},

\frac{\mathrm{d}}{\mathrm{d}t}\bm{c}(t)+\mathrm{i}\hbar/2\cdot((\bm{D}_{1}^{\mu})^{2}+\cdots+(\bm{D}_{d}^{\mu})^{2})\bm{c}(t)=\bm{0},

where the relation (3.38) is used, which gives

\bm{c}^{*}=\Big{(}\mathrm{e}^{-\mathrm{i}\hbar\Delta tD_{\mu}^{2}/2}\Big{)}^{\otimes^{d}}\bm{c}^{n},\qquad D_{\mu}={\rm diag}(\mu_{1},\cdots,\mu_{M}).

•

The updated numerical solution for the second step is $u_{\bm{j}}^{n+1}=\mathrm{e}^{-\mathrm{i}V({x}_{\bm{j}})\Delta t/\hbar}u_{\bm{j}}^{*}$ , which can be written in vector form as $\bm{u}=\mathrm{e}^{-\mathrm{i}\bm{V}\Delta t/\hbar}\bm{u}^{*}$ , where $\bm{V}$ is a diagonal matrix with

$\bm{V}_{n_{\bm{j}},n_{\bm{j}}}=V({x}_{\bm{j}}),\quad n_{\bm{j}}=j_{1}2^{d-1}+\cdots+j_{d}2^{0}.$

5.2.2 The quantum simulation of the Schrödinger equation

We only consider the 1-D case since it is straightforward to extend the arguments to high-dimensional cases by using tensor products. According to the previous discussion, we have the following algorithm (Algorithm 1) represented by the matrix-vector multiplication.

Algorithm 1 Time splitting approximations for the Schrödinger equation

1.

Given the initial data $\bm{u}^{0}$ and $n=0$ , compute the discrete Fourier coefficients $\bm{c}^{n}=\Phi^{-1}\bm{u}^{n}$ .
2.

Calculate the intermediate variables $\bm{c}^{*}=\mathrm{e}^{-\mathrm{i}\hbar D_{\mu}^{2}\Delta t/2}\bm{c}^{n}$ and $\bm{u}^{*}=\Phi\bm{c}^{*}$ .
3.

Update the numerical solution $\bm{u}^{n+1}=\mathrm{e}^{-\mathrm{i}\bm{V}\Delta t/\hbar}\bm{u}^{*}$ .

In the above algorithm, the matrix $\Phi$ plays the role of the discrete Fourier transform (DFT), where given a set of numbers $x_{0},x_{1},\cdots,x_{M-1}$ , the DFT and the inverse DFT are defined by

y_{k}=\frac{1}{\sqrt{M}}\sum\limits_{j=0}^{M-1}\mathrm{e}^{2\pi{\rm i}jk/M}x_{j},\quad k=0,\cdots,M-1

and

x_{j}=\frac{1}{\sqrt{M}}\sum\limits_{k=0}^{M-1}\mathrm{e}^{-2\pi{\rm i}jk/M}y_{k},\quad j=0,\cdots,M-1,

respectively. Denote the transformation matrix of DFT by $F$ . It is easy to find the transformation matrix in Algorithm 1 satisfies $\Phi=\sqrt{M}SF$ , where $S$ is the diagonal matrix

S={\rm diag}\Big{(}[1,-1,\cdots,1,-1]_{M\times 1}\Big{)},

which in turn gives

\bm{u}^{n+1}=\mathrm{e}^{-\mathrm{i}\bm{V}\Delta t/\hbar}SF\mathrm{e}^{-\mathrm{i}\hbar D_{\mu}^{2}\Delta t/2}F^{-1}S\bm{u}^{n}

since $S^{-1}=S$ . For convenience, the above two diagonal matrices are denoted by $D_{1}$ and $D_{2}$ , respectively, and then one has

\displaystyle\bm{u}^{n+1}=D_{1}SFD_{2}F^{-1}S\bm{u}^{n}=SD_{1}FD_{2}F^{-1}S\bm{u}^{n},

\bm{v}^{n+1}=(D_{1}FD_{2}F^{-1})\bm{v}^{n},\quad\bm{v}^{n}=S\bm{u}^{n},

where we have used the fact that $SD_{1}=D_{1}S$ .

Therefore, when preparing the variable $\bm{v}$ in the computational basis, the implementation in each iteration involves one application of an inverse quantum Fourier transform (QFT), followed by a multiplication of a diagonal unitary operator $D_{2}$ , and a QFT and another diagonal unitary operator $D_{1}$ , since the QFT is exactly the quantum version of the DFT.

According to the above discussion, the quantum simulation algorithm to find $\bm{v}^{n}:=S\bm{u}^{n}$ is described as follows (see Algorithm 2).

Algorithm 2 Quantum simulation of the Schrödinger equation

Step 0.

Initialization of the quantum state: Given $v_{j}^{0}$ encode it as

$|\psi^{0}\rangle=\frac{1}{\mathcal{N}}\sum\limits_{j=0}^{M-1}\psi_{j}^{0}|j\rangle,\quad\psi_{j}^{0}=v_{j}^{0},$

where $\mathcal{N}$ is the normalization constant. Let $n=0$ .
Step 1.

Performing inverse QFT on $|\psi^{n}\rangle$ yields $|\tilde{\psi}\rangle$ .
Step 2.

Perform a diagonal unitary operator $\Big{(}\mathrm{e}^{-\mathrm{i}\hbar\mu_{l}^{2}\Delta t/2}\Big{)}_{1\leq l\leq M}$ for $|\tilde{\psi}\rangle$ , and the resulting state is denoted as $|\psi^{\mu}\rangle$ .
Step 3.

Perform QFT on $|\psi^{\mu}\rangle$ with the output denoted by $|\psi^{*}\rangle$ .
Step 4.

Apply a diagonal unitary operator $\Big{(}\mathrm{e}^{-\mathrm{i}V(x_{j})\Delta t/\hbar}\Big{)}_{0\leq j\leq M-1}$ to $|\psi^{*}\rangle$ . The output is denoted by $|\psi^{n+1}\rangle$ .
Step 5.

Let $n\leftarrow n+1$ and go back to Step 1.

The following example is taken from Example 1 in [4].

Example 5.1.

The initial data is $u(x,0)=A_{0}(x)\mathrm{e}^{\mathrm{i}S_{0}(x)/\hbar}$ , where

A_{0}(x)=\mathrm{e}^{-25(x-0.5)^{2}},\quad S_{0}(x)=-\frac{1}{5}\ln\Big{(}\mathrm{e}^{5(x-0.5)}+\mathrm{e}^{-5(x-0.5)}\Big{)}.

We take $[a,b]=[0,1]$ and $V(x)=10$ . The position density $\rho(t,x)=|u(t,x)|^{2}$ is shown in Fig. 2. One can see that the solution is oscillatory for small $\hbar$ . For numerical descriptions, please refer to [4].

5.2.3 Gate counts for the computation of wave functions

To simplify the discussion, we set $L=b-a=1$ and the evolution time $t=1$ throughout the paper. The time-splitting scheme involves only diagonal operators and QFTs whose complexities depend on the number of qubits $m$ per dimension. Since the meshing satisfies $\Delta x=L/M$ and $M=2^{m}$ , we can determine $\Delta x$ and thus $m$ by the expected precision of the algorithm.

Theorem 5.1.

Given the error tolerance $\varepsilon$ , assume $V({x})$ is sufficiently smooth and $\hbar=\mathcal{O}(\sqrt{\varepsilon})$ . Then the Schrödinger equation (5.1) can be simulated using $m_{d}$ qubits with gate complexity $N_{\text{Gates}}$ , given respectively by,

m_{d}=\mathcal{O}\Big{(}d\log\frac{d^{1/\ell}}{\varepsilon^{1/2+5/(2\ell)}}\Big{)},\qquad N_{\text{Gates}}=\mathcal{O}\Big{(}\frac{d}{\varepsilon^{3/2}}\log\frac{d^{1/\ell}}{\varepsilon^{1/2+5/(2\ell)}}\Big{)}.

Proof.

According to Refs. [4, 37], if $\Delta x/\hbar=\mathcal{O}(1)$ and $\Delta t/\hbar=\mathcal{O}(1)$ , then for each $\ell$ , the error of the Fourier spectral method is bounded by,

\|u^{n}-u(t_{n},\cdot)\|\leq C_{\ell}n\Big{(}d\Big{(}\frac{\Delta x}{\hbar}\Big{)}^{\ell}+\frac{\Delta t^{2}}{\hbar}\Big{)}.

(5.9)

Here $C_{\ell}$ is considered as $\mathcal{O}(1)$ . The mesh strategy is

\frac{\Delta t}{\hbar}\sim\varepsilon,\qquad\frac{\Delta x}{\hbar}\sim\Big{(}\frac{\varepsilon}{nd}\Big{)}^{1/\ell},\qquad\hbar\sim\sqrt{\varepsilon},

or equivalently,

\Delta t\sim\varepsilon^{3/2},\qquad\Delta x\sim\varepsilon^{1/2+5/(2\ell)}/d^{1/\ell},

(5.10)

which is obtained by forcing both error terms to be of order $\varepsilon$ . With this choice of $\Delta x$ and $\Delta t$ , the number of qubits is

m_{d}=dm,\quad m=\log\frac{1}{\Delta x}=\mathcal{O}\Big{(}\log\frac{d^{1/\ell}}{\varepsilon^{1/2+5/(2\ell)}}\Big{)}.

The diagonal unitary operators can be implemented using $J(m_{d})=\mathcal{O}(m_{d})$ gates [45, 37]. Thus the gate complexity required to iterate to the $n$ -th step is

\displaystyle N_{\text{Gates}}\sim n(2dm\log m+2J(m_{d}))\sim 2ndm\log m=\mathcal{O}\Big{(}\frac{d}{\varepsilon^{3/2}}\log\frac{d^{1/\ell}}{\varepsilon^{1/2+5/(2\ell)}}\Big{)},

(5.11)

as required. ∎

5.3 The QLSA for the spectral discretisation

According to the discussion in Subsect. 5.2.1, one has the following ODEs

\frac{\mathrm{d}}{\mathrm{d}t}\bm{u}(t)=A\bm{u}(t),

where

A=-\mathrm{i}\Big{(}\frac{\hbar}{2}(\bm{P}_{1}^{2}+\cdots+\bm{P}_{d}^{2})+\frac{1}{\hbar}\bm{V}\Big{)}=:-\mathrm{i}\tilde{A}.

(5.12)

Obviously, $\tilde{A}$ is a real symmetric matrix.

As in Theorem 4.3, we can apply the quantum algorithm in [7] to solve the above ODEs.

Theorem 5.2.

Assume that $\|\bm{V}\|=\mathcal{O}(1)$ and $T=\mathcal{O}(1)$ . Then there exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{u}(T)/\|\bm{u}(T)\|$ with the gate complexity given by

N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{1+5/\ell}}\Big{)}.

Proof.

Let $A=V^{-1}DV$ with $D=\text{diag}(\lambda_{1},\cdots,\lambda_{N})$ , where $N=\mathcal{O}(M^{d})$ is the order of $A$ . Since the matrix $\tilde{A}$ in (5.12) is a real symmetric matrix, $\text{Re}(\lambda_{j})=0$ for any $j\in\{1,\cdots,N\}$ and $\kappa_{V}=1$ . According to Lemma 3.4, there exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{u}(T)/\|\bm{u}(T)\|$ with the gate complexity given by

N_{Gates}=\mathcal{O}\Big{(}s\|A\|T\cdot\text{Poly}(\log(sN\|A\|T/\varepsilon))\Big{)}.

It is evident that the sparsity of $A$ is $\mathcal{O}(M)=\mathcal{O}(1/\Delta x)$ . For the mesh strategy in (5.10), the norm of $A$ satisfies

\displaystyle\|A\|\leq\frac{\hbar}{2}\sum\limits_{l=1}^{d}\|\bm{P}_{l}\|+\frac{1}{\hbar}\|\bm{V}\|\lesssim\hbar dM+\frac{1}{\hbar}\lesssim\frac{d^{1+1/\ell}}{\varepsilon^{5/(2\ell)}}+\frac{1}{\varepsilon^{1/2}}\leq\frac{d^{1+1/\ell}}{\varepsilon^{1/2+5/(2\ell)}}.

Then one has

N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{1+5/\ell}}\Big{)},

where in the last equal sign we have included the additional factor $d$ arising from the matrix order (see Remark 3.8). This completes the proof. ∎

5.4 The computation of physical observables

The quantum mechanical wave function $u(t,x)$ can be considered as an auxiliary quantity used to compute physical quantities. The most basic quadratic observables [4, 42, 36] include the position density $\rho(t,x):=|u(t,x)|^{2}$ , the current density

J(t,x)=\hbar{\rm Im}(\overline{u}(t,x)\nabla u(t,x))=\frac{\hbar}{2\mathrm{i}}(\overline{u}\nabla u-u\nabla\overline{u}),

and the kinetic or total energy

E(t,x)=\frac{\hbar^{2}}{2}|\nabla u(t,x)|^{2}\qquad\mbox{or}\qquad\frac{\hbar^{2}}{2}|\nabla u(t,x)|^{2}+V(x)|\nabla u(t,x)|^{2}.

5.4.1 The expectation of observables for the quantum simulation

The observables $\rho(t_{n},x_{\bm{i}})$ , $J(t_{n},x_{\bm{i}})$ and $E(t_{n},x_{\bm{i}})$ can be expressed as the standard form of $\langle O\rangle=\langle\widetilde{\bm{u}}|O|\widetilde{\bm{u}}\rangle$ , which is the expectation of the observable $O$ . Here, $\widetilde{\bm{u}}$ is normalized to $1$ since the output of the quantum algorithms is the normalized state. For the spectral discretisation of the Schrödinger equation, one has

\|\bm{u}^{N_{t}}\|=\cdots=\|\bm{u}^{0}\|=:N_{u_{0}},

(5.13)

and $\widetilde{\bm{u}}=\bm{u}^{N_{t}}/N_{u_{0}}$ . The output state of the quantum simulation is

|\widetilde{\bm{u}}\rangle=\frac{1}{N_{u_{0}}}\sum\limits_{\bm{j}}\bm{u}^{N_{t}}|\bm{j}\rangle=\sum\limits_{\bm{j}}\widetilde{\bm{u}}|\bm{j}\rangle,\qquad\bm{j}=(j_{1},\cdots,j_{d}).

For the position density $\rho(t_{n},x_{\bm{i}})$ , it is obvious that one can measure the magnitude of the wave function using multiple shots in the computational basis $|\bm{i}\rangle$ as was done in the numerical experiments in [37], so we choose

O_{\rho}=:N_{u_{0}}^{2}|\bm{i}\rangle\langle\bm{i}|\quad\mbox{for the position density}.

Let $p=|\widetilde{\bm{u}}_{\bm{i}}|^{2}\leq 1$ . Then ${\rm Var}(O_{\rho})=N_{u_{0}}^{4}p(1-p)\leq N_{u_{0}}^{4}$ , hence the number of samples is $n_{\rho}=N_{u_{0}}^{4}/\varepsilon^{2}$ for the position density.

For the current density, we consider $d=1$ for simplicity and denote $\bm{e}_{i}=|i\rangle$ . Let $\bm{u}=\bm{u}^{N_{t}}$ for simplicity. Using the previous notations, one has

\overline{u}(x_{i})=\overline{\bm{u}}^{T}\bm{e}_{i},\qquad(\partial_{x}u)(x_{i})=(\Phi D_{\mu}\Phi^{-1}\bm{u})^{T}\bm{e}_{i}=\bm{e}_{i}^{T}(\Phi D_{\mu}\Phi^{-1}\bm{u}),

which gives

(\overline{u}\partial_{x}u)(x_{i})=\overline{\bm{u}}^{T}(\bm{e}_{i}\bm{e}_{i}^{T}\Phi D_{\mu}\Phi^{-1})\bm{u}=\bm{u}^{\dagger}(\bm{e}_{i}\bm{e}_{i}^{T}\Phi D_{\mu}\Phi^{-1})\bm{u}=:\bm{u}^{\dagger}A\bm{u},

and

(u\partial_{x}\overline{u})(x_{i})=\Big{(}(\overline{u}\partial_{x}u)(x_{i})\Big{)}^{\dagger}=\bm{u}^{\dagger}A^{\dagger}\bm{u}.

Therefore,

J(t_{n},x_{i})=\frac{\hbar}{2\mathrm{i}}\bm{u}^{\dagger}(A-A^{\dagger})\bm{u}=\frac{\hbar N_{u_{0}}^{2}}{2\mathrm{i}}\widetilde{\bm{u}}^{\dagger}(A-A^{\dagger})\widetilde{\bm{u}},

and we can choose

O_{J}:=\frac{\hbar N_{u_{0}}^{2}}{2\mathrm{i}}(A-A^{\dagger})\quad\mbox{for the current density}.

One can check that

{\rm Var}(O_{J})\leq\|O_{J}\widetilde{\bm{u}}\|^{2}\lesssim\Big{(}\frac{\hbar MN_{u_{0}}^{2}}{4}\Big{)}^{2},

where $M$ comes from $D_{\mu}={\rm diag}(-N,\cdots,N-1)$ with $N=M/2$ , which implies a multiplicative factor $n_{J}=(\hbar MN_{u_{0}}^{2})^{2}/\varepsilon^{2}=M^{2}N_{u_{0}}^{4}/\varepsilon$ for the current density since $\hbar=\mathcal{O}(\sqrt{\varepsilon})$ .

For the kinetic energy, we similarly obtain

E(t_{n},x_{i})=\frac{\hbar^{2}}{2}\bm{u}^{\dagger}(\Phi D_{\mu}\Phi^{-1}\bm{e}_{i}\bm{e}_{i}^{T}\Phi D_{\mu}\Phi^{-1})\bm{u}=:\frac{\hbar^{2}}{2}\bm{u}^{\dagger}B\bm{u}=\frac{\hbar^{2}N_{u_{0}}^{2}}{2}\widetilde{\bm{u}}^{\dagger}B\widetilde{\bm{u}}

for the one-dimensional case. We can choose

O_{E}:=\frac{\hbar^{2}N_{u_{0}}^{2}}{2}B\quad\mbox{for the kinetic energy}.

It is obvious that

{\rm Var}(O_{E})\leq\|O_{E}\widetilde{\bm{u}}\|^{2}\lesssim\Big{(}\frac{\hbar^{2}M^{2}N_{u_{0}}^{2}}{4}\Big{)}^{2},

which implies a multiplicative factor $n_{E}=(\hbar^{2}M^{2}N_{u_{0}}^{2})^{2}/\varepsilon^{2}=M^{4}N_{u_{0}}^{4}$ for the kinetic energy.

Remark 5.2.

By definition, $\rho(t,x):=|u(t,x)|^{2}$ , where

\rho(t,x)=\int w(t,x,p)\mathrm{d}p\approx\frac{1}{M^{d}}\sum_{\bm{l}}G_{\bm{l}}w^{\omega}_{\bm{j},\bm{l},n},

where $w$ is the solution to (5.6). Then

N_{u_{0}}^{2}=\|\bm{u}^{0}\|^{2}=\|\bm{\rho}^{0}\|\lesssim\|\bm{G}\|\|(\bm{w}^{\omega})^{0}\|/M^{d}\lesssim\|(\bm{w}^{\omega})^{0}\|/M^{d/2}=n_{H}.

Here, $(\bm{w}^{\omega})^{0}$ is exactly the $(\bm{\psi}^{\omega})^{0}$ in Theorem 4.2 when considering the problem (5.6).

5.4.2 The expectation of observables for the QLSA

Unlike the quantum simulation, the solution of the QLSA is a quantum state that is a superposition of the solution at all temporal and spatial points, denoted as

|\widetilde{\bm{u}}\rangle=[\widetilde{\bm{u}}^{1};\cdots;\widetilde{\bm{u}}^{N_{t}}],\qquad\widetilde{\bm{u}}^{n}=\frac{1}{N_{u}}\bm{u}^{n},

where the normalization constant is

N_{u}=\|\bm{u}\|=(\|\bm{u}^{1}\|^{2}+\cdots+\|\bm{u}^{N_{t}}\|^{2})^{1/2}=\sqrt{N_{t}}\|\bm{u}^{0}\|=:\sqrt{N_{t}}N_{u_{0}}.

Here we have used (5.13). Let $O_{\bm{i}}=|\bm{i}\rangle\langle\bm{i}|$ , $O_{n}=|n\rangle\langle n|$ and

O_{\bm{i}}^{n}=O_{n}\otimes O_{\bm{i}}=|n,\bm{i}\rangle\langle n,\bm{i}|,

where $|n\rangle$ is of size $N_{t}$ . Then the position density

\rho(t=t_{n},x_{\bm{i}})=(\bm{u}^{n})^{\dagger}O_{\bm{i}}\bm{u}^{n}=\bm{u}^{\dagger}(O_{n}\otimes O_{\bm{i}})\bm{u}=N_{t}N_{u_{0}}^{2}\cdot\langle\widetilde{\bm{u}}|O_{\bm{i}}^{n}|\widetilde{\bm{u}}\rangle.

The expectation $\langle O_{\bm{i}}^{n}\rangle:=\langle\widetilde{\bm{u}}|O_{\bm{i}}^{n}|\widetilde{\bm{u}}\rangle$ satisfies the condition that ${\rm Var}(O_{\bm{i}}^{n})$ is bounded. In this case, however, we must evaluate $\langle O_{\bm{i}}^{n}\rangle$ to precision $\mathcal{O}(\varepsilon/(N_{t}N_{u_{0}}^{2}))$ , which increases the number of samples by another factor $(N_{t}N_{u_{0}}^{2})^{2}$ when considering the general sampling law. We remark that the factor $N_{t}^{2}$ can be removed by using the dilation procedure. In this case, the multiplicative factor is still given by $n_{\rho}=N_{u_{0}}^{4}/\varepsilon^{2}$ .

For the current density, one easily finds that ( $d=1$ )

J(t_{n},x_{i})=\frac{\hbar}{2\mathrm{i}}(\bm{u}^{n})^{\dagger}(A-A^{\dagger})\bm{u}^{n}=N_{t}N_{u_{0}}^{2}\cdot\frac{\hbar}{2\mathrm{i}}\widetilde{\bm{u}}^{\dagger}(A-A^{\dagger})\widetilde{\bm{u}}.

We still need to apply the dilation procedure (3.28) to remove the unexpected multiplicative factor $N_{t}^{2}$ , and still obtain $n_{J}=M^{2}N_{u_{0}}^{4}/\varepsilon$ for the current density.

The kinetic energy can be analysed similarly, with the factor given by $n_{E}=M^{4}N_{u_{0}}^{4}$ .

5.5 Gate complexity for the computation of the observables

It is worth pointing out that the time step $\Delta t$ can be chosen independently of the small parameter $\hbar$ if one is only concerned with the computation of the physical quantities. This observation was interpreted by using the Wigner transformation approach in [4], and mathematically rigorously investigated in [31, 12]. For instance, the first-order time splitting spectral method gives [37, 12]

\|\langle O\rangle_{u^{n}}-\langle O\rangle_{u(t_{n},\cdot)}\|\leq C_{\ell}n\Big{(}d\Big{(}\frac{\Delta x}{\hbar}\Big{)}^{\ell}+\Delta t^{2}+\Delta t\hbar^{2}\Big{)},

(5.14)

where $\hbar=\mathcal{O}(\sqrt{\varepsilon})$ implies the above observation since $n\Delta t\hbar^{2}=\mathcal{O}(\varepsilon)$ . Note that the term $n\Delta t\hbar^{2}=t\hbar^{2}$ is the error between the classically evolved Wigner function and the expectation value of the Schrödinger solution [12].

5.5.1 The gate counts of the quantum simulation method

Theorem 5.3.

Given the error tolerance $\varepsilon$ , suppose that the estimate (5.14) holds with $C_{\ell}$ considered as $\mathcal{O}(1)$ . The gate complexities for the observables from the Schrödinger equation (5.1) are given by

N_{Gates}(\langle O\rangle)=\mathcal{O}\Big{(}c_{O}\frac{N_{u_{0}}^{4}d}{\varepsilon^{3}}\log\frac{d^{1/\ell}}{\varepsilon^{1/2+2/\ell}}\Big{)},

where

c_{O}=\begin{cases}1,&\text{for the position density},\\ d^{2/\ell}/\varepsilon^{4/\ell},&\text{for the current density},\\ d^{4/\ell}/\varepsilon^{8/\ell},&\text{for the kinetic energy}.\end{cases}

(5.15)

Proof.

For the observables, one can implement the mesh strategy according to (5.14), given by

\Delta t=\mathcal{O}(\varepsilon),\quad\Delta x=\mathcal{O}(\varepsilon^{1/2+2/\ell}/d^{1/\ell}).

(5.16)

Then the number of gates for outputting the quantum sate is

\mathcal{O}\Big{(}\frac{d}{\varepsilon}\log\frac{d^{1/\ell}}{\varepsilon^{1/2+2/\ell}}\Big{)}.

For computing the observables, we must add a multiplicative factor $\text{Var}(O)/\varepsilon^{2}$ if the general sampling law is used. According to the previous discussion, one has

n_{\rho}=\frac{N_{u_{0}}^{4}}{\varepsilon^{2}},\qquad n_{J}=\frac{M^{2}N_{u_{0}}^{4}}{\varepsilon}\sim\frac{N_{u_{0}}^{4}d^{2/\ell}}{\varepsilon^{2+4/\ell}},\qquad n_{E}=M^{4}N_{u_{0}}^{4}\sim\frac{N_{u_{0}}^{4}d^{4/\ell}}{\varepsilon^{2+8/\ell}}.

(5.17)

This completes the proof. ∎

5.5.2 The gate counts of the quantum linear systems algorithm

Theorem 5.4.

Given the error tolerance $\varepsilon$ , suppose that the estimate (5.14) holds with $C_{\ell}$ considered as $\mathcal{O}(1)$ . The gate complexities of the QLSA for the observables from the Schrödinger equation (5.1) are given by

N_{Gates}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}c_{O}\frac{N_{u_{0}}^{4}d^{2+2/\ell}}{\varepsilon^{3+4/\ell}}\Big{)},

where $c_{O}$ is defined by (5.15).

Proof.

For the mesh strategy in (5.16), according to the proof of Theorem 5.2, one easily finds that the number of gates for approximating the wave function is

N_{\text{Gates}}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{1+4/\ell}}\Big{)},

where in the last equal sign we have included the additional factor $d$ arising from the matrix order. For computing the observables, one just need to include the multiplicative factor $\text{Var}(O)/\varepsilon^{2}$ as given in (5.17). ∎

6 Linear representation approach for scalar nonlinear hyperbolic PDEs

We consider the linear representations for the spatially $d$ -dimensional scalar nonlinear hyperbolic PDE

\begin{cases}\partial_{t}u+F(u)\cdot\nabla_{x}u+Q(x,u)=0,\\ u(0,x)=u_{0}(x),\end{cases}

(6.1)

where $x\in\mathbb{R}^{d}$ and $u\in\mathbb{R}$ .

6.1 The Liouville representation

For this general scalar hyperbolic equation, one can still use the Liouville representation but the Schrödinger representation is not available. In fact, the Liouville representation has been considered in [38] by using the level set formalism. To this end, we first review the construction.

Let $\phi(t,x,p)$ be the level set function in $(d+1)+1=d+2$ dimensions, where $p\in\mathbb{R}$ . The zero level set of $\phi$ gives solution $u$ :

\phi(t,x,p)=0\quad\mbox{at}\quad p=u(t,x).

One easily finds that $\phi$ satisfies

\begin{cases}\partial_{t}\phi+F(p)\cdot\nabla_{x}\phi-Q(x,p)\partial_{p}\phi=0,\\ \phi(0,x,p)=p-u_{0}(x).\end{cases}

Like for the Hamilton-Jacobi PDEs, we can similarly define a function $\varphi$ such that

\begin{cases}\partial_{t}\varphi+F(p)\cdot\nabla_{x}\varphi-Q(x,p)\partial_{p}\varphi=0,\\ \varphi(0,x,p)=\delta(p-u_{0}(x)),\end{cases}

(6.2)

with the solution given by $\varphi(t,x,p)=\delta(\phi(t,x,p))$ . Eq. (6.2) is referred to as the Liouville representation of (6.1).

One can apply the quantum difference method to solve (6.2) and compute the physical observables as in Subsect. 4.1. For the spectral discretisation, one can utilize the Trotter based technique in Subsect. A.1. The similar numerical performance can be deduced, so we omit the detailed discussions in view of the length of the article and the similarity of the numerical implementation.

6.2 The KvN representation

For scalar nonlinear hyperbolic PDEs, it is not very clear how to formulate the KvN representation. Here we offer an idea, which is inspired by the evolution of the phase factor of the KvN wave function in [44].

Let $\bm{x}=(x,p)$ and denote $\bm{v}=[F(p),-Q(x,p)]$ . Then problem (6.2) can be written as

\begin{cases}\partial_{t}\varphi+\bm{v}\cdot\nabla_{\bm{x}}\varphi=0,\\ \varphi_{0}(\bm{x})=\varphi(0,\bm{x})=\varphi(0,x,p)=\delta(p-u_{0}(x)).\end{cases}

(6.3)

Let $\psi$ be the complex-valued KvN wave function and $f=\psi^{\dagger}\psi$ be the probability distribution function. Inspired by the discussion in Section II(B) of [44], we define $\psi=f^{1/2}\mathrm{e}^{\mathrm{i}\varphi}$ :

•

The amplitude $f$ satisfies the Liouville equation in conservative form:

$\begin{cases}\partial_{t}f+\nabla_{\bm{x}}\cdot(\bm{v}f)=0,\\ f(0,\bm{x})=f_{0}(\bm{x}),\end{cases}$ (6.4)

where $f_{0}(\bm{x})=\delta(\bm{x}-\bm{q}_{0})$ and $\bm{q}_{0}$ is an arbitrarily given vector.
•

The phase factor $\varphi$ satisfies the non-conservative Liouville equation in (6.3).

As in [44], one can check that the KvN wave function $\psi$ is governed by the KvN equation

\mathrm{i}\partial_{t}\psi=\mathcal{H}_{\text{KvN}}\psi=-\mathrm{i}\Big{(}\bm{v}\cdot\nabla_{\bm{x}}+\frac{1}{2}\nabla_{\bm{x}}\cdot\bm{v}\Big{)}\psi,

(6.5)

with the initial data given by

\psi_{0}(\bm{x})=\psi(0,\bm{x})=f_{0}^{1/2}(\bm{x})\mathrm{e}^{\mathrm{i}\varphi_{0}(\bm{x})}.

The original intention of [44] is to introduce the linear representation for the following nonlinear ODEs

\frac{\mathrm{d}\bm{q}}{\mathrm{d}t}=\bm{v}(t,\bm{q}).

(6.6)

Like in Sect. 3, one can first obtain the Liouville representation (6.4) corresponding to (6.6). The KvN representation is then derived by assuming the non-conservative hyperbolic equation (6.2) for the phase factor. In view of the evolution of the phase factor, we therefore propose the KvN representation for the scalar nonlinear hyperbolic PDEs.

It should be pointed out that one cannot get the phase factor $\varphi$ from $\psi=f^{1/2}\mathrm{e}^{\mathrm{i}\varphi}$ since $\mathrm{e}^{\mathrm{i}\varphi}$ is periodic with respect to $\varphi$ . For this reason, it may be impossible to define the physical observables

\langle g(t,x)\rangle=\int g(p)\varphi(t,x,p)\mathrm{d}p

as in [38]. One may instead define

\langle g_{\chi}(t,x)\rangle=\int g(p)\chi(\varphi(t,x,p))\mathrm{d}p

to remove the unexpected period, where $\chi$ is a function with period $2\pi$ . However, $\langle g_{\chi}(t,x)\rangle$ is actually very tricky to compute. It’s not clear if this is possible to do. Suppose we have $\chi(\varphi)=\sin(\varphi)$ or $\cos(\varphi)$ . If we have access to a Hamiltonian diagonal in $\varphi$ , i.e. $H=\sum_{j}\varphi_{j}|j\rangle\langle j|$ , then we can create a control unitary $U=|0\rangle\langle 0|\otimes\mathbf{1}+|1\rangle\langle 1|\otimes\mathrm{e}^{\mathrm{i}H}$ . Then a standard Hadamard test can be employed to compute $\langle g_{\chi}(t,x)\rangle=\text{Im}\langle\sqrt{g}|U|\sqrt{g}\rangle$ or $\text{Re}\langle\sqrt{g}|U|\sqrt{g}\rangle$ where $|\sqrt{g}\rangle$ has amplitudes $\sqrt{g_{j}}$ . However, by solving Eq. (6.5) alone, one cannot create access to $H$ unless every $\varphi_{j}$ is individually computed, which defeats the purpose of a quantum algorithm.

Instead we can try to create the state $|\sin(\varphi)\rangle$ or $|\cos(\varphi)\rangle$ whose amplitudes are proportional to $\sin(\varphi_{j})$ or $\cos(\varphi_{j})$ . Note that

\sin\varphi_{j}=\text{Im}(\psi_{j}/\sqrt{f_{j}})=\text{Im}(\psi_{j}/\sqrt{\psi_{j}\psi_{j}^{\dagger}})=\frac{1}{2\mathrm{i}}\frac{\psi_{j}-\psi_{j}^{\dagger}}{\sqrt{\psi_{j}\psi_{j}^{\dagger}}},\qquad\cos(\varphi_{j})=\frac{\psi_{j}+\psi_{j}^{\dagger}}{2\sqrt{\psi_{j}\psi_{j}^{\dagger}}}.

The question is if we can create a quantum state with amplitudes proportional to $\text{Im}(\psi_{j}/\sqrt{f_{j}})$ or $\text{Re}(\psi_{j}/\sqrt{f_{j}})$ . If we solve Eq. (6.5) alone to obtain $|\psi\rangle$ , then we cannot easily access these states.

An alternative to solving Eq. (6.5) alone is to include also its complex conjugate and we instead solve for $\widetilde{\psi}$ obeying

\mathrm{i}\frac{\partial}{\partial t}\widetilde{\psi}=\widetilde{\mathcal{H}}_{\text{KvN}}\widetilde{\psi}

where

\widetilde{\psi}=\begin{bmatrix}\psi-\psi^{\dagger}\\ \psi+\psi^{\dagger}\end{bmatrix},\qquad\widetilde{\mathcal{H}}_{\text{KvN}}=\begin{bmatrix}\mathcal{H}_{\text{KvN}}&0\\ 0&\mathcal{H}_{\text{KvN}}\end{bmatrix}.

We now perform quantum simulation of the state $|\widetilde{\psi}\rangle\propto|\psi-\psi^{\dagger}\rangle+|\psi+\psi^{\dagger}\rangle$ with the new Hamiltonian $\widetilde{\mathcal{H}}_{\text{KvN}}$ and we can choose to post-select either the state $|\psi-\psi^{\dagger}\rangle$ or $|\psi+\psi^{\dagger}\rangle$ . Given $|\psi\pm\psi^{\dagger}\rangle$ , we can compute the inner product $|\langle g|\psi\pm\psi^{\dagger}\rangle|^{2}$ to obtain the observable $\langle g_{\chi}(t,x)\rangle$ for $\chi(\varphi)=|\psi\pm\psi^{\dagger}|$ using standard methods. This differs from $\chi(\varphi)=\sin(\varphi)$ or $\cos(\varphi)$ up to norm factors $\sqrt{\psi\psi^{\dagger}}$ . We leave it as an open question on how to design algorithms for more general $\chi(\varphi)$ functions.

7 Summary and discussion

In this paper, we systematically studied the quantum difference methods and the quantum spectral methods for solving the linear representations of nonlinear ODEs and nonlinear PDEs. Since our studies involve many different methods, we summarize the results for computing the physical observables in Tab. 1, from which we clearly observe that the quantum simulation methods give the best performance in the computational cost. On the other hand, for the classical treatment the order of the iterative matrix in $\bm{u}^{n+1}=B\bm{u}^{n}$ scales as $\mathcal{O}(N_{x}^{d})$ , leading to the classical run time scaling as $C=\mathcal{O}(d^{d+3}/\varepsilon^{d+1})$ for the first-order hyperbolic equation discussed in [41], while the quantum cost is $Q=\mathcal{O}(d^{3}\log(1/\varepsilon))$ , where $\varepsilon=\mathcal{O}(d/N_{x})$ and the quantum difference method is employed. This implies the exponential speedup with respect to $d$ and $\varepsilon$ even for producing the observables.

Tab. 1: Time complexities for the computation of physical observables

Nonlinear Ordinary Differential Equations
Problem	Quantum simulation		Spectral QLSA		FD QLSA ( $\alpha\geq 4$ )
Problem	Subroutine	Observable	Subroutine	Observable	Subroutine	Observable
Liouville representation	$\dfrac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}^{*}$	$\dfrac{n_{L}^{4}d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}^{*}$	$\dfrac{d^{3+2/\ell}}{\varepsilon^{4+4/\ell}}^{*}$	$\dfrac{n_{L}^{4}d^{3+2/\ell}}{\varepsilon^{6+4/\ell}}^{*}$	$\dfrac{d^{\alpha}}{\varepsilon^{3}}$	$\dfrac{n_{L}^{4}d^{\alpha}}{\varepsilon^{5}}$
KvN representation	$\dfrac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}$	$\dfrac{n_{L}^{2}d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}$	$\dfrac{d^{3+2/\ell}}{\varepsilon^{2+4/\ell}}$	$\dfrac{n_{L}^{2}d^{3+2/\ell}}{\varepsilon^{4+4/\ell}}$	$\dfrac{d^{\alpha}}{\varepsilon^{3}}$	$\dfrac{n_{L}^{2}d^{\alpha}}{\varepsilon^{5}}$
Nonlinear Hamilton-Jacobi PDEs^♯
Problem	Quantum simulation		Spectral QLSA		FD QLSA ( $\alpha\geq 4$ )
Problem	Subroutine	Observable	Subroutine	Observable	Subroutine	Observable
Liouville equation	$\dfrac{d}{\varepsilon^{2}}$	$\dfrac{n_{H}^{4}d}{\varepsilon^{4}}$	$\dfrac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}$	$\dfrac{n_{H}^{4}d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}$	$\dfrac{d^{\alpha}}{\varepsilon^{3}}$	$\dfrac{n_{H}^{4}d^{\alpha}}{\varepsilon^{5}}$
Schrödinger equation	$\dfrac{d}{\varepsilon}$	$\dfrac{c_{O}N_{u_{0}}^{4}d}{\varepsilon^{3}}$	$\dfrac{d^{2+2/\ell}}{\varepsilon^{1+4/\ell}}$	$\dfrac{c_{O}N_{u_{0}}^{4}d^{2+2/\ell}}{\varepsilon^{3+4/\ell}}$

a)

The notations $\mathcal{O}$ for quantum simulations and $\widetilde{\mathcal{O}}$ for QLSA based methods are omitted. Since the dependence on the matrix order is not explicitly presented in [20], we just include the multiplicative factor $d^{\alpha-3}~{}(\alpha\geq 4)$ in the gate complexity for the finite difference discretisations (note that the order of the matrix grows exponentially with respect to the dimension).
b)

$n_{L}$ and $n_{H}$ are the sampling factors for the Liouville equation, given respectively in Theorems 3.2 and 4.2. $N_{u_{0}}$ and $c_{O}$ are defined by (5.13) and (5.15), respectively. Note that $c_{O}$ depends on $d,\varepsilon$ when computing the current density and kinetic energy. For the Hamilton-Jacobi equations, $c_{O}$ may be neglected, and $N_{u_{0}}\lesssim n_{H}^{1/2}$ (see Remark 5.2).
*)

Despite the absence of the unitary structure, we still proposed a “quantum simulation” algorithm for the Liouville representation in Appendix A, where non-unitary procedures are involved. The results are presented in Theorem A.1 and Theorem A.2 when the cost arising from multiple copies of initial quantum states is ignored.
$\sharp$ )

For the spectral discretisation of the Liouville/Schrödinger equation, we only consider a prototype case $H(x,p)=\frac{1}{2}|p|^{2}+V(x)$ .

Let us summarise the several parts of the article as follows.

(1) Motivated by the idea in [23, 38], we established the correspondence between the nonlinear dynamic system and the Liouville representation via a simple ansatz, which relates the ODE solution and the density distribution by the Dirac delta function. In this case, the ODE solution can be recast as a physical observable of the Liouville equation, which provides an efficient way to solve the nonlinear ODEs by using quantum algorithms for the resulting linear Liouville equation. We introduced the Liouville approximation and the KvN approximation from the perspective of quantum differential equations solvers for the Liouville equation with smoothed initial data, while the KvN approximation can be regarded as the “symmetrized” counterpart of the Liouville approximation (See Fig. 3 for illustration). For both linear representation approaches, we proposed the upwind difference discretisations and Fourier spectral discretisations and provided in detail the time complexity analysis for the QLSA based methods and the quantum simulation methods, including the output of the quantum states and the post-processing for the computation of physical observables. The KvN mechanism allows direct quantum Hamiltonian simulations, while the Liouville method can also be translated into the evolution of symmetric propagators with additional non-unitary procedures at every time step after an appropriate Trotter based approximation.

(2) For the nonlinear PDEs, more specifically the Hamilton-Jacobi equations, by using the level set mechanism as proposed in [38], one can map the nonlinear PDEs of $(d+1)$ -dimension to a linear $(2d+1)$ -dimensional Liouville equation, referred to as the Liouville representation for nonlinear PDEs. For a classical device, doubling the dimension of the problem may seem too costly because the cost increases exponentially with dimension. However, for quantum algorithms, the relative overhead in doubling the dimension can be up to exponentially smaller, which has been verified by the proposed quantum algorithms since no exponential terms in dimension like $d^{d}$ and $(1/\varepsilon)^{d}$ for the classical cost are included.

(3) It is well-known that the Schrödinger equation can be transformed into the quantum Liouville equation via the Wigner transform, which in turn leads to the Liouville equation when taking the semiclassical limit. In view of the close relations between their physical observables, we introduced the Schrödinger framework for solving the Liouville equation. We studied the quantum interpretation of the classical time-splitting Fourier spectral method proposed in [4] for the Schrödinger equation, and presented a comprehensive discussion in the correspondence between the time-splitting spectral method and the Trotter based Hamiltonian simulation, although this issue has been addressed (in less detail) in some earlier literature. The time-splitting spectral discretisation for the Schrödinger equation generates a discretised Hamiltonian system, which can be handled by standard Hamiltonian simulation algorithms or quantum linear systems algorithms. We analysed in detail the gate complexity of these two approaches for numerically resolving the wave function and the physical observables. Despite the advantages in terms of time complexity, it should be pointed out that the solution to the Schrödiner equatin is oscillatory, hence if one wants high-resolution (oscillation-free) numerical results the Liouville framework is preferred.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

SJ was partially supported by the NSFC grant No. 12031013, the Shanghai Municipal Science and Technology Major Project (2021SHZDZX0102), and the Innovation Program of Shanghai Municipal Education Commission (No. 2021-01-07-00-02-E00087). NL acknowledges funding from the Science and Technology Program of Shanghai, China (21JC1402900), the Shanghai Pujiang Talent Grant (no. 20PJ1408400) and the NSFC International Young Scientists Project (no. 12050410230). YY was partially supported by China Postdoctoral Science Foundation (no. 2022M712080).

Appendix A Spectral discretisation for the Liouville representation of nonlinear ODEs

A.1 The Trotter based spectral discretisation

Consider the problem (3.23). Let $u=F_{i}(x)w(t,x)$ . According to the above discussion, one has

-\mathrm{i}\frac{\partial}{\partial x_{i}}u\longrightarrow\hat{P}_{i}^{\rm d}\bm{u}=\bm{P}_{i}\bm{u}=\bm{P}_{i}\Lambda_{F_{i}}\bm{w},

where $\Lambda_{F_{i}}=\text{diag}(\bm{F}_{i})$ is a diagonal matrix, where the vector $\bm{F}_{i}=\sum_{\bm{j}}F_{i}(x_{\bm{j}})|\bm{j}\rangle$ . The resulting system of ordinary differential equations is

\begin{cases}\frac{\mathrm{d}}{\mathrm{d}t}\bm{w}(t)=-\mathrm{i}A\bm{w}(t),\\ \bm{w}(0)=\bm{w}^{0}=(w_{0}(x_{\bm{j}})),\end{cases}

(A.1)

where $A=\sum_{i=1}^{d}A_{i}$ with $A_{i}=\bm{P}_{i}\Lambda_{F_{i}}$ . One can check from (3.38) that $\bm{P}_{i}$ are Hermitian matrices. However, we note that this does not mean that $A$ is Hermitian, since in general $(\bm{P}_{i}\Lambda_{F_{i}})^{\dagger}=\Lambda_{F_{i}}\bm{P}_{i}\neq\bm{P}_{i}\Lambda_{F_{i}}$ .

The evolution of (A.1) can be formally written as

|\psi(t+\Delta t)\rangle=\mathrm{e}^{-\mathrm{i}(A_{1}+\cdots+A_{d})\Delta t}|\psi(t)\rangle,

where $|\psi(t)\rangle$ is a quantum state whose amplitudes are proportional to $\bm{w}(t)$ and the evolution operator $\exp(-\mathrm{i}A_{i}\Delta t)$ is not necessarily unitary. Let us consider the first-order product formula

U_{\Delta t}=\mathrm{e}^{-\mathrm{i}A_{d}\Delta t}\cdots\mathrm{e}^{-\mathrm{i}A_{1}\Delta t}.

(A.2)

One obtains from [62, 16] that

\mathrm{e}^{-\mathrm{i}(A_{1}+\cdots+A_{d})\Delta t}=U_{\Delta t}+C_{A}\Delta t^{2},

(A.3)

where $C_{A}$ depends on the matrix $A$ , considered as $\mathcal{O}(1)$ in the following. Therefore, the problem is reduced to the simulation of each $A_{j}$ , where $A_{j}$ is not necessarily symmetric. Take $j=1$ as an example and consider the decomposition $\Lambda_{F_{1}}=\Lambda_{F_{1}}^{+}-\Lambda_{F_{1}}^{-}$ , where $\Lambda_{F_{1}}^{\pm}=\text{diag}(d_{1}^{\pm},\cdots,d_{n}^{\pm})$ are diagonal matrices with $d_{j}^{\pm}>\alpha$ for some positive constant $\alpha$ . One can further require that $\|\Lambda_{F_{i}}^{\pm}\|\lesssim\max_{1\leq j\leq d}\|\bm{F}_{j}\|$ . Let $A_{1}^{\pm}=\bm{P}_{1}\Lambda_{F_{1}}^{\pm}$ . The Strang splitting gives

\mathrm{e}^{-\mathrm{i}A_{1}\Delta t}=\mathrm{e}^{-\mathrm{i}A_{1}^{+}\Delta t/2}\mathrm{e}^{\mathrm{i}A_{1}^{-}\Delta t}\mathrm{e}^{-\mathrm{i}A_{1}^{+}\Delta t/2}+C_{1}\Delta t^{3},

(A.4)

where $\mathrm{e}^{\mathrm{i}A_{1}^{-}\Delta t}$ and $\mathrm{e}^{-\mathrm{i}A_{1}^{+}\Delta t/2}$ can be evolved in the similar way (note that one can also use the first-order approximation). To this end, we consider the second one as an example.

The simulation of $\mathrm{e}^{-\mathrm{i}A_{1}^{+}\Delta t/2}$ is related to the following ODEs:

\frac{\mathrm{d}}{\mathrm{d}t}\bm{w}(t)=-\frac{\mathrm{i}}{2}A_{1}^{+}\bm{w}(t)=-\frac{\mathrm{i}}{2}\bm{P}_{1}\Lambda_{F_{1}}^{+}\bm{w}(t),\qquad 0\leq t\leq\Delta t.

Let $\widetilde{\bm{w}}=\sqrt{\Lambda_{F_{1}}^{+}}\bm{w}$ , where $\sqrt{\Lambda_{F_{1}}^{+}}=\text{diag}(\sqrt{d_{1}^{+}},\cdots,\sqrt{d_{n}^{+}})$ . Then the above ODEs can be reformulated as

\frac{\mathrm{d}}{\mathrm{d}t}\widetilde{\bm{w}}(t)=-\frac{\mathrm{i}}{2}\widetilde{A}_{1}^{+}\widetilde{\bm{w}}(t),\qquad 0\leq t\leq\Delta t,

(A.5)

where $\widetilde{A}_{1}^{+}=\sqrt{\Lambda_{F_{1}}^{+}}\bm{P}_{1}\sqrt{\Lambda_{F_{1}}^{+}}$ is a Hermitian matrix. The one-step simulation gives

|\widetilde{\psi}(t+\Delta t)\rangle=\mathrm{e}^{-\mathrm{i}\widetilde{A}_{1}^{+}\Delta t/2}|\widetilde{\psi}(t)\rangle,

(A.6)

where $|\widetilde{\psi}(t)\rangle$ corresponds to $\widetilde{\bm{w}}(t)$ . Since $d_{j}^{\pm}>\alpha>0$ , $\widetilde{\bm{w}}=\sqrt{\Lambda_{F_{1}}^{+}}\bm{w}$ can be viewed as a linear systems problem:

\widetilde{\bm{w}}=D_{1}^{-1}\bm{w},\qquad D_{1}=(\sqrt{\Lambda_{F_{1}}^{+}})^{-1}.

(A.7)

Similarly,

\bm{w}=D_{2}^{-1}\widetilde{\bm{w}},\qquad D_{2}=\sqrt{\Lambda_{F_{1}}^{+}}=D_{1}^{-1}.

(A.8)

Given the initial state of $\bm{w}^{0}$ , denoted by $|\psi^{0}\rangle$ . At each time step, one needs to consider the procedure

|\psi^{0}\rangle\xrightarrow{\eqref{w2tw}}|\widetilde{\psi}^{0}\rangle\xrightarrow{\eqref{tpsi1}}|\widetilde{\psi}^{1}\rangle\xrightarrow{\eqref{tw2w}}|\psi^{1}\rangle

for $\mathrm{e}^{-\mathrm{i}A_{1}^{+}\Delta t/2}$ , followed by the similar procedures for $\mathrm{e}^{\mathrm{i}A_{1}^{-}\Delta t}$ and the first $\mathrm{e}^{-\mathrm{i}A_{1}^{+}\Delta t/2}$ in (A.4), where (A.6) can be solved by quantum Hamiltonian simulations or quantum differential equations solvers.

Remark A.1.

The transition between $\widetilde{\bm{w}}$ and $\bm{w}$ in $\eqref{w2tw}$ and (A.8) may be implemented in a simpler manner, for example, the LCU method, which decomposes the diagonal (and Hermitian) matrices $D_{i}=\text{diag}(d_{i,1},\cdots,d_{i,n})$ , $i=1,2$ into a sum of two unitary operations. In fact, it is always possible to write $D_{i}=(U_{i}+V_{i})/2$ , where $U_{i}=D_{i}+\mathrm{i}\sqrt{\mathbf{1}-D^{2}_{i}}$ , and $V_{i}=D_{i}-\mathrm{i}\sqrt{\mathbf{1}-D^{2}_{i}}$ (One can assume $\|D_{i}\|\leq 1$ after an adjustment). These unitaries are also diagonal matrices with diagonal entries $\text{diag}(d_{i,1}\pm\mathrm{i}\sqrt{1-d^{2}_{i,1}},\cdots,d_{i,n}\pm\mathrm{i}\sqrt{1-d^{2}_{i,n}})$ . Applying the operation $D_{i}$ onto a quantum state can be done by a straightforward application of the LCU method (e.g. Lemma 6 in [14]). Here we can assume access to the control operation $|0\rangle\langle 0|\otimes U_{i}+|1\rangle\langle 1|\otimes V_{i}$ .

However, multiple copies are needed at every time step for both the QLSA and the LCU since they are not unitary procedures, hence the cost (i.e. number of copies needed of the initial state) will increase exponentially with $N_{t}$ . We can see the last statement more explicitly.

•

For the QLSA, when solving $D_{1}|\widetilde{w}\rangle=|w\rangle$ with the quantum state $|w\rangle$ given, one must prepare unitary operations to query the entries of $w$ . This needs post-processing or multiple uses of $|w\rangle$ , hence multiple copies of $|w\rangle$ .
•

For the LCU, to obtain the state $D_{i}|w\rangle$ for some state $|w\rangle$ , one can construct a unitary procedure acting on $|w\rangle$ and an ancilla that will output a state $|w^{\prime}\rangle=\alpha_{i}D_{i}|w\rangle+\beta_{i}|v\rangle$ , where $|v\rangle$ is some state we don’t want. One can then obtain a single copy of $D_{i}|w\rangle$ upon post-selection of $|w^{\prime}\rangle$ with $\mathcal{O}(1/\|\alpha_{i}D_{i}|w\rangle\|^{2})$ number of measurements, and hence multiple copies of $|w\rangle$ , where $\|\alpha_{i}D_{i}|w\rangle\|\leq 1$ .

We are not only interested in the cost at every time step: we want to evolve for a long time, to $N_{t}$ time-steps. Suppose we wanted to repeat this procedure $N_{t}$ times with at least $C>1$ copies at every time step. At the last $N_{t}^{\text{th}}$ step if we need only a single copy of the desired state, then $C$ copies of the state in the previous $N_{t}-1$ time-step is required, which means $C^{2}$ copies of the state in the $(N_{t}-2)^{\text{th}}$ time-step is needed. Ultimately, this requires $C^{N_{t}}$ copies of the original $|w\rangle$ state in the initial time-step.

In the following, we ignore the cost of multiple copies at each time step.

A.1.1 The QLSA for the spectral discretisation

Theorem A.1.

Assume further that $\max_{1\leq j\leq d}\|\bm{F}_{j}\|=\mathcal{O}(1)$ and $T=\mathcal{O}(1)$ .

(1)

There exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{w}(T)/\|\bm{w}(T)\|$ with the gate complexity given by

$N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)}.$
(2)

The observable of the Liouville representation can be computed with gate complexity given by

$N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{L}^{4}d^{3+2/\ell}}{\varepsilon^{6+4/\ell}}\Big{)},$

where $n_{L}=\|(\bm{\rho}^{\omega})^{0}\|/M^{d/2}$ .

Proof.

For simplicity, we omit the discussion of the cost and the error resulting from (A.7) and (A.8).

(1) Since $\widetilde{A}_{1}^{+}$ is hermitian, the eigenvalues of $\widetilde{A}_{1}^{+}$ are real and $\kappa_{V}=1$ , where $V$ is the transformation matrix associated with $\widetilde{A}_{1}^{+}$ . At each time step, by Lemma 3.4, the gate complexity of solving (A.5) within error $\eta$ is

Q_{\Delta t}=\widetilde{\mathcal{O}}(s\kappa_{V}\|\widetilde{A}_{1}^{+}\|)=\widetilde{\mathcal{O}}(M^{2})=\mathcal{O}(M^{2}\text{Polylog}(M^{d+2}/\eta)),

where $s=\mathcal{O}(M)=\mathcal{O}(1/\Delta x)$ , $\|\widetilde{A}_{1}^{+}\|=\mathcal{O}(M)$ . Therefore, $U_{\Delta t}$ defined in (A.2) can be evolved within error $\mathcal{O}(d(3\eta+\Delta t^{2}))$ [51, Proposition 1.12] with

N_{\text{Gates}}(U_{\Delta t})=\mathcal{O}(dQ_{\Delta t})=\widetilde{\mathcal{O}}(dM^{2}).

(A.9)

In the following, we choose $\eta\sim\Delta t^{2}$ , and hence $U_{\Delta t}$ can be evolved within error $\mathcal{O}(d\Delta t^{2})$ . According to the error estimate (3.17) and noting Eq. (A.3), one may have

e_{\rho}\leq C(\omega+d\Delta t/\omega+d\Delta x^{\ell}/\omega^{\ell+1}),

which is also true for the computation of the observable. The above error bounds suggest the following mesh strategy:

\omega\sim d\Delta t/\omega\sim d\Delta x^{\ell}/\omega^{\ell+1}\sim\varepsilon,

(A.10)

or equivalently,

M\sim d^{1/\ell}/\varepsilon^{1+2/\ell},\qquad\Delta t=\varepsilon^{2}/d.

(A.11)

The total number of gates required to iterate to the $n$ -th step is

\displaystyle N_{\text{Gates}}=nN_{\text{Gates}}(U_{\Delta t})=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{3+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)},

where in the last equal sign we have included the additional factor $d$ arising from the matrix order (see Remark 3.8).

(2) For the computation of the observable, according to Remark 3.4, one just needs to multiply the original gate complexity by the sampling factor $k=\mathcal{O}(n_{L}^{4}/\varepsilon^{2})$ . This completes the proof. ∎

A.1.2 Quantum simulation for the spectral discretisation

The approximate operator in (A.2) can also evolved by the quantum simulation.

Theorem A.2.

Assume further that $\max_{1\leq j\leq d}\|\bm{F}_{j}\|=\mathcal{O}(1)$ and $T=\mathcal{O}(1)$ .

(1)

There exists a quantum algorithm that produces a state $\varepsilon$ -close to $\bm{w}(T)/\|\bm{w}(T)\|$ with the gate complexity given by

$N_{Gates}=\widetilde{\mathcal{O}}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\Big{)}.$
(2)

The observable of the Liouville representation can be simulated with gate complexity given by

$N_{\text{Gates}}(\langle O\rangle)=\widetilde{\mathcal{O}}\Big{(}\frac{n_{L}^{4}d^{2+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)},$

where $n_{L}=\|(\bm{\rho}^{\omega})^{0}\|/M^{d/2}$ .

Proof.

We first quantify the number of gates used in the quantum simulation. According to Theorem 1 in [6], $\mathrm{e}^{-\mathrm{i}\widetilde{A}_{j}^{+}\Delta t/2}$ in (A.6) can be simulated within error $\eta$ with

\mathcal{O}\Big{(}\tau(m_{d}+\log^{2.5}(\tau/\eta))\frac{\log(\tau/\eta)}{\log\log(\tau/\eta)}\Big{)}

2-qubits gates, where $\tau=s\|\widetilde{A}_{j}^{+}\|_{\max}\Delta t/2$ , $s$ is the sparsity of $\widetilde{A}_{j}^{+}$ and $\|\widetilde{A}_{j}^{+}\|_{\max}$ denotes the largest entry of $\widetilde{A}_{j}^{+}$ in absolute value. This result is near-optimal by Theorem 2 therein. One can check that the sparsity of $\widetilde{A}_{j}^{+}$ is $s=\mathcal{O}(M)$ . Therefore, $U_{\Delta t}$ defined in (A.2) can be simulated within error $\mathcal{O}(d(3\eta+\Delta t^{2}))$ [51, Proposition 1.12] with

N_{\text{Gates}}(U_{\Delta t})=\mathcal{O}(d\tau m_{d}\cdot\text{polylog}),

where

\tau=M\widetilde{A}_{\max}\Delta t,\qquad\widetilde{A}_{\max}=\max_{j}\|\widetilde{A}_{j}^{\pm}\|_{\max}=\mathcal{O}(M),

\text{polylog}=\log^{2.5}(\tau/\eta)\frac{\log(\tau/\eta)}{\log\log(\tau/\eta)}.

In the following, we choose $\eta\sim\Delta t^{2}$ , and hence obtain the same mesh strategy for the QLSA.

With the mesh strategy given in (A.11), one has the number of qubits per dimension is

m=\mathcal{O}(\log M)=\mathcal{O}\Big{(}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\Big{)}.

The total number of qubits is $m_{d}=dm$ . With these settings, we obtain

\tau=\mathcal{O}\Big{(}\frac{d^{2/\ell}}{\varepsilon^{2+4/\ell}}\Delta t\Big{)},\qquad\tau/\eta=\mathcal{O}\Big{(}\frac{d^{1+2/\ell}}{\varepsilon^{4+4/\ell}}\Big{)},

and the total number of gates required to iterate to the $n$ -th step is

\displaystyle N_{\text{Gates}}=nN_{\text{Gates}}(U_{\Delta t})=\mathcal{O}\Big{(}\frac{d^{2+2/\ell}}{\varepsilon^{2+4/\ell}}\log\frac{d^{1/\ell}}{\varepsilon^{1+2/\ell}}\cdot\text{polylog}\Big{)}.

(2) For the computation of the observable, one just needs to add the multiplicative factor $\mathcal{O}(n_{L}^{4}/\varepsilon^{2})$ . This completes the proof. ∎

Remark A.2.

We emphasise that here only the total cost of the Hamiltonian simulation component at each time-step is included. The simulation protocol here is different from the traditional time-marching Hamiltonian simulation since non-unitary procedures are involved at each time step, leading to exponential increase of the cost as discussed in Remark A.1. It is worth pointing out that this obstacle can be overcome by our recently proposed technique - the Schrödingerisation approach [39, 40]. Please refer to Section 4.5 of [40] for details.

References

[1] D. Aharonov, V. Jones, and Z. Landau. A polynomial quantum algorithm for approximating the jones polynomial. Algorithmica, 55(3):395–421, 2009.
[2] K. Alexander and R. Jeffrey. The Order of Accuracy of Quadrature Formulae for Periodic Functions, pages 155–159. Birkhäuser Boston, Boston, 2009.
[3] A. Ambainis. Variable time amplitude amplification and quantum algorithms for linear algebra problems. LIPIcs. Leibniz Int. Proc. Inform., 14:636–647, 2012.
[4] W. Bao, S. Jin, and P. A. Markowich. On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime. J. Comput. Phys., 175:487–524, 2002.
[5] D. W. Berry. High-order quantum algorithm for solving linear differential equations. J. Phys. A: Math. Theor., 47(10):105301, 17 pp., 2014.
[6] D. W. Berry, A. M. Childs, and R. Kothari. Hamiltonian simulation with nearly optimal dependence on all parameters. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 792–809, 2015.
[7] D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang. Quantum algorithm for linear differential equations with exponentially improved dependence on precision. Comm. Math. Phys., 356(3):1057–1081, 2017.
[8] A. V. Bobylëv. On the Chapman-Enskog and Grad methods for solving the Boltzmann equation. Dokl. Akad. Nauk SSSR, 262(1):71–75, 1982.
[9] Y. I. Bogdanov and N. A. Bogdanova. The study of classical dynamical systems using quantum theory. In Proc. SPIE, International Conference on Micro- and Nano-Electronics, volume 9440, 2014.
[10] Y. I. Bogdanov, N. A. Bogdanova, D. V. Fastovets, and V. F. Lukichev. Quantum approach to the dynamical systems modeling. In Proc. SPIE, International Conference on Micro- and Nano-Electronics, volume 11022, 2019.
[11] Y. Cao, A. Papageorgiou, I. Petras, et al. Quantum algorithm and circuit design solving the Poisson equation. New J. Phys., 15:013021, 2013.
[12] L. Caroline and L. Christian. Computing quantum dynamics in the semiclassical regime. Acta Numer., 29:229–401, 2020.
[13] S. Chakraborty, A. Gilyén, and S. Jeffery. The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation. In 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), volume 132 of Leibniz International Proceedings in Informatics (LIPIcs), pages 33:1–33:14, 2019.
[14] A. M. Childs, R. Kothari, and R. D. Somma. Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. SIAM J. Comput., 46(6):1920–1950, 2017.
[15] A. M. Childs, J. P. Liu, and A. Ostrander. High-precision quantum algorithms for partial differential equations. Quantum, 5:574, 2021.
[16] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu. Theory of trotter error with commutator scaling. Phys. Rev. X, 11(1):011020, Feb 2021.
[17] A. W. Childs, T. Li, J. Liu, C. Wang, and R. Zhang. Quantum algorithms for sampling log-concave distributions and estimating normalizing constants. arXiv:2210.06539, 2022.
[18] A. W. Childs and J. Liu. Quantum spectral methods for differential equations. Comm. Math. Phys., 375(2):1427–1457, 2020.
[19] B. D. Clader, B. C. Jacobs, and C. R. Sprouse. Preconditioned quantum linear system algorithm. Phys. Rev. Lett., 110:250504, Jun 2013.
[20] P. C. S. Costa, D. An, Y. A. Sanders, Y. Su, R. Babbush, and D. W. Berry. Optimal scaling quantum linear systems solver via discrete adiabatic theorem. arXiv:2111.08152, 2021.
[21] P. C. S. Costa, S. Jordan, and A. Ostrander. Quantum algorithm for simulating the wave equation. Phys. Rev. A, 99:012323, 22 pp., 2019.
[22] D. Deutsch and R. Jozsa. Rapid solution of problems by quantum computation. Proc. R. Soc. London, Ser. A, 439:553, 1992.
[23] I. Y. Dodin and E. A. Startsev. On applications of quantum computing to plasma simulations. Phys. Plasmas, 28:092101, 2021.
[24] A. Dong, J. Liu, D. Wang, and Q. Zhao. A theory of quantum differential equation solvers: limitations and fast-forwarding. arXiv:2211.05246, 2022.
[25] Y. Dong, L. Lin, and Y. Tong. Ground-state preparation and energy estimation on early fault-tolerant quantum computers via quantum eigenvalue transformation of unitary matrices. arXiv:2204.05955, 2022.
[26] A. Engel, G. Smith, and S. E. Parker. Quantum algorithm for the Vlasov equation. Phys. Rev. A, 100:062315, Dec 2019.
[27] A. Engel, G. Smith, and S. E. Parker. Linear embedding of nonlinear dynamical systems and prospects for efficient quantum algorithms. Phys. Plasmas, 28:062305, 2021.
[28] Di Fang, Lin Lin, and Yu Tong. Time-marching based quantum solvers for time-dependent linear differential equations. arXiv preprint arXiv:2208.06941, 2022.
[29] P. Gérard, P. A. Markowich, N. J. Mauser, and F. Poupaud. Homogenization limits and Wigner transforms. Comm. Pure Appl. Math., 50(4):323–379, 1997.
[30] A. Gilyén, Y. Su, G. Low, and N. Wiebe. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204. STOC, 2019.
[31] F. Golse, S. Jin, and T. Paul. On the convergence of time splitting methods for quantum dynamics in the semiclassical regime. Found. Comput. Math., 21(3):613–647, 2021.
[32] Harold Grad. Asymptotic theory of the Boltzmann equation. Phys. Fluids, 6(2):147–181, 1963.
[33] A. W. Harrow, A. Hassidim, and S. Lloyd. Quantum algorithm for linear systems of equations. Phys. Rev. Lett., 103(15):150502, 4 pp., 2009.
[34] A. W. Harrow and A. Y. Wei. Adaptive quantum simulated annealing for Bayesian inference and estimating partition functions. In Proceedings of the 2020 ACM-SIAM Symposium on Discrete Algorithms, pages 193–212. SIAM, Philadelphia, PA, 2020.
[35] R. Horn and C. R. Johnson. Matrix Analysis. Second edition. Cambridge University Press, Cambridge, 2013.
[36] S. Jin. Asymptotic-preserving schemes for multiscale physical problems. Acta Numer., 31:415–489, 2022.
[37] S. Jin, X. Li, and N. Liu. Quantum simulation in the semi-classical regime. Quantum, 6:739, 2022.
[38] S. Jin and N. Liu. Quantum algorithms for computing observables of nonlinear partial differential equations. arXiv:2202.07834, 2022.
[39] S. Jin, N. Liu, and Y. Yu. Quantum simulation of partial differential equations via Schrödingerisation. arXiv:2212.13969, 2022.
[40] S. Jin, N. Liu, and Y. Yu. Quantum simulation of partial differential equations via Schrödingerisation: technical details. arXiv:2212.14703, 2022.
[41] S. Jin, N. Liu, and Y. Yu. Time complexity analysis of quantum difference methods for linear high dimensional and multiscale partial differential equations. J. Comput. Phys., 471:111641, 2022.
[42] S. Jin, P. Markowich, and C. Sparber. Mathematical and computational methods for semiclassical Schrödinger equations. Acta Numer., 20:121–209, 2011.
[43] S. Jin and S. Osher. A level set method for the computation of multivalued solutions to quasi-linear hyperbolic PDEs and Hamilton-Jacobi equations. Commun. Math. Sci., 1(3):575–591, 2003.
[44] I. Joseph. Koopman–von Neumann approach to quantum simulation of nonlinear classical dynamics. Phys. Rev. Research, 2(4):043102, 2020.
[45] I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni, and A. Aspuru-Guzik. Polynomial-time quantum algorithm for the simulation of chemical dynamics. Proceedings of the National Academy of Sciences, 105(48):18681–18686, 2008.
[46] H. Krovi. Improved quantum algorithms for linear and nonlinear differential equations. arXiv:2202.01054, 2022.
[47] C. Lasser and C. Lubich. Computing quantum dynamics in the semiclassical regime. Acta Numer., pages 229–401, 2020.
[48] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge, 2002.
[49] C. D. Levermore. Moment closure hierarchies for kinetic theories. J. Statist. Phys., 83(5-6):1021–1065, 1996.
[50] S. K. Leyton and T. J. Osborne. A quantum algorithm to solve nonlinear differential equations. arXiv:0812.4423, 2008.
[51] L. Lin. Lecture notes on quantum algorithms for scientific computation. arXiv:2201.08309, 2022.
[52] L. Lin and Y. Tong. Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems. Quantum Inf. Process., 4:361, 2020.
[53] L. Lin and Y. Tong. Heisenberg-limited ground state energy estimation for early fault-tolerant quantum computers. arXiv:2102.11340v2, 2022.
[54] Y. Lin, R. B. Lowrie, D. Aslangil, Y. Subaşi, and A. T. Sornborger. Koopman von Neumann mechanics and the Koopman representation: A perspective on solving nonlinear dynamical systems with quantum computers. arXiv:2202.02188v2, 2022.
[55] N. Linden, A. Montanaro, and C. Shao. Quantum vs. classical algorithms for solving the heat equation. arXiv:2004.06516, 2020.
[56] P.-L. Lions and T. Paul. Sur les mesures de Wigner. Rev. Mat. Iberoamericana, 9(3):553–618, 1993.
[57] R. J. Lipton and K. W. Regan. Quantum Algorithm via Linear Algebra. Cambridge, Massachusetts, 2010.
[58] J. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs. Efficient quantum algorithm for dissipative nonlinear differential equations. Proc. Natl. Acad. Sci. U. S. A., 118(35), 2021.
[59] N. Liu, T. F. Demarie, S. Tan, L. Aolita, and J. F. Fitzsimons. Client-friendly continuous-variable blind and verifiable quantum computing. Phys. Rev. A, 100:062309, 10 pp., 2019.
[60] N. Liu, J. Thompson, C. Weedbrook, S. Lloyd, V. Vedral, M. Gu, and K. Modi. Power of one qumode for quantum computation. Phys. Rev. A, 93:052304, 10 pp., 2016.
[61] S. Lloyd, G. D. Palma, C. Gokler, B. Kiani, Z. Liu, M. Marvian, F. Tennie, and T. Palmer. Quantum algorithm for nonlinear differential equations. arXiv:2011.06571, 2020.
[62] M. A. Nielsen and I. L. Chuang. Quantum Computation and Quantum Information. Cambridge, New York, 2010.
[63] L. Qi. Some simple estimates for singular values of a matrix. Linear Algebra Appl., 56:105–119, 1984.
[64] P. A. Raviart. An analysis of particle methods. In Numerical Methods in Fluid Dynamics, pages 243–324. Springer, Berlin, 1985.
[65] P. W. Shor. Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer. SIAM J. Comput., 26:1484, 1997.
[66] Y. Subasi and R. D. Somma. Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing. Phys. Rev. Lett., 122:060504, 2019.
[67] C. Xue, Y. Wu, and G. Guo. Quantum homotopy perturbation method for nonlinear dissipative ordinary differential equations. New J. Phys., 23:No. 123035, 20 pp, 2021.

	$\displaystyle I_{1}$	$\displaystyle=\|\langle G_{\rho}(t_{n})\rangle-\langle G_{\rho^{\omega}}(t_{n})\rangle\|=\Big{\|}\int_{[0,1]^{d}}G(x)(\rho(t_{n},x)-\rho^{\omega}(t_{n},x))\mathrm{d}x\Big{\|}$
		$\displaystyle=\Big{\|}\int_{[0,1]^{d}}G(x)\Big{(}\rho_{0}(X(0;x,t_{n}))-\rho_{0}^{\omega}(X(0;x,t_{n}))\Big{)}J(0;x,t_{n})\mathrm{d}x\Big{\|}$
		$\displaystyle\leq C\omega\max_{x\in[0,1]^{d}}J(0;x,t_{n})\leq C\omega\mathrm{e}^{t_{n}\\|\text{div}F\\|_{\infty}},$

	$\displaystyle e_{\psi}$	$\displaystyle=\|\langle G(t_{n})\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|$
		$\displaystyle\leq\|\langle G_{\rho}(t_{n})\rangle-\langle O_{\psi^{\omega}}(t_{n})\rangle\|+\|\langle O_{\psi^{\omega}}(t_{n})\rangle-\langle O_{\psi^{\omega},n}\rangle\|+\|\langle O_{\psi^{\omega},n}\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|$
		$\displaystyle=\|\langle G_{\rho}(t_{n})\rangle-\langle G_{\rho^{\omega}}(t_{n})\rangle\|+\|\langle O_{\psi^{\omega}}(t_{n})\rangle-\langle O_{\psi^{\omega},n}\rangle\|+\|\langle O_{\psi^{\omega},n}\rangle-\langle O_{\psi_{h}^{\omega},n}\rangle\|$
		$\displaystyle=:I_{1}+I_{2}+I_{3}.$

	$\displaystyle\\|\bm{w}^{n+1}\\|_{1}$	$\displaystyle=\sum\limits_{j}\|w_{j}^{n+1}\|\leq\sum\limits_{j}\Big{[}(1-\lambda(a_{j}^{+}-a_{j-1}^{-}))\|w^{n}_{j}\|-\lambda a_{j}^{-}\|w^{n}_{j+1}\|+\lambda a_{j-1}^{+}\|w^{n}_{j-1}\|\Big{]}$
		$\displaystyle=\sum\limits_{j}\Big{[}(1-\lambda(a_{j}^{+}-a_{j-1}^{-}))\|w^{n}_{j}\|-\lambda a_{j-1}^{-}\|w^{n}_{j}\|+\lambda a_{j}^{+}\|w^{n}_{j}\|\Big{]}$
		$\displaystyle=\sum\limits_{j}\|w_{j}^{n}\|=\\|\bm{w}^{n}\\|_{1},$

Time complexity analysis of quantum algorithms via linear representations for nonlinear ordinary and partial differential equations

Abstract

1 Introduction

2 Review of quantum nonlinear differential equation solvers

3 Linear representation methods for nonlinear ODEs: Liouville vs. KvN representations

3.1 Linear representation methods

3.1.1 The Liouville representation

3.1.2 The Koopman-von Neumann representation

3.1.3 Computation of the quantity of interest and error analysis

Lemma 3.1 (Error of the Liouville representation).

Proof.

Remark 3.1.

Lemma 3.2 (Error of KvN representation).

Proof.

Remark 3.2.

3.2 Finite difference discretisation for the Liouville representation

3.2.1 The QLSA for the finite difference discretisation

Theorem 3.1.

Proof.

Remark 3.3.

3.2.2 The algorithm for the computation of the observable

The expectation of the observable

The general sampling law

Boosting the success probability of the final state projection

3.2.3 The gate complexity of computing the observable

Lemma 3.3 (Estimation of ‖A−1​b‖\|A^{-1}b\|).

Theorem 3.2.

Proof.

Remark 3.4.

Remark 3.5.

3.3 Finite difference discretisation for the KvN representation

3.3.1 The QLSA for the finite difference discretisation

Theorem 3.3.

Proof.

Remark 3.6.

3.3.2 The gate complexity of computing the observable

Theorem 3.4.

Proof.

Remark 3.7.

3.4 Spectral discretisation for the KvN representation

3.4.1 The notations

3.4.2 The QLSA for the spectral discretisation

Lemma 3.4.

Remark 3.8.

Theorem 3.5.

Proof.

3.4.3 The quantum simulation for the spectral discretisation

Theorem 3.6.

Proof.

Remark 3.9.

Remark 3.10.

4 Liouville representation for nonlinear Hamilton-Jacobi PDEs: finite difference vs. spectral approximations

4.1 Finite difference discretisation for the Liouville equation

4.1.1 The QLSA for the finite difference discretisation

Theorem 4.1.

Proof.

4.1.2 The computation of the physical observables

Theorem 4.2.

Proof.

Remark 4.1.

4.2 Spectral discretisation for the Liouville equation

4.2.1 The QLSA for the spectral discretisation

Theorem 4.3.

Proof.

4.2.2 The quantum simulation for the spectral discretisation

Theorem 4.4.

Proof.

Remark 4.2.

5 The Schrödinger framework

5.1 The semiclassical approximation

5.2 Quantum simulations for the spectral discretisation

5.2.1 The time-splitting spectral approximations

Remark 5.1.

5.2.2 The quantum simulation of the Schrödinger equation

Example 5.1.

5.2.3 Gate counts for the computation of wave functions

Theorem 5.1.

Proof.

5.3 The QLSA for the spectral discretisation

Theorem 5.2.

Lemma 3.3 (Estimation of $\|A^{-1}b\|$ ).