# Exploiting Maximally Mixed States for Spectral Estimation by Time Evolution

Kaelyn J. Ferris  IBM Quantum, IBM T.J. Watson Research Center, Yorktown
Heights, New York 10598, USA  Zihang Wang  Department of Physics, University
of California Santa Barbara, Santa Barbara, California 93106, USA  Itay Hen
Information Sciences Institute, University of Southern California, Marina del
Rey, CA 90292, USA  Department of Physics and Astronomy, and Center for
Quantum Information Science & Technology, University of Southern California,
Los Angeles, California 90089, USA  Amir Kalev  Information Sciences
Institute, University of Southern California, Arlington, VA 22203, USA
Department of Physics and Astronomy, University of Southern California, Los
Angeles, California 90089, USA  Nicholas T. Bronn  IBM Quantum, IBM T.J.
Watson Research Center, Yorktown Heights, New York 10598, USA  Vojtěch Vlček
Department of Chemistry and Biochemistry, University of California, Santa
Barbara, CA 93106, USA  Department of Materials, University of California,
Santa Barbara, CA 93106, USA

###### Abstract

We introduce a novel approach for estimating the spectrum of quantum many-body
Hamiltonians, and more generally, of Hermitian operators, using quantum time
evolution. In our approach we are evolving a maximally mixed state under the
Hamiltonian of interest and collecting specific time-series measurements to
estimate its spectrum. We demonstrate the advantage of our technique over
currently used classical statistical sampling methods. We showcase our
approach by experimentally estimating the spectral decomposition of a 2-qubit
Heisenberg Hamiltonian on an IBM Quantum backend. For this purpose, we develop
a hardware-efficient decomposition that controls $n$-qubit Pauli rotations
against the physically closest qubit alongside expressing two-qubit rotations
in terms of the native entangling interaction. This substantially reduced the
accumulation of errors from noisy two-qubit operations in time evolution
simulation protocols. We conclude by discussing the potential impact of our
work and the future directions of research it opens.

##  I Introduction

Tackling the quantum many-body problem for any realistic system is, besides
fundamental theoretical advances, inevitably intertwined with the developments
of new computational techniques.[1] The problem complexity scales
exponentially with the system size. Yet, the physically interesting phenomena
are typically related to finding the ground or selected excited states (i.e.,
those accessible on characteristic energy scales), the relevant observables,
and measurable correlation functions. In this context, a class of random
sampling methods has a key advantage over deterministic calculations as they
allow to offset the steep computational cost by recasting the desired
quantities (e.g., correlation functions) as statistical estimators over random
realizations that explore the Hilbert space of the problem [2]. The most
conceptually straightforward problem is associated with determining the
spectrum of a Hamiltonian, $H$, by the time evolution of random vectors and
the spectral analysis of their autocorrelation. Since only information on
selected states is typically desired, this may be achieved by employing a set
of random states in a particular subspace.

In principle, a single vector constructed as a uniformly weighted linear
combination of the selected stationary states provides access to the spectrum
of interest. However, as the spectral decomposition of $H$ these stationary
states are not known a priori. As a result, each realization represents a
particular sampling of the overlap with the eigenvectors of $H$. Since every
sample has varying (and non-uniform) overlap, the calculation thus needs to be
repeated and averaged, while the standard deviation of the spectrum estimation
decreases as $1/\sqrt{N}$ with the number of samplings $N$. In principle, the
exact results are guaranteed (for unbiased sampling) when $N\to\infty$. The
prefactor of the statistical error is largely determined by the “information
redundancy” in the sampled subspace. Random sampling techniques for weakly
correlated systems generally exhibit self-averaging and the desired quantity
is well captured with only a few random vectors leading to a significant
speedup of perturbative many-body techniques.[3, 4, 5, 6] For strongly
correlated systems, however, the problem representation cannot be simply
reduced and converging the statistical estimator requires $N$ similar to the
dimensionality of the correlated Hilbert space.

In contrast to classical random sampling methods, quantum mechanics allows the
realization of uniformly distributed random vectors sampling the Hilbert
(sub)space in the form of a single maximally mixed state. As such, each such
random vector has uniform overlap with the eigenvectors of $H$. Manipulating
in parallel all the vectors that span the relevant subspace on a quantum
computer may provide a distinct advantage over a classical counterpart.
Namely, if the system is highly degenerate, a random sampling of pure states
may introduce bias in the evaluated multiplicity of the power spectrum which
can be avoided by sampling directly from a maximally mixed sate. For each
preparation, the quantum phase estimation (QPE) algorithm [7, 8] allows one in
principle to calculate a single eigenvalue of $H$ with uniform probability.
This process can be repeated until all $N$ eigenvalues have been determined.
The quantum method therefore does not suffer from the overhead of classical
random sampling, which becomes particularly demanding as applied to strongly-
correlated systems. However, QPE is not viable on near term noisy quantum
computers. As such, several modifications to this protocol have been proposed
such as robust phase estimation [9, 10], variational quantum phase estimation
[11, 12, 13, 14], exploiting spectrographic techniques [15], or by utilizing
mid-circuit measurements [16].

In this work, we outline and practically demonstrate an approach to eigenvalue
spectrum estimation – based on a time-series analysis of the Hadamard test
[17, 18] – in which a truly random (uniformly distributed) vector is prepared
from a maximally mixed state on a quantum computer. Here, the states span the
entire Hilbert space, but its projection on a selected subspace can be easily
performed by filtering [2]. We further demonstrate the advantage of the
quantum representation over the classical (statistical) sampling methods as
well as how the two approaches scale with the system size. We further provide
a detailed discussion of the practical limitations of these methods.

In Section II we generally describe the stochastic and quantum approaches to
approximating the spectrum of a given Hamiltonian by time evolution. In
Section III we discuss a simple example system on which to benchmark these two
approaches to the eigenvalue spectrum estimation. The results of these
approaches are then discussed in Section IV and conclusions and future
perspectives are given in Section V.

##  II Spectral Estimation via Time Evolution

Let $H$ be a time-independent Hamiltonian on Hilbert space of dimension $d$.
The time evolution operator under this Hamiltonian is given by

$\displaystyle U(t)=\mathrm{e}^{-iHt}\quad{\rm
taking}\quad\ket{\psi(0)}\xlongrightarrow{U}$

where $t$ denotes the evolution time and $\psi$ is an arbitrary (pure) quantum
state ($\hbar\equiv 1$ throughout). A simple yet potentially powerful
realization is that the trace of $U(t)$, which we denote by $u(t)$, has
information about the entire spectrum of $H$

$\displaystyle
u(t)=\text{Tr}\big{(}U(t)\big{)}=\sum_{j}\text{e}^{-i\lambda_{j}t},$ (1)

where the $\lambda_{j}$ are the eigenvalues of the Hamiltonian (which may be
degenerate). Consequently, the Fourier transform of $u(t)$,
$\tilde{u}(\omega)={\cal F}(u(t))$, provides the eigenvalue spectrum of $H$
(including the multiplicity of degenerate eigenvalues) in the frequency
domain:

$\displaystyle u(t)\xlongrightarrow{{\cal
F}}\tilde{u}(\omega)=2\pi\sum_{j}\delta(\omega-\lambda_{j}).$ (2)

in the limit of infinite time. This mathematical setup inspires numerous
techniques [19, 20, 21, 22] and it is also employed in the proposed quantum
protocol for estimating the spectrum of $H$. In the first case, it is employed
in a suite of random sampling algorithms that emerged as numerically
advantageous in evaluating expectation values and correlators [2]. A classical
implementation of stochastic sampling of the spectrum of the system
Hamiltonian employs a set of random vectors $\\{\ket{\psi_{\rm
stoc}^{(k)}}\\}_{k}$ (each representing a pure quantum state) for which the
stochastic resolution of identity holds:
$\lim_{K\to\infty}\frac{1}{K}\sum_{k=1}^{K}\ket{\psi^{(k)}_{\rm
stoc}}\bra{\psi^{(k)}_{\rm stoc}}={\mathbf{I}}$. In practice, the number of
stochastic vectors, $K$, is finite, and each sampling state is generated as a
linear combination:

$\ket{\psi_{\rm
stoc}^{(k)}}=\frac{1}{\sqrt{S}}\sum_{j=1}^{S}\alpha_{j}^{(k)}\ket{E_{j}^{(k)}},$
(3)

where $k\in[1,K]$ labels the state, $\ket{E_{j}^{(k)}}$ are efficiently-
computable basis states, typically chosen as eigenstates of an auxiliary
Hamiltonian $\tilde{H}$. Here, $\alpha_{j}^{(k)}\equiv e^{i\theta_{j}^{(k)}}$
is a random phase; the choice of $\tilde{H}$ is not unique, and, typically,
one employs $\tilde{H}$ that is related to the problem in question, e.g., the
corresponding mean-field approximation to $H$ is assumed in the context of
perturbation expansion [5, 6]. The corresponding spectrum for each random
state has the form, $\tilde{\psi}_{\rm
stoc}^{(k)}(\omega)=\mathcal{F}\left(\psi_{\rm stoc}^{(k)}(t)\right)$ where

$\displaystyle\psi_{\rm stoc}^{(k)}(t)=\bra{\psi_{\rm
stoc}^{(k)}}e^{-iHt}\ket{\psi_{\rm stoc}^{(k)}}.$ (4)

The statistical average of the sample spectrum should converge, with a
variance scaling $\sigma^{2}\sim 1/K^{2}$, to the spectrum of $H$,

$\tilde{\psi}_{K}(\omega)=\frac{1}{K}\sum_{k=1}^{K}\tilde{\psi}_{\rm
stoc}^{(k)}(\omega),\hskip
5.69046pt\lim_{K\to\infty}\tilde{\psi}_{K}(\omega)=\tilde{u}(\omega),$ (5)

where the last equality uses the central limit theorem.

In practice, the evolution operator, truncated to first order, has the form

$U(T)=\mathrm{exp}\left(-iHT\right)\approx\left(\mathbb{I}-iH\Delta
t\right)^{n},$ (6)

where $\Delta t=T/n$ defines the unit time step, and the total evolution time
$T$ is inversely proportional to the frequency resolution,
$\Delta\omega=2\pi/T$. The time resolution must be small compared to the
inverse of the maximum eigenvalue ($\Delta t\ll 1/\omega_{\rm max}$) in order
to perform the truncation in Eq. 6. Additionally, the frequency resolution
$\Delta\omega\ll\mathrm{Min}|\omega_{i}-\omega_{j}|$ must be small enough to
resolve the differences between eigenvalues. This sets the minimum total
evolution time: $T\gg 2\pi/\mathrm{Min}|\omega_{i}-\omega_{j}|$. The ratio
between these two time scales $T/\Delta t$ sets the operation costs
$\epsilon\sim\mathcal{O}(2\pi\omega_{\rm
max}/\mathrm{Min}|\omega_{i}-\omega_{j}|)$. As the system size grows, the
minimal difference between eigenvalues generally reduces, and the maximum
eigenvalue generally increases, resulting in an increased computational cost.

Figure 1: An illustration of a protocol to measure $\text{Tr}(U)$. (a) Circuit
diagram generating the maximally mixed state for a single qubit. Two qubits
are first prepared in a Bell state, then one is discarded while the other is
operated on further. (b) A circuit for measuring $\text{Tr}(U)$, where
$\rho_{\rm MMS}=\frac{1}{d}\mathbf{I}$. The circuit is repeated many times. In
each realization we measure the pointer qubit (top register in the figure) in
either the Pauli-$X$ basis or the Pauli-$Y$ basis (the basis of measurement is
implicit in the meter). Over many repetitions we estimate
$2\left<\sigma_{+}\right>=\left<X+i\;Y\right>$ from which we can approximate
$\frac{1}{d}\text{Tr}(U)$.

& —+⟩ 1 garbage —0⟩ ρ_MMS (b) @C=0.9em @R=0.8em &
$2\left<\sigma_{+}\right>=\frac{1}{d}\;\text{Tr}(U)$ —+⟩ 1 ρ_MMS / U /

In principle, the need for the use of random samplings can be circumvented if
each stochastic vector has a uniform overlap with the eigenstates of $H$. We
here propose and implement a protocol that employs the time evolution on
quantum computers with random vectors drawn from such a distribution,
surpassing the classical limitations of the stochastic implementation. We
argue that this protocol is amenable to implementation on near-term quantum
computers.

This protocol entails constructing a maximally mixed state $\rho_{\rm MMS}$
(shown in Fig. IIa) as the input to a computational register, in contrast to
an initial pure state (which may be prepared to have some overlap on the
relevant subspace), for usual time evolution on a quantum computer controlled
by the pointer qubit, depicted in Fig. IIb. While the protocol can provide
input for many quantum algorithms, including QPE [7, 8], here we apply time
evolution via the Hadamard test to calculate eigenvalues. We also note that
this technique can be easily extended to correlation functions in general [18,
17], because hardware requirements (e.g., error rate, circuit depth, coherence
time) are within the capabilities of current noisy quantum computers for time-
series experiments on small enough computational registers [23]. From here,
the usual $X$\- and $Y$-basis measurements as a function of time are used to
construct the time series (see Appendix)

$\displaystyle\langle
X+iY\rangle(t)=2\langle\sigma^{+}\rangle(t)=\frac{1}{d}{\rm Tr}(U(t)),$ (7)

which determines the spectrum through a fourier transformation. In our
implementation we create the maximally mixed state by entangling each qubit in
the computational register with a “garbage” qubit which remains unmeasured
(and presumably coherent) throughout the algorithm execution.

##  III Eigenvalue Spectrum Estimation

As a concrete example to contrast the stochastic sampling and quantum
evolution techniques, we consider the general one-dimensional spin-1/2
Heisenberg Hamiltonian with a $\hat{\mathbf{z}}$ direction external field
coupling, which has the form

$\hat{H}=-J\sum_{\langle
ij\rangle}\left(X_{i}X_{j}+Y_{i}Y_{j}+Z_{i}Z_{j}\right)-B\sum_{i}Z_{i},$ (8)

where $X_{i},Y_{i}$, and $Z_{i}$ are Pauli operators at site $i$, $\langle
ij\rangle$ sums nearest-neighbor (NN) sites, with a coupling strength $J$. At
each site, electrons can take spin up $\ket{\uparrow}(\equiv\ket{0})$ or spin
down $\ket{\downarrow}(\equiv\ket{1})$ configuration in the Pauli-$Z$ basis.
Physically, the Pauli operators $X_{i}$ and $Y_{i}$ generate spin flips
between NN sites, leading to the exchange interaction between spins. The
$Z_{i}$ operator, on the other hand, favors parallel alignment between NN
sites, resulting in a ferromagnetic coupling.

In this work, we asses the performance of the quantum algorithm running on a
noisy quantum hardware in comparison to stochastic methods, by implementing a
2-site Heisenberg model. In this case, the Hamiltonian can be represented in a
matrix form

$H=-\begin{bmatrix}J+2B&0&0&0\\\ 0&-J&2J&0\\\ 0&2J&-J&0\\\
0&0&0&J-2B\end{bmatrix},$ (9)

in the Pauli-$Z$ basis
$\\{\ket{\downarrow\downarrow},\ket{\uparrow\downarrow},\ket{\downarrow\uparrow},\ket{\uparrow\uparrow}\\}$.
Its eigenvalues eigenvalues are, $-J\pm 2B$ and $J\pm 2J$, where degeneracy
occurs when the external field $B=\pm 2J$.

###  III.1 Quantum Hardware Approach

####  III.1.1 Preparing a maximally mixed state

Figure 2:  Schematic depiction of preparing a maximally mixed state of 4
qubits on the heavy-hexagonal topology of IBM Quantum’s superconducting
processors. The purple qubits and dashed $X$ lines represent the garbage
qubits while the teal qubits and arrows represent the computation qubits. The
pointer qubit is represented in green with $\ket{+}$. Due to the connectivity,
three SWAP gates (represented as crossed wires in the figure) are required in
this state preparation protocol.

A maximally mixed qubit register is experimentally prepared by preparing it,
together with a garbage qubit register, in one of the Bell states, for
example, using the circuit depicted in Fig. IIa. Ideally, the garbage qubit
should be left unmeasured and protected against relaxation processes during
the computation. To understand why, consider the state
$\ket{\phi^{+}}=(\ket{0}_{c}\ket{0}_{g}+\ket{1}_{c}\ket{1}_{g})/\sqrt{2}$. If
one of the qubits is subjected to a measurement, then the other qubit is
projected onto the same state. In such a case, instead of evolving the
maximally mixed state, only a single computational basis state is evolved.
Note, that it is enough to protect the qubit against any mechanism that
corresponds to a measurement process. On the other hand if the garbage qubit
goes through a trace preserving quantum channel $\cal{E}(\cdot)$, then the
reduced state of the computation qubit is

$\rho_{c}=\frac{1}{2}\sum_{i,j=0}^{1}{\rm Tr}[{\cal
E}(\ket{i}_{g}\\!\bra{j})]\ket{i}_{c}\\!\bra{j}=\frac{1}{2}\mathbf{I}.$ (10)

Therefore, the preparation of the maximally mixed state is robust against
depolarizing channel or other trace-preserving channels acting on the garbage
qubit.

Another consideration in preparing the maximally mixed state is the limited
connectivity between qubits in planar architectures. As controlled-NOT CX
gates are the major source of error in superconducting qubit architectures,
with error rates in the range of 0.5-1%, it is important to minimize the
number of SWAP gates (consisting of up to three CX gates). A diagram depicting
how minimizing the number of CX gates may be achieved for a one-dimensional
line of qubits on IBM’s heavy-hexagonal topology [24] is displayed in Fig. 2.

Figure 3: Various decompositions of a ctrl-$R_{zz}(\theta)$ gate (note that
any other two-qubit Pauli rotation can be decomposed to $R_{zz}(\theta)$ by
placing the appropriate single-qubit gates on either side of the operation).
(a) The default decomposition using Qiskit’s standard equivalence library in
a(i) contains two Toffoli gates which decompose to six CX gates (d) is very
inefficient. Alternatively in (b), the controlled $R_{zz}(\theta)$ gate can be
written using two $R_{zz}(\theta)$ gates as shown in b(i). In (c) we show the
decomposition to a single ctrl-$R_{z}(\theta)$ gate between two CX gates c(i).
The single ctrl-$R_{z}(\theta)$ gate is then expressed by a $Z$-rotation for
half of the total angle $\theta$ followed by an $R_{zz}(-\theta/2)$ rotation
in c(ii) which realizes a circuit with a minimal amount of two-qubit gates. In
both (b) and (c), the $R_{zz}(\theta)$ gates can be expressed as a pair of CX
gates between a single $R_{z}(\theta)$ (which are realized virtually via frame
changes [25] and therefore noiseless) as shown in boxes in (b)(ii) and
(c)(iii). Here $R_{zz}(\theta)$ rotations can be constructed by pulse scaling
to reduce the effective gate time [26]

& 1 1R_zz(θ) R_zz(θ) (i) @C=1em @R=.7em & 2 2 2 2 = 1 1 R_z(θ2) R_z(-θ2) (d)
@C=0.7em @R=2.0em & 2 1 @C=0.5em @R=0.9em & 2 2 1 T 1 = 1 1 T T^† H T^† T T^†
T H

(b) @C=1em @R=.5em

& 1 1R_zz(θ) R_zz(θ) (i) @C=1em @R=.7em & 1 1 = 1R_zz(θ/2) 1R_zz(-θ/2)
R_zz(θ/2) R_zz(-θ/2) (ii) @C=1em @R=.7em & 1 1 = 1 1 1 1 R_z(θ2) R_z(-θ2)
22341.4em– 26381.4em–

(c) @C=1em @R=.7em

& 1 1R_zz(θ) R_zz(θ) (i) @C=.5em @R=.7em & 1 = R_z(θ) -1 -1 (ii) @C=0.4em
@R=.7em & 1R_zz(-θ/2) = R_z(θ/2) R_zz(-θ/2) -1 -1 (iii) @C=0.8em @R=.7em & 1 1
= R_z(θ/2) R_z(-θ/2) -1 -1 14261.2em–

####  III.1.2 Implementation of controlled time-evolution

Our protocol requires an implementation of a controlled time evolution
circuit. However, the decomposition of an arbitrary controlled unitary circuit
using Qiskit’s default equivalence library synthesizes gate sequences [27]
that cannot be realistically executed (see Fig. III.1.1a(i)) on noisy near-
term quantum hardware, even for the modest 2-qubit Hamiltonian we are
considering here. The limited connectivity between qubits potentially requires
SWAP gates, which additionally provides a challenge to minimizing the number
of CX gates.

To overcome this, we first realize that, up to single-qubit rotations, any
controlled-Pauli rotation is locally equivalent to

$cP_{R(\theta)}\sim\ket{0}\bra{0}+e^{-i\frac{\theta}{2}\bigotimes_{j}P^{i}_{j}}\ket{1}\bra{1}$
(11)

where $P^{i}\in\\{I,Z\\}$ consisting of a tensor product of identity and
Pauli-$Z$ operators, often referred to in the Pauli string notation where
position in the string indicates which qubit is operated on. These Pauli
rotations can be realized as a single-qubit $R_{z}(\theta)$ on a single choice
of qubit with a $Z$ term in the Pauli string, and CXs targeting this qubit. It
is controlled on every other qubit with a $Z$ in the Pauli string, both before
and after the $R_{z}(\theta)$. Therefore, when constructing the controlled-
rotation, one is free to place the $R_{z}(\theta)$ on the qubit with the least
distance to the pointer qubit within the connectivity graph of the quantum
computer, as this reduces the number of required SWAPs (and hence CXs).

Next, we implement the controlled $Z$-rotation via an uncontrolled rotation by
half the angle $R_{z}(\theta/2)$ followed by an $R_{zz}(-\theta/2)$ gate. The
single-qubit $Z$-rotation is implemented without error by a frame change [25].
The two-qubit $ZZ$-rotation can be implemented with scaled echoed cross
resonance gates [26]. This technique allows us to reduce the number of CX
gates by 75% (from roughly 100 to about 25), for the specific Hamiltonian we
are considering (Table 1). Performing the uncontrolled evolution first is
important as it allows the preparation of the superposition on the pointer
qubit (by Hadamard gate) as late as possible, preventing unnecessary
decoherence.

Synthesis | CX Count | CX Count(mapped) | ECR Count
---|---|---|---
Fig. III.1.1a | 50 | 67 | 104
Fig. III.1.1b | 22 | 23 | 18
Fig. III.1.1c | 16 | 16 | 15
Table 1: CX count for each circuit synthesis technique shown in Figure III.1.1
which constructs the controlled time evolution of the Heisenberg dimer. The
values shown are counted before and after mapping to the qubit connectivity of
the device. The number of echoed cross-resonance (ECR) pulses (both scaled and
unscaled) after mapping is also listed for these synthesis techniques.

The careful decomposition of the controlled time-evolved Heisenberg
Hamiltonian is highlighted in Fig. III.1.1 by considering just the $ZZ$ term.
Fig. III.1.1a shows a decomposition first using Qiskit’s standard equivalence
library for the $R_{zz}(\theta)$ gate followed by the addition of controls on
each individual gate. Fig. III.1.1b shows a decomposition by controlling the
direction of rotation of a $R_{zzz}(\theta/2)$ gate (consisting of an
$R_{zz}(-\theta/2)$ sandwiched between 2 CXs) followed by an action of the
equivalence library on the $R_{zz}$ gates in the dashed boxes. The most
hardware-efficient decomposition appears in Fig. III.1.1c where the control is
placed on the single-qubit $R_{z}(\theta)$ necessary to achieve the Pauli
rotation, which is further decomposed into an error-free $R_{z}(\theta/2)$
implemented virtually [25] and an $R_{zz}(-\theta/2)$ gate that produced by
scaling the native cross resonance interaction [26].

##  IV Results

To assess our eigenvalue estimation approach, we compare the spectral
decomposition of a Heisenberg dimer with an energy ratio $J/B=1$ between the
classical stochastic method and the quantum hardware approach. We obtain this
first by stochastically sampling the 2-site Hilbert space with $S=100$ samples
and evolve each sampled state to time $T\approx 2\pi J^{-1}$ with time steps
$dt=10^{-4}J^{-1}$ following the common implementation, as detailed, e.g., in
Ref. [2]. The resulting average expectation value of $\langle\psi_{\rm
stoc}^{(k)}(t)\rangle$ and eigenvalue spectrum $\tilde{\psi}_{K}(\omega)$ are
then normalized and plotted in green in Figs. 5a and 5b respectively – along
with the exact evolution.

Figure 4: Lifetime of maximally mixed state (a) The fidelity between the state
of the computational qubits with the maximally mixed state as a function of
time. Two garbage qubits are prepared in a maximally entangled state with two
computational qubits, and left idle for time $t$ before tomographic
measurement. The probability of obtaining the $|00\rangle$ state increases
with time due to qubit relaxation.

For the quantum approach, the time evolution was simulated using the 7-qubit
ibm_lagos device and the controlled-unitary $U_{c}(t)$ was generated using the
approach described in this work. We first prepared two Bell pairs then
executed a controlled evolution of two qubits (one from each pair) under a
2-site Heisenberg Hamiltonian. The same energy ratio of $J/B=1$ was chosen and
the system was evolved for time $T\approx 2\pi~{}J^{-1}$ and $dt\approx
0.04~{}J^{-1}$ under the first order Trotter-Suzuki approximation with a
single Trotter step. For each time step, a number of shots $N=8192$ of the
observable $\langle X+iY\rangle$ was performed. Experimentally, the system
suffers from a non-trace-preserving relaxation channel which will destroy the
maximally-mixed state between the computation and garbage qubits. To ensure
our circuits are much shorter than this decay, in Fig. 4 we show the effective
lifetime of the two-qubit maximally mixed state by preparation and tomographic
measurement of the computational register after idle time $t$. Given that our
experiments execute within 9.5$\mu s$, the fidelity of the maximally mixed
state is high enough such that our analysis will return accurate eigenvalue
spectra.

If perfect single-shot measurement is assumed, the standard deviation of the
expectation of an observable is $\Delta\langle\mathcal{O}\rangle=2/\sqrt{N}$,
where we require 2: $\langle X\rangle$ and $\langle Y\rangle$, giving standard
deviation $1/\sqrt{N}$ overall. Given a uniform sampling over the Hilbert
space $d=2^{N}$ we can select windows of time evolution to lessen the
degradation to the signal from measuring multiple eigenvalues. Additionally,
to mitigate hardware noise, we execute all $R_{zz}(\theta)$ gates via cross-
resonance pulse scaling [26]. By doing so, the error rate of these two-qubit
Pauli rotations is linearly scaled according to the angle $\theta$. The data
was then projected onto an identical time grid to the stochastic approach via
a quadratic interpolation and processed by a discrete Fourier transform. The
corresponding normalized eigenvalue spectrum is shown in blue in Fig. 5b.

The precision of unitary evolution techniques is limited by the total time
evolved, which is ultimately limited by the number of gates that can be
completed within the qubit coherence times. With this limitation, however,
results agree well compared to classical numerics for $T=6$. Further
improvements are expected as qubit coherence times increase.

(a)

(b)

Figure 5: Spectral comparison between a classical stochastic method and
controlled quantum propagation on a 2-site Heisenberg model (at $J/B\equiv
1$). (a) The time evolution amplitudes: the dash-line in black represents the
exact amplitude, e.g. $u(t)$; The lines in blue are expectation values of the
Pauli-Y operator ($\braket{Y}$), acquired from the ancilla qubit; The line in
green labels the amplitude generated by the stochastic technique, with sample
size $S=100$. The unit time step for the classical evolution is
$dt=10^{-4}~{}J^{-1}$, $T=6~{}J^{-1}$, and for easy comparison, we project the
ibm_lagos data onto an identical time and frequency grid used in the classical
algorithm via a quadratic interpolation. (b) The corresponding normalized
power spectrum generated by the discrete Fourier transform (DFT), with the
same color labeling. Both spectra generated are in close agreement with the
exact spectrum (black vertical dash-line). The harmonic peaks at even integer
frequencies are artifacts generated by the frequency interpolation scheme via
the DFT.

As further analysis, we estimate the fidelity of this controlled time
evolution for large chain lengths in Fig. 6. To obtain these estimates, we
prepare the eigenvalue spectrum estimation protocol described in Section III.1
for a given chain length of $n$. The routing of the associated circuit is then
optimized to ibm_sherbrooke, a 127-qubit device configured to a heavy-hex
connectivity graph. Finally, the calibration information of the device is
gathered to implement the weight-2 Pauli rotations (as shown in Fig. III.1.1).
An estimated fidelity was then calculated using the prepared circuit –
evolving the system to time $t\approx 2\pi~{}J^{-1}$ with the same energy
ratio and time evolution approximation – by taking a product of the error
rates obtained from the device calibration for all of the single and two-qubit
gates. We also take into account the linear scaling of the error rates for the
cross-resonance pulses for this estimate [26]. The results here demonstrate an
estimated fidelity less than $0.5$ once the chain has grown to $10$ sites. An
inset shows a sample layout of a small, 6-site chain on the ibm_sherbrooke
device with the system qubits highlighted in green and the “garbage” qubits
highlighted in pink. Additionally, we note that because the transpilation
procedure utilizes a probabilistic heuristic to insert needed SWAP gates [28],
the fidelity estimate possesses some variance in what should be a smooth
function.

Figure 6:  Estimated fidelity of time evolution as a function of chain length.
For each chain we prepare a maximally mixed stated and the controlled time
evolution protocol then estimate the fidelity by taking a produce of the error
rates obtained from the device calibration. Inset shows a connectivity graph
of the ibm_sherbrooke device with a sample layout. The computation qubits are
highlighted in green while the “garbage” qubits used to prepare the MMS are
shown in pink. The phase qubit used to measure $\langle X+iY\rangle$ is shown
in black.

##  V Perspective and Conclusions

This work explores the practical advantages of quantum realization of
randomized sampling for spectral estimation, which (similar to the classical
counterpart) enables exploiting the information redundancy. In essence, the
classical approach requires multiple re-samplings, while the quantum
realization allows constructing “ideal” (uniformly distributed) random states
even on near-term devices and yields an unbiased estimation of correlation
functions.

The quantum resources required for this sampling technique scale favorable
with system size. The time step sets the cutoff frequency and thereby largest
eigenvalue that can be observed, whereas the total time evolution sets both
the minimum eigenvalue that can be observed as well as the minimum resolution
between eigenvalues. The total time evolution is currently limited by qubit
decoherence, as circuit depth increases with larger number of Trotter steps.
Neglecting measurement error, which maybe be mitigated [29], the signal will
be limited due to shot noise to a standard deviation of $1/\sqrt{N}$ where $N$
is the number of identical experiments performed for each the $X$ and $Y$
expectation values. For example, here $N=8192$ experiments were performed
twice for each time step at a typical rate of 4 kHz. This signal is further
suppressed due to spectral weight transfer among all the eigenstates
contributing to the time evolution. Since we are typically interested only in
$\tilde{d}$ fraction of the Hilbert space explored, the signal suppression and
shot noise contribute overall as $\sim\tilde{d}/\sqrt{N}$ to the overall
error. Further improvement is possible via denoising through projecting the
measured data onto the “nearest” function satisfying the physical constraints
of the time correlators, as has been recently demonstrated for the equilibrium
and steady states [30].

Besides the hardware advances necessary to achieve longer time evolution (and
hence better frequency resolution), the near-term strategy may alternatively
employ numerical postprocessing on classical devices. In particular,
techniques based on principle mode analysis and dimension reduction technique
through low-rank approximations (such as dynamical mode decomposition [31, 32,
33, 34, 35, 36]). Alternative strategies may employ, for instance, Recursive
Neural Networks, in which the general functional forms of the time evolution
are trained on solvable model.[37]. Such a combination of multiple techniques
is likely optimal for practical applications in the near future and
circumvents the need for a huge technological leap.

In this work, we have practically demonstrated an approach for estimating the
spectrum of a quantum many-body Hamlitonian. We utilize a novel technique to
efficiently construct a controlled time evolution operator $U_{c}(t)$. Our
results show agreement with classical stochastic methods of obtaining
eigenvalue spectra and demonstrate the ability to accurately reproduce
eigenvalue spectra for a small Heisenberg system. With the introduction of
more sophisticated error mitigation techniques [38, 39, 40, 41, 42] and a
direct measurement protocol for the observable measurement [43], quantum
hardware may be able to outperform classical stochastic methods of obtaining
eigenvalues of larger systems – allowing for the accurate characterization of
condensed matter models currently inaccessible using numerical, stochastic, or
other approximating methods.

##  VI Acknowledgment

This research on random sampling (VV and ZW) was supported NSF CAREER award
DMR 1945098. IH and AK would like to acknoledge support from NSF award
#2210374. The authors thank Jeffrey Cohn for a careful reading of the
manuscript and acknowledge the use of IBM Quantum Services for this work. The
collaboration was enabled through the collaboration with Quantum Foundry at
UCSB (NSF award DMR-1906325). In addition, this material is based upon work
supported by the Defense Advanced Research Projects Agency (DARPA) under
Contract No. HR001122C0063.

## Appendix - Hadamard test

Very generally, the Hadamard test concerns with the application of a
controlled-unitary operation, that is the application of a unitary
transformation $U$ on register qubits, conditioned that an ancilla qubit is in
a $\ket{1}$ state. By applying a controlled unitary $U_{c}$ on a state $\rho$
when the control qubit is prepared in the $\ket{+}=(\ket{0}+\ket{1})/\sqrt{2}$
we obtain:

$U_{c}(\ket{+}\bra{+}\otimes\rho)U_{c}^{\dagger}=\frac{1}{2}\Big{[}\ket{0}\bra{0}\otimes\rho+(\ket{0}\bra{1}\otimes(\rho
U^{\dagger})+\ket{1}\bra{0}\otimes(U\rho)+\ket{1}\bra{1}\otimes(U\rho
U^{\dagger})\Big{]}$

For the special case where $\rho$ is the completely mixed state we have

$U_{c}(\ket{+}\bra{+}\otimes\rho)U_{c}^{\dagger}=\frac{1}{2d}\Big{[}\ket{0}\bra{0}\otimes\mathbf{I}+(\ket{0}\bra{1}\otimes
U^{\dagger}+\ket{1}\bra{0}\otimes U+\ket{1}\bra{1}\otimes\mathbf{I})\Big{]}.$

Therefore, the reduced state of the qubit after the application of the
controlled unitary operation becomes

$\rho_{\rm
qubit}=\frac{1}{2}\mathbf{I}+\frac{1}{2d}\Big{[}\ket{0}\bra{1}\text{Tr}(U^{\dagger})+\ket{1}\bra{0}\text{Tr}(U)\Big{]}$

Hence measuring $X+iY=2\ket{1}\bra{0}$ corresponds to measuring
$\text{Tr}(U)/d$.

## References

  * Martin  _et al._ [2016] R. M. Martin, L. Reining, and D. M. Ceperley,  _Interacting Electrons: Theory and Computational Approaches_ (Cambridge University Press, 2016).
  * Baer  _et al._ [2022] R. Baer, D. Neuhauser, and E. Rabani, Stochastic vector techniques in ground-state electronic structure, Annual Review of Physical Chemistry 73, 255 (2022), pMID: 35081326, https://doi.org/10.1146/annurev-physchem-090519-045916 .
  * Neuhauser  _et al._ [2014] D. Neuhauser, Y. Gao, C. Arntsen, C. Karshenas, E. Rabani, and R. Baer, Breaking the theoretical scaling limit for predicting quasiparticle energies: The stochastic $gw$ approach, Phys. Rev. Lett. 113, 076402 (2014).
  * Vlček  _et al._ [2017] V. Vlček, E. Rabani, D. Neuhauser, and R. Baer, Stochastic gw calculations for molecules, Journal of Chemical Theory and Computation 13, 4997 (2017), pMID: 28876912, https://doi.org/10.1021/acs.jctc.7b00770 .
  * Vlček  _et al._ [2018] V. Vlček, W. Li, R. Baer, E. Rabani, and D. Neuhauser, Swift g w beyond 10,000 electrons using sparse stochastic compression, Physical Review B 98, 075107 (2018).
  * Vlček [2019] V. Vlček, Stochastic vertex corrections: Linear scaling methods for accurate quasiparticle energies, Journal of Chemical Theory and Computation 15, 6254 (2019), pMID: 31557012, https://doi.org/10.1021/acs.jctc.9b00317 .
  * Abrams and Lloyd [1997] D. S. Abrams and S. Lloyd, Simulation of many-body fermi systems on a universal quantum computer, Phys. Rev. Lett. 79, 2586 (1997).
  * Lloyd [1996] S. Lloyd, Universal Quantum Simulators, Science 273, 1073 (1996).
  * Kimmel  _et al._ [2015] S. Kimmel, G. H. Low, and T. J. Yoder, Robust calibration of a universal single-qubit gate set via robust phase estimation, Phys. Rev. A 92, 062315 (2015).
  * Rudinger  _et al._ [2017] K. Rudinger, S. Kimmel, D. Lobser, and P. Maunz, Experimental demonstration of a cheap and accurate phase estimation, Phys. Rev. Lett. 118, 190502 (2017).
  * Shen  _et al._ [2023a] Y. Shen, K. Klymko, J. Sud, D. B. Williams-Young, W. A. d. Jong, and N. M. Tubman, Real-Time Krylov Theory for Quantum Computing Algorithms, Quantum 7, 1066 (2023a).
  * Klymko  _et al._ [2022] K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa, N. Wiebe, W. A. de Jong, and N. M. Tubman, Real-time evolution for ultracompact hamiltonian eigenstates on quantum hardware, PRX Quantum 3, 020323 (2022).
  * Shen  _et al._ [2023b] Y. Shen, D. Camps, A. Szasz, S. Darbha, K. Klymko, D. B. Williams-Young, N. M. Tubman, and R. V. Beeumen, Estimating eigenenergies from quantum dynamics: A unified noise-resilient measurement-driven approach (2023b), arXiv:2306.01858 [quant-ph] .
  * Cohn  _et al._ [2021] J. Cohn, M. Motta, and R. M. Parrish, Quantum filter diagonalization with compressed double-factorized hamiltonians, PRX Quantum 2, 040352 (2021).
  * Stenger  _et al._ [2022] J. P. Stenger, G. Ben-Shach, D. Pekker, and N. T. Bronn, Simulating spectroscopy experiments with a superconducting quantum computer, Physical Review Research 4, 043106 (2022).
  * Córcoles  _et al._ [2021] A. D. Córcoles, M. Takita, K. Inoue, S. Lekuch, Z. K. Minev, J. M. Chow, and J. M. Gambetta, Exploiting dynamic quantum circuits in a quantum algorithm with superconducting qubits, Physical Review Letters 127, 100501 (2021).
  * Somma  _et al._ [2002] R. Somma, G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Simulating physical phenomena by quantum networks, Physical Review A 65, 042323 (2002).
  * Somma [2019] R. D. Somma, Quantum eigenvalue estimation via time series analysis, New Journal of Physics 21, 123025 (2019).
  * Hutchinson [1989] M. Hutchinson, A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines, Communications in Statistics - Simulation and Computation 18, 1059 (1989), https://doi.org/10.1080/03610918908812806 .
  * Neuhauser [1990] D. Neuhauser, Bound state eigenfunctions from wave packets: Time→energy resolution, The Journal of Chemical Physics 93, 2611 (1990), https://pubs.aip.org/aip/jcp/article-pdf/93/4/2611/11340586/2611_1_online.pdf .
  * Wall and Neuhauser [1995] M. R. Wall and D. Neuhauser, Extraction, through filter‐diagonalization, of general quantum eigenvalues or classical normal mode frequencies from a small number of residues or a short‐time segment of a signal. I. Theory and application to a quantum‐dynamics model, The Journal of Chemical Physics 102, 8011 (1995).
  * Mandelshtam and Taylor [1997] V. A. Mandelshtam and H. S. Taylor, A low-storage filter diagonalization method for quantum eigenenergy calculation or for spectral analysis of time signals, The Journal of Chemical Physics 106, 5085 (1997), https://pubs.aip.org/aip/jcp/article-pdf/106/12/5085/10781975/5085_1_online.pdf .
  * Preskill [2018] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Chamberland  _et al._ [2020] C. Chamberland, G. Zhu, T. J. Yoder, J. B. Hertzberg, and A. W. Cross, Topological and subsystem codes on low-degree graphs with flag qubits, Phys. Rev. X 10, 011022 (2020).
  * McKay  _et al._ [2017] D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Efficient $z$ gates for quantum computing, Phys. Rev. A 96, 022330 (2017).
  * Stenger  _et al._ [2021] J. P. T. Stenger, N. T. Bronn, D. J. Egger, and D. Pekker, Simulating the dynamics of braiding of majorana zero modes using an ibm quantum computer, Phys. Rev. Res. 3, 033171 (2021).
  * Shende  _et al._ [2006] V. Shende, S. Bullock, and I. Markov, Synthesis of quantum-logic circuits, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25, 1000 (2006).
  * Li  _et al._ [2019] G. Li, Y. Ding, and Y. Xie, Tackling the qubit mapping problem for nisq-era quantum devices (2019), arXiv:1809.02573 [cs.ET] .
  * Nation  _et al._ [2021] P. D. Nation, H. Kang, N. M. Sundaresan, and J. M. Gambetta, Scalable mitigation of measurement errors on quantum computers, PRX Quantum 2, 040326 (2021).
  * Kemper  _et al._ [2023] A. F. Kemper, C. Yang, and E. Gull, On the positive definiteness of response functions in the time domain, arXiv preprint arXiv:2309.02566 (2023).
  * Schmid [2010] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech. 656, 5 (2010).
  * Schmid  _et al._ [2011] P. J. Schmid, L. Li, M. P. Juniper, and O. Pust, Applications of the dynamic mode decomposition, Theor. Comput. Fluid Dyn. 25, 249 (2011).
  * Tu  _et al._ [2014] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, On dynamic mode decomposition: Theory and applications, J. Comput. Dyn. 1, 391 (2014).
  * Reeves  _et al._ [2023] C. C. Reeves, J. Yin, Y. Zhu, K. Z. Ibrahim, C. Yang, and V. c. v. Vlček, Dynamic mode decomposition for extrapolating nonequilibrium green’s-function dynamics, Phys. Rev. B 107, 075107 (2023).
  * Yin  _et al._ [2021] J. Yin, Y.-h. Chan, F. da Jornada, D. Qiu, C. Yang, and S. G. Louie, Analyzing and predicting non-equilibrium many-body dynamics via dynamic mode decomposition, arXiv preprint arXiv:2107.09635 (2021).
  * Mejía  _et al._ [2023] L. Mejía, J. Yin, D. R. Reichman, R. Baer, C. Yang, and E. Rabani, Stochastic real-time second-order green’s function theory for neutral excitations in molecules and nanostructures (2023), arXiv:2303.06874 [physics.chem-ph] .
  * Bassi  _et al._ [2023] H. Bassi, Y. Zhu, S. Liang, J. Yin, C. C. Reeves, V. Vlcek, and C. Yang, Learning nonlinear integral operators via recurrent neural networks and its application in solving integro-differential equations (2023), arXiv:2310.09434 [cs.LG] .
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Sung  _et al._ [2022] K. J. Sung, M. J. Rancic, O. T. Lanes, and N. T. Bronn, Preparing majorana zero modes on a noisy quantum processor (2022), arXiv:arXiv:2206.00563 .
  * 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, and A. Kandala, Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * van den Berg  _et al._ [2023] E. van den Berg, Z. K. Minev, A. Kandala, and K. Temme, Probabilistic error cancellation with sparse Pauli–Lindblad models on noisy quantum processors, Nature Physics 19, 1116 (2023).
  * Majumdar  _et al._ [2023] R. Majumdar, P. Rivero, F. Metz, A. Hasan, and D. S. Wang, Best practices for quantum error mitigation with digital zero-noise extrapolation (2023), arXiv:2307.05203 .
  * Mitarai and Fujii [2019] K. Mitarai and K. Fujii, Methodology for replacing indirect measurements with direct measurements, Phys. Rev. Res. 1, 013006 (2019).