# Dynamic, Symmetry-Preserving, and Hardware-Adaptable Circuits for Quantum
Computing Many-Body States and Correlators of the Anderson Impurity Model

Eric B. Jones<EMAIL_ADDRESS>Infleqtion, Boulder CO, 80301  Cody
James Winkleblack<EMAIL_ADDRESS>National Renewable Energy
Laboratory, Golden CO, 80401  Colin Campbell  Infleqtion, Boulder CO, 80301
Caleb Rotello  National Renewable Energy Laboratory, Golden CO, 80401  Edward
D. Dahl  IonQ, College Park, MD 20740  Matthew Reynolds  National Renewable
Energy Laboratory, Golden CO, 80401  Peter Graf  National Renewable Energy
Laboratory, Golden CO, 80401  Wesley Jones  National Renewable Energy
Laboratory, Golden CO, 80401

###### Abstract

We present a hardware-reconfigurable ansatz on $N_{q}$-qubits for the
variational preparation of many-body states of the Anderson impurity model
(AIM) with $N_{\text{imp}}+N_{\text{bath}}=N_{q}/2$ sites, which conserves
total charge and spin z-component within each variational search subspace. The
many-body ground state of the AIM is determined as the minimum over all minima
of $O(N_{q}^{2})$ distinct charge-spin sectors. Hamiltonian expectation values
are shown to require $\omega(N_{q})<N_{\text{meas.}}\leq
O(N_{\text{imp}}N_{\text{bath}})$ symmetry-preserving, parallelizable
measurement circuits, each amenable to post-selection. To obtain the one-
particle impurity Green’s function we show how initial Krylov vectors can be
computed via mid-circuit measurement and how Lanczos iterations can be
computed using the symmetry-preserving ansatz. For a single-impurity Anderson
model with a number of bath sites increasing from one to six, we show using
numerical emulation that the ease of variational ground-state preparation is
suggestive of linear scaling in circuit depth and sub-quartic scaling in
optimizer complexity. We therefore expect that, combined with time-dependent
methods for Green’s function computation, our ansatz provides a useful tool to
account for electronic correlations on early fault-tolerant processors.
Finally, with a view towards computing real materials properties of interest
like magnetic susceptibilities and electron-hole propagators, we provide a
straightforward method to compute many-body, time-dependent correlation
functions using a combination of time evolution, mid-circuit measurement-
conditioned operations, and the Hadamard test.

††preprint: APS/123-QED

##  I Introduction

The electronic structure problem—solving for the ground states of collections
of interacting electrons—is famously difficult. In its full generality, the
electronic structure problem resides in the Quantum Merlin Arthur (QMA)
complexity class, which is comprised of problems whose solutions are thought
to be difficult to find even on a quantum computer [1, 2]. In spite of this, a
number of classical approximation techniques, like density functional theory
(DFT) and GW theory, are known to describe many weakly- and moderately-
correlated electronic systems well [3, 4]. Meanwhile, it remains unclear what
degree of advantage quantum computing alone will confer over classical methods
for the electronic structure problem [5]. Moreover, quantum processors remain
relatively small and noisy, with the largest processors ranging from dozens to
hundreds of qubits with typical two-qubit gate infidelities on the orders of
.1% to 1% [6, 7, 8]. As such, a compelling near-term avenue to the
demonstration of quantum utility in the electronic structure of real materials
and chemical systems is to marry the scalability of weakly-correlated
theories, like DFT and GW, with theories that can treat strong electronic
correlations on a subsystem self-consistently via embedding theory [9, 10, 11,
12].

Dynamical mean-field theory (DMFT) is one such embedding theory, which is
particularly well-suited to periodic systems [13, 14]. In DMFT, one replaces a
lattice with interactions on every site 111Generally, a “site” can refer to an
atom, a collection of atoms, or some collection of interacting electronic
orbitals. Here we consider it a single spatially localized site with one spin-
up and one spin-down degree of freedom. with two, simultaneous but
complementary descriptions of the dynamics. The first description neglects all
interactions on the extended lattice but retains hopping between sites so that
the quasi-particle band picture of the electronic dynamics is considered
valid. For real materials, this description can be treated with theories like
DFT or GW. The second description considers only a specific site, an impurity
or cluster, within the lattice and treats all on-site interactions and
transitions exactly, but approximates the coupling of the site to the rest of
the lattice as electrons hopping into and out of an effective bath. This
description is similar to the mean-field treatment of the Ising model, where
the local Green’s function plays the analogous role of the on-site
magnetization and the dynamic hybridization function plays the role of the
mean-field [16]. Self-consistency between the non-interacting lattice and
impurity descriptions is enforced by insisting that the local lattice Green’s
function and self-energy are equivalent to the impurity Green’s function and
self-energy. Self-consistency is obtained by varying the hybridization
function, termed bath fitting, and constructs a solution to the original
problem, which can simultaneously describe both delocalized, band-like
behavior, and localized dynamical correlations.

The most computationally burdensome aspect of DMFT is solving for the impurity
Green’s function and, by extension, self-energy. Classically, this is often
done either via exact diagonalization (ED) [17], which scales exponentially in
the number of impurity and bath orbitals—the number of single-electron degrees
of freedom comprising the impurity and a discrete representation of the bath,
respectively—or via continuous-time quantum Monte Carlo (CTQMC), which scales
exponentially only in the number of impurity orbitals, but which suffers from
the fermionic sign problem [18, 19, 20]. And while recent work has proposed an
algorithm for solving impurity models with quasipolynomial runtime in the
number of bath orbitals using fermionic Gaussian states, the method still
scales exponentially in the number of impurity orbitals [21]. Exact
diagonalization has been used to solve impurities of up to about 25 (spatial)
orbitals, while more specialized quantum chemical approximation schemes have
treated systems of up to 127 orbitals [22, 23]. Tensor network and
renormalization group methods show the promise to more scalably solve the
impurity problem [24, 25, 26, 27, 28, 29], but the community has naturally
looked towards quantum computers as alternative impurity solvers.

Methods to solve the impurity problem via quantum circuits target either the
frequency-dependent Green’s function or the time-dependent Green’s function,
which must then be transformed to enforce self-consistency. Initial work by
Bauer et al. and Kreula et al. proposed measuring Green’s functions of the
single-impurity Anderson model (SIAM) on the real-time axis via Trotterized
time evolution and ancilla qubits [30, 31]. Subsequently, Rungger et al.
demonstrated the computation of the impurity Green’s function directly on the
real-frequency axis using the variational quantum eigensolver (VQE) to
construct many-body states in the Lehmann representation on quantum devices
[32]. Keen et al. fit real-time Green’s functions derived from Trotterized
evolution on a quantum device to extract frequency dependent Green’s functions
[33]. Lie algebraic fast-forwarding has extended the reach of near-term
processors to calculate real-time Green’s functions as was recently
demonstrated by Steckmann et al. [34]. And the quantum equation of motion
method was used to compute Green’s functions in the Lehmann representation to
perform DMFT on 14 qubits of an IBM quantum computer as well [35]. Still, many
of the aforementioned techniques suffer from problems of scalability, either
from too-deep circuits required by Trotter error [36], a combinatorial number
of many-body states, which need to be constructed in the Lehmann
representation [32], or an exponential growth in the Hamiltonian algebra as a
function of bath sites for fast-forwarding [34].

An alternative proposed recently by Jamet et al. constructs the impurity
Green’s function on the real or imaginary frequency axis via its continued
fraction expression and a variational version of the Lanczos iteration [37].
While the associated number of many-body states that need to be prepared on a
quantum computer also scales combinatorially in the number of total orbitals,
the associated Krylov subspace can often be truncated to provide a good
approximation to the Green’s function before the full Krylov space dimension
is reached. However, experimental implementations of the Krylov variational
quantum algorithm (KVQA) have been few [38]. The main challenge in
implementing a scalable version of the KVQA is to find an ansatz that does not
suffer from barren plateaus in its gradient, as does the low-depth hardware-
efficient ansatz (HEA) [39, 40], but which also does not suffer from high
Trotterization and SWAP routing overhead as do more physically-informed
ansatzae, such as quantum alternating operator and unitary coupled cluster
ansatzae [41, 42].

Figure 1:  Topology of the charge and spin conserving
$\mathbb{Z}_{2}$-symmetric Anderson impurity model with on-site interactions
in the the impurity. There are five impurity sites and four bath sites. The
left and right subgraphs demark the two spin registers. Light green solid
lines represent hopping between and on-site energies of impurity sites (light
green circles) given by the matrix, $h$, and a complete graph topology within
each spin register. Light green dashed lines represent on-site interactions
between spin-up and spin-down sites in the impurity given by the vector, $U$.
Dark purple lines depict hybridization between the impurity and bath sites
(dark purple circles) given by the matrix $V$, and forming a complete
bipartite graph with the impurity orbitals in each spin register. The diagonal
bath energies are contained in the vector $\epsilon$. In this work we restrict
the number of impurity sites to one.

Motivated by prior work on symmetry-preserving ansatzae (SPA) for VQE [43], we
provide a variational ansatz for the preparation of many-body states of the
Anderson impurity model (AIM), which describes the impurity problem in the
context of DMFT. Our ansatz has the properties that it is (1) physics-informed
and manifestly constrained by the symmetries of the AIM, enabling efficient
trainability and high expressibility at low gate depth; and (2) adaptable to
arbitrary hardware topologies with little overhead. Our ansatz differs from
that in Ref. [43] in that it explicitly accounts for interactions within the
impurity. In Sec. II we introduce the SPA and numerically analyze its
expressibility and trainability to prepare ground states of the single-
impurity Anderson model with a varying number of bath sites. Then, in Sec. III
we apply our ansatz to the construction of one-particle Green’s functions via
the variational preparation of a sequence of Lanczos vectors. We remark that
the focus of this paper is on assessing the expressivity and trainability of
the SPA for a large number of Hamiltonian parameter instantiations rather than
closing the DMFT loop. While we expect that our ansatz will find wide utility
in the preparation of many-body states of the AIM, we ultimately expect it to
find the greatest utility in schema where ground states of the AIM are encoded
on a quantum processor either variationally or deterministically via classical
tensor network methods and parameter fixing while Green’s functions are
computed using direct time evolution [27, 44]. For this reason, we propose a
new method for computing generic $m$-point correlation functions of fermionic
theories in Sec. IV. We discuss future avenues of work such as this in Sec. V.

##  II Symmetry-Preserving Ansatz

The Hamiltonian that governs the impurity dynamics in DMFT is the Anderson
impurity model (AIM),

$\begin{split}H&=H_{imp}+H_{hyb}+H_{bath},\\\
H_{imp}&=\sum_{i,j=1}^{N_{imp}}\sum_{\sigma=\downarrow}^{\uparrow}h_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+\sum_{i=1}^{N_{imp}}U_{i}n_{i\uparrow}n_{i\downarrow}\\\
H_{hyb}&=\sum_{i=1}^{N_{imp}}\sum_{b=1}^{N_{bath}}\sum_{\sigma=\downarrow}^{\uparrow}V_{ib}c_{i\sigma}^{\dagger}c_{b\sigma}+h.c.\\\
H_{bath}&=\sum_{b=1}^{N_{bath}}\sum_{\sigma=\downarrow}^{\uparrow}\epsilon_{b}c_{b\sigma}^{\dagger}c_{b\sigma}.\end{split}$
(1)

The AIM in Eq. 1, and depicted in Fig.1, is a slight generalization of the
SIAM typically used in DMFT calculations and represents the possibility to
include more than one impurity site, which can approximate strong spatial
correlations otherwise neglected in DMFT [45, 16]. $H_{imp}$ is comprised of
hopping terms ($h_{ij}$) within, on-site energies ($h_{ii}$) of, and
interactions ($U_{i}$) on, $N_{imp}$ impurity sites. $H_{hyb}$ represents
hopping, or hybridization, between impurity and bath states with strength
$V_{ib}$. And $H_{bath}$ forms a discrete representation of the bath with
$N_{bath}$ on-site energies $\epsilon_{b}$.

Figure 2:  Symmetry-preserving ansatzae in a two-dimensional, square lattice
architecture, $U_{SPA}(\bm{\theta})$, for (a) one-impurity, one-bath, (b) one-
impurity, two-bath, (c) one-impurity, three bath, and (d) two-impurity, two-
bath site Anderson impurity models. Solid green circles (lines) represent
qubits (planar interactions) with associated $Rz(\varphi_{i})$ (Givens and
$ZZ$-rotation) gates. Dashed lines represent long-range interactions. (e)
Compilation of Givens rotation into a standard gate-set. (f) Compilation of
$ZZ$-rotation into a standard gate-set.

As is well-known, Eq. 1 commutes with the total charge, $n^{+}$, and spin
z-component, $n^{-}$, operators,

$n^{\pm}=\sum_{i}(n_{i\uparrow}\pm n_{i\downarrow})+\sum_{b}(n_{b\uparrow}\pm
n_{b\downarrow}),$ (2)

each of which generates a symmetry of the AIM that we thus refer to a charge
and spin symmetry, respectively. In addition, Eq. 1 has an overall up-down
spin degeneracy—a $\mathbb{Z}_{2}$ symmetry, which can be broken, for example,
by the presence of an external magnetic field. Each of these symmetries is
manifestly preserved under the Jordan-Wigner (J-W) mapping from fermionic
operators to qubit operators,
$c^{\dagger}_{\mu}=\frac{1}{2}(X_{\mu}-iY_{\mu})\prod_{\nu<\mu}Z_{\nu}$, where
we organize qubit indices such that $\mu\in\\{1,\ldots,N_{imp}+N_{bath}\\}$
correspond to spin up qubits and
$\mu\in\\{N_{imp}+N_{bath}+1,\ldots,2N_{imp}+2N_{bath}=N_{q}\\}$ correspond to
spin down qubits. Due to the commutativity between Eq. 1 and the operators in
Eq. 2, the AIM can be diagonalized into simultaneous eigenstates of spin,
charge and energy, indicating that we can partition a variational search over
the full $N_{q}$-qubit Hilbert space for the ground state of Eq. 1 into a
collection of searches over smaller subspaces. There are $N_{q}-1$ non-trivial
charge sectors and within each charge sector, $N\in\\{1,\ldots,N_{q}-1\\}$,
there are $\min(N,N_{q}-N)+1$ spin sectors, of which only
$\lceil(\min(N,N_{q}-N)+1)/2\rceil$ are unique due to the $\mathbb{Z}_{2}$
symmetry of Eq. 1. Hence, if one does not know, a priori, which charge-spin
sector the ground state of a specified AIM resides in, one need only perform
$O(N_{q}^{2})$ separate runs of VQE to determine the global ground state. If
instead, one has a reliable reference state, one can run a symmetry-preserving
VQE routine within a single charge-spin sector. Within each charge-spin sector
the dimension of the symmetry-constrained subspace is
$D(N,S_{z})={N_{q}/2\choose N_{\uparrow}}{N_{q}/2\choose
N_{\downarrow}}={N_{q}/2\choose(N+S_{z})/2}{N_{q}/2\choose(N-S_{z})/2}$. It
can be verified that $\sum_{N}\sum_{S_{z}(N)}D(N,S_{z})=2^{N_{q}}$ as
expected.

Figure 3:  Expressibility of the symmetry-preserving ansatz on a two-
dimensional planar qubit topology for the single impurity Anderson model. (a)
Hamiltonian seed-averaged depth, $\langle d\rangle_{\text{param.}}$, at which
the ansatz is able to approximate the exact ground state as a function of the
number of sites, $N_{\text{imp}}+N_{\text{bath}}$ for various ground state
overlap error thresholds. Error bars represent one standard error of the mean
over the 150 Hamiltonian seeds. (b) Scaling of the ansatz depth,
$\max_{\text{param.}}(d)$, required to approximate the exact ground state for
the worst-case Hamiltonian seed at each system size. The blue line indicates
the line $\max_{\text{param.}}(d)=N_{\text{imp}}+N_{\text{bath}}$ for the most
stringent error threshold $\delta=10^{-5}$.

After performing the J-W transformation, an $(N,S_{z})$-sector can be
instantiated at the start of a circuit by performing
$N_{\uparrow}=(N+S_{z})/2$ bit flips on the all-zeros, $\ket{0}^{\bigotimes
N_{q}/2}_{\uparrow}$, spin-up qubit register and $N_{\downarrow}=(N-S_{z})/2$
bit flips on the all-zeros, $\ket{0}^{\bigotimes N_{q}/2}_{\downarrow}$, spin-
down qubit register. Following Ref. [43], we place the corresponding Pauli-$X$
gates as evenly-spaced as possible within each register. Then, the
interpretation of combined charge-spin symmetry in any subsequent VQE ansatz
is that operations should preserve total Hamming weight within each qubit
register. Note that these initial excitations should technically also include
Pauli-$Z$ gates from the J-W transformation, but their effect is only to
generate a global phase on a classical product state, which we can ignore.
Fig. 2 shows the form of the SPA for various system sizes. In addition to
symmetry considerations, the structural instantiation of the SPA is also
informed by the architecture of the processor on which it is run. For
instance, it can be easily adapted to a two-dimensional, (quasi-)planar qubit
connectivity such as some neutral atom [46] and superconducting [47]
platforms.

Fig. 2a depicts the SPA for a two-site AIM where from top left to bottom right
the qubits are ordered as
$(\text{Imp.},\uparrow),\>(\text{Bath},\uparrow),\>(\text{Imp.},\downarrow),\>(\text{Bath},\downarrow)$.
Mixing of different charge configurations within each spin register is
accomplished by parameterized Givens rotations, $G(\theta_{i})$ (note that
fractional $\text{iSWAP}(\theta_{i})$ gates could also be used here).
Meanwhile, Hamming weight preservation within each spin register is
accomplished by coupling the spin up and down registers only via two-qubit
phase rotations, $ZZ(\theta_{i})=e^{-i\theta_{i}ZZ/2}$, motivated by the
Hubbard $U$ interaction in Eq. 1. Finally, an independent $R_{Z}(\varphi_{i})$
rotation is placed on each qubit at the end of the layer to represent on-site
energies and cancel unwanted phases. This single-layer structure can then be
repeated $d$ times. Figs. 2b-c demonstrate the generalization of the ansatz
for AIMs with two and three bath sites, respectively.

The SPA shares similarities with the Hamiltonian Variational Ansatz (HVA) for
the AIM [48]. In particular, the hardware embedding of the SPA attempts to
exploit the sparsity structure of Eq. 1 wherein electrons can undergo
impurity-impurity and impurity-bath hopping, but not bath-bath hopping, as the
bath is generally assumed to be diagonal. On the other hand, the SPA will
leverage hardware connections, such as $ZZ(\theta_{i})$ gates between the spin
up and spin down qubits of bath sites B1 and B2 in Fig. 2b, which would not be
present in the HVA, as it would only couple spin up and down qubits on the
impurity. The effect of these inclusions is to increase expressivity while
maintaining symmetry constraints.

Figure 4: Ground state trainability of the symmetry-preserving ansatz on a
two-dimensional planar qubit topology for the single impurity Anderson model.
(a) On a semi-log plot, $\overline{\tilde{\text{nit}}}$ is the Hamiltonian
seed-averaged number of optimizer iterations, totaled over all charge-spin
sectors, required to converge to a target overlap error, $\delta$, for a
particular site number, $N_{\text{imp}}+N_{\text{bath}}$, normalized by the
number of charge-spin sectors searched over for the variational minimization,
which is $2\times(N_{\text{imp}}+N_{\text{bath}}+1)^{2}$. Error bars represent
standard error of the mean over Hamiltonian seeds. Solid black curve is a
visual cue for worst-case scaling at the $10^{-5}$ error threshold. (b)
$\overline{\tilde{\text{nit}}}$ as a function of
$\overline{\ell(\pmb{\theta}_{s})}$, which is the “Averaged Ansatz
Parameters”– the seed-averaged length of each $\bm{\theta}_{s}$ vector, where
the length is given by $d\times(N_{e}+N_{q})$, $N_{e}$ being the number of
edges in the SPA, and where the average is taken over a Hamltonian parameter
ensemble labeled by a target error threshold and a site number,
($N_{\text{bath}}+N_{\text{imp}}$). Blue line is a power-law fit and the black
markers above the main fit are worst-case instances.

Fig. 2d demonstrates how the SPA can be expanded to include more than one
impurity site, as well as to utilize longer-range gate interactions, in this
case represented by the next-nearest neighbor (NNN) green dashed lines. Here
again, the SPA corresponds closely to the HVA with explicit hopping terms
representable between both impurity and both bath orbitals. NNN interactions
could also generate $ZZ(\theta_{i})$ interaction gates between the two
impurity orbitals for AIMs with not stricly site-local electronic
interactions. Generalizing further, six bath orbitals could be coupled to each
of the impurity orbitals in Fig. 2 given NNNN connectivity (distance-2 plus
shapes). Quasiplanar neutral atom platforms are particularly well-suited to
such generalizations, with the Rydberg interaction able to naturally generate
long-range interactions [49, 50].

Let $\ket{\bm{0}}\equiv\ket{0}^{\bigotimes N_{q}}$. We denote a variational
trial state prepared by the SPA as

$\begin{split}\ket{N,S_{z};\bm{\theta}}&=U_{(N,S_{z})}(\bm{\theta})\ket{\bm{0}},\\\
U_{(N,S_{z})}(\bm{\theta})&=U_{SPA}(\bm{\theta})\prod_{\mu\in
init.ex.}X_{\mu},\end{split}$ (3)

where
$\bm{\theta}=(\theta_{1},\ldots,\varphi_{1},\ldots,\varphi_{N_{q}},\ldots,\varphi_{d\times
N_{q}})$, $U_{SPA}(\bm{\theta})$ is a circuit of the form defined in Fig. 2
(implicitly defined to some depth, $d$), and $U_{(N,S_{z})}(\bm{\theta})$
includes the initial bit flips to initialize the ansatz in the correct charge-
spin sector. An approximation for the global ground state is then determined
via minimizing the symmetry-constrained Hamiltonian minimization over all
charge-spin sectors:

$\ket{\widetilde{GS}}=\min_{N,S_{z}}\Big{[}\min_{\bm{\theta}}\braket{N,S_{z};\bm{\theta}}{H}{N,S_{z};\bm{\theta}}\Big{]}.$
(4)

In order to perform the minimization in Eq. 4 in a fully symmetry-preserving
fashion, one needs to be able to compute Hamiltonian expectation values from
measurements, i.e., operator strings, that also commute with Eq. 2.
Fortunately, the structure of Eq. 1 also facilitates this. Ignoring terms that
are proportional to the identity, there are only two types of terms in Eq. 1
after performing the J-W transformation. The first type of term is that which
can be computed from reading out all qubits in the computational $Z$-basis.
These are the terms that come from the Hubbard-$U$ interaction(s) and diagonal
terms of $h_{ij}$ in $H_{imp}$ as well as every term in $H_{bath}$. Any string
of Pauli-$Z$ operators commutes with Eq. 2. The other type of terms are the
hopping terms, which appear in $H_{hyb}$ and include the off-diagonal elements
of $h_{ij}$ in $H_{imp}$. A hopping term has the general form
$\propto(X_{\mu}X_{\nu}+Y_{\mu}Y_{\nu})\prod_{\mu<\rho<\nu}Z_{\rho}$, which
means its expectation value can be evaluated by performing a single Givens
rotation, $G(-\pi/4)$, on the qubit pair $(\mu,\nu)$, followed by readout of
all qubits in parallel. This readout scheme corresponds to measuring the
operator
$(1/2)(X_{\mu}X_{\nu}+Y_{\mu}Y_{\nu})\prod_{\rho\neq\mu,\nu}Z_{\rho}$, which
also commutes with Eq. 2 and whose eigenvectors thus conserve charge and spin
on the whole qubit array. In fact, multiple hopping terms can be measured in
parallel in this manner as long as the qubits that the associated Givens
rotations affect do not overlap. Measuring the expectation values of Eq. 1 in
this way, with only $Z$-basis and hopping circuits, ensures symmetries are
preserved and as an added benefit enables the use of symmetry-based post-
selection as an error mitigation technique [51, 52, 53]. In the absence of
measurement circuit parallelization, the number of circuits required to
measure the expectation value of Eq. 1 is
$1+2N_{imp}N_{bath}+N_{imp}(N_{imp}-1)$. In the presence of measurement
circuit parallelization, the spin-up and spin-down terms can be measured in
parallel, reducing each of the latter two terms by a factor of two. Either
way, an upper bound of $O(N_{imp}N_{bath})$ circuits is required to compute
the Hamiltonian expectation value, assuming $N_{\text{bath}}>N_{\text{imp}}$
as is nearly always the case in DMFT. Furthermore, if the parity terms in each
measurement circuit are ignored, parallelized measurement can first be
performed on the hopping terms between the impurity and bath orbitals, each
hopping term corresponding to an edge in a complete, bipartite graph. The
chromatic index of such a graph is $\max(N_{\text{imp}},N_{\text{bath}})$– the
maximum degree of any node therein. Meanwhile, each hopping term between two
impurity orbitals constitutes an edge in a complete graph. The chromatic index
of a complete graph is either $N_{\text{imp}}$ or $N_{\text{imp}}-1$ depending
on if $N_{\text{imp}}$ is odd or even, respectively. Hence, a loose lower
bound on the number of parallelized measurement circuits that need to be run
to compute the expectation value of Eq. 1 scales as
$1+N_{\text{bath}}+N_{\text{imp}}+(-1)\sim\omega(N_{q})$. However, we note
that the time complexity of measuring these observables can be converted to
space complexity by using $\omega(N_{q})<N_{\text{meas.}}\leq
O(N_{\text{imp}}N_{\text{bath}})$ measurement circuits run trivially in
parallel on $N_{\text{meas.}}$ smaller processors. In future work, it would be
interesting to investigate the applicability of recent number conserving or
fermionic shadow tomography techniques to measuring expectation values [54,
55].

Fig. 3 summarizes the finite-size expressivity of the SPA, Eq. 3, in its
ability to solve for the ground state of Eq. 1 via the minimization condition
in Eq. 4. For AIMs with $N_{\text{imp}}+N_{\text{bath}}\in\\{2,3,4,5,6,7\\}$
(corresponding to $N_{q}\in\\{4,6,8,10,12,14\\}$) with a single impurity, we
draw 150 AIM Hamiltonians from uniform parameter distributions
$h_{ij}\in[-5,5]$, $U_{i}\in[1,10]$, $V_{ib}\in[-5,5]$, and
$\epsilon_{b}\in[-5,5]$ (all energies being in electronvolts, eV). These
distributions reflect a range of impurity parameters typically associated with
performing DMFT on real materials. We then compute the ansatz depth, $d$,
required to converge the overlap of the variational ground state estimate with
the true ground state to below some target error,
$\delta=1-|\braket{\widetilde{GS}}{GS}|$. $|\braket{\widetilde{GS}}{GS}|$ is
the Uhlmann fidelity, which has the interpretation as a distance metric
between two rays in Hilbert space [56]. Fig. 3 shows the scaling of the
Hamiltonian parameter-averaged depth, $\langle d\rangle_{\text{param.}}$, as a
function of the number of AIM sites, $N_{\text{imp}}+N_{\text{bath}}$, for
various target ground state errors
$\delta\in\\{10^{-2},10^{-3},10^{-4},10^{-5}\\}$. Hamiltonian construction and
exact diagonalization were performed using OpenFermion [57], and all quantum
circuits, compiled to common elementary gates as shown in Fig. 2e-f, were
simulated in the absence of shot or gate noise using the Qulacs circuit
simulation framework [58]. To perform the classical minimiztion within each
charge-spin sector we used the SciPy implementation of the
Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm [59]. Of the 150 sets of
Hamiltonian parameters seeded for each system size,
$\\{123,135,142,140,146,144\\}$ seeds were used, respectively, to compute data
points and error bars, which represent standard error of the mean. All
discarded seeds, other than three seeds at seven sites discarded due to the
failure of the optimizer to converge after hitting an emulation wall-time,
involved a degeneracy in the ground state calculated by exact diagonalization.
In principle, our ground state preparation method can be extended to ground
states with a degeneracy by creating superpositions of charge-spin sector
ground states with equivalent energies, but for the sake of simplicity, we
discard these classical results.

While it is difficult to make conclusions regarding asymptotic performance
based on finite-size simulations, Fig. 3 suggests that the SPA has the generic
capacity to prepare ground states on $N_{q}$ qubits in low-order polynomial
depth—roughly linear for the most stringent target ground state overlap error.
Hence it can be regarded as efficient in its expressibility. For the purposes
of variational state preparation, however, trainability of an ansatz is also
important. To assess trainability of the SPA, we define a quantity,
$\overline{\tilde{\text{nit}}}$, which is the Hamiltonian seed-averaged number
of BFGS optimizer iterations, totaled over all charge-spin sectors, required
to converge to a target overlap error, $\delta$, for a particular site number,
$N_{\text{imp}}+N_{\text{bath}}$, normalized by the number of charge-spin
sectors searched over for the variational minimization, which is
$2\times(N_{\text{imp}}+N_{\text{bath}}+1)^{2}$. Because we have already
accounted for the quadratic overhead in searching through the different
charge-spin sectors, we are interested in the ability of a classical optimizer
to determine the ground state in some typical sense within any charge-spin
sector (averaged over many AIM Hamiltonian realizations). This is what the
scaling of $\overline{\tilde{\text{nit}}}$ addresses as a function of system
size in Fig. 4a, which shows sub-exponential scaling for all target error
thresholds on a semi-log plot. To more precisely assess the scaling of
$\overline{\tilde{\text{nit}}}$ as a function of ansatz depth, we define
another quantity, $\overline{\ell(\pmb{\theta}_{s})}$, which is the “Averaged
Ansatz Parameters”—the seed-averaged length of each $\bm{\theta}_{s}$ vector,
where the length is given by $d\times(N_{e}+N_{q})$, $N_{e}$ being the number
of edges in the SPA, and where the average is taken over a Hamltonian
parameter ensemble labeled by a target error threshold and a site number,
($N_{bath}+N_{imp}$). Fig. 4b shows how $\overline{\tilde{\text{nit}}}$ scales
with $\overline{\ell(\pmb{\theta}_{s})}$ on a log-log plot. Remarkably, the
number of optimizer calls required to find the ground state as a function of
the number of ansatz parameters scales roughly as a power law on average with
fitted exponent $\sim 1.769$.

Taken together, Figs. 3-4 indicate that the SPA is efficient in both its
expressibility and trainability in preparing ground states of the single-
impurity Anderson model as the discretization of the bath becomes more fine-
grained. Namely, Fig. 3a indicates that the average-case depth to reach the
ground state is no worse than roughly linearly-scaling. In Fig. 3b we show
that for the most stringent error threshold considered, $\delta=10^{-5}$, the
worst-case scaling in depth is exactly linear after (or excluding)
$N_{\text{imp}}+N_{\text{bath}}=3$, which is consistent with near-saturation
of the Lieb-Robinson bound [60, 61]. Moreover, given the numerical evidence in
Fig. 4, we conjecture that the number of optimizer calls scales, in the
average case, as not much worse than $\sim N_{q}^{3.538}$ in the number of
sites, since the number of ansatz parameters scales quadratically in the
number of sites (one factor per layer and linear depth to represent the ground
state). Note also that worst-case scaling over the parameter set in Fig. 4
scales similarly—being shifted up by a constant factor. Hence, if one is
provided a reliable reference state—that is, a known charge-spin sector within
which to search for the ground state—as might be determined by Hartree-Fock or
a higher-fidelity mean-field theory [62, 63, 4], $O(1)$ fixed space-
parallelized measurements, and parallelized gate operation within each ansatz
layer (as is usually assumed), then the finite-size scaling observed is
consistent with $\Theta(N_{q}^{4.538})$ runtime to prepare the ground state of
the AIM. Interestingly, this heuristic, average-case scaling is somewhat worse
than Bravyi and Gossett’s result, which established an algorithm for
estimating the ground state energy of quantum impurity models with runtime
$O(N_{\text{bath}}^{3})\exp[O(N_{\text{imp}}\log^{3}(N_{\text{imp}}\gamma^{-1}))]$
[21]. For $N_{\text{imp}}=1$, like is used in this work, this scaling reduces
to $O(N_{q}^{3})$ in the number of sites (and is quasi-linear in
$\gamma^{-1}$). However, ground-state preparation is only one of the core
subroutines involved in performing the impurity component of DMFT
computations, the other being one-particle Green’s function construction.
Without an efficient method to compute Green’s functions classically, an
efficient classical ground state preparation algorithm alone is insufficient
to perform DMFT efficiently. Hence, we now discuss an application of our
ansatz to the problem of computing Green’s functions via the continued
fraction representation in frequency space. Subsequently, we provide an
efficient method for computing $m$-point correlation functions directly in the
time domain.

##  III Symmetry-Preserving Krylov Variational Quantum Algorithm

In DMFT, once a ground state of the AIM has been prepared, the single-particle
Green’s function must be computed in order to enforce self-consistency of the
impurity self-energy with the local bath self-energy [14]. The KVQA aims to
compute the impurity retarded Green’s functions,

$G^{R}_{i\sigma}(z)=||c^{\dagger}_{i\sigma}\ket{GS}||^{2}g_{\phi^{+}_{i\sigma}}(z)-||c_{i\sigma}\ket{GS}||^{2}g_{\phi^{-}_{i\sigma}}(-z),$
(5)

via their continued fraction expansion,

$g_{\phi}(z)=\bra{\phi}[z-\tilde{H}]^{-1}\ket{\phi}=\frac{1}{z-a_{0}-\frac{b_{1}^{2}}{z-a_{1}-\frac{b_{2}^{2}}{z-a_{2}-\ldots}}},$
(6)

in terms of Krylov basis coefficients, $a_{n}$ and $b_{n}$, and where $z$ is a
complex frequency and $\tilde{H}=H-E_{GS}$ [37]. Eq. 5 is taken at zero
temperature, but the following arguments can be generalized to thermal
distributions. The $a_{n}$ and $b_{n}$ are obtained from Lanczos iterations,

$\begin{split}b_{n}^{2}&=\bra{\chi_{n-1}}\tilde{H}^{2}\ket{\chi_{n-1}}-a_{n-1}^{2}-b_{n-1}^{2}\\\
\ket{\chi_{n}}&=\frac{1}{b_{n}}[(\tilde{H}-a_{n-1})\ket{\chi_{n-1}}-b_{n-1}\ket{\chi_{n-2}}]\\\
a_{n}&=\bra{\chi_{n}}\tilde{H}\ket{\chi_{n}},\end{split}$ (7)

with $b_{0}=0$, $\ket{\chi_{0}}=\ket{\phi}$, and
$a_{0}=\bra{\chi_{0}}\tilde{H}\ket{\chi_{0}}$. Two sets of Lanczos iterations
need to be performed in order to compute Eq. 5, one for each of the initial
Krylov vectors,
$\ket{\phi^{+}_{i\sigma}}=c^{\dagger}_{i\sigma}\ket{GS}/||c^{\dagger}_{i\sigma}\ket{GS}||$
and $\ket{\phi^{-}_{i\sigma}}=c_{i\sigma}\ket{GS}/||c_{i\sigma}\ket{GS}||$.

First, we provide a new method to compute initial Krylov vectors. In prior
work, $\ket{\phi^{\pm}_{i\sigma}}$ have been computed variationally, which is
unintuitive since the application of a single creation or annihilation
operator to the ground state wavefunction is a non-unitary operation [37, 38].
Many quantum processors now offer the ability to measure qubits mid-circuit
and reset them [64, 65, 66]. Meanwhile, creation and annihilation operators
can be expressed as
$c^{\dagger}_{i\sigma}\propto(\prod_{\nu<(i,\sigma)}Z_{\nu})X_{i\sigma}P_{i\sigma}(0)$
and
$c_{i\sigma}\propto(\prod_{\nu<(i,\sigma)}Z_{\nu})X_{i\sigma}P_{i\sigma}(1)$,
where $P_{i\sigma}(0)$ and $P_{i\sigma}(1)$ are projectors onto the impurity
qubit $\ket{0}$ and $\ket{1}$ state, respectively. So, in order to prepare
each initial Krylov vector, $\ket{\phi^{+}_{i\sigma}}$
($\ket{\phi^{-}_{i\sigma}}$), one prepares the ground state of Eq. 1 via VQE
with the SPA, measures the state of the relevant impurity qubit until
$\ket{0}$ ($\ket{1}$) is measured, and then applies the operator
$(\prod_{\nu<(i,\sigma)}Z_{\nu})X_{i\sigma}$. If after $M$ tries, the result
$\ket{0}$ ($\ket{1}$) is not observed, it implies that the associated norm in
Eq. 5 is zero with certainty $\sim 1/\sqrt{M}$, and the corresponding Lanczos
iterations and continued fraction do not need to be computed. Also note that
by collapsing the wavefunction of the impurity qubit via mid-circuit
measurement, the initial Krylov vectors automatically become normalized. In
order to compute the unrenormalized modulus, one can use the ground state
wavefunction:
$||c^{\dagger}_{i\sigma}\ket{GS}||^{2}=(1+\bra{GS}Z_{i\sigma}\ket{GS})/2$ and
$||c_{i\sigma}\ket{GS}||^{2}=(1-\bra{GS}Z_{i\sigma}\ket{GS})/2$ as in [37].

Figure 5:  Expressivity of the symmetry-preserving ansatz for computing
Lanczos vectors. $\langle\epsilon_{\text{rel.}}\rangle_{\text{param.}}$ is the
relative error between the classical and variationally prepared Green’s
functions, averaged over Hamiltonian parameters and
$\bar{t}_{\text{wallclock}}$ is the parameter-averaged supercomputer wallclock
time, in seconds, required to converge all of the Lanczos interations and
compute the resulting Green’s function. Error bars are one standard error of
the mean over 100 parameter sets. Purple, violet, blue, and green points
represent 2, 3, 4, and 5 sites, respectively, while circles, squares,
diamonds, and stars represent a depth of 1, 2, 3, and 4 respectively. The
black dashed line demarks the $10\%$ relative error threshold. Note that
linear SPA depth is sufficient to prepare all Green’s functions to within
$10\%$ relative error up to 4-sites.

If $(N^{GS},S_{z}^{GS})$ is the charge-spin sector in which the ground state
of Eq. 1 resides, then the two associated Krylov subspaces in which the
Lanczos iterations will be performed are related via
$(N^{\pm},S_{z}^{\pm})=(N^{GS}\pm
1,S_{z}^{GS}\pm(\delta_{\sigma\uparrow}-\delta_{\sigma\downarrow}))$. In Ref.
[37], Jamet et al. observed that a sufficient condition to compute the second
line in Eq. 7 is to find a variational state, $\ket{\chi(\bm{\theta}_{n})}$,
such that
$\braket{\chi(\bm{\theta}_{n})}{\chi(\bm{\theta}_{n-2})}=\braket{\chi(\bm{\theta}_{n})}{\chi(\bm{\theta}_{n-1})}=0$
and $\braket{\chi(\bm{\theta}_{n})}{H}{\chi(\bm{\theta}_{n-1})}=b_{n}$, where
$\ket{\chi(\bm{\theta}_{n-2})}$ and $\ket{\chi(\bm{\theta}_{n-1})}$ are the
corresponding variational states at Lanczos iterations $n-1$ and $n-2$,
respectively. The variational Lanczos state can be obtained at each iteration
by minimization of the cost function

$\begin{split}F(\bm{\theta}_{n})=&\lambda_{1}(|\braket{N^{\pm},S_{z}^{\pm};\bm{\theta}_{n}}{H}{N^{\pm},S_{z}^{\pm};\bm{\theta}_{n-1}}|-|b_{n}|)^{2}\\\
&+\lambda_{2}|\braket{N^{\pm},S_{z}^{\pm};\bm{\theta}_{n}}{N^{\pm},S_{z}^{\pm};\bm{\theta}_{n-1}}|^{2}\\\
&+\lambda_{3}|\braket{N^{\pm},S_{z}^{\pm};\bm{\theta}_{n}}{N^{\pm},S_{z}^{\pm};\bm{\theta}_{n-2}}|^{2}\end{split}$
(8)

where all parametrized unitary evolutions are performed according to Eq. 3 in
the correct $(N^{\pm},S_{z}^{\pm})$-constrained Krylov subspaces, and
$(\lambda_{1},\lambda_{2},\lambda_{3})$ are optional Lagrange multipliers.

We compute the Lanczos iterations both classically via Eq. 7 and quantum-
variationally using Eq. 8 and use Eqs. 5-6 to construct exact,
$G^{R,\text{exact}}_{I,\uparrow}(\omega+i\eta)$, and variational
approximations to, $G^{R,\text{var.}}_{I,\uparrow}(\omega+i\eta)$, the
retarded impurity Green’s function. We choose the spin-up orbital without loss
of generality due to the system’s $\mathbb{Z}_{2}$ symmetry; and $\eta=0.1$ is
a small line-broadening parameter for visual aid, as is typically used [35].
To assess the quality of the variational Lanczos method as computed by
noiseless emulation, we compute the relative error between the exact and
quantum-variationally computed Green’s functions

$\epsilon_{\text{rel.}}=\frac{|G^{R,\text{var.}}_{I,\uparrow}(\omega+i\eta)-G^{R,\text{exact}}_{I,\uparrow}(\omega+i\eta)|}{|G^{R,\text{exact}}_{I,\uparrow}(\omega+i\eta)|}.$
(9)

Fig. 5 indicates the expressivity of our ansatz in variationally preparing the
Lanczos vectors used to construct the Green’s function. Specifically, the SPA
at linear depth, i.e., $N_{\text{bath}}+N_{\text{imp}}=d$ is sufficient to
compute all Green’s functions to within a seed-averaged relative error of
$10\%$ up to $N_{\text{bath}}+N_{\text{imp}}=4$, corresponding to eight
qubits. The average is taken over Hamiltonian parameter seeds. Fig. 6
demonstrates that a $10\%$ relative error in the Green’s function is often
sufficient to capture a large majority of the quantitatively important
features of the impurity response. We remark, however, that more data is
needed at larger system sizes to draw better conclusions about whether linear
ansatz depth is in general sufficient to prepare all the Lanczos vectors
required to compute Green’s functions to reasonable accuracy. Moreover, the
trainability of our ansatz using Eq. 8 as a cost function remains an open
question, with the linear-depth, 4-site Green’s function computation taking
roughly an order of magnitude longer than the linear-depth, 3-site computation
as measured by supercomputer wallclock mean time. In particular, one can see
that in the 4-site case, even though the Krylov subspace dimension is being
kept fixed, the wallclock time increases drastically as the depth of the
ansatz is increased, indicating that trainability is degrading rather than the
number of Lanczos vectors being the cause of prolonged runtime. However, given
that the SPA is efficiently trainable for ground state preparation, it is
unclear if the form of the cost function in Eq. 8, or the ansatz itself, or
the interplay between the two, is the cause of degraded trainability. We
include the 5-site, $d=1$ point for reference, but note that the simulation
wallclock time becomes untenable for the $(5,2)$ and larger cases. Hence
expanding our results in the future would benefit from more sophisticated
circuit simulation techniques like tensor network methods run on GPUs. These
methods, however, will not help with the combinatorial proliferation of
Lanczos vectors in the Krylov subspaces of larger system size.

In Fig. 6 we examine the relationship between the qualitative and quantitative
performance of variational Green’s function preparation for a selection of
Hamiltonian seeds at $N_{\text{imp}}+N_{\text{bath}}=d=4$. Seeds 18 (Fig. 6a)
and 84 (Fig. 6b) both exhibit both very good quantitative and qualitative
accuracy, with a relative error of $\sim.68\%$ and $\sim 5.7\%$, respectively.
The relative error in seed 84 arises from small amplitude fluctuations near
the Fermi level. Seed 24 (Fig.6c) exhibits a relative error of $\sim 30\%$
arising mainly from small shifts in peak positions and bifurcation in the main
peak near $\sim 10$ eV. Finally, seed 62 (Fig. 6) represents an example of a
variationally prepared Green’s function that, despite its large relative error
of $\sim 84\%$ due to small amplitude fluctuations, peak shifts, and peak
bifurcations, retains a good qualitative description of the response function.

Figs. 5-6 suggest that at even larger system sizes, the SPA may continue to be
a useful tool to prepare the many-body states necessary to approximate one-
particle Green’s functions in the frequency domain, either via the Lanczos
method or in other representations like the Lehmann representation. However,
using a quantum processor to prepare many-body states to compute Green’s
functions in the frequency domain is ultimately less efficient than using a
quantum processor to compute Green’s functions directly in the time domain.
There is a rigorous exponential separation between quantum and classical
computers’ ability to perform time evolution, while there often at best only a
polynomial separation for quantum state preparation [5]. In the next section
we propose a new algorithm for computing arbitrary $m$-point time correlation
functions on a quantum computer.

Figure 6:  Quantitative and qualitative accuracy of variationally prepared
Green’s functions. All subplots are for $N_{\text{imp}}+N_{\text{bath}}=d=4$.
Retarded Green’s functions are computed via the continued fraction
representation. Solid lines refer to Lanczos vectors and coefficients computed
variationally using the symmetry-preserving ansatz. Dashed lines refer to
Lanczos vectors and coefficients computed exactly. (a) Hamiltonian parameter
seed 18 has $\sim.68\%$ relative error and good very good qualitative
agreement between the exact and variationally prepared Green’s function. (b)
Seed 84 has $\sim 5.7\%$ relative error and exhibits only small-amplitude
fluctuations near the Fermi level. (c) Seed 24 has a $\sim 30\%$ relative
error and contains both small shifts in peak position as well as bifurcations
in large amplitude peaks. (d) Seed 62 has a relative error of $\sim 84\%$ due
to amplitude fluctuations, peak shifts, and peak bifurcations, and yet
qualitatively reproduces many of the main features of the spectrum.

##  IV  Computation of Many-Body Correlation Functions

While the computation of the single-particle Green’s function is sufficient in
DMFT to enforce self-consistency with the weakly-correlated bath theory,
higher-order correlation functions are often needed to accurately model
properties like magnetic susceptibilities and electron-hole propagators in
realistic simulations of strongly-correlated materials [16, 67]. Prior
deterministic theorems and algorithms have been provided to perform this
computation on a quantum computer (see, for example, [68, 69]). Here we
propose another such algorithm and provide an explicit circuit construction,
which combines unitary time evolution, mid-circuit measurement-conditioned
operations, and a modified Hadamard test.

Consider a collection of real-time fermionic operators,
$\\{f_{\alpha_{j}}(t_{j})\>|\>j=1\ldots m\\}$ in the Heisenberg picture

$f_{\alpha_{j}}(t)=e^{iHt}f_{\alpha_{j}}e^{-iHt},$ (10)

where each $f_{\alpha_{j}}$ stands either for a creation
($c_{\alpha_{j}}^{\dagger}$) or annihilation ($c_{\alpha_{j}}$) operator on
the spin-orbital $\alpha_{j}$ in the Schrodinger picture. A generic $m$-point
correlation function primitive, from which more complex objects like time-
ordered correlators can be constructed, is

$\begin{split}G(m,\ldots,1)&=\langle
f_{\alpha_{m}}(t_{m})f_{\alpha_{m-1}}(t_{m-1})\ldots
f_{\alpha_{1}}(t_{\alpha_{1}})\rangle\\\
&=\bra{GS}e^{iHt_{m}}f_{\alpha_{m}}e^{-iH(t_{m}-t_{m-1})}\ldots\\\
&\>\>\>\>\>\indent\indent\>\ldots
e^{-iH(t_{2}-t_{1})}f_{\alpha_{1}}e^{-iHt_{1}}\ket{GS}.\end{split}$ (11)

We compute correlations in the ground state, leaving generalizations to
thermal states to future work. The problem with computing generic $m$-point
correlators is that the operator string in Eq. 11 is not necessarily unitary,
precluding any straightforward application of the Hadamard test. Our approach
is therefore to compute Eq. 11 by modifying the Hadamard test to incorporate
controlled operators that can be constructed from common operations available
on a quantum computer, which is to say both unitary operators and projective
measurements.

The first step in this modification is to make the following definition,

$\tilde{f}_{\alpha}\equiv\prod_{\nu<\alpha}Z_{\nu}X_{\alpha}P_{\alpha}(z),$
(12)

which was previously used in Sec. III to compute initial Krylov vectors in the
Jordan-Wigner representation. Here again, $Z_{\nu}$ and $X_{\alpha}$ are Pauli
operators and $P_{\alpha}$ is a projector onto either the $\ket{z}=\ket{0}$ or
$\ket{z}=\ket{1}$ state of orbital $\alpha$ depending on if $f_{\alpha}$
refers to a creation or annihilation operator, respectively. The action of Eq.
12 on a generic state $\ket{\psi}$ is

$\tilde{f}_{\alpha}\ket{\psi}=\begin{cases}0&\text{if \>
$||f_{\alpha}\ket{\psi}||=0$}\\\
f_{\alpha}\ket{\psi}/||f_{\alpha}\ket{\psi}||&\text{else},\end{cases}$ (13)

which can be re-written compactly as

$f_{\alpha}\ket{\psi}=||f_{\alpha}\ket{\psi}||\times\tilde{f}_{\alpha}\ket{\psi}.$
(14)

We use Eq. 14 to re-write the operator string Eq. 11 in terms of the
renormalized operators, $\\{\tilde{f}_{\alpha_{j}}\\}$.

######  Theorem 1.

Take $t_{0}=0$ and $f_{\alpha_{0}}=1$ and let
$\ket{\psi_{j}}\equiv\prod_{k=j}^{1}e^{-iH(t_{k}-t_{k-1})}f_{\alpha_{k-1}}\ket{GS}$
be an un-normalized statevector. The $m$-point correlation function, Eq. 11,
can be computed as

$G(m,\ldots,1)=\tilde{G}(m,\ldots,1)\prod_{j=m}^{1}||f_{\alpha_{j}}\ket{\psi_{j}}||,$
(15)

where all $m$ factors $||f_{\alpha_{j}}\ket{\psi_{j}}||$ can be computed
recursively using $m$ quantum circuits to measure expectation values of simple
local observables of the form $f^{\dagger}_{\alpha_{j}}f_{\alpha_{j}}$, and if
all such factors are non-zero, then the quantity
$\tilde{G}(m,\ldots,1)=\langle\prod_{j=m}^{1}\tilde{f}_{\alpha_{j}}(t_{j})\rangle$
is computable by a straightforward application of the Hadamard test, since in
this instance the operator string,
$\prod_{j=m}^{1}\tilde{f}_{\alpha_{j}}(t_{j})$, is physically implementable.

###### Proof.

First, we show that Eq. 11 can be re-written as Eq. 15. We show this
recursively via construction, beginning with Eq. 11 expressed as

$G^{(m)}=\bra{GS}e^{iHt_{m}}f_{\alpha_{m}}\ket{\psi_{m}},$ (16)

where we have assigned the shorthand $G(m,\ldots,1)\rightarrow G^{(m)}$. Using
Eq. 13, consider the action of the last fermionic operator on $\ket{\psi_{m}}$
as the base case, and noting that norms commute through operators:

$\begin{split}e^{iHt_{m}}&f_{\alpha_{m}}\ket{\psi_{m}}=||f_{\alpha_{m}}\ket{\psi_{m}}||\times
e^{iHt_{m}}\tilde{f}_{\alpha_{m}}\ket{\psi_{m}}\\\
&=||f_{\alpha_{m}}\ket{\psi_{m}}||\times\tilde{f}_{\alpha_{m}}(t_{m})e^{iHt_{m-1}}f_{\alpha_{m-1}}\ket{\psi_{m-1}}\end{split}$
(17)

Now, assume that after $m-l+1$ replacements we have

$\begin{split}G^{(m)}=\bra{GS}\prod_{j=m}^{m-l+1}&||f_{\alpha_{j}}\ket{\psi_{j}}||\tilde{f}_{\alpha_{j}}(t_{j})\times\\\
&\times e^{iHt_{m-l}}f_{\alpha_{m-l}}\ket{\psi_{m-l}}\\\
=\bra{GS}\prod_{j=m}^{m-l}&||f_{\alpha_{j}}\ket{\psi_{j}}||\tilde{f}_{\alpha_{j}}(t_{j})\times\\\
&\times e^{iHt_{m-l-1}}f_{\alpha_{m-l-1}}\ket{\psi_{m-l-1}},\end{split}$ (18)

which can be carried all the way to the $j=1$ case, resulting in Eq. 15.

Next, we prescribe how to recursively compute the norms in Eq. 15. Consider,
first, the $j=1$ base case

$||f_{\alpha_{1}}\ket{\psi_{1}}||=||f_{\alpha_{1}}e^{-iHt_{1}}\ket{GS}||$ (19)

Eq. 19 is straightforwardly-computed as an expectation value of the local
operator $f_{\alpha_{1}}^{\dagger}f_{\alpha_{1}}$ sampled from the output of
the quantum circuit $e^{-iHt_{1}}\ket{GS}$. Sampling the circuit output $M$
times will give the expectation value with precision $\sim 1/\sqrt{M}$ from
standard shot statistics. For illustration purposes consider the $j=2$ case:

$\begin{split}||f_{\alpha_{2}}\ket{\psi_{2}}||&=||f_{\alpha_{2}}e^{-iH(t_{2}-t_{1})}f_{\alpha_{1}}\ket{\psi_{1}}||\\\
&=||f_{\alpha_{2}}e^{-iH(t_{2}-t_{1})}\tilde{f}_{\alpha_{1}}\ket{\psi_{1}}||\times||f_{\alpha_{1}}\ket{\psi_{1}}||.\end{split}$
(20)

Hence, we can compute Eq. 20 with the previously-computed result of Eq. 19 and
by sampling the output of the new circuit,
$e^{-iH(t_{2}-t_{1})}\tilde{f}_{\alpha_{1}}e^{-iHt_{1}}\ket{GS}$, to compute
the expecation value of the new local observable,
$f_{\alpha_{2}}^{\dagger}f_{\alpha_{2}}$. This methodology generalizes so that
if one has already computed the factors
$\\{||f_{\alpha_{j}}\ket{\psi_{j}}||\>\text{for}\>j=1,\ldots,l-1\\}$, then the
$l^{\text{th}}$ factor can be computed as

$\begin{split}||f_{\alpha_{l}}\ket{\psi_{l}}||=||f_{\alpha_{l}}&\prod_{k=l-1}^{0}e^{-iH(t_{k+1}-t_{k})}\tilde{f}_{\alpha_{k}}\ket{GS}||\\\
\times&\prod_{j=l-1}^{1}||f_{\alpha_{j}}\ket{\psi_{j}}||.\end{split}$ (21)

Altogether, evaluation of Eq. 15 requires sampling the output of $m$ quantum
circuits a total of $M$ times each, as well as implementing the Hadamard test
twice to compute the real and imaginary parts of $\tilde{G}^{(m)}$. ∎

From a practical perspective one uses Thm. 1 in the following manner. First,
the norms in Eq. 15 are computed recursively as in Eqs. 19-21. If at any
point, one of the norms vanishes, the recursion can be stopped and the
correlation function is zero. Only if the product of norms is computed to be
non-zero to within the desired precision of the calculation should the
Hadamard test be performed to compute $\tilde{G}^{(m)}$, since it is only in
this instance that the renormalized Heisenberg operator string can be
considered as non-vanishing and physically implementable using standard
quantum circuit instructions. For concreteness we discuss a specific example
below.

###  IV.1 Single-Particle Green’s Functions

The greater and lesser Green’s functions are defined, respectively as [16]

$G^{>}(2,1)=-i\langle
c_{\alpha_{2}}(t_{2})c_{\alpha_{1}}^{\dagger}(t_{1})\rangle$ (22)

and

$G^{<}(2,1)=i\langle
c_{\alpha_{1}}^{\dagger}(t_{1})c_{\alpha_{2}}(t_{2})\rangle,$ (23)

from which other Green’s functions can be constructed. For instance, using the
Heavyside step function, the retarded Green’s function is

$G^{R}(2,1)=\Theta(t_{2}-t_{1})(G^{>}(2,1)-G^{<}(2,1)).$ (24)

Taking Eq. 22 as an example, defining the standard time-translationally
invariant quantity $t\equiv t_{2}-t_{1}$, and computing the expectation value
in the ground state via Thm. 1, we find

$G^{>}(2,1)=-i\bra{GS}\tilde{c}_{\alpha_{2}}(t)\tilde{c}^{\dagger}_{\alpha_{1}}\ket{GS}||c_{\alpha_{2}}\ket{\psi_{2}}||\>||c^{\dagger}_{\alpha_{1}}\ket{GS}||$
(25)

where
$||c^{\dagger}_{\alpha_{1}}\ket{GS}||=\sqrt{(1+\bra{GS}Z_{\alpha_{1}}\ket{GS})/2}$
is computed simply by reading out the ground state wavefuction in the
computational $z$-basis, as previously noted (Fig. 7a), and

$||c_{\alpha_{2}}\ket{\psi_{2}}||=||c^{\dagger}_{\alpha_{1}}\ket{GS}||\sqrt{(1-\bra{\tilde{\psi}_{2}}Z_{\alpha_{2}}\ket{\tilde{\psi}_{2}})/2},$
(26)

with the state
$\ket{\tilde{\psi}_{2}}=e^{-iHt}\tilde{c}^{\dagger}_{\alpha_{1}}\ket{GS}$
prepared from Eq. 12 applied to the ground state followed by Hamiltonian time
evolution for a time $t$ (Fig. 7b). Finally, if both norms are non-zero, the
real and imaginary components of the propagator in Eq. 25 are computed by two
applications of the non-unitary Hadamard test (Fig. 7c).

Figure 7:  Circuits to compute $m$-point correlation functions of fermionic
Hamiltonians. (a) Ground state preparation and readout in the $z$-basis to
compute the expectation value $||c_{\alpha_{1}}^{\dagger}\ket{GS}||$. (b)
Preparation of the state $\ket{\tilde{\psi}_{2}}$ and readout in the $z$-basis
to compute $||c_{\alpha_{2}}\ket{\psi_{2}}||$. This circuit can be generalized
to compute an arbitrary norm $||f_{\alpha_{j}}\ket{\psi_{j}}||$. (c) The
modified Hadamard test when the non-unitary Heisenberg operator string is the
specific instance
$\prod_{j=m}^{1}\tilde{f}_{\alpha_{j}}(t_{j})=\tilde{c}_{\alpha_{2}}(t)\tilde{c}^{\dagger}_{\alpha_{1}}$.
The optional implementation of the $S^{\dagger}$ gate dictates if the real or
imaginary part is computed. (d) Any controlled non-unitary Heisenberg fermion
string can be implemented by a sequence of unitary time evolutions and single
controlled Schrödinger fermion gates. (e) Example implementation of a single
controlled Schrödinger fermion. In this instance a creation operator,
$\tilde{c}_{\alpha_{1}}^{\dagger}=Z_{1}\otimes Z_{2}\otimes X_{3}P_{3}(0)$,
needs to be applied to the spin-orbital $\alpha_{1}$, which has been assigned
to the third qubit, conditioned on if the Hadamard test ancilla is in the
$\ket{1}$ state. Since the gate must be coherently controlled, we push the
projective measurement past the classically-controlled $Z$ gates using the
principle of deferred measurement and commute it past the $X$ gate. Then, the
entire sequence of controlled operations can be represented by a string of
Toffoli-equivalent gates, a CNOT, and a controlled-projector, which uses an
additional ancilla.

Fig 7d shows the decomposition of an arbitrary controlled Heisenberg fermion
string into a sequence of time evolutions and controlled Schrödinger fermion
operators. Note that if the control qubit is in the $\ket{0}$ state, none of
the fermion operators are applied and all the forward and backward time-
evolution operators multiply to the identity. Using the principle of deferred
measurement and the relationship
$X_{\alpha_{j}}P_{\alpha_{j}}(0)=P_{\alpha_{j}}(1)X_{\alpha_{j}}$ applied to
$\tilde{c}^{\dagger}_{\alpha_{1}}$ as an example, we decompose each fermion
operator into $\alpha_{j}$ Toffoli-equivalent gates (up to single-qubit
rotations), one CNOT, and a single ancilla, post-selected on the $\ket{0}$
state used to implement a controlled projection operator, C$P_{\alpha_{j}}(1)$
(Fig. 7e). The ancilla can be directly reused in the implementation of
subsequent fermionic operators. To minimize the number of Toffolis required to
compute $m$-point impurity correlation functions, one should therefore index
the impurity as low as possible and if there are multiple impurity sites, they
should be grouped together by index. In future work, Thm. 1 can be applied
directly, along with the various circuit decompositions in this section, to
the computation of quantities of typical interest in DMFT, such as impurity
electron-hole propagators and magnetic susceptibilities [70, 16, 71, 72].

##  V  Conclusion

Herein, we have provided two explicit circuit constructions for efficiently
computing the central quantities in dynamical mean-field theory. The first
construction is a hardware-adaptable and symmetry-preserving variational
ansatz that can be used to prepare many-body states of the Anderson impurity
model. We show using numerical emulation that this ansatz can prepare ground
states of the single-impurity Anderson model in depth roughly linear in the
number of bath orbitals and in sub-quartic training time for a set of small to
moderate-size models, indicating both good expressibility near the Lieb-
Robinson bound and efficient trainability. Moreover, we show that the ansatz
can be used to prepare other many-body states like the Lanczos vectors used to
compute single-particle Green’s functions in the continued fraction
representation, although the scalability of this method remains in question,
both from the perspective of expressibility and trainability and due to the
fact that an exponential number of Lanczos iterations needs to be performed to
compute exact Green’s functions in the worst case. As such, the second circuit
construction is a new method for computing arbitrary $m$-point correlation
functions of fermionic systems that uses a combination of time-evolution, mid-
circuit measurement, and a modified version of the Hadamard test. Aside from
state preparation, which our first construction addresses, and the costs
involved with time-evolving under the Anderson impurity model Hamiltonian, the
main cost of this method is in the Toffoli complexity, which in the Jordan-
Wigner representation, is linear in both $m$ and $N_{\text{imp}}$ given that
the impurity orbitals can be indexed such that they are grouped together and
assigned low indices.

Our results support the notion that electronic structure computations of
strongly correlated materials that use quantum processors as impurity solvers
alongside state-of-the-art classical mean-field theories constitute a
promising path forward to practical quantum advantage. This notion is based on
three observations made here and elsewhere [27]. First, state-of-the-art mean-
field theories like GW theory are often sufficient to capture a large majority
of the weakly-correlated physics and chemistry in describing real materials
[67]. Second, ground-state preparation of the Anderson impurity model with
$O(1)$ impurity orbitals is quasipolynomially efficient both classically and
quantumly [21]. And third, computation of time-dependent impurity response
functions, like one- and two-body Green’s functions, involves Hamiltonian time
evolution, for which there is a known exponential runtime classically [73].
Meanwhile, we show by construction, that $m$-point correlation functions can
be computed efficiently on a quantum computer given access to mid-circuit
measurement, corroborating theoretical results suggesting this to be true [68,
69]. These observations indicate that classical computers are best suited to
solving the weakly-correlated bath theory, validating and helping to prepare
the ground state of the impurity model, and enforcing DMFT self-consistency,
while the quantum computer is best suited to helping to solve for and encoding
the ground state and deterministically computing correlation functions in the
encoded ground state. Interestingly, such a workflow obviates the need for
typical quantum chemistry subroutines like quantum phase estimation (although
it may be used to increase the overlap with the true ground state). Moreover,
the fidelity of the overall computation can be systematically improved as the
fidelity of quantum computers improves, by enlarging the number of impurity
orbitals to account for more spatial correlation or by using more orbitals to
present a more fine-grained description the bath.

Within this computational paradigm, there are a number of open questions and
potential avenues for future research. First and foremost is to assess the
expressibility and trainability of our ansatz at both greater system sizes,
$N_{\text{bath}}>6$, and with more impurity sites $N_{\text{imp}}>1$,
especially when simulated and trained on a classical computer using fermionic
tensor network compression techniques [44] and in the presence of noise. Along
these lines, it would be interesting to investigate augmentations of our
ansatz that include cooling mechanisms to remove errors via auxillary qubit
reset as have been proposed and demonstrated recently in a number of spin
models [74, 75, 76]. Such augmentations may prove useful in preparing both
ground and thermal states of the Anderson impurity model, which could have
applicability in the quantum computation of high-temperature superconducting
materials. In addition, recent results in dynamic circuits for preparing
matrix product and higher-dimension tensor network states have shown promise
for state preparation in constant and sub-Lieb-Robinson depth [77, 78].
Evaluating the applicability of these methods to AIM state preparation is
likely to be fruitful. Finally, the theoretical method devised herein to
compute $m$-point correlation functions should be experimentally validated and
used to perform resource estimations on when leveraging quantum computers to
compute two-particle observables might confer quantum advantage in describing
the physics of strongly correlated materials.

###### Acknowledgements.

E. B. J. thanks François Jamet and Ivan Rungger for discussions and material
on the Krylov variational quantum algorithm; Swagata Acharya and Mark van
Schilfgaarde for discussions on the foundations of DMFT and the properties of
the Anderson impurity model; Ben Hall and David Owusu-Antwi for discussions on
symmetry-preserving ansatzae; and Rich Rines for discussions on circuit
decompositions. This work was authored in part by the National Renewable
Energy Laboratory, managed and operated by the Alliance for Sustainable
Energy, LLC for the U.S. Department of Energy (DOE) under Contract No. DE-
AC36-08G028308. A portion of this research was performed using computational
resources sponsored by the U.S. Department of Energy’s Office of Energy
Efficiency and Renewable Energy and located at the National Renewable Energy
Laboratory. The views expressed in the article do not necessarily represent
the views of the DOE or the U.S. Government. The U.S. Government retains and
the publisher, by accepting the article for publication, acknowledges that the
U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide
license to publish or reproduce the published form of this work, or allow
others to do so. for U.S. Government purposes.

## References

  * Schuch and Verstraete [2009] N. Schuch and F. Verstraete, Computational complexity of interacting electrons and fundamental limitations of density functional theory, Nature physics 5, 732 (2009).
  * O’Gorman  _et al._ [2022] B. O’Gorman, S. Irani, J. Whitfield, and B. Fefferman, Intractability of electronic structure in a fixed basis, PRX Quantum 3, 020322 (2022).
  * Jones [2015] R. O. Jones, Density functional theory: Its origins, rise to prominence, and future, Rev. Mod. Phys. 87, 897 (2015).
  * van Schilfgaarde  _et al._ [2006] M. van Schilfgaarde, T. Kotani, and S. Faleev, Quasiparticle self-consistent $gw$ theory, Phys. Rev. Lett. 96, 226402 (2006).
  * Lee  _et al._ [2023] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, _et al._ , Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry, Nature Communications 14, 1952 (2023).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Moses  _et al._ [2023] S. Moses, C. Baldwin, M. Allman, R. Ancona, L. Ascarrunz, C. Barnes, J. Bartolotta, B. Bjork, P. Blanchard, M. Bohn, _et al._ , A race track trapped-ion quantum processor, arXiv preprint arXiv:2305.03828 (2023).
  * Morvan  _et al._ [2023] A. Morvan, B. Villalonga, X. Mi, S. Mandra, A. Bengtsson, P. Klimov, Z. Chen, S. Hong, C. Erickson, I. Drozdov, _et al._ , Phase transition in random circuit sampling, arXiv preprint arXiv:2304.11119 (2023).
  * Tomczak [2015] J. M. Tomczak, Qsgw+ dmft: an electronic structure scheme for the iron pnictides and beyond, in  _Journal of Physics: Conference Series_ , Vol. 592 (IOP Publishing, 2015) p. 012055.
  * Choi  _et al._ [2016] S. Choi, A. Kutepov, K. Haule, M. van Schilfgaarde, and G. Kotliar, First-principles treatment of mott insulators: linearized qsgw+ dmft approach, npj Quantum Materials 1, 1 (2016).
  * Lee and Haule [2017] J. Lee and K. Haule, Diatomic molecule as a testbed for combining dmft with electronic structure methods such as $gw$ and dft, Phys. Rev. B 95, 155104 (2017).
  * Clinton  _et al._ [2024] L. Clinton, T. Cubitt, B. Flynn, F. M. Gambetta, J. Klassen, A. Montanaro, S. Piddock, R. A. Santos, and E. Sheridan, Towards near-term quantum simulation of materials, Nature Communications 15, 211 (2024).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. Marianetti, Electronic structure calculations with dynamical mean-field theory, Reviews of Modern Physics 78, 865 (2006).
  * Note [1] Generally, a “site” can refer to an atom, a collection of atoms, or some collection of interacting electronic orbitals. Here we consider it a single spatially localized site with one spin-up and one spin-down degree of freedom.
  * Martin  _et al._ [2016] R. M. Martin, L. Reining, and D. M. Ceperly,  _Interacting Electrons: Theory and Computational Approaches_ (Cambridge University Press, 2016).
  * Caffarel and Krauth [1994] M. Caffarel and W. Krauth, Exact diagonalization approach to correlated fermions in infinite dimensions: Mott transition and superconductivity, Phys. Rev. Lett. 72, 1545 (1994).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time monte carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011).
  * Werner and Millis [2006] P. Werner and A. J. Millis, Hybridization expansion impurity solver: General formulation and application to kondo lattice and two-orbital models, Phys. Rev. B 74, 155107 (2006).
  * Haule [2007] K. Haule, Quantum monte carlo impurity solver for cluster dynamical mean-field theory and electronic structure calculations with adjustable cluster base, Phys. Rev. B 75, 155113 (2007).
  * Bravyi and Gosset [2017] S. Bravyi and D. Gosset, Complexity of quantum impurity problems, Communications in Mathematical Physics 356, 451 (2017).
  * Sangiovanni  _et al._ [2006] G. Sangiovanni, A. Toschi, E. Koch, K. Held, M. Capone, C. Castellani, O. Gunnarsson, S.-K. Mo, J. W. Allen, H.-D. Kim, A. Sekiyama, A. Yamasaki, S. Suga, and P. Metcalf, Static versus dynamical mean-field theory of mott antiferromagnets, Phys. Rev. B 73, 205121 (2006).
  * Zhu and Chan [2021] T. Zhu and G. K.-L. Chan, Ab initio full cell $gw+\mathrm{DMFT}$ for correlated materials, Phys. Rev. X 11, 021006 (2021).
  * Wolf  _et al._ [2014] F. A. Wolf, I. P. McCulloch, O. Parcollet, and U. Schollwöck, Chebyshev matrix product state impurity solver for dynamical mean-field theory, Phys. Rev. B 90, 115124 (2014).
  * Ganahl  _et al._ [2015] M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Held, and F. Verstraete, Efficient dmft impurity solver using real-time dynamics with matrix product states, Phys. Rev. B 92, 155132 (2015).
  * Wolf  _et al._ [2015] F. A. Wolf, A. Go, I. P. McCulloch, A. J. Millis, and U. Schollwöck, Imaginary-time matrix product state impurity solver for dynamical mean-field theory, Phys. Rev. X 5, 041032 (2015).
  * Jamet  _et al._ [2023] F. Jamet, C. Lenihan, L. P. Lindoy, A. Agarwal, E. Fontana, B. A. Martin, and I. Rungger, Anderson impurity solver integrating tensor network methods with quantum computing, arXiv preprint arXiv:2304.06587 (2023).
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Numerical renormalization group method for quantum impurity systems, Reviews of Modern Physics 80, 395 (2008).
  * García  _et al._ [2004] D. J. García, K. Hallberg, and M. J. Rozenberg, Dynamical mean field theory with the density matrix renormalization group, Phys. Rev. Lett. 93, 246403 (2004).
  * Bauer  _et al._ [2016] B. Bauer, D. Wecker, A. J. Millis, M. B. Hastings, and M. Troyer, Hybrid quantum-classical approach to correlated materials, Phys. Rev. X 6, 031045 (2016).
  * Kreula  _et al._ [2016] J. M. Kreula, L. García-Álvarez, L. Lamata, S. R. Clark, E. Solano, and D. Jaksch, Few-qubit quantum-classical simulation of strongly correlated lattice fermions, EPJ Quantum Technology 3, 1 (2016).
  * Rungger  _et al._ [2019] I. Rungger, N. Fitzpatrick, H. Chen, C. Alderete, H. Apel, A. Cowtan, A. Patterson, D. M. Ramo, Y. Zhu, N. H. Nguyen, _et al._ , Dynamical mean field theory algorithm and experiment on quantum computers, arXiv preprint arXiv:1910.04735 (2019).
  * Keen  _et al._ [2020] T. Keen, T. Maier, S. Johnston, and P. Lougovski, Quantum-classical simulation of two-site dynamical mean-field theory on noisy quantum hardware, Quantum Science and Technology 5, 035001 (2020).
  * Steckmann  _et al._ [2023] T. Steckmann, T. Keen, E. Kökcü, A. F. Kemper, E. F. Dumitrescu, and Y. Wang, Mapping the metal-insulator phase diagram by algebraically fast-forwarding dynamics on a cloud quantum computer, Phys. Rev. Res. 5, 023198 (2023).
  * Selisko  _et al._ [2024] J. Selisko, M. Amsler, C. Wever, Y. Kawashima, G. Samsonidze, R. U. Haq, F. Tacchino, I. Tavernelli, and T. Eckl, Dynamical mean field theory for real materials on a quantum computer, arXiv preprint arXiv:2404.09527 (2024).
  * Jaderberg  _et al._ [2020] B. Jaderberg, A. Agarwal, K. Leonhardt, M. Kiffner, and D. Jaksch, Minimum hardware requirements for hybrid quantum–classical dmft, Quantum Science and Technology 5, 034015 (2020).
  * Jamet  _et al._ [2021] F. Jamet, A. Agarwal, C. Lupo, D. E. Browne, C. Weber, and I. Rungger, Krylov variational quantum algorithm for first principles materials simulations, arXiv preprint arXiv:2105.13298 (2021).
  * Backes  _et al._ [2023] S. Backes, Y. Murakami, S. Sakai, and R. Arita, Dynamical mean-field theory for the hubbard-holstein model on a quantum device, Phys. Rev. B 107, 165155 (2023).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature communications 9, 4812 (2018).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature communications 12, 1791 (2021).
  * Kremenetski  _et al._ [2021] V. Kremenetski, T. Hogg, S. Hadfield, S. J. Cotton, and N. M. Tubman, Quantum alternating operator ansatz (qaoa) phase diagrams and applications for quantum chemistry, arXiv preprint arXiv:2108.13056 (2021).
  * 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, A quantum computing view on unitary coupled cluster theory, Chemical Society Reviews 51, 1659 (2022).
  * Gard  _et al._ [2020] B. T. Gard, L. Zhu, G. S. Barron, N. J. Mayhall, S. E. Economou, and E. Barnes, Efficient symmetry-preserving state preparation circuits for the variational quantum eigensolver algorithm, npj Quantum Information 6, 10 (2020).
  * Provazza  _et al._ [2024] J. Provazza, K. Gunst, H. Zhai, G. K.-L. Chan, T. Shiozaki, N. C. Rubin, and A. F. White, Fast emulation of fermionic circuits with matrix product states, Journal of Chemical Theory and Computation (2024).
  * He and Lu [2014] R.-Q. He and Z.-Y. Lu, Quantum renormalization groups based on natural orbitals, Physical Review B 89, 085108 (2014).
  * Graham  _et al._ [2022] T. Graham, Y. Song, J. Scott, C. Poole, L. Phuttitarn, K. Jooya, P. Eichler, X. Jiang, A. Marra, B. Grinkemeyer, _et al._ , Multi-qubit entanglement and algorithms on a neutral-atom quantum computer, Nature 604, 457 (2022).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Libbi  _et al._ [2022] F. Libbi, J. Rizzo, F. Tacchino, N. Marzari, and I. Tavernelli, Effective calculation of the green’s function in the time domain on near-term quantum processors, Physical Review Research 4, 043038 (2022).
  * Baker  _et al._ [2021] J. M. Baker, A. Litteken, C. Duckering, H. Hoffmann, H. Bernien, and F. T. Chong, Exploiting long-distance interactions and tolerating atom loss in neutral atom quantum architectures, in  _2021 ACM/IEEE 48th Annual International Symposium on Computer Architecture (ISCA)_ (IEEE, 2021) pp. 818–831.
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Arute  _et al._ [2020] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, A. Bengtsson, S. Boixo, M. Broughton, B. B. Buckley, _et al._ , Observation of separated dynamics of charge and spin in the fermi-hubbard model, arXiv preprint arXiv:2010.07965 (2020).
  * Jones  _et al._ [2022] E. B. Jones, L. E. Hillberry, M. T. Jones, M. Fasihi, P. Roushan, Z. Jiang, A. Ho, C. Neill, E. Ostby, P. Graf, _et al._ , Small-world complex network generation on a digital quantum processor, Nature communications 13, 4483 (2022).
  * Rotello  _et al._ [2023] C. Rotello, E. B. Jones, P. Graf, and E. Kapit, Automated detection of symmetry-protected subspaces in quantum simulations, Phys. Rev. Res. 5, 033082 (2023).
  * Hearth  _et al._ [2023] S. N. Hearth, M. O. Flynn, A. Chandran, and C. R. Laumann, Efficient local classical shadow tomography with number conservation, arXiv preprint arXiv:2311.09291 (2023).
  * King  _et al._ [2024] R. King, D. Gosset, R. Kothari, and R. Babbush, Triply efficient shadow tomography, arXiv preprint arXiv:2404.19211 (2024).
  * Manenti and Motta [2023] R. Manenti and M. Motta,  _Quantum information science_ (Oxford University Press, 2023).
  * McClean  _et al._ [2020] J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, _et al._ , Openfermion: the electronic structure package for quantum computers, Quantum Science and Technology 5, 034014 (2020).
  * Suzuki  _et al._ [2021] Y. Suzuki, Y. Kawase, Y. Masumura, Y. Hiraga, M. Nakadai, J. Chen, K. M. Nakanishi, K. Mitarai, R. Imai, S. Tamiya, _et al._ , Qulacs: a fast and versatile quantum circuit simulator for research purpose, Quantum 5, 559 (2021).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17, 261 (2020).
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, The finite group velocity of quantum spin systems, Communications in mathematical physics 28, 251 (1972).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-robinson bounds and the generation of correlations and topological quantum order, Physical review letters 97, 050401 (2006).
  * Jiang  _et al._ [2018] Z. Jiang, K. J. Sung, K. Kechedzhi, V. N. Smelyanskiy, and S. Boixo, Quantum algorithms to simulate many-body physics of correlated fermions, Physical Review Applied 9, 044036 (2018).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, _et al._ , Recent developments in the pyscf program package, The Journal of chemical physics 153 (2020).
  * Graham  _et al._ [2023] T. Graham, L. Phuttitarn, R. Chinnarasu, Y. Song, C. Poole, K. Jooya, J. Scott, A. Scott, P. Eichler, and M. Saffman, Mid-circuit measurements on a neutral atom quantum processor, arXiv preprint arXiv:2303.10051 (2023).
  * Norcia  _et al._ [2023] M. A. Norcia, W. B. Cairncross, K. Barnes, P. Battaglino, A. Brown, M. O. Brown, K. Cassella, C.-A. Chen, R. Coxe, D. Crow, J. Epstein, C. Griger, A. M. W. Jones, H. Kim, J. M. Kindem, J. King, S. S. Kondov, K. Kotru, J. Lauigan, M. Li, M. Lu, E. Megidish, J. Marjanovic, M. McDonald, T. Mittiga, J. A. Muniz, S. Narayanaswami, C. Nishiguchi, R. Notermans, T. Paule, K. A. Pawlak, L. S. Peng, A. Ryou, A. Smull, D. Stack, M. Stone, A. Sucich, M. Urbanek, R. J. M. van de Veerdonk, Z. Vendeiro, T. Wilkason, T.-Y. Wu, X. Xie, X. Zhang, and B. J. Bloom, Midcircuit qubit measurement and rearrangement in a ${}^{171}\mathrm{Yb}$ atomic array, Phys. Rev. X 13, 041034 (2023).
  * Acharya  _et al._ [2023a] R. Acharya _et al._ , Suppressing quantum errors by scaling a surface code logical qubit, Nature 614, 676 (2023a).
  * Acharya  _et al._ [2023b] S. Acharya, D. Pashov, C. Weber, M. van Schilfgaarde, A. I. Lichtenstein, and M. I. Katsnelson, A theory for colors of strongly correlated electronic systems, Nature Communications 14, 5565 (2023b).
  * Pedernales  _et al._ [2014] J. Pedernales, R. Di Candia, I. Egusquiza, J. Casanova, and E. Solano, Efficient quantum algorithm for computing n-time correlation functions, Physical Review Letters 113, 020505 (2014).
  * Terhal and DiVincenzo [2000] B. M. Terhal and D. P. DiVincenzo, Problem of equilibration and the computation of correlation functions on a quantum computer, Physical Review A 61, 022301 (2000).
  * Rohringer  _et al._ [2012] G. Rohringer, A. Valli, and A. Toschi, Local electronic correlation at the two-particle level, Phys. Rev. B 86, 125114 (2012).
  * Kaufmann  _et al._ [2017] J. Kaufmann, P. Gunacker, and K. Held, Continuous-time quantum monte carlo calculation of multiorbital vertex asymptotics, Phys. Rev. B 96, 035114 (2017).
  * Kim  _et al._ [2022] A. J. Kim, J. Li, M. Eckstein, and P. Werner, Pseudoparticle vertex solver for quantum impurity models, Phys. Rev. B 106, 085124 (2022).
  * Haah  _et al._ [2021] J. Haah, M. B. Hastings, R. Kothari, and G. H. Low, Quantum algorithm for simulating real time evolution of lattice hamiltonians, SIAM Journal on Computing 52, FOCS18 (2021).
  * Matthies  _et al._ [2022] A. Matthies, M. Rudner, A. Rosch, and E. Berg, Programmable adiabatic demagnetization for systems with trivial and topological excitations, arXiv preprint arXiv:2210.17256 (2022).
  * Mi  _et al._ [2024] X. Mi, A. Michailidis, S. Shabani, K. Miao, P. Klimov, J. Lloyd, E. Rosenberg, R. Acharya, I. Aleiner, T. Andersen, _et al._ , Stable quantum-correlated many-body states through engineered dissipation, Science 383, 1332 (2024).
  * Lloyd  _et al._ [2024] J. Lloyd, A. Michailidis, X. Mi, V. Smelyanskiy, and D. A. Abanin, Quasiparticle cooling algorithms for quantum many-body state preparation, arXiv preprint arXiv:2404.12175 (2024).
  * Sahay and Verresen [2024] R. Sahay and R. Verresen, Finite-depth preparation of tensor network states from measurement, arXiv preprint arXiv:2404.17087 (2024).
  * Smith  _et al._ [2024] K. C. Smith, A. Khan, B. K. Clark, S. Girvin, and T.-C. Wei, Constant-depth preparation of matrix product states with adaptive quantum circuits, arXiv preprint arXiv:2404.16083 (2024).