Heavy quark diffusion coefficient in heavy-ion collisions via kinetic theory

We compute the heavy quark momentum diffusion coefficient $\kappa$ using QCD kinetic theory for a system going through bottom-up isotropization in the initial stages of a heavy ion collision. We find that the values of $\kappa$ are within 30% from a thermal system at the same energy density. When matching for other quantities we observe considerably larger deviations. We also observe that the diffusion coefficient in the transverse direction is larger at high occupation numbers, whereas for an underoccupied system the longitudinal diffusion coefficient dominates. The behavior of the diffusion coefficient can be understood on a qualitative level based on the Debye mass $m_D$ and the effective temperature of soft modes $T_*$. Our results for the kinetic evolution of $\kappa$ in different directions can be used in phenomenological descriptions of heavy quark diffusion and quarkonium dynamics to include the impact of pre-equilibrium stages.


I. INTRODUCTION
Several interesting signatures of the quark-gluon plasma, such as quarkonium suppression, are sensitive to the dynamics of heavy quarks in the plasma.Due to their large mass m Q ≫ T , heavy quarks are solely produced in initial hard scattering processes.Furthermore they interact strongly with the medium during the entire evolution.The typical momentum exchange of a heavy quark with the QGP is small (O(gT )) and thus a significant change in momentum requires several scatterings with the medium.This is characteristic of diffusive motion, and hence diffusion properties of heavy quarks within the QGP have been subject to extensive research [1][2][3].The interaction of the heavy quark with the medium is conventionally quantified by the momentum diffusion coefficient κ.
The diffusion coefficient is not only interesting in its own right as a property of the medium-it also plays an instrumental role in understanding and describing the evolution of quarkonia using various phenomenological frameworks [2,3], including e.g., the open quantum systems approach [4][5][6][7][8].
The aim of this paper is to bridge this gap and to study the heavy quark momentum diffusion coefficient κ during the bottom-up isotropization process [34].With this aim in mind, we use QCD effective kinetic theory with the same initial conditions and numerical setup as in [35] for the limit of an infinitely massive stationary quark.Furthermore, our aim is to understand how κ out of equilibrium compares to similar computations in thermal equilibrium.To this end, we compare values of κ to thermal systems with the same screening mass, energy density or effective infrared (IR) temperature.We establish that out of these options κ in the pre-equilibrium system is closest to the thermal value at the same energy density.We will also study the transverse (to the direction of the expansion), and longitudinal (along the beam direction) diffusion coefficients separately.We observe that there is a dynamical hierarchy between the longitudinal and transverse directions, which has a natural explanation in terms of the different stages of the bottom-up thermalization process.We will also derive a set of parametric estimates that explain the evolution of the matched quantities and κ relative to their equilibrium value in terms of the dynamical occupation number and anisotropy of the system.This paper is structured as follows.In Sec.II we discuss the kinetic theory framework, how we extract the heavy quark momentum diffusion coefficient and introduce our initial conditions and the relevant physical scales m D and T * .We continue in Sec.III by presenting our results for the diffusion coefficient during the initial stages in heavy-ion collisions and discuss the diffusion coefficient in the transverse and longitudinal directions separately.In Sec.IV we compare our results to those obtained in lattice and glasma simulations.Finally we conclude in Sec.V.
In this paper we use the same numerical setup as in Ref. [35], with the aim of reproducing the bottom-up isotropization process [34] using the numerical framework of QCD effective kinetic theory (EKT) [36].Kinetic theory is a quasiparticle description of the quark-gluon plasma in terms of its individual constituents, quarks and gluons.Here we consider only gluons for which the central dynamical quantity is the quasiparticle distribution function This is the phase space density of gluons averaged over their degrees of freedom ν g = 2d A , including d A = N 2 c −1 colors and two spins.We assume that the distributions are independent of transverse coordinate, boost invariant and do not depend on spin.The time-evolution of the distribution is given by the Boltzmann equation Here the C terms describe the number of scatterings/splittings per unit time for momentum state p.These processes involve the possibility to scatter/split into or out of a state.Effective particle 1 to 2 splittings are encoded in C 1↔2 and 2 to 2 scatterings are described by C 2↔2 .The boost invariant expansion is treated as an effective scattering term, taking the simple form [37] C These assumptions leave our distributions f dependent only on the magnitude of the momentum p and on the polar angle cos θ = ẑ • p, i.e f (p) = f (p, cos θ p ).Our spherical coordinate system is defined such that Momenta in the x and y directions are given in terms of the azimuthal angle ϕ The 2 to 2 scattering term reads [36] C where dΓ PS is a phase space integration over the incoming and outgoing momenta p, p ′ , k and k ′ .The matrix Here λ = g 2 N c is the 't Hooft coupling.In ( 9) the division by the Mandelstam variables t and u leads to infrared divergences.These can be regulated by including contributions arising from medium modifications.In the employed numerical framework, these are approximated by substituting the denominators in (9) with with q = |q| and q = p ′ − p, and symmetrically for u.Note that the limit t → 0 with s fixed is reached only when q 2 → 0. We follow the prescription of [38] and use Here m is the asymptotic mass in the gluon dispersion relation at high momentum that is perturbatively related to the Debye mass m D as The final state Bose enhancement in (8) is governed by the factors (1 + f ).The first term in ( 8) is a loss term for particles scattering out of the state and the second term is the term associated with scattering into the state.The phase space integration measure is given by with ω = q 0 , x q = cos θ q , p = dp 3 2p 0 (2π) 3 and p 0 = p.The simplification above is achieved by changing variables from p ′ to q and using delta functions to eliminate the k ′ integral.The residual three-dimensional integrals can be written as integrals over four-momenta, introducing delta functions that constrain the energy component.The integration limits are found by inspecting the energy delta functions.One also makes use of the fact that the system is azimuthally symmetric.
The collinear 1 to 2 splitting rate C 1↔2 is given by Due to collinearity all momenta are pointing in the same direction.We take the particle with momentum p to be the particle that splits, and hence p ′ = p − k ′ .Here γ parametrizes the differential rate for the splitting processes and is given by The function F is computed from Here n is a unit vector in the direction of the splitting or merging hard particle(s).The vector h is perpendicular to n, and can be related to the transverse momentum [39].In deriving (16) the Wightman correlation function is evaluated by making use of the sum rule derived in [40] using an isotropic screening approximation.The energy difference in ( 16) is given by The integral equation ( 16) is solved by using the numerical method described in [41].For a more comprehensive description of EKT we refer the reader to [36] and we discuss our discretization procedure in A.

B. Numerical implementation and discretization effects
Our numerical framework has several discretization parameters, such as minimum and maximum momenta stored for the distribution function or the number of bins in the angular discretization.In our simulations, it turns out that the dominant source of discretization effects is the infrared cutoff of the particle distribution function p min .The reason for this is that due to the longitudinal expansion, the typical momenta become smaller over time.For instance, the dimensionless ratio of the infrared cutoff and the temperature p min /T (τ ) grows with time.In practice, this means that the thermal equilibrium that the system approaches will also have discretization effects.We will take these into account by replacing the continuum thermal quantities with ones calculated from a thermal distribution with the same cutoff as the EKT simulation.This correction enables us to compare our nonequilibrium results to the thermal state that the system is actually approaching.We will apply this correction to all observables when comparing with equilibrium quantities unless stated otherwise.This effect is discussed in more detail in Appendix A 3. Our discretization framework and parameters are discussed in more detail in Appendix A 1.

C. Initial conditions: bottom-up thermalization
Our initial conditions are chosen to correspond to the bottom-up isotropization scenario, as implemented in Ref. [35].The initial distribution function at the time In this paper we consider the same two sets of initial parameters with different anisotropies ξ as in Ref. [35], given by In the following we will refer to these only by the value of the anisotropy parameter ξ.In our figures the initial condition with ξ = 10 is always represented by full lines.When ξ = 4 results are shown for comparison, we use dash-dotted transparent lines.

D. Heavy quark diffusion coefficient
The expression for the heavy quark momentum diffusion coefficient κ for pure glue QCD has been originally derived in [9].Keeping only the leading order contributions in 1 /M, where M is the mass of the heavy quark, we obtain the diffusion coefficient as The heavy quark mass is taken to be larger than any other scale in the system (e.g.M ≫ T, Q s , m D ).The inand outgoing heavy quark momenta are given by p, p ′ , the momenta of gluons in the plasma are labeled by k, k ′ , and q = p − p ′ is the momentum transfer to the heavy quark.The gluon distribution function f is obtained from solving the Boltzmann equation (2).A final state Bose-enhancement factor appears in (21) similarly as in the Boltzmann equation.We set p 0 = p ′0 = M for the heavy quark.For the gluons and To leading order in the coupling, the dominant contribution to κ is given by t-channel gluon exchange [9], while other diagrams are suppressed by inverse powers of the heavy quark mass.The corresponding matrix element is then given by where C H = (N 2 c −1) /2Nc.Note that due to the heavy quark mass, the HTL propagator reduces in this case to just Debye screening [9], unlike for massless particles in Eq. (10).
The resulting expression |M κ | 2 /M 2 entering the Boltzmann equation is independent of the heavy quark mass.In our simulations, κ is computed using Eq.(21).The discretization of (21) and a detailed description of how transverse and longitudinal diffusion coefficients are extracted, are explained in more detail in the Appendix A 2. For more details on the diffusion coefficient we refer the reader to [9] and to our previous work [30].
Perturbatively, in thermal equilibrium the screening mass (11) becomes In our kinetic simulations, we compute m D according to Eq. ( 11) using the nonequilibrium distribution f (p).The screening mass is also affected by the IR regulator p min .Thus, when comparing m D to its thermal value, we will apply a discretization correction to it.This is discussed in more detail in Appendix A 3 a.
For an isotropic distribution in the continuum, the diffusion coefficient (21) becomes This expression is very easy to evaluate in thermal equilibrium using a Bose-Einstein distribution for f (p) and the thermal result (23) for m D .We take into account the impact of a finite p min on the value of κ in a thermal system by calculating it using (24), where the screening mass now depends on the infrared cutoff p min according to Eq. (A12).We will not regulate the momentum transfers by this parameter, since the main p min dependence seems to enter κ through m D .This procedure is described in more detail in Appendix A 3 c, and our final expression for thermal κ in the presence of the IR regulator p min is given by Eq. (A17).

E. Other observables
Here we describe the observables that we are extracting during the nonequilibrium evolution in addition to the heavy quark diffusion coefficient.

Expectation values
In general, we use the following definition for the expectation value of observable X ⟨X⟩ =

Energy density
The energy density corresponds to the first moment of the distribution where ν g = 2 N 2 c − 1 for pure glue QCD.

Temperature
There is no unambiguous way to define a temperature out of equilibrium.One option is to formulate it effectively in terms of the (time-dependent) energy density which agrees with the temperature in thermal equilibrium.For ideal Bjorken hydrodynamics, the temperature scales as This scaling is also expected to hold for T ε at sufficiently late times due to the approximate τ −4/3 scaling of the energy density in the expanding system.
The effective infrared temperature T * is given by where Similarly to κ and m D , T * is also affected by the IR regulator p min .Thus we compute the corresponding equilibrium value by taking into account these effects both in m D and in the momentum integral.This is discussed in more detail in Appendix A 3 b.In thermal equilibrium in the presence of the IR regulator, T * is given by Eq. (A15).Slightly different definitions of T * , based on other integral moments, have been also introduced in the literature [42], but we do not consider them here.

Energy momentum tensor
The components of the energy momentum tensor are obtained as moments of the distribution function by The components relevant for this paper are They are connected to the longitudinal and transverse pressure by The temporal component of the energy-momentum tensor corresponds to the energy density T 00 ≡ ε.

III. RESULTS
We start by discussing the bottom-up isotropization process and how it is reproduced in our simulations in Sec.III A. We will also discuss how we highlight different stages of the isotropization process.In Sec.III B we will explain how we compare our κ results out of equilibrium FIG. 1. Occupation number as a function of anisotropy as in [35].The markers indicate different stages of the bottom up thermalization as described in the text.The color coding of the lines indicates the coupling used in the simulation, and all figures use the same color coding.Full lines correspond to the ξ = 10 initial condition (19), dash-dotted lines to the ξ = 4 initial condition (20).
to thermal values, and especially how we choose the corresponding matching scales.The comparison is then carried out in Sec.III C, first without distinguishing between directions.Then in Sec.III D we discuss the transverse and longitudinal diffusion coefficients separately.In order to better understand the observed time-evolution of κ, we will consider how the matching scales evolve during the bottom-up isotropization compared to equilibrium in Sec.III E. Finally, we will derive simple parametric estimates that can be used to explain our results in Sec.III F.

A. Different stages of the bottom up thermalization scenario
Our kinetic theory simulations follow the different phases of the bottom-up thermalization scenario [34].It consists of the following stages: during the first stage the overoccupied gluons become more dilute as the system expands, and consequently the occupation number of the hard gluons becomes of the order of unity . At this point the system is no longer describable with classical fields.In the second stage hard gluons radiate softer gluons, creating a soft thermal bath.Remarkably, at the end of this process the hard gluons become underoccupied f h ∼ α s .This occurs at the timescales Q s τ ∼ α −5 /2 s .In the final stage, the hard particles lose their energy to the soft thermal bath.For a review on the thermalization processes see e.g.[43,44].
In this scenario, the system is expected to thermalize parametrically on a timescale of the order of [34] τ BMSS = α −13 /5 s /Qs, where α s = λ 4πNc .We will therefore use this quantity to rescale the time in our figures.Note that an alternative time scale, the hydrodynamical relaxation time τ R = 4πη/s(λ) T with shear viscosity η and entropy density s, is also often used to rescale the time variable.We will investigate the universality of these time scales for different observables and couplings in a separate paper, while employing only τ BMSS in the present work.
The bottom-up thermalization process is shown in Fig. 1 in terms of the anisotropy P T /P L and the mean occupation number ⟨pλf ⟩ /⟨p⟩ for different couplings as in Ref. [35].In order to illustrate how observables behave during different stages of the thermalization process, we have placed three time markers on the curves in Fig. 1.The first marker (star) is placed during the highly occupied regime, when f ∼ 1 /λ.For smaller values of the coupling this corresponds to maximal anisotropy.However, for large couplings the first stage of the evolution proceeds differently, and the anisotropy does not increase initially.Hence, we have chosen the occupancy as the criterion for the time marker instead of maximum anisotropy.The second marker (circle) is inserted at the minimum occupancy, which in the bottom-up thermalization scenario is expected to be f ∼ α s .The third marker (triangle) is placed at P T /P L = 2.The purpose of this is to illustrate when the system is approximately close to equilibrium.We will use the same markers in other figures throughout this paper to allow the reader to connect the timeevolution of observables to the stages of the bottom-up scenario.
The curves for different values of the coupling λ also use the same color coding in all figures throughout this paper.The initial conditions with ξ = 10 from ( 19) are shown as full lines.For comparison, we will often add curves for the initial conditions with ξ = 4 in (20) as more transparent dash-dotted lines.
We see from Fig. 1 that for the extremely anisotropic initial condition (ξ = 10) the bottom-up picture is better realized at smaller values of the coupling.For intermediate couplings λ = 2, 5, 10 the system does not experience an initial growth in anisotropy.Instead, the system takes a more straightforward path to thermal equilibrium, without resolving of the different stages of the bottom-up picture in detail.However, the third stage of the scenario is still visible and emerges after the circle marker.

B. Comparing nonequilibrium to equilibrium
Our aim in this paper is to calculate the heavy quark momentum diffusion coefficient κ during the hydrodynamization process in order to eventually assess the importance and impact of the initial nonequilibrium evolu- tion on heavy quark observables.To facilitate the quantitative interpretation we will mostly present our results as ratios to the thermal equilibrium values.There is no unique method to compare equilibrium and nonequilibrium systems.A reasonable way to construct such We observe that the transverse coefficient is enhanced compared to the longitudinal one during the initial evolution.When the system becomes underoccupied, this ordering reverses.Finally when the system reaches approximate equilibrium the ratio approaches unity.
a comparison is to match a dimensionful observable to its thermal counterpart.We will consider three different quantities -the energy density ε (which leads to the temperature defined in ( 27)), the temperature of infrared modes T * , and the screening mass m D .We have chosen these observables because they are straightforward to compute both in equilibrium and out of equilibrium and are physical scales that play an important role in transport phenomena, as discussed in Sec.III E. However, we would like to emphasize that the matching could be done with respect to any dimensionful observable that can be defined in and out of equilibrium.The matching is performed by choosing the temperature of the thermal distribution to reproduce either the energy density ε, the Debye mass m D or the effective temperature of the soft modes T * calculated from the EKT simulation, and then calculating the diffusion coefficient κ with this thermal distribution.

C. Diffusion coefficient in and out of equilibrium
Figure 2 shows our result for the heavy quark diffusion coefficient during bottom-up thermalization.we observe an agreement between the equilibrium and thermal dynamics within 30% even at early times, while at late times the system thermalizes and the ratio thus approaches unity.This is one of the main results in this paper.Matching the screening mass m D (2 middle) or the effective infrared temperature T * (2 bottom), κ is much further away from the thermal system.The main observation is that at early times the non-equilibrium diffusion coefficient is considerably larger than the equilibrium coefficient for the same m D , and considerably smaller than that for the same T * .

D. Transverse and longitudinal diffusion coefficient
During most of the bottom up thermalization procedure the system is highly anisotropic.This could have a significant effect on experimental observables sensitive to the initial stage of the evolution.It is therefore interesting to study also the anisotropy of the diffusion coefficient, which we parametrize in terms of the ratio of the transverse diffusion coefficient κ T to the longitudinal one κ z .Our results for this quantity are shown in Fig. 3.In the initial overoccupied phase the transverse diffusion coefficient dominates.This can naturally be explained by the Bose enhancement, which at the early overoccupied and highly anisotropic stage benefits scatterings where only transverse momentum is exchanged.As the system becomes underoccupied, the longitudinal diffusion coefficient becomes dominant.When the system approaches equilibrium the two become equal again, as expected due to the emerging isotropization.Thus the hierarchy of the diffusion coefficients has a natural explanation in terms of the stages of the bottom-up thermalization scenario.The ordering of the coefficients κ T < κ z after the star marker, i.e., after the initially large occupancies, is also in line with what is observed using squeezed thermal distributions [31] and results from the momentum anisotropy of the system.
We can compare the pre-equilibrium system to the thermal one using the same three matching procedures as in Fig. 2 to the transverse and longitudinal diffusion coefficients separately.The results are shown in Fig. 4, and show qualitatively similar results as for the full coefficient κ.Similarly as in Fig. 2 we find that the best way to compare our nonequilibrium results to a thermal system is by matching for the same energy density ε.This way the deviations from equilibrium turn out to be not larger than approximately 40 %.For the same screening mass m D we observe very large deviations after the initially highly occupied regime, i.e., after the star time marker.In particular, we observe deviations up to a factor of 4 for the transverse and 6 for the longitudinal coefficient, respectively, for the smallest coupling.

E. Time-evolution of the relevant scales
To better understand the time dependence of κ, let us now discuss the evolution of relevant scales that it depends on.The different time dependences that are relevant for this discussion are summarized in Figs. 5, 6  and 7.As a first step, we start with the time dependence of the occupation number of the hard modes, shown in Fig. 5.It starts in an overoccupied state, then becomes underoccupied before reaching a thermal value.
Figure 6 shows the time dependence of the energy density.In the upper panel, the initial value ε ∼ 1/λ is scaled out by multiplying with the coupling λ so that all the curves start at the same value by construction.The time dependence is additionally divided by the ideal hydrodynamical behavior τ −4/3 .We see that during most of the bottom-up thermalization, the energy density decreases as ε ∼ 1/τ until it approaches the hydrodynamical behavior at the triangle marker near the end of the evolution.The bottom panel shows the temperature T ε extracted from the energy density and the temperature in ideal hydrodynamics T id , extrapolated backwards from the final state.The latter is computed by first matching it with T ε at the triangle time marker of near isotropy when the system exhibits an almost hydrodynamical behavior, and propagating backwards in time as in T id ∼ τ −1 /3 , as discussed in [42].We find that the temperature deviates from the ideal hydrodynamics temperature at most by a factor of 2 at early times for the couplings we have considered.However, when the system is underoccupied at the circle time marker, the ideal temperature is in reasonably good agreement (within 20%) with the temperature extracted from the energy density.
Let us now return to the heavy-quark diffusion coeffi- cient.Parametrically, it is given by [30] κ where Λ is the largest momentum scale of the particles the plasma; for a thermal system we have Λ ≈ T .Ignoring the logarithmic contribution, we estimate To understand the behavior of the heavy quark diffusion coefficient, we must thus look at the time dependence of m D and T * which are shown on the top and bottom, respectively, in Fig. 7.We observe that initially m D is enhanced compared to the thermal system at the same energy density.When the system becomes underoccupied, especially for weak couplings, m D gets relatively suppressed.The suppression is mainly driven by a decreasing occupation number, which is discussed in more detail in Sec.III F. For the effective temperature T * we observe a strong enhancement at the early stage, which is understood from the large occupancies encountered initially as also discussed in Sec.III F. When the system becomes underoccupied, i.e., in the vicinity of the circle time marker, T * is already very close to its equilibrium value for the same energy density.The combination of the curves in Fig. 7 explains qualitatively the behavior seen in Fig. 2; at first an enhancement compared to the equilibrium value due to large T * and m D in the overoccupied phase, followed by a suppression mostly resulting from lower values of m D .Let us now construct a parametric estimate that makes this behavior explicit.

F. Understanding the evolution of scales and κ
Let us try to qualitatively understand the time evolution of the matching scales and of the heavy-quark diffusion coefficient in terms of the occupation number and anisotropy.For the purposes of our parametric estimate, we use the following squeezed and scaled distribution where θ is the step function.We determine the typical occupation number f 0 during bottom-up thermalization by taking its value to be the typical occupation number of the hard modes, calculated as f 0 = ⟨pf ⟩ /⟨p⟩.The values of this parameter can be read off Fig. 5.The squeezing parameter δ describes the momentum anisotropy, δ ∼ ⟨pz⟩ ⟨p T ⟩ .We estimate its value during the bottom-up evolution as δ = P L /P T , where the values of the pressure ratio are visible in Fig. 1.
For the purposes of this parametric estimate, we approximate the thermal distribution with the same sim-plified form (41).For this, we obtain the value f 0 by calculating the expectation values for a thermal (Bose-Einstein) distribution (corresponding to the crosses in Fig. 1).The value is independent of λ and reads For a thermal system, we naturally have δ = 1, and we take the scale Q s in ( 41) to be Q s = T .The screening mass m D is given by (11).Performing the integral for the distribution (41), we get Thus, the ratio to the thermal case becomes This parametric estimate is shown as a gray dotted line in Fig. 7 (top) for λ = 1 (as indicated by the color coding) and ξ = 10.We will use these parameters throughout this section for other quantities as well.The time-evolution in Fig. 7 can now be understood using Fig. 1.The decrease in the value of m D when the system evolves from overoccupation to underoccupation is mainly driven by the falling occupation number f 0 .As the system evolves towards thermal equilibrium from the underoccupation regime (circle marker), the value of m D starts to increase.This process is driven by both the growing occupation number and decreasing anisotropy.As can be seen in Fig. 7, this estimate does not describe the evolution quantitatively, but offers a simple qualitative description of the evolution.Possibly our parametric estimate undershoots the actual value of m D in the underoccupied phase (the circle marker) because the modes contributing to the Debye mass are softer and not quite as anisotropic and underoccupied as the modes contributing to the pressure, leading to the factor √ δf 0 in (44) to overestimate the effect.
Similar estimates can be applied to T * given by (29).Using the result for m D , we obtain Thus T * is less sensitive to the anisotropy.Comparing with the thermal estimate, we get This result naturally explains the behavior of T * observed in Fig. 7 (bottom).The initial enhancement of T * is driven by the large occupation number, and the observed enhancement is larger than for m D since T * has a stronger dependence on the occupation number.Due to the constant term, there is no suppression from underoccupation, and T * is less sensitive to δ than m D .Hence, our parametric estimate is better than for m D , and T * does not go below the thermal value in the underoccupied regime, but is already close to unity.Using the estimates for m D and T * , we obtain for the heavy-quark diffusion coefficient Plugging in the values for f 0 and δ during the evolution (calculated from f 0 and δ in Figs. 1 and 5) leads to the gray dotted curve in Fig. 8.We observe that the general trend of enhancement followed by suppression and equilibration is somewhat exaggerated by this parametric estimate.The transparent curves in Fig. 8 show the estimate obtained using Eq. ( 40) with the actual calculated values of m D and T * from Fig. 7.They exhibit the same behavior in a less exaggerated way.

A. Comparing with glasma
For a sensible comparison between our kinetic results and κ from the glasma stage, we first reproduce the energy density of the glasma by choosing Q s = 1.4 GeV as in Ref. [21] at the initial time Q s τ = 1.The same value is also obtained in [45] in order to achieve consistency with the later hydrodynamic evolution.
Figure 9 shows the transverse and longitudinal diffusion coefficients κ Glasma T and κ Glasma z , which correspond to the situation of static quarks in the glasma [46], together with our results denoted by κ T and κ z .The main observation is that during the quasiparticle phase, κ is considerably smaller than during the glasma stage at very early times (note that the glasma results peak at O(10) GeV 2 fm ).Transverse and longitudinal diffusion coefficients behave differently -our result for κ T agrees with the glasma around the star marker signaling large occupation numbers, which is within the overlapping validity range of glasma and EKT.In contrast, the longitudinal coefficients intersect close to the circle marker, and the transition is not smooth.The circle marker indicates that the system is underoccupied, and hence the matching time falls outside of the validity range of classical-statistical simulations underlying the glasma description.
The main reason for this discrepancy is the qualitatively different behavior of κ in the glasma framework.We can always define the derivative of momentum broadening and call it "diffusion coefficient" κ.This, however, does not imply that the underlying behavior is that of diffusion, corresponding to Langevin-type behavior for ⟨p 2 ⟩.In the glasma, this manifests itself in the longitudinal direction, where for large τ momentum broadening turns into momentum narrowing (κ Glasma z becomes negative).Thus, more research is needed to understand how the 9. Comparison of our results with the glasma results from [46].Blue curves correspond to the longitudinal diffusion coefficient and red curves to the transverse diffusion coefficient.Dashed curves correspond to glasma results and solid curves to our results.The brown data point with error bars corresponds to the result of [47] at T = 1.5 Tc.The point is placed at such a value of τ that the nonequilibrium system has the same temperature Tε defined through energy density as the lattice system.transition from nondiffusive to diffusive behavior takes place between the glasma and quasiparticle pictures.
Based on the data presented in Fig. 9, we can also estimate the total effect of momentum broadening during the nonequilibrium evolution.It is given by the integrated diffusion coefficient In the Bjorken hydro limit, we have ε ∼ τ −4/3 and T ε ∼ τ −1/3 .In a rough parametric estimate κ ∼ T 3 , this leads to κ ∼ 1 /τ and ⟨p 2 ⟩ ∼ log(τ ), everything in units of Q s .Thus, we expect the pre-equilibrium kinetic phase of the evolution, roughly 0.1 − 1 fm, to have an equal contribution to heavy quark diffusion as the equilibrium phase, roughly 1 − 10 fm.Integrating over the entire EKT evolution, from τ = 0.14 fm to τ = 1 fm yields an estimate ⟨p 2 ⟩ = 0.9 GeV 2 .Fig. 9 also features a datapoint for the lattice result, which we will discuss in the following section IV B.

B. Comparing with lattice
For comparisons with lattice results, our main reference will be [47], where low-temperature (T = 1.5 T c ) and high-temperature (T = 10 4 T c ) estimates for κ /T 3 are available.Other studies have also obtained comparable results at similar temperatures [12,17,48].
The obvious approach is to compare systems with the same temperature (defined by the energy density).This, 10.Our results for λ = 2, 5, 10 together with the lattice results from [47].The star, circle and triangle markers correspond to the values of κ /T 3 ε at the time when the corresponding phase of the bottom-up evolution is reached.We only consider the static correlator without higher order (magnetic) corrections here.
however, neglects the different treatment of the QCD coupling.Since the lattice calculation is nonperturbative, it will also include effects arising from the running coupling, whereas our EKT calculation uses a fixed coupling.The second approach is to compare the systems with the same (or similar) coupling constant, which can be estimated using at the scale Q = 4πT , N f = 0, and Λ QCD /T c ≈ 2.
Figure 10 shows our EKT results for λ = 2, 5, 10 as well as lattice results at T = 1.5 T c and T = 10 4 T c in terms of the ratio κ /T 3 .For the higher temperature, the above expression gives λ ≈ 2, whereas the lower temperature -while being outside of the scope of the one-loop expression -indicates values roughly around λ ∼ 10.The EKT values are extracted at different stages of the non-equilibrium evolution, marked by the respective markers.We observe that our results for λ = 2, κ /T 3 = 0.024 ± 0.002 are in rough agreement with the lattice estimate 0.02 ≤ κ /T 3 ≤ 0.16 at the larger temperature (T = 10 4 T c ).However, at the lower temperature corresponding to larger couplings, the lattice calculation yields a larger value of κ/T 3 .This larger value is also visible in Fig. 9 where κ is shown at the time of the evolution corresponding to the same energy density.

V. CONCLUSIONS
In this paper, we have extracted the heavy quark momentum diffusion coefficient using effective kinetic the-ory simulations during the bottom-up isotropization process.In order to better quantify the importance of the pre-equilibrium evolution in relation to thermal equilibrium, we have displayed our result as the ratio κ /κ eq to a thermal distribution with the same energy density ε, Debye mass m D , or effective infrared temperature T * .Our main conclusion is that during the kinetic stages of the pre-equilibrium evolution that can be described using effective kinetic theory, i.e., after the initial glasma evolution, the diffusion coefficient is within 30% of the equilibrium value in a thermal system with the same energy density.Thus our results suggest that modeling the pre-equilibrium diffusion coefficient with the equilibrium coefficient computed for a thermal system with the same energy density (Landau matching) is a reasonably good first approximation during the hydrodynamization process.In a more detailed comparison, we find that in the early overoccupied stage, κ is larger than the equilibrium value and driven by the enhancement of both m D and T * .In the later underoccupied phase, κ is below the thermal comparison value due to the suppression of m D .We have tested initial conditions with two different anisotropies, and we have found that our conclusions are not strongly affected by the value of the initial anisotropy parameter.
We have also studied the diffusion coefficient separately for the transverse and longitudinal (beam) directions.The transverse diffusion coefficient is larger than the longitudinal one at the very early stages when the occupation numbers are large, due to the Bose enhancement available for scatterings with a momentum exchange in the transverse direction.When the system becomes underoccupied, the longitudinal diffusion coefficient dominates.When the system approaches thermal equilibrium, the transverse and longitudinal diffusion coefficients become equal as expected.In a similar manner, we have qualitatively explained the evolution of m D , T * and κ using the behavior of typical hard occupation numbers and momentum anisotropy.Moreover, we have compared our results with those obtained in the glasma and lattice frameworks.
On the phenomenological side, our results give a rough estimate for the heavy quark momentum broadening during the kinetic stage of ⟨p 2 ⟩ ≈ 1 GeV 2 and display interesting angular dependence.It would be interesting to study the implications of these results in more detail.One exciting perspective is to include them into the quantum trajectories framework [6,49] to study the impact of initial stages on quarkonium dynamics.It is also possible to address the evolution of other transport coefficients during the initial non-equilibrium dynamics: in a separate paper [50] we use the EKT setup to calculate the evolution of the jet quenching parameter q.These studies open the possibility of assessing the impact of initial stages on other phenomenological observables.innovation by the STRONG-2020 project (grant agreement No. 824093).The content of this article does not reflect the official opinion of the European Union and responsibility for the information and views expressed therein lies entirely with the authors.This work was funded in part by the Knut and Alice Wallenberg foundation, contract number 2017.0036.TL and JP have been supported by the Academy of Finland, by the Centre of Excellence in Quark Matter (project 346324) and project 321840.KB and FL would like to thank the Austrian Science Fund (FWF) for support under project P 34455, and FL is additionally supported by the Doctoral Program W1252-N27 Particles and Interactions.The authors wish to acknowledge CSC -IT Center for Science, Finland, for computational resources.We acknowledge grants of computer capacity from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533 ).The authors wish to acknowledge the Vienna Scientific Cluster (VSC) project 71444 for computational resources.Our numerical framework is the same as in [35].Instead of the distribution function, we use the quantity which represents the particle number per degree of freedom and unit (spatial) volume, as our dynamical degree of freedom.Here the w ij involve the discretization in momentum p ≡ |p| and polar angle cos θ as with the wedge functions defined as This discretization conserves particle number, energy density and ⟨p z ⟩ exactly.
Our numerical framework has in total seven discretization parameters.Two of them correspond to the number of bins used in the momentum ∈ [p min , p max ] and polar angle.There are two parameters associated to the Monte Carlo sampling procedure, one of which determines the number of samples that is needed for the time-evolution, the other one being associated to the sampling of the diffusion coefficient.For the momentum discretization we have parameters which are associated to the minimum and maximum momentum in our momentum grid p min and p max .Finally for the time-evolution we use an adaptive step-size algorithm, which needs an initial value.
We have verified that our numerical results are independent of the discretization parameter set by reproducing the results of [35].The most important dependencies on the discretization parameters are the following: Reproducing the correct behavior during the underoccupied region is very sensitive to the length of the time-step.It turns out that the minimum momentum on the grid is important to understand the approach to equilibriumas the system becomes thermal, it will approach a thermal distribution with an infrared cutoff p min .When we compare our simulations to thermal equilibrium, this has to be taken into account.This effect is described in detail in Appendix A 3. In principle we receive corrections also from the maximum momentum parameter p max (and other discretization parameters).However, in practice the dependence on p min appears to be by far the most important effect.We have also checked that our results do not depend on the accuracy of the discretization of the θ angle.

Discretization of κ
Due to cylindrical symmetry around the beam direction, the distribution function f (p) out of equilibrium is assumed to only depend on the magnitude of the momentum and the polar angle θ, and not on the azimuthal angle ϕ.We start by evaluating the p ′ integral in (21) by making use of the momentum conserving delta function.Furthermore we integrate over the radial component of k ′ using the energy conserving delta function.This leads to where we use the notation x k = cos θ k , x ′ k = cos θ k ′ for the polar angles.The azimuthal angles are denoted by ϕ k and ϕ ′ k .The cosine of the angle between k and k ′ can be written as As can be seen from (A6), the expression depends only on the difference of the azimuthal angles.Thus we can change the variables by 2π 0 dϕ kk ′ and trivially carry out one of the integrals over the azimuthal angles.Further evaluating the integral and using the results above yields 3κ (A7) The momentum transfer q 2 can be written in terms of the integration variables as We also want to study the transverse and longitudinal diffusion coefficients separately.These are given by Here q i = q T , q z , is the momentum transfer in different directions, given by The normalization is such that q 2 = 2q 2 T + q 2 z , and consequently 3κ = 2κ T + κ z .For a thermalized system this corresponds to κ T /κz → 1 and κ T ,z /κ → 1.In our numerical framework the integrals are computed using Monte-Carlo techniques.In (A7) and (A9) we have also discretized the momentum interval with IR and UV cutoffs [p min , p max ].The dominating discretization effects are discussed in more detail in Appendix A 3.

Observables and pmin dependence
Due to the discretization effects, our system does not reach the continuum thermal equilibrium.Instead, it approaches a discretized version of thermal equilibrium.This affects our observables, as illustrated in Fig. 11.We observe that the numerically obtained value for m D divided by its thermal expectation (obtained for a gluonic system with a Bose-Einstein distribution) deviates from unity at late times.However, when we take into account FIG.11.Debye mass extracted during the bottom-up evolution, divided by the value in a continuum thermal system with the same energy density.The dashed curve has been obtained by computing the screening mass in a thermal system with the IR regulator pmin given by (A12).The temperature Tε is extracted from the energy density using Eq. ( 27).The main observation is that the IR cutoff pmin decreases mD by roughly 20% from the continuum value, depending on T /pmin.
the discretization effects in terms of the momentum cutoff p min , we see a nice agreement for λ = 10.A similar agreement is observed for other curves as well, but here we show only one for clarity.In principle we could take into account also other discretization parameters such as p max and the angular discretization, but we find that the dominant discretization contribution arises from p min .We will now discuss the p min dependence in m D , T * and κ in more detail.The curve labeled as thermal in Fig. 11 is obtained using (A12).(A16) Thus in the limit x → 0 we recover the equilibrium relation λT * = λT .
In the comparisons to equilibrium distributions in the main text (e.g. in Fig. 7) these corrections are applied to m D and T * .Consequently, the ratios approach unity when the system thermalizes, contrary to the behavior observed in Fig. 11 without the corrections.

c. Diffusion coefficient κ
For κ the corrections arising from the infrared regulator p min are taken into account as follows.We start with the expression (24) and perform the q integral analytically without an infrared regulator as before.Then we replace m D with the infrared regulated expression m D (p min ) given by (A12).The resulting expression is The remaining k integral in (A17) is evaluated numerically for different values of λ and T using a small independent IR regulator for numerical stability.The results are interpolated in order to find a κ corresponding to an arbitrary temperature within the tabulated range.We have applied these corrections to the results shown in Fig. 2. We observe that the ratio to equilibrium is very close to unity when the system is approximately thermal.
We would like to emphasize that this treatment takes only into account the discretization effects arising from m D .We do not regulate the k and q integrals with the regulator p min .Thus this treatment takes only a subset of the corrections into account.Hence this treatment is not expected to maintain the ratio at unity for infinitely large times.This can be seen for instance in Figs. 2 and even more prominently 8, where the ratio starts to slightly deviate from unity at very late times.
The ratios κ /κ m D ,T * eq comparing to the thermal system at the same m D and T * are obtained using the same data.The temperature T and p min uniquely determine m D given by (A12) and T * given by (A15).

FIG. 2 .
FIG.2.Ratio of κ to the thermal value at the same energy density ε (top), the same screening mass mD (middle) and effective soft mode temperature T * (bottom).We have applied the Savitzky-Golay filter to smoothen the data.The filter is also applied to all the following figures involving κ.

FIG. 3 .
FIG.3.Ratio of the transverse and longitudinal diffusion coefficients during the bottom-up thermalization scenario.We observe that the transverse coefficient is enhanced compared to the longitudinal one during the initial evolution.When the system becomes underoccupied, this ordering reverses.Finally when the system reaches approximate equilibrium the ratio approaches unity.

FIG. 4 .
FIG. 4. Transverse (left panel) and longitudinal (right panel) diffusion coefficients compared to thermal values for the same ε (top row), mD (middle) and T * (bottom row).

FIG. 8 .
FIG.8.Ratio of the diffusion coefficient to the equilibrium value for the same energy density.The wide transparent lines show the parametric estimate given by (40) using the extracted values for mD and T * depicted in Fig.7.The gray dotted line corresponds to the parametric estimate(47).

Appendix A: Numerical framework 1 .
Details of the discretization procedure (p min )

2 DT− p min log 1 − e − p min T π 2 ,
a. Debye mass mDThe p min dependence in m D given by(11) is obtained by inserting an IR cutoff into the integral as follows m Li 2 e − p min T where the distribution is given by the Bose-Einstein distribution.The dilogarithm is defined as Li 2 (z) =

3 TTλT 2 + 2π 2 k 2  2 T 2 + 12λ 2 T 2 p
Nck 2 e k/T − 1 2 C H λ 2 e k/T − 4k 2   λT T Li 2 e − p min T − p min log 1 − e − p min T Li 2 e − p min T − p min log 1 − e − p min T Li 2 e − p min T − p min log 1 − e − p min T Li 2 e − p min T − p min log 1 − e − p min T Li 2 e − p min T − λp min T log 1 − e − p min T T Li 2 e − p min T − p min log 1 − e − p min T π min log 1 − e − p min T − T Li 2 e − p min T