Heavy Quarks in Turbulent QCD Plasmas

The quark-gluon plasma, which is produced at an early stage of ultrarelativistic heavy-ion collisions, is expected to be initially strongly populated with chromodynamic fields. We address the question how heavy quarks interact with such a turbulent plasma in comparison with an equilibrated one of the same energy density. For this purpose we derive a Fokker-Planck transport equation of heavy quarks embedded in a plasma of light quarks and gluons. We first discuss the equilibrium plasma and then the turbulent one applying the same approach, where the heavy quarks interact not with the plasma constituents but rather with the long wavelength classical fields. We first consider the three schematic models of isotropic trubulent plasma and then the simplified model of glasma with the chromodynamic fields only along the beam direction. The momentum broadening and collisional energy loss of a test heavy quark are computed and compared to those of equilibrium plasma of the same energy density.

The quark-gluon plasma, which is produced at an early stage of ultrarelativistic heavy-ion collisions, is expected to be initially strongly populated with chromodynamic fields.We address the question how heavy quarks interact with such a turbulent plasma in comparison with an equilibrated one of the same energy density.For this purpose we derive a Fokker-Planck transport equation of heavy quarks embedded in a plasma of light quarks and gluons.We first discuss the equilibrium plasma and then the turbulent one applying the same approach, where the heavy quarks interact not with the plasma constituents but rather with the long wavelength classical fields.We first consider the three schematic models of isotropic trubulent plasma and then the simplified model of glasma with the chromodynamic fields only along the beam direction.The momentum broadening and collisional energy loss of a test heavy quark are computed and compared to those of equilibrium plasma of the same energy density.

I. INTRODUCTION
The early stage of relativistic heavy-ion collisions is least known because there are hardly any experimentally accessible signals of the phase.Nevertheless one expects that the quark-gluon plasma, which is produced in the collisions, is initially strongly populated with chromodynamic fields.Within the framework of the Color Glass Condensate (CGC) approach, see e.g. the review [1], color charges of partons confined in the colliding nuclei act as sources of long wavelength chromodynamic fields which can be treated classically because of large occupation numbers of the soft modes.Since the density of color charges per transverse area of heavy nuclei is large, the corresponding momentum scale Q s is expected to be significantly bigger than the QCD scale parameter Λ QCD .Consequently, the coupling constant α s is presumably sufficiently small and perturbative methods are applicable.The system, however, is rather strongly interacting because of the high-amplitude fields present in the system.
A momentum anisotropy of the early stage quark-gluon plasma makes it unstable with respect to chromomagnetic modes which in turn cause a spontaneous generation of the fields, as explained at length in the review article [2].Therefore, the effect of strong fields is further enhanced.Following the terminology of electromagnetic plasma, we call such a nonequilibrium system of fields as the turbulent plasma meaning that numerous modes are excited in the system.In the CGC approach the non-equilibrium system of fields from the early stage of relativistic heavy-ion collisions is called glasma [1] and it can be treated as a specific realization of the turbulent QCD plasma.Leaving aside the mechanism of field generation and its detailed structure, one asks what are the transport properties of turbulent plasmas.We are specifically interested how heavy quarks -charm or beauty -behave in such a system when compared to the equilibrium plasma of the same energy density.
Heavy quarks are often treated as a probe of strongly interacting matter created in relativistic heavy-ion collisions, see e.g. the review [3].Thanks to their large masses the quarks are produced only at the earliest stage of the collision due to hard interactions of partons from incoming nuclei.Later on they propagate through a surrounding medium testing the entire history of the system.It has been long believed that the interaction of heavy quarks is significantly weaker than that of light quarks or gluons but experimental data clearly contradict the expectation.As discussed in the review [3], the behavior of mesons containing a heavy quark is rather similar to that of light mesons at both small and large transverse momenta.The problem is not fully resolved.
The medium created in relativistic heavy-ion collisions evolves fast towards the locally equilibrated quark-gluon plasma which expands hydrodynamically and ultimately is converted into a hadron gas.Final momentum spectra of heavy quarks are mostly shaped in the long-lasting equilibrium phase which is relatively well understood.An effect of a pre-equilibrium phase is often entirely ignored but this transient phase can significantly influence heavy-quark spectra because of its high density.Non-equilibrium calculations recently performed in a framework of kinetic theory [4] confirm the suggestion.However, we are interested even in the earlier phase when the medium is not described in terms of quasi-particles, as in a kinetic theory, but rather as a system dominated by classical fields which is the turbulent plasma.
A simple parametric estimate suggests that the interaction of heavy quarks in a turbulent plasma is much stronger than in the equilibrium one, if the coupling constant g is small.The momentum broadening parameter q, for example, is of order g 4 in equilibrium plasmas.Since the quark of interest actually interacts with soft gluons emitted by plasma constituents, one can think that the factor g 4 is composed of two pieces of g 2 .The first one is related to the gluon emission and the second one to the gluon absorption.If soft chromodynamic fields are present in the plasma, the interaction of the quark should be rather of order g 2 than g 4 .In Sec.VI we argue that q is indeed not of the order g 4 , not even g 2 but presumably of the order g in a turbulent plasma.
Because of their big masses, relaxation times of heavy quarks, which are produced in relativistic heavy-ion collisions, are expected to be significantly longer than that of light quarks and gluons.When an equilibrium or, more generally, a stationary state is reached by light quarks and gluons, heavy quarks need some extra time to adjust to the state of the plasma.Such a situation is naturally described in terms of the Fokker-Planck transport equation which was indeed repeatedly applied to heavy quarks in [5][6][7][8].The equation is usually derived from the Boltzmann equation by applying the so-called diffusion approximation to the collision term [9].The approximation assumes that the momentum transfer to the heavy quark in every collision is much smaller than the quark momentum.
The aim of this paper is threefold.In the first part we rederive the Fokker-Planck equation of heavy quarks which do not interact with plasma constituents but rather with soft classical fields present in the plasma.Specifically, we apply the so-called quasi-linear theory known from the electromagnetic plasma [9,10].The theory assumes that the distribution function can be decomposed into a large but slowly varying regular part and a small fluctuating or turbulent one which oscillates fast.The average over a statistical ensemble of the turbulent part is assumed to vanish and thus the average of the distribution function equals its regular part.The turbulent contribution to the distribution function obeys the collisionless transport equation while the transport equation of the regular part is determined by the fluctuation spectra which provide the collision term.The derivation presented here closely follows the procedure which was developed for QCD in [11], where, however, only the longitudinal chromoelectric field was taken into account and here the complete chromodynamic field is considered.The equilibrium correlation functions of chromodynamic fields, which are needed to obtain the quasi-linear transport equations, were derived in [12].
Our second aim is to confront the equilibrium plasma with the turbulent one.Therefore, we consider three models of isotropic turbulent plasma in the second part of the paper.Postulating a form of the correlation functions of chromodynamic fields, we derive the coefficients of the Fokker-Planck equation which can be related to the energy loss, momentum broadening and diffusion coefficient of heavy quarks in the plasma.The transport coefficients of turbulent plasma are compared to those of the equilibrium one at the same energy density.
The third aim is to study an evolution of heavy quarks at the earliest stage of relativistic have-ion collisions.Since the chromoelectric and chromomagnetic fields spanned between the receding nuclei are initially mostly parallel to the beam direction, we model the glasma with the boost invariant correlation functions of longitudinal fields.The energy loss and momentum broadening of heavy quarks are computed, assuming that all energy of the glasma is accumulated in the longitudinal fields.
At the end of the introductory remarks we note that an approach similar to ours, which was also inspired by the electromagnetic plasma studies [13], was formulated in [14], see [15,16] as well.We also mention an attempt [17] to study transport of heavy quarks in a plasma populated by strong chromodynamic fields.Unfortunately, the paper is flawed as the framework of an isotropic Langevin approach is applied to anisotropic plasmas.
Throughout the paper we use the natural system of units with c = = k B = 1; our choice of the signature of the metric tensor is (+ − −−).Lorentz indices are denoted with µ, ν = 0, 1, 2, 3 and i, j = 1, 2, 3 label the Cartesian coordinates x, y, z.The color indices of the adjoint representation of SU(N c ) gauge group are a, b = 1, 2, . . .N 2 c − 1.

II. DERIVATION OF FOKKER-PLANCK EQUATION
Our derivation of the Fokker-Planck equation of heavy quarks embedded in quark-gluon plasma starts with the transport equation of the Vlasov form where the distribution function Q(t, r, p) of heavy quarks is the N c × N c hermitian matrix which belongs to the fundamental representation of the SU(N c ) group.The distribution function depends on the time (t), position (r) and momentum (p) variables.There is no explicit dependence on the time-like component of the four-momentum p µ = (p 0 , p) as the distribution function is assumed to be non-zero only for momenta obeying the mass-shell constraint that is p 0 = E p = p 2 + m 2 .The quark velocity equals v = p/E p and being the chromodynamic potential in the fundamental representation.The mean-field term of the transport equation ( 1) is expressed through the color Lorentz force F(t, r) ≡ g E(t, r) + v × B(t, r) with the chromoelectric E(t, r) and chromomagnetic B(t, r) fields also belonging to the fundamental representation.The symbol {. . ., . . .} denotes the anticommutator.The derivation of the transport equation ( 1) is discussed in detail in the review [2].Further on we assume that the chromodynamic fields and the distribution function which enter the transport equation ( 1) can be decomposed into a regular and fluctuating or turbulent component.The distribution function is thus written down as where • • • denotes ensemble average; Q(t, r, p) is called the regular part while δQ(t, r, p) is called the fluctuating or turbulent one.It directly follows from Eq. ( 2) that δQ = 0.The regular contribution is assumed to be color neutral or white, and it is expressed as where 1 is the unit matrix in color space.Since the distribution function transforms under gauge transformations as Q → U QU † , where U is the transformation matrix, the regular contribution of the form (3) is gauge independent.We also assume that but at the same time What concerns the chromodynamic fields, we assume in accordance with Eq. ( 3) that their regular parts vanish and thus We substitute the distribution function (2) into the transport equation ( 1) and linearize the equations in the fluctuating contributions.Thus we get the equation where D ≡ ∂ ∂t + v • ∇ is the substantial or material derivative.Now we substitute the distribution functions (2) into the transport equations (1) but instead of linearizing the equation in the fluctuating contributions, we take the ensemble average of the resulting equation and trace over the color indices.Thus we get Since the regular part of distribution function is assumed to be color neutral, see Eq. ( 3), the term Tr[ F • ∇ p n ] vanishes because the fields E, B are traceless.The trace over color indices also cancels the terms originating from covariant derivatives like Tr [A µ , δQ] .We finally note that the trace Tr[ F • ∇ p δQ ] is gauge independent as the regular distribution function n(t, r, p) is.Now, we are going to write down the transport equation (8) in the Fokker-Planck form.For this purpose we observe that due to the condition (5), the space-time dependence of the regular distribution function can be neglected in the linearized transport equation (7) and then, the equation becomes easily solvable.We solve it with the initial condition using the one-sided Fourier transformation defined as The inverse transformation is where the real parameter σ > 0 is chosen in such a way that the integral over ω is taken along a straight line in the complex ω−plane, parallel to the real axis, above all singularities of f (ω, k).
The linearized transport equation (7), which is converted into the algebraic equation by means of the one-sided Fourier transformation, is solved as We stress that although we have ignored the (weak) frequency and wave number dependence of the regular distribution n, the fields E(ω, k), B(ω, k) retain their full frequency and wave number dependence in the expression (12).Inverting the one-sided Fourier transformation, one finds the solution of the linearized transport equation as where we assumed that E(ω, k) and B(ω, k) are analytic functions of ω.
With the help of the solution ( 13), the force term in the transport equation (8) becomes The second term in the r.h.s. of Eq. ( 14) can be manipulated to the form which is effectively the definition of the vector Y i (v).We also introduce the tensor and we note that, as explained in the subsequent sections, X and Y become time independent for a sufficiently long t.Then, the transport equation ( 8) can be written as the Fokker-Planck equation Since the distribution function n(t, r, p) carries no information about color degrees of freedom, the function is gauge invariant, and consequently X ij (v) and Y i (v) should be gauge invariant as well.However, one observes that X ij (v) and Y i (v) as defined by Eqs. ( 15) and ( 16) are gauge dependent because the traces are of nonlocal quantities in the definitions (15) and (16).The starting transport equation ( 1) is gauge covariant but the linearization procedure breaks the covariance because the covariant derivative is replaced by the normal one.Consequently, the solution ( 13) is not gauge covariant -the right-hand side of Eq. ( 13) transforms differently under local gauge transformations than the left-hand side.To cure the problem, one modifies the solution ( 13) by means of the link operator which is also called the gauge parallel transporter, see e.g.Sec.IIIE of the review article [2].Then, the modified solution obeys Eq. ( 8) with the covariant derivative instead of the normal one.Let us briefly discuss the procedure in a context of the Fokker-Planck equation (17).
According to Eq. ( 16), the tensor X ij (v) is determined by the traces of the field correlation functions like where chromodynamic fields are written in the adjoint representation of the SU(N c ) group which is used further on.The trace becomes gauge invariant under the replacement where Ω ab (t 1 , r 1 |t 2 , r 2 ) is the link operator defined as Here T c is the adjoint representation generator of the SU(N c ) group and P denotes the ordering along the path connecting the points (t 2 , r 2 ) and (t 1 , r 1 ).Since the fields transform as vectors under the local gauge transformation U(t, r) and the link transforms as one checks that the trace of the correlation function which includes the link is indeed gauge invariant.Consequently, the tensor X ij (v) is gauge invariant.Analogously one achieves the gauge invariance of the vector Y i (v).Further on, whenever any nonlocal correlation function shows up a presence of the link operator is implicitly assumed even so it is not explicitly written.

III. FOKKER-PLANCK EQUATION
Although this is a textbook material we briefly discuss here the Fokker-Planck equation (17).We first note that in the isotropic plasma the tensor X ij (v) and vector Y i (v) both depend on a single vector that is the heavy-quark velocity v. Therefore, they can be written as where v ≡ |v| and the coefficients X L (v), X T (v) and Y (v) are equal to The equilibrium distribution function of the form with T being the temperature of the plasma of light quarks and gluons, where heavy quarks are embedded, is expected to solve the transport equation ( 17).This is indeed the case if the coefficients X ij (v) and Y i (v) obey the condition which in the isotropic plasma reads When the plasma is isotropic and the coefficients X L (v) and X T (v) are equal to each other and independent of v, the Fokker-Planck equation reads where The quantities X ij (v) and Y i (v) have a clear physical meaning.As discussed in e.g. the classical monograph [20], the average momentum change per unit time and the correlation of momentum changes per unit time are given as Using the formulas (28) and (29), X ij (v) and Y i (v) can be related to the collisional energy loss dE dx and transverse momentum broadening q of a heavy-quark in the quark-gluon plasma, which play an important role in a theoretical description of the jet quenching phenomenon.The parameter q controls the radiative energy loss in a plasma medium [18].One easily finds that Since ∆x = v∆t, the energy loss per unit path equals which in isotropic plasmas reads The coefficient q, which is the broadening per unit path of the distribution of the test parton's momentum transverse to the initial parton's momentum, is immediately found as and in an isotropic plasma it equals When we deal with an equilibrium plasma and the coefficients X L (v), X T (v) are equal to each other and independent of v, the Fokker-Planck equation can be related to the nonrelativistic Langevin equation [20].Then, the diffusion constant D can be expressed as We are mostly interested in turbulent QCD plasmas populated with strong chromodynamic fields but we start with the equilibrium system where the fields are only at a level of thermal noise.We rederive the known Fokker-Planck equation [6] in a different way to demonstrate the reliability of our approach.

IV. EQUILIBRIUM PLASMA
We assume that the quark-gluon plasma, in which heavy quarks are embedded, is in thermodynamical equilibrium and we first derive in this section explicit expressions for the coefficients X L , X T and Y which enter the Fokker-Planck equation.

A. Computation of X and Y
As the formula (16) shows, the quantity X is given by the correlations functions E i (t, r)E j (t , r ) , B i (t, r)B j (t , r ) , E i (t, r)B j (t , r ) , and B i (t, r)E j (t , r ) which were studied in detail in [12].The explicit expressions are collected in Appendix A. Since the correlation functions are of the structure (A1), the tensor X ij is written as where the chromodynamic fields are expressed in the adjoint representation of the SU(N c ) group and F i a F j a ω, k is the fluctuation spectrum.For a translationally invariant system it is defined as A more general definition is discussed in Appendix B. Combining the equilibrium fluctuation spectra (A2), (A3) and (A4), one finds where β ≡ T −1 and ε L,T (ω, k) are chromodielectric functions which for the equilibrium plasma of massless particles are also given in Appendix A.
After performing the elementary time integration in Eq. ( 36), one is left with the integral over the four-vector k µ = (ω, k).Taking into account only the terms of the integrand which are even as a function of k µ and give nonzero contributions, one obtains where ω ≡ k • v.In the limit t → ∞, we have lim and thus we get One can show that the expression (41) properly approximates the formula (39) if the spectrum F i a F j a ω, k weakly changes as a function of ω in the interval [ω − π/t, ω + π/t].When the time grows the condition is easier and easier to fulfill.
Since the plasma under consideration is isotropic, the tensor X ij (v) is fully determined by the two functions X L (v) and X T (v) introduced in Eq. (21).Substituting the explicit form of the fluctuation spectrum (38) into Eq.( 41), one finds the coefficients X L (v) and X T (v) as , which is determined by the correlations functions E i (t, r)δQ 0 (r , p ) , and B i (t, r)δQ 0 (r , p ) can be derived directly from the formula (15).Such a derivation for a simplified case of only longitudinal electric field present in the plasma can be found in [11] where it is shown that the condition (25) or ( 26) is indeed satisfied.Here instead we refer to the condition (25) to obtain Y i (v).
Once the coefficients X ij (v) is given by Eq. ( 21) together with Eqs.(42), (43) and Y i (v) by Eqs. ( 22) and ( 26), the Fokker-Planck equation ( 17) is fully specified.We note that X ij (v) and Y i (v) depend on heavy quark momentum p and its mass m only through the velocity v = p/ m 2 + p 2 .Therefore, X L (v), X T (v) and Y (v) become independent of p when the heavy quarks of interest are truly relativistic and p 2 m 2 .Although, the coefficients X L (v), X T (v) and Y (v) are independent of the quark mass, the Fokker-Planck equation (17) does depend on m which is evident when the momentum derivatives are replaced by the velocity derivatives.
Since the coefficients X L (v), X T (v) and Y (v) depend on the quark mass only through the velocity, one might think that the corresponding Fokker-Planck equation is applicable to quarks of any mass.However, it is not true.As mentioned in the introduction, the typical momentum transfer in a single collision, which is of order gT , must be much smaller than the quark momentum that is gT mv/ √ 1 − v 2 .And here the quark mass matters.

B. Limit of small v
In the limit of small velocities of heavy quarks, the equilibrium Fokker-Planck equation gets a simpler form and the coefficients X L (v), X T (v) can be estimated analytically.Indeed, when heavy quarks are nonrelativistic (v 2  1), we have ω2 k 2 and one can use the approximate formulas (A8) and (A9).Then, Eqs. ( 42) and (43) read where the contributions originating from ε T (ω, k) appear to vanish.Introducing spherical coordinates with the axis z along the vector v, the integrals (44) and (45) are rewritten as where the trivial azimuthal integrals are performed and x ≡ cos θ with θ being the angle between v and k.
When T m D , the integrals ( 46) and (47) can be estimated as follows.One first observes that the dominant contribution comes from k ∈ [m D , T ].Assuming that βkxv 1, the integrals over x are easily computed and one obtains Approximating the integrand as k −1 , we finally get As one sees in Eq. ( 48) or (49), X L (v) and X T (v) are independent of v and equal to each other.The formula (35) is therefore applicable and the inverse diffusion constant equals which agrees with Eq. ( 12) of the study [5] when only the logarithmic term is taken into account.

C. Numerical results
Here we show some numerical results for the equilibrium QGP of N c = 3 and N f = 2.The Debye mass is computed according to the formula (A7).Fig. 1 presents the coefficients X L (v) and X T (v), which are obtained directly from Eqs. ( 42) and (43), as functions of the velocity v.The coupling constant and the temperature are α s ≡ g 2 /4π = 0.1 and T = 0.5 GeV.In Fig. 2 we show how X L (v) and X T (v) depend on the temperature T .The coupling constant is again α s = 0.1 and the velocity equals v = 0.8.
We have checked that the values of X L (v) and X T (v), which we obtained numerically, agree rather well with those computed by Svetitsky [6] except in the domain of small velocities v ≤ 0.1.The agreement is not trivial because the coefficients of the Fokker-Planck equation were derived in [6] from the matrix elements of heavy quark binary interactions with plasma constituents.To remove infrared divergences of the matrix elements, a mass parameter, corresponding to the Debye mass, was included in the gluon propagator.Since the procedure is not very accurate, it presumably explains the difference with our X L (v) and X T (v) in the domain of small velocities.On the other hand our approach does not treat properly the interactions with a momentum transfer exceeding the Debye mass.Nevertheless the results of both approaches are numerically rather similar.

V. TURBULENT QGP
In this section we consider a Fokker-Planck equation of heavy quarks in a turbulent QGP which is populated with strong chromodynamic fields.The plasma is assumed to be isotropic and translationally invariant both in time and space.The tensor X ij (v), which enters the Fokker-Planck equation (17), is given by Eq. ( 41).The method to obtain Y j (v), which is used in [11], works only for equilibrium plasmas.Therefore, we will refer to the relation (25).However, it implicitly assumes that in the long time limit the system of heavy quarks described by the Fokker-Planck equation reaches a state of thermal equilibrium with temperature T .First of all, the value of T is, in principle, unknown and one needs additional arguments to determine it.There is also a more important problem -it is unclear what are the properties of X ij (v) for which the assumption of equilibrium makes sense.If, for example, there are only magnetic fields in the plasma, X ij (v) is purely transverse and T Y i (v) = 0. Consequently, the Fokker-Planck equation reads and any stationary isotropic function n(p) solves the equation because ∇ p n(p) ∼ p.Therefore, pure magnetic fields do not drive the system to the thermal equilibrium, as expected.In spite of these concerns, we will use the relation (25) to get Y j (v).

A. Gaussian correlation functions of independent E and B fields
We start with a simple model proposed in [14] where the correlation functions of electric and magnetic fields are chosen in the following Gaussian form The correlation lengths σ t , σ r and the parameters M E , M B of dimension mass to the fourth power will be discussed later on.The fluctuation spectrum obtained from the correlation function ( 52) is also Gaussian and it equals Substituting the correlation functions (52), (53), and (54) into Eq.( 41), the tensor X ij (v) is found as which gives When v 2 σ 2 r /σ 2 t and v 2 M E /M B , the coefficients X L (v) and X T (v) are, as in the equilibrium case, equal to each other and independent of v, and The Fokker-Planck equation is then of the form (27).

B. Gaussian correlation function of vector potentials
Since the electric and magnetic fields are, in general, coupled to each other, the functions E i a E j b , E i a B j b , B i a E j b , and B i a B j b are not fully independent from each other.The electric and magnetic fields are automatically related to each other if one postulates the correlation function of the four-potential and then computes the correlation functions of the E− and B−fields.In this section we follow this path.Specifically, we assume the Gaussian correlation function of the vector potential The parameter M A of the dimension mass squared will be discussed later on.The fluctuation spectrum of the potential equals We further choose the radiation gauge and the electric and magnetic fields are obtained in the linear regime as To get the fluctuation spectra b ω,k from the spectrum (61) we refer to the relation (B7) derived in Appendix B. Using Eq. ( 63), the fluctuation spectra are found as Substituting the formulas (64), (65), and (66) into Eq.( 41), one gets after some manipulations which gives When v 2 σ 2 r /σ 2 t , we find again, as in the equilibrium case, that the coefficients X L (v), X T (v) are equal to each other, but they do depend on v. Specifically, In contrast to the equilibrium result and Eq. ( 59), X L (v) and X T (v) vanish as v → 0.

C. Stationary power spectrum of vector potential
When the momentum distribution of plasma constituents is anisotropic, the system is unstable due to the Weibel instability, see the review [2], and a strong chromomagnetic field is generated spontaneously.As shown in the numerical study [21], the fluctuation spectrum of the soft fields becomes up to the wave number k max stationary after a sufficiently long time and the spectrum decays with wave vector as k −2 .Inspired by these findings we choose, as previously, the radiation gauge (62) and consider the following spectrum where the parameters M, µ and k max , which are all of the dimension of mass, will be determined later on.We note that k max is of order of the Debye mass and that the small but nonzero parameter µ is introduced to eliminate the infrared divergence of the expression (71).In contrast to the study [21], there are only zero frequency modes in the spectrum (71), and consequently the correlators involving electric field vanish.The only non-vanishing correlator is Using the formula (36), the tensor X ij (v) is found as The coefficient X L (v) given by Eq. ( 23) vanishes because of transversality of magnetic field while X T (v) equals After taking the elementary integrals, one obtains As seen, X T (v) vanishes when v → 0 and the formula (75) simplifies to when k max µ.In Tables I and II there are collected the coefficients X L (v), X T (v) and corresponding transport coefficients dE dx , q obtained in the models of turbulent plasma discussed in the subsections V A, V B and V C. The models are called 'Gauss E & B', 'Gauss A' and 'Stationary A', respectively.We observe that dE dx is not always reliable because there is no guarantee that heavy quarks evolve in turbulent plasmas towards the state of thermodynamical equilibrium as required by the relation ( 26) that is used.We also note that the formula of q in the 'Gauss E & B' model, which holds not for a test quark with any v but for a test gluon with v = 1, was first obtained in [15].

VI. EQUILIBRIUM VS. ISOTROPIC TURBULENT PLASMA
We are going to compare the equilibrium plasma to the turbulent one at the same energy density.The energy density of the equilibrium quark-gluon plasma of N f massless flavors equals The density of energy accumulated in chromodynamic fields is expressed as If the fluctuation spectra E i a E i a ω, k and B i a B i a ω, k are of the Gaussian form (55), the energy density equals When the fluctuation spectra E i a E i a ω, k and B i a B i a ω, k are given by Eqs. ( 64) and (65), we have Finally, if the fluctuation spectrum A i a A i a ω, k is of the form (71), the purely magnetic energy density equals where we assume that k max µ.We note that in the weak coupling limit the magnitudes of electric and magnetic fields in turbulent plasmas, which are according to the above estimates of the order T 2 , are much larger than in equilibrium plasmas where the fields are typically of the order g 2 T 2 .This is the main reason why the turbulent plasmas, which are discussed here, are qualitatively different than the equilibrium one.
To reduce the number of parameters of our models of turbulent plasma we adopt the simplifying assumptions that M E = M B and σ t = σ r .Then, when the energy density of the fields in turbulent plasma (79), ( 80) and (81) equals the energy density of the equilibrium QGP (77), the parameters M E , M A and M are equal to TABLE I: The coefficients XL(v) and XT (v) in the models of turbulent plasmas FIG. 5: The energy loss as a function of v in the models of turbulent plasma.
FIG. 6: The momentum broadening q in the models of turbulent plasma.
We set N c = 3 and N f = 2.With the simplifying relations M E = M B and σ t = σ r and the equalities (82), one obtains the parametric estimates of the coefficient X T (v) which are shown in the second column of Table III where we have additionally assumed that v 2 1.Since the parameters σ t , σ r determine the field correlation length in time and space, they are of order of the screening length and thus we choose σ t = σ r = m −1 D .Similarly, k max is a soft momentum and following [21] we set k max = 5m D .Using the formula (A7), we get The parametric estimates of X T (v), which take into account the relations (83), are shown in the third column of Table III.As seen, the interaction of heavy quarks is much stronger in the turbulent plasma than in the equilibrium one if the plasma is truly weakly coupled with g 1.The effect disappears when the coupling constant α s ≡ g 2 /4π is of realistic value and gives g close to unity or even bigger.
One asks whether the strong magnetic fields discussed above are stable or may be they rapidly decay due to the Nielsen-Olesen instability [22] which has been discussed in the context of quark-gluon plasma in [23,24].Strictly speaking, the Nielsen-Olsen instability occurs in systems with a homogeneous magnetic field because the system's energy can be reduced by particles circulating in the magnetic fields.The unstable mode grows as e γt with γ = √ gB.When the field is inhomogeneous with the wave vector k, the decrement of growth is reduced to γ = gB − k 2 .The instability disappears when gB < k 2 which physically means that the length of inhomogeneity is shorter than the Larmor radius.As discussed above, B ∼ T 2 and k ∼ gT in the turbulent plasma and thus the Nielsen-Olsen instability is unavoidable as long as g 1.The lifetime of strong magnetic is then of order ( √ g T ) −1 .In Sec.VII we discuss the glasma [1] which is of our main interest.We show that the instability is absent in the glasma because the field correlation length is of the order of the inverse stauration scale Q s .In Figs. 3 and 4 we present how X L (v) and X T (v) depend on the velocity v in the plasma models under consideration.The coupling constant and temperature are α s = 0.1 and T = 0.5 GeV.Figs. 5 and 6 show the corresponding energy loss and momentum broadening as functions of v. Since the energy density of the turbulent plasma is fixed, the dx and q in the models of turbulent plasmas.
Stationary A 0 temperature T , which enters Eq. ( 25) or (26), is fixed as well.The coefficient q is computed for v ≥ 0.3 because it diverges as v → 0 in the model 'Gauss E & B'.However, q is of no physical relevance for small velocities.As seen, the coefficients X L (v) and X T (v) and consequently dE/dx and q are rather different in the equilibrium and turbulent plasmas -not only the magnitudes differ but the dependences on the heavy quark's velocity are different.The interaction of heavy quarks is particularly strong in the model 'Gauss E & B'.One wonders why the models 'Gauss E & B' and 'Gauss A', which at first glance are quite similar, give rather different results shown Figs. 3 and 4. Comparing the fluctuation spectra (64) and ( 65) to (55), one observes that the low frequency and long wavelength fields are suppressed in the 'Gauss A' model.These fields appear to contribute more effectively to the transport coefficients than the high frequency and short wavelength fields.

VII. GLASMA
When relativistic heavy ions collide color charges of partons confined in the colliding nuclei generate strong chromodynamic fields right after the collision.Since the system of infinitely contracted nuclei moving against each other with the speed of light is boost invariant, so is the configuration of generated fields.As shown in the detailed analytic study [25], the chromoelectric and chromomagnetic fields spanned between the receding nuclei are initially only parallel to the beam direction.Transverse field components start to develop later on.We focus here on the longitudinal fields which dominate the glasma's dynamics and are invariant under Lorentz transformations along the beam direction identified with the axis z.

A. Field correlation functions
Since the electric and magnetic fields are expressed in the Abelian limit through the four-potential according to Eq. ( 63), the potential generating the fields only along the axis z is of the form that is the 0 and z components of A µ depend on t and z while the x and y components on x and y.The electric field corresponding to the potential (84) depends on t and z while the magnetic field on x and y.The electric field E z a is boost invariant if it depends on t and z only through the proper time τ ≡ √ t 2 − z 2 which is the Lorentz scalar.However, we do not require the boost invariance of the fields but of the field correlators.To write down a general expression of the electric field correlator E z a (t 1 , z 1 ) E z b (t 2 , z 2 ) , we introduce the variables and we note that the proper times τ i and space-time rapidities η i are well defined only for the time-like two-vectors (t i , z i ).The boost invariant correlation function of the electric fields is assumed to be where Θ(t 2 i − z 2 i ) is the Heaviside step function and f E (τ, η) is an arbitrary function.Since τ i and η i are well defined only for t 2 i ≥ z 2 i , we require that the correlation function (86) vanish when the space-time points (t i , z i ) are localized beyond the light cone of the point (0, 0).We observe that the correlation function, which depends not only on τ 1 − τ 2 but on both τ 1 and τ 2 , is also boost invariant but we assume the translation invariance in the τ variable just for simplicity.
The boost invariant magnetic field correlator is chosen to be that is the plasma is assumed to be translationally invariant in the x and y directions.We note that the fields E z a (t, z) and B z a (x, y) are completely decoupled from each other in the Abelian limit and thus the mixed correlator E z a (t 1 , z 1 ) B z b (t 2 , z 2 ) is expected to be small or vanish.We introduce the unit vector n = (0, 0, 1) along the axis z and the correlation functions of electric and magnetic fields are written as The correlators are chosen to be of the Gaussian form with the real positive parameters ME , MB , σ τ , σ η and σ T to be determined later on.

B. Computation of X and Y
Substituting the correlation function ( 89) with (91) into Eq.( 16), the magnetic contribution to the tensor X ij (v) is found as where Since the Fokker-Planck equation is derived in the long time limit, we assume that t σ T and the integral becomes Gaussian.Thus, we get where y is the quark velocity perpendicular to n.The electric contribution to the tensor X ij (v) is where FIG. 7: The coefficient X zz E (v) divided by the approximate expression (98) for t = 10στ in the left panel and for t = 30στ in the right panel.
The integral from Eq. ( 95) with the correlation function ( 90) is difficult to compute analytically.The problem greatly simplifies in the long time limit.When t z and t σ τ one easily finds which is independent of σ η .We have checked numerically a quality of the approximation (98).In Fig. 7 we show the numerically computed coefficient X zz E (v) divided by the approximate expression (98) for t = 10σ τ and t = 30σ τ as a function of heavy-quark velocity v z .The coefficient is independent of v x and v y .The correlation length in rapidity is chosen as σ η = 1 but for a bigger σ η the approximation (98) is even more accurate.We conclude that the approximation (98) works pretty well and we write the tensor X ij (v), which includes the magnetic and electric contributions, as Using the relation (25), one finds that the vector Y i (v) equals The temperature will be estimated in the next section but, as already discussed at the beginning of Sec.V, the result (100) should be treated with a reservation.

C. Estimates of parameters
To get approximate values of the parameters, which characterize the glasma, let us estimate the density of energy released in relativistic heavy-ion collisions.When one deals with the central collisions of two nuclei of mass number A, the energy density in the center-of-mass frame of colliding nuclei is roughly where √ s is the center-of-mass energy of nucleon-nucleon collision, c inel is the inelasticity parameter which gives the fraction of all accessible energy going to particle production, R A is the radius of colliding nuclei and l is a length of the cylinder where the energy is released.Assuming that c inel = 0.5 [26] and taking A = 200, R A = 7 fm and l = 1 fm, one obtains ε coll ≈ 3.25 TeV/fm 3 for √ s = 5 TeV which is the energy of Pb-Pb collisions at LHC of the 2015 run.According to Eq. (78) the density of energy accumulated in chromodynamic fields equals where the explicit form of the correlation functions (90) and ( 91) have been used.The field energy density is controlled by the parameters ME and MB but is independent of the correlations lengths σ T , σ τ and σ η .Assuming ME = MB , as suggested [25], we get ε field = 8 ME for N c = 3. Requiring that ε field = ε coll , the parameter ME is estimated as To use the formula (25) to obtain the coefficient Y i (v), one needs a temperature of the equilibrated quark-gluon plasma of the same energy density as the glasma.Using the formula (77), the equation ε field = ε QGP provides where N c = 3 and N f = 2. Within the CGC approach the strong longitudinal chromodynamic fields are screened on transverse distance which is of the order of the inverse saturation momentum Q s .Since Q s is estimated as 2 GeV [1], we choose the transverse correlation σ T = Q −1 s = 0.5 GeV −1 .We also assume that σ τ = σ T .The remaining parameter is the coupling constant which, as previously, is chosen to be α s = 0.1.As already mention in Sec.VI, the glasma magnetic field is stable against the Nielsen-Olesen instability, because the inhomogeneity length of order Q −1 s is not much bigger than the Larmor radius (gB) −1/2 , Q s is even bigger than M 1/4 B .Using the formulas (100) and (99) combined with Eqs. ( 31) and (34), we obtain the following estimates of the energy loss and momentum broadening where θ is the angle between v and n.Because the electric field is along the beam axis, the collisional energy loss is maximal when the heavy quark moves along the axis, v n.The maximal momentum broadening occurs when v ⊥ n.Let us compare the numerical values (105) and (106) with those which are required to properly model experimental data on the charm meson suppression.The collisional energy loss and momentum broadening of a charm quark with 10 GeV momentum in the plasma of the temperature from the interval 0.35 -0.5 GeV are estimated [3] as As seen, the values (105) and (106) can be significantly larger than (107) and (108), suggesting that in spite of a short lifetime of the glasma it can provide a significant contribution to the collisional and radiative energy loss to heavy quarks from relativistic heavy-ion collisions.Consequently, the effect should be included in the phenomenology of jet quenching.

VIII. SUMMARY, CONCLUSIONS AND OUTLOOK
Applying the so-called quasi-linear theory [9,10], we have derived a general form of the Fokker-Planck equation of heavy quarks embedded in the plasma of light quarks and gluons.Since the interaction is taken into account through the correlation functions of chromodynamic fields, heavy quarks are seen as interacting not with plasma constituents but rather with fields present in the plasma.At first we have obtained the explicit form of the equation for the case of equilibrium plasma which was studied long ago [6] using the standard method where the Fokker-Planck equation simply approximates the Boltzmann one.Although our approach is noticeably different, the Fokker-Planck equation we obtained agrees with the standard one [6].
In the second part of the paper, the method developed for the equilibrium plasma has been applied to the turbulent plasma populated with strong fields.The parametric estimate shows that the interaction of heavy quarks with the turbulent plasma is much stronger than with the equilibrium one of the same energy density if the coupling constant is truly small.The effect is less prominent for a realistic value of the coupling constant and the difference depends on characteristics of the plasma fields.Within the 'Gaussian E & B' model both the energy loss and momentum broadening are significantly bigger than the equilibrium results.A dependence of dE/dx and q on heavy quark velocity also strongly depends on how the turbulent plasma is modeled.
The third part of our study is devoted to the glasma from the earliest stage of relativistic heavy-ion collisions.Assuming that there are chromoelectric and chromomagnetic fields only along the beam direction we have derived the appropriate Fokker-Planc equation.We have also shown that in spite of its short lifetime the glasma can provide a significant contribution to the collisional and radiative energy loss of heavy quarks.
Our findings clearly suggest a direction of further work.We need a more realistic model of turbulent QCD plasma from relativistic heavy-ion collisions.In contrast to the simple model discussed here, a temporal evolution of the glasma has to be taken into account and the fields cannot be purely longitudinal.The CGC studies [1] and, in particular, the analytic analysis [25] provide a very good guidance to build up such a model.

FIG. 1 :FIG. 2 :
FIG. 1: The equilibrium coefficients XL(v) and XT (v) as functions of the velocity v. FIG.2: The equilibrium coefficients XL(v) and XT (v) as functions of the temperature T .

FIG. 3 :FIG. 4 :
FIG. 3: The coefficient XL(v) as a function of the velocity v in the equilibrium plasma and the models of turbulent plasma.FIG.4: The coefficient XT (v) as a function of the velocity v in the equilibrium plasma and the models of turbulent plasma.

TABLE II :
The transport coefficients dE

TABLE III :
Parametric estimates of the coefficients XT (v) in the four plasma models under consideration.The third column takes into account that σ −1 t ∼ kmax ∼ gT .