# Practical Quantum Computing: solving the wave equation using a quantum
approach

Adrien Suau<EMAIL_ADDRESS>CERFACS, 42 Avenue G. Coriolis, 31057
Toulouse, France  LIRMM, Univ Montpellier, CNRS, Montpellier, France  Gabriel
Staffelbach  CERFACS, 42 Avenue G. Coriolis, 31057 Toulouse, France  Henri
Calandra  TOTAL SA, 2 Avenue Vignancour, 64000 Pau, France

(August 27, 2024)

###### Abstract

In the last years, several quantum algorithms that try to address the problem
of partial differential equation solving have been devised. On one side, the
“direct” quantum algorithms that aim at encoding the solution of the PDE by
executing one large quantum circuit. On the other side, variational algorithms
that approximate the solution of the PDE by executing several small quantum
circuits and making profit of classical optimisers. In this work we propose an
experimental study of the costs (in terms of gate number and execution time on
a idealised hardware created from realistic gate data) associated with one of
the “direct” quantum algorithm: the wave equation solver devised in [PCS.
Costa, S. Jordan, A. Ostrander, Phys. Rev. A 99, 012323, 2019]. We show that
our implementation of the quantum wave equation solver agrees with the
theoretical big-O complexity of the algorithm. We also explain in great
details the implementation steps and discuss some possibilities of
improvements. Finally, our implementation proves experimentally that some PDE
can be solved on a quantum computer, even if the direct quantum algorithm
chosen will require error-corrected quantum chips, which are not believed to
be available in the short-term.

wave equation solver quantum computing algorithm

##  I Introduction

Quantum computing has drawn a lot of attention in the last few years,
following the successive announcements from several world-wide companies about
the implementation of quantum hardware with an increasing number of qubits or
reduced error rates ibm (2019a); int (2019); goo (2019); Arute  _et al._
(2019); Pino  _et al._ (2020).

Along with the hardware improvement, new quantum algorithms were discovered,
yielding potential quantum speed-up and applications in various fields such as
quantum chemistry Cao  _et al._ (2018), linear algebra Harrow  _et al._
(2009); Gilyén  _et al._ (2018); Shao and Xiang (2020); Xu  _et al._ (2019);
Bravo-Prieto  _et al._ (2019); Huang  _et al._ (2019) or optimisation
Kerenidis and Prakash (2018, 2017); Farhi  _et al._ (2014). Recent works even
show that differential equations may be solved by using a quantum computer
Childs and Liu (2019); García-Ripoll (2019); Arrazola  _et al._ (2018);
Srivastava and Sundararaghavan (2018); Xin  _et al._ (2018); Berry  _et al._
(2017); Berry (2010); Leyton and Osborne (2008); Griffin  _et al._ ; Childs
_et al._ (2020); Todorova and Steijl (2020); Lubasch  _et al._ (2019). But
despite the large number of algorithms available, it is hard to find an actual
implementation of a quantum differential equation solver, Hamiltonian
simulation being the unique exception by solving the time-dependant
Schrödinger equation.

In this work, we present and analyse a quantum wave equation solver we
implemented from scratch according to the algorithm depicted in Costa  _et
al._ (2019). During the solver implementation, we had to look for a
Hamiltonian Simulation procedure. The implementations we found being too
restricted, we decided to implement our own Hamiltonian Simulation procedure,
which will also be analysed.

To the best of our knowledge, this work is the first to analyse experimentally
the characteristics of a quantum PDE solver. We checked that the practical
implementation agrees with the theoretical asymptotic complexities on several
quantities of interest such as the total gate count with respect to the number
of discretisation points used or the precision, the number of qubits required
versus the number of discretisation points used to approximate the solution or
precision of the solution when compared to a classical finite-difference
solver. Finally, we verified that the execution time of the generated quantum
circuit on today’s accessible quantum hardware was still following the
theoretical asymptotic complexities devised for the total gate count. Quantum
hardware data were extracted from IBM Q chips.

We show experimentally that it is possible to solve the 1-dimensional wave
equation on a quantum computer with a time-complexity that grows as
$\mathcal{O}\left(N_{d}^{3/2}\log\left(N_{d}\right)^{2}\right)$ where $N_{d}$
is the number of discretisation points used to approximate the solution. But
even if the asymptotic scaling is better than classical algorithms, we found
out that the constants hidden in the big-O notation were huge enough to make
the solver less performant than classical solvers for reasonable
discretisation sizes.

##  II Problem considered

We consider a simplified version of the wave equation on the 1-dimensional
line $[0,1]$ where the propagation speed $c$ is constant and equal to $1$.
This equation can be written as

$\frac{\partial^{2}}{\partial t^{2}}\phi(x,t)=\frac{\partial^{2}}{\partial
x^{2}}\phi(x,t).$ (1)

Moreover, we only consider solving Equation (1) with the Dirichlet boundary
conditions

$\frac{\partial}{\partial x}\phi(0,t)=\frac{\partial}{\partial x}\phi(1,t)=0.$
(2)

No assumption is made on initial speed $\phi(x,0)$ and initial velocity
$\frac{\partial\phi}{\partial t}(x,0)$.

The resolution of this simplified wave equation on a quantum computer is an
appealing problem for the first implementation of a PDE solver for several
reasons. First, the wave equation is a well-known and intensively studied
problem for which a lot of theoretical results have been verified. Secondly,
even-though it is a relatively simple PDE, the wave equation can be used to
solve some interesting problems such as seismic imaging Bamberger  _et al._
(1979, 1977). Finally, the theoretical implementation of a quantum wave
equation solver has already been studied in Childs  _et al._ (2018).

In this paper, we present the complete implementation of a 1-dimensional wave
equation solver using quantum technologies based on qat library. To the best
of our knowledge, this work is the first to consider the implementation of an
entire PDE solver that can run on a quantum computer. Specifically, we explain
all the implementation details of the solver from the mathematical theory to
the actual quantum circuit used. The characteristics of the solver are then
discussed and analysed, such as the estimated gate count and estimated
execution time on real quantum hardware. We show that the implementation
follows the theoretical asymptotic behaviours devised in Costa  _et al._
(2019). Moreover, the wave equation solver algorithm relies critically on an
efficient implementation of a Hamiltonian simulation algorithm, which we have
also implemented and analysed thoroughly.

##  III Implementation

The algorithm used to solve the wave equation is explained in Costa  _et al._
(2019) and uses a Hamiltonian simulation procedure. Costa  _et al._ chose the
Hamiltonian simulation algorithm described in Berry  _et al._ (2015a) for its
nearly optimal theoretical asymptotic behaviour. We privileged instead the
Hamiltonian simulation procedure explained in Ahokas (2004); Berry  _et al._
(2007) for its good experimental results based on Childs  _et al._ (2018) and
its simpler implementation (detailed in Appendix A).

The code has been written using qat, a Python library shipped with the Quantum
Learning Machine (QLM), a package developed and maintained by Atos. It has not
been extensively optimized yet, which means that there is still a large room
for improvements. Possible optimizations and improvements are discussed in
Appendix E.

All the circuits used in this paper have been generated with a subset of qat’s
gate set:

$\left\\{H,X,R_{y}\left(\theta\right),P_{h}\left(\theta\right),CP_{h}\left(\theta\right),CNOT,CCNOT\right\\}$
(3)

and have then been translated to the gate set

$\left\\{U_{1}\left(\lambda\right),U_{2}\left(\lambda,\phi\right),U_{3}\left(\lambda,\phi,\theta\right),CNOT\right\\}$
(4)

for $U_{1}$, $U_{2}$ and $U_{3}$ defined in Equation (7) of Coles  _et al._
(2018) as follow:

$U(\lambda,\phi,\theta)=\begin{pmatrix}\cos\left(\frac{\theta}{2}\right)&-e^{i\lambda}\sin\left(\frac{\theta}{2}\right)\\\
e^{i\phi}\sin\left(\frac{\theta}{2}\right)&e^{i(\lambda+\phi)}\cos\left(\frac{\theta}{2}\right)\end{pmatrix}$
(5) $U_{3}(\lambda,\phi,\theta)=U(\lambda,\phi,\theta)$ (6)
$U_{2}(\lambda,\phi)=U\left(\frac{\pi}{2},\lambda,\phi\right)$ (7)
$U_{1}(\lambda)=U(0,0,\lambda)$ (8)

######  Note 1.

The target gate set presented in Equation (4) does not correspond to the
physical gate set implemented by IBM hardware (see Equation (8) of Coles  _et
al._ (2018)). This choice is justified by the fact that IBM only provides
hardware characteristics such as gate times for the gate set of Equation (4)
and not for the real hardware gate set.

######  Note 2.

Our implementation ignores several practical details such as hardware topology
or the number of available qubits.

This implementation aims at validating in practice the theoretical asymptotic
complexities of Hamiltonian simulation algorithms and providing a proof-of-
concept showing that it is possible to solve a partial differential equation
on a quantum computer.

###  III.1 Sparse Hamiltonian simulation algorithm

######  Definition 1.

$s$-sparse matrix: A $s$-sparse matrix with $s\in\mathbb{N}^{*}$ is a matrix
that has at most $s$ non-zero entries per row and per column

######  Definition 2.

sparse matrix: A sparse matrix is a $s$-sparse matrix with
$s\in\mathcal{O}\left(\log(N)\right)$, $N$ being the size of the matrix.

In the past years, a lot of algorithms have been devised to simulate the
effect of a Hamiltonian on a quantum state Low and Chuang (2017a, 2016, b);
Low (2018); Berry and Childs (2012); Berry  _et al._ (2015a, 2007, b);
Kieferova  _et al._ (2018); Shmoys  _et al._ (2014); Childs and Wiebe (2012);
Childs and Kothari (2011). Among all these algorithms, only few have already
been implemented for specific cases qis (2019); sim (2019) but to the best of
our knowledge no implementation is currently capable of simulating a generic
sparse Hamiltonian.

The domain of application of the already existing methods being too narrow, we
decided to implement our own generic sparse Hamiltonian simulation procedure.
We based our work on the product-formula approach described in Ahokas (2004);
Berry  _et al._ (2007). One advantage of this approach is that product-formula
based algorithms have already been thoroughly analysed both theoretically
Berry  _et al._ (2007); Ahokas (2004) and practically Childs  _et al._ (2018);
Scherer  _et al._ (2017), and several implementations are publicly available,
though restricted to Hamiltonians that can be decomposed as a sum of tensor
products of Pauli matrices. Moreover, Ahokas (2004) provides a lot of
implementation details that allowed us to go straight to the development step.

Our implementation is capable of simulating an arbitrary sparse Hamiltonian
provided that it has already been decomposed into a sum of $1$-sparse
Hermitian matrices with either only real or only complex entries, each
described by an oracle. The implementation has been validated with several
automated tests and a more complex case involving the simulation of a
$2$-sparse Hamiltonian and described in Section III.2. Furthermore, it agrees
perfectly with the theoretical complexities devised in Ahokas (2004); Berry
_et al._ (2007) as studied and verified in Section IV.

###  III.2 Quantum wave equation solver

Using the Hamiltonian simulation algorithm implementation, we successfully
implemented a 1-dimensional wave equation solver using the algorithm described
in Costa  _et al._ (2019) and explained in appendices B and C.

For the specific case considered (Eq. (1) and (2)), solving the wave equation
for a time $T$ on a quantum computer boils down to simulating a $2$-sparse
Hamiltonian for a time $f(T)$, the function $f$ being thoroughly described in
Costa  _et al._ (2019) and Equation (18). The constructed quantum circuit can
then be applied to a quantum state representing the initial position
$\psi(x,0)$ and velocity $\frac{\partial\phi}{\partial t}(x,0)$, and will
evolve this state towards a quantum state representing the final position
$\phi(x,T)$ and velocity $\frac{\partial\phi}{\partial t}(x,T)$.

As for the Hamiltonian simulation procedure, the practical results we obtain
from the implementation of the quantum wave equation solver seems to match the
theoretical asymptotic complexities. See Section IV for an analysis of the
theoretical asymptotic complexities.

##  IV Results

Using a simulator instead of a real quantum computer has several advantages.
In terms of development process, a simulator allows the developer to perform
several actions that are not possible as-is on a quantum processor such as
describing a quantum gate with a unitary matrix instead of a sequence of
hardware operations. Another useful operation that is possible on a quantum
simulator and not currently achievable on a quantum processor is efficient
generic state preparation.

Our implementation uses only standard quantum gates and does not leverage any
of the simulator-only features such as quantum gates implemented from a
unitary matrix. In other words, both the Hamiltonian simulation procedure and
the quantum wave equation solver are “fully quantum” and are readily
executable on a quantum processor, provided that it has enough qubits. As a
proof, and in order to benchmark our implementation, we translated the
generated quantum circuits to IBM Q Melbourne gate-set (see Equation (4)). IBM
Q Melbourne ibm (2019b) is a quantum chip with 14 usable qubits made available
by IBM the 23th of September, 2018.

######  Note 3.

We chose IBM Q Melbourne mainly because, at the time of writing, it was the
publicly accessible quantum chip with the larger number of qubits and so was
deemed to be the closest to future quantum hardware. It is important to note
that even if IBM Q Melbourne has 14 qubits, the quantum circuits constructed
in this paper are not runnable because they require more qubits. Consequently,
because of this hardware limitation, hardware topology has also been left
apart of the study.

This allowed us to have an estimation of the number of hardware gates needed
to either solve the wave equation or simulate a specific Hamiltonian on this
specific hardware. Combining these numbers and the hardware gate execution
time published in mel (2019a), we were able to compute a rough approximation
of the time needed to solve the considered problem presented in Equations (1)
and (2) on this specific hardware.

###  IV.1 Hamiltonian simulation

As explained in Section III.1, the Hamiltonian simulation algorithm
implemented has been first devised in Ahokas (2004); Berry  _et al._ (2007). A
quick review of the algorithm along with implementation details can be found
in Appendix A. This Hamiltonian simulation procedure requires that the
Hamiltonian matrix $H$ to simulate can be decomposed as

$H=\sum_{j=1}^{m}H_{j}$ (9)

where each $H_{j}$ is an efficiently simulable Hermitian matrix.

In our benchmark, we simulated the Hamiltonian described in Equation (39).
According to Ahokas (2004), real $1$-sparse Hermitian matrices with only $1$
or $0$ entries can be simulated with $\mathcal{O}\left(n\right)$ gates and $2$
calls to the oracle, $n$ being the number of qubits the Hamiltonian $H$ acts
on.

Let $O_{i}$ be the gate complexity of the oracle implementing the $i$th
Hermitian matrix $H_{i}$ of the decomposition (9), we end up with an
asymptotic complexity of $\mathcal{O}\left(n+O_{i}\right)$ to simulate
$H_{i}$.

Applying the Trotter-Suzuki product-formula of order $k$ (see Definition 4 in
Appendix A.5 for the definition of the Trotter-Suzuki product-formula) on the
quantum circuit simulating the Hermitian matrices produces a circuit of size

$\mathcal{O}\left(5^{k}\sum_{i=1}^{m}(n+O_{i})\right).$ (10)

This circuit should finally be repeated $r$ times in order to achieve an error
of at most $\epsilon$, with

$r\in\mathcal{O}\left(5^{k}m\tau\left(\frac{m\tau}{\epsilon}\right)^{\frac{1}{2k}}\right),$
(11)

and $\tau=t\max_{i}||H_{i}||$, $t$ being the time for which we want to
simulate the given Hamiltonian and $||\cdot||$ being the spectral norm Berry
_et al._ (2007).

Merging Equations (10) and (11) gives us the complexity

$\mathcal{O}\left(5^{2k}m\tau\left(\frac{m\tau}{\epsilon}\right)^{\frac{1}{2k}}\sum_{i=1}^{m}(n+O_{i})\right).$
(12)

This generic expression of the asymptotic complexity can be specialized to our
benchmark case. The number of gates needed to implement the oracles is
$\mathcal{O}\left(n^{2}\right)$ and the chosen decomposition contains $m=2$
Hermitian matrices, each with a spectral norm of $1$. Replacing the symbols in
Equations (10) and (11) results in the asymptotic gate complexity of

$\mathcal{O}\left(5^{2k}n^{2}\right)$ (13)

for the circuit simulating $e^{-iHt/r}$ and a number

$r\in\mathcal{O}\left(5^{k}t\left(\frac{t}{\epsilon}\right)^{\frac{1}{2k}}\right)$
(14)

of repetitions, which lead to a total gate complexity of

$\mathcal{O}\left(5^{2k}n^{2}t\left(\frac{t}{\epsilon}\right)^{\frac{1}{2k}}\right).$
(15)

In order to check that our implementation follows this theoretical asymptotic
behaviour, we chose to let $k=1$ and plotted the number of gates generated
versus the three parameters that have an impact on the number of gates: the
number of discretisation points $N_{d}$ (Figure 1), the time of simulation $t$
(Figure 2) and the precision $\epsilon$ (Figure 3). The corresponding
asymptotic complexity should be

$\mathcal{O}\left(n^{2}\frac{t^{3/2}}{\sqrt{\epsilon}}\right)=\mathcal{O}\left(\log_{2}\left(N_{d}\right)^{2}\frac{t^{3/2}}{\sqrt{\epsilon}}\right).$
(16)

Figure 1: Plot of the number of quantum gates needed to simulate the
Hamiltonian described in Appendix B using the oracles implemented following
Appendix C. Graph generated with a Trotter-Suzuki product-formula order $k=1$,
a physical time $t=1$ and a precision $\epsilon=10^{-5}$. The value of the
constant $\gamma=210000$ has been chosen arbitrarily to fit the experimental
data. A small discrepancy can be observed: the data seems to be better fitted
by $\log_{2}\left(N\right)^{3/2}$ instead of the theoretical
$\log_{2}\left(N\right)^{2}$. This non-concordance may be due to the fact that
the quantity observed (number of gates) has not reached its asymptotic
behaviour yet.  Figure 2: Plot of the number of quantum gates needed to
simulate the Hamiltonian described in Appendix B using the oracles implemented
following Appendix C. Graph generated with a Trotter-Suzuki product-formula
order $k=1$, $32$ discretisation points (i.e. $n=6$ qubits) and a precision
$\epsilon=10^{-5}$. The value of $\beta=39000000$ has been chosen arbitrarily
to fit the experimental data. Figure 3: Plot of the number of quantum gates
needed to simulate the Hamiltonian described in Appendix B using the oracles
implemented following Appendix C. Graph generated with a Trotter-Suzuki
product-formula order $k=1$, $32$ discretisation points (i.e. $n=6$ qubits)
and a time $t=1$. The value of $\alpha=130000$ has been chosen arbitrarily to
fit the experimental data.

###  IV.2 Wave equation solver

Figure 4: Plot of the number of qubits needed to run the wave equation solver
for a time $t=1$, a precision $\epsilon=10^{-5}$ and a Trotter-Suzuki product-
formula of order $k=1$. The constants values $11$ and $2$ have been chosen
arbitrarily to fit the experimental data.

The first characteristic of the wave equation solver that needs to be checked
is its validity: is the quantum wave equation solver capable of solving
accurately the wave equation as described in Equations (1) and (2)?

To check the validity of the solver, we used qat simulators and Atos QLM to
simulate the quantum program generated to solve the wave equation with
different values for the number of discretisation points $N_{d}$, for the
physical time $t$ and for the precision $\epsilon$. Figure 5 shows the
classical solution versus the quantum solution and the absolute error between
the two solutions for $N_{d}=32$, $t=0.4$ and $\epsilon=10^{-3}$. The solution
obtained by the quantum solver is nearly exactly the same as the classical
solution obtained with finite differences. The error between the two solutions
is of the order of $10^{-7}$, which is $4$ orders of magnitudes smaller than
the error we asked for.

Figure 5: Comparison of the classical solver and the quantum solver. Both
solvers solved the $1$-D wave equation with $N_{d}=32$ discretisation points
and a physical time of $t=0.4$. The classical solver uses finite-differences
with a very small time-step in order to avoid as much as possible errors due
to time-discretisation. The quantum solver was instructed to solve the wave
equation with a precision of at least $\epsilon=10^{-3}$, used a Trotter-
Suzuki order of $k=1$ with the minimised bound for $r$, $r^{min}_{2}$,
described in Childs  _et al._ (2018) and re-introduced in Appendix D. The
solutions of the two solvers are too close to be able to notice a difference
(they overlap on the graph), that is why a second graph plotting the absolute
error between the two solvers is included.

Once the validity of our solver has been checked on multiple test cases, the
next interesting property we would like to verify is the asymptotic cost: does
the implemented simulator seems to agree with the theoretical asymptotic
complexities derived from Costa  _et al._ (2019) and Berry  _et al._ (2007)?

In our specific case, the Hamiltonian $H$ to simulate can be decomposed in two
$1$-sparse Hermitian matrices, both of them having a spectral norm of $1$. The
exact decomposition can be found in Section B.3. We chose to let the product-
formula order be equal to $k=1$ and reuse the asymptotic complexity found in
Equation (15) by changing the time of simulation $t$ by the time $f(t)$:

$\mathcal{O}\left(5^{2k}n^{2}f(t)\left(\frac{f(t)}{\epsilon}\right)^{\frac{1}{2k}}\right).$
(17)

Following the study performed in Costa  _et al._ (2019),

$f(t)=\frac{t}{\delta x}=t\left(N_{d}-1\right)$ (18)

where $\delta x$ is the distance between two discretisation points. Moreover,
it is possible to prove (see Appendix B.3) that

$n=\left\lceil\log_{2}(2N_{d}-1)\right\rceil$ (19)

Replacing $f(t)$ and $n$ in Equations (10) and (11) gives us a gate complexity
of

$\mathcal{O}\left(5^{k}n^{2}\right)$ (20)

to construct a circuit simulating $e^{-iHt/r}$ and a number of repetitions

$r\in\mathcal{O}\left(5^{k}tN_{d}\left(\frac{tN_{d}}{\epsilon}\right)^{\frac{1}{2k}}\right).$
(21)

Merging the two expression results in a gate complexity of

$\mathcal{O}\left(5^{2k}tN_{d}\log_{2}(N_{d})^{2}\left(\frac{tN_{d}}{\epsilon}\right)^{\frac{1}{2k}}\right).$
(22)

Choosing the Totter-Suzuki formula order $k=1$ gives us a final complexity of

$\mathcal{O}\left(N_{d}^{3/2}\log_{2}(N_{d})^{2}\frac{t^{3/2}}{\sqrt{\epsilon}}\right)$
(23)

to solve the wave equation presented in Equation (1). This theoretical result
is verified experimentally in Figure 7.

Figure 6: Estimated execution time of the wave equation solver on IBM Q
Melbourne hardware. The individual gate times have been extracted from mel
(2019b) and mel (2019a). The topology of the hardware is not accounted for,
which means that we cannot use the GF pule timings provided in mel (2019a).
Instead, we chose to use the arithmetic mean. This graph shows the estimated
time with a GF pulse that takes $347$ns to complete, GD pulse of $100$ns and a
buffer time of $20$ns.

######  Note 4.

The huge execution time and number of gates presented in Figures 6 and 7
should be interpreted with care as we may be able to reduce them in the
future. See Appendix E for an explanation on the possible optimisation.

Figure 7: Plot of the number of quantum gates needed to solve the wave
equation described in Equation (1). Graph generated with a Trotter-Suzuki
product-formula order $k=1$, a physical time $t=1$ and a precision
$\epsilon=10^{-5}$. The value of $\lambda=770000$ has been chosen arbitrarily
to fit the experimental data.

##  V Discussion

In this work, we focus on the practical cost of implementing a $1$-dimensional
quantum wave equation solver on a quantum computer. We show that a quantum
computer is able to solve partial differential equations by constructing and
simulating the quantum circuits described. We also study the scaling of the
solver with respect to several parameters of interest and show that the
theoretical asymptotic bounds are mostly verified.

In future works, one can study the possibilities of circuit optimisation
described in Appendix E. It would also be interesting to implement Neumann
boundary conditions instead of Dirichlet ones. A practical implementation
including a non-constant propagation speed $c$ has also been realised during
the writing of this paper. The results were encouraging but we decided to not
include them in the paper to gain more confidence. Finally, future works might
want to extend the wave equation solver to $2$ dimensions or more.

###### Acknowledgements.

The authors would like to thank Reims University, the ROMEO HPC center, Total,
the CCRT and Atos for their support by giving us access to Atos quantum
simulator.

## References

  * ibm (2019a) “Ibm quantum computing,” https://www.ibm.com/quantum-computing/ (2019a), accessed: 2020-03-27.
  * int (2019) “Quantum computing — intel newsroom,” https://newsroom.intel.com/press-kits/quantum-computing/ (2019), accessed: 2020-03-27.
  * goo (2019) “Quantum supremacy using a programmable superconducting processor,” https://ai.googleblog.com/2019/10/quantum-supremacy-using-programmable.html (2019), accessed: 2020-03-27.
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * Pino  _et al._ (2020) J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, C. H. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, C. Ryan-Anderson, and B. Neyenhuis,  (2020), 2003.01293v1 .
  * Cao  _et al._ (2018) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik,  (2018), 1812.09976v2 .
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Physical Review Letters 103 (2009), 10.1103/PhysRevLett.103.150502, phys. Rev. Lett. vol. 15, no. 103, pp. 150502 (2009), 0811.3171v3 .
  * Gilyén  _et al._ (2018) A. Gilyén, Y. Su, G. H. Low, and N. Wiebe,  (2018), 1806.01838v1 .
  * Shao and Xiang (2020) C. Shao and H. Xiang, Phys. Rev. A 101, 022322 (2020).
  * Xu  _et al._ (2019) X. Xu, J. Sun, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan,  (2019), 1909.03898v1 .
  * Bravo-Prieto  _et al._ (2019) C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles,  (2019), 1909.05820v1 .
  * Huang  _et al._ (2019) H.-Y. Huang, K. Bharti, and P. Rebentrost,  (2019), 1909.07344v2 .
  * Kerenidis and Prakash (2018) I. Kerenidis and A. Prakash,  (2018), 1808.09266v1 .
  * Kerenidis and Prakash (2017) I. Kerenidis and A. Prakash,  (2017), 1704.04992v3 .
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann,  (2014), 1411.4028v1 .
  * Childs and Liu (2019) A. M. Childs and J.-P. Liu,  (2019), 1901.00961v1 .
  * García-Ripoll (2019) J. J. García-Ripoll,  (2019), 1909.06619v1 .
  * Arrazola  _et al._ (2018) J. M. Arrazola, T. Kalajdzievski, C. Weedbrook, and S. Lloyd,  (2018), 1809.02622v1 .
  * Srivastava and Sundararaghavan (2018) S. Srivastava and V. Sundararaghavan,  (2018), 1812.10572v2 .
  * Xin  _et al._ (2018) T. Xin, S. Wei, J. Cui, J. Xiao, I. Arrazola, L. Lamata, X. Kong, D. Lu, E. Solano, and G. Long,  (2018), 1807.04553v1 .
  * Berry  _et al._ (2017) D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, Communications in Mathematical Physics 356, 1057 (2017), communications in Mathematical Physics 356, 1057-1081 (2017), 1701.03684v2 .
  * Berry (2010) D. W. Berry,  (2010), 10.1088/1751-8113/47/10/105301, j. Phys. A: Math. Theor. 47, 105301 (2014), 1010.2745v2 .
  * Leyton and Osborne (2008) S. K. Leyton and T. J. Osborne,  (2008), 0812.4423v1 .
  * (24) K. Griffin, S. Jain, T. Flint, and W. Chan, .
  * Childs  _et al._ (2020) A. M. Childs, J.-P. Liu, and A. Ostrander,  (2020), 2002.07868v1 .
  * Todorova and Steijl (2020) B. N. Todorova and R. Steijl, Journal of Computational Physics 409, 109347 (2020).
  * Lubasch  _et al._ (2019) M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch,  (2019), 1907.09032v2 .
  * Costa  _et al._ (2019) P. C. S. Costa, S. Jordan, and A. Ostrander, Physical Review A 99 (2019), 10.1103/PhysRevA.99.012323, phys. Rev. A 99, 012323 (2019), 1711.05394v1 .
  * Bamberger  _et al._ (1979) A. Bamberger, G. Chavent, and P. Lailly, Applied Mathematics & Optimization 5, 1 (1979).
  * Bamberger  _et al._ (1977) A. Bamberger, G. Chavent, and P. Lailly, Ann. Geophys 33, 2 (1977).
  * Childs  _et al._ (2018) A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456 (2018), proceedings of the National Academy of Sciences 115, 9456-9461 (2018), 1711.10980v1 .
  * Berry  _et al._ (2015a) D. W. Berry, A. M. Childs, and R. Kothari, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (2015) pp. 792–809, proceedings of the 56th IEEE Symposium on Foundations of Computer Science (FOCS 2015), pp. 792-809 (2015), 1501.01715v3 .
  * Ahokas (2004) G. R. Ahokas,  _Improved Algorithms for Approximate Quantum Fourier Transforms and Sparse Hamiltonian Simulations_ , mastersthesis, University of Calgary (2004).
  * Berry  _et al._ (2007) D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Communications in Mathematical Physics 270, 359 (2007), communications in Mathematical Physics 270, 359 (2007), quant-ph/0508139v2 .
  * Coles  _et al._ (2018) P. J. Coles, S. Eidenbenz, S. Pakin, A. Adedoyin, J. Ambrosiano, P. Anisimov, W. Casper, G. Chennupati, C. Coffrin, H. Djidjev, D. Gunter, S. Karra, N. Lemons, S. Lin, A. Lokhov, A. Malyzhenkov, D. Mascarenas, S. Mniszewski, B. Nadiga, D. O’Malley, D. Oyen, L. Prasad, R. Roberts, P. Romero, N. Santhi, N. Sinitsyn, P. Swart, M. Vuffray, J. Wendelberger, B. Yoon, R. Zamora, and W. Zhu,  (2018), 1804.03719v1 .
  * Low and Chuang (2017a) G. H. Low and I. L. Chuang, Physical Review Letters 118 (2017a), 10.1103/PhysRevLett.118.010501, phys. Rev. Lett. 118, 010501 (2017), 1606.02685v2 .
  * Low and Chuang (2016) G. H. Low and I. L. Chuang, “Hamiltonian simulation by qubitization,”  (2016), only available as eprint, no journal publication., 1610.06546v2 .
  * Low and Chuang (2017b) G. H. Low and I. L. Chuang, “Hamiltonian simulation by uniform spectral amplification,”  (2017b), only available as eprint. No journal publication., 1707.05391v1 .
  * Low (2018) G. H. Low,  (2018), 1807.03967v1 .
  * Berry and Childs (2012) D. W. Berry and A. M. Childs, Quantum Info. Comput. 12, 29 (2012), quantum Information and Computation 12, 29 (2012), 0910.4157v4 .
  * Berry  _et al._ (2015b) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical Review Letters 114 (2015b), 10.1103/PhysRevLett.114.090502, phys. Rev. Lett. 114, 090502 (2015), 1412.4687v1 .
  * Kieferova  _et al._ (2018) M. Kieferova, A. Scherer, and D. Berry, “Simulating the dynamics of time-dependent hamiltonians with a truncated dyson series,”  (2018), only eprint on arXiv., 1805.00582v1 .
  * Shmoys  _et al._ (2014) D. Shmoys, D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, in  _Proceedings of the 46th Annual ACM Symposium on Theory of Computing - STOC 14_ (2014) pp. 283–292, proceedings of the 46th ACM Symposium on Theory of Computing (STOC 2014), pp. 283-292 (2014), 1312.1414v2 .
  * Childs and Wiebe (2012) A. M. Childs and N. Wiebe,  (2012), 10.26421/QIC12.11-12, quantum Information and Computation 12, 901-924 (2012), 1202.5822v1 .
  * Childs and Kothari (2011) A. M. Childs and R. Kothari, in  _Theory of Quantum Computation, Communication, and Cryptography_ (Springer Berlin Heidelberg, 2011) pp. 94–103, theory of Quantum Computation, Communication, and Cryptography (TQC 2010), Lecture Notes in Computer Science 6519, pp. 94-103 (2011), 1003.3683v2 .
  * qis (2019) “Hamiltonian simulation implementation in qiskit-aqua,” https://github.com/Qiskit/qiskit-aqua/blob/master/qiskit/aqua/operators/weighted_pauli_operator.py#L837 (2019), accessed: 2020-03-27.
  * sim (2019) “Quantum algorithms for the simulation of hamiltonian dynamics,” https://github.com/njross/simcount (2019), accessed: 2020-03-27.
  * Scherer  _et al._ (2017) A. Scherer, B. Valiron, S.-C. Mau, S. Alexander, E. van den Berg, and T. E. Chapuran, Quantum Information Processing 16 (2017), 10.1007/s11128-016-1495-5, quantum Inf Process (2017) 16: 60, 1505.06552v2 .
  * ibm (2019b) “14-qubit backend: IBM Q team, ”IBM Q 16 melbourne backend specifications v1.3.0” (2019),” (2019b), retrieved from https://quantum-computing.ibm.com.
  * mel (2019a) “Melbourne hardware operation execution time,” https://github.com/Qiskit/ibmq-device-information/blob/master/backends/melbourne/V1/version_log.md#gate-specification (2019a), accessed: 2020-03-27.
  * mel (2019b) “Melbourne gate specification,” https://github.com/Qiskit/ibmq-device-information/tree/master/backends/melbourne/V1#gate-specification (2019b), accessed: 2020-03-27.
  * Vazquez (2018) A. C. Vazquez,  _Quantum Algorithm for Solving Tri-Diagonal Linear Systems of Equations_ , mastersthesis, ETH Zürich (2018).
  * Suzuki (1990) M. Suzuki, Physics Letters A 146, 319 (1990).
  * Suzuki (1986) M. Suzuki, Journal of Statistical Physics 43, 883 (1986).
  * Thapliyal and Ranganathan (2017) H. Thapliyal and N. Ranganathan,  (2017), 10.1145/2491682, j. Emerg. Technol. Comput. Syst. 9 (2013) 17:1-17:31, 1712.02630v1 .
  * Cuccaro  _et al._ (2004) S. A. Cuccaro, T. G. Draper, S. A. Kutin, and D. P. Moulton,  (2004), quant-ph/0410184v1 .
  * Draper (2000) T. G. Draper,  (2000), quant-ph/0008033v1 .
  * Barenco  _et al._ (1996) A. Barenco, A. Ekert, K.-A. Suominen, and P. Törmä,  (1996), 10.1103/PhysRevA.54.139, quant-ph/9601018v1 .
  * Cleve and Watrous (2000) R. Cleve and J. Watrous,  (2000), quant-ph/0006004v1 .
  * Häner  _et al._ (2016) T. Häner, M. Roetteler, and K. M. Svore,  (2016), quantum Information and Computation, Vol. 17, No. 7 & 8 (2017), 1611.07995v2 .
  * mul (2015) “Constructing large controlled nots,” https://algassert.com/circuits/2015/06/05/Constructing-Large-Controlled-Nots.html (2015), accessed: 2020-03-27.
  * ibm (2019c) “Ibmq backends information,” https://github.com/Qiskit/ibmq-device-information (2019c), accessed: 2020-03-27.
  * gpr (2019) “gprof2dot main page,” https://github.com/jrfonseca/gprof2dot (2019), accessed: 2020-03-27.
  * Childs  _et al._ (2019) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu,  (2019), 1912.08854v1 .
  * Nam  _et al._ (2018) Y. Nam, N. J. Ross, Y. Su, A. M. Childs, and D. Maslov, npj Quantum Information 4 (2018), 10.1038/s41534-018-0072-4.

##  Appendix A Product-formula implementation details

###  A.1 Hamiltonian simulation

Hamiltonian simulation is the problem of constructing a quantum circuit that
will evolve a quantum state according to a Hamiltonian matrix, following the
Schrödinger equation. In other words, Hamiltonian simulation algorithms
generate a quantum circuit performing the unitary transformation $U$ such that
$||U-e^{-iHt}||<\epsilon$, $H$ being a given Hamiltonian matrix, $t$ a time of
evolution and $\epsilon$ a precision with respect to $||\cdot{}||$, the
spectral norm.

Several quantum algorithms have been developed in the last few years to solve
the problem of $s$-sparse Hamiltonian simulation Low and Chuang (2017a, 2016,
b); Low (2018); Berry and Childs (2012); Berry  _et al._ (2015a, 2007, b);
Kieferova  _et al._ (2018); Shmoys  _et al._ (2014); Childs and Wiebe (2012);
Childs and Kothari (2011). Among these algorithms we decided to implement the
product-formula approach Berry  _et al._ (2007); Ahokas (2004), for the
reasons presented in Section III.1.

The product formula algorithm has three main steps: decompose, simulate,
recompose. It works by first decomposing the $s$-sparse Hamiltonian matrix $H$
that should be simulated as a sum of Hermitian matrices $H_{j}$ that are
considered easy to simulate

$H=\sum_{j=0}^{m-1}H_{j}.$ (24)

The second step is then to simulate each $H_{j}$ separately, i.e. to create
quantum circuits implementing $e^{-iH_{j}t}$ for all the $H_{j}$ in the
decomposition in Equation (24). The last step uses the simulations computed in
step two to approximate $e^{-iHt}$.

The very first questions that should be answered before starting any
implementation of the product-formula algorithm are “What is an easy to
simulate matrix?” and “What kind of Hermitian matrices are easy to simulate?”.

###  A.2 Easy to simulate matrices

One of the most desirable properties for an “easy to simulate” matrix is the
possibility to simulate it exactly, i.e. to construct a quantum circuit that
will perfectly implement $e^{-iHt}$. This property becomes a requirement when
one wants rigorous bounds on the error of the final simulation. Another
enviable property of these matrices is that they can be simulated with a low
gate number and only a few calls to the matrix oracle.

######  Definition 3 (Easy to simulate matrix).

A Hermitian matrix $H$ can be qualified as “easy to simulate” if there exist
an algorithm that takes as input a time $t$ and the matrix $H$ and outputs a
quantum circuit $C(H)_{t}$ such that

  1. 1.

The quantum circuit $C(H)_{t}$ implements exactly the unitary transformation
$e^{-iHt}$, i.e.

$\left|\left|e^{-iHt}-C(H)_{t}\right|\right|=0.$

  2. 2.

The algorithm only needs $\mathcal{O}\left(1\right)$ calls to the oracle of
$H$ and $\mathcal{O}\left(\log N\right)$ additional gates, $N$ being the
dimension of the matrix $H$.

With this definition of an “easy to simulate” matrix, we can now search for
matrices or group of matrices that satisfy this definition.

####  A.2.1 Multiples of the identity

The first and easiest matrices that fulfil the easy to simulate matrix
requirements are the multiples of the identity matrix $\left\\{\alpha
I,\,\alpha\in\mathbb{R}\right\\}$ with $I$ the identity matrix. The quantum
circuit to simulate this class of matrices can be found in Vazquez (2018) .

####  A.2.2 $1$-sparse Hermitian matrices

A larger class of matrices that can be efficiently and exactly simulated are
the $1$-sparse, integer weighted, Hermitian matrices. Quantum circuits
simulating exactly $1$-sparse matrices with integer weights can be found in
Ahokas (2004).

######  Note 5.

Procedures simulating $1$-sparse matrices with real (non-integers) weights are
also described in the paper, but these matrices do not fall in the “easy to
simulate” category because the procedures explained are exact only if all the
matrix weights can be represented exactly with a fixed-point representation,
which is not always verified.

######  Note 6.

Multiples of identity matrices presented in Section A.2.1 are a special case
of $1$-sparse matrices. The two classes have been separated because more
efficient quantum circuits exists for $\alpha I$ matrices.

###  A.3 Decomposition of $H$

Once the set of “easy to simulate” matrices has been established, the next
step of the algorithm is to decompose the $s$-sparse matrix $H$ as a sum of
matrices in this set.

There are two possible ways of performing this decomposition, each one with
its advantages and drawbacks: applying a procedure computing the decomposition
automatically, or decompose the matrix $H$ beforehand and provide the
decomposition to the algorithm.

The first solution, which is to automatically construct the oracles of the
$H_{j}$ matrices from the oracle of the $H$ matrix has been studied in Ahokas
(2004) and Childs and Kothari (2011). Thanks to this automatic decomposition
procedure, we only need to implement one oracle. This simplicity comes at the
cost of a higher gate count: each call to the automatically constructed
oracles of the matrices $H_{j}$ will require several calls to the oracle of
$H$ along with additional gates.

On the other hand, the second solution offers more control at the cost of less
abstraction and more work. The decomposition of $H$ is not automatically
computed and should be performed beforehand. Once the matrix $H$ has been
decomposed as in Equation (9), the oracles for the matrices $H_{j}$ should be
implemented. This means that we should now implement $m$ oracles instead of
only $1$ for the first solution. The main advantage of this method over the
one using automatic-decomposition is that it gives us more control, a control
that can be used to optimize even more the decomposition of Equation (24)
(less $H_{j}$ in the decomposition, $H_{j}$ matrices that can be simulated
more efficiently, …).

All the advantages and drawbacks weighted, we chose to implement the second
option for several reasons. First, the implementation of the automatic
decomposition procedure adds a non-negligible implementation complexity to the
whole Hamiltonian simulation procedure. Moreover, the automatic decomposition
procedure can be implemented afterwards and plugged effortlessly to the non-
automatic implementation. Finally, our use-case only required to simulate a
$2$-sparse Hamiltonian that can be decomposed as the sum of two $1$-sparse
Hermitian matrices, which makes the manual decomposition step manageable.

###  A.4 Simulation of the $H_{j}$

Once the matrix $H$ has been decomposed following Equation (24) with each
$H_{j}$ being an “easy to simulate” matrix, the simulation of $H_{j}$ becomes
a straightforward application of the procedures described in Section A.2.

After this step, we have access to quantum circuits implementing
$e^{-iH_{j}t}$ for $j\in\left[0,m-1\right]$ and $t\in\mathbb{R}$.

###  A.5 Re-composition of the $e^{-iH_{j}t}$

The ultimate step of the algorithm is to approximate the desired evolution
$e^{-iHt}$ with the evolutions $e^{-iH_{j}t}$. In the special case of mutually
commuting $H_{j}$, this step is trivial as it boils down to use the properties
of the exponential function on matrices and write
$e^{iHt}=e^{i\sum_{j}H_{j}t}=\prod_{j}e^{iH_{j}t}$. But in the more realistic
case where the matrices $H_{j}$ do not commute, a more sophisticated method
should be used to approximate the evolution $e^{-iHt}$. To this end, we used
the first-order Lie-Trotter-Suzuki product formula defined in Definition 4.

######  Definition 4 (Lie-Trotter-Suzuki product formula Suzuki (1990, 1986);
Childs  _et al._ (2018)).

The Lie-Trotter-Suzuki product formula approximates

$\exp\left(\lambda\sum_{j=0}^{m}\alpha_{j}H_{j}\right)$ (25)

with

$S_{2}(\lambda)=\prod_{j=0}^{m-1}e^{\alpha_{j}H_{j}\lambda/2}\prod_{j=m-1}^{0}e^{\alpha_{j}H_{j}\lambda/2}$
(26)

and can be generalized recursively to higher-orders

$\begin{split}S_{2k}\left(\lambda\right)=&\left[S_{2k-2}\left(p_{k}\lambda\right)\right]^{2}\\\
&\times
S_{2k-1}\left(\left(1-4p_{k}\right)\lambda\right)\left[S_{2k-2}\left(p_{k}\lambda\right)\right]^{2}\end{split}$
(27)

with $p_{k}=\left(4-4^{1/(2k-1)}\right)^{-1}$ for $k>1$. Using this formula,
we have the approximation

$e^{\lambda
H}=\left[S_{2k}\left(\frac{\lambda}{n}\right)\right]^{n}+\mathcal{O}\left(\frac{|\lambda|^{2k+1}}{n^{2k}}\right).$
(28)

We used the Lie-Trotter-Suzuki product formula with $\lambda=-it$ to
approximate the operator $e^{-iHt}$ up to an error of
$\epsilon\in\mathcal{O}\left(\frac{t^{2k+1}}{n^{2k}}\right)$.

##  Appendix B Hermitian matrix construction and decomposition

One of the main challenge in implementing a quantum wave equation solver lies
in the construction and implementation of the needed oracles. This appendix
describes the first step of the implementation process: the construction and
decomposition of the Hamiltonian matrix that will be simulated using the
Hamiltonian simulation procedure introduced in Appendix A.

This appendix follows the analysis performed in Costa  _et al._ (2019) and
adds details and observations that will be refereed to in Appendix C when
dealing with the actual oracle implementation.

###  B.1 Hamiltonian matrix description

In order to devise the Hamiltonian matrix that should be simulated to solve
the wave equation, the first step is to discretise Equation (1) with respect
to space. Such a discretisation can be seen as a graph $G_{\delta x}$ whose
vertices are the discretisation points and with edges between nearest
neighbour vertices. The graph $G_{\delta x}$ is depicted in Figure 8.

Figure 8: Graph $G_{\delta x}$ built from the discretisation of the
$1$-dimensional line $[0,1]$ with $N_{d}$ discretisation points (i.e. $\delta
x=\frac{1}{N_{d}-1}$).

The graph Laplacian of $G_{\delta x}$, defined as

$L(G_{\delta x})_{i,j}:=\begin{cases}\deg(v_{i})&\\!\\!\\!\\!\text{if }i=j\\\
-1&\\!\\!\\!\\!\text{if }(i\neq j)\land(v_{i}\mbox{ adjacent to }v_{j})\\\
0&\\!\\!\\!\\!\text{otherwise}\\\ \end{cases}$ (29)

can then be used to approximate the differential operator
$\frac{\partial^{2}}{\partial x^{2}}$. By using the discretisation
approximation

$\frac{\partial^{2}\phi}{\partial x^{2}}(i\delta
x,t)\approx\frac{\phi_{i-1,t}-2\phi_{i,t}+\phi_{i+1,t}}{{\delta x}^{2}}$ (30)

with $\phi_{i,t}=\phi(i\delta x,t)$, and approximating $\phi(x,t)$ with a
vector $\phi=\left[\phi_{i,t}\right]_{0\leqslant i<N_{d}}$, the matrix

$A=-\frac{1}{{\delta x}^{2}}L(G_{\delta x})$ (31)

approximates the second derivative of $\phi$ when $\delta x\to 0$ as

$\left[A\phi\right]_{i}=\frac{\phi((i-1)\delta x,t)-2\phi(i\delta
x,t)+\phi((i+1)\delta x,t)}{{\delta x}^{2}}.$ (32)

The approximation (31) is then used in Equation (1) to approximate the spatial
derivative operator:

$\frac{\partial^{2}}{\partial t^{2}}\phi=-\frac{1}{{\delta x}^{2}}L(G_{\delta
x})\phi.$ (33)

Based on this formula, Costa  _et al._ (2019) shows that simulating

$H=\begin{pmatrix}0&B\\\ B^{\dagger}&0\\\ \end{pmatrix}$ (34)

with

$BB^{\dagger}=L$ (35)

constructs a quantum circuit that will evolve a part of the quantum state it
is applied on according to the discretised wave equation (33).

A matrix $B$ satisfying (35) can be obtained directly from the graph
$G_{\delta x}$ representing the discretisation. The algorithm to construct the
matrix $B$ can be decomposed in three steps. First, the vertices
(discretisation points) should be arbitrarily ordered by assigning them a
unique index in $[0,N_{d}-1]$. Then, each edge of the graph is arbitrarily
oriented and indexed with indices in $[0,N_{d}-2]$. Finally, $B$ is computed
with the following definition

$B_{ij}=\left\\{\begin{array}[]{ll}1&\text{if edge }j\text{ is a self-loop of
vertex }i\text{,}\\\ 1&\text{if edge }j\text{ has vertex }i\text{ as
source,}\\\ -1&\text{if edge }j\text{ has vertex }i\text{ as sink,}\\\
0&\text{otherwise}\\\ \end{array}\right..$ (36)

Note that edges’ orientation and vertices/edges ordering is completely
arbitrary. Changing either the edges orientation on one of the orderings will
change the matrix $B$ but will not affect $BB^{\dagger}$ which should be equal
to $L$. This freedom in the ordering and orientation choices takes a crucial
importance in the oracle implementation as it allows us to pick the
ordering/orientation that will produce an easy-to-implement matrix $B$.

###  B.2 Dirichlet boundary conditions

Fixing boundary conditions is a requirement for most of the partial
differential equations to admit a unique well-defined solution. There exist
several boundary conditions such as Neumann, Dirichlet, Robin or Cauchy ones.
For simplicity, we restricted ourselves to the study of Equation (1) with
Dirichlet boundary condition of Equation (2).

In the case of Dirichlet boundary conditions on the $1$-dimensional line
$[0,1]$, the two boundary nodes at $x=0$ and $x=1$ can be ignored as their
value is always equal to $0$. Moreover, Costa  _et al._ (2019) shows that the
graph $G_{\delta x}^{D}$ representing the discretisation with Dirichlet
boundary conditions (2) is simply $G_{\delta x}$ with self-loops on the two
outer nodes (i.e. the ones indexed $1$ and $N_{d}-2$ as $0$ and $N_{d}-1$ are
ignored). $G_{\delta x}^{D}$ is depicted in Figure 9.

Figure 9: Graph $G_{\delta x}^{D}$ representing the discretisation of the
1-dimensional line $[0,1]$ with Dirichlet boundary conditions. The points and
edges in grey are only drawn for illustration purpose and are ignored in the
analysis because the boundary condition impose a value of $0$ on these
vertices. Loops are added to $G_{\delta x}$ to encode the fact that this graph
represents Dirichlet boundary conditions. Vertices (resp. edges) are ordered
with indices within $[0,N_{c}-1]$ (resp. $[0,N_{c}+1]$). $N_{c}$ is the number
of considered points and is equal to $N_{d}-2$ (the two extremal points are
ignored).

###  B.3 Matrices construction

All the pieces are now in place to start building the matrix
$B_{\text{d}}\in\mathbb{R}^{(N_{c}-1)\times N_{c}}$. Using the definition of
the matrix $B$ written in Equation (36) and the graph $G_{\delta x}^{D}$
depicted in Figure 9 we end up with

$B_{\text{d}}=\begin{pmatrix}1&1&0&\cdots&0\\\ 0&-1&1&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&0\\\ 0&\cdots&0&-1&1\\\ \end{pmatrix}$ (37)

We can easily check that $B_{\text{d}}B_{\text{d}}^{\dagger}$ is equal to the
well-known discretisation matrix

$B_{\text{d}}B_{\text{d}}^{\dagger}=\begin{pmatrix}2&-1&0&\cdots&0\\\
-1&2&\ddots&\ddots&\vdots\\\ 0&\ddots&\ddots&\ddots&0\\\
\vdots&\ddots&\ddots&2&-1\\\ 0&\cdots&0&-1&2\\\ \end{pmatrix},$ (38)

which validate the method of construction of $B_{\text{d}}$.

Computing $\widetilde{H_{\text{d}}}$, the Hamiltonian matrix that should be
simulated to evolve the quantum state according to the wave equation (1) with
Dirichlet boundary conditions, is now straightforward. Using Equation (34), we
directly obtain

$\widetilde{H_{d}}=\frac{1}{\delta
x}\begin{pmatrix}0&\cdots&\cdots&0&1&1&0&\cdots&0\\\
\vdots&&&\vdots&0&-1&1&\ddots&\vdots\\\
\vdots&&&\vdots&\vdots&\ddots&\ddots&\ddots&0\\\
0&\cdots&\cdots&0&0&\cdots&0&-1&1\\\ 1&0&\cdots&0&0&\cdots&\cdots&\cdots&0\\\
1&-1&\ddots&\vdots&\vdots&&&&\vdots\\\ 0&1&\ddots&0&\vdots&&&&\vdots\\\
\vdots&\ddots&\ddots&-1&\vdots&&&&\vdots\\\
0&\cdots&0&1&0&\cdots&\cdots&\cdots&0\\\ \end{pmatrix}$ (39)

As explained in Appendix A, the Hamiltonian simulation algorithm implemented
requires that the Hamiltonian to simulate is split as a sum of $1$-sparse
Hermitian matrices. There are a lot of valid decompositions for the matrix
$\widetilde{H_{d}}$ and we are free to choose the decomposition that will
simplify the most the oracle implementation or reduce the gate complexity.

We made the choice to decompose $B_{\text{d}}$ as two $1$-sparse matrices and
then reflect this decomposition on $\widetilde{H_{d}}$. Let $B_{1}$ and
$B_{-1}$ defined as

$B_{1}=\begin{pmatrix}0&1&0&\cdots&0\\\ \vdots&\ddots&1&\ddots&\vdots\\\
\vdots&&\ddots&\ddots&0\\\ 0&\cdots&\cdots&0&1\\\ \end{pmatrix}$ (40)
$B_{-1}=\begin{pmatrix}1&0&\cdots&\cdots&0\\\ 0&-1&\ddots&&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\ 0&\cdots&0&-1&0\\\ \end{pmatrix}$ (41)

we have $B_{\text{d}}=B_{1}+B_{-1}$. Let also

$\widetilde{H_{1}}=\frac{1}{\delta x}\begin{pmatrix}0&B_{1}\\\
{B_{1}}^{\dagger}&0\\\ \end{pmatrix},\;\widetilde{H_{-1}}=\frac{1}{\delta
x}\begin{pmatrix}0&B_{-1}\\\ {B_{-1}}^{\dagger}&0\\\ \end{pmatrix},$ (42)

it is easy to see that
$\widetilde{H_{d}}=\widetilde{H_{1}}+\widetilde{H_{-1}}$ and that both
$\widetilde{H_{1}}$ and $\widetilde{H_{-1}}$ are $1$-sparse Hermitian
matrices.

For convenience, we also define

$H_{1}=\begin{pmatrix}0&B_{1}\\\ {B_{1}}^{\dagger}&0\\\ \end{pmatrix},\quad
H_{-1}=\begin{pmatrix}0&B_{-1}\\\ {B_{-1}}^{\dagger}&0\\\ \end{pmatrix},$ (43)

and $H_{d}=H_{1}+H_{-1}$, the $\widetilde{H_{1}}$, $\widetilde{H_{-1}}$ and
$\widetilde{H_{d}}$ matrices rescaled to contain only integer weights. These
matrices have the interesting property that simulating $\widetilde{H_{d}}$
(resp. $\widetilde{H_{1}}$, $\widetilde{H_{-1}}$) for a time $t$ is equivalent
to simulating $H_{d}$ (resp. $H_{1}$, $H_{-1}$) for a time $\frac{t}{\delta
x}$. This property will be used in the following sections as it offers us the
opportunity to simulate the integer-weighted matrices $H_{d}$, $H_{1}$ and
$H_{-1}$ instead of the real-weighted ones $\widetilde{H_{d}}$,
$\widetilde{H_{1}}$ and $\widetilde{H_{-1}}$.

Note also that a lower bound of the number of qubits needed to solve the wave
equation for $N_{d}$ discretisation points can be computed from the dimensions
of $H_{d}$. As the non-empty upper-left block of matrix $H_{d}$ is of
dimension $\left(2N_{d}-1\right)\times\left(2N_{d}-1\right)$, we need at least

$\left\lceil\log_{2}(2N_{d}-1)\right\rceil$ (44)

qubits to simulate it. This estimation does not take into account ancilla
qubits that may be needed to implement the oracles.

##  Appendix C Oracle construction

Oracles can be seen as the interface between a quantum procedure and real-
world data. Their purpose is to encode classical data such that a quantum
algorithm can process it efficiently.

###  C.1 Oracle interface

In order to work as a bridge between the classical and the quantum worlds and
to be used by the quantum algorithm, a clear interface for the oracle should
be established.

We chose to use the interface described in (Ahokas, 2004, Eq. 4.4) with slight
modifications improving the arity of the oracle for our specific case of
$1$-sparse matrices.

More precisely, our oracles $O$ implement the following interface

$O\left|x_{0}\right\rangle\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle=\left|x_{0}\right\rangle\left|m(x_{0})\right\rangle\left|v(x_{0})\right\rangle\left|s(x_{0})\right\rangle$
(45)

with $\left|x_{0}\right\rangle$ encoding a row index as a unsigned integer,
$m(x)$ the function that returns the column index of the only non-zero element
in row $x$, $v(x)=|w(x)|$ the absolute value of the weight $w(x)$ of the first
non-zero element in row $x$ and

$s(x)=\left\\{\begin{split}0&\text{ if }w(x)\geqslant 0\\\ 1&\text{ else
}\end{split}\right..$ (46)

the sign of the first non-zero entry in row $x$.

The interface of the oracle $O$ can also be obtained with 3 separate oracles
that will each take care of computing one output:

$M\left|x_{0}\right\rangle\left|0\right\rangle=\left|x_{0}\right\rangle\left|m(x)\right\rangle$
(47)
$V\left|x_{0}\right\rangle\left|0\right\rangle=\left|x_{0}\right\rangle\left|v(x)\right\rangle$
(48)
$S\left|x_{0}\right\rangle\left|0\right\rangle=\left|x_{0}\right\rangle\left|s(x)\right\rangle$
(49)

###  C.2 Optimisation of $M$ and $S$

######  Claim 1.

The simulation algorithms provided by Ahokas (2004) have the interesting
property that if the oracle $V$ encodes a weight of zero for some inputs (i.e.
$v(x)=0$ for some $x$) then the outputs of oracles $M$ and $S$ are ignored for
those inputs.

###### Proof.

The circuit simulating a $1$-sparse $m$-bit-integer weighted Hamiltonian $H$
depicted in Figure 10 is taken from Ahokas (2004).

Figure 10: Quantum circuit re-created from (Ahokas, 2004, p. 71) that
simulates a 1-sparse integer-weighted Hamiltonian for a given time $t$.
$F_{m}$ is defined as the diagonal matrix with diagonal entries increasing
from $0$ to $2^{m}-1$.

In our special case of $1$-bit weights (i.e. $m=1$), the third quantum gate
$e^{-iZ\otimes Ft}$ can be written as

$\begin{split}e^{-iZ\otimes Ft}&=\exp\left[-i\begin{pmatrix}F&0\\\ 0&-F\\\
\end{pmatrix}t\right]\\\ &=\begin{pmatrix}e^{-iFt}&0\\\ 0&e^{iFt}\\\
\end{pmatrix}\\\ &=\begin{pmatrix}1&0&0&0\\\ 0&e^{-it}&0&0\\\ 0&0&1&0\\\
0&0&0&e^{-it}\\\ \end{pmatrix}.\end{split}$ (50)

where $F=F_{1}$ and

$F_{m}=\begin{pmatrix}0&0&\cdots&\cdots&\cdots&\cdots&0\\\
0&1&\ddots&&&&\vdots\\\ \vdots&\ddots&2&\ddots&&&\vdots\\\
\vdots&&\ddots&3&\ddots&&\vdots\\\ \vdots&&&\ddots&4&\ddots&\vdots\\\
\vdots&&&&\ddots&\ddots&0\\\ 0&\cdots&\cdots&\cdots&\cdots&0&2^{m}-1\\\
\end{pmatrix}.$ (51)

It follows from the matrix notation that if the second qubit $e^{-iZ\otimes
Ft}$ is applied on is in the state $\left|0\right\rangle$, the gate
$e^{-iZ\otimes Ft}$ is the identity transformation, i.e. the unitary operation
$e^{-iZ\otimes Ft}$ sends $\left|00\right\rangle$ (resp.
$\left|10\right\rangle$) to $\left|00\right\rangle$ (resp.
$\left|10\right\rangle$). This means that if the oracle $O$ does not set the
last qubit to $\left|1\right\rangle$ (i.e. encodes a weight of $0$ for the
$x$th row of $H$), the quantum circuit depicted in Figure 10 can be simplified
up to an identity transformation as the effects of $O$ (resp. $A$) are
reverted by $O^{\dagger}$ (resp. $A^{\dagger}$).

Rephrasing, if the $x$th row of matrix $H$ has no non-zero entries, the
effects of the oracle $O$ is ignored, which implies that the effects of the
oracles $M$ and $S$ that compose $O$ are also ignored. ∎

Using the result of Claim 1, we are free to implement any transformation that
best suits us for the set of inputs $\left|x\right\rangle$ such that the $x$th
row of the considered Hermitian matrix ($H_{1}$ or $H_{-1}$) has no non-zero
elements as long as the oracle $V$ implements the right transformation.

To illustrate clearly the _implemented_ transformations we chose to encode
with $M$ and $S$, the next sections will re-write the matrices $H_{1}$ and
$H_{-1}$ according to Equation (43) but with one $\mathbf{0}$ or $\mathbf{-0}$
in each empty row. A $\mathbf{0}$ entry at position $(i,j)$ in the matrix
means that the row $i$ was empty, the oracle $M$ will map
$\left|i\right\rangle$ to $\left|j\right\rangle$ and the oracle $S$ will
encode a positive sign, i.e. $\left|0\right\rangle$. The same reasoning
applies for $\mathbf{-0}$ entries, except that the encoded sign is now
negative, i.e. $\left|1\right\rangle$.

The following sections will explain step by step the construction of each of
the three oracles $M$, $V$ and $S$, both for the matrix $H_{1}$ ($M_{1}$,
$V_{1}$ and $S_{1}$) and the matrix $H_{-1}$ ($M_{-1}$, $V_{-1}$ and
$S_{-1}$).

###  C.3 About arithmetic and logic quantum gates

Implementing the oracles $M$, $V$ and $S$ for the matrices $H_{1}$ and
$H_{-1}$ requires several arithmetic and logic quantum gates such as or, add
or compare.

All these gates have been implemented prior to the oracle implementation and
the implementation steps are detailed in this section.

####  C.3.1 The or gate

The or gate is easily implemented using only X and CCX (or Toffoli) gates. The
implementation used is depicted in Figure 11 and uses the famous Boole algebra
formula linking not, or and and: $x\vee y=\neg(\neg x\wedge\neg y)$.

|
---
| =
---
|
---

Figure 11: Implementation of the or gate.

####  C.3.2 The add and sub gates

Most of the research papers presenting an implementation of the add or sub
gates only consider the case where the two numbers to add or subtract are
stored in quantum registers.

In our case, the oracles implementation requires an adder and subtractor that
can add or subtract to a quantum register a quantity known when the quantum
circuit is generated, i.e. not necessarily encoded on a quantum state.

######  Claim 2.

Implementing a subtractor is trivial once an adder procedure is available.

###### Proof.

A subtractor can be implemented from a generic adder by using the identity

$a-b=(a^{\prime}+b)^{\prime}$ (52)

where ′ denotes the bitwise complementation.

The circuit resulting of the application of this identity is depicted in
Figure 12 and only requires one call to the adder and $2n$ additional gates,
$n$ being the number of qubits used to represent one of the operands.

|
---
| =
---
|
---

Figure 12: Implementation of the sub gate from an add gate. The $y$ value
encoding is intentionally omitted. The subtractor will use the same encoding
as the adder (i.e. either the $y$ value is encoded on a quantum register or it
is encoded directly in the quantum circuit implementing the adder). Note that
the $y$ value is not negated.

∎

######  Note 7.

Following Claim 2 we will restrict the study to implementing an adder.
Implementing a subtractor is trivial and cheap in term of additional quantum
gates used once an adder is available.

######  Definition 5.

Generation-time value A generation-time value is a value that is known by the
programmer when generating the quantum circuit. Knowing a value at generation-
time may allow to optimise even further the generated quantum circuit. The
closest analog in classical programming would be C-like macros.

The easiest solution to overcome the problem caused by the non-compatible
input formats is to encode the quantity known at generation-time into
ancillary qubits and then use the regular adder algorithms to add to a quantum
register the value encoded in a second quantum register. Even if this solution
is trivial to implement, it has the huge downside of requiring
$\mathcal{O}\left(\log n\right)$ additional ancillary qubits to temporarily
store the generation-time value.

Another answer to the problem would be to adapt a quantum adder originally
devised to add two quantum registers to a quantum adder capable of adding a
constant value to a quantum register. Several adders Thapliyal and Ranganathan
(2017); Cuccaro  _et al._ (2004); Draper (2000) have been studied to check if
they can be modified to allow a generation-time input, i.e. if it possible to
remove completely the quantum register storing the right-hand-side (or left-
hand-side) of the addition.

The task of removing the quantum register storing one of the operands appears
to be challenging for adders based on classical arithmetic like Cuccaro  _et
al._ (2004); Thapliyal and Ranganathan (2017) but trivial for Draper’s quantum
adder introduced in Draper (2000).

######  Claim 3.

Draper’s quantum adder can be adapted into an efficient adder that takes as
right-hand side input a unsigned “generation-time” integer value and add this
value to a sufficiently large quantum register encoding another unsigned
integer.

###### Proof.

The original Drapper’s adder as introduced in Draper (2000) is illustrated in
Figure 13.

Figure 13: Original Draper’s adder example for $3$-qubit registers
$\left|a\right\rangle$ and $\left|b\right\rangle$. The round gates between the
two applications of the Quantum Fourier Transform (QFT gates) are controlled
phase gates and are defined in Draper (2000). Note that the adder wraps on
overflow, meaning that if an overflow happens, the result will be $(a+b)\mod
2^{3}$.

The only quantum gates using the quantum register $\left|b\right\rangle$ are
the controlled-phase gates. Moreover, they only use the qubits of the right-
hand-side register $\left|b\right\rangle$ as controls.

In the case of a constant value of $b$ known at generation time, we can
replace each controlled-phase gate by either a phase gate if the corresponding
bit of $b$ is $1$ or by an identity gate (or a “no-op” gate) if the bit of $b$
is $0$. Once this transformation has been performed, the quantum register
$\left|b\right\rangle$ is no longer used and can be safely removed from the
circuit. ∎

The final quantum add gate implementation is depicted in Figure 14, requires
$\mathcal{O}\left(n^{2}\right)$ gates and has a depth of
$\mathcal{O}\left(n\right)$. Following Barenco  _et al._ (1996); Draper
(2000); Cleve and Watrous (2000), the asymptotic gate count can be improved to
$\mathcal{O}\left(n\log(n)\right)$ by removing the rotation with an angle
below a given threshold that depend on hardware noise.

Figure 14: Modified Draper’s adder example for $3$-qubit register
$\left|a\right\rangle$ and $3$-bit classical constant $b$. The round gates
between the two applications of the Quantum Fourier Transform (QFT gates) are
phase gates and are defined in Draper (2000). A label $b_{i}$ above a phase
gate means that the phase gate should only be applied when the $i$th bit of
$b$ is set to $1$. Note that the adder wraps on overflow, meaning that if an
overflow happens, the result will be $(a+b)\mod 2^{3}$.

####  C.3.3 The cmp gate

For the same reasons exposed in the adder implementation in Section C.3.2, the
cmp gate cannot be implemented using the arithmetic comparator presented in
Cuccaro  _et al._ (2004) because removing the right-hand side qubits seems to
be a challenging task.

Instead, we use the idea from (Cuccaro  _et al._ , 2004, Section 4.3) that
explain how to implement a comparator only by using a quantum adder. The
comparison algorithm works by computing the high-bit of the expression $a-b$.
If this high-bit is in the state $\left|1\right\rangle$ then $a<b$.

In order to compute the high-bit of $a-b$, several options are open. The two
most promising options are described in the following paragraphs.

The first option is to use a subtractor acting on $n+1$ qubits and behaving
nicely on underflows (i.e. underflows result in cycling to the highest-value),
as illustrated in Figure 15. This approach requires $2$ calls to the
subtractor and $1$ additional $2$-qubit quantum gate.

Figure 15:  Computation of the high-bit of $a-b$ with a $n+1$-qubit
subtractor. The second quantum register is an ancilla qubit that is appended
to the quantum register storing $\left|a\right\rangle$ in order to form a
$n+1$-qubit register. The result is stored in a third quantum register as
$\left|1\right\rangle$ if $a<b$, else $\left|0\right\rangle$.

Another solution would be to use Equation (52) to change the subtraction into
an addition and then use a specialised procedure to compute the high-bit of
the addition of two numbers $a$ and $b$ ($a$ being encoded on a quantum
register and $b$ a constant). Computing the high-bit of an addition between a
quantum register and a constant can be performed with the CARRY gate
introduced in Häner  _et al._ (2016). This approach requires
$\mathcal{O}\left(n\right)$ Toffoli, CNOT and X gates.

Each of the described methods has its advantages and drawbacks.

For example, the first method crucially relies on a quantum subtractor, and
will have the same properties as the subtractor used. In our specific case, we
use the subtractor implemented with Drapper’s adder Draper (2000) as explained
in Section C.3.2, which in turn uses the quantum Fourier transform. The main
disadvantage of using the QFT when looking at practical implementation on
quantum hardware is that the QFT involves phase gates with exponentially small
angles. These gates may be implemented correctly up to a given threshold, but
very small rotation angles will inevitably not be as precise as normal
rotation angles due to the hardware limitations in precision. This problem can
be circumvented by using an approximate QFT algorithm Barenco  _et al._
(1996); Cleve and Watrous (2000) that will cut all the rotation gates that
have a rotation angle smaller than a given threshold from the generated
circuit but the algorithm will not be exact anymore (small probability of
incorrect result).

On the other hand, the CARRY gate involves only X, controlled-X and Toffoli
gates. This restriction makes this implementation more robust than the first
one to hardware errors or approximations and allows to simulate efficiently
this gate on classical hardware as it only involves classical arithmetic.

As a last word, in the future, the QFT may be implemented directly into the
hardware chips to make it more efficient because it is one of the most used
quantum procedure (and so one of the best candidate for optimisation). Taking
this possibility into account seems a little premature right now but may have
a high impact on the efficiency and precision of the first solution presented.

After summarising all the drawbacks and advantages, we decided to use the
arithmetic comparator for its linear number of gates, because it is based on
arithmetic which does not involve exponentially small rotation angles and
because the need to have $n-1$ dirty qubits to lend to the procedure is not an
issue in our implementation.

####  C.3.4 The eq gate

The last gate the oracle implementation will need is an eq gate, testing the
equality between an integer stored in a quantum register and a generation-time
constant integer.

This gate has been implemented with a multi-controlled Tofolli gate and a few
X gates before and after the control qubits of the Toffoli gates that should
be equal to $\left|0\right\rangle$. The X gates are necessary because a raw
Toffoli gate set its target qubit only when all its control are in the state
$\left|1\right\rangle$, but we want each control qubit to be equal to a
specific bit of the generation-time constant integer, which can be either
$\left|0\right\rangle$ or $\left|1\right\rangle$.

An implementation example is available in Figure 16.

|
---
| =
---
|
---

Figure 16: The eq gate implementation. X gates are applied to the second
control qubit because the only bit set to $0$ in the big-endian binary
representation of $5=101_{2}$ is at the second (middle) position.

Implementing a NOT gate controlled by $n$ qubits can be done with only one
ancilla qubit or $n-2$ garbage qubits and requires $\mathcal{O}\left(n\right)$
X, controlled-X or Toffoli gates mul (2015).

###  C.4 Oracles for $H_{1}$

As noted in Section C.2, the oracles $M_{1}$ and $S_{1}$ can be optimized by
using the fact that they can encode anything for $\left|x\right\rangle$ when
the $x$th row of $H_{1}$ is empty.

We decided to use this optimization opportunity to add regularity to the
description of the $H_{1}$ matrix. The _implemented_ matrix $H_{1}$, denoted
as $H_{1}^{\text{impl}}$, is described in Equation (53).

######  Note 8.

All indices start at $0$. The first row of a matrix has the index $0$, the
second row the index $1$ and so on. This convention is used to match Python’s
indexing that starts at $0$.

$\begin{split}H_{1}^{\text{impl}}=\begin{array}[]{r}\begin{array}[]{rrr}\overbrace{\hphantom{\begin{matrix}\vdots&0&\cdots&0\\\
\vdots&-1&\ddots&\vdots\\\ \vdots&1&\ddots&0\\\ \vdots&\ddots&\ddots&-1\\\
0&\cdots&\cdots&\cdots\\\
\end{matrix}}}^{N_{c}}&\overbrace{\hphantom{\begin{matrix}1&1&0&\cdots&0\\\
0&-1&1&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
\cdots&\cdots&\cdots&\cdots&\vdots\\\
\end{matrix}}}^{N_{c}+1}&{\overbrace{\hphantom{\begin{matrix}0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&&&\vdots\\\ \vdots&&&&&\vdots\\\ \vdots&&&&&\vdots\\\
\end{matrix}}}^{2^{q}-\left(2N_{c}+1\right)}}\\\ \end{array}\\\\[0.0pt]
\begin{array}[]{r}\rotatebox[origin={c}]{90.0}{$N_{c}$}\left\\{\vphantom{\begin{matrix}1&1&0&\cdots&0\\\
0&-1&1&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
\cdots&\cdots&\cdots&\cdots&\vdots\\\ \end{matrix}}\right.\\\\[5.0pt]
\rotatebox[origin={c}]{90.0}{$N_{c}+1$}\left\\{\vphantom{\begin{matrix}\vdots&0&\cdots&0\\\
\vdots&-1&\ddots&\vdots\\\ \vdots&1&\ddots&0\\\ \vdots&\ddots&\ddots&-1\\\
0&\cdots&\cdots&\cdots\\\ \end{matrix}}\right.\\\\[5.0pt]
\rotatebox[origin={c}]{90.0}{$2^{q}-\left(2N_{c}+1\right)$}\left\\{\vphantom{\begin{matrix}0&\cdots&\cdots&\cdots\\\
\vdots&&&\\\ \vdots&&&\\\ \vdots&&&\\\ \vdots&&&\\\ 0&\cdots&\cdots&\cdots\\\
\end{matrix}}\right.\\\
\end{array}\setcounter{MaxMatrixCols}{15}\begin{pmatrix}0&\cdots&\cdots&0&0&1&0&\cdots&0&0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&\vdots&\vdots&\ddots&\ddots&\ddots&\vdots&\vdots&&&&&\vdots\\\
\vdots&&&\vdots&\vdots&&\ddots&\ddots&0&\vdots&&&&&\vdots\\\
0&\cdots&\cdots&0&0&\cdots&\cdots&0&1&0&&&&&\vdots\\\
0&\cdots&\cdots&0&0&\cdots&\cdots&\cdots&0&\mathbf{0}&&&&&\vdots\\\
1&\ddots&&\vdots&\vdots&&&&\vdots&0&&&&&\vdots\\\
0&\ddots&\ddots&\vdots&\vdots&&&&\vdots&\vdots&&&&&\vdots\\\
\vdots&\ddots&\ddots&0&\vdots&&&&\vdots&\vdots&&&&&\vdots\\\
0&\cdots&0&1&0&\cdots&\cdots&\cdots&0&0&\cdots&\cdots&\cdots&\cdots&0\\\
0&\cdots&\cdots&0&\mathbf{0}&0&\cdots&\cdots&0&0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&&\ddots&\mathbf{0}&\ddots&&\vdots&\vdots&&&&&\vdots\\\
\vdots&&&&&\ddots&\mathbf{\ddots{}}&\ddots&\vdots&\vdots&&&&&\vdots\\\
\vdots&&&&&&\ddots&\mathbf{0}&0&0&&&&&\vdots\\\
\vdots&&&&&&&\ddots&\mathbf{0}&0&\ddots&&&&\vdots\\\
0&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&0&\mathbf{0}&0&0&\cdots&\cdots&0\\\
\end{pmatrix}\end{array}\end{split}$ (53)

According to the shape of the matrix in Equation (53), the oracle $M_{1}$
should implement the transformation

$M_{1}|x\rangle|0\rangle\mapsto\begin{cases}|x\rangle|x+(N_{c}+1)\rangle&\text{if
}x<(N_{c}+1)\\\ |x\rangle|x-(N_{c}+1)\rangle&\text{else}\\\ \end{cases}.$ (54)

$M_{1}$ can be easily implemented with the quantum circuit depicted in Figure
17.

Figure 17: Implementation of the oracle $M_{1}$. The cmp gate compare the
value of the control quantum register (interpreted as a unsigned integer) with
the parameter given (written below the cmp). If the control register is
strictly lower than the parameter, the gate set the qubit it is applied on to
$\left|1\right\rangle$. The add (resp. sub) gate used in this quantum circuit
add (resp. subtract) the value of its parameter to (resp. from) the quantum
register it is applied on only if the control qubit is in the state
$\left|1\right\rangle$.

The oracle $V$ cannot be simplified using the results from Claim 1. It should
implement the transformation written in Equation (55).

$V_{1}|x\rangle|0\rangle\mapsto\begin{cases}|x\rangle|1\rangle&\text{if
}(x<2N_{c}+1)\land(x\neq N_{c})\\\ |x\rangle|0\rangle&\text{else}\\\
\end{cases}.$ (55)

The implementation of the oracle $V_{1}$ is depicted in Figure 18.

Figure 18: Implementation of the oracle $V_{1}$. The first part, Set, sets the
qubit representing the weight to $1$ for all $\left|x\right\rangle$ such that
$x<2N_{c}+1$. As this does not correspond to the correct expression of $V$,
the second part Correct is here to set the weight register back to
$\left|0\right\rangle$ when $\left|x\right\rangle==N_{c}$. The cmp gate
compare the value of the control quantum register (interpreted as a unsigned
integer) with the parameter given (written below the cmp). If the control
register is strictly lower than the parameter, the gate set the qubit it is
applied on to $\left|1\right\rangle$. The eq gate used in this quantum circuit
test if the value of its parameters equal to the quantum register it is
applied on only if the control qubit is in the state $\left|1\right\rangle$.

The last oracle left to implement in order to be able to simulate $H_{1}$ is
$S_{1}$, the oracle encoding the signs of the non-zero entries of $H_{1}$. The
convention used to encode the sign of an entry has been taken from Ahokas
(2004) and is: a positive sign is encoded as $\left|0\right\rangle$, a
negative sign is encoded as $\left|1\right\rangle$. As shown in Equation (53),
$H_{1}$ only contains positive non-zero entries so the sign oracle $S_{1}$
should implement the simple transformation of Equation (56): the identity.

$S_{1}|x\rangle|0\rangle\mapsto|x\rangle|0\rangle$ (56)

###  C.5 Oracles for $H_{-1}$

The matrix $H_{-1}$ has less regularity than $H_{1}$, which will lead to a
more complex implementation. The _implemented_ matrix $H_{-1}$, denoted as
$H_{-1}^{\text{impl}}$, is described in Equation (57).

$\begin{split}H_{-1}^{\text{impl}}=\begin{array}[]{r}\begin{array}[]{rrr}\overbrace{\hphantom{\begin{matrix}\vdots&0&\cdots&0\\\
\vdots&-1&\ddots&\vdots\\\ \vdots&1&\ddots&0\\\ \vdots&\ddots&\ddots&-1\\\
0&\cdots&\cdots&\cdots\\\ \end{matrix}}\hskip
4.30554pt}^{N_{c}}&\overbrace{\hphantom{\begin{matrix}0&-\mathbf{0}&\ddots&&\vdots\\\
&\ddots&\ddots&\ddots&\vdots\\\ &&\ddots&\ddots&0\\\ &&&\ddots&-\mathbf{0}\\\
&&&&\ddots\\\ \cdots&\cdots&\cdots&\cdots&\cdots\\\ \end{matrix}}\hskip
2.15277pt}^{N_{c}+1}&{\overbrace{\hphantom{\begin{matrix}0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&&&\vdots\\\ \vdots&&&&&\vdots\\\ 0&&&&&\vdots\\\
\mathbf{0}&\ddots&&&&\vdots\\\ 0&-\mathbf{0}&0&\cdots&\cdots&0\\\
\end{matrix}\hskip 8.61108pt}}^{2^{q}-\left(2N_{c}+1\right)}}\\\
\end{array}\\\\[0.0pt]
\begin{array}[]{r}\rotatebox[origin={c}]{90.0}{$N_{c}$}\left\\{\vphantom{\begin{matrix}1&1&0&\cdots&0\\\
0&-1&1&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\
\cdots&\cdots&\cdots&\cdots&\vdots\\\ \end{matrix}}\right.\\\\[5.0pt]
\rotatebox[origin={c}]{90.0}{$N_{c}+1$}\left\\{\vphantom{\begin{matrix}\vdots&0&\cdots&0\\\
\vdots&-1&\ddots&\vdots\\\ \vdots&1&\ddots&0\\\ \vdots&\ddots&\ddots&-1\\\
0&\cdots&\cdots&\cdots\\\ \end{matrix}}\right.\\\\[5.0pt]
\rotatebox[origin={c}]{90.0}{$2^{q}-\left(2N_{c}+1\right)$}\left\\{\vphantom{\begin{matrix}0&-\mathbf{0}&\ddots&&\vdots\\\
&\ddots&\ddots&\ddots&\vdots\\\ &&\ddots&\ddots&0\\\ &&&\ddots&-\mathbf{0}\\\
&&&&\ddots\\\ \cdots&\cdots&\cdots&\cdots&\cdots\\\ \end{matrix}}\right.\\\
\end{array}\setcounter{MaxMatrixCols}{15}\begin{pmatrix}0&\cdots&\cdots&0&1&0&\cdots&\cdots&0&0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&\vdots&0&-1&\ddots&&\vdots&\vdots&&&&&\vdots\\\
\vdots&&&\vdots&\vdots&\ddots&\ddots&\ddots&\vdots&\vdots&&&&&\vdots\\\
0&\cdots&\cdots&0&0&\cdots&0&-1&0&\vdots&&&&&\vdots\\\
1&0&\cdots&0&0&\cdots&\cdots&\cdots&0&\vdots&&&&&\vdots\\\
0&-1&\ddots&\vdots&\vdots&&&&\vdots&\vdots&&&&&\vdots\\\
\vdots&\ddots&\ddots&0&\vdots&&&&\vdots&\vdots&&&&&\vdots\\\
\vdots&&\ddots&-1&0&&&&\vdots&\vdots&&&&&\vdots\\\
0&\cdots&\cdots&0&-\mathbf{0}&\ddots&&&\vdots&0&\cdots&\cdots&\cdots&\cdots&0\\\
0&\cdots&\cdots&\cdots&0&-\mathbf{0}&\ddots&&\vdots&0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&&&\ddots&\ddots&\ddots&\vdots&\vdots&&&&&\vdots\\\
\vdots&&&&&&\ddots&\ddots&0&\vdots&&&&&\vdots\\\
\vdots&&&&&&&\ddots&-\mathbf{0}&0&&&&&\vdots\\\
\vdots&&&&&&&&\ddots&-\mathbf{0}&\ddots&&&&\vdots\\\
0&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&0&-\mathbf{0}&0&\cdots&\cdots&0\\\
\end{pmatrix}\end{array}\end{split}$ (57)

Following the placement of the non-zero and the $\mathbf{0}$ or $-\mathbf{0}$
entries in the matrix $H_{-1}^{\text{impl}}$ of Equation (53), the oracle
$M_{-1}$ should implement the transformation

$M_{-1}|x\rangle|0\rangle\mapsto\begin{cases}|x\rangle|x+N_{c}\rangle&\text{if
}x<N_{c}\\\ |x\rangle|x-N_{c}\rangle&\text{else}\\\ \end{cases}.$ (58)

This transformation is quite similar to the one implemented by the oracle
$M_{1}$ in Equation (54): $N_{c}-1$ from the transformation of $M_{1}$ has
been replaced by $N_{c}$ in the transformation of $M_{-1}$. Thanks to this
similarity, the implementation of $M_{-1}$ will be a nearly-exact copy of the
implementation of $M_{1}$. The full implementation of the $M_{-1}$ oracle is
depicted in Figure 19.

Figure 19: Implementation of the oracle $M_{-1}$. The cmp gate compare the
value of the control quantum register (interpreted as a unsigned integer) with
the parameter given (written below the cmp). If the control register is
strictly lower than the parameter, the gate set the qubit it is applied on to
$\left|1\right\rangle$. The add (resp. sub) gate used in this quantum circuit
add (resp. subtract) the value of its parameter to (resp. from) the quantum
register it is applied on only if the control qubit is in the state
$\left|1\right\rangle$.

The weight oracle $V_{-1}$ is the simplest to implement for the matrix
$H_{-1}$, even if it cannot take advantage of the optimisation discussed in
Claim 1. The transformation that should be implemented by the oracle $V_{-1}$
is shown in Equation (59).

$V_{-1}|x\rangle|0\rangle\mapsto\begin{cases}|x\rangle|1\rangle&\text{if
}x<2N_{c}\\\ |x\rangle|0\rangle&\text{else}\\\ \end{cases}.$ (59)

The implementation of the weight oracle $V_{-1}$ is illustrated in Figure 20.

Figure 20: Implementation of the oracle $V_{-1}$. The cmp gate compare the
value of the control quantum register (interpreted as a unsigned integer) with
the parameter given (written below the cmp). If the control register is
strictly lower than the parameter, the gate set the qubit it is applied on to
$\left|1\right\rangle$.

The last oracle left to implement is $S_{-1}$, the sign oracle. Due to the
sign irregularity in the matrix $H_{-1}^{\text{impl}}$, the implementation of
$S_{-1}$ is more involved and requires several ancillary qubits. According to
the shape of the matrix $H_{-1}^{\text{impl}}$, the sign oracle $S_{-1}$
should implement the transformation defined in Equation (60).

$S_{-1}|x\rangle|0\rangle\mapsto\begin{cases}|x\rangle|0\rangle&\text{if
}(x=0)\lor(x=N_{c})\\\ |x\rangle|1\rangle&\text{else}\\\ \end{cases}.$ (60)

An implementation of the oracle $S_{-1}$ is illustrated in Figure 21.

Figure 21: Implementation of the oracle $S_{-1}$. The eq gate used in this
quantum circuit test if the value of its parameters equal to the quantum
register it is applied on only if the control qubit is in the state
$\left|1\right\rangle$. The OR gate flips the target qubit if and only if at
least one of the two control qubits is in the state $\left|1\right\rangle$.

##  Appendix D Note on precision

The Hamiltonian simulation algorithm we used to implement the wave equation
solver still suffer from loose bounds on $r$, the number of time-step we
should divide the evolution in in order to ensure that the error is at most
$\epsilon$.

Several bounds exist to determine a $r\in\mathbb{N}^{*}$ that will
analytically ensure that the maximum allowable error $\epsilon$ is not
exceeded. The definition of such bounds can be found in (Childs  _et al._ ,
2018, Appendix F) and is rewritten in the following equations.

The first bound has been devised by analytically bounding the error of
simulation due to the Trotter-Suzuki formula approximation by $\epsilon_{0}$

$\left|\left|\exp\left(-it\sum_{j=0}^{m-1}H_{j}\right)-\left[S_{2k}\left(-\frac{it}{r}\right)\right]^{r}\right|\right|\leqslant\epsilon_{0}$
(61)

and then let $\epsilon_{0}\leqslant\epsilon$ for a given desired precision
$\epsilon$. If we let $\Lambda=\max_{j}||H_{j}||$ and

$\tau=2m5^{k-1}\Lambda|t|$ (62)

then

$r^{ana}_{2k}=\left\lceil\max\left\\{\tau,\sqrt[2k]{\frac{e\tau^{2k+1}}{3\epsilon}}\right\\}\right\rceil.$
(63)

This bound is called the _analytic bound_.

A better bound called the _minimised bound_ can be devised by searching for
the smallest possible $r$ that satisfies the conditions detailed in (Childs
_et al._ , 2018, Propositions F.3 and F.4). This bound is rewritten in
Equation (64).

$r^{min}_{2k}=\min\left\\{r\in\mathbb{N}^{*}:\frac{\tau^{2k+1}}{3r^{2k}}\exp\left(\frac{\tau}{r}\right)<\epsilon\right\\}$
(64)

The main issue with the analytic and minimised bounds is their looseness: even
if the minimised bound improves over the analytic bound, both bounds can be
orders of magnitudes higher from the optimal bound for some Hamiltonian
matrices $H$. The analytic and minimised bounds are compared to the empiric
bound, i.e. the optimal value of $r$, in Figure 22 for the matrix
$H=H_{0}+H_{1}$ as defined in Equation (39). From this comparison we see that
in our specific case of $H=H_{0}+H_{1}$, the minimised bound is quite close to
the ideal empiric bound, which means that the overall performance of the
implementation will not be improved a lot by a more precise bound.

Figure 22: Plot of the repetition number $r$ computed for the matrix
$H_{0}+H_{1}$ with the analytic bound of Equation (63), the minimised bound of
Equation (64) and the empirical bound which is the optimal value of $r$,
computed by brute-force for small matrix sizes. This value of $r$ has been
computed for the product-formula introduced in Equation (26) (i.e. the order
$k=1$), a simulation time of $t=1$ and a target precision of
$\epsilon=10^{-5}$.

##  Appendix E Optimisation of the implementations

Both the Hamiltonian simulation and the wave equation solver circuits have
been presented unoptimised. In this appendix we study possible optimisations
that may be used on the implemented circuits to improve either the depth of
the circuit or the number of gates used.

######  Note 9.

As the wave equation solver uses the Hamiltonian simulation routine
underneath, all the optimisations performed on the later will affect the
former.

The methodology used to optimise a quantum program is the same as in classical
computing:

  1. 1.

Define the quantity we want to optimise. Most of the time, both in classical
and quantum computing, the quantity of interest is the execution time. In the
specific case of quantum computing, other quantities may be interesting to
optimise, such as the number of S or T gates (for gate count optimisation when
error correction codes are used) or the maximum idling time of each qubit (for
coherence error minimisation).

  2. 2.

Isolate the worst portions of code with respect to the quantity chosen in the
first step. If we are interested in the total execution time, the best
candidates for optimisation are the portions of code that take the most time
to execute.

  3. 3.

Improve the isolated portions of code with respect to the optimised quantity.

  4. 4.

Return at step 2 until the overall code reach the target goal (for example
“execution time divided by 2”) or the improvements made in the previous step
are too small when compared to the effort delivered to optimise.

We followed this scheme to determine which subroutines are good candidates for
optimisation. Because quantum program are ultimately bound to be run on
quantum hardware, we chose to optimise the execution time on a real quantum
chip. Only a few constructor give access to the characteristics of their chips
and one of them is IBM. All the characteristics of IBM Q chips are available
online on ibm (2019c) and via Qiskit API. We chose to use IBM Q 16 Melbourne
V1.3.0 because it is the largest chip in terms of number of qubits that was
publicly available at the time of writing.

Once the quantity to optimise has been clearly set, the next step is to
determine the most interesting routines to optimise. In classical computing, a
lot of tools exists to gather runtime informations and visualise this data.
One of the most famous tool is probably gprof which is able to approximate the
time spent in each function and create a call graph. Sadly, this kind of tool
is still non-existent in the field of quantum computing. That is why we
created a quantum equivalent of gprof, called qprof, that outputs exactly the
same data as gprof but for a quantum program.

The text output of qprof is enough to isolate the subroutines to optimise
first, but a visual representation can be obtained by combining qprof with a
gprof-compatible visualisation tool such as gprof2dot gpr (2019). The output
of qprof for the Hamiltonian simulation procedure implemented can be found in
Figure 23. The same graph can be obtained for the wave equation solver. As the
only difference between the wave equation solver and the Hamiltonian
simulation routine of Figure 23 is the time of evolution (and so the number of
repetitions $r$), the produced graphs are the same except for some routines
where a negligible difference of $0.1\%$ between the two call-graphs have been
observed. For this reason, we omitted the output of qprof for the quantum wave
equation solver.

###  E.1 Hamiltonian simulation procedure optimisation

Figure 23: Reduced call-graph of the Hamiltonian Simulation routine
implemented. Each node represent a quantum routine and contain the routine
name, the total runtime of the routine (time spent in subroutines included),
the self time of the routine (time spent only in the routine, not in
subroutines) and the number of calls of the routine. A routine name starting
with D- means that the routine has been inverted (D- is for dagger). All the
nodes (i.e. sub-routines) taking less than 10% of the total time have been
ignored. Instruction parallelism is not taken into account (the circuit is
flattened before analysis).

Figure 23 represent the simplified call graph of the Hamiltonian simulation
procedure. A lot of information can be extracted from this call graph:

  1. 1.

The trotter_suzuki_formula procedure is called $r_{0}=1463$ times. Dividing
the number of repetitions by a factor of $k$ would divide the number of gates
and the execution time by the same factor $k$.

  2. 2.

More than $85\%$ of the total execution time is spent in the oracles.

  3. 3.

More than $75\%$ of the total execution time is spent in the
arithmetic_compare procedure.

  4. 4.

More than half of the total execution time is spent in the high_bit_compute
procedure.

From these observations, several optimisation strategies can be devised. Some
are discussed in the following sections.

####  E.1.1 Optimising the number of repetitions $r$

One of the major cost of the algorithm comes from the number of repetitions
$r$ we need to perform in order to ensure a given precision $\epsilon$.
Several strategies can be used to try to reduce the number of repetitions $r$.

A first strategy could be to find a better upper-bound for $r$ than
$r^{ana}_{2k}$ (see Equation (63)) and $r^{min}_{2k}$ (see Equation (64)).
Such a bound can be found in Childs  _et al._ (2019). This strategy could be
interesting if the gap between the empiric and the minimised bounds was bigger
in Figure 22. As the bound $r^{min}_{2k}$ seems to tend towards $r^{emp}$, we
believe that the potential gain of this optimisation will not be significant
enough.

Another way of improving the number of repetitions $r$ would be to use a more
precise approximation, for example by using $S_{2k}(\lambda)$ with $k>1$. But
using a higher-order approximation also increase the number of calls to the
oracles by a factor $5^{k}$.

####  E.1.2 Optimising the subroutines

The call graph in Figure 23 allows us to isolate the procedures that
contribute the most to the total execution time: high_bit_compute with more
than $55\%$ of the total runtime.

The choice we made in Section C.3.3 to compare a quantum register and a
constant classical value with a procedure computing the high-bit of a
subtraction turns out to be a crucial step for the overall solver
performances.

One idea of optimisation would be to use another algorithm to compare a
quantum register with a constant value. The procedure using Drapper’s adder in
Section C.3.3 might be a good candidate on quantum hardware that can perform
the QFT efficiently for example.

####  E.1.3 LCR optimisation

In parallel of all the other optimisation methods, automatic quantum circuit
optimisers can be applied to all or parts of the quantum wave equation solver.

One of the most promising optimisation procedure for our wave equation solver
is probably the LCR method used with an automatic circuit optimisation
algorithm, both of which can be found in Nam  _et al._ (2018). The LCR method
is specifically designed to optimise quantum circuits that consist in the
repetition of the same circuit a number $N\geqslant 3$ times.