∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>e4e-mail<EMAIL_ADDRESS>
11institutetext: Atos Quantum Lab, 78340 Les Clayes-sous-Bois, France
22institutetext: Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France 33institutetext: Centre de Physique Théorique, 91120 Palaiseau, France

# Quantum computing with and for many-body physics

Thomas Ayral e1,addr1  Pauline Besserve e2, addr1, addr3  Denis Lacroix
e3,addr2  Edgar Andrés Ruiz Guzman e4,addr2

###### Abstract

Quantum computing technologies are making steady progress. This has opened new
opportunities for tackling problems whose complexity prevents their
description on classical computers. A prototypical example of these complex
problems are interacting quantum many-body systems: on the one hand, these
systems are known to become rapidly prohibitive to describe using classical
computers when their size increases. On the other hand, these systems are
precisely those which are used in the laboratory to build quantum computing
platforms. This arguably makes them one of the most promising early use cases
of quantum computing.

In this review, we explain how quantum many-body systems are used to build
quantum processors, and how, in turn, current and future quantum processors
can be used to describe large many-body systems of fermions such as electrons
and nucleons. The review includes an introduction to analog and digital
quantum devices, the mapping of Fermi systems and their Hamiltonians onto
qubit registers, as well as an overview of methods to access their static and
dynamical properties. We also highlight some aspects related to entanglement,
and touch on the description, influence and processing of decoherence in
quantum devices.

††journal: Eur. Phys. J. A

Received: date / Accepted: date

##  1 Introduction

For decades since they were envisioned by Richard Feynman in the 1980s
feynman2018 , quantum computers have been imagined as futuristic objects
overcoming the limitations of classical devices. Today, with significant
progress in the manipulations of various quantum systems Acin2018 , the
possibility of using them as a computational unit is becoming a reality
GYONGYOSI2019 ; Meter2016 . The race towards proving "quantum advantage" is
now underway Harrow2017 ; Arute2019 ; Zhong2020 . The ultimate challenge of
this race is to provide one or several reliable quantum processing units
(QPUs) with unprecedented capabilities in terms of hard memory storage or the
ability to solve specific complex problems in record time Arora2009 ;
Nielsen2010 .

Quantum computing is now at a turning point in its practical development
thanks to the growing availability of quantum machines SOEPARNO2021 . Yet,
these machines have a limited quality: The noise that degrades each operation
of a quantum circuit strongly constrains the complexity and types of
algorithms that can be used today, and requires specific denoising methods
Devitt2013 ; Gadioli2020 . This had led to the notion of Noisy Intermediate
Scale Quantum (NISQ) Preskill2018 ; Bharti2021 ; Endo2021  processors to
describe this intermediate stage of development. Despite the limitations of
the NISQ era, the possibility to experiment with actual quantum devices has
led to intensive scientific emulation to test the capacity of current
computers, propose new quantum algorithms, and ultimately prepare for the
coming second quantum revolution and surpass the limitations of classical
algorithms.

Quantum many-body systems—formed by a set of particles interacting with one
another—appear as natural test benches for quantum platforms McClean2016 ;
McArdle2018a ; Cao2019 ; Bauer2020 ; Tilly2021a ; Claudino2022 . This class of
problems is characterized by a Hilbert space size that increases steeply when
the number of one-body degrees of freedom (the number of particles or
accessible single-particle space or both) increases. This large size leads to
severe restrictions in the class of many-body systems that one can solve
exactly on classical computers.

This increase in complexity is common to all fields of physics or chemistry.
In particular, it applies to quantum devices themselves: assemblies of quantum
bits or "qubits" are also characterized by an exponential growth of the
spanned Hilbert space, as it is of size $2^{n_{q}}$ where $n_{q}$ is the
number of qubits. What is a hindrance for the classical description of many-
body systems may thus become a blessing: the fact that they contain an
exponential complexity a priori makes quantum computers suitable for tackling
many-body (and thus also exponentially complex) systems, provided such systems
can be accurately encoded and probed via a precise manipulation of qubits.
This is facilitated by the fact that particles treated in the formalism of
second quantization share many formal aspects with qubits (as will be
described in section 4). Many-body systems such as those encountered in
quantum chemistry McArdle2018a ; Cao2019 ; Bauer2020 ; Claudino2022 ;
Bassman2021 ; Fedorov2021a , condensed-matter physics Zeng2019 , atomic
physics Saffman2016 ; Bian2021 ; Graham2022 ; Bluvstein2022 , astrophysics
KerenLi2019 ; mielczarek2019 ; Czelusta2021 ; Mielczarek2021 ; Mielczarek2018
, or nuclear physics Zhang2021 ; Stetcu2022 ; Roggero2020 ; Hobday2022 ;
Romero2022 , have thus become key application domains of quantum computing.

The primary goal of the present article is to introduce quantum computing from
a double many-body physics perspective: quantum computers are many-body
systems that can help understand (among others) many-body problems. We start
by highlighting some specific aspects of selected many-body systems one might
find in nature and underline their common features (section 2) with a focus on
the complexity of treating them using classical computers. Section 3
introduces different types of quantum computing devices, namely, analog and
digital, to the nonexpert reader. This section is not only a discussion of the
basic concepts in quantum computation but also an opportunity to underline the
fact that quantum computers are built upon many-body interacting systems. In
section 4, we discuss various aspects related to the solution of quantum many-
body problems with quantum computers. We introduce selected quantum algorithms
to solve these problems, be it with post-NISQ (section 5) or NISQ (section 6)
processors. In section 7, we briefly discuss how entanglement in many-body
systems can be described on a quantum computer. Finally, in section 8, we
discuss the modeling of the noise impacting quantum hardware, how its effect
can be mitigated on current machines as well as the main principles of error
correction, which could bring about fault tolerance in the long term.

Figure 1.1: Illustration of some of the many-body systems where quantum
computing is now being explored as a disruptive technique compared to
classical computing: quantum chemistry, condensed matter, nuclear and neutrino
physics.

##  2 Quantum Many-body problems

In this section, we discuss various types of many-body systems where quantum
computing could be of help. The reader interested mainly in quantum computing
aspects can skip this section and directly go to section 3.

Many-body systems are encountered in many fields of physics. They are ruled in
all generality by a Hamiltonian that reads, in a second-quantized form:

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{\rm 1-body}+H_{\rm
2-body}+H_{\rm 3-body}+\cdots$ (2.1) $\displaystyle=$
$\displaystyle\sum_{\alpha\beta}h_{\alpha\beta}c_{\alpha}^{\dagger}c_{\beta}+\frac{1}{2}\sum_{\alpha\beta\gamma\delta}v_{\alpha\beta\gamma\delta}c^{\dagger}_{\alpha}c^{\dagger}_{\beta}c_{\gamma}c_{\delta}+
\cdots,$

where $\alpha,\beta,\gamma,\delta$ are multi-indices whose components depend
on the specific many-body system at hand. As we will see, depending on the
system, one may not need terms involving more than two particles.

The defining feature of many-body problems is that they cannot be solved
within the mean-field approximation. This approximation is obtained, at the
level of Hamiltonian (up to constant terms), by replacing all operators but
one by its average value in each term. In other words, in a many-body system,
the correlations between its key constituents—be they electrons, nucleons, or
spins—must be handled with sophisticated methods, hence their other name of
many-body systems: strongly correlated systems.

In the following subsections, we introduce a few fields—illustrated in Fig.
1.1—where many-body problems are found, with a focus on the form of the
Hamiltonians that describe them, and the typical classical methods that are
used to investigate their properties.

###  2.1 The electronic structure problem: electrons in solids and quantum
chemistry

Electrons in solids or molecules are interacting particles. In many solids,
the Coulomb interaction can be dealt with in an averaged fashion because
electrons do not come too close to each other due to the Pauli principle. Due
to this, Hartree-Fock (HF) theory, where electrons are treated as independent
particles in an average field, is often used as the building block for more
complex theoretical methods to treat many-body effects. Another powerful
technique that allows treating a priori all correlations while keeping the
independent particle or quasi-particle picture is Density Functional Theory
(DFT) Lipparini2008 .

However, in some solids and most molecules, such an averaged description of
Coulomb interactions leads to wrong predictions. For instance, in solids where
valence electrons are of $d$ or $f$ character (namely very localized), neither
HF nor DFT can predict the transition from a metallic Fermi liquid to a Mott
insulator. In such insulators, Coulomb interactions freeze the charge degree
of freedom. This type of solids has received much attention in the past forty
years since the high-temperature superconductors discovered in the mid-1980s
are believed to be doped Mott insulators. Despite prolonged efforts to crack
this problem, a complete theoretical account of the origin of high-Tc
superconductivity—and numerous other phenomena, especially when driving these
systems out of equilibrium—is still missing. This lack of explanation is due
to the complexity of dealing with such systems beyond mean field.

The simplest model describing such strongly-correlated solid-state systems is
the so-called (Fermi) Hubbard model, which, in its single-band version, reads:

$H=-t\sum_{\langle
ij\rangle}c_{i\sigma}^{\dagger}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}-\mu\sum_{i\sigma}n_{i\sigma},$
(2.2)

where $i,j=1\dots N$ denote lattice sites ($N$ is typically infinite in a
solid; $\langle ij\rangle$ denote nearest neighbors on the lattice), $t$
denotes a tunneling or hopping term, $U$ the on-site Coulomb interaction and
$\mu$ the chemical potential. The fermionic creation and annihilation
operators $c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ typically create and
annihilate electrons in localized orbitals $\phi_{i}(r)$ (sometimes called
Wannier orbitals). $n_{i\sigma}$ is a density operator which reads
$n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$. $\mu=U/2$ corresponds to half-
filling (undoped case), namely, one electron per site. Despite its apparent
simplicity, this "spherical cow" of strongly-correlated electronic systems is
difficult to solve on a classical computer in physically "interesting"
regimes, like the doped, low-temperature regime in two spatial dimensions and
the thermodynamical limit ($N\rightarrow\infty)$.

The Hubbard model can either be tackled "directly" via exact diagonalization
(ED), quantum Monte-Carlo (MC) or tensor-network methods, or "indirectly" via
embedding methods, which map the model to a smaller many-body problem as will
be described in subsection 4.2. The exponentially growing size of the Hilbert
space with $N$ quickly makes ED prohibitive, although advanced versions like
Krylov or Lanczos methods can help reach relatively large sizes Saad2011 . An
exponential sign problem typically plagues MC methods; namely, the statistical
error bar of MC methods grows exponentially with system size and decreasing
temperature, requiring an exponential number of samples and hence exponential
run time. Tensor-network methods (like Matrix Product States, MPS)
schollwock_2011 ; Ran2020 ; Cirac2021  can also be used but usually require a
space complexity (memory usage) that scales exponentially with the so-called
entanglement entropy $S$ (see section 7 for a more in-depth discussion). These
methods are potent when $S$ is constrained to small values, as happens e.g in
one dimension. However, in two and more dimensions, $S$ usually grows in a way
that makes these methods difficult to apply.

The Hubbard model proves insufficient to describe molecules in quantum
chemistry, as chemical systems are less prone to screening and thus cannot be
described by on-site interactions only. Thus, one typically deals with the
following more general Hamiltonian:

$H=\sum_{pq,\sigma}h_{pq}c_{p\sigma}^{\dagger}c_{q\sigma}+\frac{1}{2}\sum_{pqrs}\sum_{\sigma\sigma^{\prime}}v_{pqrs}c^{\dagger}_{p\sigma}c^{\dagger}_{q\sigma^{\prime}}c_{r\sigma^{\prime}}c_{s\sigma}$
(2.3)

where $p,q,r,s=1\dots N$ denote orbitals $\phi_{q}(r)$ (typically molecular
orbitals obtained after a first Hartree-Fock computation). The one- and two-
body matrix elements $h_{pq}$ and $v_{pqrs}$ are computed as integrals over
the molecular orbitals. Typically, $N$ is of the order of $10-100$.

The exact diagonalization of this Hamiltonian (ED, usually called Full
Configuration Interaction – FCI – in a quantum chemical context David1999 ) is
limited to a small number $N$ of orbitals. Variational methods that optimize
parameterized wave-functions to minimize the energy are less costly than FCI.
One simple example of such variational approaches is the HF method itself.
Another illustration of a widely used variational approach in chemistry is the
multi-configurational self-consistent field (MCSCF) method Hinze67 ;
schmidt1998 ; jensen2007 ; cramer2013 . Alternatively, non-variational methods
based on systematic expansions on top of the HF wave-function, like many-body
perturbation theory (MBPT) shavitt2009 ; Werner2009 , are also less costly
than FCI. A prominent example of such methods is the Coupled Cluster (CC)
method, based on the parametric many-body state (here limited to single and
double excitations) Bartlett2007 ; shavitt2009 ; Evangelista2018

$\displaystyle|\Psi(\vec{\theta})\rangle$
$\displaystyle=e^{T_{1}+T_{2}}\ket{\Psi_{\mathrm{HF}}}$
$\displaystyle=e^{\sum\limits_{ia,\sigma}\theta_{i}^{a}c_{i\sigma}^{\dagger}c_{a\sigma}+\sum\limits_{ijab,\sigma\sigma^{\prime}}\theta_{ij}^{ab}c^{\dagger}_{i\sigma}c^{\dagger}_{j\sigma^{\prime}}c_{a\sigma^{\prime}}c_{b\sigma}}\ket{\Psi_{\mathrm{HF}}},$
(2.4)

with $i,j$ (resp. $a,b$) empty (resp. occupied) orbitals. The non-variational
nature of the CC method comes from the non-unitary character of the so-called
excitation operator $e^{T_{1}+T_{2}}$, which implies that the variational
principle does not hold. A unitary, variational variant of the CC method,
called unitary coupled cluster (UCC), has been investigated for quantum
computing and is tackled in 6.2.5. This method is regarded as one of the most
advanced methods in chemistry (for systems with dynamical correlations, as
opposed to systems with large static correlations, where methods based on
Matrix Product States [see section 7] are among the most advanced Baiardi2020
). It can be combined with active-space methods, whose goal is to reduce the
number of relevant (or correlated or active) degrees of freedom, similar to
embedding methods for solids. This selection of degrees of freedom will be
further discussed in subsection 4.2.

###  2.2 Nuclear physics

Atomic nuclei are self-bound, strongly interacting systems with a wide range
of numbers of particles, from very few ($2$ for the deuteron) to several
hundreds for the heaviest nuclear systems existing in nature. Nucleons
organize themselves to form quantum droplets with a large variety of static
and dynamical physical phenomena Bohr1998 . These phenomena can be observed in
the laboratory through the use of accelerators. The many-body treatment of
nuclei is particularly complex due to the non-perturbative nature of the two-
body interaction with a strong repulsion at short distances between particles.

The nuclear Hamiltonian is of the form (2.1) where the multi-indices
$(\alpha,\beta,\gamma,\delta)=1,...,N$ label single-particle states
characterized by the usual quantum numbers $n,l,m,\sigma$, as well as an
isospin component $\tau$. These states can for instance be 3-dimensional
harmonic oscillators (HO) states with $\phi_{\alpha}({\bf
r})=\phi_{nlm\sigma\tau}({\bf r})$ ($\tau=-1/2$ and $+1/2$ for neutrons and
protons respectively), where the strength of the HO is optimized to reproduce
nuclei sizes Ring80 . Due to the presence of spin-orbit coupling, it is usual
to introduce the angular momentum $\vec{j}=\vec{l}+\vec{s}$, and relabel the
state with $(nljm\sigma\tau)$. Since there is no external field, the one-body
term only contains the kinetic component. The two-body interaction contains
nuclear (short-range) and Coulomb (long-range) interactions. The Coulomb part
acts only on protons, while the nuclear part, which depends on the spin of the
particles, acts on all nucleons. Note that the latter is almost the same for
all particles, a property know as isospin symmetry of the nuclear force.
Altogether, $v_{\alpha\beta\gamma\delta}$ is both spin and isospin-dependent.
In addition, the existence of 3-body and, more generally, multi-body
interactions was recognized only recently with the recent advances in
Effective-Field-Theory to construct nuclear interactions Epelbaum2009 ;
Machleidt2011 ; Hammer2013 ; Machleidt2016 ; Hammer2020 . The presence of
multi-body interactions is an extra complication compared with the electronic
structure Hamiltonian (2.3), and even in the most advanced many-body
techniques, such interactions are usually treated only approximately.

Despite this complexity, a variety of simplified Hamiltonians have been
proposed to understand specific properties of nuclei. A typical example is the
Lipkin-Meshkov-Glick model Lipkin1965 , which is often used to understand the
concept of spontaneous symmetry breaking in finite systems. Another example is
the pairing, also called Richardson Vondelft2001 , Hamiltonian that is often
used to understand superconducting effects. This Hamiltonian can be justified
by (i) assuming that only a set of particles, typically close to the Fermi
energy are active, (ii) the interaction between them is constant, and (iii) it
is non-negligible only when pairs of time-reversed states, denoted by
$(i,\bar{i})$, are involved. These states are often taken as opposite spin
particles or as particles with spin projection $j_{z}=m$ and $-m$ when a
single $j$-shell is considered as active. The Hamiltonian then reduces to:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i}\varepsilon_{i}(c^{\dagger}_{i}c_{i}+c^{\dagger}_{\bar{i}}c_{\bar{i}})-g\sum_{ij}c^{\dagger}_{i}c^{\dagger}_{\bar{i}}c_{\bar{j}}c_{j}.$
(2.5)

These simple, schematic Hamiltonians are studied today on quantum computers as
first steps towards future applications.

An overview of the microscopic approaches used to describe atomic nuclei can
be found in Ref. Lacroix2010 . The only approach able to describe the large
variety of phenomena ranging from static (nuclear structure), dynamical
(nuclear dynamics), and thermodynamical properties is nuclear Density
Functional Theory, often referred to as Energy Density Functional (EDF) theory
Bender2003 ; Nakatsukasa2016 ; Colo2020 . Another powerful approach,
restricted to studying nuclear structure properties, consists of performing a
direct CI method in a restricted subspace of single-particle states forming
the valence space Caurier2005 ; Nowacki2021 . In this approach, often referred
to as the Shell Model, the effective interaction is fine-tuned to account for
the truncation of the model space. One of the main difficulties that forces
the restriction to a set of active valence particles is the size of the many-
body Hilbert space when the number of single-particle states increases. The
current scope of restricted CI approaches is the treatment of eigenvalue
problems in spaces with $10^{11}-10^{12}$ states. These values are still far
from the requirements to treat the whole nuclear chart with all single-
particle active states.

Another breakthrough in the most recent description of interactions lies in
the possibility of getting rid of the hard core and using soft interactions
for low-energy nuclear physics problems Bogner2010 . Such interactions have
opened the way to the so-called ab-initio method that aims at treating the
nuclear many-body problem directly, starting from the bare Hamiltonian (2.1).
This advance has led to a significant boost in applying several many-body
techniques, some already mentioned in section 2.1. One can, in particular,
mention the use of the full CI technique—known, in this context, as the no-
core shell-model Quaglioni2008 ; Navratil2009 , the Green’s Function Monte-
Carlo method Pudliner1997 ; Wiringa1998 ; Wiringa2000 , the Self-Consistent
Green Function method Dickhoff2004 ; Soma2020 , the Coupled-Cluster method
Bartlett2007 ; Hagen2014 , or Many-Body Perturbation Theory Tichai2020 , among
others. One specific aspect of atomic nuclei is the necessity to generalize
some of these theories to allow for possible spontaneous symmetry breaking
like particle number or rotational symmetries (see, for instance, a few
examples of extensions in Refs. Tichai2011 ; Soma2011 ; Henderson2014 ;
Duguet2016 ; Demol2021 ). This generalization is fundamental to describing
open-shell nuclei.

Although significant progress has been made in recent years, ab-initio methods
are still applied to study the nuclear structure effects in a limited region
of the nuclear chart or to nuclear reactions involving only very light
systems. This limitation is due to the increment in the complexity of the
problem when the number of particles increases. For now, very few quantum
computing pilot applications on real devices have been made, and these have
been limited to rather simplistic nuclear Hamiltonians Dumitrescu2018 ; Lu2019
; Kiss2022 . However, the use of quantum computers for nuclear many-body
problems has recently gained momentum Roggero2020 ; Roggero2020b ; Lacroix2020
; Dimatteo2021 ; Siwach2021 ; Du2021 ; Siwach2021b ; Cervia2021 ; ruiz-
guzman2021 ; Hlatshwayo2022 ; Perezfernandez2022 ; RuizGuzman2022 ; Stetcu2022
; Siwach2022 ; Romero2022 ; Hodbay2022 ; Perez2023 ; Robin2023 ; Beck2023 .

###  2.3 Common difficulties

A defining property of many-body systems is that they are exponentially
difficult to solve on classical computers. As we will see in this section,
this exponential difficulty may trivially arise from the size of the Hilbert
space. However, most classical methods strive to circumvent this difficulty by
either using the structure of the problem to merely reduce the number of
relevant degrees of freedom or by adopting another representation—shifting the
exponential difficulty from the size of the Hilbert space to other parameters,
like, for instance, the severity of the Monte-Carlo sign problem or the
internal dimension of a tensor network.

One standard strategy to circumvent the exponential complexity is to use
parameterized variational wave function ansätze. Another strategy—based on so-
called reduced density matrices—relies on the assumption that some degrees of
freedom contain more information than others. A typical starting point is to
assume that one-body degrees of freedom are the most relevant. The information
on them is contained in the one-body reduced density matrix (1-RDM). Reducing
the information to these DoFs leads to Hartree-Fock (HF) or mean-field theory.
Usually, such simplified approaches miss significant correlation effects, and
approximations beyond the mean-field are necessary to describe many-body
systems accurately. For instance, this description can be done by truncating
the Bogolyubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy and treating 2-body
or higher DoFs explicitly through the two-body reduced density matrix (2-RDM)
or higher-order reduced density matrices Bonitz2016  (see also the discussion
of embedding or active-space methods in section 4.2).

Despite the dimensionality reduction that the aforementioned approximate
methods afford, they always, at some point, reach a computational limit on
classical computers. Quantum computers—provided fermionic problems can be
turned into spin or qubit problems (see section 4.1)—can a priori overcome
these limitations provided enough qubits can be efficiently manipulated.
Interestingly, quantum computers do not start from a blank page: many quantum
algorithms are strongly guided by the accumulated expertise and methods gained
on classical devices. An illustration of that will be discussed in section 4.

Besides the above general considerations on many-body systems, each physical
system has specificities that will render its encoding on quantum computers
more or less difficult. For instance, electrons can have two spin components
(spin up and down), while nucleons can have both spin and isospin (neutron and
proton components). Some physical systems, such as solids or atomic ones, can
be suitably described on a lattice, sometimes with only nearest neighbor two-
body interaction. This case is advantageous when encoding such a problem on a
set of qubits with limited connectivity. Some other systems might be more
complex, like atomic nuclei, with long-range two-body or, more generally,
multi-body interactions or the necessity to describe unbound states. A
prerequisite to the success of future applications is the efficient
transposition of a given many-body problem with good mapping of its
characteristics into an analog or a digital quantum computer (see section 4).
The complexity of this transposition will strongly depend on the problem
itself, but many-body problems will undoubtedly be among the first stringent
benchmarks for current and future quantum technologies.

The exponential size of the Hilbert space is not a sufficient condition for
making a problem computationally hard from classical computers. For instance,
the ground state of quadratic fermionic Hamiltonians—namely Hamiltonians that
are bilinear in the creation and annihilation operators—can be found in
polynomial time on a classical computer (see e.g Bruus2002 ). The problem
becomes complicated only when terms beyond quadratic terms are added. This
result explains why certain classes of time evolutions—quantum circuits (see
section 3.1.2 below for a definition) stemming from a quadratic Hamiltonian,
with so-called matchgates—are also simulatable classically in polynomial time.
Interestingly, other circuits, this time unrelated to "uncorrelated systems",
are efficiently simulatable. For instance, Clifford circuits, with gates
belonging only to the Clifford group Gottesman1998 , are simulatable with time
and space complexity $\mathcal{O}(n^{2})$ (with $n$ the number of qubits), a
result known as the Gottesman-Knill theorem Gottesman1998a ; Aaronson2004 .
This result holds even though these circuits can generate highly entangled
states. The key idea behind the efficiency of the simulation of such circuits
is that the states generated by Clifford circuits can be represented in a
compact (polynomial) fashion.

Leveraging a compact representation is also very common in many-body classical
methods. For instance, some Monte-Carlo methods decouple quartic (interaction)
terms in the Hamiltonian to eliminate this term at the expense of an
additional auxiliary field BAEURLE2004 . Then, the exponential difficulty is
shifted from the size of the Hilbert space to the Monte-Carlo sign problem. In
tensor network methods, like Matrix Product States schollwock_2011 , the
accuracy of the representation is tuned by the size of the internal indices
(often called the bond dimension), which is related to the degree of
entanglement of the state at stake (see section 7.3).

###  2.4 A few quantum complexity considerations

In this section, we comment on the formal expectations regarding speedups that
quantum computers can bring when solving many-body problems. For general
reviews on the topic, we refer the reader to e.g Watrous2012a ; Nielsen2010 .

Classical computational problems with a yes/no answer (decision problems) are
classified using complexity classes. For instance, finding the ground state
energy of a Hamiltonian $H$ acting on $n$ qubits (fermions) can be formulated
as a decision problem by picking numbers $a$ and $b$ such that $\epsilon\equiv
b-a>1/\mathrm{poly}(n)$ defines a region where it is promised that the ground
state energy does not lie in, and asking whether $E_{0}<a$ or $E_{0}>b$ (see
e.g., Watson2022  and references therein).

Figure 2.1: Computational complexity classes and selected problems. The three
crosses indicate three examples of complex problems. The "factoring" and
"traveling salesman" are common problems often quoted in complexity theories.
The "k-local Hamiltonian ground state estimation" is a common problem in many-
body systems (see text).

The two main classes of classical complexity are P and NP, which describe
problems that can be solved in polynomial time or whose solution can be
verified in polynomial time. NP-hard problems are problems at least as hard as
any problem in NP, and NP-complete problems are the NP-hard problems that
belong to NP. NP-complete problems are very likely (that is unless P=NP) not
to have a polynomial-time solution, i.e., they are colloquially referred to as
exponentially hard problems.

One central question is whether quantum computers can reduce the complexity of
solving many-body problems Bernstein1997 . To answer this question
systematically, two quantum complexity classes, BQP (bounded-error quantum
polynomial time) and QMA (quantum Merlin–Arthur), have been designed as
quantum counterparts to P and NP, respectively: the computer that solves the
problem or, respectively, certifies the solution of the problem is a quantum
computer. Fig 2.1 illustrates these classes and the connection between quantum
and classical complexity classifications.

A significant result is that NP-complete problems are very likely not in BQP.
In other words, it is improbable that quantum computers can solve
exponentially hard classical problems in polynomial time. The fact that the
factoring problem can be solved using Shor’s algorithm Shor1994  with an
exponential speedup is possible because factoring is not a NP-complete
problem. Importantly, this also does not mean that quantum computers are not
helpful for NP-complete problems: quantum heuristics (polynomial algorithms
with no quality guarantee) can still reach better solutions than classical
heuristics.

A natural question to know what to expect from quantum computers for many-body
problems is to which complexity class they belong. The ground state estimation
problem is QMA-complete for $k$-local (spin) Hamiltonians (Hamiltonians whose
terms act on at most $k$ qubits, see Eq. (5.5) below) as long as $k\geq 2$
Kempe2008 . QMA-completeness is retained for the ground state estimation
problem of geometrically 2-local Hamiltonians – that is to say 2-local
Hamiltonians only involving pairs of adjacent qubits – in a square lattice
qubit layout (aka a quantum spin glass) Oliveira2008 . These statements hold
for general values of the couplings of the spin Hamiltonians mentioned above.
Restrictions on the coefficients can lead to a reduction in computational
complexity. For instance, the transverse Ising model with negative transverse
field and ferromagnetic interactions (i.e., Eq. (3.4) below with $\Omega<0$
and $C<0$ [which is not the case for Rydberg atoms, where $C>0$]) is a so-
called stoquastic Hamiltonian. Its ground state energy can be approximated
polynomially in the system size $n$ and $1/\epsilon$ Bravyi2015 .

As for fermionic and bosonic models, Schuch2009  showed that the Fermi-Hubbard
model with local magnetic fields (with an additional
$\sum_{i}\bf{\sigma}_{i}\cdot\bf{B}_{i}$ term in Eq. (2.2)) is QMA-complete,
while Childs2013  showed that the Bose-Hubbard model is QMA-complete.

These formal considerations call for two comments. First, despite these
hardness results, the ground state energy of the models cited above can be
determined with great accuracy on classical computers in special regimes,
i.e., for certain classes of parameters. For instance, this is the case of the
Fermi-Hubbard model on a bipartite lattice, whose solution with quantum Monte-
Carlo does not suffer from a sign problem at half-filling (see, e.g.,
Leblanc2015 ). The existence of these special regimes means that one must be
careful to look for truly hard classical computational regimes to identify a
useful application of quantum computers. In other words, classical methods
make the most of any symmetry or structure of the problem to overcome the
underlying exponential complexity of the many-body problem so that truly hard
regimes are hard to come by. In these regimes, quantum computers ought to also
leverage these symmetries and structures in order to outperform classical
computers.

Second, even in those regimes where classical computers fail to reach an
accurate enough result, the QMA-completeness of the problem is a strong
indication that the problem will also be hard to solve on a quantum computer.
In other words, classical and quantum algorithms likely end up running into an
exponential wall (see, e.g., Lee2022  for a concrete example). This limitation
is not necessarily a showstopper: what matters is whether quantum computers
can reach regimes inaccessible to classical computers before they run into
said wall.

Lastly, complexity theory can also be used to appraise, at least formally, the
feasibility of hybrid quantum-classical algorithms like the Variational
Quantum Eigensolver (see section 6.1.1 below). For instance, the classical
optimization procedure of the energy $E(\bf{\theta})$ in VQE is generically a
challenging computational problem Bittel2021 . In practice, this does not
exclude the existence of heuristics for finding accurate enough variational
parameters for concrete (as opposed to generic) problems. What is more, the
classical counterpart—an entirely classical variational algorithm—also suffers
from the same problems. Ultimately, what matters is whether quantum processors
can accelerate parts of the computation relative to the best classical
algorithm.

##  3 Quantum computers: artificial many-body systems

This section explains how to investigate the many-body problems mentioned
above by building an artificial many-body system with a similar Hamiltonian
(aka analog quantum computers or quantum simulators) or, still starting from a
many-body system, but with individual control over the particles/degrees of
freedom, by building a gate-based (digital) quantum computer (subsection 3.1).
We then explain the basic building blocks and rules of ideal quantum
computers.

###  3.1 From analog to digital

Quantum computers are essentially synthetic many-body systems whose state can
be manipulated according to some predefined plan—aka a quantum program—and
measured to learn something.

Depending on the level of control of this quantum system, one speaks of
quantum simulators (or analog quantum computers) or quantum computers (or
digital quantum computers). While (analog) simulators offer only a limited and
specific set of controls, (digital) computers offer controls—usually called
gates—that are universal. This universality allows them to reach, in
principle, any state of the Hilbert space by performing any unitary operation.

####  3.1.1 Analog quantum computers (aka quantum simulators)

The term analog quantum computer refers to any synthetic many-body system with
a certain amount of control over its degrees of freedom. Each analog computer
is characterized (in the absence of defects) by a many-body Hamiltonian $H(t)$
whose time-dependence is "programmed" more or less at will, depending on the
experimental constraints. This Hamiltonian is usually chosen as close as
possible to the "real-life" Hamiltonians introduced in the previous section.
Thus, by measuring the properties of the analog simulator, one hopes to get
insights into the physics of real-life systems. We highlight below some
illustrations of physical systems used as analog simulators.

Ultracold atoms can be described as implementing a Fermi- or Bose-Hubbard
model, depending on the atomic isotopes used Bloch2008 . For instance, the
Bose-Hubbard model reads:

$H(t)=\frac{U(t)}{2}\sum_{i}n_{i}(n_{i}-1)-J(t)\sum_{\langle
ij\rangle}b_{i}^{\dagger}b_{j}.$ (3.1)

Here, the creation and annihilation operators $b^{\dagger}_{i}$ and $b_{i}$
create and annihilate (bosonic) atoms in orbitals $\phi_{i}(r)$ and
$n_{i}=b^{\dagger}_{i}b_{i}$; $J(t)$ is the tunneling between two neighboring
"sites" $\langle ij\rangle$ of the optical lattice, and $U(t)$ is the on-site
repulsion between two atoms. Both $U$ and $J$ can be temporally modulated by
changing the amplitudes of the lasers creating the lattice. In addition, the
interaction $U(t)$ can also be tuned by changing the background magnetic field
using a phenomenon known as the Feshbach resonance Bloch2008 .

Spin qubits, which are essentially electrons trapped in quantum dots, can also
be described by a Fermi-Hubbard model or, when neglecting charge fluctuations,
by a Heisenberg model Loss1998 :

$H(t)=J_{\mathrm{ex}}(t)\sum_{\langle
ij\rangle}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)+\sum_{i}H_{\mathrm{loc}}^{(i)},$
(3.2)

with $(X_{i},Y_{i},Z_{i})$ denoting the Pauli matrices acting on the
$i^{\mathrm{th}}$ spin, and with

$\displaystyle
H_{\mathrm{loc}}^{(i)}=\frac{\omega_{0}^{(i)}-\delta_{i}(t)}{2}Z_{i}+\Omega_{i}(t)\cos\left(\omega_{c}^{(i)}t+\phi_{i}(t)\right)X_{i}.$
(3.3)

The exchange constant $J_{\mathrm{ex}}\sim 4J(t)^{2}/U$ can be turned on and
off via the tuning of the tunneling term $J(t)$ between two dots using a gate
voltage. The local term $H_{\mathrm{loc}}$ can, for instance, come from a
magnetic field with a static ($Z_{i}$ term) and a rotating ($X_{i}$ term)
component.

Depending on which atomic levels they target, platforms of Rydberg atoms (see,
e.g., Browaeys2020 ; Henriet2020a ) may implement an Ising Hamiltonian:

$\displaystyle H(t)$ $\displaystyle=\sum_{ij,i\neq
j}\frac{C}{|r_{i}-r_{j}|^{6}}n_{i}n_{j}$
$\displaystyle+\frac{\Omega(t)}{2}\sum_{i}X_{i}-\delta(t)\sum_{i}Z_{i},$ (3.4)

with $n_{i}=(1-Z_{i})/2$, or a $XY$ Hamiltonian:

$\displaystyle H(t)$ $\displaystyle=2\sum_{ij,i\neq
j}\frac{C}{|r_{i}-r_{j}|^{3}}\left(X_{i}X_{j}+Y_{i}Y_{j}\right)$
$\displaystyle+\Omega(t)\sum_{i}X_{i}-\frac{\delta(t)}{2}\sum_{i}Z_{i}.$ (3.5)

Superconducting qubits, which are usually thought of as (digital) computers,
can also be seen as analog computers realizing a Bose-Hubbard model (see Eq.
(3.1)). There, the creation and annihilation operators refer to bosonic
excitations relative to the charge and flux variables inside Josephson
junctions.

All these Hamiltonians are of many-body nature owing to the coupling terms
(first term of each equation). Thanks to their closeness to the many-body
Hamiltonians encountered when studying quantum matter (see previous section,
2), these Hamiltonians have long been used as proxies (or simulators) for the
many-body systems one wants to understand.

By nature, quantum simulators are very specific in that (i) they implement (or
"simulate") only one class of Hamiltonians, and (ii) they usually implement
partial (often only global) control of their degrees of freedom. This
limitation is both temporal and spatial: e.g., in Eq. (3.4), the coupling term
cannot be made time-dependent as it corresponds to a van der Waals interaction
that cannot be switched off or decreased, except by moving the atoms, a very
slow operation. Also, the second (Rabi) and third (detuning) terms can most
often not be controlled at an individual site level (i.e., $\Omega$ and
$\delta$ are the same for all atoms).

Analog platforms have advantages and drawbacks. Their primary advantage is
that the limited degree of control usually allows them to work with
significantly more degrees of freedom (atoms, spins, ions, junctions, and
other building blocks). The major drawback is that their scope is limited
according to the specific resource Hamiltonian at hand, whose terms might not
be tunable at will. Experimental platforms endowed with reasonably good levels
of temporal and spatial control earn the qualification of gate-based quantum
computers.

Figure 3.1: Analog computation (a) vs. digital computation (b). In analog
computation, one directly specifies the (analog) parameters (here $\Omega(t)$
and $\delta(t)$, see, e.g., (3.3), (3.4), or (3.5)). These controls are not
necessarily local (as in (3.4) or (3.5)). In digital computations, the user
discretely describes the sought-after evolution with quantum gates, which are
usually local, i.e., act only on a few qubits (lines in the diagram).
Internally, each gate is performed using an analog description.

####  3.1.2 Digital (aka gate-based) quantum computers

Digital, or gate-based quantum computers, refer to physical setups (i) whose
description can be narrowed to an assembly of interacting $d$-level quantum
systems and (ii) whose Hamiltonian can be controlled at a local level.

In this review, we focus on the qubit case $d=2$, but note that qudit
approaches, where $d>2$, were also developed Wang2020 ; Jiang2020 ;
Ringbauer2022 .

##### Two-level quantum systems: qubits.

Let us focus on criterion (i). For instance, among the systems cited above,
Rydberg atoms or spin qubits are already naturally described as two-level
systems: two atomic levels for Rydberg atoms and two spin levels for spin
qubits. In photonic platforms, the photon’s two polarizations can play the
role of the two levels. Superconducting platforms, which are naturally
described with bosonic variables, can be restricted to a two-level subspace by
tuning their parameters so that "leakage" out of the two lowest levels—called
computational subspace—is very improbable. The two levels of the computational
subspace are usually denoted as $\ket{0}$ and $\ket{1}$. Hence, the
wavefunction of a single two-level system, or qubit, is, in general, the
superposition:

$\ket{\psi}=a_{0}\ket{0}+a_{1}\ket{1},$ (3.6)

with $a_{i}\in\mathbb{C}$ and $|a_{0}|^{2}+|a_{1}|^{2}=1$. More generally, a
$n$-qubit wavefunction $\ket{\Psi}$ is the superposition of $2^{n}$
computational basis states $|00\dots 0\rangle,|00\dots 01\rangle,\dots|11\dots
1\rangle$. We see that all states can be written as
$|b_{n-1},\cdots,b_{0}\rangle=\bigotimes_{i=0}^{n-1}|b_{i}\rangle$ where
$|b_{i}\rangle=|0_{i}\rangle$ or $|1_{i}\rangle$ refers to the state of the
$i^{\mathrm{th}}$ qubit. Below, we will use the same convention as in Eq.
(3.4), and operators that act on this qubit will be labeled by $i$ like, for
instance, the Pauli operators $(X_{i},Y_{i},Z_{i})$. Each state
$|b_{n-1},\cdots,b_{0}\rangle$ can also be labeled by a single integer
$k=\sum_{i}b_{i}2^{i}$.

##### Manipulating qubits: gates.

Criterion (ii) ensures that one can reach any state of this Hilbert space
using operations called quantum gates. Mathematically, these gates are unitary
operations $U$ acting on the wavefunction $\ket{\Psi}$:
$\ket{\Psi^{\prime}}=U\ket{\Psi}$. Such operations are performed by letting
the system evolve under a given Hamiltonian. For instance, let us consider a
$n$-qubit system described by the (non-interacting) Hamiltonian

$H=\sum_{i=1}^{n}H_{\mathrm{loc}}^{(i)},$ (3.7)

with $H_{\mathrm{loc}}^{(i)}$ defined in Eq. (3.3). $\omega_{0}^{(i)}$ is the
$i^{\mathrm{th}}$ qubit’s frequency (energy difference between the two
levels), and $\omega_{c}^{(i)}$ is the drive frequency. $\delta_{i}(t)$,
$g_{i}(t)$ and $\phi_{i}(t)$ in Eq. (3.3) are controllable fields. If one
switches off all but the $i^{\mathrm{th}}$ qubit’s field, one goes to the
frame rotating at frequency $\omega_{0}^{(i)}$, one drives at resonance
(namely $\omega_{c}^{(i)}=\omega_{0}^{(i)}$), and one neglects terms
oscillating as $2\omega_{0}^{(i)}$ (so-called rotating wave approximation),
the Hamiltonian reads, up to terms acting on the other qubits, as

$H_{\omega_{0}^{(i)}}=-\frac{\delta_{i}(t)}{2}Z_{i}+\frac{\Omega_{i}(t)}{2}\left[\cos(\phi_{i})X_{i}+\sin(\phi_{i})Y_{i}\right].$
(3.8)

If we turn off the Rabi term $\Omega_{i}$, solving the Schrödinger equation
yields a wavefunction $\ket{\psi(t)}=R_{z}^{(i)}(\theta(t))\ket{\psi(0)}$ with
$R_{z}^{(i)}(\theta)\equiv e^{-i\frac{\theta}{2}Z_{i}}$ and
$\theta(t)=-\int_{0}^{t}\delta_{i}(\tau)d\tau$. Thus, we have operated a
rotation of angle $\theta(t)$ around axis $z$ for the $i$th qubit. Similarly,
if we turn off the "detuning" term $\delta_{i}$, we are going to effect a
rotation along axes $x$ ($\phi=0$) and $y$ ($\phi=\pi/2$).

Such a time evolution is illustrated in Fig. 3.2, using the standard Bloch
sphere representation Nielsen2010 . We show the evolution of a one-qubit state
under a Rabi drive $\Omega(t)$ and a detuning drive $\delta(t)$ (see Fig.
3.1), with $\phi=0$ (green trajectory). The area under the Rabi curve
$-\Omega(t)$ is chosen to effect a $-\pi/2$ rotation, but as a consequence of
the detuning being not strictly zero, a small $z$-axis rotation is effected in
addition to the $x$-rotation. We will explain what happens when the qubit is
affected by decoherence (red trajectory) in a later section (section 8).

Figure 3.2: Bloch sphere with the North and South pole corresponding to
$|0\rangle$ and $|1\rangle$, respectively. The other poles shown in the figure
are decomposed in terms of $|0\rangle$ and $|1\rangle$ as
$|+\rangle=\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)$ and
$|i+\rangle=\frac{1}{\sqrt{2}}\left(|0\rangle+i|1\rangle\right)$. The green
and red trajectories represent the evolution of a one-qubit state, starting
from $|0\rangle$, when subject to a $R_{x}(-\pi/2)$ rotation. Green curve:
noiseless (pure state) evolution. This evolution approximately realizes a
$R_{x}(-\pi/2)$ rotation (it would be exactly such a rotation if the detuning
drive $\delta(t)$ were rigorously 0). Red curve: noisy evolution under
dephasing and relaxation noise, see section 8. Table 1: Summary of some
standard single-qubit quantum gates. The Clifford group mentioned in section
2.3 is depicted in purple and can be generated using the H,
$S=P\left(\pi/2\right)$ and CNOT gates (see Table 2). The purple and green
hatches on the phase gate $P\left(\varphi\right)$ reflect that it is Clifford
only for $\varphi=\pm\pi/2$. Table adapted from Motta2022 . Table 2: Summary
of some common two-qubit quantum gates. We use the convention where the
uppermost qubit line corresponds to the least significant bit in the binary
representation of the computational state. This convention corresponds to the
so-called "little endian" convention used for bit/qubit storing. We will use
this convention systematically throughout the article for circuit
representation. The A matrix presented from a) to d) in the Table is equal to
$\big{(}\begin{smallmatrix}a_{00}&a_{01}\\\
a_{10}&a_{11}\end{smallmatrix}\big{)}$. Some circuits are adapted from
Seki2022 .

To produce entanglement, we need a Hamiltonian with interacting spins. For
instance, if we can switch on a term $JZ_{i}Z_{j}$ in Hamiltonian (3.7) (a
term similar to the van der Waals term in Eq. (3.4)), we can perform
operations of the type $e^{-i\frac{\theta}{2}Z_{i}Z_{j}}$. Such operations can
create entanglement between qubits $i$ and $j$.

##### Universal quantum gates

It turns out that, with the one-qubit rotations $R_{x}\left(\theta\right)$,
$R_{y}\left(\theta\right)$, and $R_{z}\left(\theta\right)$ presented in Table
1, together with, for instance, a two-qubit gate called “CNOT” presented in
Table 2, one can achieve any unitary operation $U$ acting on $n$-qubits as a
finite sequence of these gates (a result known as the Solovay-Kitaev theorem
Kitaev1997 ; Nielsen2010 ). Therefore, one calls this gate set a universal
gate set. The Clifford group (see section 2.3) can become universal if a $T$
gate is added (with $T=P\left(\pi/4\right)$, see Table 1). Below, we briefly
discuss the fundamentals of digital quantum computation; for more advanced
considerations, we refer to different textbooks Nielsen2010 ; Mcmahon2007 ;
Kaye2006 ; deLima2019 .

##### Quantum circuits

A standard digital quantum computation is a sequence of simple manipulations
of the system’s Hamiltonian; each described as a quantum gate. The computation
usually starts from an initial state corresponding to all qubits in state
$\ket{0}$. In other words, the "quantum register" is in state
$|0\rangle^{\otimes n}$. One then applies gates $U_{1},U_{2},\dots,U_{m}$.
This sequence of gates is usually represented as a so-called quantum circuit,
where each line stands for a qubit (time flowing from left to right) and each
gate is pictured by a symbol that acts only on a subset of these lines. Table
1 and 2 respectively give examples of the most common gates acting on one or
two qubits.

##### Quantum measurements

After applying the gates, the register is in its final state
$\ket{\Psi}=U_{m}U_{m-1}\cdots U_{1}|0\rangle^{\otimes n}$ and one can measure
some observable. In most platforms, one can only measure the observable
$Z_{i}$ (or a tensor product $Z_{i_{1}}\otimes\cdots\otimes Z_{i_{k}}$ if one
"measures" $k$ qubits). One can translate a measurement in the $X$ or $Y$
basis into a measurement in the $Z$ basis using the insertion of one or two
gates before the measurement, as pictured in Table 3.

The outcome of a measurement of $Z_{i_{1}}\otimes\cdots\otimes Z_{i_{k}}$ is a
bitstring $b_{i_{1}}\dots b_{i_{k}}$, with $b_{i}\in\\{0,1\\}$. This bitstring
is obtained with a probability given by Born’s rule,

$p(b_{i_{1}}\dots b_{i_{k}})=\langle\Psi|P_{b_{i_{1}}\dots
b_{i_{k}}}|\Psi\rangle,$ (3.9)

with $P_{b_{i_{1}}\dots b_{i_{k}}}=|b_{i_{1}}\rangle\langle
b_{i_{1}}|\otimes\cdots\otimes|b_{i_{k}}\rangle\langle b_{i_{k}}|$ (we do not
explicitly write identities for qubits that are not being measured). The
measurement projects the register to the state $P_{b_{i_{1}}\dots
b_{i_{k}}}|\Psi\rangle/\sqrt{p(b_{i_{1}}\dots b_{i_{k}})}$, so that if one
wants to measure another observable that does not commute with
$Z_{i_{1}}\otimes\cdots\otimes Z_{i_{k}}$, one needs to rerun the circuit.

The estimation of the expectation value of an observable (hermitian operator),
like $\langle O\rangle=\bra{\Psi}{O}\ket{\Psi}$, is typically done by
measuring the given observable a number $n_{\mathrm{shots}}$ of times,
resulting in values $\\{o_{k}\\}_{k=1,n_{\mathrm{shots}}}$ that can be
averaged to yield the estimator

$\overline{O}=\frac{1}{n_{\mathrm{shots}}}\sum_{k=1}^{n_{\mathrm{shots}}}o_{k}.$
(3.10)

In the limit $n_{\mathrm{shots}}\rightarrow\infty$, this estimate converges to
$\langle O\rangle$. Due to the central limit theorem,
$\mathcal{O}(1/\varepsilon^{2})$ samples are needed to reach an accuracy
$\varepsilon$.

Quantum circuits can also be used to compute the average value of any unitary
operator, as shown in Table 4. This table also shows how to measure time-
dependent correlation functions of the form $\langle
P_{l}(t)P_{k}(t^{\prime})\rangle$ (using fermion-spin transforms, this type of
circuit can be used to compute fermionic correlation functions).

Table 3: Conversion of measurements in $X$ or $Y$ bases into measurements in
the $Z$ basis. Table 4: Here we present some standard interferometry circuits
which allow computing overlaps by making measurements on an ancillary qubit.
The first circuit evaluates quantities of the form $\langle\Psi|U|\Psi\rangle$
where $UU^{\dagger}=\mathbb{I}$, e.g., the overlap between a state and its
time-evolved counterpart. The second circuit enables to retrieve two-time
correlators of the form $\langle\Psi|A(t)B(t^{\prime})|\Psi\rangle$, with
$A(t)=U^{\dagger}(t,0)P_{k}U(t,0)$ and
$B(t^{\prime})=U^{\dagger}(t^{\prime},0)P_{l}U(t^{\prime},0)$ (where $P_{i}$
denotes a Pauli operator). On the right, $\langle{O}_{0}\rangle_{\cal C}$
indicates that the ancillary qubit (labeled by convention as the
"$0^{\mathrm{th}}$" qubit) is to be measured in the $O\in\\{X,Y,Z\\}$ basis
after the execution of the circuit $\mathcal{C}$ drawn on the left.

##### DiVincenzo criteria

The principles introduced above have been gathered in a list of five criteria
known as the DiVincenzo criteria divincenzo_physical_2000 :

  1. 1.

The ability to work with a scalable number of two-level systems (qubits)
without "leakage" out of the computational subspace.

  2. 2.

The capacity to initialize and reset qubits in a reliable (usually fast
enough) fashion.

  3. 3.

A long coherence time (compared to the typical time scales of gates,
measurements, and resets, i.e., compared to the "clock time" of the
processor).

  4. 4.

A universal set of gates (with the possibility to parallelize operations on
disjoint sets of qubits).

  5. 5.

Reliable qubit-wise measurements.

Current quantum processors are strongly impacted by decoherence effects, as
will be explicited in section 8. As a result, they do not meet all five
criteria. However, they help develop and test algorithms in real-life
conditions.

Before reviewing these algorithms, we describe how many-body problems can be
translated to forms amenable to quantum computing.

##  4 Mapping a many-body problem to a quantum computer

This section explains how to go from the many-body problem at hand to the one
that the quantum computer models. In particular, if we focus on digital (gate-
based) quantum platforms, such devices usually have constraints: (i) they have
qubits (two-level systems), not fermions/bosons; (section 4.1), (ii) they have
a limited number of qubits (section 4.2), and (iii) they have a limited
coherence (section 8.1). Thus, one needs to transform, reduce or/and map the
original problem so that the quantum computer can give insights into the
properties of the original many-body problem.

###  4.1 From fermions to qubits

The treatment of fermions on a quantum computer can be made starting from the
first or from the second quantization McArdle2018a . While we mainly focus in
this review on the second option, the two strategies, i.e., first or second
quantization as a starting point, are briefly discussed below.

####  4.1.1 Encoding fermions as qubits in second quantization

To perform any Hamiltonian simulation written in second quantized form on a
quantum computer, one must map the fermion Hamiltonian to a spin Hamiltonian.
This mapping is not unique. The most standard mapping techniques are the
Jordan-Wigner (JW) transformation JordanWigner1928 , the Bravyi-Kitaev (BK)
Bravyi2002  or the parity mapping Seeley2012  (for a comprehensive discussion
see Bauer2016 ; Fano2019 ; McArdle2018a ). We illustrate the JW case that is
often retained for many-body applications due to its relative simplicity.

Let us consider a set of fermions associated with the creation/annihilation
operators $(a^{\dagger}_{p},a_{p})$ where $p$ labels a complete basis of
single-particle states $\phi_{p}(r)$. These operators act on the many-body
vacuum by changing the occupation of orbital $p$. In the JW fermion-to-qubit
mapping, the occupation (resp. vacancy) of a state is usually encoded as the
state $|1\rangle_{p}$ (resp. $|0\rangle_{p}$). Then, the operator
$Q^{+}_{p}=|1\rangle_{p}\langle 0|_{p}=\frac{1}{2}\left(X_{p}-iY_{p}\right)$
and its hermitian conjugate $Q^{-}_{p}$ can be seen as the qubit equivalents
of the creation/annihilation operators. The difficulty is that these sets of
operators commute between each other for different qubits while they should
anticommute for fermions. One solution to this issue is to choose a specific
ordering for the one-to-one correspondence between the single-particle state
and the qubits and use the following prescription:

$\displaystyle a^{\dagger}_{p}$ $\displaystyle\longleftrightarrow$
$\displaystyle\bigotimes_{k=1}^{p-1}Z_{k}\otimes
Q_{p}^{+},\leavevmode\nobreak\
a_{p}\longleftrightarrow\bigotimes_{k=1}^{p-1}Z_{k}\otimes Q^{-}_{p}.$ (4.1)

In this transformation, the fermionic sign (which comes from the
anticommutation rules of fermions) is kept track of via the string
$\bigotimes_{k=1}^{p-1}Z_{k}$ of Pauli-$Z$ operators. At the circuit level,
this means that operations that are one-body at a fermionic level (like
$\exp(-it\\{c_{0}^{\dagger}c_{1}+\mathrm{h.c}\\})$) might become a multi-qubit
operation. For instance, $a^{\dagger}_{3}a_{1}$ leads to a term
$Q^{\dagger}_{3}Z_{2}Q_{1}$ that acts on the three qubits $(1,2,3)$.

The fermion–to–qubit mapping corresponds to the mapping of a set of
$\tilde{n}$ single-particle orbitals obeying fermionic anticommutation rules
to a set of $n$ qubits. The mapping is not unique, and depending on which
choice is made the number of qubits $n$ may not be identical to the number of
spin-orbitals $\tilde{n}$. Furthermore, the locality $d$ (a $d$-local
Hamiltonian can be expressed as the sum of Hamiltonian terms acting upon at
most $d$ qubits) is usually not conserved upon encoding. For instance, in the
Jordan-Wigner mapping, $\tilde{n}=n$ and $\tilde{d}=O(\tilde{n})$, while
another transform called the Bravyi-Kitaev transformation that also has
$\tilde{n}=n$ achieves a better locality, namely
$\tilde{d}=O(\log_{2}{\tilde{n}})$. Note that in general it is advantageous to
use the Bravyi-Kitaev encoding compared to the Jordan-Wigner encoding in terms
of circuit depth Tranter2018 . Despite this, the Jordan-Wigner mapping is
often used in practice due to its relative simplicity. Another example is the
so-called superfast fermionic encoding Setia2018 , which requires
$\tilde{n}=O(md)$ and achieves $\tilde{d}=O(d)$. Here $m$ represents the
number of Fermi modes, such that each mode can be either empty or occupied by
a fermionic particle.

Note that we have implicitly referred to $SU(2)$ fermions only, but $SU(N)$
fermions can also be handled. For instance, the JW mapping scheme presented
above readily generalizes to $SU(N)$ fermions Consiglio_2022 .

####  4.1.2 Encoding fermions in first quantization

One major drawback of the encoding based on second quantization is that the
many-body problem is encoding directly in the entire Fock space, whose size
grows exponentially with the number of considered orbitals $n$ as $2^{n}$. On
the other hand, we know that a problem of $A$ interacting particles described
on $n$ single-particle states requires at most $C^{A}_{n}$ many-body states,
just by invoking the particle number conservation. Additional symmetries can
further reduce the size of the relevant Hilbert space. Said differently,
second quantization applied to many-body problems leads to a significant dark
sector that might ultimately reduce quantum advantage. First quantization,
which was the original formulation of quantum mechanics, has been explored as
a possible solution to this problem. An exhaustive discussion of the
particles–to–qubit encoding using first quantization is out of the scope of
the present review and we only give some brief guidelines. For more details,
see for instance Ref. McArdle2018a  and the historical overview given in Fig.
5 of Ref. Cao2019 .

First quantization directly relies on the wave-function representation of the
problem of interacting particles without invoking the machinery associated
with Fock space and creation/annihilation operators. Some anticipated
advantages of this formulation are underlined below. We start from a
Hamiltonian describing a set of particles $i=1,\dots,N$ and written in first
quantization as:

$\displaystyle H=\sum_{i=1}^{N}h(i)+\sum_{i<j=1}^{N}V(i,j),$ (4.2)

where $h(i)$ is a one-body Hamiltonian acting on the $i$th particle, and
$V(i,j)$ denotes the interaction between particles $i$ and $j$. Two main
strategies have been used to encode the problem directly starting from
Hamiltonian (4.2):

  1. 1.

Grid-based method: Let us assume that the particle is described by a certain
set of coordinates $\\{\bf q\\}$, for instance, its position and spin, i.e.,
${\bf q}=(x,y,z,\sigma)$. Then, the many-body wave-function, denoted by
$|\Psi_{N}\rangle$, can be decomposed as

$\displaystyle|\Psi_{N}\rangle$ $\displaystyle=$ $\displaystyle\int
d^{3N}\Omega\Psi_{N}({\bf q}_{1},\cdots,{\bf q}_{N})|{\bf q}_{1},\cdots{\bf
q}_{N}\rangle.$ (4.3)

Here $\int d^{3N}\Omega$ denotes the integral over the $3N$ dimensional space
and the sum over the $2N$ spin components. One advantage of this
representation is that the state $|{\bf q}_{1},\cdots{\bf q}_{N}\rangle$ can
already contain the antisymmetry property Blaizot86 . Continuous variables
like position or momentum are generally not bounded and application requires
specific discretization schemes on a grid. The possibility to describe quantum
systems on a grid using quantum computers has been explored already in
pioneering articles Zalka1998 ; Lloyd1996 . The encoding can be made by
assuming that each $(x,y,z)$ component is discretized on a mesh with $L$
points. As discussed in detail in Ref. McArdle2018a , the description of a
wave-function given by (4.3) requires storing $L^{3N}\times 2^{N}$ amplitudes
and, therefore, becomes rapidly prohibitive on a classical computer. Brute-
force few-body problems treated on a classical computer are usually restricted
to small number of particles $N$. Assuming that $L=2^{m-1}$, such amplitudes
can be stored on $(3m+1)N$ qubits. Such encoding becomes more compact and
automatically treats the problem in the proper space of wave-function having
$N$ particles. Using grid-based methods, the different operators entering in
(4.2) can be encoded in the basis by rewriting them as:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle h(i)=\sum_{{\bf p},{\bf
q}}\langle{\bf p}|h|{\bf q}\rangle\left(|{\bf p}\rangle\langle{\bf
q}|\right)_{i}\\\ \displaystyle V(i,j)=\sum_{{\bf p},{\bf q},{\bf r},{\bf
s}}\langle{\bf p},{\bf r}|V|{\bf q},{\bf s}\rangle\left(|{\bf
p}\rangle\langle{\bf q}|\right)_{i}\left(|{\bf r}\rangle\langle{\bf
s}|\right)_{j},\end{array}\right.$ (4.6)

where we use the notation Babbush2019

$\displaystyle\left(|{\bf p}\rangle\langle{\bf q}|\right)_{i}\equiv|{\bf
p}\rangle_{i}\langle{\bf q}|_{i}\bigotimes_{k\neq i}I_{k}.$

With these expressions, we see that each term of the one- and two-body
components will only couple one state of the computational basis to a single
state of the same basis. Following this discretization technique, a number of
quantum algorithms eventually taking advantage of the sparsity of the
Hamiltonian have been proposed Kivlichan2017 ; Babbush2018 ; Babbush2019 ;
Su2021  (see also McArdle2018a ; Cao2019 ).

  2. 2.

Basis set method or restricted CI method: the direct solution of a many-body
problem directly on a mesh is rather challenging and requires a good control
of the precision both on the quantum algorithm and on the discretization
scheme itself Babbush2019 . As discussed in the introduction and further
elaborated in section 4.2, classical algorithms based on CI techniques are
built using a set of $M_{a}$ "active" single-particle states
$\\{|\phi_{l}\rangle_{l=1,M_{a}}\\}$ where a number $M_{p}$ of particles can
be distributed to form Slater determinants. If $M_{p}$ identifies with the
total number of particles $N$ and a sufficient number of single-particle
states are considered the CI calculation is said unrestricted; otherwise, it
is referred to as restricted CI. One can then write the Hamiltonian components
in a similar way as in Eq. (4.6) replacing the states $|{\bf p}\rangle$ by the
new single-particles states. The number of many-body states to consider is
$C^{M_{p}}_{M_{a}}$ and the Hamiltonian is sparse in this basis. As proposed
in Ref. Toloui2013 , one can label a given state as
$|m_{1},\cdots,m_{N}\rangle$ with $m_{l}=0,M_{a}-1$ denoting the orbital that
the particle $l$ is occupying. Let us assume simply that $M_{a}=2^{n_{a}}$,
then, a set of integers $m_{l}=0,M_{a}-1$ can be encoded on $n_{a}$ qubits
Abrams1997 . Accordingly, the set of states $|m_{1},\cdots,m_{N}\rangle$ can
be encoded on $N\log_{2}(M_{a})$ qubits. Illustrations of quantum algorithms
using this technique can be found in Toloui2013 ; McArdle2018a .

###  4.2 Reducing the number of degrees of freedom

NISQ devices come with a limited number of qubits and limited coherence. These
constraints limit the number of degrees of freedom (typically orbitals) of the
system one wants to study. Nevertheless, many-body condensed matter, quantum
chemistry, or nuclear problems typically comprise tens or hundreds of
orbitals. Directly tackling these large model spaces with a quantum processor
appears unfeasible, if not ill-advised. Indeed, over the last century,
numerous classical many-body methods have been devised to reduce the number of
truly correlated—sometimes called "ultra quantum"—degrees of freedom. These
classical methods then rely on advanced algorithms to solve the "reduced
model". Despite its reduced complexity, this model is usually hard to tackle
in some physically relevant regimes due to its strongly-correlated character.
This is where quantum coprocessors could be used to extend the power of
classical algorithms. We highlight below a few classical-inspired methods that
were used to reduce the number of qubits for many-body systems treated on
quantum computers. These methods are illustrated in Fig. 4.1.

Figure 4.1: Embedding methods (a) and active-space selection (b). In embedding
methods, an extended (lattice) model is self-consistently mapped to a local
(impurity or embedded) model. In active-space methods, a subset of orbitals
(usually partially filled ones) is selected to construct the active space
Hamiltonian, while the other orbitals are treated at the mean-field (Hartree-
Fock) level. Due to their smaller size, the embedded or active-space models
are better suited for a solution with today’s quantum computers.

####  4.2.1 Embedding methods

Typical condensed-matter problems are formulated on a lattice of atomic sites
(e.g., the Hubbard model introduced in Eq. (2.2)). The number of sites needed
to observe collective phenomena like phase transitions typically exceeds the
capacity of exact diagonalization or Monte-Carlo methods. Classical methods
collected under the term embedding methods have been developed to overcome
this limitation. They draw inspiration from mean-field methods in that they
self-consistently map the original, extended problem onto a smaller, more
local many-body problem (sometimes called fragment) "embedded" in a (usually)
non-interacting environment (also called bath). One can then leverage the fact
that this embedded problem has fewer correlated degrees of freedom to tackle
it with classical or, if need be, quantum methods Bauer2016 . An illustration
for these methods is provided in Figure 4.1 (a).

Examples of embedding methods include, but are not limited to, Dynamical Mean
Field Theory (DMFT)Georges1996 , the Gutzwiller or Rotationally-Invariant
Slave Boson (RISB)Lechermann2007  method, and the Density-Matrix Embedding
Theory (DMET)Knizia2012  method.

Generically, the embedded problem has the form:

$\displaystyle H$
$\displaystyle=\sum_{i=1}^{N_{c}}Un^{c}_{i\uparrow}n^{c}_{i\downarrow}-\mu\sum_{i=1}^{N_{c}}\sum_{\sigma}n^{c}_{i\sigma}$
(4.7)
$\displaystyle+\sum_{p=1}^{N_{b}}\sum_{i=1}^{N_{c}}\sum_{\sigma}\left(V_{p}c^{\dagger}_{i\sigma}a_{p\sigma}+\mathrm{h.c}\right)+\sum_{p=1}^{N_{b}}\sum_{\sigma}\varepsilon_{p}a^{\dagger}_{p\sigma}a_{p\sigma},$
(4.8)

where $a^{\dagger}_{p\sigma}$ creates electrons in the bath (of size $N_{b}$),
while $c^{\dagger}_{i\sigma}$ creates electrons in the correlated orbitals
($n_{i\sigma}^{c}=c^{\dagger}_{i\sigma}c_{i\sigma}$). Compared to the Hubbard
model, Eq. (4.8) has fewer ($N_{c}$) interacting sites (the $N_{b}$ bath sites
are uncorrelated). However, it is still a complicated many-body problem.
Typically, $N_{c}$ is adapted to the spatial resolution one wants. For regimes
with considerable correlation lengths, it can exceed the reach of advanced
Monte-Carlo methods.

The embedding methods mentioned above differ by the number of bath sites, the
observables that need to be computed (generally, Green’s functions on the
impurity or reduced density matrices), and the way the self-consistent
parameters ($\varepsilon_{p}$ and $V_{p}$) are updated. For instance, within
RISB and DMET, $N_{b}=N_{c}$; this results in the embedded model being much
smaller than the original Hubbard model.

Recent works have used quantum processors to tackle the embedded model within
an embedding method Rungger2019 ; Keen2019 ; Tilly2021 ; Yao2020 ;
besserve_ayral_2022 ; Backes2023 . They are limited so far to small sizes
($N_{c}\leq 2$) due to NISQ limitations. Until now, classical methods still
outperform quantum methods in solving these problems.

####  4.2.2 Active space methods

In quantum chemistry, like in condensed matter or nuclear physics problems,
the number of degrees of freedom can be reduced to the genuinely complicated
degrees of freedom. These are usually called active orbitals. Instead of
handling all orbitals at the same level of theory, orbitals are divided into
active ones—which require an advanced many-body method—and inactive ones—for
which mean-field (Hartree-Fock) methods will be sufficient. The active space
selection can be based on the occupation level of molecular orbitals (the
orbitals resulting from a Hartree-Fock optimization). Empty and occupied
orbitals (with occupation numbers close to $0$ or $1$, respectively) are
inactive, while partially-filled orbitals are considered active. The active
space size is adjusted according to the sought-after accuracy, available
computational capacity, or both. The so-obtained active space Hamiltonian has
the same form as the original Hamiltonian, which is given by Eq. (2.3).
However, it usually has a much smaller number of orbitals: $N$ is reduced to
$N_{a}$. Then, with a classical computer, one can tackle this reduced problem
with advanced methods like FCI (if $N_{a}$ is very small) or CC otherwise. The
simplification of a many-body problem’s solution by selecting some active
space without degrading the quality of an approach is a complex and delicate
issue where one can take advantage of quantum information concepts to make the
selection of relevant states Stein2016 ; Stein2019 .

With a quantum coprocessor, the reduction from $N$ to $N_{a}$ orbitals
directly translates, via fermion-spin transforms (see section 4.1, to a
reduced number of required qubits (namely $N_{a}$). An illustration for these
methods is provided in Figure 4.1 (b).

Many recent works use this active space selection to reduce the number of
required qubits, see, e.g., McCaskey2019 , or to explore the resource
requirements on future quantum computers Reiher2016 ; Li2018a .

##  5 Ideal algorithms

Here, we discuss some textbook methods to simulate a quantum system and solve
the eigenvalue problems on a quantum computer Nielsen2010 .

###  5.1 Quantum Phase Estimation for the eigenvalue problem

##### Description

Quantum Phase Estimation (QPE), also called Phase Estimation Algorithm, is a
generic method to shed light on the spectrum of a unitary operator with a
quantum computer Nielsen2010 ; deLima2019 . It is already well documented, and
we only give here the key ingredients of the approach.

Suppose that the operator $U$ has a set of eigenvalues $\\{e^{2\pi
i\phi^{\alpha}}\\}$. We assume that for all $\alpha$, we have
$0\leq\phi^{\alpha}<1$ and denote by $\\{|\alpha\rangle\\}$ the corresponding
eigenvectors. The system’s initial state can be decomposed as:

$\displaystyle|\Psi\rangle$ $\displaystyle=$
$\displaystyle\sum_{\alpha}c_{\alpha}|\alpha\rangle.$ (5.1)

Figure 5.1: Illustration of the QPE circuit used to get the eigenvalues of an
arbitrary unitary operator $U$. The QPE circuit uses the inverse Quantum
Fourier Transform (QFT-1) Nielsen2010 . The QPE requires $n_{a}$ additional
ancillary qubits. The precision of the method will directly depend on $n_{a}$.

The QPE method consists of the following schematic sequence:

$\displaystyle|\Psi\rangle\xrightarrow[]{\text{QPE}}\sum_{\alpha}c_{\alpha}|\alpha\rangle\otimes|\widetilde{\phi}^{\alpha}\rangle\xrightarrow[]{\text{Measure}}|\alpha\rangle\otimes|\widetilde{\phi}^{\alpha}\rangle,$
(5.2)

where the bitstring
$\widetilde{\phi}^{\alpha}=\phi^{\alpha}_{n_{a}-1}\dots\phi^{\alpha}_{0}$
($\phi^{\alpha}_{i}\in\\{0,1\\}$ being the result of the measurement of the
ancillary qubit $i$) encodes an estimation of the phase $\phi^{\alpha}$ as the
binary fraction
$0.\phi^{\alpha}_{0}\cdots\phi^{\alpha}_{n_{a}-1}\equiv\sum\limits_{j=0}^{n_{a}-1}\frac{\phi^{\alpha}_{j}}{2^{j+1}}$.
We will use the notation $0.\widetilde{\phi}^{\alpha}$ for this number,
although the order of the bits is reversed. A schematic view of the QPE
circuit is shown in Fig. 5.1.

The effect of QPE is twofold: the initial state is projected into one
eigenstate (or a set of degenerate states) having non-vanishing overlap with
the initial state, and the associated eigenvalue is retrieved with a precision
$|{\phi}^{\alpha}-0.\widetilde{\phi}^{\alpha}|\leq 1/2^{n_{a}}$, which
improves exponentially with $n_{a}$. Note that the projection is only
approximate unless the binary fraction of ${\phi}^{\alpha}$ is finite and the
number $n_{a}$ is sufficient to have
${\phi}^{\alpha}=0.\widetilde{\phi}^{\alpha}$. This projection effect was used
recently in a many-body system to restore broken symmetries (see Lacroix2020 ;
RuizGuzman2022 ).

For many-body problems, QPE can be seen as a gold standard to solve the
eigenvalue problem for the Hamiltonian $H$ in a large Hilbert space. In this
case, the operator $U$ can be chosen as the propagator itself, with

$\displaystyle U(\tau)=e^{-2\pi i\tau(H-E_{0})},$ (5.3)

where $\tau$ and $E_{0}$ are parameters chosen to map the spectrum of $H$ into
elements of $[0,1[$. The circuit to prepare the unitary $U$ starting from a
given $H$ is usually obtained via trotterization, a method summarized in
section 5.2.

Figure 5.2: Illustration of the QPE algorithm applied to the Hartree-Fock
state $|\psi\rangle=|00001111\rangle$ of a Pairing Hamiltonian, Eq. (2.5),
with one body energies at levels j, $\varepsilon_{j=0,1,\dots,n_{q}-1}=j\Delta
e$, $g=0.5\Delta e$ and $\Delta e=1$ RuizGuzman2022 . p is the probability of
measuring the energy value E in the ancillary register. The ancillary qubits
$n_{a}$ used for the QPE were 4 (a), 6 (b), and 8 (c). The vertical green and
black lines correspond to the ground and excited energies of the Hamiltonian,
respectively. The figure has been adapted from RuizGuzman2022 . Note that
panel (c) is shown in linear-log scale.

A key feature of the QPE algorithm is that its complexity for finding the
phase $\phi^{\alpha}$ associated with the eigenstate $|\alpha\rangle$ with
additive error $\epsilon$ scales as:

$\mathcal{O}\left(\frac{{\text{poly}(n_{q})}}{\epsilon S^{2}}\right).$ (5.4)

Here, $S$ denotes the overlap between the initial state $|\Psi\rangle$ and
$|\alpha\rangle$, i.e., $S=|\langle\alpha|\Psi\rangle|$. For large systems,
$S$ is generically small because it becomes exponentially difficult to attain
a high overlap value due to the exponential expansion of the Hilbert space, a
phenomenon known as the orthogonality catastrophe Lee2022 ; Louvet2023 . It is
important to note that the cost of preparing the initial state $|\Psi\rangle$
has not been included in this analysis, as we assume that a cost-efficient
approach is used: the term $1/S^{2}$ simply reflects the scaling of the number
of measurements required to identify $\phi^{\alpha}$. In the general case,
however, preparing an arbitrary initial state has exponential scaling with
respect to $n_{q}$.

Since $\epsilon\sim 1/2^{n_{a}}$, the term $1/\epsilon$ is linked both to the
number of ancilla qubits required for a specified precision and to the circuit
depth, given that there are $n_{a}$ controlled-$U$ gates and $O(n_{a}^{2})$
operations are needed for the inverse Quantum Fourier Transform. Typically,
the $U$ gates exhibit polynomial scaling with respect to the number of system
qubits $n_{q}$, i.e., $O(\text{poly}(n_{q}))$, and for electronic structure
problems, this scaling is commonly $O(n_{q}^{5})$. Furthermore, as highlighted
in a later section (5.2) and due to the no fast-forwarding theorem berry_2006
, the complexity of the $U$ gates generally scales linearly with the evolution
time $t$. This implies that the application of a sequence of
controlled-$U^{2^{k}}$ gates will usually result in exponential complexity,
given that the total evolution time needed to execute the algorithm is
$\tau(2^{n_{a}}-1)$.

These scalings illustrate the advantages and shortcomings of the “perfect” QPE
(i.e., performed on a fault-tolerant quantum computer). On the flip side, it
means that QPE requires an input state that reasonably overlaps the actual
eigenstate. One could, for instance, resort to adiabatic state preparation
Farhi2000  to rotate some simple initial state into a state exhibiting a
significant overlap with the eigenstate of interest. This requirement over the
overlap becomes problematic in high dimension for many-body problems due to
the orthogonality catastrophe referred to above. Note that this initial state
problem can be cleverly handled in some cases. A case in point is Shor’s
factoring algorithm Shor1994 , which uses QPE as the main ingredient after a
clever state preparation, providing an exponential speedup over the classical
version of the factoring algorithm.

On the bright side, the $1/\epsilon$ scaling is much better than the
$1/\epsilon^{2}$ typical of classical Monte-Carlo methods. This advantage is
used in "quantum Monte-Carlo" methods. These techniques use QPE as a critical
building block within another algorithm called quantum amplitude estimation
Brassard2000 . We point out that these techniques should not be confused with
the many-body quantum Monte-Carlo methods, which refer to purely classical
methods to solve quantum many-body problems. We also note in passing that
recent proposals have been made to hybridize "classical" quantum Monte-Carlo
methods such as Auxiliary-Field Quantum Monte-Carlo (AFQMC) or Full-
Configuration-Interaction-Quantum-Monte-Carlo (FCIQMC) with quantum algorithms
Huggins2021 ; Mazzola2022 ; Zhang2022a .

On NISQ processors, QPE is hardly applicable owing to the number of operations
required (whether the previously mentioned $2^{n_{a}}$ depth or those of the
QFT). Estimates for chemical problems yield tremendous numbers Wecker2014a ,
not compatible with the number of qubits and error rates of current and near-
term machines.

Today, intensive efforts are being made to provide less costly methods for
Hamiltonian eigenvalue problems. Some of them will be further reviewed in
section 6.

###  5.2 Trotterization

##### Description

Trotterization is a technique to implement a time evolution $U=e^{-iHt}$ on
digital quantum hardware Lloyd1996 , i.e., as a sequence of few-qubit gates.

Most many-body Hamiltonians are $k$-local, meaning they can be decomposed as a
sum of terms acting on at most $k$ qubits:

$H=\sum_{j=1}^{m}\lambda_{j}P_{j},$ (5.5)

with $P_{j}$ a product of at most $k$ Pauli operators.

Trotterization approximates the exponential of a sum $e^{-iHt}$ as the product
of the individual terms. The so-called first-order Trotter-Suzuki formula
Lloyd1996  reads

$\displaystyle e^{-iHt}$
$\displaystyle=\left(\prod_{j=1}^{m}e^{-i\lambda_{j}P_{j}\frac{t}{n_{t}}}\right)^{n_{t}}+O\left(\frac{m^{2}t^{2}}{n_{t}}\right)$
(5.6)

where $n_{t}$ is the number of Trotter steps. The rationale behind formula
(5.6) is that the whole Hamiltonian evolution is carried out in the form of
repeated sequences of step-wise evolutions
$e^{-i\lambda_{j}P_{j}\frac{t}{n_{t}}}$. Each such evolution can be simplified
as a sequence of one- and two-qubit gates.

Notably, the number $n_{t}$ of Trotter steps must be increased when $t$
increases: the circuit structure is not fixed with $t$; instead, its
complexity grows linearly with $t$. In some cases, a sublinear scaling can be
proposed, but this is not generally the case because of the no fast-forwarding
theorem berry_2006 . In the QPE described in the previous section, this
implies that a unitary operator of the form $U^{2^{n}}=\exp(-iH2^{n}t)$ has
depth $O(2^{n})$, explaining the exponential scalings of QPE.

##### Beyond standard trotterization

The possibility of reducing the significant scaling associated with the
Trotter-Suzuki methods is an active field of research. Several alternative
methods have been proposed: the Variational Fast-Forwarding Cirstoiu2020 ;
Gibbs2021 , Incremental Structure Learning Jaderberg2020 , the Adaptive
Product Formula Zhang2020 , and the Variational Time Dependent Phase
Estimation Klymko2021 , to quote some of them.

##  6 NISQ Algorithms

Quantum algorithms like the QPE presented above require, in general, a large
number of qubits or gates, or both. Because of the limitations of current
quantum platforms in qubit and gate counts, these algorithms cannot be used in
the presence of noise. Specific algorithms have been designed to circumvent
those limitations and study the applicability of today’s quantum processors to
concrete problems. These algorithms are designed, for instance, to reduce the
circuit depth by allocating only a specific task to the quantum computer or,
as when using variational methods, to allow better control of these states.

###  6.1 Variational algorithms

Variational methods are standard tools for many-body physicists using
classical computers Blaizot86 . In recent years, they have emerged as a tool
of choice for applications on quantum platforms, and are an important part of
the broader class of hybrid quantum-classical methods Bharti2021 ; Endo2021 .

####  6.1.1 Variational Quantum Eigensolver (VQE)(and VQS Kokail2018 )

The VQE algorithm, first introduced in Peruzzo2014 , aims at finding the
approximate ground-state energy and wavefunction of a Hamiltonian $H$ by
minimizing the energy over a parameterized trial space (also called ansatz),
i.e.:

$E_{\rm
VQE}={\min}_{\bm{\theta}}\left[\langle\Psi\left({\bm{\theta}}\right)|H|\Psi\left({\bm{\theta}}\right)\rangle\right]\equiv{\min}_{\bm{\theta}}\left[E(\bm{\theta})\right],$
(6.1)

where $\bm{\theta}\equiv\\{\theta_{p}\\}$ is a set of parameters that defines
the trial state $|\Psi({\bm{\theta}})\rangle$. In most applications, the
preparation of the trial state vector is made using a unitary transformation
$U({\bm{\theta}})$ of the qubit vacuum, denoted hereafter by
$|\bm{0}\rangle\equiv|0,\dots,0\rangle$ with:

$\displaystyle|\Psi({\bm{\theta}})\rangle$ $\displaystyle=$ $\displaystyle
U({\bm{\theta}})|\bm{0}\rangle.$ (6.2)

Some illustrative examples of trial state vectors and associated unitary
transformation will be given in section 6.2.

The parameterized observable $E(\bm{\theta})$ is further decomposed into
observables that are directly measurable on the quantum device. To do so,
being provided a fermion-to-qubit mapping (see paragraph 4.1), one can write
$H$ as qubit operator

$H=\sum_{k}\alpha_{k}P_{k},$ (6.3)

where each $P_{k}$ corresponds to a string of Pauli operators. The expectation
value of each operator over the trial state can be evaluated via sampling in
the computational basis ($Z$ measurements) with negligible gate overhead (as
illustrated on Table 3). $E(\bm{\theta})$ is then obtained by classically
aggregating the expectation value of each of the Pauli strings:

$E(\bm{\theta})=\sum_{k}\alpha_{k}\langle P_{k}\rangle_{\bm{\theta}}$ (6.4)

where we use the shorthand notation
$\langle.\rangle_{\bm{\theta}}\equiv\langle\Psi\left({\bm{\theta}}\right)|.|\Psi\left({\bm{\theta}}\right)\rangle$.

The expectation value $\langle P_{k}\rangle$ is computed by sampling many
instances of the parameterized quantum circuit to gather enough statistics to
curb statistical or "shot" noise. Numerous strategies to limit the sampling
overhead incurred by VQE have been put forward (like term-grouping strategies,
see, e.g., Huggins2019 ). Note that the associated time burden of the
algorithm can also be partly alleviated by running these circuits in parallel
on several quantum devices or parts of a large chip.

In VQE, the optimization of parameters is delegated to classical computers
using standard optimization methods, either gradient-free (like the Nelder-
Mead or COBYLA method) or gradient-based (like gradient descent). In order to
compute the gradients within gradient-based methods, one can apply either
finite difference methods, utilizing the generating function from Eq. (6.28)
(notably, this method might be challenging for NISQ devices due to the
requisite precision), or ’analytical’ gradients. The computation of the first
derivative needed in those methods can be done either with a Hadamard-test-
like circuit or with the parameter shift-rule technique Mitarai2018 ;
Schuld2019 ; izmaylov_analytic_2021 . The latter method tends to be more
suitable for noisy devices as it typically necessitates running the original
circuit twice while adjusting a single parameter. However, in more general
cases, the employment of a gate conditioned on an ancilla qubit may be needed.
Some optimization methods are particularly well-suited in that they are more
robust to the shot noise mentioned above, like the SPSA method or the
parameter-shift-rule-based "rotosolve" method Nakanishi2019 ; Ostaszewski2021
.

A schematic view of the VQE hybrid method is given in Fig. 6.1.

Figure 6.1: Schematic illustration of the VQE approach. A set of expectation
values of Pauli strings $\\{P_{k}\\}$ is obtained upon measurements on a
quantum processor, while the cost function reconstruction and the parameter
optimization are made through classical processing. The unitary
$U_{k}(\bm{\theta})$ comprises both the ansatz circuit instance
$U(\bm{\theta})$ as well as the additional basis rotation gates necessary to
access the expectation value of operator $P_{k}$ as described in Table 3.

The implementation of the VQE algorithm on current noisy devices faces several
challenges. To mention some important ones: the preparation of the initial
state (in the presence of noise), the accurate estimation of expectation
values (Eq. (6.1)) (despite shot noise), and the (classical) optimization of
the set of parameters $\bm{\theta}$ to find the minimum of $E_{\rm VQE}$.
VQE’s main advantage is the flexibility in the choice of the unitary
transformation.

Variational methods like VQE are not limited to digital quantum processors.
Elementary operations $\\{U_{R_{j}}(\theta_{j})\\}$ of the "resource"
Hamiltonian of an analog quantum computer can also be assembled into a
parameterized ansatz circuit $U_{R_{1}}(\theta_{1})U_{R_{2}}(\theta_{2})\cdots
U_{R_{k}}(\theta_{k})$ in order to minimize the energy of the target
Hamiltonian in the final state. This method is called generically Variational
Quantum Simulation (VQS, Kokail2018 ). It has been applied e.g to the
Schwinger model on trapped-ion processors Kokail2018  or to simple molecules
on Rydberg processors Michel2023 .

Some refinements of VQE are highlighted below.

####  6.1.2  Advanced VQE schemes

Various refined algorithmic schemes have been proposed to either extend the
scope of or overcome some limitations of plain-vanilla VQE.

##### Penalty methods

Other types of states than the ground state can be valuable to prepare, such
as excited states. VQE schemes with an alternative cost function have thus
been proposed to enforce the exploration of specific subspaces of interest.
For instance, leveraging the orthogonality of the eigenvectors of the
Hamiltonian, one can prepare a sequence of eigenstates by supplementing the
cost function with penalty terms proportional to the overlaps (inner products)
between the trial state $\Psi(\bm{\theta})\rangle$ and the eigenstates already
prepared $\\{|\Psi(\bm{\theta}^{*}_{j})\rangle\\}$:
$\mathcal{C}(\bm{\theta})=E(\bm{\theta})+\sum_{j}\lambda_{j}\langle\Psi(\bm{\theta}^{*}_{j})|\Psi(\bm{\theta})\rangle$.
This approach was dubbed the Variational Quantum Deflation scheme
higgott_wang_brierley_2019 . More generally, penalty terms chosen accordingly
can focus the variational search on, e.g., a specific spin sector
ryabinkin_genin_izmaylov_2018 .

##### Subspace-search VQE

A significant inconvenience of the penalty methods listed above is the need to
evaluate inner products. Hence, the proposal in nakanishi_mitarai_fujii_2019
to leverage the preservation of the inner product under unitary transformation
to look for the circuit that best maps a set of orthogonal states to the
Hamiltonian’s eigenstates.

##### ctrl-VQE

VQE can be applied more fundamentally by optimizing the control pulses which
underlie quantum operations (see Fig. 3.1) meitei_2021 . This low-level
optimization allows to address the time-limited character of coherence in NISQ
devices.

##### Orbital-rotating VQE schemes

To use shorter circuits, one had better work in a suitable one-particle
orbital basis, depending on the target state. Some VQE schemes were proposed
to tailor the basis to the problem. We can distinguish two different
approaches: (i) a "classical dressing" of the Hamiltonian observable with a
general orbital rotation whose parameters must be determined along with the
circuit’s parameters (Orbital Optimized-VQE sokolov_2020 ; koridon_2021 ), and
(ii) iterative basis updates: a converged variational state in the current
basis (starting from the usual site-spin basis in lattice models for instance)
is leveraged to extract ground state features, setting forth advantageous
updates to the single-particle basis. The terms of the Hamiltonian are
transformed accordingly before a new VQE optimization is run. The latter
strategy was applied in two different settings. The first one—dubbed permVQE
tkachenko_2021 —consists in mere single-particle basis permutations guided by
the form of the current converged state’s mutual information matrix. This
matrix measures the information shared by pairs of qubits. For a nearest-
neighbor qubit topology, one aims at concentrating its high magnitude
coefficients around the diagonal so as to lower the count of entangling gates
required. The second one—called NOization besserve_ayral_2022 —performs
general basis updates aimed at iteratively reaching a specific, state-
dependent basis known as the Natural Orbitals basis. This basis is associated
with a compact state representation, namely the one encompassing the lowest
number of computational basis states (Slater determinants). In this case, the
quantity being monitored is the 1-RDM, whose eigenbasis provides the new basis
for the subsequent VQE run.

##### Projective Quantum Eigensolver (PQE)

The PQE method Stair2021a  minimizes the residuals (that measure the non-
orthogonality of excited states to the ground state manifold) instead of the
energy, yielding accuracies on a par with VQE’s using fewer resources, and
with less size-dependence.

##### State-Average VQE

: A possible way to access excited states as well as to handle degeneracies is
to use the quantum equivalent of the state-averaged complete active space
self-consistent field used in classical computers Siegbahn1981 . This method,
dubbed State-Average-VQE (SA-VQE) Yalouz2021  provides a way to obtain excited
states in a reduced Hilbert space that can be used even in presence of
degeneracies.

####  6.1.3 On the use of variational principles in quantum computers

Although it is slightly out of the scope of the present review, we want to
mention broader applications of variational principles in quantum computing.

One can, for instance, use McLachlan’s variational principle (MVP) McLachlan64
to obtain approximate quantum systems’ unitary evolutions. In that case, the
MVP takes the form:

$\displaystyle\delta\left\|{(i\hbar\partial/\partial
t-H)\left|{\Psi(t)}\right\rangle}\right\|=0,$ (6.5)

where $\left\||\Psi\rangle\right\|\equiv\sqrt{\langle\Psi|\Psi\rangle}$. This
variational principle can be interpreted as a cost function which, given a
parametric form for the trial state, minimizes the deviation of the
approximate evolution $i\hbar\partial_{t}|\Psi(t)\rangle$ from the true
evolution $H|\Psi\rangle$. This principle can be connected with other
variational principles generally used in many-body problems Blaizot86 . A
complete discussion of the MVP and its manipulation is out of the scope of the
present article but can be found in Ref. yuan_theory_2019  in the quantum
computing context. Still, it is interesting to mention that the MVP is not
restricted to real-time unitary evolution but can be adapted to other problems
involving non-unitary motion or mixed-state evolution.

Mixed state evolution: The approximate solution of a density matrix evolution
that could differ from that of a pure state density can be obtained from the
MVP:

$\displaystyle\delta\left\|i\hbar d\rho/dt-{\cal L}(\rho)\right\|^{2}=0,$
(6.6)

where $\rho$ is the density matrix (see section 8.1 for a definition), and the
Liouvillian ${\cal L}(\rho)$ is a general functional of the density. This
variational principle can be used to find an approximation of
$i\hbar\dot{\rho}={\cal L}(\rho)$. The pure state Hamiltonian evolution can be
recovered by setting $\rho=|\Psi\rangle\langle\Psi|$ and ${\cal
L}(\rho)=[H,\rho]$. Besides this case, it can also be used to simulate
dissipative processes using a Lindblad-type equation (see section 8) for
${\cal L}(\rho)$ Endo2020 .

Imaginary-time propagation When the real-time evolution is replaced by an
imaginary-time evolution ($t\rightarrow i\tau$), an initial state with a good
enough overlap with the exact ground state is projected to the exact ground
state by the evolution operator $U(\tau)=e^{-\tau H}$. This method is a
standard practical one to obtain the lowest energies and associated
eigenstates of the Hamiltonian on a classical computer. The operator $U(\tau)$
is non-unitary and cannot a priori be directly implemented on a quantum
computer. This problem was overcome in Ref. mcardle_variational_2019 ;
McArdle2018a  using the following variant of the MVP:

$\displaystyle\delta\left\|{(\partial_{\tau}+\left[H-\langle
H\rangle_{\tau}\right])\left|{\Psi(\tau)}\right\rangle}\right\|=0.$ (6.7)

This approach is called Quantum Imaginary-Time Evolution (QITE). It can be
variational if $\ket{\Psi}$ is parameterized by variational parameters, or not
Motta2019 .

It is interesting to mention that variational techniques, that are fundamental
tools in many-body systems, have also been exported to general learning
problems. Such techniques are a growing field of interest today Mitarai2018 ;
Schuld2018 ; Benedetti2019 ; Biamonte2021 ; Schuld2021 ; Goto2021 ; Cerezo2021
.

###  6.2 Variational ansätze for many-body problems

####  6.2.1 Adiabatic state preparation vs. variational state preparation

Figure 6.2: Ground state preparation of $\ket{\Psi_{0}^{(1)}}$ from a simple-
to-prepare ground state $\ket{\Psi_{0}^{(0)}}$, either via adiabatic state
preparation or variational state preparation. The optimal variational state is
represented with a dot in the variational manifold. It is the one with the
highest overlap with $\ket{\Psi_{0}^{(1)}}$. Its distance to the target state
$\ket{\Psi_{0}^{(1)}}$ measures the expressivity of the ansatz. The orange
sphere materializes the variability on the state obtained upon optimization of
the variational ansatz, due to (i) the optimization being done with regards to
an energetic criterion rather than the overlap (ii) the presence of local
minima in the optimization landscape, (iii) shot noise (all of which incur
solely errors making the state move in the variational manifold), and (iv)
hardware noise (possibly also incurring errors kicking the state out of the
variational manifold).

The adiabatic state preparation method is an essential inspiration to many
variational state preparation techniques that we will describe. Its goal is
reach the ground state $\ket{\psi_{0}^{(1)}}$ of a Hamiltonian $H_{1}$ by
starting from the (easy-to-prepare) ground state $\ket{\psi_{0}^{(0)}}$ of an
initial Hamiltonian $H_{0}$.

This is achieved by designing a time-dependent Hamiltonian:

$H(t)=(1-s(t))H_{0}+s(t)H_{1},$ (6.8)

where $s$ denotes a "schedule" such that $s(t=0)=0$ and $s(t=T)=1$. Typically,
$\ket{\psi_{0}^{(0)}}$ can be a Slater determinant.

The adiabatic theorem guarantees that provided one proceeds slowly enough
(adiabatically) with regards to the spectral gap along the path, the system
remains in the ground state of the instantaneous Hamiltonian $H(t)$. This
property is linked to the so-called Gell-Mann and Low theorem Gell-Mann1951 ;
Fetter2012 . Upon trotterizing the evolution under perturbed Hamiltonians, one
thus has a general—but costly—recipe for ground state preparation: adiabatic
state preparation (ASP) Farhi2000 . This method was used in, e.g., Bauer2016 ,
which furthermore makes use of QPE (see paragraph 5.1 for a discussion on the
limitations of QPE for state preparation) to pin the state into the ground
state of the perturbed Hamiltonian.

By construction, ASP gives rise to long time evolutions and hence deep
circuits, which is a problem for NISQ processors. It can nevertheless be used
as a formal inspiration to design variational states to be used in the
variational methods described above: these can be regarded as a way to find
unitary "shortcuts" (sometimes called diabatic evolution) to go from a simple
initial state $\ket{\psi_{0}^{(0)}}$ to a target ground state
$\ket{\psi_{0}^{(1)}}$, as pictured on Fig. 6.2.

The general strategy to borrow from ASP to design VQE states is the following:
Hamiltonian (6.8) induces a unitary evolution that can be trotterized as
$\prod_{k=1}^{n_{t}}e^{-i(1-s(t_{k}))H_{0}\delta t}e^{-is(t_{k})H_{1}\delta
t}$, with $\delta_{t}=T/n_{t}$ and $t_{k}=k\delta_{t}$. This unitary evolution
can be transformed into a variational circuit

$U(\bm{\theta})=\prod_{k}e^{-i\theta_{2k}H_{0}}e^{-i\theta_{2k+1}H_{1}}.$
(6.9)

In the following subsections, we consider ansätze with increasing complexity
starting from a rather academic methodology (Hartree-Fock approach, Bogolyubov
transformation) and then present the recent efforts to craft more minimal
ansätze so that they avoid running into the usual limitations of VQE (limited
coherence time of NISQ devices, or the barren plateaus encountered in
numerous-parameter ansatz optimizations Romero2018 ).

####  6.2.2 Uncorrelated ansätze

##### Hartree-Fock theory and Slater determinants

In many-body problems, a complete basis of the Fock space is given by the set
of Slater determinants:

$\displaystyle|\delta_{n_{q}-1},\dots,\delta_{0}\rangle=\prod_{k}\left[a^{\dagger}_{k}\right]^{\delta_{k}}|\bm{0}\rangle,$
(6.10)

where $\\{a^{\dagger}_{k}\\}_{k=0,n_{q}-1}$ correspond to the creation
operators associated with a complete set of single-particle states
$\\{|k\rangle\\}$. Here it is assumed that the one-body Hilbert space is
finite with dimension $n_{q}$. For a set of $A$ particles, the Hartree-Fock
procedure consists in estimating the ground state energy of a (possibly
correlated) Hamiltonian ${H}$ as $E_{\rm HF}=\langle\Psi_{\rm HF}|H|\Psi_{\rm
HF}\rangle$ assuming that the trial wave-function is a Slater determinant
given by

$\displaystyle|\Psi_{\rm HF}\rangle$ $\displaystyle=$
$\displaystyle\prod_{\alpha}\left[b^{\dagger}_{\alpha}\right]^{\gamma_{\alpha}}|\bm{0}\rangle.$
(6.11)

In this equation, only $A$ coefficients $\gamma_{\alpha}$ are equal to $1$
(corresponding to particles below the Fermi level) while $n_{q}-A$ of them are
equal to zero (particle states). The creations operators
$\\{b^{\dagger}_{\alpha}\\}$ are associated to a complete basis
$\\{|\alpha\rangle\\}$, with the relationship:

$\displaystyle|\alpha\rangle$ $\displaystyle=$
$\displaystyle\sum_{k}|k\rangle\langle k|\alpha\rangle\leavevmode\nobreak\
\longrightarrow b^{\dagger}_{\alpha}=\sum_{k}a^{\dagger}_{k}U_{k\alpha},$
(6.12)

with $U_{k\alpha}=\langle k|\alpha\rangle$ a unitary transformation to be
variationally determined. In other words, the HF procedure consists in
variationally finding the best Slater determinant approximation to the target
ground state. It is a very simple mean-field approach to the original problem:
it replaces the many-body problem with a set of particles influencing one
another through a self-consistent, average one-body potential (or mean field),
instead of actual interactions. Practical aspects of finding the HF solution,
which consists in minimizing the energy with respect to the variations of
unitary matrix $U$, are standard.

Below, we focus on the implementation of HF on a quantum computer because the
preparation of arbitrary Slater determinants is the starting point of many
more advanced methods (indeed, implementing HF on a quantum computer is per se
not useful as HF can be implemented efficiently—i.e in polynomial time—on a
classical computer).

When the Jordan-Wigner transformation is used to map the Fock space to qubit
space (see section 4.1), the set of Slater determinants, given by Eq. (6.10),
directly identifies with the computational basis
$\\{|\delta_{n_{q}-1},\dots,\delta_{0}\rangle\\}$ with
$\\{\delta_{k}=0,1\\}_{k=0,n_{q}-1}$, that we can rewrite as:

$\displaystyle|\delta_{n_{q}-1},\dots,\delta_{0}\rangle=\prod_{k}\left[X_{k}\right]^{\delta_{k}}|\bm{0}\rangle.$
(6.13)

The key ingredient to implementing the HF theory on a quantum computer is to
be able to realize a general unitary transformation $U$ on the fermionic
modes, as defined in Equation (6.12), using a parametric circuit and the
possibility to prepare a state equivalent to (6.11) with exactly $A$
particles. The technique that has been employed, for instance, in Ref.
Dumitrescu2018 ; Arute2020 , is based on the use of the Thouless theorem
Thouless1960  (see also Appendix E of Ref. Ring80 ). Starting from one of the
Slater determinants given in (6.10) and denoted generically as
$|\Psi_{0}\rangle$, one can generate an ensemble of new Slater determinants
given by

$\displaystyle|\Psi({Z})\rangle$ $\displaystyle\equiv$ $\displaystyle
e^{i\sum_{ij}Z_{ij}a^{\dagger}_{i}a_{j}}|\Psi_{0}\rangle,$ (6.14)

where $Z$ is supposed to be hermitian. These states identify with the form
(6.11) with $b^{\dagger}_{\alpha}({Z})=\sum_{k}(e^{i{Z}})_{\alpha
k}a^{\dagger}_{k}$, i.e. $U^{*}=e^{iZ}$. The circuit that performs the
Thouless transformation starting from a product state as given by Eq. (6.13)
is discussed in detail in Ref. Arute2020 . Additional discussion on the
application of HF can be found in Refs. Wecker2014a ; McArdle2018a .

##### General quasi-particle vacuum

The transformation (6.12) and the Thouless method can be generalized to a
larger class of trial states known as quasi-particle vacua (also known as
Gaussian states in other contexts). Below, such vacuum states are written
generically as $|\Psi_{\beta}\rangle\propto\prod\beta_{\alpha}|\bm{0}\rangle$,
where $\\{\beta_{\alpha},\beta^{\dagger}_{\alpha}\\}$ denotes a complete set
of quasi-particles creation/annihilation operators. These operators can be
connected through a generalization of Eq. (6.12) given by Ring80 ; Blaizot86 :

$\displaystyle\beta^{\dagger}_{k}$ $\displaystyle=$
$\displaystyle\sum_{k}\left[a^{\dagger}_{k}U_{k\alpha}+a_{k}V_{k\alpha}\right].$
(6.15)

Using general quasi-particle vacua instead of restricting to Slater
determinants leads to the Hartree-Fock-Bogolyubov (HFB) theory, where the
$U(1)$ symmetry, associated with particle number conservation, is broken. The
advantage of breaking this symmetry is the possibility to describe superfluid
systems brink_broglia_2005 .

The quantum state preparation of a general quasi-particle vacuum using the
Thouless transformation, as done at the HF level in Arute2020 , was addressed
in, e.g., Dallaire-Demers2018  and relies on two main arguments. First, the
mapping between Thouless’ transformation
$\mathcal{U}(Z)=e^{i\sum_{ij}Z_{ij}\gamma_{i}\gamma_{j}}$ (where the
$\gamma_{i}$ denote Majorana modes, i.e., $\gamma_{2k}=a_{k}+a_{k}^{\dagger}$
and $\gamma_{2k+1}=-i(a_{k}-a_{k}^{\dagger})$) and quantum gates can be found
by leveraging the decomposition of $R=e^{iZ}$ as a product of elementary
Givens rotations hoffman_1972 . Let $M$ be the number of fermionic modes (and
hence, qubits), then $R=\prod\limits_{k=1}^{2M\choose 2}r_{k}(\theta_{k})$,
where the Givens rotations consist of $M$ local phase rotations which can be
implemented as $R_{z}$ gates and $2M(M-1)$ $SO(4)$ rotations acting non-
trivially only on two modes. The $SO(4)$ rotations can be rendered in a
quantum circuit by means of matchgates valiant_2001 ; jozsa_miyake_2008  and
SWAP gates. All in all, Thouless’ transformation can be implemented as a
nearest-neighbour matchgate circuit with depth $O(M)$.

Several applications have also been explored in quantum computers using a
Bardeen–Cooper–Schrieffer (BCS)-like ansatz given by:

$\displaystyle|\Psi\rangle$ $\displaystyle=$
$\displaystyle\prod_{k}\left(u_{k}+v_{k}a^{\dagger}_{k}a^{\dagger}_{\bar{k}}\right)|\bm{0}\rangle,$
(6.16)

where $(k,\bar{k})$ refers to two single-particle states forming a pair of
time-reversed states, and where $u^{2}_{k}+v^{2}_{k}=1$. The encoding of such
a state on a qubit register is not unique. If the brute force JWT is used to
make a direct mapping between single-particle states and qubits, the trial
state can be written as:

$\displaystyle|\Psi(\bm{\theta})\rangle$ $\displaystyle=$
$\displaystyle\bigotimes_{k}\left[\sin(\theta_{k})|00\rangle_{k}+\cos(\theta_{k})|11\rangle_{k}\right],$
(6.17)

where we made the identification $u_{k}=\sin(\theta_{k})$ and
$v_{k}=\cos(\theta_{k})$. To map Eq. (6.16) into Eq. (6.17), time-reversed
states are assumed to be represented by adjacent qubits, and $|.\rangle_{k}$
denotes the two qubits associated with these states. We recognize a
generalized Bell state that can be obtained by performing a
$R_{y}(\theta_{k})$ rotation on one of the qubits followed by a CNOT operation
with the second qubit. Such encoding was used, for instance, in Verstraete2000
; ovrum2007quantum ; Jiang2018 ; Lacroix2020 ; Khamoshi2020 ; Khamoshi2022 .

This encoding is general and allows treatment of systems where one or several
pairs are broken (usually referred to as nonzero seniority brink_broglia_2005
) as illustrated in Ref. Lacroix2020  for instance, to treat odd systems. If
we restrict to the situation with seniority $0$, i.e., when no pairs are
broken, one can reduce the number of qubits by directly encoding the
occupation of the two adjacent time-reversed states onto one qubit. In this
case, the state $|1\rangle_{k}$ or $|0\rangle_{k}$ represents the simultaneous
occupation or not of the two time-reversed particles $(k,\bar{k})$. This
technique, used in Refs. Khamoshi2020 ; ruiz-guzman2021 ; RuizGuzman2022 , has
the advantage of reducing by a factor of $2$ the required number of qubits
compared to the case where one particle is encoded on one qubit. It also
avoids the use of controlled operations since we have, for this encoding
scheme:

$\displaystyle|\Psi(\bm{\theta})\rangle$ $\displaystyle=$
$\displaystyle\bigotimes_{k}\left[\sin{\theta_{k}}|0\rangle_{k}+\cos{\theta_{k}}|1\rangle_{k}\right]$
(6.18) $\displaystyle=$
$\displaystyle\prod_{k}R^{(k)}_{y}\left(\pi-2\theta_{k}\right)|\bm{0}\rangle.$

Quasiparticle-like states have been extensively explored in quantum computers
Wecker2015a ; Kivlichan2018  (see also Dallaire-Demers2018 ), as well as their
experimental implementation, Arute2020 .

We focused here on relatively standard quasi-particle vacuum states that play
a particular role in many-body systems and lead to the HF and HFB frameworks.
Below is a selection of other ansätze that are widely discussed today in the
literature.

####  6.2.3 Hamiltonian Variational Ansatz

Inspired both from adiabatic state preparation (ASP, see Fig. 6.2) and the
Quantum Approximate Optimization Algorithm circuit (QAOA Farhi_2014 ), which
is ubiquitous in quantum combinatorial optimization, the Hamiltonian
Variational Ansatz (HVA wecker_2015 ) state reads

$\ket{\Psi(\bm{\theta})}=\prod_{l=1}^{L}\left(\prod_{k}e^{-i\theta_{k}^{l}H_{k}}\right)\ket{\Psi_{0}},$
(6.19)

where the terms $H_{k}$ come from decomposing $H$ as $H=\sum_{k}H_{k}$, with
$[H_{k},H_{k^{\prime}}]\neq 0$ for $k\neq k^{\prime}$. The dimension $L$ of
index $l$ is referred to as the depth of the ansatz. The initial state
$\ket{\Psi_{0}}$ is the ground state of some $H_{k_{0}}$ that is taken not to
act first on $\ket{\Psi_{0}}$. Optimizing the HVA parameters amounts to
optimizing the Hamiltonian schedule $s$ of ASP (as introduced in Eq. (6.8)).
HVA was applied, e.g. to the study of the 1-D Hubbard model in
anselme_martin_2022 .

####  6.2.4 Hardware-Efficient Ansatz (HEA)

Today’s quantum computers are not uniformly accurate in performing different
operations. Knowing the strength or weaknesses of a given platform, one might
adapt the ansatz to the operations most efficiently realized. The HEA
technique consists in writing the trial state from a set of operations that
are "native" in the quantum processor. This heuristic approach can optimize
the trial state construction with respect to the specific hardware, but also
restricts the type of trial states that can be constructed Kandala2017 ;
Barkoutsos2018 . One issue with the HEA is that the classical optimization of
the variational parameters can become difficult due to gradients that vanish
exponentially with the number of qubits, a phenomenon dubbed the barren
plateau problem McClean2018 . Interestingly, it was observed on the other hand
that beyond a certain depth at a given number of qubits, gradients increased
again for e.g. the transverse-field Ising model as well as for the Sachdev-Ye-
Kitaev model Kim_2021 . This phenomenon is general, and is easier to
understand with the help of the concept of overparametrization Larocca_2023 .
Moreover, HEA are typically defined with gates taken from the native gateset
of the specific technology that is used, and are therefore strongly hardware-
driven. This might in practice lead to quantum circuits that do not respect
some of the symmetries of the physical problem of interest, leading to
specific additional difficulties.

####  6.2.5 Unitary coupled cluster (UCC)

The UCC offers a framework that naturally extends the HF method based on the
Thouless approach described in section 6.2.2 (see also section 2.1). The trial
wavefunction is written in a generalized form Romero2018 :

$\displaystyle|\Psi_{\rm UCC}(\bm{\theta})\rangle$ $\displaystyle=$
$\displaystyle
e^{T(\bm{\theta})-T^{\dagger}(\bm{\theta})}|\Psi_{0}\rangle\equiv
U(\bm{\theta})|\Psi_{0}\rangle,$ (6.20)

where $T$ can be expanded as a set of operators of increasing complexity with
$T=T_{1}+T_{2}+\dots$. Here, $T_{1}$, $T_{2}$, … stands for single, double, …
particle-hole excitation operators with respect to the state
$|\Psi_{0}\rangle$, with

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\sum_{i,j}T^{(1)}_{ij}a^{\dagger}_{i}a_{j},$ (6.21)
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\sum_{i,j}T^{(2)}_{ij,kl}a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k},$
$\displaystyle\cdots$

We see, in particular, that the HF is recovered by using a Slater determinant
and restricting $T$ to single excitations. After truncation, the state
prepared using Eq. (6.20) is used as a trial state in the VQE approach
discussed in section 6.1.1. This technique is currently widely applied in
quantum chemistry (see the recent review Abhinav2022  and references therein)
, but mostly in noise-free simulations. The reason is that whereas it has the
capacity to yield highly accurate results, the very large depth of the UCC
ansatz circuit leads to degraded performances on real quantum platforms. Ways
around this prohibitive gate count include circuit compilation—tailoring the
implementation to the subspace under scrutiny, as in Ref. Dumitrescu2018 , as
well as adaptive approaches, for instance as applied to atomic nuclei in Ref.
Kiss2022 . Finally, the possibility of combining such an approach with $U(1)$
symmetry that is relevant for strongly interacting systems like nuclei is of
current interest Henderson2014 ; Qiu2019  and is currently explored on quantum
computers too Khamoshi2020 ; Khamoshi2022 .

####  6.2.6 The Low-Depth Circuit Ansatz (LDCA)

Elaborating on the general quasiparticle vacua preparation routine reviewed in
paragraph 6.2.2, the LDCA circuit Dallaire-Demers2018  (standing for Low-Depth
Circuit Ansatz) possibly allows reaching any state due to the insertion of
$R_{zz}$ rotation gates into to the Hartree-Fock-Bogolyubov circuit. It also
allows the replication of similar layers in the ansatz to increase its
representability systematically. Intuitively, $R_{zz}(\theta)\equiv
e^{-i\theta/2{Z}_{p}{Z}_{q}}$ gates generate correlated states as these gates
correspond to density-density interactions within a Jordan-Wigner encoding:
the density operator of orbital $p$, $n_{p}\equiv a^{\dagger}_{p}a_{p}$, maps
to the qubit operator $\frac{I-{Z}_{p}}{2}$; therefore $n_{q}n_{p}$
interactions translate into $Z_{p}Z_{q}$ terms.

The LDCA is competitive compared to the UCC approach in terms of circuit
depth, especially when double or higher excitations are considered. Despite
its gentle scaling in terms of the number of quantum operations, it still
incurs a prohibitive gate count with respect to NISQ capacities, e.g.,
besserve_ayral_2022 .

####  6.2.7 Projected Ansätze

An important cornerstone for future applications, especially in nuclear
systems, is the possibility of making, for instance, symmetry restoration
after symmetry breaking. A discussion for the particle number symmetry was
made in section 6.2.2. Assume, for instance, that a state
$|\Psi(\bm{\theta})\rangle$ can be prepared on a quantum computer and that
such state does not respect the symmetry of the physical problem that is
encoded in the Hamiltonian $H$. Instead of using $|\Psi(\bm{\theta})\rangle$
in the variational principle, one can use the projected wavefunction:

$\displaystyle|\Psi_{\rm P}^{\prime}(\bm{\theta})\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\langle\Psi(\bm{\theta})|P_{S}|\Psi(\bm{\theta})\rangle}}P_{S}|\Psi(\bm{\theta})\rangle.$
(6.23)

Here, $P_{S}$ is a projector onto the subspace of the Hilbert space containing
states with the desired property (for instance, the proper symmetries). Such a
strategy is used in many-body physics to grasp specific correlations between
particles or strongly entangled states that are hard to describe otherwise
Ring80 ; Blaizot86 . One difficulty is that the projector is a non-unitary
operation and cannot be directly implemented on a quantum computer. Several
methods have been proposed recently to construct projected states Lacroix2020
; Lacroix2022 . These methods have been combined with VQE in Ref.
RuizGuzman2022 , leading to the Quantum-Variation After Projection (Q-VAP)
framework.

####  6.2.8 Adapt-VQE

Plain-vanilla VQE takes a fixed variational state as an input to the
computation. This incurs the risk of overfitting the target state if the
variational manifold is too large. The ADAPT-VQE method iteratively constructs
this ansatz instead. It relies on a predefined operator pool from which
operators are drawn adaptively along the optimization procedure. Typically one
selects the operator maximizing the gradient at the current step so that its
addition to the circuit has the most significant effect on the variational
energy. The aim is to reduce the number of parameters of the ansatz at the
expense of an increased measurement overhead due to the gradients. In the
initial proposal (ADAPT-VQE grimsley_economou_barnes_mayhall_2019 )), the
operators in the pool were fermionic, but the large gate overheads resulting
from long Jordan-Wigner strings can be avoided by directly using qubit
operators instead (qubit-ADAPT VQE tang_shkolnikov_2021 ). This method can
reach chemical accuracy at a relatively low gate count, at least with
noiseless computers (see, e.g., Haidar2022  for an example with large
molecules).

####  6.2.9 Tensor-network-inspired quantum circuits

Tensor-network states refer to widely used representations of the wave
function $|\Psi\rangle$ of a many-body problem. Instead of storing the
information contained in a $n$-qubit wavefunction in a multi-array
$a_{b_{1},b_{2},\dots b_{n}}$ (with storage cost $2^{n}$ complex floating-
point numbers), one assumes a particular factorization of this multi-array,
with the hope that storing the different factors will be less costly. The
Matrix Product State (MPS) class is a widespread subclass of tensor networks.
It consists in factorizing the multi-array as a product of matrices
$[A^{(k)}]^{b_{k}}$:

$a_{b_{1},b_{2},\dots
b_{n}}=\sum_{\alpha_{1},\dots\alpha_{n-1}}[A^{(1)}]^{b_{1}}_{\alpha_{1}}[A^{(2)}]^{b_{2}}_{\alpha_{1},\alpha_{2}}\cdots[A^{(n)}]^{b_{n}}_{\alpha_{n-1}}$
(6.24)

The internal indices $\alpha_{k}$ have a dimension called the bond dimension,
which is usually denoted as $\chi$. As we will see in section 7.3, this
parameter is closely connected to the degree of entanglement one can access
with such a state.

MPS can be used as an inspiration or starting point for quantum computations.
Specifically, methods have been proposed to convert a given MPS to a quantum
circuit Ran2020 ; Rudolph2022 . Paragraph 7.3.3 presents an outline of the
quantum circuitry involved. An important outcome is that MPS with a given bond
dimension requires a circuit with a logarithmic depth in $\chi$, resulting in
a complexity gain and thus the possibility to use a quantum computer to
generate states with an entanglement level inaccessible to classical computers
due to too large a bond dimension. The conversion from MPS to a circuit can
also be used to warm-start a variational quantum computation Rudolph2022a .

Other tensor networks can also be used as a comparison point or inspiration to
quantum circuit design, like Tree Tensor Networks Yuan2021a , the Multiscale
Entanglement Renormalization Ansatz Miao2021 , or Projected Entangled Pair
States Cirac2021 .

###  6.3 Beyond variational methods

As described in section 6.1.1, variational methods usually target the search
of approximate ground states. Here, we discuss several methods that can give
access to excited states too. Previously (section 5.1), we saw that the QPE
could be a tool of choice to obtain energy eigenvalues and associated
eigenvectors. Unfortunately, it cannot be used in current devices, and
probably, it will take some time before the fidelity of quantum machines
becomes sufficiently high to apply it. The search for alternative methods,
less costly in terms of quantum resources, is therefore an intensive domain of
activity today. We report below some of the methods that have been proposed to
access excited state properties. Note that these methods sometimes only give
access to the energies, not the states associated with them.

####  6.3.1 Quantum Subspace Expansion methods

A common strategy used in classical computers to obtain the approximate
solution of a diagonalization problem when a complete CI solution is
prohibitive consists of iteratively constructing subspaces of the total
Hilbert space ${\cal H}$ of increasing complexity. In many cases, at a given
level $M$ of complexity, a subspace is generated by a set of states
$\\{|\Psi_{0}\rangle,\dots,|\Psi_{M-1}\rangle\\}$ that spans a subspace
denoted as ${\cal H}_{M}$. The method to obtain the states is usually
iterative in the sense that $|\Psi_{k+1}\rangle$ is constructed from
$|\Psi_{k}\rangle$ using specific operations. Arnoldi or Lanczos methods are
famous examples widely used on classical computers Saad2011 . The generic
strategy of using an increasing number of states to form a subspace of the
Hilbert space will be called hereafter Quantum Subspace Expansion (QSE)
McClean2017a . A schematic view of the QSE strategy is shown in Fig. 6.3. This
strategy’s success depends on its capacity to construct the relevant subspace
for a given problem. Key ingredients are the seed state $|\Psi_{0}\rangle$ and
the rules for the iterative generation of states.

Figure 6.3: Illustration of QSE philosophy where the eigenvalue problem is
considered in a subspace with increasing dimension.

The generated states are usually not orthogonal with one another. Any state
belonging to the reduced space ${\cal H}_{M}$ can be written as
$|\Psi\rangle=\sum\limits_{K=0}^{M-1}c_{K}|\Psi_{K}\rangle$. An approximate
eigenvalue $E$ in this space can be obtained by solving the generalized
eigenvalue problem written as a set of $K=0,M-1$ equations given by:

$\displaystyle\sum\limits_{K^{\prime}=0}^{M-1}c_{K^{\prime}}H_{KK^{\prime}}$
$\displaystyle=$ $\displaystyle
E\sum\limits_{K^{\prime}=0}^{M-1}c_{K^{\prime}}O_{KK^{\prime}},$ (6.25)

with $H_{KK^{\prime}}=\langle\Psi_{K}|H|\Psi_{K^{\prime}}\rangle$ and
$O_{KK^{\prime}}=\langle\Psi_{K}|\Psi_{K^{\prime}}\rangle$ the overlap matrix.
Such equations can be solved in a two-step process by first diagonalizing the
overlap matrix prior to the Hamiltonian diagonalization. Potentially, in
${\cal H}_{M}$, $M$ approximate eigenstates can be obtained, which makes the
method quite attractive.

On classical computers, a typical choice for the states is the Krylov basis
where $|\Psi_{k+1}\rangle=H|\Psi_{k}\rangle$. Great effort is currently
devoted to the possibility of extending the Krylov space technique to quantum
computers. Then, the strategy is to compute the matrix elements of the two
matrices $O$ and $H$ in the reduced space using the quantum computer.
Subsequently, this matrix is diagonalized via classical methods Seki2021 . A
discussion on the possibility of obtaining the strict equivalent of the Krylov
basis using derivatives of the generating function $F(t)$ introduced below in
Eq. (6.28) was scrutinized in Ref. ruiz-guzman2021 . This analysis was made
using the fact that
$H_{KK^{\prime}}=\langle\Psi_{0}|H^{K+K^{\prime}+1}|\Psi_{0}\rangle$ and that
$F(t)$ is the generating function of the Hamiltonian moments. However, this
method is exceptionally susceptible to numerical noise.

Alternatively, one can generate the states using unitary transformation of the
seed state such that $|\Psi_{K}\rangle=U_{K}|\Psi_{0}\rangle$. This approach
is particularly well adapted to quantum computing, where circuits are
automatically unitary. In the Quantum Krylov technique, the hamiltonian
propagator itself is used such that $U_{K}\equiv e^{-iH\tau_{K}}$ where a set
of times $\\{\tau_{K}\\}_{K=0,\dots,M}$ has been assumed with the convention
$\tau_{0}=0$. The possibility of using Krylov-inspired techniques on quantum
computers has more generally attracted much attention in recent years
McClean2017a ; Parrish2019 ; Stair2020 ; Seki2021 ; Bharti2021b ; Bharti2021c
; Bespalova2021 ; Jamet2021 ; Cortes2022 ; Seki2022 ; Haug2022 ; Wei2022 ;
RuizGuzman2022  (see also the survey Aulicino2022 ).

####  6.3.2 QPE-inspired quantum algorithms

Different methods inspired by the QPE algorithm have been proposed to reduce
the quantum resources in ancilla qubits or the number of operations in the
quantum circuit. For instance, the methods we will discuss subsequently use
only one ancilla qubit and, contrary to the QPE algorithm, do not require the
inverse quantum Fourier Transform. Ref. McArdle2018a  discusses a
comprehensive list of these methods.

As an alternative to the standard phase estimation, Kitaev’s algorithm
Kitaev1996  and the iterative QPE algorithm based on the semiclassical quantum
Fourier transform Griffiths1996 ; Dobsicek2007  (see also Svore2013 ; Lin2022
) were proposed to find the eigenvalue of a single eigenstate. More recently,
further progress has been made with the Rodeo algorithm Choi2021  that appears
as a practical tool in the NISQ context Qian2021 ; Bee-Lindgren2022 . We
briefly describe below how these iterative techniques can be implemented.

Figure 6.4: Circuit used for the Rodeo algorithm Choi2021 . Given that there
are no entangling gates between the ancillary qubits, it is possible to use a
procedure where only one is used and measured multiple times.

We consider again an initial state $|\Psi\rangle$ that decomposes onto the
Hamiltonian eigenstates $\\{|\alpha\rangle\\}$ (associated to a set of
eigenvalues $\\{E_{\alpha}\\}$) as
$|\Psi\rangle=\sum_{\alpha}c_{\alpha}|\alpha\rangle$. Provided that the
initial state is prepared on a quantum register, the circuit depicted in Fig.
6.4 is applied, and the ensemble of ancilla qubits is measured. The parameter
$E$ appearing as a scaling factor in the phase rotations shown in Fig. 6.4 can
be freely varied. Given that there are no entangling gates between the ancilla
wires in Fig. 6.4, it is possible to see this circuit as a consecutive series
of measurements on a single qubit, similar to the iterative QPE procedure
Dobsicek2007 . To be more specific, assuming $n_{a}$ indirect measurements, a
set of times $(\tau_{1},\dots,\tau_{n_{a}})$ are considered and the controlled
operation of the $j^{\mathrm{th}}$ measurement is made using
$U(\tau_{j})=e^{-i{H}\tau_{j}}$ with ${H}$ the Hamiltonian. It can then be
shown Choi2021 ; Qian2021  that the probability to obtain only the $|0\rangle$
state in all of the consecutive $n_{a}$ measurements of the ancilla qubits is:

$\displaystyle p_{0^{n_{a}}}\left(E,\\{\tau_{i}\\}\right)$ $\displaystyle=$
$\displaystyle\sum_{\alpha}|c_{\alpha}|^{2}\prod_{i=1}^{n_{a}}\cos^{2}\left(\left(E_{\alpha}-E\right)\frac{\tau_{i}}{2}\right).$
(6.26)

As the number of repetitions $n_{a}$ increases, the above function of $E$
peaks around the $E_{\alpha}$ values. The flexibility in choosing the
$\\{\tau_{i}\\}$ values can be further used to improve the convergence. Below,
we discuss two main options:

(i) Fixed times prescription: We can assume that we have
$\tau_{i}=\tau/2^{i-1}$ and $\tau=\frac{\pi 2^{n_{a}-2}}{|E_{up}-E_{low}|}$
where $E_{up}$ (resp. $E_{low}$) is an upper bound (resp. lower bound) on the
spectrum of the Hamiltonian. An example of the resulting probability
$p_{0^{n_{a}}}$ given by Eq. (6.26) at various $E$ is shown in Fig. 6.5. From
the positions of the peaks in the distribution, as well as their amplitude, we
can extract approximate eigenenergies and the weights of the associated
eigenstates in the decomposition of $\ket{\Psi}$, $|c_{\alpha}|^{2}$. (ii)
Rodeo prescription: The key idea behind the Rodeo method is to assume a
Gaussian statistical ensemble of times $\\{\tau_{i}\\}$ with an adjustable
Gaussian width $\sigma$. Averaging over the statistical ensemble gives the
probability:

$\centering
p_{0^{n_{a}}}\left(E\right)=\sum_{\alpha}|c_{\alpha}|^{2}\left[\frac{1+e^{-\left(E_{\alpha}-E\right)^{2}\sigma^{2}/2}}{2}\right]^{n_{a}},\@add@centering$
(6.27)

that is also strongly peaked around the eigenenergies. An illustration of the
Rodeo prescription is also given in Fig. 6.5. The Rodeo method has two
advantages compared to the fixed times’ (original Kitaev method) approach.
First, the probabilities are flattened away from the energies, which helps to
identify peaks. Second, the extra parameter $\sigma$ can be used as a
resolution to rapidly scan a given energy range (see Qian2021 ).

Figure 6.5: Illustration of the iterative methods using the fixed times (black