# Utilizing Quantum Processor for the Analysis of Strongly Correlated
Materials

Hengyue Li, Yusheng Yang, Pin Lv, Jinglong Qu, Zhe-Hui Wang, Jian Sun,
Shenggang Ying

###### Abstract

This study introduces a systematic approach for analyzing strongly correlated
systems by adapting the conventional quantum cluster method to a quantum
circuit model. We have developed a more concise formula for calculating the
cluster’s Green’s function, requiring only real-number computations on the
quantum circuit instead of complex ones. This approach is inherently more
suited to quantum circuits, which primarily yield statistical probabilities.
As an illustrative example, we explored the Hubbard model on a 2D lattice. The
ground state is determined utilizing Xiaohong, a superconducting quantum
processor equipped with 66 qubits, supplied by QuantumCTek Co., Ltd.
Subsequently, we employed the circuit model to compute the real-time retarded
Green’s function for the cluster, which is then used to determine the lattice
Green’s function. We conducted an examination of the band structure in the
insulator phase of the lattice system. This preliminary investigation lays the
groundwork for exploring a wealth of innovative physics within the field of
condensed matter physics.

##  1 Introduction

Quantum computing garners interest for its potential to solve complex problems
beyond classical computing’s reach. Pioneering algorithms like Shor’s
algorithm [1, 2] and Grover’s algorithm [3] have marked significant
milestones. Additionally, emerging algorithms such as the HHL algorithm [4]
are in development, expanding the quantum computing repertoire. At the same
time, quantum processor hardware is also undergoing significant development.
Unlike classical computers, where manufacturing techniques have already
converged on electronic circuits based on silicon semiconductor materials,
quantum computing presents a unique challenge. Various technical approaches
are being explored for quantum computer realization, with superconducting (SC)
systems emerging as a frontrunner. Over the recent decades, there has been
notable progress not only in the number of qubits but also in their quality in
SC systems [5]. However, achieving practical fault-tolerant quantum computing
(FTQC) compatible with standard quantum algorithms requires logical qubits of
exceptionally high quality. Among various error-correcting codes, the surface
code [6] has shown promise due to its high error threshold. Although some
demonstration experiments [7, 8] have shown a modest improvement in the
quality of qubits encoded using surface code, these advancements are still
insufficient for the implementation of large-scale quantum computing systems.

A pragmatic strategy in the current stage of quantum computing is the
development of noisy intermediate-scale quantum (NISQ) applications. Recent
advancements in this area include the quantum variational method, like
variational quantum eigensolver (VQE) [9, 10] and quantum machine learning
algorithms [11, 12]. Representative algorithms of NISQ era, such as VQE, have
found applications in diverse fields, including quantum chemistry [9, 13, 14]
and high-energy physics [15]. It represents a hybrid quantum-classical
approach, utilizing quantum circuits to produce states that approximate the
ground state of a Hamiltonian, and employing a classical optimizer to optimize
across these circuits. However, applications of NISQ era in other research
fields, such as material science, are comparatively less explored.

Materials exhibiting strong correlation effects, such as heavy-fermion systems
[16] and non-Fermi-liquid systems [17], have long fascinated researchers. In
recent decades, high-temperature superconductors, perhaps one of the most
attractive materials, have also started to draw a lot of attention, promising
to change the way we use technology. Other interesting phenomena from strongly
correlated materials, like Mott metal-insulator transitions [18], fractional
quantum Hall effect [19], and spin-state transitions [20, 21], inspire the
pursuit of new metamaterials. All these materials don’t behave as expected
according to traditional band theory [22]. The excitement around these
materials comes from discovering unusual behaviors, such as novel quantum
phase transitions, and unexpected electric and magnetic properties. However,
simulating strongly correlated materials presents significant challenges at
the present time because of the complex interplay of electron-electron
interactions within these materials. Current simulation methods, predominantly
first-principles techniques like density functional theory (DFT), face hurdles
with systems containing strongly correlated electrons. Concurrently, the
development of non-perturbation theories has offered alternative perspectives
in quantum physics to simulate strongly correlated materials. Dynamic mean
field theory (DMFT) [23, 24] has notably achieved early success by
innovatively mapping the orbitals in a lattice system to a local impurity
model, based on the assumption that the lattice self-energy is local. To fully
consider the local correlations, cluster perturbation theory (CPT) [25] has
subsequently been developed. Additionally, considering broken-symmetry states,
other cluster methods with a meanfield approach have been developed [26, 27,
28]. Collectively, these approaches for considering local correlations are
called quantum cluster methods (QCM) [29].

Intuitively, leveraging NISQ circuits for quantum simulation of strongly
correlated systems appears appropriate and potentially sufficient [30].
Various encoding methods [31, 32] enable the mapping of finite-sized fermion
clusters onto a circuit model. The core of QCM involves iteratively
calculating the Green’s function of a cluster system. Several methods [33, 34]
for computing Green’s functions on a circuit model have been developed, making
it theoretically possible to implement QCM calculations on a quantum circuit
model. It is the purpose of this paper to introduce a systematic procedure for
implementing each critical step of QCM on the circuit model. As an
illustrative example, we investigate the paramagnetic insulator phase of a
standard Hubbard model on a 2-dimensional (2D) square lattice. The results
obtained solely through the circuit model will be compared with exact results
obtained from classical simulations.

In the remainder of this paper, we organize the content as follows. Standard
formulas and steps for QCM are revisited in Sec. 2. The approach to
determining the ground state via VQE is presented in Sec. 3.1. Section 3.2
introduces a streamlined process for calculating the real-time retarded
Green’s function. Section 4 briefly summarizes the entire calculation
procedure. The results of the ground state energy, computed on the quantum
processor Xiaohong provided by QuantumCTek, are discussed in Sec. 5.1. The
findings regarding the cluster Green’s function are presented and analyzed in
Sec. 5.2. Using QCM, the one-particle excitation spectra are explored in Sec.
5.3. Finally, Sec. 6 offers a concise summary of our findings.

##  2  Reflective overview of the quantum cluster method

We revisit fundamental theories of QCM for strongly correlated systems. A
general description of a fermion system can be expressed through a
Hamiltonian, as given by

$\displaystyle H$
$\displaystyle=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}+\sum_{ijkl}X_{ijkl}c_{i}^{\dagger}c_{j}^{\dagger}c_{k}c_{l},$
(1)

where the operator $c_{i}^{\dagger}$ ($c_{j}$ ) creates (annihilates) an
electron at site $i$ ($j$). Without loss of generality, we suppress the index
of spin (and also other possible orbital degrees of freedom). The hopping
energy, denoted as $t_{ij}$, and the interacting energy, $X_{ijkl}$ , as
outlined in Eq. 1, are typically computed using ab initio methods with proper
basis sets. To illustrate, the interacting energy $X_{ijkl}$ is defined by the
integral

$\displaystyle X_{ijkl}$ $\displaystyle=\int\int
d\boldsymbol{r}d\boldsymbol{r}^{\prime}\phi_{i}^{*}(\boldsymbol{r})\phi_{j}^{*}(\boldsymbol{r}^{\prime})V(\boldsymbol{r}-\boldsymbol{r}^{\prime})\phi_{k}(\boldsymbol{r}^{\prime})\phi_{l}(\boldsymbol{r}),$
(2)

where $\\{\phi_{i}(\boldsymbol{r})\\}$ represents a complete set of states and
$V(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ is the electron-electron
interacting energy [35]. This model is versatile and can be adapted to
represent many well-known models. For instance, by focusing solely on the
interacting term $X_{ijkl}\delta_{ik}\delta_{jl}$ and assigning $i$ and $j$ to
different spin indices on the same site, Eq. 1 transforms into the standard
Hubbard model [36, 37], described as:

$\displaystyle H$
$\displaystyle=-\gamma\sum_{ij\sigma}\left(c_{i\sigma}^{\dagger}c_{j\sigma}+h.c.\right)+Un_{i\uparrow}n_{i\downarrow}-\mu\sum_{i\sigma}n_{i\sigma},$
(3)

where $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ is the particle number
for spin $\sigma$ on the $i$-th site, $\mu$ is the chemical potential,
$-\gamma$ is the hopping integral, and $U$ represents the interaction term of
electron repulsion.

In QCM, the original lattice is tiled through super-lattice expanded by
$\\{\boldsymbol{e}_{i}\\}$, with super-lattice-cell of size $L$, as shown in
Fig. 1. For each operator $c_{i}$, a distinct cluster index $\boldsymbol{R}$
and a corresponding in-site index $\alpha$ are identified. This allows for the
operator to be re-expressed in a more refined form, transforming

$\displaystyle c_{i}$ $\displaystyle\to\psi_{\boldsymbol{R}\alpha}.$ (4)

With this notation, Eq. 1 is decoupled by tiling clusters and can be rewritten
as
$H=\sum_{\boldsymbol{R}}H_{c}\left(\\{\psi_{\boldsymbol{R}\alpha}\\}\right)+\sum_{\boldsymbol{R}\boldsymbol{R}^{\prime}}T(\boldsymbol{R}-\boldsymbol{R}^{\prime})$.
Here, $H_{c}\left(\\{\psi_{\boldsymbol{R}\alpha}\\}\right)$ represents a
cluster-specific Hamiltonian, while
$T(\boldsymbol{R}-\boldsymbol{R}^{\prime})$ denotes the Hamiltonian terms that
facilitate connections between the cluster $\boldsymbol{R}$ and
$\boldsymbol{R}^{\prime}$. Each cluster can be considered as a reference
system for abstracting the self-energy [25]. The lattice Green’s function,
$\boldsymbol{\mathcal{G}}(\boldsymbol{q},\omega)$, is given by

$\displaystyle\boldsymbol{\mathcal{G}}(\boldsymbol{q},\omega)$
$\displaystyle=\left[\boldsymbol{G}^{-1}(\omega)-\tau_{\boldsymbol{q}}\right]^{-1},$
(5)

where $\boldsymbol{q}$ is the reciprocal vector in the reduced Brillouin zone,
$\boldsymbol{G}(\omega)$ (a $L\times L$ matrix) is the Green’s function of a
local cluster, and $\tau_{\boldsymbol{q}}$ is a linear combination of
$T(\boldsymbol{0}-\boldsymbol{r})$ around a center cluster $\boldsymbol{0}.$
Details of these derivations can be found in Appendix Appendix: Details of
QCM.

We recall some basic knowledge of Green’s function theory for calculating each
matrix element $G_{ij}$ of $\boldsymbol{G}$. The definition of the real-time
retarded Green’s function [38] is given by

$\displaystyle G_{ij}^{R}(t-t^{\prime})$
$\displaystyle=-i\Theta(t-t^{\prime})\langle\\{c_{i}(t),c_{j}^{\dagger}(t^{\prime})\\}\rangle$
(6)

where $\Theta$ is the Heaviside step function, the curly bracket $\\{*,*\\}$
represents the anti-commutator, $c_{i}(t)=e^{iHt}c_{i}e^{-iHt}$ is the
Heisenberg representation of the operator $c_{i}$, the bracket
$\langle...\rangle=\text{Tr}\left(e^{-H/T}...\right)/\text{Tr}\left(e^{-H/T}\right)$
represents the thermal average. Throughout this paper, we set $\hbar=1$. In
zero temperature case $(T=0)$, the calculation of the thermal average value is
specialized by $\langle g|...|g\rangle$, where $|g\rangle$ is the ground state
of the system (assuming a non-degenerate scenario). With relation
$\Theta(t-t^{\prime})=\frac{-1}{2\pi
i}\int_{-\infty}^{\infty}d\omega\frac{e^{-i(t-t^{\prime})\omega}}{\omega+i\eta},$
where $\eta\to 0^{+}$ is an infinite small real number, we rewrite the real-
time retarded Green’s function as

$\displaystyle G_{ij}^{R}(t)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega
G_{ij}(\omega+i\eta)e^{-i\omega t},$ (7)

where $G_{ij}(\omega+i\eta)$ is the Fourier transform of the real-time
retarded Green’s function. With analytic continuation to full complex plane
$\omega+i\eta\to z$ and Kramers–Kronig relations [39], the Green’s function is
given by

$\displaystyle G_{ij}(z)$ $\displaystyle=\int
dz^{\prime}\frac{\rho_{ij}(z^{\prime})}{z-z^{\prime}},$ (8)

where

$\displaystyle\rho_{ij}(\omega)$
$\displaystyle=-\frac{1}{\pi}\mathcal{I}\left[G_{ij}(\omega+i\eta)\right]$ (9)

is the spectral function. Many useful Green’s functions can be derived from
this form. For instance, the Matsubara Green’s function can be calculated by
$G_{ij}(i\omega_{n})$, where $\omega_{n}$ is the Matsubara frequency [35]. In
Lehmann representation, at zero temperature case the spectral function is
given by

$\displaystyle\rho_{ij}(\omega)$ $\displaystyle=\sum_{n}\langle
g|c_{i}|n\rangle\langle n|c_{j}^{\dagger}|g\rangle\delta\left(\omega-
E_{n}+E_{0}\right)$ $\displaystyle+\sum_{n}\langle n|c_{i}|g\rangle\langle
g|c_{j}^{\dagger}|n\rangle\delta\left(\omega-E_{0}+E_{n}\right),$ (10)

where $\\{|n\rangle\\}$ and $\\{E_{n}\\}$ are eigenstates and eigenvalues of
the system respectively. When $i=j$, $\rho_{i}(\omega)\equiv\rho_{ij}(\omega)$
represents the density of states at orbital $i$. It obeys the summation rule
$\int\rho_{i}(\omega)=1$. Practically, we calculate the spectral function by
Eq. (9). And the Green’s function in frequency domain, $G_{ij}(\omega+i\eta),$
is calculated by the inverse of Eq. (7) and given by

$\displaystyle G_{ij}(\omega+i\eta)$
$\displaystyle=\int_{0}^{\infty}dte^{i\omega t}\left[e^{-\eta
t}G_{ij}^{R}(t)\right].$ (11)

In actual computations, we rely on an artificial broadening parameter $\eta$
to facilitate the calculation, instead of employing an impractical value of
$0^{+}$. Henceforth, we disregard the subscript “$R$” in $G_{ij}^{R}(t)$.

Figure 1:  The original lattice is divided into distinct clusters and inter-
cluster segments. Each cluster (red areas), referred to as a reference system,
is independently analyzed using various methods, such as exact diagonalization
(ED). The inter-cluster section is treated as a perturbation affecting the
reference system. To address symmetry breaking, we introduce additional mean-
fields. Meanwhile, achieving the metallic phase requires integrating
additional environmental elements, typically bath sites, which are hybridized
with the cluster. These additional parameters, mean-fields, and hybridization
functions are determined through a self-consistent variational process.

##  3 Calculation on circuits

###  3.1  Ground state

The QCM effectively segments the original lattice into distinct clusters.
These clusters, each characterized by a specific size labeled $L$, are
regarded as individual cluster systems. The VQE method is employed to
investigate the ground state of the system using quantum processor. We
employed a parameterized unitary transformation, denoted as
$U_{G}(\boldsymbol{\theta})$, to construct the parameterized state from
scratch as
$|g(\boldsymbol{\theta})\rangle=U_{G}(\boldsymbol{\theta})|\boldsymbol{0}\rangle.$
In VQE method, the expectation value of the Hamiltonian serves as an upper
limit for the ground state energy. Thus, the ground state energy $E_{0}$ is
bounded by $E_{0}\leq\langle
g(\boldsymbol{\theta})|H_{c}|g(\boldsymbol{\theta})\rangle$, where $H_{c}$ is
the Hamiltonian of the system.

In order to performed VQE by quantum circuits, it is essential to transform
the Hamiltonian into a format compatible with quantum processor measurements.
We adopt the standard Jordan-Wigner (JW) encoding [31]:

$\displaystyle c_{i}$
$\displaystyle=\frac{1}{2}\left(X_{i}-iY_{i}\right)\prod_{\alpha<i}\left(-Z_{\alpha}\right),$
(12)

where $\sigma_{\alpha}$ represents the $\sigma$ gate ($\sigma=I,X,Y$ and $Z$)
on the $\alpha$-th qubit. With this encoding, the Hamiltonian of the cluster
$H_{c}\left(\\{c_{i}\\}\right)$ is decomposed by

$\displaystyle H_{c}$ $\displaystyle=\sum_{i}\xi_{i}P_{i},$ (13)

where the coefficient $\xi_{i}$ is a real number and
$P_{i}=\sigma_{0}\otimes\sigma_{1}\otimes...$ is the tensor product of Pauli
matrices. Hence the energy of the system is parameterized as
$E(\boldsymbol{\theta})=\sum_{i}\xi_{i}\langle P_{i}\rangle_{\theta}$, where
$\langle P_{i}\rangle_{\theta}=\langle
g(\boldsymbol{\theta})|P_{i}|g(\boldsymbol{\theta})\rangle.$ The ground state
can be ascertained by minimizing $E(\boldsymbol{\theta})$. The optimal point,
$\boldsymbol{\theta}_{0}$, yields the minimal value
$E(\boldsymbol{\theta}_{0})$. Consequently, the approximated ground state is
represented as
$|\tilde{g}\rangle=U_{G}(\boldsymbol{\theta}=\boldsymbol{\theta}_{0})|\boldsymbol{0}\rangle$.

###  3.2 Green’s function

The Green’s function can be calculated in both time and frequency domains. In
the time domain, one can compute the real-time retarded Green’s function [40,
41] or imaginary-time Green’s function [42]. To obtain the spectral function,
these methods require additional processes, such as Fourier transform or
numerical analytic continuation. Based on Lehmann representation [43], the
Green’s function can be calculated directly in the frequency domain using the
variational quantum eigensolver [44] to find excitation state energies. In
this work, we calculate the real-time retarded Green’s function, as we believe
simulating a time-evolution problem on a quantum circuit is more natural.

The calculation of a retarded Green’s function, based on the given ground
state $|g\rangle$, is defined by Eq. 6. In the previous calculation,
additional parameter controls were needed to consider the real and imaginary
parts of each component of the Green’s function separately [40]. In this
context, we have restructured the formulas so that only the computation of
real values is necessary. Using JW encoding, Eq. 6 is reformulated as:

$\displaystyle G_{ij}(t)$
$\displaystyle=\frac{\left(-1\right)^{i+j}}{4}\left(W_{ij}^{-}-iW_{ij}^{+}\right),$
(14)

where $W_{ij}^{-}$ and $W_{ij}^{+}$ represent the real and imaginary parts,
respectively, and

$\displaystyle W_{ij}^{\pm}$
$\displaystyle=\mathcal{F}\left(\bar{X_{i}},\bar{Y_{j}}\right)\pm\mathcal{F}\left(\bar{Y_{i}},\bar{X_{j}}\right).$
(15)

In above expression, we have defined
$\bar{\sigma_{i}}=\prod_{\alpha<i}Z_{\alpha}\sigma_{i}$. And
$\mathcal{F}(\bar{\sigma_{i}},\bar{\sigma_{j}})$ is given by

$\displaystyle\mathcal{F}\left(\bar{\sigma_{i}},\bar{\sigma_{j}}\right)=$
$\displaystyle\langle
g|e^{iHt}\bar{\sigma_{i}}e^{-iHt}\bar{\sigma_{j}}|g\rangle$
$\displaystyle+\langle
g|\bar{\sigma_{j}}e^{iHt}\bar{\sigma_{i}}e^{-iHt}|g\rangle.$ (16)

The advantage of basing all calculations on Eq. 16 is that
$\mathcal{F}\left(\bar{\sigma_{i}},\bar{\sigma_{j}}\right)$ is a real number,
and Eq. 16 can be directly computed using the circuit shown in Fig. 2. This
circuit, which includes parameters $i,j,$$t$ and $N_{\tau}$, is repeatedly
called as a subroutine by the program. At the end, the attached qubit on the
top line in Fig. 2 is measured. The probability of measuring 0, denoted as
$p_{+}$, is obtained through sampling from the quantum circuit. The final
value of Eq. 16, is given by
$\mathcal{F}\left(\bar{\sigma_{i}},\bar{\sigma_{j}}\right)=2(2p_{+}-1)$, with
a detailed derivation available in Appendix Appendix: Circuit for
$\mathcal{F}\left(P_{i},P_{j}\right)$.

In this approach, to simulate the time-evolution of $U_{t}=e^{-iHt}$, usually
Suzuki-trotter [45] decomposition is used. It is given by

$\displaystyle e^{-iHt}$
$\displaystyle\approx\underbrace{e^{-iH\tau}...e^{-iH\tau}}_{N_{\tau}}.$ (17)

To achieve the necessary precision, it is imperative to use rather small time
slices $\tau=t/N_{\tau}$, which results in a very long circuit depth. However,
the current state of quantum hardware, with its noise levels, poses a
challenge in simulating such long circuits with the required accuracy. A
potential workaround is to employ variational principles, such as McLachlan’s
[46, 47], to approximate the original decomposition circuit with a shorter one
[40, 41]. Nevertheless, this introduces computational overhead, as fitting an
ansatz circuit to replicate the time evolution of the Hamiltonian requires
solving an additional system of differential equations [47]. These challenges
raise concerns about the scalability of this method as the system size
increases. In this work, to determine the noise threshold necessary for the
direct application of the Suzuki decomposition, we persist in employing Eq. 17
to calculate the Green’s function on a noisy simulator with a comprehensive
depolarization error [48] rate of 0.01%, affecting each gate and channel. To
maintain consistency in the subsequent calculations, we utilize an exact
ground state $|g\rangle$ as input, despite having previously obtained a ground
state $|\tilde{g}\rangle$ through VQE on quantum hardware.

Figure 2:  Following the fundamental principles of the Hadamard test [49], the
circuit is devised to compute $\mathcal{F}\left(P_{i},P_{j}\right)$ as
outlined in Eq. 16. The topmost qubit serves as an ancilla qubit and is
configured for the final measurement. The result of this measurement, denoted
as $p_{+}$ , determines the value of Eq. 16.

##  4  Computation procedures

All of our numerical calculations are based on isQ [50]. In the preceding
sections, we delved into methodologies for determining the ground state and
computing the Green’s function. We summarize below the standard steps for the
whole calculation:

  1. 1.

Initialization of the ansatz: Set the ansatz $U_{G}(\boldsymbol{\theta})$, and
based on the variational principle by calling quantum hardware iteratively, a
set of optimized parameters $\boldsymbol{\theta}_{0}$ (therefore the optimized
ground state
$|\tilde{g}\rangle=U_{G}(\boldsymbol{\theta}_{0})|\boldsymbol{0}\rangle$) is
found.

  2. 2.

Calculation of the retarded Green’s function: For the specified input
$|g\rangle$, we employ the circuit illustrated in Fig. 2 to compute Eq. 16.
Subsequently, this allows us to calculate the real-time retarded Green’s
function as defined in Eq. 14.

  3. 3.

Fourier transformation for frequency domain: The cluster Green’s function in
the frequency domain, $G_{ij}(\omega)$, is derived via a Fourier transform, as
outlined in Eq. 11. We adopt the Gauss–Legendre quadrature method [51] for
implementing an exceptionally precise and efficient integration.

  4. 4.

Computation of the lattice Green’s function: Following the acquisition of each
matrix element $G_{ij}$ of $\boldsymbol{G}$, proceed to calculate the lattice
Green’s function $\boldsymbol{\mathcal{G}}$ by applying Eq. 5.

##  5 Results and discussions

###  5.1  Ground state energy

Figure 3:  The ground state energies as a function of interaction $U$. Inset
details illustrate the mitigation within Hamiltonian $H_{c}$ at $U=3$, showing
the contribution of each term.

For clarity and simplicity, we omit the consideration of symmetry-breaking
cases, thereby reverting the calculation to a standard CPT calculation.
Meanwhile, since the reference system lacks bath sites, our study is limited
to the insulator phase. Despite these simplifications, we still possess
sufficient computational resources to consider the anti-ferromagnetic phase
[52], having selected a $2\times 1$ cluster as the reference system, as
illustrated in Fig. 1.

With consideration of the spin degree of freedom, the quantum state of the
cluster is encoded by four qubits. In our calculations, we set the hopping
energy $\gamma=1$ as the energy unit, and $\mu=\frac{U}{2}$ to maintain the
half-filling condition. In our case, in pursuit of high efficiency, we employ
a well-designed ansatz that consists solely of a single variational parameter,
denoted as $\varphi$ . The variational process is further simplified by
directly minimizing analytical sinusoidal function which is obtained by
fitting sparse sampling data from scanning $\varphi$ [15]. More details of the
variational ansatz can be found in Appendix Appendix: Ansatz.

The final total energies of the cluster system as a function of $U$ are shown
in Fig. 3. The raw data of energy (red solid line), which is calculated from
the original circuit on the quantum hardware without mitigation, decreases
with increasing U. The exact result by classical simulation is calculated and
shown as the green dashed line in Fig. 3. We see that the raw data is
consistent with the exact result quantitatively despite some discrepancies. We
recorrect this result by post-condition of digital zero noise extrapolation
(DZNE) with circuit folding [53]. For the $U$>0 case, there are 6 terms, in
addition to a constant term with $P=I$, in the Hamiltonian, see Eq. 13. The
result of each term at $U=3$ are shown in the insertion of Fig. 3, from which
we see that these results have been improved to varying degrees. The corrected
energy is shown as the black solid line in Fig. 3. Compared with the raw data,
the result is enhanced significantly.

###  5.2 Green’s function and spectral function

Figure 4:  Comparison of retarded Green’s function calculated by a quantum
circuit (cyan solid line) and the exact result (orange dashed line). (a) The
imaginary part of the real-time retarded Green’s function. (b) The Green’s
function after being multiplied by a factor $e^{-\eta t}$ which effectively
reduces the error at large $t$.

In our subsequent analyses, we maintain $U=3$. The computation of the real-
time retarded Green’s function necessitates a delicate balance in the Suzuki
approximation’s precision. Specifically, the time slice $\tau$ in the
decomposition of $e^{-iHt}$ (refer to Eq. 17) must be carefully calibrated. If
$\tau$ is excessively large, precision is compromised, while if it is too
small, the cumulative errors from a circuit with long depth would overshadow
any meaningful data obtained from the results. To address this, we
standardized $\tau$ with $N_{\tau}=60$ for all simulations.

As an illustrative case, we computed the on-site Green’s function at the first
site $(i=j=1)$ for a single spin. Due to symmetry, the results for spin up and
spin down are the same. Given the particle-hole symmetry in the system, the
real part of this Green’s function vanishes, and we have exclusively
calculated and presented the imaginary part, $\text{Im}[G(t)]$, represented by
the solid cyan lines in Fig. 4. For comparative analysis, we also calculated
the exact result, represented as orange dashed lines in the same figure. As
shown in Fig. 4 (a), in scenarios with small $t$, the circuit-derived results
closely align with the exact results, maintaining satisfactory precision. This
alignment is attributed to the smaller value of $\tau$ (due to the fixed value
of $N_{\tau}$), ensuring the accuracy of the Suzuki decomposition. Minor
discrepancies in the small $t$ region primarily stem from circuit noise.
However, as $t$ increases, $\tau$ consequently becomes larger, leading to a
decline in the accuracy of the Suzuki approximation. This trend is observable
in Fig. 4 (a), where beyond $t\approx 22$, the circuit results significantly
diverge from the exact ones. Despite this, we demonstrate that such large
discrepancies at higher $t$ values do not substantially affect subsequent
calculations. Recall Eq. 11, where the integrand includes the factor $e^{-\eta
t}$. For practical computations, we have set $\eta=0.2$. The exponential decay
induced by $e^{-\eta t}$ effectively mitigates the impact of discrepancies
resulting from larger $\tau$ values at higher $t$ regions, as illustrated in
Fig. 4 (b). This property can serve as a valuable guide for devising more
refined computational strategies, allowing us to balance precision with the
allocation of computational resources effectively.

Figure 5:  The quantum circuit calculations (cyan solid line) and exact
results (orange dashed line). (a) The real part of the Green’s function in the
frequency domain. (b) The corresponding spectral function of the cluster.

To perform the integration specified in Eq. 11, we employ a total of $100$
sampling points and set the upper limit to $t_{\infty}=30$. We have
ascertained the sufficiency of this integration range through our prior
discussion, as illustrated in Fig. 4 (b). The real part of $G(\omega+i\eta)$
and the corresponding spectral function are presented in Fig. 5 (a) and (b),
respectively. We depict the results obtained through the circuit model as
solid cyan lines and compare them with the exact results represented by dashed
orange lines. It is noteworthy that, despite the presence of notable
discrepancies arising from circuit noise in the calculation of the real-time
retarded Green’s function, the ultimate outcome—$G$ in the frequency domain—is
within reasonable accuracy. The positions of the poles of the Green’s function
remain accurate, albeit with a decrease in amplitudes attributable to the
effects of noisy errors. Our examination of the summation rule reveals that
$\int d\omega\rho(\omega)\sim 0.88$, indicating an approximate $12\%$ loss.
This diminution is also visually evident in Fig. 5 (b), where the inner peaks
appear marginally reduced. This loss predominantly stems from quantum noise.

###  5.3 One-particle excitation spectra

Figure 6:  Intensity plot of one-particle excitation spectra for the lattice
system revealing the distribution and intensity of excitations across various
energy levels and momentum states.

Once we obtain the cluster Green’s function, we can proceed with the
calculation of physical properties relevant to the lattice model. The lattice
Green’s function in the original Brillouin zone is defined as

$\displaystyle g_{\sigma}(\boldsymbol{k},\omega)$
$\displaystyle=\frac{1}{L}\sum_{ij}e^{-i\boldsymbol{k}\cdot(\boldsymbol{r}_{i}-\boldsymbol{r}_{j})}\mathcal{G}_{i\sigma,j\sigma}(\boldsymbol{k},\omega),$
(18)

where $\sigma=\uparrow$ or $\downarrow$ represents the spin index, and $i,j$
span all sites within the cluster. The definition of the single-particle
excitation spectra $\rho_{\sigma}(\boldsymbol{k},\omega)$ remains the same as
Eq. 9. Due to the symmetry, we calculate
$\rho_{\uparrow}(\boldsymbol{k},\omega)$, and the results are shown in Fig. 6.
The obtained data captures information around the typical critical points,
$\Gamma$, X, and M, which are consistent with the classical simulation of a
Hubbard model on a 2D square lattice by conventional methods. A key
characteristic of this model, the Mott-gap, is evident around the Fermion
surface due to interaction $U$, as depicted in Fig. 6. Despite some loss in
spectral function weight, clear band structures are still observable. This
analysis of one-particle excitation spectra paves the way for computing other
physical quantities, such as optical conductivity [54]. Such calculations also
hold significant practical relevance in experimental contexts, as this spectra
can be extracted through techniques like angle-resolved photoemission
spectroscopy (ARPES).

##  6  Conclusion

In summary, we have thoroughly demonstrated the feasibility of studying
strongly correlated systems using quantum circuits. To account for local
correlations using the QCM, we segment the lattice model into reference
clusters to abstract self-energy. Employing the VQE with an appropriate ansatz
implemented on quantum hardware, we successfully obtained the ground state of
a reference cluster, further refining the results through the error mitigation
of DZNE method. Formulas for calculating Green’s functions within the circuit
model were derived. We observed that the real-time retarded Green’s function
results are consistent with exact results within an acceptable error in the
small $t$ region, and the influence of any discrepancies in the large $t$
region on calculating the lattice Green’s function has been suppressed.
Utilizing the lattice Green’s function, we investigated the one-particle
excitation spectra. Our computations successfully reproduced the typical
properties of the band structure of a Hubbard model on a 2D square lattice.

##  7 Author infomation

### Corresponding authors

#### Hengyue Li

Arclight Quantum Computing Inc., Beijing 100191, China; Email:
<EMAIL_ADDRESS>
#### Yusheng Yang

Arclight Quantum Computing Inc., Beijing 100191, China; Email:
<EMAIL_ADDRESS>
#### Shenggang Ying

Arclight Quantum Computing Inc., Beijing 100191, China; Institute of Software,
Chinese Academy of Sciences, Beijing 100190, China; Email<EMAIL_ADDRESS>
### Authors

#### Pin Lv

China Telecom Quantum Information Technology Group Co., Ltd., Hefei 230031,
China

#### Jinglong Qu

High-Temperature Materials Institute, Central Iron and Steel Research
Institute, Beijing, 100081, China; Beijing GAONA Materials & Technology Co.,
Ltd., Beijing, 100081, China; Sichuan CISRI-Gaona Forging Co., Ltd., Deyang,
Sichuan 618099 China

#### Zhe-Hui Wang

QuantumCTek (Shanghai) Co., Ltd., Shanghai 200120, China

#### Jian Sun

Arclight Quantum Computing Inc., Beijing 100191, China

##  8 Acknowledgments

We express our sincere gratitude to the quantum cloud platform of QuantumCTek
Co., Ltd. (https://quantumctek-cloud.com) for their invaluable support in
providing quantum computing resources, which were instrumental in conducting
the computational experiments for this study. Without access to their quantum
hardware infrastructure, this work would not have been feasible. Additionally,
we extend our heartfelt appreciation to Guolong Cui from Arclight Quantum
Computing Inc. for his guidance and support throughout the research process.

## Appendix: Details of QCM

Utilizing the index transformation outlined in Eq. 4, we can reformulate the
original Hamiltonian (Eq. 1) as

$\displaystyle H$
$\displaystyle=\sum_{\boldsymbol{R}\alpha\boldsymbol{R}^{\prime}\beta}\mathcal{T}_{\alpha,\beta}^{\boldsymbol{R}-\boldsymbol{R}^{\prime}}\psi_{\boldsymbol{R}\alpha}^{\dagger}\psi_{\boldsymbol{R}^{\prime}\beta}$
$\displaystyle+\sum_{\boldsymbol{R}}\sum_{\alpha\beta\gamma\delta}X_{\boldsymbol{R}\alpha,\boldsymbol{R}\beta,\boldsymbol{R}\gamma,\boldsymbol{R}\delta}\psi_{\boldsymbol{R}\alpha}^{\dagger}\psi_{\boldsymbol{R}\beta}^{\dagger}\psi_{\boldsymbol{R}\gamma}\psi_{\boldsymbol{R}\delta}.$
(19)

In Eq. 19, $\mathcal{T}^{\boldsymbol{R}-\boldsymbol{R}^{\prime}}$ is defined
as a $L\times L$ matrix, where $L$ denotes the site of the cluster. We can
write
$t_{\boldsymbol{R}\alpha,\boldsymbol{R}^{\prime}\beta}=\mathcal{T}_{\alpha,\beta}^{\boldsymbol{R},\boldsymbol{R}^{\prime}}$
as $\mathcal{T}_{\alpha,\beta}^{\boldsymbol{R}-\boldsymbol{R}^{\prime}}$
because of the translation symmetry. We assume all interactions to be local
within clusters. If inter-cluster interactions exist, clusters need to be
decoupled using a mean-field approximation. The inter-site Coulomb interaction
$Vn_{i}n_{j}$, for instance, can be approximated at the Hartree-Fock level as
: $Vn_{i}n_{j}\approx V(n_{i}\langle n_{j}\rangle+\langle n_{i}\rangle
n_{j}-\langle n_{i}\rangle\langle n_{j}\rangle)$. The fluctuation term
$(n_{i}-\langle n_{i}\rangle)(n_{j}-\langle n_{j}\rangle)$ is neglected in
this approximation, introducing additional mean fields $\\{\langle
n_{i}\rangle\\}$ into the expression. These mean fields can be determined
either through iterative calculation or by minimizing the self-energy
functional $\Omega_{t}[\Sigma]$ [28, 54, 29]. Translation symmetry implies
that the Hamiltonian parameters for each cluster $\boldsymbol{R}$ are
identical, leading to
$X_{\boldsymbol{R}\alpha,\boldsymbol{R}\beta,\boldsymbol{R}\gamma,\boldsymbol{R}\delta}=\chi_{\alpha\beta\gamma\delta}$.

Further, the first term in Eq. 19 is divided into inter-cluster and intra-
cluster parts. By incorporating the identity
$1=\delta_{\boldsymbol{R}\boldsymbol{R^{\prime}}}+(1-\delta_{\boldsymbol{R}\boldsymbol{R^{\prime}}})$
into the first term of Eq. 19, we obtain:
$\sum_{\boldsymbol{R}}\sum_{\alpha\beta}\psi_{\boldsymbol{R}\alpha}^{\dagger}\mathcal{T}_{\alpha,\beta}^{0}\psi_{\boldsymbol{R}\beta}+\sum_{\boldsymbol{R}\neq\boldsymbol{R}^{\prime}}\sum_{\alpha\beta}\psi_{\boldsymbol{R}\alpha}^{\dagger}\mathcal{T}_{\alpha,\beta}^{\boldsymbol{R}-\boldsymbol{R}^{\prime}}\psi_{\boldsymbol{R}^{\prime}\beta}$.
Therefore, the Hamiltonian (Eq. 19) is decoupled by clusters and can be
rewritten as:

$\displaystyle H$
$\displaystyle=\sum_{\boldsymbol{R}}H_{c}\left(\\{\psi_{\boldsymbol{R}\alpha}\\}\right)+\sum_{\boldsymbol{R}\boldsymbol{R}^{\prime}}T(\boldsymbol{R}-\boldsymbol{R}^{\prime}),$
(20)

where

$\displaystyle H_{c}\left(\\{\psi_{\boldsymbol{R}\alpha}\\}\right)$
$\displaystyle=\sum_{\alpha\beta}\psi_{\boldsymbol{R}\alpha}^{\dagger}\mathcal{T}_{\alpha,\beta}^{0}\psi_{\boldsymbol{R}\beta}$
$\displaystyle+\sum_{\alpha\beta\gamma\delta}\chi_{\alpha\beta\gamma\delta}\psi_{\boldsymbol{R}\alpha}^{\dagger}\psi_{\boldsymbol{R}\beta}^{\dagger}\psi_{\boldsymbol{R}\gamma}\psi_{\boldsymbol{R}\delta},$
(21)

and

$\displaystyle T(\boldsymbol{R}-\boldsymbol{R}^{\prime})$
$\displaystyle=\sum_{\alpha\beta}\psi_{\boldsymbol{R}\alpha}^{\dagger}\mathcal{T}_{\alpha,\beta}^{\boldsymbol{R}-\boldsymbol{R}^{\prime}}\psi_{\boldsymbol{R}^{\prime}\beta}.$
(22)

The subscript “c” in $H_{c}$ signifies that it represents the Hamiltonian
corresponding to a local cluster. For the non-interacting case where
$\chi_{\alpha\beta\gamma\delta}=0$, the Green’s function of the lattice
Hamiltonian is computed by performing the Fourier transform:

$\displaystyle\psi_{\boldsymbol{q}\alpha}$
$\displaystyle=\frac{1}{\sqrt{M}}\sum_{\boldsymbol{R}}\psi_{\boldsymbol{R}\alpha}e^{-i\boldsymbol{q}\cdot\boldsymbol{R}},$
(23)

where $M=\frac{N}{L}$ is the number of clusters in the super-lattice ($N$
being the total number of sites in the original lattice), and the summation
over $\boldsymbol{R}$ covers all super-lattice cells, with $\boldsymbol{q}$ as
the reciprocal vector in the reduced Brillouin zone. Using the inverse
transform of Eq. 23, the non-interacting Hamiltonian is given by
$H_{0}=\sum_{\boldsymbol{q}}\psi_{\boldsymbol{q}}^{\dagger}\left(\mathcal{T}^{0}+\tau_{\boldsymbol{q}}\right)\psi_{\boldsymbol{q}}$,
where $\psi_{\boldsymbol{q}}=\left(\begin{array}[]{c}\psi_{\boldsymbol{q}1}\\\
...\\\ \psi_{\boldsymbol{q}L}\end{array}\right),$ and
$\tau_{\boldsymbol{q}}=\sum_{\boldsymbol{r}}^{\prime}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}\mathcal{T}^{\boldsymbol{r}}$,
with the summation $\sum_{\boldsymbol{r}}^{\prime}$ representing the summation
over all connected clusters. In the non-interacting case, the Green’s function
$\boldsymbol{\mathcal{G}}_{0}(\boldsymbol{q},\omega)$ is given by:

$\displaystyle\boldsymbol{\mathcal{G}}_{0}^{-1}(\boldsymbol{q},\omega)$
$\displaystyle=\omega-\mathcal{T}^{0}-\tau_{\boldsymbol{q}}.$ (24)

For the interacting case, according to the CPT the Green’s function
$\boldsymbol{\mathcal{G}}(\boldsymbol{q},\omega)$ is:

$\displaystyle\boldsymbol{\mathcal{G}}^{-1}(\boldsymbol{q},\omega)$
$\displaystyle=\boldsymbol{\mathcal{G}}_{0}^{-1}(\boldsymbol{q},\omega)-\boldsymbol{\Sigma}(\omega),$
(25)

where $\boldsymbol{\Sigma}(\omega)$ represents the self-energy abstracted from
a reference system $H_{c}\left(\\{f_{i}\\}\right)$. In cases of symmetry
breaking, mean fields $\Delta$ are added to the reference system [28]. The
reference system’s Green’s function for the interacting case is therefore
given as:

$\displaystyle\boldsymbol{G}^{-1}(\omega)$
$\displaystyle=\omega-\mathcal{T}^{0}-\boldsymbol{\Sigma}(\omega).$ (26)

Following Eqs. 2425, and 26, we arrive at the relation:

$\displaystyle\boldsymbol{\mathcal{G}}^{-1}(\boldsymbol{q},\omega)$
$\displaystyle=\boldsymbol{G}^{-1}(\omega)-\tau_{\boldsymbol{q}}.$ (27)

## Appendix: Circuit for $\mathcal{F}\left(P_{i},P_{j}\right)$

In Figure 2, the quantum circuit is initialized with the input state
$|\psi_{in}\rangle=|0\rangle|\boldsymbol{0}\rangle$. The output state, before
the final measurement, is represented as:

$\displaystyle|\psi_{out}\rangle$
$\displaystyle=\frac{1}{2}|0\rangle\left(P_{i}U_{t}+U_{t}P_{j}\right)|g\rangle$
$\displaystyle+\frac{1}{2}|1\rangle\left(U_{t}P_{j}-P_{i}U_{t}\right)|g\rangle.$
(28)

When measuring the ancilla qubit, the probabilities of observing outcomes 0
and 1 are denoted as $p_{+}$ and $p_{-}$, respectively. These probabilities
are calculated as follows:

$\displaystyle p_{\pm}$ $\displaystyle=\frac{1}{4}\langle
g|\left[2\pm\left(P_{j}U_{t}^{\dagger}P_{i}U_{t}+U_{t}^{\dagger}P_{i}U_{t}P_{j}\right)\right]|g\rangle.$
(29)

Consequently, the difference between $p_{+}$ and $p_{-}$ is given by:

$\displaystyle p_{+}-p_{-}$ $\displaystyle=\frac{1}{2}\langle
g|\left(P_{j}U_{t}^{\dagger}P_{i}U_{t}+U_{t}^{\dagger}P_{i}U_{t}P_{j}\right)|g\rangle,$
(30)

This leads to the derivation of
$\mathcal{F}\left(P_{i},P_{j}\right)=2(p_{+}-p_{-})\xrightarrow{p_{+}+p_{-}=1}2(2p_{+}-1).$

## Appendix: Ansatz

Figure 7:  The quantum circuit layout employed in our experiment. Selected
qubits are interconnected and notable for their high readout fidelity, as
highlighted in the middle section of the diagram. The circuit features a
custom-designed ansatz, adapted specifically for our search problem. PM= Pauli
Measurement and $Y2M=R_{y}(-\frac{\pi}{2})$ is a dedicated gate in Xiaohong.

In our experiments, we configured the circuit sampling to consist of 12,000
shots. Figure 7 presents the qubits selected for our study. These qubits,
interconnected with adjacent ones, exhibited outstanding quality during the
experiment, as indicated by the high readout fidelity depicted in the central
part of the figure. We developed a specialized ansatz tailored to our search
problem, as shown in Fig. 7. Given the symmetry inherent in our problem, we
incorporated a single variable parameter, implemented using an $R_{y}$ gate.
To generate the necessary entanglement, two arrays of CNOT gates were
employed. Our design permits experimentation with various single-qubit gate
configurations interspersed between these CNOT gate arrays. For instance,
substituting the $X$ gates with $R_{y}(\pi)$ and $R_{y}(-\pi)$, we derived
another effective ansatz. We believe that a more simplified ansatz is
achievable.

## References

  * [1] P.W. Shor.  Algorithms for quantum computation: discrete logarithms and factoring.  In Proceedings 35th Annual Symposium on Foundations of Computer Science, SFCS-94. IEEE Comput. Soc. Press.
  * [2] Peter W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM Journal on Computing, 26(5):1484–1509, 1997.
  * [3] Lov K. Grover.  A fast quantum mechanical algorithm for database search.  In Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing, STOC ’96, pages 212–219, New York, NY, USA, 1996. Association for Computing Machinery.
  * [4] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Phys. Rev. Lett., 103:150502, Oct 2009.
  * [5] Morten Kjaergaard, Mollie E. Schwartz, Jochen Braumüller, Philip Krantz, Joel I.-J. Wang, Simon Gustavsson, and William D. Oliver.  Superconducting qubits: Current state of play.  Annual Review of Condensed Matter Physics, 11(1):369–395, March 2020.
  * [6] Austin G. Fowler, Matteo Mariantoni, John M. Martinis, and Andrew N. Cleland.  Surface codes: Towards practical large-scale quantum computation.  Phys. Rev. A, 86:032324, Sep 2012.
  * [7] Sebastian Krinner, Nathan Lacroix, Ants Remm, Agustin Di Paolo, Elie Genois, Catherine Leroux, Christoph Hellings, Stefania Lazar, Francois Swiadek, Johannes Herrmann, Graham J. Norris, Christian Kraglund Andersen, Markus Müller, Alexandre Blais, Christopher Eichler, and Andreas Wallraff.  Realizing repeated quantum error correction in a distance-three surface code.  Nature, 605(7911):669–674, May 2022.
  * [8] Youwei Zhao, Yangsen Ye, He-Liang Huang, Yiming Zhang, Dachao Wu, Huijie Guan, Qingling Zhu, Zuolin Wei, Tan He, Sirui Cao, Fusheng Chen, Tung-Hsun Chung, Hui Deng, Daojin Fan, Ming Gong, Cheng Guo, Shaojun Guo, Lianchen Han, Na Li, Shaowei Li, Yuan Li, Futian Liang, Jin Lin, Haoran Qian, Hao Rong, Hong Su, Lihua Sun, Shiyu Wang, Yulin Wu, Yu Xu, Chong Ying, Jiale Yu, Chen Zha, Kaili Zhang, Yong-Heng Huo, Chao-Yang Lu, Cheng-Zhi Peng, Xiaobo Zhu, and Jian-Wei Pan.  Realization of an error-correcting surface code with superconducting qubits.  Phys. Rev. Lett., 129:030501, Jul 2022.
  * [9] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, AlÃ¡n Aspuru-Guzik, and Jeremy L. O’Brien.  A variational eigenvalue solver on a photonic quantum processor.  Nature Communications, 5(1), July 2014.
  * [10] Jules Tilly, Hongxiang Chen, Shuxiang Cao, Dario Picozzi, Kanav Setia, Ying Li, Edward Grant, Leonard Wossnig, Ivan Rungger, George H. Booth, and Jonathan Tennyson.  The variational quantum eigensolver: A review of methods and best practices.  Physics Reports, 986:1–128, 2022.  The Variational Quantum Eigensolver: a review of methods and best practices.
  * [11] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  Nature, 549(7671):195–202, September 2017.
  * [12] Maria Schuld and Francesco Petruccione.  Supervised Learning with Quantum Computers.  Springer International Publishing, 2018.
  * [13] Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink, Jerry M. Chow, and Jay M. Gambetta.  Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.  Nature, 549(7671):242–246, September 2017.
  * [14] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Benjamin Chiaro, Roberto Collins, William Courtney, Sean Demura, Andrew Dunsworth, Edward Farhi, Austin Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Steve Habegger, Matthew P. Harrigan, Alan Ho, Sabrina Hong, Trent Huang, William J. Huggins, Lev Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Erik Lucero, Orion Martin, John M. Martinis, Jarrod R. McClean, Matt McEwen, Anthony Megrant, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Hartmut Neven, Murphy Yuezhen Niu, Thomas E. O’Brien, Eric Ostby, Andre Petukhov, Harald Putterman, Chris Quintana, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Doug Strain, Kevin J. Sung, Marco Szalay, Tyler Y. Takeshita, Amit Vainsencher, Theodore White, Nathan Wiebe, Z. Jamie Yao, Ping Yeh, and Adam Zalcman.  Hartree-fock on a superconducting qubit quantum computer.  Science, 369(6507):1084–1089, August 2020.
  * [15] E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski.  Cloud quantum computing of an atomic nucleus.  Phys. Rev. Lett., 120:210501, May 2018.
  * [16] P. Coleman.  Heavy fermions: electrons at the edge of magnetism.  2006\.
  * [17] Sung-Sik Lee.  Recent developments in non-fermi liquid theory.  Annual Review of Condensed Matter Physics, 9(1):227–244, March 2018\.
  * [18] N. F. MOTT.  Metal-insulator transition.  Rev. Mod. Phys., 40:677–683, Oct 1968.
  * [19] T. H. Hansson, M. Hermanns, S. H. Simon, and S. F. Viefers.  Quantum hall physics: Hierarchies and conformal field theory techniques.  Rev. Mod. Phys., 89:025005, May 2017.
  * [20] Y. Tokura, Y. Okimoto, S. Yamaguchi, H. Taniguchi, T. Kimura, and H. Takagi.  Thermally induced insulator-metal transition in ${\mathrm{lacoo}}_{3}:$ a view based on the mott transition.  Phys. Rev. B, 58:R1699–R1702, Jul 1998.
  * [21] Kichizo Asai, Atsuro Yoneda, Osamu Yokokura, J. M. Tranquada, G. Shirane, and Key Kohn.  Two spin-state transitions in lacoo3.  Journal of the Physical Society of Japan, 67(1):290–296, January 1998.
  * [22] J.C. Slater.  Band theory.  Journal of Physics and Chemistry of Solids, 8:21–25, January 1959\.
  * [23] Walter Metzner and Dieter Vollhardt.  Correlated lattice fermions in $d=\infty$ dimensions.  Phys. Rev. Lett., 62:324–327, Jan 1989.
  * [24] Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg.  Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions.  Rev. Mod. Phys., 68:13–125, Jan 1996.
  * [25] D. Sénéchal, D. Perez, and M. Pioro-Ladrière.  Spectral weight of the hubbard model through cluster perturbation theory.  Phys. Rev. Lett., 84:522–525, Jan 2000.
  * [26] M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Pruschke, and H. R. Krishnamurthy.  Nonlocal dynamical correlations of strongly interacting electron systems.  Phys. Rev. B, 58:R7475–R7479, Sep 1998.
  * [27] A. I. Lichtenstein and M. I. Katsnelson.  Antiferromagnetism and d-wave superconductivity in cuprates: A cluster dynamical mean-field theory.  Phys. Rev. B, 62:R9283–R9286, Oct 2000.
  * [28] M. Potthoff, M. Aichhorn, and C. Dahnken.  Variational cluster approach to correlated electron systems in low dimensions.  Phys. Rev. Lett., 91:206402, Nov 2003.
  * [29] David SÃ©nÃ©chal.  An introduction to quantum cluster methods, 2010.
  * [30] Richard P. Feynman.  Simulating physics with computers.  International Journal of Theoretical Physics, 21(6-7):467–488, June 1982.
  * [31] P. Jordan and E. Wigner.  Z. Physik, 47(9-10):631–651, September 1928.
  * [32] Sergey B. Bravyi and Alexei Yu. Kitaev.  Fermionic quantum computation.  Annals of Phys., 298(1):210–226, May 2002.
  * [33] Suguru Endo, Iori Kurata, and Yuya O. Nakagawa.  Calculation of the green’s function on near-term quantum computers.  Phys. Rev. Res., 2:033281, Aug 2020.
  * [34] Taichi Kosugi and Yu-ichiro Matsushita.  Construction of green’s functions on a quantum computer: Quasiparticle spectra of molecules.  Phys. Rev. A, 101:012330, Jan 2020.
  * [35] Gerald D. Mahan.  Many-Particle Physics.  Springer US, 2000.
  * [36] Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 276(1365):238–257, November 1963.
  * [37] Nature Physics, 9(9):523–523, September 2013.
  * [38] Eleftherios N Economou.  Green’s Functions in Quantum Physics.  Springer Series in Solid-State Sciences. Springer, Berlin, Germany, 3 edition, June 2006.
  * [39] John S. Toll.  Causality and the dispersion relation: Logical foundations.  Phys. Rev., 104:1760–1770, Dec 1956.
  * [40] Suguru Endo, Iori Kurata, and Yuya O. Nakagawa.  Calculation of the green’s function on near-term quantum computers.  Phys. Rev. Res., 2:033281, Aug 2020.
  * [41] Francesco Libbi, Jacopo Rizzo, Francesco Tacchino, Nicola Marzari, and Ivano Tavernelli.  Effective calculation of the green’s function in the time domain on near-term quantum processors.  Phys. Rev. Res., 4:043038, Oct 2022.
  * [42] Rihito Sakurai, Wataru Mizukami, and Hiroshi Shinaoka.  Hybrid quantum-classical algorithm for computing imaginary-time correlation functions.  Phys. Rev. Res., 4:023219, Jun 2022.
  * [43] Jacopo Rizzo, Francesco Libbi, Francesco Tacchino, Pauline J. Ollitrault, Nicola Marzari, and Ivano Tavernelli.  One-particle green’s functions from the quantum equation of motion algorithm.  Phys. Rev. Res., 4:043011, Oct 2022.
  * [44] Ken M. Nakanishi, Kosuke Mitarai, and Keisuke Fujii.  Subspace-search variational quantum eigensolver for excited states.  Phys. Rev. Res., 1:033062, Oct 2019.
  * [45] Masuo Suzuki.  Generalized trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems.  Communications in Mathematical Physics, 51(2):183–190, June 1976\.
  * [46] A.D. McLachlan.  A variational solution of the time-dependent schrodinger equation.  Molecular Physics, 8(1):39–44, January 1964.
  * [47] Ying Li and Simon C. Benjamin.  Efficient variational quantum simulator incorporating active error minimization.  Phys. Rev. X, 7:021050, Jun 2017.
  * [48] Michael A Nielsen and Isaac L Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, Cambridge, England, December 2010.
  * [49] Lin Lin.  Lecture notes on quantum algorithms for scientific computation, 2022.
  * [50] Jingzhe Guo, Huazhe Lou, Jintao Yu, Riling Li, Wang Fang, Junyi Liu, Peixun Long, Shenggang Ying, and Mingsheng Ying.  isq: An integrated software stack for quantum programming.  IEEE Transactions on Quantum Engineering, 4:1–16, 2023.
  * [51] William H. Press, Brian P Flannery, Saul A Teukolsky, and William T Vetterling.  Numerical recipes in FORTRAN 77: Volume 1 of FORTRAN numerical recipes volume 1.  Cambridge University Press, Cambridge, England, 2 edition, September 1992\.
  * [52] M. Louis NÃ©el.  PropriÃ©tÃ©s magnÃ©tiques des ferrites ; ferrimagnÃ©tisme et antiferromagnÃ©tisme.  Annales de Physique, 12(3):137–198, 1948.
  * [53] Tudor Giurgica-Tiron, Yousef Hindy, Ryan LaRose, Andrea Mari, and William J. Zeng.  Digital zero noise extrapolation for quantum error mitigation.  In 2020 IEEE International Conference on Quantum Computing and Engineering (QCE). IEEE, October 2020.
  * [54] Hengyue Li, Junya Otsuki, Makoto Naka, and Sumio Ishihara.  Optical response in the excitonic insulating state: Variational cluster approach.  Phys. Rev. B, 101:125117, Mar 2020.