# Subspace methods for electronic structure simulations on quantum computers

Mario Motta  IBM Quantum, IBM Research - Almaden, San Jose, CA 95120, USA
William Kirby  IBM Quantum, IBM Research - Cambridge, Cambridge, MA 02142, USA
Ieva Liepuoniute  IBM Quantum, IBM Research - Almaden, San Jose, CA 95120, USA
Kevin J. Sung  IBM Quantum, T. J. Watson Research Center, Yorktown Heights, NY
10598, USA  Jeffrey Cohn  IBM Quantum, IBM Research - Almaden, San Jose, CA
95120, USA  Antonio Mezzacapo  IBM Quantum, T. J. Watson Research Center,
Yorktown Heights, NY 10598, USA  Katherine Klymko  National Energy Research
Scientific Computing Center (NERSC), Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  Nam Nguyen  Integrated Vehicle Systems, Applied
Mathematics, Boeing Research & Technology, Huntington Beach, CA 92647, USA
Nobuyuki Yoshioka  Department of Applied Physics, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo, Japan  Julia E. Rice  IBM Quantum, IBM Research
- Almaden, San Jose, CA 95120, USA

###### Abstract

Quantum subspace methods (QSMs) are a class of quantum computing algorithms
where the time-independent Schrödinger equation for a quantum system is
projected onto a subspace of the underlying Hilbert space. This projection
transforms the Schrödinger equation into an eigenvalue problem determined by
measurements carried out on a quantum device. The eigenvalue problem is then
solved on a classical computer, yielding approximations to ground- and
excited-state energies and wavefunctions.

QSMs are examples of hybrid quantum-classical methods, where a quantum device
supported by classical computational resources is employed to tackle a
problem. QSMs are rapidly gaining traction as a strategy to simulate
electronic wavefunctions on quantum computers, and thus their design,
development, and application is a key research field at the interface between
quantum computation and electronic structure.

In this review, we provide a self-contained introduction to QSMs, with
emphasis on their application to the electronic structure of molecules. We
present the theoretical foundations and applications of QSMs, and we discuss
their implementation on quantum hardware, illustrating the impact of noise on
their performance.

The simulation of ground- and excited-state properties of electronic systems
is an important application of quantum computing algorithms [1, 2, 3, 4]. In
recent years, new approaches based on the projection of the Schrödinger
equation onto a subspace of the many-electron Hilbert space have been proposed
to extract Hamiltonian eigenpairs [5, 6, 7, 8, 9, 10, 11]. These methods,
which we will collectively refer to as quantum subspace methods (QSMs), are
rapidly emerging as a prominent and promising class of quantum computing
algorithms for near-term and fault-tolerant quantum devices.

This review provides a survey of the state of the art of QSMs. It is aimed at
practitioners of electronic structure (ES) interested in familiarising
themselves with quantum computation and QSMs in particular, as well as at
practitioners of quantum computation (QC) interested in the simulation of
electronic structure.

Figure 1: Structure of this review. Abbreviations indicate configuration
interaction (CI), Hartree-Fock (HF), equation of motion (EOM), multireference
CI with singles and doubles (MRCISD), quantum subspace expansion (QSE),
selected CI (SCI).

This work begins in Section I, with a review of basic concepts of ES and QC,
that readers can choose to read based on their background and expertise.
Afterward, Section II describes classical subspace methods, grouping them in
algorithms based on the notions of Krylov space, configuration interaction
(CI), and equation-of-motion (EOM), as shown in Fig. 1. This Section is
primarily destined for practitioners of QC interested in surveying subspace
methods for ES simulations on classical computers. It also provides the
notation and defines the concepts used in the remainder of the review. In
Section III, QSMs are reviewed. These are divided into algorithms based on
$k$-body fermionic operators, the Hamiltonian Krylov space, and quantum
subroutines implementing real- and imaginary-time evolution. Sections II and
III are interconnected, and discuss the computational cost (especially the
scaling of simulation on classical/quantum computers) and accuracy (especially
the convergence properties) of classical and quantum subspace methods. Section
IV presents some applications of QSMs. It illustrates the concerted use of
classical and quantum computers, for example to account for static and
dynamical electronic correlation and simulate spectral functions. Section V
focuses on the implementation of QSMs on quantum computing devices. It
emphasises the importance and the challenges posed by finite statistical
samples (shot noise) and decoherence, the economisation of quantum circuits
and measurements to contain the computational cost of QSMs, and opportunities
for the mitigation of errors affecting quantum computers based on the
formalism of QSMs. Finally, conclusions and perspectives are given in Section
VI.

##  I Basic concepts

In this Section, we provide readers with basic concepts of ES (Subsection I.1)
and QC (Subsection I.2). Readers may consult these Subsections at their
discretion, based on their background and expertise. Readers interested in ES
methods for classical computers may consult textbooks Refs. [12, 13, 14, 15,
16, 17] and reviews Refs. [18, 19, 20]. Readers interested in QC applied to ES
may consult textbooks Refs. [21, 22, 23] and reviews Refs. [1, 2, 3, 4, 24]
respectively.

###  I.1 The electronic structure problem

Our starting point is the molecular Born-Oppenheimer Hamiltonian written in
second quantisation and in a basis of orthonormal one-electron orbitals
$\\{\varphi_{p}\\}_{p=1}^{M}$,

$\hat{H}=E_{\mathrm{nuc}}+\sum_{\begin{subarray}{c}pr\\\
\sigma\end{subarray}}h_{pr}\hat{c}_{p\sigma}^{\dagger}\hat{c}_{r\sigma}^{\phantom{\dagger}}+\sum_{\begin{subarray}{c}prqs\\\
\sigma\tau\end{subarray}}\frac{(pr|qs)}{2}\hat{c}_{p\sigma}^{\dagger}\hat{c}_{q\tau}^{\dagger}\hat{c}_{s\tau}^{\phantom{\dagger}}\hat{c}_{r\sigma}^{\phantom{\dagger}}\;,$
(1)

where indices $p,r,q,s$ label spatial orbitals in a finite orthonormal basis
set, $\sigma,\tau\in\\{\uparrow,\downarrow\\}$ are spin indices, and
$\hat{c}_{p\sigma}^{\dagger}/\hat{c}_{r\sigma}^{\phantom{\dagger}}$
creates/destroys an electron in orbital $p/r$ with spin $\sigma/\tau$. The
internuclear electrostatic interaction energy is

$E_{\mathrm{nuc}}=\sum_{a<b}^{N_{nuc}}\frac{Z_{a}Z_{b}}{\|{\bf{R}}_{a}-{\bf{R}}_{b}\|}\;,$
(2)

where ${\bf{R}}_{a}$ and $Z_{a}$ are the position and atomic number of nucleus
$a$, and the coefficients

$\begin{split}h_{pr}&=\int
d{\bf{r}}\,\varphi^{*}_{p}({\bf{r}})\,\left[-\frac{1}{2}\,\frac{\partial^{2}}{\partial{\bf{r}}^{2}}-\sum_{a=1}^{N_{nuc}}\frac{Z_{a}}{\|{\bf{r}}-{\bf{R}}_{a}\|}\right]\,\varphi_{r}({\bf{r}})\;,\\\
(pr|qs)&=\int d{\bf{r}}_{1}\int
d{\bf{r}}_{2}\,\frac{\varphi^{*}_{p}({\bf{r}}_{1})\varphi_{r}({\bf{r}}_{1})\,\varphi^{*}_{q}({\bf{r}}_{2})\varphi_{s}({\bf{r}}_{2})}{\|{\bf{r}}_{1}-{\bf{r}}_{2}\|}\;,\end{split}$
(3)

describe the one-electron part of the Hamiltonian and the electron-electron
electrostatic interaction respectively. Atomic units are used throughout and
the numbers of spin-up and spin-down electrons, nuclei, and orbitals are
denoted by $N_{\uparrow}$, $N_{\downarrow}$, $N_{nuc}$ and $M$ respectively.
The electronic structure problem consists in solving for the ground (i.e. the
lowest-energy) and low-lying excited states of the Hamiltonian Eq. (1),

$\hat{H}|\Psi_{\mu}\rangle=E_{\mu}|\Psi_{\mu}\rangle\;,\;\mu=0\dots D-1\;,$
(4)

where $D$ is the dimension of the Hilbert space on which $\hat{H}$ acts (for a
system of $N_{\uparrow}$ spin-up and $N_{\downarrow}$ spin-down electrons in
$M$ orbitals, $D=\binom{M}{N_{\uparrow}}\binom{M}{N_{\downarrow}}$). In this
work, we focus on subspace methods, a family of linear variational methods to
approximately solve Eq. (4) for a quantum-chemical many-electron system.

####  I.1.1 Hartree-Fock method

The exact solution of the Schrödinger equation (4) for a molecule with more
than two electrons is a formidable problem. Within the Born-Oppenheiner
approximation, the difficulty originates from the two-electron terms of (1),
which introduce correlations in the motion of electrons under the potential
generated by nuclei. For this reason, the Schrödinger equation is solved
approximately. A simple approximation is the Hartree-Fock (HF) method, in
which the wave function is a Slater determinant

$|\Phi_{C}\rangle=\prod_{i=1}^{N_{\downarrow}}\hat{c}_{c_{i}\downarrow}^{\dagger}\prod_{i=1}^{N_{\uparrow}}\hat{c}_{c_{i}\uparrow}^{\dagger}|\mbox{\O}\rangle\;,$
(5)

where $|\mbox{\O}\rangle$ is the vacuum state, and orbitals
$|c_{i}\rangle=\sum_{p}C_{pi}|\varphi_{p}\rangle$ are determined by minimising
the energy

$E_{C}=\frac{\langle\Phi_{C}|\hat{H}|\Phi_{C}\rangle}{\langle\Phi_{C}|\Phi_{C}\rangle}\;,$
(6)

yielding the Hartree-Fock energy, $E_{\mathrm{HF}}=\min_{C}E_{C}$, and
wavefunction $|\Phi_{\mathrm{HF}}\rangle$. This procedure defines a set of
orthonormal molecular orbitals, $|c_{p}\rangle$, divided into occupied
($p\leq\max(N_{\uparrow},N_{\downarrow})$, denoted by indices $i,j,k,l$) and
virtual ($p>\max(N_{\uparrow},N_{\downarrow})$, denoted by indices $a,b,c,d$),
and a set of Slater determinants of the form

$|{\bf{x}}\rangle=\prod_{p\sigma}\left(\hat{c}_{c_{p}\sigma}^{\dagger}\right)^{x_{p\sigma}}|\mbox{\O}\rangle\;,\;x_{p\sigma}\in\\{0,1\\}^{2M}\;,\;\sum_{p}x_{p\sigma}=N_{\sigma}\;.$
(7)

We will call these Slater determinants configurations, a term that should not
be confused with configuration state functions (i.e. eigenfunctions of the
total spin operator, which are generally linear combinations of Slater
determinants). Configurations span the Hilbert space of
$(N_{\uparrow},N_{\downarrow})$ electrons in $M$ spatial orbitals, and
approximate electronic ground and excited states.

####  I.1.2 Electronic correlation

The eigenstates of Eq. (4) cannot be expressed as single determinants.
Recognising this fact, Löwdin [25] introduced the concept of electron
correlation energy, $E_{\mathrm{corr}}$, as the difference between the exact
non-relativistic ground-state energy of the molecule and $E_{\mathrm{HF}}$.
Conventionally, correlation energy is divided into static and dynamical, as
initially proposed by Sinanoglu [26]. Dynamical correlation arises when the
ground-state wavefunction,
$|\Psi_{0}\rangle=\sum_{{\bf{x}}}c_{{\bf{x}}}|{\bf{x}}\rangle$, can be
qualitatively approximated by the HF configuration,
$|c_{{\bf{x}}_{\mathrm{HF}}}|\gg|c_{{\bf{x}}}|$ for all configurations
${\bf{x}}\neq{\bf{x}}_{\mathrm{HF}}$, in presence of smaller corrections from
other configurations. Wavefunctions with this property are called single-
reference, signalling that no individual configuration mixes significantly
with the HF configuration when representing the ground-state wavefunction as a
linear combination of the form
$|\Psi_{0}\rangle=\sum_{{\bf{x}}}c_{{\bf{x}}}|{\bf{x}}\rangle$.

Static correlation, on the other hand, occurs when the ground state
wavefunction cannot be qualitatively approximated by a single Slater
determinant. Wavefunctions with this property are called multi-reference,
signalling that the Hartree-Fock configuration interacts significantly with
other low-energy configurations. Accounting for this near-degeneracy effect
requires diagonalising an appropriate secular matrix
$\langle{\bf{x}}|\hat{H}|{\bf{y}}\rangle$, where the configurations ${\bf{x}}$
and ${\bf{y}}$ arise from all possible occupations of a set of active orbitals
by a set of active electrons, i.e. an active space (see Fig. 2). For molecules
containing first- and second-row atoms, it is desirable to include all valence
orbitals (including bonding, nonbonding, and antibonding orbitals) in the
active space, and to optimize active-space orbitals to self-consistency [27,
28]. Static correlation is essential for the accurate representation of e.g. a
molecule’s dissociation into its constituent atoms. Dynamical correlation is
essential for the accurate determination of chemical properties at any
molecular geometry. This effect is important at both long range, where it
describes dispersion, and at short range. In the latter case, it is associated
with the behaviour of the wavefunction as two electrons approach each other
[29]. Properly accounting for dynamical correlation requires wavefunctions
that explicitly incorporate electronic distances or the transfer of electronic
correlation from wavefunctions to operators through canonical transformations
[30].

It should be noted that today, dynamical correlation energy usually refers to
the difference between the exact energy and a reference energy, often known as
the energy of a zeroth order reference wavefunction (within a given one-
particle basis set). For situations dominated by a single configuration, the
reference energy is thus the Hartree-Fock energy. However, for situations with
static correlation, the zeroth order reference energy corresponds to that of
the multi-reference wavefunction required to describe the low-lying electronic
states.

Figure 2: Schematic representation of an active space of 5 electrons in 4
orbitals. Spin-up/down electrons are represented by up/down-pointing arrows.
Active, inactive occupied, and inactive virtual orbitals are shown in green,
red, and blue respectively.

###  I.2 Basic quantum computing concepts

The building blocks of quantum computers are qubits: a qubit is a physical
system whose states are described by unit vectors in a two-dimensional Hilbert
space $\mathcal{H}\simeq\mathbbm{C}^{2}$. A system of $N_{q}$ qubits, also
called a register, has states described by unit vectors in the Hilbert space
$\mathcal{H}_{N_{q}}=\mathcal{H}^{\otimes N_{q}}$. An orthonormal basis of the
Hilbert space $\mathcal{H}_{N_{q}}$ is given by the following vectors, called
computational basis states,

$|{\bf{z}}\rangle=\bigotimes_{\ell=0}^{N_{q}-1}|z_{\ell}\rangle=|z_{n-1}\dots
z_{0}\rangle=|z\rangle\;,\;{\bf{z}}\in\\{0,1\\}^{N_{q}}\;,\;z=\sum_{\ell=0}^{N_{q}-1}z_{\ell}\,2^{\ell}\in\\{0\dots
2^{N_{q}-1}\\}\;.$ (8)

Starting from a register of $N_{q}$ qubits initialised in the state
$|{\bf{0}}\rangle\in\mathcal{H}_{N_{q}}$, an $N_{q}$-qubit state
$|\Phi\rangle$ can be prepared applying single- and multi-qubit unitary
transformations, or gates.

####  I.2.1 Single-qubit gates

Single-qubit Pauli operators

$\mathbbm{1}=\left(\begin{array}[]{cc}1&0\\\ 0&1\\\
\end{array}\right)\;,\;X=\left(\begin{array}[]{cc}0&1\\\ 1&0\\\
\end{array}\right)\;,\;Y=\left(\begin{array}[]{rr}0&-i\\\ i&0\\\
\end{array}\right)\;,\;Z=\left(\begin{array}[]{rr}1&0\\\ 0&-1\\\
\end{array}\right)\;,$ (9)

are very important in quantum computation. A single qubit can be prepared in a
generic state $|\Phi\rangle$ by initialisation in $|{\bf{0}}\rangle$ and
application of single-qubit Pauli rotations, i.e. unitary transformations of
the form

$R_{\sigma}\left(\theta\right)=e^{-i\frac{\theta}{2}\sigma}\;,\;\sigma\in\\{X,Y,Z\\}\;,\;\theta\in[0,2\pi)\;.$
(10)

In matrix form, with $\alpha=\cos(\theta/2)$, $\beta=\sin(\theta/2)$, and
$\gamma=e^{i\theta/2}$,

$R_{x}\left(\theta\right)=\left(\begin{array}[]{rr}\alpha&-i\beta\\\
-i\beta&\alpha\\\
\end{array}\right)\;,\;R_{y}\left(\theta\right)=\left(\begin{array}[]{rr}\alpha&-\beta\\\
\beta&\alpha\\\
\end{array}\right)\;,\;R_{z}\left(\theta\right)=\left(\begin{array}[]{rr}\gamma^{*}&0\\\
0&\gamma\\\ \end{array}\right)\;.$ (11)

It is also useful to recall that
$R_{\sigma}\left(\theta\right)=\cos(\theta/2)\mathbbm{1}-i\sin(\theta/2)\sigma$.
A generic single-qubit gate with unit determinant, $U$, can be written as

$U=R_{z}\left(\theta_{2}\right)R_{y}\left(\theta_{1}\right)R_{z}\left(\theta_{0}\right)=\left(\begin{array}[]{rr}e^{-i\frac{\theta_{0}+\theta_{2}}{2}}\cos(\theta_{1}/2)&-e^{i\frac{\theta_{0}-\theta_{2}}{2}}\sin(\theta_{1}/2)\\\
e^{-i\frac{\theta_{0}-\theta_{2}}{2}}\sin(\theta_{1}/2)&e^{i\frac{\theta_{0}+\theta_{2}}{2}}\cos(\theta_{1}/2)\\\
\end{array}\right)$ (12)

for suitable angles $\theta_{0},\theta_{1},\theta_{2}$, in what is called a
ZYZ decomposition [31] (one can think of these angles as equivalent to Euler
angles up to the double-cover of $SO(3)$ by $SU(2)$, see e.g. [32]).

For some superconducting qubit architectures, the available 1-qubit gates are
the “square-root of X” and phase gates (the latter are applied in a virtual
way by changing the phase of subsequent electromagnetic pulses [33], which
reduces both the number of gates in the circuit and the errors, since virtual
phase gates are implemented on classical software and not as physical
operations), respectively

$\sqrt{X}=e^{i\pi/4}R_{x}\left(\pi/2\right)=\left(\begin{array}[]{rr}\frac{1+i}{2}&\frac{1-i}{2}\\\
\frac{1-i}{2}&\frac{1+i}{2}\\\
\end{array}\right)\;,\;\mathsf{P}(\lambda)=e^{-i\lambda/2}R_{z}\left(\lambda\right)=\left(\begin{array}[]{rr}1&0\\\
0&e^{i\lambda}\\\ \end{array}\right)\;.$ (13)

Combining $\sqrt{X}$ and phase gates one can construct $X$ and $Y$ rotations,
e.g.
$R_{y}\left(\theta\right)=\sqrt{X}R_{z}\left(\theta\right)\sqrt{X}^{\dagger}$.
Similarly, a generic single-qubit gate with unit determinant can be written in
the following VZ decomposition, as

$\begin{split}\mathsf{U}_{3}(\theta,\phi,\lambda)&=\left(\begin{array}[]{rr}\cos(\theta/2)&-e^{i\lambda}\sin(\theta/2)\\\
e^{i\phi}\cos(\theta/2)&e^{i(\phi+\lambda)}\cos(\theta/2)\\\
\end{array}\right)\\\
&=e^{-i\frac{\pi+\theta}{2}}\mathsf{P}(\phi+\pi)\sqrt{X}\mathsf{P}(\theta+\pi)\sqrt{X}\mathsf{P}(\lambda)\;.\end{split}$
(14)

Important single-qubit gates are the Hadamard, S, and T gates, respectively

$\mathrm{Had}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{rr}1&1\\\ 1&-1\\\
\end{array}\right)\;,\;S=\mathsf{P}(\pi/2)\;,\;T=\mathsf{P}(\pi/4)\;.$ (15)

####  I.2.2 Two-qubit gates

Single-qubit gates are not sufficient for universal quantum computation, since
they cannot entangle qubits. This goal is accomplished by combining single-
and two-qubit gates. An important example is the canonical gate

$\mathsf{Can}(t_{x},t_{y},t_{z})=e^{-i\frac{\pi}{2}(t_{x}X\otimes
X+t_{y}Y\otimes Y+t_{z}Z\otimes Z)}=R_{xx}\left(\pi
t_{x}\right)R_{yy}\left(\pi t_{y}\right)R_{zz}\left(\pi t_{z}\right)\;,$ (16)

where the parameters $t_{x},t_{y},t_{z}$ lie in the so-called [34] Weyl
chamber $1/2\geq t_{x}\geq t_{y}\geq t_{z}\geq 0$ $\cup$ $1/2\geq(1-t_{x})\geq
t_{y}\geq t_{z}\geq 0$, shown in Fig. 3. Any 2-qubit gate can be decomposed
into a canonical gate and single-qubit gates by the so-called KAK
decomposition [34, 35, 36, 37, 38]. The following well-known two-qubit gates,
for example, are given (up to a global phase) by

$\begin{split}\mathsf{cNOT}&=\left(R_{z}\left(-\frac{\pi}{2}\right)R_{y}\left(-\frac{\pi}{2}\right)\otimes
R_{x}\left(-\frac{\pi}{2}\right)\right)\mathsf{Can}(1/2,0,0)\left(R_{y}\left(\frac{\pi}{2}\right)\otimes\mathbbm{1}\right)\;,\\\
\mathsf{SWAP}&=\mathsf{Can}(1/2,1/2,1/2)\;,\\\ \mathsf{iSWAP}&=(Z\otimes
Z)\mathsf{Can}(1/2,1/2,0)\;,\\\
\sqrt{\mathsf{SWAP}}&=\mathsf{Can}(1/4,1/4,1/4)\;,\\\ \end{split}$ (17)

and correspond to precise points in the Weyl chamber, shown in Fig. 3. The
cross-resonance and $\mathsf{fSim}$ gates, native to superconducting devices
with fixed and tunable frequencies respectively [39, 40, 41, 42], are given by

$\begin{split}R_{zx}\left(\theta\right)&=(\mathbbm{1}\otimes\mathrm{Had})\mathsf{Can}(\theta/\pi,0,0)(\mathbbm{1}\otimes\mathrm{Had})\;,\\\
\mathsf{fSim}(\theta,\phi)&=\left(R_{z}\left(-\frac{\phi}{2}\right)\otimes
R_{z}\left(-\frac{\phi}{2}\right)\right)\mathsf{Can}(\theta/\pi,\theta/\pi,\phi/(2\pi))\;.\\\
\end{split}$ (18)

The following $\mathsf{XX+YY}$, Givens, and controlled-phase unitaries, very
important in the simulation of fermionic systems (see Subsection I.2.6), are
given by

$\begin{split}\mathsf{V}_{\mathsf{XX+YY}}(\theta,\beta)&=(\mathbbm{1}\otimes
R_{z}\left(\beta\right))\mathsf{Can}(\theta/\pi,\theta/\pi,0)(\mathbbm{1}\otimes
R_{z}\left(-\beta\right))\;,\\\ \mathsf{G}(\theta)&=(T\otimes
T^{\dagger})\mathsf{Can}(\theta/\pi,\theta/\pi,0)(T^{\dagger}\otimes T)\;,\\\
\mathsf{cP}(\theta)&=\left(R_{z}\left(-\frac{\theta}{2}\right)\otimes
R_{z}\left(-\frac{\theta}{2}\right)\mathrm{Had}\right)\mathsf{Can}(\theta/(2\pi),0,0)(\mathbbm{1}\otimes\mathrm{Had})\;.\\\
\end{split}$ (19)

Figure 3: Location of some important 2-qubit gates in the Weyl chamber. Gates
are defined in the main text, black points indicate parameter-free gates, and
colored lines/surfaces indicate gates with one/two parameters.

####  I.2.3 Multi-qubit Pauli operators

These operators are defined as

$\hat{\sigma}_{{\bf{m}}}=\hat{\sigma}_{m_{N_{q}-1}}\otimes\dots\otimes\hat{\sigma}_{m_{0}}\;,\;\hat{\sigma}_{m}\in\\{\mathbbm{1},X,Y,Z\\}\;.$
(20)

Pauli operators are conceptually useful and important since they form a basis
for the space of linear operators on $\mathcal{H}_{N_{q}}$. This fact allows
us to represent, either exactly or approximately, multi-qubit unitaries as
products of exponentials of Pauli operators. The exponential of a Pauli
operator $P=\hat{\sigma}_{\bf{m}}$ such that $\hat{\sigma}_{m}\neq\mathbbm{1}$
for all $m$ can be applied to a register of qubits [43] by first mapping it
onto the exponential of a Pauli-$Z$ operator,

$R_{P}\left(\theta\right)=e^{-\frac{i\theta}{2}P}=\hat{V}^{\dagger}e^{-\frac{i\theta}{2}Z\otimes\dots\otimes
Z}\hat{V}\;,\;\hat{V}=\bigotimes_{\ell=0}^{N_{q}-1}\hat{A}_{m_{\ell}}\;,\;\hat{A}_{m}^{\dagger}\hat{\sigma}_{m}\hat{A}_{m}=Z\;,$
(21)

and then by introducing a ladder of $\mathsf{cNOT}$ gates that maps
$Z\otimes\dots\otimes Z$ onto $Z_{N_{q}-1}$,

$R_{P}\left(\theta\right)=\hat{V}^{\dagger}\hat{W}^{\dagger}\left(R_{z}\left(\theta\right)\right)_{N_{q}-1}\hat{W}\hat{V}\;,\;\hat{W}=\prod_{\ell=0}^{N_{q}-2}\mathsf{cNOT}_{\ell,\ell+1}\;.$
(22)

####  I.2.4 Quantum measurements

Quantum computers allow for measurement of one or more qubits in the
computational basis, Eq. (8). Measuring a register of $N_{q}$ qubits prepared
in a state $|\Phi\rangle$ a number $N_{s}$ of times yields a collection of
samples, or “shots”, $\\{{\bf{z}}_{i}\\}_{i=0}^{N_{s}-1}$, which are binary
strings of length $N_{q}$, ${\bf{z}}_{i}\in\\{0,1\\}^{N_{q}}$, statistically
independent and distributed according to
$p({\bf{z}})=|\langle{\bf{z}}|\Phi\rangle|^{2}$. The expectation value of a
diagonal operator, i.e. any operator of the form
$\hat{B}=\sum_{\bf{z}}f({\bf{z}})|{\bf{z}}\rangle\langle{\bf{z}}|$ can be
estimated as $\langle\psi|\hat{B}|\psi\rangle=\mu\pm\sigma/\sqrt{N_{s}}$ with

$\mu=\frac{1}{N_{s}}\sum_{i=0}^{N_{s}-1}f({\bf{z}}_{i})\;,\;\sigma^{2}=\frac{1}{N_{s}-1}\sum_{i=0}^{N_{s}-1}\left(f({\bf{z}}_{i})-\mu\right)^{2}\;.$
(23)

For example, $\hat{B}=Z\otimes\dots\otimes Z$ has
$f({\bf{z}})=\prod_{\ell=0}^{N_{q}-1}(-1)^{z_{\ell}}$. Non-diagonal operators
can be measured by prepending a unitary to a computational basis measurement.
For example, a generic Pauli operator can be expressed as

$\hat{\sigma}_{{\bf{m}}}=\hat{V}^{\dagger}{(Z\otimes\dots\otimes
Z)}\hat{V}=\hat{V}^{\dagger}\left[\sum_{{\bf{z}}}f({\bf{z}})\,|{\bf{z}}\rangle\langle{\bf{z}}|\right]\hat{V}\;,\;f({\bf{z}})=\prod_{\ell=0}^{N_{q}-1}(-1)^{z_{\ell}}\;,$
(24)

and its expectation value can be computed by applying Eq. (23) to measurement
outcomes drawn from the probability distribution
$p({\bf{z}})=|\langle{\bf{z}}|\hat{V}|\Phi\rangle|^{2}$. Since a generic
Hermitian operator can be written as a linear combination of Pauli operators,
its expectation value can be computed by measuring Pauli operators only.

Statistical uncertainties in measurement is a crucial aspect of quantum
computation. Quantum algorithms should be understood and formulated in terms
of random variables, and their outcomes should be accompanied by carefully
estimated statistical uncertainties. These aspects cannot be overlooked in the
implementation and design of quantum algorithms, including QSMs.

####  I.2.5 Noisy quantum devices

Noisy quantum devices are subject to decoherence (i.e., unwanted interaction
with the environment) and imperfect implementation of quantum operations
(i.e., gates and measurements). Qubits undergo relaxation and dephasing on
timescales known as $T_{1}$ and $T_{2}$ respectively which, for
superconducting qubits, are on the scale of $10^{2}$ to $10^{3}\,\mu s$.
Single-qubit gates are fast and precise operations, with duration
$T_{1q}\simeq 20\,ns$ and error $\varepsilon_{1q}\simeq 0.1\%$, whereas two-
qubit operations and measurements have duration $T_{2q}\simeq 200\,ns$,
$T_{m}\simeq 700\,ns$ and error $\varepsilon_{2q}\simeq 1\%$,
$\varepsilon_{m}\simeq 1\%$ respectively [44, 45]. The qubit coherence times
define the timescales over which they lose quantum information to decoherence,
and the gate times and errors define the duration and accuracy of a
computation: in the absence of error correction, the duration of a computation
should not exceed the qubit coherence time, and the number of gates should be
such that the accumulation of error does not prevent algorithms from yielding
accurate results.

It should also be noted that many quantum computer architectures have limited
qubit connectivity, i.e. two-qubit gates can only be applied across certain
pairs of qubits. This fact conflicts with the quantum circuit model which
allows general 2-qubit interactions and thus implicitly assumes a completely
connected network of qubits. To remedy this, Beals et al [46] demonstrated the
use of $\mathsf{SWAP}$ gates, which can be inserted into the quantum circuit
to enable it to be executed on a quantum computer with limited qubit
connectivity (see Fig. 5b for an example). An important optimization problem
is how to insert these gates to run a quantum circuit while minimising the
adverse effect of these additional gates on the performance of the algorithm
(runtime and error) [47, 48, 49, 50, 51].

Two parameters, respectively called width and depth, are often used to
characterise the cost of a quantum circuit. Width refers to the number of
qubits that comprise the circuit, $N_{q}$. Depth is the maximal length $d$ of
a path from the input (qubit initialization) to the output (measurement
operation) of a circuit. Width and depth are both limiting factors in the
execution of quantum algorithms (large width corresponds to many qubits, and
large depth to many operations carried out sequentially), and depth is an
important computational bottleneck when the physical duration of a circuit,
roughly approximated by $dT_{2q}$, is comparable with the coherence time of a
qubit.

Techniques for the mitigation of readout [52, 53] and gate [54, 55, 56] errors
exist, which alleviate these bottlenecks and, along with continued improvement
in device manufacturing and operation and research in algorithm design and
refinement, are fundamental activities towards practical quantum computation.

####  I.2.6 Qubit mappings for fermions

The Fock space of electrons occupying $M$ spatial orbitals has the same
dimension, $4^{M}$, as the Hilbert space of $N_{q}=2M$ qubits. Therefore, it
is possible to construct a one-to-one mapping between the two spaces. There
are combinatorially many ways to map a quantum system to a set of qubits [57,
58] and, since fermions exhibit non-locality of their state space, due to
their antisymmetric exchange statistics, any representation of fermionic
systems on collections of qubits must introduce non-local structures [59].
Here, we consider the Jordan-Wigner (JW) transformation [60, 61, 62, 63], that
maps electronic configurations with generic particle number, Eq. (7), onto
computational basis states, Eq. (8),

$\big{(}\hat{c}_{M-1,\downarrow}^{\dagger}\big{)}^{x_{M-1,\downarrow}}\dots\big{(}\hat{c}_{0,\uparrow}^{\dagger}\big{)}^{x_{0,\uparrow}}|\mbox{\O}\rangle\mapsto|{\bf{x}}\rangle\;,$
(25)

and fermionic creation and annihilation operators ($\hat{c}_{p}^{\dagger}$ and
$\hat{c}_{p}^{\phantom{\dagger}}$ respectively) onto non-local qubit operators
of the form

$\begin{split}\hat{c}_{p\sigma}^{\dagger}\mapsto\frac{X_{p\sigma}-iY_{p\sigma}}{2}\otimes
Z_{p\sigma-1}\otimes\dots\otimes Z_{0}\equiv
S^{(+)}_{p\sigma}\,Z^{p\sigma-1}_{0}\;,\\\
\hat{c}_{p\sigma}^{\phantom{\dagger}}\mapsto\frac{X_{p\sigma}+iY_{p\sigma}}{2}\otimes
Z_{p\sigma-1}\otimes\dots\otimes Z_{0}\equiv
S^{(-)}_{p\sigma}\,Z^{p\sigma-1}_{0}\;.\end{split}$ (26)

In Eqs. (25) and (26) we used $Z^{p}_{r}=\otimes_{l=r}^{p}Z_{l}$ to denote a
tensor product of Pauli $Z$ operators acting on qubits $r$ to $p$, and we used
the following notation for qubit indices,

$p\sigma=\left\\{\begin{array}[]{ll}p&\sigma\,=\,\uparrow\\\
p+M&\sigma\,=\,\downarrow\\\ \end{array}\right.$ (27)

The non-locality of these operators is required to preserve canonical
anticommutation relations between creation and destruction operators. The main
limitation of the JW transformation is that the qubit representation of
$\hat{c}_{p\sigma}^{\dagger}$ has $\mathcal{O}(M)$ non-locality [64, 65],
which immediately translates to $k$-body fermionic operators. Another
limitation is that, since the JW transformation operates in the Fock space
(i.e. it allows the description of states with any particle number, spin, and
point group symmetry), on noisy quantum devices the conservation of particle
number, spin, and point group symmetry (an important requirement in typical ES
calculations) is not guaranteed.

The qubit representation of Eq. (1) can be derived by using Eq. (26) to
represent each product of creation and destruction operators as a linear
combination of Pauli operators. For the one-body part
$\hat{H}_{1}=\sum_{pr,\sigma}h_{pr}\hat{c}_{p\sigma}^{\dagger}\hat{c}_{r\sigma}^{\phantom{\dagger}}$,
assuming real-valued coefficients, one has

$\hat{H}_{1}\to\sum_{p\sigma}\frac{h_{pp}}{2}(1-Z_{p\sigma})+\sum_{p<r,\sigma}\frac{h_{pr}}{2}Z^{p\sigma-1}_{r\sigma+1}(X_{p\sigma}X_{r\sigma}+Y_{p\sigma}Y_{r\sigma})=\sum_{\alpha}c_{\alpha}\sigma_{{\bf{m}}_{\alpha}}\;,$
(28)

where we used the notation $(p\uparrow)=p$, $(p\downarrow)=p+M$ for qubit
indices. This expression, combined with Eq. (23), immediately shows that
$\hat{H}_{1}$ can be measured Pauli-by-Pauli. A step of time evolution under
$\hat{H}_{1}$ can similarly be approximated with e.g. a Trotter product
formula,

$e^{-i\Delta t\hat{H}_{1}}\to\prod_{\alpha}e^{-i\Delta
tP_{\alpha}}+\mathcal{O}(\Delta t^{2})\;.$ (29)

For the two-body part, one has a similar but more complicated expression [65].

Figure 4: Left (a): Implementation of a change-of-basis unitary $e^{\hat{K}}$
as a circuit $W_{\kappa}$ comprising $2M$ $R_{z}$ gates and $M(M-1)$
$\mathsf{V_{XX+YY}}$ gates arranged in $M$ layers (marked as
$R_{z}/\mathsf{V}$ and omitting parameters to avoid clutter), illustrated for
a system of $M=4$ spatial orbitals. Center (b): measurement of $\hat{H}_{1}$
using the circuit $W_{\kappa}$ (teal block) and a computational basis
measurement (red meter symbols). Right (c): time evolution under $\hat{H}_{1}$
using the circuit $W_{\kappa}$ and a single layer of single-qubit $Z$
rotations (purple blocks) with ${\bf{x}}$ defined in Eq. (32).

Alternatively, $\hat{H}_{1}$ can be diagonalised by the exponential of a one-
body operator,

$\hat{H}_{1}=e^{-\hat{K}}\left(\sum_{p\sigma}\eta_{p}\,\hat{c}_{p\sigma}^{\dagger}\hat{c}_{p\sigma}^{\phantom{\dagger}}\right)e^{\hat{K}}\;,\;\hat{K}=\sum_{p<r\sigma}\kappa_{pr}(\hat{c}_{p\sigma}^{\dagger}\hat{c}_{r\sigma}^{\phantom{\dagger}}-\hat{c}_{r\sigma}^{\dagger}\hat{c}_{p\sigma}^{\phantom{\dagger}})\;,$
(30)

where $\eta_{p}$ are the eigenvalues of the matrix with elements $h_{pr}$. The
qubit representation of $\hat{H}_{1}$ is

$\hat{H}_{1}\to
W_{\kappa}^{\dagger}\left(\sum_{p\sigma}\frac{\eta_{p}}{2}(1-Z_{p\sigma})\right)W_{\kappa}\;,$
(31)

where $W_{\kappa}$ is a unitary that can be implemented with a quantum circuit
(see Fig. 4a) containing $\mathcal{O}(M^{2})$ $\mathsf{XX+YY}$ gates,
$\mathcal{O}(M)$ depth, and requiring linear qubit connectivity only [66, 67,
68, 69, 70]. As a consequence, $\hat{H}_{1}$ can be measured as shown in Fig.
4b and time evolution can then be exactly implemented up to a global phase
(see Fig. 4c) as

$e^{-it\hat{H}_{1}}\to
W_{\kappa}^{\dagger}\,U_{1}(t\boldsymbol{\eta})\,W_{\kappa}\;,\;U_{1}({\bf{x}})=\bigotimes_{p\sigma}\big{(}R_{z}\left(-x_{p}\right)\big{)}_{p\sigma}\;.$
(32)

Figure 5: Top (a): quantum circuit implementing a step of time evolution under
the ES Hamiltonian, with a Trotter product formula based on a low-rank
decomposition of the two-body part, Eq. (37). The circuits $W_{\kappa}$ and
$U_{1}\big{(}{\bf{x}}\big{)}$ are defined as in Fig. 4, and the circuit
$U_{2}\big{(}{\bf{x}}\big{)}$ is shown in the bottom panel (b) for $M=4$
spatial orbitals, with purple blocks labeling single-qubit $Z$ rotations, and
green/blue/red blocks connected by vertical black lines labeling controlled-
phase rotations implementing
$\uparrow,\uparrow$/$\downarrow,\downarrow$/$\uparrow,\downarrow$ terms in Eq.
(37). Note the use of a $\mathsf{SWAP}$ network to implement two-qubit gates
acting on distant qubits assuming linear device connectivity only, and the
fact that the $\mathsf{SWAP}$ network inverts the qubit order.

Similarly, the two-body part of the Hamiltonian can be mapped on a qubit
operator using a low-rank approximation [69]. One can write the electron
repulsion integral (ERI) with a density fitting [71, 72, 73, 74, 75, 76, 77,
78, 79] or Cholesky decomposition [80, 81, 82, 83, 84, 85, 86] as

$\frac{(pr|qs)}{2}\simeq\sum_{\gamma=1}^{N_{\gamma}}L^{\gamma}_{pr}L^{\gamma}_{qs}\;,\;N_{\gamma}=\mathcal{O}(M)\;,$
(33)

and obtain

$\hat{H}=E_{\mathrm{nuc}}+\hat{J}_{1}+\sum_{\gamma=1}^{N_{\gamma}}\hat{L}_{\gamma}^{2}\;,\;\hat{L}_{\gamma}=\sum_{pr,\sigma}L^{\gamma}_{pr}\hat{c}_{p\sigma}^{\dagger}\hat{c}_{r\sigma}^{\phantom{\dagger}}\;,$
(34)

where
$\hat{J}_{1}=\sum_{pr,\sigma}\big{(}h_{pr}-\sum_{q\gamma}L^{\gamma}_{pq}L^{\gamma}_{qr}\big{)}\hat{c}_{p\sigma}^{\dagger}\hat{c}_{r\sigma}^{\phantom{\dagger}}$
and each $\hat{L}_{\gamma}$ is real-valued and symmetric due to the 8-fold
symmetry of the ERI. The qubit representation of Eq. (34) is

$\hat{H}\to
E_{\mathrm{nuc}}+W^{\dagger}_{\kappa_{0}}\Big{(}\sum_{p\sigma}\zeta_{p}\,n_{p\sigma}\Big{)}W_{\kappa_{0}}+\sum_{\gamma=1}^{N_{\gamma}}W_{\kappa_{\gamma}}^{\dagger}\Big{(}\sum_{\begin{subarray}{c}pq\\\
\sigma\tau\end{subarray}}\lambda_{p}^{\gamma}\lambda^{\gamma}_{q}\,n_{p\sigma}n_{q\tau}\Big{)}W_{\kappa_{\gamma}}\;,\;$
(35)

where $n_{p\sigma}=(1-Z_{p\sigma})/2$, $\zeta_{p}$ ($\lambda_{p}^{\gamma}$)
are the eigenvalues of $\hat{J}_{1}$ ($\hat{L}_{\gamma}$) and $W_{\kappa_{0}}$
($W_{\kappa_{\gamma}}$) is a circuit representation of the unitary
diagonalising $\hat{J}_{1}$ ($\hat{L}_{\gamma}$). The operator Eq. (35) can be
measured term-by-term, and a step of time evolution under $\hat{H}$ can be
approximated with e.g. a Trotter product formula,

$e^{-i\Delta
t\hat{H}}\to\prod_{\gamma=1}^{N_{\gamma}}\left(W^{\dagger}_{\kappa_{\gamma}}U_{2}(\sqrt{\Delta
t}\,\boldsymbol{\lambda}^{\gamma})W_{\kappa_{\gamma}}\right)W^{\dagger}_{\kappa_{0}}U_{1}(\Delta
t\boldsymbol{\zeta})W_{\kappa_{0}}+\mathcal{O}(\Delta t^{2})\;.$ (36)

The circuit in Eq. (36) comprises $N_{\gamma}+2$ change-of-basis circuits, a
layer of single-qubit $Z$ rotations, and $N_{\gamma}$ circuits of the form

$U_{2}(\boldsymbol{x})=\prod_{pq}\left(\mathsf{cP}_{2x_{p}x_{q}}\right)_{p\uparrow,q\downarrow}\prod_{\begin{subarray}{c}p<q\\\
\sigma\end{subarray}}\left(\mathsf{cP}_{2x_{p}x_{q}}\right)_{p\sigma,q\sigma}U_{1}(\boldsymbol{x}^{2})\;.$
(37)

The circuit $U_{2}(\boldsymbol{x})$ has $M^{2}+M(M-1)$ $\mathsf{cP}$ gates and
$2M$ single-qubit $Z$ rotations, depth $2M+1$, and requires all-to-all qubit
connectivity, see Fig. 5.

##  II Classical subspace methods

Subspace methods construct finite-dimensional approximations to the eigenvalue
equation Eq. (4), an approach called Galerkin’s method [87]. Specifically,
given a set of many-electron wavefunctions ${\bf{v}}_{0}\dots\bf{v}_{n-1}$, we
seek a function of the form
$|\tilde{\Psi}^{(n)}_{\mu}\rangle=\sum_{\beta=0}^{n-1}C_{\beta\mu}|{\bf{v}}_{\beta}\rangle$
that solves the projection of the Schrödinger equation Eq. (4) onto the
subspace $V_{n}=\mathrm{Span}({\bf{v}}_{0}\dots\bf{v}_{n-1})$, i.e., for each
$\alpha,\mu\in[n]\coloneqq\\{0,1,...,n-1\\}$,

$\sum_{\beta=0}^{n-1}\langle{\bf{v}}_{\alpha}|\hat{H}|{\bf{v}}_{\beta}\rangle\;C_{\beta\mu}=\tilde{E}^{(n)}_{\mu}\sum_{\beta=0}^{n-1}\langle{\bf{v}}_{\alpha}|{\bf{v}}_{\beta}\rangle\;C_{\beta\mu}\;.$
(38)

The generalised eigenvalue equation (GEEV) (38) will be compactly written as
$HC=SC\tilde{E}^{(n)}$, where $H$ and $S$ denote the Hamiltonian and overlap
matrices, respectively, whose entries are

$H_{\alpha\beta}\coloneqq\langle{\bf{v}}_{\alpha}|\hat{H}|{\bf{v}}_{\beta}\rangle,\quad
S_{\alpha\beta}\coloneqq\langle{\bf{v}}_{\alpha}|{\bf{v}}_{\beta}\rangle.$
(39)

The approximate eigenvalues $\tilde{E}^{(n)}_{\mu}$ and eigenvectors
$|\tilde{\Psi}^{(n)}_{\mu}\rangle$ differ from the exact ones since the
approximate eigenvectors are restricted to the subspace $V_{n}$. The quality
of a subspace method depends on several factors: (i) the nature of the basis
vectors $|{\bf{v}}_{\beta}\rangle$, (ii) the dimension $n$ of the employed
subspace, (iii) the cost (and associated error) of computing the matrices $H$,
$S$ in Eq. (39), and (iv) the numeric stability of Eq. (38), connected with
the condition number of the overlap matrix $S$.

###  II.1 Krylov space methods

A Krylov space is associated to a matrix-vector pair $(A,{\bf{v}}_{0})$, and
is the subspace spanned by powers of $A$ applied to ${\bf{v}}_{0}$,
$\mathrm{K}_{n}(A,{\bf{v}}_{0})=\mathrm{Span}({\bf{v}}_{0}\dots
A^{n-1}{\bf{v}}_{0})$. In the context of classical numerical eigensolvers, $A$
is often chosen to be the matrix whose spectrum one wishes to study [88]. To
study the low-energy eigenpairs of a Hamiltonian, one applies powers of the
Hamiltonian to a trial vector [89],
$|{\bf{v}}_{\alpha}\rangle=\hat{H}^{\alpha}|{\bf{v}}_{0}\rangle$,
$\alpha=0\dots n-1$. We will often call this ubspace the “Hamiltonian Krylov
space”. The overlap and Hamiltonian matrices are
${S_{\alpha\beta}=\langle{\bf{v}}_{0}|\hat{H}^{\alpha+\beta}|{\bf{v}}_{0}\rangle=f_{\alpha+\beta}}$
and
${H_{\alpha\beta}=\langle{\bf{v}}_{0}|\hat{H}^{\alpha+\beta+1}|{\bf{v}}_{0}\rangle=f_{\alpha+\beta+1}}$,
where $f_{\ell}=\langle{\bf{v}}_{0}|\hat{H}^{\ell}|{\bf{v}}_{0}\rangle$ and
$\ell=0\dots 2n-1$.

An intuitive justification of the Krylov space comes from the notion of
imaginary-time evolution (ITE): for a Hamiltonian with a non-degenerate ground
state $|\Psi_{0}\rangle$, any trial state $|{\bf{v}}_{0}\rangle$ non-
orthogonal to $|\Psi_{0}\rangle$ is mapped onto the ground state by ITE for a
sufficiently long time $\tau>0$,

$|\Phi_{\tau}\rangle=\frac{e^{-\tau\hat{H}}|{\bf{v}}_{0}\rangle}{\|e^{-\tau\hat{H}}{\bf{v}}_{0}\|}=\frac{\sum_{\mu}e^{-\tau(E_{\mu}-E_{0})}|\Psi_{\mu}\rangle\langle\Psi_{\mu}|{\bf{v}}_{0}\rangle}{\sqrt{\sum_{\mu}e^{-2\tau(E_{\mu}-E_{0})}|\langle\Psi_{\mu}|{\bf{v}}_{0}\rangle|^{2}}}\;,\;\lim_{\tau\to\infty}|\Phi_{\tau}\rangle=|\Psi_{0}\rangle\;.$
(40)

As $e^{-\tau\hat{H}}$ can be approximated by a truncated Taylor series,
${e^{-\tau\hat{H}}\simeq\sum_{\alpha=0}^{n-1}\frac{(-\tau)^{\alpha}}{\alpha!}\hat{H}^{\alpha}}$,
the state $|\Phi_{\tau}\rangle$ can be approximated by a linear combination of
vectors $\hat{H}^{\alpha}|{\bf{v}}_{0}\rangle$, i.e. by a vector in the
Hamiltonian Krylov space.

Eq. (40) suggests that the quality of the $n$-dimensional Hamiltonian Krylov
space is affected by several factors, such as the overlap between the initial
state ${\bf{v}}_{0}$ and the ground state $\Psi_{0}$,
$\langle\Psi_{0}|{\bf{v}}_{0}\rangle$, and the spectral gap, $E_{1}-E_{0}$,
whose inverse is the decay rate of excited-state contributions in Eq. (40). A
more formal result is the Kaniel-Paige inequality [90, 91],

$0\leq\tilde{E}_{0}^{(n)}-E_{0}\leq(E_{D-1}-E_{0})\left[\frac{\tan\theta(\Psi_{0},{\bf{v}}_{0})}{T_{n-1}(\gamma_{0})}\right]^{2}\;,$
(41)

where $E_{0}<\dots<E_{D-1}$ and $\Psi_{0}\dots\Psi_{D-1}$ are the eigenpairs
of $\hat{H}$, $\theta(\Psi_{0},{\bf{v}}_{0})$ is the angle between $\Psi_{0}$
and ${\bf{v}}_{0}$, $T_{n-1}$ is the $(n-1)$-th Chebyshev polynomial, and
$\gamma_{0}=1+2\frac{E_{1}-E_{0}}{E_{D-1}-E_{1}}$ is related to the spectral
gap of $\hat{H}$. The left-hand side of Eq. (41) shows that the Krylov method
produces an upper bound for the ground-state eigenvalue. On the right-hand
side, $\tan$ shows that the bound is tighter when ${\bf{v}}_{0}$ is closer to
$\Psi_{0}$. Furthermore, as $\lim_{n\to\infty}T_{n-1}(1+x)=\infty$ for all
$x>0$, increasing the dimension of the Krylov space results in a progressively
more accurate estimation of $E_{0}$. This property holds provided that the
spectral gap is positive. Expanding this analysis, the Saad inequality [92]
quantifies the accuracy of the eigenvalue approximations yielded by a Krylov
space for excited states. For a generic eigenvector $\Psi_{\mu}$ non-
orthogonal to ${\bf{v}}_{0}$,

$0\leq\tilde{E}_{\mu}^{(n)}-E_{\mu}\leq(E_{D-1}-E_{\mu})\left[\frac{L_{\mu}^{(n)}\tan\theta(\Psi_{\mu},{\bf{v}}_{0})}{T_{n-1-\mu}(\gamma_{\mu})}\right]^{2}\;,$
(42)

where $\gamma_{\mu}=1+\frac{E_{\mu+1}-E_{\mu}}{E_{D-1}-E_{\mu+1}}$, and
$L_{\mu}^{(n)}=\prod_{\nu<\mu}\frac{E_{D-1}-\tilde{E}^{(n)}_{\nu}}{E_{\mu}-\tilde{E}^{(n)}_{\nu}}$.
Not only do these results show that the best approximate eigenvalues
represented in the Krylov space converge, they converge extremely quickly:
since ${T_{n}(\gamma)\geq\frac{1}{2}(\gamma+\sqrt{\gamma^{2}-1})^{n}}$ for
$\gamma\geq 1$, the upper bounds in Eqs. (41) and (42) converge exponentially
with the Krylov space dimension provided the corresponding gap condition
holds.

From a numerical standpoint, the Krylov method has an important limitation:
the condition number of the overlap matrix grows exponentially with subspace
dimension $n$. This fact is expressed by the Beckermann-Townsend inequality
[93] stating that, for any two integers $j,k$ such that $0\leq j+2k\leq n-1$,

$\frac{\sigma_{j}(S)}{\sigma_{j+2k}(S)}\geq\frac{1}{4}\exp\left[\frac{\pi^{2}}{4\log\left(\frac{4(n-n\%2)}{\pi}\right)}\right]^{k-n\%2}\;,$
(43)

where $\sigma_{\ell}(S)$ is the $\ell$-th singular value of the overlap matrix
$S$ and $\%$ denotes modulo 2. Choosing $n$ odd, $j=0$, and $2k=n-1$, shows
that the condition number is lower-bounded by a quantity,

$\mathrm{cond}(S)=\frac{\sigma_{0}(S)}{\sigma_{n-1}(S)}\geq\frac{1}{4}\exp\left[\frac{\pi^{2}}{4\log\left(\frac{4(n-1)}{\pi}\right)}\right]^{\frac{n-1}{2}-1}\;,$
(44)

that diverges very rapidly with $n$. This may be intuitively understood as a
trade-off for achieving the exponential convergence of the lowest energy in
the subspace toward the true ground-state energy. However, as we will see in
Section II.1.1, this ill-conditioning can be largely mitigated by constructing
an orthonormal basis for the Krylov space.

####  II.1.1 Lanczos method

Lanczos is a specific variant of Krylov subspace methods. As $n$ increases,
Krylov vectors tend to become almost linearly dependent, see the divergence of
$\mathrm{cond}(S)$ in Eq. (44), leading to ill-conditioning. In the Lanczos
method [94], the objective is to create an orthonormal basis
$|{\bf{q}}_{\alpha}\rangle$ for the Hamiltonian Krylov space, such that
$H_{\alpha\beta}=\langle{\bf{q}}_{\alpha}|\hat{H}|{\bf{q}}_{\beta}\rangle$ is
tridiagonal. More specifically, one constructs Lanczos vectors
$|{\bf{q}}_{\alpha}\rangle$ according to

$\begin{split}|{\bf{q}}_{0}\rangle&=|{\bf{v}}_{0}\rangle\;,\\\
|{\bf{q}}_{1}\rangle&\propto\Big{(}\mathbbm{1}-|{\bf{q}}_{0}\rangle\langle{\bf{q}}_{0}|\Big{)}\hat{H}|{\bf{q}}_{0}\rangle\;,\\\
|{\bf{q}}_{\alpha+1}\rangle&\propto\Big{(}\mathbbm{1}-|{\bf{q}}_{\alpha}\rangle\langle{\bf{q}}_{\alpha}|-|{\bf{q}}_{\alpha-1}\rangle\langle{\bf{q}}_{\alpha-1}|\Big{)}\hat{H}|{\bf{q}}_{\alpha}\rangle\;,\\\
\end{split}$ (45)

where the proportionality symbols indicate that the
$|{\bf{q}}_{\alpha}\rangle$ should be normalized. We can see that the
resulting $H$ is tridiagonal, and the basis is orthogonal, as follows. By
construction, $\hat{H}|{\bf{q}}_{\alpha}\rangle$ is in the span of
$\\{|{\bf{q}}_{\alpha+1}\rangle,|{\bf{q}}_{\alpha}\rangle,|{\bf{q}}_{\alpha-1}\rangle,...,|{\bf{q}}_{0}\rangle\\}$.
Assuming as an inductive hypothesis that
$\\{|{\bf{q}}_{\alpha}\rangle,|{\bf{q}}_{\alpha-1}\rangle,...,|{\bf{q}}_{0}\rangle\\}$
are orthogonal, an equivalent statement is that

$H_{\beta\gamma}=\langle{\bf{q}}_{\beta}|\hat{H}|{\bf{q}}_{\gamma}\rangle$
(46)

is upper-Hessenberg up to dimension $\alpha+1$, i.e. $H_{\beta\gamma}=0$
whenever $\beta>\gamma+1$. However, $H$ is also Hermitian (since $\hat{H}$
is), which implies that it is tridiagonal in the basis of
$|{\bf{q}}_{\alpha}\rangle$. Hence, $\hat{H}|{\bf{q}}_{\alpha}\rangle$ is
actually guaranteed to lie in the span of
$\\{|{\bf{q}}_{\alpha+1}\rangle,|{\bf{q}}_{\alpha}\rangle,|{\bf{q}}_{\alpha-1}\rangle\\}$,
so the partial orthogonalization in the last line of Eq. (45) is in fact
sufficient to orthogonalize with respect to all previous basis vectors.

Finally, it can be proved by induction over $n$ that the Lanczos and Krylov
bases span identical subspaces. While the vectors $|{\bf{q}}_{\alpha}\rangle$
are orthonormal assuming arithmetic operations are carried out exactly,
computational simulations use floating-point arithmetic, which results in loss
of orthonormality and spurious eigenpairs [95, 88]. Practical implementations
of the Lanczos algorithm mitigate its numerical instability by preventing
orthogonality loss through repeated re-orthogonalisation of each newly
generated vector with all the previously generated ones [96].

####  II.1.2 Davidson method

The standard diagonalisation algorithm in classical electronic structure is
the Davidson method [97, 98, 99]. Unlike the Krylov subspace method, the
Davidson method iteratively extends an $n$-dimensional subspace
$V_{n}=\mathrm{Span}[{\bf{v}}_{0}\dots\bf{v}_{n-1}]$ by adding a vector
$|{\bf{w}}\rangle$ to its basis. This vector is chosen such that
$\langle{\bf{w}}|{\bf{v}}_{\alpha}\rangle=0$ and the approximation for the
ground state of $\hat{H}$, restricted to the subspace
$V_{n+1}=\mathrm{Span}[{\bf{v}}_{0}\dots\bf{v}_{n-1}{\bf{w}}]$, is as accurate
as possible. The extension of $V_{n}$ is guided by a linear transformation
$\hat{P}$ called a preconditioner. A simple and widespread example is the
Jacobi or diagonal preconditioner [100], which is efficient for diagonally-
dominant matrices. Given the best approximation
$|\tilde{\Psi}^{(n)}_{0}\rangle$ to the ground state in $V_{n}$ and the
corresponding eigenvalue $\tilde{E}^{(n)}_{0}$, let us assume that one seeks
an eigenstate of the form
$|{\bf{w}}\rangle=\alpha|\tilde{\Psi}^{(n)}_{0}\rangle+\beta|{\bf{e}}_{\ell}\rangle$,
where $|{\bf{e}}_{\ell}\rangle$ is the $\ell$-th element of the canonical
basis. The coefficients $\alpha,\beta$, which are the solution of a simple
$2\times 2$ eigenvalue equation, can be expanded to first order in
$\langle{\bf{e}}_{\ell}|\tilde{\Psi}^{(n)}_{0}\rangle$ to give

$|{\bf{w}}\rangle=|\tilde{\Psi}^{(n)}_{0}\rangle+|{\bf{e}}_{\ell}\rangle\langle{\bf{e}}_{\ell}|\hat{P}|{\bf{r}}\rangle\;,\;|{\bf{r}}\rangle=\left[\hat{H}-\tilde{E}^{(n)}_{0}\right]|\tilde{\Psi}^{(n)}_{0}\rangle\;,\;\hat{P}=\sum_{\ell}\frac{|{\bf{e}}_{\ell}\rangle\langle{\bf{e}}_{\ell}|}{\langle{\bf{e}}_{\ell}|\hat{H}|{\bf{e}}_{\ell}\rangle-\tilde{E}^{(n)}_{0}}\;.$
(47)

The vector $|{\bf{r}}\rangle$ and the operator $\hat{P}$ are called the
residue and the diagonal preconditioner respectively. In the Davidson method
with the Jacobi preconditioner, one applies the preconditioner to the residual
vector, yielding
$|{\bf{w}}\rangle=|\tilde{\Psi}^{(n)}_{0}\rangle+\hat{P}|{\bf{r}}\rangle$, and
then expands the subspace by orthonormalising $|{\bf{w}}\rangle$ against
${\bf{v}}_{0}\dots\bf{v}_{n-1}$. There exist alternatives to the diagonal
preconditioner illustrated here, which are suited to non-diagonally-dominant
matrices, such as the “pspace” preconditioner of Olsen et al [101].

###  II.2 Configuration interaction (CI) methods

CI methods are a type of subspace method for ES calculations. CI methods
operate within a subspace of the Hilbert space spanned by a collection of
electronic configurations, i.e. Slater determinants of the form Eq. (7).
Interaction means constructing linear combinations,
$|\tilde{\Psi}_{\mu}\rangle=\sum_{\alpha\in
A}C_{\alpha\mu}|{\bf{x}}_{\alpha}\rangle$, of electronic configurations
${\bf{x}}_{\alpha}$ drawn from a set $A$, which may be predefined (as in
truncated CI approaches) or adaptively constructed (as in selected CI
approaches).

####  II.2.1 Full and truncated CI

If the set $A$ includes all possible configurations of the appropriate
symmetry, the resulting method is called full CI (FCI) and exactly solves the
electronic Schrödinger equation, Eq. (4). Exact diagonalizaton in an active
space is called complete active space (CAS) CI, and CASCI with variational
optimization of the active-space orbitals is called CASSCF. The FCI method has
many attractive characteristics, especially that it exactly solves the
Schrödinger equation. However, the computational cost of FCI grows
combinatorially with the numbers of electrons and orbitals,
$|A|=D=\binom{M}{N_{\uparrow}}\binom{M}{N_{\downarrow}}$ and thus,
notwithstanding considerable progress, an exact FCI approach is feasible only
for relatively small basis sets [99, 102, 103, 104, 105]. Motivated by this
observation, attempts were made to search for approximate FCI approaches by a
judicious choice of the configurations [106, 107, 108].

The first term in the expansion of the exact ground state onto configurations,
$|\Psi_{\mu}\rangle=\sum_{{\bf{x}}}C_{{\bf{x}}\mu}|{\bf{x}}\rangle$ where
$|{\bf{x}}\rangle$ denotes a configuration as in Eq. (7), is often the HF
state. The other terms can be characterised by the number of electronic
excitations from occupied to virtual orbitals in the HF state (singles,
doubles, triples, quadruples, etc). The CI space may be truncated by retaining
configurations with a finite number of excitations, to save computational
resources. Well-established examples of truncated CI methods are CIS (single
excitations), and CISD (singles and doubles excitations) [109]. An important
limitation of truncated CI methods is their lack of size-consistency, i.e. the
energy $E(A+B)$ of two infinitely-separated systems $A$ and $B$ is not equal
to $E(A)+E(B)$ [110]. Size consistency is generally regarded as being more
important than the provision of an energy upper-bound. Therefore, CISD results
are often modified using various corrections, such as the one introduced by
Langhoff and Davidson [111], to make the energies approximately size-
consistent, although this adjustment sacrifices variationality.

####  II.2.2 Selected CI

Selected CI methods rely on the same principle as the usual CI approaches.
However, in selected CI methods, determinants are not chosen solely based on
the number of electronic excitations. Instead, they are adaptively selected
from the entire set of determinants based on their estimated contribution to
the FCI wavefunction. This is because, even inside a predefined set of
determinants, only a fraction significantly contributes to the wave function
[112, 113].

The first multireference selected-CI algorithm going beyond singles and
doubles was the CIPSI (perturbatively selected configuration interaction
scheme) method of Huron et al [114, 115, 116], who proposed to iteratively
select external determinants (i.e. determinants which are not present in the
variational set) using a perturbative criterion. Recent years have witnessed
the development of various selected-CI approaches. Though based on the
original CIPSI method, these approaches feature significant variations in the
core idea as well as in the algorithm, in part motivated by the availability
of more powerful computational hardware. These include: (i) stochastic and
semi-stochastic approaches such as the FCI quantum Monte Carlo [117, 118, 119]
and the heat-bath CI approach [120, 121, 122]; (ii) purely variational
approaches such as the iterative and static-dynamic-static CI [123, 124];
(iii) size-consistent approximate FCI approaches such as full coupled-cluster
reduction (FCCR) [125] and the many-body expansion FCI (MBE-FCI) [126, 127,
128].

###  II.3 The equation-of-motion approach

The equation of motion (EOM) approach [129] is an alternative subspace-
construction technique, that remedies some of the limitations of truncated CI
methods (e.g. the lack of size-extensitivity and size-consistency of ground-
state energies and lack of size-intensivity of excitation energies). Starting
from the Schrödinger equation, Eq. (4), EOM introduces a set of excitation
operators $\hat{Q}^{\dagger}_{\mu}=|\Psi_{\mu}\rangle\langle\Psi_{0}|$ with
the properties: (i)
$\hat{Q}^{\dagger}_{\mu}|\Psi_{0}\rangle=|\Psi_{\mu}\rangle$, (ii)
$\hat{Q}_{\mu}|\Psi_{0}\rangle=0$, and (iii)
$\langle\Psi_{0}|\hat{Q}^{\dagger}_{\mu}=0$. Eq. (4) is then rewritten as

$[\hat{H},\hat{Q}^{\dagger}_{\mu}]|\Psi_{0}\rangle=\Delta
E_{\mu}\hat{Q}^{\dagger}_{\mu}|\Psi_{0}\rangle\;,$ (48)

where $\Delta E_{\mu}=E_{\mu}-E_{0}$ are excitation energies. Taking the
overlap with a generic state $\langle\Psi_{0}|\delta\hat{Q}$ leads to the
equation

$\langle\Psi_{0}|[\delta\hat{Q},[\hat{H},\hat{Q}^{\dagger}_{\mu}]]|\Psi_{0}\rangle=\Delta
E_{\mu}\langle\Psi_{0}|[\delta\hat{Q},\hat{Q}^{\dagger}_{\mu}]|\Psi_{0}\rangle\;,$
(49)

using the fact that $\langle\Psi_{0}|[\hat{H},\hat{Q}^{\dagger}_{\mu}]$ and
$\langle\Psi_{0}|\hat{Q}^{\dagger}_{\mu}$ are both zero. Eq. (49) is called an
EOM for the excitation operators $\hat{Q}^{\dagger}_{\mu}$ and energies
$\Delta E_{\mu}$. It can be turned into a GEEV, and then numerically solved,
by expanding the excitation operators as
$\hat{Q}^{\dagger}_{\mu}=\sum_{J}X^{\mu}_{J}\hat{F}_{J}-Y^{\mu}_{J}\hat{F}^{\dagger}_{J}$
where the $\hat{F}_{J}$ are some set of many-electron operators and
$X^{\mu}_{J},Y^{\mu}_{J}$ are coefficients. Then combining the equations
associated to each assignment of $\delta\hat{Q}=\hat{F}_{I}$ or
$\delta\hat{Q}=\hat{F}_{I}^{\dagger}$ for each $I$, and replacing the exact
unknown ground state $\Psi_{0}$ with an approximation $\Phi$, yields the
matrix equation

$\left(\begin{array}[]{cc}M&Q\\\ Q^{*}&M^{*}\\\
\end{array}\right)\left(\begin{array}[]{cc}X^{\mu}\\\ Y^{\mu}\\\
\end{array}\right)=\Delta\tilde{E}_{\mu}\left(\begin{array}[]{cc}V&W\\\
-W^{*}&-V^{*}\\\ \end{array}\right)\left(\begin{array}[]{cc}X^{\mu}\\\
Y^{\mu}\\\ \end{array}\right)$ (50)

with

$\begin{array}[]{ll}V_{IJ}=\phantom{-}\langle\Phi|[\hat{F}^{\dagger}_{I},\hat{F}_{J}^{\phantom{{\dagger}}}]|\Phi\rangle\;,&M_{IJ}=\phantom{-}\langle\Phi|[\hat{F}^{\dagger}_{I},[\hat{H},\hat{F}_{J}^{\phantom{{\dagger}}}]]|\Phi\rangle\;,\\\
\\\
W_{IJ}=-\langle\Phi|[\hat{F}^{\dagger}_{I},\hat{F}^{\dagger}_{J}]|\Phi\rangle\;,&Q_{IJ}=-\langle\Phi|[\hat{F}^{\dagger}_{I},[\hat{H},\hat{F}^{\dagger}_{J}]]|\Phi\rangle\;.\\\
\end{array}$ (51)

The Tamm-Dancoff (TDA) approximation to Eq. (50) is obtained by forcing
$Y^{\mu}=0$, and leads to the simplified eigenvalue equation
$MX^{\mu}=\Delta\tilde{E}_{\mu}VX^{\mu}$. The cost and accuracy of an EOM
calculation depends on two factors: the state $\Phi$, which may not be an
accurate approximation for the ground state in statically correlated systems,
and the truncation of electronic excitations connecting ground and excited
states (e.g. singles and doubles), which may not always be sufficient for
multireference excited states.

####  II.3.1 Time-dependent Hartree-Fock and equation-of-motion coupled-
cluster

Notable examples, summarised in Table 1, are the random phase approximation
(RPA) or time-dependent HF (TD-HF) method [130, 131, 132, 133, 134], the CIS
method (which is the TDA of TD-HF) [135, 136], and the EOM-CCSD method for
electronic excitations [137]. EOM-CCSD can also be formulated for electron-
removing excitations (IP) [138, 139, 140, 141], and electron-attaching
excitations (EA) [142]. EOM-CCSD is a widespread method to compute response
and excited-state properties, from energy differences like singlet-triplet
gaps, ionisation potentials, and electron affinities, to quasiparticle
properties, Green’s functions, the density of states, and spectral functions
[143, 144, 145]. Some of the advantages associated with the EOM-CCSD formalism
are its theoretical rigour, the accuracy and correct scaling behavior of
energy differences, and the ability to systematically improve results.
However, standard quantum chemistry methods such as EOM-CCSD sometimes face
challenges in a quantitative determination of excited states and their
properties, notably for same-symmetry conical intersections [146, 147, 148,
149] and when the ground state has a prominent multi-reference character [150,
151, 152, 153]. Since quantum algorithms are projected to provide accurate
ground-state wavefunctions, even in the case of statically-correlated systems,
they can improve these features of EOM-CCSD with practical computational
expenses.

name | $\Phi$ | $\hat{F}_{I}$ | TDA
---|---|---|---
RPA (TD-HF) | $\Phi_{\mathrm{HF}}$ | $\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{ai\sigma}$ | no
CIS | $\Phi_{\mathrm{HF}}$ | $\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{ai\sigma}$ | yes
IP-EOM-CCSD | $\Phi_{\mathrm{CCSD}}$ | $\\{\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{a\sigma}\cup\\{\hat{c}_{i\sigma}^{\phantom{\dagger}}\hat{c}_{a\tau}^{\dagger}\hat{c}_{j\tau}^{\phantom{\dagger}}\\}_{aij\sigma\tau}$ | no
EA-EOM-CCSD | $\Phi_{\mathrm{CCSD}}$ | $\\{\hat{c}_{a\sigma}^{\dagger}\\}_{a\sigma}\cup\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{i\tau}^{\phantom{\dagger}}\\}_{abi\sigma\tau}$ | no
EE-EOM-CCSD | $\Phi_{\mathrm{CCSD}}$ | $\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{ai\sigma}\cup\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{j\tau}^{\phantom{\dagger}}\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{abij\sigma\tau}$ | no
Table 1: Examples of EOM-based classical subspace methods.
$\Phi_{\mathrm{HF}}/\Phi_{\mathrm{CCSD}}$ denotes the Hartree-Fock/coupled-
cluster with singles and doubles ground-state, $ij$/$ab$ label
occupied/unoccupied spatial orbitals in the Hartree-Fock state, and
$\sigma\tau$ label spin polarisations.

###  II.4 Nonlinear subspace methods

Although they are somewhat less closely tied to existing quantum subspace
methods, it is worth mentioning the existence of subspace methods for
nonlinear eigenvalue problems [154, 155]. These methods find natural
applications in HF or other SCF (self-consistent field) calculations. The
high-level idea of these algorithms is to track a subspace that is iteratively
evolved along with the SCF updates. By doing so, the outer SCF iteration and
the inner subspace iteration are combined into a single loop, which can in
principle lead to substantial computational savings, although these methods
are quite modern so practical evidence is somewhat limited. For a recent
summary of these methods and their theoretical analysis, see [156].

##  III Quantum subspace methods

Figure 6: Schematics of quantum and classical subspace methods (QSM, CSM). The
search for Hamiltonian eigenpairs is restricted from a full Hilbert space
$\mathcal{H}$ to a subspace $V$ spanned by vectors
$\\{|{\bf{v}}_{\alpha}\rangle\\}_{\alpha=0}^{n-1}$ (left). Overlap and
Hamiltonian matrices are computed using a classical or quantum computer
(middle, enclosed in a rounded rectangle, with QSM/CSM abbreviating
quantum/classical subspace method). A generalised eigenvalue equation (GEEV)
is solved on a classical computer, yielding approximate Hamiltonian eigenpairs
(right).

The quantum subspace expansion (QSE) method [5, 6, 7, 8, 10, 11, 157, 158,
159, 160, 161, 162, 163, 164, 165] constructs a variational subspace spanned
by a set of states $\\{|{\bf{v}}_{\alpha}\rangle\\}_{\alpha=0}^{n-1}$, as
sketched in Fig. 6. Hamiltonian eigenstates are approximated by a linear
combination of basis states,
$|\tilde{\Psi}^{(n)}_{\mu}\rangle=\sum_{\alpha}C_{\alpha\mu}|{\bf{v}}_{\alpha}\rangle$.
To obtain the expansion coefficients $C_{\alpha\mu}$, one computes matrix
elements $S_{\alpha\beta}=\langle{\bf{v}}_{\alpha}|{\bf{v}}_{\beta}\rangle$
and $H_{\alpha\beta}=\langle{\bf{v}}_{\alpha}|\hat{H}|{\bf{v}}_{\beta}\rangle$
in the chosen basis and solves the GEEV $HC=SC\tilde{E}^{(n)}$. Within QSE, a
quantum device is used to compute the matrix elements $H_{\alpha\beta}$ and
$S_{\alpha\beta}$, and a classical computer to solve the GEEV and obtain
approximate eigenpairs
$\\{\tilde{E}^{(n)}_{\mu},|\tilde{\Psi}^{(n)}_{\mu}\rangle\\}_{\mu}$. The
quantum circuits required to compute $H_{\alpha\beta}$ and $S_{\alpha\beta}$
depend on the choice of the basis
$\\{|{\bf{v}}_{\alpha}\rangle\\}_{\alpha=0}^{n-1}$. Important families of QSMs
will be illustrated in Subsections III.1, III.3, and III.2, and alternative
approaches will be discussed in Subsection III.4.

The main advantage of QSE (and other subspace methods) over the well-
established variational quantum eigensolver (VQE) is that QSE does not require
a non-linear parameter optimisation, which is an NP-hard problem [166], as
part of the eigenvalue approximation. While both QSE and VQE employ quantum
and classical computers in synergy, the nature of their interaction is
profoundly different: in VQE, each update of the quantum circuit parameters
requires a new call of the quantum computer, while in QSE all the circuits
required to measure $H_{\alpha\beta}$ and $S_{\alpha\beta}$ can be sent to the
quantum computer in just one call. Furthermore, QSE circuits can be naturally
parallelised over multiple quantum computers (or groups of qubits within a
quantum computer). A limitation of QSE is that the eigenstates
$|\tilde{\Psi}^{(n)}_{\mu}\rangle$ are never actually stored on a quantum
device, so computing properties after a QSE calculation generally requires
additional measurements.

As mentioned above, the accuracy and computational cost of a QSE calculation
depends on the choice of the basis states
$\\{|{\bf{v}}_{\alpha}\rangle\\}_{\alpha=0}^{n-1}$. However, no specific
prescription is provided for this selection, similar to the original
Galerkin’s method. In this sense, QSE can be regarded as a family of QSMs,
differentiated by the choice of the basis states. We will now discuss
representative QSMs, starting from those based on $k$-body fermionic operators
(e.g., single and double electronic excitations), then moving to those based
on real- and imaginary-time evolution under the Hamiltonian, and polynomials
of the Hamiltonian.

###  III.1 QSMs based on k-body fermionic operators

The simplest example of a QSM is QSE based on $k$-body fermionic operators,
e.g. single- and double-excitation operators applied to an initial state
$|\Phi\rangle$,

$\begin{split}|\tilde{\Psi}^{(n)}_{\mu}\rangle&=\Big{[}C_{\mu}+\sum_{ai,\sigma}C^{ai,\sigma}_{\mu}\,\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}+\sum_{aibj,\sigma\tau}C_{\mu}^{aibj,\sigma\tau}\,\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{j\tau}^{\phantom{\dagger}}\hat{c}_{i\sigma}^{\phantom{\dagger}}\Big{]}|\Phi\rangle=\sum_{\alpha}C_{\alpha\mu}\,\hat{O}_{\alpha}|\Phi\rangle\;,\end{split}$
(52)

with
$\hat{O}_{\alpha}\in\\{\mathbbm{1}\\}\cup\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{ai\sigma}\cup\\{\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{j\tau}^{\phantom{\dagger}}\hat{c}_{i\sigma}^{\phantom{\dagger}}\\}_{abij\sigma\tau}$.
This flavor of QSE can be regarded to as a multi-reference CISD method
(MRCISD) where the wavefunction $|\Phi\rangle$, prepared on a quantum device,
is not a single Slater determinant but a correlated electronic state.

QSE based on single- and double-excitation operators has had significant
success in the last few years, see e.g. Section V for a list of
implementations on quantum hardware. As sketched in Fig. 7, it requires
measuring a set of operators on a register of qubits prepared in the state
$|\Phi\rangle$. Importantly, it does not increase the depth of the quantum
circuit required to prepare $|\Phi\rangle$ and measure the target operators.
This feature is especially beneficial on near-term hardware limited by qubit
coherence times and errors of two-qubit gates. Furthermore, it typically leads
to well-conditioned eigenvalue equations.

However, it is important to note that the overhead of measurements is
substantial in quantum MRCISD. This method requires estimating $k$-body
reduced density matrices (RDMs) of order higher than $k=2$, which poses a
significant computational bottleneck. We will discuss this issue further in
Section V.4.

Furthermore, a significant drawback of the quantum MRCISD approach, which it
shares with the classical CISD and MRCISD methods, is the lack of size-
intensivity in the computed excitation energies. This limitation can result in
quantitative errors and non-physical predictions. This may become a severe
limitation when QSE is applied to larger systems and/or the underlying ground-
state wavefunction is imprecise. We emphasise the existence of alternatives to
$k$-body fermionic operators in the construction of a subspace, for example,
Pauli operators [6] and elements of the Hamiltonian [167, 168].

Figure 7: Quantum circuits used in the multi-reference CISD (a) and self-
consistent quantum EOM (b) methods. Within MRCISD a register of qubits is
prepared in a state $|\Phi\rangle$, and operators defined in Eq. (52) are
measured (red meter blocks). Within quantum EOM, a different set of operators
is measured, with $\hat{F}_{I}$ defined as in Table 1. In the case of self-
consistent quantum EOM, a unitary operator $e^{\hat{X}}$ is applied to the
register before the measurement (teal block). The expectation values of the
operators in the top/bottom row define the overlap/Hamiltonian matrices of the
method, respectively.

####  III.1.1 Quantum equation-of-motion

In search of a size-intensive alternative to MRCISD, the quantum EOM method
(qEOM) was proposed by Ollitrault et al [9] for electronic excitation energies
(EEs). qEOM is essentially a transposition of the EOM-CC approach to quantum
computers. A quantum algorithm is used to produce an approximation
$|\Phi\rangle$ for the ground-state of the Hamiltonian, and the matrix
elements Eq. (51) are measured on the quantum computer over the state
$|\Phi\rangle$. While qEOM provides good agreement for EEs with the exact
results obtained by exact diagonalisation [9], it does not necessarily satisfy
the vacuum annihilation conditions $\hat{Q}^{\dagger}_{\mu}|\Phi\rangle=0$
ensuring that the ground-state wavefunction cannot be de-excited. This may
result in the appearance of large errors when the formalism is extended to
calculate properties such as IPs and EAs. Moreover, the qEOM method, like QSE
based on single- and double-excitation operators, requires high-body RDMs
which significantly increases the measurement cost. To remedy this limitation,
Asthana et al [169] proposed a self-consistent generalisation of qEOM, where
excitation operators have the form
$\hat{G}_{I}=e^{-\hat{X}}\hat{F}_{I}e^{\hat{X}}$ where $\hat{X}$ is a linear
combination of single- and double-excitations,

$\hat{X}=\sum_{ai,\sigma}x^{a}_{i}\,\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}+\sum_{abij,\sigma\tau}x^{ab}_{ij}\,\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{j\tau}^{\phantom{\dagger}}\hat{c}_{i\sigma}^{\phantom{\dagger}}-\mathrm{h.c.}\;,$
(53)

and $\hat{F}_{I}$ as in Table 1. The self-consistent qEOM method satisfies the
vacuum annihilation conditions, produces size-intensive and real-valued energy
differences between ground and excited/charges states, and requires the
measurement of 1- and 2-body RDMs. Notwithstanding these desirable features,
qEOM shares with quantum MRCISD a significant measurement overhead and, as
sketched in Fig. 7, deeper quantum circuits, due to the need of applying
$e^{\hat{X}}$ to $|\Phi\rangle$.

###  III.2 Hamiltonian Krylov-space methods

In the previous Subsection, we explored QSMs based on the application of
$k$-body fermionic operators to an initial state. These methods are compelling
for near-term devices as they do not require deep quantum circuits, but they
suffer from high measurement costs and accuracy limitations akin to those of
classical truncated CI calculations. Alternatively, there are QSMs based on
Krylov spaces, and a specific example is a Chebyshev quantum Krylov. However,
constructing a Krylov space generated by powers of the Hamiltonian on a
quantum computer is non-trivial. This is because the Hamiltonian and its
powers are non-unitary operations, and thus cannot be trivially mapped onto a
quantum circuit. The use of block-encoding unitaries was proposed to
circumvent this limitation.

####  III.2.1 Chebyshev quantum Krylov

The Chebyshev quantum Krylov method is a quantum implementation of the
classical Krylov space [165]. This method is exact in the sense that it has no
algorithmic error in the construction of the Krylov space, i.e.,
approximations such as Trotterization or other approximate time evolution
algorithms are not required. The method exploits the fact that the classical
Krylov space generated by powers of the Hamiltonian is mathematically
equivalent to the subspace generated by any basis for polynomials of the
Hamiltonian (powers may be viewed as the monomial basis).

The method is based on the notion of block-encoding [170] of a Hamiltonian
$\hat{H}$, i.e. a unitary operator $\hat{U}_{b}$ acting on an extended Hilbert
space and such that

$\langle
G,\cdot|\hat{U}_{b}|G,\cdot\rangle=\langle\cdot|\hat{H}|\cdot\rangle\;,$ (54)

for some state $|G\rangle$ (in the remainder of this paragraph we will omit
operator hats to avoid clutter). An example of block-encoding is the
following: for a Hamiltonian operator
$H=\sum_{\ell}h_{\ell}\,\sigma_{\boldsymbol{m}_{\ell}}$ where
$\sum_{\ell}|h_{\ell}|=\sum_{\ell}h_{\ell}=1$ and
$\sigma_{\boldsymbol{m}_{\ell}}$ are Pauli operators, the unitary
$U_{b}=\sum_{\ell}|\ell\rangle\langle\ell|\otimes\sigma_{\boldsymbol{m}_{\ell}}$,
in conjunction with the state
$|G\rangle=\sum_{\ell}\sqrt{h_{\ell}}|\ell\rangle=G|\boldsymbol{0}\rangle$,
defines a block encoding for $\hat{H}$. Under the assumption that the block
encoding is self-inverse, $U_{b}^{2}=\mathbbm{1}$, one has that

$\langle
G,\cdot|(RU_{b})^{\alpha}|G,\cdot\rangle=\langle\cdot|T_{\alpha}(H)|\cdot\rangle\;,$
(55)

where $R=\big{[}2|G\rangle\langle G|-\mathbbm{1}\big{]}\otimes\mathbbm{1}$ is
a reflection operator and $T_{\alpha}$ the $\alpha$-th Chebyshev polynomial of
the first kind. In other words, $(RU_{b})^{\alpha}$ is a block encoding of
$T_{\alpha}(H)$. Since the Chebyshev polynomials $T_{\alpha}$ with
$\alpha=0\dots n-1$ are a basis for polynomials of degree less than $n$, the
subspace spanned by the states
$|{\bf{v}}_{\alpha}\rangle=T_{\alpha}(H)|{\bf{v}}_{0}\rangle$, $k=0\dots n-1$,
coincides with the $n$-dimensional Krylov space generated by powers of the
Hamiltonian. The overlap and Hamiltonian matrix elements are given by

$S_{\alpha\beta}=\frac{f_{\alpha+\beta}+f_{|\alpha-\beta|}}{2}\;,\;H_{\alpha\beta}=\frac{f_{\alpha+\beta+1}+f_{|\alpha+\beta-1|}+f_{|\alpha-\beta|+1}+f_{|\alpha-\beta-1|}}{4}\;,$
(56)

where the function
$f_{\alpha}=\langle{\bf{v}}_{0}|T_{\alpha}(H)|{\bf{v}}_{0}\rangle$ can be
expressed as

$f_{\alpha}=\begin{cases}\langle{\bf{v}}_{0}|(U_{b}R)^{\alpha/2}R(RU_{b})^{\alpha/2}|{\bf{v}}_{0}\rangle&\alpha\;\mathrm{even}\\\
\langle{\bf{v}}_{0}|(U_{b}R)^{\lfloor\alpha/2\rfloor}U_{b}(RU_{b})^{\lfloor\alpha/2\rfloor}|{\bf{v}}_{0}\rangle&\alpha\;\mathrm{odd}\\\
\end{cases}$ (57)

and measured with the circuits in Fig. 8.

Figure 8: Quantum circuits used in the quantum Krylov method. For
$\alpha=0\dots n-1$ even/odd, we apply the left/right circuit. As defined in
the text, $U_{b}$ block-encodes the Hamiltonian $\hat{H}$, $G$ prepares the
state that identifies the block containing $\hat{H}$, and $R$ reflects around
that state. In panel (b), when the ancillae collapse in the state
$|\ell\rangle$ upon measurement, the Pauli operator
$\sigma_{\boldsymbol{m}_{\ell}}$ is measured on the main qubit register.

This quantum Krylov method is a technique to construct a Krylov space on a
quantum computer. It addresses the challenge posed by the non-unitarity of
Hamiltonian powers using a block-encoding unitary transformation. Due to its
relationship with the classical Lanczos method, it is provably convergent to
the ground state (see Section II.1) and can benefit from heuristic quantum
computing algorithms providing initial states having high overlap with the
ground state. Furthermore, it offers an important advantage in terms of the
memory required to store the Krylov space information, i.e. it bypasses the
need to store combinatorially large CI strings in memory. On the other hand,
like the classical Lanczos method, it may be numerically ill-conditioned (see
Section II.1) and sensitive to shot noise (see Section V.2 for a detailed
discussion). Furthermore, the use of a block-encoding unitary requires
ancillae and deep quantum circuits, making it more suited for future quantum
devices [165].

####  III.2.2 Gaussian-power quantum Krylov

This method constructs a subspace using the conventional power function acting
on the Hamiltonian, as described in Section II.1. However, the reference state
$|{\bf{v}}_{0}\rangle$ is now replaced with the state
$e^{-\frac{1}{2}(H-E_{0})^{2}\tau^{2}}|{\bf{v}}_{0}\rangle$. In other words,
each of the basis states $|{\bf{v}}_{\alpha}\rangle$ is generated as follows:

$|{\bf{v}}_{\alpha}\rangle=(\hat{H}-E_{0})^{\alpha}e^{-\frac{1}{2}(H-E_{0})^{2}\tau^{2}}|{\bf{v}}_{0}\rangle,\;\alpha=0\dots
n-1$ (58)

where $E_{0}$ is an arbitrary constant of choice. This is known as the
Gaussian-power basis [171]. The spectral norm for each basis vector
$|{\bf{v}}_{\alpha}\rangle$ can be bounded as
$|||{\bf{v}}_{\alpha}\rangle||_{2}\leq\big{(}\dfrac{n-1}{e\tau^{2}}\big{)}^{\frac{n-1}{2}}$.
Therefore, when $e\tau^{2}\geq n-1$, the spectral norm of
$|{\bf{v}}_{\alpha}\rangle$ decreases exponentially with respect to $\alpha$,
and the rate is parameterized by the parameter $\tau$. This is useful since
the exact projected matrices $H$ and $S$ are unknown due to the statistical
errors that arise from finite sampling. The exponential suppression of the
spectral norm of $|{\bf{v}}_{\alpha}\rangle$ is the main factor leading to a
significant decrease in statistical errors and minimizes the number of
measurements required in this approach.

One way to realise Eq. (58) on a quantum computer is through the method of
Linear Combination of Unitaries (LCU), as proposed by [171]. In this approach,
the author showed that the generating function
$(\hat{H}-E_{0})^{\alpha}e^{-\frac{1}{2}(\hat{H}-E_{0})^{2}\tau^{2}}$ in Eq.
(58) can be rewritten as

$(\hat{H}-E_{0})^{\alpha}e^{-\frac{1}{2}(\hat{H}-E_{0})^{2}\tau^{2}}=\left(\frac{i}{\sqrt{2\tau^{2}}}\right)^{\alpha}\int_{-\infty}^{\infty}dt\,H_{\alpha}\Big{(}\dfrac{t}{\sqrt{2\tau^{2}}}\Big{)}\,g_{\tau}(t)\,e^{-it(H-E_{0})}\;,$
(59)

where $H_{\alpha}(u)$ are Hermite polynomials, and
$g_{\tau}(t)=e^{-t^{2}/2\tau^{2}}/\sqrt{2\pi\tau^{2}}$. In this
representation, the only function that depends on the Hamiltonian is the real-
time evolution operator. Several techniques, e.g. Trotterisation and
particularly those discussed in Section III.3.1, can be employed for this
purpose. Another alternative realization of Eq. (58) on a quantum computer is
through the technique of block encoding.

###  III.3 QSMs based on time evolution

Here we present a third family of QSMs, that make use of real- or imaginary-
time evolution operations to construct subspaces. These methods are the
Quantum Filter Diagonalisation (QFD) and Quantum Lanczos (QLanczos).

####  III.3.1 Quantum Filter Diagonalisation

This method constructs a subspace by applying the time evolution operator to
an initial state [10],

$|{\bf{v}}_{\alpha}\rangle=e^{-it_{\alpha}\hat{H}}|{\bf{v}}_{0}\rangle\;,\;\alpha=0\dots
n-1$ (60)

where the times $t_{\alpha}$ are a set of times often, though not necessarily,
given by $t_{\alpha}=\alpha\Delta t$ for some time step $\Delta t>0$. Note
that technically the subspace is a Krylov space only if
$t_{\alpha}=\alpha\Delta t$, since then the operators generating the subspace
are powers of the unitary $e^{-i\Delta t\hat{H}}$. In the remainder of the
Section, we will assume this construction is used. The same basic idea is
known elsewhere in the literature as quantum subspace diagonalization [172]
and variational quantum phase estimation [162]; herein we will use QFD since
that was the name conferred in the original paper proposing this method [10].
A natural and compelling generalisation of QFD is the multireference selected
quantum Krylov (MRSQK) algorithm [11], where time evolution is applied to a
set of initial states,

$|{\bf{v}}_{\alpha
k}\rangle=e^{-it_{\alpha}\hat{H}}|{\bf{v}}_{0k}\rangle\;,\;\alpha=0\dots
n-1\;,\;k=0\dots n_{k}-1\;.$ (61)

QFD is a particularly compelling method because it possesses a formal error
analysis [172], even accounting for noise on the quantum device.

Before discussing that analysis, we mention a common pitfall in thinking about
the error of QFD. For small $\Delta t$, the QFD Krylov space is approximately
equal to the classical Krylov space (which we use to refer to the Krylov space
generated by powers of the Hamiltonian) up to high-order terms,

$\sum_{\beta}c_{\beta}|{\bf{v}}_{\beta}\rangle=\sum_{\alpha=0}^{n-1}\left[\sum_{\beta=0}^{n-1}M_{\alpha\beta}c_{\beta}\right]\hat{H}^{\alpha}|{\bf{v}}_{0}\rangle+\mathcal{O}(\Delta
t^{n})\;,\;M_{\alpha\beta}=\frac{(-i\beta\Delta t)^{\alpha}}{\alpha!}\;,$ (62)

as the matrix $M$ is invertible (it is the product of a Vandermonde and a
diagonal matrix). Since the residual terms in $\mathcal{O}(\Delta t^{n})$
vanish as $\Delta t$ is taken to zero, one might be tempted to run the
algorithm for some very small $\Delta t$ since the errors resulting from
classical Krylov methods provably converge [90, 91, 92], see Section II.1.
However, for very small $\Delta t$, the QFD subspace approaches linear
dependence (time-evolution operators approximate the identity operator
closely). Such a scenario induces ill-conditioning and thus complicated
tradeoffs in the choice of $\Delta t$.

However, those tradeoffs turn out to be illusory, essentially because the
residual term in (62) does not actually reduce the quality of the lowest-
energy state in the QFD subspace, so $\Delta t$ does not have to be pushed as
close to zero as possible. This fact was demonstrated in Epperly et al [172]
(in particular Theorem 3.1). Choosing $\Delta t=\frac{\pi}{\Delta E_{L}}$ and
$t_{\alpha}=\frac{\pi\alpha}{\Delta E_{L}}$, $\alpha=-n/2\dots n/2$, where
$\Delta E_{L}$ is the energy difference between the ground state and some
excited state with index $0<L<D-1$, one has

$0\leq\tilde{E}^{(n)}_{0}-E_{0}\leq\left(1+\frac{\pi\Delta E_{1}}{\Delta
E_{L}}\right)^{-n}8\sum_{\mu=1}^{L}\Delta
E_{\mu}\frac{w_{\mu}}{w_{0}}+2\sum_{\mu=L+1}^{D-1}\Delta
E_{\mu}\frac{w_{\mu}}{w_{0}}\;,$ (63)

$w_{\mu}=|\langle\Psi_{\mu}|{\bf{v}}_{0}\rangle|^{2}$. The error
$\tilde{E}^{(n)}_{0}-E_{0}$ contains two terms, arising from the components of
$|{\bf{v}}_{0}\rangle$ on eigenstates with energies respectively below and
above $E_{L}$. The first term decreases exponentially with subspace dimension
$n$, albeit at a rate determined by the ratio $\frac{\Delta E_{1}}{\Delta
E_{L}}$, meaning that larger QFD subspaces remove unwanted low-energy
components of $|{\bf{v}}_{0}\rangle$ as long as the system has a non-zero gap.
The second is independent of $n$, meaning that high-energy components of
$|{\bf{v}}_{0}\rangle$ cannot be eliminated by increasing the dimension of the
QFD subspace, but only by choosing a finer mesh of times $t_{\alpha}$
(equivalent to increasing $L$) or a trial state with support on the low-energy
subspace ($w_{\mu}\simeq 0$ for $\mu>L$).

By choosing $L=D-1$ we get an upper bound on the ground state energy
approximation that is free from conditions on the support of the trial state,
obtaining

$0\leq\tilde{E}^{(n)}_{0}-E_{0}\leq\left(1+\frac{\pi\Delta E_{1}}{\Delta
E_{D-1}}\right)^{-n}8\sum_{\mu=1}^{D-1}\Delta E_{\mu}\frac{w_{\mu}}{w_{0}}.$
(64)

This shows that the smallest $\Delta t$ should ever be chosen is ${\Delta
t=\frac{\pi}{\Delta E_{D-1}}\geq\frac{\pi}{\|H\|}}$, where $\|\cdot\|$ is the
spectral norm (but could in practice be replaced by more easily computed upper
bounds). The convergence of the method is thus not improved by making $\Delta
t$ smaller past this point, even if ill-conditioning were not a concern. We
note that, in the presence of noise and imperfect time-evolution simulation,
the above choice of $\Delta t$ may no longer be optimal, and it should instead
be chosen heuristically. Roughly, the relation of the analysis of Epperly et
al [172] to the classical analysis of Kaniel, Paige, and Saad [90, 91, 92] is
that while the latter is based on finding approximate projectors among
polynomials of the Hamiltonian (which are elements of the classical Krylov
space), the former is based on finding approximate projectors among
trigonometric polynomials of the Hamiltonian (which are elements of the QFD
space), without any need to approximate to the classical Krylov space along
the way. This analysis is discussed further in Subsection V.2.

Figure 9: Quantum circuits used in the QFD method to compute
$O_{\alpha\beta}=\langle{\bf{v}}_{\alpha}|\hat{O}|{\bf{v}}_{\beta}\rangle$. An
ancilla is prepared in the state $|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$,
and the unitary transformation $|0\rangle\langle
0|\otimes\hat{U}_{\beta}+|1\rangle\langle 1|\otimes\hat{U}_{\alpha}$ is
applied, where $\hat{U}_{\alpha}=e^{-it_{\alpha}\hat{H}}$ and
$\hat{U}_{\beta}=e^{-it_{\beta}\hat{H}}$ are time-evolution unitaries (or
approximations thereof) controlled by the state of the ancilla. The operator
$2|1\rangle\langle 0|\otimes\hat{O}$ is then measured, yielding an estimate of
$O_{\alpha\beta}$. For $\hat{O}=\mathbbm{1}/\hat{H}$ one obtains
overlap/Hamiltonian matrix elements respectively.

From a quantum computing perspective, QFD is appealing because it is based on
real-time evolution, an operation that can be simulated on a quantum computer
with polynomial cost and controllable accuracy (more formally, it lies in the
complexity class BQP). The method lies between the variational quantum
eigensolver (VQE) and the phase estimation algorithm (PEA) in terms of
required quantum circuit resources and conceptual simplicity. Compared to the
Chebyshev quantum Krylov method discussed in Sec. III.2.1, QFD is appealing
because time-evolutions admit relatively near-term-friendly approximate
implementations.

On the other hand, QFD leads to an eigenvalue equation that is in general ill-
conditioned in practice, much like in classical Krylov methods (see Subsection
II.1), although the Beckermann-Townsend inequality [93] does not apply since
in QFD the Krylov space is not generated by powers of a Hermitian matrix. This
ill-conditioning is aggravated by the presence of decoherence and shot noise,
although it can be ameliorated by a thresholding procedure [172], as discussed
in Subsection V.2. Furthermore, time evolution is an expensive operation for
the electronic structure Hamiltonian, as noted in Subsection I.2.6, and the
implementation of QFD in general requires controlled application of a time-
evolution circuit, as sketched in Fig. 9.

In addition to the theoretical analysis in [172], the accuracy of QFD and the
dependence of QFD energies on time step, subspace dimension, and shot noise,
were extensively analyzed numerically by Klymko et aland Shen et al [162,
173]. These papers also introduced and analyzed the so-called unitary version
of QFD, where the GEEV can be written in a unitary form, requiring only
measurement of overlap matrix elements, and restoring the Toeplitz structure
of the GEEV in the case of Trotterized time evolution with a uniform grid.
Additionally, QFD’s asymptotic computational cost was reduced by means of low-
rank approximations by Cohn et al [158] and stochastic compilation of the time
evolution subroutine was proposed by Stair et al [174]. In Subsection V.4 we
examine the issue of controlled time evolution, and describe procedures to
bypass this expensive operation under suitable assumptions.

####  III.3.2 Quantum Lanczos based on imaginary-time evolution

The quantum Lanczos (QLanczos) method constructs a subspace of vectors by
applying an imaginary-time evolution (ITE) operator to an initial state,

$|{\bf{v}}_{\alpha}\rangle=\frac{e^{-\alpha\Delta\tau\hat{H}}|{\bf{v}}_{0}\rangle}{\|e^{-\alpha\Delta\tau\hat{H}}|{\bf{v}}_{0}\rangle\|}\;,\;\alpha=0\dots
n-1\;,$ (65)

where $\Delta\tau>0$ is an imaginary-time step. Like QFD, QLanczos is formally
connected with the classical Krylov space in the sense of Eq. (62). However,
an important difference between QFD and QLanczos is that the latter is based
on an operation, the ITE, that is neither unitary nor linear. On a fault-
tolerant quantum device, non-linear operations can be implemented using
techniques like quantum signal processing [175, 176, 177, 178]. On a near-term
device, ITE can be approximated by variational ansatz-based simulations [179],
or by the quantum ITE (QITE) algorithm [8]. The latter approximates a step of
ITE on a wavefunction $\Phi$ by the action of a quantum circuit
$\hat{U}(\theta)$,

$\frac{e^{-\Delta\tau\hat{H}}|\Phi\rangle}{\|e^{-\Delta\tau\hat{H}}\Phi\|}\simeq\hat{U}(\boldsymbol{\theta})|\Phi\rangle\;.$
(66)

Both members of the equation above are expanded to the first order around the
identity,

$\Big{[}\mathbbm{1}-\Delta\tau\big{(}\hat{H}-E_{\Phi}\big{)}\Big{]}|\Phi\rangle\simeq\Big{[}\mathbbm{1}+i\sum_{m}\theta_{m}\hat{G}_{m}\Big{]}|\Phi\rangle\;,$
(67)

where $E_{\Phi}=\langle\Phi|\hat{H}|\Phi\rangle$. The parameters $\theta_{m}$
are determined solving the equation $A\boldsymbol{\theta}=\boldsymbol{b}$,
where $A_{mn}=\langle\Phi|\\{\hat{G}_{m},\hat{G}_{n}\\}|\Phi\rangle$ and
$b_{m}=-i\Delta\tau\langle\Phi|[\hat{H},\hat{G}_{m}]|\Phi\rangle$. For a
$k$-local spin Hamiltonian [8], a natural choice is
$\hat{U}(\boldsymbol{\theta})=\prod_{m}\exp(ix_{m}\sigma_{m})$, where
$\sigma_{m}$ are spatially local Pauli operators. For electronic systems
[180], a compelling alternative is the use of one- and two-body fermionic
operators,

$\hat{U}(\boldsymbol{\theta})=\exp\left(\sum_{ai}\theta^{a}_{i}\hat{c}_{a\sigma}^{\dagger}\hat{c}_{i\sigma}^{\phantom{\dagger}}+\sum_{aibj}\theta^{ab}_{ij}\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{j\tau}^{\phantom{\dagger}}\hat{c}_{i\sigma}^{\phantom{\dagger}}-\mathrm{h.c.}\right)\;.$
(68)

A limitation of QITE is the growth of circuit depth with number of imaginary-
time steps, which can be ameliorated e.g. using the step-merged approach of
Gomes et al [180], where the approximation
$\hat{U}(\boldsymbol{\theta}_{1})\hat{U}(\boldsymbol{\theta}_{2})\simeq\hat{U}(\boldsymbol{\theta}_{1}+\boldsymbol{\theta}_{2})$
is made, in order to implement QITE with a fixed-depth circuit, albeit with a
possible loss of accuracy.

On the other hand, QLanczos does not require ancillae and controlled
operations, because matrix elements can be written in terms of norms and
expectation values,

$\langle{\bf{v}}_{\alpha}|{\bf{v}}_{\beta}\rangle=\frac{n^{2}_{\frac{\alpha+\beta}{2}}}{n_{\alpha}n_{\beta}}\;,\;\langle{\bf{v}}_{\alpha}|\hat{H}|{\bf{v}}_{\beta}\rangle=\frac{n^{2}_{\frac{\alpha+\beta}{2}}}{n_{\alpha}n_{\beta}}h_{\frac{\alpha+\beta}{2}}\;,\;$
(69)

with $n_{\alpha}=\|e^{-\alpha\Delta\tau\hat{H}}|{\bf{v}}_{0}\rangle\|$ and
$h_{\alpha}=\langle{\bf{v}}_{\alpha}|\hat{H}|{\bf{v}}_{\alpha}\rangle$, and
these quantities can be evaluated without ancillae and controlled operations
as the ITE unfolds [8].

###  III.4 Alternative approaches

In the previous Subsections, we discussed representative examples of QSMs.
Research in the design, implementation, and refinement of QSMs extends beyond
these classes of methods. In this Subsection, we describe alternative QSMs.

####  III.4.1 Subspace-based variational quantum simulations

Electronic ground- and excited-states can be approximated by variational
quantum simulations, exemplified by the variational quantum eigensolver (VQE)
[181], wherein the ground-state wavefunction, $|\Psi_{0}\rangle$, and energy,
$E_{0}$, are approximated by variationally optimising a parameterised
wavefunction ansatz $|\Phi(\boldsymbol{\theta})\rangle$,

$E_{\mathrm{VQE}}=\min_{\boldsymbol{\theta}}\,\langle\Phi(\boldsymbol{\theta})|H|\Phi(\boldsymbol{\theta})\rangle=\min_{\boldsymbol{\theta}}\,E(\boldsymbol{\theta})\;.$
(70)

The energy $E(\boldsymbol{\theta})$ is evaluated on a quantum computer, and
parameters $\boldsymbol{\theta}$ are optimised on a classical computer.
However, the quality of a VQE calculation depends on the ansatz and the
convergence of the optimisation procedure. Literature [9, 182, 183, 184, 185,
186] indicates that VQE applied to small active spaces can yield energies
close to those of CASCI, but also that it can require a computational cost
prohibitive for near-term devices [187] and feature symmetry-breaking, non-
differentiable potential-energy curves, and exponentially expensive parameter
optimisation when hardware-efficient ansatzes are used [188]. QSMs can be
constructed starting from a wavefunction generated by a VQE simulation, making
them a natural and compelling approach to enhance the quality of variational
simulations.

An example of this improvement is the multistate-contracted variational
quantum eigensolver (MC-VQE) [189], where a parametrised Ansatz
$\hat{U}(\boldsymbol{\theta})$ is applied to a linear combination of states
$|{\bf{e}}_{\alpha}\rangle$ that are e.g. qubit representations of electronic
configurations,

$|\Phi(\boldsymbol{\theta},{\bf{c}})\rangle=\hat{U}(\boldsymbol{\theta})\sum_{\alpha}c_{\alpha}|{\bf{e}}_{\alpha}\rangle\;.$
(71)

The parameters $\boldsymbol{\theta},{\bf{c}}$ are jointly optimised. In
particular, optimisation of coefficients ${\bf{c}}$ for a given
$\boldsymbol{\theta}$ is a ground-state search in the subspace spanned by the
states $\hat{U}(\boldsymbol{\theta})|{\bf{e}}_{\alpha}\rangle$. In a similar
vein, the non-orthogonal VQE (NO-VQE) method [7] proposes to construct the
variational ansatz

$|\Phi(\boldsymbol{\theta},{\bf{c}})\rangle=\sum_{\alpha}c_{\alpha}|\Phi_{\alpha}(\boldsymbol{\theta}_{\alpha})\rangle\;,\;|\Phi_{\alpha}(\boldsymbol{\theta}_{\alpha})\rangle=\hat{U}_{\alpha}(\boldsymbol{\theta}_{\alpha})|{\bf{v}}_{0}\rangle\;,$
(72)

where the parameters $\boldsymbol{\theta}_{\alpha}$ are optimised in an outer
loop and the coefficients $c_{\alpha}$ are determined by a ground-state search
in the subspace spanned by the states
$|\Phi_{\alpha}(\boldsymbol{\theta}_{\alpha})\rangle$. NO-VQE can be
considered a generalisation of MC-VQE; on the other hand, it requires a
Hadamard test (circuit in Fig. 10a), whereas the quantum circuits required by
MC-VQE are of the form in Fig. 10b (no Hadamard test). The non-orthogonal
quantum eigensolver (NOQE) is a modification [163] of NO-VQE that relies on a
quantum computer to synthesise wavefunctions of the form
$|\Phi_{\alpha}(\boldsymbol{\theta}_{\alpha})\rangle$, where parameters are
suggested by classical perturbative calculations on top of spin-unrestricted
Slater determinants, and to compute overlap and Hamiltonian matrices to solve
for Hamiltonian eigenstates in the subspace spanned by the states
$|\Phi_{\alpha}(\boldsymbol{\theta}_{\alpha})\rangle$, without further
parameter optimisation. NOQE can be considered a quantum transposition of the
classical non-orthogonal configuration interaction (NOCI) method [190, 191].

Figure 10: Quantum circuits used in the (a) NO-VQE and (b) MC-VQE methods. NO-
VQE requires an ancilla and controlled unitaries, whereas MC-VQE requires the
preparation of superposition states of the form
$|{\bf{e}}^{p}_{\alpha\beta}\rangle=(|{\bf{e}}_{\alpha}\rangle+i^{p}|{\bf{e}}_{\beta}\rangle)/\sqrt{2}$
with $p=0,1,2,3$. Measuring the operator $\hat{O}=\mathbbm{1}/\hat{H}$ in NO-
VQE yields overlap/Hamiltonian matrix elements respectively. Within MC-VQE,
the overlap matrix is equal to the identity due to the orthonormality of the
states $|{\bf{e}}_{\alpha}\rangle$, so that only Hamiltonian measurements are
required.

####  III.4.2 Quantum selected CI

Quantum computing algorithms can be used to inform classical subspace
calculations. An interesting example is the quantum selected CI method by
Kanno et al [192], where an approximation $\tilde{\Phi}$ for the ground state
of a quantum system is prepared on a quantum computer, for example using a
variational method. Measuring a register of qubits prepared in
$|\Phi\rangle=\sum_{{\bf{x}}}c_{{\bf{x}}}|{{\bf{x}}}\rangle$ yields a
bitstring ${{\bf{x}}}$ with probability ${p_{{\bf{x}}}=|c_{{\bf{x}}}|^{2}}$.
Executing the measurement $N$ times yields bitstrings
${{{\bf{x}}}_{m},m=1\dots N}$. Having sampled those bitstrings, one can
construct a matrix $H_{mn}=\langle{{\bf{x}}}_{m}|H|{{\bf{x}}}_{n}\rangle$ on a
classical computer using the Slater-Condon rules, and diagonalise it on a
classical computer, to find another approximation
$|\tilde{\Phi}\rangle=\sum_{m=1}^{N}\tilde{c}_{{\bf{x}}_{m}}|{{\bf{x}}}_{m}\rangle$
for the ground state. This algorithm is variational, and intrinsically
resistant to decoherence. While the probability distribution $p_{{\bf{x}}}$ of
an actual simulation may differ from $|c_{{\bf{x}}}|^{2}$ due to decoherence,
$|\tilde{\Phi}\rangle$ is prepared on a classical computer. This algorithm
presents several challenges and opportunities for refinement. First, ensuring
the generation of continuous potential energy curves: when independent finite
samples are drawn from the probability distribution $p_{{\bf{x}}}$ along a
potential energy curve, discrepancies in the corresponding bitstrings can
result in irregular potential energy curves. Second, addressing sampling
inefficiency: in cases involving dynamically correlated wavefunctions,
$p_{{\bf{x}}}$ tends to be concentrated around specific bitstrings, leading to
repetitive sampling and reduced efficiency. Third, developing cost-effective
and systematic approaches to modify the probability distribution
$p_{{\bf{x}}}$ in such a way as to prepare low-energy states
$|\tilde{\Phi}\rangle$.

####  III.4.3 Eigenvalue continuation

Eigenvector Continuation (EC) as a QSM method that considers a family of
Hamiltonians, $\hat{H}_{\boldsymbol{\lambda}}$, that depend on one or more
parameters ${\boldsymbol{\lambda}}$ [193, 194]. This situation is encountered
in ES when ${\boldsymbol{\lambda}}$ is e.g. a set of nuclear coordinates or an
external field applied to a molecule. A set of approximations
$\Phi_{{\boldsymbol{\lambda}}_{\alpha}}$ for low-energy states of the
Hamiltonian at different points ${{\boldsymbol{\lambda}}_{\alpha}}$ of the
parameter space are chosen as the subspace basis. In other words, the ground
state of $\hat{H}_{\boldsymbol{\lambda}}$ is approximated as a linear
combinaton of the form
$\sum_{\alpha}c_{\alpha}|\Phi_{{\boldsymbol{\lambda}}_{\alpha}}\rangle$, and
the coefficients $c_{\alpha}$ are determined using the standard QSE method. EC
allows to construct a subspace using physics-informed basis states, thereby
allowing accurate calculations of low-energy eigenvalues. In general, EC
requires a Hadamard test to measure overlap and Hamiltonan matrices (see
Section V.4 for further discussion). For accurate results, the Hamiltonian and
overlap matrices need to be measured accurately, and the overlap matrix needs
to be well-conditioned thanks to a suitable choice of the
${\boldsymbol{\lambda}}$ parameters, or regularized via e.g. a threshold
procedure (see Section V.2 for a discussion).

##  IV Applications of quantum subspace methods

Some of the most natural and widespread applications of QSMs are the
calculations of ground and low-lying excited electronic states. QSMs can also
be used to study various phenomena, assuming that these can be accurately
approximated using electronic eigenfunctions from the chosen subspace. In this
Section, we illustrate the description of dynamical electron correlation, the
fast-forwarding of time evolution, and the computation of frequency-dependent
response functions.

###  IV.1 Dynamical electronic correlation

Near-term quantum devices are predominantly used to carry out active-space
simulations in the context of electronic structure (see Section I.1.2 and Fig.
2). While active-space simulations can describe forms of static electronic
correlation, they cannot recover the dynamical correlation (described in
Section I.1.2). A similar situation is encountered in classical ES methods
when an active-space calculation is conducted, e.g. with CASCI or CASSCF.
Retrieval of dynamical correlation is then carried out by combining the
active-space method with a low-level calculation on the full basis set, e.g.
second-order perturbation theory. While a straightforward possibility is to
simulate a realistic basis set with larger quantum simulations, this strategy
is not always compatible with near-term devices [187].

The possibility of extending QSMs with post-processing on classical computers
to account for dynamical electronic correlation has been recognised and is an
active research area. In this section, we will present illustrative studies.
First, starting from the Hamiltonian Eq. (1), we divide orbitals into core
(indices $i,j,k\in C$), active (indices $t,u,v,w\in A$), and virtual (indices
$a,b,c\in V$). For simplicity, we consider core orbitals to be frozen, i.e.,
they are doubly occupied and electrons in them are never excited into other
orbitals. Since the Hamiltonian with frozen-core orbitals can be transformed
into one containing only active and virtual orbitals, we will now exclude the
core space. Active orbitals are considered crucial because electrons in these
orbitals exhibit static correlation, making them the focus of treatment on
quantum computers. Virtual orbitals contribute to additional dynamical
correlation.

A hybrid quantum-classical technique for dynamical correlation is virtual QSE
[157]. This method starts from an approximation for the ground state of the
system as $|\Phi,\mbox{\O}\rangle$, where $|\Phi\rangle$ is an active-space
wavefunction and $|\mbox{\O}\rangle$ the vacuum state for virtual orbitals.
Then, it defines a subspace spanned by states of the form
$\hat{O}_{\alpha}|\Phi,\mbox{\O}\rangle$, where $\hat{O}_{\alpha}$ belongs to
the set

$D=\Big{\\{}\hat{c}_{a\sigma}^{\dagger}\hat{c}_{t\sigma}^{\phantom{\dagger}},\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\tau}^{\dagger}\hat{c}_{u\tau}^{\phantom{\dagger}}\hat{c}_{t\sigma}^{\phantom{\dagger}}\;,\;uv\in
A,ab\in V,\tau\sigma\in\\{\uparrow,\downarrow\\}\Big{\\}}\;.$ (73)

Overlap and Hamiltonian matrices are formed as in standard QSE. Evaluating
these matrices requires tracing out virtual degrees of freedom using Wick’s
theorem [195], e.g.

$S_{bu\tau,at\sigma}=\langle\Phi,\mbox{\O}|\hat{c}_{u\tau}^{\dagger}\hat{c}_{b\tau}^{\phantom{\dagger}}\hat{c}_{a\sigma}^{\dagger}\hat{c}_{t\sigma}^{\phantom{\dagger}}|\Phi,\mbox{\O}\rangle=\delta_{ab}\delta_{\sigma\tau}\langle\Phi|\hat{c}_{u\sigma}^{\dagger}\hat{c}_{t\sigma}^{\phantom{\dagger}}|\Phi\rangle\;,$
(74)

where
$\langle\Phi|\hat{c}_{u\sigma}^{\dagger}\hat{c}_{t\sigma}^{\phantom{\dagger}}|\Phi\rangle$
is an element of the active-space one-body density matrix. Detailed formulas
and representative applications are given in Refs. [196, 157]. Virtual QSE
allows accounting for dynamical correlation by allowing single and double
excitations in the virtual space. However, the quality of the resulting
energies is comparable to that of an MRCISD calculation, and it requires
computing high-order RDMs.

An alternative to virtual QSE is second-order perturbation theory [197]. This
formalism is a quantum transposition of classical $N$-electron valence
perturbation theory (NEVPT2) [198, 199, 200, 201]. The Hamiltonian is
partitioned as $\hat{H}=\hat{H}_{D}+\hat{V}$, where $\hat{H}_{D}$ is the Dyall
Hamiltonian [202], i.e. the sum between the active space Born-Oppenheimer
Hamiltonian and the restriction of the Fock operator to the external space,
and $\hat{V}=\hat{H}-\hat{H}_{D}$ is treated as a perturbation. The second-
order energy contribution can be written as

$-\Delta E_{\mathrm{PT2}}=\sum_{\nu\neq
0}\frac{|\langle\Psi_{\nu}|\hat{V}|\Psi_{0}\rangle|^{2}}{E_{\nu}-E_{0}}\;,$
(75)

where $(\Psi_{\nu},E_{\nu})$ are the eigenpairs of the Dyall Hamiltonian, and
$\nu=0$ labels the ground state. Eq. (75) yields the exact (or uncontracted)
NEVPT2. Implementing uncontracted NEVPT2 has a combinatorial cost with active-
space size, due to the summation over excited states. This limitation can be
remedied using strongly-contracted NEVPT2 [203], which requires high-order
ground-state RDMs, or partially-contracted NEVPT2 [197], which approximates
the sum over excited states. The latter is very naturally interfaced with
QSMs. One can observe that the action of the perturbation over the ground
state,

$\hat{V}|\Phi\rangle=\Bigg{[}\sum_{\begin{subarray}{c}a\\\
\sigma\end{subarray}}\hat{c}_{a\sigma}^{\dagger}\hat{O}^{(1)}_{a,\sigma}+\sum_{\begin{subarray}{c}a<b\\\
\sigma\end{subarray}}\hat{c}_{a\sigma}^{\dagger}\hat{c}_{b\sigma}^{\dagger}\hat{O}^{(2)}_{ab,\sigma}+\sum_{ab}\hat{c}_{a\uparrow}^{\dagger}\hat{c}_{b\downarrow}^{\dagger}\hat{O}^{(3)}_{ab}\Bigg{]}|\Psi_{0}\rangle\;,$
(76)

removes a particle with spin $\sigma$, or two particles with identical spins
$\sigma$, or two particles with opposite spin from the active space, exciting
the remaining electrons through the action of suitable operators
$\hat{O}^{(1)}_{a,\sigma}$, $\hat{O}^{(2)}_{ab,\sigma}$, $\hat{O}^{(3)}_{ab}$.
Therefore, one can restrict the summation in Eq. (75) to excited states with
$(N_{\uparrow}-\Delta N_{\uparrow},N_{\downarrow}-\Delta N_{\downarrow})$
particles, where $(N_{\uparrow},N_{\downarrow})$ is the number of electrons in
the active-space ground-state wavefunction and $\Delta N_{\uparrow},\Delta
N_{\downarrow}=0,1,2$ and $\Delta N_{\uparrow}+\Delta N_{\downarrow}=1,2$. The
latter can be evaluated, for example, with active-space QSE calculation [197],
allowing for recovery of dynamical correlation with accuracy between strongly-
contracted and uncontracted NEVPT2.

###  IV.2 Response functions

An important application of QSMs is the computation of frequency-dependent
response functions,

$C_{AB}(\omega)=\int_{-\infty}^{\infty}\frac{dt}{2\pi}\,e^{i\omega
t}\,\langle\Psi_{0}|\hat{A}e^{-it(\hat{H}-E_{0})}\hat{B}|\Psi_{0}\rangle\;,$
(77)

where $\hat{A},\hat{B}$ are two operators.

When $\hat{A}=\hat{B}={\bf{n}}\cdot\hat{\boldsymbol{\mu}}$ is the component of
the dipole moment along direction ${\bf{n}}$, Eq. (77) is called the dipole
spectral function, and characterises the absorption of ultraviolet and visible
light by a molecule [204, 205, 206, 207, 208, 209, 210]. In solid-state
systems, when $\hat{A}=\hat{c}_{{\bf{p}}\sigma}^{\phantom{\dagger}}$ and
$\hat{B}=\hat{c}_{{\bf{p}}\sigma}^{\dagger}$ for a plane-wave with momentum
${\bf{p}}$ and spin $\sigma$, Eq. (77) is called the quasiparticle spectral
function and it is used to compute the cross-section of angular-resolved
photoemission spectroscopy (ARPES) experiments [211, 212, 213, 214, 215, 216].
Eq. (77) is therefore physically relevant, and challenging to compute as it
involves a time-evolution operator. The Lehmann representation of Eq. (77),

$C_{AB}(\omega)=\sum_{\mu}\delta(\omega-\Delta
E_{\mu})\langle\Psi_{0}|\hat{A}|\Psi_{\mu}\rangle\langle\Psi_{\mu}|\hat{B}|\Psi_{0}\rangle\;,$
(78)

indicates that QSMs are natural strategies to approximate frequency-dependent
response functions, by simply replacing the exact Hamiltonian eigenfunctions
in Eq. (78) with the approximations yielded by a quantum subspace calculation
[6, 217, 218, 219].

###  IV.3 Additional applications

####  IV.3.1 Fast-forwarding time evolution

An important goal of quantum simulations is to understand how the properties
of a physical system evolve over time. In the context of ES, the non-
equilibrium character of time-dependent Hamiltonians (e.g., in the presence of
an external, oscillating electromagnetic field) requires solving the time-
dependent Schrödinger equation

$i\hbar\frac{d}{dt}|\Phi(t)\rangle=\hat{H}(t)|\Phi(t)\rangle$ (79)

to access frequency-dependent polarisabilities and other optical properties
[208, 209, 210, 220, 221, 222, 223].

The simulation of Eq. (79) is a particularly compelling application for a
quantum computer, as this problem lies in the complexity class BQP, meaning
that quantum computers can approximate the state $|\Phi(t)\rangle$ with
accuracy $\varepsilon$ at a polynomial cost in $t,\varepsilon^{-1}$, and
system size. While the complexity of simulating $|\Phi(t)\rangle$ on a quantum
computer is linear in time for a generic quantum system [224, 225], there
exist exceptions to this lower bound: certain Hamiltonians that can be
efficiently diagonalized, e.g. frustration-free Hamiltonians [226] and free
fermions (see Section I.2.6). Time evolution under those Hamiltonians can be
simulated exactly at a cost independent of $t$.

For a generic quantum system, one can approximately solve for
$|\Phi(t)\rangle$ by fast-forwarding time evolution. Within the fast-
forwarding approach, one approximately diagonalizes a short-time propagator
$\hat{U}_{\Delta t}=e^{-i\Delta t\hat{H}}$ and uses knowledge of its
eigenpairs to approximate a generic-time propagator $\hat{U}_{t}$ [227, 228].
While the error of the fast-forwarding procedure scales linearly with time
[228, 229], accurate approximations may be obtained in specific situations:
for example, if $|\Phi(0)\rangle$ has support over $n$ Hamiltonian
eigenstates, $|\Phi(t)\rangle$ can be simulated exactly in the $n$-dimensional
Hamiltonian Krylov space.

QSMs offer a natural avenue to study fast-forwarded time evolution [230, 162].
Examples of the use of QSMs to fast-forward time evolution are the Subspace
Variational Quantum Simulator (SVQS) [231], fixed-state Variational Fast
Forwarding (fs-VFF) [232], and Classical-Quantum Fast Forwarding (CQFF) [168].

##  V Implementation

In this Section, we list some implementations on quantum hardware, outline
important challenges posed by the hardware implementation of QSMs, and discuss
recent research aimed at addressing these challenges.

###  V.1 Hardware implementations

Recent hardware simulations of QSMs for electronic structure applications are
listed in Table 2. For each simulation, we provide details on the studied
system and properties the QSM used, the number of qubits, the depth and number
of $\mathsf{cNOT}$ gates in the simulated circuits, along with the hardware
used. The simulations listed in the table predominantly employed JW mapping as
described in Subsection I.2.6, except for Ref. [6, 197, 233], which used
parity mapping and qubit-tapering techniques [234, 235], and Ref. [236], which
used a first-quantisation mapping. Refs. [219, 237] used a qubit-reduction
technique called entanglement forging [238].

Simulations display improving trends in qubit number, circuit depth, and
number of $\mathsf{cNOT}$ gates. To make substantial progress in this domain
requires (in addition to the continuous refinement of hardware manufacturing
and control) a deeper understanding and effective resolution of the challenges
posed by the implementation of QSMs on quantum hardware. The remainder of this
section addresses some of these challenges and the corresponding research
efforts.

Reference | System | Properties | Algorithm | $N_{q}$ | ($\mathsf{depth},\mathsf{cNOT}$) | Hardware
---|---|---|---|---|---|---
[6] |  H2 | gs/es | MRCISD | 2 | (3,1) | sc
[233] | OLED | gs/es | qEOM | 2 | (3,2) | sc
[239] |  H2 | gs/es | MRCISD | 4 | (9,6) | sc
[158] | stilbene | gs/es | QFD | 5 | (102,78) | sc
[197] |  OH/OH- | gs | MRCISD | 6 | (20,15) | sc
[240] | NV | gs/es | MRCISD | 4 | (7,2) | sc
[236] | NV | gs/es | MRCISD | 4 | (20,14) | sc
[219] |  H3S+ | gs/es/rf | MRCISD | 6 | (21,19) | sc
[241] |  CH4 | gs/es/rf | MRCIS | 6 | (12,7) | ti
[237] |  C4H4N2 | gs/es | MRCISD | 8 | (32,31) | sc
Table 2: List of hardware simulations of QSMs for electronic structure
applications. The abbreviations “gs/es/rf” stand for ground state/excited
states/response functions, the abbreviations “sc/ti” stand for
superconducting/trapped-ion, and the abbreviations “OLED/NV” stand for organic
light-emitting diode and nitrogen-vacancy. Algorithm abbreviations are defined
in the main text.

###  V.2 Effect of shot noise and decoherence

QSMs involve extracting matrix elements of the Hamiltonian between pairs of
subspace basis states, as well as inner products of subspace basis states (see
Fig. 6). As discussed in Subsection V.4.2, due to the probabilistic nature of
quantum operations, repeating quantum measurements on multiple copies of the
quantum circuit is necessary for attaining accurate outcomes and
probabilities. Therefore, unlike classical SMs, where numerical errors only
arise due to machine precision, QSMs are inherently accompanied by errors
originating from a finite number of samples. Moreover, as mentioned above,
QSMs often do not employ orthogonal subspace bases, leading to ill-
conditioning in the classical post-processing. Although this also appears in
some variants of classical SMs (see Subsection II.1), in the noisier context
of QSMs it can have the effect of amplifying errors in the generalised
eigenvalue problem unless handled carefully. The relationship between sampling
noise (or noise from other sources in a quantum algorithm) and its effects on
eigenvalue estimates from QSMs is a subject of intense research [172, 242].

Let us first consider the simplified situation of a QSM targeting an
orthogonal subspace, i.e. $S=\mathbbm{1}$. In that situation, assuming the
computed Hamiltonian matrix $\tilde{H}$ differs from the exact one $H$ by a
perturbation $\Delta_{H}$, the Bauer-Fike theorem [243] provides an upper
bound for the difference $\tilde{E}_{\mu}-E_{\mu}$ between the computed and
exact Hamiltonian eigenvalues,

$\min_{\nu}|\tilde{E}_{\mu}-E_{\nu}|\leq\mathrm{cond}_{p}(V)\|\Delta_{H}\|_{p}\;,$
(80)

where $V$ is the eigenvector matrix that diagonalises $H$,
$\mathrm{cond}_{p}(V)=\|V\|_{p}\|V^{-1}\|_{p}$, and $\|V\|_{p}$ is the
$p$-norm of $V$. This result quantifies the intuitive fact that eigenvalue
perturbations are related to the condition number of the Hamiltonian and the
magnitude of the perturbation.

Compared with the eigenvalue problem, the generalised eigenvalue of a matrix
pair $(H,S)$ tends to be more sensitive to noise. In a generic QSM, the
matrices $(H,S)$ may be estimated incorrectly due to finite sampling, quantum
hardware decoherence, and algorithm-dependent errors. These errors result in
perturbations $\Delta_{H}$, $\Delta_{S}$ of $H,S$. As demonstrated in
literature [242], these perturbations can affect the solution of the GEEV,

$\tilde{H}\tilde{C}_{\mu}=\tilde{S}\tilde{C}_{\mu}\tilde{\tilde{E}}_{\mu}^{(n)}$
(81)

where $\tilde{H}=H+\Delta_{H}$ and $\tilde{S}=S+\Delta_{S}$ are perturbed
Hamiltonian and overlap matrices, respectively,
$\tilde{\tilde{E}}_{\mu}^{(n)}$ is the $\mu$-th perturbed eigenvalue and
$\tilde{C}_{\mu}$ the corresponding eigenvector. Because solving the GEEV
involves the calculation of $\tilde{S}^{-1/2}$, small singular values of
$\tilde{S}$ amplify the noise in the matrix pair significantly. Such cases in
which $\tilde{S}$ has small singular values are called ill-conditioned
problems.

A previous study by Mathias and Li [244] reported an improved perturbation
theory for GEEVs, using a geometrical approach on the complex plane describing
the quadratic form of the problem. Subsequently, Epperly et al [172] used
perturbation theory to describe QFD perturbation with a real-time evolution
ansatz. They also proposed regularizing the GEEV, i.e., reducing the condition
number of $\tilde{S}$, by a thresholding procedure in which the least
significant eigenvectors of $\tilde{S}$ are projected out of both $\tilde{H}$
and $\tilde{S}$, and the corresponding dimensions removed. This yields the
matrices $\tilde{A}=V^{\dagger}_{\varepsilon}\tilde{H}V_{\varepsilon}$ and
$\tilde{B}=V^{\dagger}_{\varepsilon}\tilde{S}V_{\varepsilon}$ where
$V_{\varepsilon}$ is the matrix whose columns are the $n_{\varepsilon}\leq n$
eigenvectors of $\tilde{S}$ with eigenvalues above a threshold
$\varepsilon>0$. Mathias and Li characterised the relationship between
$\Delta_{H}$, $\Delta_{S}$ and $\tilde{\tilde{E}}_{\mu}^{(n)}$. Defining
$(\tilde{A},\tilde{B})=(A+\Delta_{A},B+\Delta_{B})$ and
$\chi^{2}=\|\Delta_{A}\|^{2}+\|\Delta_{B}\|^{2}$, and assuming that

$\begin{split}2n_{\varepsilon}^{2}\chi^{2}&\leq\lambda_{\varepsilon}^{2}\;,\\\
|\tan^{-1}\tilde{E}_{1}^{(n_{\varepsilon})}-\tan^{-1}\tilde{E}_{0}^{(n_{\varepsilon})}|&\geq\sin^{-1}\frac{n_{\varepsilon}\chi}{\lambda_{\varepsilon}}\;\end{split}$
(82)

where $\tilde{E}_{\mu}^{(n_{\varepsilon})}$ are the eigenvalues of the exact
pair $(A,B)$ and $\lambda_{\varepsilon}$ the smallest singular value of $B$,
it follows that the lowest eigenvalue of the perturbed pair
$(\tilde{A},\tilde{B})$ satisfies

$|\tan^{-1}\tilde{\tilde{E}}_{0}^{(n_{\varepsilon})}-\tan^{-1}\tilde{E}_{0}^{(n_{\varepsilon})}|\leq\sin^{-1}\frac{\sqrt{2}n_{\varepsilon}\chi}{d_{0}}\;,$
(83)

where $\tilde{E}_{0}^{(n_{\varepsilon})}$ is the lowest energy of the
unperturbed subspace subject to the same thresholding procedure. Here
$d_{0}^{-1}$ is the condition number of
$\tan^{-1}\tilde{E}_{0}^{(n_{\varepsilon})}$, given by

$d_{0}=|{\bf{x_{0}}}\cdot(A+iB){\bf{x_{0}}}|,$ (84)

where ${\bf{x_{0}}}$ is the unit-norm eigenvector of $(A,B)$ with the lowest
eigenvalue. Epperly et al [172] obtained a different bound, based on the
perturbations $\Delta_{H}$, $\Delta_{S}$,

$|\tan^{-1}\tilde{\tilde{E}}_{0}^{(n_{\varepsilon})}-\tan^{-1}\tilde{E}_{0}^{(n_{\varepsilon})}|\leq\mathcal{O}\left(\eta^{1/(1+\alpha)}d_{0}^{-1}\right)\;,$
(85)

where $\alpha$ is a constant ranging from 0 to $1/2$, and
$\eta^{2}=\|\Delta_{H}\|^{2}+\|\Delta_{S}\|^{2}$ is related to $\varepsilon$
and $\chi$ by $\chi\leq\mathcal{O}\left(\eta^{1/(1+\alpha)}/n\right)$ and
$\varepsilon=\mathcal{O}\left(\eta^{1/(1+\alpha)}\right)$. In summary, the
perturbation bound in Eq. (85) indicates that the perturbation error is
sublinear to the error matrix norms and condition number after the truncation
of the basis. Thus, with additional information about the error matrix norms
$\eta$, one can establish a sampling error analysis for QSMs.

Alternative approaches for post-processing the quantum data have been
suggested which circumvent solving the GEEV entirely and thus don’t require
dealing with the ill-conditioning of the overlap matrix. Inspiring by signal
processing techniques, the quantum exponential least squares routine (QCELS)
[245, 246] solves a nonlinear least-squares problem to approximate the best
amplitude and phase parameters fitting the measured overlap matrix data. QCELS
has been shown to converge quickly for ground state estimation when the
initial state has large overlap with the ground state. Another approach,
observable dynamic mode decomposition (ODMD) [247], builds off of standard
dynamic mode decomposition (DMD) [248, 249], which was originally developing
for approximating classical dynamical systems. ODMD requires measuring only
the real or imaginary part of the overlap matrix and embedding this data into
a pair of time shifted Hankel matrices. These Hankel matrices can be used to
construct the DMD matrix through a least squares formulation, whose
eigenvalues approximate the ground state energy, avoiding the GEEV and thus
improving conditioning and stability. ODMD has shown to have better
convergence than many competing methods in the case of low initial overlap
with the ground state.

###  V.3 Error mitigation and generalised QSE

In the previous Subsection, we discussed how decoherence alters the
expectation values and variances of random variables sampled on a quantum
computer, which in turn impacts the accuracy and precision of a simulation. To
remedy this limitation, techniques for the mitigation of readout [52, 53] and
gate [54, 55, 56] errors were introduced.

Remarkably, QSMs are known to possess error-mitigating properties [250]. In
particular, QSE was predicted and experimentally confirmed to approximate
excited states and reduce errors by performing additional measurements and
solving an eigenvalue problem [6]. The observation that QSE measurement can be
used to mitigate errors lies at the core of the generalised QSE approach
[160]. Consider a quantum circuit that would prepare a pure state
$|\Phi\rangle$ in the absence of quantum noise, but instead prepares a density
operator $\hat{\rho}$. Then, by measuring a set of operators
$\\{\hat{\sigma}_{\alpha}\\}_{\alpha}$, one can construct a state of the form

$\rho_{\mathrm{EM}}=\frac{\hat{A}^{\dagger}\hat{P}\hat{A}}{\mbox{Tr}[\hat{A}^{\dagger}\hat{P}\hat{A}]}\;,\;\hat{A}=\sum_{\alpha}C_{\alpha}\hat{\sigma}_{\alpha}\;,$
(86)

where $\hat{P}$ is a positive operator that is practically taken as
$\mathbbm{1}$ or $\hat{\rho}$, and the coefficients
$\\{C_{\alpha}\\}_{\alpha}$ are determined, if the goal is to approximate the
ground state of a Hamiltonian, by solving the eigenvalue equation

$HC=SC\tilde{E}\;,\;H_{\alpha\beta}=\mbox{Tr}[\hat{\sigma}_{\alpha}^{\dagger}\hat{P}\hat{\sigma}_{\beta}\hat{H}]\;,\;S_{\alpha\beta}=\mbox{Tr}[\hat{\sigma}_{\alpha}^{\dagger}\hat{P}\hat{\sigma}_{\beta}]\;.$
(87)

An alternative to the measurement of Pauli operators is the use of powers
[160, 251] of the density operators,
$\rho_{\mathrm{EM}}=\sum_{\alpha,\beta=0}^{n-1}C_{\alpha}^{*}C_{\beta}\,\hat{P}^{\alpha+\beta}$,
where $\hat{P}\in\\{\mathbbm{1},\hat{\rho}\\}$ is chosen according to
$(\alpha+\beta)\%2$. The reason behind this choice is the fact that raising
$\hat{\rho}$ to powers suppresses non-dominant eigenvalues, thereby mitigating
errors. Another alternative is the fault subspace, where a target density
operator $\hat{\rho}$ can be prepared approximately, i.e. one can produce a
set of density operators $\hat{\rho}(\eta_{\ell})$ characterised by variable
error rates $\\{\eta_{\ell}\\}_{\ell}$ corresponding to the amplification of
achievable error rates (e.g. by gate repetition, probabilistic error
amplification, decoherence amplification, or cross-talk boost). One can then
produce a state of the form
$\rho_{\mathrm{EM}}\propto\hat{\sigma}^{\dagger}\hat{\sigma}$ where
$\hat{\sigma}=\sum_{\ell}C_{\ell}\,\hat{\rho}(\eta_{\ell})$ for suitable
coefficients $C_{\ell}$ [160, 252], and $\hat{P}=\mathbbm{1}$ is assumed in
Eq. (86). It should be noted that, when $\rho_{\mathrm{EM}}$ involves the
$n$-th power of $\hat{\rho}$, an ancilla is required to prepare
$\hat{\rho}^{n}$ via controlled $\mathsf{SWAP}$ operations over $n$ copies of
$\hat{\rho}$ [253]. While it is challenging to require multiple copies of a
state $\hat{\rho}$, this issue can be addressed by running deeper quantum
circuits, i.e., of depth increased by a factor of $n$ [254].

However, numerical evidence suggests generalised QSE is capable of mitigating
various stochastic, coherent, and algorithmic errors [160], making it a
compelling error mitigation technique for quantum devices. Of particular
relevance is the possibility to measure, in addition to the Hamiltonian,
symmetry operators, e.g. the total spin, $\hat{S}^{2}$, and its component
along the $z$-axis, $\hat{S}_{z}$. This is because the measurement of symmetry
operators allows modifying the coefficients $\\{C_{\alpha}\\}_{\alpha}$, thus
mitigating errors arising from symmetry breaking.

###  V.4 Optimisation of quantum circuits

As sketched in Figs. 8, 9, and 10, implementing a QSM may require ancillae and
deep circuits containing controlled unitary operations. Furthermore, a
considerable overhead of quantum measurements may be needed, see e.g. Fig. 7.
The economisation of these operations stands to impact the successful
application of QSMs to ES problems. In this Subsection, we briefly describe
some recent research aimed at economising quantum circuits and measurements.

Figure 11: Quantum circuit implementing the generalised $\mathsf{SWAP}$
network (a), and quantum circuits to measure matrix elements without Hadamard
test (b).

####  V.4.1 Circumventing the Hadamard test

QSMs may require the measurement of quantities of the form
$O_{\alpha\beta}=\langle{\bf{v}}_{0}|\hat{U}^{\dagger}_{\alpha}\hat{O}\hat{U}_{\beta}|{\bf{v}}_{0}\rangle$,
where $\hat{U}_{\alpha}$, $\hat{U}_{\beta}$ are unitary operators, and
$\hat{O}$ is a Hermitian operator (e.g. the Hamiltonian or the identity).
Here, without loss of generality, we will assume that
$\hat{O}=\hat{V}^{\dagger}\hat{\Delta}\hat{V}$ for a unitary $\hat{V}$ and a
diagonal Pauli operator $\hat{\Delta}$.

In general, the measurement of $O_{\alpha\beta}$ requires an ancilla qubit and
the application of the controlled versions of $\hat{U}_{\alpha}$ and
$\hat{U}_{\beta}$ (see e.g. Figs. 9 and 10), a protocol known as Hadamard test
in QC literature [63, 255]. The controlled $\hat{U}_{\alpha}$ and
$\hat{U}_{\beta}$ operations are particularly expensive: every single-qubit
and $\mathsf{cNOT}$ gate in a circuit representation of $\hat{U}_{\alpha}$ and
$\hat{U}_{\beta}$ is replaced by a controlled-single-qubit operation
(requiring 2 $\mathsf{cNOT}$ gates) and a Toffoli gate (requiring 3
$\mathsf{cNOT}$ gates) [31] respectively. Furthermore, unless the device has
all-to-all connectivity, implementing the controlled unitaries involves a
considerable overhead of $\mathsf{SWAP}$ gates, as discussed in Subsection
I.2.5.

As a way to bypass the Hadamard test, Huggins et alproposed [7] a generalised
$\mathsf{SWAP}$ network, sketched in Fig. 11a. This method requires the
existence of a state $|\mbox{\O}\rangle$ that is an eigenstate of
$\hat{U}_{\ell}$ with $\ell=\alpha,\beta$ (a typical example in ES is the
vacuum state). Under such assumption, the application of controlled
$\hat{U}_{\alpha}$, $\hat{U}_{\beta}$ can be replaced with two controlled
multi-qubit $\mathsf{SWAP}$ gates. While this method does not remove the need
for an ancilla and controlled operations, it makes the computational cost of
the computation of $O_{\alpha\beta}$ independent of the structure of
$\hat{U}_{\alpha}$, $\hat{U}_{\beta}$. A further improvement, sketched in Fig.
11b, is due to Cortes and Gray [161]. Under the same assumption, and the
requirement that a superposition of the form
$|{\bf{v}}_{+}\rangle=\left(|{\bf{v}}_{0}\rangle+|\mbox{\O}\rangle\right)/\sqrt{2}$
can be prepared efficiently, one can evaluate $O_{\alpha\beta}$ from the
quantities

$\begin{split}f_{0}&=|\langle\mbox{\O}|\hat{U}^{\dagger}_{\alpha}\hat{O}\hat{U}_{\beta}|\mbox{\O}\rangle|^{2}\;,\\\
f_{1}&=|\langle{\bf{v}}_{0}|\hat{U}^{\dagger}_{\alpha}\hat{O}\hat{U}_{\beta}|{\bf{v}}_{0}\rangle|^{2}\;,\\\
f_{2}&=|\langle{\bf{v}}_{+}|\hat{U}^{\dagger}_{\alpha}\hat{O}\hat{U}_{\beta}|{\bf{v}}_{+}\rangle|^{2}\;.\\\
\end{split}$ (88)

Indeed, writing $f_{2}$ in terms of $|{\bf{v}}_{0}\rangle$ and
$|\mbox{\O}\rangle$, and writing in polar form the complex numbers
$\langle\mbox{\O}|\hat{U}^{\dagger}_{\alpha}\hat{O}\hat{U}_{\beta}|\mbox{\O}\rangle=\sqrt{f_{0}}e^{i\theta_{0}}$
and $O_{\alpha\beta}=\sqrt{f_{1}}e^{i\theta_{1}}$ (where the absolute values
$\sqrt{f_{0}}$ and $\sqrt{f_{1}}$ are measurable and the phases $\theta_{0}$
and $\theta_{1}$ are respectively known and unknown), one obtains the relation
$4f_{2}=f_{1}+f_{0}+2\sqrt{f_{0}f_{1}}\cos(\theta_{1}-\theta_{0})$.

####  V.4.2 Quantum measurements

An important challenge to overcome for the practicality of QSMs is improving
the computational cost of quantum measurements. This is particularly important
because the measurement requirements in QSMs are typically much higher than in
other quantum computing algorithms (e.g. VQE) since one must measure the
expectation values of multiple operators in addition to that of the electronic
Hamiltonian. Numerous quantum measurement techniques have been proposed
recently, and their adoption and refinement stand to benefit the practicality
of QSMs.

In the case of MRCISD, Takeshita et al [157] noted that the expectation values
that one needs to measure in order to construct the subspace spanned by
singles and doubles on top of a trial wavefunction are equal to contractions
between $k$-body density matrices with $k=1\dots 6$ and suitable coefficients
defining the electronic Hamiltonian. This high computational cost prompted
them to propose the use of a cumulant approximation [256] to express $k$-body
density matrices as linear combinations of antisymmetrised products of
cumulant operators, e.g.

$\begin{split}\rho^{(1)}&=\Delta^{(1)}\\\
\rho^{(2)}&=\Delta^{(2)}+\Delta^{(1)}\wedge\Delta^{(1)}\\\
\rho^{(3)}&=\Delta^{(3)}+\Delta^{(1)}\wedge\Delta^{(1)}\wedge\Delta^{(1)}+3\Delta^{(2)}\wedge\Delta^{(1)}\\\
\dots{}\end{split}$ (89)

where $\rho^{(k)}$ and $\Delta^{(k)}$ are the $k$-body density matrix and
cumulant respectively, and “$\wedge$” denotes an antisymmetrised tensor
product. Measuring density matrices of order up to $l$ allows computation of
cumulants of order up to $l$. By then requiring that cumulants of order $k>l$
vanish, one can approximate density matrices of order $k>l$ using cumulants of
order up to $k$ only. The cumulant approximation reduces the computational
cost of MRCISD calculations but introduces approximations, that affect
physical properties computed with MRCISD. Furthermore, it is not generalizable
to other QSMs.

A more general study was carried out by Choi et al [257]. Their starting point
is the observation that the operators to be measured in a QSE calculation,
$\hat{A}_{d}\in\\{\hat{O}^{\dagger}_{\alpha}\hat{H}\hat{O}_{\beta},\hat{O}^{\dagger}_{\alpha}\hat{O}_{\beta}\\}_{\alpha\beta}$,
can be written as linear combinations of Pauli operators,
$\hat{A}_{d}=\sum_{{\bf{m}}\in A_{d}}c_{{\bf{m}}d}\,\sigma_{\bf{m}}$, where
$A_{d}$ is the set of indices labeling Pauli operators that appear in
$\hat{A}_{d}$. The standard approach to economising such measurement is “Pauli
grouping”, wherein Pauli operators with indices $A=\cup_{d}A_{d}$ are divided
into sets $D_{f}$ of mutually commuting (and therefore jointly measurable)
Pauli operators, $A=\cup_{f}D_{f}$ with $[\sigma_{\bf{m}},\sigma_{\bf{n}}]=0$
for all ${\bf{m}},{\bf{n}}\in D_{f}$. Expectation values can then be
determined by writing operators as $\hat{A}_{d}=\sum_{f}\sum_{{\bf{m}}\in
A_{d}\cap
D_{f}}c^{(f)}_{{\bf{m}}d}\,\sigma_{\bf{m}}=\sum_{f}\overline{A}^{(f)}_{d}$
(note that coefficients $c^{(f)}_{{\bf{m}}d}$ and $c_{{\bf{m}}d}$ may differ
when the sets $D_{f}$ have non-zero intersection. Each set $D_{f}$ has an
associated unitary transformation $\hat{W}_{f}$ that turns Pauli operators in
the set into $Z$ (diagonal) Pauli operators, which can be measured jointly,
obtaining sample averages $\overline{\sigma}^{(f)}_{\bf{m}}$. By classical
post-processing, one can obtain unbiased estimates of expectation values,

$\overline{A}_{d}=\sum_{f}\sum_{{\bf{m}}\in A_{d}\cap