Shear viscosity of classical Yang-Mills field

We investigate the shear viscosity $\eta$ of the classical Yang-Mills (CYM) field on a lattice by using the Green-Kubo formula, where the shear viscosity is calculated from the time-correlation function of the energy-momentum tensor in equilibrium. Dependence of the shear viscosity $\eta(g,T)$ on the coupling $g$ and temperature $T$ is represented by a scaling function $f_\eta(g^2T)$ as $\eta(g,T)=Tf_\eta(g^2T)$ due to the scaling-invariant property of the CYM. The explicit functional form of $f_\eta(g^2T)$ is successfully determined from the calculated shear viscosity: It turns out that $\eta(g,T)$ of the CYM field is proportional to $1/g^{1.10-1.88}$ at weak coupling, which is a weaker dependence on $g$ than that in the leading-order perturbation theory but consistent with that of the"anomalous viscosity"$\eta\propto 1/g^{1.5}$ under the strong disordered field. The obtained shear viscosity is also found to be roughly consistent with that estimated through the analysis of the anisotropy of the pressure of the CYM dynamics in the expanding geometry with recourse to a hydrodynamic equation.


I. INTRODUCTION
It is widely believed that the initial stage of relativistic heavy-ion collisions is well described by the classical Yang-Mills (CYM) field [1][2][3]. The real-time lattice simulations of the CYM field have provided important insights into the non-equilibrium dynamics of this dense gluon matter, which is known as glasma. In the current understanding of the early dynamics of the heavy-ion collisions, the instabilities [2][3][4][5][6][7] and the chaoticity [8] of the CYM field lead to a rapid thermalization of the glasma, so that a hydrodynamic description makes sense [9]. The hydrodynamic analyses showed that the created matter may be an almost perfect fluid whose shear viscosity is close to the lower bound predicted by superstring theories with an Einsteinian classical limit [10]. The mechanism that leads such a near-minimal shear viscosity has been studied from various points of view [11][12][13][14][15][16].
Then a natural question that arises here is whether the CYM-field description of the early stage can be further extended to account for the transition region from the glasma to the low-viscosity, near-thermal gluonic fluid. To answer this question, the equilibrium transport properties of the CYM field must be investigated. It was shown in Ref. [17] that the system reaches a quasistationary state and the entropy increases slowly after a rapid isotropization of pressure. In Ref. [3], the shear viscosity is deduced through the fitting of the time dependence of the energy density of the CYM field with that expected as a solution of the viscous hydrodynamics equation, while the pressure anisotropy still remains. These analyses suggest that the CYM field may be relaxed by the intrinsic dynamics to a quasi-local equilibrium state. Hence it should make sense to explore its transport properties such as the viscosities.
One of the ways to extract the shear viscosity of the CYM field is to fit the hydrodynamic parameters to the non-equilibrium time evolution of the CYM field as done in Ref. [3]. Another, more standard way is to use the Green-Kubo formula obtained in the linear response theory. The latter approach has been successfully applied to the classical scalar theory [15,16].
In this article, we investigate the real-time dynamics of CYM field close to thermal equilibrium and extract the shear viscosity of the CYM field at thermal equilibrium. We employ the Green-Kubo formula in order to investigate the shear viscosity arising from the CYM field's dynamics. In the Green-Kubo approach, the shear viscosity can be extracted from the thermal expectation value of the time-correlation function of the energy-momentum tensor. We evaluate the expectation value as the ensemble average measured with thermally equilibrated CYM field configurations. We also make a detailed analysis of the time-correlation function and the power spectrum (the Fourier transform of the time-correlation function) in order to gain deeper insight into the relaxation of the energy-momentum tensor and the origin of the shear viscosity.
It is worth noting that our analysis is based on the scaling property of the CYM theory. The classical equation of motion of the CYM theory is invariant under the scale transformation, g → γg, A → A/γ and E → E/γ. As the thermal classical gauge theory must be regulated on a lattice with spacing a, this scale invariance implies that intensive quantities in the thermal CYM theory, such as the shear viscosity, are functions of the dimensionless quantity g 2 T a. Setting a = 1, any intensive quantity thus becomes a function of the product of the squared coupling g 2 and the temperature T measured in lattice units. One of the main results of the present work thus is to determine the scaling function of the shear viscosity, f η (g 2 T ), over a wide range of the scaling variable g 2 T .
This article is organized as follows: In Sect. II, we introduce the lattice formulation of the CYM field, the Green-Kubo formula for the shear viscosity, the classical field ensemble, and some properties of the CYM field on a lattice. In Sect. III, we show the numerical results of the time-correlation function of the energy-momentum tensor, the power spectrum, and the shear viscosity. We discuss the g and T dependence of the obtained shear viscosity in terms of the lattice unit by using the scaling functions. In Sect. IV, we summarize our work.

II. CLASSICAL YANG-MILLS FIELD AND ITS SHEAR VISCOSITY
A. Classical Yang-Mills Field Theory on Lattice We consider the CYM theory on a L 3 lattice in the temporal gauge, A a 0 (x) = 0. Its Hamiltonian in the noncompact formalism is given as (using a = 1): Here the electric and magnetic gauge fields, (E a i (x), B a i (x)), are given by ) are the canonical variables and ∂ F i is the forward difference operator. The time evolution of the CYM field is obtained by solving the Hamilton equation of motion, By solving this equation of motion, we calculate the time-evolution of the space-averaged off-diagonal matrix elements of the energy-momentum tensor defined as on the lattice. Here (E a i (x), B a i (x)) are the electric and magnetic fields defined at a shifted spatial point It should be noted that the electric and magnetic fields (E a i (x), B a i (x)) are defined at the same spatial position for different i because A a i (x) resides on the link between x and x +î. Therefore, (E a i (x), B a i (x)) make the matrix elements more symmetric in the space directions than (E a i (x), B a i (x)).

B. Green-Kubo Relation
The shear viscosity in any systems close to equilibrium can be evaluated by means of the Green-Kubo formula [18]: where T and V denote temperature and volume in lattice units, respectively. τ ij (t) is the space-averaged offdiagonal matrix element of the energy-momentum tensor, and · · · denotes the expectation value in equilibrium. C(t) is the direction-averaged time-correlation function of τ ij (t).
While the correlation disappears in the long timeseparation limit t → ∞, and thus C(t) is expected to approach V τ ij (0) 2 eq = 0 in this limit, the integral in (9) converges only when the correlation function falls off faster than 1/t or exhibits oscillatory behavior with a decreasing amplitude. We will check the convergence in Sect. III C. We will also compute the Fourier transform of C(|t|), and discuss the structure of ρ(ω) in comparison with that in the scalar theory in Appendix B.

C. Classical Field Ensemble
For the calculation of the shear viscosity via the Green-Kubo formula, we need the thermal expectation value of the time correlation of the energy-momentum tensor. We evaluate the expectation value from the ensemble average of the classical field configurations in thermal equilibrium as where N conf is the total number of configurations used in the evaluation, and O i is the observable measured with the i-th configuration. Statistical errors of the expectation values, which are of the order of 1/ √ N conf , are estimated and taken into account. We prepare the ensemble of thermally equilibrated CYM field configurations by the long time evolution of the classical CYM field starting from the random initial configurations. The procedure we have used is explained in detail in Sect. III A.

D. Equipartition of Electric Field Energy
The electric field strength can be utilized to examine thermalization and also to measure the temperature. In equilibrium, the distribution of canonical variables (A, E) is given by the canonical partition function Z, where Z B = DA exp(−H B /T ) is the magnetic field part of the partition function: From the second line in (13), the average squared electric field is found to be equal to the temperature, (E a i (x)) 2 = T at each spatial position x. From the third line, we find that the same relation applies to the mean square Fourier component of the electric field for each lattice momentum k = 2π(n x , n y , n z )/L (n i = 0, 1, · · · L − 1). Note that E a i (k) and E a i (−k) are not independent but related as E a i (−k) = (E a i (k)) * , since E a i (x) is real. Note that, while E a i (k) itself is not gauge invariant, the property (15) does not depend on the chosen gauge. In Sect. III A, we will examine thermalization from the momentum independence of |E a i (k)| 2 and obtain the temperature from |E a i (k)| 2 .

E. Scaling Property
Here we explain the scaling property of the CYM theory, which plays an important role in our analysis of the shear viscosity. The equation of motion of the CYM field is invariant under the following scaling transformation, Then, the lattice temperature defined in Eq. (15) scales as where we have omitted the color and vector indices. This relation can be rewritten as (γg) 2 T γg = g 2 T g . Thus, once we prepare the ensemble at a given value of g 2 T , we can obtain the thermal average of observables at various values of (g, T ) having the same g 2 T by simply rescaling the observables. For example, the time-correlation function C(t), its power spectrum (the Fourier transform of C(t)) ρ(ω), and the shear viscosity η are given as where we have used the scaling property of the energymomentum tensor, These relations show that there exist scaling functions of These scaling functions are related with each other as follows, When multiplied by appropriate powers of T , these scaling functions give the values of C(t; g, T ), ρ(ω; g, T ), and η(g, T ). In this study, we perform calculations at several values of g with T ≈ 1. Results will shown in terms of the above scaling functions.

A. Calculation Setup and Initial Configuration
Our calculations were performed on 16 3 , 24 3 and 32 3 lattices with periodic boundary conditions for the SU (2) Yang-Mills field keeping the average energy per degree of freedom at unity, which corresponds to the CYM system with the temperature of T ≈ 1. The coupling constant was varied in the wide range g = 0.15 − 20. The equation of motion is solved by leapfrog integration with ∆t = 0.01. The thermal expectation value of an observable is evaluated from the ensemble average using N conf = 1000 thermal configurations.
Here we explain the procedure for generating a set of thermally equilibrated gauge field configurations. First, we generate totally N conf field configurations, (A (i) , E (i) ) (i = 1, 2, · · · , N conf ), at initial time t = t ini , as where ∂ B i is the backward difference operator, R (i),a j (x) is a Gaussian random number with unit variance: and N E is a normalization factor chosen to make the average energy per each degree of freedom equal to unity. This initial configuration satisfies the lattice Gauss law: Next, we evolve each configuration with the classical equation of motion until it equilibrates. We confirm the equilibration by checking for equipartition of the energy given in Eq. (15). In Fig. 1, we show the time evolution of the reduced chi-square of the fit |E 1 1 (k)| 2 = T fit with T fit as fitting parameter using the N conf = 1000 configurations with g = 1 and L = 32. The reduced chi-squared value reaches unity for t − t ini > t eq ∼ 200, which indicates that Eq. (15) holds after this time. We can consider the configurations at t − t ini > t eq ∼ 200 to be in thermal equilibrium for t − t ini > t eq . In Appendix A, we confirm the equilibration using another quantity, the local electric energy distribution. Finally, we extract the temperature of the thus prepared configurations via Eq. (15). In Fig. 2, we show |E 1 1 (k)| 2 as a function of the lattice frequency ω k , at g = 1 on the 32 3 lattice. The expectation value has been evaluated from the configurations at t − t ini = 200. It appears that |E 1 1 (k)| 2 is around unity at any ω k as we expected. By fitting a constant T fit to the result, we obtain T fit = 1.16 (χ 2 r = 1.00).
The reduced chi-squared is around unity (χ 2 r ∼ 1), and thus we can consider the fitted value of T fit = 1.16 as the temperature of the system with (g, L) = (1, 32). In Table I further below, we summarize the values of g 2 T used in our calculations; they cover the range 0.0256 ≤ g 2 T ≤ 499.

B. Time correlation function of energy-momentum tensor
We shall now discuss the time-correlation function of the energy-momentum tensor, C(t), evaluated in the thermally equilibrated CYM field.
In Fig. 3, we show the normalized time-correlation functions, f C (t) = C(t)/T 2 , in the short time range, t < 5. The time-correlation function is seen to decrease rapidly over the interval 0 < t < 1.5. In this range it can be well described by the sum of two damped oscillators , superimposed on a constant, as shown by the dashed lines. Therefore, the rapid decrease in the short time range is found to come from the two damped oscillatory motions.
In Fig. 4, we show f C (t) in the long-time range, 5 < t < 50. The remaining correlation is found to decay more slowly in this later time period, and the decay time scale significantly depends on the value of g 2 T . Before discussing the shear viscosity, we consider the Fourier transform of the time-correlation function C(|t|), referred to as the power spectrum, whose low-frequency limit gives the shear viscosity. The power spectrum of C(|t|) is given by where ∆t is the time step size and N t is the total number of steps in the time evolution, and is half of the data points in summation. In the last line, we represent the integration over the interval of [−∞, ∞] by the discrete Fourier transformation with the period 2N t , in which the power spectrum ρ(ω) has entries at each discrete frequency ω = 2πn/(2N t ∆t) (n = −N t , −N t +1, · · · N t −1).
In actual calculations, the upper bound of the integration, t cut = N t ∆t, is finite but must be large enough so that the integral (summation) in (36) converges. As an example of the convergence check, we show in Fig. 5 the integral of the normalized time-correlation function C(t)/T 2 over the interval 0 ≤ t ≤ t cut = N t ∆t, which is directly related to the shear viscosity η via (9) and corresponds to ρ(0) in (36). The integral at each value of g 2 T shows convergent behavior at large t cut , while small oscillations around a constant are still visible. Thus, the integral in Eq.(36) can be approximated by the integral over the finite interval as long as N t is large enough. This is in contrast with the classical φ 4 theory case, in which a very long tail appears in the time-correlation function and must be subtracted in advance [16]. oscillator part corresponds to the broad bump of the spectral function, which has large strength in the high frequency region, 4 ≤ ω ≤ 6. The remaining part of the power spectrum corresponds to the long-time decay part of the time-correlation function and show the peaks at ω ≤ 1. The number of peaks decreases and the peak heights become smaller with increasing g 2 T . The peak at ω = 0 remains even at large g 2 T , and is the highest peak in the g 2 T region shown here. The ω = 0 peak is sharp at small g 2 T but seems to be a continuous function of ω. This supports the convergence of the integral in (36) at ω → 0. It would be interesting to discuss origin of the peaks other that at ω = 0, which may be due collective modes in the low frequency region. However this is beyond the scope of this work, which is concerned with the shear viscosity. At larger values of g 2 T , the power spectrum becomes smoother and is increasingly dominated by the damped oscillator part. In Appendix B, we show f ρ (ω) in the extremely large g 2 T region. The broad bump in the high frequency region caused by the damped oscillatory behavior in the timecorrelation function was also found in Ref. [16], where it was conjectured that the damped oscillatory behavior is caused by the propagation of two modes having momenta of k and −k (two-momentum mode), where φ k denotes the Fourier transform of the classical scalar field and · · · denotes the thermal expectation value. The frequency difference between different twomomentum modes causes the phase decoherence and the damping of the time-correlation function. Since this can take place even in the free field case and the damped oscillations in the CYM field have shapes similar to those in the φ 4 theory, the damped oscillation in the early stage and the bump around ω ∼ 5 may be also caused by the interference of the two-momentum modes.

D. Shear viscosity
We now discuss the shear viscosity. The shear viscosity η is obtained by taking the low frequency limit of the power spectrum, η = lim ω→0 ρ(ω)/T . In Table I, we summarize the g 2 T dependence of the normalized shear viscosity, f η (g 2 T ) = η/T on the 16 3 , 24 3 and 32 3 lattices. As shown in (25), η/T is invariant under the scaling transformation given in (17) and a function of of the single parameter g 2 T . In Fig. 7, we show f η as a function of g 2 T . The results on the 16 3 , 24 3 and 32 3 lattices are consistent with each other. We find that f η (g 2 T ) is a rapidly falling function at small values of g 2 , but flattens out in the region g 2 T > 1. This weaker fall-off is readily understood from the behavior of f C (t) and f ρ (ω). The damped oscillation in f C (t) is insensitive to g 2 T as seen in Fig. 3, and its contribution to f ρ (ω) in the low frequency region, which is responsible for the viscosity, becomes more dominant and remains even in the larger g 2 T region.
Next we analyze the functional form of f η (g 2 T ) and understand the g 2 T dependence of f η (g 2 T ) in the range, 0.0256 < g 2 T < 499, by fitting our results to a polynomial function with parameters α, β, γ and δ, where x = g 2 T . The form of F (x) is motivated by the results shown in Fig. 7. We expect that the first term of F (x) describes the contribution from the long-time tail in f C (t), which seems to increase like an inverse power of g 2 T and dominates at small g 2 T . The second term of F (x) is expected to describe the contribution from the damped oscillation in f C (t), which is less sensitive to g 2 T . The fit to the f η (g 2 T ) obtained from the Green-Kubo formula in Eq. (9) on the 32 3 lattice yields the following parameters: F (x) : α = 0.09 ± 0.07, β = 1.49 ± 0.39, γ = 0.33 ± 0.06, δ = 0.35 ± 0.07.
In Fig. 8, we compare F (g 2 T ) and f η (g 2 T ) on the 32 3 lattice. The fit function F (x) describes f η (g 2 T ) in this g 2 T range well. These results imply that the shear viscosity of the CYM field is proportional to 1/g 1.10−1.88 at weak coupling, which indicates that the shear viscosity in the classical limit depends more weakly on g compared with the leading-order perturbative theory result, η ∝ 1/[g 4 ln (g −1 )]. Such a weak dependence on g is consistent with the g-dependence of the "anomalous viscosity" under the strong disordered field, η ∝ 1/g 1.5 [12].

E. Comparison with previous estimates
We now compare our results with those of previous estimates obtained by analyzing the CYM field evolution in terms of the viscous hydrodynamic equation [3]. Since we have evaluated the shear viscosity using the Green-Kubo formula based on the linear response theory, the so-obtained shear viscosity characterizes the relaxation process around the equilibrium and is not sensitive to the initial conditions. In Ref. [3], Epelbaum and Gelis deduced the shear viscosity of the CYM field from the anisotropy of the pressure by comparing the energy density evolution in the Analytical fit to the shear viscosity normalized by 1/T , fη = η/T , on the 32 3 lattice with the fit function expanding geometry to first-order viscous hydrodynamics. The shear viscosity in Ref. [3] was obtained for a lattice spacing a −1 = Q s = 2 GeV and g = 0.5 using the phenomenological glasma initial condition on a 64 × 64 × 128 lattice. The authors of Ref. [3] obtained the result η/ε 3/4 0.3, where ε is the energy density. Shear viscosity normalized by 1/T , fη = η/T . The red, blue and green points show the results with T phys = T /a = 0.2 GeV, 0.4 GeV and 0.8 GeV, respectively, at a = (2 GeV) −1 . The black point shows the value of the shear viscosity in the expanding and anisotropic CYM field, which is estimated in [3] In Fig. 9, we show the shear viscosity normalized by ε 3/4 with ε ≈ 3(N 2 c − 1)T = 9T at lattice spacing a = (2 GeV) −1 for the physical temperatures T phys = T /a = 0.2 GeV, 0.4 GeV, and 0.8 GeV as a function of g. These temperatures in physical units are in the range of temperatures reached in relativistic heavy-ion collisions at RHIC and LHC at the start of the hydrodynamical evolution. The ratio η/ε 3/4 in the present work is smaller than the value obtained in Ref. [3] at g = 0.5 by around 30-40 %. These differences may reflect the effect of a partial persistence of the initial condition and incomplete thermalization in the approach of Ref. [3]. However, we find it encouraging that our results obtained in equilibrium lie in the same numerical range as the shear viscosity value deduced from nonequilibrium dynamics in Ref. [3].
It is worth noting that these values of the shear viscosity are much smaller than the leading-order perturbative result η pert /T 3 ≈ 1200 for g = 0.5 [11]. This difference is conjectured [3] to be the manifestation of the anomalously small viscosity for systems made of strong disordered fields [12], an effect that is not included in the perturbative calculation.

IV. SUMMARY
We calculated the shear viscosity of the classical Yang-Mills (CYM) field on a lattice by applying the Green-Kubo formula to the time-correlation function C(t) of the energy-momentum tensor in equilibrium. The time evolution of the CYM field was calculated for several values of the coupling constant g starting from equilibrated configurations prepared at temperature T ≈ 1, and C(t) was evaluated as the ensemble average. The shear viscosity as a function of the coupling and temperature η(g, T ) was obtained from the low-frequency limit of the power spectrum of C(t). The dependence of the shear viscosity on g and T was represented as η(g, T )/T = f η (g 2 T ) with a scaling function f η in accordance with the scaling property of the CYM theory. The functional form of the scaling function f η (g 2 T ) was extracted as a polynomial fit function.
The time-correlation function C(t) was found to exhibit damped oscillatory behavior at early times followed by a slow decay. The damped oscillatory behavior is reflected in a broad bump around ω ≈ 5 in the power spectrum and may be caused by the decoherence among the two-momentum modes, which is attributed to a universal early-stage dynamics showing a damped-oscillatory behavior as observed in the case of the classical φ 4 theory [16].
The slowly decaying part of C(t) produces a sharp peak at ω = 0 in the power spectrum. With decreasing g 2 T , the height of the peak at ω = 0 increases and its width narrows. Thus the ω = 0 peak suggests the appearance of a collective mode with long life-time, which is most naturally identified with the hydrodynamic mode.
The scaling function f η rapidly increases with decreasing g 2 T at small g 2 T and its dependence on the coupling, f η ∝ 1/g 1.10−1.88 , is found to be much weaker than the perturbative estimate. This weaker dependence may show the realization of anomalous viscosity under the strong disordered field [12]. The value of the shear viscosity in thermal equilibrium obtained in our analysis is found to be roughly consistent with that obtained from the glasma energy density evolution in the boostinvariant expanding geometry [3]. Thus the validity of the estimate of the shear viscosity from the energy density is confirmed.
While the transport properties of the CYM field around equilibrium may be relevant to the small shear viscosity in high-energy heavy-ion collisions, there still remain problems to be considered in order to discuss the shear viscosity around equilibrium in a more rigorous way. The behavior of high-dimensional observables such as the time-correlation function of the energy-momentum tensor, C(t), is dominated by hard thermal particles with momenta of order temperature and is sensitive to the ultraviolet cutoff in the classical field theory. In addition, the equipartition property of classical fields leads to the Rayleigh-Jeans divergence in the continuum limit, a → 0. It should be remembered, however, that the scaling function allows us to evaluate the shear viscosity in the low momentum region by using the classical Yang-Mills field simulation at large lattice temperatures, where the classical field description is justified.
There are several sophisticated ways to circumvent these problems on the contribution of high-momentum contribution, such as renormalizing the theory by introducing counterterms in the classical field [19], applying the effective dynamics obtained by integrating hard particles [20], taking account of the explicit coupling of fields and particles [21], and switching to dynamics described by the kinetic theory [22]. Applying the classical field theory with quantum statistical nature [23] would be another interesting direction to be studied. This thermal distribution function can be obtained from the canonical partition function in Eq. (13), ( a dE a 1 ) exp(−ε E /T ) = √ 2ε E dε E dΩ exp(−ε E /T ), where Ω is the solid angle in the 3 dimensional color space. In Fig. 10, we show the distribution of ε E (x) in the configurations at t − t ini = 200 with g = 1 on the 32 3 lattice. This numerical distribution agrees with the thermal distribution whose temperature is 1.16(g 2 T = 1.16). From this agreement, we conclude that the configurations are sufficiently thermalized at t − t ini = 200. The power spectra in the very large g 2 T region also show interesting features. In Fig. 11, we show the normalized power spectrum, f ρ (ω) = ρ(ω)/T 2 at g = 1, 3, 10 (g 2 T = 1. 16, 10.6, 122) by solid curves. In this g 2 T 1 region, f ρ (ω) consists of the peak at ω = 0 and the bump at ω = (5 − 7). This structure is well described by the Fourier transform of the exponential decay and damped oscillations,

10
Dashed curves show the Fourier transform of F exp+DO (t), which almost completely agrees with the numerical results and thus is hard to discern in the figure.
In particular, the result obtained at the strongest coupling, g = 10 (g 2 T = 116), agrees well with the fit. This agreement of the functional form may be due to the dominance of the four-point interaction both in the large g 2 T region of the CYM theory and in the classical massless φ 4 theory.