Heavy quark drag and diffusion coefficients in the pre-hydrodynamic QCD plasma

Kinetic and chemical equilibrations play important roles in the formation of the quark-gluon plasma (QGP) in relativistic heavy-ion collisions (HICs). These processes further influence the production of hard and electromagnetic probes in HICs, in particular, the thermalization of heavy quarks, which are produced at an extremely early time before the formation of the QGP. We calculate the drag and diffusion coefficients of heavy quarks in the pre-hydrodynamic quantum chromodynamic (QCD) plasma with the state-of-the-art QCD effective kinetic theory (EKT) solver. We present the time, momentum, and angular dependencies of these coefficients for gluon and quark contributions separately, showing the effects of isotropization and chemical equilibration from the QCD plasma. We also provide a simple formula to estimate the heavy quark drag and diffusion coefficients, as well as its energy loss, within the pre-hydrodynamic plasma at different coupling strengths based on the attractor theory. We then discuss the validity of these estimations with leading-order calculations and leading-logarithmic rescaling factors.


Introduction
Thermalization is omnipresent and there are two main categories of thermalization systems that are being widely studied due to the simplification of degrees of freedom at certain scales.One is the thermalization of manybody closed quantum systems, where the microscopic dynamics of single particles can be coarse-grained and the emergent behavior is more interesting.The other one is the thermalization of open quantum systems, where the objects have distinguished degrees of freedom or scales from the background medium environment, and tracing out the environment leaves a simple dynamical description of the system in the medium.Relativistic heavy-ion collisions (HICs) are such experiments that both categories are present for us to understand the fundamental strong interaction.
On the other hand, heavy quarks have distinguished mass scales from light partons in the QGP, and are produced as open quantum systems due to their large mass thresholds.Heavy quark thermalization is contributed by energy loss and diffusion.Most of the heavy quarks are produced within the momentum range p ≲ m HQ , where the radiative energy loss can be neglected and the collisional energy loss dominates [33].They are produced at a time scale τ ≃ 1 mHQ before the hydrodynamization of the QGP at τ h ≃ 4πη T s , and relax at a much later time τ R ≃ mHQ T τ h .Thus, most of the heavy quark thermalization simulations [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] are focusing on the hydrodynamic stage when the QGP is nearly thermalized.There are some efforts in addressing the heavy quark thermalization in the pre-hydrodynamic glasma or QGP [50][51][52] taking care of the anisotropy or chemical effects.The EKT allows us to calculate the heavy quark diffusion [53] as well as jet momentum broadening [54] dynamically during the pre-hydrodynamic stage from the first principle.Furthermore, the recent developments of the EKT to a full QCD level will complete this picture.
In this article, we will show the first-principle calculations of heavy quark drag and diffusion coefficients in the QCD plasma, from the state-of-the-art QCD effective kinetic theory (QCD EKT) solver [13] including both gluon and quark dynamics at fixed coupling strengths.Since the pre-hydrodynamic QCD plasma undergoes a rapid drop in temperature and a transition from weakly coupled to strongly coupled, the interpolation from the weakly coupled plasma to the strongly coupled plasma as a realistic medium profile would be favored.Although the EKT breaks down at strong coupling, the universal attractor can perform such an interpolation regardless of the coupling strength.It also avoids running complex EKT simulations at various coupling strengths by provid-ing the rescaling of key medium profiles like the time τ and the temperature T which are essential ingredients for heavy quark thermalization.Augmented with a proper rescaling factor for the heavy quark diffusion, one would arrive at a realistic heavy quark energy loss estimation in the pre-hydrodynamic stage.In this paper, we perform the first step towards this goal by justifying the validity of the rescaling at weak coupling and demonstrating the breakdown of the rescaling at strong coupling with leading-order (LO) perturbative QCD (pQCD) calculations and leading-logarithmic (LL) rescaling factors.

Pre-hydrodynamic QCD plasma & attractor
The QCD plasma out-of-equilibrium before the formation of the hydrodynamic state can be described by the QCD EKT, with a Bjorken expansion at the early stage of HICs.Within the EKT, the evolution of gluon and light quark/anti-quarks as constitutes of the QCD plasma is formulated as a set of coupled Boltzmann equations [55] with a = g, q, q and flavor number The expansion term with prefactor 1/τ anisotropizes the plasma at an early time, while the collision terms C a [f ] take over the evolution at a later time and drive the plasma to reach hydrodynamic equilibrium, in both kinetic and chemical sense.Details of the QCD EKT and its numerical implementations can be found in our previous paper [13].
The early-time expansion and the later-on hydrodynamization are independent of the microscopic interactions in the kinetic theory, resulting in a universal attractor solution.This solution connects the energy density of the plasma at any time e(τ ) to its initial value e 0 in a simple and universal way [27,28] The function E(ω) is called the energy attractor in terms of the universal and dimensionless time scale ω = τ T s 4πη .The effective temperature can be evaluated by Landau matching T = 30e(τ ) π 2 ν eff 1 4 .The degeneracy factor ν eff = ν g + 7  4 ν q N f = 47.5 and C ∞ = 0.87 are both constants.The shear viscosity over entropy density ratio η/s directly reflects how strong the interaction is and how quick the equilibration can be.Indeed, at any ω, there is a universal energy attractor E(ω) that characterizes the degree of thermalization, valued from E(ω → 0) = 0 to E(ω → ∞) = 1.One can evaluate the corresponding  (τ 0 e 0 ) This means that a more strongly coupled plasma with a smaller η/s requires a shorter time to reach a certain degree of thermalization, while a more weakly coupled plasma with a larger η/s requires a longer time.With a fixed initial energy density, the shorter thermalization time in strongly coupled plasma also results in higher initial temperatures for the following hydrodynamics.As a consequence, the universality of the attractor gives the approximate relation in the pre-hydrodynamic stage for varying coupling strengths The corresponding η/s can be extracted from fitting to hydrodynamic constitutive relation as it was done in [13].Fig. 1 shows comparisons of the rescaled time τ (η/s) −4/3 and temperature T (η/s) 1/3 at various 'tHooft coupling λ = g 2 N c from the QCD EKT simulations.One observes the rescaling relation of Eq. ( 4), even up to large coupling λ = 60 where the η/s ≃ 0.084 is already close to the lower bound from holography η/s = 1/4π [56], although slight deviation is observed since the kinetic theory should break down at such a large coupling.
There is only one scale in a conformal theory and one can fix it by matching experimental data at the end of the QGP evolution.In equilibrium, E(ω ≫ 1) = 1 and sT = e + p at zero net-baryon density.The equation of state e = 3p always holds in a conformal theory.One has the entropy density in equilibrium (τ s) eq = 4 3 (τ 4 3 e) eq (τ eq π 2 ν eff 30 which is related to the charged particle multiplicity via dN ch /dη = N ch S (τ s) eq S ⊥ with S/N ch = 8.36 [57] and S ⊥ the transverse area of the collision.As a consequence, one can constrain the initial condition τ 0 e 0 with the following relation [27] dN ch dη = 4 3 The average η/s represents a measure of the thermalization speed for both the pre-hydrodynamic stage and the following hydrodynamic stage until freeze-out.It is, however, dominated by the hydrodynamic stage due to a much shorter lifetime of the pre-hydrodynamic period compared to the hydrodynamic period.A strongly coupled fluid dynamic features a small η/s close to a holographic bound η/s = 1/4π [56].Matching the LHC 5.02 TeV Pb+Pb collision data [58], dN ch /dη = 1942 and S ⊥ = 138fm 2 to Eq. ( 6), one arrives at τ 0 e 0 = 1.961GeV 3 in this specific collision system.One has further the rescaling formula to evaluate the initial energy density in other circumstances The initial energy is mainly deposited by an overoccupied, anisotropic, and gluon-saturated state; and described by the CGC-inspired distribution [11,59] with the momentum decomposition in transverse and longitudinal directions ⃗ p = (⃗ p ⊥ , p ∥ ).The anisotropic parameter is typically chosen as ξ = 10.A typical value for the 'tHooft coupling λ 0 = g 2 0 N c chosen in simulating weakly coupled gauge field in the CGC effective theory is λ 0 = 10, which is reasonable as well at the initial time for

Equilibrium limit
Characteristic scales ω ~1st-order hydro p L /e Kinetic p L /e Kinetic e q /(2N f )/e g FIG.2: Typical characteristic scales of isotropization pL/e (red) and chemical equilibration eq/eg (blue) in terms of universal time ω, compared to the hydrodynamic limit (black) and their equilibrium limits (dashed).These curves and related discussions can also be found in our previous article on QGP thermalization [13].
our QCD kinetic simulation when the system has a high temperature.The general 'tHooft coupling λ = g 2 N c entering into the QCD kinetic simulation is controlling the thermalization speed, which can be further reflected in the macroscopic coefficient η/s.We keep the weak coupling λ = 10 for the QCD plasma as the default throughout the QCD kinetic simulation and perform rescaling to evaluate strongly coupled plasma, where the validity of both the kinetic theory and the perturbation theory breaks down.General rescaling can be achieved due to the universality of the attractor solutions, from the basic principles of energy conservation and conformality, regardless of the coupling or modeling.To extend the rescaling from the QCD plasma to the heavy quark, we will discuss the validity of rescaling for the transport coefficients in a standalone section later.The initial time for the QCD kinetic evolution is approximately the inverse of the saturation scale for the gauge fields τ 0 ≃ 1/Q s before the formation of the quasi-particles in the kinetic theory.Without losing generality, by choosing τ 0 = Q −1 s we can calculate the initial energy density and one gets τ 0 e 0 = 0.5858Q 3 s = 1.961GeV 3 .Now one can estimate that Q s = 1.496GeV and τ 0 = Q −1 s = 0.134fm, smaller than the typical hydrodynamization time τ h ≃ 0.2−0.6fmbut at the same order of magnitude.
By solving the QCD EKT, one gets the time evolution of the distributions f g,q,q (⃗ p, τ ).Certain quantities can characterize the equilibration of the QCD plasma, such as the energy-momentum tensor The longitudinal pressure over energy density ratio p L /e = T zz /T τ τ characterizes the isotropization of the plasma with an equilibrium limit 1/3.The quark over gluon energy density ratio e q /e g characterizes the chem-ical equilibration of the plasma with an equilibrium limit (7ν q N f )/(4ν g ).We show these characteristic scales in Fig. 2 in terms of the universal time scale ω for the QCD plasma at constant coupling λ = 10.The anisotropy of the plasma approaches the hydrodynamic limit at around ω ≃ 1 − 2 while its equilibrium limit has to be reached after a much longer time.The chemical equilibration roughly finishes later than ω ≃ 2 − 3 where the quark over gluon density ratio tends to be a plateau.These non-equilibrium parton distributions will deviate the heavy quark transport coefficients from thermal cases, from the initial time τ 0 ≃ 1/Q s to the hydrodynamization time τ h ( ω ≃ 1 − 2), before the formation of the hydrodynamic plasma.This opens an opportunity to extend heavy quark simulations to the pre-hydrodynamic stage of HICs.

Heavy quark thermalization
The heavy quark thermalization with soft collisions from the background QCD plasma can be described by a stochastic differential equation (Langevin) in phase space (⃗ x, ⃗ p) with a Wiener process dW j ∼ N (0, dτ ) correlated as ⟨dW i dW j ⟩ = δ ij dτ .Applying Ito's lemma to Eq. ( 10) up to order O(dτ ), the Kolmogorov equation (Fokker-Planck) reads There are two evolving chemical contributions for the drag coefficients A i (⃗ p, τ ) and diffusion coefficients B ij (⃗ p, τ ) = 1 2 σ ik (⃗ p, τ )σ jk (⃗ p, τ ), from gluon collisions gQ → gQ and from quark/antiquark collisions qQ → qQ (including the factor 2N f in the quark sector) The drag and diffusion coefficients for heavy quark Q collided by a parton a = g, q, q in the QCD plasma can be calculated as [34] A Due to the small occupation of heavy quarks f Q (⃗ p ′ , τ ) ≪ 1, the Fermi-blocking factor for the heavy quark can be neglected.
We calculate the amplitude squares |M aQ→aQ | 2 of heavy quark scattering in LO pQCD, with a dynamical and isotropic screening mass fitting to the Hard Thermal Loop (HTL) calculation, the same as it is implemented in the QCD EKT simulation for QCD plasma [9].For details, see our previous work [13].
To calculate the drag and diffusion coefficients, we assume a charm quark mass m HQ = 1.5 GeV, and the coupling in the collisional amplitude squares the same as the background QCD plasma α s = g 2 4π = λ 4πNc with λ = 10 if not claimed otherwise.
Heavy quark drag and diffusion coefficients Due to the rotation symmetry in the transverse plane, without loss of generality, we define the momentum direction of the heavy quark in the transverse plane as the x-axis.That is ⃗ p = (p, cos(θ), ϕ) = (p, 0, 0) in cylindrical coordinate.The symmetry of the integration makes the vector ⃗ A simply along ⃗ p in transverse plane, giving trivial values of A y , A z , as well as off-diagonal terms in the B ij matrix, even if the plasma is anisotropic.
We plot the time evolution of coefficients A x , B xx , B yy , B zz in Fig. 3 including their gluon and quark compo- nents.The time evolution of the coefficients in terms of ω roughly features a pow-law behavior.The lower panel of Fig. 3 shows the isotropization, that B yy is deviated from B zz at an early time and is approaching B zz at a late time.The increasing trends in quark contribution show in both panels, that the drag and diffusion coefficients contributed by quarks become comparable to gluon at a late time.However, due to the Bose-enhancement and Fermi-blocking factors from quantum statistics, the quark contribution in the coefficients is not as significant as it is in the energy density of the QCD plasma as we see in Fig. 2, where e q ≃ 2e g at the later time.
The momentum dependencies of the coefficients are shown in Fig. 4 for various times, where we have normalized everything by either A x (p = m HQ ) or B xx (p = m HQ ) so that these two coefficients are fixed at the point (1, 1).One finds a linear dependence of A x (p)/A x (m HQ ) for p/m HQ ≲ 2. Isotropization is shown in the lower panel that B yy gradually approaches B zz at a later time.One also finds B yy ≃ B xx for p ≪ m HQ for all time.The increasing trend of the quark contribution is also presented.Now we release our constraints for the heavy quark momentum direction in the transverse plane.Still, look at the typical momentum p = m HQ , but in the x-z plane so that ⃗ p = (p, cos(θ), ϕ) = (p, cos(θ), 0) in cylindrical coordinate.We therefore present the coefficients as a function of cos(θ) = p z /p in Fig. 5. Now the breaking symmetries in the integrals makes the vector ⃗ A not necessarily along ⃗ p, resulting in many more nontrivial coefficients A x,z , B xx,xz,yy,zz , but they have vanishing points at certain angles.For example, at cos(θ) = 0, the heavy quark momentum is in the transverse plane and A z vanishes; while at cos(θ) = ±1, the heavy quark momentum is in the longitudinal plane, and A x vanishes.The isotropization can be clearly seen, for example, the ratio of A z (cos(θ) = 1)/A x (cos(θ) = 0) and B zz (cos(θ) = 1)/B xx (cos(θ) = 0) only goes to 1 when the time is large and the medium becomes isotropic.The off-diagonal diffusion coefficient B xz features a cos(θ) dependence as ⟩⟩ which vanishes at cos(θ) = 0, ±1 for all time and peaks at cos(θ) = 1/ √ 2 when approaches the equilibrium.

Rescaling of transport coefficients
With the kinetic theory simulated weakly coupled plasma and the universality of attractor valid even at strongly coupled regime, we may estimate the certain physical processes in the pre-hydrodynamic QCD plasma with different coupling strengths even at strong couplings.Indeed, for any time convolution of a physical  16) from [33].The values of the parameters for different velocities are taken from [33] as well.
Strongly coupled plasmas featured by different coupling strengths also result in different values of the transport coefficients for heavy quarks.One may generically consider the rescaling of coefficients in plasma from the original coupling strength λ o to coupling strength λ, and from the original temperature T o to temperature T as Phenomenological studies of heavy flavor energy loss suggest a large overall rescaling factor K ≃ 5 compared to the pQCD calculation [60] at α s ≃ 0.3 (comparable to λ = 4πα s N c ≃ 10), empirically representing the nonperturbative effect.This empirical value of the rescaling factor is however dominated by the contributions from the hydrodynamic period.The pre-hydrodynamic rescaling at higher temperatures would favor a smaller factor K ≤ 5.Although the rescaling for the conformal plasma with attractor theory can be safely rescaled to the strongly coupled regime, the rescalings for the drag and diffusion coefficients are theoretically non-trivial.For instance, a rescaling from LO pQCD calculation would suggest an LL factor bo ln(1/λo)+co , but it can only be safely restricted to the weakly-coupled regime.
More specifically, LL calculations for a heavy quark in a thermal background would suggest rescaling factors for the drag and diffusion coefficients as [33] with ≃ −0.647 [33].These LL factors break down at large couplings with a negative coefficient from the logarithmic term and deviate quite a bit from numerical results [61].The appearance of the negativity when increasing the coupling strength is, however, delayed by a larger heavy quark ve- Thermal QCD plasma at T=0 FIG.7: Fitting the K □ factors (dashed) in the square bracket term regularized as Eq. ( 17) with numerical LO pQCD calculations (solid) in a thermal QCD plasma at fixed T and various velocities.The fitting restricts itself at weak couplings ranging from λ = 0.01 to λ = 10.The LL K □ factors without regulators are also presented (dotted).
locity.As a consequence, a positive drag coefficient A i for λ o = 10 appears when v ≳ 0.60 from the LL factor, and a monotonously increasing K A LL factor for λ < λ 0 = 10 occurs when v ≳ 0.85.Life is much easier for a weaker coupling λ o = 5, such that the positivity appears when v ≳ 0.10 and a monotonously increasing K A LL factor requires v ≳ 0.50.
Due to the anisotropy in the pre-hydrodynamic plasma, we will only discuss the corresponding logarithmic terms in the square-bracket of Eq. ( 16).In Fig. 6, we present the time-dependent rescaled coefficients A x (T (ω), λ)/K Ax LL (T (ω); T o (ω), λ; λ o ) and B xx (T (ω), λ)/K Bxx LL (T (ω); T o (ω), λ; λ o ) calculated at various couplings, in comparison to our default coupling λ 0 = 10 (and K = 1 for λ = 10) for v = 0.70, 0.80, 0.90.It shows that the LL factors give close and presumably convergent rescaling results at weaker couplings, while the rescaling results diverge quickly at stronger couplings.A next-leading-order (NLO) correction cannot amend this due to poor convergence of the perturbative expansion [61] and it clearly presents the breakdown of pertur- bation theory at strong couplings.It also appears that the LL factors rescale better at larger velocities.For example, at v = 0.90, the diffusion coefficients have almost perfect rescalings.It is conceivable that at the limit of v → 1, the heavy quark becomes a high-energy jet and the perturbation theory is well valid.However, at such a large velocity, radiational processes dominate over collisional processes.One may notice that some rescaled curves for large couplings are missing at low velocity due to negative values of the LL factors attributed by the ln(1/µ) term.
The negativity from the logarithmic term of the LL factors in Eq. ( 16) is due to an infrared cut of the parton momentum down to m D , which is not a consistent treatment at strong couplings.If one does not impose the parton momentum to be much larger than the screening mass k ≃ T ≫ m D ≃ gT or g ≪ 1, one would expect an effective regulator in the logarithmic term.Instead of the LL factors, we can employ a parameterization of the K factors with regulators ãb , ãf , bb , bf , cb , cf Fitting the numerical LO pQCD calculations with these regularized factors K □ presented as Fig. 7, we may have a better handle of the rescaling for larger couplings.The rescaling results with the regularized factor are presented for p = m HQ , v ≃ 0.71 in Fig. 8.A nice rescaling is shown even at strong couplings.
Fast thermalization of the QGP requires the plasma to be strongly coupled when close to the hydrodynamic limit, while the perturbative calculations fail at large couplings.Although the regularized LL factors from the LO pQCD suggested by Eq. ( 17) converge the large coupling rescalings, it is still an artificial treatment dropping higher-order corrections.Rescaling of heavy quark drag and diffusion coefficients at large 'tHooft couplings λ = g 2 N c ≫ 1 might be achieved by the AdS/CFT correspondence, which suggests B ij ∼ √ λT 3 [62,63].Indeed, the non-polynomial rescaling from the AdS/CFT correspondence clearly indicates a non-perturbative effect.
However, the smooth transition of the rescaling form from the weakly coupled pQCD polynomial results to the strongly coupled AdS/CFT non-polynomial results may not be possible to construct since their picture of heavy quark diffusion is very different.This aspect is out of the scope of the current study and we use the generic factors Phenomenological consequences Since we rescale from λ 0 = 10 which has (η/s) 0 = 1, we have accordingly, the energy loss and diffusion of heavy quarks at different couplings λ Focusing on the transverse plane and assuming the heavy quark initial momentum to be p 0 ≲ 2m HQ in the xdirection, one has roughly a linear momentum dependence on the drag coefficient A The infinitesimal form of Eq. ( 18) is before the plasma reaches the hydrodynamic stage.With a proper rescaling factor K □ presumably from some nonperturbative numerical methods, the above formula is useful to simply estimate the heavy quark energy loss in the pre-hydrodynamic stage in HICs.

Conclusions & Outlook
In this article, we calculate the heavy quark drag and diffusion coefficients in a weakly coupled prehydrodynamic QCD plasma from a first principle and the state-of-the-art QCD EKT solver.We present the time, momentum, and angular dependencies of the coefficients A i , B ij with all indices.With arguments from the attractor theory, we provide a simple formula to evaluate the heavy quark energy loss in pre-hydrodynamic plasma with different coupling strengths.As a first step towards this goal, we study the rescaling of transport coefficients with the LL factors.Although a trend of convergence shows at weak couplings, the rescalings augmented with LL factors diverge at strong couplings due to the failure of perturbation theory.Releasing the restriction of the infrared cut for parton momentum down to the screening mass, one may employ a regulator in the logarithmic terms of LL factors and fit the regulator to the LO pQCD calculations, which avoids the negativity from native LL factors and leads to convergent rescalings for large couplings.
One needs to keep in mind that the attractor rescaling of the QCD plasma is valid at both weak and strong couplings.The QCD EKT simulations of the time and temperature profiles of the plasma appear to satisfy this rescaling even up to a large coupling λ = 60 whose corresponding η/s ≃ 0.084 is already close to the holographic lower bound, regardless of the failure of the semiclassical kinetic theory at strong couplings.In a realistic pre-hydrodynamic stage in HICs, the rapid drop in temperature and fast thermalization of the plasma favors a quick transition from weakly coupled to strongly coupled, which may not be easy to simulate by kinetic theory or other theories.Attractor theory provides a simple way without any complex simulation to perform a rescaling in the weakly/strongly coupled transition for the QCD plasma.
However, the difficulty in calculating the heavy quark transport coefficients may come from the invalidity of different theories at different scales and coupling strengths.A generic rescaling from a weakly coupled regime to a strongly coupled regime is theoretically non-trivial.Still, it is conceivable upon incorporating non-perturbative approaches like the T-matrix calculation [64] or functional renormalization group (FRG) for calculating the transport coefficients, which we leave for future studies.

2 √ 2
an isotropic coefficient fitted to the HTL calculations.The coefficients for the boson and the fermion sectors a b , a f , b b , b f , c b , c f are velocity dependent.The zerovelocity limit for a static heavy quark (v = 0) gives the coefficients for diffusion b b = c b = ln 2 + ξ and b