Turbulence in nonabelian gauge theory
Kolmogorov wave turbulence plays an important role for the thermalization process following plasma instabilities in nonabelian gauge theories. We show that classical-statistical simulations in SU(2) gauge theory indicate a Kolmogorov scaling exponent known from scalar models. In the range of validity of resummed perturbation theory this result is shown to agree with analytical estimates. We study the effect of classical-statistical versus quantum corrections and demonstrate that the latter lead to the absence of turbulence in the far ultraviolet.
Heavy-ion collision experiments explore strongly interacting matter starting from a transient far-from-equilibrium initial state. Available data from the Relativistic Heavy Ion Collider reveals remarkable, unexpected properties such as rapid apparent thermalization with robust collective phenomena. The explanation of these findings from quantum chromodynamics (QCD) provides a challenge for theory. There are theoretical indications that nonequilibrium dynamics related to plasma instabilities could be crucial for the understanding of the plasma’s temporal evolution [1, 2, 3, 4]. These instabilities lead to exponential growth of occupation numbers in long wavelength modes on time scales much shorter than the asymptotic thermal equilibration time. After the fast initial period of exponential growth the dynamics slows down considerably. The system is still far from equilibrium at this stage and the subsequent slow evolution may be characterized by Kolmogorov wave turbulence .
In the turbulent range characteristic properties such as values for scaling exponents can be insensitive to the details of the underlying microscopic theory. The associated universality may lead to quantitative agreements for very different theories or energy scales. Wave turbulence has been studied in detail in weakly-coupled scalar theories in the context of early-universe reheating [6, 7]. Though the underlying mechanisms are very different for QCD and for scalar inflaton dynamics, the subsequent evolution after an instability follows similar patterns: There a tachyonic or parametric resonance instability leads to the exponential growth of occupation numbers, followed again by a slow cascade [8, 9, 10, 11, 12]. Wave turbulence in scalar theories is described by nonthermal distributions such as with for cubic self-interactions in the presence of a field expectation value, or for quartic interactions () describing particle (energy) cascades . These distributions are distinct from a thermal high-temperature distribution . The nonthermal power-law behavior can be obtained from stationary solutions of perturbative Boltzmann equations neglecting quantum corrections . A new class of nonperturbative scaling solutions, with strongly enhanced low-momentum fluctuations, has recently been found . A very detailed picture of the different scaling solutions emerges using nonequilibrium renormalization group techniques .
In contrast to the detailed picture for scalar quantum theories, our knowledge about gauge theories is still in its infancies. Even at sufficiently high momenta, where perturbative estimates are expected to be applicable, conflicting statements can be found about the nature of power cascades following plasma instabilities [14, 15]. It has been noted that the value for the exponents characterizing scaling behavior for gauge theories could deviate from those known for scalar field dynamics. For instance, a detailed numerical study in Ref.  for SU(2) gauge theory in Coulomb gauge suggests the value two for the scaling exponent, which is not fully understood perturbatively and seems to have no corresponding value in scalar models. The estimate employs a separation of scales between suitably defined ’soft’ and ’hard’ momenta for sufficiently small characteristic running gauge coupling. The respective hard-loop effective theory of soft excitations is a nonabelian version of the linearized Vlasov equations of traditional plasma physics, which are based on collisionless kinetic theory for hard particles coupled to a soft classical field . Vlasov equations can also be derived starting from classical-statistical field theory . Therefore, they may also be considered as an approximation of the classical-statistical field theory limit of the respective quantum gauge theory [3, 4]. Classical-statistical lattice gauge theory provides a quantitative description in the presence of sufficiently large energy density or occupation numbers per mode.
In this work we show that classical-statistical simulations following plasma instabilities in SU(2) gauge theory indicate a Kolmogorov scaling exponent known from scalar models (Sec. 2). More precisely, we demonstrate that the numerical simulation is consistent with characterizing particle cascades for scattering associated to the phenomenon of weak wave turbulence. We show that this matches resummed perturbative estimates based on the two-particle irreducible (2PI) effective action (Sec. 3). In this perturbative regime we study the effect of quantum corrections and demonstrate that turbulence is absent at sufficiently short distances.
2 Classical-statistical lattice gauge theory
Collisions of heavy nuclei leave behind a plasma of quarks and mostly gluons in a nearly flat region of space because of Lorentz contraction along the beam- or -axis. If the quarks and gluons streamed freely into the surrounding space, then the distribution of particles would become locally highly anisotropic. In this case the stress tensor quickly acquires an oblate form with . The fastest growing mode of plasma instabilities then has its wave vector along the normal direction and generates a prolate contribution to the stress, which pushes the system towards greater isotropy [1, 2, 3, 4].
We describe this physics using classical-statistical Yang-Mills theory in Minkowski space-time following closely Ref. , to which we refer for further technical details. It employs the Wilsonian lattice action for SU(N) gauge theory in Minkowski space-time:
with and spatial Lorentz indices . It is given in terms of the plaquette variable , where . Here is the parallel transporter associated with the link from the neighboring lattice point to the point in the direction of the lattice axis . The definitions and contain the lattice parameter , where denotes the spatial and the temporal lattice spacings, and we will consider .
We specify to SU(2) as the gauge group. The dynamics is solved in temporal axial gauge. Varying the action (1) w.r.t. the spatial link variables yields the classical equations of motion. Variation w.r.t. to a temporal link gives the Gauss constraint. The coupling can be scaled out of the equations of motion and we will set for the simulations. We define the gauge fields as
where are the Pauli matrices with color index . Correlation functions are obtained by repeated numerical integration of the classical lattice equations of motion and Monte Carlo sampling of initial conditions. To be specific, we consider the extreme anisotropy case described by an effectively -like initial Gaussian distribution. The initial time derivatives are set to zero, which implements the Gauss constraint at all times. Shown results are from computations on lattices, where gauge fixing to Coulomb gauge is employed using a stochastic overrelaxation algorithm . We comment on the observed stability of the results on and lattices below.
Fig. 1 shows the equal-time correlation function
as a function of momentum for various instances of time in units of the energy density . The direction of the gauge fields is in the transverse -plane and the momenta are parallel to the axis. This correlation function is typically associated to an occupation number distribution divided by frequency. Because of the oblate initial conditions the distribution along the axis is characterized by small values at early times. Exponential growth due to plasma instabilities leads to a rapid population of long-wavelength modes, while fluctuation effects induced by the growth in the lower momentum modes lead to amplification in a broad range of momenta. The exponential growth stops around , and we concentrate on the subsequent evolution. The behavior at earlier times is analyzed in detail in Ref. .
After the exponential growth period the dynamics slows down considerably. As time proceeds, more and more ultraviolet modes approach a power-law described by
with the ”occupation number exponent” . Fig. 2 shows the simulation results at time for transverse (square) and longitudinal momenta (circle). Given that longitudinal and transverse momentum modes were populated very differently at initial time, the differences reflect the comparably small remaining anisotropy at this time. The solid line represents a power-law fit (4) with
which matches the data in a significant momentum range rather well.
In order to obtain the fit value for from the simulations, we first employ a least-square fit to the transverse and longitudinal results separately at a given time. Subsequently, we average the results for the transverse and longitudinal momenta. For the longitudinal momenta we choose the fit range . For the transverse-momentum results, which display power-law behavior in a larger domain at earlier times, we use the range . Convergence in the infrared proceeds rather slowly. However, changing the lower bound for the fit range from down to leads only to changes in which are comparable to the statistical error. Fig. 3 displays the fit values for as a function of time. The comparably large error for comes from the insufficient isotropization at these times. Subsequently, the values for change very little as a function of time and the errors become small. The dashed lines correspond to the results for perturbative Kolmogorov exponents for energy () and particle () cascades, which is explained in Sec. 3. The late-time behavior for is consistent with a slow approach to a classical thermal equilibrium distribution with exponent one (solid line). A remarkable agreement between simulation results and is observed in the quasi-stationary period around . We obtain the best fit value (5) by taking the time-average in this quasi-stationary period.
We emphasize that the errors quoted are purely statistical. The different dynamics for longitudinal and transverse modes leads to the fact that the longitudinal distribution corresponds to somewhat larger estimates for than what is inferred for transverse modes at not too late times. If this difference is employed to estimate the systematic error for , the error will be less than ten per cent.
Our results seem to be incompatible with the value two for reported in Refs. [14, 18] based on Vlasov equations. Unfortunately, the origin of that value is not fully understood so far. In particular, a perturbative analysis in terms of Kolmogorov wave turbulence similar to the one in Sec. 3 seems not to be able to explain it. There is still no consensus about the underlying reasons for the different numerical estimates of using Vlasov and classical-statistical approaches, respectively. It has been pointed out that in the Vlasov treatment the hard modes represent static sources, which do not isotropize. In contrast, total energy is conserved in the classical-statistical lattice gauge theory simulation and all modes isotropize at sufficiently late times. However, for weak coupling and assuming a sufficiently large separation of scales both approaches are expected to agree. One could conclude that simulations on very much larger lattices to achieve the clear separation of scales underlying the Vlasov approach might be required to compare the two. We verified that simulations on lattices as big as , with much less statistics, show no indication for quantitatively relevant changes of the value. We found that already simulations on lattices accurately reproduce (5), which indicates that the results are rather robust.
Clearly, none of these approximations are sufficient to quantitatively address the asymptotic late-time behavior, which should finally be characterized by a Bose-Einstein distribution for the gluons. While for low momenta the employed classical-statistical simulations are expected to give an accurate description for sufficiently high energy densities or occupation numbers, the high-momentum behavior will be altered by quantum corrections. In the following, we address aspects of this question using resummed perturbation theory at high momenta.
3 Resummed perturbative quantum theory
For quantum fields one can define two independent two-point correlation functions out of equilibrium, which may be associated to the anti-commutator and the commutator
respectively. Loosely speaking, the spectral function determines which states are available, while the statistical propagator contains the information about how often a state is occupied. The spectral function is related to the retarded propagator and the advanced one as . A tremendous simplification of thermal equilibrium is that the spectral and statistical functions are related by the fluctuation-dissipation relation, which is not assumed here .
In the following we will show how a power-law behavior (4) can be explained for sufficiently high momenta in resummed perturbation theory, if quantum corrections are neglected. Since the scaling solution is time and space translation invariant, the correlators in (6) can be Fourier transformed to and with four-momentum .
which can be directly verified by plugging in the above definitions. This equation is well-known in nonequilibrium physics and will be the starting point for our calculation. In the language of Boltzmann dynamics it states that ”gain terms” equal ”loss terms” for which stationarity is achieved . Thermal equilibrium trivially solves (8), which we do not consider in the following. Instead, we will look for possible non-thermal stationary solutions in perturbation theory relevant at sufficiently high momenta.
To exemplify the evaluation of the self-energies, we consider the one-loop diagram of Fig. 4 appearing in the two-particle irreducible (2PI) effective action scheme, where self-energies are expressed in terms of self-consistently dressed propagators while vertices remain undressed. Using and at one-loop order the gauge field contributions to the color-diagonal self-energies can be written as 
with the notation . The Lorentz and momentum structure of the two three-vertices appearing in the one-loop expression for SU(2) gauge theory is contained in
In accordance with (4), we are looking for scaling solutions which behave as
under rescaling with the real parameter . The spectral function is antisymmetric, , while and we assume symmetry in Lorentz indices. As in Ref. , scaling solutions can be efficiently identified by integrating (8) over external spatial momentum and suitable scaling transformations for coordinates. This is similar to what is done in the context of weak Kolmogorov wave turbulence for scalar models using a Boltzmann equation , and in this way the problem is reduced to simple algebraic conditions for exponents. Here, the analysis is complicated by the fact that for the perturbative result to order the two terms in (8) vanish separately on-shell, i.e. for . The non-vanishing contributions of the diagram in Fig. 4 to order can be conveniently analyzed by writing (7) as
This is used to replace one propagator line in the one-loop graph of Fig. 4. To order , the self-energies appearing in (12) can then be taken as (9). Following a well-trodden path [7, 13, 6, 5], we obtain by integrating over spatial momentum from the stationarity condition (8):
with . In (13) the first term in brackets, , vanishes for , and for in the on-shell limit for and such that only scattering contributes. With (14) one observes that these correspond to the well-known Kolmogorov exponents for energy and for particle cascades,
respectively. In (13) the second term in brackets, , does not vanish with (15). Doing the equivalent calculation starting from the classical-statistical gauge theory instead of the quantum theory leads to the same stationarity equation (13), however, without this term . Correspondingly, neglecting quantum corrections one observes rather good agreement of the perturbative solution with the value (5) indicated by the full numerical solution of the classical-statistical gauge theory. A similar analysis can be performed for the two-loop gluon diagrams.
For sufficiently high momenta, where characteristic occupation numbers are of order one, the quantum corrections appearing in (13) become relevant, which preempts a scaling solution (15) in the far UV. On the other hand, parametrically the perturbative estimate certainly breaks down for occupancies as large as , where order-one corrections to the self-energies occur at any order in a loop expansion. As a consequence, one expects a window of momenta for which the scaling solution (15) can describe quantitatively the dynamics. In view of this, it is remarkable that the simulation results shown in Fig. 2 exhibit for low momenta only small deviations from the perturbative behavior. Given that the numerical results are rather insensitive to going to lattices, simulations on much larger lattices might be required to settle the important question of whether a new nonperturbative infrared fixed point exists in nonequilibrium QCD, similar to what is observed in scalar theories .
Despite important differences the discussion of the thermalization process in heavy-ion collisions and cosmology after inflation shows remarkable similarities. In both cases nonequilibrium instabilities can lead to a fast period of exponential growth of occupation numbers. This is followed by a slow period, where the quantitative agreement concerns characteristic exponents or scaling functions. Our numerical as well as analytical results indicate that the nonabelian gauge theory can lead to perturbative Kolmogorov scaling exponents known from scalar models.
In QCD simulations for large volumes an ambitious analysis in terms of gauge invariant quantities such as stress-tensor correlation functions would be desirable for a quantitative understanding of the nonperturbative infrared dynamics. One can hope that this leads to an improved theoretical understanding of the plasma’s temporal evolution, which is crucial for the successful outcome of the experimental heavy-ion program.
This work is supported in part by the BMBF grant 06DA267, and by the DFG under contract SFB634. Some of the ideas were initiated during the program on ”Nonequilibrium Dynamics in Particle Physics and Cosmology” (2008) at the Kavli Institute for Theoretical Physics in Santa Barbara, supported by the NSF under grant PHY05-51164.
- S. Mrówczyński, Phys. Lett. B 214 (1988) 587; ibid. Phys. Lett. B 314 (1993) 118; ibid. Phys. Rev. C 49 (1994) 2191; ibid. Phys. Lett. B 393 (1997) 26; P. Arnold, J. Lenaghan and G. D. Moore, JHEP 08 (2003) 002; P. Romatschke and M. Strickland, Phys. Rev. D 68 (2003) 036004; S. Mrówczyński, A. Rebhan and M. Strickland, Phys. Rev. D 70 (2004) 025004; A. Rebhan, P. Romatschke and M. Strickland, Phys. Rev. Lett. 94 (2005) 102303; P. Arnold, G. D. Moore and L. G. Yaffe, Phys. Rev. D 72 (2005) 054003; A. Rebhan, P. Romatschke and M. Strickland, JHEP 09 (2005) 041; P. Romatschke and A. Rebhan, Phys. Rev. Lett. 97 (2006) 252301; B. Schenke and M. Strickland, Phys. Rev. D 74 (2006) 065004; B. Schenke, M. Strickland, C. Greiner and M. H. Thoma, Phys. Rev. D 73 (2006) 125004; D. Bödeker and K. Rummukainen, JHEP 07 (2007) 022; P. Arnold and G. D. Moore, Phys. Rev. D 76 (2007) 045009; A. Rebhan, M. Strickland and M. Attems, Phys. Rev. D 78 (2008) 045023.
- A. Dumitru and Y. Nara, Phys. Lett. B 621 (2005) 89; A. Dumitru, Y. Nara and M. Strickland, Phys. Rev. D 75 (2007) 025016; A. Dumitru, Y. Nara, B. Schenke and M. Strickland, Phys. Rev. C 78 (2008) 024909.
- P. Romatschke and R. Venugopalan, Phys. Rev. Lett. 96 (2006) 062302; ibid. Phys. Rev. D 74 (2006) 045011.
- J. Berges, S. Scheffler and D. Sexty, Phys. Rev. D 77 (2008) 034504.
- V. Zakharov, V. Lvov, and G. Falkovich, Kolmogorov Spectra of Turbulence, Wave Turbulence (Springer-Verlag, Berlin Heidelberg New York, 1992).
- For a review see R. Micha and I. I. Tkachev, Phys. Rev. D 70 (2004) 043538; ibid. Phys. Rev. Lett. 90 (2003) 121301.
- J. Berges, A. Rothkopf and J. Schmidt, Phys. Rev. Lett. 101 (2008) 041603.
- L. Kofman, A.D. Linde, A.A. Starobinsky, Phys. Rev. Lett. 73 (1994) 3195; ibid. Phys. Rev. D 56 (1997) 3258.
- J.H. Traschen, R.H. Brandenberger, Phys. Rev. D 42 (1990) 2491; Y. Shtanov, J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 51 (1995) 5438.
- S.Yu. Khlebnikov, I.I. Tkachev, Phys. Rev. Lett. 77 (1996) 219; ibid. Phys. Rev. Lett. 79 (1997) 1607; T. Prokopec, T.G. Roos, Phys. Rev. D 55 (1997) 3768.
- G. N. Felder, J. Garcia Bellido, P. B. Greene, L. Kofman, A. D. Linde and I. Tkachev, Phys. Rev. Lett. 87 (2001) 011601; G. N. Felder, L. Kofman and A. D. Linde, Phys. Rev. D 64 (2001) 123517; J. Garcia Bellido, M. Garcia Perez and A. Gonzalez Arroyo, Phys. Rev. D 67 (2003) 103501.
- J. Berges and J. Serreau, Phys. Rev. Lett. 91 (2003) 111601; A. Arrizabalaga, J. Smit and A. Tranberg, JHEP 0410 (2004) 017.
- J. Berges and G. Hoffmeister, Nucl. Phys. B 813 (2009) 383.
- P. Arnold and G. D. Moore, Phys. Rev. D 73 (2006) 025006; ibid. Phys. Rev. D 73 (2006) 025013.
- A. H. Mueller, A. I. Shoshi and S. M. H. Wong, Nucl. Phys. B 760 (2007) 145.
- For a review see J. Berges, AIP Conf. Proc. 739 (2005) 3; arXiv:hep-ph/0409233.
- A. Cucchieri and T. Mendes, Nucl. Phys. B 471 (1996) 263.
- M. Strickland, J. Phys. G 34 (2007) S429.
- J. Berges, Phys. Rev. D 70 (2004) 105010.