#  Quantum algorithms for generator coordinate methods

Muqing Zheng  Lehigh University, Bethlehem, Pennsylvania, 18015, USA  Pacific
Northwest National Laboratory, Richland, Washington, 99354, USA  Bo Peng
Pacific Northwest National Laboratory, Richland, Washington, 99354, USA
Nathan Wiebe  University of Toronto, Toronto, Ontario, M5G 1Z8, Canada
Pacific Northwest National Laboratory, Richland, Washington, 99354, USA  Ang
Li  Pacific Northwest National Laboratory, Richland, Washington, 99354, USA
Xiu Yang  Lehigh University, Bethlehem, Pennsylvania, 18015, USA  Karol
Kowalski<EMAIL_ADDRESS>Pacific Northwest National Laboratory,
Richland, Washington, 99354, USA

###### Abstract

This paper discusses quantum algorithms for the generator coordinate method
(GCM) that can be used to benchmark molecular systems. The GCM formalism
defined by exponential operators with exponents defined through generators of
the Fermionic $U(N)$ Lie algebra (Thouless theorem) offers a possibility of
probing large sub-spaces using low-depth quantum circuits. In the present
studies, we illustrate the performance of the quantum algorithm for
constructing a discretized form of the Hill-Wheeler equation for ground and
excited state energies. We also generalize the standard GCM formulation to
multi-product extension that when collective paths are properly probed, can
systematically introduce higher rank effects and provide elementary mechanisms
for symmetry purification when generator states break the spatial or spin
symmetries. The GCM quantum algorithms also can be viewed as an alternative to
existing variational quantum eigensolvers, where multi-step classical
optimization algorithms are replaced by a single-step procedure for solving
the Hill-Wheeler eigenvalue problem.

##  I Introduction

The rapid development of quantum technologies and quantum algorithms addresses
long-standing computational challenges of many-body physics and quantum
chemistry. While the primary target is to overcome exponential growth in
complexity associated with approaching the exact limit in the simulations, the
possibility of identifying physically meaningful solutions to problems of
interest is equally important. Although Quantum Phase Estimation algorithms
Luis and Peřina (1996); Cleve  _et al._ (1998); Berry  _et al._ (2007); Childs
(2010); Wecker, Hastings, and Troyer (2015a); Häner  _et al._ (2016); Poulin
_et al._ (2017) have been designed in a way that adequately addresses both of
these issues, their applicability is currently limited by the necessity of
using complex quantum circuits with corresponding depths that preclude its
practical applications on existing quantum computing platforms, dominated by
Noisy Intermediate-Scale Quantum (NISQ) devices. Instead, hybrid algorithms
such as various Variational Quantum Eigensolver (VQE) Peruzzo  _et al._
(2014); McClean  _et al._ (2016); Romero  _et al._ (2018); Shen  _et al._
(2017); Kandala  _et al._ (2017, 2019); Colless  _et al._ (2018); Huggins  _et
al._ (2020); Ryabinkin  _et al._ (2018); Cao  _et al._ (2019); Ryabinkin  _et
al._ (2020); Izmaylov  _et al._ (2019); Lang, Ryabinkin, and Izmaylov (2021);
Grimsley  _et al._ (2019a, b); Cerezo  _et al._ (2021); McArdle  _et al._
(2020); Bharti  _et al._ (2022); Anand  _et al._ (2022) are currently being
intensively tested on NISQ quantum computers to characterize the properties of
correlated quantum systems. In this effort, in many aspects, the VQE
formulations, for example, based on the unitary coupled-cluster (UCC)Bartlett,
Kucharski, and Noga (1989); Taube and Bartlett (2006); Hoffmann and Simons
(1988); Kutzelnigg (1991); Evangelista, Chan, and Scuseria (2019)
representation of the wave function, mirror the standard conventional
formulations of CC theory, where a large number of excitations are included in
the cluster operator and simultaneously optimized in the iterative process.
This algorithm leads to another set of challenges associated with the
potential problems with the convergence of iterative process (commonly
referred to as the barren minimum problem) and representation of UCC Ansatz on
quantum registers, which may result in long Trotter-like products of
exponential operators defined by multi-qubit gates. Although several
strategies mitigating these problems have recently been proposed, utilizing
VQE-UCC formulations and extending these methods beyond specific system-size
limits may be challenging and require reformulation of the quantum many-
problem into recently introduced quantum flow equations,Kowalski (2021) where
large subspaces of Hilbert space can be sampled through the constant-depth
small-dimensionality coupled eigenproblems.

Instead, in this paper, we explore the applicability of the Generator
Coordinate Method (GCM) Hill and Wheeler (1953); Rodríguez-Guzmán, Egido, and
Robledo (2002); Bender, Heenen, and Reinhard (2003); Ring and Schuck (2004);
Yao  _et al._ (2010); Egido (2016); Hizawa, Hagino, and Yoshida (2021) as an
alternative to popular VQE formulations. The main difference with the VQE
method is that the GCM method avoids highly-nonlinear parametrization of the
wave function and provides an efficient mean for direct extension of the
probed subspaces. Additionally, it offers an efficient utilization of Ansatzes
represented by low-depth quantum circuits.

The GCM method was one of the first attempts to combine two distinct aspects
of many-body theories: independent-particle models and theories describing
collective phenomena, where the approximate eigenstates $|\Psi_{\rm
GCM}\rangle$ of the Hamiltonian $H$ are expressed using a family of $N$-body
wave functions $|\Phi({\bf q})\rangle$,

$|\Psi_{\rm GCM}\rangle=\int d{\bf q}|\Phi({\bf q})\rangle f({\bf q})\;,$ (1)

where ${\bf q}$ is a set of collective variables that describe correlation
effects in many-body systems, and usually, the corresponding $|\Phi({\bf
q})\rangle$ is represented a complicated linear combination of Slater
determinants. The scalar $f({\bf q})$ is referred to as the weight function.
The advantage of the GCM approach is obtaining ground states and classes of
excited states described by the chosen set of generator coordinates. In
general, in the analogy to the coherent state representation,Klauder and
Skagerstam (1985); Zhang, Gilmore  _et al._ (1990) the family $|\Phi({\bf
q})\rangle$ forms an over-complete basis. Upon substituting (1) into
Schödinger equation one gets the so-called Hill-Wheeler integral equations for
unknown $f({\bf q})$ coefficients

$\int d{\bf q}’[{\bf H}({\bf q},{\bf q}’)–E{\bf S}({\bf q},{\bf q}’)]f({\bf
q}’)=0\;,$ (2)

where the integral kernels ${\bf H}$ and ${\bf S}$ are defined as

$\displaystyle{\bf H}({\bf q},{\bf q}’)$ $\displaystyle=$
$\displaystyle\langle\Phi({\bf q})|H|\Phi({\bf q}^{\prime})\rangle\;,$ (3)
$\displaystyle{\bf S}({\bf q},{\bf q}’)$ $\displaystyle=$
$\displaystyle\langle\Phi({\bf q})|\Phi({\bf q}^{\prime})\rangle\;.$ (4)

In typical applications, the Hill-Wheeler equation is usually solved
numerically by discretization, which transforms integral equation (3) into an
algebraic eigenvalue problem.

There are two categories of GCM formulations in applications to many-body
quantum systems. The first category’s purpose is to restore broken symmetries
of the $|\Phi({\bf q})\rangle$ states. For example, Bardeen-Cooper-Schrieffer
states are not eigenstates of the particle number operator $N$. To project the
Bardeen-Cooper-Schrieffer states onto the Hilbert space with the desired
number of particles $N_{0}$, one uses the projection operator expressed in
terms of integral over gauge angle $\phi_{N_{0}}$,

$P_{N_{0}}=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi_{N_{0}}e^{i\phi_{N_{0}}(N-N_{0})}\;.$
(5)

In this case, the properties of $f_{k}({\bf q})$ coefficients are determined
by the properties of the symmetry (projection) operators (e.g., particle
number projection operator). In the second category of GCM formulations, the
unknown weight function is optimized to capture correlation effects encoded in
generator coordinates. However, designing an adequate grid or path to probe
the Hill-Wheeler equation is a rather empirical procedure that requires much
intuition and prior knowledge of the sought-after many-body system.

The appealing feature of the GCM method, especially from the point of possible
quantum computing applications, is the possibility of combining low-depth
representations of $|\Phi({\bf q})\rangle$ functions with simple, one-step
optimization conditions for weight function. In this approach, the role of
quantum computing is to map a discrete number of states $\\{|\Phi({\bf
q}_{p})\rangle\\}_{p=1}^{M}$ to quantum register and to evaluate matrix
elements for Hamiltonian and overlap matrices (Eqs. (3) and (4)). In contrast,
solving a generalized eigenvalue problem in a discrete basis representation
takes place only once on a classical computer, avoiding multiple instances of
quantum-classical communication as in the VQE formalism. In the discussed
formalism, we follow an “algebraic” GCM formulation discussed by Fukutome in
Ref.Fukutome, 1981a, extend the standard GCM to multi-product exponential
formulas, and provide an algorithm for sampling coordinate space in a way that
provides selective approaching classes of excited Slater determinants. In this
context, the multi-product GCM formulation alleviates some problems associated
with the usage of high-fidelity of Trotter-type expansions.

We demonstrate the performance of the quantum GCM algorithm on the example of
the H4 benchmark system in various configurations. We show that by a judicious
choice of the $|\Phi({\bf q})\rangle$, one can recover high-level of accuracy
both in weakly and strongly correlated regimes using $|\Phi({\bf q})\rangle$s
that use the manifold of single excitations (or $U(N)$ Lie algebra
generators). The obtained level of accuracy is similar to the one obtained
with the advanced VQE formulations.

##  II Theory

The GCM was initially introduced to describe collective effects in nuclei.Hill
and Wheeler (1953); Goeke and Reinhard (1980); Ring and Schuck (2004); Hizawa,
Hagino, and Yoshida (2021) Fukutome in his seminal paper Fukutome (1981a)
considered the GCM from the Lie-algebra theoretical standpoint. We use that
language throughout this paper. Let us assume that the Fermion system is
described by annihilation and creation operators $a_{p}$ and $a_{q}^{\dagger}$
that satisfy the following set of anticommutation relations:

$[a_{p},a_{q}]_{+}=[a_{p}^{\dagger},a_{q}^{\dagger}]_{+}=0\;,\;\;[a_{p},a_{q}^{\dagger}]_{+}=\delta_{pq}$
(6)

As discussed in Ref.Fukutome, 1981a there are several Fermion algebras
including $U(N)$, $SO(2N)$, and $SO(2N+1)$ Lie algebras, and Clifford algebras
that can be used to characterized approximate many-body wave functions (here,
$N$ stands for the number of single particle states). These algebras can be
defined by the following set of operators

$E^{p}_{q}=a_{p}^{\dagger}a_{q}-\frac{1}{2}\delta_{pq}\;,\;E_{pq}=a_{p}a_{q}\;,\;E^{pq}=a_{p}^{\dagger}a_{q}^{\dagger}\;.$
(7)

For example,

  * •

$U(N)$: $\\{E^{p}_{q}\\}$ ,

  * •

$SO(2N)$: $\\{E^{p}_{q},E_{pq},E^{pq}\\}$ ,

  * •

$SO(2N+1)$: $\\{a_{p},a_{q}^{\dagger},E^{p}_{q},E_{pq},E^{pq}\\}$ .

The $U(N)$ algebra is the only algebra where particle number operator commutes
with all operators belonging to $U(N)$ algebra.

Let us focus attention on the $U(N)$ algebra. By $\Gamma({\bf Z})$ we
designate an anti-Hermitian operator defined as

$\Gamma({\bf
Z})=\sum_{p,q}z_{pq}E^{p}_{q}=\sum_{p,q}z_{pq}a_{p}^{\dagger}a_{q}\;\;,\;\;z_{qp}^{\star}=-z_{pq}$
(8)

where ${\bf Z}$ can be viewed as an anti-Hermitian matrix $[z_{pq}]$ where
$z_{qp}^{\star}=-z_{pq}$ or set of indexed parameter $\\{z_{pq}\\}$ satisfying
$z_{qp}^{\star}=-z_{pq}$. Canonical transformation $U({\bf Z})$ generated by
$\Gamma({\bf Z})$ takes the form

$U({\bf Z})=e^{\Gamma({\bf Z})}\;.$ (9)

Standard canonical Thouless transformation of the Hartree-Fock (HF)
determinant $|\Phi\rangle$ can be obtained by acting with $U({\bf Z})$ onto
$|\Phi\rangle$, i.e.,

$|\Phi({\bf Z})\rangle=U({\bf Z})|\Phi\rangle=e^{\Gamma({\bf Z})}|\Phi\rangle$
(10)

The standard unitary CC model with singles (UCCS) is a special case of (10)
where $z_{ij}=z_{ab}=0$, where $i,j,\ldots$ and $a,b,\ldots$ correspond to
spin-orbital indices occupied and unoccupied in the Slater determinant
$|\Phi\rangle$ (where ${\bf Z}$ can be identified with the set of elements
$\\{z_{ia}\\}$). In comparison to $|\Phi\rangle$, $|\Phi({\bf Z})\rangle$ can
contain certain classes of correlation effects (that are ${\bf Z}$ dependent;
it can contain elements of symmetry breaking and effects responsible for the
appearance of HF instabilities).

The parametrized states (10) can be viewed as a non-orthogonal (in general
overcomplete) basis in the Hilbert space. The generator coordinate method
utilizes this fact by representing the wave function in the form

$|\Psi_{\rm GCM}\rangle=\int d{\bf Z}|\Phi({\bf Z})\rangle f({\bf Z})\;,$ (11)

where collective variables ${\bf q}$ of Eq.(1) are now identified with ${\bf
Z}$ matrix, i.e.,

${\bf q}\rightarrow{\bf Z}\;.$ (12)

A typical way for solving Hill-Wheeler equations (2) is through discretization
of the ${\bf Z}$ domain. Let us introduce a set of the ${\bf Z}$-points -
$Q=\\{{\bf Z}_{i}\\}_{i=1}^{M}$ then equations (2) takes the form of non-
orthogonal eigenvalue problem

${\bf H}{\bf f}=E{\bf S}{\bf f},$ (13)

where ${\bf H}$ and ${\bf S}$ ($M\times M$ matrices) and ${\bf f}$
($M$-dimensional vector) are defined as follows:

$\displaystyle{\bf H}_{pq}$ $\displaystyle=$ $\displaystyle\langle\Phi({\bf
Z}_{p})|H|\Phi({\bf Z}_{q})\rangle\;,$ (14) $\displaystyle{\bf S}_{pq}$
$\displaystyle=$ $\displaystyle\langle\Phi({\bf Z}_{p})|\Phi({\bf
Z}_{q})\rangle\;,$ (15) $\displaystyle{\bf f}_{p}$ $\displaystyle=$
$\displaystyle f({\bf Z}_{p})\;.$ (16)

Using this form of discretization, the optimal form of wave function (we
assume the ground-state wave function in this paper) is given by the expansion

$|\Psi_{GCM}\rangle\simeq\sum_{p=1}^{M}{\bf f}_{p}e^{\Gamma({\bf
Z}_{p})}|\Phi\rangle$ (17)

which is reminiscent of recently discussed non-orthogonal variational
approaches discussed in the context of quantum computing.Huggins  _et al._
(2020) If $M$ is equal for a given basis set to a dimension of the full
configuration interaction (FCI) problem (for a given spin and spatial
symmetry) and $e^{\Gamma({\bf Z}_{p})}|\Phi\rangle$ are linearly independent,
then the expansion (17) with optimized $f_{i}$ coefficients describes the
exact electronic wave function.

##  III Multi-product extension of the Generator Coordinate formalism

A possible extension of the GCM expansion given by Eq.(11) can be provided by
the expansion involving multiple products (the product GCM formalism,
abbreviated as PGCM(k)) of $k$ exponential operator (which is inspired by a
recent progress in the development of dynamical GCM methods Goeke and Reinhard
(1978); Reinhard and Goeke (1978a, b); Goeke and Reinhard (1980); Hizawa,
Hagino, and Yoshida (2021)). For example, one can introduce the following
expansion

$|\Psi_{\rm PGCM}^{(k)}\rangle=\int d{\bf Z}(1)\ldots d{\bf
Z}(k)|\Phi^{(k)}({\bf Z}(1),\ldots,{\bf Z}(k))\rangle f({\bf Z}(1),\ldots,{\bf
Z}(k))$ (18)

where

$|\Phi^{(k)}({\bf Z}(1),\ldots,{\bf Z}(k))\rangle=e^{\Gamma({\bf Z}(k))}\ldots
e^{\Gamma({\bf Z}(1))}|\Phi\rangle\;.$ (19)

In the above representation, $\Gamma({\bf Z}(i))$ can belong to various Lie
algebras. In the following we will focus on the case where all $\Gamma({\bf
Z}(i))$’s ($i=1,\ldots,k)$ belong to the same $U(N)$ Lie algebra; that is,

$\Gamma({\bf
Z}(i))=\sum_{pq}z_{pq}(i)E^{p}_{q}\;,\;z_{qp}(i)^{\star}=-z_{pq}(i)$ (20)

In fact, the PGCM formula may be viewed as a special case of the GCM where

${\bf q}\rightarrow{\bf Z}(1)\times\ldots\times{\bf Z}(k)\;.$ (21)

We will demonstrate that the $U(N)$ case of PGCM lends itself for an efficient
way of representing higher-rank excitations in quantum computing. In
particular, this goal can be achieved by using a simple algorithm for
discretization of ${\bf Z}(i)$ domains in the GCM method. In particular, we
can show that PGCM(k) can be used to approximate $2k$-tuple excitations in the
configuration-interaction-type expansion.

Let us focus on the specific case when $k=2$ (i.e., the PGCM(2)formalism). In
this case we will represent the $|\Psi^{(2)}_{\rm PGCM}\rangle$ wave function,
given by the formulas:

$|\Psi_{\rm PGFM}^{(2)}\rangle=\int d{\bf Z}(1)d{\bf Z}(2)|\Phi^{(2)}({\bf
Z}(1),{\bf Z}(2))\rangle f({\bf Z}(1),{\bf Z}(2))$ (22)

and

$|\Phi^{(2)}({\bf Z}(1),{\bf Z}(2))\rangle=e^{\Gamma({\bf
Z}(2))}e^{\Gamma({\bf Z}(1))}|\Phi\rangle\;.$ (23)

The GCM algorithm can be adapted easily for the product representation of the
trial wave functions. Now the discretization of the problem, in analogy to
Eq.(13), involves a $Q$ set defined as $Q=\\{{\bf Z}_{I}\\}_{I}$, where ${\bf
Z}_{I}={\bf Z}(1)_{p}\times{\bf Z}(2)_{q}$, where we use composite index
$I=(p,q)$. Controllable sampling algorithms of the sub-spaces of the Hilbert
space corresponding to higher-rank excitations with the PGCM formalism based
on the $U(N)$ algebras require careful selection of the sampling points, which
will be discussed in the following section in the context of PGCM(2) method
applications to the H4 system.

##  IV H4 model: the choice of the GCM sampling points

Figure 1: Schematic representation of the orbital energies and the enumeration
scheme for the corresponding spin-orbitals for the H4 model in STO-3G basis
set.Hehre, Stewart, and Pople (1969) Spin-orbitals 1,2 and 3,4 correspond to
occupied $\alpha$ and $\beta$ spin-orbitals, respectively, while spin-orbitals
5,6 and 7,8 correspond to virtual $\alpha$ and $\beta$ spin-orbitals.

We use the H4 model of Ref.Jankowski and Paldus, 1980 as a benchmark system.
The system’s geometry can be defined by a single parameter $\alpha$. For
$\alpha=0.5$, the H4 model corresponds to a linear chain of hydrogen atoms
with distances between adjacent hydrogen atoms equal to 2.0 a.u.; for
$\alpha=0.005$, the system is almost in a square configuration. The main
difference between $\alpha=0.5$ and $\alpha=0.005$ H4 models is in the
structure of the corresponding ground-state wave function. While for
$\alpha=0.5$, the ground-state wave function is dominated by the restricted
Hartree-Fock (RHF) determinant, for $\alpha=0.005$, the ground state is quasi-
degenerate. The contribution of the RHF Slater determinant is almost the same
as the doubly excited configuration where two electrons from the highest
occupied orbital are promoted to the lowest lying virtual orbital. The spin-
orbital numbering scheme is shown in Fig.1. In the PGCM(2) formalism for H4,
we will adopt the following sampling scheme, where the choice of the cluster
operators is consistent with analyzing HF equation stability conditions based
on the Thouless theorem Thouless (1960); Čížek and Paldus (1967); Seeger and
Pople (1977); Fukutome (1981b):

$\displaystyle|R_{0}\rangle$ $\displaystyle=$ $\displaystyle|\Phi\rangle\;,$
(24) $\displaystyle|R_{1}^{(\pm)}(t_{1})\rangle$ $\displaystyle=$
$\displaystyle e^{\pm t_{1}R_{1}}|\Phi\rangle=e^{\pm
t_{1}\\{(a_{5}^{\dagger}a_{2}+a_{7}^{\dagger}a_{4})-(a_{2}^{\dagger}a_{5}+a_{4}^{\dagger}a_{7})\\}}|\Phi\rangle\;,$
(25) $\displaystyle|R_{2}^{(\pm)}(t_{2})\rangle$ $\displaystyle=$
$\displaystyle e^{\pm t_{2}R_{2}}|\Phi\rangle=e^{\pm
t_{2}\\{(a_{6}^{\dagger}a_{1}+a_{8}^{\dagger}a_{3})-(a_{1}^{\dagger}a_{6}+a_{3}^{\dagger}a_{8})\\}}|\Phi\rangle\;,$
(26) $\displaystyle|R_{3}^{(\pm)}(t_{3})\rangle$ $\displaystyle=$
$\displaystyle e^{\pm t_{3}R_{3}}|\Phi\rangle=e^{\pm
t_{3}\\{(a_{6}^{\dagger}a_{2}+a_{8}^{\dagger}a_{4})-(a_{2}^{\dagger}a_{6}+a_{4}^{\dagger}a_{8})\\}}|\Phi\rangle\;,$
(27) $\displaystyle|R_{4}^{(\pm)}(t_{4})\rangle$ $\displaystyle=$
$\displaystyle e^{\pm t_{4}R_{4}}|\Phi\rangle=e^{\pm
t_{4}\\{(a_{5}^{\dagger}a_{1}+a_{7}^{\dagger}a_{3})-(a_{1}^{\dagger}a_{5}+a_{3}^{\dagger}a_{7})\\}}|\Phi\rangle\;,$
(28) $\displaystyle|R_{5}(t_{5})\rangle$ $\displaystyle=$ $\displaystyle
e^{t_{5}R_{3}}e^{t_{5}R_{4}}|\Phi\rangle\;,$ (29)
$\displaystyle|R_{6}(t_{6})\rangle$ $\displaystyle=$ $\displaystyle
e^{t_{6}R_{4}}e^{t_{6}R_{3}}|\Phi\rangle\;,$ (30)
$\displaystyle|R_{2}^{(\pm)}R_{1}^{(\pm)}(t_{7})\rangle$ $\displaystyle=$
$\displaystyle e^{\pm t_{7}R_{2}}e^{\pm t_{7}R_{1}}|\Phi\rangle=e^{\pm
t_{7}\\{(a_{6}^{\dagger}a_{1}+a_{8}^{\dagger}a_{3})-(a_{1}^{\dagger}a_{6}+a_{3}^{\dagger}a_{8})\\}}e^{\pm
t_{7}\\{(a_{5}^{\dagger}a_{2}+a_{7}^{\dagger}a_{4})-(a_{2}^{\dagger}a_{5}+a_{4}^{\dagger}a_{7})\\}}|\Phi\rangle\;.$
(31)

These sampling vectors can be naturally tied to the general form of the
$|\Phi^{(2)}({\bf Z}(1),{\bf Z}(2))\rangle$ basis given by Eq.(23). For
example, $|R_{0}\rangle$ corresponds to ${\bf Z}(1)={\bf Z}(2)=0$. For
$|R_{1}^{(+)}(t_{1})\rangle$, ${\bf Z}(2)=0$, and
$z(1)_{52}=z(1)_{74}=-z(1)_{25}=-z(1)_{47}=t_{1}$ while remaining matrix
elements are equal to zero, etc. Another advantage of using this form of
sampling vectors is that their combinations provide a rudimentary (yet not
exact) mechanism for eliminating symmetry impurities when $R_{i}$ operators
break the symmetry of the reference state $|\Phi\rangle$. For example, if in
general $R_{1}^{(\pm)}(t_{1})$ the operator is expressed in terms of
excitations that produce a triplet state when acting on the reference
function, then this singly-excited impurity (or instability) is eliminated by
taking a combination of $|R_{1}^{(+)}\rangle+|R_{1}^{(-)}\rangle$ states; that
is,

$\displaystyle|R_{1}^{(+)}(t_{1})\rangle+|R_{1}^{(-)}(t_{1})\rangle=(e^{t_{1}R_{1}}+e^{-t_{1}R_{1}})|\Phi\rangle$
$\displaystyle=(2+t_{1}^{2}R_{1}^{2}+\frac{2}{4!}t_{1}^{4}R_{1}^{4}+\ldots)|\Phi\rangle\;,$
(32)

where linear triplet “impurities” are eliminated. These combinations also can
allow us to selectively approach doubly excited configurations, or using
combinations of $|R_{2}^{(\pm)}R_{1}^{(\pm)}(t_{7})\rangle$, quadruply excited
ones using manifold of single excitations and very simple quantum circuits to
represent $e^{\pm tR_{i}}$ operators. This analysis can be extended to higher-
order PGCM(k) formulations to include higher-rank excitations. It should be
stressed that $t_{i}$ parameters in Eqs.(24)-(31) can be chosen at random. The
effects of random $t_{i}$ parameters on ground-state FCI energies are
illustrated in Appendix A.

##  V General outline of quantum algorithm and post-processing on conventional
computers

After selecting GCM sampling points, the remaining work is involves computing
matrices ${\bf H}$ and ${\bf S}$ through Eq.(14) and (15) using quantum
computers. To accommodate gate-based quantum computers, the expectations are
computed in the form of

$\displaystyle{\bf H}_{pq}$
$\displaystyle=\sum_{j}h_{j}\langle\Phi|\left(e^{-\Gamma({\bf
Z}_{p})}P_{j}e^{\Gamma({\bf Z}_{q})}\right)|\Phi\rangle$ (33)
$\displaystyle\text{and }{\bf S}_{pq}$
$\displaystyle=\langle\Phi|\left(e^{-\Gamma({\bf Z}_{p})}e^{\Gamma({\bf
Z}_{q})}\right)|\Phi\rangle,$ (34)

where Hamiltonian matrix $H$ is transformed to the linear combination of Pauli
strings $H=\sum_{j}h_{j}P_{j}$ under Jordan-Wigner (JW) transformation and
$|\Phi\rangle$ is the HF state. The expectations can be easily evaluated using
algorithms like the Hadamard test on fault-tolerant quantum computers. A pure
quantum algorithm proposed for GCM and its complexity analysis is given in
Appendix B, where we propose to exploit a block encoding for the operation
$\mathbf{S}^{-1}$ and use the phase estimation to compute the eigenvalues for
the non-orthogonal eigenproblem (13). Note that the pure quantum algorithm
gives a favorable scaling but relies on the use of the controlled unitary
circuits that makes the quantum simulations on NISQ devices challenging.

For near-term devices, we focus on hybrid quantum-classical approach. there is
a trade-off between the depth of the quantum circuit and classical computation
time. Note that, following the JW transformation, each of the operators
$R_{1}$ to $R_{4}$ in Eq.(25) to (28) can be transformed to four commuting
Pauli strings. That is to say, the matrix exponential $e^{\pm t_{i}R_{i}}$ for
$i\in\\{1,2,3,4\\}$ are exactly the linear combinations of Pauli strings. So,
for both H4 models, the calculations of expectation in Eq.(33) and (34) are
essentially equivalent to $\langle\Phi|P|\Phi\rangle$ for some Pauli string
$P$, while qubit-wise commuting (QWC) terms are grouped to perform a
simultaneous measurement. The only implemented form of the circuit is
illustrated in Fig.2 and the whole process is summarized in Alg.1. The entire
simulated quantum computation is achieved using Qiskit.tA v  _et al._ (2021)
We put a brief discussion about trotterization under Qiskit in Appendix C,
along with other implementation details related to duplicated Pauli strings
and effects of random parameters in Appendix A.

Algorithm 1 Quantum GCM (QuGCM) for near-term devices

1:Hamiltonian matrix $H=\sum_{j}h_{j}P_{j}$, HF state $|\Phi\rangle$, and a
set $\\{\Gamma({\bf Z}_{i})\\}_{i=1}^{M}$ where the index $i$ could be a
composite up to $k$ terms as in Eq.(23)

2:Transform all $\\{\Gamma({\bf Z}_{i})\\}_{i=1}^{M}$ using JW transformation

3:Generate unitaries $\\{V_{i}\\}_{i=1}^{M}$ for $V_{i}:=e^{\Gamma({\bf
Z}_{i})}$ with Eq.(9) and Eq.(20)

4:Trotterize each element in $\\{V_{i}\\}_{i=1}^{M}$ to a linear combination
of Pauli strings

5:for each $V_{q}$ in $\\{V_{i}\\}_{i=1}^{M}$ do

6: Compute $\sum_{j}h_{j}P_{j}V_{q}$ classically

7: for each $V_{p}$ in $\\{V_{i}\\}_{i=1}^{M}$ do

8: Compute $\sum_{j}h_{j}V_{p}^{\dagger}(P_{j}V_{q})$ classically

9: Compute $V_{p}^{\dagger}V_{q}$ classically

10: Evaluate ${\bf
H}_{pq}:=\sum_{j}h_{j}\langle\Phi|V_{p}^{\dagger}P_{j}V_{q}|\Phi\rangle$ in a
quantum device

11: Evaluate ${\bf S}_{pq}:=\langle\Phi|V_{p}^{\dagger}V_{q}|\Phi\rangle$ in a
quantum device

12: end for

13:end for

14:Solve the general eigenvalue problem ${\bf H}{\bf f}=E{\bf S}{\bf f}$
classically

15:return interested eigenvalues and eigenvectors

Figure 2: Shallow circuit for Alg.1. State-preparation gates initialize the
ground state to the HF state. Pauli gates are a single layer of 1-qubit Pauli
gates. QWC rotation gates consist of 1-qubit rotation gates that allow us to
measure each wire on Pauli-$X,Y,$ or $Z$ bases with the final Pauli-$Z$-basis
measurement. So all gates in the circuit are 1-qubit gates and the total
number of 1-qubit gates is $O(n)$.

##  VI Complexity Analysis and effects of finite samplings

Recall that $M$ is the number of selected ${\bf Z}$-points, which is
equivalent to the number of sampling vectors. Let $n$ be the number of spin
orbitals. We know the Hamiltonian matrix, $H$, consist of at most $O(n^{4})$
Pauli strings, so operations for multiplying matrices and solving the
generalized eigenvalue problem in Alg.1 dominate the classical part of the
operation complexity. If an exponent of a cluster operator has up to $c$
number of Pauli strings after trotterization and the multiplication between
two length-$n$ Pauli strings takes at most $O(n)$ operations, then the matrix
multiplication in line 10 and 11 of Alg.1 takes $O(n\cdot c^{2}\cdot n^{4})$
and $O(n\cdot c^{2})$ operations, respectively, because all matrices are
decomposed into the linear combinations of Pauli terms. As we have $M^{2}$
iterations in Alg.1 and the final generalized eigenvalue problem of $M\times
M$ matrices takes $O(M^{3})$ operations, Hybrid PGCM(2) has the overall
classical part of the worst-case scaling of $O(c^{2}n^{5}M^{2}+M^{3})$
operations.

The specific value of $c$ highly depends on the how users want to approximate
the molecular models. It is affected by the number of creation and
annihilation operators pairs in each single-excitation cluster operators, the
value of $k$, and number of trotterization steps if it is necessary. As we
discussed earlier, we can properly choose two pairs of creation and
annihilation operators in each of single-excitation cluster operators without
the requirements of trotterization. The results in Section VII are promising
enough when only single and double excitations are considered, i.e., when
$k=2$.

Regarding the quantum part of the operation complexity, it includes the number
of circuits and the number of measurements in every circuit. Naively, line 10
and line 11 in Alg.1 indicate there are $O(n\cdot c^{2}\cdot n^{4})$ number of
circuits and every circuit contains only $O(n)$ number of 1-qubit gates. Many
existing methods can be applied to reduce the number of terms associated with
$V_{p}^{\dagger}(HV_{q})$ and $V_{p}^{\dagger}V_{q}$. For example, to reduce
the order $O(n\cdot c^{2}\cdot n^{4})$ to $O(n\cdot c^{2}\cdot n^{2\sim 3})$
for the number of circuits in every iteration, the linear combination of
unitaries technique,Childs and Wiebe (2012) amplitude amplification
approach,Brassard  _et al._ (2002) Hamiltonian simulation,Berry  _et al._
(2015) qubitization,Low and Chuang (2017) or the direct block-encodingGilyén
_et al._ (2019) methods can be typically employed at the cost of introducing
deeper circuits and implementing controlled unitary operations. Nevertheless,
these approaches come with a probability of failure and require advanced
circuit and error mitigation that might go beyond the capability of the
current NISQ devices. Toward a more feasible NISQ approach, in light of the
unitary partitioning scheme proposed by Izmaylov et al.Izmaylov  _et al._
(2020), Peng and Kowalski recently proposed a more efficient unitary
partitioning approach guided by the single-reference trial state used in the
simulation.Peng and Kowalski (2022) In particular, through numerical tests
over a wide range of molecules in different bases, they found that the non-
unitary wave operators when acting on single-reference trial wave function
(such as $e^{\Gamma(\mathbf{Z}_{i})}|\Phi\rangle$ and
$He^{\Gamma(\mathbf{Z}_{i})}|\Phi\rangle$ in the present discussion) can be
represented by a much more compact unitary basis, thus providing a more
efficient route for performing the general non-unitary quantum simulations.

It is worth mentioning that the above discussion is focused on the number of
terms/operations that can be efficiently reduced through groupings that
feature the commutativity or anti-commutativity of Pauli strings, while the
total number of measurements required from the number of groups also
critically depends on the covariance between the contributing terms,
$\text{cov}\big{(}P_{i},P_{j}\big{)}$, and the desired precision $\epsilon$.
Given that we can always write a matrix operator as a linear combination of
Pauli strings, the total number of measurements for evaluating the expectation
value of an operator with respect to the trial wave function can then be
expressed as Gonthier  _et al._ (2020); Rubin, Babbush, and J. (2018); Wecker,
Hastings, and Troyer (2015b)

$\displaystyle\text{\\# of Measurements}=$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \Bigg{(}\displaystyle\frac{\sum_{G}\sqrt{\sum_{i,j,\in
G}h_{i}h_{j}\text{cov}\big{(}P_{i},P_{j}\big{)}}}{\epsilon}\Bigg{)}^{2}$ (35)

with $G$ indexing the groups. Therefore, it is likely that the number of
groups decreases at the cost of introducing larger covariances that could
essentially increase the total number of measurements required to achieve a
desired precision. The variance of an individual Pauli string $P_{i}$ can be
bounded by

$\displaystyle\text{var}(P_{i})=\braket{P_{i}^{2}}-\braket{P_{i}}^{2}=1-\braket{P_{i}}^{2}\leq
1,$ (36)

which provides bounds to the covariance of any contributing terms

$\displaystyle\left|\text{cov}(P_{i},P_{j})\right|\leq\text{var}(P_{i})\text{var}(P_{j})\leq
1.$ (37)

Thus, according to Eq.(35), the bounds for the standard deviations of the
matrix entry ${\bf H}_{pq}$ and ${\bf S}_{pq}$ under a finite number of
measurements are

$\displaystyle\epsilon_{{\bf H}_{pq}}$ $\displaystyle\leq\frac{\sum_{G_{{\bf
H}_{pq}}}\sqrt{\sum_{i,j,\in G_{{\bf H}_{pq}}}|h_{i}h_{j}|}}{\sqrt{\text{\\#
of Measurements}}}$ (38) $\displaystyle\epsilon_{{\bf S}_{pq}}$
$\displaystyle\leq\frac{\sum_{G_{{\bf S}_{pq}}}\sqrt{\sum_{i,j,\in G_{{\bf
S}_{pq}}}|h_{i}h_{j}|}}{\sqrt{\text{\\# of Measurements}}}$ (39)

Let $\tilde{{\bf H}}_{pq}$ and $\tilde{{\bf S}}_{pq}$ be the entries computed
from finite number of measurements. We can empirically estimate the effects of
the uncertainty by considering $\tilde{{\bf H}}_{pq}$ and $\tilde{{\bf
S}}_{pq}$ as normal random variables

$\displaystyle\tilde{{\bf H}}_{pq}$ $\displaystyle\sim N\left({\bf
H}_{pq},\epsilon^{2}_{{\bf H}_{pq}}\right)\text{ and }\tilde{{\bf S}}_{pq}$
$\displaystyle\sim N\left({\bf S}_{pq},\epsilon^{2}_{{\bf S}_{pq}}\right).$
(40)

To illustrate the maximum possible fluctuations brought by finite-sampling
errors, we set each $\epsilon_{{\bf H}_{pq}}$ and $\epsilon_{{\bf S}_{pq}}$ to
their maximum as in Eq.(38) and Eq.(39), and summarize the results in Fig.3.
It is clear that to generally reach chemical accuracy, the $\alpha=0.005$ case
requires about two orders of magnitude more measurements than the
$\alpha=0.500$ case using Alg.1.

Figure 3: Differences between the ground-state energy estimations from PGCM(2)
and FCI formalism in milli-Hartrees for $\alpha=0.005$ (up) and $\alpha=0.500$
(down) H4 models, respectively. The variance of each random sampling is set to
the maximum according to Eq.(38) and Eq.(39). The red-shaded region is the
range of chemical accuracy $\pm 1.5936\text{mHa}$. The blue-shaded region is
the $95\%$ confidence interval estimated from 100 simulations with the red
bars marking the rough number of measurement per entry at which the energy
difference would be within the chemical accuracy.

We also notice that there are also many other advanced measurement schemes
proposed recently. One example is to simultaneously obtain expectation values
of multiple observations by randomly measuring and projecting the quantum
state into classical shadows Huang, Kueng, and Preskill (2020, 2021); Aaronson
(2018); Struchalin  _et al._ (2021); Chen  _et al._ (2021); Zhao, Rubin, and
Miyake (2021); Acharya, Saha, and Sengupta (2021); Hadfield (2021); Hillmich
_et al._ (2021); Zhang  _et al._ (2021). In principle, the algorithm enables
measurements of $m$ low-weight observations using only $O(\log_{2}m)$ samples.
The practical performance of the algorithm for model and molecular
Hamiltonians on NISQ devices, in terms of accuracy and efficiency, is still
under intense study.

##  VII Results

The GCM results for the ground-state energies and excitation energies
corresponding to low-lying states of the symmetry of the reference function
(singlet $A_{1}$ states) are collated in Tables 1 and 2, respectively. To
evaluate the accuracy of ground-state simulations, we compared GCM results
with the results obtained with the RHF, multi-configurational self-consistent
field (MCSCF) formalism for active space defined by four electrons and three
active orbitals (MCSCF(4e,3o)), configuration interaction method with singles
and doubles (CISD), CC method with singles and doubles (CCSD), VQE formalism,
and FCI formalism. We used the equation-of-motion CC approach (EOMCCSD) and
FCI formalism for excited states. Note that any quantum computation in VQE and
GCM assume infinite number of measurements.

Inspection of the results in Table 1 indicates that for both geometries, the
GCM results are in very good agreement with the FCI result despite the
simplicity of the expansions given by Eqs.(24)-(31). Interestingly, the GCM
energies are significantly better in both cases than the VQE ones. The
excellent performance of the GCM formalism is well illustrated by the strongly
correlated $\alpha=0.005$ case, where CISD and CCSD methods struggle to
capture needed correlation effects. For the weakly correlated variant of H4
($\alpha=0.500$), the GCM formalism reproduces nearly FCI-level accuracy with
the energy error of 0.022 milli-Hartree.

In Table 2 we juxtaposed the excitation energies ($\omega_{1}$, $\omega_{2}$,
and $\omega_{3}$) obtained with the GCM approaches for three lowest-lying
${}^{1}A_{1}$ symmetry states, of H4 model for $\alpha=0.005$ and
$\alpha=0.500$. For the $\omega_{1}$ excitation energies, the GCM approach
provides consistently better estimates of their exact (FCI) values than the
ubiquitous EOMCCSD approach. While for $\omega_{2}$ GCM prediction is better
than the EOMCCSD one only for the $\alpha=0.005$, for non-degenerate case
($\alpha=0.500$), the GCM prediction is by 0.6 eV off the FCI value. For the
$\omega_{3}$ excitation energies, the GCM is capable of providing estimates
within 0.160 eV ($\alpha=0.005$) and 0.33 eV ($\alpha=0.500$) of error. Again,
given the simplicity of the GCM formulations, one should view the GCM
estimates of the excitation energies as quite satisfactory.

FCI | RHF | MCSCF(4e,3o) | CISD | CCSD | VQE | GCM
---|---|---|---|---|---|---
H4 $\alpha=0.005$
-1.942993 | 151.407 | 82.645 | 5.507 | 3.331 | 0.905 | 0.147
H4 $\alpha=0.500$
-2.151007 | 75.764 | 43.276 | 1.866 | 0.003 | 0.058 | 0.022
Table 1: Differences (in milliHartrees) with total ground-state FCI energies. Method | $\omega_{1}$ | $\omega_{2}$ | $\omega_{3}$
---|---|---|---
H4 $\alpha=0.005$
FCI | 4.183 | 6.040 | 18.484
EOMCCSD | 4.275 | 6.079 | 18.350
GCM | 4.179 | 6.038 | 18.635
H4 $\alpha=0.500$
FCI | 12.565 | 14.214 | 21.293
EOMCCSD | 12.738 | 14.179 | 21.255
GCM | 12.589 | 14.840 | 21.622
Table 2: Singlet excitation energies in eV

##  VIII Conclusion

In this work, we explored use of the GCM in the context of quantum computing.
For this purpose, we introduced the multi-product extension of the GCM
formalism that enables one to construct state vectors in Hilbert space using
various types of the Fermion Lie algebra and general quantum algorithms that
allow one to perform GCM calculations in a way that can be viewed as a
specific case of the quantum algorithms for configuration interaction
formalisms involving non-orthogonal basis. In the present studies, we focused
entirely on the $U(N)$ algebra, where resulting state vectors can be
interpreted in terms of the Thouless theorem. This analogy is essential in the
sampling process of the parametrized unitary canonical transformations. It
enables one to construct corresponding state vectors and corresponding linear
space where higher-order excitations (e.g., double, triple, quadruple, etc.)
excitations can be selectively approached using the language of single
excitations to define generators ($\Gamma({\bf Z})$) of the canonical
transformations. The discussed procedures can be easily related to the
searches of various type instabilities in independent particle formulations
with the HF method as a specific example.

Using the H4 system as a benchmark, we showed that the quantum GCM algorithm
could provide ground-state energies competitive to the VQE simulations
involving the explicit form of the double excitations. In contrast to standard
VQE algorithms, the GCM formalism also yields the values of excited-state
energies. We showed that GCM excitation energies corresponding to the low-
lying excited states could be competitive with the excitation energies
obtained with the popular EOMCCSD approach.

Although in the present form, the GCM formalism falls into the category of
hybrid formulations, unlike the VQE method, it avoids multiple quantum-
classical machine data transfer and unstable and tedious iterative processes.
Also, our near-term implementation of the quantum GCM algorithm can run in
parallel easily using multiple quantum and classical nodes.

In future studies, we plan to use the full potential of the quantum algorithms
based on the quantum CI formulations and explore the possibility of using
other types of Fermionic algebras, i.e., $SO(2N)$ and $SO(2N+1)$ Lie algebras.

##  IX ACKNOWLEDGMENTS

Muqing Zheng, Bo Peng, Nathan Wiebe, Ang Li, and Karol Kowalski were supported
by the “Embedding QC into Many-body Frameworks for Strongly Correlated
Molecular and Materials Systems” project, which is funded by the U.S.
Department of Energy, Office of Science, Office of Basic Energy Sciences, the
Division of Chemical Sciences, Geosciences, and Biosciences. The Pacific
Northwest National Laboratory is operated by Battelle for the U.S. Department
of Energy under Contract DE-AC05-76RL01830. Xiu Yang was supported by National
Science Foundation CAREER DMS-2143915. Muqing Zheng and Xiu Yang both also
were supported by Defense Advanced Research Projects Agency as part of the
project W911NF2010022: The Quantum Computing Revolution and Optimization:
Challenges and Opportunities.

##  X Data availability

The code and data in the paper are openly available in the GitHub
repository.Zheng  _et al._ (2022)

##  Appendix A Duplicated Pauli strings and selection of parameters

For each of H4 models, while $M=15$, the symmetries of ${\bf H}$ and ${\bf S}$
allows us to only compute 120 iterations instead of 225 iterations. Among
those 120 iterations, as shown in Table 3, we only need to measure 4,000 Pauli
strings because more than $97\%$ of them are duplicated. This brings the two
orders of magnitude reduction on the total number of measurements in practice.

| Num. of unique terms | Total num. of terms | Ratio ($\%$)
---|---|---|---
H4 $\alpha=0.005$
${\bf S}$ | $1444$ | $38904$ | $3.71\%$
${\bf H}$ | $3423$ | $118426$ | $2.89\%$
All | $4180$ | $157330$ | $2.66\%$
H4 $\alpha=0.500$
${\bf S}$ | $1766$ | $38888$ | $4.54\%$
${\bf H}$ | $3454$ | $118621$ | $2.91\%$
All | $4370$ | $157509$ | $2.77\%$
Table 3: Number of unique Pauli terms and their ratios after grouping among
all iterations.

Meanwhile, we conducted the following experiments for $\alpha=0.005$ and
$\alpha=0.500$ H4 models to demonstrate the influence of the random choices of
$t_{i}$ parameters on the estimations of the ground-state FCI energies. For
each molecular model, we generated 50 sets of $\\{t_{i}\\}_{i=1}^{7}$ for
$t_{i}$ in $[0,1)$, $[0,100)$, and $[100,100)$, respectively. Each of 300 sets
of $\\{t_{i}\\}_{i=1}^{7}$ produced an estimation of the ground-state FCI
energy from Alg. 1 (assume infinite number of shots). The distributions of
differences between the estimations from PGCM(2) and FCI formalism under
various settings are illustrated in Fig.4 using box-plots. It worth noting
that when we produce $t_{i}$’s from $[0,1)$, about $25\%$ random generations
can provide an estimation of ground-state energy within the range of chemical
accuracy.

Figure 4: Differences between the ground-state energy estimations from PGCM(2)
in various random $\\{t_{i}\\}$’s and FCI formalism in milli-Hartrees for
$\alpha=0.005$ (up) and $\alpha=0.500$ (down) H4 models, respectively. The
horizontal red-shaded region shows the range of chemical accuracy $\pm
1.5936\text{ mHa}$. The fact that minimum value and the 25th percentile are in
the shaded region when $t_{i}\in[0,1)$ in both H4 models indicates that random
generations of ${t_{i}}$ in that range is relatively likely to provide a good
approximation in our algorithm.

##  Appendix B Linear-System Inspired Algorithms for GCM

Perhaps the most direct way to use quantum computing to solve the electronic
structure problem using Generator Coordinate Methods is by simply solving the
non-orthogonal eigenvalue problem by dilating it to a square matrix in a
higher dimensional space. In this case, it is most convenient to express our
eigenvalue problem as

${\bf S}^{-1}{\bf H}{\bf f}=E{\bf f},$ (41)

where we have assumed here that ${\bf S}$ is an invertible matrix. Now let us
define an isometric extension of our original space. We do this to exploit a
block encoding for the operation ${\bf S}^{-1}$, which allows us to express it
as a unitary operation in a higher dimensional space. Specifically, let

$\displaystyle{\bf U}$ $\displaystyle=\begin{bmatrix}{\bf S}^{-1}/\|{\bf
S}^{-1}\|&{\Box}\\\ {\Box}&{\Box}\end{bmatrix},$ (42)

be a unitary matrix (i.e., ${\bf U}^{\dagger}={\bf U}^{-1}$) for arbitrary
matrices ${\Box}$. Also let

$\displaystyle{\bf J}$ $\displaystyle=\begin{bmatrix}{\bf H}/\alpha&0\\\
0&0\end{bmatrix},\leavevmode\nobreak\ \leavevmode\nobreak\ {\bf
g}=\begin{bmatrix}{\bf f}\\\ 0\end{bmatrix},\leavevmode\nobreak\
\leavevmode\nobreak\ \bf{Z}=\begin{bmatrix}\bf{I}&0\\\
0&-\bf{I}\end{bmatrix},\leavevmode\nobreak\ \leavevmode\nobreak\ {\bf
P}=\frac{1}{2}\left(\bf{I}+\bf{Z}\right).$ (43)

Inside this enlarged space, the eigenvalue equation reads for eigenvalue $E$

${\bf P}{\bf U}{\bf J}{\bf P}{\bf g}=\frac{E\bf{g}}{\|{\bf S}^{-1}\|\alpha}.$
(44)

Now let ${\bf H}/\alpha=\sum_{j}h_{j}{\bf U}^{(H)}_{j}$ for a set of unitary
${\bf U}_{j}^{(H)}$. We then have that our Hamiltonian in the enlarged space
can be expressed as a similar linear combination of unitaries.

${\bf J}=\sum_{j}h_{j}|0\rangle\\!\langle 0|\otimes{\bf
U}^{(H)}_{j}=\sum_{j}\frac{h_{j}}{2}\left({\bf I}\otimes{\bf
U}^{(H)}_{j}+\mathbf{Z}\otimes{\bf U}^{(H)}_{j}\right),$ (45)

Thus the entire product can be expressed as

${\bf P}{\bf U}{\bf J}{\bf P}=\sum_{j}\frac{h_{j}}{2}{\bf
P}\left({\bf{U}}\left({\bf I}\otimes{\bf
U}^{(H)}_{j}\right)+{\bf{U}}\left(\mathbf{Z}\otimes{\bf
U}^{(H)}_{j}\right)\right){\bf P},$ (46)

Thus the enlarged Hamiltonian ${\bf J}$ can be expressed as a linear
combination of unitary matrices that further has the exact same value of
$\alpha=\sum_{j}|h_{j}|$.

For us to use phase estimation to compute these eigenvalues, we also need a
further modification. Note that $e^{-i{\bf UJ}}$ is not necessarily unitary
because ${\bf UJ}$ is not necessarily Hermitian. We address this issue by
considering a further embedding:

$\displaystyle{\bf K}$ $\displaystyle=\begin{bmatrix}0&{\bf PUJP}\\\ {\bf
PJ}{\bf U}^{\dagger}{\bf P}&0\end{bmatrix}$ $\displaystyle=|0\rangle\\!\langle
1|\otimes{\bf PUJP}+|1\rangle\\!\langle 0|\otimes{\bf PJ}{\bf U}^{\dagger}{\bf
P}$ $\displaystyle=\frac{1}{2}\left({\bf X}\otimes{\bf PUJP}+(i{\bf
Y})\otimes{\bf PUJP}+{\bf X}\otimes{\bf PJ}{\bf U}^{\dagger}{\bf P}\right.$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \left.+(-i{\bf
Y})\otimes{\bf PJ}{\bf U}^{\dagger}{\bf P}\right)$ (47)

This expression also takes the form of a linear combination of unitary
matrices; however, the coefficient sum now obeys
$\sum_{j}|h_{j}^{\prime}|=2\alpha$, which does not change the overall
complexity. Finally note that ${\bf K}$ is anti-diagonal block matrix. We can
therefore search for an eigenvector of the form ${\bf h}=\begin{bmatrix}0&{\bf
g}\end{bmatrix}^{\dagger}$ as the eigenvectors of ${\bf K}$ can be taken to
have support only on one of the two blocks. The eigenvalues of $\mathbf{K}$
corresponding the the eigenvector $\mathbf{g}$ of ${\bf PUJP}$ must be $\pm E$
from this construction. To see this note that

${\rm Tr}(\mathbf{K}^{2}\mathbf{h}\mathbf{h}^{\dagger})={\rm
Tr}(\mathbf{h^{\dagger}K^{\mathrm{2}}h})=(\mathbf{Kh})^{\dagger}(\mathbf{Kh})=E^{2}.$
(48)

Thus $\mathbf{h}$ is an eigenvector of $\mathbf{K}$ and the eigenvalues of
$\mathbf{K}$ in the support of $\mathbf{h}$ must be $\pm E$. Thus, we can
estimate the absolute value of $E$ by using phase estimation on ${\bf K}$.

To perform this simulation using qubitization ideas, we will need to propose
prepare-and-select circuits for the coefficients. Specifically,

$U_{\rm
prep}\ket{0}=\sum_{j}\sum_{\mu=0,1}\sum_{\nu=0,1}\sum_{\omega=0,1}\frac{\sqrt{h_{j}/4}}{\sqrt{2\alpha}}\ket{j}\ket{\mu}\ket{\nu}\ket{\omega}$
(49)

and further:

$\displaystyle U_{\rm sel}\ket{j}\ket{\mu}\ket{\nu}\ket{\omega}\ket{\psi}$
$\displaystyle=\ket{j}\ket{\mu}\ket{\nu}\ket{\omega}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\otimes(-1)^{\mu\omega}{\bf X}^{\mu}(i{\bf Y})^{1-\mu}\otimes{\bf Z}^{\nu}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \otimes{{\bf
U}^{\dagger}}^{\omega}{\bf U}_{j}^{(H)}{\bf U}^{1-\omega}\ket{\psi}$ (50)

This then forms a block encoding of our operator for the generalized
eigenvalue problem:

$(\bra{0}\otimes\mathbf{I})U_{\rm prep}^{\dagger}U_{\rm sel}U_{\rm
prep}(\ket{0}\otimes\mathbf{I})={\bf K}/(2\alpha)$ (51)

Using qubitization, we can convert this into a unitary with result Gilyén  _et
al._ (2019), we can perform phase estimation to learn $E$ with a variance of
at most $\epsilon^{2}$ using a number of queries to $U_{\rm sel}$ and $U_{\rm
prep}$ that is in $O\left(\alpha\|{\bf S}^{-1}\|/\epsilon\right)$. However,
each application of $U_{\rm sel}$ requires $O(1)$ queries to ${\bf U}$, which
is a block encoding of the inverse of ${\bf S}$. A deeper analysis of the cost
is possible; however, to do this, we first need to have consensus on the input
model used for the simulation.

The easiest way to compute a matrix element for the overlap matrix is through
the use of a controlled unitary. The overlap that, the estimate takes the
form:

$\bra{\Phi}V^{\dagger}_{p}V_{q}\ket{\Phi}={\bf S}_{p,q}.$ (52)

where $V_{p}$ is the basis transform operation such that for the reference
state $\ket{\Phi}$, $V_{p}\ket{\Phi}=\ket{\Phi({\bf Z}_{p})}$. To use the
Hadamard test to reconstruct this circuit we require a single application of
$V^{\dagger}_{p}$ and a single application of $V_{q}$. The probability that
the control qubit that governs this circuit is $0$ is::

$P(0|p,q)=\frac{1+{\rm Re}(\langle{\Phi({\bf Z}_{p})}|{\Phi({\bf
Z}_{q})}\rangle)}{2}.$ (53)

If needed, the imaginary part can be similarly found by applying an
$S^{\dagger}$ gate to the control. If we apply amplitude amplification to the
result, then we can construct a matrix $\mathbf{AA(P)}$ with eigenvalues
inside the sector:

$\lambda(\mathbf{AA(P)})=e^{\pm i\cos^{-1}\left(\sqrt{\frac{1+{\rm
Re}(\langle{\Phi(Z_{p})}|{\Phi(Z_{q})}\rangle)}{2}}\right)},$ (54)

and is equivalent to the following matrix in the two dimensional space spanned
by the initial state and the marked state:

$\begin{bmatrix}\sqrt{\frac{1+{\rm
Re}(\langle{\Phi(Z_{p})}|{\Phi(Z_{q})}\rangle)}{2}}&-\sqrt{\frac{1-{\rm
Re}(\langle{\Phi(Z_{p})}|{\Phi(Z_{q})}\rangle)}{2}}\\\ \sqrt{\frac{1-{\rm
Re}(\langle{\Phi(Z_{p})}|{\Phi(Z_{q})}\rangle)}{2}}&\sqrt{\frac{1+{\rm
Re}(\langle{\Phi(Z_{p})}|{\Phi(Z_{q})}\rangle)}{2}}\end{bmatrix}$ (55)

Applying quantum signal processing we can apply a transformation $u\mapsto
2u^{2}-1$. Note that this is 1) an even degree polynomial and 2 it is between
$[-1,1]$ for $u$ in a similar range. This means that quantum signal processing
can be used to apply this transformation on the top block of the matrix to
prepare a block encoding of the form Gilyén  _et al._ (2019):

$U_{\bf S}=\begin{bmatrix}{{{\bf
Re}(\langle{\Phi(Z_{p})}|{\Phi(Z_{q})}\rangle)}}&\Box\\\
\Box&\Box\end{bmatrix}$ (56)

Where specifically we have that, the initial state used in the amplitude
amplification routine block-encodes the overlap matrix $\mathbf{S}$. In all,
this process costs $O({\rm polylog}(1/\epsilon))$ queries to the state
preparation oracle to produce this block encoding.

Next we need to prepare a block encoding of ${\mathbf{S}}^{-1}$. Using the
results of Ref. 75 we can prepare such a block encoding using $O(\|{\bf
S}\|\|{\bf S}^{-1}\|{\rm polylog}(1/\epsilon))$ queries where $\epsilon$ is
our target accuracy. Thus, the overall query complexity is the number of
queries made to $\mathbf{K}$ multiplied by the number of queries per $U$ or
$U_{j}^{(H)}$. The former result then shows that the final cost of the
simulation using this approach (in terms of queries to the prepare and select
oracles of ${\mathbf{H}}$ and the queries to $Z_{p},Z_{q}$ scale as:

$\tilde{O}\left(\frac{\alpha\|{\bf{S}}\|\|{\bf
S}^{-1}\|^{2}}{\epsilon}\right).$ (57)

This shows that the query complexity of the simulation using this approach
does not necessarily scale with the dimension of the space. It does, however,
depend strongly on the one-norm of the coefficients of the Hamiltonian and the
norm of the inverse of the overlap matrix. Thus, the worse conditioned the
matrix is, the worse we expect the performance of the algorithm to be. In
contrast, learning the matrix $\mathbf{S}$ in high dimensional spaces to
perform the inverse can be expensive as noted in the main body. This approach
gives a theoretical alternative in such cases that has favorable scaling
asymptotically at the price of the algorithm requiring a large number of
qubits.

##  Appendix C Trotterization of a matrix exponential

This section aims to transform a matrix exponential into a linear combination
of Pauli group matrices. Because the standard form of a matrix exponential in
Qiskit is $e^{-itA}$ instead of $e^{tA}$ as we have in the main body for some
parameter $t$ and matrix $A$ and we implemented the quantum part of our
algorithm in Qiskit, we keep $i=\sqrt{-1}$ in the power. Because we will use
dot product between matrices, we do not omit the operator $\otimes$ when we do
the Kronecker product between Pauli matrices for clarity (e.g., $XY$ means the
dot product of $X$ and $Y$ instead of $X\otimes Y$). Let $n$ be the number of
qubits, $I_{n}$ the identity matrix in space
$\mathbb{C}^{2^{n}}\times\mathbb{C}^{2^{n}}$, and $P\in\\{X,Y,Z,I\\}^{\otimes
n}$ an $n$-qubit Pauli group matrix. Then, for any $n$-qubit Pauli group
matrix, we have:

$P^{2}=\left(\bigotimes_{l=1}^{2^{n}}\sigma_{l}\right)\left(\bigotimes_{l=1}^{2^{n}}\sigma_{l}\right)=\left(\bigotimes_{l=1}^{2^{n}}\sigma_{l}^{2}\right)=I_{n}$
(58)

given $\sigma_{l}^{2}=I$ , for all $\sigma_{l}\in\\{X,Y,Z,I\\}$ and the
property of Kronecker product

$(A\otimes B)(C\otimes D)=(AC)\otimes(BD)$ (59)

for some matrices $A,B,C,D$ in the appropriate dimensions.

Now, by doing Taylor expansion of $e^{itP}$ at $t=0$, we obtain the following
exact conversion:

$\displaystyle e^{itP}$ $\displaystyle=$ $\displaystyle
I_{n}+(itP)+\frac{(itP)^{2}}{2!}+\frac{(itP)^{3}}{3!}+\frac{(itP)^{4}}{4!}+\frac{(itP)^{5}}{5!}\cdots$
(60) $\displaystyle=$ $\displaystyle I_{n}+itP-
I_{n}\frac{t^{2}}{2}-iP\frac{t^{3}}{3!}+I_{n}\frac{t^{4}}{4}+iP\frac{t^{5}}{5!}\cdots\qquad$
$\displaystyle=$ $\displaystyle
I_{n}\left(1-\frac{t^{2}}{2}+\frac{t^{4}}{4}+\cdots\right)+iP\left(t-\frac{t^{3}}{3!}+\frac{t^{5}}{5!}+\cdots\right)$
$\displaystyle=$ $\displaystyle\cos(t)I_{n}+i\sin(t)P.$

where $t$ is a scalar parameter. If we deal with more than a single Pauli
group matrix, then we need to use Suzuki trotterization. For
$P_{l}\in\\{X,Y,Z,I\\}^{\otimes n}$, the first-order Suzuki trotterization is:

$e^{-i\sum_{l=1}^{m}s_{l}P_{l}}=\prod_{l=1}^{m}e^{-is_{l}P_{l}}+O(m^{2}s^{2}).$
(61)

where $s:=\max_{l}\,s_{l}$. So for large $t$, to control error, we need to
separate the evolution into multiple steps:

$e^{-i\sum_{l=1}^{m}s_{l}P_{l}}=\left(\prod_{k=1}^{m}e^{-i(s_{l}/r)P_{l}}\right)^{r}+O(m^{2}s^{2}/r).$
(62)

where $r$ is the number of evolution steps. There is also a second-order
formula for smaller error

$e^{-i\sum_{l=1}^{m}s_{l}P_{l}}=\left(\prod_{l=1}^{m}e^{-i\frac{s_{l}}{2r}P_{l}}\prod_{l=m}^{1}e^{-i\frac{s_{l}}{2r}P_{l}}\right)^{r}+O(m^{3}s^{3}/r^{2}).$
(63)

So, if we set the trotterization error level at $\epsilon$, then we need to
split the evolution into $r\in
O\left(\frac{m^{3/2}s^{3/2}}{\sqrt{\epsilon}}\right)$ steps with the second-
order formula. By choosing $t=-\frac{s}{2r}$, Eq.(60) gives:

$e^{-i\frac{s}{2r}P}=\cos\left(\frac{s}{2r}\right)I_{n}-i\sin\left(\frac{s}{2r}\right)P.$
(64)

As a result, Eq.(63) becomes:

$e^{-i\sum_{l=1}^{m}s_{l}P_{l}}\approx\left(\prod_{l=1}^{m}\left[\cos\left(\frac{s_{l}}{2r}\right)I_{n}-i\sin\left(\frac{s_{l}}{2r}\right)P_{l}\right]\prod_{l=m}^{1}\left[\cos\left(\frac{s_{l}}{2r}\right)I_{n}-i\sin\left(\frac{s_{l}}{2r}\right)P_{l}\right]\right)^{r}.$
(65)

In this case, we approximate the exponential of the linear combination of the
Pauli group matrices by another linear combination of Pauli group matrices,
where the expectation of the latter one can be evaluated in a gate-based
quantum computer easily. Note that in the $r=1$ case, if we have $m$ terms in
the power, we end up with at most $O(m^{2})$ terms after the transformation,
which is exactly the scenario we had in H4 examples in the main text.

## References

  * Luis and Peřina (1996) A. Luis and J. Peřina, Phys. Rev. A 54, 4564 (1996).
  * Cleve  _et al._ (1998) R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proc. R. Soc. Lond. A 454, 339 (1998).
  * Berry  _et al._ (2007) D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Comm. Math. Phys. 270, 359 (2007).
  * Childs (2010) A. M. Childs, Comm. Math. Phys. 294, 581 (2010).
  * Wecker, Hastings, and Troyer (2015a) D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev. A 92, 042303 (2015a).
  * Häner  _et al._ (2016) T. Häner, D. S. Steiger, M. Smelyanskiy, and M. Troyer, in  _SC ’16: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ (2016) pp. 866–874.
  * Poulin  _et al._ (2017) D. Poulin, A. Kitaev, D. S. Steiger, M. B. Hastings, and M. Troyer, arXiv preprint arXiv:1711.11025 (2017).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. R. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nat. Commun. 5, 4213 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Romero  _et al._ (2018) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Quantum Sci. Technol. 4, 014008 (2018).
  * Shen  _et al._ (2017) Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung, and K. Kim, Phys. Rev. A 95, 020501 (2017).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Colless  _et al._ (2018) J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Phys. Rev. X 8, 011021 (2018).
  * Huggins  _et al._ (2020) W. J. Huggins, J. Lee, U. Baek, B. O’Gorman, and K. B. Whaley, New J. Phys. 22, 073009 (2020).
  * Ryabinkin  _et al._ (2018) I. G. Ryabinkin, T.-C. Yen, S. N. Genin, and A. F. Izmaylov, Journal of chemical theory and computation 14, 6317 (2018).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Chemical reviews 119, 10856 (2019).
  * Ryabinkin  _et al._ (2020) I. G. Ryabinkin, R. A. Lang, S. N. Genin, and A. F. Izmaylov, J. Chem. Theory Comput. 16, 1055 (2020).
  * Izmaylov  _et al._ (2019) A. F. Izmaylov, T.-C. Yen, R. A. Lang, and V. Verteletskyi, J. Chem. Theory Comput. 16, 190 (2019).
  * Lang, Ryabinkin, and Izmaylov (2021) R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov, JJ. Chem. Theory Comput. 17, 66 (2021).
  * Grimsley  _et al._ (2019a) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nature communications 10, 1 (2019a).
  * Grimsley  _et al._ (2019b) H. R. Grimsley, D. Claudino, S. E. Economou, E. Barnes, and N. J. Mayhall, Journal of chemical theory and computation 16, 1 (2019b).
  * Cerezo  _et al._ (2021) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Nature Reviews Physics 3, 625 (2021).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Reviews of Modern Physics 92, 015003 (2020).
  * Bharti  _et al._ (2022) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, _et al._ , Reviews of Modern Physics 94, 015004 (2022).
  * Anand  _et al._ (2022) A. Anand, P. Schleich, S. Alperin-Lea, P. W. Jensen, S. Sim, M. Díaz-Tinoco, J. S. Kottmann, M. Degroote, A. F. Izmaylov, and A. Aspuru-Guzik, Chemical Society Reviews (2022).
  * Bartlett, Kucharski, and Noga (1989) R. J. Bartlett, S. A. Kucharski, and J. Noga, Chem. Phys. Lett. 155, 133 (1989).
  * Taube and Bartlett (2006) A. G. Taube and R. J. Bartlett, Int. J. Quantum Chem. 106, 3393 (2006).
  * Hoffmann and Simons (1988) M. R. Hoffmann and J. Simons, J. Chem. Phys. 88, 993 (1988).
  * Kutzelnigg (1991) W. Kutzelnigg, Theor. Chim. Acta 80, 349 (1991).
  * Evangelista, Chan, and Scuseria (2019) F. A. Evangelista, G. K.-L. Chan, and G. E. Scuseria, J. Chem. Phys. 151, 244112 (2019).
  * Kowalski (2021) K. Kowalski, Physical Review A 104, 032804 (2021).
  * Hill and Wheeler (1953) D. L. Hill and J. A. Wheeler, Physical Review 89, 1102 (1953).
  * Rodríguez-Guzmán, Egido, and Robledo (2002) R. Rodríguez-Guzmán, J. Egido, and L. Robledo, Nuclear Physics A 709, 201 (2002).
  * Bender, Heenen, and Reinhard (2003) M. Bender, P.-H. Heenen, and P.-G. Reinhard, Reviews of Modern Physics 75, 121 (2003).
  * Ring and Schuck (2004) P. Ring and P. Schuck,  _The nuclear many-body problem_ (Springer Science & Business Media, 2004).
  * Yao  _et al._ (2010) J. Yao, J. Meng, P. Ring, and D. Vretenar, Physical Review C 81, 044311 (2010).
  * Egido (2016) J. L. Egido, Physica Scripta 91, 073003 (2016).
  * Hizawa, Hagino, and Yoshida (2021) N. Hizawa, K. Hagino, and K. Yoshida, Physical Review C 103, 034313 (2021).
  * Klauder and Skagerstam (1985) J. R. Klauder and B.-S. Skagerstam,  _Coherent states: applications in physics and mathematical physics_ (World scientific, 1985).
  * Zhang, Gilmore  _et al._ (1990) W.-M. Zhang, R. Gilmore, _et al._ , Reviews of Modern Physics 62, 867 (1990).
  * Fukutome (1981a) H. Fukutome, Progress of Theoretical Physics 65, 809 (1981a).
  * Goeke and Reinhard (1980) K. Goeke and P.-G. Reinhard, Annals of Physics 124, 249 (1980).
  * Goeke and Reinhard (1978) K. Goeke and P.-G. Reinhard, Annals of Physics 112, 328 (1978).
  * Reinhard and Goeke (1978a) P.-G. Reinhard and K. Goeke, Physical Review C 17, 1249 (1978a).
  * Reinhard and Goeke (1978b) P.-G. Reinhard and K. Goeke, Nuclear Physics A 312, 121 (1978b).
  * Hehre, Stewart, and Pople (1969) W. J. Hehre, R. F. Stewart, and J. A. Pople, The Journal of Chemical Physics 51, 2657 (1969).
  * Jankowski and Paldus (1980) K. Jankowski and J. Paldus, International Journal of Quantum Chemistry 18, 1243 (1980).
  * Thouless (1960) D. J. Thouless, Nuclear Physics 21, 225 (1960).
  * Čížek and Paldus (1967) J. Čížek and J. Paldus, The Journal of Chemical Physics 47, 3976 (1967).
  * Seeger and Pople (1977) R. Seeger and J. A. Pople, The Journal of Chemical Physics 66, 3045 (1977).
  * Fukutome (1981b) H. Fukutome, International Journal of Quantum Chemistry 20, 955 (1981b).
  * tA v  _et al._ (2021) A. tA v _et al._ , “Qiskit: An open-source framework for quantum computing,”  (2021).
  * Childs and Wiebe (2012) A. M. Childs and N. Wiebe, Quantum Inf. Comput. 12 (2012), 10.26421/qic12.11-12.
  * Brassard  _et al._ (2002) G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, “Quantum amplitude amplification and estimation,” in  _Quantum Computation and Quantum Information_, Vol. 305, edited by J. S. J. Lomonaco (AMS Contemporary Mathematics, 2002) pp. 53–74.
  * Berry  _et al._ (2015) D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Phys. Rev. Lett. 114, 090502 (2015).
  * Low and Chuang (2017) G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118, 010501 (2017).
  * Gilyén  _et al._ (2019) A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) p. 193–204.
  * Izmaylov  _et al._ (2020) A. F. Izmaylov, T.-C. Yen, R. A. Lang, and V. Verteletskyi, J. Chem. Theory Comput. 16, 190 (2020).
  * Peng and Kowalski (2022) B. Peng and K. Kowalski, arXiv preprint arXiv:2206.09005 (2022).
  * Gonthier  _et al._ (2020) J. F. Gonthier, M. D. Radin, C. Buda, E. J. Doskocil, C. M. Abuan, and J. Romero, “Identifying challenges towards practical quantum advantage through resource estimation: the measurement roadblock in the variational quantum eigensolver,” (2020), arXiv:2012.04001 [quant-ph] .
  * Rubin, Babbush, and J. (2018) N. C. Rubin, R. Babbush, and M. J., New J. Phys. 20, 053020 (2018).
  * Wecker, Hastings, and Troyer (2015b) D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev. A 92, 042303 (2015b).
  * Huang, Kueng, and Preskill (2020) H. Y. Huang, R. Kueng, and J. Preskill, Nat. Phys. 16, 1050 (2020).
  * Huang, Kueng, and Preskill (2021) H.-Y. Huang, R. Kueng, and J. Preskill, Phys. Rev. Lett. 127, 030503 (2021).
  * Aaronson (2018) S. Aaronson, in  _Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2018 (Association for Computing Machinery, New York, NY, USA, 2018) pp. 325–338.
  * Struchalin  _et al._ (2021) G. I. Struchalin, Y. A. Zagorovskii, E. V. Kovlakov, S. S. Straupe, and S. P. Kulik, PRX Quantum 2, 010307 (2021).
  * Chen  _et al._ (2021) S. Chen, W. Yu, P. Zeng, and S. T. Flammia, arXiv preprint arXiv:2011.09636 (2021).
  * Zhao, Rubin, and Miyake (2021) A. Zhao, N. C. Rubin, and A. Miyake, Phys. Rev. Lett. 127, 110504 (2021).
  * Acharya, Saha, and Sengupta (2021) A. Acharya, S. Saha, and A. M. Sengupta, arXiv preprint arXiv:2105.05992 (2021).
  * Hadfield (2021) C. Hadfield, arXiv preprint arXiv:2105.12207 (2021).
  * Hillmich  _et al._ (2021) S. Hillmich, C. Hadfield, R. Raymond, A. Mezzacapo, and R. Wille, arXiv preprint arXiv:2105.06932 (2021).
  * Zhang  _et al._ (2021) T. Zhang, J. Sun, X.-X. Fang, X. Zhang, X. Yuan, and H. Lu, arXiv preprint arXiv:2106.10190 (2021).
  * Zheng  _et al._ (2022) M. Zheng _et al._ , “QuGCM,” https://github.com/pnnl/QuGCM (2022).
  * An and Lin (2022) D. An and L. Lin, ACM Transactions on Quantum Computing 3, 1 (2022).