# Evolving Glasma and Kolmogorov Spectrum^{1}

^{1}

## Abstract

We present a pedagogical introduction to the theoretical framework of the Color Glass Condensate (CGC) and the McLerran-Venugopalan (MV) model. We discuss the application of the MV model to describe the early-time dynamics of the relativistic heavy-ion collision. Without longitudinal fluctuations the classical time evolution maintains boost invariance, while an instability develops once fluctuations that break boost invariance are included. We show that this “Glasma instability” enhances rapidity-dependent variations as long as self-interactions among unstable modes stay weak and the system resides in the linear regime. Eventually the amplitude of unstable modes becomes so large that the growth of instability gets saturated. In this non-linear regime the numerical simulations of the Glasma lead to turbulent energy flow from low-frequency modes to higher-frequency modes, which results in a characteristic power-law spectrum. The power found in numerical simulation of the expanding Glasma system turns out to be consistent with Kolmogorov’s scaling.

## 1 Introduction

Relativistic heavy-ion collision experiments have aimed to create a new state of matter out of color-deconfined particles, i.e. a quark-gluon plasma (QGP) in extreme environments in the laboratory. Presumably it was only up to after the Big Bang when the Early Universe was still hot enough to realize the QGP in nature. Experimental data on intrinsic properties of the QGP suggest that this new state of QCD matter found in the heavy-ion collision is not a weakly-coupled plasma but rather a strongly-coupled fluid. The hydrodynamic description of the time evolution has successfully reproduced the measured particle distributions, in particular, the azimuthal distribution of emitted particles in non-central collisions. The great success of the hydrodynamic model is a strong evidence for thermalization. After thermal equilibrium is achieved, the time evolution of the QGP is somehow under theoretical control, though some uncertainties remain in the determination of the equation of state and the implementation of dissipative effects.

In contrast to the hydrodynamic regime after thermalization, our understanding is still quite limited about the early-time dynamics toward thermalization. This kind of problem is in general one of the most difficult physics challenges. One would naturally anticipate that the system may form a turbulent fluid during transient stages right after the collision. In fact, turbulence is a common phenomenon in our daily life whenever the Reynolds number exceeds a certain threshold. Nevertheless, it poses a very difficult theory problem even today. It is widely known that a renowned physicist, Richard Feynman described turbulence as “the last great unsolved problem of classical physics” and it is indeed so also in the context of the QGP study.

A good news is that we already have a powerful theoretical tool to investigate the very early-time dynamics of the high-energy heavy-ion collision. One may well think at a first glance that any microscopic description of nucleus-nucleus collisions is too complicated to handle in terms of the QCD first principle. This is truly so unless the collisional energy is sufficiently high. Actually microscopic information of nucleus before collision is far from simple on its own. The high-energy limit, however, allows for drastic simplification that makes the calculation feasible. In the Regge limit, precisely speaking, the strong interactions exhibit totally different characteristics from low-energy hadron physics. In this particular case the c.m. energy scale is infinitely larger than other energy scales such as the transferred momentum squared . Then, although the strong coupling constant is small due to asymptotic freedom with large , a resummation is required for a series of non-small terms where Bjorken’s . In terms of the Feynman diagram this resummation represents quantum processes to emit softer gluons successively. Intuitively, as one goes to larger and thus goes to down from , one should consider more radiated gluons within a bin of to . Naturally the wave function of nucleus is -dependent and we should expect more and more gluons inside at smaller .

Eventually it should be the most suitable to treat gluons as coherent fields rather than particles once the gluon density is high enough. This is reminiscent of photons in the Weizsäcker-Williams approximation. In the first approximation in the high-energy limit of QCD, therefore, the classical fields are appropriate ingredients in theoretical computations, which is as a consequence of quantum radiations. Such a classical description of high-energy QCD is called the Color Glass Condensate (CGC) [1]. In this way, the initial condition for the relativistic heavy-ion collision should be formulated by means of the CGC theory. As long as the gluon distribution function stays large, the CGC picture holds, and it is finally superseded by a particle picture of plasma. Some unnamed physical state between the CGC and the QGP was given a name of Glasma as a mixture of “glass” and “plasma” [2]. The Glasma time evolution thus follows the classical equations of motion and hopefully it should be transformed smoothly into a hydrodynamic regime. In this sense the Glasma should play a central role to figure out what the initial condition for the hydrodynamic model should look like.

Because QCD or the Yang-Mills theory involves gluonic self-interactions, it is generally hard to find an analytical solution of the classical equations of motion except for some simple situations. Hence, one has to resort to numerical methods to go into quantitative estimates, and the numerical implementation has been established. In view of numerical outputs, however, there is no indication seen in the Glasma evolution toward thermalization. If all quantum fluctuations are completely frozen in particular, the classical dynamics respects boost invariance. The coordinate rapidity is simply shifted under a boost on the system, that means that the boost-invariant system is insensitive to . The pressure is then highly anisotropic depending on the beam-axis direction or on the transverse plane. This makes a sharp contrast to thermalized matter in which the pressure is isotropic.

Quantum (or structural – see Sec. 5) fluctuations including -dependence break boost invariance explicitly. Interestingly enough, it was discovered that a small modulation along the longitudinal direction grows up exponentially as a function of time and -dependent modes show instability behavior, which is sometimes referred to as the “Glasma instability” [3]. It is of paramount importance to understand the nature of the Glasma instability to fill in a missing link between the CGC initial state and the initial condition for hydrodynamics. There are several theoretical attempts to account for qualitative features of the Glasma instability [4, 5, 6, 7, 8, 9]. Here, instead of doing so, we will think about subsequent phenomena at later stages; turbulence may be formed by instability growth, and then it is sensible to anticipate the decay of turbulence and the associated scaling law in the energy spectrum [10].

This article contains lectures on the basic facts of the MV model for those who are not necessarily familiar with QCD at high energy. In Sec. 2 a stationary-point approximation on the functional integration is introduced, which leads to a classical treatment of the high-energy QCD problems. The classical equations of motion in the pure Yang-Mills theory are further discussed in Sec. 3. Then, some numerical results for the -independent case are presented next in Sec. 4, and those for the case with -dependent fluctuations in Sec. 5. Some evidence is presented for the realization of the power-law scaling in the energy spectrum in the non-linear regime where the instability stops. Section 6 is devoted to outlooks.

## 2 Scattering amplitude and the Eikonal approximation

We will see the essence of the Eikonal approximation and the scattering amplitudes of our interest in QCD physics. We will first consider the case of light projectile and dense target, and proceed to the case of dense projectile and dense target next.

### 2.1 Eikonal approximation

Before addressing QCD application, we shall first consider a scattering problem in non-relativistic Quantum Mechanics. If we want to solve a problem of potential scattering, we should treat the Schrödinger equation,

(1) |

with the following boundary condition,

(2) |

at large distance (). The term involving represents the scattered wave and one can obtain the scattering amplitude from . The lowest-order estimate immediately gives an expression in the Born approximation. In the high-energy limit, however, the scattering angle is small and the incident and the scattered waves interfere strongly there. In this situation the following Ansatz is more convenient,

(3) |

which is called the Eikonal approximation in analogy with the terminology in Optics. The wave length of incident wave, , is shorter than the potential range when is large enough, and the differential equation for is

(4) |

with . Then, in the high-energy limit, the second term is negligible as compared to the first term, which is easily integrated out to give,

(5) |

In the gauge theory is replaced by . Let us consider the scattering of the target particles moving at the speed of light in the positive- direction and the projectile particles in the negative- direction. Then, in general, in the Eikonal approximation, the scattering matrix is

(6) |

with the Wilson lines corresponding to the Eikonal phase (5),

(7) |

in the light-cone coordinates; and , as sketched in Fig. 1.

The momenta conjugate to and are the energy and the longitudinal momentum . In Eq. (6) the weights, and , represent the wave functions of the target and the projectile at the Bjorken variables and , respectively. Since the hard particles ( with being the total longitudinal momentum of the target) are included in the wave function, the functional integration in should contain the softer gauge fields with . In the following subsections let us discuss how to make an approximation on this scattering amplitude.

### 2.2 Light projectile and dense target

For the simplest example let us take the projectile as a color-dipole, i.e. the scattering amplitude is then,

(8) |

The product of ’s can be re-expressed on the exponential [11, 12] as

(9) |

where the last expressions is an approximation valid for sufficiently large . Then, with large , the functional integration in Eq. (8) can be estimated by means of the stationary approximation at the solution of

(10) |

The solution of the above classical equations of motion thus represents the contribution from soft gluons with . Then, finally, the dipole scattering amplitude is

(11) |

This expression is easily generalized to an arbitrary operator and

(12) |

Once the dependence in is known [13, 14], small- evolution is deduced for in general and in this way the BFKL equation up to the quadratic-order of , and besides, the BK [15] and JIMWLK equations including full-order of are derived.

### 2.3 Dense projectile and dense target

The discussions so far are quite generic but the above stationary-point approximation needs slight modifications when not only the target parton density but also the projectile is large as in the situation of the relativistic heavy-ion collision. Then, the scattering amplitude reads

(13) |

with the source action given approximately as

(14) |

This time, the stationary-point is shifted also by the effect of the presence of and it is determined by the classical equations of motion with two sources,

(15) |

Using the solution of the above classical equations of motion, one can obtain the general formula similar to previous Eq. (12) as

(16) |

In what follows we discuss how to solve this above Eq. (16) to investigate the early-time dynamics of the heavy-ion collision.

Then, there are two ingredients necessary to estimate physical observables using Eq. (16). One is the solution of Eq. (15), which is impossible to find in an analytical way unfortunately. The other is to figure out the wave functions and , which is again impossible to do so by solving QCD exactly. Given an initial condition at a certain , in principle, the evolution equation such as the JIMWLK equation leads to the wave function at any . The theoretical framework of such a description of scattering processes with non-linearity of abundant gluons is called the Color Glass Condensate (CGC). The CGC theory is not a phenomenological model but an extension of conventional perturbative QCD with resummation in a form of background fields. In any case, however, the initial condition for small- evolution is necessary for actual computations. This part needs an Ansatz as explained in the next section.

## 3 Equations of motion and the MV model

Here we will see that the analytical solution is written down for one-source problem. Although it is impossible to give an analytical formula for the solution of two-source problem, the initial condition on the light-cone can be specified. Also we will introduce a Gaussian approximation for the wave functions that defines the MV model.

### 3.1 One-source problem

Let us first consider how to solve the equations of
motion (10) for the light-dense scattering. It is actually
easy to find a special solution in the same way as the classical
electromagnetism. The important point is that the source
is independent of because of time dilation of
particles moving at the speed of light^{2}

(17) |

In later discussions it is more convenient to adopt the light-cone gauge, , which is achieved by a gauge rotation by that solves . That is,

(18) |

After the gauge transformation by , only the transverse components are non-vanishing,

(19) |

and . In the same way, for the projectile moving in the opposite direction to the target, the equations of motion have a solution as

(20) |

with

(21) |

in gauge. In this manner the one-source problem is readily solvable. It is, however, impossible to solve the two-source problem in Eq. (15) as simply as above.

### 3.2 Two-source problem

In the presence of two sources it is most suitable to make use of the Bjorken coordinates that represent an expanding system. The time and the longitudinal variables are replaced by the proper time and the space-time rapidity, respectively, as

(22) |

The temporal gauge in this coordinate, , has a close connection to the light-cone gauge discussed in the previous subsection because leads to a condition on and on . Therefore, the experience in the one-source problem turns out to be useful here.

The equations of motion can be expressed in the Bjorken coordinates as

(23) |

and the conjugate momenta are defined as

(24) |

These equations (23) and (24) determine the time evolution uniquely once the initial condition at an initial time is specified.

The solutions of the one-source problems are consistent with the gauge-fixing condition since we found for the target and the projectile both. Then solves Eq. (23) too if there is no interference from the projectile source. This means that should be a solution in the region outside of the light-cone, , and in (see Fig. 2 for illustration).

From this fact, it is naturally understood that the initial condition at or may be a superposition of these two solutions,

(25) |

Here, though this is a very simple Ansatz by a superposition, the consequence is quite non-trivial. The field strength associated with these gauge fields is

(26) |

and thus, the longitudinal component of the chromo-magnetic fields appears from the non-Abelian interactions. By solving the equations of motion, one can also find similar expressions for the chromo-electric fields [16]

(27) |

These field strengths stand for characteristic properties of the initial condition of the relativistic heavy-ion collisions in the CGC or the so-called Glasma picture, an intuitive illustration for which is displayed in Fig. 3. It is important to point out that the initial conditions at are independent of , namely, boost invariant. Because there is no explicit in the equations of motion, boost invariance is kept during the time evolution.

### 3.3 McLerran-Venugopalan model

One of the simplest and reasonable Ansätze for the wave function is the Gaussian approximation, that is given by

(28) |

which defines the McLerran-Venugopalan (MV) model. Here represents either or . In this model setup characterizes the typical energy scale. Indeed, once the parton saturation manifests itself, any details in the structure are lost and only the transverse parton density should be a relevant scale. In principle in the MV model is to be interpreted as in the parton saturation. With the Gaussian wave function (28), the expectation value is obtained by a decomposition into the two-point function that is read as

(29) |

In other words, the Gaussian approximation (28) assumes no correlation at all between spatially distinct sites. Evaluating the Gaussian average with various functionals of is an interesting mathematical excercise [17].

Especially it is feasible to evaluate the initial energy density,

(30) |

at using the initial conditions (25) and (27) and the Gaussian weight (28). The results are a bit disappointing because it involves both the UV and the IR divergences, which can be regularized by (momentum cutoff) and (gluon mass). Then, the initial energy density is found to be

(31) |

after some calculations [4, 18, 19]. In the numerical simulation with discretized grid in a finite-volume box [20, 21] there are natural UV and IR cutoffs incorporated from the beginning. The IR cutoff is, however, not included as a gluon mass but originates from the system size .

## 4 Numerical method and the boost-invariant results

The model parameters should be fixed first. In the case of the gold-gold collision at , the empirical choice is and , which correspond to , , and . Now that the model definition and the model parameters are given, we can calculate and numerically. Then, in the high-energy limit, it is a common assumption that the nucleus source is infinitesimally thin, i.e.

(32) |

Then, the Wilson lines are replaced, respectively, by

(33) |

with the solution of the Poisson equation,

(34) |

One should be very careful about this replacement because this is not an approximation on Eqs. (18) and (21). Even though the longitudinal extent in the color source is infinitesimally thin, it should not be legitimate to drop the path ordering [22]. Therefore, we have to think that the numerical MV model is something distinct from the original MV model in the continuum variables.

In the MV model the distribution of the color source is random and there is no correlation between different sites. Figure 4 illustrates an example of the initial as a solution of the Poisson equation. Because the operator involves spatial average, we see that the spatial distribution of is rather smooth even though the source has random fluctuations. This smoothness is not physical, however, and the gauge fields are furiously fluctuating as shown in the right panel of Fig. 4. It is worth noting that the color-flux tube picture as sketched in Fig. 3 is not the case in the MV model and the JIMWLK evolution is indispensable to take account of the flux tube structure.

The physical observables are measured by taking an ensemble average of results with different initial ’s. It is useful to compute not only the energy density (31) but also other combinations of the energy-momentum tensor. In particular the following pressures are important in order to judge how anisotropic the system is;

(35) | |||

(36) |

where the longitudinal and transverse chromo-electric and chromo-magnetic fields are defined as

(37) | |||

(38) |

The numerical results from the numerical Glasma simulation are presented in Fig. 5. From this figure it is clear that there are only longitudinal fields and right at the collision () as explained with Fig. 3. The transverse fields are developing as increases, and eventually the longitudinal and the transverse fields approach each other at . This does not mean isotropization, however, because there are two ( and ) components in the transverse direction. One can understand what is happening by taking a careful look at the pressures. The longitudinal pressure goes to zero for and this means that particles move on the free stream along the longitudinal direction. Thus, the system remains far from thermal equilibrium.

Here we encounter a perplexing situation. We know the the CGC should be a correct description of the initial dynamics of the heavy-ion collision. Even if the CGC cannot reach isotropization and thermal equilibrium, it should be a natural anticipation that the CGC can at least capture a correct tendency toward thermalization. This anticipation is not the case at all, however, as seen in Fig. 5. Then, is there anything missing in our discussions so far?

## 5 Glasma instability and the scaling spectrum

Yes, there is. We have neglected fluctuations on top of boost-invariant CGC-background field and thus dependence at all. Such a treatment is not always justified. As a matter of fact, longitudinal structures in a finite extent by not using but treating correctly in the path ordering would give rise to -dependent fluctuations. Also, quantum fluctuations have dependence as well [8, 23].

One might have thought that small disturbances in the longitudinal direction could make only a slight difference. But, the fact is that there is a tremendous difference between results with and without -dependent fluctuations. With random fluctuations in space, the Glasma simulation would lead to significant decay from the CGC background fields with -independent zero-mode into -fluctuating non-zero modes [3]. An example is shown in Fig. 6.

To obtain the results as presented in Fig. 6, only modes are disturbed,

(39) |

where is the longitudinal extent that we took as in the simulation. Once is given, the transverse fluctuations, , are chosen in such a way to satisfy the Gauss law. Then, physical observables of our interest are Fourier transformed from space to space. It is now obvious that the CGC background fields at are relatively larger than fluctuations by small and the second dominant mode should sit at . In this article we limit ourselves to the simplest choice of .

Usually unstable modes grow up exponentially, but in the expanding geometry, the instability implies a slightly weaker growth according to . The horizontal axis in the left panel of Fig. 6 is, thus, not the dimensionless time itself, but . We can surely confirm that the longitudinal pressure component at increases linearly in Fig. 6 when plotted with the logarithm of the pressure as a function of .

The detailed structure of the instability in the spectral pattern is interesting to see. The right panel of Fig. 6 is the spectrum corresponding to the simulation results shown in the left.

From now on, let us consider the fate of the instability; it is simply impossible for the instability to last for ever. The saturation of instability growth can be observed in two ways. The first case is that the initial magnitude (appearing in Eq. (39)) is taken to be large enough to accommodate non-linear effects. The second is the large-time behavior – simply we wait until the unstable modes grow up considerably.

As shown in the left panel in Fig. 7 the instability for stops and the spectrum is flattened after all. In the saturated regime at we see that the pressure at slightly decreases rather than increasing. This behavior can be interpreted as follows; As long as the non-zero modes are small (i.e. in the linear regime), the energy decay from the dominant CGC background at makes non-zero modes enhanced exponentially. The injected energy is much bigger than the escaped energy toward higher- modes in the linear regime. This balance changes gradually with increasing amplitude of unstable modes, and eventually a steady energy flow is expected when non-Abelian nature of non-zero modes becomes substantially large (i.e. in the non-linear regime). At even larger time, as hinted from Fig. 7, the spectrum is flattened and the UV-cutoff effects at large should be appreciable. Then, an intriguing question is; what is going on in the non-linear regime before the UV-cutoff effects contaminate the simulation?

To address this question, we shall take the latter case, namely the long-time run of the simulation with small for better numerical stability. The qualitative features in the results for shown in Fig. 8 are just the same as the top curve in Fig. 7 that represents the results for . We are interested in the energy spectrum in the non-linear regime that can be immediately identified on Fig. 8.

Figure 9 is the energy spectrum (not the Fourier-transformed longitudinal pressure but the energy contribution from respective modes). That is, what is shown in Fig. 9 is given as

(40) |

We see that there is a clear tendency to approach a scaling form in the non-linear regime (as in the right panel of Fig. 9). With some rescaling we find that this power is exactly consistent with the Kolmogorov value, . Generally speaking, in non-expanding systems, it is not a surprise that the Kolmogorov power emerges because this value can be guessed by dimensional analysis. Such a dimensional argument may work even for the expanding system. In the Bjorken coordinates is dimensionless, but the physical scale in the longitudinal direction is to be interpreted as , and thus the corresponding momentum should be . Then, this quantity gives a dimension of length (with multiplied appropriately). Because of an expected energy flow in space, its rate is also a relevant quantity. Then, the characteristic length and time scales of the system are expressed by two quantities with the following dimensions;

(41) |

The energy spectrum, on the other hand, has the dimension,

(42) |

that is reproduced exactly by a unique combination of . Therefore, it is concluded that should exhibit the power-law scaling in terms of whose power is identical to the Kolmogorov’s value, as long as the system stays in the non-linear regime. From the left panel of Fig. 9, we can see only the scaling region or the so-called inertial region realized. The dissipative region at high is not found, probably because of the UV-cutoff effects.

This nice scaling is lost at further later time. We can understand from the right panel of Fig. 9 that the energy flow is stuck at the UV edge and the energy spectrum is artificially pushed up by the UV-cutoff contamination.

## 6 Outlooks

It was a surprise that the Kolmogorov’s scaling law could be realized in the Glasma simulation in the expanding geometry. The dimensional argument is not so strong to constrain the shape of the energy spectrum uniquely. It should be an important test whether the Kolmogorov’s behavior can be confirmed or not in other simulations of the pure Yang-Mills theory. It has been established that the non-Abelian plasma generally has an instability associated with anisotropy that grows up until the non-Abelianization occurs [24]. Similar phenomena of instability tamed by non-Abelian interactions are found in many other examples. Then, presumably, there must appear the power-law scaling in the region around the non-Abelianization.

The turbulent decay and the associated Kolmogorov’s power-law are interesting discoveries from the long-run simulation of the Glasma. But, it cannot answer anything about the realistic thermalization mechanism from the Glasma. The turbulence is certainly a tendency into thermalized matter, but the energy flow is a slow process and the relevant time scale, , is far outside of the validity region of the CGC description. There must be still something missing that can accelerate the thermalization speed.

If this missing piece were finally set in, the transient Glasma would provide us with the initial input for the hydrodynamic equations. Even in this case, the analysis we have seen should be useful. The energy spectrum in space should carry important information. Then, the power could agree with or deviate from as suggested in strong-coupling studies [25, 26].

## Acknowledgments

The author thanks the organizers of the 51st Cracow School of Theoretical Physics for a kind invitation. This article summarizes lectures given there. He also thanks Hiro Fujii, Francois Gelis, and Yoshimasa Hidaka for collaborations. The central parts in these lectures are based on works done in collaborations with them.

### Footnotes

- thanks: Presented at the 51st Cracow School of Theoretical Physics
- Even though the time dependence is frozen at the speed of light, non-commutativity of color charges needs distinction by the order of interactions along the -axis, which introduces a label that plays the role of time. Such “-dependence” is dropped off in the large limit only.

### References

- For reviews, see: E. Iancu, A. Leonidov and L. McLerran, arXiv:hep-ph/0202270; E. Iancu and R. Venugopalan, arXiv:hep-ph/0303204; L. McLerran, arXiv:hep-ph/0311028; F. Gelis, T. Lappi and R. Venugopalan, Int. J. Mod. Phys. E 16, 2595 (2007).
- T. Lappi and L. McLerran, Nucl. Phys. A 772, 200 (2006).
- P. Romatschke and R. Venugopalan, Phys. Rev. Lett. 96, 062302 (2006); Phys. Rev. D 74, 045011 (2006).
- K. Fukushima, Phys. Rev. C76, 021902 (2007).
- H. Fujii, K. Itakura, Nucl. Phys. A809, 88-109 (2008). [arXiv:0803.0410 [hep-ph]].
- A. Iwazaki, Prog. Theor. Phys. 121, 809-822 (2009).
- K. Dusling, T. Epelbaum, F. Gelis, R. Venugopalan, Nucl. Phys. A850, 69-109 (2011).
- K. Dusling, F. Gelis, R. Venugopalan, [arXiv:1106.3927 [nucl-th]].
- T. Epelbaum, F. Gelis, [arXiv:1107.0668 [hep-ph]].
- K. Fukushima, F. Gelis, [arXiv:1106.1396 [hep-ph]].
- J. Jalilian-Marian, S. Jeon and R. Venugopalan, Phys. Rev. D 63, 036004 (2001).
- K. Fukushima, Nucl. Phys. A770, 71-83 (2006).
- E. Iancu, A. Leonidov and L. D. McLerran, Nucl. Phys. A 692, 583 (2001).
- E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran, Nucl. Phys. A 703, 489 (2002).
- Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999).
- A. Kovner, L. D. McLerran and H. Weigert, Phys. Rev. D 52, 6231 (1995); ibid. D 52, 3809 (1995).
- K. Fukushima, Y. Hidaka, JHEP 0706, 040 (2007).
- T. Lappi, Phys. Lett. B643, 11-16 (2006).
- H. Fujii, K. Fukushima, Y. Hidaka, Phys. Rev. C79, 024909 (2009).
- A. Krasnitz and R. Venugopalan, Nucl. Phys. B 557, 237 (1999); Phys. Rev. Lett. 84, 4309 (2000); ibid. 86, 1717 (2001).
- A. Krasnitz, Y. Nara and R. Venugopalan, Phys. Rev. Lett. 87, 192302 (2001); Nucl. Phys. A 717, 268 (2003); ibid. 727, 427 (2003).
- K. Fukushima, Phys. Rev. D77, 074005 (2008).
- K. Fukushima, F. Gelis, L. McLerran, Nucl. Phys. A786, 107-130 (2007).
- P. B. Arnold, G. D. Moore, and L. G. Yaffe, Phys. Rev. D72, 054003 (2005).
- J. Berges, D. Gelfand, S. Scheffler, and D. Sexty, Phys. Lett. B677, 210–213 (2009).
- M. Carrington and A. Rebhan, [arXiv:1011.0393 [hep-ph]].