,

  *

# Quantum algorithms for scientific applications

R. Au-Yeung1, B. Camino2, O. Rathore3 and V. Kendon1 1 Department of Physics,
University of Strathclyde, Glasgow G4 0NG, United Kingdom  2 Department of
Chemistry, UCL, London WC1E 6BT, United Kingdom  3 Department of Physics,
Durham University, Durham DH1 3LE, United Kingdom  rhonda.au-
<EMAIL_ADDRESS><EMAIL_ADDRESS>0000-0002-0082-5382
0000-0002-9569-2219 0000-0002-6551-3056

###### Abstract

Quantum computing promises to provide the next step up in computational power
for diverse application areas. In this review, we examine the science behind
the quantum hype and breakthroughs required to achieve true quantum advantage
in real world applications. Areas that are likely to have the greatest impact
on high performance computing (HPC) include simulation of quantum systems,
optimisation, and machine learning. We draw our examples from materials
simulations and computational fluid dynamics which account for a large
fraction of current scientific and engineering use of HPC. Potential
challenges include encoding and decoding classical data for quantum devices,
and mismatched clock speeds between classical and quantum processors. Even a
modest quantum enhancement to current classical techniques would have far-
reaching impacts in areas such as weather forecasting, engineering, aerospace,
drug design, and realising “green” materials for sustainable development. This
requires significant effort from the computational science, engineering and
quantum computing communities working together.

###### type:

Review Article

Keywords: quantum algorithms, quantum computing, scientific computing

###### Contents

  1. 1 Introduction
    1. 1.1 Scope of this review
  2. 2 Designing quantum algorithms
    1. 2.1 Universal quantum computers
    2. 2.2 Quantum computing stack
    3. 2.3 Foundational quantum algorithms
      1. 2.3.1 Quantum Fourier transform
      2. 2.3.2 Quantum phase estimation
      3. 2.3.3 Shor’s algorithm
      4. 2.3.4 Grover’s algorithm
      5. 2.3.5 Quantum amplitude amplification
      6. 2.3.6 Quantum walks
      7. 2.3.7 Harrow-Hassidim-Lloyd algorithm
    4. 2.4 Quantum optimisation
    5. 2.5 Variational quantum algorithms
      1. 2.5.1 Quantum approximate optimisation algorithm
      2. 2.5.2 Variational quantum eigensolvers
      3. 2.5.3 Variational quantum simulation
      4. 2.5.4 Quantum neural networks and parameterised quantum circuits
      5. 2.5.5 Limitations of VQAs
    6. 2.6 Performance benchmarks: quality, speed, scale
      1. 2.6.1 Low- and high-level benchmarks
      2. 2.6.2 Quantifying quantum device performance
      3. 2.6.3 Quantum circuit simulation methods
      4. 2.6.4 Measuring quantum advantage
  3. 3 Quantum simulation of quantum systems
    1. 3.1 Hamiltonian evolution
    2. 3.2 Quantum algorithms for Monte Carlo
    3. 3.3 Quantum computational chemistry
      1. 3.3.1 Quantum chemistry’s many-electron problem
      2. 3.3.2 HF and the electron correlation problem
      3. 3.3.3 Density functional theory (DFT)
      4. 3.3.4 Post Hartree-Fock methods
      5. 3.3.5 Post Hartree-Fock quantum algorithms
      6. 3.3.6 Quantum optimisation for chemistry and materials science
    4. 3.4 High-energy physics quantum simulations
    5. 3.5 Important applications in industry
      1. 3.5.1 High-temperature superconductivity
      2. 3.5.2 Fertilisers
  4. 4 Quantum algorithms for classical simulation
    1. 4.1 Solving the Navier-Stokes equations
    2. 4.2 Lattice Boltzmann method
    3. 4.3 Quantum lattice Boltzmann
    4. 4.4 Smoothed particle hydrodynamics
    5. 4.5 Quantum SPH algorithm
    6. 4.6 Other quantum CFD algorithms
  5. 5 Interfacing with high performance computing
    1. 5.1 Clock speed mismatch
    2. 5.2 Quantum data encoding
    3. 5.3 Data interconnections
  6. 6 Outlook
  7. ORCID

##  1 Introduction

The most advanced quantum computers are now approaching exascale classical
computing capabilities. For most scientific and engineering applications,
“exascale” means the system’s aggregate performance crosses a threshold of
$10^{18}$ IEEE 754 double precision (64-bit) operations (multiplications
and/or additions) per second (exaFLOPS). On the other hand, machines just
below exascale are used to verify quantum supremacy experiments. For example,
random circuit sampling and Gaussian boson sampling methods respectively used
the Summit ($\sim$150 petaFLOPs) Arute et al. [2019] and Sunway TaihuLight
($\sim$93 petaFLOPs) Li et al. [2022a] supercomputers. To compare, the Sony
PlayStation 5 digital edition is listed as having a peak of $\sim$10 teraFLOPS
($\sim 10^{12}$ FLOPS). These supremacy tasks are not yet known to have any
industrial or other useful applications. However, they do indicate that a
quantum co-processor running a well-chosen sub-task could provide a huge
increase in exascale hardware capabilities.

Figure 1: Microprocessor trend data, 1980-2021. Data up to year 2010 collected
by M. Horowitz, F. Labonte, O. Shacham, K. Olukotun, L. Hammond, and C.
Batten, data for 2010-2021 by K. Rupp. Both datasets taken from Rupp [2022].

Reaching for exascale computing calls for a paradigm shift from traditional
computing advances Mann [2020], Betcke et al. [2022]. Historical advances in
computing were achieved via increasing clock speed and memory performance,
adding more cores and parallelisation. Moore’s law Moore [1965] predicts that
the exponential growth of silicon-based transistors on computer chips is about
to hit physical limits Khan et al. [2018], Leiserson et al. [2020]. The growth
predicted was borne out by major chip and computer system vendors until about
2005, when CPU clock speeds reached around 4 GHz (Figure 1). Significantly
increasing the clock speed beyond this would require enormous effort to cool
the processor to prevent malfunction or permanent hardware damage from
overheating.

Launched in 1993, the TOP500 project assembles and maintains a list of the
world’s 500 most powerful computer systems Strohmaier et al. [2023a]. As of
November 2023, Oak Ridge National Laboratory’s Frontier system is the only
true exascale machine worldwide. In recent years, GPUs have enabled a step
change in processing power – all of the GREEN500 Strohmaier et al. [2023b] use
them. Only Fugaku in the top ten of the current TOP500 list Strohmaier et al.
[2023a] does not use GPUs. Physical constraints and power requirements are
forcing computer vendors to develop new strategies to achieve more compute
power. Novel hardware systems, including quantum computing, will be required
to significantly extend global computational capacity beyond current
capabilities.

The idea of using quantum systems to process information more efficiently than
classical von Neumann computers was introduced some 40 years ago Feynman
[1982], Deutsch [1985]. Since then there has been steady progress in quantum
computing, building on the initial achievements of Shor’s algorithm for
factorisation Shor [1994] that could break current methods of encryption, and
Grover’s algorithm Grover [1996, 1997] for searching unsorted data. Quantum
simulation is another promising field of quantum information processing
Georgescu et al. [2014], Altman et al. [2021]. Many near-term applications of
quantum computers fall under this umbrella. Quantum simulation involves
modelling the quantum properties of systems that are directly relevant to
understanding modern materials science, high-energy physics and quantum
chemistry. As the practical relevance of quantum computing becomes clearer,
interest has grown beyond the confines of academia and many countries now have
national strategies to develop quantum computing and quantum technology more
generally, e.g., Department for Science, Innovation and Technology [2023],
Acín et al. [2018], Monroe et al. [2019].

The consensus is that quantum devices will not simply replace HPC. It makes
little sense (in terms of engineering and economics) to use quantum methods
for problems that classical computers already do well. Hence, a hybrid
solution is the most efficient, cost-effective and productive way to approach
quantum computation. Instead of simply sending jobs to one or more CPUs,
hybrid architectures delegate different parts of the problem to different
types of processors, the most common currently being GPUs. Conceptually,
quantum processor units (QPUs) are a natural extension to enhance the
processing power of HPC. However, increasing the number and diversity of
processors compounds the challenges for efficient programming and scheduling.
Quantum computing hardware is not yet mature enough to seamlessly integrate
with HPC.

Figure 2: Quantum computing commercial landscape. Some examples of companies
(blue) partnered with quantum tech companies (orange) to advance QC research.
Cambridge Quantum Computing (CQC) and Honeywell now merged as Quantinuum
(red). We include CQC, Honeywell and Quantinuum separately as some
collaborations were formed before the merger in Dec 2021. Universities,
research institutes and government organisations not shown. Compiled by R. Au-
Yeung from press releases and publicly disclosed interactions with end users.

After several decades of largely academic-based theoretical and experimental
development, growing commercial interest has resulted in venture capital
MacQuarrie et al. [2020], Gibney [2019], Bova et al. [2021] providing funding
for many existing and start-up companies. The current landscape shows growing
collaborations between quantum start-ups, established big players, and with
well-established financial, pharmaceutical and automotive companies who wish
to integrate quantum technologies into their development and production
processes (Figure 2). Notable players include Google, IBM and Rigetti, who are
developing superconducting qubit systems, and Intel who are building systems
based on quantum dots, with chips that can be made in its existing foundries.
Cloud providers such as Amazon Web Services (AWS), are now selling time on
different quantum systems, including QuEra, which uses neutral atoms in
tweezer arrays, IonQ ion trap quantum processors, and superconducting qubit
systems from Rigetti. Current systems are small and imperfect, known as noisy,
intermediate scale quantum computers (NISQ) Preskill [2018], Bharti et al.
[2022]. They are not yet powerful enough to solve useful problems, but do
allow researchers to test and develop quantum algorithms and applications to
prepare for more advanced quantum hardware.

Several commercial partnerships are focused on integrating quantum processors
with HPC. Finnish quantum start-up IQM is mapping quantum applications and
algorithms directly to quantum processors to develop application-specific
superconducting computers IQM [2023]. The result is a quantum system optimised
to run particular applications such as HPC workloads. The Atos-Pasqal
collaboration are working to incorporate Pasqal’s neutral-atom quantum
processors into HPC environments Atos [2020]. NVIDIA’s cuQuantum software
development kit (SDK) provides the tools needed to integrate and run quantum
simulations in HPC environments using GPUs McCaskey [2022], with the goal of
adding QPUs to the SDK. Quantum software vendor Zapata anticipates the
convergence of quantum computing, HPC, and machine learning, and has created
the Orquestra platform to develop and deploy generative AI applications Zapata
AI [2023].

###  1.1 Scope of this review

Given the growing government and commercial investment in the development of
quantum computing, it is timely to consider how it might accelerate
computational science and engineering. Now considered the third pillar of
scientific research alongside theoretical and experimental science, scientific
computing has become ubiquitous in scientific and engineering applications. It
leverages the exponential growth in computing power of the last decades, and
receives even more government investment than quantum technology, given its
essential role in science, engineering, and innovation.

In the UK, HPC investments by the Engineering and Physical Sciences Research
Council (EPSRC) total £466 million over 10 years up to 2019 Herr et al.
[2019], Department for Science, Innovation and Technology [2022]. This is
anticipated to add between £3 billion and £9.1 billion to the UK economy. The
European High-Performance Computing Joint Undertaking (EuroHPC JU) is a €7
billion joint initiative between the EU, European countries and private
partners to develop a HPC ecosystem in Europe European Commission [2023].
These examples show that governments are willing to invest significant
resources into developing exascale technologies that can provide numerous
societal and economic benefits.

Quantum computing has a broad reach across applications, so this review
necessarily has to select which areas to focus on. We choose two important
applications, namely, materials science and fluid dynamics simulations. In
computational chemistry, the hard problems are fully quantum and therefore
natural for quantum computers Feynman [1982]. In fluid dynamics, the problem
is purely classical. The hardness stems from the need to solve nonlinear
differential equations. It is less clear that there will be substantial
quantum advantages here, although there are multiple interesting proposals for
possible quantum algorithms. Our chosen applications make up a significant
fraction of current HPC use. Between them they are diverse enough to offer
insights for other areas of computational science and engineering simulations.
Although we discuss quantum machine learning techniques relevant for chemistry
and fluid simulations, the wider field of “big data” and artificial
intelligence (AI) are outside the scope of this review. There are indeed
potential quantum enhancements in AI, but they merit their own dedicated
reviews (see, e.g., Dunjko and Briegel Dunjko and Briegel [2018]).

In this review, we begin by describing the concepts behind quantum computers
and giving examples of important foundational quantum algorithms in section 2.
For conciseness, we minimise the amount of mathematical detail where possible.
Interested readers are encouraged to consult references such as Nielsen and
Chuang [2010], Abhijith et al. [2022] for details and refer to the original
papers for proofs. There are already many excellent reviews that cover
different classes of quantum algorithms (see e.g., Montanaro [2016], Albash
and Lidar [2018], Childs and van Dam [2010]), including a “Quantum Algorithm
Zoo” website Jordan [2021] which cites over 400 papers. In this review, we aim
to emphasise the physical mechanisms and qualitative insights into quantum
algorithms and their real-world applications. We will not discuss abstract
topics such as complexity analysis Bernstein and Vazirani [1993], Vazirani
[2002] in depth. We also mention alternative quantum computing methods, such
as quantum annealing (section 2.4) for chemistry and materials science
simulations (section 3.3.6) and variational quantum algorithms (section 2.5).
We outline some of the many proposed performance benchmarks for quantum
computers in section 2.6.

Next we focus on our two example fields: quantum simulation of quantum
chemical systems (section 3), and quantum algorithms for classical fluids
simulation (section 4). We summarise the current (classical) methods and
shortcomings, then explore how quantum computers may improve performance. With
the theoretical potential established, we turn to the practical challenges of
how to combine quantum components into HPC architectures in section 5. This
includes the clock speed mismatches. and how to encode classical data into
quantum states. Finally, we summarise and set out future research directions
in section 6.

There are other important topics on the road to practical quantum computing
that already have their own reviews. These include: quantum error correction
Terhal [2015], Roffe [2019]; verification and testing Eisert et al. [2020],
Carrasco et al. [2021] (we only cover benchmarks in section 2.6); and
comparisons between different hardware and software platforms available on the
market MacQuarrie et al. [2020]. Quantum annealing has been well-reviewed
recently: Hauke et al. Hauke et al. [2020] cover methods and implementations
while Yarkoni et al. Yarkoni et al. [2022] focus on industrial applications of
optimisation problems.

Despite numerous publications hyping the quantum revolution (see, e.g., Arute
et al. [2019], Wu et al. [2021], Kim et al. [2023]), it pays to read the fine
print Leymann and Barzen [2020], Coveney and Highfield [2020], Aaronson
[2015]. The methodology in Chancellor et al. [2020] is an instructive example
of how to fairly and objectively evaluate different use cases. This helps to
avoid evaluations that focus on the positives and ignore weaknesses.

A key message from the IQM-Atos 2021 white paper Perini and Ciarletta [2021]
is that HPC centres must create a long-term strategic plan to successfully
integrate quantum computing into their workflow. In particular, they recommend
establishing a three-step roadmap:

  1. 1.

Gap analysis and quantum solution identification (now). Where are the
bottlenecks in classical HPC? How can we improve the solution accuracy?

  2. 2.

Quantum solution design and integration (mid-term). What quantum methods can
we use to resolve the issues identified in gap analysis? How do we integrate
this into HPC architectures?

  3. 3.

Quantum computing use-case development and implementation (long-term). Which
applications and problems can get the most benefit from quantum speedup?

Our review provides an introduction to the tools and topics needed to make
such plans.

##  2 Designing quantum algorithms

Developing useful quantum algorithms is extremely challenging, especially
before quantum computing hardware is widely available. It is like building a
car without having a road or fuel to take it for a test drive. Since quantum
computing has a different logical basis from classical computing, it is not as
straightforward as taking a successful classical algorithm and mapping it to
the quantum domain. For example, a classical CNOT gate makes a copy of a
classical bit. In quantum theory, the no-cloning theorem Wootters and Zurek
[1982] states that it is impossible to make copies of unknown quantum states.

Table 1: Notation in our paper. Notation |  Description
---|---
$a$, $A$ |  scalar, can be real or complex
$a(t)$, $A(t)$ |  scalar function defined in terms of variable $t$
$\vec{a}$ |  vector, contains elements $\vec{a}_{0}$, $\vec{a}_{1}$, …
$\mathbf{A}$ |  matrix
$\ket{a}$, $\ket{A}$ |  quantum state vector
$\hat{A}$ |  quantum operator
$j$, $k$ |  index

In this section, we outline the basic concepts underpinning quantum computing.
Then, we describe some existing quantum algorithms that can potentially be
adapted for practical scientific applications. To avoid confusion, we use the
notation in Table 1 to differentiate between different quantities. These
include operators (hats), matrices (bold font), vectors (arrows), scalars and
functions (italics), and bra-ket notation for quantum states.

###  2.1 Universal quantum computers

The most widely used model for quantum computing is the gate or circuit model.
The basic unit of quantum information is the quantum bit (qubit), a two-state
system. Qubits are the quantum analogue of classical bits. They are two-level
quantum systems that represent linear combinations of two basis states

$\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$ (1)

where $\alpha,\beta\in\mathbb{C}$ and $|\alpha|^{2}+|\beta|^{2}=1$. In vector
notation, we usually choose

$\ket{0}=\begin{bmatrix}1\\\ 0\end{bmatrix},\quad\ket{1}=\begin{bmatrix}0\\\
1\end{bmatrix}.$ (2)

Multiple qubits can be entangled, forming systems with an exponentially large
state space of size $N=2^{n}$ for $n$ qubits. This corresponds to the tensor
product of matrices, for example,

$\ket{0}\otimes\ket{0}=\begin{bmatrix}1\\\ 0\\\ 0\\\ 0\end{bmatrix}.$ (3)

The superposition property means that an $N$-qubit system can represent
$2^{N}$ states whereas classical bits represent only one state at a time (a
single $N$-bit string).

Figure 3: Bloch sphere representation of qubit state $\ket{\psi}$ with
examples of gate operations.

We visualise single qubits using a Bloch sphere (Figure 3). To perform
calculations, we manipulate the qubit state using (reversible) gate
operations. These are analogous to the classical logic gates, such as AND, OR
and NOT.

Table 2: Examples of commonly used quantum gates. Middle column uses circuit notation where horizontal lines represent qubits. Boxes, crosses and dots indicate the gate applied. Time flows from left to right. Right-hand column shows the matrix representation. One-qubit gate |  Circuit symbol |  Matrix
---|---|---
Pauli-X or NOT (bit flip) |  |  $\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}$
Pauli-Z (phase flip) |  |  $\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix}$
Hadamard (create superposition) |  |  $\dfrac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\ 1&-1\end{bmatrix}$
T gate |  |  $\begin{bmatrix}1&0\\\ 0&e^{i\pi/4}\end{bmatrix}$
Phase shift |  |  $\begin{bmatrix}1&0\\\ 0&e^{i\vartheta}\end{bmatrix}$
Z-rotation |  |  $\begin{bmatrix}e^{-i\varphi/2}&0\\\ 0&e^{i\varphi/2}\end{bmatrix}$
Two-qubit gates |  Circuit symbol |  Matrix
CNOT (flip target qubit if control qubit is $\ket{1}$) |  |  $\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&1\\\ 0&0&1&0\end{bmatrix}$
SWAP (swap two qubit states) |  |  $\begin{bmatrix}1&0&0&0\\\ 0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\end{bmatrix}$

Gate operations are unitary operators applied to the qubits. Applying such an
operator $\hat{U}$ to a qubit state, then its conjugate transpose
$\hat{U}^{\dagger}$ brings the qubit back to its original state
($\hat{U}\hat{U}^{\dagger}=\mathbb{1}$). Gate operations can be written as
matrices (see examples in Table 2. In quantum computing, we use a convenient
graphical way to depict the unitary transformations (gate operations) on a one
or $n$-qubit system. The building blocks of the transformations are $2\times
2$ and $4\times 4$ unitary matrices on single qubits and two qubits
respectively. As in classical computing, any computation can be built up from
a small set of universal gates. A universal gate set includes at least one
entangling gate acting on two or more qubits, and one or more single qubit
rotations.

###  2.2 Quantum computing stack

Figure 4: Quantum computer architecture among some sub-fields of quantum
computation. Reproduced from van Meter and Horsman [2013]

As with classical computers, we can decompose the quantum computer
architecture in terms of its various levels of abstraction, called the
“stack”. In gate-based models, this extends from the user interface and
compiler down to the low-level gate operations on the physical hardware itself
(Figure 4).

At the highest levels, quantum algorithms give instructions for solving the
numerical problem. Compiling and translating quantum gates in the middle
levels (software stack) compress the algorithms to accelerate performance,
while error-correction techniques mitigate quantum computing errors. The lower
levels comprise the quantum computer’s building blocks (hardware) for
manipulating and storing qubits Alexeev et al. [2021]. Note that there are two
classes of qubits: physical and logical qubits Kjaergaard et al. [2020].

Inside a quantum computer, physical qubits are physical devices that behave as
a two-level quantum system, whereas logical qubits are groups of physical
qubits that acts as a single qubit in a quantum algorithm. In reality,
physical qubits often operate in noisy environments. Many physical qubits are
needed for quantum error correction procedures Roffe [2019] to produce one
logical qubit to perform useful computations. Engineering qubits in real
physical hardware is highly challenging. The requirements for scalable quantum
computing laid out by DiVincenzo DiVincenzo [2000] are still relevant today:
scalable, well-characterised qubits; the ability to initialise the qubits;
long decoherence times; a “universal” set of quantum gates; and the ability to
measure the qubits.

For example, Reiher et al. Reiher et al. [2017] examine the cost of performing
quantum simulations to find the potential energy landscape of a quantum
chemical structure (FeMoco, more details in section 3.5.2). They estimated 313
physical qubits are needed for each logical qubit, if there is one error for
every $10^{9}$ quantum gates. By extension, we would likely need thousands of
physical qubits per logical qubit, or much better error rates, for complicated
real-world quantum chemistry simulations.

In the rest of our review, we will discuss in terms of logical qubits unless
otherwise stated.

###  2.3 Foundational quantum algorithms

In this section, we introduce some foundational quantum algorithms which can
be used as components in quantum algorithms for a wide variety of
applications. They also illustrate how quantum algorithms can provide powerful
speed ups over classical algorithms for the same problems. For readers seeking
a comprehensive graduate-level introduction, we recommend Abhijith et al.
[2022], Portugal [2022], Nielsen and Chuang [2010] plus IBM’s Qiskit software
documentation Qiskit [2023].

####  2.3.1 Quantum Fourier transform

The quantum Fourier transform (QFT) is a key subroutine in many quantum
algorithms, most prominently Shor’s factoring algorithm Shor [1994]. The QFT
operation changes the quantum state from the computational (Z) basis to the
Fourier basis. It is the quantum analogue of the discrete Fourier transform
and has exponential speed-up (complexity $O((\log N)^{2})$) compared to the
fast Fourier transform’s $O(N\log N)$ complexity for problem size $N$. The QFT
performs a discrete Fourier transform on a list of complex numbers encoded as
the amplitudes of a quantum state vector. It stores the result “in place” as
amplitudes of the updated quantum state vector. Measurements performed on the
quantum state identify individual Fourier components – the QFT is not directly
useful for determining the Fourier-transformed coefficients of the original
list of numbers, since these are stored as amplitudes (see section 5.2 for
more details on quantum encoding).

When applying the QFT to an arbitrary multi-qubit input state, we express the
operation as

$\hat{U}_{\textrm{QFT}}=\frac{1}{\sqrt{N}}\sum_{j,k=0}^{N-1}e^{2\pi
ijk/N}\ket{j}\\!\bra{k}$ (4)

which acts on a quantum state in $N$-dimensional Hilbert space
$\ket{\tilde{x}}=\sum_{k=0}^{N-1}x_{k}\ket{k}$ and maps it to
$\ket{\tilde{y}}=\sum_{j=0}^{N-1}y_{j}\ket{j}$. This is a unitary operation on
$n$ qubits that can also be expressed as the complex $N\times N$ matrix

$\displaystyle\mathbf{F}_{N}=$
$\displaystyle\frac{1}{\sqrt{N}}\begin{bmatrix}1&1&1&1&...&1\\\
1&\omega&\omega^{2}&\omega^{3}&...&\omega^{N-1}\\\
1&\omega^{2}&\omega^{4}&\omega^{6}&...&\omega^{2(N-1)}\\\
1&\omega^{3}&\omega^{6}&\omega^{9}&...&\omega^{3(N-1)}\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
1&\omega^{N-1}&\omega^{2(N-1)}&\omega^{3(N-1)}&...&\omega^{(N-1)^{2}}\end{bmatrix}$
(5)

where $\omega=e^{2\pi i/N}$ and $N=2^{n}$. Note that the Hadamard gate (Table
2) acts as the single-qubit QFT. It transforms between the Z-basis states
$\ket{0}$ and $\ket{1}$ to the X-basis states

$\ket{\pm}=\frac{1}{\sqrt{2}}(\ket{0}\pm\ket{1}).$ (6)

Figure 5: Quantum Fourier transform operation on $n$ qubits. QFT =
$\mathbf{F}_{N}$, equation 5. Adapted from Camps et al. [2020].

Figure 5(a) shows the quantum circuit for the QFT operation. Each multi-qubit
gate block can be decomposed into a sequence of single or two-qubit gates. For
example, the $\mathbf{B}_{k+1}$ block is written in terms of controlled
operations involving two qubits (Figure 5(b)). The $\mathbf{B}_{k+1}$ gate is
applied to the last $k+1$ qubits. All controlled-$R$ gates are diagonal
matrices and commute,

$\mathbf{R}_{N}=\begin{bmatrix}1&0\\\ 0&e^{-2\pi i/N}\end{bmatrix}$ (7)

so the order in which they are applied does not change the outcome. The bit-
reversal permutation matrix $\mathbf{P}_{n}$ can first be written in terms of
$\lfloor n/2\rfloor$ swap gates (Figure 5(c)). Since each swap can be
implemented with three CNOT gates, $\mathbf{P}_{n}$ requires $\lfloor
3n/2\rfloor$ CNOT gates. Hence the QFT requires $O(n^{2})$ elementary gates in
total.

The semi-classical nature of the QFT circuit means it can be efficiently
simulated on classical computers Griffiths and Niu [1996], Browne [2007] when
used at either the start or end of a quantum computation. The idea is to
measure a qubit, then use the result to produce a classical signal which
controls a one-qubit transformation on the next qubit before it is measured,
and so forth.

####  2.3.2 Quantum phase estimation

Kitaev’s quantum phase estimation (QPE) algorithm Kitaev [1996] is a core
component of many quantum algorithms and an important technique in algorithm
design. While the standard procedure is based on the QFT circuit, various
improvements can reduce the dependence on QFT while retaining its accuracy by
using, for example, approximations to relax the QFT constraints Ahmadi and
Chiang [2012] or a modified algorithm from quantum metrology Ni et al. [2023].
We outline the QPE algorithm below.

The basic phase estimation problem is, given a unitary operator

$\hat{U}\ket{\psi}=\lambda\ket{\psi},$ (8)

we want to find an eigenvalue $\lambda$ of the eigenvector $\ket{\psi}$.
Because $\hat{U}$ is unitary, the eigenvalue can be expressed as
$\lambda=e^{2\pi i\theta}$ and the real aim is to estimate the eigenvalue
phase $\theta$.

Figure 6: Structure of quantum phase estimation circuit, where iQFT is the
inverse of the QFT operation.

The phase estimation algorithm contains two main steps (Figure 6). In the
circuit, the top register contains $n$ “counting” qubits in state $|0\rangle$.
The second register contains as many qubits as needed to store $\ket{\psi}$,
the eigenstate of $\hat{U}$. We apply a set of Hadamard gates to the first
register (this is equivalent to a QFT on the all zeros state) and
controlled-$\hat{U}$ operations on the second register with $\hat{U}$ raised
to successive powers of two. The final state of the first register is

$\displaystyle\frac{1}{\sqrt{2^{n}}}\sum_{k=0}^{2^{n}-1}e^{2k\pi
i\theta}|k\rangle.$ (9)

Applying the inverse QFT (iQFT) gives

$\frac{1}{\sqrt{2^{n}}}\sum_{k=0}^{2^{n}-1}e^{2k\pi
i\theta}|k\rangle|\psi\rangle\to|\tilde{\theta}\rangle|\psi\rangle$ (10)

where $\tilde{\theta}$ is an estimate for $\theta$. The final stage is to read
out the state of the first register by measuring in the computational basis.
Hence, we estimate $\theta$ to $n$ bits of precision.

Quantum phase estimation is equivalent to simulating quantum system dynamics.
Given an eigenstate $\ket{\psi}$, we can estimate the eigenvalue for
Hamiltonian $\hat{H}$ by applying $\hat{H}$ via the time-evolution operator
$\hat{U}=e^{-i\hat{H}t}$. When the algorithm implements $\hat{U}^{k}$ for
increasing powers $k$, this corresponds to evolving the system for increasing
times $kt$ Aspuru-Guzik et al. [2005]. The phases $\theta$ are eigenvalues of
time evolution operator $\hat{U}$ which we can map back on to the
eigenenergies $E$ of $\hat{H}$ Aspuru-Guzik et al. [2005],

$\displaystyle\hat{U}\ket{\psi}=e^{i\hat{H}t}\ket{\psi}=e^{2\pi
i\theta}\ket{\psi},$ (11) $\displaystyle E=2\pi\theta/t.$ (12)

Computational chemistry often requires finding the eigenenergies, as they can
help deduce many chemical properties like ionisation potential and equilibrium
constants. We discuss some of these applications in section 3.

The phase estimation circuit outputs an eigenphase as the measurement outcomes
of the top register. It also prepares the corresponding eigenstate
$\ket{\psi}$ in the lower register. This is a useful starting point for
calculating other observables besides energy. Once we have the eigenstates and
eigenenergies of interest, we can understand how a (closed) system evolves in
time by decomposing the initial state into a sum of eigenstates and evolving
each eigenstate according to the phase found in QPE. The final state is an
interferometric sum of all components.

There are two important performance metrics for the QPE: maximum run time, and
total run time over all repetitions from each circuit in the algorithm.
Maximum and total run time approximately measure the circuit depth and total
cost of the algorithm respectively. The QPE algorithm in Figure 6 is
unsuitable for early fault-tolerant quantum computers since these devices can
only implement low circuit depths and limited numbers of logical qubits Ni et
al. [2023].

####  2.3.3 Shor’s algorithm

While Feynman’s 1982 arguments Feynman [1982] formalised the idea of using
quantum computers to solve quantum problems, what really caught people’s
attention was Shor’s prime factoring algorithm in 1994 Shor [1994], Ekert and
Jozsa [1996]. It provided a way for large enough (future) quantum computers to
break RSA encryption Rivest et al. [1978], then considered unbreakable due to
the exponential computational cost of factoring large numbers. RSA encryption
works by multiplying two large prime numbers $p$ and $q$. The product $C=pq$
is so large that the computational cost of factorising – and hence recovering
the encrypted message – is impractical. Shor’s algorithm showed that factoring
could be solved by quantum computers in polynomial time.

Shor’s factoring algorithm reduces the factorisation problem into an order-
finding problem. Given the modular exponential function

$f(x)=a^{x}\text{mod}\,C,\;0<a<C$ (13)

where $a$ is a random number coprime to $C$, the order-finding problem
involves finding the period $r$ of function $f(x)$. For an even $r$ and if
$a^{r/2}\neq-1\text{mod}C$, at least one prime factor of $C$ is given by
$\text{gcd}(a^{r/2}\pm 1,C)$.

Figure 7: Quantum circuit schematics for Shor’s algorithm.

There are two main components of Shor’s algorithm: modular exponentiation
(calculating $a^{x}\text{mod}C$) and inverse quantum Fourier transform (iQFT).
Figure 7 shows a high-level circuit diagram of the algorithm. First we
initialise the registers: we use $n=\lceil\log_{2}(C+1)\rceil$ qubits to store
the number $C$. To extract the period $r$ of $f(x)$, we need one register
initialised as $\ket{0}$ and containing $2n$ qubits to store values $x$, and
another register initialised as $\ket{1}$ with $n$ qubits to store $f(x)$.

Applying the Hadamard gate on the qubits in $\ket{0}$ produces a superposition
state of all integers, such that
$\ket{0}\to\sum_{j=0}^{2^{n}-1}\ket{j}/\sqrt{2^{n}}$. Then the unitary
$\hat{U}$ associates each input $j$ with the value $a^{j}\text{mod}\,C$
(equation 13), for some random $a$ which has no common factors with $C$, to
get

$\frac{1}{\sqrt{2^{n}}}\sum_{j=1}^{2^{n}-1}\ket{j}\ket{a^{j}\text{mod}C}.$
(14)

The first register is in a superposition of $2^{n}$ terms $\ket{j}$ and the
circuit calculates the modular exponentiation for $2^{n}$ values of $j$ in
parallel. Next, inverse QFT on the first register gives

$\frac{1}{2^{n}}\sum_{k=0}^{2^{n}-1}\sum_{j=0}^{2^{n}-1}e^{2\pi
ijk/2^{n}}\ket{k}\ket{a^{j}\text{mod}C},$ (15)

where interference causes only $k$ terms with integer multiples of $2^{n}/r$
to have a substantial amplitude. This is analogous to the quantum phase
estimation in equations 9 and 10. Finally, measuring the first register gives
the period $r$. Ideally, the measurement outcome is an integer multiple of
$2^{n}/r$ with high probability. We can then deduce $r$ classically, with
continued fractions.

Note that Kitaev’s QPE algorithm Kitaev [1996] gives an alternate derivation
of Shor’s algorithm. It reduces the order-finding problem to one involving
phase estimation by essentially replacing the unitary $\hat{U}$ (Figure 7)
with a QPE subroutine.

Classical computers have to date factored semi-primes up to 829 bits long
Zimmermann [2019]. Each extra bit doubles the classical computational cost.
However, quantum computers large enough to beat this are still some way into
the future. This is due to the error correction overheads required to carry
out the quantum computation accurately enough van Meter and Devitt [2016].

####  2.3.4 Grover’s algorithm

Grover’s algorithm Grover [1996, 1997] is a technique for searching
unstructured data. The canonical example is finding the name corresponding to
a given phone number. It can be used as a building block in other
applications, such as database search and constraint satisfaction problems
Ambainis and Špalek [2006].

Figure 8: Grover’s algorithm. (a) Initial state vector $\ket{d}$ with angle
$\vartheta/2$. (b) Vector $\ket{d}$ after inversion and reflection. (c) Final
state vector $\ket{\psi}$.

Suppose there are $N$ items labelled by bit strings $x=0$ to $N-1$, and there
is a unique element $x_{0}$ such that $f(x)=1$ if and only if $x=x_{0}$. The
quantum algorithm uses $n$-qubit basis states $\ket{x}$ corresponding to the
labels of the items, with $n=\lceil\log_{2}N\rceil$.

Grover’s algorithm uses an oracle, which in computer science is a mathematical
device for providing part of the computation that is not included in the
analysis. Consequently, Grover’s algorithm is only of practical use in
situations where we can perform the oracle function efficiently, i.e., $f(x)$
can be evaluated for a superposition of all possible inputs $\ket{x}$. This
can happen when the algorithm is used as a subroutine, but is not generally
possible when the data are classical. The oracle
$\hat{\mathcal{O}}=(-1)^{f(x)}$ has the effect

$\hat{\mathcal{O}}\ket{x}=\begin{cases}-\ket{x_{0}},&x=x_{0}\\\
\ket{x},&\textrm{otherwise.}\\\ \end{cases}$ (16)

Figure 8 illustrates in two dimensions the amplitude amplification part of
Grover’s algorithm, through inversion about the mean. For large $N$, the
initial state

$\ket{d}=\frac{1}{\sqrt{N}}\sum_{x}\ket{x}$ (17)

is an equal superposition of all $\ket{x}$. Note that $\ket{d}$ is almost
orthogonal to the marked state $\ket{x_{0}}$. In Figure 8(a), the small angle
$\vartheta/2\simeq 1/\sqrt{N}$. The oracle $\hat{\mathcal{O}}$ acts on only
$\ket{x_{0}}$ by inverting the sign of its amplitude. Geometrically,
$\hat{\mathcal{O}}$ reflects vector $\ket{d}$ about the horizontal axis. Then
applying the Grover operator

$\hat{U}_{G}=2\ket{d}\bra{d}-\mathbbm{1}$ (18)

reflects $\hat{\mathcal{O}}\ket{d}$ about the axis defined by $\ket{d}$.
Figure 8(b) shows the action of $\hat{\mathcal{O}}$ and $\hat{U}_{G}$. Since
$\vartheta$ is a small angle, this achievement is modest. Solving for
$r\vartheta=\pi/2$ gives an approximate number of required iterations $r$,

$r=\Big{\lfloor}\frac{\pi}{2\vartheta}\Big{\rfloor}=\Big{\lfloor}\frac{\pi}{4}\sqrt{N}\Big{\rfloor}$
(19)

to reach the final state within $\vartheta/2$ of $\ket{x_{0}}$. The
probability of measuring $\ket{x_{0}}$ is better than $1-1/N$.

Since the number of iterations $r$ scales as $\sqrt{N}$, Grover’s algorithm
generally provides a square root speed up over a classical search, which would
need to check on average $N/2$ items to find $x_{0}$. Among the proposed
applications, a Grover search-based method for matched filtering can help
detect gravitational waves Gao et al. [2022] by analysing noisy data collected
by the Laser Interferometer Gravitational Observatory (LIGO). Matched
filtering has detected numerous gravitational waves since LIGO picked up its
first signal in September 2015, but it is time-consuming and resource-
intensive. The quantum method provides a square-root speedup in the
computational cost of searching through a large template bank, assuming the
oracle function can be efficiently implemented.

There are alternative quantum algorithms that solve the search problem, such
as quantum walks Shenvi et al. [2003] (section 2.3.6). These are easily
generalised to find multiple items Childs and Eisenberg [2005] and are
intuitive for graph-based search problems.

####  2.3.5 Quantum amplitude amplification

Grover’s algorithm was subsequently generalised to the framework of quantum
amplitude amplification (QAA) Brassard et al. [2002]. We can understand QAA as
an iterative process that starts with a uniform superposition of all states.
Then the algorithm increases the probability amplitude of some target state
while reducing all other probability amplitudes at each iteration. Amplitude
amplification is a valuable subroutine that appears in many quantum
algorithms. In the Harrow-Hassidim-Lloyd (HHL) algorithm for solving linear
systems (section 2.3.7), the variable time QAA algorithm reduces the number of
repetitions needed to obtain the correct solution which reduces the algorithm
run-time Ambainis [2012]. QAA also appears in applications that involve
optimisation, determining graph connectivity, matching patterns, quantum
counting, and searching for crypto-keys.

####  2.3.6 Quantum walks

In the last few decades, classical random walks and Markov chains have
generated powerful new algorithms in computer science and mathematics Motwani
and Raghavan [1995]. They were thus a natural place to look for quantum
equivalents. Early in the development of quantum algorithms both discrete-time
Aharonov et al. [2001], Ambainis et al. [2001] and continuous-time quantum
walks Farhi and Gutmann [1998] were introduced with algorithmic applications
in mind. The first proven quantum speed ups for quantum walk algorithms came
soon after, with a search algorithm Shenvi et al. [2003] equivalent to
Grover’s algorithm Grover [1996], and an algorithm for transport across a
particular type of disordered graph Childs et al. [2003] that provides an
exponential speed up with respect to an oracle. There are many comprehensive
reviews of quantum walks and their applications, including Venegas-Andraca
[2012] and Kadian et al. [2021].

Quantum walks are widely studied for their mathematical properties, and for
providing simple models of physical phenomena. In an algorithmic context, they
provide quantum speed up for sampling problems (see section 3.2 on quantum
Monte Carlo). The quantum walk dynamics provide both faster spreading Ambainis
et al. [2001] and faster mixing Aharonov et al. [2001], but also localisation.
These are exploited in the search algorithm Shenvi et al. [2003]. The fast
mixing dynamics can also be exploited in a quantum annealing (section 2.4)
context, see, for example Callison et al. [2019]. Quantum walks also provide
models of universal quantum computing Childs [2009], Lovett et al. [2010], but
they are not suitable for physical implementation. Instead these are mostly
useful for complexity proofs.

Quantum walks with multiple non-interacting walkers and quantum particle
statistics include the boson sampling model Aaronson and Arkhipov [2011], a
computational model intermediate between classical and quantum computing.
Quantum walks with multiple interacting walkers Childs et al. [2013] are a
full universal model of quantum computing. They are a special case of quantum
cellular automata (see, for example, Schumacher and Werner [2004]) and can
provide an alternative physical architecture for quantum computing, most
suitable for neutral atoms in optical lattices Karski et al. [2009].

####  2.3.7 Harrow-Hassidim-Lloyd algorithm

In 2009, Harrow, Hassadim and Lloyd (HHL) Harrow et al. [2009] presented their
algorithm to solve one of the most basic classical problems in scientific
computing. The HHL algorithm efficiently calculates the solution
$\vec{x}=\mathbf{A}^{-1}\vec{b}$ of the linear system
$\mathbf{A}\vec{x}=\vec{b}$ for a sparse, regular $N\times N$ Hermitian matrix
$\mathbf{A}$. The hard part of this problem is to invert $\mathbf{A}$ and
obtain $\mathbf{A}^{-1}$. To see how a quantum method might help, we present a
high level description of the HHL algorithm. For further details and analysis,
see Dervovic et al. [2018], Abhijith et al. [2022].

Figure 9: High level HHL circuit. Adapted from Dervovic et al. [2018].

We start with the ancilla, clock, and input registers (Figure 9). The HHL
algorithm represents $\vec{b}$ as a quantum state in the computational basis

$\ket{\vec{b}}=\sum_{j}b_{j}\ket{j}.$ (20)

and stores it in the input register with the amplitude encoding method
(section 5.2). Then using Hamiltonian simulation, we apply

$\hat{U}=e^{i\hat{A}t}$ (21)

to $\ket{b}$ for a superposition of times $t$ determined by the clock
register. This step corresponds to quantum phase estimation (QPE) which
expresses $\ket{b}$ in terms of the eigenbasis $\ket{u_{j}}$ of $\mathbf{A}$
and finds the corresponding eigenvalues $\lambda_{j}$ which are stored in the
clock register. Now we can express the system state as

$\sum_{j}\beta_{j}\ket{u_{j}}\ket{\lambda_{j}}$ (22)

with $\beta_{j}=\langle u_{j}|\vec{b}\rangle$ and

$\ket{\vec{b}}=\sum_{j}\beta_{j}\ket{u_{j}}.$ (23)

Note that we do not know what the values of $\ket{u_{j}}$.

In the next step, we perform a controlled rotation on the ancilla register
with clock register as the control. The procedure extracts the eigenvalues of
$\mathbf{A}^{-1}$ by performing a $\sigma_{y}$ rotation, conditioned on
$\lambda_{j}$. It is equivalent to the linear map,
$\ket{\lambda_{j}}\to\lambda_{j}^{-1}\ket{\lambda_{j}}$. This transforms the
system to

$\sum_{j}\beta_{j}\ket{u_{j}}\ket{\lambda_{j}}\left(\sqrt{1-\frac{c^{2}}{\lambda_{j}^{2}}\ket{0}+\frac{c}{\lambda_{j}}\ket{1}}\right)$
(24)

with normalisation constant $c$.

The final step uses inverse QPE to uncompute the $\ket{\lambda_{j}}$ register.
This leaves the remaining state as

$\sum_{j}\beta_{j}\ket{u_{j}}\ket{\lambda}\left(\sqrt{1-\frac{c^{2}}{\lambda_{j}^{2}}}\ket{0}+\frac{c}{\lambda_{j}}\ket{1}\right).$
(25)

Measuring and post-selecting the ancilla qubit of $\ket{1}$ outputs the state

$\ket{x}=A^{-1}\ket{b}\approx
c\sum_{j}\frac{\beta_{j}}{\lambda_{j}}\ket{u_{j}}.$ (26)

Note that if the ancilla measurement outputs $\ket{0}$, we rerun the
algorithm.

There are many limitations of the HHL algorithm Aaronson [2015], Dervovic et
al. [2018]. In short, they are:

  * •

$\mathbf{A}$ must be Hermitian.

  * •

Eigenvalues of $\mathbf{A}$ must be in $[0,1)$.

  * •

We must be able to efficiently implement QPE with $e^{i\hat{A}t}$.

  * •

We must be able to efficiently prepare $\ket{b}$.

  * •

We cannot recover the full solution $\ket{x}$. Instead we may recover the
expectation value of some operator $\hat{M}$, such as $\bra{x}\hat{M}\ket{x}$.

The importance of HHL for practical applications has led to many improvements
on the original algorithm, such as using time variable amplitude amplification
to increase the success probability during ancilla measurement Ambainis
[2012]. This is a more general case of amplitude amplification (section
2.3.5), better suited to quantum algorithms where individual branches finish
at different times.

While the overall algorithm has seen improvements, so have its individual
components such as Hamiltonian simulation in the QPE step (equation 21). QPE
remains a bottleneck by requiring $\Omega(1/\epsilon)$ applications of a
unitary operation to estimate its eigenvalues to precision $\epsilon$. Hence,
even if we use the fastest Hamiltonian simulation methods, the overall
improvement would likely be modest with total complexity still being
polynomial in $(1/\epsilon)$. One remedy involves directly applying the matrix
inverse to bypass the phase estimation subroutine Childs et al. [2017], thus
exponentially improving dependence on the precision parameter. This approach
has subsequently found good use in applications ranging from accelerating
finite element methods (FEMs) Montanaro and Pallister [2016] to approximating
the hitting times of Markov chains Chowdhury and Somma [2017].

Alternatively, HHL can be “hybridised” into a quantum-classical method. Given
the resource limitations of NISQ-era hardware, hybrid approaches will likely
pave the way for meaningful applications in the near future. For example, an
iterative-HHL approach Saito et al. [2021] contains a classical iterative
process to improve accuracy beyond the limit imposed by number of qubits in
QPE. However, the convergence rate for more general cases remains uncertain.
Other approaches include streamlining the algorithm, such as using a classical
information feed forward step after the initial QPE to reduce the circuit
depth for subsequent steps in HHL Lee et al. [2019a]. Of course, if the
reduced circuit is not applicable, it would then require repeating QPE to
determine an appropriate reduction. A hybrid approach can achieve comparable
precision to conventional HHL with fewer qubits and multiple phase estimation
modules Gao et al. [2023a]. It can also improve phase estimation itself to
mitigate the long coherence times Angara et al. [2020].

Solving dense matrices, particularly relevant in machine learning and kernel
methods Wilson et al. [2015], has been tackled by using quantum singular value
estimation (QSVE) to gain a polynomial advantage over traditional HHL
approaches Wossnig et al. [2018], Bravo-Prieto et al. [2020]. There are many
other QML algorithms where HHL is an important subroutine, such as data
classification using quantum support vector machines Rebentrost et al. [2014],
quantum principal component analysis Lloyd et al. [2014], or quantum linear
discriminant analysis Cong and Duan [2016]; quantum ordinary linear regression
(QOLR) Wiebe et al. [2012], Wang [2017] and QOLR for prediction Schuld et al.
[2016]; quantum ridge regression (QRR) Yu et al. [2021]; quantum
recommendation systems Kerenidis and Prakash [2017]; quantum singular value
thresholding Duan et al. [2018]; and quantum Hopfield neural networks
Rebentrost et al. [2018]. More advanced HHL-based methods include solving
linear differential systems Berry [2014], Berry et al. [2017]; using
techniques from adiabatic optimisation (section 2.4) to solve linear equations
Subaşı et al. [2019]; and reformulating least-squares regression as a QUBO
problem and solving with adiabatic quantum computation Date and Potok [2021].

###  2.4 Quantum optimisation

Figure 10: The full quantum annealing process from problem formulation to
solution. The graphs illustrate how the initially uniform distribution evolves
under the driver Hamiltonian to concentrate in the global minimum representing
the solution. Adapted from Au-Yeung et al. [2023].

Combinatorial optimisation problems are widely encountered across industry
Yarkoni et al. [2022], academia Hauke et al. [2020], and medical science Emani
et al. [2021]. Typically, the problem is to find the best of an exponentially
large set of solutions. Mathematically, a basic optimisation problem consists
of the objective (cost) function which is the output we want to maximise or
minimise, the input variables, and optional constraints.

Many optimisation problems can be expressed as quadratic unconstrained binary
optimisation (QUBO) problems. These are encoded using an upper-triangular
$N\times N$ real matrix $\mathbf{Q}$ and vector of binary variables $\vec{x}$.
The aim is to minimise the cost function

$f(\mathbf{x})=\sum_{j}\mathbf{Q}_{jj}\vec{x}_{j}+\sum_{j<k}\mathbf{Q}_{jk}\vec{x}_{j}\vec{x}_{k}.$
(27)

Many cost functions contain a large number of false local minima, making it
difficult for classical algorithms to find the true global minimum.

Classical simulated annealing Kirkpatrick et al. [1983] is based on the idea
that thermal fluctuations move the system out of local minima and toward the
lowest potential energy state as temperature decreases. The simulated
annealing algorithm simulates a random walker that travels through the search
space or optimisation landscape. The rate at which the temperature is reduced
determines how fast the system evolves and how likely it is to avoid becoming
stuck in local minima. Quantum annealing Finilla et al. [1994], Kadowaki and
Nishimori [1998] was conceived as an alternative that adds quantum tunnelling
to the system evolution, allowing it to escape from local minima more easily.
Figure 10 illustrates the quantum annealing process.

Conveniently, classical optimisation problems can be mapped to finding the
ground state of a classical Ising Hamiltonian Lucas [2014], Choi [2010]. This
leads naturally to a quantum version using the transverse field Ising
Hamiltonian,

$\displaystyle\hat{H}_{\text{Ising}}=A(t)\underbrace{\sum_{j}\hat{\sigma}^{x}_{j}}_{\hat{H}_{0}}+B(t)\underbrace{\sum_{j,k}\left(h_{j}\hat{\sigma}^{z}_{j}+J_{jk}\hat{\sigma}^{z}_{j}\hat{\sigma}^{z}_{k}\right)}_{\hat{H}_{p}}$
(28)

with are the Pauli $x$\- and $z$-matrices $\hat{\sigma}^{x,z}$, symmetric
interaction strength $J_{jk}=J_{kj}$ of qubit spins $q_{j}$ and $q_{k}$, and
on-site energy $h_{j}$. Note that the qubit spin basis states are $\ket{\pm
1}$. There is a straightforward translation from QUBO $\vec{x}$ to Ising
$\vec{q}$, using $\vec{x}_{j}=(\vec{q}_{j}+1)/2$. The resulting values of the
$\\{h_{j},J_{jk}\\}$ variables then encode the problem into the Ising
Hamiltonian $\hat{H}_{p}$.

The transverse field in the $x$-direction $\hat{H}_{0}$ provides the dynamics
that rotate the qubits from their initial state to the final target state,
ideally the ground state of problem Hamiltonian $\hat{H}_{p}$. The control
functions $A(t)$ and $B(t)$ must be specified. Usually $A(t)$ varies from 1 to
0, and $B(t)$ from 0 to 1. The choice of how they vary determines the method
used to find the ground state and hence the problem solution.

Adiabatic quantum computing (AQC), introduced by Farhi et al. Farhi et al.
[2001] in 2001, uses the quantum adiabatic theorem to guarantee that the
system remains in the ground state. It assumes that the control functions
$A(t)$ and $B(t)$ are varied slowly enough to keep the system in the
instantaneous ground state. AQC has several useful properties: it is
inherently robust against noise Childs et al. [2001]. It is computationally
equivalent to gate-based quantum computers Aharonov et al. [2004], making it a
universal quantum computing paradigm Aharonov et al. [2008].

While AQC provides a sound theoretical underpinning for quantum annealing, in
reality, the adiabatic conditions imposed by AQC are rarely met. Meeting them
would require long run times that are inefficient and impractical. At the
other extreme, continuous-time quantum walk algorithms Farhi and Gutmann
[1998] can efficiently locate ground states for certain problems Callison et
al. [2019] where it is viable to make many short repeats (section 2.3.6).
Further diabatic methods are reviewed in Crosson and Lidar [2021]. Many of
these, including the quantum approximate optimisation algorithm (QAOA) or
quantum alternating operator ansatz Farhi et al. [2014], Hadfield et al.
[2019] use classical optimisation in the controls to produce an efficient
quantum process. We describe QAOA in subsection 2.5.1 after introducing
variational quantum algorithms more generally. Exploiting the non-adiabatic
and open system effects in quantum annealing is the topic of many recent
efforts to derive efficient controls Banks et al. [2023], Schulz et al.
[2023]. For more details, Albash and Lidar Albash and Lidar [2018] provide a
comprehensive introduction to adiabatic quantum computing and quantum
annealing. Hardware currently available for quantum annealing includes
superconducting systems, trapped ions, and Rydberg atoms (see Hauke et al.
[2020] for an overview).

###  2.5 Variational quantum algorithms

Figure 11: High level diagram of hybrid quantum-classical training algorithm
for variational circuit. Quantum device (yellow) calculates terms of a cost
function, then classical device calculates better circuit parameters
$\vec{\theta}$ (blue, green). Repeat cycle until desired accuracy is achieved.

Variational quantum algorithms (VQAs) distribute the job of solving a problem
between a parameterised quantum circuit and conventional classical optimiser
Cerezo et al. [2021]. Figure 11 shows schematically how it works. The quantum
circuit is defined by a set of parameters $\vec{\theta}$. We evaluate a
classical cost function $C(\vec{\theta})$ with the measured outputs
$\vec{\theta}$ of the quantum circuit. Classical optimisation then updates the
parameters $\vec{\theta}$, guided by the cost function. This general framework
can solve a wide range of problems, including quantum chemistry and materials
science Kokail et al. [2019], Kandala et al. [2017] (further details in
section 3).

Most VQAs can be grouped into one of two categories: variational quantum
optimisation and variational quantum simulation. The optimisation algorithm
involves optimising parameters under a static target cost function. For
example, when the cost function corresponds to the expectation value of some
Hamiltonian, minimising the cost function gives the ground state energy of the
Hamiltonian. If the Hamiltonian encodes an optimisation problem, as in section
2.4, the state that minimises the cost function yields the solution to the
original problem. If the Hamiltonian corresponds to a quantum system, such as
a molecule, the cost function gives the molecule’s ground state energy. On the
other hand, variational quantum simulation aims to simulate dynamical
processes such as the Schrödinger time evolution of a quantum state Yuan et
al. [2019].

The past years have seen rapid developments in the variational approach, and
it is clear there are beneficial properties for NISQ-era devices Cerezo et al.
[2021]. It is possible to use many separate repeats of a shallow (short)
quantum circuit. This avoids the need for long coherence times, and the
variational parameter tuning compensates for several types of errors. Suitable
error mitigation strategies have been developed for NISQ devices Botelho et
al. [2022], Cai et al. [2023]. VQAs can be applied to a wide range of linear
optimisation problems and extended to solve nonlinear variational problems
Lubasch et al. [2020].

Variational algorithms may potentially demonstrate a quantum advantage on NISQ
devices. However, they are highly demanding on quantum hardware, and require
large numbers of sequential calls to quantum devices. The low-throughput of
quantum hardware is a severe limitation. It is a major barrier towards
studying larger and more complex problems in industry.

There are two important algorithms that make use of the variational principle
for optimisation: the quantum approximate optimisation algorithm (QAOA) and
the variational quantum eigensolver (VQE).

####  2.5.1 Quantum approximate optimisation algorithm

The quantum approximate optimisation algorithm (QAOA) was originally
introduced to approximately solve combinatorial optimisation problems Farhi et
al. [2014] on gate-based architecture. This was generalised as the quantum
alternating operator ansatz Hadfield et al. [2019]. Instead of applying the
full Hamiltonian $\hat{H}=\hat{H}_{M}+\hat{H}_{P}$ continuously in time, the
aim is to map some input state to the ground state of a problem Hamiltonian
$\hat{H}_{P}$ by sequentially applying the ansatz

$\hat{U}(\gamma,\beta)=\prod_{k}e^{-i\beta_{k}\hat{H}_{M}}e^{-i\gamma_{k}\hat{H}_{P}},$
(29)

with mixer Hamiltonian $\hat{H}_{M}$. The ansatz is applied in short discrete
steps, that approximate the continuous time evolution, with parameters
$\vec{\theta}=(\vec{\gamma},\vec{\beta})$ to be optimised. In-depth studies of
QAOA performance Zhou et al. [2020a], Wang et al. [2020] developed an
efficient parameter-optimisation procedure to find high-quality parameters. It
has been shown to converge to quantum annealing in the transverse Ising
setting Brady et al. [2021].

####  2.5.2 Variational quantum eigensolvers

Variational quantum eigensolvers (VQEs) were originally developed to find the
ground state energy of molecules McClean et al. [2016], Kandala et al. [2017].
They have since been adapted to solve optimisation problems Cerezo et al.
[2021], Amaro et al. [2022]. VQEs are similar to QAOA, but analyse the
Hamiltonian differently. The cost function is defined as the expectation value
of the quantum state energy, $\bra{\psi(\theta)}\hat{H}\ket{\psi(\theta)}$.
The aim is to minimise the expectation value of Hamiltonian $\hat{H}$ over a
trial state $\ket{\psi(\theta)}=\hat{U}(\theta)\ket{\psi_{0}}$ for ansatz
$\hat{U}(\vec{\theta})$ and initial state $\ket{\psi_{0}}$.

The Rayleigh-Ritz principle states that the lowest eigenvalue of a Hermitian
operator $\hat{H}$ is upper-bounded by the minimum expectation value found by
varying the state $\ket{\psi}$. In other words, the statement

$\bra{\psi(\theta)}\hat{H}\ket{\psi(\theta)}\geq E_{0}$ (30)

is valid with the equality holding if $\ket{\psi(\theta)}$ is the ground state
$\ket{E_{0}}$ of $\hat{H}$. Hence, the optimised vector $\vec{\theta}$ that
approximates the ground state (or eigenvector corresponding to the lowest
eigenvalue) is the choice that minimises $\braket{\hat{H}}$. Using this
principle, we can break down the VQE methodology into three distinct parts:

  1. 1.

Prepare a parameterised initial state as the starting ansatz
$\ket{\psi_{0}(\vec{\theta})}$.

  2. 2.

Estimate expectation value $\braket{\hat{H}}(\vec{\theta})$ in the quantum
processor (this may require multiple runs using quantum phase estimation).

  3. 3.

Use a classical optimiser to find a new set of $\vec{\theta}$ values that
decreases $\braket{\hat{H}}(\vec{\theta})$.

  4. 4.

Repeat above procedure to achieve convergence in $\langle\hat{H}\rangle$.
Parameters $\vec{\theta}$ at convergence define the desired state.

####  2.5.3 Variational quantum simulation

Simulating the dynamics of quantum systems is expected to be one of the first
useful applications of quantum computers. Variational approaches make best use
of limited quantum hardware by optimising the parameters classically. There
are two basic algorithms for variational quantum simulation, specifically for
solving real- and imaginary-time evolution of quantum systems. These can be
described respectively by the time-dependent Schrödinger equation

$\frac{d\ket{\psi(t)}}{dt}=-i\hat{H}\ket{\psi(t)}$ (31)

with Hamiltonian $\hat{H}$ and state $\ket{\psi(t)}$, and normalised Wick-
rotated Schrödinger equation

$\frac{d\ket{\psi(\tau)}}{d\tau}=-(\hat{H}-\braket{\hat{H}})\ket{\psi(\tau)},$
(32)

where $\tau=-it$ and
$\braket{\hat{H}}=\bra{\psi(\tau)}\hat{H}\ket{\psi(\tau)}$ preserve the norm
of state $\ket{\psi(\tau)}$.

The direct approach to quantum simulation is to apply a unitary circuit
$e^{-i\hat{H}t}$ to the initial state. However this increases the circuit
depth polynomially with respect to the evolution time $t$ Low and Chuang
[2019] which is difficult to achieve in hardware with short coherence times.
On the other hand, variational quantum simulations use an ansatz quantum
circuit

$\ket{\phi(\vec{\theta}(t))}=\hat{U}(\vec{\theta}(t))\ket{\phi_{0}}$ (33)

applied to an initial state $\ket{\phi_{0}}$ to represent $\ket{\psi(t)}$.
Therefore, we must choose the variational parameters to map the time evolution
of the Schrödinger equation for $\ket{\psi(t)}$ on to the evolution of
parameters $\vec{\theta}(t)$. The aim is to obtain a shorter circuit than
would be needed for direct simulation. How the parameters evolve and simulate
the time evolution depends on the variational principle used Yuan et al.
[2019].

The variational method for generalised time evolution can perform matrix
multiplication and solving linear systems of equations Endo et al. [2020]. It
is an alternative method to those presented in Bravo-Prieto et al. [2020], Xu
et al. [2021], or the HHL algorithm (section 2.3.7).

####  2.5.4 Quantum neural networks and parameterised quantum circuits

Classical neural networks (NNs) LeCun et al. [2015] allow us to classify and
cluster large datasets by training the parameters associated with the neural
connections. This is analogous to the learning process in the brain. Quantum
neural networks (QNNs) combine artificial neural networks with quantum
principles Abbas et al. [2021]. NNs are highly nonlinear models, whereas QNNs
obey the laws of quantum mechanics and contain linear operators acting on
quantum states. QNNs can offload nonlinearities to classical computers, or use
quantum kernels. The latter option encodes classical data into Hilbert space
using a nonlinear transformation, the quantum feature map Schuld and Killoran
[2019].

There are many different proposals for QNN architectures, such as the fully
quantum, feed-forward neural networks Farhi and Neven [2018], and
convolutional neural networks Cong et al. [2019]. Most QNNs follow similar
steps Beer et al. [2020]: initialise a neural network architecture; specify a
learning task; implement a training algorithm; and simulate the learning task.

Figure 12: Example of quantum neural network using variational quantum
circuit.

Quantum neural networks are a subclass of variational quantum algorithms Abbas
et al. [2021]. QNNs are quantum circuits that contain parameterised gate
operations Schuld et al. [2020] (Figure 12). They offer a concrete way to
implement quantum reinforcement learning algorithms in the NISQ era. The
approach requires redefining problems as variational optimisation problems
then using hybridised quantum-classical devices to find approximate solutions.

There are typically three steps. Given a dataset of pairs $(x,y)$, we want to
learn a model function $f$ that maps each $x$ to its corresponding target $y$.
We encode the input data (usually classical) into a quantum state using a
quantum feature map Schuld et al. [2021] $\phi:V\to V^{\prime}$, from input
space $V$ to enhanced feature space $V^{\prime}$. The choice of feature map
aims to enhance the QNN performance. Then the kernel is an inner product
$K(\vec{x}_{j},\vec{x}_{k})\equiv\braket{\phi(\vec{x}_{j}),\phi(\vec{x}_{k})}$.
This corresponds to a state preparation circuit $\hat{U}_{\phi}(x)$ that maps
classical datapoint $x$ to an $n$-qubit quantum state $\ket{\phi(x)}$,

$\hat{U}_{\phi}(x)\ket{0}^{\otimes n}=\ket{\phi(x)}.$ (34)

Next, we use a variational model to build the quantum state

$\ket{\psi(x,\vec{\theta})}=V(\vec{\theta})\ket{\phi(x)}$ (35)

which is parameterised by vector $\vec{\theta}$. It contains parameterised
gate operations that are optimised for a particular task, analogous to
classical ML techniques Dunjko and Briegel [2018]. The aim is to choose
parameters $\vec{\theta}$ to minimise a certain loss function. Lastly, we
measure an observable $\hat{M}$. Classical post-processing procedures on the
measurement outcomes $\vec{z}=[z_{1},z_{2},...,z_{n}]$ extracts the output of
model

$f(\vec{z})=\bra{\psi(x,\vec{\theta})}\hat{M}\ket{\psi(x,\vec{\theta})}.$ (36)

While the approach is straightforward, the state space of a polynomial-sized
quantum circuit is much smaller than the full Hilbert space. Hence, the main
advantage of variational circuits is the feature map: when $\hat{U}_{\phi}(x)$
cannot be simulated efficiently on classical computers, the overall quantum
speedup can be up to exponential.

In quantum chemistry, PQCs such as the unitary coupled-cluster (UCC) circuits
can produce states that approximate the molecular ground states Lee et al.
[2019b]. UCCs are based on the coupled-cluster method in computational
chemistry. They are too deep to be easily implemented on today’s quantum
computers. Instead we use the hardware-efficient ansatz, a shallower PQC
composed of gates that are natural operations on the quantum hardware Kandala
et al. [2017].

The dissipative QNN is analogous to the classical feedforward neural network
Schuld et al. [2014], Beer et al. [2020], a widely-used framework in classical
ML. Many other QNNs have been proposed for different types of data analysis,
such as quantum Boltzmann machines Amin et al. [2018], Kieferová and Wiebe
[2017] which are a natural fit for quantum annealing Benedetti et al. [2016],
Xu and Oates [2021] applications. Constructing QNN architectures is currently
an active area of research. For a concise review, see Cerezo et al. [2022] and
references within.

####  2.5.5 Limitations of VQAs

Variational algorithms have a range of problems (see Cerezo et al. [2021,
2022] for an overview). Here we note two common problems.

Figure 13: Ideal solution landscape (blue) corrupted by (a) noise and (b)
barren plateau, both orange.

Hardware noise is a defining characteristic of NISQ computing. Noise affects
all aspects of a quantum model, corrupting information as it forward-
propagates in a quantum circuit. It particularly affects deeper circuits with
longer runtimes. Noise can hinder the trainability of hybrid algorithms as it
leads to noise-induced barren plateaus Wang et al. [2021a], ultimately
deforming the landscape features and shifting the minima positions (Figure
13(a)). Suggested antidotes include partial quantum error correction Bultrini
et al. [2023], building relatively shallow circuits (depth grows sublinearly
in the problem size) Wang et al. [2021a], and error mitigation techniques to
partially alleviate noise-induced trainability issues Wang et al. [2021b],
Endo et al. [2021]. Alternatively, we can create models with noise-resilient
properties Fontana et al. [2021] that do not alter the minima positions.

Quantum machine learning and VQAs contain expensive classical optimisation
loops. They minimise a cost function by navigating a typically non-convex loss
function landscape to find its global minimum. The barren plateau phenomenon
occurs when the optimisation landscape becomes exponentially flat (on average)
as qubit number increases McClean et al. [2018], Cerezo et al. [2021]. At the
same time, the valley containing the global minimum also shrinks exponentially
with problem size, leading to a narrow gorge Arrasmith et al. [2022] (Figure
13(b)). The cost function gradient is exponentially suppressed except in an
exponentially small region. Hence finding the minimum requires computational
resources that scale exponentially with the number of qubits used, which can
potentially destroy any quantum speedup.

###  2.6 Performance benchmarks: quality, speed, scale

Benchmarking for high performance computers is well established, and essential
for managing HPC facilities and optimising their throughput. For example, the
High Performance Linpack (HPL) benchmark Petitet et al. [2018] solves a dense
linear system to measure the throughput of a computing system. Running the
same computational tasks on different HPC hardware can provide some level of
side-by-side comparison of different facilities. However, the results from
standard benchmarks cannot fully predict performance for a specific
application, and bespoke tests are often run to tune performance on a per-
application basis.

Benchmarks for quantum computers are under active development. Given the
diverse types of qubits under development, fair comparisons of quantum
hardware are even more challenging than for classical hardware. There are five
common physical qubit realisations for building quantum processors. Hardware
based on superconducting qubits Krantz et al. [2019], Blais et al. [2021] is
arguably the most popular choice Bravyi et al. [2022], Blais et al. [2020],
followed by ion-trap-based systems Bruzewicz et al. [2019], Malinowski et al.
[2023], Moses et al. [2023]. Semiconductor-based spin qubits Chatterjee et al.
[2021], Burkard et al. [2023] are less advanced but potentially offer more
scalability in future. Photonic platforms Pelucchi et al. [2022], Moody et al.
[2022] require less stringent cooling regimes but other aspects of the
engineering (e.g., single photon sources and photon loss) are challenging.
Neutral atoms, especially using highly excited Rydberg states Adams et al.
[2020], are versatile and already showing impressive results for directly
simulating quantum many-body systems Surace et al. [2020], Ebadi et al.
[2021].

The criteria for good quantum computer hardware are fairly straightforward
DiVincenzo [2000]. It must have as many high quality qubits as possible.
Qubits must be individually controlled to generate complex (entangled) states.
We must be able to apply a large number of sequential operations before the
qubits lose coherence. The control and gate operations are characterised by
the quality of single- and two-qubit gates. Scaling up to obtain good
performance on large numbers of qubits requires error correction Terhal [2015]
to be implemented in the hardware. This needs mid-circuit measurements to be
possible, as well as measuring the qubit states at the end of the computation.

####  2.6.1 Low- and high-level benchmarks

Low-level benchmarks measure the performance of individual qubits, single- and
two-qubit gate fidelities, state preparation and measurement operations. We
can often directly measure the $T_{1,2}$ decoherence times, gate speed or gate
operation reliability (gate fidelity). However these characteristics do not
provide an accurate benchmark of quantum processor performance. Using low-
level benchmarks alone can label a slow high-fidelity quantum computer as
equivalent to a fast low-fidelity quantum computer Campbell and Deffner
[2017].

For more thorough testing of gate operations, there is randomised benchmarking
Knill et al. [2008], Magesan et al. [2011], cycle benchmarking Erhard et al.
[2019], or gate set tomography Nielsen et al. [2021]. Randomised benchmarking
protocols typically characterise the quality of Clifford operations which are
classically simulable. For non-Clifford operations, random circuit sampling
and direct randomised benchmarking can be used Liu et al. [2022], Proctor et
al. [2019]. These methods may not be scalable to large numbers of qubits as it
is difficult to classically simulate non-Clifford operations to check the test
outcomes. A class of methods that aims to circumvent this issue is mirror
benchmarking Proctor et al. [2022a], Mayer et al. [2021]. The trivial ideal
output of the circuits involved in these benchmarks makes the results easily
verifiable and scalable to large numbers of qubits.

At a slightly higher level, specific quantum circuits can be used to test
hardware performance for gate sequences. The circuits can be random,
application-agnostic circuits Cross et al. [2019], Proctor et al. [2022b], or
circuits run as part of certain algorithms Mills et al. [2021], Quetschlich et
al. [2023], Li et al. [2022b]. Other types of gate sequences can be used
including identity operations and entanglement generation circuits Michielsen
et al. [2017]. For non-gate-based hardware, a generalised form of randomised
benchmarking can be used, for example in programmable analogue quantum
simulators Derbyshire et al. [2020].

For benchmarking at the algorithmic level, qubit layout and connectivity
affect algorithm efficiency. Some types of hardware (e.g., neutral atoms Adams
et al. [2020]) can easily be reconfigured between runs. Other platforms, such
as superconducting or semiconductor qubits have connectivity baked in at the
manufacturing stage.

Algorithm- Georgopoulos et al. [2021] or application-level Lubinski et al.
[2023], Donkers et al. [2022] benchmarks require knowing the ideal output in
order to evaluate the solution quality. For running entire algorithms, there
are higher-level benchmarks for how well a quantum computer performs. For
example, the one-dimensional Fermi-Hubbard model is exactly solvable using the
Bethe ansatz Dallaire-Demers et al. [2020]. There are benchmarks comprising a
suite of many different applications (e.g., amd J. Bausch and Gilyén [2021],
Tomesh et al. [2022], Murali et al. [2019]). Most include well-known
algorithms such as the quantum Fourier transform (section 2.3.1), quantum
phase estimation (section 2.3.2), Grover’s search algorithm (section 2.3.4),
VQE (section 2.5.2), the Bernstein-Vazirani algorithm Bernstein and Vazirani
[1993], and the hidden shift algorithm van Dam et al. [2006]. In general,
device performance on simple problems or circuits may not scale Lubinski et
al. [2023]. It is also necessary to define suitable high-level success
criteria, e.g., for quality Blume-Kohout and Young [2020], Mills et al.
[2021].

It is instructive to test at both low and high levels because extrapolating
from low level performance to higher level is unreliable. Cross-talk errors
can dominate and these are not captured by lower level tests using only one-
or two- qubit gates. Even at a higher level, there may be large differences
between performance on random and structured circuits Lubinski et al. [2023],
Proctor et al. [2022b].

Determining whether an actual computation has run correctly requires
verification. Classical verification methods are becoming increasingly
sophisticated for larger, more complex problems (see, for example, Groen et
al. [2022]). Verifying quantum computations Carrasco et al. [2021], especially
for NISQ hardware, requires some extra tools, such the accreditation scheme of
Ferracin et al. Ferracin et al. [2019] for gate-based computations, which has
been experimentally implemented Ferracin et al. [2021]. Accreditation makes
use of similar tools to low-level benchmarking, like mirror circuits, to
quantify the errors and hence assign a confidence level to the correctness of
the output of the computation, which has been run alongside the test circuits.
Accreditation has been extended beyond the circuit model to certain types of
quantum simulation and quantum annealing Jackson et al. [2023a].

####  2.6.2 Quantifying quantum device performance

Various metrics have been proposed to capture the performance of quantum
computers in a single number or a few parameters. Most of these are fairly low
level benchmarks. For example, quantum volume is a benchmark developed by IBM
Cross et al. [2019] to test hardware performance for gate sequences. Quantum
volume is calculated by running randomised square circuits with different
numbers of qubits. The quantum volume is given by the largest square circuit
that can be run on the device which passes an acceptance criterion that
quantifies the output quality. It has been generalised to more realistic non-
square circuits Blume-Kohout and Young [2020], Miller et al. [2022]. Metrics
such as Circuit Layer Operations Per Second (CLOPS) Wack et al. [2021], layer
fidelity McKay et al. [2023], and qubit number have been proposed to
collectively quantify the quality, speed, and scale of quantum computers.
Single-number metrics are useful to make cross-device comparisons
straightforward but similar to classical metrics, they may not accurately
describe quantum computer performance for all kinds of algorithms. McGeoch
McGeoch [2019] proposes guidelines that go beyond the circuit model for
reporting and analysing quantum annealer performance.

It is also important to measure how quickly the quantum processor can execute
layers of a parameterised model circuit similar to those used to measure
quantum volume. Increased quantum processor speed is critical to support near-
term variational algorithms, which require thousands of iterations. Some
benchmarks evaluate speed as a sub-score and combine it with other factors to
determine a final metric Donkers et al. [2022]. There also exist standalone
metrics for speed Wack et al. [2021]. Lubinski et al. Lubinski et al. [2023]
note that evaluating measures of speed for comparing different platforms is
difficult. Timing information (when available at all) is reported differently
by different hardware vendors, and gate speeds naturally vary by several
orders of magnitude between hardware types (see section 5). Another time-
related factor is consistency, or how the performance varies over time. This
determines how often hardware must be recalibrated. Consistency can be
captured by repeating the same tests, and is an important consideration for
current hardware which requires significant downtime for regular
recalibration.

The above methods (section 2.6.1) and metrics (section 2.6.2) are useful for
benchmarking hardware dependent characteristics. On top of this, the compilers
and software stack efficiency can significantly affect performance, in much
the same way as classical HPC requires highly tuned software to optimise
performance. Quantum software benchmarking is still at an early stage of
development.

####  2.6.3 Quantum circuit simulation methods

In the current NISQ era, simulating quantum circuits is important for multiple
purposes. Simulations of ideal quantum circuits allow quantum algorithm
development and testing in a more reliable environment than quantum hardware
currently provides, albeit limited to about 45 qubits. Several benchmarking
tools require knowledge of the ideal circuit outputs produced by such
simulations, to verify the outputs of the quantum hardare. It also establishes
a classical computational bar that quantum computers must pass to demonstrate
a quantum advantage. One of the main tests used to claim a quantum advantage
from a NISQ device involves sampling the bit-strings from a random quantum
circuit (RQC) Boixo et al. [2018], Villalonga et al. [2019]. The result must
be within some variational distance of the output distribution defined by the
circuit Bouland et al. [2019].

It is important to carry out classical simulations as efficiently as possible,
so that reports of quantum advantage are not exaggerated. The tensor network
contraction method is one of the best classical methods for simulating RQCs
for sizes close to the quantum advantage regime Wahl and Strelchuk [2023],
Huang et al. [2021] Approximating a quantum state using tensor networks is
equivalent to applying lossy data compression that preserves the most
important properties of the quantum state (see, e.g., Orús [2014], Huggins et
al. [2019]). They are thus especially effective for simulating low-depth noisy
quantum circuits Zhou et al. [2020b], Orús [2019], Ran et al. [2021]. Tensor
network states mathematically represent quantum many-body states based on
their entanglement structure. Because RQCs (like in Google’s work Arute et al.
[2019]) contain noisy quantum gates that create a small amount of
entanglement, applying tensor-network algorithms is adequate for their
simulation Zhou et al. [2020b]. Contracting a whole tensor network involves a
contraction tree which describes a series of pairwise merges that turn the
network into a single tensor. It is analogous to performing a series of sums
and multiplications of the tensor entries to obtain an expectation value.
Examples of tensor network-based benchmarking methods include the holoQUADS
algorithm Foss-Feig et al. [2021, 2022] and qFlex Villalonga et al. [2019].

####  2.6.4 Measuring quantum advantage

It is not straightforward to determine exactly when a quantum computer
provides a practical advantage over a classical computer. The latter are
benchmarked using FLOPS, which measure the number of floating-point operations
per second. The equivalent speed benchmark for quantum processors is CLOPS
which measures the number of circuit layer operations per second. It does not
make sense to convert between the two metrics, since they do not correspond to
the same unit of computation.

Demonstrating a quantum advantage at the level of whole algorithms also
highlights the difference between theoretical computational complexity and
practical implementation. Many quantum machine learning (QML) algorithms make
different assumptions from the classical algorithms they are compared with
Aaronson [2015]. In some cases, it is possible to find “dequantised” versions
of QML algorithms Tang [2022]. These are fully classical algorithms that
process classical data and perform only polynomially slower than their quantum
counterparts. If a dequantised algorithm exists, then its quantum counterpart
cannot give an exponential speedup on classical data. On the other hand,
dequantisation shows how replacing quantum linear algebra algorithms with
classical sampling techniques can potentially create a classical algorithm
that runs exponentially faster than any other known classical algorithm Tang
[2019]. Quantum-inspired classical algorithms that provide a significant
improvement on current best methods are a useful benefit from quantum
algorithm research, if they turn out to be practical to implement on classical
HPC.

Many QML algorithms rely on pre-processed data via quantum encoding.
Breakthroughs that claim a quantum speed-up often neglect the quantum encoding
step in the algorithm runtime, while making comparisons to classical
algorithms that operate on raw, unprocessed data. This means it is often
possible to achieve equivalent speed-ups using only classical resources, if we
perform the same pre-processing for classical algorithms. For example, the
pre-processing in Tang [2019] allows us to apply “dequantised” classical
sampling techniques to the data, and achieve a significant speedup over
previous classical methods. We discuss the issues around encoding classical
data into quantum systems in more depth in section 5.2.

Comparing state-of-the-art quantum computers with the best classical algorithm
for the same problem is not enough to establish a quantum advantage. For the
RQC used in Google’s claim of quantum supremacy Arute et al. [2019], classical
algorithms had not been optimised for the problem: it is not a useful problem
and there was little research on it. Subsequent work brought down the
classical runtime by some orders of magnitude. Similar advances have occurred
with boson sampling Zhong et al. [2020], Bulmer et al. [2022], Oh et al.
[2023].

Current research predicts a narrow window in terms of the number of qubits
where circuit model NISQ hardware can potentially outperform classical
hardware for certain problems Bouland et al. [2022]. Larger numbers of qubits
need longer gate sequences for the same algorithm. Errors inevitably build up
and render the computation unreliable without active error correction. Recent
proof-of-concept studies of quantum error correction Ryan-Anderson et al.
[2021], Sivak et al. [2023] begin to show errors suppressed for extended gate
sequences, heralding a new era of fault tolerant quantum computing.

##  3 Quantum simulation of quantum systems

Simulating quantum systems on classical computers has long been known for the
exponential explosion problem, where the size of the quantum space increases
exponentially with the number of qubits. This may be overcome using quantum
simulation Feynman [1982]: we use some controllable quantum system to study
another less controllable or accessible system. The method is especially
relevant for studying many-body problems Weimer et al. [2021] in condensed
matter physics, high-energy physics, atomic physics, quantum chemistry, and
cosmology. For authoritative reviews, see Georgescu et al. [2014], Altman et
al. [2021], Daley et al. [2022].

Table 3: Compare digital vs analogue quantum simulation. |  |  Analogue |  |  Digital
---|---|---|---|---
Hamiltonian |  |  Directly implement system Hamiltonian by experimentally building a controllable quantum system (simulator). |  |  Express Hamiltonian in second quantised form, then construct Hamiltonian out of quantum logic gates.
System vs simulator |  |  System may be completely different from simulator, e.g., using atoms in optical lattices to simulate electrons in a crystal. |  |  Implements circuits containing quantum gates, therefore has same level of technical difficulty as quantum computer.
Physical model |  |  Needs to match the simulator Hamiltonian within the available experimental controls. |  |  In principle, any quantum system can be simulated using digital (universal) quantum simulation.
Hardware |  |  Easier to implement for certain problems well matched to the simulator hardware. Likely to become practicable before universal quantum computers become available. |  |  Requires high depth circuits (many quantum gates) therefore unlikely to be practical before fully error corrected quantum computers are available.

Quantum simulation can be implemented using quantum computers (digital quantum
simulation) or with analogue devices (analogue quantum simulation) that
require less control and are easier to construct. In the latter, ultracold
atomic systems are often used to simulate condensed matter physics Trotzky et
al. [2012], Choi et al. [2016] and quantum chemistry Argüello-Luengo et al.
[2019] problems. We summarise the relative merits of the different platforms
in Table 3. In the next subsections, we describe the simulation techniques and
applications in more detail.

###  3.1 Hamiltonian evolution

The basic problem in most quantum simulations requires time evolution under
the system Hamiltonian. For Hamiltonian $\hat{H}$, the unitary operator
$\hat{U}(t)$ evolves the dynamics for a time $t$,

$\hat{U}(t)=\exp\left(-i\int_{0}^{t}\hat{H}(t)dt\right),$ (37)

using units where $\hbar=1$. For many physical systems, the system Hamiltonian
$\hat{H}$ can be written as a sum of terms that describe local interactions,

$\hat{H}=\sum_{j=1}^{n}\hat{H}_{j}.$ (38)

In analogue quantum simulation, the system Hamiltonian $\hat{H}_{\text{sys}}$
is directly mapped Somaroo et al. [1999] on to the simulator Hamiltonian
$\hat{H}_{\text{sys}}\leftrightarrow\hat{H}_{\text{sim}}$.

Implementing $\hat{U}(t)$ on a digital quantum computer usually relies on
numerical integration, often called Trotterisation Hatano and Suzuki [2005],
i.e., breaking up the evolution time into small steps of duration $\Delta t$.
The first-order Trotter-Suzuki formula gives

$\hat{U}(\Delta t)=e^{-i\sum_{j}\hat{H}_{j}\Delta
t}=\prod_{j}e^{-i\hat{H}_{j}\Delta t}+O((\Delta t)^{2}).$ (39)

Here, acceptable precision comes at the expense of small $\Delta t$, which
requires many quantum gates, and hence extremely deep circuits, well beyond
NISQ capabilities for useful problem sizes.

Figure 14: Conceptual diagram of universal simulator on digital quantum
computing device Tacchino et al. [2020].

Figure 14 gives a conceptual illustration of how to use Trotterisation to
perform quantum simulation on a quantum computer. A physical model describes
the quantum state evolution in a physical space, $\Psi(t)$; this evolution can
be approximated to arbitrary precision by mapping the given model onto a spin-
type Hamiltonian which is easily encoded on to a qubit-based register. We use
the Suzuki-Trotter formula (equation (39)) to perform the time evolution. The
sequence of unitary operations can be programmed through a quantum circuit and
directly run on a quantum computer, giving the approximated evolved state as
output $\psi(t)$ starting from a given input state $\ket{\psi(0)}$. Although
time evolution can be readily realised using Trotter product formulas, long
evolution times can lead to deep circuits. Recent work has demonstrated that
variational approaches can speed up the calculation of Trotter terms Benedetti
et al. [2021]. Alternatively, tensor network methods can be more efficient
than Trotterisation. This is particularly relevant when optimising quantum
circuits for Schrödinger time evolution McKeever and Lubasch [2023a] and, by
extension, solve combinatorial optimisation problems in adiabatic quantum
computing McKeever and Lubasch [2023b].

###  3.2 Quantum algorithms for Monte Carlo

Monte Carlo methods have a long scientific history with some of the earliest
known applications in the Los Alamos project Eckardt [1987]. Monte Carlo
simulation uses random sampling and statistical modelling to estimate
mathematical functions and mimic the operations of complex systems. All Monte
Carlo methods follow the same framework: model a system as a series of
probability density functions, repeatedly sample from the functions, then
calculate the statistics of interest. The method requires large numbers of
samples to achieve reasonable accuracy, often on the order of millions of
samples which can take hours or even days. When applied to quantum systems
(quantum Monte Carlo), the classical computation requires less memory than a
full quantum simulation, trading memory for the time needed to perform the
samples.

Quantum algorithms that can speed up Monte Carlo methods have been discussed
since the early 2000s. If sampling is done using a Markov chain, i.e., a
random walk, quantum versions of random walks 2.3.6 can generally provide a
quadratic improvement in the sampling rate .Szegedy [2004]. This faster mixing
of quantum processes is at the core of many polynomial quantum speed ups, and
can be applied to a wide range of classical sampling problems Montanaro
[2015].

When the problem involves simulating a quantum system, further advantage can
be gained with a fully quantum algorithm. For example, the time-evolution
operator can be expressed in the form of a summation, where each term in the
sum corresponds to a deterministic quantum-circuit configuration Berry et al.
[2015], Childs and Wiebe [2012] or by sampling random circuits Campbell
[2019].

Methods for quantum systems involve calculating the transition amplitude
$\bra{\psi_{f}}e^{i\hat{H}t^{*}}\hat{O}e^{-i\hat{H}t}\ket{\psi_{i}}$ with
initial state $\ket{\psi_{i}}$, final state $\ket{\psi_{f}}$, operator
$\hat{O}$, Hamiltonian $\hat{H}$. This is used to estimate the expectation
value of the operator $\hat{O}$. Both real ($t$) or imaginary ($it$) time
evolution can be used. Using imaginary time involves decaying rather than
oscillating functions, which can converge faster. In the following example
using imaginary time, the ground-state energy of an interacting Hamiltonian is

$\lim_{t\to\infty}\frac{\bra{\psi_{0}}e^{\hat{H}t/2}\hat{H}e^{-\hat{H}t/2}\ket{\psi_{0}}}{\bra{\psi_{0}}e^{\hat{H}t}\ket{\psi_{0}}}$
(40)

where the trial ground state $\ket{\psi_{0}}$ overlaps with the true ground
state.

The many diverse applications of Monte Carlo to quantum systems provide
corresponding opportunities for more efficient computation using quantum
computers (for reviews, see for example Ceperley and Alder [1986], Foulkes et
al. [2001]).

###  3.3 Quantum computational chemistry

One of the original motivations for developing quantum computers was to be
able to accurately simulate and characterise systems of interacting fermions
Abrams and Lloyd [1997], Ortiz et al. [2001]. Obtaining accurate solutions of
the electronic structure of many-body systems is a major challenge in
computational chemistry. Research has been carried out to develop algorithms
to solve the electronic structure problem using quantum computers Lin and Tong
[2022], Huggins et al. [2022], Somma [2019]. Fault-tolerant quantum computers
are predicted to solve the electronic structure problem for many-body systems
Cao et al. [2019] in polynomial time, for example, by using the quantum phase
estimation algorithm Aspuru-Guzik et al. [2005], Abrams and Lloyd [1997],
Kitaev [1996] (section 2.3.2). It is possible that a quantum advantage exists
even in the absence of a clear polynomial scaling. A detailed overview of the
state-of-the-art methods developed in this field can be found in Cao et al.
[2019] and McArdle et al. [2020].

In practice, the possibility of achieving exponential quantum advantage for
quantum chemistry problems remains controversial Lee et al. [2023]. It is
challenging for quantum simulations to achieve quantum advantage over the best
classical computations, given the decades of work already done on classical
methods for quantum systems. As classical algorithms continue to improve,
quantum information processing has been useful in developing advanced
“quantum-inspired” classical techniques Arrazola et al. [2019]. In this
regard, we review new techniques for both classical and quantum hardware from
different perspectives. We cover a range of methods applied to problems such
as ab initio simulations, Hubbard models, spin Hamiltonians, interacting Fermi
liquids, solid-state systems, and molecular systems.

Modern computational chemistry aims to simulate molecules and materials within
1 kcal/mol accuracy, which is required to predict chemical phenomena at room
temperature. Despite the tremendous progress in quantum chemistry algorithms
over the past century, the precision of 1 kcal/mol can only be achieved for
small systems due to the polynomial scaling (higher than O($x^{4}$) with
system size $x$) of advanced computational chemistry methods Grimme and
Schreiner [2018]. Hence, current research has a strong emphasis on developing
hybrid quantum-classical and quantum algorithms; estimating quantum resources;
and designing quantum hardware architectures for simulations of practical
interest.

Even with noisy hardware, the long-term impact of using quantum simulation for
quantum chemistry and materials is promising, given the recent efforts using
quantum Monte Carlo methods (see Section 3.2). For example, Google’s Sycamore
quantum processor used a hybrid quantum-classical Monte Carlo method Arute et
al. [2019] to calculate the atomisation energy of the strongly correlated
square H${}_{\text{4}}$ molecule. The results are competitive with those of
state-of-the-art classical methods. The algorithm relies on being able to
efficiently prepare a good initial guess of the ground state on the quantum
hardware. This is a non-trivial step.

Understanding the strengths and weaknesses of classical algorithms can
indicate where we may see a quantum advantage. It also indicates where issues
may arise in quantum simulation approaches in the near and long term. For
example, Tubman et al. Tubman et al. [2018] used classical simulations to
provide resource estimates for crucially important but often neglected aspects
of quantum chemistry simulations. These include efficient state preparation
for quantum methods such as phase estimation.

One promising research direction is to improve existing quantum methods for
estimating electronic ground-state energies. Two of the most popular methods,
quantum phase estimation (QPE) Aspuru-Guzik et al. [2005] (Section 2.3.2) and
VQEs Kandala et al. [2017] (Section 2.5), require significant quantum
computing resources Gonthier et al. [2022]. This motivates research to find
less expensive alternatives. For example, while VQEs do not require quantum
error correction protocols, it is challenging to choose a suitable ansatz Tang
et al. [2021], optimise potentially thousands of parameters McClean et al.
[2018], and have suitable noise resilience protocols Fontana et al. [2021],
Wang et al. [2021a].

Quantum phase estimation can help accurately simulate long-time dynamics and
determine eigenstate energies Kitaev [1996], but requires a large number of
qubits and quantum gates Wecker et al. [2014], Babbush et al. [2018]. It has
been experimentally implemented for small molecules on a variety of quantum
computing architectures, including superconducting qubits O’Malley et al.
[2016]. However, larger molecules of practical interest require a gate depth
several orders of magnitude beyond what is currently feasible Reiher et al.
[2017]. Therefore, chemical researchers are keen to improve techniques for
these types of simulations, including considerable efforts to estimate the
resources required to perform QPE simulations Babbush et al. [2018].

The partition function is a central quantity in statistical mechanics and
materials science. All thermodynamic quantities of interest can be derived
from the partition function, which then help us understand the behaviour of
many-body systems. Calculating partition functions of physical systems is
generally a #P-hard complexity class problem Lidar [2004], Jaeger et al.
[1990] – both classical and quantum methods are unlikely to provide exact
solutions. Quantum algorithmic attempts to estimate partition functions
involve quantum phase estimation Temme et al. [2011] and quantum walks Wocjan
et al. [2009], Arunachalam et al. [2022]. Hybrid quantum-classical algorithms
Wu and Wang [2021] and quantum-inspired algorithms Jackson et al. [2023b] have
also been proposed, with quantum algorithms potentially able to do better than
classical.

####  3.3.1 Quantum chemistry’s many-electron problem

The core objective of quantum chemistry is to solve the stationary Schrödinger
equation for an interacting electronic system in the static external field of
the nuclei, which describes isolated molecules in the Born-Oppenheimer
approximation. The problem can be defined as follows. There are $N$ electrons
in the field of stationary $M$ point charges (nuclei). Electrons interact only
through electrostatic Coulomb terms and are collectively described by a wave
function with $4N$ coordinates $\vec{x}_{i}=(\vec{r}_{i},\vec{\omega}_{i})$,
i.e., $3N$ spatial coordinates $\vec{r}_{i}$ and $N$ spin coordinates
$\omega_{i}$. The antisymmetric electronic wave functions are solutions of the
time-independent Schrödinger equation

$\displaystyle\hat{H}\Psi(\vec{r})=E\Psi(\vec{r})$ (41)

where we describe the motion of the $N$ electrons ($i$, $j$ indices) in the
field of the $M$ nuclei ($k$ index) with Hamiltonian McArdle et al. [2020]

$\displaystyle\hat{H}=\underbrace{-\sum_{i}\frac{\nabla_{i}^{2}}{2}}_{=\hat{H}_{K}}\underbrace{-\sum_{i,k}\frac{z_{k}}{|\vec{r}_{i}-\vec{r}_{k}|}}_{=\hat{H}_{e\text{-}n}}+\underbrace{\frac{1}{2}\sum_{j\neq
i}\frac{1}{|\vec{r}_{i}-\vec{r}_{j}|}}_{=\hat{H}_{e\text{-}e}}$ (42)

in atomic units $z$. The Hamiltonian contains the kinetic term $\hat{H}_{K}$,
the electron-electron interaction $\hat{H}_{e\text{-}e}$, and the electron-
nuclei interaction $\hat{H}_{e\text{-}n}$. When there are two or more
electrons, we cannot solve the Schrödinger equation exactly. This is the
fundamental challenge of computational chemistry.

There are two main families of approximations to solve this problem: Hartree-
Fock (HF) Slater [1951] methods with refinements and density functional theory
(DFT) Hohenberg and Kohn [1964], Kohn and Sham [1965]. Each has their own
opportunities for quantum enhancements. We summarise the main methods and
their limitations before discussing potential quantum enhancements.

####  3.3.2 HF and the electron correlation problem

The Hartree-Fock (mean-field) method Echenique and Alonso [2007] describes the
potential experienced by each electron resulting from the average field due to
all electrons. This means that the $\hat{H}_{e\text{-}e}$ term in equation
(42) is replaced by a mean-field potential to make the equation solvable.
However, this approximation does not correctly describe the electron-electron
correlation. Electron correlation accounts for the fact that electrons do not
move independently in the average field of other electrons rather, their
movements are influenced by the instantaneous positions and motions of other
electrons. The particle exchange energy resulting from the anti-symmetric wave
function (i.e., Pauli principle)

$\Psi(\vec{x}_{1},...,\vec{x}_{i},\vec{x}_{j},...,\vec{x}_{N})=-\Psi(\vec{x}_{1},...,\vec{x}_{j},\vec{x}_{i},...,\vec{x}_{N}),$
(43)

means we cannot express the wave function as a separable product of single-
particle wave functions,

$\Psi(\vec{x}_{1},\vec{x}_{2},...,\vec{x}_{N})\neq\psi_{1}(\vec{x}_{1})\psi_{2}(\vec{x}_{2})...\psi_{N}(\vec{x}_{N})$
(44)

where the spatial and spin parts are defined as
$\psi_{i}(\vec{x}_{i})=\phi_{i}(\vec{r}_{i})\chi_{i}(\vec{\omega}_{i})$.
Instead, the HF method approximates the many-body wave function with a Slater
determinant

$\Psi(\vec{x}_{1},\vec{x}_{2},...,\vec{x}_{N})=\\\
\frac{1}{\sqrt{N!}}\begin{vmatrix}\psi_{1}(\vec{x}_{1})&\psi_{2}(\vec{x}_{1})&...&\psi_{N}(\vec{x}_{1})\\\
\psi_{1}(\vec{x}_{2})&\psi_{2}(\vec{x}_{2})&...&\psi_{N}(\vec{x}_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}(\vec{x}_{N})&\psi_{2}(\vec{x}_{N})&...&\psi_{N}(\vec{x}_{N})\end{vmatrix},$
(45)

where the set of $\\{\psi_{i}(\vec{x}_{j})\\}$ are chosen to minimise the
ground state Hartree-Fock energy

$E_{0}=\bra{\Psi_{0}}\hat{H}\ket{\Psi_{0}}.$ (46)

The wave function depends on the Hamiltonian, and, in turn, the mean-field
part of the Hamiltonian depends on the wave function. This means that the HF
method must be solved iteratively using the self-consistent field (SCF)
method. We start from an initial guess for the wave function, build an initial
Hamiltonian, and then calculate the initial energy. We change the wave
function, update the Hamiltonian, and calculate the energy until the
difference between cycles is lower than a defined threshold. The SCF method
relies on the variational principle, which states that the energy of any
approximate wave function is higher than or equal to the exact energy. The
final energy is the exact HF energy $E_{0}$. This value is always higher than
the true ground-state energy because the formalism cannot fully account for
electron correlation effects. Recovering the correlation energy is the main
goal of post-HF methods.

For small systems, such as atoms, the unknown $\psi$ functions in equation
(45) can be solved numerically on a three-dimensional grid of points. Larger
systems require grids that are too dense to reach convergence. Hence we
express the $\psi$ functions as a linear combination of a set of known
functions, the basis set:

$f(x)=\sum_{i=1}^{M}c_{i}\chi_{i}(x).$ (47)

In practice, the basis set could be built from any type of function $\chi_{i}$
following two criteria: they should reproduce the behaviour of the electron
density (the highest density should be located at the atom positions, and they
should monotonically decrease to zero at an infinite distance from the nuclei)
and they should be easy to integrate. Among the most widely used basis sets,
there are Gaussian and Slater-type basis sets for molecules Molecular Sciences
Software Institute and Pacific Northwest National Lab/Environmental Molecular
Sciences Laboratory (PNNL/EMSL)(2023) [MolSSI] and plane waves for periodic
systems Singh and Nordström [2006].

Selecting the appropriate basis set is a critical consideration. It is a
balance between achieving the desired accuracy and managing computational
resources. Classical computer simulations often employ extensive basis sets,
which can provide high levels of accuracy, including up to 30 functions per
atom Peterson [2007], Jensen et al. [2016]. The computational capacity of
classical computers typically means that the basis set size is not the primary
limitation on accuracy. On the contrary, NISQ devices impose limitations on
the basis set size due to hardware constraints Chien and Klassen [2022],
making the basis set size a more significant factor affecting accuracy as
McArdle et al. discuss in this comprehensive review McArdle et al. [2020].

####  3.3.3 Density functional theory (DFT)

Density functional theory (DFT) Capelle [2006] focuses on the electronic
density

$\displaystyle n(\vec{r})=N\int d^{3}r_{2}\int d^{3}r_{3}...\int
d^{3}r_{N}|\Psi(\vec{r},\vec{r}_{2},...,\vec{r}_{N})|^{2}$ (48)

rather than electron coordinates. DFT systematically maps the many-body
problem to a single-body problem through the particle density $n(\vec{r})$
from which we calculate all other observables.

The Hohenberg-Kohn (HK) theorem lies at the heart of DFT. It states that the
system’s ground-state energy is uniquely determined by the ground-state
density $n_{0}(\vec{r})$. The ground-state wave function
$\Psi_{0}(\vec{r}_{1},\vec{r}_{2},...,\vec{r}_{N})$ must both reproduce the
ground-state density $n_{0}(\vec{r})$ and minimise the energy. Hence the aim
is to minimise the sum of kinetic, interaction and potential energies

$\displaystyle E_{0}[n(\vec{r})]$ $\displaystyle=\min_{\Psi\to
n}\bra{\Psi}\hat{H}_{K}+\hat{H}_{e-e}+\hat{H}_{e-n}\ket{\Psi}$ (49)
$\displaystyle=\underbrace{\min_{\Psi\to
n}\bra{\Psi}\hat{H}_{K}+\hat{H}_{e-e}\ket{\Psi}}_{=:F[n]}+\int
d^{3}r\;n(\vec{r})v(\vec{r})$ (50)

with Hohenberg-Kohn density functional $F[n]$. We then minimise $E_{0}$ with a
variational method using Lagrange multipliers to obtain the Kohn-Sham
equations Baseden and Tye [2014],

$-\frac{\nabla^{2}}{2}\psi_{i}(\vec{r})+v_{\text{eff}}(\vec{r})\psi_{i}(\vec{r})=\epsilon_{i}\psi_{i}(\vec{r}).$
(51)

The effective potential $v_{\text{eff}}(\vec{r})$ is defined as the sum of the
external potential (the interaction with the nuclei)
$v_{\text{ext}}(\vec{r})$, the Coulomb interaction among the electrons
$v_{\text{H}}(\vec{r})$, and the exchange-correlation potential
$v_{\text{xc}}(\vec{r})$:

$v_{\text{eff}}(\vec{r})=v_{\text{ext}}(\vec{r})+v_{\text{H}}(\vec{r})+v_{\text{xc}}(\vec{r}).$
(52)

The exchange-correlation potential takes into account non-classical effects
such as the Pauli principle and the electron correlation. This functional,
which must be parameterised, determines the quality of the result.

DFT is system-specific and does not reach the desired chemical accuracy.
Similarly to the Hartree-Fock method, it does not take into account direct
electron correlation, although this can be compensated for by the functional.
Additionally, this method has the disadvantage of self-interaction, which
occurs when each electron interacts with the electron density generated by all
the electrons, including itself. This interaction is exactly cancelled out in
the Hartree-Fock approach.

Because DFT is based on the three-dimensional distribution of the electron
density, it is possible to directly calculate it on a three-dimensional grid
around the atoms and bonds that represent the electron density at different
points in space. However, in practical implementations, DFT frequently uses
auxiliary functions to discretise the electron density on the grid. This
approach ensures precise numerical integration, computational efficiency,
management of exchange-correlation functionals, and allows us to apply linear
scaling methods for larger systems. Consequently, DFT uses basis set functions
similar to those in HF theory.

####  3.3.4 Post Hartree-Fock methods

Since neither HF nor DFT directly take into account the electron-electron
correlations, the repulsive force between electrons needs to be added a
posteriori. Post-Hartree-Fock (post-HF) methods are more accurate than HF, but
at greater computational cost. For example, the full configuration interaction
(FCI) method considers all possible electron configurations (excitations)
within the chosen basis set and gives the exact solution of the Schrödinger
equation within a basis set expansion Eriksen and Gauss [2021]. The number of
all possible occupancy combinations scales as $O(N^{M})$ with $N$ electrons
and $M$ basis functions. It can only be performed on classical computers for
small systems. Another example is Møller-Plesset perturbation theory Møller
and Plesset [1934] which takes the core terms $\hat{H}_{K}+\hat{H}_{e-n}$ as
an unperturbed Hamiltonian and treats the electron repulsion terms
$\hat{H}_{e-e}$ as a perturbation. Applying the Rayleigh-Schrödinger
perturbation theory gives an approximate correction of the correlation energy
to the HF energy. In comparison, coupled cluster techniques Čížek [1966] apply
an exponential cluster operator $\hat{T}=\sum_{k}\hat{T}_{k}$ containing
$k$-fold excitation operators $\hat{T}_{k}$ to the HF wave function
$\ket{\Psi_{0}}$. We obtain the exact solution $e^{\hat{T}}\ket{\Psi_{0}}$
through perturbation theory. These methods can be used for larger systems in
their truncated form. For example, instead of the FCI approach, we could
consider a single excitation to reduce the computational cost. However, this
will also decrease the accuracy of the final energy. Another approach is to
reduce the number of electrons and orbitals considered in the CI expansion by
selecting a complete active space (CAS) where all excitations are considered.
This approach relies on the premise that chemical reactions primarily involve
higher-energy electrons. Although we may not obtain the exact total energy,
energy differences (such as those between reactants and products) should
account for the electron correlation.

The discussion above focuses on the simulation of molecules. The same concepts
apply when simulating crystalline materials when the basis set consists of
plane waves Singh and Nordström [2006] or Bloch functions built from local
orbitals Hoffmann [1987]. Crystalline orbitals, which are the periodic
equivalent of molecular orbitals, are delocalised over the cell and,
therefore, are difficult to use in a CAS approach. One method to solve this
problem is to use the Wannier single-particle electron basis Wannier [1937]. A
detailed analysis of how these functions can be used in the context of quantum
algorithms can be found in Clinton et al. [2022].

####  3.3.5 Post Hartree-Fock quantum algorithms

Quantum computing holds significant promise for enhancing two primary
categories of chemistry challenges: post-Hartree-Fock (post-HF) methods and
optimisation/global minimum search tasks. Both categories deal with the
challenges of combinatorial explosion. The former arises when encompassing all
possible excitations, as seen in the FCI approach, while the latter emerges in
scenarios where a multitude of configurations need to be evaluated, such as in
determining the thermodymamic properties of complex materials or the optimum
folded structure of particular proteins. Efforts have been made to map
chemical problems to quantum computing-based simulations Cao et al. [2019],
Kassal et al. [2011]; understand the possible benefits of using quantum
hardware Bauer et al. [2020]; and explore the possibility of QML Bernal et al.
[2022]. It is worth noting that the FCI method, while being a significant
improvement over HF or DFT for a given basis set, is still constrained in its
accuracy by the choice of basis set.

Another difficulty arises when simulating electrons due to their fermionic
properties (half-integer spin). When mapping the wave function to the qubit
system, one normally starts from the second quantised form of the wave
function, where the computational basis state is a Slater determinant in the
Fock space. Then the Jordan-Wigner or Bravyi-Kitaev mappings Tranter et al.
[2018] are used to represent the correspondence between each fermionic state
and qubit. The latter step is the bottleneck: the proposed methods have a
scaling with respect to $N$ basis functions of between $O(N^{4})$ Xia et al.
[2018] and $O(N^{2})$ Babbush et al. [2014] depending on the basis set used.

So far, gate-based quantum hardware has been the technology of choice for
performing post-HF simulations. On NISQ devices, FCI is suitable for small
molecules such as H2 and LiH McArdle et al. [2020], while, for larger chemical
systems, an embedding scheme is generally preferred. Embedding techniques are
well-established in classical computing for quantum chemical systems, see for
example Catlow et al. [2017] for an overview of common approaches. The model
system is typically divided into two regions, an active region that is treated
at a high level of theory, using DFT or methods beyond DFT, embedded within an
environment or inactive (frozen) region. Different schemes are available for
embedding via systems of point charges, molecular mechanical models Catlow et
al. [2017], electronic density Berger et al. [2014], a density matrix Knizia
and Chan [2012] or Green’s function-based method Anisimov et al. [1997]. This
approach has high potential to be used with quantum computing systems, where
quantum algorithms are applied to the active region in the embedding model.
For example, Gao et al. Gao et al. [2021] selected an orbital active space
derived from an initial Hartree-Fock (HF) wavefunction to study the lithium
superoxide dimer rearrangement. The method incorporates relevant orbitals
present in reactants, products, and those involved in the transition state. In
another instance, Gujarati et al. Gujarati et al. [2023] chose active space
orbitals from DFT orbitals to study the water splitting reaction on magnesium
surfaces.

Accurately calculating electronic excited-state properties is another
challenging area of research Westermayr and Marquetand [2021]. This is
relevant when we want to simulate spectroscopic experiments involving high-
energy transitions above the ionisation potential, non-equilibrium dynamics,
light-matter interactions and other excited-state based phenomena. In these
cases, explicitly solving the time-dependent Schrödinger equation is the most
natural route Goings et al. [2018], Li et al. [2020]. Such a method requires
going beyond DFT and HF to properly account for electron correlation in real
time. The time-dependent density matrix renormalisation group method Baiardi
[2021] is one possibility but it may not perform well for the dynamics of
highly entangled states. This is an ongoing field of research.

####  3.3.6 Quantum optimisation for chemistry and materials science

Quantum annealing (section 2.4) has been demonstrated to be a valuable tool
for solving optimisation problems in chemistry, biochemistry Fedorov and
Gelfand [2021] and materials science Camino et al. [2023]. Given a set of
variables (for example, atom position in the solid solution lattice or bond
angle in protein folding), the aim is to find the configurations of those
variables that correspond to the lowest energy. The configuration space to
explore is proportional to $\binom{n}{k}$ with $n$ variables and $k$ possible
values that each variable can take. This configurational space rapidly becomes
difficult to sample using classical methods for large values of $n$. Quantum
annealing can offer an alternative means of finding the lowest energy
configurations by exploiting quantum superposition and tunnelling. See Camino
et al.’s Camino et al. [2023] tutorial for more details.

Quantum annealers have been used to solve the electronic Hamiltonian
eigenvalue-eigenvector problem Teplukhin et al. [2020] and calculate the
excited electronic states Teplukhin et al. [2021] for biatomic and triatomic
molecules. In industry, researchers at Volkswagen used D-Wave to calculate the
ground state energies of H2 and LiH as a proof-of-principle Streif et al.
[2019], where their mixed discrete-continuous optimisation algorithm found the
lowest eigenstate of the qubit coupled cluster method and used a quantum
annealer to solve the discrete part of the problem Genin et al. [2019]. Menten
AI used D-Wave to address the rotamer optimisation problem for protein design
without major simplifications or a decrease in accuracy Mulligan et al.
[2019]. GlaxoSmithKline studied the mRNA codon optimisation problem to help
drug discovery Fox et al. [2021] using quantum annealers. Fujitsu developed a
quantum-inspired annealing machine to screen chemicals for materials discovery
Hatakeyama-Sato et al. [2021].

Optimisation can also be performed on gate-based hardware when the problem is
formulated as a QAOA problem (section 2.5.1). For example, the design of
deuterated molecules for OLED applications has been studied using an approach
that combines DFT-derived energies and VQE Gao et al. [2023b].

###  3.4 High-energy physics quantum simulations

In 2016, the first full quantum simulation of a high-energy physics (HEP)
experiment demonstrated the creation of electron-positron pairs from energy
Martinez et al. [2016]. Gauge field theories are the underlying formalism
describing interactions among elementary particles in the Standard Model, and
can also be extended to physics beyond the Standard Model. Quantum simulators
can be used to study gauge field theories and investigate non-perturbative
dynamics at strong coupling. It is generally difficult to investigate the
latter using conventional Monte Carlo techniques because of the “sign problem”
Alexandru et al. [2022]. Quantum field theories are also used to describe
condensed matter many-body systems. Work on quantum simulation in this
setting, e.g., quantum algorithms Byrnes and Yamamoto [2006], and analogue
Zohar and Reznik [2011] and digital simulators Büchler et al. [2005], has
evolved into an experimental field that has attracted the attention of
particle and nuclear physicists Cloët et al. [2019], Bauer et al. [2023]. In
particular, after the first quantum degenerate Bose and Fermi gases were
experimentally created Anderson et al. [1995], DeMarco and Jin [1999], these
became a versatile analogue simulator to study the behaviour of condensed-
matter systems Schäfer et al. [2020], Daley [2023]. The quantum gas microscope
is another significant advance, allowing research into quantum-degenerate gas
phenomena where particle indistinguishability plays a key role Bakr et al.
[2009], Cheuk et al. [2015]. Classical numerical methods cannot simulate these
quantum phenomena, particularly for strongly-correlated fermionic systems and
(non-Landau-Ginzburg) topological matter. Rapid advances in atomic, optical,
molecular, and solid-state platforms Altman et al. [2021] have given us
analogue quantum simulators to explore HEP models. Better HEP simulations may
eventually help us understand difficult problems like quark confinement and
the properties of dense nuclear matter, e.g., neutron stars. Bauer et al.
Bauer et al. [2023] provide a roadmap of quantum simulations for HEP.

The Large Hadron Collider (LHC) experiments produce enormous amounts of data
as they search for new physics. This requires highly sensitive signal
detection methods. To discover the Higgs boson particle, researchers analysed
data from $10^{15}$ proton-proton collisions per experiment Wu and Yoo [2022].
The high-luminosity LHC (HL-LHC) programme will produce over 300 times as much
data, which will need sophisticated automated methods to extract interesting
signals. Quantum machine learning techniques may offer new methods for real
time big data analysis Wu and Yoo [2022].

While there have been proof-of-concept studies of quantum error correction
Ryan-Anderson et al. [2021], Sivak et al. [2023], digital quantum simulators
currently do not have sufficient resources to solve HEP problems at a
significant scale. In the NISQ era, VQE and QAOA platforms produce meaningful
results even with shallow circuits and without active error-correction
procedures (section 2.5). There has been work on simulating field theories,
including finding the lowest lying energy spectra of low-dimensional Abelian
and non-Abelian lattice gauge theories Kokail et al. [2019], Atas et al.
[2021], Irmejs et al. [2023]. It is important to note that many lattice gauge
theory simulations on NISQ hardware require fully or partially removing the
redundant degrees of freedom or by imposing symmetries. Simplifications like
these are essential in the resource-limited NISQ era.

###  3.5 Important applications in industry

Quantum simulation will potentially solve many industrial problems and bring
about significant societal and environmental benefits. We give two examples in
this section. Researchers are using the Fermi-Hubbard model as an analogue
quantum simulator for high-temperature superconductors, specifically copper
oxide compounds. On the other hand, digital quantum simulation may help us
understand how to produce ammonia fertilisers using less energy intensive
processes.

####  3.5.1 High-temperature superconductivity

Theoretical high-energy and condensed matter physics share common fundamental
concepts such as symmetry breaking, renormalisation group, and Feynman
diagrams. Interesting phenomena occur in strongly correlated electronic
systems where several physical interactions (spin, charge, lattice, and/or
orbital) are simultaneously active. In particular, the discovery of high-
temperature superconductors (HTS) in the 1980s Bednorz and Müller [1986],
Sheng and Hermann [1988] launched decades of diligent efforts to understand
and use these compounds Zhou et al. [2021], Lilia et al. [2022].

Figure 15: Crystal structure of cuprate superconductor YBa2Cu3O6+x Dagotto
[1994].

The Bardeen-Cooper-Schrieffer (BCS) theory Bardeen et al. [1957a, b]
introduced the idea of electron-hole Cooper pairs to explain conventional low-
temperature superconductivity. However, the pairing mechanism in
‘unconventional’ superconductors is still unknown. These materials do not
conform to conventional BCS theory or its extensions. Conventional
superconductors are used in many application areas such as MRI machines and
high speed trains, but they can only operate at relatively low temperatures.
This severely limits their wider use. One of the holy grails of materials
science is to develop superconducting materials that operate at room
temperature. This would revolutionise many technologies Lilia et al. [2022].
In particular, room-temperature superconductors would decrease the heat wasted
from electronic devices and allow them to run more efficiently. On a larger
scale, HTS can help achieve the International Energy Agency’s roadmap to
carbon-free economies International Energy Agency [2021] via nuclear fusion-
generated electricity Ball [2021].

Researchers have focused on the two-dimensional Fermi-Hubbard model as an
analogue quantum simulator. This model is believed to capture the important
behaviour of HTS Dagotto [1994], Fradkin et al. [2015], specifically the
cuprate supercondustors (a popular HTS) in the copper-oxygen planes (“CuO2
planes” in Figure 15). Despite the Fermi-Hubbard model’s apparent simplicity,
its theoretical properties are far from fully understood. It is difficult to
solve accurately as the model exhibits competing orders in its phase diagram
where it is most relevant to cuprates. Experimental work has involved building
analogue quantum simulators of the Fermi-Hubbard model using quantum dots
Hensgens et al. [2017], Wang et al. [2022] and ultracold atoms in optical
lattices Mazurenko et al. [2017], Hirthe et al. [2023].

The Fermi-Hubbard model is widely used to benchmark numerical methods for
strongly correlated systems LeBlanc et al. [2015]. Recent efforts explored
variational approaches Cade et al. [2020], Cai [2020], Stanisic et al. [2022],
while Google simulated the Fermi-Hubbard model using digital superconducting
quantum processors Google Quantum AI and collaborators [2020].

A state-of-the-art circuit-model algorithm for simulating the two-dimensional
Fermi-Hubbard time dynamics on an $8\times 8$ lattice requires roughly
$10^{7}$ Toffoli (CNOT) gates Kivlichan et al. [2020]. This includes the
overhead for fault tolerance, given the gate fidelities of current and near-
term hardware. A significant contribution to the gate count comes from the
phase estimation procedure. Clinton et al. Clinton et al. [2021] estimate the
resources required to simulate the Fermi-Hubbard time dynamics on a $5\times
5$ square lattice to $\leq 10$ % accuracy using 50 qubits and 1,243,586
standard two-qubit gates. This (optimistic) estimate assumes the effects of
decoherence and errors in the circuit can be neglected.

There is much work to be done before we have a complete microscopic theory of
HTS Lilia et al. [2022], Zhou et al. [2021]. The combination of theory,
simulation, materials synthesis and experiment has been crucial to progress in
the last two decades. For example, the breakthrough discovery of hydrogen-rich
superconductors Drozdov et al. [2015] may not have happened for another
century Castelvecchi [2023a] had there not been significant advances in
simulation and algorithms for chemical structure prediction Li et al. [2014].

Under high pressure, super-hydrides are hydrogen-rich compounds that exhibit
superconductivity at high temperatures. The high pressure pushes atoms so
close togather that they begin to interact through both their outer and inner
electrons. Understanding this electronic behaviour is considered one of the
top problems in the HTS field. In the search for other super-hydrides,
researchers are simulating combinations of different elements that combine