Classical and Quantum Computing of Shear Viscosity for 2 + 1 D SU(2) Gauge Theory

,


I. INTRODUCTION
The scientific goal of relativistic heavy ion collisions is to study the deconfined phase of nuclear matter at finite temperature and/or density, known as the quarkgluon plasma (QGP).The most striking property of the QGP created in current heavy ion collision experiments is its small shear viscosity, as shown by the good agreement between the experimental data on various particles' yields and azimuthal distributions and a description that is mainly based on relativistic hydrodynamic equations with small shear viscosity [1,2].The smallness of the shear viscosity indicates the QGP created in current heavy ion collision experiments is strongly coupled.Interestingly, the current value of the ratio between the shear viscosity and the entropy density η s extracted from experimental data [3,4] is consistent with the AdS/CFT calculation result for a N = 4 supersymmetric Yang-Mills plasma in the strong coupling limit, which is 1 4π [5].Theoretically, shear viscosity can be calculated from real time two-point correlation functions of stress-energy tensors via the Kubo formula.However, this computation is hard for both perturbative and nonperturbative approaches in QCD [6].Perturbatively, certain diagrams have to be resummed due to the existence of "pinching poles" [7][8][9] and the convergence of the perturbative series is poor when the temperature of the QGP is below 1 GeV [10], which is the temperature range of most interest in current collision experiments.Nonperturbatively, Euclidean lattice QCD methods have been applied to calculate the relevant two-point correlation functions in imaginary time [11][12][13].However, extraction of the shear viscosity from the imaginary time correlation function involves an ill-defined "inverse problem," and thus is not under good theoretical control.Different frequency dependence of the real time correlation function can give the same Euclidean correlation function in imaginary time.These limitations of current theoretical studies urgently demand a new technique for transport coefficients calculations, since a fully theoretical determination of the shear viscosity provides an independent check of the experimental extraction and is thus of great value.These calculations can test the current hydrodynamical framework for heavy ion collisions and deepen our understanding of nuclear matter in conditions that are unachievable in experiments at the moment.
In this paper, we consider the Hamiltonian formulation of lattice gauge theory and investigate the calculation of shear viscosity from the retarded correlation function obtained from real time Hamiltonian evolution.Our study is motivated by recent developments in quantum computing for lattice gauge theories [14][15][16][17][18][19], which follow Feynman's idea [20] to use quantum computers and hopefully will be able to tame the exponential growth of the Hilbert space and perform more efficient quantum simulations than classical devices.A previous work studied the construction of lattice operators for stress-energy tensors and quantum algorithms for thermal state preparation in the same context [21].Here, more specifically, we will carry out detailed calculations for the SU(2) pure gauge theory in 2 + 1 dimensions (2 + 1D) on a small lattice via both classical and quantum methods and analyze various systematics.The paper is organized as follows: In Sec.II, we will briefly review the Kubo formula for the shear viscosity calculation in the context of the SU(2) pure gauge theory in 2 + 1D.Then, in Sec.III, we will explain the lattice Hamiltonian formulation for the calculation, followed by an introduction of a quantum algorithm in Sec.IV.Various systematics of the calculation will be discussed in Sec.V. We will show both classical and quantum results in Sec.VI and draw conclusions in Sec.VII.

II. SHEAR VISCOSITY IN 2 + 1D SU(2) GAUGE THEORY
The Lagrangian density of the continuum 2+1D SU(2) gauge theory can be written as where g is the coupling constant, and ν is the non-Abelian field strength tensor with A a µ as the gauge field.In particular, F a 0i denotes the electric field along the i-th spatial direction and F a ij is related to the non-Abelian magnetic field.The alphabet indices a, b, c ∈ [1, 2, 3] label the SU(2) adjoint indices.
Stress-energy tensors of the theory are given as Standard linear response analysis and gradient expansion of the stress-energy tensor for relativistic hydrodynamics in the Minkowski spacetime with a small metric perturbation h xy (t) lead to [22,23] where G xy r (ω) can be expressed in terms of the retarded Green's function of T xy2 G xy r (ω) = dt e iωt G xy r (t) ≡ dt d 2 x e iωt G xy r (t, x) The density matrix in the definition of G xy r (t, x) describes the thermal state at temperature By utilizing translation invariance, the retarded correlation function can also be written as where A denotes the area of the system.
Combining everything together gives We have two ways to evaluate G xy r (t) as shown in Eqs. ( 4) and (6).If we know all the eigenstates |n⟩ of the system and their corresponding eigenenergies E n , we can write When the system has an infinite number of states, as continuum quantum field theories do, the symbol Rigorously speaking, Eq. ( 3) is the tree-level matching condition between the hydrodynamic effective theory and the full theory.In this sense, η can be thought of as a Wilson coefficient of the hydrodynamic effective theory.The matching condition between η and the retarded Green's function becomes more complicated once nonlinear terms are taken into account in T xy of the hydrodynamic effective theory.These nonlinear terms contribute at one-loop level and generate the so-called long-time tails [25], which lead to a logarithmic divergence in two spatial dimensions [26,27].As a result, η becomes scale dependent; it is some function of η(ω) that equals G xy r (ω) when ω is small {see, e.g., Eq. ( 52) in Ref. [27] as the one-loop matching condition}.In this work, we will only consider the tree-level matching and use Eq. ( 8) to calculate the shear viscosity for simplicity.Classically, Eq. ( 8) can be evaluated on a lattice by solving all eigenenergies and eigenstates.We will also take the continuum limit along the renormalization group flow of the coupling.

A. General Setup
The Kogut-Susskind Hamiltonian [28] of the 2+1D SU(2) gauge theory can be discretized on a honeycomb lattice as shown in Fig. 1 where a in the denominator is the side length of the honeycomb, and we have shifted the energy reference point.The honeycomb plaquette operator is defined as the trace of the product of the six Wilson lines U (x, î) on the edges of one honeycomb.The two-vector x = (i, j) labels the position of a honeycomb on the lattice plane along the directions specified as in Fig. 1.The electric field is projected along three unit directions E a i ≡ êi • E a , where the three unit vectors are defined as in Fig. 1.On each link, only one type of projected electric field lives, i.e., i is 1, 2 or 3.More details can be found in Ref. [29].Physical states satisfy Gauss's law at each vertex for every a.
We use the electric basis that labels each link by the quantum number j.In this basis, the electric energy is diagonal [30][31][32] The matrix element of the plaquette term (magnetic energy) has been worked out to be [29,33,34] (see Refs. [35][36][37][38][39] for the square plaquette case): where {j} ({J}) labels the states on the six links of the honeycomb plaquette before (after) the action of the operator.The product is over all the vertices V of the honeycomb plaquette, attached to which are two internal links labeled by the subscripts a and b and an external link labeled by x.
For the 2+1D SU(2) lattice gauge theory, gauge invariant states can be uniquely represented by the j values if Honeycomb lattice on which the central position of each plaquette is labeled by (i, j) along the two axes shown, starting from i = j = 0.For example, the plaquette marked in blue is located at (i = 1, j = 2).The red dots represent the positions at which the stress-energy tensors are evaluated.The location of a red dot can be identified by the coordinates of the three plaquettes sharing the same red dot.The six plaquettes labeled by the K values around the blue-colored plaquette at (i = 1, j = 2) are used when the magnetic term at (i = 1, j = 2) is written out explicitly, as in Eq. (18).
each vertex has at most three links joined.On a square lattice, one has to introduce extra labels besides the j values in order to represent gauge invariant states, which leads to additional computational cost.This is an advantage of using the honeycomb lattice [29].
From Eqs. ( 2) and ( 10), we find Using the electric field projection, we find x and E a 2 = −E a y .Combining with Gauss's law E a 1 + E a 2 + E a 3 = 0, we can express T xy as In this expression, we need to specify the position where T xy is defined, since the two electric fields E a 1 and E a 3 are defined on different links.We use the convention that the vertex joining the two electric fields represents the position of T xy .On a 3 × 4 lattice as shown in Fig. 1, we should specify 12 positions for different T xy 's.We choose the 12 red points in Fig. 1 as our convention, which can be easily generalized to bigger lattices.Summing over all red points gives where N plaq is the total number of honeycomb plaquettes on the lattice and is equal to the number of red dots, as shown in Fig. 1.Using Eqs. ( 4) and ( 6) leads to B. Truncation at jmax = 1 2 For quantum computation discussed later, we need to decompose the Hamiltonian and T xy in terms of tensor products of Pauli matrices, for which quantum circuits of implementation are known.This decomposition has been done for the case with the local Hilbert space truncated at j max = 1 2 .Under this truncation, the Hamiltonian can be represented as a 2D Ising-like model [29] aH = H el + H mag (18) where Π ± i,j = (1 ± σ z i,j )/2 are the projection operators onto the spin-up and spin-down states that represent the plaquette state at (i, j) as labeled in Fig. 1.H el represents the electric part of the Hamiltonian and H mag stands for the magnetic part.We have multiplied the Hamiltonian by the lattice spacing a such that every quantity is unitless The index K comes from a periodic (K mod 6) chain {K = 0 : (i, j + 1), K = 1 : (i + 1, j), K = 2 : (i + 1, j − 1), K = 3 : (i, j − 1), K = 4 : (i − 1, j), K = 5 : (i − 1, j + 1)}, as shown in Fig. 1.
The xy component of the stress-energy tensor for the honeycomb plaquette located at (i, j) is calculated from Eq. ( 15) where the two electric fields are those attaching the red vertex, which is at the upper right corner of the honeycomb, as shown in Fig. 1 T

C. Closed Boundary Condition
For a finite lattice, we use a closed boundary condition in which all the links outside the lattice boundary are in j = 0 states.In other words, no electric fluxes go out of the lattice.For the case with j max = 1 2 truncation, the imposed closed boundary condition is equivalent to setting all the spins outside the boundary to be pointing down.
We choose the closed boundary condition since it makes the quantum circuit construction more convenient for the j max = 1 2 case.The periodic boundary condition results in an overall spin-flipping degeneracy in the spin representation of physical states [29].Lifting up the degeneracy will distort the expressions of the Hamiltonian and T xy away from the spin representations in Eqs.(18) and (20).Then the corresponding quantum circuits are unknown and thus need explicit constructing, which can be computationally expensive.This concludes our discussion of the calculation setup.Classical computing results of the shear viscosity that are obtained from Eqs. ( 8) and (20) will be shown in Sec.VI.In the next section, we will introduce a quantum algorithm to evaluate the retarded correlation function.
In this section, we present a quantum algorithm to calculate the retarded Green's function G xy r (t).A schematic diagram of the quantum circuit for a four-qubit system (e.g., a 2 × 2 lattice with j max = 1 2 ) is shown in Fig. 2, which consists of four parts: thermal state preparation, unitary transformation driven by T xy for the evaluation of the commutator in G xy r (t), real time evolution and measurements.

A. Thermal State Preparation
We first discuss the quantum circuit for thermal state preparation that uses the algorithm of Refs.[40,41].It is based on imaginary time propagation.The imaginary time propagation techniques are well-known methods for classically preparing ground or thermal states.However, for a quantum algorithm, we have to deal with the nonunitary nature of the imaginary time operator.In the algorithm we will use, this is overcome with a diluted operator using ancilla qubits [40].
The algorithm is schematically shown in Fig. 3.The physical qubits representing a state of the lattice gauge theory are the first, third, fifth and seventh qubits from the top of the figure.We first apply a Hadamard gate to each of them and then apply controlled-NOT (CNOT) gates with these four physical qubits as controls and four auxiliary qubits (second, fourth, sixth and eighth from the top) as targets.After the CNOT gates, the four auxiliary qubits are measured.The measurement outcomes are not needed for the remaining circuit.So effectively, the measurements serve as partial trace and the resulting part implements the algorithm of Ref. [41] to prepare the thermal state, which is shown explicitly in Fig. 3.The e ±i π 4 Σα part is used to calculate the commutator, as will be explained in Sec.IV B.
physical state is a maximally mixed state where n s is the number of qubits in the physical system and equals to four in the example shown in the figure.
Then one applies the quantum imaginary time propagation (QITP th ) to the physical qubits (whose number is four in the figure) plus an additional ancilla qubit, which is the bottom qubit in the figure: where p is a parameter to be tuned (in the limit τ → 0, p is equal to the success probability), and E T is another parameter chosen to be smaller than or equal to the ground state energy.The operators e −τ H and √ 1 − e −2τ H act on the physical system qubits.If p = 1 and τ = β 2 are chosen and the measurement of the ancilla qubit returns |0⟩, one can show the physical system is in the Gibbs state where p s indicates the success probability of measuring the ancilla in |0⟩ state.
The thermal state preparation algorithm used here is efficient for preparing high temperature thermal states.At low temperature, it becomes less efficient as the Hilbert space dimension increases.In particular, the success probability, i.e., the probability of measuring the ancilla qubit in |0⟩, decreases with the Hilbert space dimension at low temperature.The lower bound (that corresponds to T = 0) is determined by the degeneracy of the ground state divided by the dimension of Hilbert space; in the worst scenario, it is given by 1  2 ns , with n s the number of system qubits.For future applications in QCD, we only need temperatures above the confinementdeconfinement crossover temperature, which is around 150 MeV.So the reduced efficiency of the thermal state 3. Quantum circuit for preparing the thermal state in quantum processors.The dashed box shows the gates for initializing the system qubits to be the maximally mixed state.The first, third, fifth and seventh qubits from the top represent the system while the second, fourth, sixth and eighth are used to initialize the maximally mixed state.The bottom qubit is the ancilla qubit used to implement the imaginary time evolution in a unitary way.
preparation algorithm at lower temperature may not be a severe problem.In Sec.V E, we will perform a more detailed analysis of the thermal state preparation success probability as a function of the system size.
In practice, to explicitly construct the quantum circuit for the QITP th from the Hamiltonian in Eq. ( 18), we apply Trotter decomposition at first order Nτ iτ =1 (i,j) where ∆τ = τ Nτ is the imaginary time Trotter step, and (i, j) denotes a plaquette position as in Eq. ( 18).In the Trotter decomposition, we need to add one additional ancilla qubit for each non-commuting Hamiltonian term, unless the quantum hardware allows ancilla qubit reset in the middle of the circuit.More concretely, the QITP th can be written as where σ y anc denotes the Pauli-y matrix acting on the ancilla qubit.The arccos function is best computed in the diagonal basis of H el/mag ij .The electric part is already diagonal in the computational basis.For the magnetic part at (i, j), one needs to apply a Hadamard gate before and another one after the QITP th circuit to the qubit representing the state at (i, j) in order to convert to a basis where H mag ij is diagonal.Once the arccos part is made diagonal, one can decompose it into sums of tensor products of identity and σ z matrices.The procedure to construct quantum circuits for tensor products of Pauli matrices is well known [42].Appendix B shows a quantum circuit for the magnetic Hamiltonian term on a 2 × 2 lattice with j max = 1 2 .As we will show later in Sec.V D, we only need to use a single Trotter step in the imaginary time because its Trotter error is very small, i.e., ∆τ = β 2 .

B. Circuit for Commutator of T xy
Using Eq. ( 20), the commutator [T xy sum (t), ] can be rewritten as the sum of two commutators: where we have decomposed T xy ij (0) operator in terms of two Pauli strings as in Eq. (20).To simplify the notations, we have introduced Σ α as Σ 0 = σ z i,j+1 σ z i+1,j , and Σ 1 = σ z i,j+1 σ z i+1,j .The commutator between a Pauli string A and a (generic) unitary operator B can be easily evaluated in quantum processor using the following relation [43,44] Applying Eq. ( 27) to Eq. ( 26) with A and B identified as Σ α and T xy sum (t), respectively, we find where we have used O(t) = e iHt O(0)e −iHt and omitted the argument t = 0 of T xy sum .

C. Real Time Evolution
As just discussed, the quantum computation procedure requires real time evolution.Here, we briefly discuss its implementation by applying the Trotter decomposition at first order.We divide the total time length t into N t steps with the step size ∆t = t Nt .The full real time propagator becomes a product of the short-time electric and magnetic Hamiltonian evolution operators as the following: where H mag ij (H el ij ) indicates the magnetic (electric) part of the Hamiltonian at the (i, j) plaquette position.Each Hamiltonian piece just contains commuting strings of Pauli matrices, for which the construction of quantum circuits is well known [42].For example, the circuit for the magnetic part is discussed in Ref. [29].A discussion of the Trotter errors in the real time evolution for the calculation of G xy r (t) can be found in Sec.V D. After implementing the real time evolution gates, we measure the quantum state in the computational basis in which the operator T xy sum (0) = k,l T xy kl (0) is diagonal.

D. Post-Measurement Processing
The last step to obtain the retarded correlation function is to perform a post-measurement analysis.We note that the stress-energy tensor operator T xy sum is a sum of Pauli-z tensor products; so it is diagonal in the computational basis.Therefore, the thermal expectation value of the commutator [T xy sum (t), Σ α ] is given by We now summarize the quantum algorithm: After preparing the thermal state, we first apply the gates for e ±i π 4 Σα , as shown in Fig. 2 [+ for the first line and − for the second line in Eq. ( 28)].Then we evolve the resulting state in real time e −iHt .Finally, we perform projective measurements in the computational basis.The measurement results allow us to reconstruct the thermal expectation value of the commutator in Eq. ( 28) by postmeasurement processing and taking the difference of the results obtained from the two different circuits that differ in the sign of the Σ α term.Because we have two different Pauli strings in order to evaluate Eq. ( 26), we have to run four different quantum circuits that differ in the Pauli operators (different α) and the signs of the term i π 4 Σ α .
To conclude this section, we prove that the action of the proposed quantum circuit gives the correct value of G xy r (t).The first set of gates prepares the thermal state density matrix, ρ T = e −βH Z .Then, after applying the gates for e ±i π 4 Σα and the real time evolution gates U t , the density matrix of the quantum processor becomes Measuring the quantum processors in the computational basis, we obtain the diagonal part of the density matrix.
Using the cyclic property of the trace, and the fact that T xy (0) is diagonal in the computational basis, we can write: Using Eq. ( 27), the thermal expectation value of the commutator [T xy sum (t), Σ α ] is obtained from A final usage of Eq. ( 26) leads to the retarded Green's function of T xy .

V. CALCULATION SYSTEMATICS
Before presenting results, we discuss various systematics of the calculation.The analysis we will present in this section is important to understand whether one can obtain the physical quantity (i.e., shear viscosity) at a given accuracy with a given amount of quantum computing resource.This analysis will be useful for large scale quantum computation of transport coefficients in the future.

A. Continuum Limit and Renormalization
For physical limits, one needs to take the continuum a → 0 and the infinite volume limits and remove the local Hilbert space truncation by setting j max → ∞.In the continuum limit a → 0, the coupling of the theory needs proper renormalization.In 2+1D SU(2) pure gauge theory, the mass dimension of g is 0.5, so a unitless quantity for the coupling is ag 2 .What needs to be done is to tune ag 2 as a → 0 such that a physical observable is invariant.An example of physical observables is the correlation length of the ground state, which can be extracted from the subsystem size dependence of the entanglement entropy [45].In a renormalization scheme where both the pressure and the trace of the stress-energy tensor T µ µ do not require any additional renormalization other than the running coupling, the renormalization of the coupling is given by [46] which means the rescaled coupling ag 2 ∝ a.In the Hamiltonian approach, the lattice gauge theory also has a truncation in the local Hilbert space, labeled by j max .
Whether the running coupling has a non-trivial dependence on j max should be studied analytically and tested against numerical calculations, which are left for future work.
Besides the running coupling, the stress-energy tensor component T xy may need additional renormalization since the lattice breaks the Lorentz invariance of the continuum theory.Generally, one can write where R and B represent renormalized and bare quantities.The bare quantity is the direct numerical result obtained by using the truncated lattice gauge theory with a lattice spacing a and local Hilbert space truncation j max .Z(µ, a, j max ) is the additional operator renormalization factor needed when taking the continuum limit, which should be distinguished from the partition function introduced in Eq. ( 5) and depends on the final renormalization scheme in which we want to obtain the quantity η R , e.g., the MS scheme with the renormalization scale µ.One way to obtain Z(µ, a, j max ) is to perform a lattice perturbation theory calculation.Alternatively, one can develop gradient flow methods for the Hamiltonian approach to regularize the stress-energy tensor operator, as done in Euclidean lattice calculations [47,48].Compared with Euclidean lattice calculations, the potential dependence on j max is new and needs systematic understanding.
As can be seen, taking the continuum limit and renormalizing the quantity properly can be complicated.Since our current study is on a small lattice with a low-j max truncation, we will only take into account the coupling renormalization here and leave studies of additional operator renormalization to future work.

B. Local Truncation Effect
We then consider local Hilbert space (i.e., j max ) truncation effect with a given lattice spacing a, which means the coupling ag 2 is fixed, and a fixed lattice size.The truncation effect can lead to artifacts in thermodynamic description of the system, since they constrain the energy spectrum of the system.This has been considered on a SU(2) plaquette chain [49], and here, we generalize it to the honeycomb lattice case.
Specifically, we first consider the internal energy density ε and the entropy density s as a function of j max on a 2 × 2 honeycomb lattice with ag 2 = 1, which are defined as where γ and α s T 2 , respectively, which are expectations from the continuum theory.The fits use the results with j max = 2 and in the temperature range T < 2.5 (in lattice units).The fitted parameters are α ε = 0.0191, γ = −0.0466and α s = 0.0234.We see that the 2 × 2 lattice with j max = 1 already shows continuum behavior for ag 2 = 1 in the temperature range T ≲ 3.In the region 3 < T < 5, the curves of ε and s already converge with j max = 1.5.Further increasing j max is of no help for continuum physics, and one has to use a bigger lattice.
At last, we provide an analytic estimate of the j max needed.At a given lattice size, if we want to describe all the states below a fixed energy E [corresponding to the Hamiltonian in Eq. ( 9)] with an accuracy 1 − ϵ, the minimum j max needed is at most where N l is the number of links on the lattice, and E = FIG. 5. Internal energy and entropy densities in lattice units as functions of temperature for several lattice sizes with jmax = 1 2 and ag 2 = 1.Ni and Nj specify the size of the honeycomb lattice along the i and j axes respectively, as shown in Fig. 1.
9g 2 a 2 N p , where N p denotes the number of plaquettes on the lattice.We provide a proof of this formula in Appendix A.

C. Finite Volume Effect
Next we study the lattice size dependence of the internal energy and entropy densities with ag 2 = 1 and j max = 1  2 .The results are shown in Fig. 5, where it can be seen that these two densities change little with the lattice size.Increasing the lattice size leads to a decrease of ε at high temperature, which indicates that one needs to increase j max to maintain the same energy density at a given temperature.Both energy and entropy densities saturate at high temperature, which is a finite size effect.Figure 6 shows the energy eigenstate density on a 4 × 4 lattice with j max = 1 2 and ag 2 = 1.The total number of states is 2 16 = 65536.When the energy is below 35 (in lattice units), the state density ρ(E) keeps increasing with E. This is qualitatively similar to the continuum theory.However, once the energy exceeds 35, ρ(E) starts to drop with E, which is an artifact originated from the finite size of the Hilbert space.Another way of inspecting the finite size effect is to compare η(t f ) calculated via different methods on different lattices.In Fig. 7a, we compare the result of η(t f ) obtained from the commutator [ T xy (t), T xy (0)] as in the third line of Eq. ( 8), where T xy (0) is located at (i, j) = (1, 1), with that from the commutator [ T xy (t), T xy (0)]/A as in the second-to-last line of Eq. ( 8) on a 4 × 4 lattice with the cutoff j max = 1 2 and ag 2 = 1 at β = 0.3 (in lattice units).We see that η(t f ) oscillates in both results with the [ T xy (t), T xy (0)] case oscillating more severely.The oscillation is caused by the finite state density.Inspecting Eq. ( 8), we can see the oscillating factor f (t f ) only smooths out if the energy levels are dense enough.Thus we expect at lower temperatures that the oscillating is more severe when the Hilbert space is of fixed size.The reason why the calculation with the [ T xy (t), T xy (0)] commutator oscillates more severely is the additional non-positive definite term in the double sums over n, m, i.e., ⟨n| T xy |m⟩⟨m|T xy |n⟩, which can be positive or negative.On the other hand, the |⟨n| T xy |m⟩| 2 term is positive semi-definite so it leads to a smoother result.If we had used periodic boundary conditions, spatial translation invariance would have been preserved and as a result, the two methods would have given the same result.We use a closed boundary condition here since it leads to a more straightforward construction of the relevant quantum circuits due to the absence of an overall spin-flipping degeneracy explained earlier and long-range qubit interaction in the quantum hardware.In the following, we will mainly use the method with the [ T xy (t), T xy (0)]/A commutator, unless explicitly mentioned otherwise.
In Fig. 7b we compare the results of η(t f ) calculated     on two lattices of different sizes at the same conditions: j max = 1 2 , β = 0.3.We see the smaller lattice leads to a bigger finite size effect due to the much smaller state density.It is expected that the physical state density will increase exponentially with the lattice size (e.g.N p plaquettes) and j max as roughly (2jmax+1) 3Np 2 2Np [50], where the denominator is a rough estimate of the Gauss's law constraint effect.
We also find that the oscillation becomes less severe at higher temperatures due to the increase of the density of active states, as shown in Fig. 7c.

D. Trotter Error
In this subsection, we show the effect of implementing the Trotter decomposition in calculating the transport coefficients.In our case, we have two possible contributions: the error deriving from applying the Trotter decomposition in the real time evolution and that when preparing the thermal density matrix.
In order to quantify the different Trotter errors, we compute the G xy r (t), defined in Eq. ( 17), as a function of time t on a 3 × 3 lattice with j max = 1  2 for different values of the real time and imaginary time Trotter steps.We set ag 2 = 1 and β = 0.2 in lattice units.Panel (a) of Fig. 8 depicts the behavior of G xy r (t) as a function of time using the exact time evolution operators (without applying the Trotter decomposition).We constrain the y axis range in order to have a better visual comparison with the uncertainty plots.Panel (b) of Fig. 8 shows the absolute difference between the exact G xy r results3 and those when we apply the Trotter decomposition for real time evolution.The thermal density matrix is computed exactly.We observe that for ∆t > 0.1, we get a significant error: The magnitude of the green line at late time in (b) is of the order of 5 × 10 −4 , which is comparable with the black line in (a).Using smaller time steps is needed to improve accuracy, which may result in a deeper quantum circuit.Panel (c) of Fig. 8 shows the absolute error in G xy r as a function of time, when the Trotter decomposition is implemented for preparing the thermal state.The real time propagator is computed exactly.We observe that we have a negligible uncertainty.Panel (d) of Fig. 8 illustrates the results of combining the Trotter decomposition for the real time evolution and that for preparing the thermal state.As these plots show, the major Trotter error comes from the real time evolution.

E. Thermal State Preparation Success Probability
Earlier in Sec.IV A, we mentioned that in practical application for QCD, only the deconfined temperature regime roughly above 150 MeV is of interest.So the exponential decay of success probability in the thermal state preparation at low temperature may not be a problem.However, one may worry that even at high temperature, the success probability still decreases exponentially with the system size, just with a small prefactor in the exponent.Here, we provide a detailed analysis of this dependence.
The success probability of the thermal state preparation is given by where we have chosen E T = E 0 , and d H denotes the total Hilbert space dimension.We plot the success probability as a function of the lattice size on a periodic4 honeycomb lattice with j max = 1 2 in Fig. 9, where the lattice spacing a is fixed such that the coupling is ag 2 = 1.At this fixed coupling, the first excited state has an energy of about 6.2 above the ground state, i.e., E 1 − E 0 ≈ 6.2, which does not change much as the lattice size increases.High temperature in this case should at least be T ≳ 1/6.2 ≈ 0.16.Three temperatures are shown for comparison in the plot.
We see that the success probability decays exponentially with the system size, but the high temperature case has a very small prefactor in the exponent.
For physical application, one must make sure the energy density ε does not change with the lattice size at fixed a and β, since it captures the local equilibrium physics.Figure 5 and Eq.(37) imply that one has to increase j max as the lattice becomes bigger, in order to keep the energy density fixed.Our current computing resources forbid us to analyze this case so we leave it to future work.

F. Numerical Integration Error
The uncertainty control in numerical integration is well known.We take it as an individual issue to discuss because it plays an important role in practical calculations.
In quantum computing, one will calculate G xy r (t) at many time points and then integrate to obtain η as in Eq. ( 8).Under a given accuracy, if one can reduce the number of time points at which to compute G xy r (t) on quantum computers, one will reduce the amount of quantum resources needed to achieve the accuracy.
To demonstrate the integration error, we define the Riemann sum version of η(t f ) as where N t = t f /∆t.Its difference from the exact result |η sum (t f ) − η(t f )| is shown in Fig. 10 for a 4 × 4 lattice with ag 2 = 1 , j max = 1 2 and β = 0.3 (in lattice units).Comparing with Fig. 7a, we find the relative error up to t f = 300 is less than 5% for dt = 0.5 and is roughly 1% for dt = 0.2.
It is clear that the error grows with ∆t, and its magnitude scales as t 2 f /N t = t f ∆t.One can easily improve this by using the midpoint value in the Riemann sum, i.e., 2k−1 2 G xy r ( 2k−1 2 ∆t), whose error magnitude scales as

G. Fitting Uncertainty
Next we discuss the uncertainty when one extracts the bare value of the shear viscosity from the infinite time limit of η(t f ).As we have seen in Secs.V B and V C, the result of η(t f ) oscillates at late time due to the finite volume and local truncation effects, and thus, the infinite time limit is not well defined.To overcome this issue, we fit the time dependence of η(t f ) via some function that becomes constant at large t f .We consider two fitting functions where a 1,2 , b 1,2 , c 1,2 and d 2 are parameters.
Figure 11 shows the fitting results of η(t f ) as a function of t f on a 4 × 4 lattice with j max = 1 2 , ag 2 = 0.6 and β = 0.2 (in lattice units), where we can observe the exponential behavior.We use all the 500 time points up to t f = 500 (two neighboring time points are separated by ∆t f = 1) in the fit.The fitting is implemented through the python scipy curve fit function.From Eq. ( 8), the bare value of η is given by the plateaus coefficient (a 1,2 ) value.We see that the two functions are very similar.Indeed, the a i values fitted in Fig. 11 are identical up to fitting uncertainties: For f 1 , we find a 1 = 0.0587(3), while for f 2 , we have a 2 = 0.0590 (4).
At late time, due to the finite volume and local truncation effects, the signal becomes noisy.In order to reflect this real time fluctuation in our fitting uncertainty for the plateau value, we choose different ranges of t f used in the fit and estimate the average and uncertainty associated with them.For example, if we choose t f ∈ [0, 400], [0, 500], [0, 600], [0, 700] or [0, 1000] in the fit, we find the values shown in Table I when using the function f 1 .In the result section, we will use the function f 1 and apply this procedure to estimate the uncertainties 5 .11. η(t f ) as a of time, fitted by two different functional forms on a 4 × 4 lattice with jmax = 1 2 , ag 2 = 0.6 and β = 0.2 in lattice units.

H. Number of Shots and CNOT Gates
Finally, we estimate the number of times (shots) that one needs to repeat simulating the same quantum circuit and performing measurements, as well as the number of CNOT gates needed.It turns out that a huge number of shots is necessary to evaluate correctly the retarded Green's function from the commutator [T xy sum (t), T xy ij (0)] because its thermal expectation value is very small, around 10 −4 − 10 −2 .We assume that we prepare the thermal density matrix with an efficient quantum algorithm and will neglect the thermal state preparation part in the following estimate of the number of shots.As we mentioned previously, the implemented QITP th becomes inefficient when studying big systems at low temperature due to the large drop of the success probability.If we include the proposed QITP th algorithm for the thermal state preparation, the number of shots estimated below must be divided by the QITP th success probability.
Given a number of shots n shot and a probability of measuring the wavefunction of the quantum circuit containing the Pauli string Σ α of the sign s at time t to be in the computational basis state |b⟩, i.e., P s α (b), the measurement uncertainty can be computed by the binomial . The uncertainty of the measured retarded Green's function can be obtained from Eq. ( 30) The upper bound of this equation can be obtained by setting P s α (b) = 1 d H , ∀b where d H is the dimension of the system Hilbert space.Hence, we obtain the upper bound at a fixed t as where, in the second line, we use that T xy sum (0) is given by a sum of Pauli-z strings, so |⟨b|T 2 is used for all cases.A rough estimation gives that d T is smaller than 2N x N y .
Given the retarded Green's function G xy r (t), the required number of shots to achieve a relative error ϵ = is given by where we use that G xy r (t) ∼ 10 −3 .The number of shots needed increases polynomially with the lattice size because of the d 2 T factor in Eq. ( 43).The maximum number of CNOT gates that need to be implemented for studying the real time evolution is given by 70 per time step and lattice point: 64 for implementing the evolution driven by the magnetic part of the Hamiltonian and six for the electric part.Implementing the π 4 Σ α gate requires two CNOT gates for one α and one sign.Hence, on a N x × N y honeycomb lattice with j max = 1 2 , to implement the π 4 Σ α gate and real time evolution with a number of time steps N t , the number of CNOT gates we need is

I. Classical Computational Resource Estimate
In Appendix C, we study the computational time to classically simulate the quantum circuit for the real time evolution via the matrix product state method.The computational time grows exponentially with the Trotter steps when the bond dimension is large, which is expected to be necessary to describe states with volumelaw entanglement, as is usually the case when the system thermalizes at late time.Performing such real time calculations by exact diagonalization on a large lattice with high j max truncation also requires exponentially growing resources.For example, the Hilbert space on a 3 × 3 lattice with j max = 1 is 519233 and that on a 5 × 5 lattice with j max = 1 2 is 2 25 , when all external links are in j = 0 states.The needed space to store a single T xy ij matrix in the energy eigenbasis is about 1.1 TB and 4.5 PB, respectively (single-precision float).Quantum computing may be able to overcome this difficulty.

η/s in the continuum
We will show results of the ratio between the shear viscosity η and the entropy density s as a function of temperature in the continuum limit.The continuum limit is taken via extrapolating toward a = 0.In this work, we will only consider the renormalization group equation of the coupling as a varies as written in Eq. (34).Additional operator renormalization is left for future work which may lead to a 20% change roughly [51].
We start this procedure by fixing a "physical" temperature, which is kept invariant as lattice spacing a changes.We set this "physical" temperature T 0 (β 0 = 1 T0 ) to be the temperature when ag 2 = 1, which is a number in the lattice unit corresponding to ag 2 = 1.Then the temperatures at other couplings (i.e., other lattice spacings) are T = ag 2 T 0 (β = β0 ag 2 ).It would also be useful to fix the temperature scale by the confinement-deconfinement crossover temperature, which is left for future work.
The black circles of Fig. 12 represent the obtained results for the ratios η(β0,ag 2 ) s(β0,ag 2 ) at six different couplings ag 2 = {0.4,0.5, 0.55, 0.6, 0.65, 0.7} for β 0 = 0.2.We stop at ag 2 = 0.4 rather than going to smaller couplings since at such a small coupling, the low-j max truncation effect is very large, signaling large oscillation of η(t f ) at late time, which deteriorates the fitting of the plateau value.The vertical uncertainty bars associated with the black points describe the fitting uncertainties explained in Sec.V G.The fitting uncertainties are very small except for ag 2 = 0.4.It can be seen that an exponential function can describe the trend of the black points, where c 0 , c 1 , c 2 are fitting parameters.
To quantify the systematic uncertainty of using the function f (ag 2 ), we choose three different data sets to perform the fit: {ag 2 } 1 = {0.5, 0.55, 0.6, 0.65}, {ag 2 } 2 = {0.4,0.5, 0.55, 0.6, 0.65} and {ag 2 } 3 = {0.4,0.5, 0.55, 0.6, 0.65, 0.7}, which are shown in Fig. 12 in green, orange and blue, respectively.The fitted parameter values are listed in Table II.The band with the same color of the line represents the relative uncertainty at one sigma of the fits.The band grows at large ag 2 since the fitting function grows exponentially there and a small uncertainty in either c 1 or c 2 leads to a big uncertainty after proper error propagation.
The continuum limit for η s at ag 2 = 0 is obtained as c 0 + c 1 .The results for β 0 = 0.2 are shown in the last column of Table II for the three different fitting ranges of ag 2 .The three values are compatible within two-sigma error bars.However, we observe a change in the continuum limit value when the data point at ag 2 = 0.7 is included in the fitting.We attribute this to the potentially larger lattice discretization effect at bigger lattice spacing.
Iterating this procedure for different "physical" temperatures T 0 , we obtain the temperature dependence of η s in the continuum.We show in Fig. 13 results from all the three fitting data sets.The blue and orange dots are slightly shifted horizontally for better visualization.We do not study temperatures higher than β 0 = 0.15 since higher j max truncation is needed to accurately describe highly excited states, as discussed in Sec.V B. The uncertainty grows rapidly at lower temperatures (e.g., β 0 = 0.225), since not many states of the theory are effectively contributing to the retarded Green's function (suppressed by e −β0E ) and then the density of contributing states is not large enough to suppress the real time fluctuation in η(t f ) due to our small lattice and local Hilbert space truncation, as seen in Sec.V C. 2 .We slightly shift the data horizontally for better visualization of the three fittings using different data sets.
Our results of η s are consistent with the holographic result 1 4π within uncertainties, which is shown as the dashed line in Fig. 13.We also observe a trend of decrease in η s as temperature increases from the blue dots.However, this trend is not obvious from the green and orange dots.All these should be further studied on bigger lattices with higher j max truncation in the future, to better understand the finite volume and local Hilbert space truncation effects.

Structure of Spectral Function
We also study the off-diagonal matrix elements of T xy in small frequency ω ranges, which are related to the spectral function that is defined as When ω is small, |⟨n| T xy |m⟩| 2 is closely related to ρ xy (ω) ω : Our results of |⟨n| T xy |m⟩| on a 4 × 4 lattice with j = 1 2 are shown in Fig. 14 for two values of ag 2 : 0.6 and 1.0, where we use eigenstates in the energy windows 15 < E n , E m < 17 and 26 < E n , E m < 28, respectively.Previous calculations showed no structure in ρ xy (ω) ω when ω is small (which means it is flat in ω) for strongly coupled N = 4 supersymmetric Yang-Mills theory while there is a peak structure in perturbative QCD results [52,53].Our results exhibit peak structures at small ω at the two couplings studied.The peak is broader when the coupling is smaller (note the x-axis scale is different in the two plots in Fig. 14).Whether these peak structures persist with higher j max truncation should be studied in the future.

B. Quantum Computing Results
We test the quantum circuit proposed in Sec.IV to evaluate G xy r (t) on a 2 × 2 lattice with j max = 1 2 and ag 2 = 1 for β = 0.15 (in lattice units).Parts of the thermal state preparation and real time evolution quantum circuits can be found in Appendix B.
We will show results of the imaginary part of the retarded Green's function (it is purely imaginary) at different times, which are calculated from the commutator [T xy sum (t), T xy 10 (0)].We do not evaluate the commutator by using T xy sum (0) at t = 0 because it requires 16 quantum circuits per time step.Indeed, for a N i × N j lattice, the decomposition of T xy sum contains 2N i N j Pauli strings [see Eqs.(20) and (26)].Since the commutator formula in Eq. ( 27) is valid only for a single Pauli string, we need to prepare different circuits for each of the 2N i N j terms in T xy sum ; for each commutator, we have two quantum circuits with different signs in the exponent.
We first run the quantum circuits on the Quantinuum H2-1E emulator [54] to compute the first two time points with a real time Trotter step of ∆t = 0.025 and a single imaginary time Trotter step of ∆τ = β 2 = 0.075.For each quantum circuit, we measure the evolved circuit, which collapses the wavefunction onto some basis state, and repeat for a total number of n shot = 500 times 6 .We expect that running the same quantum circuits on the real machine would give us very close results because the Quantinuum emulator is very close to the real hardware [55,56].For the given parameters, the success probability of obtaining the thermal density matrix is 0.189, which means only 18.9% of the total 500 shots give useful measurement outcomes.The required number of CNOT gates to prepare the thermal state and perform the π 4 Σ α gate is 76.For each real time step, the number of CNOT gates is 34.So the total number of CNOT gates grow linearly as 76 + 34N t with the total real time steps N t .The number of implemented CNOT gates differs from Eq. ( 44) that describes the general case when we have a seven-body term in the Hamiltonian.Here, for a 2 × 2 lattice, we only have four-body interactions.
The results obtained from the Quantinuum H2-1E are shown in blue squares in Fig. 15, where the exact diagonalization results are also shown for comparison, as well as the results obtained from the classically evolved quantum circuits with the same real and imaginary time Trotter steps, which we call "exact Trotter".The perfect agreement between the exact diagonalization results and the classically evolved quantum circuit results show the Trotter errors are negligible.The associated statistical uncertainties with 500 shots are huge since the small expectation values of the commutator (which are 10  15.Imaginary part of the retarded Green's function as a function of time, obtained from exact diagonalization (black solid line), classically evolved quantum circuits (red dashed line), Quantinuum H2-1E emulator with n shots = 500 shots (blue dots) and noiseless IBM simulator using 10 7 shots (green dots).The red dashed line is obtained by connecting nearest points at intervals of ∆t = 0.025 with straight lines.In the right panel, we zoom in on the y-axis so the statistical uncertainty associated with the green dots can be seen.
shots (more than a million) to reconstruct the observable from the projective samples of the state.To confirm this, we run the same quantum circuits on the noiseless ibmq qasm simulator IBM simulator [57,58] with n shot = 10 7 and obtain results shown as green circles in Fig. 15.We observe that the IBM simulator results obtained from 10 7 shots have much smaller statistical uncertainties and are compatible with the exact results, demonstrating the validity of the proposed quantum circuit in computing the retarded Green's function to obtain the shear viscosity.

VII. CONCLUSIONS
In this work, we considered calculating the shear viscosity nonperturbatively by using the Hamiltonian lattice approach for the 2+1-dimensional SU(2) Yang-Mills theory.The shear viscosity is obtained from the real time retarded Green's function of the stress-energy tensor via the Kubo formula.We included the renormalization of the coupling when taking the continuum limit but not additional operator renormalization.By exactly diagonalizing the theory on a 4 × 4 honeycomb lattice with j max = 1 2 , we found the ratio of the shear viscosity and the entropy density is consistent with the well-known holographic result 1 4π at several temperatures.On the other hand, our results showed a peak structure in the spectral function divided by the frequency (i.e., ρ xy (ω) ω ) when the frequency is small, which is qualitatively different from the holographic result but similar to the perturbative one.
The finite volume and local Hilbert space truncation effects are probably large in our current studies, which motivate future calculations with higher j max values in the Hilbert space truncation on bigger lattices.These calculations probably can neither be easily done by exact diagonalization due to the memory limitation nor by the matrix product state classical simulation method due to the exponential growth of the computational time.
So we propose a quantum algorithm to evaluate the shear viscosity on quantum devices.We analyzed various systematics of the calculation.We tested the reliability of the quantum algorithm on a 2×2 lattice with j max = 1 2 and found the quantum results agree with the classical ones, despite the huge number of shots needed to accurately evaluate the retarded Green's function.Moreover, the success probability of the QITP th algorithm used in preparing the thermal state decreases with the lattice size exponentially, with a very small prefactor in the exponent at high temperature.Therefore, it would be desirable to upgrade or develop more efficient methods to prepare the thermal state for the evaluation of the retarded Green's function.Future physical calculations of the shear viscosity beyond our current study would probably require robust large-scale quantum computers that are capable of performing a large number of shots.

Appendix C: Matrix Product State Classical Simulation
Classically, one can simulate the quantum circuit by implementing the matrix product state (MPS) algorithms [75,76].We simulate just the real time evolution part of the quantum circuit on the Aer qiskit emulator with the MPS flag activated, starting from the default initial state of the quantum processor |000 . . .0⟩ on a 5 × 5 lattice with j max = 1 2 , ag 2 = 1 and ∆t = 0.05.We investigate the dependence of the computational time on the bond dimension.Its value corresponds to a trade- a For more accurate results, we have to enhance the number of shots with a further increase of the computational time.
off between faster calculations (lower value) and better accuracy (higher value).
Figure 18 shows the computational time as a function of the number of Trotter steps (t s ) for different MPS qiskit emulators [57,58]  In both cases, we impose that at zero time step, the computational time is 0. The obtained fitting parameters are reported in the first two rows of Table III.
In the same table, using the fitted parameters, we estimate the order of magnitude of the needed computational time to perform t s time steps.We have to evolve our system for long time in order to evaluate the shear viscosity by integrating over time.The total computational time would be given by the sum of computational times for all the time steps.Moreover, we have to add the computational time to prepare the thermal state in realistic analyses.The exponential trend and the numbers shown in Table III suggest that it is almost impossible to use the classical MPS simulation of the quantum circuit to accurately calculate the shear viscosity (imposing higher B D values).
Indeed, in the quick calculations with a lower bond dimension in the MPS algorithm, the obtained results become less accurate after a long evolution time.We implement the real time evolution using the MPS qiskit function fixing the bond dimension B D = 20 and B D = 100 on a 5 × 5 honeycomb lattice with j max = 1 2 , ag 2 = 1 and ∆t = 0.05. Figure 19

FIG. 2 .
FIG. 2. Schematic diagram of the implemented quantum circuit to calculate the retarded Green's function of T xy on a 2 × 2 lattice with jmax = 1 2 .The e − βH 2 where |b⟩ denotes the computational basis states, which, e.g., on a 2 × 2 lattice with j max = 1 2 , are |0000⟩, |0001⟩, . . .|1111⟩ in the spin representation.The symbol P ± α (b) indicates the measured probability of the |b⟩ state for the circuit with Σ α = e ±i π 4 Σα evolved from time 0 to t.The time dependence of P ± α (b) is omitted for notational simplicity.

FIG. 4 .
FIG.4.Internal energy and entropy densities in lattice units as functions of temperature for several values of jmax on a 2 × 2 lattice with ag 2 = 1.The dashed lines are expectations from the continuum theory.

FIG. 7 .
FIG. 7. η defined in Eq. (8) as a function of time.(a) Results on a 4 × 4 lattice with jmax = 1 2 and ag 2 = 1 at β = 0.3 in lattice units, obtained from two ways that are equivalent if the system obeys translational invariance.(b) Results on two lattices with different sizes but the same jmax = 1 2 and ag 2 = 1 at β = 0.3.(c) Results on a 4×4 lattice with jmax = 1 2

FIG. 8 .
FIG. 8. Trotter error in computing G xy r (t) as a function of time on a 3×3 lattice with jmax = 1 2 and ag 2 = 1 for β = 0.2 in lattice units.(a) Exact G xy r (t) results obtained from Eq. (17).(b) Uncertainty when Trotter decomposition is only applied in real time evolution.(c) Uncertainty when Trotter decomposition is solely applied in preparing the thermal state.(d) Uncertainty when Trotter decomposition is applied in both real time evolution and thermal state preparation.

2 FIG. 9 .
FIG. 9. Success probability of thermal state preparation as a function of the honeycomb lattice size for several temperatures.The lattice spacing a is fixed corresponding to ag 2 = 1 and jmax = 1 2 is used.

5 FIG. 10 .
FIG.10.Numerical integration uncertainty as a function of final time for three different finite elements of time.

ag 2
FIG. 12. Results of the coupling dependence of η s for β0 = 0.2 on a 4 × 4 lattice with jmax = 1 2 .Black points represent the calculated η s at different couplings, lines indicate the fitting results and bands describe one sigma uncertainty of each fitting.

2 )
FIG. 13.Obtained ηs results as a function of β0 on a 4 × 4 lattice with jmax = 1 2 .We slightly shift the data horizontally for better visualization of the three fittings using different data sets.
Using the fitting result of Fig.18, we estimate the order of magnitude of the required computational time to implement the real time evolution for specific time steps and 1024 shots a on a laptop for different bond dimension values on a 5 × 5 lattice with ag 2 = 1 and jmax = 1 2 .

FIG. 19 .
Figure 18  shows the computational time as a function of the number of Trotter steps (t s ) for different MPS qiskit emulators[57,58] with bond dimensions B D = 20, 100 and without truncation, namely B D = ∞, which indicates the exact result.All the simulations are −4 − 10 −3 ) require a significant number of FIG.
shows the obtained relative error defined as B D (t) indicates the probability of the initial state |i⟩ at time t using B D bond dimension, i.e., p B D (t) = ⟨i|U t |i⟩.Different lines correspond to the obtained results starting from a random bit state, generated by random applications of X gate on the default initial state.We can observe a rapid increase of the relative error with the time step.Therefore, the shear viscosity calculations that require long time evolution would be extremely imprecise using MPS algorithms with low bond dimensions.