Heavy Quark Momentum Broadening in a Non-Abelian Plasma away from Thermal Equilibrium

We perform classical-statistical real-time lattice simulations to compute real-time spectral functions and momentum broadening of quarks in the presence of strongly populated non-Abelian gauge fields. Based on a novel methodology to extract the momentum broadening for relativistic quarks, we find that the momentum distribution of quarks exhibit interesting non-perturbative features as a function of time due to correlated momentum kicks it receives from the medium, eventually going over to a diffusive regime. We extract the momentum diffusion coefficient for a mass range describing charm and bottom quarks and find sizeable discrepancies from the heavy quark limit.

Introduction Heavy quarks are excellent probes of the Quark Gluon Plasma (QGP) formed in collisions (HIC) of heavy nuclei, performed in large-scale collider experiments at RHIC and LHC [1].Heavy charm (c) and bottom (b) quarks are particularly interesting since these are almost exclusively formed in the very early stages of ultra-relativistic HICs through perturbative hard scattering process of the colliding nuclei.Subsequently, these heavy quarks interact with the QGP formed in HICs, thus charting out its entire time history until hadronization.Historically, the in-medium dynamics of heavy quarks has primarily been modeled as a Brownian motion in a medium consisting of colored degrees of freedom through the Langevin equations [2,3].More recently, the effect of color interactions [4] have been implemented including color memory effects [5,6] in the Langevin dynamics.Despite the existence of successful phenomenological models, major theoretical challenges in the description of heavy quark dynamics still remain to be addressed.Firstly the QCD medium traversed by the heavy quark both in the pre-and near-equilibrium stages interacts with it non-perturbatively [7][8][9][10], and gluon (color) exchange processes cannot be described semi-classically except in the limit of infinite number of colors [11].Because of this the dynamics of heavy quarks in a colored medium can be described using Langevin equation only under specific circumstances [11].Secondly, in view of the experimental data on the elliptic flow v 2 for charm quarks [12][13][14][15] -the non-relativistic approximation often applied for the heavy quarks may not be well justified for charm quarks.Since the experimentally observed v 2 of charm quarks is consistent with that of light quarks, this may hint towards the need for a full relativistic treatment of charm quark dynamics beyond the existing studies at or close to the infinite-mass limit [16].
To understand the interactions of the heavy quarks with the QCD medium, a quantity of interest is the heavy quark diffusion coefficient κ, which quantifies the autocorrelation function of the noise term that model the random momentum kicks that it receives in the medium during its evolution [17,18].In a plasma in local ther-mal equilibrium, the fluctuation-dissipation theorem relates κ to the drag coefficient which quantifies the rate at which kinetic equilibration occurs.Since its reconstruction at the level of spectral function is difficult [19], at leading order in 1/m, where the mass m of the heavy quarks which is assumed to be much larger than all other scales, the coefficient κ can be related to the spectral function of the color-electric field correlator [20,21].In this approximation, κ has been calculated within perturbative QCD, showing a poor convergence [22] emphasizing the need for its non-perturbative calculation.Such a non-perturbative calculation in the continuum limit using lattice techniques have been performed both in the quenched case [22][23][24][25][26][27][28][29] as well as in QCD medium with dynamical fermions [30].These results on κ are significantly lower than the perturbative estimates [31,32] moving towards the upper edge of the theoretical systematic uncertainty band of weak-coupling calculations at temperatures of ∼ 1 GeV [29].Very recently 1/m 2 corrections to κ have been calculated, which receive different contributions, significant among which is the color magnetic part of the Lorentz force [33].The κ that arises out of the color magnetic correlator has non-trivial renormalization effects which are being investigated only in recent years [34][35][36][37][38].The magnetic part of κ is comparable to its electric counterpart [37], which implies a 20 (34)% correction to κ for bottom (charm) quarks respectively using lattice results of their mean-squared velocities [39].
Recent phenomenological studies also demonstrate that a significant amount of v 2 for heavy quarks already starts to build up during the pre-equilibrium stages of HIC [40,41], hence a precise determination of κ in this regime is important.At early times this medium has non-perturbatively high occupation numbers of lowmomentum gluons, hence the precise determination of κ is challenging.The κ has been estimated within transport models [1, [42][43][44] as well as in classical-statistical lattice gauge theory simulations in the infinite quark mass limit [45,46] showing a sizeable difference from perturbative kinetic theory estimates [47].It will be thus desirable to develop a more fundamental theoretical for-malism to measure κ without resorting to an expansion in powers of 1/m.We thus, for the first time report the value of the heavy quark diffusion coefficient in a nonequilibrium plasma, treating them as relativistic particles.We choose a specific non-equilibrium condition of a highly occupied non-Abelian plasma in a self-similar scaling regime at its classical fixed point and perform a time evolution of quarks inside such a medium.We monitor the momentum broadening as a function of time and measure its κ from the width of the distribution.Our formalism is very general which enables us to measure κ for a wide range of bare quark masses.
Details of the classical-statistical lattice simulations We follow previous works [48][49][50][51] and first generate gauge field configurations using a classical-statistical evolution of the color fields using the classical Yang-Mills equations of motion defined using the (Minkowski) Wilson gauge action on a lattice of spatial volume N 3 and spacing a. Color electric fields E a i and gauge links U i are timeevolved up to a reference time of Qt = 500, 1000, 1500 using a leap-frog integrator, where Q is the characteristic hard scale in the problem.By these times the gauge field evolution reaches a self-similar regime where the momentum distribution function f g attains a fixed point scaling with time, of the form Q with a scaling function f s [49,50].We perform our simulations with two different initial conditions for the SU (2) gauge group for performing statistical averages, and if not stated otherwise we consider lattices of size N = 256 with lattice spacing a s = 0.5/Q and a time step a t = 0.025 Q in our work.
Over the course of the self-similar evolution, the nonequilibrium plasma establishes a separation of scales [52,53], which is analogous to the separation of scales in a high-temperature equilibrium plasma.The hardest scale is set by the highest momenta Λ(t) ∼ Q(Qt) 1/7 , the electric screening (Debye) scale is m D (t) ∼ Q(Qt) −1/7 and the magnetic screening scale is simply σ(t) ∼ Q(Qt) −3/10 .While at initial times all these scales are ∼ Q, they start to separate dynamically culminating in the scaling regime where σ(t) < m D (t) ≪ Λ(t), such that at Qt = 1500 (1750), where we typically start (end) our measurements, the scales are Λ = 2.1 Q (2.14 Q), m D = 0.21 Q, (0.206 Q) and √ σ = 0.03 Q, (0.029 Q) respectively [54].
We then study the dynamics of relativistic quarks in this scaling regime, as described in detail in Appendix A. Starting from a reference time Qt, we initialize the fermionic wavefunction as ϕ u/v λ,s (t, x) = u/v λ,s (P)e ±iP.x . 1 , representing a free relativistic quark/anti-quark with a fixed momentum P, spin index s = 1, 2 and color in-1 We globally employ temporal axial (A 0 = 0) gauge, and fix the residual gauge freedom by further implementing the Coulomb dex λ = 1, • • • , N c at the initial time Qt.Subsequently, the quark wave functions ϕ u,v λ,s (t, x) are then evolved quantum-mechanically as a function of time, t ′ , in the background of classical color gauge fields using the O(a 2 )improved Wilson-Dirac Hamiltonian on the lattice, where D i (i = 1, 2, 3) represents the spatial part of the Wilson-Dirac operator in presence of the gauge fields.By performing a spatial Fourier transform of the timeevolved quark wavefunctions at time t ′ > t denoted as φu,v λ,s,p (t ′ , p), the quark spectral function can be calculated by projection on the free Dirac spinors corresponding to (anti) particles (v λ,s (p)) u λ,s (p) as [55] which can then be decomposed into scalar (ρ S ), pseudoscalar (ρ P ), vector (ρ µ V ), axial-vector (ρ µ A ) and tensor (ρ µν T ) components using the Clifford basis.
Similarly, the momentum distribution dN/d 3 q at time t ′ = t + ∆t is calculated by projecting the evolved quark wave functions φu λ,s (t ′ , x) (measured in the Coulomb gauge) onto the non-interacting eigenstates where the sums over λ, λ ′ , s, s ′ correspond to averaging over initial and final spin and color states.Since for initial momentum P ≡ (0, 0, 0), the momentum distribution dN/d 3 q measures the momentum acquired by the quark through interactions with the medium, the second moment of the distribution can then used to calculate the diffusion coefficient κ as discussed further in the penultimate section (also see Appendix B).
Quark spectral functions -Effective masses and lifetimes We first investigate the behavior of the spectral function of a quark inside the medium, which contains information about its dynamics in an interacting quantum field theory.We will be considering a wide range of bare quark masses m bare = m in units of the initial hard scale Q and study the corresponding spectral functions inside the highly-occupied non-Abelian plasma described in the previous section.We anticipate that in the nonequilibrium Glasma created in the initial stages of highenergy heavy-ion collisions the characteristic hard momentum scale Λ ≳ 1 GeV, such that the values m/Q = 1-2 and m/Q = 4-8 represent quarks with masses close to charm and bottom quarks respectively.gauge condition (∂ i A i = 0) at the time of initialization and measurement of gauge dependent observables [49,55].
We recall that in a previous study [55], it was shown that the behavior of the spectral function for massless quarks can be well described by a modification of the hard-thermal loop (HTL) result [56] introducing an effective damping rate γ, such that for vanishing spatial momentum |p| = 0, the vector component ρ V of the spectral function can be parametrized as [55], Re ρ 0 V (t) = e −γ.t cos(m eff .t).We find that this ansatz also describes the spectral functions of heavy flavor quarks, when the medium induced quark mass m eff and damping rate γ are considered to depend on the bare quark mass m bare .Only for intermediate quark masses, (m eff − m bare ) ∼ m bare , we find that the spectral function is not well described by this rather simple form (see Appendix C).By performing a fit (c.f.Appendix C for details) of the above ansatz to the real-time spectral functions in a time interval Q∆t = 500(100) for bare quark masses m bare /Q < 0.01 (≥ 0.6), we obtain the medium induced masses and the widths shown in Fig. 1.The medium induced mass m eff − m bare is ∼ 0.08 Q for a range of bare masses for which m bare < (m eff − m bare ) beyond which it decreases rapidly falling to less than a tenth of that value for m bare = 2.1 Q ∼ Λ(t).This suggests that for the typical mass range near the bottom quarks the mediuminduced effects are negligibly tiny whereas for the charm quarks, it can still be significant but less than that of the light quarks.As a naive comparison, for 2 + 1 flavor QCD in thermal equilibrium, the thermal mass for light quarks comes out to be m th /T = 0.725( 14) at around T = 3T c ≃ 470 MeV [57], though the spectral functions in both these regimes can be very different and needs a dedicated study [55,58].From the plot of the damping factor γ as a function of the bare quark mass in Fig. 1, we conclude that the width of the effective mass peak in the spectral function is quite large ∼ 20% for the range of m bare corresponding to the light quarks, falling to less than ∼ 10% for heavier quarks.Importantly, if the width is too broad, a heavy quark in the medium will correspond to a collective excitation and cannot be simply treated as a weakly interacting quasi-particle.Since the typical spectral width of heavy quasi-particles is significantly smaller than the lighter counterparts, these can be treated as well-defined quasi-particles and their dynamics can be modelled in terms of a Langevin equation [2].Conversely, for light quarks one would require a more involved treatment of the medium-particle interactions, as their large width of the spectral function indicates a significant modification due to the medium.Interestingly, particles with bare mass of the same order as the charm quark, have intermediate widths between the very heavy and light quarks and thus require a careful treatment, beyond the infinitemass approximation.
Momentum broadening of heavy quarks We next study the scenario where a heavy quark with some fixed initial momentum propagates through the medium.It is expected that if the quark behaves like a Brownian particle and experiences uncorrelated kicks from the gluons, it will diffuse with time and its momentum distribution will show a spread about its initial momentum value which can be quantified though the heavy quark momentum diffusion coefficient κ.By analyzing the momentum spectrum dN/d 3 q of a quark evolving in the dense gluonic medium for different values of the bare quark mass m, we can then investigate to what extent the dynamics of a heavy quark can be approximated in terms of the motion of a non-relativistic Brownian particle.
We present a compact summary of our findings in Fig. 2, where we show the momentum distribution along the q x -direction (obtained by integrating over spatial momenta q y and q z ) at different times Q∆t as a function of the O(a 2 ) improved lattice momentum qx = −4/3 sin(2πn x /N ) + 1/6 sin(4πn x /N ) for quark masses m/Q = 1.2, 4.2.Starting with an initial distribution which is a δ-function at momentum q ≡ (0, 0, 0), we observe that the momentum broadening is clearly dependent on the quark mass, and larger for the lighter quark mass i.e 1.2 Q compared to 4.2 Q. Significant differences emerge in the tails of the momentum distribution, which corresponds to momentum transfer due to hard scattering with the medium constituents, while the broadening of the central peak is quite similar for the different masses.
We further quantify the momentum broadening through the second moment of the momentum distribution, which measures the effective width of the distribution and hence the variance of the momentum along q x : Our results presented in Fig. 3 show distinct features as a function of time; first a transient oscillatory behavior due to the initial drift in the correlated gluonic medium and then a subsequent monotonic increase with time.In the same Fig. 3, the inset shows a comparative study of the momentum broadening at very early times, for a range of quark masses and also in the static i.e. m ≫ Λ(t) limit, shown as a gray line.The ⟨q 2 ⟩ calculated in the static limit as obtained from the color-electric field two-point correlator shows a very rapid early-time growth ∝ (∆t) 2 , before settling to an approximately linear rise ∝ ∆t at later times.By following the arguments of [47], the initial quadratic rise occurs due to the fact that the high energy (∼ Λ(t)) classical gauge field modes act coherently for very initial times ∆t ≲ Λ(t) −1 .Even though a classical particle would feel the this coherent force, the quarks being inherently quantum mechanical excitations, do not feature this coherent early-time behavior and instead only exhibit the approximately linear rise, at later times where they can effectively be treated as classical colored particles.By comparing the late time momentum broadening ⟨q 2 ⟩ for different masses m in Fig. 3, we also observe that for m/Q ≥ 6, the results converge towards the infinite mass limit given by the color electric field correlator.However, in the mass range of charm and bottom quarks there are still significant differences which highlights the importance of inclusion of relativistic effects, especially for the lighter charm quarks.
Since at Qt = 1500 the evolution of the nonequilibrium plasma is slow compared to the time scale of the measurement, a fit to the momentum broadening data at late times Q∆t > 250 for quark masses m/Q ≥ 1.2 shows a linear growth as a function of time.Such a linear growth then signifies the onset of a diffusive regime, and we can therefore extract the momentum diffusion coefficient κ from the slope of the linear fit to ⟨q 2 ⟩ at late times.For m/Q = 1.2 (2.1), we FIG.3: Evolution of the second moment ⟨q 2 ⟩ of the momentum distribution dN/d 3 q as a function of the propagation time ∆t in the medium.Colored curves indicate the results for different bare quarks masses, while the gray curve shows the result obtained in the infinite mass limit from the color-electric field two-point correlation function, as in [45].Differences between the relativistic quantum mechanical treatment and the effective description at early times, are highlighted in the inset.
extract a κ = 37.13(5) × 10 −5 Q 3 (16.54(4)× 10 −5 Q 3 ) which shows a significant departure when compared to the earlier estimate of κ = 5.9 × 10 −5 Q 3 obtained in the infinite mass limit [45].On the other hand κ = 8.64(2) × 10 −5 Q 3 (7.24(1)× 10 −5 Q 3 ) for the bare quark mass m/Q = 4.2 (6.0) similar to the bottom quark is much closer to the static estimate.The results for the variation of the momentum diffusion coefficient as a function of the inverse of bare quark mass in units of Q, at different times of initiation of quark motion Qt = 500-1500, are summarized in Fig. 4. The lines in the bottom left corner of the plot are the values of κ obtained in the infinitely massive limit and agree well with the calculation in [45].To quantify whether this deviation in the range of masses between charm and bottom quarks can be understood in terms of 1/m 2 correction to the static limit estimate, κ(∞), we perform a fit of the values of 2 where this ansatz is appropriate, and then extrapolate to lower values of m.If charm, bottom quarks also follow this trend then the κ(m ≤ 4.2) will lie on the extrapolated graphs, which is not observed in our data, re-emphasizing the need for a relativistic treatment of heavy-flavor dynamics.
Summary & outlook Based on microscopic ab-initio lattice simulations we show for the first time the detailed dynamics of a relativistic quark in a highly occupied non-Abelian plasma.We propose a new method to measure the momentum distribution of a quark traversing in the plasma, and based on this method report a significant FIG.4: The momentum diffusion coefficient κ as a function of inverse of bare quark mass, at different times Qt = 500, 1000, 1500.Solid lines at the bottom left corner correspond to the values obtained in the infinite-mass limit using the color-electric correlator.
Shaded bands represent 1/m 2 correction to the infinite-mass limit.
excess of the momentum diffusion of quarks in the mass range from charm to bottom, than previously calculated in the static limit, thus highlighting the need for treating their dynamics through a relativistic Dirac equation.
Beyond the corrections to the infinite mass limit, we also observe significant correlations through the interaction of quarks with the medium, which manifest themselves in oscillations of the ⟨q 2 ⟩(∆t), indicating that a Markovian description is only applicable on relatively large time scales.Since a highly occupied gluonic plasma considered in our analysis is believed to be produced in the initial stages of heavy-ion collisions, our study re-emphasizes the fact that a significant momentum broadening and hence a build-up of elliptic flow of charm and bottom quarks may already occur in the initial stages of heavyion collisions [40,41,46,47].However, to unambiguously confirm this, our calculation needs to be extended to an expanding space-time geometry to properly capture the dynamics of the correlated gluon fields in the Glasma.Since our formalism for investigating the dynamics of quark probes is quite general, one can also study the dynamics of light-quark jets through a highly-occupied non-Abelian plasma, which we would like to address in a future study.
Acknowledgements.This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 'Stronginteraction matter under extreme conditions'-project number 315477589 -TRR 211.Sa.S. gratefully acknowledges partial support from the Department of Science and Technology, Govt. of India through a Ramanujan fellowship when a major portion of this work was performed and from the Institute of Mathematical Sciences.The authors acknowledge computing time provided by the Paderborn Center for Parallel Computing (PC2).
When investigating the dynamics of heavy and light flavor quarks, we start with a single quark in a fixed momentum mode labelled by momentum P, spin polarization s and color λ and let it evolve in the background of gauge fields in the self-similar regime up to some time t ′ .By following the methodology of [59][60][61] the associated fermionic field operator Ψ(t ′ , x) can be defined in terms of the time-independent operators b † λ,s (t = 0, p) and d † λ,s (t = 0, p) and time dependent wavefunctions ϕ u,v λ,s as, where for our initial conditions the creation/annihilation operators satisfy Now in order to extract the momentum distribution of quarks inside the medium, we make use of the fact that at any time t, the fermion field operator Ψ(t, x) can be decomposed in terms of time-evolved creation and annihilation operators, labeled by the momentum p, the color index λ and spin index s , as By inverting the above expression using the orthonormality of Dirac spinors, we can obtain then the timedependent creation operators b λ ′ ,s ′ (t, q) in terms of Fourier transform of the time-evolved fermion fields as, Now, performing Fourier transform of the quark fields and substituting them back in Eq. 7, we get the timeevolved annihilation operator b λ,s (t, q) entirely in terms of time-dependent fermionic wavefunctions ϕ u,v λ,s (t ′ , x).
The momentum mode occupancy distribution of fermions in a particular momentum state labelled by q is defined as the expectation value of the number density operator at that particular value of the momentum, Now, instead of calculating this observable in terms of b λ ′ ,s ′ (t, q) and b † λ ′ ,s ′ (t, q) we can use Eq. 7 and substitute b's in terms of the Fourier modes u λ ′ ,s ′ (q) and φu,v λ,s (t ′ , q), the later quantity representing the fermion wavefunctions in the Fourier space, given as φu,v λ,s (t Thus the momentum distribution can now be rewritten as, where we have used the initial conditions given in Eq. 5 to sum over various spin polarizations and color quantum numbers.
Appendix B: Calculating the second moment of momentum distribution function through fermionic wavefunctions and electric field correlators Once we obtain the momentum occupancy distribution at time t from fermionic wavefunctions following the procedure described in the earlier subsection, for the whole set of lattice momenta (with n x , n y , n z lattice indices), we calculate the second moment (as a function of ∆t) of the distribution and define it as a measure of momentum broadening, i.e., where, qi is the lattice momenta defined as, where n i = 0, ..., N − 1 is the lattice momentum mode index and 0, 1, .., C k−1 are coefficients used for O(a k ) improvement of the Wilson-Dirac operator.Specifically, for the NLO i.e.O(a 2 ) improvement, the coefficients are C 1 = 4/3 and C 2 = −1/6.Now, from the definition of the second moment of the momentum distribution, we can extract ⟨q 2 i ⟩/Q 2 as defined in Eq. 10.
Now, to compare the results to those obtained in the infinite-mass limit, we need to calculate ⟨q 2 i ⟩/Q 2 using color-electric field correlators [45] which is given as, where the repeated index a implies a sum over number of color generators a = 1, ..., N 2 c − 1.On the lattice, the second moment of the momentum distribution can be defined as, where, Êa i (t n ) = ga 2 s E a a (t) denote the components of lattice electric field defined at time t n and a t is the lattice spacing in the temporal direction.
In the diffusive regime, when the time variation of is linear, we can calculate momentum diffusion coefficient κ from the second moment of momentum distribution from the following relation, The factor of 1/3 in the expression of κ/Q 3 comes from the fact that we have isotropy in the space directions and the momentum diffusion coefficient measures the amount of broadening along any of the spatial dimensions.We have calculated the slope using Q∆t = 350 which is long enough for the diffusive and hence a linear dependence in time to set in.Values of κ/Q 3 thus obtained for various quark masses and from the color-electric field correlator i.e. in the infinite massive limit are shown in table I.
Quark mass, m/Q κ/Q 3 (×10   Re ρ 0 V (Q∆t) = exp(−γ • ∆t) cos(m eff • ∆t).In Fig. 5 we plot the spectral function measured in our simulations, the fit obtained using our ansatz and the residual difference that is obtained subtracting the fit function from the data.It is evident from the figure that our ansatz for the vector spectral function corresponding to light quark masses m bare /Q < 0.01 as well as heavy quark masses m bare /Q ≥ 0.6 provides an excellent fit to the data, yielding a quite small residual.However, for intermediate quark masses for which m bare ∼ (m eff − m bare ) the residual, i.e. the difference between the fit-function and data, for the spectral function is noticeably larger and we conclude that in this mass range the spectral function is not well described by the aforementioned ansatz.
FIG. 1: Dependence of the medium-induced mass m eff − m bare (triangles) and the damping factor γ (circles) on the bare quark mass m bare at Qt = 1500.Vertical lines indicate magnetic ( √ σ) and electric (m D ) screening scales, as well as the hard scale Λ.

FIG. 2 :
FIG.2: Evolution of the momentum distribution dN/dq x as a function of time Q∆t for quark masses m/Q = 1.2 and m/Q = 4.2.Starting from vanishing momentum dN/dq x = δ(q x ) at initial time Qt = 1500, the quarks acquire momentum through interaction with the non-Abelian plasma.

Appendix C : FIG. 5 :
FIG. 5: HTL-inspired fit function (green lines) to the Re ρ 0 V component of the spectral function, shown along with the lattice data (blue points) for a wide range of bare quark mass.In each plot the light-blue points represent the residual i.e. the difference between the fit function and the data.

TABLE I :
The momentum diffusion coefficient as a function of quark mass.