Quantum algorithms for transport coefficients in gauge theories

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.


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.
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 ∼ 10 fm. To simulate the internal dynamics of each nucleus, a resolution smaller than the nuclear radius is required, thus the lattice spacing should be ∼ 0.1 fm. Together, these two scales imply 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 ∼ 10 fm [53]. In order to keep the trotterization errors small, one must use a small time step, again ∼ 0.1 fm. This suggests ∼ 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 ∼ 1 fm box. This reduces qubit costs by a factor of ∼ 10 3 .
Heavy-ion collisions allow for a similar decomposition. In particular, the transport coefficients of the QGP should require L ∼ 1 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 ∼ 200 MeV matter) of a ∼ 0.1 fm. Further, theoretical models [56,57] and experimental data [58][59][60] suggest that thermalization happens rapidly -on the scale of 1 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 wellknown 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 ∼ 200 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 µν (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 ϕ, for instance, is governed by Fick's law: 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 η and ζ, which are defined by their role in the Navier-Stokes equations: where the mass density ρ 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 for the mass density ρ, 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 η 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 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 F and small wavenumbers k, the equilibrium "displacement" is 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 η. To see this, consider the linearization of Eq. (2) about the equilibrium solution ρ = ρ 0 , u i = 0. The mode U cos(kx) obeys Thus a shear wave u 1 ( x) ∝ cos(kx 2 ) is seen to decay as e − η ρ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 µν 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 long 1 times: In the context of a classical molecular dynamics simulation, the calculation of C(t) is straightforward: dx sin(kx)T 01 (x) is measured at various times, and the autocorrelations appear in the resulting stochastic timeseries. At first blush, the appropriate procedure on a quantum computer is rather different. In particular, because measurements cannot be performed without altering 1 We 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 η can.
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 state 2 . 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 ω, k → 0 limit of the two-point correlator of T 12 : 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 η 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 ζ typically proceed from the Green-Kubo relation, which is in terms of the diagonal components of the T µν [82]:

III. ENERGY-MOMENTUM TENSOR
Transport coefficients in quantum field theory are defined in terms of the n-point correlators of T µν . Thus, in order to extract them from lattice calculations, one must define a lattice T µν . 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 µν 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 I. They are constructed from link operatorsÛ and their conjugate momentaπ, defined in Sec. III.2 below. The plaquetteP and cloverĈ are defined in Sec. III.1.

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 lowdimension 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 where the mostly minus convention is used. Following the notation of [83], we normalize the group generators λ a to and functions such as F µν are defined with this normalization. Before we discretize Eq. (8), we must establish some lattice notation. Links are denoted as U n,µ where n is the site the link starts at, and µ 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: On a spacetime lattice, the terms in the EMT are polynomials of plaquettes. Note that, as terms in Eq. (8) contain products of F µν 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 where g s is the coupling constant. Here a µ are the lattice spacings in the µ 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 → 0 before performing continuum extrapolations. If instead we considered Tr F µν (n + 1/2[μ +ν]) 2 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 µν plane: In cases of F 0i , the O(a 0 ) error induced in the time direction is acceptable as we take the limit a 0 a in the Hamiltonian formalism. This implies that we don't need to average over four plaquettes in both directiont andî. Instead, we need to average over only two plaquettes inî direction around the site: The second term of Eq. (8), F µα F να , requires us to construct a discretization for F µν itself. The naive discretization for the field strength tensor is Note again that F N µν (n) is evaluated on a lattice site. The RHS of Eq. (14) approximates the value at the center Tr πn,jImĈij(n) + Tr Û † n−ĵ,jπ n−ĵ,jÛn−ĵ,j ImĈij(n) 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: From these clover operators we construct the improved clover discretization of the field strength: Here we have used P −µ,ν (n), for instance, to denote the plaquette in which the first link begins at site n and ends at site n −μ. In an Abelian theory, we of course have P −µ,ν (n) = P µ,ν (n − µ) † , 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 2 0 ) regardless of orientation. However, as mentioned, working in the Hamiltonian formalism implies a 0 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. (13) -averaged over two plaquettes as in Fig. 1: This is enough to implement T µν correctly up to O(a 2 , a 0 ):

III.2. Transfer matrix
With a choice of discretization of F µν and F 2 µν , one can compute a lattice T µν in the action formulation. In the Hamiltonian formulation, T µν (which are functions of the space of field configurations) must be replaced by operatorsT µν 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 → 0 limit) yields a Hamiltonian usable for quantum simulations. Now, consider a perturbation to that action by a term proportional to O. The free energy, differentiated with respect to the perturbation, yields the O . Connecting the perturbed system to one in the Hamiltonian formalism via the transfer matrix, we obtain a perturbed Hamiltonian, revealing what operatorÔ 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 action 3 : 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 piecesT = e −iHa0 (this is the transfer matrix) and inserting a complete set of states: Here |U 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 matrixT = e −ia0Ĥ should be given by 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Û that are diagonal in the position basis, and their conjugate operatorsπ defined below. As spatial plaquettes are already written only with U s, they remain in the same form and all arguments U become operatorsÛ . 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, |g U m,j when m = n and i = j |U m,j otherwise (24) which can be written with Hermitian operatorsπ aŝ In short,π generate rotationsR of each link. Here, in the case of SU (N ), the x a are N 2 − 1 real parameters parameterizing the group element g ∈ SU (N ), andπ 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 Exchanging g for real parameters x a , one obtainŝ dx a n,i e ix a n,iπ a n,i In the limit a 0 → 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 ofT is gauge invariant, we take the saddle-point approximation around x = 0 to obtain For brevity we have written combined three indices into ρ ≡ (n, i, a). Together with the spatial plaquette terms from V , the Hamiltonian iŝ ReTr P ij (n) . (29) The first term in the Hamiltonian can also be written as Tr[π 2 n,i ] by defining the N × N matrixπ n,i =π a n,i λ a . Under gauge transformations, we havê and this implicitly defines the transformation law for the operatorsπ 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 O -in the sections below this perturbation will be T µν . In the action formulation, one way to measure the expectation value of O at time t 0 in a system governed by the action S is to perturb the action by the function O(U ) and define Differentiating Z -or rather its logarithm, the free energy -with respect to yields: A perturbed Hamiltonian H , yet to be determined, corresponds to this perturbed path integral. The corresponding transfer matrix has matrix elements given by Differentiating the transfer matrix with respect to yields i times the desired operator; therefore we see the Hamiltonian should be perturbed tô So, by perturbing the Lagrangian with a component of T µν and finding the corresponding perturbed Hamiltonian H , one can read off the operatorsT µν . To be precise, if the perturbing parameter is , then the corresponding operator in the Hamiltonian formalism is the coefficient of ( ).

III.3. Tµµ 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 To simplify notation, S , H , and T 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 to give the same dynamics as by the action S , they should be connected via transfer matrix T = e −ia0H as with 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 where we have introduced shorthand Performing the saddle-point approximation gives: with Here we have abbreviated ρ ≡ (n, i, a) and σ ≡ (m, j, b).
Finally, we can construct T µµ from these operators: 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) givesT µµ up to O(a). To improve theT µµ operators up to O(a 2 ), we use Eq. (12) and Eq. (13) and take the average around the site n 0 : These operators enable us to constructT µµ up to discretization errors that are O(a 2 , a 0 ).

III.4. Tij in the Hamiltonian formulation
Let us now move to deriving the operatorsT ij , that is, the off-diagonal spatial parts of the EMT: 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 thatT is given by Eq. (36) with: 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π. UsingR(g) operators,T can be written (53) Evaluating the integral via the saddle point x = 0 (exact in the limit a 0 → 0) gives: (54) which at O( ) yieldsĤ : Tr[π n0,1πn0,2 ] − a 3 Tr F N ik (n 0 )F N jk (n 0 ) (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 Thus the naiveT ij (n 0 ) in the Hamiltonian formulation iŝ The clover approximations are obtained from F C ij in Eq. (16) and F B i0 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, We use the definitions of Fig. 2 for the links around n 0 . For the example of U n0,1 , we denote operators and functions on them as U 1 ,Û 1 ,π 1 , and g 1 = e ix a 1 λ 2 . ThenT is After the saddle-point approximation around x = 0,T simplifies and becomeŝ Matrix elements of (M 1 ) ρσ are zero other than where sites n 1 , n 2 are as labeled in Fig. 2. Now by expanding Eq. (60) to linear order in , we obtain the perturbed Hamiltonian Tr [π 1π3 ] + Tr π 1Û † 2π 2Û2 Thus F i0 F j0 is generally implemented as 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 µν are fully constructed following Eq. (50).

III.5. T0i 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: 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 from which we find the followingT : Tr[(gn 0 ,2−g † n 0 ,2 )(P12(n0)−P † 12 (n0))]+iV .
From the perturbed Hamiltonian, we read off Therefore, the operator to give T 0i with naive discretization isT The naive discretization induces O(a) error in T 0i . To improve up to O(a 2 ), we use the clovers, F C ij in Eq. (16) and F B i0 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 In the following we introduce notation for links around the site n 0 as in Fig.3. The transfer matrix is then whereĈ is the clover operator. In the limit of a 0 → 0, one can evaluate the integral via the saddle point around which yields to O( ), Tr π 1 ImĈ 12 (n 0 ) FromĤ , we read off, for general F 0j F ij (n 0 ), 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 welldescribed by hydrodynamics. The equilibration time is believed to be ∼ 1 fm; comparable to the characteristic scale of QCD. It is highly plausible that such rapid equilibration is a generic feature of strongly-coupled theories [12,97]. Accordingly, one expects that if an out-ofequilibrium 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 ∼ 200 MeV. Very cold temperatures, near the true ground state, may be difficult to prepare due to frustration. Very high temperature states (T Λ QCD ) likely have difficulty thermalizing, as the fluid is closer to a free gas. At T ∼ 200 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 timeevolution (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 O = Z −1 Tr e −βH 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 where a subscript of sys denotes a state of the physical system being studied, and one of 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 so that for an operator O 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 descriptions 4 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 |Ψ of the system Hamiltonian H sys sampled from the set of all eigenstates with energy density is near , expectation values Ψ|O|Ψ 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 , 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 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 couple , which is initially switched off. The coupling Hamiltonian satisfies 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 couple must be small in the sense that at all stages in the evolution the absolute value of its expectation value: where H sys 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 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 formally requires a heat bath at least as large as the system; otherwise the condition HB ψ k |ψ k HB = δ k ,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 HB ψ k |ψ k HB = δ k ,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 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 H sys = 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 At that stage, the system is decoupled from the first heat bath and coupled to the second. We perform more time evolution such that H sys = E 2 , with E 2 < E 1 . Continuing the process, we see that the energy is monotonically decreasing: 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 pump = |0, 0, 0, · · · 0, 0, 0 and then to choose H pump so its ground state is the fiducual state. One can reinitialize by simply measuring the qubits using σ z to determine if the qubit is in |1 or |0 and if the result is |1 act on it with the unitary operator σ 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 latticescale 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 sys that acts on the system, H pump , that acts in the pump space and H couple that couples the two. The coupling dynamics is constructed to be gauge invariant, H couple must have non-zero commutators with H sys and H 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 H sys = E I relative to its ground state, the pump subspace is then initialized to its ground state, with H pump = 0. Since H couple satisfies inequality (81), the total energy (relative to the ground state) of the combined system H sys +H pump +H couple , 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 H sys + H pump + H couple remains at E I . In step III the combined system undergoes time evolution under H sys +H pump +H 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, H sys reaches some value E III where E III < E I . The time evolution then is switched off suddenly leaving H sys = E III < E I . At that point, the cycle repeats. Thus each time through the cycle, H sys drops. One could continue iterating the cycle and lowering H sys until the inequality (81) ceases to hold.
Consider the reduced density matrix of Eq. (79) for the system as after tracing over the pump. By construction step I preserves ρ 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,ρ I is given bŷ pump b| sys a| (83) where the states |a sys and |b 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 bŷ Thus the full density matrix at step II is given bŷ where |g pump is the ground state of the pump. In the remainder of the cycleρ evolves toρ III = U † III ρ II U III with U III = exp −i(H sys + H pump + H couple )τ 3 where τ 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 couple , and H pump , and τ 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 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 couple is negligible. Thus a sensible approach would be to make H 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 pump . It is straightforward to see that N cycles , the number of cycles needed to go from H sys = E i (presumably with energy density at the lattice scale) to a final configuration with H sys = E f scales logarithmically with the ratio of E i to E f : 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 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 pump .
A simple compromise would be to choose the overall strength of H 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 H sys . 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 × f with 1 > f c > 0. One might, for example, take f c = 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 sys = E i to one with H sys = E f is thus given by Eq. (86) with . (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 H sys 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 ρ sys which, to good approximation, are thermal.
The basic approach is to start from preceding approach and produce ρ sys with H sys 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 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 pump to be match the expectation value of H pump in thermal equilibrium with a heat bath at temperature T .
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 ρ sys to become approximately thermal with temperature T .

IV.4. Adiabatic state preparation
An entirely different approach to preparing a lowtemperature 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 ∼ 10 fm with fine resolutions ∼ 0.1 fm evolved for long times ∼ 10 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 ∼ 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 Z N [70,110] and D 4 [69,95] could be computed with a more modest ∼ 10 2 and ∼ 10 3 qubits respectively.