Entropy production and isotropization in Yang-Mills theory with use of quantum distribution function

We investigate thermalization process in relativistic heavy ion collisions in terms of the Husimi-Wehrl (HW) entropy defined with the Husimi function, a quantum distribution function in a phase space. We calculate the semiclassical time evolution of the HW entropy in Yang-Mills field theory with the phenomenological initial field configuration known as the McLerran-Venugopalan model in a non-expanding geometry, which has instabilty triggered by initial field fluctuations. HW-entropy production implies the thermalization of the system and it reflects the underlying dynamics such as chaoticity and instability. By comparing the production rate with the Kolmogorov-Sina\"i rate, we find that the HW entropy production rate is significantly larger than that expected from chaoticity. We also show that the HW entropy is finally saturated when the system reaches a quasi-stationary state. The saturation time of the HW entropy is comparable with that of pressure isotropization, which is around $1$ fm/c in the present calculation in the non-expanding geometry.


Introduction
A new form of matter consisting of deconfined quarks and gluons is formed in high-energy heavy-ion collisions at Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) [1,2,3,4,5]. The created matter is opaque for colored particles, shows hydrodynamical behavior and collectivity of quarks, and finally decays into hadrons. Then it is considered to be a quark gluon plasma (QGP) [6,7]. While quantitative studies on the QGP properties are in progress using hydrodynamical models combined with jet and hadronic transport, its formation process is not yet clear. For example, the early thermalization problem remains as one of the serious problems in highenergy heavy-ion collisions [8]. Hydrodynamical-model analyses suggest that the created matter becomes close to local equilibrium at τ th = 0.6-1.0 fm/c after the contact, and this thermalization time is significantly shorter than the perturbative QCD estimate [9,10].
In tackling the early thermalization problem, the classical Yang-Mills field plays an important role. The created matter in the initial stage is described well by the classical Yang-Mills field, and is often called "glasma" [11]. In the glasma, both the color-electronic and -magnetic fields are parallel to the collision axis, the pressure is anisotropic, and the anisotropy leads to instabilities triggered by initial field fluctuations [12,13,14,15,16,17,18,19,20]. Fluctuations of classical fields may be regarded as particles, then the glasma instability is expected to produce many particles and cause early thermalization. Actually, recent studies [21,22,23,24] have successfully shown early-time isotropization of the pressure required by the hydrodynamical model analyses. For the detailed understanding of the thermalization process, however, we need to evaluate the entropy of the system, which is a very important key concept that characterises thermalization. In Ref. [25], the required amount of entropy produced in the glasma is estimated to be 3000 per rapidity or the 55 % of the total entropy. Nevertheless, in many of previous works [12,13,14,15,16,17,18,21,22,23,24,26,27], the entropy production itself has not been discussed, and the relation between the isotropization and thermalization remains unclear.
Entropy production in the classical Yang-Mills field theory has been discussed based on the Kolmogorov-Sinaï entropy production rate (KS rate) [28,29,30,31] and the Husimi-Wehrl entropy [32,33]: Entropy of classical systems is obtained by the Wehrl entropy [34,35], S W = −Trf log f , where f is the phase space distribution function and Tr denotes the integral over the phase space. In quantum systems, the Wigner function [36,37,38,39] (f W ) is a candidate of the distribution function, since it is defined through a mere Weyl transformation [40] of the density matrix and should contain full information equivalent to the density matrix. However the Wigner function is not appropriate as the distribution function to discuss entropy production; it is not semi-positive definite and cannot be regarded as the phase space probability distribution. In addition, even if f W ≥ 0 is satisfied everywhere, the Wehrl entropy does not increase in the semiclassical time evolution due to the Liouville theorem.
One possible solution for the phase space distribution function in calculating the Wehrl entropy is the Husimi function [41] (f H ). The Husimi function is obtained from the Wigner function by smearing in the phase space within the allowance of the uncertainty principle, and it is shown to be semi-positive definite. In fact, Husimi function is an expectation value of the density matrix with respect to the wave packet with the minimal uncertainty, which is nothing but a coherent state [42,43]. We call the Wehrl entropy defined with the Husimi function the Husimi-Wherl (HW) entropy [32,33,34,35,44]. The HW entropy is shown to be approximately the same as the von Neumann entropy at high temperatures [44]. In inverse harmonic oscillators, the HW entropy is found to increase in time and the growth rate agrees with the KS rate, the sum of the positive Lyapunov exponents [44]. The increase of the HW entropy implies information loss caused by instabilities and/or chaoticities combined with the coarse-graining in the phase space, and it is expected to play a crucial role in thermalization.
We can obtain the HW entropy in field theories by regarding the field strength and its canonical conjugate momentum as the phase space variables. In Ref. [33], the present authors have calculated the semiclassical time evolution of the HW entropy of the classical Yang-Mills fields with a random initial condition, and have confirmed that the HW entropy growth rate is consistent with the KS rate [30]. This agreement suggests the entropy production is caused by the chaoticity of the classical Yang-Mills fields, since the KS rate characterizes the chaoticity of the system.
In this article, we discuss entropy production in Yang-Mills field theory starting from the glasmalike configuration given by the McLerran-Venugapalan(MV) model [45,46] in the non-expanding geometry [47] based on the framework developed in [32,33]. Quantum fluctuations are incorporated around the initial glasma-like field configuration, and we compare the time scales of the entropy production with that of other quantities such as the pressure isotropization and the equilibration of the local energy distribution.
This paper is organized as follows. In Sec. 2, we introduce the quantum distribution functions and entropy in field theories as well as the initial condition in the MV model in the non-expanding geometry. In Sec. 3, we explain the numerical method to calculate the semiclassical time evolution of the HW entropy and pressure. We show the results in the Sec. 4. Section 5 is devoted to the summary of our work.

Quantum distribution functions and entropy in Yang-Mills theory
The Husimi-Wehrl entropy of the Yang-Mills field is obtained as a natural extension of that in quantum mechanics by regarding (A(x), E(x)) as canonical variables. We define the Wigner and Husimi functions on the lattice, as a straightforward extension of those in quantum mechanics [38,39]. The semiclassical time-evolution of the Wigner function is given by the classical equation of motion (see Eq.(6) below), then we can obtain the Husimi function from thus constructed Wigner function at each time.
In the SU(N c ) Yang-Mills field theory on a L 3 lattice in the temporal gauge, the Hamiltonian in the non-compact formalism is given by where is the field strength tensor, and N D = 3L 3 (N 2 c − 1) is the total degrees of freedom (DOF). We take the dimensionless gauge field A and conjugate momentum E and spacetime variables x normalized by the lattice spacing a throughout this article. Then the Wigner function f W [A, E; t] is defined by a Weyl transform of the density matrixρ as where A · E = i,a,x A a i (x)E a i (x) denotes the inner product. It should be noted that the coupling constant g appears in the denominator in the integral measure, since g is included in the definitions of A and E. The expectation value of a physical quantity X is given by integrating the product of f W [A, E; t] and X in the Weyl representation denoted by X W [A, E] as; where DΓ = DA DE/(2π g 2 ) ND . For instance, let X be transverse (longitudinal) pressure P T,L . The pressure is given by the diagonal part of the energy-momentum tensor T µν . The expectation value is then given by where the E a ⊥ (B a ⊥ ) is the transverse component of the color electric (magnetic) field, and E a The color magnetic field is defined as B ai = − ǫ ijk 2 F a jk , and the ǫ ijk is a completely antisymmetric (Levi-Civita) tensor (ǫ 123 = 1). The time evolution of the Wigner function is derived from the von Neumann equation, In the semiclassical approximation in which we ignore O( 2 ) terms, f W is found to be constant along the classical trajectory given by the classical equation of motion (EOM) [48], The Husimi function is defined as the smeared Wigner function with the minimal Gaussian packet, where ∆ = a∆ phys is the dimensionless parameter corresponding to the Gaussian-smearing range.
It should be noted that the Husimi function is also obtained as the expectation value of the density matrix in the coherent state as in quantum mechanics [42]. Then the Husimi function is semi-positive definite, f H [A, E; t] ≥ 0, while the Wigner function is not. We finally define the Husimi-Wehrl entropy as the Boltzmann's entropy or the Wehrl's classical entropy [34] by adopting the Husimi function for the phase space distribution, The HW entropy is gauge invariant, and the semiclassical time evolution does not break the gauge invariance as shown in Appendix A.

Initial condition
We consider two nuclei moving at the velocity of light along the z axis. These nuclei collide at time t = 0 and z = 0, and glasma is formed between the two nuclei. In the framework of the color glass condensate (CGC), the gluons with small Bjorken x are described by the classical field and those with large Bjorken x and quarks are regarded as color sources. The color-source distribution is assumed to be Gaussian in the McLerran-Venugopalan (MV) model [45,46].
We adopt a glasma-like initial condition which mimics the MV model in the non-expanding geometry [47]. As in the MV model, we generate the Gaussian random color sources for a target nucleus ρ (t) and a projectile ρ (p) , where x ⊥ ≡ (x, y) and a, b are the color indices. On the lattice, the delta function δ (2) (x ⊥ − y ⊥ ) is replaced by the Kronecker delta δ x⊥,y⊥ /a 2 , and Eq.
where µ = aµ phys and ρ (i)a is given in the lattice unit. Gauge fields are given by α , which are created by the color sources; Gauge fields, electric fields and magnetic fields are then given by The above gauge fields are classical and uniform in the z direction, and there is no quantum fluctuations for a given source. In order to make the initial Wigner function f W [A, E, t = 0] taking into account quantum fluctuations, the uncertainty relation between A and E, the initial Wigner function is set to be a glasma-like field configuration with a Gaussian fluctuation around it. With A MV and E MV being the solutions of Eqs. (13) and (14), the initial Wigner function is obtained as where ω = aω phys is the parameter of the Gaussian width.

Physical scale
We have two dimensionful parameters, g 2 µ phys and ω phys , in the initial condition, one dimensionful parameter, ∆ phys , in the calculation of the HW entropy, and one dimensionless parameter, g 2 , in addition to the lattice spacing a. The factor g 2 appears from the field redefinition, gA → A and gE → E, then the uncertainty relation is modified as (∆A) 2 (∆E) 2 ≥ ( g 2 /2) 2 for each component of A and E. This relation is consistent with the classical field dominance in the weak coupling regime. Since the saturation scale Q s is the fundamental scale in the color glass condensate, we take µ phys ≃ Q s and ω phys ≃ Q s . We now set the physical scale. We consider heavy-ion collisions at RHIC and LHC energies, then g = 1(α s = 0.15) and Q s ≃ 2 GeV may be reasonable. We also assume that the total lattice area in the xy plane is equal to the transverse area of the colliding nuclei. Then the parameters are fixed as From these equations, we get The lattice spacing a is inversely proportional to the lattice size L.

Numerical methods
It is not an easy task to perform numerical calculation of the HW entropy, Eq. (10), especially in field theories. We need 2N D dimensional integral of a function f H , which additionally requires 2N D dimensional integral to obtain, where N D is very large in field theories. The logarithmic term − log f H takes a large value when f H is small and the integrand exhibits an acute peak. The Monte-Carlo method is then effective and necessary for a large-dimensional integral.
We have developed numerical methods to calculate the time evolution of the Husimi-Wehrl entropy in semiclassical approximation in quantum mechanical systems [32] and the Yang-Mills field theory [33]. In this section, we recapitulate our formalism. We introduce two methods based on the test particle (TP) method to calculate the HW entropy, which were applied to Yang-Mills field theory in Ref. [33]. The test particle method is applied also to calculate other physical quantities such as pressure.

Test particle method and Husimi-Wehrl entropy
In the TP method, we express the Wigner function by a sum of the delta functions, where N TP is the total number of the test particles, the number of the delta functions used to express the Wigner function. The variables It is noteworthy that the Husimi function here is a smooth function in contrast to the corresponding Wigner function in Eq. (20). With the Wigner function Eq. (21), the HW entropy in the test-particle method Eq. (10) is now obtained as, Note here that the integral over (A, E) has a support only around the positions of the test particles (A α (t), E α (t)) due to the Gaussian function for each α, and we can effectively perform the Monte-Carlo integration. We generate random numbers (A α,k , E α,k ) (k = 1, · · · , N MC ) with zero mean and standard deviations of ( g 2 /2∆, g 2 ∆/2), with N MC being the total number of Monte-Carlo samples. Then we obtain the HW entropy as shown in the second line of Eq. (22). The width parameter ∆ needed to define the Husimi function f H is set to be ∆/ω = 1. At present, ∆ is treated merely as an input parameter. We have checked the dependence of results on ∆ and confirm that main conclusions remain unchanged.
The TP method has a following problem. In the case where α = β in Eq. (22), the Husimi function, the argument of the logarithm, tends to take a large value, which generally leads to an underestimate of the HW entropy. Since this underestimate arises from the f W sampled with a finite number of delta functions (test particles), the HW entropy in the TP method is essentially underestimated, though this artifact vanishes when N T P → ∞. In order to evade the problem, we also introduce a parallel test particle (pTP) method, where we prepare independent sets of test particles (A α , E α ) and (A β , E β ) for in and out of the logarithm in Eq. (22). In the pTP method, the HW entropy tends to be overestimated. The phase space distance of test particles grows exponentially in chaotic or unstable systems, then we may not have any test particle (A β , E β ) inside the logarithm in the vicinity of the test particle (A α , E α ) prepared outside the logarithm. In this case, the argument of the logarithm becomes very small, and − log f H is overestimated. While both the TP and pTP methods have problems stemming from the formalism, the results should converge at large N TP from below and above in the TP and pTP methods, respectively, and the converged value of S HW exists between the TP and pTP results at a finite N TP .

Product ansatz
While the TP and pTP methods can be, in principle, applied to the field theory on the lattice, the DOF is large and numerical-cost is demanding. For example, we need to adopt very large number of test particles, N TP , to make the Monte-Carlo integration converge. Since the Husimi function is equivalent to the expectation value of the density matrix in the coherent state, it has a value in the range of 0 ≤ f H ≤ 1. The Gaussian Eq. (9) take the maximal value 2 ND , then the required number of test particles is N TP > 2 ND in order to respect the f H range. Thus we need to invoke some approximation scheme in practical calculations.
We here adopt a product ansatz to avoid this difficulty. In the ansatz, we assume that the total Husimi function is given as a product of that for each degree of freedom, where I = (i, a) denotes the direction (i = x, y, z) and color indices (a = 1, 2, 3), and f is obtained by integrating out other degrees of freedom than I. By substituting this ansatz into Eq. (10), we obtain the HW entropy as a sum of the HW entropy for each degree of freedom; It is found that the HW entropy obtained with product ansatz is found to overestimate the entropy by 10-20 % in a few-dimensional quantum mechanical system [33]. Secondly, the maximum value of the HW entropy in the TP method is shifted with the product ansatz, while the minimum value remains unchanged. For a one-dimensional case, there is a minimum of S HW = 1 [35,49]. When the Wigner function is a Gaussian, f W (A, E) = G(A, E; ω), the Husimi function is also a Gaus- , and the HW entropy is found to be S HW = N D (1 − log[2 √ ∆ω/(∆ + ω)]) ≥ N D . The equality holds when we take ∆ = ω. The HW entropy will have an upper bound in the TP method, when all the test particles are separated from each other, and we find S HW ≤ N D + log(N TP /2 ND ) [32]. In the TP method with the product ansatz, the HW entropy for each DOF has the above upper bound for N D = 1, S (I) HW ≤ 1 + log(N TP /2). Thus the upper bound of the HW entropy with the product ansatz becomes larger than that without the ansatz, Thirdly, the HW entropy in the product ansatz S (PA) HW is not gauge invariant. Nevertheless we might expect that the gauge dependence does not cause serious problems in entropy production because gauge degrees of freedom dose not significantly contribute to chaoticity and instability [19,31], and that the production rate of the HW entropy from random initial condition in the product ansatz agrees with the gauge invariant KS rate [33].

Vacuum subtraction
When we calculate observables in field theories, it is generally necessary to subtract vacuum expectation values. It also applies to the present semiclassical treatment. Let X be a physical quantity and X MV be the expectation value calculated by using the Winger function, as given in Eqs. (3) and (20). When we calculate an expectation value of X, we subtract the vacuum contribution X vac arising from quantum fluctuations. We have evaluated the vacuum expectation value by using the fluctuation part of (A, E), with A α (0) = A MV + δA α (0) and E α (0) = E MV + δE α (0). For example, in the case of X = P T,L , the expectation values are given by

Results
We shall now discuss the numerical results of the time evolution of the HW entropy and the pressure based on the numerical methods explained in Sec. 3. We mainly show the results on the 64 3 lattice, and also show some of the results on the 16 3 and 32 3 lattices for comparison. The 64 3 lattice may be a reasonable choice to discuss heavy-ion collisions at RHIC and LHC based on the classical Yang-Mills fields. The classical Yang-Mills field theory is a low-energy effective theory and has a ultraviolet cut off. At L = 64, the lattice spacing is a ≃ 2Q −1 s , which corresponds to the diameter of one color flux tube.

Husimi-Wehrl entropy production
In Fig. 1, we show the time evolution of the HW entropy on the 32 3 and 64 3 lattices obtained by the TP and pTP methods with the product ansatz. The HW entropy per DOF starts from the minimum value, S HW /N D = 1, then increases rapidly and almost linearly until g 2 µt = 3 at almost a common rate on the 32 3 and 64 3 lattices, and shows slow increase in the later stage. In the later stage, e.g. g 2 µt = 10, the HW entropy takes a smaller value on the 64 3 lattice. The pTP method gives the upper bound of the HW entropy and the TP methods gives the lower bounds, then we can guess that the converged value in the limit of N TP → ∞ exists between the results of the two methods as discussed in Ref. [33].  The growth rate of the HW entropy in the linearly increasing stage may be characterised by the Kolmogorov-Sinaï (KS) rate, which is the sum of positive Lyapunov exponents and reflects the underlying dynamics. Due to the scale invariance of classical Yang-Mills, the KS rate scales as λ KS /L 3 = c KS × ε 1/4 [30], where ε = H /L 3 is the energy density. Then the HW entropy is expected to increase as We While the local KS rate obtained from the second derivative of the Hamiltonian at initial time is sensitive to the gluon-field configuration itself, the intermediate KS rate represents the intrinsic property of the chaotic system that dose not depend on initial conditions. In Fig. 1, we compare the numerically obtained HW entropy and that expected from the KS rates. The black straight lines in Fig. 1 show the entropy increase expected from the local and intermediate KS rates. In the present calculation on the 64 3 lattice, the total energy (energy density) amounts to H = 6.5 × 10 5 (ε = 2.48), then the slope from the local (intermediate) KS rates, c LLE(ILE) g 2 µ , is evaluated to be 0.14(0.074). As seen in Fig. 1, the growth rate of the HW entropy in the early time is around dS HW /d(g 2 µt)/N D ≃ 0.14, which is close to the local KS rate and significantly larger than the intermediate KS rate. This comparison implies that we cannot explain the entropy production from the MV model initial condition only by intrinsic chaoticity, and that some instability may be the trigger of the entropy production. In fact, the initial field configuration of the MV model has strong instabilities and the HW entropy is considered to saturate even without showing the intermediate KS rate.
The HW entropy production rate per degrees of freedom in the early stage is almost independent of the lattice size. In addition to 32 3 and 64 3 lattices shown in Fig. 1, similar production rate is found on smaller lattices, 4 3 , 8 3 and 16 3 . At least, this lattice-size independence dose not come from the chaoticity of the system because the KS rates depend on the lattice size. The energy density on the 32 3 lattice is ε = 3.53, and the slopes from the local and intermediate KS rates are evaluated as 0.08 and 0.04, respectively, which are smaller than the KS rates on the 64 3 lattice. This fact suggests that another possible mechanism exists to create the HW entropy such as the initial instability.
The amount of the produced entropy on the 64 3 lattice is ∆S/N D ≃ 0.4 and may be in the same order of the expected entropy production. The longitudinal thickness of glasma at the initial stage should be in the order of Q −1 s , and the present calculation in the nonexpanding geometry corresponds to a very thick nuclei, aL = 120 Q −1 s . The produced entropy per unit rapidity for color SU(3) is expected to be This value is around half of the expected entropy, ∆S/∆Y ≃ 4500 [25], but several systematic uncertainties in the present setup could easily account for a factor of two. Calculation of the entropy production in an expanding geometry is desired.

Classical equilibration
The present calculation shows that the HW entropy approximately saturates at g 2 µt ≃ 7(5) on the 32 3 (64 3 ) lattice, then some kind of quasi-stationary state is formed. In Fig. 2, we show the distribution of the electric and magnetic local energies, In the thermal equilibrium in the classical regime, the distribution of (A, E) would be described by the Boltzmann distribution, For the electric energy distribution, we can rewrite the measure as d 3 E = √ H E dH E dΩ with dΩ being the solid angle in the color space, and the distribution function can be given as √ H E exp(−H E /T ). Actually, the electric energy distribution in the later stage is described well by this distribution except for the high energy region as shown by the solid line in Fig. 2. The magnetic energy distribution is also found to follow the same function but with a different temperature. Similar Boltzmann distribution of the energy is found in Ref. [29]. Thus the saturation of the HW entropy seems to be related to the quasi-statinary state, where approximate equilibrium is reached among the electric energies and among the magnetic energies but with a different temperature.
The saturation time and saturated value of the HW entropy per DOF decrease with increasing lattice size, as shown in Fig. 1. It should be noted that the above quasi-statinary state is, however, different from the true equilibrium of gluons: In addition that the electric and magnetic temperatures are different, the long-term evolution with the classical Yang-Mills equation does not reach the Bose-Einstein distribution of the high-momentum modes but reach the classical statistical distribution.
Since the classical statistical distribution in field theories does not have a well-defined continuum limit, it is reasonable to find the lattice size dependence of the saturation time and saturated value of the HW entropy.

Isotropization of pressure
In Fig. 3, we show the time evolution of the ratio of the pressure to the energy density ratio in the longitudinal and transverse directions, P L,T /ε, on the 16 3 , 32 3 and 64 3 lattices. Because the energy-momentum tensor is traceless, the relation 2P L /ε + P T /ε = 1 is satisfied. While the classical configuration of the MV model has the completely anisotropic pressure P L = −P T at initial time, the quantum fluctuations modifies this relation and the initial value P L /ε (P T /ε) is not equal to 1.0 (−1.0).
The isotropization of the pressure can be found to occur in Fig. 3. The lattice size dependence of the isotropization time is strong in the smaller lattices, L < 32, and For larger lattices (L ≥ 32), the isotropization time almost converges g 2 µt ≃ 10, as seen from the L = 32 and L = 64 results, This isotropization time roughly agrees with the time of the HW entropy saturation. It also happens to agree with the isotropization time obtained in the expanding geometry with fluctuation effects from the finite coupling [21].

Summary and conclusion
The aim of this paper is to understand the thermalization in the relativistic heavy ion collisions by focusing on the entropy production. We calculate the Husimi-Wherl (HW) entropy in the Yang-Mills field theory with the phenomenological initial condition given by the McLerran-Venugopalan (MV) model but in the non-expanding geometry. The HW entropy constructed from the Husimi function plays an important role in thermodynamics of quantum systems and its production implies the thermalization of the system. We calculate the semiclassical time evolution of the Wigner function by solving the classical equation of motion which keeps the gauge invariance of the HW entropy. In actual calculations, we use the product ansatz to reduce the numerical cost at the cost of breaking the gauge invariance of the HW entropy. Nevertheless the HW entropy in the product ansatz agrees with the result without the product ansatz within 10-20% in a few-dimensional quantum mechanical system [33] and production rate of the HW entropy in the product ansatz agrees with the Kolmogolov-Sinaï (KS) rate.
We have found that the HW entropy increases linearly in early time and saturates at later times. The growth rate of the HW entropy is independent of the lattice size and is significantly larger than the intermediate KS rate, defined as the sum of the positive intermediate Lyapunov exponents. This implies that we cannot explain the entropy production from MV initial condition only by intrinsic chaoticity. It also suggests that the large amount of the entropy may be produced by the initial instability. When the HW entropy saturates, the electric and magnetic local energy distributions reach the classical statistical equilibrium except for the high energy regions. The saturation time agrees with the equilibrium time of the local energy distribution and the isotropization time of the pressure, which suggests the thermalization of the gluon field is realized in the sense of the HW entropy production. The saturation time is around 1 fm/c. In order to reach more quantitative and realistic conclusions, the evaluation the HW entropy in the expanding geometry is desired, which is under progress.