# Analog quantum simulation of partial differential equations

Shi Jin  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
School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai,
200240, China  Shanghai Artificial Intelligence Laboratory, Shanghai, China
Nana Liu<EMAIL_ADDRESS>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  Shanghai Artificial Intelligence Laboratory, Shanghai,
China  University of Michigan-Shanghai Jiao Tong University Joint Institute,
Shanghai 200240, China.

###### Abstract

Quantum simulators were originally proposed for simulating one partial
differential equation (PDE) in particular – Schrödinger’s equation. Can
quantum simulators also efficiently simulate other PDEs? While most
computational methods for PDEs – both classical and quantum – are digital
(PDEs must be discretised first), PDEs have continuous degrees of freedom.
This suggests that an analog representation can be more natural. While digital
quantum degrees of freedom are usually described by qubits, the analog or
continuous quantum degrees of freedom can be captured by qumodes. Based on a
method called Schrödingerisation, we show how to directly map $D$-dimensional
linear PDEs onto a $(D+1)$-qumode quantum system where analog or continuous-
variable Hamiltonian simulation on $D+1$ qumodes can be used. This very simple
methodology does not require one to discretise PDEs first, and it is not only
applicable to linear PDEs but also to some nonlinear PDEs and systems of
nonlinear ODEs. We show some examples using this method, including the
Liouville equation, heat equation, Fokker-Planck equation, Black-Scholes
equations, wave equation and Maxwell’s equations. We also devise new protocols
for linear PDEs with random coefficients, important in uncertainty
quantification, where it is clear how the analog or continuous-variable
framework is most natural. This also raises the possibility that some PDEs may
be simulated directly on analog quantum systems by using Hamiltonians natural
for those quantum systems.

##  I Introduction

One of the oldest and currently most promising application areas for quantum
devices is quantum simulation. Popularised by Feynman in the early 1980s, it
is important for more efficient simulation – compared to its classical
counterpart – of one special partial differential equation (PDE):
Schrödinger’s equation. If quantum simulators can be helpful for simulating
Schrödinger’s equation, it is hoped that they may also be helpful for
simulating other PDEs. As with large-scale quantum systems, classical methods
for other high-dimensional and large-scale PDEs often suffer from the curse-
of-dimensionality, which a quantum treatment might in certain cases be able to
mitigate. To enable simulation of PDEs on quantum devices that obey
Schrödinger’s equations, it is crucially important to develop good methods for
mapping other PDEs onto Schrödinger’s equations.

Before tackling the question of how to realise this mapping, we first need to
consider how PDEs are represented. Representation of information can be
broadly classed into digital, analog or their mixture. Digital information
uses discrete and finite degrees of freedom while analog uses continuous
and/or infinite degrees of freedom. The simulation of PDEs using any digital
representation first requires a discretisation of the PDE. Methods include
finite difference, finite element, finite volume, spectral methods, and
particle and Monte-Carlo methods. The implementation details and costs are
dependent on the discretisation chosen. The majority of classical and quantum
methods for PDEs today are based on this paradigm. The digital representation
of information in the quantum regime is typically via qubits. Many quantum
algorithms for simulating linear PDEs [1, 2, 3, 4, 5] and systems of nonlinear
ODEs and nonlinear PDEs [6, 7, 8, 9, 10, 11, 12, 13] based on qubits have been
theoretically proposed. For PDEs the algorithms can have the complexity
$O(\text{poly}(D,\log(1/N)))$ in the best case scenario – as a quantum
subroutine – where $D$ is the dimension and $N$ is the size of the
discretisation needed.

However, PDEs are intrinsically analog because they typically describe
dynamics in space and time which have continuous degrees of freedom. A simple
example of solving PDEs in the analog manner, by direct analogy and in common
use today, is wind tunnels, for solving Navier-Stokes equations. It is
purpose-built for simulating airflow over different geometries and is
sufficiently accurate for its purpose. A downside is the difficulty in
adapting this, and other specialised devices, to solve more general PDEs and a
more universal (though not necessarily fully universal) approach would be
preferred. Classical methods for solving ODEs in the analog manner, by
indirect analogy, include the use of electronic or mechanical components
described by classical physics. Solving PDEs can be difficult with these
devices, however, unless they are first discretised into a finite system of
ODEs, and the analog machines solve the resulting ODEs [14, 15]. To retain the
analog nature of PDEs, without first discretising, with the possibility of a
more universal approach to simulate large classes of PDEs and with improved
efficiency, we can look beyond classical analog systems to quantum systems.
See Figure 1 for a simplified categorisation.

Analog information in the quantum regime is commonly known as continuous-
variable (CV) quantum information [16, 17]. Many quantum systems, like quantum
optical systems and superconducting circuits, are naturally continuous-
variable. Indeed, quantum mechanics in the form of wave mechanics was
introduced by Schrödinger using wave functions, which is a continuous
description of information in the state, while the qubit representation was
introduced much later. Past methods proposed for CV quantum systems to
simulate linear PDEs focus on a CV version of a quantum linear systems solver
[18]. Methods for ODEs include approximating ODEs through a variational
quantum circuit [19] or using Hamiltonian simulation [20].

We focus on a new method by directly mapping general linear PDEs evolving in
time with $D$ spatial dimensions to a Schrödinger-like equation in $D+1$
spatial dimensions. This then enables us to use Hamiltonian simulation in
continuous time. This mapping – which is exact – is called Schrödingerisation,
recently introduced by the authors in [21, 22]. This method, in contrast to
methods based on qubits, enables one to easily identify the necessary
Hamiltonian to simulate any given linear PDE, see Figure 2. This simplicity
also enables an easier way of identifying possible quantum simulators that
could be used. Quantum simulators provide the most plausible near-term device
for efficiently simulating a wide range of Hamiltonians [23, 24]. The size of
the necessary quantum systems are small, requiring only $(D+1)$ qumodes, which
are the CV counterparts to qubits. Simulating the same PDE using qubits in the
near-term, on the other hand, even for low-dimensional PDEs, would be much
more difficult as much larger systems are required when $N$ is large, as
system sizes required is $\sim O(D\log N)$. Furthermore, algorithms based on
qubits require unitaries that need to be decomposed into smaller gates, rather
than relying on a single Hamiltonian that may have the possibility of being
naturally realised. While for CV the most general PDEs still require such a
decomposition into smaller CV gates, there is still the possibility to find
cases where such a decomposition is not necessary. It is also clear from our
CV formalism which parts of a given PDE quantum exploit entangling gates and
when non-Gaussian operations are required. Both of these are necessary to go
beyond efficient classical simulability [25].

In Section II, we briefly introduce PDEs and CV quantum systems. In Section
III, we begin with the simplest cases of PDEs to illustrate analog quantum
simulation for PDEs without the need of Schrödingerisation. For general linear
PDEs, we use the Schrödingerisation method in Section IV for homogeneous PDEs,
Section IV.1 for inhomogeneous PDEs and Section IV.2 for PDEs with more than
first-order time derivatives. We follow with some examples of important PDEs
in Section V, including the Liouville equation, heat equation, the Fokker-
Planck equation, the Black-Scholes equation, wave equation and Maxwell’s
equations (which is a system of PDEs). We then examine nonlinear PDEs and
systems of nonlinear ODEs in Section VII. In Section VI, we propose a new
quantum algorithm for uncertain linear PDEs, based on stochastic Galerkin
methods. Uncertain PDEs play an important role in the subject of uncertainty
quantification. For these algorithms it is particularly clear that the CV or
analog representation is the most natural and elegant.

We summarise our findings in Table 2.

Figure 1: Table summarising basic possible approaches to the simulation of
PDEs. The column ‘analog’ (A) refers to schemes representing the solution with
continuous-variables. The column ‘digital’ (D) refers to schemes representing
the solution with finite and discrete degrees of freedom. The (QA) regime
(orange box) refers to algorithms that simulate the PDEs with the use of
continuous-variable quantum systems. There is no need to discretise the PDEs
before mapping to the quantum system. This is the main focus of this paper, by
using a technique called Schrödingerisation, where the evolution in time is
also continuous and the corresponding Hamiltonian has a very simple-to-derive
form. The unitary generated by the Hamiltonian can either be realised directly
by mapping onto an existing quantum system (analogue quantum computing -
different notion to analog), or it can be decomposed into a sequence of gates
(universal continuous-variable quantum computation). The majority of proposed
quantum algorithms for PDEs reside in QD, where discretisation of the PDE is
necessary, before it can be mapped onto qubit dynamics. There are also methods
in CA, which are primarily for ODEs. To use these for PDEs, one must
discretise the PDEs spatially first into a system of ODEs, hence solving PDEs
require simultaneously solving many ODEs. The regime CD comprises the majority
of algorithms in classical scientific computing, like finite difference
methods, etc. Hybrid analog-digital scenarios in this framework are also
possible, for example see Sections IV.1, IV.2 and V.2.5. There are also hybrid
classical-quantum schemes, which we will not discuss in this paper.

.

##  II Background

###  II.1 Partial differential equations

The most general PDEs for $u(t,x_{1},...,x_{D})$ in $D$ spatial dimensions and
first-order in time has the form

$\displaystyle\frac{\partial u}{\partial t}+F(u,\nabla
u,\nabla^{2}u,...,\nabla^{n}u)=0,\qquad t\geq 0,\qquad
x_{j}\in\mathbb{R},\qquad j=1,...,D,$ (1)

where if $F$ is a linear function, then the above is a linear PDE. If $F$ is a
nonlinear function of any of its arguments $u,\nabla u,$…etc, then the above
is a nonlinear PDE. If all the terms in Eq. (1) involve $u$, then it is called
a homogeneous PDE. If there are also terms of the kind $f(x_{1},...,x_{D})\neq
0$, independent of $u$, then it is an inhomogeneous PDE.

For most PDEs, whether linear or nonlinear, analytical solutions seldom exist.
This means that most of the time one resorts to numerical algorithms that are
implemented on classical (digital) computers. These algorithms are among the
most important achievements in scientific and engineering computing. Different
types of partial differential equations have drastically different behavior in
their solutions. For example, they could differ in regularity (smoothness) or
have different physical properties (e.g. positivity, conservation, invariance,
entropy conditions,… etc). Thus different PDEs require drastically different
numerical strategies.

However, all these methods generally involve some kind of discretisation of
spatial degrees of freedom. For general PDEs, one can write it as a system of
ODEs by first discretising in space. There are two classes of discretisation
methods. One is the so-called Eulerian framework, in which one discretises the
spatial variables by a finite difference, finite element, finite volume, or
spectral method, on a fixed mesh. The advantage of Eulerian methods is their
high order accuracy but they suffer from the curse-of-dimensionality in high
spatial dimensions. The other is the Lagrangian method, including the so-
called particle, and mesh-freemethods. This is popular for high dimensional
problems, for instance particle or Monte-Carlo methods for kinetic equations
(the Boltzmann equation and Vlasov type equations), and vortex methods for
incompressible Euler and Navier-Stokes equations in fluid dynamics. These
Lagrangian methods do not suffer from the curse-of-dimensionality but they are
of lower order (typically only half to first order) methods, and also become
inaccurate in areas where there are not enough particles to accurately
approximate, for example, field equations (like in a Particle-in-Cell
simulation of Vlasov-Poisson or Vlasov-Maxwell equations in plasma).

The curse-of-dimensionality presents a major obstacle for numerical methods
for PDEs, and it arises chiefly out of the need to perform spatial
discretisation in the Eulerian framework. This is problematic especially for
systems with large dimensions (like $n$-body Schrödinger equations, kinetic
equations such as the Boltzmann equation, Landau and radiative transfer
equations, and PDEs with uncertainties) or with requirement of high
resolutions (for example problems with small scales and high frequencies).

To resolve this issue, quantum algorithms have been proposed on an analogous
problem, which is the simulation of quantum states whose amplitudes contain
the classical solutions. There is a large body of literature on this topic
where the bulk of these quantum methods focus on representing the discretised
solution in the form of qubits. For instance, consider a problem in $D$
dimensions to be solved with $N$ mesh points per space dimension, then there
are $N^{D}$ points to be used, which can be embedded into $d\log_{2}N$ qubits.
Then, when certain conditions hold (see references in the introduction), it
can be shown that the quantum cost can be polynomial in $D$ by using methods
mostly based on quantum linear algebra, thus greatly alleviating the curse-of-
dimensionality problem.

However, very few methods to date have taken advantage of existing continuous-
variable degrees of freedom in physical systems to simulate PDEs. By using
continuous-variable systems to perform these simulations, no discretisation is
actually necessary, at least in the ideal scenario. We will briefly introduce
quantum continuous-variable systems in the next section and later embed both
linear and nonlinear PDEs in these systems.

###  II.2 Continuous-variable quantum systems

Quantum information processing protocols are often introduced in terms of
qubits, which are the quantum counterparts to classical bit-strings
$\\{0,1\\}^{\otimes m}$. A qubit uses the eigenbasis
$\\{|0\rangle,|1\rangle\\}$ for the computational basis. An $m$-qubit system
is a tensor product of $m$ qubits which lives in $2^{m}$-dimensional Hilbert
space.

A continuous-variable (CV) quantum state, or ‘qumode’, on the other hand,
spans an infinite-dimensional Hilbert space. A qumode is the quantum analogue
of a continuous classical degree of freedom, like position, momentum or energy
before being quantised. A qumode is equipped with observables with a
continuous spectrum, such as the position $\hat{x}$ and momentum $\hat{p}$
observables of a quantum particle. Its eigenbasis can be chosen to be for
instance $\\{|x\rangle\\}_{x\in\mathbb{R}}$, which are the eigenstates of
$\hat{x}$. In this basis, a qumode can be expressed as
$|u\rangle=(1/\|\boldsymbol{u}\|)\int u(x)|x\rangle dx$, where
$\|\boldsymbol{u}\|^{2}=\int dx|u(x)|^{2}$ is the normalisation constant. A
system of $m$-qumodes is a tensor product of $m$ qumodes. See [26] for a good
introduction on the CV formalism and how it can be used in quantum
information.

Despite the relatively less well-known terminology of qumodes, physicists will
no doubt observe that $u(x)/\|\boldsymbol{u}\|$ is a familiar quantity as it
is precisely the wavefunction of the quantum state in the position basis.
Thus, despite CV quantum information processing and computation being less
well-known compared to processing with qubits, it is clear that these are the
natural states of quantum systems. Extra effort is generally required to
reduce the wavefunction or qumode description of quantum systems in the
laboratory to that of qubits.

The quadrature operators of a qumode are $\hat{x}$ and $\hat{p}$, where
$[\hat{x},\hat{p}]=i$. If we let $|x\rangle$ and $|p\rangle$ denote the
eigenvectors of $\hat{x}$ and $\hat{p}$ respectively, then $\langle
x|p\rangle=\exp(ixp)/\sqrt{2\pi}$. The position and momentum eigenstates each
form a complete eigenbasis so $\int dx|x\rangle\langle x|=I=\int
dp|p\rangle\langle p|$. Here the quantised momentum operator $\hat{p}$ is also
associated with the spatial derivative
$-i\hat{p}\leftrightarrow\partial/\partial x$. One can see this easily from
the same commutation relation being obeyed $[x,i\partial/\partial x]u=iu$.
Suppose we define the state $|\partial u/\partial x\rangle=\int(\partial
u/\partial x)|x\rangle dx$. Then it is straightforward to show
$-i\hat{p}|u\rangle=|\partial u/\partial x\rangle$. For simplicity, we ignore
normalisation constants and it is simple to prove the following. Let
$\tilde{u}(p)$ denote the Fourier transform of $u(x)$. Then
$-i\hat{p}|u\rangle=-i\int u(x)\hat{p}|x\rangle dx=-i\iint
u(x)\hat{p}|p\rangle\langle p|x\rangle dxdp=-i\iint u(x)p\exp(-ixp)|p\rangle
dxdp=-i\int\tilde{u}(p)p|p\rangle
dp=-i\iint\tilde{u}(p)p|x^{\prime}\rangle\langle x^{\prime}|p\rangle
dpdx^{\prime}=-i\iint\tilde{u}(p)p\exp(ix^{\prime}p)|x^{\prime}\rangle
dx^{\prime}dp=-\iint(\partial/\partial
x^{\prime})(\tilde{u}(p)\exp(ix^{\prime}p))|x^{\prime}\rangle
dx^{\prime}dp=\int\partial u(x^{\prime})/\partial x^{\prime}|x^{\prime}\rangle
dx^{\prime}=|\partial u/\partial x\rangle$. Similarly
$(-i\hat{p})^{n}|u\rangle=|\partial^{n}u/\partial x^{n}\rangle$. This simple
observation will be important for our simulation of PDEs later.

These quadrature operators can also be written in terms of creation and
annihiliation operators $\hat{a},\hat{a}^{\dagger}$. These are related to the
quadrature operators by $\hat{x}=(\hat{a}+\hat{a}^{\dagger})/\sqrt{2}$ and
$\hat{p}=(\hat{a}-\hat{a}^{\dagger})/(\sqrt{2}i)$. The number operator is
defined as $\hat{n}=\hat{a}^{\dagger}\hat{a}$ with corresponding eigenvector
$|n\rangle$ and eigenvalue $n$, where the latter is interpreted as particle
number. The eigenbasis $\\{|n\rangle\\}_{n=0}^{\infty}$ forms a so-called
number state or Fock state basis set. In this notation, $|0\rangle$ is the
vacuum state.

Consider a system of $m$-qumodes, then we define the position/momentum
operator only acting on the $j^{\text{th}}$ mode as

$\displaystyle\hat{p}_{j}=I^{\otimes j-1}\otimes\hat{p}\otimes I^{\otimes
m-j},\qquad\hat{x}_{j}=I^{\otimes j-1}\otimes\hat{x}\otimes I^{\otimes
m-j},\qquad[\hat{x}_{j},\hat{p}_{k}]=i\delta_{jk}.$ (2)

Measurement of the quadrature operators is also routinely performed, for
instance through homodyne detection. Mathematically, this corresponds to a
projective measurement of a quadrature operator. In fact, the expectation
value of any linear combination of quadratures
$\begin{pmatrix}\hat{x}_{\theta}\\\
\hat{p}_{\theta}\end{pmatrix}=\begin{pmatrix}\cos(\theta)&\sin(\theta)\\\
-\sin(\theta)&\cos(\theta)\end{pmatrix}\begin{pmatrix}\hat{x}_{0}\\\
\hat{p}_{0}\end{pmatrix}$, $\theta\in[0,\pi/2)$ can be measured using homodyne
detection in a rotated basis.

Operations on qumodes fall into two classes: Gaussian and non-Gaussian. A
universal gate set for qumodes can be formed from a set of Gaussian operations
and only a single kind of non-Gaussian operation. This is analogous to how a
universal gate set for qubits consist of Clifford and non-Clifford gates,
where Clifford gates are the analogy to Gaussian gates and non-Clifford gates
are analogous to non-Gaussian operations. A pure qumode Gaussian state can be
expressed as $U_{G}|0\rangle$ where $U_{G}$ is the Gaussian unitary operation
characterised by the affine map
$(\boldsymbol{S},\boldsymbol{d}):(\hat{x},\hat{p})^{T}\rightarrow\boldsymbol{S}(\hat{x},\hat{p})^{T}+\boldsymbol{d}$
for $\boldsymbol{S}\in\mathbb{R}^{2\times k}$,
$\boldsymbol{d}\in\mathbb{R}^{2}$. For example, a coherent state
$|\alpha\rangle$ is a Gaussian state prepared by a Gaussian operation
$U_{G}(\alpha)=\exp(\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a})=\exp(i\sqrt{2}(\hat{x}\text{Im}(\alpha)-\hat{p}\text{Re}(\alpha)))$.
The Hamiltonian generating this $U_{G}$ only contains terms that are linear in
the quadratures. Gaussian operations also contain operations whose generating
Hamiltonians are quadratic in the quadrature operations. For quadratic
operations on a single mode, examples include single-mode squeezing operators.
For quadratic operations across two-modes, an example is the so-called CZ gate
$U_{G}=\exp(i\hat{x}\otimes\hat{x})$. For operations across two-modes like the
CZ gate, we note that there is entanglement generation. For any order higher
than two, we refer to these as non-Gaussian operations.

While Gaussian operations are generally quite accessible in laboratory
settings, non-Gaussian operations are less so. However, they are known to be
necessary for universal CV computation [17] and also necessary to go beyond
efficient classical simulability [25]. While concatenations of Gaussian
operations only produce other Gaussian operations, only a single kind of non-
Gaussian operation and Gaussian operations are sufficient to produce any other
non-Gaussian operation. Thus, Gaussian operations with a single kind of non-
Gaussian gate forms a universal gate set for qumodes. A typical non-Gaussian
gate has the form $U_{nG}=\exp(it\hat{x}^{r})$, $r\geq 3$. Using methods
similar to Trotterisation – employing commutation relations of CV operators
$[\hat{x},\hat{p}]=i$ – any CV unitary can be decomposed into a finite number
of CV gates from the universal set [27].

#### Physical platforms

In this paper, we show how any PDE can be simulated using qumodes through
finding the appropriate generating Hamiltonian $\boldsymbol{H}$. We will find
that for some PDEs, only Gaussian operations are necessary while for others
non-Gaussian operations are essential. We also show how entangling gates like
CZ gates are necessary for most PDEs. These unitaries can be generated in one
of two main ways. The first way is by breaking the unitary up into the set of
gates in the CV universal gate set (with a set of Gaussian gates and a non-
Gaussian gate), like we mentioned above. This is simplest in the case of time-
independent Hamiltonians. The second way is through the so-called analogue
quantum simulation, meaning that one looks for Hamiltonians of systems that
naturally have the form one wants, up to the tuning of parameters in the
Hamiltonian. Here time-dependence of Hamiltonians is handled automatically.
The ease of using this approach depends on the individual experimental
platform. For both of these cases, there are many possible physical platforms.
See [28] for good discussions on the framework of analogue quantum simulation.

Quantum mechanical systems are naturally continuous-variable. Examples of
candidates for simulation of qumodes include superconducting circuits [29],
optical systems [30, 31], optomechanical oscillators [32], atoms in optical
lattices [33] and vibrational modes of trapped ions [34]. These can also be
used for simulation of hybrid continuous-variable discrete-variable systems
[35]. These are systems consisting of both qumodes and qubits, with coupling
between qumode and qubit degrees of freedom. We will see these hybrid
Hamiltonians for the simulation of inhomogeneous PDEs, PDEs with derivative in
time of more than first-order, and a finite system of PDEs.

##  III Simple examples: Linear convection and higher-order PDEs

We illustrate first with the linear convection equation – where a more general
approach is not necessary – to show how simply the continuous-variable
framework can be employed.

Consider the solutions $u(x_{1},...,x_{D},t)$ to a $D$-dimensional initial
value problem of the linear convection equation

$\displaystyle\frac{\partial u}{\partial t}-\sum_{j=1}^{D}a_{j}\frac{\partial
u}{\partial x_{j}}=0,\qquad u(0,x)=u_{0}(x),\qquad
x_{j}\in\mathbb{R}^{D},\qquad a_{j}\in\mathbb{R}$ (3)

with constant $a_{j}$. We can define an infinite-dimensional vector
$\boldsymbol{u}$, whose entries are $u(x,t)$ for different
$x=(x_{1},...,x_{D})$ for a specified $t$, so we write

$\displaystyle\boldsymbol{u}\equiv\int
u(x_{1},...,x_{D},t)|x_{1},x_{2},...,x_{D}\rangle dx_{1}...dx_{D}$ (4)

with $l_{2}$ norm $\|\boldsymbol{u}(t)\|$. The integral sign
$\int=\int_{-\infty}^{\infty}$ unless otherwise specified. By simulating the
solutions (in a quantum way), we mean to prepare the normalised quantum state
with $D$ continuous-variable modes

$\displaystyle|u(t)\rangle=\frac{1}{\|\boldsymbol{u}(t)\|}\int
u(x_{1},...,x_{D},t)|x_{1},x_{2},...,x_{D}\rangle
dx_{1}...dx_{D}=\frac{\boldsymbol{u}(t)}{\|\boldsymbol{u}(t)\|}.$ (5)

We use the observation that $-i\hat{p}_{j}\leftrightarrow\partial/\partial
x_{j}$ here (see Section II.2) and rewrite the $D$-dimensional convection
equation as

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\mathbf{H}\boldsymbol{u},\qquad\mathbf{H}=\sum_{j=1}^{D}\mathbf{H}_{j},\qquad\mathbf{H}_{j}=a_{j}\hat{p}_{j},\qquad\boldsymbol{H}_{j}=\boldsymbol{H}^{\dagger}_{j}.$
(6)

We note that this can be considered now as an (infinitely large) ODE system
since we have now converted differentiation with respect to $x$ by an operator
applied to $\boldsymbol{u}$. Since $a_{j}$ is a constant and $\boldsymbol{H}$
is time independent, this equation has the solution

$\displaystyle\boldsymbol{u}(t)=e^{-i\boldsymbol{H}t}\boldsymbol{u}(0),\qquad\text{similarly}\qquad|u(t)\rangle=e^{-i\boldsymbol{H}t}|u(0)\rangle$
(7)

where the latter equality holds since
$\|\boldsymbol{u}(t)\|=\|\boldsymbol{u}(0)\|$ because the convection equation
is a conservative equation. This means that preparing $|u(t)\rangle$ only
requires a sequence of displacement operators $U=\prod_{j=1}^{D}U_{j}$ where
$U_{j}=\exp(-i\mathbf{H}_{j}t)$. This only requires $D$ single-mode unitary
operations $U=\exp(-ia_{j}\hat{p}_{j}t)$ on mode $j$ and identity operations
on every other mode.

That this transformation works is not surprising, since the displacement
operator essentially performs a quantum Fourier transform, which is exactly
what turns a convection equation into a Schrödinger type equation. From this
decomposition, it looks like the solution associated with each dimension can
evolve separately, without interfering with the behaviour in other dimensions.
Thus no entangling gates are necessary. That only single-mode operations are
necessary, without mode coupling, is not surprising, since when $a_{j}$ is a
constant, evolution in each dimension is not coupled with other dimensions.
Therefore, the equation for each dimension can be solved separately, i.e.
$u(t,x_{1},...,x_{D})=u(t,x_{1}-a_{1}t,...,x_{D}-a_{D}t)$.

As a example, suppose we want to prepare $|u(t)\rangle$ for $D=2$, given
initial conditions $u(0,x)$. Suppose
$u(0,x_{1},x_{2})\propto\exp(-x_{1}^{2}/\sigma_{1}^{2}-x_{2}^{2}/\sigma_{2}^{2})$.
Then we can prepare $|u(0)\rangle$ by just preparing a product state of two
Gaussian states $|u(0)\rangle\propto\int
dx_{1}\exp(-x_{1}^{2}/\sigma_{1}^{2})|x_{1}\rangle\int
dx_{2}\exp(-x_{2}^{2}/\sigma_{2}^{2})|x_{2}\rangle$. To prepare the evolution,
one just needs to apply the displacement operator $\exp(-ic_{1}\hat{p}_{1}t)$
onto the first mode and $\exp(-ic_{2}\hat{p}_{2}t)$ onto the second mode of
the initial state $|u(0)\rangle$. The resulting state is thus $|u(t)\rangle$.
This is simply a product state of the independent evolution of the two initial
modes. Note that this simplicity comes from the fact that the initial state
was also in a product state.

One sees that this linear convection equation is particularly simple for two
reasons:

  1. (1)

Behaviour in each dimension is independent of each other, meaning that one can
deduce the final $D$-dimensional solution by propagating the time-dependence
of each dimension independently. As an example, if the initial data is
factorisable (e.g., can use the separation of variables), then the final state
is also factorisable. In this case, the solution in $D$ dimensions can be
obtained by solving for $D$ differential equations independently, where each
equation only depends on one of the original $D$ variables.

  2. (2)

$\boldsymbol{u}$ evolves via a unitary operation $\exp(-i\boldsymbol{H}t)$,
which means a quantum simulation can be used directly. This is not true for
most PDEs.

Next we give examples of PDEs which are simple in the sense of point (2)
above, where the solution $\boldsymbol{u}$ directly evolves through a unitary
operation. However, they are not necessarily simple in the sense of point (1).
In the latter case, we can have $x$-dependence in the coefficients and
couplings between the modes corresponding to different dimensions are
necessary. We will see that this coupling can be achieved with quantum
entangling operations.

Consider the following linear PDEs with odd powers of derivatives in $x$ and
we can obtain a dynamical equation for $\boldsymbol{u}$ of the form in Eq.
(6). For instance, consider the conservative $D$-dimensional PDEs of the form

$\displaystyle\frac{\partial u}{\partial
t}-\sum_{k=1}^{K}\sum_{j=1}^{D}a_{k,j}(x_{1},...,x_{j-1},x_{j+1},...,x_{D})\frac{\partial^{2k-1}u}{\partial
x_{j}^{2k-1}}=0,\qquad u(0,x)=u_{0}(x),\qquad x_{j}\in\mathbb{R}^{D},\qquad
a_{j}\in\mathbb{R}$ (8)

with highest derivative order $2K-1$. The coefficients $a_{k,j}$ can be any
$x$-dependent function, except without dependence on $x_{j}$. Note this
equation preserves the $l_{2}$ norm of $u$, a necessary condition for direct
quantum simulation. This is not necessarily the most physically-motivated
variable coefficient PDE example, but we are just using it to illustrate that
it is possible for a system to be simple in the sense of point (2) but not
point (1). Using the same $-i\hat{p}_{j}\leftrightarrow\partial/\partial
x_{j}$ correspondence, one can rewrite this $D$-dimensional equation as a
system similar to Eq. (6)

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{H}\boldsymbol{u},\qquad\boldsymbol{H}=\sum_{k=1}^{K}(-1)^{k+1}\sum_{j=1}^{D}\boldsymbol{H}_{kj}=\boldsymbol{H}^{\dagger},$
$\displaystyle\boldsymbol{H}_{kj}=a_{k,j}(\hat{x}_{1},...,\hat{x}_{j-1},\hat{x}_{j+1},...,\hat{x}_{D})\hat{p}_{j}^{2k-1}.$
(9)

As an example, if $a_{2,1}=x_{2}x_{5}$, then the corresponding operator
$\boldsymbol{H}_{21}=\hat{p}_{1}^{3}\hat{x}_{2}\hat{x}_{5}=\boldsymbol{H}_{21}^{\dagger}$.
Here we note that entangling gates are necessary, so it is no longer simple in
the sense of point (1).

Although it is not simple in the sense of point (1), this setup in simple in
the sense of point (2) due to the hermiticity of $\boldsymbol{H}$. The
hermitian property comes from the fact that each quadrature is hermitian and
they each act on a different mode. If $a_{k,j}$ would depend on $\hat{x}_{j}$,
then $a_{k,j}\hat{p}_{j}^{2k-1}$ is no longer hermitian. Thus for more general
PDEs, we do not have the same equation for $\boldsymbol{u}$ in Eq. (6), which
is already of a form directly simulable through a Schrödinger-like equation,
with dynamics determined by a unitary operator $\exp(-i\boldsymbol{H}t)$ when
$\boldsymbol{H}$ is time independent.

For more general linear PDEs, we can instead transform them into a Schrödinger
form through adding a single ancilla qumode. This is a simple transformation
called Schrödingerisation, recently introduced in [21, 22] and analysed for
qubit-based systems. We will study this procedure in the following Section IV
– for PDEs which are generally not simple with respect to either point (1) or
point (2) – and see how it can be naturally used in the continuous-variable
setting. In Section VIII, we also have an extended Schrödingerisation
formalism when $D$ ancilla modes are added instead.

##  IV Method for general linear PDEs: Schrödingerisation for continuous-
variable quantum systems

A linear homogeneous PDE for $u(x,t)$ with first derivative in $t$ (and any
order of derivative in $x$) can be written as

$\displaystyle\frac{\partial u}{\partial
t}+\sum_{k=1}^{K}\sum_{j=1}^{D}a_{k,j}(x_{1},...,x_{D})\frac{\partial^{k}u}{\partial
x_{j}^{k}}+b(x_{1},...,x_{D})u=0.$ (10)

We remark that the above coefficients can also depend on time $t$ and our
formalism below works in the same way. The coefficients corresponding even $k$
should have the sign $(-1)^{k/2}$ for stability. Since it will not affect the
demonstration of our method, we will not label this dependence unless when
relevant.

By embedding $u$ in a continuous-variable state $\boldsymbol{u}$, where
$\boldsymbol{u}(t)=\int_{-\infty}^{\infty}u(t,x)|x\rangle dx$ and we use
$x=x_{1},...,x_{D}$, one obtains the following equation for $\boldsymbol{u}$

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}(\hat{x}_{1},...,\hat{x}_{D},\hat{p}_{1},...,\hat{p}_{D})\boldsymbol{u},\qquad\boldsymbol{u}(0)$
(11)

where
$\boldsymbol{A}(\hat{x}_{1},...,\hat{x}_{D},\hat{p}_{1},...,\hat{p}_{D})$,
acting on $\boldsymbol{u}$, is a linear operator consisting of factors of
$\hat{x},\hat{p}$ depending on the form of the original PDE. This is easy to
see, since any $\partial^{l}/\partial x_{j}^{k}$ gives a contribution of
$(-i\hat{p}_{j})^{k}$, and any $x$-dependent coefficient $a(x_{1},...,x_{D})$
gives a contribution $a(\hat{x}_{1},...,\hat{x}_{D})$. Thus, for Eq. (10), one
has

$\displaystyle\boldsymbol{A}=\sum_{k=1}^{K}\sum_{j=1}^{D}a_{k,j}(\hat{x}_{1},...,\hat{x}_{D})(-i)^{k+1}\hat{p}_{j}^{k}-ib(\hat{x}_{1},...,\hat{x}_{D})$
(12)

We will examine $\boldsymbol{A}$ explicitly for different classes of PDEs
later.

Figure 2: Simple overview of the Schrödingerisation procedure

For the convection equation and certain higher-order PDEs in Section III,
$\boldsymbol{A}=\boldsymbol{H}$. This means one could directly simulate the
evolution of $\boldsymbol{u}$ with the unitary operator
$\exp(-i\boldsymbol{A}t)=\exp(-i\boldsymbol{H}t)$. For more general
$\boldsymbol{A}\neq\boldsymbol{A}^{\dagger}$, however, one cannot directly use
Hamiltonian simulation. For more general cases, we will use a procedure called
Schrödingerisation [21, 22]. We first observe we can always decompose any
operator into its hermitian and anihermitian components

$\displaystyle\boldsymbol{A}=\boldsymbol{A}_{1}-i\boldsymbol{A}_{2},\qquad\boldsymbol{A}_{1}=(\boldsymbol{A}+\boldsymbol{A}^{\dagger})/2=\boldsymbol{A}_{1}^{\dagger},\qquad\boldsymbol{A}_{2}=i(\boldsymbol{A}-\boldsymbol{A}^{\dagger})/2=\boldsymbol{A}_{2}^{\dagger}.$
(13)

We introduce a new parameter $\xi\in\mathbb{R}$ where for $\xi>0$, and define
$w(t,x,\xi)=\exp(-\xi)u(t,x)$ and
$\boldsymbol{w}(t,\xi)=\exp(-\xi)\boldsymbol{u}(t)$, which is called a warped
phase transformation. Then one can rewrite Eq. (11) as

$\displaystyle\frac{\partial\boldsymbol{w}(t,\xi)}{\partial
t}=\boldsymbol{A}_{2}\frac{\partial\boldsymbol{w}(t,\xi)}{\partial\xi}-i\boldsymbol{A}_{1}\boldsymbol{w}(t,\xi).$
(14)

We use even extension to the domain $\xi<0$ (so we can perform the Fourier
transform) by extending the initial condition evenly (see [21, 22]) to

$\displaystyle\boldsymbol{w}(0,\xi)=\boldsymbol{w}(t=0,\xi)=\exp(-|\xi|)\boldsymbol{u}(t=0)=\exp(-|\xi|)\boldsymbol{u}(0).$
(15)

Let $\tilde{\boldsymbol{w}}$ denote the Fourier transform of $\boldsymbol{w}$
with respect to $\xi$ with the corresponding Fourier mode $\eta$, i.e.
$\tilde{\boldsymbol{w}}(t,\eta)=\int
d\xi\exp(-i\xi\eta)\boldsymbol{w}(t,\xi)$. Then Eq. (14) becomes the following
set of equations for all $\eta$

$\displaystyle\frac{d\tilde{\boldsymbol{w}}(t,\eta)}{dt}=-i(\eta\boldsymbol{A}_{2}+\boldsymbol{A}_{1})\tilde{\boldsymbol{w}}(t,\eta),\qquad\tilde{\boldsymbol{w}}(0,\eta)=\frac{2}{1+\eta^{2}}\boldsymbol{u}(0).$
(16)

Now we can define the infinite-dimensional $\tilde{\boldsymbol{w}}(t,\eta)$
for all values of $\eta$ by $\tilde{\boldsymbol{v}}(t)$, so we can write
$\tilde{\boldsymbol{v}}(t)=\iint\tilde{w}(t,x,\eta)|x\rangle|\eta\rangle
dxd\eta$. Here $\hat{\eta}$ can be any quadrature, e.g.
$\cos(\theta)\hat{x}+\sin(\theta)\hat{p}$ for any $\theta\in[0,2\pi)$. We can
define an operator $\hat{\eta}$ such that
$\hat{\eta}|\eta\rangle=\eta|\eta\rangle$. This means
$(I\otimes\hat{\eta})\tilde{\boldsymbol{v}}(t)=\iint\eta\tilde{w}(t,x,\eta)|x\rangle|\eta\rangle
dxd\eta$. Then $\tilde{\boldsymbol{v}}(t)$ obeys

$\displaystyle\frac{d\tilde{\boldsymbol{v}}(t)}{dt}=-i(\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I)\tilde{\boldsymbol{v}}(t)=-i\boldsymbol{H}\tilde{\boldsymbol{v}}(t),\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger},$
$\displaystyle\tilde{\boldsymbol{v}}(0)=\int\frac{2}{1+\eta^{2}}|\eta\rangle
d\eta\boldsymbol{u}(0)=\int e^{-|\xi|}|\xi\rangle
d\xi\boldsymbol{u}(0)=|\Xi\rangle\boldsymbol{u}(0),\qquad|\Xi\rangle=\int
e^{-|\xi|}|\xi\rangle d\xi.$ (17)

This equation is now of the form where we can directly use Hamiltonian
simulation, where
$\boldsymbol{H}=\boldsymbol{H}(\hat{x}_{1},...,\hat{x}_{D},\hat{p}_{1},...,\hat{p}_{D})$.
In the case the coefficients of the PDE are time-independent, then
$\boldsymbol{H}$ is also time-independent and the solution is
$\tilde{\boldsymbol{v}}(t)=\exp(-i\boldsymbol{H}t)\tilde{\boldsymbol{v}}(0)$.
Now
$|v(t)\rangle=(1/\|\boldsymbol{v}\|)\iint_{-\infty}^{\infty}w(t,x,\xi)|x\rangle|\xi\rangle
dxd\xi$ with
$\|\boldsymbol{v}(t)\|=\|\boldsymbol{v}(0)\|\equiv\|\boldsymbol{v}\|$. This
equality holds since $w(t,x,\xi)$ obeys a conservative equation and hence the
$l_{2}$ norm is preserved. Thus we can write

$\displaystyle|\tilde{v}(t)\rangle=e^{-i\boldsymbol{H}t}|\tilde{v}(0)\rangle,\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I,\qquad|\tilde{v}(0)\rangle=|\Xi\rangle|u(0)\rangle.$ (18)

Now, for the PDE in Eq. (10) for $a_{k,j},b\in\mathbb{R}$ with
$\boldsymbol{A}$ given in Eq. (12), we can easily write the Hamiltonian as

$\displaystyle\boldsymbol{H}=b\otimes\hat{\eta}+\frac{i}{2}\sum_{k=1}^{K}\sum_{j=1}^{D}(-1)^{k+1}[a_{k,j},\hat{p}_{j}^{k}]\otimes\hat{\eta}+\frac{1}{2}\sum_{k=1}^{K}\sum_{j=1}^{D}(-i)^{k+1}(a_{k,j}\hat{p}_{j}^{k}+\hat{p}_{j}^{k}a_{k,j})\otimes
I.$ (19)

This form can be greatly simplified as we deal with specific PDEs, as we will
see in the following sections. This simple procedure is summarised in Figure
2. From this it is also simple to see how, when given any $\boldsymbol{H}$,
one can work backwards to find the corresponding PDE.

Given $|\tilde{v}(t)\rangle$, we can perform an inverse quantum Fourier
transform (continuous) on $|\tilde{v}(t)\rangle$ to retrieve $|v(t)\rangle$.
This is straightforward to do for instance on optical systems and is
essentially just a change in basis from one quadrature to its conjugate. We
can take a projection
$\hat{P}_{>0}=I\otimes\int_{0}^{\infty}|\xi^{\prime}\rangle\langle\xi^{\prime}|d\xi^{\prime}$
to select only the $\xi^{\prime}>0$ parts of the state $|v(t)\rangle$. Here
$\xi^{\prime}$ can represent any quadrature, e.g. momentum. Since
$w(t,x,\xi)=e^{-\xi}u(t,x)$ for $\xi>0$ [8], then one can see

$\displaystyle\hat{P}_{>0}|v(t)\rangle=\frac{1}{\|\boldsymbol{v}\|}\int_{-\infty}^{\infty}\int_{0}^{\infty}w(t,x,\xi)|x\rangle|\xi\rangle
dxd\xi=\frac{1}{\|\boldsymbol{v}\|}\int_{0}^{\infty}e^{-\xi}|\xi\rangle
d\xi\int_{-\infty}^{\infty}u(t,x)|x\rangle dx$
$\displaystyle=\frac{\|\boldsymbol{u}(t)\|}{\sqrt{2}\|\boldsymbol{v}\|}|\Xi^{\prime}\rangle|u(t)\rangle=\frac{\|\boldsymbol{u}(t)\|}{\sqrt{2}\|\boldsymbol{u}(0)\|}|\Xi^{\prime}\rangle|u(t)\rangle$
(20)

where $|\Xi^{\prime}\rangle=\sqrt{2}\int_{0}^{\infty}\exp(-\xi)|\xi\rangle
d\xi$. The second equality uses $\int_{0}^{\infty}w(t,x,\xi)d\xi=u(t,x)$. Thus
the probability of retrieving the state $|\Xi\rangle|u(t)\rangle$ after the
projective measurement $\hat{P}_{>0}$ is given by
$\left(\|\boldsymbol{u}(t)\|/(\sqrt{2}\|\boldsymbol{v}\|)\right)^{2}=(\|\boldsymbol{u}(t)\|^{2}/2)\iint_{-\infty}^{\infty}w^{2}(0,x,\xi)dxd\xi$.
Using the initial condition $w(0,x,\xi)=\exp(-|\xi|)u(0,x)$ for all $\xi$,
then one can write
$\|\boldsymbol{v}\|^{2}=\int_{-\infty}^{\infty}e^{-2|\xi|}d\xi\int_{-\infty}^{\infty}u^{2}(0,x)dx=\|\boldsymbol{u}(0)\|^{2}$.
Thus the probability of retrieving $|u(t)\rangle$ from $|v(t)\rangle$ through
$\hat{P}_{>0}$ is $\|\boldsymbol{u}(t)\|^{2}/(2\|\boldsymbol{u}(0)\|^{2})$.
This factor is an inevitable consequence of mapping any PDE into the
Schrödinger-like conservative PDE. Any method that maps linear PDEs into a
Schrödinger-like conservative PDE will have this factor.

We see that in all these cases, we require the preparation of the initial
state $|\tilde{v}(0)\rangle=|\Xi\rangle|u(0)\rangle$ where $|\Xi\rangle=\int
e^{-|\xi|}|\xi\rangle d\xi$. However, $|\Xi\rangle$ is an unusual state which
is not necessarily easily attainable in a laboratory setting. However, it is
possible to approximate this state for instance by the experimentally
accessible Gaussian state
$|G\rangle=\int(1/(\sqrt{s}\pi^{1/4}))\exp(-\xi^{2}/(2s^{2}))|\xi\rangle d\xi$
with squeezing parameter $s$. Then the quantum fidelity between $|G\rangle$
and $|\Xi\rangle$ is
$|\langle\Xi|G\rangle|=\sqrt{2s}\exp(s^{2}/2)\pi^{1/4}\text{Erfc}(s/\sqrt{2})$.
It’s clear that there is some optimal $s$ in which we can attain maximum
overlap. In fact, the maximum fidelity $|\langle\Xi|G\rangle|\approx 0.986$ is
near $1$ at $s\approx 0.925$ which is almost near unity. Now if the initial
condition $u(0,x)$ is a Gaussian function with respect to $x$, then
$|u(0)\rangle$ can also be efficiently prepared. For a discussion on the error
propagation due to using a different ancilla state, see Appendix A.

We can also generalise the method in this section to deal with inhomogeneous
PDEs and higher-order derivatives in $t$ in the following sections.

Figure 3: Circuit for homogeneous linear PDE with first-order derivatives in
time. Here $\hat{P}_{>0}$ requires only projection onto $\xi>0$, and not any
particular value of $\xi$.

###  IV.1 Inhomogeneous PDE

A linear inhomogeneous PDE for $u(x,t)$ with first derivative in $t$ when
$f\neq 0$ can be written as

$\displaystyle\frac{\partial u}{\partial
t}+\sum_{k=1}^{K}\sum_{j=1}^{D}a_{k,j}(x_{1},...,x_{D})\frac{\partial^{k}u}{\partial
x_{j}^{k}}+b(x_{1},...,x_{D})u=f(x_{1},...,x_{D}).$ (21)

We denote that vector giving the value $f(x_{1},...,x_{D})$ at the point
$\boldsymbol{x}=(x_{1},...,x_{D})$ by $\boldsymbol{f}$. One can then turn the
inhomogeneous equation

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}(\hat{x}_{1},...\hat{x}_{D},\hat{p}_{1},...,\hat{p}_{D})\boldsymbol{u}+\boldsymbol{f}(\hat{x}_{1},...,\hat{x}_{D})$
(22)

into a homogeneous one by absorbing the inhomogeneous term $\boldsymbol{f}$
[22]. By dilating the system
$\boldsymbol{u}\rightarrow\boldsymbol{y}=\boldsymbol{u}\otimes|0\rangle+\boldsymbol{f}\otimes|1\rangle$,
from Eq. (30) one can derive

$\displaystyle\frac{d\boldsymbol{y}}{dt}=\frac{d}{dt}\begin{pmatrix}\boldsymbol{u}\\\
\boldsymbol{f}\end{pmatrix}=-i\boldsymbol{B}\boldsymbol{y},\qquad\boldsymbol{B}=\begin{pmatrix}\boldsymbol{A}&iI\\\
0&0\end{pmatrix},\qquad\boldsymbol{A}=\boldsymbol{A}_{1}-i\boldsymbol{A}_{2},\qquad\boldsymbol{B}=\boldsymbol{B}_{1}-i\boldsymbol{B}_{2},$
$\displaystyle\boldsymbol{B}_{1}=(\boldsymbol{B}+\boldsymbol{B}^{\dagger})/2=\boldsymbol{B}_{1}^{\dagger},\qquad\boldsymbol{B}_{2}=i(\boldsymbol{B}-\boldsymbol{B}^{\dagger})/2=\boldsymbol{B}_{2}^{\dagger},$
$\displaystyle\boldsymbol{B}_{1}=\begin{pmatrix}\boldsymbol{A}_{1}&iI/2\\\
-iI/2&0\end{pmatrix}=\boldsymbol{A}_{1}\otimes\frac{1}{2}(I+\sigma_{z})+\frac{I}{2}\otimes\sigma_{y},\qquad\boldsymbol{B}_{2}=\begin{pmatrix}\boldsymbol{A}_{2}&-I/2\\\
-I/2&0\end{pmatrix}=\boldsymbol{A}_{2}\otimes\frac{1}{2}(I+\sigma_{z})-\frac{I}{2}\otimes\sigma_{x}.$
(23)

We can Schrödingerise this new system (homogeneous with respect to
$\boldsymbol{y}$) where we transform
$\boldsymbol{y}\rightarrow\boldsymbol{v}\rightarrow\tilde{\boldsymbol{v}}$ to
obtain

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{B}_{2}\otimes\hat{\eta}+\boldsymbol{B}_{1}\otimes
I=\boldsymbol{H}^{\dagger}.$ (24)

This means that one just needs to include one extra ancilla qubit to the
$(D+1)$-qumodes. To recover $|u\rangle$, we first recover $|y\rangle$ using
the same method before for homogeneous PDEs and then project onto the
$|0\rangle$ qubit. The generating Hamiltonian $\boldsymbol{H}$ thus requires
an extra coupling to $\sigma_{x},\sigma_{y},\sigma_{z}$.

Figure 4: Circuit for inhomogeneous linear PDE with first-order derivatives in
time. Here $\hat{P}_{>0}$ requires only projection onto $\xi>0$, and not any
particular value of $\xi$. In the central register, $|0\rangle\langle 0|$
denotes projection onto the qubit $|0\rangle$.

###  IV.2 PDE with higher-order derivatives in $t$

We illustrate how our formalism can be applied when the PDE has second-order
derivatives in $t$. We look at the following homogeneous linear PDE (includes
for example the wave equation as a special case)

$\displaystyle\frac{\partial^{2}u}{\partial
t^{2}}+c_{0}(x_{1},...,x_{D})\frac{\partial u}{\partial
t}+\sum_{j=1}^{d}c_{j}(x_{1},...,x_{D})\frac{\partial^{2}u}{\partial
x_{j}\partial
t}-\sum_{l=1}^{L}\sum_{j=1}^{D}a_{j,l}(x_{1},...,x_{D})\frac{\partial^{l}u}{\partial
x_{j}^{l}}+b(x_{1},...,x_{D})u=0.$ (25)

One can rewrite Eq. (25) into the following equation for $\boldsymbol{u}$

$\displaystyle\frac{d^{2}\boldsymbol{u}}{dt^{2}}+\boldsymbol{\Gamma}\frac{d\boldsymbol{u}}{dt}+i\boldsymbol{A}\boldsymbol{u}=0,\qquad\boldsymbol{\Gamma}\equiv\boldsymbol{c}_{0}-i\sum_{j=1}^{D}\boldsymbol{c}_{j}(I^{\otimes
j-1}\otimes\hat{p}_{j}\otimes I^{\otimes
D-j}),\qquad\boldsymbol{\Gamma}=\boldsymbol{\Gamma}_{1}-i\boldsymbol{\Gamma}_{2}$
$\displaystyle\boldsymbol{\Gamma}_{1}=(\boldsymbol{\Gamma}+\boldsymbol{\Gamma}^{\dagger})/2=\boldsymbol{\Gamma}_{1}^{\dagger},\boldsymbol{\Gamma}_{2}=i(\boldsymbol{\Gamma}-\boldsymbol{\Gamma}^{\dagger})/2=\boldsymbol{\Gamma}_{1}^{\dagger},$
$\displaystyle\boldsymbol{A}=-i\left(\sum_{l=1}^{L}\sum_{j=1}^{D}a_{j,l}(\hat{x_{1}},...,\hat{x}_{D})(i\hat{p}_{j})^{l}+b(\hat{x}_{1},...,\hat{x}_{D})\right)$
(26)

where $\boldsymbol{c}_{i}=c_{i}(\hat{x_{1}},...,\hat{x}_{D})$. The goal is to
transform this into an equation of the form Eq. (11).

We can dilate the system
$\boldsymbol{u}\rightarrow\boldsymbol{y}=\boldsymbol{u}\otimes|0\rangle+(d\boldsymbol{u}/dt)\otimes|1\rangle$.
Then the following equation for $\boldsymbol{y}$ can be obtained

$\displaystyle\frac{d\boldsymbol{y}}{dt}=\begin{pmatrix}d\boldsymbol{u}/dt\\\
d^{2}\boldsymbol{u}/dt^{2}\end{pmatrix}=-i\boldsymbol{V}\boldsymbol{y},\qquad\boldsymbol{V}=\begin{pmatrix}0&iI\\\
\boldsymbol{A}&-i\boldsymbol{\Gamma}\end{pmatrix},\qquad\boldsymbol{V}=\boldsymbol{V}_{1}-i\boldsymbol{V}_{2},\qquad\boldsymbol{A}=\boldsymbol{A}_{1}-i\boldsymbol{A}_{2}$
(27)
$\displaystyle\boldsymbol{V}_{1}=(\boldsymbol{V}+\boldsymbol{V}^{\dagger})/2=\boldsymbol{V}_{1}^{\dagger},\boldsymbol{V}_{2}=i(\boldsymbol{V}-\boldsymbol{V}^{\dagger})/2=\boldsymbol{V}_{2}^{\dagger},\qquad\boldsymbol{A}_{1}=(\boldsymbol{A}+\boldsymbol{A}^{\dagger})/2=\boldsymbol{A}_{1}^{\dagger},\boldsymbol{A}_{2}=i(\boldsymbol{A}-\boldsymbol{A}^{\dagger})/2=\boldsymbol{A}_{2}^{\dagger},$
$\displaystyle\boldsymbol{V}_{1}=\frac{1}{2}\begin{pmatrix}0&\boldsymbol{A}^{\dagger}+iI\\\
\boldsymbol{A}-iI&-i(\boldsymbol{\Gamma}-\boldsymbol{\Gamma}^{\dagger})\end{pmatrix}=\boldsymbol{A}_{1}\otimes\sigma_{x}+\boldsymbol{A}_{2}\otimes\sigma_{y}-\frac{I}{2}\otimes\sigma_{y}+\boldsymbol{\Gamma}_{2}\otimes\frac{1}{2}(I-\sigma_{z}),$
$\displaystyle\boldsymbol{V}_{2}=\frac{1}{2i}\begin{pmatrix}0&-\boldsymbol{A}^{\dagger}+iI\\\
\boldsymbol{A}+iI&-i(\boldsymbol{\Gamma}+\boldsymbol{\Gamma}^{\dagger})\end{pmatrix}=\boldsymbol{A}_{2}\otimes\sigma_{x}-\boldsymbol{A}_{1}\otimes\sigma_{y}+\frac{I}{2}\otimes\sigma_{x}-\boldsymbol{\Gamma}_{1}\otimes\frac{1}{2}(I-\sigma_{z}).$
(28)

We can then proceed to Schrödingerise the system
$d\boldsymbol{y}/dt=-i\boldsymbol{V}\boldsymbol{y}$ as before. By applying the
transformation $\boldsymbol{y}\rightarrow\tilde{\boldsymbol{v}}$ in
Schröingerisation one has

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}+i(\boldsymbol{V}_{2}\otimes\hat{\eta}+\boldsymbol{V}_{1}\otimes
I)\tilde{\boldsymbol{v}}=\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{V_{2}}\otimes\hat{\eta}+\boldsymbol{V}_{1}\otimes
I=\boldsymbol{H}^{\dagger}.$ (29)

It is straightforward to generalise to higher-order derivatives in $t$. For
instance, it is simple to see that if there is at most an
$n^{\text{th}}$-order derivative in $t$, then we need $\log_{2}(n)$ auxiliary
qubits.

Figure 5: Circuit for homogeneous linear PDE with up to second-order
derivatives in time. $\hat{P}_{>0}$ requires only projection onto $\xi>0$, and
not any particular value of $\xi$. In the central register, $|0\rangle\langle
0|$ denotes projection onto the qubit $|0\rangle$.

##  V Examples of Schrödingerisation with linear PDEs

In this section, we show examples of different linear PDEs and the
Hamiltonians $\boldsymbol{H}$ of their corresponding Schrödingerised systems.
For PDEs with $D$ spatial dimensions, we see that no more than $D+1$ qumodes
are necessary. In some cases, only two-mode Gaussian interactions are
sufficient, while other PDEs require non-Gaussian interactions. This depends
primarily on the $x$-dependence of the coefficients in the PDE and the
highest-order spatial derivative in the PDE.

###  V.1 Linear first-order PDEs

Here we focus on linear first-order homogeneous PDEs for
$u(t,x_{1},...,x_{D})$

$\displaystyle\frac{\partial u}{\partial
t}+\sum_{j=1}^{D}a_{j}(x_{1},...,x_{D})\frac{\partial u}{\partial
x_{j}}+b(x_{1},...,x_{D})u=0,\qquad a_{j}\in\mathbb{R},b\in\mathbb{C}.$ (30)

The inhomogeneous $f\neq 0$ case is treated in Section IV.1 and is
straightforward to apply here. Here $a_{j}$ must be real otherwise the
equation is unstable. We can rewrite

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=-\sum_{j=1}^{D}a_{j}(\hat{x}_{1},...,\hat{x}_{D})\hat{p}_{j}-ib(\hat{x}_{1},...,\hat{x}_{D}),\qquad\boldsymbol{A}=\boldsymbol{A}_{1}-i\boldsymbol{A}_{2}$
(31)

where $\boldsymbol{A}_{1}$ and $\boldsymbol{A}_{2}$ are hermitian operators
defined in Eq. (13). Following Section IV, we can Schrödingerise the system

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger}$
$\displaystyle\boldsymbol{A}_{1}=\frac{1}{2}\sum_{j=1}^{D}\left(a_{j}(\hat{x}_{1},...,\hat{x}_{D})\hat{p}_{j}+\hat{p}_{j}a_{j}(\hat{x}_{1},...,\hat{x}_{D})\right)+\frac{i}{2}(b(\hat{x}_{1},...,\hat{x}_{D})-b^{*}(\hat{x}_{1},...,\hat{x}_{D})),$
$\displaystyle\boldsymbol{A}_{2}=\frac{i}{2}\sum_{j=1}^{D}[\hat{p}_{j},a_{j}(\hat{x}_{1},...,\hat{x}_{D})]+\frac{1}{2}(b(\hat{x}_{1},...,\hat{x}_{D})+b^{*}(\hat{x}_{1},...,\hat{x}_{D}))$
(32)

We can summarise the necessary $\boldsymbol{H}$ in Table 1, where for brevity
we only write out the highest-order terms. In all the examples in the table,
there are no parts of the Hamiltonian more difficult than single-mode and two-
mode Gaussian operations, except for when $\real(\gamma_{jk})\neq 0$ in the
last line.

Table 1: Highest-order Hamiltonian terms $\boldsymbol{H}$ for simulating
linear first-order homogeneous PDEs with mostly Gaussian operations. Below
$\hat{\Theta}\equiv(1/2)\sum_{jk}\alpha_{jk}(\hat{x}_{j}\hat{p}_{k}+\hat{p}_{k}\hat{x}_{j})\otimes
I$, $\alpha_{jk}\in\mathbb{R}$ and $\beta_{j}$, $\gamma_{jk}\in\mathbb{C}$.

|  $\qquad a_{j}=$ constant | $\qquad a_{j}=\sum_{k=1}^{D}\alpha_{jk}x_{k}$
---|---|---
$b=$ constant | $\real(b)I\otimes\hat{\eta}+\sum_{j}a_{j}\hat{p}_{j}\otimes I$ | $\hat{\Theta}$
$b=\sum_{j=1}^{D}\beta_{j}x_{j}$ | $\sum_{j}\real(\beta_{j})\hat{x}_{j}\otimes\hat{\eta}$ | $\sum_{j}\real(\beta_{j})\hat{x}_{j}\otimes\hat{\eta}+\hat{\Theta}$
$b=\sum_{j,k=1}^{D}\gamma_{jk}x_{j}x_{k}$ |  $\sum_{j,k}\real(\gamma_{jk})\hat{x}_{j}\hat{x}_{k}\otimes\hat{\eta}-\sum_{jk}\imaginary(\gamma_{jk})\hat{x}_{j}\hat{x}_{k}\otimes I$ |  $\sum_{jk}\real(\gamma_{jk})\hat{x}_{k}\hat{p}_{j}\otimes\hat{\eta}-\sum_{jk}\imaginary(\gamma_{jk})\hat{x}_{j}\hat{x}_{k}\otimes I+\hat{\Theta}$

Note that if $a$, $b$ are of any higher order polynomials, one would require
non-Gaussian gates. For instance, quadratic functions could be made possible
with Kerr nonlinearity. For the examples in the table, we see that certain
PDEs require only single-mode operations, with generating Hamiltonians like
$\hat{\eta},\hat{x}\hat{p}$, $\hat{x}^{2}$, while others require pairwise
entangling operations like $\hat{\eta}\otimes\hat{x}$. Those that only require
single-mode operations are simple in the sense of point (1) in Section III,
like the linear convection equation.

A well-known first-order linear PDE for example is the Liouville equation
which is Eq. (30) with $b=0$. In the case of linear dependence with
$a_{j}=\sum_{k=1}^{D}\alpha_{kj}x_{k}$,
$\boldsymbol{H}=(1/2)\sum_{j=1}^{D}\alpha_{jj}I\otimes\hat{\eta}+\hat{\Theta}\otimes
I$, which requires pairwise interactions of the form
$\hat{x}_{j}\hat{p}_{k}+\hat{p}_{k}\hat{x}_{j}$ when $j\neq k$. So this is not
simple in the sense of point (1) in Section III, but we require at most
Gaussian operations.

###  V.2 Examples of important linear PDEs

Here we give a few physically important examples of linear PDEs. Most of these
have second-order spatial derivatives with the exception of Maxwell’s
equation. There do exist physically important PDEs with spatial derivatives of
orders higher than two, for example the KdV equation and the quantum
hydrodynamic equations (with Bohm potential) which contain third order
derivatives, and the Kuramoto-Sivashinsky equation and the Cahn-Hilliard
equations which have fourth order derivatives. However, these are typically
nonlinear PDEs. We concentrate on linear PDEs in this section with second
order spatial derivatives and discuss nonlinear PDEs in Section VII.

We will start with problems with first-order derivatives in $t$, including the
heat equation, Fokker-Planck equation and Black-Scholes equation. Then we
proceed to the wave equation, with second-order in time and space derivatives.
For simplicity in these examples we focus on the homogeneous case and the
extension to $f\neq 0$ can be done straightforwardly using the formulation in
Section IV.1. Finally, we show Maxwell’s equations which is a system of
inhomogeneous linear first-order PDEs.

####  V.2.1 The heat equation

Consider the $D$-dimensional general heat equation with a forcing term

$\displaystyle\frac{\partial u}{\partial
t}+\sum_{i=1}^{D}\frac{\partial}{\partial
x_{i}}\left(\sum_{j=1}^{D}D_{ij}(x_{1},...,x_{D})\frac{\partial u}{\partial
x_{j}}\right)+V(x_{1},...,x_{D})u=0,\qquad
D_{ij}(x_{1},...,x_{D}),\,V(x_{1},...,x_{D})\in\mathbb{R},$ (33)

where $D_{ij}$ are the components of the diffusion matrix. This heat equation
can be transformed into

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=i\sum_{i,j=1}\hat{p}_{i}D_{ij}(\hat{x}_{1},...,\hat{x}_{D})\hat{p}_{j}-iV(\hat{x}_{1},...,\hat{x}_{D}).$
(34)

Applying the results of Section IV for the homogeneous case, the
Schrödingerised system becomes

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger},$
$\displaystyle\boldsymbol{A}_{1}=\frac{1}{2}\sum_{i,j=1}^{D}[\hat{p}_{j}D_{ij},\hat{p}_{i}]$
$\displaystyle\boldsymbol{A}_{2}=\frac{1}{2}\sum_{i,j=1}^{D}(\hat{p}_{i}D_{ij}\hat{p}_{j}+\hat{p}_{j}D_{ij}\hat{p}_{i})+V(\hat{x}_{1},...,\hat{x}_{D}).$
(35)

Suppose one has a constant diffusion coefficient
$D_{ij}(x)=a\delta_{ij}\in\mathbb{R}$ and a linear driving term
$V(x_{1},...,x_{D})=\sum_{j=1}^{D}k_{j}x_{j}$, both of which are time-
independent. Then $\boldsymbol{A}_{1}=0$ and we only have the
$\boldsymbol{A}_{2}$ term. The Hamiltonian $\boldsymbol{H}$ in this case
generates an entangling gate where
$\boldsymbol{H}=\sum_{j=1}^{D}(a\hat{p}_{j}^{2}+k_{j}\hat{x}_{j})\otimes\hat{\eta}$.
Thus, the most difficult gate requires a pairwise third-order operator of the
form $\exp(-i\hat{p}^{2}\otimes\hat{\eta}t)$. One can either identify an
analog quantum system to naturally realise this Hamiltonian, or easily
decompose this into more elementary continuous-variable gates. Observe that
$\boldsymbol{H}=\sum_{j=1}^{D}\boldsymbol{H}_{j}$ such that
$[\boldsymbol{H}_{i},\boldsymbol{H}_{j}]=0,$ $\forall\,i,j=1,...,D$, so one
can factorise
$\exp(-i\boldsymbol{H}t)=\exp(-i\boldsymbol{H}_{1}t)...\exp(-i\boldsymbol{H}_{D}t)$.
Here $\boldsymbol{H}_{j}=\boldsymbol{h}_{j,p}+\boldsymbol{h}_{j,x}$,
$\boldsymbol{h}_{j,p}\equiv a_{j}\hat{p}_{j}^{2}\otimes\hat{\eta}$,
$\boldsymbol{h}_{j,x}\equiv k_{j}\hat{x}_{j}\otimes\hat{\eta}$. Using the
Campbell-Baker-Hausdorff relation with $[\hat{x},\hat{p}^{2}]=2i\hat{p}$, one
obtains the commutation relation
$[\boldsymbol{h}_{j,x},\boldsymbol{h}_{j,p}]=2ia_{j}k_{j}\hat{p}_{j}\otimes\hat{\eta}=(2i/a_{j})\boldsymbol{h}_{j,p}$,
then one can decompose
$\exp(-i\boldsymbol{H}_{j}t)\propto\exp(i\boldsymbol{h}_{j,x}t)\exp(i\boldsymbol{h}_{j,p}t/a_{j})$.
Thus in $D$ dimensions, we can decompose $\exp(-i\boldsymbol{H}t)$ into $2D$
elementary continuous-variable gates. These are two-mode entangling gates: $D$
two-mode Gaussian gates of the form $\exp(-i\hat{x}\otimes\hat{\eta}t)$ and
$D$ nonlinear gates of the form $\exp(-i\hat{p}^{2}\otimes\hat{\eta}t)$.

####  V.2.2 The Fokker-Planck equation

The Fokker-Planck equation gives the time evolution of the probability density
function $u(t,x_{1},...,x_{D})$ of the velocity of a particle under the impact
of drag and random forces. In the Fokker-Planck equation below, $\mu_{j}$ are
the components of the drift vector and $D_{j}$ are the components of the
diffusion vector

$\displaystyle\frac{\partial u}{\partial
t}+\sum_{j=1}^{D}\frac{\partial}{\partial
x_{j}}(\mu_{j}(x_{1},...,x_{D})u)+\sum_{j=1}^{D}\frac{\partial^{2}}{\partial
x_{j}^{2}}(D_{j}(x_{1},...,x_{D})u)=0,\qquad\mu_{j},D_{ij}\in\mathbb{R}.$ (36)

The drift vector and diffusion coefficients $D_{j}$ only have real-valued
components, and $D_{j}>0$. Then it is clear

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=-\sum_{j=1}^{D}\hat{p}_{j}\mu_{j}(\hat{x}_{1},...,\hat{x}_{D})+i\sum_{j=1}^{D}\hat{p}^{2}_{j}D_{j}(\hat{x}_{1},...,\hat{x}_{D}).$
(37)

Applying Schrödingerisation we see

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger},$
$\displaystyle\boldsymbol{A}_{1}=\frac{1}{2}\sum_{j=1}^{D}(\hat{p}_{j}\mu_{j}+\mu_{j}\hat{p}_{j})+\frac{i}{2}\sum_{j=1}^{D}[\hat{p}^{2}_{j},D_{j}]$
$\displaystyle\boldsymbol{A}_{2}=-\frac{i}{2}\sum_{j=1}^{D}[\mu_{j},\hat{p}_{j}]-\frac{1}{2}\sum_{j=1}^{D}(\hat{p}_{j}^{2}D_{j}+D_{j}\hat{p}^{2}_{j})$
(38)

For example, $\mu_{j}=c_{j}\hat{x}_{j}$ with constant $c_{j}\in\mathbb{R}$ in
the case of linear drift. For additive noise $D_{j}=a_{j}$ where $a_{j}>0$ is
a positive constant and for multiplicative noise, in the case of geometric
Brownian motion, $D_{j}\sim d_{j}x_{j}^{2}$ with real constants $d_{j}$. In
the limit of linear drift and additive noise,
$\boldsymbol{A}=-\sum_{j=1}^{D}c_{j}\hat{p}_{j}\hat{x}_{j}-ia\sum_{j=1}^{D}\hat{p}_{j}^{2}$
and for Schrödingerisation we require the following Hamiltonian acting on
$D+1$ modes

$\displaystyle\boldsymbol{H}=\frac{1}{2}\sum_{j=1}^{D}\left(c_{j}I-a\hat{p}_{j}^{2}\right)\otimes\hat{\eta}+\frac{1}{2}\sum_{j=1}^{D}c_{j}(\hat{x}_{j}\hat{p_{j}}+\hat{p}_{j}\hat{x}_{j})\otimes
I.$ (39)

When we have linear drift with geometric Brownian motion, the Hamiltonian
becomes

$\displaystyle\boldsymbol{H}=-\frac{1}{2}\sum_{j=1}^{D}(c_{j}I-d_{j}(\hat{x}_{j}^{2}\hat{p}_{j}^{2}+\hat{p}_{j}^{2}\hat{x}^{2}_{j}))\otimes\hat{\eta}+\sum_{j=1}^{D}\left(\frac{c_{j}}{2}+d_{j}\right)(\hat{x}_{j}\hat{p}_{j}+\hat{p}_{j}\hat{x}_{j})\otimes
I$ (40)

where for both examples, we don’t require more than pairwise interactions.

####  V.2.3 The Black-Scholes equation

We can also illustrate our approach with the Black-Scholes equation [36]

$\displaystyle\frac{\partial u}{\partial
t}+\frac{1}{2}\sigma^{2}x^{2}\frac{\partial^{2}u}{\partial
x^{2}}+rx\frac{\partial u}{\partial x}-ru=0,\qquad\sigma,r\in\mathbb{R}.$ (41)

This is a PDE that evaluates the price of a financial derivative, where $r$
and $\sigma$ are constants. Then $\boldsymbol{u}$ satisfies

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=\frac{1}{2}\sigma^{2}\hat{x}^{2}\hat{p}^{2}+ir\hat{x}\hat{p}+rI.$
(42)

We can Schrödingerise the system

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger}$
$\displaystyle\boldsymbol{A}_{1}=\sigma^{2}\left(\frac{\hat{p}}{2}+\hat{x}\hat{p}\right),\qquad\boldsymbol{A}_{2}=-\frac{\sigma^{2}}{4}(\hat{x}^{2}\hat{p}^{2}+\hat{p}^{2}\hat{x}^{2})-\frac{r}{2}.$
(43)

Since this is a one-dimensional PDE evolving in time, only two qumodes are
necessary to simulate $|u\rangle$. The highest-order terms in this Hamiltonian
correspond to nonlinear gates of the form
$\exp(-i\hat{x}^{2}\hat{p}^{2}\otimes\hat{\eta}t)$.

####  V.2.4 The wave equation

Consider a $D$-dimensional wave equation of the form

$\displaystyle\frac{\partial^{2}u}{\partial
t^{2}}-\sum_{j=1}^{D}a_{j}(x_{1},...,x_{D})\frac{\partial^{2}u}{\partial
x_{j}^{2}}=-V(x_{1},...,x_{D})u,\qquad u(0,x),\qquad
x_{j}\in\mathbb{R}^{D},a_{j}\in\mathbb{R}^{+},V\in\mathbb{R}.$ (44)

This is just a special case of Eq. (25) with $a_{j,2}=a_{j}$, $b=V$ and all
the other coefficients in Eq. (25) are equal to zero. We can then directly use
the results in Section IV.2. Then we see that $\boldsymbol{u}$ obeys

$\displaystyle\frac{d^{2}\boldsymbol{u}}{dt^{2}}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=-i\left(\sum_{j=1}^{D}a_{j}(\hat{x}_{1},...,\hat{x}_{D})\hat{p}_{j}^{2}+V(\hat{x}_{1},...,\hat{x}_{D})\right)$
$\displaystyle\boldsymbol{A}_{1}=\frac{i}{2}\sum_{j=1}^{D}[\hat{p}^{2}_{j},a_{j}],\qquad\boldsymbol{A}_{2}=\frac{1}{2}\sum_{j=1}^{D}(a_{j}\hat{p}^{2}_{j}+\hat{p}_{j}^{2}a_{j})+V(\hat{x}_{1},...,\hat{x}_{D}).$
(45)

Following Section IV.2 we see that Schrödingerisation can be applied to the
dilated system
$\boldsymbol{u}\rightarrow\boldsymbol{y}=\boldsymbol{u}\otimes|0\rangle+(d\boldsymbol{u}/dt)\otimes|1\rangle$,
and we find the following system of $D+1$-qumodes and one qubit

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{V}_{2}\otimes\hat{\eta}+\boldsymbol{V}_{1}\otimes
I,$ (46)
$\displaystyle\boldsymbol{V}_{1}=\boldsymbol{A}_{1}\otimes\sigma_{x}+\left(\boldsymbol{A}_{2}-\frac{I}{2}\right)\otimes\sigma_{y},\qquad\boldsymbol{V}_{2}=-\boldsymbol{A}_{1}\otimes\sigma_{y}+\left(\boldsymbol{A}_{2}+\frac{I}{2}\right)\otimes\sigma_{x}.$
(47)

The scenarios where $V\neq 0$ could include equations like the Klein-Gordon
equation where we use $V=m^{2}I$, which is easily achieved by including this
into the $\boldsymbol{A}_{2}$ term. Then when $a_{j}$ is a constant
$\boldsymbol{A}_{1}=0$ and the highest-order term in $\boldsymbol{H}$ would be
of the form $\hat{p}^{2}\otimes\hat{\eta}\otimes\hat{\sigma}_{x,y}$, which is
a coupling between two qumodes and one qubit.

####  V.2.5 Maxwell’s equations

We can also consider not just a scalar PDE, but rather a system of linear PDEs
like Maxwell’s equations. Assume a linear medium with permittivity
$\varepsilon$ and permeability $\mu$. The speed of the electromagnetic wave in
this medium is then $v=1/\sqrt{\mu\varepsilon}$. We define in Eq. (V.2.5) the
quantities $\mathcal{F}$, $\mathcal{J}$ and $\mathcal{D}$, where $E_{x,y,z}$
and $B_{x,y,z}$ are the $x,y,z$ components of the electric and magnetic fields
respectively and $J_{x,y,z}$ are the $x,y,z$ components of the current density
and $\rho$ is the charge density. The vectors $\boldsymbol{F}$ in bold like
$\boldsymbol{E}_{x,y,z}$, $\boldsymbol{B}_{x,y,z}$, $\boldsymbol{J}_{x,y,z}$
and $\boldsymbol{\rho}$ represents $\boldsymbol{F}=\iiint
F(x,y,z)|x\rangle|y\rangle|z\rangle dxdydz$. Then

$\displaystyle\mathcal{F}=\frac{1}{\sqrt{2}}\begin{pmatrix}\sqrt{\varepsilon}E_{x}\\\
\sqrt{\varepsilon}E_{y}\\\ \sqrt{\varepsilon}E_{z}\\\ 0\\\
(\sqrt{\mu})^{-1}B_{x}\\\ (\sqrt{\mu})^{-1}B_{y}\\\ (\sqrt{\mu})^{-1}B_{z}\\\
0\end{pmatrix},\qquad\boldsymbol{\mathcal{F}}=\sqrt{\varepsilon}\boldsymbol{E}_{x}|0\rangle+\sqrt{\varepsilon}\boldsymbol{E}_{y}|1\rangle+\sqrt{\varepsilon}\boldsymbol{E}_{z}|2\rangle+(\sqrt{\mu})^{-1}\boldsymbol{B}_{x}|4\rangle+(\sqrt{\mu})^{-1}\boldsymbol{B}_{y}|5\rangle+(\sqrt{\mu})^{-1}\boldsymbol{B}_{z}|6\rangle,$
$\displaystyle\mathcal{J}=\frac{1}{\sqrt{2\varepsilon}}\begin{pmatrix}J_{x}\\\
J_{y}\\\ J_{z}\\\ 0\\\ 0\\\ 0\\\ 0\\\
-v\rho\end{pmatrix},\qquad\boldsymbol{\mathcal{J}}=\boldsymbol{J}_{x}|0\rangle+\boldsymbol{J}_{y}|1\rangle+\boldsymbol{J}_{z}|2\rangle-\boldsymbol{v}\boldsymbol{\rho}|8\rangle,$
$\displaystyle\mathcal{D}=\begin{pmatrix}0&-\partial_{z}&\partial_{y}&-\partial_{x}\\\
\partial_{z}&0&-\partial_{x}&-\partial_{y}\\\
-\partial_{y}&\partial_{x}&0&-\partial_{z}\\\
\partial_{x}&\partial_{y}&\partial_{z}&0\end{pmatrix}=\begin{pmatrix}0&i\hat{p}_{z}&-i\hat{p}_{y}&i\hat{p}_{x}\\\
-i\hat{p}_{z}&0&i\hat{p}_{x}&i\hat{p}_{y}\\\
i\hat{p}_{y}&-i\hat{p}_{x}&0&i\hat{p}_{z}\\\
-i\hat{p}_{x}&-i\hat{p}_{y}&-i\hat{p}_{z}&0\end{pmatrix}=-\hat{p}_{z}\otimes
I\otimes\sigma_{y}-\hat{p}_{x}\otimes\sigma_{y}\otimes\sigma_{x}+\hat{p}_{y}\otimes\sigma_{y}\otimes\sigma_{z}.$
(48)

We consider $\boldsymbol{\mathcal{F}}$ as a state of a 3 qumodes and 3 qubits.
Then Maxwell’s equations can be written as

$\frac{\partial\mathcal{F}}{\partial
t}=v\begin{pmatrix}\textbf{0}&\mathcal{D}(1-\bar{\mu})\\\
-\mathcal{D}(1-\bar{\varepsilon})&\textbf{0}\end{pmatrix}\mathcal{F}-\mathcal{J}=v\left(\frac{1}{2}\mathcal{D}(\bar{\varepsilon}-\bar{\mu})\otimes\sigma_{x}+\frac{i}{2}\mathcal{D}(2-\bar{\varepsilon}-\bar{\mu})\otimes\sigma_{y}\right)\mathcal{F}-\mathcal{J}$
(49)

where $\bar{\varepsilon}=\ln\varepsilon/2$, $\bar{\mu}=\ln\mu/2$. Then the
following equation governing $\boldsymbol{\mathcal{F}}$ with an inhomogeneous
term holds: $\boldsymbol{\mathcal{J}}$

$\displaystyle\frac{d\boldsymbol{\mathcal{F}}}{dt}=-i\boldsymbol{A}\boldsymbol{\mathcal{F}}-\boldsymbol{\mathcal{J}}$
$\displaystyle\boldsymbol{A}=\frac{v(\hat{x},\hat{y},\hat{z})}{2}\mathcal{D}(2I-\bar{\varepsilon}(\hat{x},\hat{y},\hat{z})-\bar{\mu}(\hat{x},\hat{y},\hat{z}))\otimes\sigma_{y}-i\frac{v(\hat{x},\hat{y},\hat{z})}{2}\mathcal{D}(\bar{\varepsilon}(\hat{x},\hat{y},\hat{z})-\bar{\mu}(\hat{x},\hat{y},\hat{z}))\otimes\sigma_{x}$
(50)

One can then apply our formulation to the inhomogeneous case in Section IV.1.
We Schrödingerise the system by dilating
$\boldsymbol{\mathcal{F}}\rightarrow\boldsymbol{y}=\boldsymbol{\mathcal{F}}\otimes|0\rangle+\boldsymbol{\mathcal{J}}\otimes|1\rangle$
and transform $\boldsymbol{y}\rightarrow\tilde{\boldsymbol{v}}$ to obtain

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\frac{1}{2}(I+\sigma_{z})\otimes\hat{\eta}-\frac{I}{2}\otimes\sigma_{x}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes\frac{1}{2}(I+\sigma_{z})\otimes
I+\frac{I}{2}\otimes\sigma_{y}\otimes I,$ (51)

where $\boldsymbol{A}_{1}=(\boldsymbol{A}+\boldsymbol{A}^{\dagger})/2$,
$\boldsymbol{A}_{2}=i(\boldsymbol{A}-\boldsymbol{A}^{\dagger})/2$. This is
quantum simulation on a system of 4 qumodes and 4 qubits.

##  VI PDEs with uncertainty

Often PDEs have coefficients which depend on $L$ stochastic variables, where
$L$ might be large. This defines an uncertain PDE, the study of which is the
subject of ”Uncertainty Quantification” [37]. We illustrate our idea in a
simple first-order convection equation with uncertain coefficients of the
following form

$\displaystyle\frac{\partial u}{\partial
t}-\sum_{j=1}^{D}c_{j}(z_{1},...,z_{L},x_{1},...,x_{D})\frac{\partial
u}{\partial x_{j}}=0,\qquad z_{1},...,z_{L}\sim
W(z_{1},...,z_{L})=\prod_{l=1}^{L}w(z_{l})=\prod_{l=1}^{L}\frac{\exp(-z_{l}^{2})}{\sqrt{\pi}},\qquad
c_{j}\in\mathbb{R}$ (52)

where $z_{1},...,z_{L}$ are random variables each sampled from the Gaussian
probability distribution function $w(z_{l})$. Typically these equations become
very expensive to solve in classical computers as $L$ increases, due to the
curse of dimensionality also with respect to $L$. In the Monte-Carlo or other
sample based methods such as the stochastic collocation method, one needs to
sample $z_{1},...,z_{D}$ multiple times and solve a new PDE – due to different
$c$ – for each sample.

Our goal is instead to simulate ensemble-averaged quantities of $u$, without
solving Eq. (52) multiple times, one for each sample of $z=(z_{1},...,z_{L})$.
We use the polynomial chaos based stochastic Galerkin approach [38] and call
our method the quantum stochastic Galerkin method. Let
$\boldsymbol{n}=(n_{1},...,n_{L})\in\mathbb{N}_{0}^{L}$ be a multi-index. Then
any function $u$ (in $l_{2}$ space with respect to variables
$z_{1},...,z_{D}$) can be expressed as

$\displaystyle
u(t,z_{1},...,z_{L},x_{1},...,x_{D})=\sum_{\boldsymbol{n}\in\mathbb{N}_{0}^{L}}u_{\boldsymbol{n}}(t,x_{1},\cdots,x_{D})P_{\boldsymbol{n}}(z_{1},...,z_{L}).$
(53)

where

$\displaystyle u_{\boldsymbol{n}}(t,x_{1},...,x_{D})=\int
u(t,z_{1},...,z_{L},x_{1},...,x_{D})P_{\boldsymbol{n}}(z_{1},...,z_{L})W(z_{1},...,z_{L})dz_{1}...dz_{L}$
(54)

with the normalisation

$\displaystyle\int
P_{\boldsymbol{n}}(z_{1},...,z_{L})P_{\boldsymbol{m}}(z_{1},...,z_{L})W(z_{1},...,z_{L})dz_{1}...dz_{L}=\delta_{\boldsymbol{n}\boldsymbol{m}}.$
(55)

$P_{\boldsymbol{n}}$ is just the normalized multivariate Hermite polynomials
of degree $|\boldsymbol{n}|=n_{1}+\cdots+n_{L}$. Since
$\boldsymbol{n}\in\mathbb{N}_{0}^{L}$, one can identify $\boldsymbol{n}$ to
particle numbers, so we can define
$\\{|\boldsymbol{n}\rangle\\}_{\boldsymbol{n}\in\mathbb{N}_{0}^{L}}$ as the
Fock or number basis, so
$\sum_{\boldsymbol{n}\in\mathbb{N}_{0}^{L}}|\boldsymbol{n}\rangle\langle\boldsymbol{n}|=I$.
Since $z_{1},...,z_{D}\in\mathbb{R}$, we can associate $z$ with a quadrature
value, like position eigenvalues. Thus $\int|z_{1},...,z_{L}\rangle\langle
z_{1},...,z_{L}|dz_{1}...dz_{L}=I$. This implies
$\int\langle\boldsymbol{n}|z_{1},...,z_{L}\rangle\langle
z_{1},...,z_{L}|\boldsymbol{m}\rangle
dz_{1}...dz_{L}=\delta_{\boldsymbol{n}\boldsymbol{m}}$. Combining this with
the normalisation in Eq. (55) and $W=\prod_{l=1}^{L}w(z_{l})$ in Eq. (52)
allows one to identify

$\displaystyle P_{\boldsymbol{n}}(z)=\frac{\langle
z_{1}...z_{L}|\boldsymbol{n}\rangle}{W^{1/2}(z_{1},...,z_{L})}=\prod_{l=1}^{L}P_{n_{l}}=\prod_{l=1}^{L}\frac{\langle
z_{l}|n_{l}\rangle}{w^{1/2}(z_{l})}=\prod_{l=1}^{L}\frac{H_{n_{l}}(z_{l})}{\sqrt{2^{n_{l}}n_{l}!}},\qquad
w(z_{l})=\frac{\exp(-z_{l}^{2})}{\sqrt{\pi}},$ (56)

where $H_{n_{l}}(z_{l})$ are the one-dimensional Hermite polynomials and
$P_{n_{l}}=H_{n_{l}}(z_{l})/\sqrt{2^{n_{l}}n_{l}!}$. The last equality is
well-known in quantum optics when $|n_{l}\rangle$ is the Fock state and
$|z_{l}\rangle$ is the position eigenstate [39].

To recover the ensemble-averaged quantities of $u$, we observe that
$P_{0}(z_{l})=H_{0}(z_{l})=1$. From the expansion in Eq.(54), one can recover
the expectation value of $u$

$\displaystyle u_{\boldsymbol{0}}(t,x_{1},...,x_{D})=\int
u(t,z_{1},...,z_{L},x_{1},...,x_{D})w(z_{1})...w(z_{L})dz_{1}...dz_{L}\equiv\mathbb{E}_{z_{1},...,z_{L}}(u).$
(57)

Similarly, to recover the variance of $u$, using the expansion of $u$ in Eq.
(53) and the normalisation condition in Eq. (55), it is simple to see that

$\displaystyle\text{Var}_{z}(u)\equiv\int
u^{2}(t,z_{1},...,z_{L},x_{1},...,x_{D})w(z_{1})...w(z_{l})dz_{1}...dz_{L}-\mathbb{E}_{z_{1},...,z_{L}}^{2}(u)=\sum_{\boldsymbol{n}\in\mathbb{N}_{0}^{L},|n|\geq
1}u_{\boldsymbol{n}}^{2}(t,x_{1},...,x_{D}).$ (58)

This means it is of more physical relevance to solve for
$u_{\boldsymbol{n}}(t,x_{1},...,x_{D})$ than the individual
$u(t,z_{1},...,z_{L},x_{1},...,x_{D})$ values for every random sample of $z$.

If $u(t,x_{1},...,x_{D},z_{1},...,z_{L})$ solves the PDE in Eq. (52), then
using the expansion in Eq. (53) it is simple to show

$\displaystyle\frac{\partial}{\partial
t}u_{\boldsymbol{n}}(t,x_{1},...,x_{D})-\sum_{j=1}^{D}\sum_{\boldsymbol{m}\in\mathbb{N}_{0}^{L}}\Lambda_{j\boldsymbol{n}\boldsymbol{m}}(x_{1},...,x_{D})\frac{\partial}{\partial
x_{j}}u_{\boldsymbol{n}}(t,x_{1},...,x_{D})=0,$
$\displaystyle\Lambda_{j\boldsymbol{n}\boldsymbol{m}}(x_{1},...,x_{D})=\int
c_{j}(z_{1},...,z_{L},x_{1},...,x_{D})w(z_{1})...w(z_{L})P_{\boldsymbol{n}}(z_{1},...,z_{L})P_{\boldsymbol{m}}(z_{1},...,z_{L})dz_{1},...,dz_{L}$
(59)

where the coefficients of the PDE are now completely independent of $z$! From
definition of $\Lambda_{j\boldsymbol{n}\boldsymbol{m}}$ and Eq. (56), one can
rewrite

$\displaystyle\Lambda_{j\boldsymbol{n}\boldsymbol{m}}=\int
c_{j}(z_{1},...,z_{L},x_{1},...,x_{D})\langle\boldsymbol{m}|z_{1},...,z_{L}\rangle\langle
z_{1},...,z_{L}|\boldsymbol{n}\rangle dz_{1},...,dz_{L}$
$\displaystyle=\langle\boldsymbol{m}|\int
c_{j}(x_{1},...,x_{D},\hat{z}_{1},...,\hat{z}_{L})|z_{1},...,z_{L}\rangle\langle
z_{1},...,z_{L}|\boldsymbol{n}\rangle
dz_{1}...dz_{L}=\langle\boldsymbol{m}|c_{j}(x_{1},...,x_{D},\hat{z}_{1},...,\hat{z}_{L})|\boldsymbol{n}\rangle.$
(60)

where $\hat{z}_{l}|z_{l}\rangle=z_{l}|z_{l}\rangle$ is the position operator.

Thus instead of solving the $z$-dependent equation multiple times in Eq. (52)
for every sample of $z_{1},...,z_{L}$, we can solve the $z$-independent
equation in Eq. (VI). We can define $\boldsymbol{u}(t)$ as the infinite-
dimensional vector in both the Fock and the $x$ bases

$\displaystyle\boldsymbol{u}(t)\equiv\int
dx\sum_{\boldsymbol{n}\in\mathbb{Z}_{0}^{L}}u_{\boldsymbol{n}}(t,x_{1},...,x_{D})|\boldsymbol{n}\rangle|x_{1},...,x_{D}\rangle$
(61)

and the corresponding $(d+L)$-qumode quantum state
$|u(t)\rangle=(1/\mathcal{N}_{u})\boldsymbol{u}_{\boldsymbol{n}}(t)$ where we
assume the normalisation $\mathcal{N}^{2}_{u}=\int
dx\sum_{\boldsymbol{n}\in\mathbb{N}_{0}^{L}}|u_{\boldsymbol{n}}(t,x_{1},...,x_{D})|^{2}<\infty$.

Then it is simple to see from Eq. (VI) and Eq. (VI) that $\boldsymbol{u}$
obeys

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=\sum_{j=1}^{D}c_{j}(\hat{z}_{1},...,\hat{z}_{L},\hat{x}_{1},...,\hat{x}_{D})\hat{p}_{j},\qquad\boldsymbol{u}(0)=\int\sum_{\boldsymbol{n}\in\mathbb{N}_{0}^{L}}u_{\boldsymbol{n}}(0,x_{1},...,x_{D})|\boldsymbol{n}\rangle|x_{1},...,x_{D}\rangle
dx_{1}...dx_{D}.$ (62)

We can easily Schrödingerise this equation

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I,\qquad\tilde{\boldsymbol{v}}(0)=|\Xi\rangle\boldsymbol{u}(0),$
$\displaystyle\boldsymbol{A}_{1}=\frac{1}{2}\sum_{j=1}^{D}(c_{j}(\hat{z}_{1},...,\hat{z}_{L},\hat{x}_{1},...,\hat{x}_{D})\hat{p}_{j}+\hat{p}_{j}c_{j}(\hat{z}_{1},...,\hat{z}_{L},\hat{x}_{1},...,\hat{x}_{D}))$
$\displaystyle\boldsymbol{A}_{2}=\frac{i}{2}\sum_{j=1}^{D}[c_{j}(\hat{z}_{1},...,\hat{z}_{L},\hat{x}_{1},...,\hat{x}_{D}),\hat{p}_{j}].$
(63)

Note that when $c_{j}$ doesn’t depend on $\hat{x}_{j}$, then
$[c_{j},\hat{p}_{j}]=0=\boldsymbol{A}_{2}$, thus
$\boldsymbol{A}=\boldsymbol{A}^{\dagger}$ and we can directly use Hamiltonian
simulation and no Schrödingerisation is needed.

We emphasise here that although $\boldsymbol{A}_{1},\boldsymbol{A}_{2}$
depends on the operators $\hat{z}_{1},...,\hat{z}_{L}$, they are only
considered as position operators – not random variables – that act on the
number basis $|\boldsymbol{n}\rangle$. Eq. (VI) is therefore independent of
random coefficients and
$c_{j}(\hat{z}_{1},...\hat{z}_{L},\hat{x}_{1},...,\hat{x}_{D})$ is a
deterministic operator. However, randomness and information about the
distribution $w$ do enter, but only through the initial conditions. To prepare
the initial state, we first rewrite the initial state in the
$|z_{1},...,z_{L}\rangle,|x_{1},...,x_{D}\rangle$ bases

$\displaystyle\boldsymbol{u}(0)=\iint
u(0,z_{1},...,z_{L},x_{1},...,x_{D})w^{1/2}(z_{1})...w^{1/2}(z_{L})|z_{1},...,z_{L}\rangle|x_{1},...,x_{D}\rangle
dz_{1}...dz_{L}dx_{1}...dx_{D},\qquad
w(z_{l})=\frac{e^{-z_{l}^{2}}}{\sqrt{\pi}}.$ (64)

This means that all the randomness of the PDE only resides in the initial
state preparation with the inclusion of $w^{1/2}(z_{l})$, without needing to
solve the PDE in Eq. (52) multiple times to estimate ensemble averaged values
like the expectation and variance of $u$. For example, if the $x$ and $z$
dependence in the initial state $u(0,z_{1},...,z_{L},x_{1},...,x_{D})$ are all
Gaussian, then $|u(0)\rangle$ is a Gaussian state since $w$ is also Gaussian.
This makes the initial state efficient to prepare. Furthermore, if this
initial condition holds, this is the unique initial state to prepare for
solving any uncertain linear PDE, and is not dependent on any particular
features of the original PDE itself.

Another interesting observation is that here it is the probability amplitude
of the distribution of $z_{l}$, i.e., $w^{1/2}(z_{l})$ rather than $w(z_{l})$,
that is embedded in the amplitude of the pure quantum state $|u(0)\rangle$,
which is a uniquely quantum way of embedding the probability distribution.
This is a continuous-variable version of the qubit-based Grover-Rudolph
embedding [40], which is also considered in quantum walks. A classical
sampling of $z_{l}$, on the other hand, would be equivalent to preparing a
density matrix, instead of a pure state, in $z_{l}$. We leave the discussion
of the significance of this embedding for now [41].

Given the output $\boldsymbol{u}(t)$, if one can measure in the number basis
$|\boldsymbol{n}\rangle\langle\boldsymbol{n}|$, one is able to retrieve a
state $\propto\int
u_{\boldsymbol{n}}(t,x_{1},...,x_{D})|x_{1},...,x_{D}\rangle dx_{1}...dx_{D}$,
whose amplitudes correspond to already ensemble-averaged physically important
quantities. The choice of $\boldsymbol{n}=\boldsymbol{0}$ for example can give
the state whose amplitudes approximate $\mathbb{E}_{z_{1},...,z_{L}}(u)$.

Typically, for smooth, say $C^{\infty}$, functions $u$, the coefficients in
the above expansion in $n_{l}$ decays exponentially in $n_{l}$, so only small
$n_{l}$ values need to be kept in the expansion, which can be used to retrieve
accurately physically meaningful quantities. Suppose we keep terms in the
expansion to $n_{l}=n_{max}$ for each of the $L$ modes. Then it is known that
the error in approximating $u$ with an expansion in $u_{n_{l}}$ scales like
$\sim 1/n_{max}^{m}$, where $m$ is the regularity of $u$, i.e., all
derivatives up to $\partial^{m}u/\partial z_{l}^{m}$ exist and are bounded,
but no higher derivatives exist. Thus, to capture most of the ensemble
information about $u$, we only need to meausure $n_{l}$ up to $n_{max}$ which
can be small, typically $n_{max}\sim 3-5$ will be enough, for smooth functions
$u$.

For a simple example, it is known that the following approximation is often
used in the Uncertainty Quantification community (either through the Karhunen-
Loeve expansion, or the stochastic Galerkin approximation [42, 38])

$\displaystyle
c_{j}(z_{1},...,z_{L},x_{1},...,x_{D})\sim\sum_{l=1}^{L}a_{jl}(z_{1},...,z_{L})b_{jl}(x_{1},...,x_{D}).$
(65)

Often the Karhunen-Loeve expansion is used and we can assume the linear
approximation

$\displaystyle
c_{j}(z_{1},...,z_{L},x_{1},...,x_{D})=c_{j1}(x_{1},...,x_{D})+\sum_{l=1}^{L}c_{j2l}(x_{1},...,x_{D})z_{l},\qquad
c_{j1},c_{j2l}\in\mathbb{R}.$ (66)

If $c_{j1}(x_{1},...,x_{D})=\sum_{k=1}^{D}c_{j1k}x_{k}$,
$c_{j2l}(x_{1},...,x_{D})=\sum_{k=1}^{D}c_{j2lk}x_{k}$ are also linear with
constants $c_{j1k},c_{j2lk}$, this means one can consider Hamiltonians of the
general form

$\displaystyle\boldsymbol{H}=\frac{1}{2}\sum_{j,k=1}^{D}\left(c_{j1k}+\sum_{l=1}^{L}c_{j2lk}\hat{z}_{l}\right)(\hat{x}_{k}\hat{p}_{j}+\hat{p}_{j}\hat{x}_{k})\otimes
I-\frac{1}{2}\sum_{j=1}^{D}\left(c_{j1j}+\sum_{l=1}^{L}c_{j2lj}\hat{z}_{l}\right)\otimes\hat{\eta}.$
(67)

Finally, we remark how continuous-variable quantum systems are naturally
suited for this problem through the use of the continuous-variable number and
position eigenstates. The inner product of these naturally output the
normalised Hermite polynomials. An analogous quantum algorithm with qubits to
approximate the effect of the normalised Hermite polynomials, on the other
hand, would be very complicated.

Other classes of linear PDEs and distributions can also be considered
following a similar approach and using Schrödingerisation. For more details
and more general PDEs, see our upcoming paper on quantum stochastic Galerkin
methods [41].

##  VII Nonlinear PDEs

The most natural way to simulate nonlinear PDEs with quantum simulation is to
find a way to represent the nonlinear problem in a linear way, since quantum
mechanics is fundamentally linear. We must distinguish two types of approaches
that converts a nonlinear PDE into a linear PDE. One class of methods
necessarily involve approximations (e.g. either through linearisation of the
nonlinearity or through discretisation). They include the linear approximation
(e.g. Carlemann linearisation) and the linear representation of nonlinear ODEs
methods (e.g. through Koopman-von Neumann), where in the latter case the
nonlinear ODE is approximated as a finite system of nonlinear ODEs. Typically,
these methods are not useful beyond the regime of weak nonlinearity.

The second approach is the linear representation for nonlinear PDEs, where no
approximations are required to map between the nonlinear and linear PDEs. For
our continuous-variable approach, we choose this second method of applying an
exact formulation without approximations.

While this second approach is the most desirable, it is not simple to find
such an exact mapping for all nonlinear PDEs, or know if such mappings even
exist. However, we do know some important classes, which include nonlinear
Hamilton-Jacobi PDEs (Hamiltonian system) and scalar hyperbolic PDEs (non-
Hamiltonian system). Quantum algorithms in the qubit setting and without using
Schrödingerisation were found for these PDEs in [8, 9]. We now employ our
continuous-variable formalism based on Schrödingerisation to simulate these
nonlinear PDEs.

###  VII.1 The level-set formalism

We briefly describe the level set formalism for the nonlinear scalar
hyperbolic and nonlinear Hamilton-Jacobi PDEs. The reader can refer to the
references [8, 9] for more details on the methodology in the qubit setting. We
begin with a $(D+1)$-dimensional scalar nonlinear hyperbolic PDE with $M$
initial conditions

$\displaystyle\frac{\partial u^{[k]}}{\partial
t}-\sum_{j=1}^{D}F_{j}(u^{[k]})\frac{\partial u^{[k]}}{\partial
x_{j}}+Q(x_{1},...,x_{D},u^{[k]})=0,\qquad u^{[k]}(0,x)=u^{[k]}_{0}(x),\quad
k=1,...,M$ (68)

We introduce a level set function $\phi(t,x_{1},...,x_{D},\chi)$, where
$\chi\in\mathbb{R}$. It is defined so its zero level set is the solution $u$:

$\phi^{[k]}(t,x_{1},...,x_{D},\chi)=0\quad{\text{at}}\quad\chi=u^{[k]}(t,x_{1},...,x_{D})$
(69)

Then defining $X=(x_{1},...,x_{D},\chi)$, $a_{j}=F_{j}(\chi)$ for $j=1,...,D$
and $a_{D+1}=-Q(x_{1},...,x_{D},\chi)$

$\displaystyle\Psi(t,x_{1},...,x_{D},\chi)=\frac{1}{M}\sum_{k=1}^{M}\delta(\phi^{[k]}(t,x_{1},...,x_{D},\chi))$
(70)

Then $\Psi$ evolves according to the linear PDE in one extra dimension

$\displaystyle\frac{\partial\Psi^{[k]}}{\partial
t}-\sum_{j=1}^{D+1}\frac{\partial}{\partial X_{j}}(a_{j}(X)\Psi^{[k]})=0$ (71)
$\displaystyle\Psi^{[k]}(0,x_{1},...,x_{D},\chi)=\frac{1}{M}\sum_{k=1}^{M}\delta(\chi-u^{[k]}_{0}(x_{1},...,x_{D})).$
(72)

If $\boldsymbol{\Psi}(t)\equiv\int dX\Psi(t,X)|X\rangle$, then it satisfies
the following linear ODE

$\displaystyle\frac{d\boldsymbol{\Psi}}{dt}=-i\boldsymbol{A}\boldsymbol{\Psi},\qquad\boldsymbol{A}=\sum_{j=1}^{D+1}\hat{P}_{j}a_{j}(\hat{X})=\sum_{j=1}^{D}\hat{p}_{j}F_{j}(\hat{\chi})-\hat{\zeta}Q(\hat{x}_{1},...,\hat{x}_{D},\hat{\chi})$
(73)

where $-i\hat{P}_{j}=\partial/\partial X_{j}$,
$-i\hat{p}_{j}=\partial/\partial x_{j}$,
$-i\hat{\zeta}=\partial/\partial\chi$. Since $\hat{p}_{j}$ and $\hat{\chi}$
operate on different modes and if $F_{j}(u^{[k]})$ is a real-valued function,
then the first term in $\boldsymbol{A}$ is hermitian. Thus when $Q=0$ and
$F\in\mathbb{R}$, one has $\boldsymbol{A}=\boldsymbol{A}^{\dagger}$, so one
can directly apply Hamiltonian simulation to Eq. (73). Futhermore, one only
needs at most pairwise entangling quantum operations of the form
$F(\hat{\eta})\hat{p}$. The nonlinearity of the original PDE in Eq. (68) is in
the function $F_{j}$, so we see that due to the term
$F_{j}(\hat{\chi})\hat{p}_{j}$ in $\boldsymbol{H}$, then when $F(u)\propto u$,
two-mode Gaussian operations are in fact sufficient. For higher order
nonlinearity in $F(u)$, we need to go beyond Gaussian operations to simulate
any nonlinearity.

However, since $\hat{\chi}$ and $\hat{\zeta}$ operate on the same mode, then
the second term depending on $Q$ in $\boldsymbol{A}$ is not necessarily
hermitian. To deal with general $Q\neq 0$, we can apply Schrödingerisation in
Section IV to this system $\boldsymbol{\Psi}\rightarrow\tilde{\boldsymbol{v}}$
to get

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger}$
$\displaystyle\boldsymbol{A}_{1}=\sum_{j=1}^{D}\hat{p}_{j}F_{j}(\hat{\chi})-\frac{1}{2}\sum_{j=1}^{D}(\hat{\zeta}Q(\hat{x}_{1},...,\hat{x}_{D})+Q(\hat{x}_{1},...,\hat{x}_{D})\hat{\zeta}),\qquad\boldsymbol{A}_{2}=\frac{i}{2}\sum_{j=1}^{D}[Q(\hat{x}_{1},...,\hat{x}_{D}),\hat{\zeta}].$
(74)

This Hamiltonian operators on $D+1$ modes.

The process for the nonlinear Hamilton-Jacobi equation proceeds similarly. For
simplicity, we choose $M=1$. If $S$ obeys the nonlinear Hamilton-Jacobi PDE
and $u=\gradient S$, then $u$ solves the nonlinear hyperbolic system of
conservation equations in gradient form

$\displaystyle\frac{\partial u_{j}}{\partial t}+\frac{\partial}{\partial
x_{j}}H(x_{1},...,x_{D},u_{1},...,u_{D})=0,\qquad\forall j=1,...,D,\qquad
H\in\mathbb{R}$ (75)

where the function $H=H(x_{1},...,x_{D},p_{1},...,p_{2})$ is a Hamiltonian. We
now introduce $D$ variables $\chi_{1},...,\chi_{D}\in\mathbb{R}$. Then one can
define a system of level set functions $\phi=(\phi_{1},...,\phi_{D})$ where
each $\phi_{j}(t,x_{1},...,x_{D},\chi_{1},...,\chi_{D})=0$ at
$\chi_{j}=u_{j}$. The system $\phi$ solves the linear Liouville equation

$\displaystyle\frac{\partial\phi}{\partial t}+\sum_{j=1}^{D}\frac{\partial
H}{\partial\chi_{j}}\frac{\partial\phi}{\partial
x_{j}}-\sum_{j=1}^{D}\frac{\partial H}{\partial
x_{j}}\frac{\partial\phi}{\partial\chi_{j}}=0$ (76)

and the initial data can be chosen as $\phi_{j}(t=0)=\chi_{j}-u_{j}(t=0)$. We
can similarly define

$\displaystyle\Psi(t,x_{1},...,x_{D},\chi_{1},...,\chi_{D})=\prod_{j=1}^{D}\delta(\phi_{j}(t,x_{1},...,x_{D},\chi_{1},...,\chi_{D})),\qquad\Psi(t=0)=\prod_{j=1}^{D}\delta(\chi_{j}-u_{j}(t=0)).$
(77)

which obeys the following linear transport PDE

$\displaystyle\frac{\partial\Psi}{\partial t}+\sum_{j=1}^{D}\frac{\partial
H}{\partial\chi_{j}}\frac{\partial\Psi}{\partial
x_{j}}-\sum_{j=1}^{D}\frac{\partial H}{\partial
x_{j}}\frac{\partial\Psi}{\partial\chi_{j}}=0.$ (78)

Then

$\displaystyle\frac{d\boldsymbol{\Psi}}{dt}=-i\boldsymbol{H}\boldsymbol{\Psi},\qquad\boldsymbol{H}=i\sum_{j=1}^{D}[\hat{\zeta}_{j}H(\hat{x}_{1},...,\hat{x}_{D},\hat{p}_{1},...,\hat{p}_{D}),\hat{p}_{j}]=\boldsymbol{H}^{\dagger}.$
(79)

which acts on $2D+1$ qumodes. $\boldsymbol{H}=\boldsymbol{H}^{\dagger}$ since
$H$ is interpreted as a Hamiltonian and each of the quadratures is hermitian.
Here we can directly use Hamiltonian simulation and no Schrödingerisation is
actually required.

We remark that, since we are solving the Liouville equation, the solution
obtained by the level set formalism generates the multivalued solutions [43,
44, 45, 46], not the viscosity solution [47, 48], for scalar nonlinear
hyperbolic PDEs and Hamilton-Jacobi equations.

###  VII.2 Systems of nonlinear ODEs

Another way of representing nonlinear PDEs exactly is as a infinite system of
nonlinear ODEs. Discretisation methods lead to the approximation of solutions
of nonlinear PDEs by a finite system of $N$ nonlinear ODEs. Linearisation
methods also lead to solving only a finite number of nonlinear ODEs.
Discretisation methods fall largely into two classes: Eulerian and Lagrangian.
For example, grid-based Eulerian methods solve the PDEs on a fixed grid while
Lagrangian methods are mesh-free methods, like particle methods. Linearisation
methods include moment-closure methods and Carlemann linearisation. $N$ can be
larger or smaller depending on the discretisation or linearisation technique.

This means we first need to know how to simulate a system of $N$ nonlinear
ODEs for $\gamma_{n}(t)$ where each $\gamma_{n}(t)$ obeys a nonlinear ODE

$\displaystyle\frac{d\gamma_{n}(t)}{dt}=-F_{n}(\gamma_{0},...,\gamma_{N-1},t),\qquad
n=0,...,N-1$ (80)

and $F_{n}$ is a nonlinear function of its arguments. We can use a similar
mapping as the level set mapping to solve nonlinear ODEs [8, 9]. This is also
related to the Koopman-von Neumann methods, see [9] for a discussion on the
relationship. Here we introduce $N$ auxiliary variables $q_{0},...,q_{N-1}$
and a function $\Phi(t,q_{0},...,q_{N-1})$ defined by

$\displaystyle\Phi(t,q_{0},...,q_{N-1})=\prod_{n=0}^{N-1}\delta(q_{n}-\gamma_{n}(t)),\qquad
q_{n}\in\mathbb{R}.$ (81)

Then it is simple to check that $\Phi(t,q_{0},...,q_{N-1})$ satisfies, in the
weak sense [8], the linear $N+1$-dimensional PDE

$\displaystyle\frac{\partial\Phi(t,q_{0},...,q_{N-1})}{\partial
t}-\sum_{n=0}^{N-1}\frac{\partial}{\partial
q_{n}}(F_{n}(q_{0},...,q_{N-1})\Phi(t,q_{0},...,q_{N-1}))=0.$ (82)

We can treat this linear PDE with the same methods as before. We can define
$\hat{q}_{n}$ to be a quadrature operator whose eigenstate is $|q_{n}\rangle$
and its conjugate quadrature operator as $\hat{Q}_{n}$. Thus
$[\hat{q}_{n},\hat{Q}_{n}]=i$ and we can use the correspondence
$-i\hat{Q}_{n}\leftrightarrow\partial/\partial q_{n}$. Then defining
$\boldsymbol{\Phi}(t)\equiv\int\Phi(t,q_{0},...,q_{N-1})|q_{0},...,q_{N-1}\rangle
dq_{0}...dq_{N-1}$, we see that Eq. (82) becomes

$\displaystyle\frac{d\boldsymbol{\Phi}(t)}{dt}=-i\boldsymbol{A}\boldsymbol{\Phi}(t),\qquad\boldsymbol{A}=\sum_{n=0}^{N-1}\hat{Q}_{n}F_{n}(\hat{q}_{0},...,\hat{q}_{N-1}).$
(83)

We can easily use the same Schrödingerisation methods in Section (IV)
$\boldsymbol{\Phi}\rightarrow\tilde{\boldsymbol{v}}$ and

$\displaystyle\frac{d\tilde{\boldsymbol{v}}}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\boldsymbol{A}_{2}\otimes\hat{\eta}+\boldsymbol{A}_{1}\otimes
I=\boldsymbol{H}^{\dagger}$
$\displaystyle\boldsymbol{A}_{1}=\frac{1}{2}\sum_{n=0}^{N-1}(\hat{Q}_{n}F_{n}(\hat{q}_{0},...,\hat{q}_{N-1})+F_{n}(\hat{q}_{0},...,\hat{q}_{N-1})\hat{Q}_{n})$
$\displaystyle\boldsymbol{A}_{2}=\frac{i}{2}\sum_{n=0}^{N-1}[\hat{Q}_{n},F_{n}(\hat{q}_{0},...,\hat{q}_{N-1})]$
(84)

which is a Hamiltonian simulation with $N+1$ qumodes. Given that $F_{n}$ must
be a nonlinear function of its arguments, from the form of $\boldsymbol{A}$,
it is clear we must necessarily go beyond Gaussian operations, and the order
of nonlinearity of the Hamiltonian is either one or two orders larger
(depending on if Schrödingerisation is needed) than the order of the
nonlinearity of the original ODE system.

Table 2: Summary of resources for simulating PDEs with qumodes, Below $D$ denotes the spatial dimension of the PDE. $K$ is the maximum order of the spatial derivative of the PDE and $L$ denotes the maximum number of stochastic variables in an uncertain PDE. The Hamiltonian $\boldsymbol{H}$ refers to the Hamiltonian used in Schrödingerisation in the most general scenario. Quantities like $a,b,L,\text{D},V,\mu,c,Q,H,F$ are $\hat{x}$-dependent operators. They represent the $x$-dependent coefficients of the PDEs. The coefficients $c$ also depends on $\hat{z}$. Here $\hat{\eta}$, $\hat{\zeta}$ can be any quadrature, and $\hat{x}$, $\hat{p}$ can also be any quadrature with the requirement $[\hat{x},\hat{p}]=i$. Below $\\{A,B\\}\equiv AB+BA$ is the Poisson bracket and $\sigma$ refers to the Pauli matrices. The exact number of terms in $\boldsymbol{H}$ in the table below depends on the form of the coefficients of the PDE and must be worked out on an individual basis. If these coefficients can have an expansion with up to $G$ terms in $x$, then there could be $G$ times more terms in $\boldsymbol{H}$.  Equation | Number of qumodes | Max order term in $\boldsymbol{H}$ | Number of terms in $\boldsymbol{H}$
---|---|---|---
General linear PDEs |  |  |
Homogeneous IV | $\displaystyle D+1\par$ | $\displaystyle[a_{K,j},\hat{p}_{j}^{K}]\otimes\hat{\eta},\qquad b\otimes\hat{\eta}$ | $\displaystyle O(DK)$
$1^{\text{st}}$-order time derivative |  |  |
Inhomogeneous IV.1 | Add 1 more qubit | $\displaystyle[a_{K,j},\hat{p}_{j}^{K}]\otimes\hat{\eta}\otimes\sigma,\qquad b\otimes\hat{\eta}\otimes\sigma$ | $\displaystyle O(DK)$
$n$-th order time derivative IV.2 | $\displaystyle\quad\text{Add}\log_{2}(n)\text{ more qubits}$ | $\displaystyle\text{Similar terms to $1^{\text{st}}$-order}\otimes\sigma^{\otimes\log_{2}n}$ | $\displaystyle O(DK)$
Examples |  |  |
Liouville V.1 | $\displaystyle D+1$ | $\displaystyle i[\hat{p},L]\otimes\hat{\eta},\qquad\\{\hat{p},L\\}\otimes I$ | $\displaystyle O(D)$
Heat V.2.1 | $\displaystyle D+1$ | $\displaystyle\\{\hat{p}\text{D},\hat{p}\\}\otimes\hat{\eta},\qquad V\otimes\hat{\eta}$ | $\displaystyle O(D+1)$
Fokker-Planck V.2.2 | $\displaystyle D+1\par$ | $\displaystyle\\{\hat{p}^{2},\text{D}\\}\otimes\hat{\eta},i[\mu,\hat{p}]\otimes\hat{\eta},\\{\hat{p},\mu\\}\otimes I,i[\hat{p}^{2},\text{D}]\otimes I$ | $\displaystyle O(D)$
Black-Scholes V.2.3 | $\displaystyle 2$ | $\displaystyle(\hat{x}^{2}\hat{p}^{2}+\hat{p}^{2}\hat{x}^{2})\otimes\hat{\eta}$ | $\displaystyle 5$
Wave V.2.4 | $\displaystyle D+1+\text{one qubit}$ | $\displaystyle\\{a,\hat{p}^{2}\\}\otimes\hat{\eta}\otimes\sigma,\quad V\otimes\eta\otimes\sigma,\quad i[a,\hat{p}^{2}]\otimes I\otimes\sigma$ | $\displaystyle O(D)$
Maxwell ($D=3$) V.2.5 | $\displaystyle 4\text{ qumodes}+4\text{ qubits}$ | See Section V.2.5 | $\displaystyle 30$
Uncertain linear PDE VI |  |  |
Convection | $\displaystyle D+L+1$ | $\displaystyle i[c,\hat{p}]\otimes\hat{\eta},\qquad\\{c,\hat{p}\\}\otimes I$ | $\displaystyle O(DL)$
Nonlinear |  |  |
Scalar hyperbolic VII.1 | $\displaystyle D+1$ | $\displaystyle i[Q,\hat{\zeta}]\otimes\hat{\eta},\qquad\\{Q,\hat{\zeta}\\}\otimes I$ | $\displaystyle O(D)$
Hamilton-Jacobi VII.1 | $\displaystyle 2D$ | $\displaystyle[\hat{\zeta}H,\hat{p}]$ | $\displaystyle O(D)$
$N$ nonlinear ODEs VII.2 | $\displaystyle N+1$ | $\displaystyle i[Q,F]\otimes\hat{\eta},\qquad\\{Q,F\\}\otimes I$ | $\displaystyle O(N)$

##  VIII Extended Schrödingerisation

In the previous sections using Schrödingerisation, we only required a single
auxiliary qumode with quadrature operators $\hat{\eta},\hat{\xi}$. Thus a
linear PDE with $D$ spatial dimensions is simulable on a $D+1$-qumode quantum
device. This is the smallest number of auxiliary qumodes required, so from the
standpoint of complexity for universal analog quantum simulation, this is the
most favourable. However, it may not be suitable for direct simulation on some
quantum systems for large $D$. This is because the Schrödingerisation in
Section IV in general requires $O(D)$ couplings from other modes onto this
single mode. For analog quantum systems where multiple couplings to one
specific mode is difficult, we can instead modify our Schrödingerisation
procedure so we have instead $O(D)$ coupling terms with $D$ distinct auxiliary
modes. We call this the extended Schrödingerisation procedure.

Here, instead of introducing a single $\xi\in\mathbb{R}$, we introduce $D$
auxiliary variables $\xi_{1},...,\xi_{D}\in\mathbb{R}$ and we can define the
extended warped phase transformation as
$\boldsymbol{w}(t,\xi_{1},...,\xi_{D})=\exp(-\xi_{1})...\exp(-\xi_{D})\boldsymbol{u}(t)$
for $\xi_{1},...,\xi_{D}>0$. We follow a similar derivation as before. From
Section IV, for a general homogeneous linear PDE with first-order time
derivative and maximum $K^{\text{th}}$-order spatial derivative, we can write

$\displaystyle\frac{d\boldsymbol{u}}{dt}=-i\boldsymbol{A}\boldsymbol{u},\qquad\boldsymbol{A}=\sum_{j=1}^{D}\boldsymbol{a}_{j},\qquad\boldsymbol{a}_{j}=\sum_{k=1}^{K}\left(a_{k,j}(\hat{x}_{1},...,\hat{x}_{D})(-i)^{k+1}\hat{p}^{k}_{j}-\frac{i}{K}b(\hat{x}_{1},...,\hat{x}_{D})\right).$
(85)

We can decompose each
$\boldsymbol{a}_{j}=\boldsymbol{a}_{j1}-i\boldsymbol{a}_{j2}$ into its
hermitian and anti-hermitian components, where
$\boldsymbol{a}_{j1}=(1/2)(\boldsymbol{a}_{j}+\boldsymbol{a}_{j}^{\dagger})$,
$\boldsymbol{a}_{j2}=(i/2)(\boldsymbol{a}_{j}-\boldsymbol{a}_{j}^{\dagger})$.
Note here $\boldsymbol{a}_{j1}$ corresponds to odd-derivative terms, while
$\boldsymbol{a}_{j2}$ corresponds to the even-derivative terms. We observe
that for every $\xi_{j}$, $j=1,...,D$, we can write
$\boldsymbol{w}=-\partial\boldsymbol{w}/\partial\xi_{j}$. Thus we can rewrite
Eq. (85) as

$\displaystyle\frac{\partial\boldsymbol{w}(t,\xi_{1},...,\xi_{D})}{\partial
t}=\sum_{j=1}^{D}\boldsymbol{a}_{j2}\frac{\partial\boldsymbol{w}(t,\xi_{1},...,\xi_{D})}{\partial\xi_{j}}-i\boldsymbol{a}_{j1}\boldsymbol{w}(t,\xi_{1},...,\xi_{D}).$
(86)

We can similarly perform an even extension to the domain
$\xi_{1},...,\xi_{D}<0$ by extending the initial condition evenly

$\displaystyle\boldsymbol{w}(0,\xi_{1},...,\xi_{D})=\exp(-|\xi_{1}|)...\exp(-|\xi_{D}|)\boldsymbol{u}(0).$
(87)

If $\tilde{\boldsymbol{w}}$ denotes the Fourier transform of $\boldsymbol{w}$
with respect to $\xi_{1},...,\xi_{D}$ with corresponding Fourier modes
$\eta_{1},...,\eta_{D}$, then for all $\eta_{1},...,\eta_{D}$

$\displaystyle\frac{d\tilde{w}(t,\eta_{1},...,\eta_{D})}{dt}=-i\sum_{j=1}^{D}(\eta_{j}\boldsymbol{a}_{j2}+\boldsymbol{a}_{j1})\tilde{\boldsymbol{w}}(t,\eta_{1},...,\eta_{D}),\qquad\tilde{\boldsymbol{w}}(0,\eta_{1},...,\eta_{D})=\left(\frac{2}{1+\eta^{2}_{1}}\right)...\left(\frac{2}{1+\eta^{2}_{D}}\right)\boldsymbol{u}(0).$
(88)

Now we can also define an extended
$\tilde{\boldsymbol{v}}(t)\equiv\int\tilde{w}(t,x_{1},...,x_{D},\eta_{1},...,\eta_{D})|x_{1}...x_{D}\rangle|\eta_{1}...\eta_{D}\rangle$.
Thus $\tilde{\boldsymbol{v}}(t)$ now obeys a Schrödinger-like equation

$\displaystyle\frac{d\tilde{\boldsymbol{v}}(t)}{dt}=-i\boldsymbol{H}\tilde{\boldsymbol{v}},\qquad\boldsymbol{H}=\sum_{j=1}^{D}\boldsymbol{h}_{j}=\boldsymbol{H}^{\dagger},\qquad\boldsymbol{h}_{j}=\boldsymbol{a}_{j2}\otimes\hat{\eta}_{j}+\boldsymbol{a}_{j1}\otimes
I,\qquad\tilde{\boldsymbol{v}}(0)=|\Xi\rangle^{\otimes D}\boldsymbol{u}(0)$
(89)

which now operates on $2D$ qumodes instead of $D+1$ qumodes. However, the
advantage with this approach is that each $\boldsymbol{h}_{j}$ corresponds to
coupling to a different auxiliary mode. For example, for those PDEs that
require at most pairwise couplings, we only need to simulate $O(D)$
interactions between $D$ completely different pairs of modes instead of $O(D)$
interactions between one particular auxiliary mode and $D$ other modes. The
extension to various PDEs is straightforward. If the Hamiltonian of the
original Schrödingerised system has terms including
$\sum_{j=1}^{D}F_{j}(\hat{x})\hat{p}_{j}^{K}\otimes\hat{\eta}$, then the
corresponding modified terms in the extended Schrödingerised scheme would be
of the form
$\sum_{j=1}^{D}F_{j}(\hat{x})\hat{p}_{j}^{K}\otimes\hat{\eta}_{j}$.

The projection back onto $|u(t)\rangle$ through
$\hat{P}_{>0}=I\otimes\int_{0}^{\infty}|\xi_{1}...\xi_{D}\rangle\langle\xi_{1}...\xi_{D}|d\xi_{1}...d\xi_{D}$
then works in the same way, except now the probability of success in
retrieving $|u(t)\rangle$ becomes
$\|\boldsymbol{u}(t)\|^{2}/(2^{D}\|\boldsymbol{u}(0)\|)$. The initial state
preparation for $|\Xi\rangle^{\otimes D}$ can still be done using Gaussian
states $|G\rangle^{\otimes D}$ for example. However, in this case, the new
fidelity becomes $|\Xi|G\rangle|^{D}$, which decays with $D$. If
$|\langle\Xi|G\rangle|\approx 0.986$ is possible at the optimal $s$, then for
the preparation of $|\Xi\rangle^{\otimes D}$ for $D=5$, we can reach a
fidelity of $\approx 0.932$. Thus, for very high $D$, this extended approach
is much more costly and the original Schrödingerisation scheme in Section IV
is preferred. However, for some physical systems where the extended
Schrödingerisation might be easier to physically implement in the near term,
low $D$ demonstrations are still possible.

##  IX Discussion

In our formulation, we considered primarily the cost in preparation of the
state $|\tilde{v}(t)\rangle$ through Hamiltonian $\boldsymbol{H}$ simulation
when given the initial state $|\Xi\rangle|u(0)\rangle$. We have also outlined
a procedure of retrieving the state $|u(t)\rangle$ whose amplitudes are
proportional to the solutions of the PDE. In cases where $\boldsymbol{H}$ can
be realised directly in analogue quantum simulation, we don’t need to consider
gate count. The procedure is efficient in the number of qumodes $D+1$. When
extended Schrödingerisation is used, there are only $2D$ qumodes. However,
when direct quantum simulation is difficult, then a concatenation of
continuous-variable quantum gates is required. The exact number needs to
computed on an individual basis since it also depends on the form of the
coefficients of the PDE. Given that we know the number of terms of the
Hamiltonian scales like $O(D)$, and if the terms in the Hamiltonian are
pairwise interactions and almost commute, then we expect the gate count to
also be of order $O(D)$, where each gate consists of no more than pairwise
interactions. A detailed analysis will be left to future work and to see in
particular which classes of PDEs are more efficient than others to simulate
and which are not efficient. Note that for classical finite difference methods
(performed on classical digital devices), the scaling is exponential in $D$.
It is also interesting to consider that, in continuous-variable quantum
computation, non-Gaussian operations that are necessary for potential quantum
speedup [25], otherwise efficient classical simulation can suffice. We leave
discussions on this and its connections to efficient simulability of
differential PDEs to future analysis.

We emphasise that, like in the qubit-based approaches to PDEs, our algorithm
is currently also in the form of a quantum subroutine: that we are preparing
quantum states who amplitudes are related to the solution of the PDE. The next
stage is to consider efficient measurements on these final states that can
capture for us observables. Recent methods based on entangled CV measurements
for example [49] be used a step forward, which we leave to the next stage of
investigation.

We observe that when there is time-dependence in the coefficients of the PDE,
our Schrödingerisation procedure works in the same way and in this case we
find $d\tilde{\boldsymbol{v}}(t)/dt=-i\boldsymbol{H}(t)\tilde{\boldsymbol{v}}$
where now $\boldsymbol{H}$ is time-dependent coefficient. When there is an
existing physical system with a corresponding time-dependent Hamiltonian
$\boldsymbol{H}(t)$, then the simulation is straightforward. However, if we
require a decomposition of the corresponding unitary into smaller unitaries
time-ordered operators are needed and it is not so simple. It is less studied
in the continuous-variable scenario and we leave this for future analysis.

An important aspect of solving PDEs is the consideration of boundary
conditions. Using Schrödingerisation to deal with both physical [50] and
artificial [51] boundary conditions have been investigated by the authors in
the case of qubits. How this can be extended to analog systems can be left for
future work when we study more specific systems. It is also important to
remark that if one wishes to simulate analog quantum systems under artificial
boundary conditions, then using another analog quantum system via
Schrödingerisation is necessary. This is because a classical analog (i.e., not
digital) system even in low dimensions could have difficulty in efficiently
simulating the system along with reproducing uniquely quantum features.

An important consideration for near term realisation is to ask if existing
analog quantum systems can be repurposed to simulate linear PDEs, without
decomposing the corresponding unitary into smaller gates. Our formalism can be
very easily used to identify, presented by a given $\boldsymbol{H}$, the
corresponding linear PDE, if it exists. This can be done by following the
arrows backwards in the procedure in Fig. 2. We leave discussions of realistic
errors beyond the ideal limit and noise and any possible truncation in the
physical degrees of freedom to future work.

As a very simple example, suppose we have a quantum optomechanical system with
driving $\boldsymbol{H}_{OM}=\Delta(\hat{x}_{a}^{2}\otimes
I+\hat{p}_{a}^{2}\otimes
I)+\Omega(I\otimes\hat{x}_{b}^{2}+I\otimes\hat{p}_{b}^{2})+g_{1}\hat{x}_{a}\otimes\hat{x}_{b}+g_{2}(\hat{p}_{a}\otimes
I),\Delta,\Omega,g_{1},g_{2}\in\mathbb{R}$, where modes $a,b$ are the photonic
and mechanical modes respectively. Tuning $\Omega\rightarrow 0$, this is able
to simulate a Schrödinger equation with a complex potential and a convection
term $\partial u/\partial t=i\partial^{2}u/\partial t^{2}+g_{2}\partial
u/\partial x-(g_{1}x+i\Delta x^{2})u$. Complex potentials can model absorbing
boundary conditions and is a technique in artificial boundary conditions and
can also be used to model dissipative effects like optical phonon scattering.
Coupled superconducting circuits have a similar Hamiltonian and can include
more nonlinearities. If only independent pairwise couplings are considered in
$\boldsymbol{H}$, then the extended Schrödingerisation method can be used.

There are also many other platforms to explore and we leave this as a
challenge for experimentalists to consider how their own quantum systems –
which only needs a system size of $D+1$ continuous degrees of freedom – can be
repurposed for the simulation of linear PDEs with $D$ spatial dimensions or a
system of $D$ nonlinear ODEs.

Finally, we remark that, historically, analog devices dominated computation
before digital computation, and digital devices didn’t completely dominate
analog computation until the 1980s, and ideas about quantum computation only
emerged after the analog era faded. Today there is a resurgent interest in
classical analog proposals as we move beyond the von Neumann architecture, for
example see [52]. In quantum computation, it is reasonable to consider perhaps
that, just like in classical computation, the first useful devices would also
be analog. Since analog and digital models are fundamentally different, it is
not sensible to compare them directly – either classical or quantum algorithms
– but rather it is important to cultivate both kinds of algorithms to grow our
toolbox for computational methods. We encourage the reader to consider a more
concerted effort toward analog quantum algorithms to enrich methods for
scientific computing.

## Acknowledgements

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). Both authors are also
supported by the Fundamental Research Funds for the Central Universities.

## References

  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm, Physical review letters 110, 250504 (2013).
  * Childs  _et al._ [2021] A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021).
  * Costa  _et al._ [2019] P. C. Costa, S. Jordan, and A. Ostrander, Quantum algorithm for simulating the wave equation, Physical Review A 99, 012323 (2019).
  * Gonzalez-Conde  _et al._ [2022] J. Gonzalez-Conde, A. Rodriguez-Rozas, E. Solano, and M. Sanz, Simulating option price dynamics with exponential quantum speedup, arXiv:2101.04023v3 (2022).
  * Linden  _et al._ [2020] N. Linden, A. Montanaro, and C. Shao, Quantum vs. classical algorithms for solving the heat equation, arXiv preprint arXiv:2004.06516 (2020).
  * Dodin and Startsev [2021] I. Y. Dodin and E. A. Startsev, On applications of quantum computing to plasma simulations, Physics of Plasmas 28, 092101 (2021).
  * Joseph [2020] I. Joseph, Koopman–von neumann approach to quantum simulation of nonlinear classical dynamics, Physical Review Research 2, 043102 (2020).
  * Jin and Liu [2022] S. Jin and N. Liu, Quantum algorithms for computing observables of nonlinear partial differential equations, arXiv preprint arXiv:2202.07834 (2022).
  * Jin  _et al._ [2023a] S. Jin, N. Liu, and Y. Yu, Time complexity analysis of quantum algorithms via linear representations for nonlinear ordinary and partial differential equations, Journal of Computational Physics 487, 112149 (2023a).
  * Liu  _et al._ [2021] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proceedings of the National Academy of Sciences 118 (2021).
  * Krovi [2023] H. Krovi, Improved quantum algorithms for linear and nonlinear differential equations, Quantum 7, 913 (2023).
  * Lloyd  _et al._ [2020] S. Lloyd, G. De Palma, C. Gokler, B. Kiani, Z.-W. Liu, M. Marvian, F. Tennie, and T. Palmer, Quantum algorithm for nonlinear differential equations, arXiv preprint arXiv:2011.06571 (2020).
  * Leyton and Osborne [2008] S. K. Leyton and T. J. Osborne, A quantum algorithm to solve nonlinear differential equations, arXiv preprint arXiv:0812.4423 (2008).
  * Karplus and Soroka [1959] W. J. Karplus and W. W. Soroka,  _Analog methods: Computation and simulation_ (1959).
  * Jackson [1960] A. S. Jackson,  _Analog computation_ (1960).
  * Braunstein and Van Loock [2005] S. L. Braunstein and P. Van Loock, Quantum information with continuous variables, Reviews of modern physics 77, 513 (2005).
  * Lloyd and Braunstein [1999] S. Lloyd and S. L. Braunstein, Quantum computation over continuous variables, Physical Review Letters 82, 1784 (1999).
  * Arrazola  _et al._ [2019] J. M. Arrazola, T. Kalajdzievski, C. Weedbrook, and S. Lloyd, Quantum algorithm for nonhomogeneous linear partial differential equations, Physical Review A 100, 032306 (2019).
  * Knudsen and Mendl [2020] M. Knudsen and C. B. Mendl, Solving differential equations via continuous-variable quantum computers, arXiv preprint arXiv:2012.12220 (2020).
  * [20] A. Barthe, M. Grossi, J. Tura, and V. Dunjko, Continuous variables quantum algorithm for solving ordinary differential equations, In preparation .
  * Jin  _et al._ [2022a] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations via schrodingerisation, arXiv preprint arXiv:2212.13969 (2022a).
  * Jin  _et al._ [2022b] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations via schrodingerisation: technical details, arXiv preprint arXiv:2212.14703 (2022b).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Reviews of Modern Physics 86, 153 (2014).
  * Altman [2021] e. a. Altman, Ehud, Quantum simulators: Architectures and opportunities, PRX Quantum 2, 017003 (2021).
  * Bartlett  _et al._ [2002] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Efficient classical simulation of continuous variable quantum information processes, Physical Review Letters 88, 097904 (2002).
  * Navarrete-Benlloch [2016] C. Navarrete-Benlloch,  _An introduction to the formalism of quantum information with continuous variables_ (Morgan & Claypool Publishers, 2016).
  * Sefi and Van Loock [2011] S. Sefi and P. Van Loock, How to decompose arbitrary continuous-variable quantum operations, Physical review letters 107, 170501 (2011).
  * Hangleiter  _et al._ [2022] D. Hangleiter, J. Carolan, and K. Thebault,  _Analogue Quantum Simulation: A New Instrument for Scientific Understanding_ (Springer, 2022).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Applied physics reviews 6 (2019).
  * Pfister [2019] O. Pfister, Continuous-variable quantum computing in the quantum optical frequency comb, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 012001 (2019).
  * Fukui and Takeda [2022] K. Fukui and S. Takeda, Building a large-scale quantum computer with continuous-variable optical technologies, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 012001 (2022).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Hammerer  _et al._ [2010] K. Hammerer, A. S. Sørensen, and E. S. Polzik, Quantum interface between light and atomic ensembles, Reviews of Modern Physics 82, 1041 (2010).
  * Ortiz-Gutiérrez  _et al._ [2017] L. Ortiz-Gutiérrez, B. Gabrielly, L. F. Muñoz, K. T. Pereira, J. G. Filgueiras, and A. S. Villar, Continuous variables quantum computation over the vibrational modes of a single trapped ion, Optics Communications 397, 166 (2017).
  * Andersen  _et al._ [2015] U. L. Andersen, J. S. Neergaard-Nielsen, P. Van Loock, and A. Furusawa, Hybrid discrete-and continuous-variable quantum information, Nature Physics 11, 713 (2015).
  * Black and Scholes [1973] F. Black and M. Scholes, The pricing of options and corporate liabilities, Journal of political economy 81, 637 (1973).
  * Le Maître and Knio [2010] O. Le Maître and O. M. Knio,  _Spectral methods for uncertainty quantification: with applications to computational fluid dynamics_ (Springer Science & Business Media, 2010).
  * Xiu and Karniadakis [2002] D. Xiu and G. E. Karniadakis, The wiener–askey polynomial chaos for stochastic differential equations, SIAM journal on scientific computing 24, 619 (2002).
  * Gerry and Knight [2005] C. Gerry and P. L. Knight,  _Introductory quantum optics_ (Cambridge university press, 2005).
  * Grover and Rudolph [2002] L. Grover and T. Rudolph, Creating superpositions that correspond to efficiently integrable probability distributions, arXiv preprint quant-ph/0208112 (2002).
  * Jin and Liu [2023] S. Jin and N. Liu, Quantum stochastic galerkin methods, In preparation (2023).
  * Smith [2013] R. C. Smith,  _Uncertainty quantification: theory, implementation, and applications_ , Vol. 12 (Siam, 2013).
  * Engquist and Runborg [1996] B. Engquist and O. Runborg, Multi-phase computations in geometrical optics, Journal of Computational and Applied Mathematics 74, 175 (1996).
  * Sparber  _et al._ [2003] C. Sparber, P. A. Markowich, and N. J. Mauser, Wigner functions versus wkb-methods in multivalued geometrical optics, Asymptotic Analysis 33, 153 (2003).
  * Jin and Li [2003] S. Jin and X. Li, Multi-phase computations of the semiclassical limit of the schrödinger equation and related problems: Whitham vs wigner, Physica D: Nonlinear Phenomena 182, 46 (2003).
  * Jin and Osher [2003] 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, 575 (2003).
  * Lax [1973] P. D. Lax,  _Hyperbolic systems of conservation laws and the mathematical theory of shock waves_ (SIAM, 1973).
  * Crandall and Lions [1983] M. G. Crandall and P.-L. Lions, Viscosity solutions of hamilton-jacobi equations, Transactions of the American mathematical society 277, 1 (1983).
  * Wu  _et al._ [2023] Y.-D. Wu, G. Chiribella, and N. Liu, Quantum-enhanced learning of continuous-variable quantum states, arXiv preprint arXiv:2303.05097 (2023).
  * Jin  _et al._ [2023b] S. Jin, X. Li, N. Liu, and Y. Yu, Quantum simulation for partial differential equations with physical boundary or interface conditions, arXiv preprint arXiv:2305.02710 (2023b).
  * Jin  _et al._ [2023c] S. Jin, N. Liu, X. Li, and Y. Yu, Quantum simulation for quantum dynamics with artificial boundary conditions, arXiv preprint arXiv:2304.00667 (2023c).
  * Bournez and Pouly [2021] O. Bournez and A. Pouly, A survey on analog models of computation, in  _Handbook of Computability and Complexity in Analysis_ (Springer, 2021) pp. 173–226.

##  Appendix A Error estimates for initial state preparation

From section IV, by using $|G\rangle$ instead of $|\Xi\rangle$ as the ancilla
state for example and with the initial state well-approximated with fidelity
$\geq 1-\delta$, then the final state $|\tilde{v}(t)\rangle$ is also well-
approximated with fidelity $\geq 1-\delta$. When there are no errors in the
implementation of the inverse Fourier transform, $|v(t)\rangle$ can similarly
be well-approximated with fidelity $\geq 1-\delta$.

For the error in retrieving $|u(t)\rangle$, we take the heat equation as an
example to demonstrate the error estimation method. Here $\partial u/\partial
t=\nabla^{2}_{x}u$, which means the corresponding $w$ obeys

$\displaystyle\frac{\partial w}{\partial t}=-\nabla^{2}_{x}\frac{\partial
w}{\partial\xi}$ (90)

where $u(t,x)=\int_{0}^{\infty}w(t,x,\xi)d\xi$. In the Schrödingerisation
framework, we choose $w(0,x,\xi)=\exp(-|\xi|)u(0,x)$, so we need the state
$|\Xi\rangle$ for the quantum implementation. However, in the case where we
start with a Gaussian state $|G\rangle$, we need a different initial condition

$\displaystyle
w^{\prime}(0,x,\xi)=\frac{1}{\sqrt{s}\pi^{1/4}}e^{-\xi^{2}/(2s^{2})}u(0,x)$
(91)

where $u(0,x)$ is still identical to the original scenario and only the
ancilla mode changes. Here $w^{\prime}(t,x,\xi)$ still obeys the same PDE as
$w(t,x,\xi)$ in Eq. (90). Then we can define
$u^{\prime}(t,x)=\int_{0}^{\infty}w^{\prime}(t,x,\xi)d\xi$. Now we want to
estimate

$\displaystyle\|\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t)\|^{2}=\int_{-\infty}^{\infty}\left(\int_{0}^{\infty}(w(t,x,\xi)-w^{\prime}(t,x,\xi))d\xi\right)^{2}dx.$
(92)

Both $w(0,x,\xi)$ and $w^{\prime}(0,x,\xi)$ have exponential decay in $\xi>0$
and are waves moving to the left of the $\xi$ domain. This means there exists
$\xi^{*}>0$, $\xi^{*}=O(1)$ such that $w,w\sim 0$ for $\xi>\xi^{*}$. This
means we can replace
$\int_{0}^{\infty}(w-w^{\prime})^{2}d\xi\rightarrow\int_{0}^{\xi^{*}}(w-w^{\prime})^{2}d\xi$.
This implies

$\displaystyle\|\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t)\|^{2}\lesssim\int_{-\infty}^{\infty}\left(\int_{0}^{\xi^{*}}(w(t,x,\xi)-w^{\prime}(t,x,\xi))d\xi\right)^{2}dx$
$\displaystyle\leq\int_{-\infty}^{\infty}\int_{0}^{\xi^{*}}d\xi^{\prime}\int_{0}^{\xi^{*}}(w(t,x,\xi)-w^{\prime}(t,x,\xi))^{2}d\xi
dx\leq\xi^{*}\|\boldsymbol{v}-\boldsymbol{v}^{\prime}\|^{2}$ (93)

where the Cauchy-Schwartz inequality is used in the second line, and
$\boldsymbol{v}^{\prime}(t)=\iint_{-\infty}^{\infty}w^{\prime}(t,x,\xi)|x\rangle|\xi\rangle
dxd\xi$. We showed before $\|\boldsymbol{v}\|=\|\boldsymbol{u}(0)\|$. Using
the normalised Gaussian ancilla state it is simple to see
$\|\boldsymbol{v}^{\prime}\|=\iint_{-\infty}^{\infty}w^{\prime}(0,x,\xi)^{2}d\xi
dx=\|\boldsymbol{u}(0)\|=\|\boldsymbol{v}\|$. Therefore we can write

$\displaystyle 1-\delta\leq|\langle
v|v^{\prime}\rangle|=1-\frac{1}{2\|\boldsymbol{u}(0)\|^{2}}\|\boldsymbol{v}-\boldsymbol{v}^{\prime}\|^{2}.$
(94)

Inserting into Eq. (A) we have the bound

$\displaystyle\|\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t)\|^{2}\lesssim\xi^{*}\delta\|\boldsymbol{u}(0)\|^{2}.$
(95)

The quantum fidelity between $|u(t)\rangle$ and $|u^{\prime}(t)\rangle$ can be
written

$\displaystyle|\langle
u(t)|u^{\prime}(t)\rangle|=1-\frac{1}{2}\left\|\frac{\boldsymbol{u}(t)}{\|\boldsymbol{u}(t)\|}-\frac{\boldsymbol{u}^{\prime}(t)}{\|\boldsymbol{u}^{\prime}(t)\|}\right\|^{2}.$
(96)

Now

$\displaystyle\left\|\frac{\boldsymbol{u}(t)}{\|\boldsymbol{u}(t)\|}-\frac{\boldsymbol{u}^{\prime}(t)}{\|\boldsymbol{u}^{\prime}(t)\|}\right\|=\left\|\frac{\|\boldsymbol{u}^{\prime}(t)\|\boldsymbol{u}(t)-\|\boldsymbol{u}(t)\|\boldsymbol{u}^{\prime}(t)}{\|\boldsymbol{u}(t)\|\|\boldsymbol{u}^{\prime}(t)\|}\right\|=\left\|\frac{\|\boldsymbol{u}^{\prime}(t)\|(\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t))+(\|\boldsymbol{u}^{\prime}(t)\|-\|\boldsymbol{u}(t)\|)\boldsymbol{u}^{\prime}(t)}{\|\boldsymbol{u}(t)\|\|\boldsymbol{u}^{\prime}(t)\|}\right\|$
$\displaystyle\leq\frac{\|\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t)\|}{\|\boldsymbol{u}(t)\|}+\frac{\|\|\boldsymbol{u}^{\prime}(t)\|-\|\boldsymbol{u}(t)\|\|}{\|\boldsymbol{u}(t)\|}\leq
2\frac{\|\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t)\|}{\|\boldsymbol{u}(t)\|}.$
(97)

Therefore

$\displaystyle|\langle u(t)|u^{\prime}(t)\rangle|\geq
1-2\left(\frac{\|\boldsymbol{u}(t)-\boldsymbol{u}^{\prime}(t)\|}{\|\boldsymbol{u}(t)\|}\right)^{2}\gtrsim
1-2\xi^{*}\delta\frac{\|\boldsymbol{u}(0)\|^{2}}{\|\boldsymbol{u}(t)\|^{2}}.$
(98)