QCD Plasma Instabilities and Isotropization

We study the coupled Wong Yang-Mills equations for both U(1) and SU(2) gauge groups and anisotropic particle momentum distributions numerically on a lattice. For weak fields with initial energy density much smaller than that of the particles we confirm the existence of plasma instabilities and exponential growth of the fields which has been discussed previously. Also, the SU(2) case is qualitatively similar to U(1), and we do find strong ``abelianization'' of the non-Abelian fields during the period of exponential growth. However, the effect nearly disappears when the fields are strong. This is because of the very rapid isotropization of the particle momenta by deflection in a strong field on time scales comparable to that for the development of Yang-Mills instabilities. This happens without entropy production which makes it very interesting for the phenomenology of high-energy heavy-ion collisions.

High-energy heavy-ion collisions release a large amount of partons from the wavefunctions of the colliding nuclei. Partons with large transverse momenta originate from high-Q 2 hard interactions which can be computed from perturbative QCD [1]. On the other hand, partons with "small" transverse momenta on the order of the so-called saturation momentum Q s (given by the square root of the total color charge density per unit rapidity and unit area in the incoming nuclei) are much more abundant if Q s ≫ Λ QCD and are better viewed as a classical non-Abelian field [2].
If the presence of the soft classical field is neglected, which amounts to assuming that Q s ∼ Λ QCD , the timeevolution of the hard partons after they come on-shell can be studied by means of the Boltzmann equation with a collision kernel, which is the so-called parton-cascade approach [3,4]. The collision kernel could be truncated at the level of elastic binary collisions (perhaps with a summation of time-like and space-like parton showers in the leading logarithmic approximation [3]); recently, an attempt to fully include 2 ↔ 3 processes beyond the relaxation time and leading-logarithmic approximations has also been made [5].
On the other hand, for large nuclei and at high energies the saturation scale Q s is expected to grow much larger than Λ QCD [2,6] and so the presence of the classical field can no longer be neglected. The "bottom-up scenario" [7] generalizes the parton cascade description of the time-evolution after the collision to include the soft classical modes, too. Soft gluon radiation is found to be the dominant process leading to equilibration [5,7,8] (see also papers by Wong in [4]).
Recently, it has been argued that collective processes due to the soft gauge field should be taken into account. Specifically, QCD plasma instabilities may develop due to anisotropic distributions of released hard partons [9] and modify the "bottom-up scenario" significantly [10]. The hard loop effective action for anisotropic hard modes was formulated in [11] and unstable soft modes were analyzed in [12]. Numerical studies of its static limit [13] revealed the interesting tendency of the non-Abelian gauge fields to "abelianize" during the stage of instability in the sense that locally commutators become much smaller than the fields themselves (see below). The "abelianization" has also been seen in solutions of the full non-linear hard loop effective action [14]. It is argued that, because of abelianization, non-Abelian effects should not cause instabilities to saturate; rather, similarly to the Abelian case, the fields should continue to grow until their energy density becomes comparable to that of the hard modes [13,14,15], i.e. until the growing fields begin to have a significant effect on the dynamics of the particles.
It is interesting to note the following difference between isotropization by propagation of particles in a strong random field versus that via scattering and gluon multiplication. Namely, in the absence of a collision kernel the entropy of any specific initial condition is conserved, while the standard parton cascade approach produces additional entropy [4,16]. An ensemble average over sufficiently random initial field configurations can nevertheless increase the entropy of the soft modes by a moderate (logarithmic) amount; this follows from the equivalence of the averaged classical field description to a Boltzmann equation to leading and subleading orders in the occupation number [17].
In heavy-ion collisions, it might not be necessary to achieve "true" thermalization in the sense of maximizing the entropy during the first few fm/c of the reaction; isotropization could be sufficient [15]. In fact, data [18] from RHIC indicate that the number of charged particles per participant in the final state is only ∼ 30% lower in central d+Au collisions than it is in central Au+Au. This perhaps indicates that the equilibration process expected to occur in Au+Au (but not in d+Au) does not produce a large amount of entropy [19]. Hence, the mechanism of isotropization of particles via strong fields could be very interesting for the phenomenology of heavy-ion collisions.
In this letter we solve the classical transport equation for hard gluons with non-Abelian color charge Q a in the collisionless approximation [20], It is coupled to the Yang-Mills equation where f (x, p, Q) denotes the one-particle phase space distribution function [20]. These equations were shown to reproduce the "hard thermal loop" effective action [20] near equilibrium. If fluctuations on top of the mean fields are not neglected, one obtains a collision term from their moments [21]. The same set of transport equations were also derived within the worldline formalism for the one loop effective action of QCD [22]; the emergence of classical transport from a quantum kinetic equation derived within the closed-time-path formalism was discussed in ref. [23]. For recent reviews of semi-classical transport theory for non-Abelian plasmas see ref. [24]. Furthermore, we refer to ref.
[25] for a study of particle production and propagation in Abelian fields, including back-reaction and collisions in the relaxation time approximation. The specific point of the present letter, however, is to study possible non-Abelian plasma instabilities due to anisotropic particle distributions [9,10,11,12,13,14,15]. We employ the test-particle method [26], replacing the continuous distribution f (x, p, Q) by a large number of test particles: where r i and p i are the coordinates of an individual test particle. This Ansatz leads to Wong's equations [20,27] for the i-th test particle [28]. The time evolution of the Yang-Mills field can be followed by the standard Hamiltonian method [29]. Numerical techniques to solve the classical field equations coupled to particles have been developed in ref. [30]. Our update algorithm is closely related to the one explained there which generalizes the Abelian version of the charge conservation method in particle simulations [31].
In the following, we assume that the fields only depend on time and on one spatial coordinate, x, which reduces the Yang-Mills equations to 1+1 dimensions. The hard modes represented by classical particles are allowed to propagate in three spatial dimensions. For simplicity, we also restrict ourselves to the case without expansion here; the more realistic case with longitudinal expansion [32] will be addressed in the future.
The initial anisotropic phase-space distribution of hard gluons is taken to be This represents a quasi-thermal distribution in two dimensions, with "temperature" p hard which now takes over the role of the saturation momentum mentioned above.
We have checked explicitly that no instability occurs when the particle distribution is taken to be isotropic. The initial field amplitudes are sampled from a Gaussian distribution with a width tuned to a given initial energy density. We solve the Yang-Mills equations in A 0 = 0 gauge and also set A = 0 (i.e. all gauge links =1) at time t = 0; the initial electric field is taken to be polarized in a random direction transverse to the x-axis. Gauss' law is then used to obtain the initial charge distribution. All results shown below were obtained using a lattice with N = 512 sites; we have checked the numerical accuracy by comparing to N = 256, 1024 lattices (for the same physical parameters) and by monitoring conservation of the total energy and of Gauss' law. The total energy was conserved to within 5 × 10 −4 (5 × 10 −3 ) over the course of the simulations for the weak (strong) field initial conditions, and the maximal violation of Gauss' law was 10 −9 for SU (2) and 10 −33 for U (1) (in lattice units).
Before coming to our results, we also comment on the occurence of "anomalous Cherenkov radiation". This corresponds to anomalous hard radiation from soft modes which may occur for simulations on a discrete lattice, as the dispersion relation of the fields may contain real space-like modes. For example, taking ω(k) = k for free fields in the continuum leads to the dispersion relation ω(k) = 2| sin(ka/2)|/a on a one-dimensional lattice (|k| < π/a). Consequently, hard field modes with k ∼ 1/a would then get populated on the lattice because their "mass" √ ω 2 − k 2 is imaginary. The situation could perhaps be improved by employing higher-order discretization schemes for the Yang-Mills action or by damping hard modes exponentially at t = 0. However, we have not done so at present. Our tests with different lattices do not indicate a significant dependence of either the growth or the saturation of the instability on the lattice spacing. Also, our solutions for isotropic particle momentum distributions appear to be stable when the number of test particles is taken to infinity. A possible physical reason for this observation could be that interactions among soft modes with k ≪ 1/a and with large occupation numbers, and interactions of those modes with the particles, which give the largest contribution to the total energy (see below), actually dominate. Nevertheless, this effect may deserve a more careful numerical study in the future.
We first show results for a large separation of initial particle and field energy densities which should qualitatively resemble the conditions studied in [13,14,15]. The results shown in Fig. 1 corresponds to a lattice of physical size L = 40 fm, a hard scale p hard = 10 GeV and a particle density of n/g 2 = 10 fm −3 . For the U (1) case we observe a rapid exponential growth of the magnetic field energy density starting at about t/L ≈ 0.1, turning into a slower growth at t/L ≈ 0.5; at this point the magnetic fields have grown sufficiently to affect the particles which visibly start loosing energy. The electric field grows less rapidly and equipartitioning is not achieved within the depicted time interval. This indicates that the field strengths are still too high for linear response to apply. In the non-Abelian case the growth of the magnetic field saturates earlier, and the electric field has comparable strength by the end of the simulation. Also, it appears that the saturation of the magnetic instability occurs before it has a noticeable effect on the particles since their energy density is nearly constant. Nevertheless, at a purely qualitative level the U (1) and SU (2) simulations are not vastly different, as anticipated in [13]. This is analyzed further in Fig. 2 of the rms average and the average of the relative size of the field commutator defined by [13] .
The behavior of φ rms is similar to that of the field energy density shown above. Initially,C is constant but then starts dropping exponentially when the magnetic instability sets in, indicating the partial "abelianization" of the fields [13,14]. The rate by whichC drops in the intermediate stage is roughly comparable to the growth rate of φ rms ; also, the abelianization appears to stop after C dropped by about one order of magnitude, at about the same time when the exponential growth of the fields and of φ rms saturates. Finally, in Fig. 3 we show the time evolution of the longitudinal and transverse components of the energymomentum tensor of the particles, i.e. the kinetic pressure. For both U (1) and SU (2) we observe a rapid growth of the longitudinal pressure, which is zero initially. Again, the rate is somewhat smaller for the non-Abelian case. The approach to "isotropization" of the kinetic pressure is clearly correlated to the stage of exponential growth of the soft fields seen before [15]. However, for both cases T xx remains significantly smaller than the transverse component for times < ∼ L. The initial conditions above were chosen such as to verify qualitatively the picture emerging in the hard loop approximation, where the field energy density is (and remains) much smaller than that of the particles and so the back-reaction can be neglected [13,14,15] color glass condensate model of high-energy collisions one does not expect such a strong separation of energy densities, however. Since our numerical solution includes the back-reaction of the fields on the particles, we study the situation with stronger fields next. Specifically, the simulations below were performed with the following set of physical parameters: length L = 10 fm, hard scale p hard = 1 GeV, particle density n/g 2 = 500 fm −3 and an initial field energy density of about 20 GeV/fm 3 .   4 shows the time evolution of the energy densities for these initial conditions. This case clearly differs from the weak-field limit shown before. Over the time interval shown, the electric field energy density is practically constant for both U (1) and SU (2). The Abelian magnetic field does exhibit a slow growth, draining some energy from the particle reservoir. For SU (2), however, after a short initial growth the magnetic field energy decreases again, to saturate pretty much at its initial value. Therefore, the particle energy density is also more or less constant over the depicted time interval.   5 confirms this observation via the φ rms observable: the initial growth saturates much earlier than before. Similarly, the average commutatorC stays constant for some time then drops by about a factor 5 (during the period where the magentic field drops !) and saturates at ≈ 10%, which is an order of magnitude larger than for the weak field case from Fig. 2. This indicates a much smaller degree of "abelianization" for strong fields.
Perhaps surprisingly, Fig. 6 nevertheless shows a very rapid isotropization of the kinetic pressure for both U (1) and SU (2) (note that t/L = 0.1 corresponds to t=1 fm in physical units for this lattice). Moreover, the degree of isotropization is much higher, i.e. the transverse and longitudinal pressures are closer than in Fig. 3. The very fast and nearly complete isotropization is, of course, the reason why field instabilities can not be sustained over a significant period of time in this case. It is caused by the bending of the particle trajectories in the strong field which is very different from conventional parton cascade transport with small-angle perturbative scattering (and no field). The random initial fields then cause a rapid isotropization of the particle momenta via eq. (5). Note that this does not require hard modes in the fields, which indeed would violate the assumed separation of momentum scales, but large field amplitudes.
In summary, we have studied instabilities in the coupled Wong Yang-Mills equations for strongly anisotropic initial particle momentum distributions. For both U (1) and SU (2) gauge groups we do observe a period of exponential growth of the fields when their initial energy density is far less than that of the hard modes (particles). This, in turn, leads to partial isotropization of the particle momentum distributions and of the kinetic pressure. Although we find somewhat smaller field growth and isotropization rates for the non-Abelian case, we nevertheless qualitatively confirm the picture developed in [13,14] in that the non-Abelian fields "abelianize" efficiently during the period of exponential growth.
For large initial field amplitudes, corresponding to a smaller ratio of initial particle to field energy densities, our results are qualitatively different. We observe a very rapid isotropization of the particle momentum distributions which is due to bending of their trajectories in the strong fields on a time-scale that is relevant for the physics of high-energy collisions. This, however, inhibits the development of instabilities of the Yang-Mills fields. Nevertheless, these results, too, suggest that the presence of the strong non-Abelian fields should be taken into account to understand the process of isotropization in the early stages of high-energy collisions.
We thank F. Gelis, M. Gyulassy, A. Mueller, R. Pisarski and M. Strickland for helpful discussions and J. Lenaghan and M. Strickland for drawing our attention to this subject.
Note added: After this manuscript was submitted for publication, a paper appeared [33] which presents an analytical discussion of a possible effective potential for anisotropic QCD plasmas beyond the hard loop approximation. Also, a modified "bottom-up" scenario for gluon thermalization in high-energy heavy-ion collisions [34] and 3 + 1d simulations of the full hard-loop effective theory [35] appeared.