NuQS Collaboration

# Quantum algorithms for transport coefficients in gauge theories

Thomas D. Cohen<EMAIL_ADDRESS>Department of Physics, University of Maryland,
College Park, MD 20742, USA  Henry Lamm<EMAIL_ADDRESS>Fermi National
Accelerator Laboratory, Batavia, Illinois, 60510, USA  Scott Lawrence
<EMAIL_ADDRESS>Department of Physics, University of Colorado,
Boulder, CO 80309, USA  Yukari Yamauchi<EMAIL_ADDRESS>Department of
Physics, University of Maryland, College Park, MD 20742, USA

###### Abstract

In the future, _ab initio_ quantum simulations of heavy ion collisions may
become possible with large-scale fault-tolerant quantum computers. We propose
a quantum algorithm for studying these collisions by looking at a class of
observables requiring dramatically smaller volumes: transport coefficients.
These form nonperturbative inputs into theoretical models of heavy ions; thus,
their calculation reduces theoretical uncertainties without the need for a
full-scale simulation of the collision. We derive the necessary lattice
operators in the Hamiltonian formulation and describe how to obtain them on
quantum computers. Additionally, we discuss ways to efficiently prepare the
relevant thermal state of a gauge theory.

††preprint: FERMILAB-PUB-21-091-T

##  I Introduction

The ultimate promise of quantum computers is that physical properties
effectively uncalculable by classical algorithms can be obtained [1]. Within
high-energy physics, it is believed that quantum algorithms will shed light
upon topics where nonperturbative nonequilibrium dynamics play a role [2]. The
limitations of current classical methods are particularly acute in heavy-ion
collisions.

The standard picture of heavy-ion collisions is that two nuclei collide at
high energy, and through multiple scatterings form an expanding fireball of
quark-gluon plasma (QGP). As the plasma cools, the quarks and gluons
rehadronize into mesons and baryons that are then measured in the detector.
The current theoretical framework divides the collision into overlapping steps
which can each be approximated by a semiclassical model [3]. In particular,
the incredible effectiveness of hydrodynamics to model the QGP is notable [4,
5]. While these calculations can account for some of the nonperturbative and
nonequilibrium dynamics of QCD at the LHC and RHIC [6, 7, 8, 9, 10, 11, 12,
13, 14, 15], the large number of free parameters and uncontrolled systematics
used in theoretical modeling are unpalatable [16, 17, 18, 19]. Of interest to
this work are the transport coefficients: diffusivity [20, 21, 22],
conductivity [21, 22, 23], and viscosity [24, 25, 26, 27, 28, 29, 30, 30, 19].

In principle, a complete nonperturbative calculation would be possible from
lattice field theory (LFT). Unfortunately, all known formulations of LFT
suitable for classical computations use the Euclidean metric in order to allow
stochastic methods to sample the path integral. While this allows for a broad
range of observables to be computed, its capabilities at finite-density or for
Minkowksi observables are practically limited by _sign problems_ [31] or the
ill-posedness of analytic continuation [32]. Despite these issues, efforts
have been undertaken to extract nonperturbative inputs [29] — e.g. parton
distrubtion functions (PDFs) [33, 34, 35, 36, 37, 38], phase diagrams [39, 40,
41, 42], transport coefficients [43, 44, 45, 46, 47, 48] — for the models from
lattice QCD, albeit with large uncertainties.

Quantum computers provide a natural facility for studying real-time dynamics
of quantum systems. With sufficiently large quantum resources, one should be
able to access arbitrary Minkowski matrix elements in LFT; however, such
calculations are forbiddingly expensive. For example, consider the task of
simulating heavy ion collisions [49]. This corresponds to an extension of the
scattering calculations considered in [50, 51, 52]. To prepare well separated
wave-packets for the ions, the spatial lattice must be many times the diameter
of a heavy nucleus, which is $\sim 10\;\mathrm{fm}$. To simulate the internal
dynamics of each nucleus, a resolution smaller than the nuclear radius is
required, thus the lattice spacing should be $\sim 0.1\;\mathrm{fm}$.
Together, these two scales imply $\gtrsim 10^{6}$ lattice sites are required,
and so are at least that many qubits. To estimate the circuit depth, one
should evolve the calculation for sufficient time that the two nuclei can
collide and the final states become sufficiently separated, again $\sim
10\;\mathrm{fm}$ [53]. In order to keep the trotterization errors small, one
must use a small time step, again $\sim 0.1\;\mathrm{fm}$. This suggests $\sim
10^{2}$ time steps, ignoring the circuits required for state preparation and
measurement. These estimates are certainly beyond near-term prospects. From
this example, we can see that the need to represent multiple scales accurately
drives the large resource requirements.

The authors of [54] pointed out that by decomposing scattering simulations
into a nonpertubative input (in that case the hadronic tensor) convolved with
perturbative expressions, one reduces the resource requirements. In that
example, instead of preparing two well-separated protons, one prepares a
single proton in a $L\sim 1\;\mathrm{fm}$ box. This reduces qubit costs by a
factor of $\sim 10^{3}$.

Heavy-ion collisions allow for a similar decomposition. In particular, the
transport coefficients of the QGP should require $L\sim 1\;\mathrm{fm}$. Since
the transport coefficients represent the hydrodynamics of the theory, they
should require reduced resolution compared to partonic observables. Using
thermodynamic observables of QCD as a guide [55], one might anticipate a
required lattice spacing (for studying $\sim 200\;\mathrm{MeV}$ matter) of
$a\sim 0.1\;\mathrm{fm}$. Further, theoretical models [56, 57] and
experimental data [58, 59, 60] suggest that thermalization happens rapidly –
on the scale of $~{}1\;\mathrm{fm}$ – which would reduce the circuit depth.
One must further emphasize that the current theoretical uncertainties on the
transport coefficients of the QGP plasma are $O(1)$. Together, these arguments
strongly suggest that the transport coefficients of gauge theories represent
serious targets for practical quantum advantage in particle physics [2].

There is another point that weights heavily in favor of focusing upon
transport coefficients. For many problems, it is not _a priori_ obvious that
factoring the physical system into separate regimes is valid. In the case of
PDFs, for instance, factorization theorems have been proven only in specific
kinematic regimes [61]. It is well-known that using PDFs outside of these
regimes leads to issues [62, 63]. In heavy-ion collisions, the validity of the
hydrodynamic approximation has been studied at length (see [4] for a review)
and the transport coefficients can be nonperturbatively defined on the lattice
[44].

Past work has extensively discussed real-time evolution in a gauge theory via
digital quantum simulation. Here, we detail novel aspects of quantum
simulation in the hydrodynamic regime of gauge theories. In particular, we
construct implementations of the stress-energy tensor in the Hamiltonian
formulation, allowing hydrodynamic correlators to be measured.

The most involved part of a quantum simulation is typically the preparation of
the initial state [64, 65, 66, 67, 68, 69, 70, 71]. For transport
coefficients, the desired initial state is in thermal equilibrium. When
studying Yang-Mills or QCD, the most interesting temperatures are $\sim
200\;\mathrm{MeV}$. The fact that we are interested in relatively high-
temperature states is a key advantage to studying QGP transport over lower
temperature processes like scattering: thermal states are markedly easier to
prepare. In addition to descriptions of appropriate lattice measurements, this
paper details several viable methods for thermal state preparation.

This paper is organized as follows. Sec. II is devoted to describing different
methods of extracting transport coefficients from a quantum simulation — these
approaches are largely independent of the simulation scheme and the
construction of lattice observables. Since hydrodynamic transport coefficients
are derived from correlators of the energy-momentum tensor $T_{\mu\nu}$ (EMT),
in Sec. III we derive lattice operators for the energy-momentum tensor within
the Hamiltonian formalism. Methods for preparing quantum thermal states are
elucidated in Sec. IV. We conclude in Sec. V with some general points and
discuss where future theoretical work is required.

##  II Estimating Transport

Dissipative processes cannot be seen in perfect thermal equilibrium, and so
are not accessible to Euclidean lattice calculations. Near equilibrium (and in
the limit of small gradients), these processes are characterized by a small
number of low-energy constants, termed transport coefficients. The rate of
diffusion of a conserved quantity $\varphi$, for instance, is governed by
_Fick’s law_ :

$\vec{J}=-D\frac{d\varphi}{d\vec{x}}$ (1)

where $J$ is the diffusive flux, and $D$, the diffusion constant, is the low-
energy constant of interest. In this section, we will describe how to
determine such transport coefficients in a quantum simulation.

Many other transport coefficients can be defined, typically related to the
behavior of locally conserved quantities. In the case of hydrodynamics, the
most interesting are the shear and bulk viscosities $\eta$ and $\zeta$, which
are defined by their role in the Navier-Stokes equations:

$\rho\left(\frac{\partial u_{i}}{\partial
t}+u_{j}\nabla_{j}u_{i}\right)+\nabla
p=\eta\partial_{j}\left(\partial_{i}u_{j}+\partial_{j}u_{i}\right)+\left(\zeta-\frac{2}{3}\eta\right)\partial_{i}\partial_{j}u_{j}\text{,}$
(2)

where the mass density $\rho$ and velocity $u$ are functions of both space and
time, and the pressure $p$ is, in equilibrium, determined from the energy
density by the equation of state. Strictly speaking, in the case of field
theory we should use the relativistic Navier-Stokes equations [4], but up to a
substitution of the energy density $\epsilon$ for the mass density $\rho$, the
nonrelativistic Eq. (2) yields the same results with easier intuition. In
actual heavy ion simulations, it is necessary to include higher-order terms in
the gradient expansion, introducing new transport coefficients [72, 73, 74].
We do not discuss the determination of these second-order coefficients here.

For these transport coefficients, more or less the same set of methods are
available for their computation. For concreteness, in this section we will
focus on the shear viscosity. The methods described below generalize easily to
bulk viscosity, diffusion, and other transport properties.

Several well-established methods exist for computing the shear viscosity
$\eta$ of fluids in molecular dynamics simulations [75]. Three are worth
summarizing here; of these, two have natural analogues in quantum simulation
of gauge theories.

In the periodic perturbation (PP) method [75, 76], one imposes a shearing
force $\mathcal{F}$ on the system. The shear viscosity acts to resist this
force, so that in equilibrium, a small shear wave is created. In the limit of
small $\mathcal{F}$ and small wavenumbers $k$, the equilibrium “displacement”
is

$\mathcal{V}=\mathcal{F}\frac{\rho}{\eta k^{2}}\text{.}$ (3)

This method works well in molecular dynamics simulations. Unfortunately, it
intrinsically involves the imposition of a non-conservative force. Such a
force, by definition, cannot be imposed by the addition of any term to the
Hamiltonian; equivalently, any simulation of this system must be nonunitary.
This renders it inapplicable for calculations of quantum systems.

The transverse current autocorrelation function (TCAF) method [75, 76]
proceeds from the observation that a shear wave, once created, decays
exponentially with a decay constant proportional to $\eta$. To see this,
consider the linearization of Eq. (2) about the equilibrium solution
$\rho=\rho_{0}$, $u_{i}=0$. The mode $\mathcal{U}\cos(kx)$ obeys

$\rho_{0}\frac{\partial\mathcal{U}}{\partial
t}=-\eta\left(\vec{k}\cdot\vec{k}\mathcal{U}+\vec{k}(\vec{k}\cdot\mathcal{U})\right)-\left(\zeta-\frac{2}{3}\eta\right)\vec{k}(\vec{k}\cdot\mathcal{U})\text{.}$
(4)

Thus a shear wave $u_{1}(\vec{x})\propto\cos(kx_{2})$ is seen to decay as
$e^{-\frac{\eta}{\rho k^{2}}t}$. Fitting this exponential decay gives the
shear viscosity. Conveniently, the fact that the shear viscosity is encoded in
the decay constant means that the normalization of the $T_{\mu\nu}$ need not
be correct. (In fact, any operator that couples to the appropriate shear mode
is likely to show the correct exponential decay.)

In principle, then, we imagine creating a small shear wave in an equilibrated
fluid and watching it decay. This requires carefully taking the limit of a
small perturbative force. Instead, we can allow “Maxwell’s angels”, _i.e._
thermodynamic and quantum fluctuations, to create the shear wave. The shear
viscosity is now provided by fitting the TCAF at long111We assume the quantum
fluid is described at long wavelengths and times by pure Navier-Stokes, and
therefore neglect the phenomenon of long-time tails [77, 78]. For realistic
fluids, the exponential decay will eventually be cut off by a power-law. This
power-law decay can be extracted just as $\eta$ can. times:

$\displaystyle C(t)=$ $\displaystyle\int dx\,dy\,\sin(ky)\sin(kx)\langle
T_{01}(x,t)T_{01}(y,0)\rangle$ $\displaystyle\propto e^{-\frac{\eta}{\rho
k^{2}}t}\text{.}$ (5)

In the context of a classical molecular dynamics simulation, the calculation
of $C(t)$ is straightforward: $\int dx\,\sin(kx)T_{01}(x)$ is measured at
various times, and the autocorrelations appear in the resulting stochastic
time-series. At first blush, the appropriate procedure on a quantum computer
is rather different. In particular, because measurements cannot be performed
without altering the wavefunction, one might expect to be required to perform
one set of measurements for each time at which one wants to measure $C(t)$.
This is an expensive prospect.

Happily, we are explicitly interested in the behavior of _hydrodynamic_
fluctuations, whose apparent behavior is classical. The measurement of the
amplitude of the shear wave contains very little information relative to the
whole wavefunction. This measurement therefore constitutes a sort of “weak
measurement,” which can be made without unduly disturbing the quantum
state222The notion of a weak masurement — see e.g. [79] — is here simply used
to denote a measurement from which very little information is obtained about
the state.. Therefore, we can measure the amplitude of the shear wave at one
time, continue evolving, and subsequently measure the amplitude at many later
times.

Similarly, the energy density may be measured in the simulated system prior to
any time evolution is done, without significantly disturbing the wavefunction.
In fact, as long as the simulated system can be kept coherent, the TCAF can be
determined to arbitrary precision with one long time evolution — in the large-
volume limit. Away from the large-volume limit, the number of measurements
that can be made can be heuristically expected to be linear in the volume of
the system.

The third method follows from the Green-Kubo relation [80] connecting shear
viscosity to the $\omega,k\rightarrow 0$ limit of the two-point correlator of
$T_{12}$:

$\eta=\frac{V}{T}\int_{0}^{\infty}\langle T_{12}(t)T_{12}(0)\rangle\text{.}$
(6)

This method is reported to converge slowly in molecular dynamics simulations
[75], but has nevertheless proven useful in practice, _e.g._ [81]. Note that
this method, unlike the TCAF method, does require the operator for the stress
tensor to be correct, including the normalization.

These methods for computing $\eta$ have natural analogs for most other
transport coefficients. However, in the case of bulk viscosity, which is the
next-most relevant hydrodynamic property, there is no simple analog of the
TCAF method. Instead, molecular dynamics calculations of the $\zeta$ typically
proceed from the Green-Kubo relation, which is in terms of the diagonal
components of the $T_{\mu\nu}$ [82]:

$\zeta=\frac{V}{9T}\int_{0}^{\infty}dt\;\langle
T_{ii}(t)T_{jj}(0)\rangle\text{.}$ (7)

##  III Energy-Momentum Tensor

Transport coefficients in quantum field theory are defined in terms of the
$n$-point correlators of $T_{\mu\nu}$. Thus, in order to extract them from
lattice calculations, one must define a lattice $T_{\mu\nu}$. Here, we derive
such operators in the lattice Hamiltonian formalism natural for use in quantum
simulations. In Sec. III.1, we summarize the discretization of the EMT on a
spacetime lattice (as distinct from the spatial-only Hamiltonian lattice). In
Sec. III.2 we review the derivation of Kogut-Susskind Hamiltonian from the
transfer matrix [83] which provides a framework for deriving the Hamiltonian
formulation of the EMT. In the rest of the section, we derive the necessary
components of $T_{\mu\nu}$ for quantum simulation: diagonal components in
III.3, spatial components $T_{ij}$ in III.4, and time-like components $T_{0i}$
in III.5.

We provide both naive operators (which have $O(a)$ corrections) and tree-level
improved operators analogous to the ‘clover’ of the spacetime lattice (these
have $O(a^{2})$ corrections). The operators derived are summarized in Table 1.
They are constructed from link operators $\hat{U}$ and their conjugate momenta
$\hat{\pi}$, defined in Sec. III.2 below. The plaquette $\hat{P}$ and clover
$\hat{C}$ are defined in Sec. III.1.

Table 1: Gauge-invariant lattice operators in the Hamiltonian formalism in $3+1d$ dimensions: naive operators with $O(a)$ errors and improved operators with errors that are $O(a^{2})$. Components of the energy-momentum tensor $T_{\mu\nu}$ are constructed as linear combinations of these operators according to Eq. (8). The plaquette $\hat{P}$ and clover $\hat{C}$ are defined in Eq. (10) and Eq. (15), respectively. Spatial indices are $i\neq j\neq k$. Operator | $O(a)$ | $O(a^{2})$
---|---|---
$\mbox{Tr}F_{0i}F_{0i}(n)$ | $\frac{g_{s}^{2}}{a^{4}}\mbox{Tr}\left[\pi_{n,i}^{2}\right]$ | $\sum_{{x=0,1}}\frac{g_{s}^{2}}{2a^{4}}\mbox{Tr}\left[\pi_{n-x\hat{i},i}^{2}\right]$
$\mbox{Tr}F_{0i}F_{0j}(n)$ | $\frac{g_{s}^{2}}{a^{4}}\mbox{Tr}\left[\pi_{n,i}\pi_{n,j}\right]$ |  $\frac{g_{s}^{2}}{4a^{4}}\left(\mbox{Tr}\left[\hat{\pi}_{n,i}\hat{\pi}_{n,j}\right]+\mbox{Tr}\left[\hat{\pi}_{n,i}\hat{U}_{n-\hat{j},j}^{\dagger}\hat{\pi}_{n-\hat{j},j}\hat{U}_{n-\hat{j},j}\right]+\mbox{Tr}\left[\hat{U}_{n-\hat{i},i}^{\dagger}\hat{\pi}_{n-\hat{i}}\hat{U}_{n-\hat{i},i}\hat{\pi}_{n,j}\right]\right.$ $\left.+\mbox{Tr}\left[\hat{U}_{n-\hat{i},i}^{\dagger}\hat{\pi}_{n-\hat{i},i}\hat{U}_{n-\hat{i},i}\hat{U}_{n-\hat{j},j}^{\dagger}\hat{\pi}_{n-\hat{j},j}\hat{U}_{n-\hat{j},j}\right]\right)$
$\mbox{Tr}F_{0j}F_{ij}(n)$ | $-\frac{1}{a^{4}}\mbox{Tr}\left[\hat{\pi}_{n,j}\mbox{Im}\hat{P}_{ij}(n)\right]$ | $-\frac{1}{2a^{4}}\left(\mbox{Tr}\left[\hat{\pi}_{n,j}\mbox{Im}\hat{C}_{ij}(n)\right]+\mbox{Tr}\left[\hat{U}_{n-\hat{j},j}^{\dagger}\hat{\pi}_{n-\hat{j},j}\hat{U}_{n-\hat{j},j}\mbox{Im}\hat{C}_{ij}(n)\right]\right)$
$\mbox{Tr}F_{ij}F_{ij}(n)$ | $\frac{2}{g_{s}^{2}a^{4}}\mbox{Re}\mbox{Tr}\left[1-\hat{P}_{ij}(n)\right]$ | $\sum_{x=0,1}\sum_{y=0,1}\frac{1}{2g_{s}^{2}a^{4}}\mbox{Re}\mbox{Tr}\left[1-\hat{P}_{ij}(n-x\hat{i}-y\hat{j})\right]$
$\mbox{Tr}F_{ij}F_{kj}(n)$ | $\mbox{Tr}[\hat{F}_{ij}^{N}(n)\hat{F}_{kj}^{N}(n)]$ | $\mbox{Tr}[\hat{F}_{ij}^{C}(n)\hat{F}_{kj}^{C}(n)]$

###  III.1 Energy-momentum tensor on a lattice

The energy-momentum tensor is the Noether’s current of spacetime translational
symmetry. In LFT, this symmetry is explicitly broken to a discrete subgroup
and thus naive lattice currents may not be conserved. We expect to restore the
symmetry (i.e. the Ward identity) in the continuum limit, albeit
renormalization is required. Another tactic instead constructs combinations of
the low-dimension operators that mix under the discrete rotations and
translations [84, 85, 86] such that a lattice Ward identity is satisfied. A
final method extracts a UV-finite version of the EMT from gradient flow [87,
88]. In this work, we study two lattice EMTs: the naive EMT accurate to $O(a)$
and the clover EMT with $O(a^{2})$ errors. There exists an $O(a^{4})$ EMT
[89], but it is left for future work, as it likely first requires developing
an appropriately improved Hamiltonian [90, 91, 92].

In the continuum, the EMT of a gauge theory is

$T_{\mu\nu}=\frac{1}{4}g_{\mu\nu}\mbox{Tr}\left[F_{\alpha\beta}F^{\alpha\beta}\right]-\mbox{Tr}\left[F_{\mu\alpha}F_{\nu}^{\alpha}\right].$
(8)

where the mostly minus convention is used. Following the notation of [83], we
normalize the group generators $\lambda^{a}$ to

$\mbox{Tr}[\lambda^{a}\lambda^{b}]=\delta_{ab}$ (9)

and functions such as $F_{\mu\nu}$ are defined with this normalization. Before
we discretize Eq. (8), we must establish some lattice notation. Links are
denoted as $U_{n,\mu}$ where $n$ is the site the link starts at, and $\mu$ is
the direction of the link. The fundamental gauge-invariant object, the
plaquette, is defined as the product of four links around a closed loop as in
Fig. 1:

$P_{\mu\nu}(n)=U_{n,\mu}U_{n+\hat{\mu},\nu}U_{n+\hat{\nu},\mu}^{\dagger}U_{n,\nu}^{\dagger}.$
(10)

Figure 1: Schematic view of the plaquette, the clover, and the half-clover.

On a spacetime lattice, the terms in the EMT are polynomials of plaquettes.
Note that, as terms in Eq. (8) contain products of $F_{\mu\nu}$ in different
directions, they are most sensibly computed either at lattice sites or at the
center of spacetime volumes. For the sake of simplicity, we will focus on
discretizing the EMT on lattice sites.

The first term in Eq. (8) is the Lagrangian. On the lattice, it can be written

$\mbox{Tr}\left[F_{\mu\nu}(n)^{2}\right]=\frac{2}{g_{s}^{2}a_{\mu}^{2}a_{\nu}^{2}}\mbox{Re}\mbox{Tr}[1-P_{\mu\nu}(n)]+O(a)\text{,}$
(11)

where $g_{s}$ is the coupling constant. Here $a_{\mu}$ are the lattice
spacings in the $\mu$ direction. We denote the temporal spacing $a_{0}$, and
assume all other spacings to be the same and denote them as $a$. When working
in the Hamiltonian formalism, we take $a_{0}\rightarrow 0$ before performing
continuum extrapolations. If instead we considered
$\mbox{Tr}\left[F_{\mu\nu}(n+1/2[\hat{\mu}+\hat{\nu}])^{2}\right]$ located at
the center of the plaquette, we would find Eq. (11) is accurate to $O(a^{2})$.

To improve the discretization on a site up to $O(a^{2})$, we can simply
average over the RHS of Eq. (11) for four plaquettes around the site $n$ in
$\mu\nu$ plane:

$\displaystyle\mbox{Tr}\left[F_{\mu\nu}(n)^{2}\right]$ $\displaystyle=$
$\displaystyle\sum_{x={0,1}}\sum_{y={0,1}}\frac{1}{2g_{s}^{2}a_{\mu}^{2}a_{\nu}^{2}}$
$\displaystyle\mbox{Re}\mbox{Tr}[1$ $\displaystyle-$ $\displaystyle
P_{\mu\nu}(n-x\hat{\mu}-y\hat{\nu})]+O(a^{2})\text{.}$ (12)

In cases of $F_{0i}$, the $O(a_{0})$ error induced in the time direction is
acceptable as we take the limit $a_{0}\ll a$ in the Hamiltonian formalism.
This implies that we don’t need to average over four plaquettes in both
direction $\hat{t}$ and $\hat{i}$. Instead, we need to average over only two
plaquettes in $\hat{i}$ direction around the site:

$\displaystyle\mbox{Tr}\left[F_{oi}(n)^{2}\right]=$
$\displaystyle\sum_{x={0,1}}\frac{1}{g_{s}^{2}a_{\mu}^{2}a_{\nu}^{2}}$ ReTr
$\displaystyle[1-P_{\mu\nu}(n-x\hat{i})]+O(a^{2},a_{0})\text{.}$ (13)

The second term of Eq. (8), $F_{\mu\alpha}F_{\nu\alpha}$, requires us to
construct a discretization for $F_{\mu\nu}$ itself. The naive discretization
for the field strength tensor is

$F_{\mu\nu}^{N}(n)=-\frac{i}{2g_{s}a_{\mu}a_{\nu}}\left(P_{\mu\nu}(n)-P_{\mu\nu}^{\dagger}(n)\right)+O(a)\text{.}$
(14)

Note again that $F_{\mu\nu}^{N}(n)$ is evaluated on a lattice site. The RHS of
Eq. (14) approximates the value at the center of the plaquette up to
$O(a^{2})$. To improve the on-site discretization scheme up to $O(a^{2})$, one
can use so-called ‘clover’ operators as shown in Fig. 1:

$C_{\mu\nu}(n)=\frac{1}{4}\left[P_{\mu,\nu}+P_{\nu,-\mu}+P_{-\mu,-\nu}+P_{-\nu,\mu}\right](n)\text{.}$
(15)

From these clover operators we construct the improved clover discretization of
the field strength:

$F_{\mu\nu}^{C}(n)=-\frac{i}{2g_{s}a_{\mu}a_{\nu}}\left(C_{\mu\nu}(n)-C_{\mu\nu}^{\dagger}(n)\right)+O(a^{2})\text{.}$
(16)

Here we have used $P_{-\mu,\nu}(n)$, for instance, to denote the plaquette in
which the first link begins at site $n$ and ends at site $n-\hat{\mu}$. In an
Abelian theory, we of course have
$P_{-\mu,\nu}(n)=P_{\mu,\nu}(n-\mu)^{\dagger}$, as the starting site does not
matter. This is not the case in general.

The clover improvement ensures that the leading discretization errors are
$O(a^{2},a_{0}^{2})$ regardless of orientation. However, as mentioned, working
in the Hamiltonian formalism implies $a_{0}\ll a$. As a result, $O(a_{0})$
corrections are acceptable while $O(a)$ are not. This, combined with the
difficulty of implementing time-nonlocal operators, motivates the “half-
clover” operator — analogous to Eq. (III.1) — averaged over two plaquettes as
in Fig. 1:

$B_{i0}(n)=\frac{1}{2}\left(P_{i0}(n)+P_{0(-i)}(n)\right)\text{.}$ (17)

This is enough to implement $T_{\mu\nu}$ correctly up to $O(a^{2},a_{0})$:

$F^{B}_{i0}(n)=\frac{-i}{2g_{s}a_{0}a}\left(B_{i0}(n)-B_{i0}^{\dagger}(n)\right).$
(18)

###  III.2 Transfer matrix

With a choice of discretization of $F_{\mu\nu}$ and $F_{\mu\nu}^{2}$, one can
compute a lattice $T_{\mu\nu}$ in the action formulation. In the Hamiltonian
formulation, $T_{\mu\nu}$ (which are functions of the space of field
configurations) must be replaced by operators $\hat{T}_{\mu\nu}$ which act on
the Hilbert space of gauge links on a spatial lattice. A few options exist for
performing this transformation. One is to use the Legendre transform; however,
this is complicated by gauge invariance. The transfer matrix formalism has the
advantage of being manifestly gauge-invariant.

Given a path integral, we can construct a Hilbert space and a transfer matrix
such that the transfer matrix completely characterizes the path integral.
Taking the logarithm of the transfer matrix (in our case, in the
$a_{0}\rightarrow 0$ limit) yields a Hamiltonian usable for quantum
simulations.

Now, consider a perturbation to that action by a term proportional to
$\mathcal{O}$. The free energy, differentiated with respect to the
perturbation, yields the $\langle\mathcal{O}\rangle$. Connecting the perturbed
system to one in the Hamiltonian formalism via the transfer matrix, we obtain
a perturbed Hamiltonian, revealing what operator $\hat{\mathcal{O}}$ to use on
a quantum computer.

The Hamiltonian which gives the same dynamics as the Wilson action is derived
in [83]; this is the unperturbed case of the above procedure. In this section,
we first summarise the derivation of the Kogut-Susskind Hamiltonian $H_{KS}$
[93] via the transfer matrix, and then discuss the perturbed case which yields
specific operators of interest. Our starting point is the Minkowski path
integral corresponding to the Wilson action333Strictly speaking, the Wilson
action is in Euclidean space, with the sign of the second term of Eq. (19)
flipped. This section can be formulated in metric with no difference provided
the Hamiltonian limit is taken [94]; we have chosen Minkowski to preserve a
more straightforward correspondence with the quantum simulation.:

$\displaystyle S_{W}$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{N_{t}}K(t)+V(t)$ (19) $\displaystyle K(t)$
$\displaystyle=$
$\displaystyle\sum_{n}\sum_{i}\frac{a}{g_{s}^{2}a_{0}}\mbox{Re}\mbox{Tr}\left[1-P_{0i}(n,t)\right]$
(20) $\displaystyle V(t)$ $\displaystyle=$
$\displaystyle\sum_{n}\sum_{i<j}\frac{a_{0}}{g_{s}^{2}a}\mbox{Re}\mbox{Tr}\left[P_{ij}(n,t)\right]$
(21)

Here $n$ denotes a spatial site and $i,j$ denote spatial directions. The
corresponding path integral is of course the integral over all field
configurations of $e^{iS}$. Given a Hamiltonian $H$, we could also construct a
path integral by splitting the time evolution operator $e^{-iHt}$ into a
product of many nearly-identity pieces $\hat{T}=e^{-iHa_{0}}$ (this is the
transfer matrix) and inserting a complete set of states:

$Z=\int\mathcal{D}U\langle
U_{t}|\hat{P}\hat{T}|U_{t-a_{0}}\rangle\cdots\langle
U_{a_{0}}|\hat{P}\hat{T}|U_{0}\rangle=\int\mathcal{D}U\;e^{iS}\text{.}$ (22)

Here $|U\rangle$ is a basis state in “position basis”: an element of the gauge
group is specified at each link. The gauge projection operator $P$ has been
inserted between every pair of states in order to obtain a time-translation-
invariant action.

As a result, we see that for the Hamiltonian $H$ to give the same dynamics as
the action of Eq. (19), the matrix elements of the transfer matrix
$\hat{T}=e^{-ia_{0}\hat{H}}$ should be given by

$\langle
U^{\prime}|\hat{T}|U\rangle=e^{i(K(t)+V(t))}=e^{i\sum_{n}\left(\frac{a}{g_{s}^{2}a_{0}}\sum_{i}\mbox{Re}\mbox{Tr}\left[1-U_{n,i}U^{{}^{\prime}\dagger}_{n,i}\right]+\frac{a_{0}}{g_{s}^{2}a}\sum_{i,j}\mbox{Re}\mbox{Tr}\left[U_{n,i}U_{n+\hat{i},j}U_{n+\hat{j},i}^{\dagger}U_{n,j}^{\dagger}\right]\right)}\text{.}$
(23)

We are working in temporal gauge: $U_{n,0}=1$ for all $n$. For our derivations
this is an irrelevant technicality, but see [83, 95] for detailed expositions
of the relationship between timelike links and the gauge projection operator.

Eq. (23) suffices to define the transfer matrix, but we would like to express
it in terms of more natural objects in the Hamiltonian formulation: in
particular, the $\hat{U}$ that are diagonal in the position basis, and their
conjugate operators $\hat{\pi}$ defined below. As spatial plaquettes are
already written only with $U$s, they remain in the same form and all arguments
$U$ become operators $\hat{U}$. On the other hand $K$ results in operators not
diagonal in this basis, and so one cannot read off the kinetic part of $H$
from Eq. (23) directly. To find an operator which satisfies Eq. (23), we
introduce unitary operators defined on each link,

$\hat{R}_{n,i}(g)|U_{m,j}\rangle=\begin{cases}|g\,U_{m,j}\rangle\mbox{\,\,\,when\,\,}m=n\mbox{\,\,and\,\,}i=j\\\
|U_{m,j}\rangle\mbox{\,\,\,\,\,\,\,otherwise}\end{cases}$ (24)

which can be written with Hermitian operators $\hat{\pi}$ as

$\hat{R}_{n,i}(g)=e^{ix^{a}\hat{\pi}^{a}_{n,i}}.$ (25)

In short, $\hat{\pi}$ generate rotations $\hat{R}$ of each link. Here, in the
case of $SU(N)$, the $x^{a}$ are $N^{2}-1$ real parameters parameterizing the
group element $g\in SU(N)$, and $\hat{\pi}^{a}_{n,i}$ are eight Hermitian
operators (differential operators on the group $SU(N)$) associated to $x^{a}$
defined on each link. Using these operators, we rewrite the transfer matrix

$\displaystyle\hat{T}=\int_{g\in G}$
$\displaystyle\prod_{n,i}\left[dg_{n,i}\;\hat{R}_{n,i}(g_{n,i})\right]$
$\displaystyle\times
e^{i\frac{a}{g_{s}^{2}a_{0}}\sum_{n,i}\mbox{Re}\mbox{Tr}\left[1-g^{\dagger}_{n,i}\right]+i\hat{V}}\text{.}$
(26)

Exchanging $g$ for real parameters $x^{a}$, one obtains

$\displaystyle\hat{T}=\int_{-\pi}^{\pi}$
$\displaystyle\prod_{n,i,a}\left[dx_{n,i}^{a}\;e^{ix_{n,i}^{a}\hat{\pi}_{n,i}^{a}}\right]$
(27) $\displaystyle\times$ $\displaystyle
e^{i\frac{a}{g_{s}^{2}a_{0}}\sum_{n,i}\mbox{Re}\mbox{Tr}\left[1-e^{ix_{n,i}^{a}\lambda^{a}}\right]+i\hat{V}}\text{.}$

In the limit $a_{0}\rightarrow 0$, the integral can be evaluated via the
saddle point method. The saddle point in this case is degenerate, being the
gauge orbit of the point $x=0$. As the whole expression of $\hat{T}$ is gauge
invariant, we take the saddle-point approximation around $x=0$ to obtain

$\displaystyle\hat{T}$ $\displaystyle\sim$ $\displaystyle\int
dx\;e^{ix_{\rho}\hat{\pi}_{\rho}+\frac{ia}{2g_{s}^{2}a_{0}}x_{\rho}x_{\rho}+i\hat{V}}=Ne^{-ia_{0}\frac{g_{s}^{2}}{2a}\hat{\pi}_{\rho}\hat{\pi}_{\rho}+i\hat{V}}\text{.}$
(28)

For brevity we have written combined three indices into $\rho\equiv(n,i,a)$.
Together with the spatial plaquette terms from $V$, the Hamiltonian is

$\hat{H}_{KS}=\frac{g_{s}^{2}}{2a}\sum_{n,i,a}\hat{\pi}_{n,i}^{a}\hat{\pi}_{n,i}^{a}-\frac{1}{g_{s}^{2}a}\sum_{\begin{subarray}{c}n\\\
i<j\end{subarray}}\mbox{Re}\mbox{Tr}\left[\hat{P}_{ij}(n)\right]\text{.}$ (29)

The first term in the Hamiltonian can also be written as
$\mbox{Tr}[\hat{\pi}^{2}_{n,i}]$ by defining the $N\times N$ matrix
$\hat{\pi}_{n,i}=\hat{\pi}^{a}_{n,i}\lambda^{a}$. Under gauge transformations,
we have

$\hat{\pi}_{n,i}\rightarrow g_{n}^{-1}\hat{\pi}_{n,i}g_{n}\text{,}$ (30)

and this implicitly defines the transformation law for the operators
$\hat{\pi}^{a}_{n,i}$ as well. Note that the Hamiltonian in either form is
Hermitian and manifestly gauge-invariant.

In the rest of the section, we repeat this procedure in the presence of a
perturbation $\mathcal{O}$ — in the sections below this perturbation will be
$T_{\mu\nu}$. In the action formulation, one way to measure the expectation
value of $\mathcal{O}$ at time $t_{0}$ in a system governed by the action $S$
is to perturb the action by the function $\mathcal{O}(U)$ and define

$Z_{\epsilon}=\int\mathcal{D}U\;e^{i(S_{0}+\epsilon\mathcal{O}(t_{0}))}\text{.}$
(31)

Differentiating $Z_{\epsilon}$ — or rather its logarithm, the free energy —
with respect to $\epsilon$ yields:

$i\langle\mathcal{O}(t_{0})\rangle=Z_{0}^{-1}\left(\frac{\partial
Z_{\epsilon}}{\partial\epsilon}\right)_{\epsilon\rightarrow 0}\text{.}$ (32)

A perturbed Hamiltonian $H_{\epsilon}$, yet to be determined, corresponds to
this perturbed path integral. The corresponding transfer matrix has matrix
elements given by

$\langle
U^{\prime}|e^{-ia_{0}H_{\epsilon}}|U\rangle=e^{i(S+\epsilon\mathcal{O})}\text{.}$
(33)

Differentiating the transfer matrix with respect to $\epsilon$ yields $i$
times the desired operator; therefore we see the Hamiltonian should be
perturbed to

$\hat{H}_{\epsilon}=\hat{H}_{0}-\epsilon\hat{O}/a_{0}\text{.}$ (34)

So, by perturbing the Lagrangian with a component of $T_{\mu\nu}$ and finding
the corresponding perturbed Hamiltonian $\hat{H}_{\epsilon}$, one can read off
the operators $\hat{T}_{\mu\nu}$. To be precise, if the perturbing parameter
is $\epsilon$, then the corresponding operator in the Hamiltonian formalism is
the coefficient of $(\epsilon)$.

###  III.3 $T_{\mu\mu}$ in the Hamiltonian formulation

The diagonal components of the EMT consist of two kinds of terms:
$F_{0i}F_{0i}$ and $F_{ij}F_{ij}$. Without loss of generality, let us consider
$T_{11}$ in two spatial dimensions. Perturbing the action by terms
proportional to $F_{01}(n_{0})^{2}$ and $F_{12}(n_{0})^{2}$, we have

$S_{\epsilon}=S_{W}+\epsilon
a_{0}a^{3}\mbox{Tr}\left[F_{01}(n_{0})^{2}+F_{12}(n_{0})^{2}\right]\text{.}$
(35)

To simplify notation, $S_{\epsilon}$, $H_{\epsilon}$, and $T_{\epsilon}$ will
have different meanings in this and each subsequent subsection, corresponding
to the different types of perturbations being considered.

The perturbing terms are defined on the spacetime lattice via Eq. (11). For
$H_{\epsilon}$ to give the same dynamics as by the action $S_{\epsilon}$, they
should be connected via transfer matrix $T_{\epsilon}=e^{-ia_{0}H_{\epsilon}}$
as

$\langle
U^{\prime}|T_{\epsilon}|U\rangle=e^{iS_{\epsilon}}=e^{i(K_{\epsilon}+V_{\epsilon})}$
(36)

with

$\displaystyle K_{\epsilon}(U^{\prime},U)$ $\displaystyle=$ $\displaystyle
K+\epsilon\frac{2a}{g_{s}^{2}a_{0}}\mbox{Re}\mbox{Tr}\left[1-U_{n_{0},1}U^{{}^{\prime}\dagger}_{n_{0},1}\right]$
(37) $\displaystyle V_{\epsilon}(U^{\prime},U)$ $\displaystyle=$
$\displaystyle
V+\epsilon\frac{2a_{0}}{g_{s}^{2}a}\mbox{Re}\mbox{Tr}\left[1-P_{12}(n_{0})\right]$
(38)

The spatial plaquettes in Eq. (38) correspond to diagonal operators; it is
only for Eq. (37) that the transfer matrix formalism is useful. A little
algebra verifies that this Lagrangian corresponds to the transfer matrix

$\hat{T}_{\epsilon}=\int\mathcal{D}g\;e^{iK(g)+i\epsilon\frac{2a}{g_{s}^{2}a_{0}}\text{ReTr}[1-g_{n_{0},1}^{\dagger}]+i\hat{V}_{\epsilon}}$
(39)

where we have introduced shorthand

$\displaystyle\int\mathcal{D}g$ $\displaystyle\equiv$ $\displaystyle\int_{g\in
G}\prod_{n,i}dg_{n,i}\hat{R}_{n,i}(g_{n,i})$ (40) $\displaystyle K(g)$
$\displaystyle\equiv$
$\displaystyle\frac{a}{g_{s}^{2}a_{0}}\sum_{n,i}\mbox{Re}\left[1-g^{\dagger}_{n,i}\right].$
(41)

Performing the saddle-point approximation gives:

$\displaystyle\hat{T}_{\epsilon}$ $\displaystyle\sim$ $\displaystyle\int
dx\;e^{ix_{\rho}\hat{\pi}_{\rho}-x_{\rho}M_{\rho\sigma}x_{\sigma}+i\hat{V}_{\epsilon}}$
$\displaystyle\rm{with\,\,}$ $\displaystyle
M_{\rho\sigma}=\frac{-ia}{2g_{s}^{2}a_{0}}\delta_{\rho\sigma}-\epsilon\frac{ia}{g_{s}^{2}a_{0}}\delta_{nn_{0}}\delta_{mn_{0}}\delta_{i1}\delta_{j1}\delta_{ab}\text{.}$

Here we have abbreviated $\rho\equiv(n,i,a)$ and $\sigma\equiv(m,j,b)$. The
integral gives to $O(\epsilon)$,

$\hat{H}_{\epsilon}=\hat{H}_{KS}-\epsilon\left(\frac{g_{s}^{2}}{a}\mbox{Tr}\left[\hat{\pi}_{n_{0},1}^{2}\right]+\frac{2}{g_{s}^{2}a}\mbox{Re}\mbox{Tr}\left[1-\hat{P}_{12}(n_{0})\right]\right)$
(43)

We read off the operators for $F_{\mu\nu}F_{\mu\nu}$:

$\displaystyle\mbox{Tr}\left[\hat{F}_{0i}(n_{0})^{2}\right]$ $\displaystyle=$
$\displaystyle\frac{g_{s}^{2}}{a^{4}}\mbox{Tr}\left[\hat{\pi}_{n_{0},i}^{2}\right]$
(44) $\displaystyle\mbox{Tr}\left[\hat{F}_{ij}(n_{0})^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{2}{g_{s}^{2}a^{4}}\mbox{Re}\mbox{Tr}\left[1-\hat{P}_{ij}(n_{0})\right]\text{.}$
(45)

Finally, we can construct $T_{\mu\mu}$ from these operators:

$\displaystyle\hat{T}_{00}(n_{0})$ $\displaystyle=$
$\displaystyle\frac{g_{s}^{2}}{2a^{4}}\sum_{i}\mbox{Tr}[\hat{\pi}_{n_{0},i}^{2}]$
(46) $\displaystyle+$
$\displaystyle\frac{1}{g_{s}^{2}a^{4}}\sum_{i<j}\mbox{Re}\mbox{Tr}[1-\hat{P}_{ij}(n_{0})]$
$\displaystyle\hat{T}_{ii}(n_{0})$ $\displaystyle=$
$\displaystyle\frac{g_{s}^{2}}{2a^{4}}\mbox{Tr}\left[-\hat{\pi}_{n_{0},i}^{2}+\hat{\pi}_{n_{0},j}^{2}+\hat{\pi}_{n_{0},k}^{2}\right]$
$\displaystyle+\frac{1}{g_{s}^{2}a^{4}}$ Re
$\displaystyle\mbox{Tr}\left[1-\hat{P}_{ij}(n_{0})-\hat{P}_{ik}(n_{0})+\hat{P}_{jk}(n_{0})\right]$
(47)

Note that $T_{00}(n)$ is the density of the Kogut-Suskind Hamiltonian $H_{KS}$
up to a constant term. Because we have worked only at tree level, the trace of
the EMT operator vanishes as is the case in Eq. (8).

Eq. (46) gives $\hat{T}_{\mu\mu}$ up to $O(a)$. To improve the
$\hat{T}_{\mu\mu}$ operators up to $O(a^{2})$, we use Eq. (III.1) and Eq.
(III.1) and take the average around the site $n_{0}$:

$\displaystyle\mbox{Tr}\left[\hat{F}_{0i}(n_{0})^{2}\right]$ $\displaystyle=$
$\displaystyle\sum_{x=0,1}\frac{g_{s}^{2}}{2a^{4}}\mbox{Tr}\left[\hat{\pi}_{n_{0}-x\hat{i},i}^{2}\right]$
(48) $\displaystyle\mbox{Tr}\left[\hat{F}_{ij}(n_{0})^{2}\right]$
$\displaystyle=$
$\displaystyle\sum_{x=0,1}\sum_{y=0,1}\frac{1}{2g_{s}^{2}a^{4}}$ (49)
$\displaystyle\mbox{Re}\mbox{Tr}\left[1-\hat{P}_{ij}(n_{0}-x\hat{i}-y\hat{j})\right]\text{.}$

These operators enable us to construct $\hat{T}_{\mu\mu}$ up to discretization
errors that are $O(a^{2},a_{0})$.

###  III.4 $T_{ij}$ in the Hamiltonian formulation

Let us now move to deriving the operators $\hat{T}_{ij}$, that is, the off-
diagonal spatial parts of the EMT:

$T_{ij}=\mbox{Tr}\left[-F_{i0}F_{j0}+F_{ik}F_{jk}\right]\text{.}$ (50)

This definition holds both on the spacetime lattice and as an operator
equation on the Hamiltonian lattice. We first work with the naive
discretization, and then with the clover discretization. Without loss of
generality, let us take $T_{12}$ as an example and perturb the Wilson action
with terms in Eq. (50). We find that $\hat{T}_{\epsilon}$ is given by Eq. (36)
with:

$\displaystyle K_{\epsilon}$ $\displaystyle=$ $\displaystyle K+\epsilon
a_{0}a^{3}\mbox{Tr}\left[F^{N}_{10}F^{N}_{20}\right]$ (51) $\displaystyle
V_{\epsilon}$ $\displaystyle=$ $\displaystyle V+\epsilon
a_{0}a^{3}\mbox{Tr}\left[F^{N}_{13}F^{N}_{23}\right]\text{.}$ (52)

As before, the spatial plaquettes in Eq. (52) can be directly converted to
operators. The time-like plaquettes in Eq. (51) will ultimately appear as
various $\hat{\pi}$. Using $\hat{R}(g)$ operators, $\hat{T}_{\epsilon}$ can be
written

$\hat{T}_{\epsilon}=\int\mathcal{D}g\;e^{iK(g)-i\epsilon\frac{a}{4g_{s}^{2}a_{0}}\text{Tr}\left[(g_{n_{0},1}^{\dagger}-g_{n_{0},1})(g_{n_{0},2}^{\dagger}-g_{n_{0},2})\right]+i\hat{V}_{\epsilon}}\text{.}$
(53)

Evaluating the integral via the saddle point $x=0$ (exact in the limit
$a_{0}\rightarrow 0$) gives:

$\displaystyle\hat{T}_{\epsilon}$ $\displaystyle\sim$ $\displaystyle\int
dx\;e^{ix_{\rho}\hat{\pi}_{\rho}-x_{\rho}M_{\rho\sigma}x_{\sigma}+i\hat{V}_{\epsilon}}=Ae^{-\frac{1}{4}\hat{\pi}_{\rho}M^{-1}_{\rho\sigma}\hat{\pi}_{\sigma}+i\hat{V}_{\epsilon}}$
(54) $\displaystyle
M_{\rho\sigma}=\frac{-ia}{2g_{s}^{2}a_{0}}\delta_{\rho\sigma}+\epsilon\frac{ia}{g_{s}^{4}a_{0}}\delta_{nn_{0}}\delta_{mn_{0}}\delta_{i1}\delta_{j2}\delta_{ab}$

which at $O(\epsilon)$ yields $\hat{H}_{\epsilon}$:

$\hat{H}_{\epsilon}=\hat{H}_{KS}-\epsilon\frac{g_{s}^{2}}{a}\mbox{Tr}[\hat{\pi}_{n_{0},1}\hat{\pi}_{n_{0},2}]-\epsilon
a^{3}\mbox{Tr}\left[\hat{F}^{N}_{ik}(n_{0})\hat{F}_{jk}^{N}(n_{0})\right]$
(55)

where the second term in RHS correspond to $a_{0}a^{3}F_{10}F_{20}$. More
generally, operators for $F_{i0}F_{j0}$ are

$\mbox{Tr}\left[\hat{F}^{N}_{i0}\hat{F}^{N}_{j0}(n_{0})\right]=\frac{g_{s}^{2}}{a^{4}}\mbox{Tr}[\hat{\pi}_{n_{0},i}\hat{\pi}_{n_{0},j}]$
(56)

Thus the naive $\hat{T}_{ij}(n_{0})$ in the Hamiltonian formulation is

$\hat{T}^{N}_{ij}(n_{0})=-\frac{g_{s}^{2}}{a^{4}}\mbox{Tr}\left[\hat{\pi}_{n_{0},i}\hat{\pi}_{n_{0},j}\right]+\mbox{Tr}\left[\hat{F}^{N}_{ik}(n_{0})\hat{F}^{N}_{jk}(n_{0})\right]\text{.}$
(57)

The clover approximations are obtained from $F_{ij}^{C}$ in Eq. (16) and
$F_{i0}^{B}$ in Eq. (18). As before, the transition from the action formalism
to the Hamiltonian is straightforward for $F_{ij}$, so we focus only on the
$F_{10}F_{20}$ term. For these,

$\displaystyle K_{\epsilon}(U^{\prime},U)$ $\displaystyle=$ $\displaystyle
K+\epsilon
a_{0}a^{3}\mbox{Tr}\left[F^{B}_{10}(n_{0})F^{B}_{20}(n_{0})\right].$ (58)

We use the definitions of Fig. 2 for the links around $n_{0}$. For the example
of $U_{n_{0},1}$, we denote operators and functions on them as
$U_{1},\hat{U}_{1},\hat{\pi}_{1}$, and $g_{1}=e^{ix_{1}^{a}\lambda^{2}}$.

Figure 2: Half-clovers $B_{10}(n_{0})$ and $B_{20}(n_{0})$ at site $n_{0}$.

Then $\hat{T}_{\epsilon}$ is

$\hat{T}_{\epsilon}=\int\mathcal{D}g\,e^{iK(g)-i\frac{\epsilon
a}{16g_{s}^{2}a_{0}}\text{Tr}\left[(g_{1}^{\dagger}-g_{1}+\hat{U}_{0}^{\dagger}(g_{0}^{\dagger}-g_{0})\hat{U}_{0})(g_{3}^{\dagger}-g_{3}+\hat{U}_{2}^{\dagger}(g_{2}^{\dagger}-g_{2})\hat{U}_{2})\right]+i\hat{V}}$
(59)

After the saddle-point approximation around $x=0$, $\hat{T}_{\epsilon}$
simplifies and becomes

$\displaystyle\hat{T}_{\epsilon}$ $\displaystyle\sim$ $\displaystyle\int
dx\;e^{ix_{\rho}\hat{\pi}_{\rho}-x_{\rho}M_{\rho\sigma}x_{\sigma}+i\hat{V}}=Ae^{-\frac{1}{4}\hat{\pi}_{\rho}M^{-1}_{\rho\sigma}\hat{\pi}_{\sigma}+i\hat{V}}$
(60)
$\displaystyle\mbox{\,\,with\,\,}M_{\rho\sigma}=-\frac{ia}{2g_{s}^{2}a_{0}}\delta_{\rho\sigma}-\frac{i\epsilon
a}{4g_{s}^{2}a_{0}}(M_{1})_{\rho\sigma}$

Matrix elements of $(M_{1})_{\rho\sigma}$ are zero other than

$\displaystyle\left(M_{1}\right)_{(n_{0},1,a)(n_{0},2,b)}$
$\displaystyle=\delta_{ab}$ (61)
$\displaystyle\left(M_{1}\right)_{(n_{0},1,a)(n_{2},2,b)}$
$\displaystyle=\mbox{Tr}\left[\lambda^{a}\hat{U}_{2}^{\dagger}\lambda^{b}\hat{U}_{2}\right]$
$\displaystyle\left(M_{1}\right)_{(n_{1},1,a)(n_{0},2,b)}$
$\displaystyle=\mbox{Tr}\left[\hat{U}_{0}^{\dagger}\lambda^{a}\hat{U}_{0}\lambda^{b}\right]$
$\displaystyle\left(M_{1}\right)_{(n_{1},1,a)(n_{2},2,b)}$
$\displaystyle=\mbox{Tr}\left[\hat{U}_{0}^{\dagger}\lambda^{a}\hat{U}_{0}\hat{U}_{2}^{\dagger}\lambda^{b}\hat{U}_{2}\right]$

where sites $n_{1},n_{2}$ are as labeled in Fig. 2. Now by expanding Eq. (60)
to linear order in $\epsilon$, we obtain the perturbed Hamiltonian

$\displaystyle H_{\epsilon}=$ $\displaystyle
H_{K,S}-\epsilon\frac{g_{s}^{2}}{4a}\bigg{(}\mbox{Tr}\left[\hat{\pi}_{1}\hat{\pi}_{3}\right]+\mbox{Tr}\left[\hat{\pi}_{1}\hat{U}_{2}^{\dagger}\hat{\pi}_{2}\hat{U}_{2}\right]$
$\displaystyle+\mbox{Tr}\left[\hat{U}_{0}^{\dagger}\hat{\pi}_{0}U_{0}\hat{\pi}_{3}\right]+\mbox{Tr}\left[\hat{U}_{0}^{\dagger}\hat{\pi}_{0}\hat{U}_{0}\hat{U}_{2}^{\dagger}\hat{\pi}_{2}\hat{U}_{2}\right]\bigg{)}.$
(62)

Thus $F_{i0}F_{j0}$ is generally implemented as

Tr
$\displaystyle\left[\hat{F}^{B}_{i0}\hat{F}^{B}_{j0}(n_{0})\right]=\frac{g_{s}^{2}}{4a^{4}}\bigg{(}\mbox{Tr}\left[\hat{\pi}_{n_{0},i}\hat{\pi}_{n_{0},j}\right]$
$\displaystyle+\mbox{Tr}\left[\hat{\pi}_{n_{0},i}\hat{U}_{n_{0}-\hat{j},j}^{\dagger}\hat{\pi}_{n_{0}-\hat{j},j}\hat{U}_{n_{0}-\hat{j},j}\right]$
$\displaystyle+\mbox{Tr}\left[\hat{U}_{n_{0}-\hat{i},i}^{\dagger}\hat{\pi}_{n_{0}-\hat{i}}\hat{U}_{n_{0}-\hat{i},i}\hat{\pi}_{n_{0},j}\right]$
$\displaystyle+\mbox{Tr}\left[\hat{U}_{n_{0}-\hat{i},i}^{\dagger}\hat{\pi}_{n_{0}-\hat{i},i}\hat{U}_{n_{0}-\hat{i},i}\hat{U}_{n_{0}-\hat{j},j}^{\dagger}\hat{\pi}_{n_{0}-\hat{j},j}\hat{U}_{n_{0}-\hat{j},j}\right]\bigg{)}.$
(63)

The spatial $F_{ik}F_{ik}$ via clovers can be directly converted to operators.
With these operators, operators for measuring the spatial off-diagonal
components of $T_{\mu\nu}$ are fully constructed following Eq. (50).

###  III.5 $T_{0i}$ in the Hamiltonian formulation

In this subsection we derive Hermitian operators for $T_{0i}$, whose all terms
contain time-like plaquettes and thus need to be appropriately converted for
quantum simulations. As an example, $T_{01}$, via naive discretization, is
written as:

$\displaystyle T_{01}$ $\displaystyle=$
$\displaystyle\mbox{Tr}\left[F_{02}F_{12}+F_{03}F_{13}\right]$ (64)
$\displaystyle=$
$\displaystyle\mbox{Tr}\left[F^{N}_{02}F^{N}_{12}+F^{N}_{03}F^{N}_{13}\right]+O(a)$
(65)

To find operators for $T_{01}$, As two terms are in the same form
$F^{N}_{0j}F^{N}_{ij}$, let us perturb the Wilson action with only
$F^{N}_{02}F^{N}_{12}$ at particular site $n_{0}$ to obtain

$\displaystyle K_{\epsilon}(U^{\prime},U)$ $\displaystyle=$ $\displaystyle
K+\epsilon a^{3}a_{0}\mbox{Tr}\left[F^{N}_{02}F^{N}_{12}\right]$ (66)

from which we find the following $\hat{T}_{\epsilon}$:

$\hat{T}_{\epsilon}=\int\mathcal{D}g\,\,e^{iK(g)-\frac{i\epsilon}{4g_{s}^{2}}\text{Tr}\left[(g_{n_{0},2}-g_{n_{0},2}^{\dagger})(\hat{P}_{12}(n_{0})-\hat{P}_{12}^{\dagger}(n_{0}))\right]+i\hat{V}}.$
(67)

In the limit of $a_{0}\rightarrow 0$, one can approximate the integral via
saddle-point approximation around $x=0$:

$\displaystyle\hat{T}_{\epsilon}$ $\displaystyle\sim$ $\displaystyle\int
dx\;e^{\frac{ia}{2g_{s}^{2}a_{0}}x_{\rho}\delta_{\rho\sigma}x_{\sigma}+ix_{\rho}\hat{\pi}^{\prime}_{\rho}+i\hat{V}}$
(68) $\displaystyle\text{where }\hat{\pi}^{\prime}_{\rho}$ $\displaystyle=$
$\displaystyle\hat{\pi}_{\rho}+\frac{\epsilon}{g_{s}^{2}}\delta_{nn_{0}}\delta_{i2}\mbox{Tr}\left[\lambda^{b}\mbox{Im}\hat{P}_{12}(n_{0})\right]\text{.}$

Evaluating the Gaussian integral gives, at $O(\epsilon)$,

$\hat{H}_{\epsilon}=\hat{H}_{KS}+\epsilon\frac{1}{a}\mbox{Tr}\left[\hat{\pi}_{n_{0},2}\mbox{Im}\hat{P}_{12}(n_{0})\right]\text{.}$
(69)

From the perturbed Hamiltonian, we read off

$\mbox{Tr}\left[\hat{F}^{N}_{02}\hat{F}^{N}_{12}(n_{0})\right]=-\frac{1}{a^{4}}\mbox{Tr}\left[\hat{\pi}_{n_{0},2}\mbox{Im}\hat{P}_{12}(n_{0})\right]\text{.}$
(70)

Therefore, the operator to give $T_{0i}$ with naive discretization is

$\hat{T}_{0i}(n_{0})=-\sum_{j\neq
i}\frac{1}{a^{4}}\mbox{Tr}\left[\hat{\pi}_{n_{0},j}\mbox{Im}\hat{P}_{ij}(n_{0})\right]\text{.}$
(71)

The naive discretization induces $O(a)$ error in $T_{0i}$. To improve up to
$O(a^{2})$, we use the clovers, $F_{ij}^{C}$ in Eq. (16) and $F_{i0}^{B}$ in
Eq. (18) instead. Again as an example, we add $F_{02}F_{12}$ to the Wilson
action to derive the transfer matrix from a perturbation given via

$\displaystyle K_{\epsilon}(U^{\prime},U)$ $\displaystyle=$ $\displaystyle
K+\epsilon a^{3}a_{0}\mbox{Tr}\left[F^{B}_{02}(n_{0})F^{C}_{12}(n_{0})\right]$
(72)

In the following we introduce notation for links around the site $n_{0}$ as in
Fig.3.

Figure 3: Clover $C_{12}(n_{0})$ and half-clover $B_{02}(n_{0})$ at site
$n_{0}$.

The transfer matrix is then

$\hat{T}_{\epsilon}=\int\mathcal{D}g\;e^{iK(g)-i\frac{\epsilon}{8g_{s}^{2}}\text{Tr}\left[(g_{1}-g_{1}^{\dagger}+\hat{U}_{0}^{\dagger}(g_{0}-g_{0}^{\dagger})\hat{U}_{0})(\hat{C}_{12}(n_{0})-\hat{C}_{12}^{\dagger}(n_{0}))\right]+i\hat{V}}$
(73)

where $\hat{C}$ is the clover operator. In the limit of $a_{0}\rightarrow 0$,
one can evaluate the integral via the saddle point around $x=0$:

$\displaystyle\hat{T}_{\epsilon}\sim\int$ $\displaystyle dx$
$\displaystyle\;e^{\frac{ia}{2g^{2}a_{0}}x_{\rho}\delta_{\rho\sigma}x_{\sigma}+ix_{\rho}\hat{\pi}^{\prime}_{\rho}+i\hat{V}}$
(74) $\displaystyle\hat{\pi}^{\prime}_{\rho}=\hat{\pi}_{\rho}$
$\displaystyle+$
$\displaystyle\frac{\epsilon}{2g_{s}^{2}}\delta_{nn_{0}}\delta_{i2}\left(\mbox{Tr}\left[\lambda^{a}\mbox{Im}\hat{C}_{12}(n_{0})\right]\right.$
(75) $\displaystyle+$
$\displaystyle\left.\mbox{Tr}\left[\hat{U}_{0}^{\dagger}\lambda^{a}\hat{U}_{0}\mbox{Im}\hat{C}_{12}(n_{0})\right]\right)$

which yields to $O(\epsilon)$,

$\displaystyle\hat{H}_{\epsilon}=\hat{H}_{KS}$ $\displaystyle+$
$\displaystyle\epsilon\frac{1}{2a}\left(\mbox{Tr}\left[\hat{\pi}_{1}\mbox{Im}\hat{C}_{12}(n_{0})\right]\right.$
(76) $\displaystyle+$
$\displaystyle\left.\mbox{Tr}\left[\hat{U}_{0}^{\dagger}\hat{\pi}_{0}\hat{U}_{0}\mbox{Im}\hat{C}_{12}(n_{0})\right]\right)$

From $\hat{H}_{\epsilon}$, we read off, for general $F_{0j}F_{ij}(n_{0})$,

$\displaystyle\mbox{Tr}\left[\hat{F}^{B}_{0j}\hat{F}^{C}_{ij}(n_{0})\right]=-\frac{1}{2a^{4}}\left(\mbox{Tr}\left[\hat{\pi}_{n_{0},j}\mbox{Im}\hat{C}_{ij}(n_{0})\right]\right.$
$\displaystyle+\left.\mbox{Tr}\left[\hat{U}_{n_{0}-\hat{j},j}^{\dagger}\hat{\pi}_{n_{0}-\hat{j},j}\hat{U}_{n_{0}-\hat{j},j}\mbox{Im}\hat{C}_{ij}(n_{0})\right]\right)\text{.}$
(77)

From this operator, $T_{0i}(n)$ can be fully constructed.

##  IV QGP state preparation

We have described the operators, acting on the Hamiltonian lattice, that
provide the shear viscosity transport coefficients. In order to evaluate these
expectation values, we need circuits corresponding to time evolution under the
physical Hamiltonian; such circuits are described for a general gauge theory
in [95]. In addition we need to prepare states for them to act on which will
enable one to reproduce the appropriate thermal expectation values. The
purpose of this section is to describe several practical methods for this
thermal state preparation.

Before describing the methods, it is worth examining whether there are _in
principle_ difficulties preparing a thermal state, such that the preparation
necessarily requires time exponential in the volume. This is famously the case
for certain frustrated spin systems [96] at low temperatures. If a similar
barrier exists for thermal state preparation in gauge theories, our methods
will be of no practical use.

To see there is unlikely to be a problem with the preparation of a QCD (or
Yang-Mills) thermal state, it is useful to appeal to experiment. A key feature
of heavy-ion phenomenology is that after the initial hard interactions, the
system rapidly equilibrates into a system that is well-described by
hydrodynamics. The equilibration time is believed to be $\sim 1\;\mathrm{fm}$;
comparable to the characteristic scale of QCD. It is highly plausible that
such rapid equilibration is a generic feature of strongly-coupled theories
[97, 12]. Accordingly, one expects that if an out-of-equilibrium state is
prepared, it equilibrates on roughly the natural time scales of the theory.
Thus, it seems implausible that preparing a thermal state would take an
exponentially long time.

It is particularly convenient that the thermal states of interest to us have
$T\sim 200\;\mathrm{MeV}$. Very cold temperatures, near the true ground state,
may be difficult to prepare due to frustration. Very high temperature states
($T\gg\Lambda_{\mathrm{QCD}}$) likely have difficulty thermalizing, as the
fluid is closer to a free gas. At $T\sim 200\;\mathrm{MeV}$, the strongly
interacting fluid is at a temperature comparable to other physical scales.
Even without experimental evidence, we might expect such systems to
equilibrate quickly.

Given this discussion it is natural to assume that polynomial-time thermal
state preparation is possible — even lacking a formal proof. With this
assumption, we now turn to the task of finding a practical way to construct
such a state.

In this section several approaches are explored to prepare thermal states. We
explore multiple paths since at this stage the relative merits with respect to
future quantum devices is unknown. One general note applies to all methods
discussed here. Common digitization schemes result in many possible states on
the quantum computer that correspond to no physical state. An important
example is gauge invariance: the physical states are those that are unchanged
by any gauge transformation. Naively, this is a serious difficulty, as we need
to ensure that the prepared state is one of the rare gauge-invariant ones.
However, as discussed in [95, 98, 99] and elsewhere, because time-evolution
(even under a Suzuki-Trotter approximation) commutes with gauge
transformations, any state preparation method that builds a gauge-invariant
state and then performs time-evolution will automatically respect the gauge
symmetry up to quantum noise.

###  IV.1 Thermal states

Before looking at detailed schemes, it behooves us to first consider what we
mean by a thermal state. Ideally, we want a single quantum state in which
expectation values match those given by the canonical ensemble
$\langle\mathcal{O}\rangle=Z^{-1}\mbox{Tr}\;e^{-\beta H}\mathcal{O}$. In such
a state the physical system — in this case of the appropriately truncated
version of the lattice gauge theory — is modeled by a subspace of a larger
system. We can define a thermal state as

$|T\rangle=\sum_{k}c_{k}|\psi_{k}\rangle_{\rm sys}\otimes|\psi_{k}\rangle_{\rm
comp}$ (78)

where a subscript of $\mathrm{sys}$ denotes a state of the physical system
being studied, and one of $\mathrm{comp}$ denotes the state of the
complementary system. The coefficients $c_{k}$ are fixed by considering what
happens when we trace out the complementary system; in particular, we require

$\rho_{T}^{\rm sys}\equiv{\rm Tr^{\rm comp}}\left[|T\rangle\langle
T|\right]=\frac{\sum_{k}e^{-\beta E_{k}}|\psi_{k}\rangle_{\rm sys}\,{}_{\rm
sys}\langle\psi_{k}|}{\sum_{j}e^{-\beta E_{j}}}$ (79)

so that for an operator ${\cal O}^{\rm sys}$ that acts entirely in the system
subspace, the quantum expectation values match the desired thermal ones. Thus,
the density matrix appears as the classically uncertain state of a quantum
system once the complementary system has been traced out.

Obtaining a true thermal state is difficult, but generally unnecessary: a good
approximation is sufficient. In fact, there are a wide variety of possible
ensembles, all of which agree in the thermodynamic limit — see [100] for an
elegant example. By exploiting the equivalence between microcanonical and
canonical descriptions444Here, we are not considering other conserved
quantities in the system so that the grand canonical ensemble is not relevant.
in the thermodynamic limit, we should be able to learn about systems even when
we cannot obtain a distribution of energies with probabilities given by
Boltzmann weights.

In the microcanonical approach to statistical mechanics, one computes
properties at a fixed energy rather than at a fixed temperature. Given the
standard assumptions relating statistical mechanics to thermodynamics, in the
thermodynamic limit of large systems, quantities computed in the
microcanonical and canonical descriptions will agree [101]. For ‘typical’
quantum field theories (presumably including both Yang-Mills and QCD), given
an eigenstate $|\Psi\rangle$ of the system Hamiltonian $H^{\mathrm{sys}}$
sampled from the set of all eigenstates with energy density is near
$\epsilon$, expectation values $\langle\Psi|\mathcal{O}|\Psi\rangle$ will
agree in the thermodynamic limit with the canonical values.

This suggests a cheap approach to obtain thermodynamic expectation values.
Rather than preparing a mixed state that exactly reproduces the canonical
ensemble, we prepare a typical pure state with the desired energy density.
This can be done, in practice, by preparing _any_ pure state with the desired
energy density, and then time evolving to allow the state to thermalize.

###  IV.2 The heat bath approach

While there is no practical way of constructing a perfect thermal state
$|T\rangle$, it is quite straightforward to find an algorithm to create a
reasonable approximation. The method exploits the physical principle of a heat
bath. The key idea is that one considers a very large total system — much
larger than the system subspace — with the complementary subspace serving as
heat bath. The heat bath is subject to dynamics under the Hamiltonian $H^{\rm
HB}$, which allows explicit construction of its eigenstates.

One starts with initial conditions in which the heat bath and the system are
decoupled; the heat bath is prepared in a known initial state of low energy
and the system is in some high-energy state. This physical state is chosen
primarily for simplicity — the details of the state will not matter provided
that the heat bath is sufficiently large. The heat bath and the system of
interest are dynamically coupled in a gauge invariant way via a Hamiltonian,
$H^{\rm couple}$, which is initially switched off. The coupling Hamiltonian
satisfies

$[H^{\rm couple},H^{\rm sys}]\neq 0\;,\;\;[H^{\rm couple},H^{\rm HB}]\neq
0\text{,}$ (80)

so that when it is switched on the coupling allows energy to flow between the
system and the heat bath. If one waits sufficiently long, one might reasonably
expect the system to approximately thermalize. The $H^{\rm couple}$ must be
small in the sense that at all stages in the evolution the absolute value of
its expectation value:

$\langle H^{\rm couple}\rangle\ll\langle H^{\rm sys}\rangle-\langle H^{\rm
sys}\rangle_{\rm vac}\;.$ (81)

where $\langle H^{\rm sys}\rangle_{\rm vac}$ is the expecation value of the
vacuum state. This condition ensures that the details of the coupling between
the heat bath and the system has a negligible effect on the final result.

Once the system has thermalized, one can switch off $H^{\rm couple}$. This is
essentially modeling the physical process by which physical systems
thermalize; and the basic assumption underlying statistical mechanics is that
the details of how the system thermalizes should not matter. By choosing
various initial configurations of the heat bath, one can evolve the system
into (approximately) thermalized systems at various temperatures.

While this method should work as a matter of principle, it has a strong
practical disadvantage: it requires a very large number of qubits to encode
the heat bath. Of course, the general notion of a heat bath in thermodynamics
is that it should be large. Moreover, in the present context there are some
particular challenges requiring this.

Recall that $|T\rangle$ formally requires a heat bath at least as large as the
system; otherwise the condition ${}_{\rm
HB}\langle\psi_{k}^{\prime}|\psi_{k}\rangle_{\rm HB}=\delta_{k^{\prime},k}$
cannot be met. However as a practical matter the heat bath is coupled to the
system in an passive way; it does not couple system states to heat bath states
in a manner that naturally pushes the system towards configurations satisfying
${}_{\rm HB}\langle\psi_{k}^{\prime}|\psi_{k}\rangle_{\rm
HB}=\delta_{k^{\prime},k}$. Rather to achieve some approximation to that
condition, one requires a large heat bath subspace in which approximate
orthogonality is likely to emerge naturally.

There is a second practical issue that suggests the need for large heat baths.
For a generic interacting theory, the only states we know how to write
explicitly in a practical way are states with energies at the scale of the
lattice spacing. Thus, to cool the system to a temperatures of physical
interest one needs to transfer a substantial energy from the system to the
heat bath. This in turn means that the heat capacity of the heat bath must be
large.

While heat bath must be large, it is clear that its size scales polynomially
in the physical size and lattice spacing of the system. Thus, the heat bath is
a useful demonstration of an explicit method with polynomial scaling. However,
given the large size of heat bath that the approach requires, it is unlikely
to be the optimal way to pursue thermal physics on quantum computers.

###  IV.3 Active cooling

A principal problem with the heat bath approach is that it is essentially
passive. The only active step is the connecting the system with the heat bath
through $H^{\rm couple}$. One reason the heat bath needs to be large in such a
passive scheme is that it must absorb all of the excess energy of the system,
which will start at an energy density of lattice scales.

Let us consider a slight variation on the heat bath approach. Suppose that
instead of a single large heat bath one had $N$ smaller heat baths each
starting in its ground state configuration. Moreover, let the system begin in
a state with $\langle H^{\rm sys}\rangle=E_{0}$. We connect the system to the
first heat bath alone, via a small coupling Hamiltonian. The coupled system
evolves dynamically for some time after which this coupling is switched off.
Assuming that the coupling Hamiltonian is small in sense of Eq. (81), the
dynamics must reduce the energy of the system: the energy of the coupled
system is conserved and the energy of the complementary system can only
increase since it starts in the ground state. Thus at the end of the time
which the first heat bath is coupled to the system $\langle H^{\rm
sys}\rangle=E_{1}$ with $E_{1}<E_{0}$.

At that stage, the system is decoupled from the first heat bath and coupled to
the second. We perform more time evolution such that $\langle H^{\rm
sys}\rangle=E_{2}$, with $E_{2}<E_{1}$. Continuing the process, we see that
the energy is monotonically decreasing:

$E_{N}<E_{N-1}\cdots<E_{2}<E_{1}<E_{0}\;.$ (82)

Including enough smaller heat bath allows the system to lose as much energy as
one wishes provided the coupling Hamiltonian is sufficiently small.

Superficially, this approach — like the heat bath approach — is essentially
passive. The only active step appears to be the coupling and decoupling of the
system with the various heat baths. This is misleading, however, there was
another active step: preparing the various heat baths to be in their ground
state. If one was not able to do that, the approach would not be viable.
However, if one can initialize these into their ground states, there is no
need to use multiple different “heat baths”; rather the same “heat bath” can
be reused and reinitialized to its ground state between cycles. The system
will yield precisely the same the results as it would have had one used
multiple “heat bath”. The virtue of this is that qubit costs needed obtain a
typical energy density for the system is greatly reduced compared to the heat
bath method.

In what follows, we refer not to a “heat bath” but suggestively to a “pump”
since it does not keep the system in thermal equilibrium. The pump can be of
similar size to the physical system — or smaller — and the key point is that
it absorbs energy from the system, moves the energy into the environment and
is reinitialized to its ground state. This step of reintitalizing the
complementary system is an active one. In effect the pumps is acting as
quantum heat pump or refrigerator [102] (for a review of quantum refrigerators
see [103]), thus the name.

Of course, there is no way for such a device to act entirely via unitary
evolution within the Hilbert space of qubits. The act of the pump dropping
back into its ground state is clearly not unitary. It is worth recalling that
quantum computing requires more than just control of the unitary behavior of
qubits: one of the DiVincenzo criteria [104] is the ability to initialize the
qubits to some fiducial state. This initialization is not unitary and
necessarily involves entangling the states of the Hilbert space of the quantum
computer with the environment. This method requires that the one can
separately initialize the system and the pump, and that the initialization
(which as we will see requires measuring the pump qubits) of the pump can be
done much faster than the physical system can decohere – a nontrival device
specification.

It is useful to have a concrete model in mind for the purposes of modeling.
Conceptually, if not practically, the simplest way to think of this is to
choose the fiducial state of the pump as $|g\rangle^{\rm pump}=|0,0,0,\cdots
0,0,0\rangle$ and then to choose $H^{\rm pump}$ so its ground state is the
fiducual state. One can reinitialize by simply measuring the qubits using
$\sigma_{z}$ to determine if the qubit is in $|1\rangle$ or $|0\rangle$ and if
the result is $|1\rangle$ act on it with the unitary operator $\sigma_{x}$. In
effect, this can be viewed as a qubit version of Maxwell’s demon, and a cycle
that employed such a device might be referred to as a quantum demon
refrigerator.

Let us provide a simple model of a cycle which will lower the system energy
which started from a state with a lattice-scale energy density. The Hilbert
space is spanned by outer product states of a space characterizing the system
and the pump, which we will now assume to be comparable in size. Again we have
three Hamiltonian, $H^{\rm sys}$ that acts on the system, $H^{\rm pump}$, that
acts in the pump space and $H^{\rm couple}$ that couples the two. The coupling
dynamics is constructed to be gauge invariant, $H^{\rm couple}$ must have non-
zero commutators with $H^{\rm sys}$ and $H^{\rm pump}$ and must be small in
the sense of satisfying inequality (81).

A schematic description of an active cooling cycle is given in Fig. 4. In step
I the time evolution has been switched off, the system begins in some pure or
mixed state with $\langle H^{\rm sys}\rangle=E_{I}$ relative to its ground
state, the pump subspace is then initialized to its ground state, with
$\langle H^{\rm pump}\rangle=0$. Since $H^{\rm couple}$ satisfies inequality
(81), the total energy (relative to the ground state) of the combined system
$\langle H^{\rm sys}+H^{\rm pump}+H^{\rm couple}\rangle$ , is also $E_{I}$ up
to a small correction. In step II, the time evolution associated with the
combined Hamiltonian is switched on suddenly, since this is sudden $\langle
H^{\rm sys}+H^{\rm pump}+H^{\rm couple}\rangle$ remains at $E_{I}$. In step
III the combined system undergoes time evolution under $H^{\rm sys}+H^{\rm
pump}+H^{\rm couple}$ for a fixed time. During this time evolution, the total
energy is conserved and any net energy flow goes from the system to the pump.
After some time, $\langle H^{sys}\rangle$ reaches some value $E_{III}$ where
$E_{III}<E_{I}$. The time evolution then is switched off suddenly leaving
$\langle H^{sys}\rangle=E_{III}<E_{I}$. At that point, the cycle repeats. Thus
each time through the cycle, $\langle H^{sys}\rangle$ drops. One could
continue iterating the cycle and lowering $\langle H^{sys}\rangle$ until the
inequality (81) ceases to hold.

Figure 4: An active cooling cycle.

Consider the reduced density matrix of Eq. (79) for the system as after
tracing over the pump. By construction step I preserves $\rho^{\rm sys}_{T}$.
This is obvious since during there is no coupling between the system and the
pump.

Suppose at step I the full density matrix, $\hat{\rho}^{I}$ is given by

$\hat{\rho}^{I}=\sum_{a,a^{\prime},b,b^{\prime}}|a^{\prime}\rangle^{\rm
sys}|b^{\prime}\rangle^{\rm
pump}\,\hat{\rho}_{a,b;a^{\prime}b^{\prime}}\,{}^{\rm pump}\langle b|{}^{\rm
sys}\langle a|$ (83)

where the states $|a\rangle^{\rm sys}$ and $|b\rangle^{\rm pump}$ are states
in an orthonormal basis for the system space and complementary space
respectively. Thus matrix elements of the reduced matrix elements are given by

$\hat{\rho}^{I\,\rm
sys}_{a,a^{\prime}}=\sum_{b,b^{\prime},b^{\prime\prime}}{}^{\rm pump}\langle
b^{\prime\prime}||b\rangle^{\rm
pump}\hat{\rho}_{a,b;a^{\prime}b^{\prime}}{}^{\rm pump}\langle
b^{\prime\prime}||b\rangle^{\rm pump}$ (84)

Thus the full density matrix at step II is given by

$\hat{\rho}^{II}=\sum_{a,a^{\prime}}|a^{\prime}\rangle^{\rm sys}|g\rangle^{\rm
pump}\,\left(\sum_{b,b^{\prime},b^{\prime\prime}}{}^{\rm pump}\langle
b^{\prime\prime}||b\rangle^{\rm
pump}\hat{\rho}_{a,b;a^{\prime}b^{\prime}}^{I}{}^{\rm pump}\langle
b^{\prime\prime}||b\rangle^{\rm pump}\right){}^{\rm pump}\langle g|{}^{\rm
sys}\langle a|$ (85)

where $|g\rangle^{\rm pump}$ is the ground state of the pump. In the remainder
of the cycle $\hat{\rho}$ evolves to
$\hat{\rho}^{III}=U_{III}^{\dagger}\rho_{II}U_{III}$ with
$U_{III}=\exp\left({-i(H^{\rm sys}+H^{\rm pump}+H^{\rm
couple})\tau_{3}}\right)$ where $\tau_{3}$ is the time the system evolves for
in step III.

In the active cycle, energy is pumped out of the combined system (system plus
pump) and into the environment by the act of initializing the pump. In the
process the entropy of the pump drops and thus the entropy of the environment
increases. For that reason we label this approach an active cooling cycle. One
might quibble that this is a bit of a misnomer since the system is not
thermally equilibrated and thus, it is not clear that the energy pumped out
can be accurately described as heat. But the essence of this active cycle is
very much the same as the cooling in a quantum refrigerator. Moreover as
described below one can use a variation on this approach to achieve a good
approximation to a true thermal equilibrium.

Clearly, this active cycle approach is rather general and variations on this
theme can be developed. As formulated, the approach requires explicit choices
for the size of the pump as well as for the form and strength of $H^{\rm
couple}$, and $H^{\rm pump}$, and $\tau_{3}$. It is clear that to get high
performance with this method one must choose these well. It is an open
question as to what optimal choices are for these.

One obvious approach is to tailor the overall strength of $H^{\rm couple}$ to
the iteration. There is a trade off: strong coupling leads to rapid
transmission of energy from the system to the pump and reduce the overall time
for reducing the energy. However, this comes at cost; large coupling limits
the lowest energy density of the system one can achieve. The cycle can only be
shown to remove energy from the system when the energy in $H^{\rm couple}$ is
negligible. Thus a sensible approach would be to make $H^{\rm couple}$ large
during early cycles in order to facilitate rapid energy transfer and in later
cycles reduce it to allow reduction to lower energy densities.

A similar approach might be taken with regard to the pump. There is a freedom
to set the overall energy scale of $H^{\rm pump}$. It is straightforward to
see that $N^{\rm cycles}$, the number of cycles needed to go from $\langle
H^{\rm sys}\rangle=E_{i}$ (presumably with energy density at the lattice
scale) to a final configuration with $\langle H^{\rm sys}\rangle=E_{f}$ scales
logarithmically with the ratio of $E_{i}$ to $E_{f}$:

$N^{\rm cycles}={\rm A}\log\left(\frac{E_{i}}{E_{f}}\right)$ (86)

where $A$ is a numerical coefficient of order unity, provided that one tailors
the value of the overall strength of the cycle to the cycle in an appropriate
way.

To see how this comes about, consider the trade offs involved in setting the
scale of $H^{\rm pump}$. If it is set too large then it is difficult to induce
transition in the pump and this energy flow will be very slow. On the other
hand if it is too small the maximum amount of energy that can be absorbed in a
cycle is limited. This is clearly true since the system is finite. Moreover,
at some point the energy flow from the system of interest to the pump becomes
negligibly small (either because the system and pump are equilibrating towards
zero net flow). The amount of energy transferred before the energy flow
becomes negligible will clearly depends sensitively on $H^{\rm pump}$.

A simple compromise would be to choose the overall strength of $H^{\rm comp}$
to be large enough so that some modest fixed fraction, $f$ of the system
energy at the beginning of the cycle is transferred before the energy flow
becomes negligible. However, the exact value depends on the initial
configuration of system with the strength increasing with $\langle H^{\rm
sys}\rangle$. Thus one might change the strength of each cycle to keep $f$
approximately the same in each cycle. It would be natural to end each cycle
well before the fraction of the energy, $f$, is transferred, since to the
extent the system equilbrates the energy transfer slows down as the fraction
approaches $f$. For simplicity, let us assume that each cycle stops when the
fraction of energy transferred is $f_{c}\times f$ with $1>f_{c}>0$. One might,
for example, take $f_{c}=\frac{1}{2}$. Thus each cycle reduces by a factor of
$1-f_{c}f$ and the number of cycles it takes to go from a configuration with
$H^{\rm sys}=E_{\rm i}$ to one with $\langle H^{\rm sys}\rangle=E_{\rm f}$ is
thus given by Eq. (86) with

$A=\frac{1}{\log(1-ff_{c})}\;.$ (87)

Of course, the method outlined above is undoubted not the optimal way to
reduce the energy of system of interest from $E_{i}$ to $E_{f}$ given various
resource constraints. But it clearly demonstrates that the minimum number of
cycles needed can be quite modest since the optimal choice will scale no worse
than logarithmically in $E_{i}/E_{f}$.

As given so far this approach can produce density matrices with $\langle
H^{\rm sys}\rangle$ at energy densities of interest and this should allow for
a microcanonical extraction of transport coefficients. While this is
sufficient for our purposes it is useful to note that a variation on this
method is likely to produce density matrices $\rho^{\rm sys}$ which, to good
approximation, are thermal.

The basic approach is to start from preceding approach and produce $\rho^{\rm
sys}$ with $\langle H^{\rm sys}\rangle$ that is at the correct general scale
for the energy of the thermal ensemble at temperature $T$ that we wish to
study. This can be done in comparatively few cycles.

Next one continues to cycle but instead of reinitializing the pump to its
ground state, each time reinitialize the pump to some well-defined state with=
$E^{\rm pump}$. This can be done by initializing to the ground state and then
making unitary transformations to bring the pump to that state. By hypothesis,
the dynamics of the pump is simple enough to do so. The value of $E_{\rm
pump}$ to be match the expectation value of $H^{\rm pump}$ in thermal
equilibrium with a heat bath at temperature $T$.

$\begin{split}E^{\rm pump}&=-T^{2}\frac{Z^{\prime}(T)}{Z(T)}\;\;{\rm with}\\\
Z(T)&=\sum_{i}e^{-iE_{i}/T}\end{split}$ (88)

If in each cycle, step III is allowed to last long enough for the system and
the pump to equilibrate, then after multiple cycles one would expect that
$\rho^{\rm sys}$ to become approximately thermal with temperature $T$.

###  IV.4 Adiabatic state preparation

An entirely different approach to preparing a low-temperature thermal state
begins from the adiabatic theorem [105]. Adiabatic state preparation is a
method for obtaining the _ground_ state of a Hamiltonian. Given the ground
state of the physical Hamiltonian, it is an easy matter to add energy.
Allowing the resulting system to thermalize, we can obtain a thermal state at
any desired energy density. A great benefit to adiabatic methods is that they
require no ancillary qubits to perform the preparation. The primary difficulty
is in obtaining the ground state.

Adiabatic state preparation begins with a Hamiltonian whose ground state is
known, and can be prepared directly on the quantum computer. One allows
quantum time-evolution to progress while slowly changing the Hamiltonian from
the initial one to the Hamiltonian of interest. In principle, such a scheme is
guaranteed to yield the ground state of the system of interest to good
approximation provided that the Hamiltonian is varied slowly. For gauge
theories there is an additional constraint. One wishes to find the ground
state in the physical subspace. As a result, it is natural to restrict all
Hamiltonians in the path from the initial one to the final one to be invariant
under spatial gauge transformations.

There is a generic, practical issue for any approach of this sort: it requires
a slow evolution and hence long times. These constraints are polynomial in all
parameters, but nevertheless may result in significant hurdles for the
forseeable future. Quantifying the time evolution needed is model-dependent
and remains for future work.

An obvious state to start with is the weak-coupling limit of the lattice gauge
theory. In the case of pure Yang-Mills theory, the ground state is the gauge-
invariant projection of the ground state of $N^{2}-1$ free vector fields. The
preparation of this state is described in [106]. However, several difficulties
may prevent this limit from being practical. In the weak coupling limit, the
spatial volume being simulated is much smaller than the confinement scale. In
order to reach the physical regime, this boundary must be crossed, and it is
likely that a phase transition (associated with an exponentially small or even
vanishing gap) is encountered. Even if there is no phase transition along the
path to the desired coupling, the weak-coupling limit is inaccessible to many
proposed truncation schemes for the gluon fields, including approximation by a
finite subgroup [107, 108] and momentum-space truncation [109].

The ground state of the strong-coupling limit is simpler to prepare [106] and
is natural in both the subgroup approximation and the character expansion.
Furthermore, the strong-coupling limit is connected to realistic physical
couplings without a phase transition in the infinite volume limit; therefore,
the desired coupling can be reached without encountering a small gap, even one
polynomially small in lattice units.

Other possible adiabatic paths may exist. One tempting possibility is to add a
Higgs field, allowing a gap to be created in the weak coupling limit and
avoiding the confinement transition.

##  V Discussion

The simulation of heavy ion collisions requires large volumes $\sim
10\;\mathrm{fm}$ with fine resolutions $\sim 0.1\;\mathrm{fm}$ evolved for
long times $\sim 10\;\mathrm{fm}$ — corresponding to enormous quantum resource
requirements. In contrast, the transport coefficients can be extracted from
the hydrodynamic limit of thermal states. Thus, such observables allow for
simplified state preparation and fewer resources. In this work, we propose an
algorithm for extracting the transport coefficients in lattice gauge theory
from thermal states. In order to implement this, it was necessary to derive a
lattice version of the stress-energy tensor in the Hamiltonian formalism.
Further, this algorithm requires the preparation of a thermal state. Here,
multiple viable methods have been suggested. Future work should investigate
the relative merits of each state preparation method in specific theories.
Naive estimates would suggest that with $\sim 10^{4}$ qubits, one could
produce a rough calculation of the viscosity of 3+1d pure glue $SU(3)$ [107,
108]. Such a calculation could account for all systematic errors from lattice
computations. In the near-term, 2+1d $\mathbb{Z}_{N}$ [110, 70] and $D_{4}$
[95, 69] could be computed with a more modest $\sim 10^{2}$ and $\sim 10^{3}$
qubits respectively.

###### Acknowledgements.

We are indebted to Guy Moore and Paul Romatschke for many insightful comments
on hydrodynamics. T.D.C. and Y.Y. are supported by the U.S. Department of
Energy under Contract No. DE-FG02-93ER-40762. Y.Y. is additionally supported
by the Jefferson Science Associates 2020-2021 graduate fellowship program.
H.L. is supported by a Department of Energy QuantiSED grant. Fermilab is
operated by Fermi Research Alliance, LLC under contract number DE-
AC02-07CH11359 with the United States Department of Energy. S.L. is supported
by the U.S. Department of Energy under Contract No. DE-SC0017905.

## References

  * Feynman [1982] R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
  * Carena  _et al._ [2020] M. Carena, H. Lamm, S. Lawrence, Y.-Y. Li, J. D. Lykken, L.-T. Wang, and Y. Yamauchi, Snowmass 2021 LOI TF10-077 (2020).
  * Gelis [2021] F. Gelis, (2021), arXiv:2102.07604 [hep-ph] .
  * Romatschke and Romatschke [2019] P. Romatschke and U. Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Shen and Yan [2020] C. Shen and L. Yan, Nucl. Sci. Tech. 31, 122 (2020), arXiv:2010.12377 [nucl-th] .
  * Stephanov  _et al._ [1998] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998), arXiv:hep-ph/9806219 .
  * Stephanov  _et al._ [1999] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999), arXiv:hep-ph/9903292 .
  * Auvinen and Petersen [2013] J. Auvinen and H. Petersen, Phys. Rev. C 88, 064908 (2013), arXiv:1310.1764 [nucl-th] .
  * Ryu  _et al._ [2015] S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * Alba  _et al._ [2018] P. Alba, V. Mantovani Sarti, J. Noronha, J. Noronha-Hostler, P. Parotto, I. Portillo Vazquez, and C. Ratti, Phys. Rev. C 98, 034909 (2018), arXiv:1711.05207 [nucl-th] .
  * Denicol  _et al._ [2018] G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 98, 034916 (2018), arXiv:1804.10557 [nucl-th] .
  * Dore  _et al._ [2020] T. Dore, J. Noronha-Hostler, and E. McLaughlin, Phys. Rev. D 102, 074017 (2020), arXiv:2007.15083 [nucl-th] .
  * Nunes da Silva  _et al._ [2020] T. Nunes da Silva, D. Chinellato, M. Hippert, W. Serenone, J. Takahashi, G. S. Denicol, M. Luzum, and J. Noronha, (2020), arXiv:2006.02324 [nucl-th] .
  * Reichert  _et al._ [2020] T. Reichert, G. Inghirami, and M. Bleicher, (2020), arXiv:2011.04546 [nucl-th] .
  * Parotto  _et al._ [2020] P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, and M. Stephanov, Phys. Rev. C 101, 034901 (2020), arXiv:1805.05249 [hep-ph] .
  * Cheng and Shen [2021] C. Cheng and C. Shen, (2021), arXiv:2103.09848 [nucl-th] .
  * Rocha  _et al._ [2021] G. S. Rocha, G. S. Denicol, and J. Noronha, (2021), arXiv:2103.07489 [nucl-th] .
  * Plumberg  _et al._ [2021] C. Plumberg, D. Almaalol, T. Dore, J. Noronha, and J. Noronha-Hostler, (2021), arXiv:2103.15889 [nucl-th] .
  * McLaughlin  _et al._ [2021] E. McLaughlin, J. Rose, T. Dore, P. Parotto, C. Ratti, and J. Noronha-Hostler, (2021), arXiv:2103.02090 [nucl-th] .
  * Son and Stephanov [2004] D. T. Son and M. A. Stephanov, Phys. Rev. D 70, 056001 (2004).
  * Mitra and Chandra [2016] S. Mitra and V. Chandra, Phys. Rev. D 94, 034025 (2016).
  * Rougemont  _et al._ [2017] R. Rougemont, R. Critelli, J. Noronha-Hostler, J. Noronha, and C. Ratti, Phys. Rev. D 96, 014032 (2017).
  * Cassing  _et al._ [2013] W. Cassing, O. Linnyk, T. Steinert, and V. Ozvenchuk, Phys. Rev. Lett. 110, 182301 (2013).
  * Noronha-Hostler  _et al._ [2012] J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. C 86, 024913 (2012).
  * Vujanovic  _et al._ [2018] G. Vujanovic, G. S. Denicol, M. Luzum, S. Jeon, and C. Gale, Phys. Rev. C 98, 014902 (2018).
  * Denicol  _et al._ [2016] G. Denicol, A. Monnai, and B. Schenke, Phys. Rev. Lett. 116, 212301 (2016).
  * Niemi  _et al._ [2016] H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016).
  * Paquet  _et al._ [2017] J.-F. Paquet, C. Shen, G. Denicol, S. Jeon, and C. Gale, Nucl. Phys. A 967, 429 (2017).
  * Okamoto and Nonaka [2018] K. Okamoto and C. Nonaka, Phys. Rev. C 98, 054906 (2018).
  * Vujanovic  _et al._ [2020] G. Vujanovic, J.-F. m. c. Paquet, C. Shen, G. S. Denicol, S. Jeon, C. Gale, and U. Heinz, Phys. Rev. C 101, 044904 (2020).
  * de Forcrand [2009] P. de Forcrand, PoS LAT2009, 010 (2009), arXiv:1005.0539 [hep-lat] .
  * Tripolt  _et al._ [2019] R.-A. Tripolt, P. Gubler, M. Ulybyshev, and L. Von Smekal, Comput. Phys. Commun. 237, 129 (2019), arXiv:1801.10348 [hep-ph] .
  * Karpie  _et al._ [2019] J. Karpie, K. Orginos, A. Rothkopf, and S. Zafeiropoulos, JHEP 04, 057, arXiv:1901.05408 [hep-lat] .
  * Cichy  _et al._ [2019] K. Cichy, L. Del Debbio, and T. Giani, (2019), arXiv:1907.06037 [hep-ph] .
  * Sufian  _et al._ [2019] R. S. Sufian, J. Karpie, C. Egerer, K. Orginos, J.-W. Qiu, and D. G. Richards, Phys. Rev. D99, 074507 (2019), arXiv:1901.03921 [hep-lat] .
  * Alexandrou  _et al._ [2020] C. Alexandrou, K. Cichy, M. Constantinou, J. R. Green, K. Hadjiyiannakou, K. Jansen, F. Manigrasso, A. Scapellato, and F. Steffens, (2020), arXiv:2011.00964 [hep-lat] .
  * Zhang  _et al._ [2020] K. Zhang, Y.-Y. Li, Y.-K. Huo, P. Sun, and Y.-B. Yang, (2020), arXiv:2012.05448 [hep-lat] .
  * Huo  _et al._ [2021] Y.-K. Huo _et al._ (Lattice Parton), (2021), arXiv:2103.02965 [hep-lat] .
  * Bazavov  _et al._ [2017] A. Bazavov _et al._ , Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
  * Borsanyi  _et al._ [2020] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 125, 052001 (2020), arXiv:2002.02821 [hep-lat] .
  * Giordano  _et al._ [2020] M. Giordano, K. Kapas, S. D. Katz, D. Nogradi, and A. Pasztor, JHEP 05, 088, arXiv:2004.10800 [hep-lat] .
  * Bellwied  _et al._ [2021] R. Bellwied, S. Borsanyi, Z. Fodor, J. N. Guenther, S. D. Katz, P. Parotto, A. Pasztor, D. Pesznyak, C. Ratti, and K. K. Szabo, (2021), arXiv:2102.06625 [hep-lat] .
  * Aarts and Martinez Resco [2002] G. Aarts and J. M. Martinez Resco, JHEP 04, 053, arXiv:hep-ph/0203177 .
  * Meyer [2011] H. B. Meyer, Eur. Phys. J. A 47, 86 (2011), arXiv:1104.3708 [hep-lat] .
  * Kohno  _et al._ [2014] Y. Kohno, M. Asakawa, and M. Kitazawa, Phys. Rev. D 89, 054508 (2014), arXiv:1112.1508 [hep-lat] .
  * Aarts  _et al._ [2015] G. Aarts, C. Allton, A. Amato, P. Giudice, S. Hands, and J.-I. Skullerud, JHEP 02, 186, arXiv:1412.6411 [hep-lat] .
  * Itou and Nagai [2020] E. Itou and Y. Nagai, JHEP 07, 007, arXiv:2004.02426 [hep-lat] .
  * Brambilla  _et al._ [2020] N. Brambilla, V. Leino, P. Petreczky, and A. Vairo, Phys. Rev. D 102, 074503 (2020), arXiv:2007.10078 [hep-lat] .
  * De Jong  _et al._ [2020] W. A. De Jong, M. Metcalf, J. Mulligan, M. Płoskoń, F. Ringer, and X. Yao, (2020), arXiv:2010.03571 [hep-ph] .
  * Jordan  _et al._ [2011] S. P. Jordan, K. S. M. Lee, and J. Preskill, Quant. Inf. Comput. (2011), arXiv:1112.4833 [hep-th] .
  * Jordan  _et al._ [2014] S. P. Jordan, K. S. M. Lee, and J. Preskill, arXiv:1404.7115 [hep-th] (2014).
  * Jordan  _et al._ [2018] S. P. Jordan, H. Krovi, K. S. M. Lee, and J. Preskill, Quantum 2, 44 (2018), arXiv:1703.00454 [quant-ph] .
  * Trainor [2014] T. A. Trainor, Int. J. Mod. Phys. E 23, 1430011 (2014), arXiv:1303.4774 [hep-ph] .
  * Lamm  _et al._ [2020] H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS), Phys. Rev. Res. 2, 013272 (2020), arXiv:1908.10439 [hep-lat] .
  * Sharma [2013] S. Sharma, Adv. High Energy Phys. 2013, 452978 (2013), arXiv:1403.2102 [hep-lat] .
  * Stachel  _et al._ [2014] J. Stachel, A. Andronic, P. Braun-Munzinger, and K. Redlich, J. Phys. Conf. Ser. 509, 012019 (2014), arXiv:1311.4662 [nucl-th] .
  * Schlichting and Teaney [2019] S. Schlichting and D. Teaney, Ann. Rev. Nucl. Part. Sci. 69, 447 (2019), arXiv:1908.02113 [nucl-th] .
  * Abelev  _et al._ [2013a] B. Abelev _et al._ (ALICE), Phys. Rev. C 88, 044910 (2013a), arXiv:1303.0737 [hep-ex] .
  * Abelev  _et al._ [2013b] B. B. Abelev _et al._ (ALICE), Phys. Rev. Lett. 111, 222301 (2013b), arXiv:1307.5530 [nucl-ex] .
  * Abelev  _et al._ [2014] B. B. Abelev _et al._ (ALICE), Phys. Lett. B 728, 216 (2014), [Erratum: Phys.Lett.B 734, 409–410 (2014)], arXiv:1307.5543 [nucl-ex] .
  * Collins  _et al._ [1989] J. C. Collins, D. E. Soper, and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. 5, 1 (1989), arXiv:hep-ph/0409313 .
  * Ma  _et al._ [2014] J. P. Ma, J. X. Wang, and S. Zhao, Phys. Lett. B 737, 103 (2014), arXiv:1405.3373 [hep-ph] .
  * Baumgart  _et al._ [2019] M. Baumgart, O. Erdoğan, I. Z. Rothstein, and V. Vaidya, Phys. Rev. D 100, 096008 (2019), arXiv:1811.04120 [hep-ph] .
  * Bilgin and Boixo [2010] E. Bilgin and S. Boixo, Physical Review Letters 105, 170405 (2010), arXiv:1008.4162 [quant-ph] .
  * Riera  _et al._ [2012] A. Riera, C. Gogolin, and J. Eisert, Phys. Rev. Lett. 108, 080402 (2012).
  * Lamm and Lawrence [2018] H. Lamm and S. Lawrence, Phys. Rev. Lett. 121, 170501 (2018), arXiv:1806.06649 [quant-ph] .
  * Brandão and Kastoryano [2019] F. G. Brandão and M. J. Kastoryano, Communications in Mathematical Physics 365, 1 (2019).
  * Clemente  _et al._ [2020] G. Clemente _et al._ (QuBiPF), Phys. Rev. D 101, 074510 (2020), arXiv:2001.05328 [hep-lat] .
  * Harmalkar  _et al._ [2020] S. Harmalkar, H. Lamm, and S. Lawrence (NuQS), arXiv:2001.11490 [hep-lat] (2020).
  * Gustafson and Lamm [2020] E. J. Gustafson and H. Lamm, arXiv:2011.11677 [hep-lat] (2020).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. Brandao, and G. K.-L. Chan, Nature Physics 16, 205 (2020).
  * Muller [1967] I. Muller, Z. Phys. 198, 329 (1967).
  * Israel [1976] W. Israel, Annals Phys. 100, 310 (1976).
  * Hiscock and Lindblom [1985] W. A. Hiscock and L. Lindblom, Phys. Rev. D 31, 725 (1985).
  * Hess [2002] B. Hess, The Journal of chemical physics 116, 209 (2002).
  * Pousaneh and de Wijn [2020] F. Pousaneh and A. S. de Wijn, Molecular Physics 118, 1622050 (2020).
  * Alder and Wainwright [1970] B. Alder and T. Wainwright, Physical review A 1, 18 (1970).
  * Kovtun and Yaffe [2003] P. Kovtun and L. G. Yaffe, Physical Review D 68, 025007 (2003).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010), arXiv:0810.4729 [cond-mat.mes-hall] .
  * Zwanzig [1965] R. Zwanzig, Annual Review of Physical Chemistry 16, 67 (1965).
  * Fernandez  _et al._ [2004] G. Fernandez, J. Vrabec, and H. Hasse, Fluid phase equilibria 221, 157 (2004).
  * Jaeger  _et al._ [2018] F. Jaeger, O. K. Matar, and E. A. Müller, The Journal of chemical physics 148, 174504 (2018).
  * Creutz [1977] M. Creutz, Phys. Rev. D15, 1128 (1977), [,132(1976)].
  * Caracciolo  _et al._ [1992] S. Caracciolo, P. Menotti, and A. Pelissetto, Nucl. Phys. B 375, 195 (1992).
  * Caracciolo  _et al._ [1990] S. Caracciolo, G. Curci, P. Menotti, and A. Pelissetto, Annals Phys. 197, 119 (1990).
  * Giusti and Pepe [2015] L. Giusti and M. Pepe, Phys. Rev. D 91, 114504 (2015), arXiv:1503.07042 [hep-lat] .
  * Suzuki [2013] H. Suzuki, PTEP 2013, 083B03 (2013), [Erratum: PTEP 2015, 079201 (2015)], arXiv:1304.0533 [hep-lat] .
  * Del Debbio  _et al._ [2013] L. Del Debbio, A. Patella, and A. Rago, JHEP 11, 212, arXiv:1306.1173 [hep-th] .
  * Bilson-Thompson  _et al._ [2003] S. O. Bilson-Thompson, D. B. Leinweber, and A. G. Williams, Annals Phys. 304, 1 (2003), arXiv:hep-lat/0203008 .
  * Luo  _et al._ [1999] X.-Q. Luo, S.-H. Guo, H. Kroger, and D. Schutte, Phys. Rev. D 59, 034503 (1999), arXiv:hep-lat/9804029 .
  * Carlsson and McKellar [2001] J. Carlsson and B. H. J. McKellar, Phys. Rev. D 64, 094503 (2001), arXiv:hep-lat/0105018 .
  * Luo  _et al._ [1994] X.-Q. Luo, Q.-Z. Chen, G.-c. Xu, and J.-Q. Jiang, Phys. Rev. D 50, 501 (1994).
  * Kogut and Susskind [1975] J. Kogut and L. Susskind, Phys. Rev. D 11, 395 (1975).
  * Kanwar and Wagman [2021] G. Kanwar and M. L. Wagman, (2021), arXiv:2103.02602 [hep-lat] .
  * Lamm  _et al._ [2019] H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS), Phys. Rev. D 100, 034518 (2019), arXiv:1903.08807 [hep-lat] .
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Physical review letters 94, 170201 (2005).
  * Kurkela  _et al._ [2020] A. Kurkela, W. van der Schee, U. A. Wiedemann, and B. Wu, Phys. Rev. Lett. 124, 102301 (2020).
  * Halimeh  _et al._ [2020] J. C. Halimeh, H. Lang, J. Mildenberger, Z. Jiang, and P. Hauke, arXiv:2007.00668 [quant-ph] (2020).
  * Tran  _et al._ [2021] M. C. Tran, Y. Su, D. Carney, and J. M. Taylor, Quantum 2, 010323 (2021), arXiv:2006.16248 [quant-ph] .
  * Sugiura and Shimizu [2013] S. Sugiura and A. Shimizu, Phys. Rev. Lett. 111, 010401 (2013), arXiv:1302.3138 [cond-mat.stat-mech] .
  * Pathria and Beale [2011] R. Pathria and P. Beale,  _Statistical Mechanics_ (Academic Press, Boston, 2011).
  * Roßnagel  _et al._ [2016] J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, Science 352, 325 (2016).
  * Kosloff and Levy [2014] R. Kosloff and A. Levy, Annual Review of Physical Chemistry 65, 365 (2014).
  * [104] D. P. DiVincenzo, Fortschritte der Physik 48, 771.
  * Messiah [1969] A. Messiah,  _Quantum mechanics, volume II_ (North-Holland Publishing Company, 1969).
  * Lawrence [2020] S. Lawrence,  _Sign Problems in Quantum Field Theory: Classical and Quantum Approaches_ , Ph.D. thesis, Maryland U. (2020), arXiv:2006.03683 [hep-lat] .
  * Alexandru  _et al._ [2019] A. Alexandru, P. F. Bedaque, S. Harmalkar, H. Lamm, S. Lawrence, and N. C. Warrington (NuQS), Phys. Rev. D 100, 114501 (2019), arXiv:1906.11213 [hep-lat] .
  * Ji  _et al._ [2020] Y. Ji, H. Lamm, and S. Zhu (NuQS), Phys. Rev. D 102, 114513 (2020), arXiv:2005.14221 [hep-lat] .
  * Ciavarella  _et al._ [2021] A. Ciavarella, N. Klco, and M. J. Savage, (2021), arXiv:2101.10227 [quant-ph] .
  * Yamamoto [2021] A. Yamamoto, PTEP 2021, 013B06 (2021), arXiv:2008.11395 [hep-lat] .