††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work††thanks: These authors contributed
equally to this work

# Walking through Hilbert Space with Quantum Computers

Tong Jiang  Department of Chemistry and Chemical Biology, Harvard University,
Cambridge, MA, USA  Jinghong Zhang  Department of Chemistry and Chemical
Biology, Harvard University, Cambridge, MA, USA  Moritz K. A. Baumgarten
Department of Chemistry and Chemical Biology, Harvard University, Cambridge,
MA, USA  Meng-Fu Chen  Department of Chemistry and Chemical Biology, Harvard
University, Cambridge, MA, USA  Hieu Q. Dinh  Department of Chemistry and
Chemical Biology, Harvard University, Cambridge, MA, USA  Aadithya Ganeshram
Department of Chemistry and Chemical Biology, Harvard University, Cambridge,
MA, USA  Nishad Maskara  Department of Physics, Harvard University, Cambridge,
MA, USA  Anton Ni  Department of Chemistry and Chemical Biology, Harvard
University, Cambridge, MA, USA  Joonho Lee<EMAIL_ADDRESS>Department
of Chemistry and Chemical Biology, Harvard University, Cambridge, MA, USA
Google Research, Venice, CA 90291, United States

###### Abstract

Computations of chemical systems’ equilibrium properties and non-equilibrium
dynamics have been suspected of being a "killer app" for quantum computers.
This review highlights the recent advancements of quantum algorithms tackling
complex sampling tasks in the key areas of computational chemistry: ground
state, thermal state properties, and quantum dynamics calculations. We review
a broad range of quantum algorithms, from hybrid quantum-classical to fully
quantum, focusing on the traditional Monte Carlo family, including Markov
chain Monte Carlo, variational Monte Carlo, projector Monte Carlo, path
integral Monte Carlo, etc. We also cover other relevant techniques involving
complex sampling tasks, such as quantum-selected configuration interaction,
minimally entangled typical thermal states, entanglement forging, and Monte
Carlo-flavored Lindbladian dynamics. We provide a comprehensive overview of
these algorithms’ classical and quantum counterparts, detailing their
theoretical frameworks and discussing the potentials and challenges in
achieving quantum computational advantages.

##  I Introduction

Figure 1:  Quantum Algorithms for Sampling Problems in Computational Chemistry
Covered in this Review. This figure shows the integration of quantum computing
for ground state, thermal state calculations, and quantum dynamics. Ground
state algorithms encompass hybrid quantum-classical algorithms, including
auxiliary-field quantum Monte Carlo (AFQMC), full configuration interaction
quantum Monte Carlo (FCIQMC), unitary coupled-cluster Monte Carlo (UCCMC),
variational Monte Carlo (VMC), selected configuration interaction (Selected
CI), and entanglement forging. Cooling with Monte Carlo-flavored Lindbladians
is covered as a fully quantum algorithm for ground state calculation. In the
thermal state calculation category, the algorithms include quantum-enhanced
classical Monte Carlo, stochastic series expansion (SSE), quantum Metropolis,
and minimally entangled typical thermal states (METTS) as hybrid methods,
complemented by quantum-quantum Metropolis and Monte Carlo-flavored
Lindbladians as fully quantum algorithms. For quantum dynamics, we covered
quantum circuit Monte Carlo, hybrid path-integral Monte Carlo (PIMC), and
Boson sampling. Each hybrid method connects back to the central unit,
illustrating the intermediate quantities quantum computing contributes to this
methodology.

Quantum chemistry strives to obtain accurate, approximate solutions to
Schrödinger’s equation as efficiently as possible [1, 2]. The primary aim is
to extend the capabilities of computational methods to handle calculations
over larger spatial and longer temporal scales through algorithmic advances.
Despite the naïve exponential (i.e., computationally intractable) complexity,
progress in quantum chemistry has never ceased. The range of applications
enabled by quantum chemistry has been ever growing from drug discovery to
materials design [3, 4, 5, 6, 7].

Owing to its inherent relevance to practical applications, quantum chemistry
has recently received significant interest from quantum computing communities,
both industry and academia [8, 9, 6, 10, 11]. A key goal is finding algorithms
and applications to achieve a “practical” quantum advantage. Practical quantum
advantages can be assessed in two ways. One way is to see if the end-to-end
run time is shorter for a fixed error with access to the quantum computer.
Another way is to see if we can perform more accurate quantum chemistry
calculations for a fixed runtime with access to the quantum computer. Even
beyond the near-term intermediate scale (NISQ) architectures, we do not have a
definitive example with practical quantum advantages known in the field.

One prominent example of NISQ and beyond-NISQ algorithms is variational
algorithms stemming from the variational quantum eigensolver (VQE) [12, 13,
14]. The central idea of VQE is to optimize a relevant objective function for
given problems by varying quantum circuit parameters [12]. The limitations of
this approach stem from the limited expressivity of finite-depth quantum
circuits and the optimization difficulties involving many non-linear
parameters accompanied by vanishing gradients exhibiting barren plateau [15].
In other words, one may fail to represent the true solution with a given
circuit ansatz or to find the ground state even with an expressive circuit due
to optimization difficulties [16]. Nonetheless, developing different types of
VQE algorithms has become one of the most active research directions in
quantum chemistry over the last few years. Interested readers are referred to
review articles on this subject [14].

Another commonly discussed quantum algorithm is called quantum phase
estimation (QPE) [17], geared towards fault-tolerant quantum computers. The
goal of QPE is to estimate an eigenvalue (usually the lowest one) of a
hermitian operator (Hamiltonian), $\hat{H}$. In QPE, an initial state with
nonzero overlap with the state of interest (typically the ground state) is
prepared and then time-evolved under the Hamiltonian. The target eigenvalue is
subsequently obtained through the inverse quantum Fourier transform. The
probability of finding the target eigenvalue is inversely proportional to the
overlap between the exact eigenstate and the initial state. As one grows
system size, this overlap decays exponentially with system size if one uses an
initial state with a fixed energy per unit volume [18]. While many
advancements in recent years economized the implementation of the necessary
time evolution with sophisticated quantum chemistry Hamiltonians [19, 20], the
cost of state preparation associated with this overlap makes it unclear
whether QPE can demonstrate quantum advantages in computational chemistry [18,
21, 22].

This Review summarizes a relatively less discussed class of quantum algorithms
developed to sample probability distributions relevant to computational
chemistry. On a classical computer, such sampling is most commonly performed
by Monte Carlo (MC) algorithms. Depending on the problem setup, classical MC
sampling may face significant computational challenges, either due to the
sample complexity (i.e., the number of samples required to obtain
statistically reliable answers) or the computational cost per statistical
sample. Despite their prevalence in computational chemistry, efforts to
enhance MC algorithms and develop sampling techniques using quantum computing
have been relatively limited in the community. On the one hand, our goal is to
provide computational chemists with an overview of quantum algorithms that may
help improve routine computational chemistry tasks. On the other hand, we want
to provide quantum information scientists with some background on the existing
computational chemistry methods and connections between classical and quantum
algorithms. For the fundamental principles and basic concepts in quantum
computing, we refer the readers to Nielsen and Chuang [23] and to the review
papers [6, 9] that are more specialized for computational chemistry. We also
reference another review paper that discusses sampling problems in quantum
computing [24]. We hope to provide complementary viewpoints and more
background for common computational chemistry sampling problems.

This review is organized as follows: Section II discusses sampling problems
and approaches in computational chemistry, including ground-state, finite-
temperature, and real-time dynamics problems. This section sets the foundation
for the subsequent exploration of quantum algorithms that aim to enhance these
classical approaches. Section III delves into quantum algorithms for ground-
state calculations, exploring various quantum-classical hybrid and quantum
algorithms. Section IV introduces quantum algorithms for thermal state
properties calculations. Section V focuses on quantum algorithms for quantum
dynamics. Finally, Section VI provides a summary and outlook.

##  II Sampling problems in computational chemistry

In computational chemistry, computational kernels often provide a sample from
an underlying distribution of interest. Classically, the most common approach
is Monte Carlo sampling on the classical computer to directly estimate local
observables (e.g., energy and correlation functions) from the distribution
[25, 26, 27]. Here, we provide several contexts where sampling problems arise
in computational chemistry.

###  II.1 Electronic ground state methods

A common task in computational chemistry is to determine the ground state of
an electronic Hamiltonian within the Born–Oppenheimer approximation for
$N_{\text{el}}$ electrons and $N_{\text{at}}$ atoms (in atomic units),

$\displaystyle\hat{H}_{\text{el}}$
$\displaystyle=\hat{H}_{0}-\sum_{i}^{N_{\text{el}}}\frac{\nabla_{i}^{2}}{2}-\sum_{i}^{N_{\text{el}}}\sum_{A}^{N_{\text{at}}}\frac{Z_{A}}{|\mathbf{r}_{i}-\mathbf{R}_{A}|}+\sum_{i>j}^{N_{\text{el}}}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|},$
(1)

where $\\{\mathbf{r}_{i}\\}$ are the positions of electrons,
$\\{\mathbf{R}_{A}\\}$ are the positions of atoms, $\\{\mathbf{Z}_{A}\\}$ are
the charges of atoms, and
$H_{0}=\sum_{A<B}\frac{Z_{A}Z_{B}}{|\mathbf{R}_{A}-\mathbf{R}_{B}|}$ is the
nuclear repulsion interaction. The second-quantized form of the one-body and
two-body interactions is

$\hat{H}_{1}=\sum_{pq}h_{pq}a_{p}^{\dagger}a_{q}$ (2)
$\hat{H}_{2}=\frac{1}{2}\sum_{pqrs}v_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{s}a_{r},$
(3)

where $\\{a_{p}^{(\dagger)}\\}$ denotes the fermionic second-quantized
operator, $\hat{H}_{1}$ represents the one-body terms, and $\hat{H}_{2}$
encapsulates the two-body terms. The coefficients $h_{pq}$ and $v_{pqrs}$ are
the one- and two-electron integrals, respectively, defined as:

$h_{pq}=\int\mathrm{d}\mathbf{r}\>\phi_{p}^{*}(\mathbf{r})\left(-\frac{1}{2}\nabla^{2}-\sum_{A}\frac{Z_{A}}{|\mathbf{r}-\mathbf{R}_{A}|}\right)\phi_{q}(\mathbf{r})$
(4) $\displaystyle
v_{pqrs}=\int\int\mathrm{d}\mathbf{r}_{1}\mathrm{d}\mathbf{r}_{2}\>\phi_{p}^{*}(\mathbf{r}_{1})\phi_{q}^{*}(\mathbf{r}_{2})\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\phi_{r}(\mathbf{r}_{1})\phi_{s}(\mathbf{r}_{2}),$
(5)

where $v_{pqrs}=(pr|qs)$ follows chemist’s notation (or the Mulliken
notation), and $\\{\phi_{p}\\}$ is a single-particle basis set chosen a
priori. With this Hamiltonian in mind, one aims to find the ground state
corresponding to its lowest eigenstate (i.e., the ground state energy). We
discuss several classical heuristics that approximately compute the ground
state energy and other properties.

####  II.1.1 Variational Monte Carlo

Variational Monte Carlo (VMC) is a quantum Monte Carlo method commonly used in
both ab initio and model settings [28, 29, 30, 31, 32]. VMC optimizes a set of
variational parameters, $\bm{\theta}$, to minimize the variational energy
expression

$E(\bm{\theta})=\min_{\bm{\theta}}\frac{\langle\Psi({\bm{\theta}})|\hat{H}|\Psi({\bm{\theta}})\rangle}{\langle\Psi({\bm{\theta}})|\Psi({\bm{\theta}})\rangle}.$
(6)

Certain wavefunction ansätze, $|\Psi(\bm{\theta})\rangle$, cannot be
efficiently stored or directly applied with $\hat{H}$ in the computational
basis. Therefore, MC sampling is employed as a practical alternative to
deterministic evaluations of the expectation value of the Hamiltonian. A
stochastic sampling of the variational energy is achieved by inserting a
resolution-of-the-identity in a chosen orthonormal computational basis,
$\\{|n\rangle\\}$,

$\displaystyle E(\bm{\theta})=\sum_{n}E_{L}(\bm{\theta},n)P(\bm{\theta},n),$
(7)

with the local energy

$\displaystyle E_{L}(\bm{\theta},n)$
$\displaystyle=\frac{\langle\Psi({\bm{\theta}})|\hat{H}|n\rangle}{\langle\Psi({\bm{\theta}})|n\rangle}=\frac{\sum_{m}\langle\Psi({\bm{\theta}})|m\rangle\langle
m|\hat{H}|n\rangle}{\langle\Psi({\bm{\theta}})|n\rangle}$ (8)

and the probability distribution

$P(\bm{\theta},n)=\frac{|\langle\Psi({\bm{\theta}})|n\rangle|^{2}}{\langle\Psi({\bm{\theta}})|\Psi({\bm{\theta}})\rangle}.$
(9)

By drawing samples from $P(\bm{\theta},n)$, the variational energy can be
approximated as

$E(\bm{\theta})=\frac{1}{N_{\text{samp}}}\sum_{i}^{N_{\text{samp}}}E_{L}(\bm{\theta},n_{i}).$
(10)

Similarly, one can obtain the gradient of the energy with respect to the
parameters,

$\frac{\partial
E(\bm{\theta})}{\partial\bm{\theta}}=\frac{1}{N_{\text{samp}}}\sum_{i}^{N_{\text{samp}}}\frac{\partial
E_{L}(\bm{\theta},n_{i})}{\partial\bm{\theta}}.$ (11)

Markov chain Monte Carlo samples the normalized distribution in Eq. 9, most
commonly via the Metropolis-Hastings algorithm [33]. We explain details of the
Markov chain Monte Carlo algorithm in the context of finite-temperature
methods in Section II.2.1.

Regarding the cost of VMC, there are two primary considerations for its
efficiency. First, one should consider the cost of sampling $|n\rangle$ from
the distribution, $|\langle\Psi|n\rangle|^{2}$. Within the Metropolis-Hastings
algorithm, this sampling cost arises from two factors. The first is the cost
of evaluating the overlap, $\langle\Psi|n\rangle$, which determines the
transition probability to a new state. Another is the mixing time, which
quantifies how quickly the Markov chain converges to the stationary
distribution. The second consideration is the cost of the local energy
evaluation in Eq. 7. Based on the overlap between $|\Psi\rangle$ and
$|n\rangle$ and the Hamiltonian matrix element in the computational basis,
$\langle m|\hat{H}|n\rangle$, one can evaluate the local energy. One should
account for the sparsity of $\hat{H}$ in the cost model given that the
Hamiltonian $\hat{H}$ in Eq. (1) is sparse.

For a VMC algorithm to be efficient, with costs scaling polynomially with
system size $N$, three key elements must meet this criterion: the overlap
evaluation, the number of non-zero elements in any given column of $\hat{H}$,
and the mixing time of the algorithm. Each of these components should scale as
$\mathcal{O}(\textrm{poly}(N))$ [34]. A more detailed cost analysis can be
provided once specific variational wavefunctions and Hamiltonians are
determined. The goal of quantum algorithms in VMC’s context is to accelerate
these computational aspects.

####  II.1.2 Projector Monte Carlo

Projector Monte Carlo (PMC) methods [27, 35, 36, 37, 38, 39, 40, 41, 26]
stochastically perform imaginary-time evolution (ITE) to sample the ground
state of many-body systems. In the infinite imaginary time limit, one can
obtain the ground state from any initial state $\ket{\Psi_{I}}$ with non-zero
overlap with the ground state $\ket{\Psi_{0}}$,

$\left|\Psi_{0}\right\rangle\propto\lim_{\tau\rightarrow\infty}\exp(-\tau\hat{H})\left|\Psi_{I}\right\rangle.$
(12)

which will be projected to the ground state in the long time limit, as
graphically shown in Fig. 5(a).

Figure 2: The stochastic propagation of walkers $\psi_{i}$ of projector
quantum Monte Carlo methods. The ground state wavefunction $\Psi_{0}$ is
approximated by sampling weight-averaged walkers in the long imaginary time
limit.

A deterministic implementation of imaginary time evolution is generally
exponentially expensive. PMC performs this ground state projection
stochastically to remove this exponential scaling bottleneck.

PMC samples the action of the propagator using an ensemble of statistical
samples with weights, $\\{w_{i}\\}$, and their states,
$\\{|\psi_{i}\rangle\\}$. This is graphically shown in Fig. 2. The global
wavefunction at a given imaginary time $\tau$ is written as

$|\Psi(\tau)\rangle=\sum_{i}w_{i}(\tau)|\psi_{i}(\tau)\rangle.$ (13)

The energy estimate at $\tau$ is

$E(\tau)=\frac{\langle\Psi_{\mathrm{T}}|\hat{H}|\Psi(\tau)\rangle}{\langle\Psi_{\mathrm{T}}|\Psi(\tau)\rangle},$
(14)

where $|\Psi_{\mathrm{T}}\rangle$ is a trial wavefunction used to evaluate the
energy estimate and often to importance-sample (see below).

Different PMC approaches vary in the details of the implementation of a
projection operator in Eq. 12. It is most common to consider a discrete
propagator,

$\ket{\Psi_{0}}=\lim_{n\to\infty}\hat{P}^{n}(\Delta\tau)\ket{\Psi_{I}},$ (15)

where $\Delta\tau$ is the time step to discretize the imaginary time.
Depending on the details of the discrete-time propagator and the computational
basis state, one may use diffusion Monte Carlo (DMC) [42, 43, 44, 45, 46],
Green’s function Monte Carlo [47, 41, 48], full configuration interaction
quantum Monte Carlo [39, 49, 50], and auxiliary field quantum Monte Carlo [35,
36, 51].

Table 1: Classifications of the commonly used PMC approaches [52, 53] Method | Projector | Computational basis | Approximation
---|---|---|---
DMC | $\mathrm{e}^{\Delta\tau(E_{\mathrm{T}}-\hat{H})}$ | first quantization | Fixed-node
GFMC | $1-\Delta\tau(\hat{H}-E_{\mathrm{T}})$ | first/second quantization | Fixed-node
FCIQMC | $1-\Delta\tau(\hat{H}-E_{\mathrm{T}})$ | second quantization | None
AFQMC | $\mathrm{e}^{\Delta\tau(E_{\mathrm{T}}-\hat{H})}$ | second quantization | Phaseless

In nearly all fermionic simulations, except a few special cases [54, 55, 56,
57, 58, 59], PMC methods suffer from the fermionic sign problem. This can be
accredited to the non-stoquasticity of the general fermionic Hamiltonians on a
local basis [60], which generally poses no constraint on the sign of the
walker weights. Thus, the expectation values have a large variance due to many
samples canceling each other to produce a final signal. These cancellations
will lead to a vanishing signal-to-noise ratio of the observable estimates.
Consequently, the computational scaling (i.e., the sample complexity) becomes
exponentially large when one seeks to compute observables up to a certain
additive error. It was proven that solving the sign problem is NP-hard using
classical algorithms [60]. Approximations to PMC methods aim to control the
sign problem and produce an overall polynomial-scaling algorithm. However,
they inevitably introduce biases in the final statistical estimates.
Controlling the sign problem is most commonly done by removing walkers that
cross the node of the trial wavefunction as shown in Fig. 3, resulting in
consistent signs in the weights. Practical PMC methods are classified based on
whether they control the sign problem, the computational basis (first versus
second quantizations), and the form of the projector, as shown in Table 1. We
provide more details of PMC methods based on second quantization as they have
been used in quantum algorithms in recent years.

$\mathcal{S}(D)$$\langle\Psi_{\mathrm{T}}|\phi_{\mathrm{w}}\rangle$ =
0$\phi_{1}(\tau)$$\phi_{1}(\tau+\Delta\tau)$$\phi_{2}(\tau)$$\phi_{2}(\tau+\Delta\tau)$$\phi_{3}(\tau)$$\phi_{3}(\tau+\Delta\tau)$
Figure 3: Propagation of walkers $\phi_{\mathrm{w}}$ in the manifold of Slater
determinants $\mathcal{S}(D)$. Walkers $\phi_{1}$ and $\phi_{2}$ remain in the
allowed hyperspace of similar signed overlap with the trial
$\Psi_{\mathrm{T}}$, whereas walker $\phi_{3}$ will be removed from the
simulation after crossing the overlap node.

##### Full Configuration Interaction Quantum Monte Carlo

Specifically for solving fermionic problems, the full configuration
interaction quantum Monte Carlo (FCIQMC) method offers a memory- and cost-
efficient framework [39, 49, 50]. The imaginary time propagation of an initial
state

$\ket{\Psi_{I}}=\sum_{\mathbf{i}}C_{\mathbf{i}}|\Phi_{\mathbf{i}}\rangle$ (16)

which can be written in terms of the evolution of its coefficients,

$\displaystyle C_{\mathbf{i}}(\tau)$
$\displaystyle=\sum_{\mathbf{j}}(\mathrm{e}^{-\tau
H})_{\mathbf{ij}}C_{\mathbf{j}}(0)=\sum_{\mathbf{j}}\left[\left(\mathrm{e}^{-\Delta\tau
H_{\mathbf{ij}}}\right)^{n}\right]_{\mathbf{ij}}C_{\mathbf{j}}(0)$
$\displaystyle\approx\sum_{\\{\mathbf{j_{1},...,j_{n-1},j}\\}}(1-\Delta\tau
H_{\mathbf{ij_{1}}})...(1-\Delta\tau
H_{\mathbf{j_{n-1}j}})\>C_{\mathbf{j}}(0).$ (17)

It suffices to take the first-order approximation to the imaginary time
propagator, as the eigenspaces, as well as the order of eigenvalues, are
preserved independently of the truncation order. FCIQMC works in a
configuration interaction (CI) basis $\\{|\Phi_{\mathbf{i}}\rangle\\}$,
including all determinants generated by successive excitations from a
reference determinant $|\Phi_{0}\rangle$. Indices $\mathbf{i}$ enumerate
excitation tuples, e.g. a single excitation $(i,a)$, which corresponds to an
excitation from orbital $i$ to orbital $a$. This reference is generally
obtained from a Hartree-Fock calculation. The matrix elements
$H_{\mathbf{ij}}$ are defined in the same CI basis,
$\langle\Phi_{\mathbf{i}}|\hat{H}|\Phi_{\mathbf{j}}\rangle$. During the
propagation of walkers across the CI space via a stochastic application of the
imaginary time propagator, one needs to evaluate the transition matrix
elements, $H_{\mathbf{ij}}$ for each walker. These elements need not be stored
and can be evaluated on the fly.

FCIQMC discretizes the coefficients $C_{\mathbf{i}}(\tau)$ by representing
them as a finite number of walkers $N_{\mathbf{i}}(\tau)$ on determinants
$|\Phi_{\mathbf{i}}\rangle$. The non-unitary nature of FCIQMC may lead to an
exponential growth in the number of walkers
($N_{\mathrm{w}}=\sum_{\mathbf{i}}|N_{\mathbf{i}}|$) as commonly seen in other
PMC methods. To manage this, a dynamical shift $S(\tau)$ of the spectrum of
the Hamiltonian is introduced. This shift is adjusted such that the total
number of walkers remains approximately constant throughout the simulation,
i.e. approximately setting the highest eigenvalue of the propagator to 1. The
walkers are updated according to

$\displaystyle C_{\mathbf{i}}(\tau+\Delta\tau)$
$\displaystyle=C_{\mathbf{i}}(\tau)\underbrace{-\Delta\tau(H_{\mathbf{ii}}-S(\tau))C_{\mathbf{i}}(\tau)}_{\text{Kill/Clone}}$
$\displaystyle\underbrace{-\Delta\tau\sum_{\mathbf{j}\neq\mathbf{i}}H_{\mathbf{ij}}C_{\mathbf{j}}(\tau)}_{\text{Spawn}}.$
(18)

A detailed description of the resulting population dynamics and implementation
can be found in Ref. 39. Each time step one can evaluate the projected energy,

$E(\tau)=\frac{\bra{\Psi_{\mathrm{T}}}H|\Psi(\tau)\rangle}{\bra{\Psi_{\mathrm{T}}}\Psi(\tau)\rangle}=\frac{\sum_{\mathbf{i}}\bra{\Psi_{\mathrm{T}}}H|\Phi_{\mathbf{i}}\rangle
N_{\mathbf{i}}(\tau)}{N_{\mathrm{T}}(\tau)},$ (19)

where $N_{\mathrm{T}}(\tau)$ is the total number of walkers located on
determinants contributing to $\bra{\Psi_{\mathrm{T}}}$. One can show that for
$\tau\to\infty$, the projected energy converges to the exact ground state if
$\ket{\Psi(\tau)}$ faithfully represents the true ground state. This is the
case in the limit of an infinite number of walkers present in the simulation
[39, 49].

While being exact, FCIQMC does not control the sign problem and its
application is limited to relatively small systems due to the exponential-
scaling sample complexity. Some efforts have been put forward to reduce the
sign problem via the initiator approximation, but they have not changed the
asymptotic scaling. One could hope to use quantum computers as a part of
FCIQMC to reduce the sign problem, an approach we will briefly review later in
Section III.2.3.

##### Green’s Function Monte Carlo

Green’s function Monte Carlo (GFMC) was first proposed by Kalos in 1962 [61].
It was then used in solving the ground state of various systems, including
Heisenberg Model [62], lattice fermions [63], molecular systems [64], and
bosonic systems such as Helium [65, 66]. In GFMC, the ground state is also
reached by performing imaginary time evolution successively, similar to all
other PMC methods. Suppose $\ket{\Psi_{I}}$ is the initial wavefunction that
has a non-zero overlap with the exact ground state. We define

$\ket{\Psi^{(n)}}:=(1-\Delta\tau(\hat{H}-E_{\mathrm{T}}))^{n}\ket{\Psi_{I}},$
(20)

and then we expand the projection operator in a complete basis set
$\\{\ket{x}\\}$, which yields

$\Psi^{(n+1)}(x^{\prime})=\sum_{x}\mathcal{G}_{x^{\prime},x}\Psi^{(n)}(x)$
(21)

where $\Psi(x)=\langle x|\Psi\rangle$, and $\mathcal{G}_{x^{\prime},x}$ is the
Green’s function defined as

$\mathcal{G}_{x^{\prime},x}=\langle
x^{\prime}|(1-\Delta\tau(\hat{H}-E_{\mathrm{T}}))|x\rangle.$ (22)

It is then natural to view Eq. 21 as a master equation with $\Psi^{(n)}(x)$
being the probability distribution and $\mathcal{G}_{x^{\prime},x}$ being the
transition probability if all of the entries of the Green’s function and the
wavefunction are non-negative. However, this is not always the case. If a
system under some local basis cannot guarantee the positivity of all the
matrix elements in $\mathcal{G}_{x^{\prime},x}$, then the sign problem will
arise. For simplicity, we will only consider systems that do not exhibit sign
problems in this section. We further decompose Green’s function into the
normalization factor and the stochastic matrix describing the transition
probability:

$\mathcal{G}_{x^{\prime},x}=b_{x}p_{x^{\prime},x}$ (23)

where $b_{x}$ is the normalization factor given by

$b_{x}=\sum_{x^{\prime}}\mathcal{G}_{x^{\prime},x}.$ (24)

For each configuration $x^{(n)}$ at time step $n$, we assign a weight
$w^{(n)}$ to take account of the normalization factor, and hence the update
rule of the weight is given by

$w^{(n+1)}=b_{x^{(n)}}w^{(n)}.$ (25)

Then, after a sufficiently long time, the equilibrium probability distribution
over $\\{\ket{x}\\}$ is the ground state expanded on this local basis.

For most of the fermionic Hamiltonians, there is no guarantee that the entries
of Green’s function are all positive. Hence, we need to control the sign
problem by using fixed-node approximation. The fixed node approximation
removes the sign problem in sampling but introduces a bias. Recently, progress
has been made in enhancing GFMC with quantum computers, thereby reducing the
bias introduced by the fixed-node approximation [67]. Interested readers can
refer to Ref. 27 for further details on the theory of GFMC.

##### Auxiliary Field Quantum Monte Carlo

In the past decades, AFQMC [35, 36] has established itself as an exceptionally
accurate and efficient method for tackling many-body systems. In recent
developments, AFQMC has proven to be a powerful tool for solving ab initio
problems, with extensive use across a broad spectrum of quantum chemical
problems [51, 68, 69, 70, 71, 72, 73, 74, 75, 76]. AFQMC performs an open-
ended random walk to propagate an initial representation of the wave function
$\ket{\Psi_{I}}$ in imaginary time,

$\ket{\Psi_{0}}=\lim_{N\to\infty}\left(\mathrm{e}^{-\Delta\tau(\hat{H}-E_{\mathrm{T}})}\right)^{N}\ket{\Psi_{I}}.$
(26)

Applying the Trotter decomposition [77] and the subsequent
Hubbard–Stratonovich transformation [78, 79], we obtain a formula of the
short-time propagator expressed as an integral over auxiliary fields

$\mathrm{e}^{-\Delta\tau\hat{H}}=\int\mathrm{d}^{N_{\text{aux}}}\mathbf{x}\
p(\mathbf{x})\hat{B}(\mathbf{x})+\mathcal{O}(\Delta\tau^{2}),$ (27)

where $\mathbf{x}$ is a vector of $N_{\text{aux}}$ auxiliary fields subjected
to Gaussian distribution, $p(\mathbf{x})$ is the standard normal distribution
and $\hat{B}(\mathbf{x})$ is a one-body propagator. The integral is carried
out via Monte Carlo by writing the global wavefunction at each time step $n$
as a weighted sum of $N_{\text{w}}$ so-called random walkers,

$\ket{\Psi^{(n)}}=\sum_{i=1}^{N_{\text{w}}}w_{i}^{(n)}\frac{\ket{\phi^{(n)}_{i}}}{\langle\Psi_{\mathrm{T}}|\phi^{(n)}_{i}\rangle},$
(28)

where $\ket{\phi_{i}^{(n)}}$ are single determinant states. The weights of the
walkers are updated according to the overlap ratio

$w_{i}^{(n+1)}=w_{i}^{(n)}\times\frac{\langle\Psi_{\mathrm{T}}|\hat{B}(\mathbf{x}_{i}^{(n)})|\phi^{(n)}_{i}\rangle}{\langle\Psi_{\mathrm{T}}|\phi^{(n)}_{i}\rangle}.$
(29)

Following this, the weights are complex numbers in general, and the average of
any observable will suffer from an exponentially large variance, known as the
phase problem [51]. One could control the phase problem by imposing
constraints on the weights of the walkers. The phaseless approximation [36] is
employed in practical AFQMC calculations to achieve this. It ensures the
positivity of the weights during propagation (also depicted in Fig. 3) by
applying

$w_{i}^{(n+1)}=w_{i}^{(n)}\times\left|\frac{\langle\Psi_{\mathrm{T}}|\hat{B}(\mathbf{x}_{i}^{(n)})|\phi^{(n)}_{i}\rangle}{\langle\Psi_{\mathrm{T}}|\phi^{(n)}_{i}\rangle}\right|\times\max\\{0,\cos\theta_{i}^{(n)}\\},$
(30)

where $\theta_{i}^{(n)}$ is the argument of the complex overlap ratio. This
approximation introduces biases to the calculation, however, one can mitigate
the bias by employing more accurate trial wavefunctions [80, 81, 72, 52, 75].
In the limit of the trial wavefunction being the exact ground state of the
Hamiltonian, the bias can be completely removed. The estimate of the ground
state energy is given by the mixed estimator

$E(n\Delta\tau)=\frac{\sum_{i}w_{i}^{(n)}E_{L,i}^{(n)}}{\sum_{i}w_{i}^{(n)}}$
(31)

where $E_{L,i}^{(n)}$ denotes the local energy of the $i$-th walker at time
step $n$,

$E_{L,i}^{(n)}=\frac{\langle\Psi_{\mathrm{T}}|\hat{H}|\phi_{i}^{(n)}\rangle}{\langle\Psi_{\mathrm{T}}|\phi_{i}^{(n)}\rangle}.$
(32)

Similarly to GFMC, one could hope to use a more sophisticated trial
wavefunction implemented through quantum computers to reduce the phaseless
bias. This is the central intuition behind the hybrid quantum-classical AFQMC
that we will discuss later in Section III.2.1 We only briefly reviewed the
AFQMC method as needed to understand its quantum counterpart. For a detailed
review of AFQMC, we refer the readers to Ref 51.

####  II.1.3 Coupled Cluster Monte Carlo

The coupled cluster (CC) wavefunction is parametrized by a cluster operator
$\hat{T}=\sum_{\mathbf{i}\in\mathcal{S}}t_{\mathbf{i}}\hat{E}_{\mathbf{i}}$,
giving

$|\Psi_{\mathrm{CC}}\rangle=\mathrm{exp}(\hat{T})|\Phi_{\mathrm{ref}}\rangle,$
(33)

where $\mathcal{S}$ defines the set of excitation operators,
$\hat{E}_{\mathbf{i}}$ is the $i$-th excitation operator that creates an
excited determinant via
$\hat{E}_{\mathbf{i}}|\Phi_{\text{ref}}\rangle=\ket{\Phi_{\mathbf{i}}}$, and
$t_{\mathbf{i}}$ is the cluster amplitude for the $i$-th excitation. Here,
$|\Phi_{\text{ref}}\rangle$ is some reference determinant, usually chosen from
a previous Hartree-Fock calculation. This exponential parametrization can
express the exact solution if all possible excitations are included in the
definition of $\hat{T}$. In practice, however, the cluster operator is
truncated after a certain excitation order to keep computations tractable.

Similarly to FCIQMC, coupled cluster Monte Carlo (CCMC)[82, 83] utilizes a CI
expansion to stochastically solve the CC projection,

$\displaystyle\langle\Phi_{\mathbf{i}}|\hat{H}-E|\Psi_{\mathrm{CC}}\rangle=0,$
(34)

for all excitations $\mathbf{i}$. As a result of the projection equations in
Eq. (34) one can write a first order imaginary-time propagator,

$\langle\Phi_{\mathbf{i}}|1-\tau(\hat{H}-E)|\Psi_{\mathrm{CC}}\rangle=\langle\Phi_{\mathbf{i}}|\Psi_{\mathrm{CC}}\rangle.$
(35)

Unlike FCIQMC, the propagation prescribed by Eq. (35) is not straightforward
due to non-linearity, as generally

$\langle\Phi_{\mathbf{i}}|\Psi_{\mathrm{CC}}\rangle=t_{\mathbf{i}}+\mathcal{O}(t^{2}).$
(36)

This can be seen when considering a double excitation $\mathbf{i}=(i,j,a,b)$,

$\langle\Phi_{\mathbf{i}}|\Psi_{\mathrm{CC}}\rangle=t^{ab}_{ij}+t^{a}_{i}t^{b}_{j}+t^{a}_{j}t^{b}_{i}.$
(37)

The projection of the CC wave function onto a doubly excited determinant has
contributions from products of cluster amplitudes. These products correspond
to composite excitors $\hat{E}_{\mathbf{i}}\hat{E}_{\mathbf{j}}$ in the
cluster operator. From Eq. (35), it is not clear how to formulate population
dynamics, including composite clusters. The CCMC evolution circumvents this
problem by requiring to spawn on non-composite clusters. Hence, one only needs
to store walkers on these clusters instead of storing walkers on all possible
excitors and their products. This, however, also limits CCMC to systems where
composite cluster amplitudes are expected to be sufficiently small such that
Eq. (36) is a valid approximation. The CCMC iterative procedure is then
captured by

$\displaystyle t_{\mathbf{i}}(\tau+\Delta\tau)$
$\displaystyle=t_{\mathbf{i}}(\tau)-\Delta\tau\sum_{\mathbf{j}}H_{\mathbf{i}\mathbf{j}}\langle\Phi_{\mathbf{j}}|\Psi(\tau)\rangle.$
(38)

Unfortunately, CCMC suffers from the sign problem [83], which is closely
related to that in FCIQMC. Due to the exponential-scaling sample complexity,
it is not always clear whether there is a significant computational benefit to
use CCMC over standard deterministic CC algorithms. Extending this method, one
could hope to sample a different variant of CC wavefunctions (called unitary
CC) using quantum computers or reduce the sign problem with quantum computers.
We will review an attempt along this line in Section III.3.

####  II.1.4 Selected Configuration Interaction

Configuration interaction (CI) is one of the conceptually simplest post-HF
methods for calculating correlation energy. The basic idea is to diagonalize
the Hamiltonian in a selected basis of $N_{\textrm{el}}$-electron Slater
determinants, or in other words, construct a ground state wavefunction as a
linear combination of Slater determinants with the coefficients determined
variationally [84]. Similarly to CC,

$|\Psi_{\text{CI}}\rangle=\hat{C}|\Phi_{\text{ref}}\rangle,$ (39)

where
$\hat{C}=\sum_{\mathbf{i}\in\mathcal{S}}c_{\mathbf{i}}\hat{E}_{\mathbf{i}}$.
If one defines the set $\mathcal{S}$ to include all possible excitations, one
recovers the exact FCI approach with an exponential scaling cost. Therefore,
the most commonly seen polynomial scaling approximations truncate the space of
determinants based on excitation level. These approximations include well-
known methods such as CI-singles (CIS), CI-singles and doubles (CISD), etc. It
is worth noting that truncated CI methods are well known to lack size
consistency, limiting their usefulness [1], unlike the CC wavefunction with
the same truncation.

Going beyond the excitation-based truncation in CI, selected CI (SCI) has seen
revived interests in quantum chemitry [85, 86, 87, 88, 89, 90]. In essence,
SCI methods define the excitation subspace $\mathcal{S}$ via a sampling
protocol (not necessarily based on Monte Carlo) intending to find important
determinants for electron correlation. One of the earliest SCI methods was
developed roughly 50 years ago, named configuration interaction perturbatively
selected iteratively (CIPSI), where an iterative procedure based on a
perturbative estimate is used to improve the choice of subspace [91]. CIPSI
has been shown to yield near-exact energies for systems much beyond the reach
of the limit of FCI [92, 93, 76]. Furthermore, new selection schemes exploring
the tradeoff between the selection efficiency and accuracy have been proposed,
example being adaptive sampling configuration interaction (ASCI) [94, 95],
adaptive configuration interaction (ACI)[96, 97], heat-bath configuration
interaction (HCI) [98, 99, 100], and iterative CI (iCI) [101, 102, 103, 104].

A Monte Carlo alternative to these SCI methods is the Monte Carlo
configuration interaction (MCCI) [89, 90]. In the MCCI approach, the FCI space
is divided into three parts: the trial CI vector, CI vectors that interact
with the trial (only allow single or double excitations from the trial CI),
and those that do not interact. Each cycle commences with branching processes
where random single and double excitations are uniformly sampled to generate
new vectors, thus forming a defined subspace. Within this subspace, the
Hamiltonian matrix is diagonalized to obtain the lowest eigenstate, which then
serves as the trial CI vector for the next cycle. Following diagonalization, a
selection mechanism evaluates the components of the vector to decide which
should be retained or discarded based on their contributions to the system’s
energy and wavefunction. The importance of a state is determined by the
magnitude of the energy reduction estimated through second-order perturbation
theory. This procedure is iteratively repeated until the vector’s expansion
stabilizes in length or a predefined energy criterion is met. Typically, the
Hartree-Fock wavefunction is employed as the initial trial CI vector. The
quantum version of MCCI, further discussed in Section III.5, extends these
principles to quantum computational frameworks.

Despite significant progress in SCI-related methods and impressive results,
these methods all scale exponentially with system size due to the
(exponentially) growing number of important determinants. This ultimately
limits the scope of these methods to small systems. Nevertheless, one could
hope to sample relevant determinants more efficiently using quantum computers,
which is an idea we will review later in Section III.5.

###  II.2 Finite-temperature methods

In addition to the ground state calculations in the previous section,
understanding finite-temperature properties is also common in computational
chemistry and materials science, as they provide insights into the equilibrium
properties and responses to external perturbations of systems at finite
temperatures. The classical Monte Carlo framework gives rise to a multitude of
finite temperature methods such as finite temperature AFQMC [105, 106, 107,
70], density matrix QMC [108], path-integral MC [109, 110, 111, 112],
stochastic series expansion [113, 114, 115, 116], thermo-field Monte Carlo
[117], etc. Not all of these methods have a quantum algorithm counterpart, so
we will only review those with the quantum counterpart.

####  II.2.1 Classical Monte Carlo

Statistical physics provides the theoretical framework to study the thermal
properties of macroscopic systems composed of many microscopic local degrees
of freedom. The Ising model serves as a fundamental model for exploring phase
transitions in magnetic systems. The demonstration of the order-to-disorder
phase transition within the two-dimensional and three-dimensional Ising model
with Monte Carlo methods is considered one of the significant achievements of
20th-century statistical physics [118, 119]. The studies profoundly influenced
understanding critical phenomena and the nature of phase changes in various
physical systems. Most importantly, such studies are commonly performed by
Monte Carlo sampling.

We briefly recapitulate the Markov chain Monte Carlo (MCMC) with the classical
Ising model as an example. The classical Ising model consists of discrete
magnetic spins on a lattice that can be in one of two states (up or down) and
interact with their nearest neighbors. The Hamiltonian gives the energy of a
particular state of the system:

$E(\textbf{s})=-\sum^{n}_{j<k}J_{jk}s_{j}s_{k}-\sum_{j}^{n}h_{j}s_{j}$ (40)

where each spin configuration $\textbf{s}\in\\{-1,+1\\}^{n}$. The central
object of interest that allows us to study the properties at varying
temperatures is the canonical partition function,

$Z=\sum_{\mathbf{s}}\mathrm{e}^{-\beta E(\mathbf{s})}$ (41)

at temperature $1/\beta$. Systems that obey Boltzmann statistics are
ubiquitous in physics, and sampling these distributions efficiently is an
important subroutine in several algorithms.

One could design a Markov chain that converges to a stationary point of the
Boltzmann distribution. To achieve the desired stationary condition, the
Markov chain must be irreducible and aperiodic. A Markov chain can be defined
by a transition matrix $P$, where
$P_{\textbf{s's}}=P(\textbf{s}^{\prime}|\textbf{s})$. A single step on this
Markov chain from a current state distribution $p$ can be expressed as
$Pp=p^{\prime}$. The detailed balance condition (a form of Bayes’ theorem) of
the Markov chain needs to be satisfied to converge to the desired stationary
distribution $\pi$:

$P(\textbf{s}^{\prime}|\textbf{s})\pi(\textbf{s})=P(\textbf{s}|\textbf{s}^{\prime})\pi(\textbf{s}^{\prime})$
(42)
$\sum_{s}P(\textbf{s}^{\prime}|\textbf{s})\pi(\textbf{s})=\pi(\textbf{s}^{\prime})$
(43)

with the target probability being the Boltzmann distribution,

$\pi(\mathbf{s})=\frac{\mathrm{e}^{-\beta E(\mathbf{s})}}{Z}.$ (44)

If we know a function $f(\textbf{s})\propto\pi(\textbf{s})$, the most common
way to construct a Markov chain that converges to $\pi$ is through an
acceptance-rejection based algorithm. In the proposal step, from the current
state s, the next state s’ is proposed with probability
$Q(\textbf{s}^{\prime}|\textbf{s})$. In the accept/reject step, the proposed
step is accepted with probability $A(\textbf{s}^{\prime}|\textbf{s})$. The
samples stay in the same configuration if the proposal is rejected. There are
multiple ways of choosing $A$: the most popular are the Metropolis-Hastings
algorithm, similar to that of VMC, and the Gibbs sampling algorithm. The
Metropolis–Hastings choice for $A$ is given by

$A(\textbf{s}^{\prime}|\textbf{s})=\min\left(1,\frac{\pi(\textbf{s}^{\prime})}{\pi(\textbf{s})}\frac{Q(\textbf{s}|\textbf{s}^{\prime})}{Q(\textbf{s}^{\prime}|\textbf{s})}\right).$
(45)

The mixing time, defined as the expected number of iterations required for the
Markov chain’s distribution to get within $\epsilon$ of the true distribution
in total variation distance is given by

$t_{\epsilon}\leq\delta^{-1}\ln\left(\frac{1}{\epsilon\min\pi(\textbf{s})}\right),$
(46)

where $\delta=1-\max_{|\lambda|<1}|\lambda|$ and $\lambda$ is an eigenvalue of
the transition matrix $P$ [33].

One potential problem in classical MCMC is that the ratio
$\pi(\textbf{s}^{\prime})/\pi(\textbf{s})$ can be very small, and moves are
accepted with very low probability, resulting in a significantly long
convergence time. Many efforts in quantum counter parts of the MCMC algorithms
are targeted toward accelerating the convergence time of MCMC. We will
introduce the enhancement of MCMC using quantum computers in Section IV.1.2
and the quantum version of Metropolis in Section IV.1.3 by utilizing the
quantum phase estimation.

####  II.2.2 Lindbladian dynamics

The Markov chain in Eq. (42) can also be expressed as [120]

$P_{s^{\prime}s}=\underbrace{\sum_{a\in
A}p(a)\gamma\left(E(\textbf{s}^{\prime})-E(\textbf{s})\right)\bm{A}_{\textbf{s}^{\prime}\textbf{s}}^{a}}_{\text{"Accept"
}}+\underbrace{\vphantom{\sum_{a\in
A}}\bm{R}_{\textbf{s}^{\prime}\textbf{s}}\delta_{\textbf{s}^{\prime}\textbf{s}}}_{\text{"Reject"
}},$ (47)

where $\gamma(\omega):=\min\left(1,\mathrm{e}^{-\beta\omega}\right)$ and
$\omega$ is the energy difference of the transition.
$\bm{A}_{\textbf{s}^{\prime}\textbf{s}}^{a}$ are the stochastic matrices
corresponding to the transitions (e.g., the spin-flip in Ising model), as
weighted by the Metropolis factor $\gamma(\omega)$. The rejection part
$\bm{R}_{\textbf{s}^{\prime}\textbf{s}}$ is a diagonal matrix determined by
the probability preserving constraints. A continuous-time Markov chain
generator can be defined as

$\displaystyle\bm{L}_{s^{\prime}s}$ $\displaystyle=\sum_{a\in
A}p(a)(\underbrace{\vphantom{\sum_{a\in
A}}\gamma\left(E(\textbf{s}^{\prime})-E(\textbf{s})\right)\bm{A}_{\textbf{s}^{\prime}\textbf{s}}^{a}}_{\text{"transition"
}}$ (48)
$\displaystyle-\underbrace{\delta_{\textbf{s}^{\prime}\textbf{s}}\sum_{\textbf{s''}}\gamma\left(E(\textbf{s''})-E(\textbf{s})\right)\bm{A}_{\textbf{s''}\textbf{s}}^{a}}_{\text{"decay"
}})$

In section IV.2, we will introduce the quantum analog of MCMC for sampling the
thermal states of any many-body Hamiltonian $\hat{H}$ using quantum computer
with a similar form as Eq. (48) while construct a Lindbladian,

$\mathrm{e}^{\mathcal{L}_{\beta}t}\left[\bm{\rho}_{\beta}\right]=\bm{\rho}_{\beta}\quad\text{
where
}\quad\bm{\rho}_{\beta}=\mathrm{e}^{-\beta\bm{H}}/\operatorname{Tr}\left(\mathrm{e}^{-\beta\bm{H}}\right)$
(49)
$\mathcal{L}_{\beta}[\bm{\rho}]=\underbrace{-\mathrm{i}[\bm{B},\bm{\rho}]}_{\text{"coherent"
}}+(\text{ "dissipative") }$ (50)

Lindbladians have been combined with Monte Carlo sampling in the context of
both ground and thermal state preparation, as well as real-time dynamics, as a
branch of full quantum algorithms, which we will cover in Section III.6 and
Section IV.2.

####  II.2.3 Stochastic Series Expansion

Another method to evaluate thermal averages of many-body systems is the
stochastic series expansion (SSE) [113, 114]. In essence, SSE samples the
trace of the $L$-th order Taylor expansion of $\exp{(-\beta\hat{H})}$. This
idea goes back to the 1960s when it was first introduced by Handscomb [115].
Unlike most Monte Carlo schemes, SSE introduces no systematic errors stemming
from the artificial time-discretization [121]. While SSE has been primarily
applied to spin chains [122, 123, 124, 125], its framework is universal and
may, in principle, be applied to generic quantum systems. However,
applications to such systems are challenging due to the severity of the sign
problem in the method [126].

Considering an $N$-body Hamiltonian $\hat{H}=-\sum_{i=1}^{N}\hat{H}_{b_{i}}$,
the partition function can be computed by

$\displaystyle Z$
$\displaystyle=\mathrm{Tr}(\mathrm{e}^{-\beta\hat{H}})=\sum_{a}\langle
a|\mathrm{e}^{-\beta\hat{H}}|a\rangle$
$\displaystyle=\sum_{a}\sum_{n}\frac{(-\beta)^{n}}{n!}\langle
a|\left(-\sum_{i=1}^{N}\hat{H}_{b_{i}}\right)^{n}|a\rangle$
$\displaystyle\approx\sum_{a}\sum_{n}^{L}\sum_{\\{b_{i}\\}}\frac{\beta^{n}}{n!}\langle
a|\hat{H}_{b_{n}}...\hat{H}_{b_{1}}|a\rangle.$ (51)

SSE generates a random walk in a joint state- and index-space by sampling
tuples $(a,n,\\{b_{i}\\})$ and evaluating corresponding matrix elements
contributing to $Z$. This sampling procedure is commonly implemented via the
Metropolis-Hastings algorithm, ultimately sampling relative weights adhering
to the probability distribution

$\displaystyle\pi(a,n,\\{b_{i}\\})$ $\displaystyle=\frac{\beta^{n}}{n!}\langle
a|\hat{H}_{b_{n}}\cdots\hat{H}_{b_{1}}|a\rangle$ (52) $\displaystyle\geq 0,$
(53)

where the inequality is a requirement for $\pi(a,n,\\{b_{i}\\})$ to define a
probability distribution. This condition is straightforwardly satisfied by
shifting the system Hamiltonian by a sufficiently large constant value [127].
For large $N$ this procedure is computationally tractable only if one imposes
an additional _no-branching_ requirement,

$\forall\>a,b\>\exists\>a^{\prime}:\hat{H}_{b}|a\rangle\propto|a^{\prime}\rangle,$
(54)

alleviating the exponential cost of evaluating and storing all diagonal and
off-diagonal operator string matrix elements. This condition limits the choice
of basis to only a small subset depending on the system Hamiltonian.

With the no-branching requirement, one can no longer shift the spectrum of the
system arbitrarily, giving rise to samples with arbitrary signs. This sign
problem directly manifests in,

$\displaystyle\langle\hat{A}\rangle_{\pi}$
$\displaystyle=\frac{\sum_{C}\pi(C)A(C)}{\sum_{C}\pi(C)}$
$\displaystyle=\frac{\sum_{C}\mathrm{sgn}(\pi(C))|\pi(C)|A(C)}{\sum_{C}\mathrm{sgn}(\pi(C))|\pi(C)|}$
$\displaystyle=\frac{\langle\mathrm{sgn}(\pi(C))A(C)\rangle_{|\pi|}}{\langle\mathrm{sgn}(\pi(C))\rangle_{|\pi|}},$
(55)

where we combined all summation indices into a configuration index,
$C=(a,n,\\{b_{i}\\})$. The average sign
$\langle\mathrm{sgn}(\pi(C))\rangle_{|\pi|}$ generally decays exponentially
with inverse temperature $\beta$ and system size [128], leading to a severe
sign problem for large systems and systems at low temperatures. Because of the
sign problem, SSE simulations generally require exponential computational
costs, making them unfeasible for more complicated systems. As in Sec. IV.3,
quantum computers can naturally represent superpositions of arbitrary basis
states, allowing for dispensing with the no-branching condition.

####  II.2.4 Minimally Entangled Typical Thermal States

Another approach to compute the thermal density matrices and properties is by
sampling the minimally entangled typical thermal states (METTS) [129, 130]
obtained from the ITE of product states. METTS is defined by

$|\phi(i)\rangle=P(i)^{-1/2}\mathrm{e}^{-\beta\hat{H}/2}|i\rangle$ (56)

with the unnormalized probability

$P(i)=\langle i|\exp(-\beta\hat{H})|i\rangle.$ (57)

$|\phi(i)\rangle$ is expected to have minimal entanglement entropy if the
typical states $|i\rangle$ are chosen as non-interacting classical product
states (CPS). The von Neumann entanglement entropy of CPS after bipartite
decomposition is strictly zero,

$|i\rangle=\left|i_{1}i_{2}i_{3}\ldots i_{N}\right\rangle$ (58)

where $|i_{j}\rangle$ is the arbitrarily chosen local basis for site $j$.

Sampling with probability $P(i)/Z$ where $Z=\text{Tr}{e^{-\beta\hat{H}}}$, one
can compute the thermal density matrix and observable expectation from an
ensemble of the pure states,

$\rho_{\beta}=Z^{-1}\sum_{i}P(i)|\phi(i)\rangle\langle\phi(i)|$ (59)
$\langle\hat{A}\rangle=Z^{-1}\sum_{i}P(i)\langle\phi(i)|\hat{A}|\phi(i)\rangle$
(60)

The rejection-free Markov chain sampling that propses moves from one CPS
$|i\rangle$ to another $|j\rangle$ is similar to the Monte Carlo step
described in Section II.2.1. The Markov chain is based on the transition
probability $p\left(i\rightarrow j\right)=\left|\left\langle
j|\phi(i)\right\rangle\right|^{2}$, leading to the detailed balance

$\frac{P(i)}{Z}p(i\rightarrow j)=\frac{P(j)}{Z}p(j\rightarrow i),$ (61)

recovering the correct thermal distribution. This approach is commonly used
with matrix product state (MPS) based methods to carry out the imaginary time
evolution and sample from the transition probability.

As the MPS ansatz is primarily effective for low-dimensional systems thanks to
the area law [131, *PhysRevLett.90.227902, *RevModPhys.82.277,
*orus2014practical], making extensions to 2D and 3D quantum many-body systems
challenging in terms of both storage and ITE demands. One may hope to lift
these limitations using quantum computers within METTS, which we will review
in Section IV.4.

###  II.3 Real-time dynamics

The non-equilibrium dynamics of quantum systems is studied by solving the
time-dependent Schrödinger’s equation:

$\mathrm{i}\partial_{t}|\psi(t)\rangle=\hat{H}(t)|\psi(t)\rangle.$ (62)

If one is interested in pure-state dynamics at $T=0$, this equation of motion
contains all the necessary dynamical information. For instance, the dynamics
of a local observable, $\hat{O}(t)$, is given by

$\langle O(t)\rangle=\langle\psi(t)|\hat{H}(t)|\psi(t)\rangle,$ (63)

where $|\psi(t)\rangle$ is determined by Eq. 62. To approximately solve this
problem, “real-time” wavefunction methods have been developed, such as real-
time time-dependent density functional theory [135], real-time coupled cluster
methods [136], etc.

In another context, one may consider a bipartite Hamiltonian,

$\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{H}_{A-B}$ (64)

where $\hat{H}_{\text{A}}$, $\hat{H}_{\text{B}}$, and $\hat{H}_{A-B}$ are
Hamiltonians for system $A$, system $B$, and the interaction between the two,
respectively. In these problems, one is only interested in the dynamics of
degrees of freedom of $A$, which can be completely characterized by a reduced
density operator,

$\hat{\rho}_{A}=\text{Tr}_{B}(|\psi(t)\rangle\langle\psi(t)|\mathrm{e}^{-\beta\hat{H}})$
(65)

where $\beta$ is the inverse temperature. In this case, one is also often
interested in thermal correlation functions given as

$\langle
A(0)B(t)\rangle=Z^{-1}\text{Tr}(\hat{A}(0)\hat{B}(t)\exp(-\beta\hat{H}))$ (66)

where $\hat{A}$ and $\hat{B}$ act only on the system $A$. These computational
tasks are usually considered part of open quantum system dynamics or impurity
problems but can also be considered in mixed quantum-classical dynamics with
additional approximations.

Directly calculating dynamical properties using Monte Carlo sampling is more
challenging than the computation of the static properties. The additional
challenge is due to so-called the dynamical sign problem, which becomes
apparent after inspecting the path integral sampling of the real-time
propagator, $\exp(\mathrm{i}\hat{H}t)$. There is an indirect way to address
this problem via analytic continuation, where one takes imaginary-time
correlation functions and analytically continues to the real-time domain.
While the recent advances made are quite encouraging [137, 138, 139, 140, 141,
142, 143], analytic continuation is a formally ill-conditioned inverse problem
that could miss some of the key features in the target real-time observables.
This section will focus on sampling methods that work directly in real time.

####  II.3.1 Classical and Mixed Quantum–Classical Methods

For a broad range of problems found in computational chemistry, treating the
entire system classically or a smaller part of it quantum mechanically while
handling the rest classically can be well-suited. Prominent examples include
studying the dynamics and structure of liquid and chemical reactions in the
condensed phase. In essence, one can sample from an equilibrium distribution,

$Z=\text{Tr}(\mathrm{e}^{-\beta\hat{H}})=\int\mathrm{d}\mathbf{z}\exp(-\beta
H(\mathbf{z}))$ (67)

where $\mathbf{z}$ is a classical variable that completely determines the
system’s energy. One could also perform the dynamics of these classical
variables following Newton’s equations to calculate thermal equilibrium
correlation functions. This approach is known as the molecular dynamics
approach [144].

One can also partition the sampling into quantum (system $A$) and classical
variables (system $B$),

$\displaystyle Z$ $\displaystyle=\text{Tr}(\mathrm{e}^{-\beta\hat{H}})$ (68)
$\displaystyle=\int\mathrm{d}\mathbf{z}\exp(-\beta
H_{B}(\mathbf{z}))\text{Tr}_{A}(\mathrm{e}^{-\beta(\hat{H}_{A}+\hat{H}_{A-B}(\mathbf{z}))}).$

The Born–Oppenheimer ab initio dynamics fit into this form where $A$ is the
electronic and $B$ is the nuclear degrees of freedom, respectively. Suppose
the electronic energy levels of $A$ are close, so multiple states contribute
to the trace. In that case, one may consider non-adiabatic effects and use
mixed quantum-classical approaches such as Ehrenfest dynamics [145] and
surface hopping [146].

Computational bottlenecks in this area include the cost of evaluating the
forces on the nuclei and the duration of the molecular dynamics simulations
necessary to probe the physical and chemical questions of interest. One may
hope that quantum computers may be able to improve the accuracy or efficiency
of force evaluation (i.e., the same goal as the ground state algorithms) [19,
147] or the efficiency of sampling the relevant regions of the phase space
such that the simulation time can be shortened (similar to what enhanced
sampling approaches accomplish.)

####  II.3.2 Diagrammatic Monte Carlo

A relatively widely used QMC method for both finite temperature (imaginary-
time) and real-time problems is diagrammatic Monte Carlo (DiagMC). DiagMC
assumes a random variable $y$ to be distributed according to

$Q(y)=\sum_{m=0}^{\infty}\sum_{\xi_{m}}\int\cdots\int\mathrm{d}x_{m}\cdots\mathrm{d}x_{1}F(\xi_{m},y,x_{1},\cdots,x_{m}).$
(69)

This random variable is then sampled via a Metropolis-Hastings algorithm from
a configuration space of diagram order $m$, topology $\xi_{m}$, and internal
variables $\\{x_{i}\\}$. The algorithm is similar to sampling from
multidimensional integrals but additionally accounts for the variable
integration multiplicity of individual terms of the sum [148, 149]. Eq. 69 can
be used to sample time-dependent observables expanded in the Dyson series
[150]:

$\displaystyle\langle\hat{O}(t)\rangle$
$\displaystyle=\sum_{m=0}^{\infty}\sum_{n=0}^{m}\int\mathrm{d}s_{m}\cdots\int\mathrm{d}s_{1}(-1)^{n}\mathrm{i}^{m}$
(70)
$\displaystyle\times\left\langle\hat{H}^{\prime}(s_{m})\cdots\hat{H}^{\prime}(s_{n+1})\hat{O}^{\prime}(t_{\text{max}})\hat{H}^{\prime}(s_{n})\cdots\hat{H}^{\prime}(s_{1})\right\rangle$

where $\hat{A}^{\prime}(t)$ denotes the time-dependent operator in the
interaction picture
$\hat{A}^{\prime}(t)=\mathrm{e}^{\mathrm{i}\hat{H}_{0}t}\hat{A}\mathrm{e}^{-\mathrm{i}\hat{H}_{0}t}$,
and $s_{i}$ is the time index ordered according to the Keldysh contour
causality. Due to the $(-1)^{n}\mathrm{i}^{m}$ factor in the expansion, the
average sign of the diagrams decays exponentially as the expansion order
increases, which is termed as the dynamical sign problem [151].

Recently, Cohen et al. [152, 153] proposed an inchworm modification to DiagMC,
effectively dealing with the dynamical sign-problem and enabling real-time
propagation of the Anderson impurity model with computational cost scales
quadratically with time. To propagate for $t_{i}$ to $t_{f}$ one starts from
an exact short-time propagator for a small time-segment
$[t_{i},t_{\scriptscriptstyle\uparrow}]$ and
$t_{\scriptscriptstyle\uparrow}<t_{f}$. The algorithms subsequently inches
forward by $\Delta\tau$ to obtain an approximately exact propagator for the
interval $[t_{i},t_{\scriptscriptstyle\uparrow}+\Delta\tau]$ by sampling
diagrams either completely contained in
$[t_{\scriptscriptstyle\uparrow},t_{f}]$, or connecting the previous
approximation of the exact propagator on
$[t_{i},t_{\scriptscriptstyle\uparrow}]$ to the interval
$[t_{\scriptscriptstyle\uparrow},t_{f}]$. While the results of the Anderson
impurity model are impressive, more general problems are still prohibitive to
handle due to the dynamical sign problem. More details can be found in Ref.
152. One may hope to reduce the dynamical sign problem by utilizing quantum
computers. This idea is explored in a quantum algorithm we will review in
Section V.1.

$Q(y)=$+$\propto\xi$++$\propto\xi^{2}$+$\cdots$ Figure 4: Pertubrative
expansion of $Q(y)$ in powers of the perturbation strength $\xi$. All diagrams
displayed here correspond to diagonal elements of 0-boson diagrams.

##  III Quantum Algorithms for Ground State

###  III.1 Quantum Algorithms for Variational Monte Carlo

####  III.1.1 Quantum-enhanced Variational Monte Carlo

As mentioned in the introduction, VQE is a class of hybrid quantum-classical
algorithms proposed as a near-term alternative to the QPE for finding the
ground state of a Hamiltonian [14]. The main argument for the utility of VQE
relies on its compatibility with the resource constraint on the NISQ devices
[154]. The accuracy of the VQE is highly dependent on the choice of quantum
circuit ansätze, many of which were inspired by the classical CC theory [155,
156, 157].

Leveraging the VQE framework, Montanaro et al. proposed using the output from
VQE to accelerate the MCMC convergence of VMC simulation and termed this
quantum-enhanced VMC (QEVMC) [158]. Specifically, the output state from VQE is
used for initial walker samples in the MCMC algorithm that generates the
desired distribution $P(\bm{\theta},n)$ (Eq. 9). Since both VQE and VMC are
performed variationally, the VMC is guaranteed to find lower energy than the
original VQE input. The hope is that if VQE output can capture the ground
state to a good extent, or equivalently, the initial and final distributions
have more resemblances, then the number of steps in MCMC is expected to be
less. The other implicit assumption in this workflow is the efficient
classical calculation of the overlap $\langle n|\Psi(\bm{\theta})\rangle$,
which, however, could be hard for some $|\Psi(\bm{\theta})\rangle$ [159], for
instance, the unitary coupled-cluster quantum state [160].

Using the VQE distribution from Ref. 161, Montanaro et al. demonstrated that
MCMC converges to a target energy more rapidly than a single determinant
distribution. For Fermi-Hubbard and transverse field Ising models, Montanaro
et al. further numerically showed that the VQE distribution obtained for
modest lattice sites can be used to speed up MCMC in larger site calculations.

####  III.1.2 Non-unitary Post-Processing via Quantum Monte Carlo

Given the limited quantum resources in the NISQ era, classical post-processing
of quantum measurements could be used to reduce the circuit depth and improve
the accuracy of quantum-implemented variational wavefunction ansatz
$|\psi(\bm{\theta})\rangle=U(\bm{\theta})|\mathbf{0}\rangle$ [162]. For
example, inspired by classical QMC, Mazzola et al. proposed using the Jastrow
factor as a non-unitary post-processing projector. Specifically, they defined
the state
$|\psi(\bm{\theta},\bm{\phi})\rangle=\mathcal{P}{\bm{\phi}}|\psi(\bm{\theta})\rangle$,
where $\mathcal{P}{\bm{\phi}}$ serves to filter out higher energy components
from the suboptimal variational quantum heuristic ansatz. This approach was
applied successfully to several lattice models [163].

In general, the non-unitary classical processing operator $\hat{f}$ takes in
bit string $|n\rangle$ sampled from the quantum circuit and outputs f($n$).
With this, we then consider the wavefunction ansatz written as

$|\psi_{f}(\bm{\theta},\bm{\phi})\rangle=\hat{f}(\bm{\phi})U(\bm{\theta})|\mathbf{0}\rangle,$
(71)

and in the computational basis, the classical post-processing unit $\hat{f}$
reads

$\hat{f}(\bm{\phi})=\sum_{n}f_{\phi}(n)|n\rangle\langle n|.$ (72)

In this framework, the optimization is carried out simultaneously for the
quantum gate parameters $\bm{\theta}$ and the classical post-processing
parameters $\bm{\phi}$. Using efficient implementation with classical post-
processing, the energy expectation value

$E(\bm{\theta},\bm{\phi})=\frac{\langle\psi_{f}(\bm{\theta},\bm{\phi})|\hat{H}|\psi_{f}(\bm{\theta},\bm{\phi})\rangle}{\langle\psi_{f}(\bm{\theta},\bm{\phi})|\psi_{f}(\bm{\theta},\bm{\phi})\rangle}$
(73)

and its gradient (e.g., through parameter shift rule [164, 165]) requires no
more than the extra polynomial number of gates or measurements in addition to
the original circuit $U(\bm{\theta})$. Since the final wavefunction is not
normalized, one needs to account for the overlap explicitly. Following Eq. 71
and Eq. 72, the overlap can be written as

$\langle\psi_{f}(\bm{\theta},\bm{\phi})|\psi_{f}(\bm{\theta},\bm{\phi})\rangle=\sum_{n}|f_{\bm{\phi}}(n)|^{2}|\langle
n|U(\bm{\theta})|\mathbf{0}\rangle|^{2}.$ (74)

$|n\rangle$ can be sampled by measuring $U(\bm{\theta})|\mathbf{0}\rangle$ in
the computational basis, and the output bit string is sent to the classical
unit. This requires no modification to the existing implementation of
$U(\bm{\theta})$ on the quantum circuit.

The evaluation of the energy expectation value can be done in several ways.
One way is to transform $\hat{H}$ as
$\hat{f}(\bm{\phi})\hat{H}\hat{f}(\bm{\phi})$ and extra observables need to be
measured from the quantum computers [166]. However, this transformation
introduces an exponential number of Pauli strings that must be implemented on
the quantum hardware. The other method relies on an entangled copy to
reconstruct the measurement probability $|\langle
n|U(\bm{\theta})|\mathbf{0}\rangle|^{2}$ [163]. However, to achieve a target
precision $\epsilon$, the number of samples/measurements required scales
exponentially with the system size, or equivalently, the number of qubits
[167].

To overcome these exponential overheads, Zhang et al. [167] proposed an
efficient scheme for evaluating $\langle\psi_{f}|H|\psi_{f}\rangle$ by
appending a measurement circuit $M$ to the original $U(\bm{\theta})$. The
following discussion assumes that $\hat{H}$ consists of only one Pauli string.
In the case where $H$ has only $I$ and $Z$ operators $M$ is identity because
$H$ is diagonal in the computational basis. If the Hamiltonian has at least
one spin-flip operator (e.g $X$ or $Y$), then one can isolate one qubit that
is flipped and write the Hamiltonian as follows,

$\hat{H}=\sum_{\mathbf{s}}S(0\mathbf{s})|0\mathbf{s}\rangle\langle
1\tilde{\mathbf{s}}|+S(1\tilde{\mathbf{s}})|1\tilde{\mathbf{s}}\rangle\langle
0\mathbf{s}|,$ (75)

where $\mathbf{s}$ and $\tilde{\mathbf{s}}$ indicate all qubits excluding the
isolated qubit, and they are related by
$H|0\mathbf{s}\rangle=S(1\tilde{\mathbf{s}})|1\tilde{\mathbf{s}}\rangle$, with
$S(\mathbf{s})$ be the coefficients. The coefficients have an additional
constraint $S(0\mathbf{s})S(1\tilde{\mathbf{s}})=1$ because $\hat{H}$ is a
single Pauli string (e.g $\hat{H}^{2}=\mathrm{1}$). It is straightforward to
construct eigenvectors of the Hamiltonian Eq. 75, denoting as
$|+,\mathbf{s}\rangle$ and $|-,\mathbf{s}\rangle$ for each unique bit string
$\mathbf{s}$ and are spanned by
$\\{|0\mathbf{s}\rangle,|1\tilde{\mathbf{s}}\rangle\\}$. The numerator in Eq.
73 is then written as

$\displaystyle\langle\psi_{f}|H|\psi_{f}\rangle=\sum_{\mathbf{s}}|\langle+,\mathbf{s}|U(\bm{\theta})|\mathbf{0}\rangle|^{2}f(0\mathbf{s})f(1\tilde{\mathbf{s}})$
(76)
$\displaystyle-|\langle-,\mathbf{s}|U(\bm{\theta})|\mathbf{0}\rangle|^{2}f(0\mathbf{s})f(1\tilde{\mathbf{s}}).$

The measurement circuit $M$ to obtain
$|\langle+,\mathbf{s}|U(\bm{\theta})|\mathbf{0}\rangle|^{2}$ can then be
constructed by appending a controlled rotation gate with the source being the
isolated qubit and the target being all other qubits on which spin-flip
operators from the Hamiltonian ($X/Y$) acted. The isolated qubit is then
measured in the $X$/$Y$ basis, and all other qubits are measured in the $Z$
basis. The controlled gates ensure the bit string on all other qubits is
flipped properly, conditioned on the outcome of the isolated qubits. The
number of extra two-qubit gates depends on the number of $X$/$Y$ operators in
the Hamiltonian.

The implementation presented in the paper is efficient if the number of these
operators does not scale exponentially with the system size, or equivalently,
the Hamiltonian has local interactions only. Furthermore, Zhang et al. [167]
showed that the number of measurements when introducing the measurement
circuit $M$ scales as $\mathcal{O}(\text{Poly}(r)/\text{Poly}(\epsilon))$,
where $\epsilon$ indicates the precision (standard deviation from the
measurement outcome)and $r$ is the range of the function $f_{\bm{\phi}}(n)$.
Since the classical processing unit can be designed within a specific range,
the number of measurements has polynomial scaling similar to that of the
initial VQE with the same target precision.

Moreover, the use of quantum circuits like those presented here can be viewed
as an improvement to VMC. Specifically, given that the energy expectation
value in Eq. 73 can be evaluated efficiently up to an additive error with
polynomial quantum resources, one can perform VMC without resorting to Monte-
Carlo algorithms to sample the local energy following distribution in Eq. 9.
With this insight, recent algorithmic developments of hybrid algorithms have
introduced methods to efficiently sample the expectation values
$\langle\Psi(\bm{\theta})|\hat{H}|\Psi(\bm{\theta})\rangle$ without relying on
MCMC. Interested reads are referred to Ref. [168, 169, 170, 171] for more
technical details.

####  III.1.3 Quantum Data-enhanced Variational Monte Carlo

The accuracy of VMC heavily relies on the parametrization of the variational
ansätze. Employing more expressive ansätze could enhance the performance of
VMC calculations, an avenue for which where neural networks (NN) have been
valuable. Recent years have seen the development of different NN ansätze in
VMC calculation applying to various problems, including the ground state
energies of various spin models [172], ground state properties of real solids
[173], ground and excited state energies of Heisenberg models for frustrated
magnets [174, 175, 176] and steady states of open quantum systems [177, 178],
etc.

Recently, with the emergence of programmable quantum devices composed of
Rydberg atom arrays [179, 180, 181, 182], researchers have begun to explore
the ground state and phase transitions in these platforms. A Rydberg array
consisting of $N$ qubits can be described using the following Rydberg
Hamiltonian parametrized by Rabi frequency $\Omega$, detuning $\Delta$, and
the long-range van der Waals interaction between Rydberg states $V_{ij}$:

$\hat{H}=\frac{\Omega}{2}\sum_{i=1}^{N}\hat{\sigma}_{i}^{x}-\Delta\sum_{i=1}^{N}\hat{n}_{i}+\sum_{i<j}V_{ij}\hat{n}_{i}\hat{n}_{j}.$
(77)

Here $\hat{\sigma}_{i}^{x}=\ket{g}_{i}\bra{r}_{i}+\ket{r}_{i}\bra{g}_{i}$,
where $\ket{g}_{i}$ denotes the ground state of the atom at site $i$, and
$\ket{r}_{i}$ is the Rydberg state. The interaction between Rydberg states
$V_{ij}=V_{0}/|\mathbf{x}_{i}-\mathbf{x}_{j}|^{6}$ comes from the so-called
Rydberg blockade mechanism [183], which can be tuned by the effective blockage
range $R_{\mathrm{b}}/a$, where $R_{\mathrm{b}}\equiv(V_{0}/\Omega)^{1/6}$ is
the blockade radius, and $a$ is the inter-atomic spacing in the array.

NN ansätze for VMC turns out to be a suitable tool for studying this
Hamiltonian, as shown in Ref. [184]. Prior to this work, generative models
such as the restricted Boltzmann machine (RBM) have been proven effective for
state reconstruction from the measurements using a quantum circuit [185]. They
compared few-body observables and R’enyi entropy derived from the
reconstructed wavefunction and found strong agreement with both experimental
outcomes and simulation results obtained using the Lindblad master equation.
Combining those efforts, Moss et. al. proposed a two-phase optimization of a
neural network ansatz to target the ground state of the Rydberg Hamiltonian
[186], where they first trained the neural network using the raw, unprocessed
experimental data obtained from the measurements of a Rydberg atom array
undergoing the phase transition between disordered and checkerboard phases
[182], and performed subsequent variational optimization with the resulting NN
wavefunction. They employed both one-dimensional (1D) and two-dimensional (2D)
recurrent neural network (RNN) to encode the ground state wavefunction, a type
of universal function approximators [187] that can encode arbitrarily complex
functions by increasing the number of parameters.

A notable improvement was shown using a hybrid data-driven and Hamiltonian-
driven approach. Naive variational optimization of 1D RNN wavefunction
struggles to reach the ground state, nevertheless, the pretraining stage
significantly improves the quality of the resulting variational wavefunction,
as shown by the fact that the energy calculated using the hybrid strategy
aligns more closely with the true ground state energy. Although variationally
optimized 2D RNN wavefunction itself can capture most of the characteristics
of the ground state, the pretraining step still accelerates the convergence
towards the ground state. Subsequent calculations on other observables further
demonstrate the significance of the pretraining step. These results highlight
the possibility of using analog quantum simulations in QMC. Namely, one can
use the data from quantum computations to accelerate convergence and enhance
accuracy in classical simulations.

###  III.2 Hybrid Projector Quantum Monte Carlo

####  III.2.1 Hybrid Quantum–Classical AFQMC

The quantum-classical hybrid QMC (QC-QMC) method has opened new avenues for
projector QMC methods, such as AFQMC, in simulating chemical systems [159].
Subsequent studies in refining the framework and assessing the performance
[171, 188, 189, 190] have been sharpening the path towards applying quantum-
classical hybrid AFQMC (QC-AFQMC) towards more precise and experimentally
relevant quantum chemical calculations. In QC-AFQMC [159], the trial
wavefunction for the phaseless approximation is prepared on a quantum device.
At the same time, imaginary time propagation and energy estimation are
performed on a classical computer. In their original work, they prepared a CC-
type trial wave function on a Sycamore processor [191]. The QC-AFQMC algorithm
was shown to be robust to some of the common limitations imposed by circuit
noise and hardware imperfections in the NISQ era. The imaginary-time evolution
and sign problem in PMC methods and the workflow of QC-AFQMC are illustrated
in Fig. 5.

Figure 5: a. Schematic illustration of the imaginary-time evolution, which
converges to the ground state in the long-time limit. b. Demonstration of the
sign problem. The noise-to-signal ratio will grow exponentially if no
constraint is imposed (top). Constrained QMC controls the sign problem at the
expense of introducing a bias, and using quantum constraint reduces the bias.
c. Summary of the QC-AFQMC algorithm. The walker wavefunctions, denoted as
$\\{\ket{\phi_{i}}\\}_{\tau}$, are evolving in time according to the
propagator $\mathrm{e}^{-\Delta\tau\hat{H}}$, which is achieved by updating
the gate parameters each time step. The resulting observables are calculated
on a classical computer with overlap queries to the quantum processor. Figure
adapted from Ref. 159.

Using this algorithm, they first calculated the atomization energy of H4, a
typical testing ground for electronic structure methods for strongly
correlated systems. There are two equally important mean-field states for the
ground state of this system, rendering AFQMC with a single Slater determinant
trial inaccurate. However, AFQMC with a generalized valance bond perfect
pairing [192] (GVB-PP) trial wavefunction prepared by a quantum circuit is
able to produce results with chemical accuracy (1 kcal/$\mathrm{mol}$). It is
also noteworthy that although the trial prepared with the quantum circuit is
far from ideal, i.e., the trial energy significantly deviates from the
expected outcomes from a noiseless quantum computer. Despite this, QC-AFQMC
still achieves chemical accuracy in the atomization energy. Calculating the
potential energy surface of N2 and the cold curve of diamond using a minimal
unit cell to nearly chemical accuracy also showcased the improvement of the
AFQMC energy using a more quantum circuit trial state as well as the noise-
resilience of the QC-AFQMC algorithm.

As mentioned in Section. II.1.2, the key quantity in (QC-)AFQMC is the overlap
between the trial wavefunction and the walkers. Shadow tomography [193, 169]
was used to efficiently measure the overlap between the quantum trial
wavefunction and the walker wavefunction, enabling the prediction of
properties of a quantum system with relatively few measurements. More
specifically, classical shadows were utilized to obtain the representation of
the trial state on a computational basis, thus bypassing repeated overlap
measurements for each time step and walker. A detailed discussion of classical
shadows techniques can be found in Sec. III.2.2.

Two computational considerations are needed to scale up QC-AFQMC. The first
consideration is the cost of classical post-processing. The original QC-AFQMC
work used Clifford-based classical shadows, causing an exponential scaling
postprocessing cost. This exponential overhead has been removed using a
different channel in classical shadows [171]. Even with the new shadows
technique, the classical post-processing per sample was found to be too high
($\mathcal{O}(N^{8})$), but this cost was brought down to $\mathcal{O}(N^{4})$
in a recent work [75]. We will discuss more details in Section III.2.2.

The second consideration is the measurement sample complexity associated with
the overlap estimation. The overlap of the walker wavefunction with a quantum
trial wavefunction is expected to decay exponentially with the system size
[159]. This means for a fixed relative error one needs to refine the overlap
estimation exponentially more when increasing the system size. This results
into an exponential-scaling measurement complexity in QC-QMC. On the classical
computer, one can completely address this issue by directly estimating the
logarithm of the overlap, which is currently not possible to achieve on the
quantum computer. To postpone this orthogonality catastrophe, it was pointed
out that only consider a small physical space by leveraging the virtual
correlation technique [159]. This way, the overlap will decay, in the worst
case, exponentially with respect to the correlation length within the active
space rather than the system size. Only the overlap within the active space is
then computed with a quantum computer.

To further investigate the exponential challenge in overlap evaluation,
Mazzola et al. [194] conducted numerical tests on a one-dimensional transverse
Ising model (TFIM), described by the Hamiltonian

$\hat{H}=-J\sum_{k}^{L}\sigma_{k}^{z}\sigma_{k+1}^{z}-\Gamma\sum_{k}^{L}\sigma_{k}^{x}.$
(78)

Even though this model can be mapped to a free fermion problem where AFQMC is
exact without sampling, the authors studied this with GFMC. They discovered
that even with an exponentially increasing number of states using a Jastrow-
type trial wavefunction, the statistical error in the local energy reaches
approximately 100% when $\Gamma/J=1$. Additionally, they observed that the
error in the energy per site increases exponentially relative to the number of
sites, even when using the exact ground state wavefunction as the trial
wavefunction. However, as Lee et al. showed in Ref. 195, the error in the
local energy depends on not only the form of the trial wavefunction but also
the interaction strength $\Gamma/J$. They found that when $\Gamma/J=0.5$, the
error in the local energy per site has weak system size dependence, whereas
they obtained similar results for $\Gamma/J=1$. Furthermore, they noted a
significant reduction in error when employing more sophisticated trial
wavefunctions and walker wavefunctions, specifically
$\ket{\Psi_{\mathrm{T}}}=\mathrm{e}^{-0.05\hat{H}}\ket{\Psi_{\mathrm{MC}}}$,
where $\ket{\Psi_{\mathrm{MC}}}$ is the trial wavefunction used in Ref. 194.
The walkers are generated by spin flips from $\ket{\Psi_{\mathrm{MC}}}$. The
comparison is shown in Fig. 6. This demonstrates the effect of using a more
accurate trial wavefunction. Notably, using a more accurate trial wavefunction
does not eliminate the exponential overhead. Instead, the exponentially
vanishing overlap persists, albeit with a significantly slower decay rate.
Consequently, the Monte Carlo simulation may not be bottlenecked by the
overlap measurement.

Figure 6: Distribution of overlap and local energy per site across states in
12 sites 1-dimensional TFIM for different trial and walker wavefunction as
well as the number of measurements. For (a) and (b),
$\ket{\Psi_{\mathrm{T}}}=\ket{\Psi_{\mathrm{MC}}}$ and single computational
basis state walkers were used, while for (c) and (d), the trial wavefunction
is given by
$\ket{\Psi_{\mathrm{T}}}=\mathrm{e}^{-0.05\hat{H}}\ket{\Psi_{\mathrm{MC}}}$
and the walkers are generated by spin flips from $\ket{\Psi_{\mathrm{MC}}}$.
16 independent simulations were run with $M$ samples each time for the
measurement of overlap and local energy. Figure replotted with the data from
Ref. 195.

M. Amsler et al. further tested the efficacy of the QC-AFQMC algorithm in
various systems, including molecules such as H4, ozone, and oxygen, as well as
model systems such as the quasi-1D Fermi–Hubbard model derived from CuBr2
[188]. For H4, they compared the AFQMC energy using VQE and RHF as the trial
wavefunction, where they used a unitary coupled-cluster single and double
(UCCSD) ansatz for the VQE calculation and found that both the VQE energy and
AFQMC energy using VQE trial are close to FCI result. In addition, AFQMC
energies using a VQE trial stay within chemical accuracy compared to FCI
results except at intermediate bond lengths, and it is considerably closer to
FCI results than the plain VQE approach. They also found that the relative
energies of ozone and singlet oxygen compared to ground-state oxygen can be
accurately predicted using AFQMC with a VQE trial wavefunction. However, this
calculation yields a slightly larger error than AFQMC with a complete active
space self-consistent field (CASSCF) trial wavefunction. They also demonstrate
the capability of AFQMC to accurately compute the ground state energy of
extended systems, as evidenced by calculations on the Fermi–Hubbard model. In
conclusion, quantum trial wavefunctions, such as those obtained from VQE,
generally enhance the quality of AFQMC calculations. However, one needs to be
cautious about the exponential scaling measurement overhead.

####  III.2.2 Classical Shadows for QCQMC

We use $\rho$ to denote the density matrix of the $N$-qubit prepared quantum
state. Classical shadows provide a convenient and powerful framework for
efficiently estimating the expectation value of general observables $O$ from
measurements [159]. In QC-AFQMC, we want to compute expectation values of
operators like $O=|\phi\rangle\langle 0|$ where $|\phi\rangle$ is a Slater
determinant state, and $|0\rangle$ is the vacuum state. These observables
contain many non-trivial matrix elements and are not easy to directly measure
on any single basis.

The key idea of classical shadows is to combine information-dense snapshots
(i.e., projective measurements of each qubit) with measurements in different
bases to estimate arbitrary observables with rigorous convergence guarantees.
The shadow procedure constructs an efficient classical representation
$\hat{\rho}$ of the quantum state $\rho$ from measurements, such that
expectation values of generic observables can be approximately evaluated by

$\mathrm{tr}[O\hat{\rho}]\approx\mathrm{tr}[O\rho]$ (79)

up to an additive error $\epsilon$. Before measurement, a unitary $U$ is
applied to the state to select the measurement basis. Then the state $U\rho
U^{\dagger}$ is measured in the computational basis
($\left.\\{|b\rangle\\}_{b\in\\{0,1\\}^{N}}\right)$). Thus the state is mapped
from $\rho\rightarrow U^{\dagger}|b\rangle\langle b|U$ where $U^{\dagger}$
inverts the unitary applied pre-measurement. In the typical shadows procedure,
$U$ is sampled from a distribution $\mathcal{U}$, so we denote averaging over
unitaries by $\mathbb{E}_{U\sim\mathcal{U}}$. Then, for a sampled $U$, the
probability of measuring any particular $b$ is given by

$p_{U}(b)=\langle b|U\rho U^{\dagger}|b\rangle,$ (80)

The map from the input state to a measured state is known as the measurement
channel and can be written explicitly as

$\mathcal{M}(\rho)=\mathbb{E}_{U\sim\mathcal{U}}\sum_{b}\langle b|U\rho
U^{\dagger}|b\rangle U^{\dagger}|b\rangle\langle b|U$

By inverting $\mathcal{M}(\rho)$, an estimate for $\rho$ can be constructed
from a measured state $U_{m}^{\dagger}|b_{m}\rangle\langle b_{m}|U_{m}$, where
$U_{m}$ and $b_{m}$ are the unitary and bitstring sampled during the $m$-th
run of the experiment.

The average over measured states is known as a classical shadow and is written
as

$\displaystyle\hat{\rho}$ $\displaystyle=\frac{1}{M}\sum_{m}\hat{\rho}_{m}$
(81) $\displaystyle\hat{\rho}_{m}$
$\displaystyle=\mathcal{M}^{-1}(U_{m}^{\dagger}|b_{m}\rangle\langle
b_{m}|U_{m})$ (82)

In general, $\mathcal{M}(\rho)$ is only invertible if the ensemble of
rotations $\mathcal{U}$ is tomographically complete [159]. In that case, it is
easy to show that $\mathbb{E}[\hat{\rho}_{m}]=\rho$, the true state is
recovered on average. More generally, if $\mathcal{M}$ is not invertible, a
pseudo-inverse can be used [196], in which case $\mathbb{E}[\hat{\rho}_{m}]$
recovers a projection of $\rho$ onto the measurable (i.e. "visible")
operators. Then, the expectation value of (visible) observables $O$ can be
estimated accurately by evaluating $\mathrm{tr}[O\hat{\rho}]$.

Different shadow procedures are characterized by different choices of
ensembles $\mathcal{U}$. There are two key characteristics that determine the
performance of a shadow procedure. The first is the sample complexity, which
depends on the variance of the estimator and quantifies how many measurements
are needed to estimate observables $O$ to within an additive error $\epsilon$.
The second is the classical post-processing complexity. For certain special
ensembles, like random single-qubit rotations or random Clifford circuits, the
states $U^{\dagger}|b\rangle\langle b|U$ and inverse channel
$\mathcal{M}^{-1}$ take a simple classically tractable form. The complexity
also depends on the type of observable being measured. For example, Pauli
observables can be efficiently classically computed with respect to Clifford
(i.e. stabilizer) states.

In the original QC-AFQMC proposal, the shadow tomography is performed using
random unitaries sampled from the Clifford group of $N$ qubits (henceforth
referred to as Clifford shadows) [159], but requires exponential classical
post-processing cost. Wan et al. removed the exponential overhead by replacing
the Clifford group with the ensemble of random matchgate circuits [171].
Matchgate circuits are qubit representations of fermionic Gaussian unitaries
$U_{Q}$, which are unitaries that satisfy

$U_{Q}^{\dagger}\gamma_{\mu}U_{Q}=\sum_{\nu=1}^{2n}Q_{\mu\nu}\gamma_{\nu},$
(83)

where $\gamma_{\mu}$ is the $\mu$-th Majorana operator obtained by
$\gamma_{2j-1}=a_{j}+a_{j}^{\dagger}$,
$\gamma_{2j}=-i(a_{j}-a_{j}^{\dagger})$. The key part of their results is the
$\mathcal{O}(N^{4})$ complexity for obtaining the overlap from a single
snapshot. Additionally, if there are $M$ observables, the number of classical
shadow samples required to estimate every $\text{Tr}(O_{i}\rho),i=1,\ldots,M$
to within additive error $\varepsilon$ with probability at least $1-\delta$ is
[171]

$N_{\text{samp}}=\mathcal{O}\left(\frac{\log(M/\delta)}{\varepsilon^{2}}\max_{1\leq
i\leq M}\mathrm{Var}[\hat{o}_{i}]\right),$ (84)

where $\hat{o}_{i}$ is defined as

$\hat{o}_{i}=\text{Tr}\left(O_{i}\mathcal{M}^{-1}(\hat{U}^{\dagger}|\hat{b}\rangle\langle\hat{b}|\hat{U})\right).$
(85)

Here $\hat{U}$ is a random variable sampled from the distribution
$\mathcal{U}$ and $\mathbb{P}[\ket{\hat{b}}=\ket{b}|\hat{U}=U]=\langle b|U\rho
U^{\dagger}|b\rangle$. It is also proved that the variance bound of
$\hat{o}_{i}$ is $\mathcal{O}(\sqrt{N}\log N)$.

Combining these facts, the overall complexity for classical processing of the
overlap obtained from matchgate shadows scales as
$\mathcal{O}(\frac{\log(M/\delta)}{\varepsilon^{2}}MN^{4.5}\log N)$, and the
post-processing for the energy estimates scales as
$\mathcal{O}(\frac{\log(M/\delta)}{\varepsilon^{2}}MN^{8.5}\log N)$ [190], as
opposed to the exponential scaling using Clifford shadows. Nevertheless, as
Huang et al. showed in Ref. 190, even though the classical post-processing for
matchgate shadows is asymptotically efficient, it is still practically
computationally intensive. They carried out an experiment and found noise
resilience in the overlap ratio obtained by matchgate shadows. However, it
would require up to $10^{11}$ CPU hours to perform the classical post-
processing of a (54e, 50o) calculation according to the prediction in Ref.
190. Kiser et al. also investigated the comparison of different protocols of
overlap estimation [189]. Taking account of the $\mathcal{O}(N^{2})$ number of
samples suffice to maintain the fixed additive accuracy for the energy
estimator [68], they found that the total cost for evaluating all overlaps
scales as $\mathcal{O}(N^{9})$ for matchgate shadows, and
$\mathcal{O}(N^{10})$ for orbital rotation shadows, which is almost
intractable in real applications.

This computational cost in classical post-processing presents a significant
challenge to the practical application of QC-AFQMC using classical shadows for
overlap estimation, and addressing this issue is a crucial task for future
research. Nonetheless, there is promise in reducing the scaling to
$\mathcal{O}(N^{4\text{-}5})$ with a differential programming approach [75].

The simulation of classical shadows circuit is exponentially hard for
classical computers as the random unitatires maximize the entanglement growth.
One could hope to find a suitable application of classical shadows in
computational chemistry to demonstrate quantum advantage.

####  III.2.3 Hybrid Quantum–Classical FCIQMC

Figure 7: (a) Sketch of the spawning, death/cloning, and annihilation steps in
FCIQMC with the determinants transformed with ADAPT-VQE quantum wavefunction.
(b) The potential energy surface for the N2 with different methods under the
STO-3G basis set. Figure modified from Ref. 197.

Unlike the polynomially scaling ph-AFQMC method, FCIQMC preserves the
exponential scaling of FCI with exact on-average accuracy [198]. The
exponential scaling sample complexity of FCIQMC is due to the sign problem. A
simple choice to reduce the sign problem in FCIQMC will be performing a
suitable unitary transformation of the Hamiltonian by
$\tilde{H}=U^{\dagger}HU$ to reduce the “non-stoquaticity" [199]. One can
observe that an exact diagonalization leads to a fully diagonal Hamiltonian
matrix free of the sign problem. This diagonalization, however, corresponds to
a brute-force FCI calculation and renders any subsequent computations
obsolete, as the problem is already solved.

Zhang et al. [197] prepared the quantum circuit $U$ from an ADAPT-VQE
simulation to transform the basis of original Slater determinants,
$U|\psi_{i}\rangle\rightarrow|\tilde{\psi}_{i}\rangle$, as shown in Fig. 7(a).
This combination was found to significantly reduce the sign problem in FCIQMC,
yielding results that surpass those obtained using ADAPT-VQE alone, as shown
in Fig. 7(b).

Kanno et al. [200] also used the solution from VQE with a hybrid tensor
network (HTN) ansatz to prepare the trial wavefunction for QC-FCIQMC
calculations. Subsequently, they employed quantum wavefunction exclusively for
mixed energy estimations during propagation steps, a technique referred to as
quantum-assisted energy evaluation (QAEE) in Ref.201. Such an implementation
is possible because the FCIQMC approach works without applying constraints in
the presence of a sign problem.

A critical step in FCIQMC is the calculation of Hamiltonian matrix elements
$H_{ij}=\langle\psi_{i}|\hat{H}|\psi_{j}\rangle$ to determine the walker
spawning probability. Classical FCIQMC can easily compute the non-zero
elements for a given state $|\psi_{i}\rangle$ using the Slater-Condon rules.
However, calculating $\tilde{H}_{ij}$ in a transformed basis as needed in QC-
FCIQMC is not straightforward without prior knowledge of which indices $j$
give non-zero $\tilde{H}_{ji}$. As for $\hat{H}=\sum_{k}h_{k}\hat{P}_{k}$
where $h_{k}$ is the coefficient for Pauli string $\hat{P}_{k}$, a sampling
approach was introduced to measure $\tilde{H}_{ji}$ to determine the spawning
probability in QC-FCIQMC on the quantum computer [197]:

$|\tilde{H}_{ji}|^{2}=\sum_{kk^{\prime}}h_{k}h_{k^{\prime}}p_{kk^{\prime}}^{i}(j)$
(86)

with

$p_{kk^{\prime}}^{i}(j)=\operatorname{Re}\left\langle
i\left|U^{\dagger}P_{k}U\Pi_{j}U^{\dagger}P_{k^{\prime}}U\right|i\right\rangle$
(87)

which, for fixed $i,k,k^{\prime}$, can be viewed as the distribution over
different $j$s and can be sampled with quantum circuit measurements and post-
processing [197]. The sign of $\tilde{H}_{ji}$ can be further estimated with
another circuit. Although all the elements of $\tilde{H}_{ji}$ can be
precomputed and stored prior to FCIQMC calculations, this method faces
scalability challenges as system size increases. It is important to find an
efficient sampling approach for the spawning probability to scale up the
capability of QC-FCIQMC beyond classical FCIQMC.

###  III.3 Hybrid Quantum–Classical Unitary CCMC

Recently, Filip proposed a present a Monte Carlo informed approach to
projective quantum eigensolver (PQE) algorithm [202], termed MC-PQE [203].
This method was designed to accurately estimate ground state energies and used
the intrinsic noise of measurements to justify an FCIQMC-type approach. In
this algorithm, the quantum device is used to sample Hilbert space throughout
imaginary time propagation. To encode the UCC ansatz on a quantum computer,
the Jordan–Wigner decomposition of fermionic creation and annihilation
operators is used, followed by the application of a Pauli gadget [204] to
exponentiate the resulting Pauli strings. The algorithm evaluates unlinked
residuals

$r_{\mathbf{i}}^{u}=\langle\Phi_{\mathbf{i}}|\hat{H}-S(\tau)|\Psi(\tau)\rangle$
(88)

and overlaps

$s_{\mathbf{i}}=\langle\Phi_{\mathbf{i}}|\Psi(\tau)\rangle$ (89)

on a quantum device. As introduced in Sec. II.1.2, $\mathbf{i}$ runs over all
possible excitations, and $|\Phi_{\mathbf{i}}\rangle$ are corresponding
determinants. The shift $S(\tau)$ is the same as in Eq. (II.1.2) and is
dynamically adjusted to stabilize the total walker number around some target
number. Instead of solving the equations defined in Eq. (88) on a classical
computer as proposed by Stair et al. [202], $r_{\mathbf{i}}^{u}$ is assumed to
be a random variable. Each time step, all residuals are evaluated, and
populations on determinants $|\Phi_{\mathbf{i}}\rangle$ are updated according
to

$N_{\mathbf{i}}(\tau+\Delta\tau)=N_{\mathbf{i}}(\tau)-\Delta\tau
N_{0}r_{\mathbf{i}}^{u},$ (90)

where $N_{0}$ is the initial number of walkers on the reference determinant.
To evaluate the projected energy, one can compute
$E_{\mathrm{proj}}(\tau)=r_{0}^{u}+S(\tau)s_{0}$. Repeating the evolution
prescription in Eq. (90) eventually leads to an estimate of the ground state
as predicted by the UCC ansatz. While this quantum version of the UCCMC
algorithm circumvents the expensive sampling of excitation and de-excitation,
the encoding of UCC ansätze on quantum devices usually involves deep quantum
circuits. Furthermore, this does not necessarily alleviate the sign problem
inherent to CCMC. As it is in its early exploration, it is currently unclear
what advantage this would offer beyond UCC and classical CCMC.

###  III.4 Entanglement Forging with Monte Carlo Sampling

The VQE algorithm, when combined with entanglement forging (EF) [205],
represents a promising strategy that effectively simulates a $2N$-qubit system
using quantum circuits with only $N$-qubit after classical postprocessing. The
Schmidt decomposition of the wavefunction is the starting point of EF for a
$N+N$ qubits bipartite quantum system $\hat{H}_{A}\otimes\hat{H}_{B}$.

$|\psi\rangle=U\otimes
V\sum_{n=1}^{2^{N}}\lambda_{n}|\sigma_{n}\rangle_{A}|\sigma_{n}\rangle_{B},$
(91)

where $U$ and $V$ are unitary quantum circuits. Fig. 8(a) graphically shows
the decomposition of $2N$-qubit system. The VQE optimization contains the
parameters in the $U$ and $V$ circuits, as well as the Schmidt coefficients.
In this approach, evaluating expectation values requires the products of
Schmidt coefficients, $\lambda_{m}\lambda_{n}$ (for $1\leq m,n\leq 2^{N}$). We
note that the number of these combinations grows exponentially with the system
size.

With EF, the energy evaluation involves the expectation value of operators
$O_{A}\otimes O_{B}$ with Pauli strings $O_{A},O_{B}\in\\{I,X,Y,Z\\}^{\otimes
N}$,

$\displaystyle\langle O_{A}\otimes
O_{B}\rangle=\sum_{n=1}^{2^{N}}(\lambda_{n}^{2}\langle\sigma_{n}|\tilde{O}_{A}|\sigma_{n}\rangle\langle\sigma_{n}|\tilde{O}_{B}|\sigma_{n}\rangle+\sum_{m=1}^{n-1}\lambda_{n}\lambda_{m}\times$
(92)
$\displaystyle\sum_{p\in\mathbb{Z}_{4}}(-1)^{p}\langle\phi_{\sigma_{n}\sigma_{m}}^{p}|U^{\dagger}O_{A}U|\phi_{\sigma_{m}\sigma_{m}}^{p}\rangle\langle\phi_{\sigma_{n}\sigma_{m}}^{p}|V^{\dagger}O_{B}V|\phi_{\sigma_{n}\sigma_{m}}^{p}\rangle)$

where
$\left|\phi_{xy}^{p}\right\rangle=\left(|x\rangle+\mathrm{i}^{p}|y\rangle\right)/\sqrt{2}$
with $p\in\mathbb{Z}_{4}=\\{0,1,2,3\\}$. Fig. 8(b) shows the
$|\phi_{01}^{p}\rangle$ for 2 qubit case with minimal basis. Eq. (92) can be
decomposed into measuring expectation values with independent quantum circuits
of $N$ qubits.

$\left\langle\psi\left|O_{A}\otimes
O_{B}\right|\psi\right\rangle=\sum_{a}\mu_{a}\operatorname{Tr}\left(O_{A}\rho_{a}\right)\operatorname{Tr}\left(O_{B}\rho_{a}\right)$
(93)

where $\rho_{a}$ is $N-$qubit pure state quantum circuit, and $\mu_{a}$ is
proportional to the products of two independent Schmidt coefficients. We
further define

$O_{A/B}=V_{A/B}^{\dagger}\left(\sum_{x\in\\{0,1\\}^{N}}O_{A/B}(x)|x\rangle\langle
x|\right)V_{A/B},$ (94)

where $V_{A/B}$ includes only Clifford gates that map the Pauli strings $X$
and $Y$ appearing in $O_{A/B}$ to Pauli string $Z$ and $O_{A/B}(x)$ will take
the values 1 or -1 depending on the parity of the bitstring $x$. This converts
the problem of measuring the $2N-$qubit observables (Eq. (92)) to measuring
the observables on a series of $N$ qubits system. Namely, each experiment
prepares a state $V_{A/B}\rho_{a}V_{A/B}^{\dagger}$ and measures it in a
computational basis.

The computation of Eq. (93) appears to be not scalable, which involves
exponentially many terms An importance sampling approach was proposed where an
unbiased estimator for Eq. (93) is obtained by $M$ samples [205],

$f=\frac{\|\mu\|_{1}}{M}\sum_{j=1}^{M}\operatorname{sgn}\left(\mu_{a_{j}}\right)O_{A}\left(x^{j}\right)O_{B}\left(y^{j}\right)$
(95)

where $\|\mu\|_{1}\equiv\sum_{a=1}^{\ell}\left|\mu_{a}\right|$. $O_{A}(x^{i})$
and $O_{B}(x^{j})$ are measured by preparing the state
$V_{A}\rho_{a_{i}}V_{A}^{\dagger}$ and $V_{B}\rho_{a_{j}}V_{B}^{\dagger}$ and
measuring in the computational basis. The sampling probability is given by

$\pi_{a}=\left|\mu_{a}\right|/\|\mu\|_{1},$ (96)

which is proportional to the products of Schmidt coefficients
$\lambda_{m}\lambda_{n}$. The measurements of $2N$-qubit quantum circuit are
thereby mapped to a series of $N-$qubit quantum circuit measurements by
importance sampling with probability proportional to $\lambda_{n}\lambda_{m}$,
as shown in Fig. 8(c). It was proven that $S$ measurement experiments ($S=2M$
for measuring both $V_{A}\rho_{a_{i}}V_{A}^{\dagger}$ and
$V_{B}\rho_{a_{i}}V_{B}^{\dagger}$) are needed for a target precision
$\epsilon$:

$S\sim\left(\frac{1}{\epsilon}\sum_{n,m}\left|\lambda_{n}\lambda_{m}\right|\right)^{2}=\frac{\|\vec{\lambda}\|_{1}^{4}}{\epsilon^{2}}$
(97)

Since the one-norm decreases toward 1 in the limit of weak entanglement, the
overhead of entanglement forging is smaller for simulations of states
divisible into weakly entangled fragments (i.e., halves). As an example, Fig.
8(d) shows the decay of leading Schmidt coefficients for molecular ground
states with bipartition between spin-up and spin-down particles. Apart from
this Monte Carlo approach, other approaches using generative neural networks
[206, 207] were also proposed to address the sampling process. The precise
assumptions and conditions in which this method could demonstrate a practical
quantum advantage in computational chemistry seem unclear as it is still in
its early development stage.

Figure 8: Entanglement forging of a bipartite system. (a) Quantum state
$|\Psi\rangle$ obtained by a combination of bitstring states $\ket{b_{n}}$.
(b) Example of a two-qubit state rewritten using superpositions of single-
qubit states. (c) Reconstruction of the entangled state $\ket{\Psi}$ from
$N$-qubit circuits initialized in bitstrings and pairwise superpositions. (d)
The rapid or slow decay of Schmidt coefficients indicates weak or strong
entanglement in the ground state of different molecules. Figure adapted from
Ref. 205.

###  III.5 Quantum Selected Configuration Interaction

Recently, Kanno et al. [208] proposed a quantum version of Monte-Carlo
configuration interaction [89, 90] by employing a quantum device to select
determinants that are used to construct the effective Hamiltonian. In this
quantum selected configuration interaction (QSCI) algorithm, an input state
$\ket{\psi_{in}}$ and a set of $n$-qubit string basis states $\ket{i}$ are
prepared on a quantum computer by VQE or other algorithms. One then
iteratively projects $\ket{\psi_{in}}$ onto the computational basis for a
total of $N_{\text{shot}}$ times to obtain

$\displaystyle\ket{\psi}=\sum_{i=0}^{2^{n}-1}\alpha_{i}\ket{i}.$ (98)

Then, one chooses the $R$ most frequent configurations to define the CI
subspace $\mathcal{S}$. Finally, the diagonalization of the Hamiltonian in the
subspace $\mathcal{S}$ is carried out on a classical machine, and the ground
state energy is obtained. One could straightforwardly extend this approach to
calculating excited states by repeating this procedure on multiple input
states.

An important technical detail in QSCI is how to prepare the input state on a
quantum computer. This could be done using a quantum circuit optimized using
VQE, such that QSCI can essentially be viewed as a post-processing method to
refine VQE results. Related to this approach, the adaptive derivative-
assembled pseudo-Trotter ansatz (ADAPT) [209] was also used in a method termed
ADAPT-QSCI [210]. This algorithm provides a way of systematically preparing
viable input states for QSCI with reduced gate count.

Recently, Robledo-Moreno et al. reported a similar hybrid quantum-classical
selected CI approach with additional error mitigation techniques [211]. In
this algorithm, a truncated version of the local unitary cluster Jastrow
(LUCJ) ansatz $\ket{\Psi}$ [212, 159],

$\displaystyle\ket{\Psi}=\prod_{\mu=1}^{L}\mathrm{e}^{\hat{K}_{\mu}}\mathrm{e}^{\mathrm{i}\hat{J}_{\mu}}\mathrm{e}^{-\hat{K}_{\mu}}\ket{\Psi_{\mathrm{RHF}}},$
(99)

is used to generate a set of configurations $\\{\chi\\}$. Here,
$\hat{K}_{\mu}$ is a generic one-body operator, and $\hat{J}_{\mu}$ is a
density-density operator. This set of configurations is then subject to a
classical post-processing routine where an iterative configuration recovery
loop is run until convergence, and the approximated ground state energy is
obtained. While they showcased that this algorithm can be used for calculating
challenging chemical systems such as iron-sulfur clusters using up to 77
qubits and 6400 classical nodes, the resulting energy shows an error orders of
maginitude greater than chemical accuracy (1 kcal/mol.)

QSCI could be useful for problems where sampling from a quantum state is
classically hard or for reducing the subspace size needed in CI to reach a
certain level of accuracy. It is worth noting here that since the final
diagonalization step is carried out classically, the subspace size $R$ needs
to be within the limits of classical computation. Furthermore, its advantage
over classical SCI methods such as CIPSI and HCI is still unclear.

###  III.6 Quantum Cooling via Lindbladians

A heuristic approach for the ground state problems is to implement Lindbladian
time-evolution to simulate the cooling of an open quantum system. For a set of
jump operators $L^{a}$, and "coherent" Hamiltonian term $H$, the quantum state
is evolved by the differential equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\rho=\mathcal{L}[\rho]=-\mathrm{i}[H,\rho]+\sum_{a}L^{a}\rho(L^{a})^{\dagger}-\frac{1}{2}\\{(L^{a})^{\dagger}L^{a},\rho\\}.$
(100)

In Ref. [213] a general prescription is given, for implementing jump operators
$L^{a}$, such that the ground-state is the unique steady-state of
$\mathcal{L}$.

Jump operators $L^{a}$ are constructed as linear combinations of Heisenberg-
evolved operators. Let $A^{a}$ be proposed jump operators, and $A^{a}(s)$ be
their Heisenberg-evolution,

$\displaystyle
A^{a}(s)=\mathrm{e}^{\mathrm{i}Hs}A^{a}\mathrm{e}^{-\mathrm{i}Hs}.$ (101)

Then, $L^{a}$ is given by,

$\displaystyle{}L^{a}=\int_{-\infty}^{\infty}\mathrm{d}s\>f(s)A(s)$ (102)

for a filter-function $f(s)$, defined by

$f(s)=\frac{1}{2\pi}\int\mathrm{d}\omega\hat{f}(\omega)\mathrm{e}^{-\mathrm{i}\omega
s},$

where $\hat{f}(\omega)$ is written in frequency space. The key intuition is
that $\hat{f}(\omega)$ should be chosen to only allow cooling transitions
while preventing heating transitions. This resembles what one obtains from
ITE. Unlike ITE, this algorithm does not need to assume the non-zero overlap
between the initial state and the exact ground state. This also offers a clear
advantage over QPE for problems where preparing a state with good overlap is
challenging. As such, it should be real, non-negative, and satisfy

$\hat{f}(\omega)=0$

for $\omega\geq 0$. To ensure $L^{a}$ can be implemented efficiently, it is
crucial to choose $f(s)$ such that it is smooth and decays rapidly with the
magnitude of $s$, ensuring the filter function can be approximated by a
discretized and truncated version of the integral in Eq. (102). Then, the
linear combination of unitary techniques [214], which use ancilla qubits, can
be applied to engineer $\mathcal{L}$.

The simplest architecture requires only a single ancilla coupled to the system
(see Fig. 9) and may be ideal for near-term devices. By alternating a
controlled $A^{a}$ operation with Hamiltonian time-evolution
$e^{-iH\tau_{s}}$, the circuit implements a Trotterized approximation to the
Linbladian evolution with jump operator $L^{a}$. Ultimately, the performance
of a cooling scheme is mainly determined by the mixing time of the Linbladian,
which depends on the choice of filter function $f(s)$ and proposed jump
operators $A$. Determining fast-mixing schemes from problems of interest
remains an open problem.

Figure 9: Single-ancilla scheme for quantum cooling algorithm from Ref. 213. A
linear combination of Heisenberg-evolved operators is engineered by
alternating controlled perturbations $C_{i}$ with Hamiltonian time-evolution
$\mathrm{e}^{-\mathrm{i}Ht}$. The ancilla is periodically discarded and reset,
realizing an effective jump operator $L^{a}$ on the system, which generates
energy-lowering transitions between system eigenstates.

##  IV Quantum Algorithms for Thermal State Properties

###  IV.1 Quantum Metropolis algorithms

####  IV.1.1 Quantum Analog of Metropolis Algorithm

Temme et al. [215] first proposed a quantum analog of the classical Metropolis
sampling algorithm discussed in Section II.2.1. To generalize Metropolis
sampling to the quantum setting, one needs to consider the following two
aspects: performing the rejection step and proving that the fixed point of the
random walk is indeed the Gibbs state. For the former, the no-cloning theorem
prevents the classical solution of storing a copy of the original state and
reverting back to it in this case. This issue can be handled using the
Mariott–Watrous algorithm to rewind the state to the original [216]. For the
latter, the fixed point was proven in an idealized case, where perfect phase
estimation can be performed. When energy measurements are only given to a
finite resolution, phase estimation does not have the same performance
guarantees. A variant using a boosted shift-invariant phase estimation had
certain performance guarantees, although the shift-invariant boosting has
recently been proved impossible [120]. More specifically, it was proved that
there is no continuous family of “boosted shiftinvariant in place phase
estimation” unitaries. Although this paper marks an important advancement in
the field and served as an inspiration to many works that followed, such as
those discussed in Section IV.1.2 and Section IV.1.3, the algorithm’s
practical feasibility remains an open question.

####  IV.1.2 Quantum Enhancement of Classical Markov Chain Monte Carlo

In Section II.2.1, we briefly discussed the computational challenges
associated with classical MCMC and identified aspects where quantum computing
could potentially accelerate these processes. In this section, we review the
quantum enhancements of classical MCMC. One of the first ways of quantum
enhancement was proposed by Szegedy in 2004 to accelerate the convergence of
MCMC is incorporating a quantum walk instead of using classical state
propagation [217]. Then, in 2007, Richter reported the potential for a quantum
advantage using quantum MCMC with decoherent quantum walks by demonstrating
potential quadratic speedups in mixing times [218]. In 2008, Somma and
coworkers could apply quantum walks in MCMC to solve combinatorial
optimization problems with quadratic speedup [219]. The same year, Wocjan and
Abeyesinghe proposed a new quantum sampling method to prepare quantum states
[220]. While these works all offer quadratic speedups, more recent work in
2020 by Lemieux et al. presented an efficient implementation of Szegedy’s
quantum walk, which could offer better than quadratic speedups [221].

A different instance in which quantum computers can accelerate the convergence
of MCMC is finding the ground state (i.e., zero temperature) spin
configuration in the Ising model. For the Ising model, each spin configuration
$\textbf{s}\in\\{-1,+1\\}^{n}$ has an associated energy (Eq. (40)). The
probability of each spin configuration can also be described using Boltzmann
statistics. Sampling from this distribution is widely used, and it is often a
computational bottleneck in various algorithms, especially in the low-
temperature limit, because sampling from the Boltzmann distribution then
approximates the challenge of finding the ground state or low energy
configurations.

Recently, there have been instances of quantum algorithms for sampling and
estimating partition functions demonstrating quadratic speedup in convergence
times [222, 223]. In a work by scientists at IBM in 2023, the classical
proposal step is replaced with a quantum proposal step [224]. The spin of each
electron for any given state can be represented by a single qubit. A unitary
$U$ is applied to state s on the quantum computer under the constraint

$|\langle\textbf{s}^{\prime}|U|\textbf{s}\rangle|=|\langle\textbf{s}|U|\textbf{s}^{\prime}\rangle|,$
(103)

and each qubit is measured. Then, the Metropolis-Hastings acceptance
probability is used to accept or reject this proposal. The choice of unitary
here guarantees that
$\frac{Q(\textbf{s}|\textbf{s}^{\prime})}{Q(\textbf{s}^{\prime}|\textbf{s})}=1$,
and for a Boltzmann distribution, computation of
$\frac{\pi(\textbf{s}^{\prime})}{\pi(\textbf{s})}$ is straightforward. This
Markov chain provably converges to the Boltzmann distribution but is hard to
mimic classically. This is the source of potential quantum advantage. The
unitary $U$ is the time-evolution operator corresponding to a Hamiltonian
constructed as a sum of the classical Ising model and a term enabling quantum
transitions.

The convergence of this algorithm was probed experimentally, and they observed
approximately cubic-to-quartic improvements in $\delta$ as a function of the
number of qubits, where $\delta$ is defined in section II.2.1. While there is
no precise mathematical bound on the improvement that this algorithm provides
over its classical counterparts, the experimental demonstration seems
promising. Because only portions of this algorithm are performed on a quantum
computer with classical intermediate steps, this algorithm has lower circuit
depths than purely quantum MCMC algorithms like those discussed in Section
IV.1.3. This makes its implementation appealing in the near future. Finding
such a quantum speed-up for computational chemistry problems using the quantum
MCMC algorithm remains an open question.

####  IV.1.3 Quantum–Quantum Metropolis Algorithm

Yung and Aspuru-Guzik [225] reported a Quantum–Quantum Metropolis Algorithm
(Q2MA), which is a quantum generalization of the classical Metropolis sampling
algorithm [226], and can be used for preparing thermal states of quantum
systems with quantum speedup. The Q2MA method builds on the Markov chain
quantization method by Szegedy [217], which was formulated for classical
Hamiltonians. However, a naive implementation of Szegedy’s quantum walk for
quantum Hamiltonians is hindered by the no-cloning theorem. [227]

Yung and Aspuru-Guzik were able to circumvent this problem by preparing states
of the form $\ket{i}\equiv\ket{\psi_{i}}\ket{\tilde{\psi}_{i}}\ket{0}$, where
$\ket{\psi_{i}}$ is the eigenstate of the quantum Hamiltonian $H$ and
$\ket{\tilde{\psi}_{i}}$ is its time-reversal counterpart. Hence, both states
are eigenstates of $H$ with the same eigenvalue $E_{i}$. Leveraging this
property, a generalized Szegedy operator $W$ is constructed and used to
implement the projective measurement through quantum phase estimation. This
routine is termed quantum simulated annealing (QSA) [219]. Finally, the
thermal state can be prepared by applying QSA on an initial maximally
entangled state to yield the coherent encoding of the thermal state:

$\displaystyle\ket{\alpha_{0}}=\sum_{i=0}^{N-1}\sqrt{\mathrm{e}^{-\beta
E_{i}}/\mathcal{Z}}\ket{i}$ (104)

This formalism also falls within the thermofield double framework, as
discussed in Section II.2.4. After tracing out the ancilla qubits, this form
is equivalent to the thermal state
$\rho_{th}=(1/\mathcal{Z})\sum_{i}\mathrm{e}^{-\beta
E_{i}}\ket{\psi_{i}}\bra{\psi_{i}}$. Here,
$\mathcal{Z}=\text{Tr}(\mathrm{e}^{-\beta H})$ is the partition function.

The quadratic speedup of the Q2MA can be traced back to the
$O(1/\sqrt{\delta})$ scaling when applying $W$ for the projective measurement
in QSA, compared to the $O(1/\delta)$ scaling for classical Markov-chain
algorithms [228]. Here, $\delta$ is the gap of the transition matrix in the
Markov chain, as we mentioned in Section II.2.1. This quadratic speedup is the
main distinction between the Q2MA method and the quantum Metropolis sampling
method mentioned in Section IV.1.1. However, some limitations of this method
include that (1) the algorithm only works when the spectrum of $H$ is non-
degenerate [229] and (2) perfect quantum phase estimation is needed [120].

###  IV.2 Full Quantum Gibbs Sampler

A series of recent works have developed systematic techniques for preparing
Gibbs states on quantum computers [230, 231, 232, 120, 233]. For a Hamiltonian
$H$, which is in general composed of non-commuting terms, the Gibbs state is
defined by

$\rho_{\beta}=\mathrm{e}^{-\beta H}/\>\mathrm{Tr}[\mathrm{e}^{-\beta H}].$

One approach for preparing the Gibbs state is to engineer a Lindbladian,
$\mathcal{L}_{\beta}$ which has $\rho_{\beta}$ as the unique stationary state,

$\displaystyle{}\mathrm{e}^{\mathcal{L}_{\beta}t}[\rho_{\beta}]=\rho_{\beta}.$
(105)

In nature, Lindbladians naturally occur in the dynamics of a reduced system
with weak interaction between a small system and a large Markovian bath at
inverse temperature $\beta$. We briefly discussed this in the context of
classical Markov chain MC in Section II.2.2 and the quantum version for ground
state calculation in Section III.6. The quantum algorithm aims to use ancilla
qubits to mimic the action of the bath. However, since quantum resources are
limited, it is desirable to keep the ancilla overhead as small as possible.

It turns out that engineering a Lindbladian to satisfy Eq. (105) is
challenging. As in the design of classical Monte Carlo algorithms, if care is
not taken, then the quantum analog of detailed balance will not be achieved.
An explicit Lindbladian which satisfies quantum detailed balance for the Gibbs
state can be constructed from the Davies’ generators [234, *Davies1976Jun,
*Rivas2011Apr, *Mozgunov2020Feb]. Given a Hamiltonian $H$, with spectrum
$E_{i},|\psi_{i}\rangle$, a set of jump operators $\\{A^{a}\\}_{a\in A}$, and
Bohr frequencies $B=\\{\nu|\nu=E_{i}-E_{j}\\}$, the Davies generators are

$\mathcal{L}_{\mathrm{Davies}}(\rho)=\sum_{a\in A}\sum_{v\in
B}\gamma(v)\left(A_{v}^{a}\rho(A_{v}^{a})^{\dagger}-\frac{1}{2}\\{(A_{v}^{a})^{\dagger
A_{v}^{a}},\rho\\}\right),$ (106)

where

$A_{v}^{a}\propto\int_{-\infty}^{\infty}\mathrm{d}t\>\mathrm{e}^{-\mathrm{i}\nu
t}\mathrm{e}^{\mathrm{i}Ht}A^{a}\mathrm{e}^{-\mathrm{i}Ht}.$ (107)

Intuitively, $A^{a}_{v}$ generates transitions between pairs of eigenstates,
separated exactly by energy $\nu$. By considering an input state diagonal in
the energy basis, it becomes clear that $\mathcal{L}_{\mathrm{Davies}}$
generates a classical Markov chain on the space of energy eigenstates. The
appropriate detailed balance condition is satisfied if the function
$\gamma(\omega)$ satisfies
$\gamma(\omega)=\mathrm{e}^{-\omega\beta}\gamma(-\omega)$. This ensures that
heating transitions are penalized relative to cooling transitions, with the
ratio being determined by the target inverse temperature.

However, applying $A^{a}_{\nu}$ in a generic quantum system is challenging due
to the energy-time uncertainty principle. Conceptually, the Fourier transform
picks out transitions between $\omega,\omega^{\prime}$ that change the energy
by exactly $\nu$. As such, to approximate the Fourier transform by finite-time
evolution, the evolution time must scale inversely with energy spacings. Such
spacings are, in general, exponentially small in system size. Instead, one may
hope to use shorter-time approximations to the Davies’ generators. However,
typically, these shorter time approximations actually break the detailed
balance condition, meaning that they do not have the Gibbs state as their
stationary state. Nevertheless, one may still hope that these approximations
prepare approximate Gibbs states.

This question was analyzed in Ref. [120]. As in the ground-state cooling case,
a filter function is introduced to regularize the integral.

$\displaystyle
A^{a}_{\nu}=\int_{-\infty}^{\infty}\mathrm{d}tf(t)\mathrm{e}^{-\mathrm{i}\nu
t}\mathrm{e}^{\mathrm{i}Ht}A^{a}\mathrm{e}^{-\mathrm{i}Ht}$ (108)

Then, they bound the distance between the stationary state
$\rho_{\mathrm{ss}}$ of finite-time approximations to the Davies’ Lindbladian
and the target Gibbs state. For a Gaussian filter function, they find that the
approximation error scales as

$||\rho_{\mathrm{ss}}(\mathcal{L}_{\beta})-\rho_{\beta}||_{1}=\tilde{O}\left(\frac{\beta}{\sigma_{t}}t_{\mathrm{mix}}(\mathcal{L}_{\beta})\right).$
(109)

Therefore, to achieve a target error rate $\epsilon$, the typical width of the
Gaussian should scale as $\sigma_{t}\sim\beta t_{\mathrm{mix}}/\epsilon$. This
sets the amount of Hamiltonian simulation needed for each Lindbladian time
step. The overall scaling to prepare the Gibbs state using this approach is,
therefore, $\tilde{O}\left(\frac{\beta t_{\mathrm{mix}}}{\epsilon}\cdot
t_{\mathrm{mix}}\right).$

A subsequent breakthrough was achieved in Ref. [238], which showed how to
construct a Lindbladian that is efficiently implementable and has the quantum
Gibbs state as its steady state. The key insight was to realize that
additional coherent evolution could be used to counteract some of the unwanted
contributions coming from the dissipative part and restore detailed balance.
Since the detailed balance is exact, this reduces the need for long-time
evolution in the implementation of a single Lindbladian time-step. As such,
the resources per time-step now scale as $\tilde{O}(\beta)$. Ref. [233]
proposed a further generalization and improvement of the algorithm, by
constructing an exact Gibbs sampler using a finite number of jump operators,
simplifying implementation. Since the Linbladians used are constructed from
linear combinations of Heisenberg-evolved operators, they can be implemented
on a quantum processor using a linear combination of unitary techniques [214].
Finding a quantum advantage using this and other related algorithms is an
active research question [232].

###  IV.3 Quantum Stochastic Series Expansion

Figure 10: (a) Periodic 3-site Ising model. (b) Circuit realizing the
corresponding quantum SSE algorithm. The circuit is partitioned into $N=3$
system qubits and $n=3$ ancilla qubits. In step I, ancilla qubits are prepared
in states $|-_{q_{i}}\rangle$ via the application of Pauli X and Hadamard
gates. In step II the system qubits are prepared in a sampled state
$|a\rangle_{A}$. Step III consists of the sequential application of controlled
unitaries defined in Eq. (IV.3). Finally, in step IV, qubits are rotated back
and measured. Adapted from Ref. 239.

Similar to classical SSE (Sec. II.2.3), quantum SSE [239] evaluates a thermal
expectation value of observable $\hat{O}$,

$\langle\hat{O}\rangle_{\beta}=\frac{1}{Z}\sum_{n,a,\\{b_{i}\\}}\frac{\beta^{n}}{n!}\langle
a|\hat{H}_{b_{n}}\cdots\hat{H}_{b_{1}}|a\rangle\langle a|\hat{O}|a\rangle,$
(110)

where the tuple of the perturbation order, state, and index
$\\{(n,a,\\{b_{i}\\})\\}$, as detailed in Section II.2.3, is sampled. On
quantum computers, superpositions of states are represented naturally,
potentially enabling efficient evaluation of the pertinent matrix elements in
arbitrary bases. In particular, Ref. 239 introduces a quantum circuit (Fig.
10) evaluating these matrix elements while sampling the state-index space
classically. Therefore, on the quantum device, the no-branching condition (see
Section II.2.3) is not needed, and one may proceed to add a sufficiently large
constant to the Hamiltonian, evading negative weights.

Let us consider the quantum SSE procedure for a simple Ising model,

$\hat{H}=-\sum_{\langle ij\rangle}J_{jk}\sigma_{x}^{(j)}\sigma_{x}^{(k)}.$
(111)

where the spectrum of the Hamiltonian is shifted by a sufficiently large
number $K=\sum_{\langle jk\rangle}|J_{jk}|$. One can define an initial state
$|\Psi_{\mathrm{in}}\rangle$ on a ($N+n$)-qubit circuit,

$|\Psi_{\mathrm{in}}\rangle=|a_{A}\rangle|+_{B_{1}}\rangle\cdots|+_{B_{n}}\rangle,$
(112)

with $n$ being the order of the SSE expansion, $N$ being the number of
particles and$|+_{B_{i}}\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$. Accommodating
the additional shift of the Hamiltonian, one defines a controlled unitary

$\displaystyle U_{A,B_{i}}|a_{A}\rangle|X_{B_{i}}\rangle=$
$\displaystyle\begin{cases}\qquad\qquad\qquad\quad
1_{A}\>\>|a_{A}\rangle|X_{B_{i}}\rangle\quad&\text{if}\quad|X_{B_{i}}\rangle=|0_{B_{i}}\rangle\\\
\left[\mathrm{sgn}(h_{b})\otimes_{j=1}^{N}\sigma_{b^{j}}^{A_{j}}\right]|a_{A}\rangle|X_{B_{i}}\rangle\quad&\text{if}\quad|X_{B_{i}}\rangle=|1_{B_{i}}\rangle,\end{cases}$
(113)

which extracts the relative weight of a state-index configuration by
evaluating the estimator

$|\langle\Psi_{\mathrm{in}}|U_{A,B_{i}}|\Psi_{\mathrm{in}}\rangle|^{2}=\left|\frac{\langle
a_{A}|\hat{H}_{b_{n}}\cdots\hat{H}_{b_{1}}|a_{A}\rangle}{2^{n}h_{b_{n}}\cdots
h_{b_{1}}}\right|^{2}.$ (114)

To evaluate the matrix elements in Eq. (114) up to a target additive precision
$\epsilon$, the circuit needs to be measured a total number of
$\mathcal{O}(\frac{1}{\epsilon^{2}})$ times. The full circuit to compute the
expectation values is provided in Fig. 10. Although quantum SSE is able to
alleviate the sign problem for arbitrary Hamiltonians successfully, its
circuit depth depends on the inverse temperature $\beta$. Low temperatures
lead to non-negligible contributions from increasingly long operator strings,
limiting the applicability of quantum SSE. It is yet to be seen whether
quantum SSE offers a substantial advantage over classical SSE, as it does not
control the sign problem either.

###  IV.4 Quantum Minimally Entangled Typical Thermal States

QMETTS [240], the quantum analog of METTS (see Section II.2.4) was realized
with quantum imaginary time evolution (QITE), as shown in Fig. 11. The METTS
algorithm was originally integrated into the finite temperature density matrix
renormalization group (DMRG), which is the de facto method for low-dimensional
systems. One aims to extend the applicability of METTS beyond low-dimensional
many-body systems using QMETTS. The basic procedure of QMETTS is summarised in
Algorithm 1. A related stochastic approach was also developed for finite-
temperature dynamical observables [241].

Algorithm 1 QMETTS Algorithm

1:Start from a product state and perform imaginary time evolution (QITE) up to
time $\beta$.

2:Measure the expectation value of an operator $\hat{O}$ to obtain its thermal
average.

3:Measure a product operator, such as
$\hat{Z}^{1}\hat{Z}^{2}\cdots\hat{Z}^{N}$, to collapse back to a random
product state.

4:Repeat from step 1.

The main bottleneck in QMETTS is the QITE process, which requires quantum
resources with the quantum circuit depth exponentially increasing with system
correlation length and linearly growing with $\beta$. QITE was improved by
several algorithmic improvements: exploiting symmetries to minimize quantum
resources, optimizing circuits to decrease the depth, and implementing error-
mitigation techniques to enhance the quality of raw hardware data [241].
Recently, an adaptive variational version of QITE, using the McLachlan
variational principle [242, 243, 244, 245] for time propagation, has been
incorporated into QMETTS. This approach dynamically generates compact and
problem-specific quantum circuits, making them suitable for noisy
intermediate-scale quantum (NISQ) hardware [246]. To demonstrate a quantum
advantage of QMETTS, one needs to find a chemical example with a limited
correlation length that a classical tensor network method cannot efficiently
simulate. This question appears to be open at the time of writing.

Figure 11: Classical and quantum simulations of QMETTS for a single qubit
model with Rigetti’s QVM and Aspen-1 QPUs. Adapted from Ref. 240.

##  V Quantum Algorithms for Quantum dynamics

###  V.1 Quantum Circuit Monte Carlo

Much like most of the methods mentioned earlier, QMC suffers from the sign
problem when applied to real-time dynamics. Yang et al. proposed a new way of
mitigating the sign problem by using a hybrid quantum-classical algorithm
termed quantum-circuit Monte Carlo (QCMC) [247]. Quantum computers are used to
tackle the sign-problem-inducing parts of the time evolution, while the rest
of the algorithm is still carried out on a classical device. On the quantum
computer, the time evolution operator is implemented as a summation of unitary
operators $U(s)$ over auxiliary fields $s$:

$\displaystyle\mathrm{e}^{-\mathrm{i}H\Delta t}=\sum_{s}c(s)U(s)$ (115)

Specifically, in this work, they proposed summation formulae based on
extensions of the Lie–Trotter–Suzuki product [248, 249], including expansions
with Pauli operators and replacing the leading order with rotation operator,
with the goal of minimizing the number of gates. These formalisms can be
viewed as perturbative expansions of the Lie-Trotter-Suzuki product, which
parallels the idea of DiagMC (Sec. II.3.2). On a fault-tolerant quantum
computer, the circuit depth of the time evolution simulation scales
polynomially using this algorithm.

More specifically, they partition the real-time propagator into a product of
its first-order Trotter-Suzuki $S_{1}(\Delta\tau)$ approximation and a
correction term $V_{1}(\Delta\tau)$ according to the first-order product
formula,

$\mathrm{e}^{-\mathrm{i}H\Delta t}=V_{1}(\Delta t)S_{1}(\Delta t),$ (116)

with

$\displaystyle S_{1}(\Delta\tau)=\mathrm{e}^{-\mathrm{i}H_{M}\Delta
t}\cdots\mathrm{e}^{-\mathrm{i}H_{1}\Delta t}$ (117)

and $H=\sum_{i=1}^{M}H_{i}$, where $H_{i}$ are Hermitian. From the expression
for $S_{1}(\Delta t)$ we can readily obtain

$V_{1}(\Delta t)=\mathrm{e}^{-\mathrm{i}H\Delta
t}\mathrm{e}^{\mathrm{i}H_{M}\Delta t}\cdots\mathrm{e}^{\mathrm{i}H_{1}\Delta
t}.$ (118)

Taylor expanding each exponential in Eq. (118) and regrouping the resulting
terms leads to

$V_{1}(\Delta t)=\sum_{k=0}^{\infty}F^{(k)}_{1}(\Delta t).$ (119)

These terms $F^{(k)}_{1}(\Delta t)$ admit to the summation

$\displaystyle