# Quantum Dynamical Emulation

Jacob M. Leamer  0000-0002-2125-7636 Tulane University, New Orleans, LA 70118,
USA  Quantum Algorithms and Applications Collaboratory, Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA  Denys I. Bondar
0000-0002-3626-4804 Tulane University, New Orleans, LA 70118, USA  Gerard
McCaul  0000-0001-7972-456X Tulane University, New Orleans, LA 70118, USA
<EMAIL_ADDRESS>
(February 2024)

###### Abstract

We introduce the concept of Quantum Dynamical Emulation, a constructive method
for mapping the solutions of non-unitary dynamics to a weighted set of unitary
operations. This allows us to derive a new correspondence between real and
imaginary time, which we term Imaginary Time Quantum Dynamical Emulation
(ITQDE). This enables an imaginary time evolution to be constructed from the
overlaps of states evolved in opposite directions in real time. We show that a
single trajectory evolved using ITQDE can be used not only to infer ground and
thermal states, but also to resolve information about the complete Hamiltonian
spectrum. We further employ ITQDE to derive novel thermodynamic results,
including a generalisation of the Hubbard-Stratonovich transform. We go on to
develop a quantum algorithm for computing the spectra of quantum systems that
is based on this premise. We demonstrate the utility of this method through
numerical simulation, as well as quantum hardware implementations.

##  I Introduction

Few mathematical concepts have furthered our understanding of the real to a
greater degree than the imaginary. From its first use in the solution of cubic
polynomials, the introduction of imaginary numbers dramatically expands the
scope of both the solvable and conceivable. In the domain of physics, the
alchemy of $i$ has an enduring capacity to surprise. For example, the Wick
rotation [1] is celebrated as a tangible link between statistical mechanics
and quantum dynamics, where oscillation is transformed to decay. In this
sense, evolution in imaginary time might rightly be thought of as the dynamics
of equilibration.

This equilibrium information, and methods to obtain it, are perhaps the
overriding preoccupation of physics. Beyond the ground state, the spectral gap
is also of fundamental interest due to its role in phase transitions [2, 3, 4,
5, 6], as well as in quantum information [7]. As such, there has been a great
deal of interest in developing methods for determining excited states [8, 9,
10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Moreover, with the discovery of
dynamical phase transitions an open question has arisen as to how these relate
to the system equilibrium. Furthermore, in the regime of quantum chaos, a
parallel question arises as to whether the signatures of chaotic dynamics can
be found in a system’s spectral characteristics [20, 21].

Imaginary time evolution is a powerful method for obtaining equilibrium
information about a quantum system, and is used extensively in computational
methods such as quantum Monte Carlo. The premise is that, as a quantum state
is evolved in imaginary time $\tau$, higher energy states decay exponentially
faster than lower energy states. As long as the system has some initial
support on the ground state, then in the limit of large $\tau$, the system
converges to the ground state.

Ground state calculations are of course central to numerous applications [22,
23, 24], e.g., in solid state physics and quantum chemistry [25, 26, 27, 28],
but with suitable interpretation of the Hamiltonian, can also represent
solutions to computational problems. Beyond ground states, imaginary time
evolution can also be used to calculate thermal states, also known as Gibbs
states, at finite temperatures [29, 30, 31]. Such states are crucial in a wide
range of inquiry, ranging (naturally) from the thermalization of quantum
systems [32, 33, 34, 35, 36] to efficient optimization procedures [37, 38].

While imaginary time evolution is a celebrated technique across multiple
computational methods, it does not map naturally to quantum computers, which
are widely hoped to advance our ability to perform difficult simulations of
quantum systems [39, 40]. This is because algorithms for typical, circuit-
model quantum computers are constructed using quantum logic gates, which are
unitary operations, while imaginary time evolution is not a unitary process.
Nevertheless, for the reasons outlined above, there has been intense interest
in adapting this powerful method to quantum algorithms [41, 42, 43, 44, 45,
46, 47]. Thus far, these methods rely either on variational methods to
approximate the non-unitary operation of imaginary time [48], or probabilistic
approaches [49]. In both cases, the non-constructive nature of these
techniques manifests itself in undesirable features. For example Quantum
Imaginary Time Evolution (QITE) requires unitary operations on a scaling
domain of qubits [48], while Probabilistic Imaginary Time Evolution (PITE)
[49] exhibits an exponential decrease in success probability as one evolves
further into imaginary time.

The apparent incompatibility of imaginary time simulations with the building
blocks of circuit-model quantum computation also speaks to a broader
challenge: our ability to manipulate physical systems is stranded in real
time, while the ultimate answers to the questions we seek to ask can be most
transparently sought with the use of non-unitary dynamics. Any result able to
constructively represent such non-unitary dynamics in terms of unitary
operations would constitute a novel tool for the exploration of the topics
outlined above.

Here, we propose a new approach to this problem, and others, which we term
Quantum Dynamical Emulation (QDE). As outlined above, for reasons both of
aesthetics and practical importance, we choose to apply it first to the
emulation of imaginary time evolutions, terming this specification of QDE
Imaginary Time Quantum Dynamical Emulation (ITQDE). As a first application we
employ it in the implementation of spectral calculations in quantum circuits,
but emphasise that the methods pioneered by QDE are likely to find useful
application in a wide variety of contexts. Using ITQDE, we specifically
demonstrate that expectation values of a state evolved in imaginary time can
be expressed in terms of the overlaps of states that are propagated forwards
and backwards in real time.

Unlike the aforementioned quantum algorithms, which aim at generating
imaginary time in a quantum circuit, ITQDE is a constructive, deterministic
approach that is consequently applicable to a much broader context, and which
establishes novel links between unitary evolutions and Gaussian \- rather than
Gibbs - states. It is established without requiring any reference to
algorithmic results, and yields a number of surprising links between physical
and thermalising dynamics. In particular ITQDE recovers in its continuum limit
a result fundamental to stochastic and quantum thermodynamics - the Hubbard-
Stratonovich transformation, offering a new interpretation of its origin in
the process. Not only this, ITQDE is also of a form manifestly implementable
in the current circuit-model paradigm for quantum computation, possessing the
significant conceptual advantage that it relies on no variational or
probabilistic arguments, and - at most - only a single ancilla qubit. We use
this fact to demonstrate ITQDE can be leveraged to perform sampling-based
spectral calculations of correlated materials in quantum computers.

The remainder of the paper is organised as follows. We begin by expounding the
underlying philosophy of QDE in Sec. II, before applying it to derive the
ITQDE correspondence for representing imaginary time-like dynamics in terms of
unitary operations. The analytical and numerical properties of the technique
are also discussed, together with its continuous limit. We then extend the
derived correspondence in Sec. III to develop a technique for resolving
spectral information from dynamical overlaps. Following this, in Sec. IV we
sketch several potential applications of ITQDE in the realm of thermodynamics.
Sec. V then provides numerical illustrations of ITQDE resolving the spectra of
correlated materials. A first experimental result utilising this
correspondence is then presented in Sec. VI, where a quantum algorithm for
implementing ITQDE is detailed, together with the results obtained from its
execution in both simulation and hardware. We then close the paper with a
discussion of further future applications and generalisations of this work in
Sec. VII.

##  II Linking real and imaginary evolutions with quantum dynamical emulation

At root, QDE exploits a freedom in the interpretation of a set of unitary
operations applied to a state $\rho$. We first define a superoperator

$\mathcal{L}[\rho]=\sum_{ij}k_{ij}\mathcal{U}_{i}\rho\mathcal{U}_{j}$ (1)

where $k_{ij}$ are weights on the application of some set of unitary
operations $\\{\mathcal{U}_{j}\\}$. We note immediately that in general this
need not conform to the form of a Kraus map, and in such a case is not trace-
preserving and therefore not directly physical. Nevertheless this operation is
composed purely of unitary operations on states, which are individually
implementable in a physical context. As we shall see, this allows one to
ultimately interpret the action of this operation as being equivalent to an
averaging over physically realisable unitary evolutions. This is the heart of
QDE, enabling constructive unitary representations of non-unitary dynamics.

To make this connection, we first stipulate our set of unitaries are
parametrised by some increment $\Delta\tau$ (or fractional powers of that
increment). This allows for the interpretation of $\mathcal{L}[\rho]$, as a
infinitesimal generator for a differential equation. Specifically,
$\mathcal{L}[\rho]$ can be Taylor expanded to the form:

$\mathcal{L}[\rho]=\rho+\Delta\tau\mathcal{G}[\rho]+\mathcal{O}(\Delta\tau^{2}),$
(2)

such that $\mathcal{L}^{m}[\rho]$ corresponds to a state that is the solution
to the differential equation

$\frac{{\rm d}}{{\rm
d}\tau}\tilde{\rho}(\tau)=\mathcal{G}[\tilde{\rho}(\tau)]$ (3)

for $\tau=m\Delta\tau$.

Surprisingly, it is relatively easy to choose $\\{\mathcal{U}_{j}\\}$ such
that this differential equation is explicitly of a non-unitary character.
Indeed, if Eq. (1) is not of the Kraus form then the evolution is not trace
preserving and the resulting state definitionally cannot be obtained
unitarily. The result of this differential equation is approximately
equivalent to $m$ repeated applications of this superoperator
$\mathcal{L}^{m}[\rho]$, i.e., $\rho_{m}=\mathcal{L}^{m}[\rho]$. While
$\mathcal{L}^{m}[\rho]$ is not directly implementable, it can after
normalisation be interpreted as anaverage of physically realisable unitary
operations of states. In this regard, $\mathcal{L}^{m}[\rho]$ is realisable as
a set of unitary operations which, in combination, correspond to a state that
is not otherwise directly accessible by such operations. The approximate
equating of the two states $\rho_{m}$ and $\tilde{\rho}(\tau)$ establishes QDE
as a constructive method for unitary representations of non-unitary dynamics
and is schematically represented in Fig. 1.

QDE may rightly be thought of as a spiritual successor to Ensemble Rank
Truncation [50] and corner space techniques [51, 52, 53], where it was
conceived of as an efficient deterministic method for the simulation of
Lindblad equations. QDE greatly extends this idea, however, to capture more
general dynamics. This is most easily demonstrated by outlining QDE’s
application to the problem of unitary representations of imaginary-time
evolved states. We begin with the unitary

Figure 1: Schematic illustration of Quantum Dynamical Emulation (QDE). The
action of a set of unitary operations on a state has an alternative
interpretation as a generator for differential equation. Depending on the
chosen form of $\\{\mathcal{U}_{j}\\}$, this may be a non-unitary evolution.
QDE equates the result of such a non-unitary evolution to the unitary
operation $\mathcal{L}^{m}[\rho]$.

$\mathcal{U}=e^{-i\sqrt{\frac{\Delta\tau}{2}}H},$ (4)

which represents evolution in real time under a Hamiltonian $H$ for a time
$\sqrt{\frac{\Delta\tau}{2}}$. With this, we define the superoperator

$\displaystyle\mathcal{L}[\rho]$
$\displaystyle=\frac{1}{2}\left({\mathcal{U}}\rho(\tau){\mathcal{U}}+{\mathcal{U}}^{\dagger}\rho(\tau){\mathcal{U}}^{\dagger}\right),$
(5)

whose effect at first order in $\Delta\tau$ is given by

$\displaystyle\mathcal{L}[\rho]$
$\displaystyle=\rho-\frac{\Delta\tau}{4}\left(\big{\\{}\rho,H^{2}\big{\\}}+2H\rho
H\right)+\mathcal{O}(\Delta\tau^{2}),$ (6)

where $\\{\cdot,\cdot\\}$ denotes the anticommutator. From Eq. 6, we observe
that for $\Delta\tau\rightarrow 0$, $\mathcal{L}[\rho]$ corresponds to an
infinitesimal evolution of $\rho$ according to

$\displaystyle\frac{d\rho}{d\tau}=-\frac{1}{4}\left(\\{\rho,H^{2}\\}+2H\rho
H\right),$ (7)

whose (unnormalized) solution is given by

$\displaystyle\rho(\tau)=\sum_{ij}e^{-\frac{1}{4}\tau(E_{i}+E_{j})^{2}}p_{ij}(0)|E_{i}\rangle\langle
E_{j}|.$ (8)

In Eq. 8, $\rho(\tau)$ is expressed in the eigenbasis of the Hamiltonian
$H=\sum_{i}E_{i}|E_{i}\rangle\langle E_{i}|$ with general initial condition
$\rho(0)=\sum_{ij}p_{ij}(0)|E_{i}\rangle\langle E_{j}|$.

In the limit of large $\tau$, the dominant contribution in $\rho(\tau)$ will
be the eigenstate of $H$ with the smallest squared eigenvalue, i.e.,

$\rho(\tau)\rightarrow e^{-\tau
E_{\ell}^{2}}p_{\ell\ell}(0)|E_{\ell}\rangle\langle E_{\ell}|$ (9)

for $\tau\rightarrow\infty$, where $E_{\ell}^{2}=\min_{j}E_{j}^{2}$ is the
ground state energy of $H^{2}$. Other matrix elements of $\rho(\tau)$ will be
exponentially suppressed in $\tau$ according to Eq. 8.

The result in Eq. 9 is directly analogous to the outcome of evolving the state
$\rho(0)$ in _imaginary time_ under the Hamiltonian $H^{2}$, and could
accordingly be used to study ground states of $H^{2}$ via Eq. (9). However, we
emphasize that we have obtained Eq. 9 utilizing purely unitary, real-time
evolutions per Eq. 5.

Having obtained Eq. (8), we now consider the repeated application of Eq. 5,
denoting $m$ applications of $\mathcal{L}$ by $\mathcal{L}^{m}[\rho]$. We
consider an initial condition that is given by the pure state

$\displaystyle\rho(0)=|\psi_{0}\rangle\langle\psi_{0}|,$ (10)

noting that the generalization to mixed states is straightforward. Introducing
the following notation for forwards and backwards evolved states,

$\displaystyle|\psi_{j+1}\rangle$
$\displaystyle={\mathcal{U}}|\psi_{j}\rangle,$
$\displaystyle|\psi_{j-1}\rangle$
$\displaystyle={\mathcal{U}}^{\dagger}|\psi_{j}\rangle,$ (11)

a single application of $\mathcal{L}$ yields

$\rho(\Delta\tau)=\mathcal{L}^{1}[\rho(0)]=\frac{1}{2}(|\psi_{1}\rangle\langle\psi_{-1}|+|\psi_{-1}\rangle\langle\psi_{1}|).$
(12)

The effect of repeated applications of $\mathcal{L}$ can then be ascertained
inductively. Namely, after $m$ steps, $\rho(m\Delta\tau)$ may be written as

$\rho(m\Delta\tau)=\mathcal{L}^{m}[\rho(0)]=\sum_{j=0}^{2m}a^{(m)}_{j}|\psi_{j-m}\rangle\langle\psi_{m-j}|,$
(13)

where the $a^{(m)}_{j}$ denote coefficients that remain to be calculated.
These coefficients have two important properties: first,
$a^{(m)}_{0}=a^{(m)}_{2m}=1$, and second, the lack of a static term in
$\mathcal{L}$ guarantees $a^{(m)}_{j}=0$ for odd $j$.

Applying $\mathcal{L}$ again to Eq. (13), we obtain

$\displaystyle
2\mathcal{L}^{m+1}[\rho(0)]=\sum_{j=0}^{j=2m}a^{(m)}_{j}(|\psi_{j-(m+1)}\rangle\langle\psi_{m+1-j}|$
(14) $\displaystyle+|\psi_{j+1-m}\rangle\langle\psi_{m-1-j}|).$

Rearranging indices, we find

$\displaystyle
2\mathcal{L}^{m+1}[\rho(0)]=|\psi_{-m-1}\rangle\langle\psi_{m+1}|+|\psi_{m+1}\rangle\langle\psi_{-m-1}|$
(15)
$\displaystyle+\sum_{j=2}^{j=2m}\left(a_{j}^{(m)}+a_{j-2}^{(m)}\right)|\psi_{j-(m+1)}\rangle\langle\psi_{m+1-j}|.$

This yields the following recursion relation

$\displaystyle 2a_{j}^{(m+1)}=a_{j}^{(m)}+a_{j-2}^{(m)},$ (16)

whose solution is given by Pascal’s identity

$\displaystyle{{m}\choose{k}}={{m-1}\choose{k}}+{{m-1}\choose{k-1}},$ (17)

yielding

$a^{(m)}_{j}=\begin{cases}\frac{1}{2^{m}}{m\choose\frac{j}{2}}&j\
\text{even},\\\ 0&j\ \text{odd}.\end{cases}$ (18)

Using $\tau=m\Delta\tau$, we substitute Eq. 18 into Eq. 13 to finally obtain

$\rho(\tau)=\frac{1}{2^{m}}\sum_{j=0}^{j=m}{{m}\choose{j}}|\psi_{2j-m}\rangle\langle\psi_{m-2j}|.$
(19)

Eq. (19) can be used to simulate the evolution of an initial state $\rho(0)$
in imaginary time under $H^{2}$ using only unitary operations, and thus
enables what we term Imaginary Time Quantum Dynamical Emulation (ITQDE). The
error between Eq. (19) and the continuous-time limit given previously in Eq.
(8) is $\mathcal{O}(m\Delta\tau^{2})$. In the interests of concision, we do
not explicitly retain this error term in the sections that follow, but stress
that that ITQDE as presented is contingent on this discrete approximation. The
continuous limit of ITQDE of this approximation is addressed in Sec. IV.

In many cases, it is desirable to additionally enforce that the state is
normalized at all (imaginary) times. This can be accomplished by introducing
the partition function, which we denote by
$Z(\tau)=\text{Tr}\left({\rho}(\tau)\right)$. We denote the normalized state
by $\overline{\rho}(\tau)\equiv\frac{{\rho}(\tau)}{Z(\tau)}$. With this
normalization, we can express the expectation value of an observable, $O$,
under the state $\overline{\rho}(\tau)$ as

$\langle
O\rangle(\tau)=\frac{1}{2^{m}Z(\tau)}\sum_{j=0}^{j=m}{m\choose{j}}\langle\psi_{2j-m}|O|\psi_{m-2j}\rangle,$
(20)

where

$Z(\tau)=\frac{1}{2^{m}}\sum_{j=0}^{j=m}{m\choose{j}}\langle\psi_{2j-m}|\psi_{m-2j}\rangle,$
(21)

noting that when substituting Eq. (21) into Eq. (20), the factors of $2^{m}$
cancel.

Equations (20) and (21) illustrate that an observable expectation value with
respect to the imaginary time evolved state $\overline{\rho}(\tau)$ is
expressible as a weighted sum of overlaps between states that have been
evolved unitarily in opposite directions in real time. In order to estimate
the expectation value of an observable under the ground state,
$|E_{\ell}\rangle$, of $H^{2}$, one could utilize Eq. (20) carried out to
large $\tau$, in line with Eq. (9), where $\rho(0)$ can be any state with
nonzero support on $|E_{\ell}\rangle$.

### Large $m$ limit

In settings where many time steps are needed to achieve convergence (i.e., the
large $m$ limit), evaluating Eqs. (20) and (21) can be computationally
challenging due to the presence of the binomial coefficients ${{m}\choose j}$,
whose magnitudes increase very quickly with $m$. As a practical matter, this
challenge can be addressed by first shifting the summation index to make the
symmetric nature of the sum explicit, i.e., replacing
$\sum_{j=0}^{j=m}{m\choose{j}}$ by
$\sum_{j=-\frac{m}{2}}^{j=\frac{m}{2}}{m\choose{\frac{m}{2}}+j}$ in Eqs. (20)
and (21), and then considering the large $m$ asymptotic for the binomial
coefficient

${m\choose m/2+j}\sim\frac{2^{m}}{\sqrt{m\pi/2}}e^{-2j^{2}/m}.$ (22)

Henceforth, we assume $m$ to be even. This expression is readily derivable
from the Stirling approximation [54], with an error on the order
$\mathcal{O}(j^{3}m^{-2})$, although as a practical matter this error becomes
negligible even for relatively small values of $m$. Substituting Eq. (22) into
Eqs. (20) and (21) yields

$\displaystyle\tilde{\langle O\rangle}(\tau)$
$\displaystyle=\frac{1}{\tilde{Z}(\tau)}\sum_{j=0}^{j=\frac{m}{2}}\frac{{e}^{-\frac{2j^{2}}{m}}}{\sqrt{m\pi/2}}\langle\psi_{2j}|O|\psi_{-2j}\rangle+\text{c.c.}$
(23)

and

$\displaystyle\tilde{Z}(\tau)$
$\displaystyle=\frac{1}{\sqrt{m\pi/2}}\sum_{j=0}^{j=\frac{m}{2}}e^{-\frac{2j^{2}}{m}}\langle\psi_{2j}|\psi_{-2j}\rangle+\text{c.c.},$
(24)

where we use the tilde to denote that the large $m$ approximation has been
used. Equations (23) and (24) can be used as alternative equations when $m$
becomes large.

### Maximally mixed initial condition

Finally, in the case where the initial condition is given by the identity
$\rho(0)={I}$, a number of further results can be inferred from ITQDE. This
initial condition corresponds to the maximally mixed or infinite temperature
state. Beginning from this condition means that at all $\tau$ $H^{2}$ and
$\rho$ necessarily commute, simplifying Eq. (7) such that we obtain at time
$\tau$ the Gaussian state

$\displaystyle{\rho}(\tau)=e^{-\tau H^{2}},$ (25)

which is equivalent to a thermal (Gibbs) state with an effective Hamiltonian
$H^{2}$ at inverse temperature $\beta=\tau$ [55, 56]. This illustrates that in
this setting, evolving $\rho(0)\propto{1}$ under Eq. (7) is functionally
equivalent to imaginary time propagation under $H^{2}$. This suggests that it
is desirable to consider this initial condition, given its natural
relationship to thermal states.

Furthermore, given $\rho(0)=I$, the final (unnormalized) state in Eq. (19)
becomes

$\rho(\tau)=\frac{1}{2^{m}}\sum_{j=0}^{j=m}{{m}\choose{j}}e^{-i(2j-m)\sqrt{2\Delta\tau}H}$
(26)

where the approximation error between Eq. (26) and the Gaussian state solution
given in Eq. (25) is also $\mathcal{O}(m\Delta\tau^{2})$. In this setting, Eq.
(20) can be used to estimate the expectation values of observables under
thermal states associated with inverse temperature $\beta=\tau$.

##  III ITQDE for spectral calculations

We now detail how the ITQDE correspondence can be used to perform spectral
calculations. We first observe that when $O=H^{2}$, evaluating Eq. 20 for
large $\tau$ would output an estimate of its ground state energy
$E_{\ell}^{2}$. Similarly, eigenenergies can be estimated with respect to the
thermal state described by Eq. (25). To calculate these at some desired
inverse temperature $\beta$ one need only implement Eq. (20) to $\tau=\beta$
from an initial condition of $\rho(0)=\frac{1}{Z(0)}I$.

It is however possible to go further, and exploit ITQDE to perform full
spectral calculations that additionally resolve excited states. To demonstrate
this, we first note that for a $d$-dimensional system, there are $2^{d}$
possible Hamiltonians $H_{\kappa}$ that square to the same $H^{2}$. In other
words, for all Hamiltonians $H_{\kappa}$, parameterized by length-$d$ binary
strings $\kappa$ with elements $\kappa_{j}\in\\{0,1\\}$ such that

$H_{\kappa}=\sum_{j=0}^{d-1}(-1)^{\kappa_{j}}E_{j}|E_{j}\rangle\langle
E_{j}|,$ (27)

we have that the squares of all these $H_{\kappa}$ are identically given by

$H^{2}=\sum_{j=0}^{d-1}E_{j}^{2}|E_{j}\rangle\langle E_{j}|.$ (28)

This reinforces that as $\tau\rightarrow\infty$, $\rho(\tau)$ will select the
eigenstate of $H$ with the smallest squared eigenvalue per Eq. (9). We can
take advantage of this behavior by introducing a shift $\lambda$ into the
Hamiltonian such that

$H^{(\lambda)}=H+\lambda I.$ (29)

This shift allows one to directly tune which eigenstate has the lowest squared
eigenvalue, and accordingly, which eigenstate $\rho^{(\lambda)}(\tau)$
approaches as the system is evolved to higher $\tau$ [57]. Importantly, the
effect of this shift on $\mathcal{U}$ is simply to introduce a global phase of
$e^{-i\sqrt{\frac{\Delta\tau}{2}}\lambda}$, i.e., such that
$\mathcal{U}^{(\lambda)}\equiv
e^{-i\sqrt{\frac{\Delta\tau}{2}}H^{(\lambda)}}=e^{-i\sqrt{\frac{\Delta\tau}{2}}\lambda}e^{-i\sqrt{\frac{\Delta\tau}{2}}H}$.
Given this, we may describe the expectation value of an observable $O$ with
respect to $\rho^{(\lambda)}(\tau)$ as

$\displaystyle\langle O^{(\lambda)}\rangle(\tau)$
$\displaystyle=\frac{1}{Z^{(\lambda)}(\tau)}\sum_{j=0}^{\frac{m}{2}}\frac{e^{2i\lambda
j\sqrt{\frac{\Delta\tau}{2}}}}{2^{m}}{{m}\choose{\frac{m}{2}+j}}\langle\psi_{2j}|O|\psi_{-2j}\rangle$
(30) $\displaystyle\quad+\text{c.c.},$

where

$Z^{(\lambda)}(\tau)=\frac{1}{2^{m}}\sum_{j=0}^{\frac{m}{2}}e^{2i\lambda
j\sqrt{\frac{\Delta\tau}{2}}}{{m}\choose{\frac{m}{2}+j}}\langle\psi_{2j}|\psi_{-2j}\rangle+\text{c.c.},$
(31)

and we have exploited the fact that each overlap has a conjugate partner to
reduce the number of terms in the sums from $m$ to $\frac{m}{2}$.

Meanwhile, in the large $m$ limit we obtain

$\displaystyle\langle\tilde{O}^{(\lambda)}\rangle(\tau)$
$\displaystyle=\frac{1}{\tilde{Z}^{(\lambda)}(\tau)}\sum_{j=0}^{j=\frac{m}{2}}\frac{{e}^{-\frac{2j^{2}}{m}+2i\lambda
j\sqrt{\frac{\Delta\tau}{2}}}}{\sqrt{m\pi/2}}\langle\psi_{2j}|O|\psi_{-2j}\rangle$
(32) $\displaystyle\quad+\text{c.c.}$

and

$\displaystyle\tilde{Z}^{(\lambda)}(\tau)$
$\displaystyle=\frac{1}{\sqrt{m\pi/2}}\sum_{j=0}^{j=\frac{m}{2}}e^{-\frac{2j^{2}}{m}+2i\lambda
j\sqrt{\frac{\Delta\tau}{2}}}\langle\psi_{2j}|\psi_{-2j}\rangle+\text{c.c.},$
(33)

by substituting Eq. (22) into Eqs. (30) and (31).

Crucially, when $O=H$, varying the value of $\lambda$ allows us to calculate
the entire energy spectrum by simply measuring the overlaps
$\langle\psi_{2j}|H|\psi_{-2j}\rangle$ and
$\langle\psi_{2j}|\psi_{-2j}\rangle$, $j=0,\cdots,m/2$ corresponding to a
single trajectory in $\tau$, and then evaluating Eq. (30) for different values
of $\lambda$ as a numerically trivial post-processing step. We also note
parenthetically here that for the ground state specifically $\lambda$ can also
be used to accelerate convergence. More details regarding this use are located
in App. A.

Somewhat amusingly (and certainly unintentionally), ITQDE implements what has
been a long-standing proposal in the context of Monte-Carlo methods [58].
Namely the spectral sweeping enabled by sampling $H^{2}$. Previous efforts in
this context have however foundered on the fact that $H^{2}$ is not only far
more singular than $H$ for any realistic Hamiltonian, but cannot guarantee
non-negative weights on paths [59]. The fact that ITQDE never actually employs
$H^{2}$ in its calculation of spectra raises the possibility that such issues
may also be avoided in a Monte-Carlo context.

Finally, we note that ITQDE’s ability to resolve spectra is a direct
consequence of the imaginary time evolution it implements generating a
Gaussian, rather than thermal (or Gibbs) state. In the case that one seeks to
obtain the thermal state directly via ITQDE, it is necessary to obtain the
‘square root’ of the Hamiltonian $\sqrt{H}$. While we do not address this in
the present work, such an operator will always exist for positive semi-
definite $H$ (which can itself be guaranteed with a suitable shift of the
spectrum), and quantum algorithmic techniques have already been proposed for
the calculation of $\sqrt{H}$ [60].

##  IV Thermodynamic applications of ITQDE

Beyond spectral calculations, the character of the ITQDE correspondence has a
number of interesting properties that are worth commenting on. In particular,
the presence of forward and reverse trajectories in the correspondence ought
not to be ignored, given the ubiquity of this feature in both quantum [61] and
stochastic [62] thermodynamics. In the latter case fluctuation theorems are
derived from considerations of the probability of trajectories against their
reversed counterpart [63], suggesting similar results may be obtained via the
methods outlined in previous sections. Indeed, if one were to employ the large
$m$ approximation according to Eq. (22), the result

$e^{-\tau{H}^{2}}=\frac{1}{\sqrt{\pi}}\sum_{j=-\frac{m}{2}}^{j=\frac{m}{2}}\sqrt{\frac{2}{m}}e^{-\frac{2j^{2}}{m}}e^{-2ij\sqrt{2\Delta\tau}{H}}$
(34)

can be used to take the expectation of the operator ${\rm e}^{-\tau H^{2}}$
with respect to $\ket{\psi_{0}}$. From this we directly obtain (reiterating
the implicit $\mathcal{O}(m\Delta\tau^{2})$ error)

$\langle{\rm e}^{-\tau
H^{2}}\rangle=\sqrt{\frac{2}{m\pi}}\sum_{j=0}^{j=\frac{m}{2}}e^{-\frac{2j^{2}}{m}}\langle\psi_{2j}|\psi_{-2j}\rangle+\text{c.c.},$
(35)

a result that bears a remarkable similarity to the Crooks fluctuation theorem
[64], where the work distribution function of a stochastic process can be
expressed as a ratio of the probabilities for forward and reversed processes.
In the current instance, we find the expected Hamiltonian Gaussian
distribution of a state is expressible in terms of an exponentially weighted
sum of overlaps between forward and backward evolutions from that state.
Indeed, the ITQDE derived expression is a a manifestly quantum one, insofar as
it utilises overlaps rather than probabilistic ratios, and that it links to
distributions over $H$, rather than the non-observable work variable.

More broadly, it is interesting to observe the tension between microscopic
reversibility and thermal equilibration manifest itself here, in the sense
that ITQDE suggests that equilibrium properties are the product of both
forwards and reversed trajectories, rendering the question of directionality
redundant. The direction of the arrow of time hardly matters when its
consequences are manifestly symmetric. Moreover, given the only sense in which
a time increment appears is ${\Delta\tau}$ \- which in the real time
propagation occurs under a square root, the notion of time reversal is most
naturally captured by the equivalent transformation of $H\to-H$. Note however
that when we consider this transformation in the case of Gibbs states - the
physical thermal distribution - the necessity of employing $\sqrt{H}$ in the
propagator means this reversal symmetry is broken. This fact, to the
philosophically inclined, may be of some interest.

It is also possible to invert the ITQDE correspondence. First, taking the same
large $m$ approximation of Eq. (34) and inserting $\lambda$, we have:

${\rm
e}^{-\tau(H-\lambda)^{2}}=\frac{1}{\sqrt{\pi}}\sum_{j=-\frac{m}{2}}^{j=\frac{m}{2}}\sqrt{\frac{2}{m}}{\rm
e}^{-\frac{2j^{2}}{m}}{\rm e}^{-2ij\sqrt{2{\rm d}\tau}H}{\rm
e}^{-2ij\sqrt{2{\rm d}\tau}\lambda}.$ (36)

Following this, we note the orthogonality relationship for complex
exponentials is given by:

$\sum_{\lambda=0}^{\lambda=K-1}{\rm
e}^{\frac{2i\pi\lambda}{K}(j-k)}=K\delta_{jk}.$ (37)

Let us say that $\frac{m}{2}$ is a square number, such that
$K=\sqrt{\frac{m}{2}}$ is an integer. This also means that $K\sqrt{2{\rm
d}\tau}=\sqrt{\tau}$. Then, applying the complex exponential and summing over
$\lambda$ on both sides, we have

$\sum^{\lambda=K-1}_{\lambda=0}{\rm
e}^{-\tau(H-\frac{\pi}{\sqrt{\tau}}\lambda)^{2}}{\rm e}^{\frac{2i\pi\lambda
k}{K}}=\frac{1}{\sqrt{\pi}}{\rm e}^{-\frac{2k^{2}}{m}}{\rm e}^{-2ik\sqrt{2{\rm
d}\tau}H}.$ (38)

By defining $t=2k\sqrt{2{\rm d}\tau}$ and rearranging, we obtain an expression
for the propagator:

${\rm e}^{-itH}=\sqrt{\pi}{\rm
e}^{\frac{t^{2}}{4\tau}}\sum^{\lambda=K-1}_{\lambda=0}{\rm
e}^{-\tau(H-\frac{\pi}{\sqrt{\tau}}\lambda)^{2}}{\rm e}^{\frac{2i\pi\lambda
t}{\sqrt{\tau}}}.$ (39)

The fundamental motivation for such a representation is that it allows one to
cast dynamical processes purely in terms of equilibrium information, with time
acting only as a weighted phase on shifted Gaussian distributions. This
expression may therefore likely be useful in obtaining alternate formulations
of dynamical quantities such as a system’s dynamic response [65] or Out-Of-
Time-Order Correlators (OTOCs) [66]. This latter quantity is a measure of
information scrambling, and recent proposals have demonstrated its intrinsic
relationship to thermodynamics [67]. In this sense Eq.(39) offers a natural
route to further exploration of this connection. Moreover, in the context of
both OTOCs and Lyapunov exponents, a $t^{2}$ exponent often emerges, identical
to that obtained in the inverted ITQDE correspondence above. Such coincidences
are suggestive of deeper links that might be uncovered with a systematic
application of ITQDE to these topics.

Lastly, if we consider the continuous limit of Eq. (26), a further connection
to open systems theory is revealed. Returning to Eq. (34) and defining
$\delta=\sqrt{\frac{2}{m}}$ and $x=j\delta$, the right-hand-side of Eq. (34)
may be expressed as

$e^{-\tau{H}^{2}}=\frac{1}{\sqrt{\pi}}\sum_{j=-\frac{m}{2}}^{j=\frac{m}{2}}\delta
e^{-x^{2}}e^{-2ix\sqrt{\tau}{H}}+\mathcal{O}(\delta^{2}),$ (40)

such that in the limit $\delta\to 0$ we recover the celebrated Hubbard-
Stratonovich (HS) transformation [68, 69, 70, 71]:

$\displaystyle e^{-\tau
H^{2}}=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-x^{2}}e^{-2ix\sqrt{\tau}H}{\rm
d}x.$ (41)

This result is usually understood as a correspondence linking a deterministic
quadratic potential to a stochastic linear one. Here, however, it emerges in a
generative manner from ITQDE as the continuous limit of infinitesimal dynamics
where the state space and its dual are oppositely evolved. Based on this
limit, it is in fact possible to infer a generalisation of the HS
transformation by repeating the same limiting procedure, but beginning from a
general initial state. Setting $\ket{\psi_{0}}=\sum_{j}a_{j}\ket{E_{j}}$, we
find

$\displaystyle{\rm e}^{-\tau(E_{k}+E_{l})^{2}/4}\ket{E_{k}}\bra{E_{l}}$
$\displaystyle\quad=\frac{1}{2^{n}}\sum_{j=-\frac{n}{2}}^{j=\frac{n}{2}}{n\choose\frac{n}{2}+{j}}{\rm
e}^{-2ij\sqrt{2{\rm d}\tau}(E_{k}+E_{l})}\ket{E_{k}}\bra{E_{l}}.$ (42)

If we define the operator

$\rho_{\rm mix}=\sum_{kl}a_{k}a_{l}{\rm
e}^{-\tau(E_{k}E_{l})/2}\ket{E_{k}}\bra{E_{l}},$ (43)

then the previous $\delta\to 0$ limit then yields

$\displaystyle e^{-\tau H^{2}/2}\rho_{\rm mix}e^{-\tau H^{2}/2}$
$\displaystyle\quad=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}e^{-x^{2}}e^{-ix\sqrt{\tau}H}\ket{\psi_{0}}\bra{\psi_{0}}e^{-ix\sqrt{\tau}H}{\rm
d}x.$ (44)

##  V Numerical illustrations

Figure 2: Results of spectral calculations performed for an eight site Ising
Hamiltonian using $\Delta\tau=0.4\times 10^{-5}$ and $m=1000$. Dashed lines
indicate the eigenenergies of $H_{\mathrm{TFIM}}$. Solid lines indicate the
convergence of the spectral calculation to these eigenenergies. For a), we
plot the calculated energy as $\lambda$ against propagation in imaginary time
$\tau$ and find that our method converges to successively lower eigenstates of
$H_{\mathrm{TFIM}}$ before reaching the ground state $E_{0}$. For b), we plot
the change in the steady state expectation value of $H_{\mathrm{TFIM}}$ as
$\lambda$ is varied, illustrating that a complete course-grained spectrum is
obtainable using only expectation values estimated from a single trajectory
only.

To demonstrate the capability of Eq. (30) to resolve the spectrum of $H$, we
now illustrate an application of this procedure to the 1D transverse field
Ising model, which is a simple model composed of nearest-neighbor coupled
spins on a lattice in the presence of an external, uniform magnetic field. It
has been studied extensively over the years in the contexts of quantum phase
transitions [72, 73], quantum spin glasses [74, 75], and the quantum annealing
process [76, 77, 78, 79] among others. It is also a useful benchmark problem
for developing and testing quantum algorithms, as the model is exactly
solvable [80, 72], and spins map naturally to qubits. The Hamiltonian is given
by

$\displaystyle
H_{\mathrm{TFIM}}=-J\sum_{l=0}^{L-1}Z_{l}Z_{l+1}-h\sum_{l=0}^{L-1}X_{l},$ (45)

where $Z_{l}$ and $X_{l}$ denote Pauli-$Z$ and Pauli-$X$ matrices,
respectively, that act on the spin occupying lattice site $l$, $J$
characterizes the strength of the nearest-neighbor interactions, $h$ is the
magnetic field strength, and $L$ is the total number of lattice sites. Here,
we consider a model with $J=1$, $h=2$, $L=8$, and open boundary conditions. We
use Eqs. (30) and (31) to calculate the expectation value $\langle
H_{\mathrm{TFIM}}^{(\lambda)}\rangle(\tau)$ for different values of $\lambda$,
which are plotted as a function of $\tau$ in Fig. 2(a). These results show
that changing the value of $\lambda$ causes $\langle
H_{\mathrm{TFIM}}^{(\lambda)}\rangle(\tau)$ to converge to different energies,
corresponding in this case to the energies associated with each of the
different bands of eigenstates. This premise is further reinforced in Fig.
2(b), where we plot the steady state expectation value $\langle
H_{\mathrm{TFIM}}^{(\lambda)}\rangle(\tau\rightarrow\infty)=\langle
H_{\mathrm{TFIM}}\rangle(\tau\rightarrow\infty)+\lambda$ against $\lambda$,
from which the complete spectrum is inferred.

Figure 3: Relative error calculated according to Eq. (46) in spectral
calculations performed using the large $m$ approximation in Eqs. (32) and (33)
compared to the results obtained using Eqs. (30) and (31) for the 8 site 1D
transverse field Ising model. In both cases, $\Delta\tau=0.4\times 10^{-5}$
and $m=1000$ Figure 4: Results of spectral calculations performed for a four
site, 2D Fermi-Hubbard model at half-filling with $d\tau=0.003$ and $m=1500$.
Dashed lines indicate the eigenenergies of $H_{\mathrm{FH}}$, and solid lines
show the convergence of $\langle H_{\mathrm{FH}}^{(\lambda)}\rangle(\tau)$ to
these eigenenergies. For a), we plot the calculated energy for various
$\lambda$ against $\tau$, and observe good convergence to successively lower
eigenstates of $H_{\mathrm{FH}}$ before reaching the ground state, $E_{0}$.
For b), we plot the change in the steady state expectation value, i.e.,
$\langle H_{\mathrm{FH}}^{(\lambda)}\rangle(\tau\rightarrow\infty)$ as
$\lambda$ is varied, illustrating that a coarse-grained spectrum can be
obtained. Figure 5: Relative error in spectral calculations of the steady
state expectation values for the eigenenergies of the 2D Fermi-Hubbard model
at half-filling for 2, 4, 6, and 8 sites. Each eigenenergy error is plotted
symmetrically about the spectral centre, which is aligned on the figure with
the number of sites that spectrum corresponds to. All calculations were
performed using a step size of $\Delta\tau=0.003$ for $m=1500$ steps to a
final imaginary time $\tau=5$ with the initial condition
$\rho_{\mathrm{FD}}(0)=I$.

We demonstrate the use of the large $m$ approximation by comparing the results
obtained using Eqs. (32) and (33) to the ones obtained using Eqs. (30) and
(31) for this 8 site transverse field Ising model. The relative difference

$\displaystyle\zeta=\frac{|\langle O^{(\lambda)}\rangle-\tilde{\langle
O^{(\lambda)}\rangle}|}{|\langle O_{g}\rangle|},$ (46)

for each eigenenergy of $H_{\mathrm{TFIM}}$ is shown in Fig. 3. Here we use
$\langle O_{g}\rangle$ as the minimum eigenenergy of $H_{\mathrm{TFIM}}$ to
avoid dividing by small numbers when the calculated eigenvalue is near 0. This
demonstrates that the large $m$ approximation obtains results that are similar
to the direct approach using the binomial coefficients.

We also apply the large $m$ approximation to the 2D Fermi-Hubbard model. The
Fermi-Hubbard model is a simple lattice model that aims to capture key aspects
of strongly correlated fermionic systems, and has been used to study phase
transitions [81], superconductivity [82], and magnetism [83]. The Hamiltonian
is given by

$H_{\mathrm{FH}}=\sum_{\langle
i,j\rangle,\sigma}t_{ij}(c_{i\sigma}^{\dagger}c_{j\sigma}+\text{h.c.})+U\sum_{j=0}^{L-1}n_{j\uparrow}n_{j\downarrow}-\mu\sum_{j=0}^{L-1}n_{j\sigma},$
(47)

where $t_{ij}$ denotes the tunneling amplitude between lattice sites $i$ and
$j$, $c_{j\sigma}^{\dagger}$ and $c_{j\sigma}$ are the creation and
annihilation operators associated with a fermion of spin $\sigma$ at lattice
site $j$, respectively, $U$ characterizes the on-site interaction of fermions,
$n_{j\sigma}=c^{\dagger}_{j\sigma}c_{j\sigma}$ is the occupation number
operator, $\mu$ is the chemical potential, and $L$ is the number of lattice
sites.

In Fig. 4(a), we plot results showing how the expectation value $\langle
H_{\mathrm{FH}}^{(\lambda)}\rangle(\tau)$ varies with $\tau$ and $\lambda$
when the model parameters are set according to $t_{ij}=t=-1$ for all nearest-
neighbor pairs, $U=2$, $\mu=0.5$, $L=4$ (corresponding to a $2\times 2$
lattice), and periodic boundary conditions. These results again demonstrate
that changing $\lambda$ causes $\langle
H_{\mathrm{FH}}^{(\lambda)}\rangle(\tau)$ to converge to different energies in
the spectrum of $H_{\mathrm{FH}}$ and that reasonable results are obtained
even when using the large $m$ approximation. The steady state expectation
value $\langle H_{\mathrm{FM}}^{(\lambda)}\rangle(\tau\rightarrow\infty)$ of
this model for varying $\lambda$ is plotted in Fig. 4(b) and further
demonstrates that changing $\lambda$ can be used to infer the energetic
spectrum.

To examine the performance of our method as the number of sites in the 2D
Fermi-Hubbard model increases, we calculate the relative error in the
calculation of each eigenenergy compared to the exact values obtained by
diagonalizing $H_{\mathrm{FM}}$. This is plotted against the number of sites
in the lattice in Fig. 5. We observe from these results that our method is
able to accurately calculate the energy eigenvalues of the 2D Fermi-Hubbard
model up to 8 sites and that the accuracy is not dependent on the number of
sites in the lattice. The exception is the 2 site case, where much more
accurate results are obtained. This is most likely due to the fact that in the
case of 2 sites the model is effectively its much simpler one dimensional
equivalent.

All of the results discussed in this section were obtained using the initial
condition $\rho(0)=I$. In App. B, we discuss the effects of choosing a single
pure state as the initial condition on the spectral calculation results.

##  VI Quantum algorithm for spectral calculations using ITQDE

Figure 6: Quantum circuits for performing Hadamard tests to estimate a)
$\mathrm{Re}(\langle\psi_{-2j}|P_{\alpha}|\psi_{2j}\rangle)$ and b)
$\mathrm{Im}(\langle\psi_{-2j}|P_{\alpha}|\psi_{2j}\rangle)$ for use in Eqs.
(32) and (33). We denote by ${P}_{\alpha}$ the $\alpha$-th Pauli string in the
decomposition of the observable $O$, per Eq. (51).
$H=\frac{1}{\sqrt{2}}\bigl{(}\begin{smallmatrix}1&1\\\
1&-1\end{smallmatrix}\bigr{)}$ and $S=\bigl{(}\begin{smallmatrix}1&0\\\
0&i\end{smallmatrix}\bigr{)}$ denote the Hadamard and phase gate,
respectively.

We now present a proof-of-concept application of ITQDE to quantum computation.
In particular, we implement ITQDE as a quantum algorithm, where a quantum
computer is utilized to estimate the overlaps (i.e., of the form
$\langle\psi_{-2j}|\psi_{2j}\rangle$ and
$\langle\psi_{-2j}|O|\psi_{2j}\rangle$) that appear in Eqs. (30)-(33). We
focus this section on how to construct a quantum algorithm for evaluating Eqs.
(30)-(33) that is independent of the initial state $|\psi_{0}\rangle$, since
in practice, $|\psi_{0}\rangle$ can be selected in a problem-dependent manner.

There are a number of methods by which the relevant overlaps might be
calculated in a quantum computer (see for example App. D), but here we
consider the most conceptually direct route provided by the Hadamard test
[84]. This requires only a single ancilla qubit, and is detailed fully in App.
C. The Hadamaed test functions by initializing the ancilla qubit in the state
$|0\rangle$ and the primary system register in a state $|\varphi\rangle$,
followed by Hadamard gates on the ancilla qubit, interleaved with a unitary
operation $\mathcal{W}$ that is applied to the primary system register and
controlled on the state of the ancilla qubit, and culminating with a final
ancilla qubit measurement [see Fig. 6(a)]. Repeating this procedure many times
and averaging the results of the final ancilla measurement then outputs an
estimate for the real part of the overlap between $|\varphi\rangle$ and
$\mathcal{W}|\varphi\rangle$ in terms of the probability of measuring 0,
$p_{0}$ on the ancilla qubit, such that

$\text{Re}\left(\langle\varphi|\mathcal{W}|\varphi\rangle\right)=2p_{0}-1.$
(48)

A modified version of the Hadamard test, i.e., one that incorporates an
additional phase gate on the ancilla qubit, can be used to obtain the
imaginary part of the overlap in an analogous manner [see Fig. 6(b)].

In order to utilize the Hadamard test in our setting, we now discuss
approaches for using it to estimate $\langle\psi_{-2j}|\psi_{2j}\rangle$ and
$\langle\psi_{-2j}|O|\psi_{2j}\rangle$, beginning with the former, which can
be expressed as

$\langle\psi_{-2j}|\psi_{2j}\rangle=\langle\psi_{0}|\mathcal{U}^{4j}|\psi_{0}\rangle.$
(49)

Hadamard tests can be performed to estimate the real and imaginary parts of
the right-hand-side of Eq. (49) by associating
$\mathcal{W}=\mathcal{U}^{4j}=e^{-i\sqrt{2j\Delta\tau}H}$ and
$|\varphi\rangle=|\psi_{0}\rangle$ in Eq. (48).

In order to estimate overlaps of the second form,
$\langle\psi_{-2j}|O|\psi_{2j}\rangle$, we first expand $O$ in the Pauli
operator basis according to

$O=\sum_{\alpha=1}^{M}c_{\alpha}P_{\alpha},$ (50)

where $c_{\alpha}$ and $P_{\alpha}$ denote a real coefficient and Pauli basis
operator, respectively. We assume that $M=\text{poly}(n)$, noting that when
$O$ is taken to be a $k$-local observable, this assumption holds. Importantly,
because the Pauli operators are unitary, this expansion allows us to then
estimate $\langle\psi_{-2j}|O|\psi_{2j}\rangle$ by taking the sum of a linear
combination of constituent overlaps, each evaluated using separate Hadamard
tests, according to

$\langle\psi_{-2j}|O|\psi_{2j}\rangle=\sum_{\alpha=1}^{M}c_{\alpha}\langle\psi_{0}|\mathcal{U}^{2j}P_{\alpha}\mathcal{U}^{2j}|\psi_{0}\rangle,$
(51)

where for the $\alpha$-th Hadamard test, we have
$\mathcal{W}=\mathcal{U}^{2j}P_{\alpha}\mathcal{U}^{2j}$, and $M$ Hadamard
tests are needed in total to evaluate Eq. (51). Using this approach to
estimate the real and imaginary parts of all $\frac{m}{2}+1$ overlaps in Eqs.
(32) and (33) then requires $(m+2)(M+1)=\mathcal{O}(mM)$ separate circuits
(each repeated sufficiently many times to ensure convergence of the estimate).

Given this Hadamard test formulation, a method for implementing the (real)
time evolution $\mathcal{U}^{4j}=e^{-i\sqrt{2j\Delta\tau}H}$ on the primary
system register is then all that is additionally required to complete the
specification of the quantum algorithm. It is worth noting that in order to
simulate imaginary time evolution by an imaginary time $\tau$ under a
Hamiltonian $H^{2}$, we only need real time evolution under $H$ by time
$\sqrt{2\Delta\tau}$. For physically-realistic Hamiltonians, it is known that
$\mathcal{U}^{4j}$ can be implemented _efficiently_ on a quantum computer [85]
via a variety of different Hamiltonian simulation algorithms that proceed in a
time that scales as $\text{poly}(n)$ on $n$ qubits.

### Quantum device implementations

Figure 7: Spectral calculation results for a 2-site transverse field Ising
model with $J=1$ and $h=2$, obtained via implementation of ITQDE on an IBM
superconducting qubit processor (ibm_brisbane). Here we use $\Delta\tau=0.01$
and $m=100$. In a), the expectation values $\langle
H_{TFIM}^{(\lambda)}\rangle(\tau)$ output from Eqs. (32) and (33), with
overlaps obtained from ibm_brisbane, are shown (dashed curves) against the
true eigenenergies (solid lines of corresponding color). Each overlap is
obtained using 4000 samples. In b), the change in the steady state expectation
of $H_{TFIM}$ as $\lambda$ changes is shown, illustrating that the coarse-
grained spectrum can still be inferred when performing ITQDE on quantum
hardware. In both figures, the errorbars denote 95% confidence.

Finally, we present the results from implementing ITQDE directly on quantum
devices. For this, we utilize the IBM superconducting qubit processor
ibm_brisbane, along with IBM’s sampler primitive [86]. Our first
implementation considers the transverse field Ising model with $J=1$, $h=2$,
$L=2$, and the results are presented in Fig. 7. In this demonstration, the
quantum dynamical emulation procedure is repeated for 104 different initial
states, prepared by applying random Clifford operations to the initial
$|00\rangle$ state, and averaging the overlaps at each step. This averaging
procedure (over initial Clifford realizations) aims to produce statistics that
converge, with increasing realizations, to those that would be obtained with
an initialization in the maximally mixed state. We expect that the convergence
towards the middle energies could be improved by sampling initial states from
additional random Cliffords.

##  VII  Outlook

In this work we have proposed a method for constructing the solutions of non-
unitary dynamics from purely unitary operations, which we term Quantum
Dynamical Emulation. This was applied in the context of imaginary time to
derive the ITQDE correspondence between real and imaginary time evolution.
This in effect provides a unitary representation of a non-unitary evolution,
which has been exploited here to develop a quantum algorithm for spectral
calculations. The results presented here represent a first application of
quantum dynamical emulation; however, this technique is ripe for further
refinement in terms of approximations to improve its efficiency, in the scope
of its application, and in understanding its limitations. ITQDE is after all
an approximate equivalence based on a discretisation, and it is clear that the
success of this approximation will depend upon the interplay of free
parameters such as $\Delta\tau$ and $\lambda$, as well as the type of
Hamiltonian being studied. There is after all no such thing as a free lunch,
and a likely approach to refining this technique will be understanding where
the computational complexity implied by the ability to calculate spectra lies.

For example, the error in the present approximation is of order $\Delta\tau$,
and while in the cases considered thus far this has not proved problematic, an
immediate avenue for future investigation will be how to improve upon the
order of this error. This could be achieved either by the application of
RK4-like techniques, of from instead working backwards from the continuous
limit of the H-S transform. In such a case, the use of techniques such as
Gaussian quadrature may be employed to obtain highly accurate discretisations
while requiring relatively few points. This procedure is sketched in App. E,
together with an error bound for such an approach. This would not only have
the potential to improve accuracy, but do so while reducing the required
number of overlaps for effective implementation of ITQDE.

While the majority of the calculations presented in this work use the large
$m$ approximations in Eqs. (32) and (33) to avoid calculating very large
factorials, the exponential weighting means that most of the weight in the
calculation is focused around the initial state and short-time overlaps.
Embracing this feature may allow for truncated expressions with the potential
to extrapolate out to much later $\tau$ than would otherwise be possible. This
would allow for calculations with much shallower circuit depths than direct
implementation would require.

As a proof-of-concept, ITQDE has been successfully applied to the calculation
of spectral calculations on current quantum hardware. Prioritising conceptual
clarity over performance, this has been achieved via the Hadamard test
methodology for overlap measurement. This method’s use of controlled unitaries
however means that it is a technique best suited to fault-tolerant
computation. At present however quantum computing remains a NISQy business,
and for this reason interferometric protocols that forego such controlled
operations [87, 88, 89] are likely to be more suitable for the implementation
of ITQDE in present hardware. An immediate application of such an approach
would be in order to perform topological data analysis using recently proposed
imaginary time methods [90].

Beyond the near term, the asymptotic scaling of ITQDE when combined with post-
Trotter methods of implementation remains to be investigated. For example,
while in a naive Trotterisation of the propagator, the circuit depth will
scale with time, techniques such as Cartan decomposition [91] fix this depth
for all times. Furthermore proposals for the estimation of thermal states [60]
employ techniques such as block encodings based on Linear Combination of
Unitaries and the TTS algorithm [92], and these same techniques may also be
applied to future implementations of ITQDE. Given ITQDE’s ready adaptability
to both fault-tolerant and NISQ friendly methods of implementation, it is
likely to find application both in the present, and in a future where large-
scale, fault tolerant quantum computers become a reality. In fact, we expect
the question of efficient implementation in both the near and far term -
together with a comparative study to existing methods for imaginary time and
spectral calculations - to be addressed in forthcoming work.

More broadly, the fact that ITQDE reconstructs spectra using only measurement
overlaps means the method is platform agnostic, and can be applied not only to
a quantum circuit context, but any scenario in which such overlaps can be
measured. A future research question is therefore whether ITQDE can be
efficiently combined with other techniques (e.g. quantum state tomography),
such that these overlaps are explicitly recoverable from measurements in as
broad a context as possible. Moreover, as alluded to in Sec. V, the indirect
manner in which $H^{2}$ is introduced into calculations potentially
circumvents the issues which usually prevent the spectral sweeps demonstrated
here from being implemented in a Monte-Carlo setting. The fusion of these two
approaches to imaginary time simulation could therefore constitute a further
source of inspiration for extending ITQDE.

We have also provided a preliminary sketch of potential uses for ITQDE in the
context of stochastic and quantum thermodynamics. Given these results, it is
difficult to resist the impression of deep links between ITQDE and
thermodynamic theory. We expect that continued study of these connections have
the potential to yield both surprising and useful results. In fact, while the
work performed here has restricted itself to a quantum setting, Hilbert space
representations of classical dynamics [93, 94, 95, 96] \- including the
recently developed waveoperator representation [97] \- represent excellent
prospects for applying both the current work and its future extensions to the
domain of classical dynamics.

Finally, we reiterate that ITQDE is a proof of concept for the more general
procedure of QDE. The objective of unitary representations of non-unitary
dynamics has recently produced a number of proposals [98, 99], but while all
these approaches are animated by the same motive, the methods used to achieve
this stand in stark contrast. A natural question then for future work is to
understand the commonalities in these approaches, perhaps with an eventual
view towards their amalgamation. In the near term, the next implementation of
QDE is likely a form of Lindbladian dynamics based on Ref.[50], which could in
turn be applied to the realisation of recent proposals for error-mitigation in
quantum simulation [100]. More fundamentally, an open question for QDE is both
the limit of what exotic states might be expressible via an operator-valued
differential equation, and the degree to which QDE renders them realisable and
computable.

In closing this paper, we have attempted to highlight the many opportunities
and challenges promised by the present work. It should be considered as an
invitation to the broader physics community, with its diversity of interests
and expertise, to seize this opportunity to attack questions old and new with
a fresh perspective. We therefore close on a note we hope reflects the sense
of adventure we have felt, and continue to feel in this ongoing work. With
QDE, the imaginary has become reality.

###### Acknowledgements.

The authors would like to thank Alicia Magann and Christian Arenz for their
vital contributions, without which the present work would have been
impossible. J.M.L. is supported by Sandia National Laboratories’ Laboratory
Directed Research and Development Program. G.M. and D.I.B. were supported by
the Army Research Office (grant W911NF-23-1-0288; program manager Dr. James
Joseph). Sandia National Laboratories is a multimission laboratory managed and
operated by National Technology & Engineering Solutions of Sandia, LLC, a
wholly owned subsidiary of Honeywell International Inc., for the U.S.
Department of Energy’s National Nuclear Security Administration under contract
DE-NA0003525. This article has been authored by an employee of National
Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-
NA0003525 with the U.S. Department of Energy (DOE). The employee owns all
right, title and interest in and to the article and is solely responsible for
its contents. The United States Government retains and the publisher, by
accepting the article for publication, acknowledges that the United States
Government retains a non-exclusive, paid-up, irrevocable, world-wide license
to publish or reproduce the published form of this article or allow others to
do so, for United States Government purposes. The DOE will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-
plan. This paper describes objective technical results and analysis. Any
subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy, Army
Research Office, or the United States Government. SAND2024-02496O.

## References

  * Wick [1954] G. C. Wick, Properties of bethe-salpeter wave functions, Phys. Rev. 96, 1124 (1954).
  * Yang  _et al._ [2013] Y. Yang, H. Li, L. Sheng, R. Shen, D. N. Sheng, and D. Y. Xing, Topological phase transitions with and without energy gap closing, New Journal of Physics 15, 083042 (2013), publisher: IOP Publishing.
  * Hastings and Koma [2006] M. B. Hastings and T. Koma, Spectral Gap and Exponential Decay of Correlations, Communications in Mathematical Physics 265, 781 (2006).
  * Cubitt  _et al._ [2015] T. S. Cubitt, D. Perez-Garcia, and M. M. Wolf, Undecidability of the spectral gap, Nature 528, 207 (2015).
  * Han  _et al._ [2012] T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Fractionalized excitations in the spin-liquid state of a kagome-lattice antiferromagnet, Nature 492, 406 (2012), number: 7429 Publisher: Nature Publishing Group.
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010), number: 7286 Publisher: Nature Publishing Group.
  * Hartmann  _et al._ [2006] M. J. Hartmann, M. E. Reuter, and M. B. Plenio, Excitation and entanglement transfer versus spectral gap, New Journal of Physics 8, 94 (2006).
  * Jacquemin  _et al._ [2011] D. Jacquemin, B. Mennucci, and C. Adamo, Excited-state calculations with TD-DFT: from benchmarks to simulations in complex environments, Physical Chemistry Chemical Physics 13, 16987 (2011), publisher: Royal Society of Chemistry.
  * Ramos and Pavanello [2018] P. Ramos and M. Pavanello, Low-lying excited states by constrained DFT, The Journal of Chemical Physics 148, 144103 (2018), publisher: American Institute of Physics.
  * Ferré  _et al._ [2016] N. Ferré, M. Filatov, and M. Huix-Rotllant, eds.,  _Density-Functional Methods for Excited States_, Topics in Current Chemistry, Vol. 368 (Springer International Publishing, Cham, 2016).
  * Kaldor [1975] U. Kaldor, Degenerate many‐body perturbation theory: Excited states of H ${}_{\textrm{2}}$, The Journal of Chemical Physics 63, 2199 (1975).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Physical Review B 34, 5390 (1986), publisher: American Physical Society.
  * Hybertsen and Louie [1985] M. S. Hybertsen and S. G. Louie, First-Principles Theory of Quasiparticles: Calculation of Band Gaps in Semiconductors and Insulators, Physical Review Letters 55, 1418 (1985), publisher: American Physical Society.
  * Godby  _et al._ [1986] R. W. Godby, M. Schlüter, and L. J. Sham, Accurate Exchange-Correlation Potential for Silicon and Its Discontinuity on Addition of an Electron, Physical Review Letters 56, 2415 (1986), publisher: American Physical Society.
  * Godby  _et al._ [1988] R. W. Godby, M. Schlüter, and L. J. Sham, Self-energy operators and exchange-correlation potentials in semiconductors, Physical Review B 37, 10159 (1988), publisher: American Physical Society.
  * Hunt  _et al._ [2018] R. J. Hunt, M. Szyniszewski, G. I. Prayogo, R. Maezono, and N. D. Drummond, Quantum Monte Carlo calculations of energy gaps from first principles, Physical Review B 98, 075122 (2018).
  * Blume  _et al._ [1997] D. Blume, M. Lewerenz, P. Niyaz, and K. B. Whaley, Excited states by quantum Monte Carlo methods: Imaginary time evolution with projection operators, Physical Review E 55, 3664 (1997), publisher: American Physical Society.
  * Lüchow  _et al._ [2003] A. Lüchow, D. Neuhauser, J. Ka, R. Baer, J. Chen, and V. A. Mandelshtam, Computing Energy Levels by Inversion of Imaginary-Time Cross-Correlation Functions, The Journal of Physical Chemistry A 107, 7175 (2003).
  * Chou and Kouri [2013] C.-C. Chou and D. J. Kouri, Multidimensional Supersymmetric Quantum Mechanics: A Scalar Hamiltonian Approach to Excited States by the Imaginary Time Propagation Method, The Journal of Physical Chemistry A 117, 3449 (2013), publisher: American Chemical Society.
  * Joshi  _et al._ [2022] L. K. Joshi, A. Elben, A. Vikram, B. Vermersch, V. Galitski, and P. Zoller, Probing Many-Body Quantum Chaos with Quantum Simulators, Phys. Rev. X 12, 011018 (2022), publisher: American Physical Society.
  * Luukko and Räsänen [2013] P. J. J. Luukko and E. Räsänen, Imaginary time propagation code for large-scale two-dimensional eigenvalue problems in magnetic fields, Computer Physics Communications 184, 769 (2013).
  * Needs  _et al._ [2010] R. J. Needs, M. D. Towler, N. D. Drummond, and P. López Ríos, Continuum variational and diffusion quantum Monte Carlo calculations, Journal of Physics: Condensed Matter 22, 023201 (2010).
  * Lehtovaara  _et al._ [2007] L. Lehtovaara, J. Toivanen, and J. Eloranta, Solution of time-independent Schrödinger equation by the imaginary time propagation method, Journal of Computational Physics 221, 148 (2007).
  * Chin  _et al._ [2009] S. A. Chin, S. Janecek, and E. Krotscheck, Any order imaginary time propagation method for solving the Schrödinger equation, Chemical Physics Letters 470, 342 (2009).
  * Tasaki [1998] H. Tasaki, The Hubbard model - An introduction and selected rigorous results, J. Phys. Cond. Mat. 10, 4353 (1998).
  * Ren  _et al._ [2023] W. Ren, W. Fu, X. Wu, and J. Chen, Towards the ground state of molecules via diffusion Monte Carlo on neural networks, Nature Communications 14, 10.1038/s41467-023-37609-3 (2023), number: 1 Publisher: Nature Publishing Group.
  * Köster  _et al._ [2011] A. M. Köster, P. Calaminici, E. Orgaz, D. R. Roy, J. U. Reveles, and S. N. Khanna, On the Ground State of Pd13, Journal of the American Chemical Society 133, 12192 (2011), publisher: American Chemical Society.
  * Ceperley and Alder [1980] D. M. Ceperley and B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, Physical Review Letters 45, 566 (1980), publisher: American Physical Society.
  * Chen and Stoudenmire [2020] J. Chen and E. M. Stoudenmire, Hybrid purification and sampling approach for thermal quantum systems, Physical Review B 101, 195119 (2020), publisher: American Physical Society.
  * Sinha  _et al._ [2024] A. Sinha, M. M. Rams, and J. Dziarmaga, Efficient representation of minimally entangled typical thermal states in two dimensions via projected entangled pair states, Physical Review B 109, 045136 (2024), publisher: American Physical Society.
  * Kadow  _et al._ [2023] W. Kadow, F. Pollmann, and M. Knap, Isometric tensor network representations of two-dimensional thermal states, Physical Review B 107, 205106 (2023), publisher: American Physical Society.
  * Gogolin  _et al._ [2011] C. Gogolin, M. P. Müller, and J. Eisert, Absence of thermalization in nonintegrable systems, Phys. Rev. Lett. 106, 040401 (2011).
  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems, Reports on Progress in Physics 79, 056001 (2016).
  * Cramer [2012] M. Cramer, Thermalization under randomized local hamiltonians, New Journal of Physics 14, 053051 (2012).
  * Reichental  _et al._ [2018] I. Reichental, A. Klempner, Y. Kafri, and D. Podolsky, Thermalization in open quantum systems, Phys. Rev. B 97, 134301 (2018).
  * Shirai  _et al._ [2018] T. Shirai, T. Mori, and S. Miyashita, Floquet–Gibbs state in open quantum systems, Eur. Phys. J. Spec. Top. 227, 323 (2018).
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671 (1983), https://www.science.org/doi/pdf/10.1126/science.220.4598.671 .
  * Krzakala  _et al._ [2007] F. Krzakala, A. Montanari, F. Ricci-Tersenghi, G. Semerjian, and L. Zdeborova, Gibbs states and the set of solutions of random constraint satisfaction problems, Proceedings of the National Academy of Sciences 104, 10318 (2007), https://www.pnas.org/doi/pdf/10.1073/pnas.0703685104 .
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Optimal hamiltonian simulation by quantum signal processing, Phys. Rev. Lett. 118, 010501 (2017).
  * Jiang  _et al._ [2018] Z. Jiang, K. J. Sung, K. Kechedzhi, V. N. Smelyanskiy, and S. Boixo, Quantum algorithms to simulate many-body physics of correlated fermions, Physical Review Applied 9, 044036 (2018).
  * Wright  _et al._ [2021] L. Wright, F. Barratt, J. Dborin, G. H. Booth, and A. G. Green, Automatic post-selection by ancillae thermalization, Phys. Rev. Res. 3, 033151 (2021).
  * Turro  _et al._ [2022] F. Turro, A. Roggero, V. Amitrano, P. Luchi, K. A. Wendt, J. L. Dubois, S. Quaglioni, and F. Pederiva, Imaginary-time propagation on a quantum chip, Physical Review A 105, 022440 (2022).
  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Quantum computational chemistry, Rev. Mod. Phys. 92, 015003 (2020).
  * Kondappan  _et al._ [2023] M. Kondappan, M. Chaudhary, E. O. Ilo-Okeke, V. Ivannikov, and T. Byrnes, Imaginary-time evolution with quantum nondemolition measurements: Multiqubit interactions via measurement nonlinearities, Phys. Rev. A 107, 042616 (2023), publisher: American Physical Society.
  * Kosugi  _et al._ [2023] T. Kosugi, H. Nishi, and Y.-i. Matsushita, Exhaustive search for optimal molecular geometries using imaginary-time evolution on a quantum computer, npj Quantum Inf 9, 1 (2023), number: 1 Publisher: Nature Publishing Group.
  * Jouzdani  _et al._ [2022] P. Jouzdani, C. W. Johnson, E. R. Mucciolo, and I. Stetcu, Alternative approach to quantum imaginary time evolution, Phys. Rev. A 106, 062435 (2022).
  * Hejazi  _et al._ [2023] K. Hejazi, M. Motta, and G. K.-L. Chan, Adiabatic quantum imaginary time evolution (2023), arXiv:2308.03292 [quant-ph].
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nature Physics 16, 205 (2020).
  * Kosugi  _et al._ [2022] T. Kosugi, Y. Nishiya, H. Nishi, and Y.-i. Matsushita, Imaginary-time evolution using forward and backward real-time evolution with a single ancilla: First-quantized eigensolver algorithm for quantum chemistry, Phys. Rev. Res. 4, 033121 (2022).
  * G. McCaul  _et al._ [2021] G. McCaul, K. Jacobs, and D. I. Bondar, Fast computation of dissipative quantum systems with ensemble rank truncation, Phys. Rev. Research 3, 013017 (2021).
  * Le Bris and Rouchon [2013] C. Le Bris and P. Rouchon, Low-rank numerical approximations for high-dimensional Lindblad equations, Phys. Rev. A 87, 022125 (2013).
  * Finazzi  _et al._ [2015] S. Finazzi, A. Le Boité, F. Storme, A. Baksic, and C. Ciuti, Corner-Space Renormalization Method for Driven-Dissipative Two-Dimensional Correlated Systems, Phys. Rev. Lett. 115, 080604 (2015).
  * Donatella  _et al._ [2021] K. Donatella, Z. Denis, A. Le Boité, and C. Ciuti, Continuous-time dynamics and error scaling of noisy highly entangling quantum circuits, Phys. Rev. A 104, 062407 (2021).
  * Spencer [2014] J. Spencer,  _Asymptopia_ , Student mathematical library (American Mathematical Society, Providence, RI, 2014).
  * Leamer  _et al._ [2023] J. M. Leamer, W. Dawson, and D. I. Bondar, Positivity preserving density matrix minimization at finite temperatures via square root (2023), arXiv:2103.07078 .
  * Bondar  _et al._ [2016] D. I. Bondar, A. G. Campos, R. Cabrera, and H. A. Rabitz, Efficient computations of quantum canonical Gibbs state in phase space, Physical Review E 93, 10.1103/PhysRevE.93.063304 (2016).
  * Hirsch and Schrieffer [1983a] J. E. Hirsch and J. R. Schrieffer, Dynamic correlation functions in quantum systems: A Monte Carlo algorithm, Physical Review B 28, 5353 (1983a), publisher: American Physical Society.
  * Hirsch and Schrieffer [1983b] J. E. Hirsch and J. R. Schrieffer, Dynamic correlation functions in quantum systems: A monte carlo algorithm, Phys. Rev. B 28, 5353 (1983b).
  * Ceperley and Bernu [1988] D. M. Ceperley and B. Bernu, The calculation of excited state properties with quantum monte carlo, The Journal of Chemical Physics 89, 6316–6328 (1988).
  * Chowdhury and Somma [2016] A. N. Chowdhury and R. D. Somma, Quantum algorithms for Gibbs sampling and hitting-time estimation (2016), arXiv:1603.02940 [quant-ph].
  * McCaul and Bondar [2021a] G. McCaul and D. I. Bondar, How to win friends and influence functionals: deducing stochasticity from deterministic dynamics, The European Physical Journal Special Topics 230, 733 (2021a).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Repo. Prog. Phys. 75, 126001 (2012).
  * Spinney and Ford [2012] R. E. Spinney and I. J. Ford,  _Fluctuation relations: a pedagogical overview_ (2012) pp. 1–58, arXiv:1201.6381 .
  * Crooks [1999] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Phys. Rev. E 60, 2721 (1999).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Reports on Progress in Physics 29, 255 (1966).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, Journal of High Energy Physics 2014, 10.1007/jhep03(2014)067 (2014).
  * Touil and Deffner [2024] A. Touil and S. Deffner, Information scrambling – a quantum thermodynamic perspective (2024), arXiv:2401.05305 .
  * Hubbard [1959] J. Hubbard, Calculation of partition functions, Phys. Rev. Lett. 3, 77 (1959).
  * McCaul  _et al._ [2017] G. M. G. McCaul, C. D. Lorenz, and L. Kantorovich, Partition-free approach to open quantum systems in harmonic environments: An exact stochastic liouville equation, PRB 95, 125124 (2017).
  * McCaul  _et al._ [2018] G. M. G. McCaul, C. D. Lorenz, and L. Kantorovich, Driving spin-boson models from equilibrium using exact quantum dynamics, PRB 97, 224310 (2018).
  * McCaul and Bondar [2021b] G. McCaul and D. I. Bondar, How to win friends and influence functionals: deducing stochasticity from deterministic dynamics, EPJST 230, 733 (2021b).
  * Dziarmaga [2005] J. Dziarmaga, Dynamics of a Quantum Phase Transition: Exact Solution of the Quantum Ising Model, Physical Review Letters 95, 245701 (2005).
  * Sachdev [1999] S. Sachdev,  _Quantum phase transitions_ (Cambridge University Press, Cambridge ; New York, 1999).
  * Kopeć  _et al._ [1989] T. K. Kopeć, K. D. Usadel, and G. Büttner, Instabilities in the quantum Sherrington-Kirkpatrick Ising spin glass in transverse and longitudinal fields, Physical Review B 39, 12418 (1989), publisher: American Physical Society.
  * Laumann  _et al._ [2008] C. Laumann, A. Scardicchio, and S. L. Sondhi, Cavity method for quantum spin glasses on the Bethe lattice, Physical Review B 78, 134424 (2008), publisher: American Physical Society.
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum Computation by Adiabatic Evolution (2000), arXiv:quant-ph/0001106.
  * Rønnow  _et al._ [2014] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Defining and detecting quantum speedup, Science 345, 420 (2014), https://www.science.org/doi/pdf/10.1126/science.1252319 .
  * Shin  _et al._ [2014] S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani, How “Quantum” is the D-Wave Machine? (2014), arXiv:1401.7087 [quant-ph].
  * Boixo  _et al._ [2014] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nature Physics 10, 218 (2014), number: 3 Publisher: Nature Publishing Group.
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, QuSpin: a Python package for dynamics and exact diagonalisation of quantum many body systems. Part II: bosons, fermions and higher spins, SciPost Physics 7, 020 (2019).
  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * Anderson [1987] P. W. Anderson, The resonating valence bond state in La2CuO4 and superconductivity, Science 235, 1196 (1987).
  * Hofrichter  _et al._ [2016] C. Hofrichter, L. Riegger, F. Scazza, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Direct probing of the mott crossover in the $\mathrm{SU}(n)$ fermi-hubbard model, Phys. Rev. X 6, 021030 (2016).
  * Cleve  _et al._ [1998] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Quantum algorithms revisited, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339 (1998).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996).
  * Qiskit contributors [2021] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2021).
  * Cortes and Gray [2022] C. L. Cortes and S. K. Gray, Quantum krylov subspace algorithms for ground- and excited-state energy estimation, Phys. Rev. A 105, 022417 (2022).
  * Kyriienko [2020] O. Kyriienko, Quantum inverse iteration algorithm for programmable quantum simulators, npj Quantum Information 6, 10.1038/s41534-019-0239-7 (2020).
  * Scali  _et al._ [2023] S. Scali, C. Umeano, and O. Kyriienko, Quantum topological data analysis via the estimation of the density of states (2023), arXiv:2312.07115 .
  * Scali  _et al._ [2024] S. Scali, C. Umeano, and O. Kyriienko, The topology of data hides in quantum thermal states (2024), arXiv:2402.15633 .
  * Khaneja and Glaser [2001] N. Khaneja and S. J. Glaser, Cartan decomposition of su(2n) and control of spin systems, Chemical Physics 267, 11 (2001).
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating hamiltonian dynamics with a truncated taylor series, Phys. Rev. Lett. 114, 090502 (2015).
  * Chruściński [2006] D. Chruściński, Koopman’s approach to dissipation, Rep. Math. Phys. 57, 319 (2006).
  * McCaul  _et al._ [2019] G. McCaul, A. Pechen, and D. I. Bondar, Entropy nonconservation and boundary conditions for hamiltonian dynamical systems, Phys. Rev. E 99, 062121 (2019).
  * McCaul and Bondar [2022] G. McCaul and D. I. Bondar, Free to harmonic unitary transformations in quantum and Koopman dynamics*, Journal of Physics A: Mathematical and Theoretical 55, 434003 (2022).
  * Bondar  _et al._ [2013] D. I. Bondar, R. Cabrera, D. V. Zhdanov, and H. A. Rabitz, Wigner phase-space distribution as a wave function, Phys. Rev. A 88, 052108 (2013), arXiv:1202.3628 .
  * McCaul  _et al._ [2023] G. McCaul, D. V. Zhdanov, and D. I. Bondar, Wave operator representation of quantum and classical dynamics, Phys. Rev. A 108, 052208 (2023).
  * An  _et al._ [2023] D. An, J.-P. Liu, and L. Lin, Linear combination of hamiltonian simulation for nonunitary dynamics with optimal state preparation cost, Phys. Rev. Lett. 131, 150603 (2023).
  * Jin  _et al._ [2023] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations: Applications and detailed analysis, Phys. Rev. A 108, 032603 (2023).
  * Rossini  _et al._ [2023] M. Rossini, D. Maile, J. Ankerhold, and B. I. C. Donvil, Single-qubit error mitigation by simulating non-markovian dynamics, Phys. Rev. Lett. 131, 110603 (2023).
  * [101] DLMF, NIST Digital Library of Mathematical Functions, https://dlmf.nist.gov/, Release 1.1.10 of 2023-06-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.

##  Appendix A Ground State Acceleration via $\lambda$

Figure 8: Illustration of the convergence to the ground state in a transverse
field Ising model. In all cases, convergence is achieved significantly faster
than evolution with the Bloch equation. The bottom panel demonstrates how
$\tau_{c}$ \- the temperature at which $|\langle{H}\rangle-E_{0}|<10^{-4}$
($E_{0}$ being the ground state energy) - has the expected exponential
dependence on $\lambda$.

Here, we discuss how tuning $\lambda$ can markedly improve the rate of
convergence to the ground state with respect to $\tau$. This phenomenon can be
most easily illustrated with the case of a two level Hamiltonian. We denote
the energies of the two levels in terms of the average energy $\bar{E}$ and
difference $2\Delta$, i.e. $E_{0}=\bar{E}-\Delta$ and $E_{1}=\bar{E}+\Delta$.
If one were to evolve via the Bloch equation i.e.

$\frac{d\rho}{d\tau}=-\frac{1}{2}\\{\rho,H\\}$ (52)

one would obtain (after normalisation) the Gibbs state

$\frac{1}{Z}e^{-\tau\left(H-\lambda\right)}=\frac{1}{Z}\left(e^{\tau\Delta}|E_{0}\rangle\langle
E_{0}|+e^{-\tau\Delta}|E_{1}\rangle\langle E_{1}|\right).$ (53)

This implies that the ratio of probabilities for occupying the excited state
$p(E_{1})$ versus the ground state $p(E_{0})$ is
$\frac{p(E_{1})}{p(E_{0})}=e^{-2\tau\Delta}$ and that the shift parameter has
no effect. The value of $\tau$ to which one needs to evolve in order to obtain
a good approximation to the ground state is determined by $\Delta$, with the
ground state only being well-approximated when
$\tau\sim\mathcal{O}(\frac{1}{\Delta})$. For many complex systems, the low
level energy states may be nearly degenerate and thus require an impractically
large $\tau$ to obtain a good estimate for the ground state. If one uses any
$\lambda<E_{0}$, however, we find

$\displaystyle\frac{1}{Z}{\rm
e}^{-\tau\left(H+\lambda\right)^{2}}=\frac{1}{Z}$ $\displaystyle\bigg{(}{\rm
e}^{2\tau\Delta(\bar{E}+\lambda)}\outerproduct{E_{0}}{E_{0}}\bigg{.}$
$\displaystyle\bigg{.}+$ $\displaystyle{\rm
e}^{-2\tau\Delta(\bar{E}+\lambda)}\outerproduct{E_{1}}{E_{1}}\bigg{)}$ (54)
$\displaystyle\implies\frac{p(E_{1})}{p(E_{0})}={\rm
e}^{-4\beta\Delta(\bar{E}+\lambda)}.$ (55)

The rate of convergence to the ground state with respect to $\tau$ can
therefore be tuned by $\lambda$, such that the required final value is reduced
by a factor of $2\left(\bar{E}+\lambda\right)$. This phenomenon is illustrated
in Fig. 8, where decreasing $\lambda$ has the predicted exponential effect in
the imaginary time $\tau_{c}$ at which $\rho^{(\lambda)}$ reaches the ground
state, i.e. $|\langle H\rangle-E_{0}|<10^{-4}$.

##  Appendix B Additional Numerical Analyses

Figure 9: Results of spectral calculations performed for an eight site Ising
Hamiltonian. Dashed lines indicate the eigenenergies of $H_{TFIM}$. Solid
lines indicate the convergence of the spectral calculation to these
eigenenergies. For a), we plot the calculated energy as $\lambda$ against
propagation in imaginary time $\tau$ and find that our method converges to
successively lower eigenstates of $H_{TFIM}$ before reaching the ground state
$E_{0}$. For b), we plot the change in the steady state expectation value of
$H_{TFIM}$ as $\lambda$ is varied, illustrating that a complete course-grained
spectrum is obtainable using only expectation values estimated from a single
trajectory only.

Here, we apply our method to compute the spectrum of a transverse-field Ising
model and initialize the system in a single pure state
$|\psi_{0}\rangle=\sum_{j}a_{j}|E_{j}\rangle$. Results are plotted in Fig. 9.

In this case, we specifically initialize in the all-zero state given by
$|\psi_{0}\rangle=|0\rangle^{\otimes 8}$ is used. This implies that
$\rho^{(\lambda)}(\tau)$ has the form

$\rho^{(\lambda)}(\tau)=\sum_{ij}e^{-\frac{1}{4}\tau(E_{i}+E_{j}+2\lambda)^{2}}a_{i}a_{j}|E_{i}\rangle\langle
E_{j}|.$ (56)

Consequently, when $\lambda$ is chosen such that $\langle
H^{(\lambda)}\rangle(\tau\rightarrow\infty)$ converges to a band of closely
spaced levels, it tends to converge towards the eigenstate corresponding to
the lowest energy level in the band (i.e., rather than selecting out each
level within the band separately). This occurs because of additional
coherences in the energetic basis that are not present when
$\rho^{(\lambda)}(0)\propto{1}$, thus reinforcing the notion that it is
desirable to initialize the method accordingly in the maximally mixed state.

##  Appendix C Hadamard Test

Qubits are initialised in the state $\ket{0}$ before a Clifford gate is
applied, in order to randomise the initial state such that over statistical
sampling it mimics ${\rho}_{S}={1}$. A Hadamard gate is also applied to the
ancilla qubit, such that the overall state of the system at the first barrier
is described by:

$\displaystyle\ket{\Psi_{1}}=\frac{1}{\sqrt{2}}(\ket{\psi_{0}}\ket{0}+\ket{\psi_{0}}\ket{1}),$
(57)

where $\ket{\psi_{0}}$ denotes the randomised initial state of our system. We
then apply $2j$ copies of the unitary evolution operator ${\mathcal{U}}$ to
$\ket{\psi_{0}}$ to simulate $2j$ steps of the evolution. Each of these
operators is controlled by the ancilla qubit. Then we have

$\displaystyle\ket{\Psi_{2}}$
$\displaystyle=\frac{1}{\sqrt{2}}({\ket{\psi_{0}}\ket{0}}+{\mathcal{U}}^{2j}\ket{\psi_{0}}{\ket{1}}),$
(58)
$\displaystyle=\frac{1}{\sqrt{2}}({\ket{\psi_{0}}\ket{0}}+\ket{\psi_{2j}}{\ket{1}}).$
(59)

In general, one needs to decompose an observable into a series of Pauli
strings so that the expectation value can be evaluated efficiently. Applying
the Pauli string ${P}_{\alpha}$ controlled on the ancilla qubit yields

$\displaystyle\ket{\Psi_{3}}$
$\displaystyle=\frac{1}{\sqrt{2}}({\ket{\psi_{0}}\ket{0}}+{P}_{\alpha}\ket{\psi_{2j}}{\ket{1}}).$
(60)

Applying another $2j$ copies of ${\mathcal{U}}$ controlled on the ancilla then
yields

$\displaystyle\ket{\Psi_{4}}$
$\displaystyle=\frac{1}{\sqrt{2}}({\ket{\psi_{0}}\ket{0}}+\mathcal{U}^{2j}{P}_{\alpha}\ket{\psi_{2j}}{\ket{1}}).$
(61)

The final Hadamard acting on the ancilla mixes the superposition again to give

$\displaystyle\ket{\Psi_{5}}=$
$\displaystyle\frac{1}{2}\bigg{(}|0\rangle\big{[}|\psi_{0}\rangle+\mathcal{U}^{2j}P_{\alpha}|\psi_{2j}\rangle\big{]}$
(62)
$\displaystyle+|1\rangle\big{[}|\psi_{0}\rangle-\mathcal{U}^{2j}P_{\alpha}|\psi_{2j}\rangle\big{]}\bigg{)}$
(63)

Now the probability of measuring $\ket{0}$ on the ancilla is

$\displaystyle p_{0}$
$\displaystyle=\frac{1}{4}\big{(}\langle\psi_{0}|+\langle\psi_{2j}|P_{\alpha}(\mathcal{U}^{2j})^{\dagger}\big{)}\big{(}|\psi_{0}\rangle+\mathcal{U}^{2j}P_{\alpha}|\psi_{2j}\rangle\big{)},$
(64)
$\displaystyle=\frac{2+2\mathrm{Re}(\bra{\psi_{2j}}{P}_{k}\ket{\psi_{-2j}})}{4}.$
(65)

From this, it is possible to obtain the real part of the overlap via a
measurement of the ancilla, as

$\displaystyle\mathrm{Re}(\bra{\psi_{2j}}{P}_{k}\ket{\psi_{-2j}})=2p_{0}-1.$
(66)

Obtaining the imaginary part of this overlap can be achieved by modifying the
circuit by applying a $S$ gate to the ancilla between the first Hadamard gate
and $\mathcal{U}^{2j}$ gate.

##  Appendix D Estimator Method

Figure 10: Energy expectation values for a 3 site transverse field Ising model
with $J=1$, $h=2$ as the system propagates in imaginary time obtained. The
overlaps required for calculating Eqs. (32) and (33) are obtained using IBM’s
estimator primitive with the ibm qasm simulator to measure the expectation
values of ${S}_{\mathrm{Re}}^{j}$ and ${S}_{\mathrm{Im}}^{j}$ as in Eq. (69)

Another method of obtaining results for larger systems is to calculate the
necessary overlaps directly from expectations. Explicitly, the real and
imaginary parts of the overlap $\langle\psi_{2j}|O|\psi_{-2j}\rangle$ can be
expressed as expectations of the Hermitian operators $S_{\mathrm{Re}}^{j}$ and
$S_{\mathrm{Im}}^{j}$ i.e.

$\displaystyle S_{\mathrm{Re}}^{j}$
$\displaystyle=\frac{1}{2}\left({\mathcal{U}}^{2j}+({\mathcal{U}}^{2j})^{\dagger}\right),$
(67) $\displaystyle S_{\mathrm{Im}}^{j}$
$\displaystyle=\frac{1}{2i}\left({\mathcal{U}}^{2j}-({\mathcal{U}}^{2j})^{\dagger}\right),$
(68)

such that

$\displaystyle\mathrm{Re}/\mathrm{Im}\left(\langle\psi_{-2j}|O|\psi_{2j}\rangle\right)=\langle\psi_{0}|O{S}_{\mathrm{Re}/\mathrm{Im}}^{j}|\psi_{0}\rangle.$
(69)

This alternative approach reduces the number of qubits used in the simulation,
which reduces the errors in the calculation due to noisy hardware. However,
this approach does not scale well in comparison to the Hadamard test approach
discussed above. This is because measuring the expectation values of
${S}_{\mathrm{Re}/\mathrm{Im}}^{j}$ in the circuit requires decomposing
${\mathcal{U}}^{2j}$ and $({\mathcal{U}}^{2j})^{\dagger}$ into series of Pauli
strings for every step. With the Hadamard test approach, we only needed to
perform the Pauli string decomposition $O$ once and then we could use it each
time we evaluated the circuit. We demonstrate the utility of this approach for
a 1D transverse field Ising model with $J=1$, $h=2$, and $L=3$ by plotting the
energy expectations against imaginary time in Fig. 10. These results are
obtained using IBM’s qasm simulator with the estimator primitive backend
rather than a real machine. To mimic starting with the maximally mixed
entangled state, we repeat the calculations with 50 different random initial
states and averaging calculated overlaps at each step. These results
demonstrate that the alternative approach of measuring the expectations of
$S_{\mathrm{Re}}^{j}$ and $S_{\mathrm{Im}}^{j}$ can be used to calculate the
energy spectrum of the system for larger systems than the Hadamard test
approach for near term quantum computers.

##  Appendix E Error Bounds

One may obtain an error for ITQDE by considering the continuous limit of the
correspondence, namely the Hubbard-Stratonovich transformation as obtained in
Sec. III. Let us evaluate the integral via Gauss-Hermite quadrature (see Eqs.
(3.5.15), (3.5.19), and (3.5.28) of Ref. [101]),

$\displaystyle e^{-\tau H^{2}}$
$\displaystyle=\frac{1}{\sqrt{\pi}}\sum_{m}\int_{-\infty}^{\infty}e^{-x^{2}}e^{-2i\sqrt{\tau}xE_{m}}\ket{E_{m}}\bra{E_{m}}dx$
$\displaystyle=\frac{1}{\sqrt{\pi}}\sum_{m}\sum_{k=1}^{n}w_{k}e^{-2i\sqrt{\tau}x_{k}E_{m}}\ket{E_{m}}\bra{E_{m}}+{\mathcal{E}}$
$\displaystyle=\frac{1}{\sqrt{\pi}}\sum_{k=1}^{n}w_{k}e^{-2i\sqrt{\tau}x_{k}H}+{\mathcal{E}},$
(70)

where the coordinate grid $\\{x_{k}\\}$ and weights $\\{w_{k}\\}$ are
tabulated in Ref. [101]. Moreover, there exists reals $\\{\xi_{m}\\}$ such
that the error term is

$\displaystyle{\mathcal{E}}$
$\displaystyle=\frac{n!2^{n}}{(2n)!}(-i\sqrt{\tau})^{2n}\sum_{m}E_{m}^{2n}e^{-2i\sqrt{\tau}\xi_{m}E_{m}}\ket{E_{m}}\bra{E_{m}}$
$\displaystyle=\frac{n!2^{n}}{(2n)!}(-i\sqrt{\tau})^{2n}\sum_{m,m^{\prime}}E_{m^{\prime}}^{2n}e^{-2i\sqrt{\tau}\xi_{m}E_{m}}\ket{E_{m}}\bra{E_{m}}\delta_{m,m^{\prime}}$
$\displaystyle=\frac{n!2^{n}}{(2n)!}(-i\sqrt{\tau})^{2n}\sum_{m,m^{\prime}}E_{m^{\prime}}^{2n}\ket{E_{m^{\prime}}}\bra{E_{m^{\prime}}}e^{-2i\sqrt{\tau}\xi_{m}E_{m}}\ket{E_{m}}\bra{E_{m}}$
$\displaystyle=\frac{n!2^{n}}{(2n)!}(-i\sqrt{\tau})^{2n}H^{2n}\sum_{m}e^{-2i\sqrt{\tau}\xi_{m}E_{m}}\ket{E_{m}}\bra{E_{m}}.$

To write it more compactly, we can say that there is a unitary operator
${\mathcal{U}}$ such that the error ${\mathcal{E}}$ is

$\displaystyle{\mathcal{E}}=\frac{n!(2\tau)^{n}}{(2n)!}H^{2n}{\mathcal{U}}.$
(71)

The error is written in the polar operator form.

From Eq. (5.6.1) of Ref. [101], we get

$\displaystyle\sqrt{\frac{2\pi}{n+1}}\left(\frac{n+1}{e}\right)^{n+1}<n!<\sqrt{\frac{2\pi}{n+1}}\left(\frac{n+1}{e}\right)^{n+1}e^{\frac{1}{12n+12}}.$
(72)

Hence, for a unitary invariant matrix norm $\|\cdot\|$, we obtain the bound

$\displaystyle\|{\mathcal{E}}\|<\frac{(n+1)^{n+1/2}}{(2n+1)^{2n+1/2}}e^{\frac{1}{12n+12}}(2\tau
e)^{n}\|H^{2n}\|,$ (73)

which leads to

$\displaystyle\|{\mathcal{E}}\|=\mathcal{O}\left(\left(\frac{\tau
e}{2n}\right)^{n}\|H^{2n}\|\right),\quad n\to\infty.$ (74)