Heavy quark diffusion in an overoccupied gluon plasma

We extract the heavy-quark diffusion coefficient \kappa and the resulting momentum broadeningin a far-from-equilibrium non-Abelian plasma. We find several features in the time dependence of the momentum broadening: a short initial rapid growth of, followed by linear growth with time due to Langevin-type dynamics and damped oscillations around this growth at the plasmon frequency. We show that these novel oscillations are not easily explained using perturbative techniques but result from an excess of gluons at low momenta. These oscillation are therefore a gauge invariant confirmation of the infrared enhancement we had previously observed in gauge-fixed correlation functions. We argue that the kinetic theory description of such systems becomes less reliable in the presence of this IR enhancement.


I. INTRODUCTION
Transport coefficients, such as viscosities, diffusion coefficients and conductivities contain information about microscopic properties of the medium. In the framework of QCD matter produced in ultrarelativistic heavy-ion collisions, the evaluation of such transport coefficients has been a longstanding problem. Perturbative evaluations at Leading Order (LO) have been available for a long time [1][2][3]. More recently perturbative calculations have been pushed to next-to-leading order (NLO) accuracy [4][5][6][7][8]. In equilibrium, there have been attempts to extract transport coefficients also using nonperturbative lattice QCD methods [9][10][11].
Heavy quarks are unique probes of the transport properties of the quark gluon plasma (QGP) because of their large mass compared to the other scales of the medium. Pair production and annihilation processes are negligible, and all the heavy quarks within the medium are created in the hard processes preceding the formation of the QGP. Heavy quark observables carry information about the entire history of the medium.
In conventional transport approaches to heavy-ion collisions, the effects of early-time, nonequilibrium evolution are usually ignored. Only very recently studies have addressed the importance of the nonequilibrium evolution. For heavy quark diffusion specifically, a Fokker-Planck approach to the evolution of heavy quarks in a non-equilibrium gluon plasma or "glasma" present in the early stages of the evolution was used in [12,13]. The authors find that the glasma phase can have a sizable contribution to momentum broadening and energy loss of heavy quarks. At later stages of the non-equilibrium evolution when the quasiparticle description is valid, recent studies have indicated that the pre-equilibrium effects can be important [14,15]. In [16] jet momentum broadening in the glasma was investigated. The main result is that a colored particle can accumulate sizable momentum broadening during the glasma phase ( p 2 ⊥ = 1−4GeV 2 ). One might thus expect the pre-equilibrium phase to be important also for heavy quarks.
The heavy quark momentum diffusion coefficient κ can be studied in multiple ways. In thermal equilibrium, it has been calculated with perturbative methods [1,2,[17][18][19][20] and studied with a standard lattice approach [21][22][23][24][25][26]. Another possibility is to use lattice gauge theory in the classical approximation. This technique has been applied to the heavy quark diffusion coefficient κ ∞ and the jet quenching coefficientq in thermal equilibrium systems [27][28][29]. However, one of the benefits of the classical approach is that one can also study nonperturbative systems out of equilibrium, as we will do here. Once the heavy quark diffusion coefficient is known, one can use it to understand heavy quark flow and spectra by incorporating the diffusion process in a simulation of the heavy ion collision [18,[30][31][32].
The heavy quark diffusion coefficient κ is not only important for momentum broadening of heavy quarks, but it also has applications for quarkonia. Quarkonia can be modelled using an open quantum system approach [33][34][35][36][37], and their time-evolution is governed by the Lindblad equation [38,39]. The equation of motion needs two transport coefficients as an input, one of which is the heavy-quark diffusion coefficient.
Our aim in this paper is to understand momentum broadening p 2 and the evolution of the momentum diffusion coefficient κ of heavy quarks in a far-fromequilibrium overoccupied system, with the main motivation coming from initial stages in ultrarelativistic heavy ion collisions. After the collision, occupation numbers of gluonic fields at the characteristic momentum scale Q are non-perturbatively large ∼ 1/g 2 [40,41] during initial stages in a weak-coupling thermalization picture. In this case, classical-statistical simulations are applicable and have been widely used  to understand the pre-equilibrium dynamics in the collision.
In this paper we simulate a highly occupied plasma in SU (2) Yang-Mills theory in a three dimensional fixed box in a self-similar regime [48,56,[65][66][67]. We extract the heavy-quark diffusion coefficient κ and momentum broadening p 2 in this far-from-equilibrium system using suitable definitions of these gauge-invariant observables for out-of-equilibrium dynamics. We explain how these generalized definitions can be used in quarkonium and diffusion studies and compare our results with perturbative calculations. Unlike e.g. in [68][69][70][71] we do not explicitly follow the motion of the quarks or quarkonia in the color field. Instead, we work in the infinite quark mass limit, where the quark is stationary, and measure the force acting on the quark from the chromoelectric field. More specifically, we will extract for the first time the heavy-quark diffusion coefficient in a classical nonequilibrium system. In addition to this, we find new features in the time-dependence of momentum broadening of a heavy quark. While we consider a self-similar regime, the explanation of the different features in terms of a perturbative calculation is more general. Therefore similar features could be found in gluonic plasmas with other initial conditions, e.g., in the Glasma state in the initial stage of a heavy ion collision. In particular, we observe modulations of the growing p 2 with the plasmon frequency ω pl , which we are able to attribute to an excess of gluons at low momenta as compared to perturbative predictions. To our knowledge, this is the first time that such oscillations are connected to the gluonic IR enhancement that has been observed earlier in gauge fixed correlators [72]. Our result complements the observation of gauge-invariant condensation of Ref. [73] in the same non-Abelian systems as studied here.
This paper is structured as follows. First, in Sec. II we will discuss the isotropic overoccupied gluonic system that we are studying. This system has been extensively analyzed in previous works, so we will be brief and concentrate on collecting the relevant numerical results and parametric time dependences that we will need in our subsequent analysis. We will then, in Sec. III discuss how the motion of heavy quarks in a dense gluonic system is related to the unequal time correlator of chromoelectric fields, how it is connected to quarkonium and diffusion studies, and present our numerical calculation of this correlator in the overoccupied gluonic system. In Sec. IV we will construct two microscopic models for calculating these correlators from a momentum distribution of gluonic quasiparticles in the system, and compare our numerical results to these models in Sec. V. We will then conclude with a brief discussion of our results in Sec. VI.

II. HIGHLY OCCUPIED NON-ABELIAN PLASMA
The system we are studying is described by an SU (2) pure gauge theory. Its classical evolution starts from an initial condition that is characterised by a single-particle occupation number distribution The overall properties of this system are rather well understood from several studies [48, 56-58, 65, 74, 75]. We use the same initial conditions and numerical methods as in our previous paper [72]. Thus we will only briefly summarize the physical properties and the calculational methods of this system here, referring the reader to the references for more details. The momentum scale Q controls the typical hard momentum in the initial distribution and we will write dimensionful quantities scaled with a suitable power of Q. For low momenta p Q, occupation numbers are large f (t, p) 1 and after a transient time, the system approaches a universal far-from-equilibrium attractor state, characterized by self-similar dynamics The scaling exponent β and scaling function f s (p) of this state are universal and insensitive to details of the initial conditions or to the precise value of the (weak) coupling. For the considered d = 3 spatial dimensions, one has β = −1/7. The momentum scale that dominates the energy density Λ, i.e., the momentum of hard excitations, grows with time as More precisely, we define the hard scale as the value of p for which the integrand of the energy density ∝ p 2 ω(p)f (t, p) is maximal, where the dispersion relation can be approximated by a relativistic dispersion ω(p) ≈ p 2 + m 2 . The (non-)thermal mass m that gluons obtain in a medium can be computed perturbatively (and self-consistently) as where we have also included a mass correction in the denominator. The scale separation between this mass scale and the hard (temperature for a thermal system) forms the basis for the perturbative Hard-Loop (HTL) framework. The mass is connected to the Debye mass m D and to the plasmon frequency ω pl via The plasmon scale will turn out to be essential in our study of heavy quark diffusion far from equilibrium. Eq. (4) implies that ω pl and m D decrease with time in the self-similar regime as Physically, the plasmon frequency is the lowest energy that quasiparticle excitations can have and has been computed in Ref. [72] as the peak position of the spectral function at vanishing momentum ρ(ω, p = 0). The quasi-particle peak of ρ(ω, p = 0) is of Lorentzian form and has a width γ pl ≡ γ(p = 0), which corresponds to the damping rate of plasmon excitations. Its value for Qt = 1500 and the same initial conditions as employed here has been extracted numerically from fits to ρ(ω, p = 0) in Ref. [72]. Using the HTL formalism, it can also be computed as γ HTL pl ≈ 6.64N c g 2 T * /(24π) [76]. Here T * is the effective temperature of the soft field modes, which is given by Thus, in the scaling regime, we would expect the decay rate γ pl to decrease with time like g 2 T * , i.e., as From Eqs. (3), (6) and (8) we obviously have, at late enough times, the hierarchy γ pl ω pl Λ, which is reminiscent of the hierarchy of scales g 2 T gT T in thermal equilibrium with temperature T . The extracted values in Ref. [72] at Qt = 1500 are Λ (Qt = 1500) = 2.1 Q m D (Qt = 1500) = 0.21 Q g 2 T * (Qt = 1500) = 0.03 Q γ pl (Qt = 1500) = 0.003 Q , which indeed shows the expected separation of scales. At momenta p ω pl we observed in our earlier work that the actual occupation number distribution [72] displays a feature that we refer to as an "IR enhancement," a feature also seen in earlier studies (see, e.g., [48,57]). By this term we mean that the occupation number is significantly larger than the behavior f (p) ∼ T * /p expected from perturbation theory. A gauge theory with a non-conserved number of particles is not expected to exhibit actual condensation (see, e.g., [48,75,77]), and we do not interpret this excess of gluons as an indication of condensation in the proper sense of the word. We will discuss this feature more quantitatively in Sec. IV A.
Due to the large occupation numbers f ∼ 1/g 2 1, the non-perturbative quantum problem can be accurately mapped onto a classical-statistical lattice gauge theory with lattice spacing a s and lattice size N 3 s . Its farfrom-equilibrium evolution can be then studied using computer simulations solving classical equations of motion in temporal axial gauge A 0 = 0. The equations are formulated in a gauge-invariant way using link fields U j (t, x) = exp(ig a s A j (t, x), that replace the usual gauge fields A j (t, x) in the numerics, and chromo-electric fields E j (t, x). Since we use the same initial conditions and numerical method as in our previous paper [72], we refer the reader there and to references therein for details of our numerical approach.
Data has been averaged over typically 10-15 configurations and error bars correspond to the standard error of the mean. If not stated otherwise, we use the lattice spacing Qa s = 0.5 and lattice sizes ranging from N 3 s = 128 3 to 264 3 . We also vary a s and N s to check for possible lattice artifacts (see Appendix A).

III. HEAVY QUARK DIFFUSION AND MOMENTUM BROADENING
A. Heavy quark motion in a color field We consider a heavy quark with mass M in the highly occupied non-Abelian plasma far from equilibrium described above. We take the quark mass to be the largest momentum scale in our system, i.e., larger than a typical hard momentum scale M Λ. Consequently the formation time of the heavy quark is much shorter than any other timescale in the system. We can then assume that the production mechanism of the heavy quark factorizes from its interaction with the color field degrees of freedom. Thus we do not need to specify a particular production mechanism here, but just concentrate on the subsequent interactions of the quark with the medium. The large mass of the quark also implies that the wave function is sufficiently localized so that we can use the classical equation of motion for its momentum: Hereṗ i ≡ dpi dt and the force F i is proportional to the chromo-electric field E i acting on the heavy quark. Averaging over color states of the quark and the ensemble of color field configurations gives a zero mean force acting on the quark ṗ = 0. On the other hand, the variance of the force is given by the force-force correlator [16] The electric field correlator is evaluated at the same spatial location because in the limit of large M , the velocity of the heavy quark is negligible. In the last line, we have defined the (statistical) correlation function EE (t, t ).
The trace is taken in the fundamental representation and in the last line we used temporal gauge with U 0 (t , t ) = 1 as well as a summation over repeated color indices of the adjoint representation a = 1, . . . , N 2 c − 1. The correlation function (11) is shown in Fig. 1 for the time Qt = 1500 as a function of the time difference t − t. The signal starts with a large initial oscillation (upper panel) that quickly fades away on a time scale ∼ 1/Λ. As shown in the inset, the correlator oscillates subsequently with a small frequency ∼ ω pl . The amplitude of these oscillations is roughly a factor 10 3 times smaller than the initial quick oscillations.
It is useful to consider the Fourier transform of the correlator (11) with respect to relative time ∆t = t − t. Normally we would define this in a symmetric way as for a fixed (central) timet ≡ (t+t )/2, where we exploited that the statistical correlator is an even function. For the purpose of numerically extracting the frequency space The measured momentum broadening given by (11). We observe that after a rapid initial rise the momentum broadening increases roughly linearly in p. The linear rise is attributable to the interactions and the "steps" can be understood in the SR framework below. The initial rapid rise is shown in the inset with higher resolution. It corresponds to a decoherence effect whose early and later ∆t behavior of Eqs. (20) and (21) is shown as black-dashed and blue horizontal lines, respectively. correlation function for Fig. 1, we approximate it by for a fixed lower limit t. This is a good approximation as long as the time difference |t − t | is small compared to the rate of change as a function of the central time. In a power law cascade this latter can be estimated as the lifetime of the system, and thus the approximation is a good one when t ≈ t |t − t |, i.e., ω 1/t, which is the case here.
The frequency space signal is shown in the lower panel of Fig. 1. While it has a broad peak around ω ∼ Λ, the relevant part for heavy quark diffusion is located at low frequencies (see inset). This also illustrates a practical challenge related to the measurement: to obtain the lowfrequency behavior correctly, high accuracy is required. The observed structure of the low-frequency part can be easily understood with the help of the spectral function ρ(ω, p), which will be discussed below in Sec. IV. Here we note that the finite piece at ω = 0 stems from Landau damping of longitudinally polarized gluonic fields. On the other hand, the steep rise at ω ≈ ω pl results from quasiparticle excitations, which can only contribute for frequencies ω ω pl .

B. Momentum broadening
So far, we have discussed the force-force correlation, which corresponds to ṗ i (t )ṗ i (t ) of a heavy quark traversing the non-equilibrium gluon plasma (see Eq. (11)). Integrating it, we arrive at the momentum broadening of the heavy quark after its creation at time t as This gauge-invariant physical observable is shown in Fig. 2. Its evolution shows three important features that are associated with different time scales for ∆t: (i) rapid growth at a short time scale of the order of the inverse hard scale ∆t ≈ 2π/Λ ; (ii) damped oscillations with period ∆t ≈ 2π/ω pl ; (iii) overall approximately linear growth ∼ ∆t for 1/Λ ∆t t.
Each of these properties of momentum broadening has a different physical explanation and we will elaborate on them in this work. The modulations with frequency ω pl in (ii) are a new feature that is, to our knowledge, observed in this work for the first time. We will discuss the relation of these oscillations to the quasiparticle properties of the plasma in Sec. IV B and show (see Fig. 6) that they are related to an enhancement of infrared modes in the system over the perturbative expectation. While we will study these features (i) -(iii) in detail below, we argue here that the different time scales result from the structure of EE (t, ω). With the approximation that this correlator depends weakly on the central timet = (t + t )/2 within the integration intervals, we can write With the frequency space correlator EE (t, ω) shown in Fig. 1, we can distinguish the features based on different regimes of ∆t. When ∆t 1/Λ, basically all frequencies are included in the integral. Then the broad peak of EE (t, ω) dominates the integration, which corresponds to frequencies with ω ∼ Λ and the rapid growth observed under (i). For larger time scales 1/Λ ∆t ∼ 1/ω pl , the broad peak provides a constant shift in p 2 and the evolution of the integral is instead dominated by lower frequencies like those depicted in the inset of Fig. 1 (bottom). This corresponds to the properties (ii) and (iii).
We can compute the early rapid growth of (i) analytically. For that, we remind ourselves that the fields in the correlator EE are evaluated at the same location x = x . Due to spatial translation invariance, we can hence write 1 In the second line we used a generalized fluctuationdissipation relation to express the frequency space EE correlator in terms of the equal time EE correlator and the spectral function in frequency space. In our earlier paper [72] we have verified that such a relation holds very well in the system considered here. Since at small ∆t the integral is dominated by high frequencies, we then used the spectral function of free quasi-particles and expressed the equal time correlator in terms of the distribution f (t, p) as Here 2(N 2 c −1) counts the degrees of freedom of transverse gluons and ω p is their dispersion relation.
The last line of Eq. (17) already explains our observations of the initial rise. The integral over momenta is dominated by the hard scale and hence, ω p ∼ p ∼ Λ. At early times ∆t 1/Λ, one can approximate with the energy density ε = 2(N 2 3 . This is shown as the black dashed curve in the inset in Fig. 2. At later times ∆t 1/Λ, the momentum integral involves rapid oscillations such that the approximation sin 2 (ω p ∆t/2) ≈ 1/2 can be used, leading to This is shown as the blue horizontal line in the inset in Fig. 2. Therefore, the initial fast growth stops due to decoherence after ∆t ∼ 1/Λ at a value ∼ ω 2 pl .

C. Heavy quark diffusion
We can define the heavy-quark diffusion coefficient as the time derivative of the accumulated squared momentum (15) at late times. For that purpose, let us define the function κ (t, ∆t) as follows where we again approximated the dependence on the central time in the correlator EE (t, ω). As we will see in the following, the limit ∆t → ∞ of κ (t, ∆t) gives the quantity that is commonly known as the heavy quark diffusion coefficient.
The expression (22) has also an interpretation for quarkonium evolution and decay in a non-Abelian plasma. As detailed in Ref. [33], it is related to the real part of the color-singlet self-energy (note that we have a different definition of the time arguments) where r is the distance between the heavy quark and anti-quark. It is also proportional to the decay width Γ in thermal equilibrium if ∆t is larger than any other time scale of the system. We emphasize, however, that the relation in the last line of (24) is more general and also holds for finite ∆t. It is beneficial to understand how the definition of κ (t, ∆t) in (22) is related to the usual definition of the heavy-quark diffusion coefficient κ therm in thermal equilibrium (see, e.g., [19]). Since thermal equilibrium is time translation invariant, there is no dependence on t and we can simply set t = 0. Thermal equilibrium is also time reversal invariant and therefore the correlator EE is an even function of the time difference. Then we have where the superscript "therm" is a reminder of the thermal state considered here while the subscript ∞ corresponds to the limit ∆t → ∞. We note that in the literature, this coefficient is usually referred to as κ. In the second line of Eq. (25), we used temporal gauge (U 0 = 1). Obviously, we lose exact time translation invariance once we consider far-from-equilibrium systems as in the present work. However, we can still compute a time-dependent heavy-quark diffusion coefficient κ ∞ (t) in analogy to Eq. (25) as Here the limit ∆t → ∞ is replaced by the condition that ∆t is larger than the longest life-time of quasi-particles in the plasma, which is given by the inverse damping rate of the zero mode γ pl . This is sufficient to assure that no contributions from quasi-particles enter the definition of κ ∞ (t). We cannot, however, formally take the infinite time difference limit, but require ∆t t. This is done so that the dependence of the correlator EE (t, t ) on the central time is weak compared to its dependence on the time difference. We note that, as for the thermal case [78,79], 3κ ∞ (t) is the zero-frequency value of the Fourier transform of the force-force correlator shown in Fig. 1 The transport coefficient κ ∞ (t) enables us to formulate a Langevin equation for heavy quarks in analogy to thermal equilibrium. Since the integral of the forceforce correlator over a long but not infinite time interval t ∆t 1/γ pl is a constant, we can approximate the dynamics of the heavy quark at this timescale by a random momentum kick with the same normalization This leads to a physical picture of the dynamics of a heavy quark in the medium as a Langevin process, which is commonly used in phenomenological applications. Finally, we explain here how we extract κ (t, ∆t) from our simulations. In this paper we only consider sufficiently late times well within the self-similar regime so that we can stay in the limit ∆t t. Then one can, for the purpose of convenience in the numerical evaluation, replace our definition (22) with a version where one of the electric fields is always evaluated at the lower, instead of the upper, time Also averaging over the (lattice) volume leads us to the equation which we employ in our numerical extraction of κ (t, ∆t). Comparison of the plasmon mass scale to the frequency extracted from the unequal time electric field correlation function κ (t, ∆t) given by (30). The plasmon mass scale is extracted using Eq. (4). The oscillation frequency of κ (t, ∆t) is extracted by fitting it to the damped oscillator in Eq. (31). The main observation is that the frequencies are closely related.

D. Self-similar behavior of κ (t, ∆t)
In Sec. III B we have seen that after a quick initial growth, p 2 (t, ∆t) grows more slowly with time ∆t, approximately linearly, and involves damped oscillations. We summarized these observations under the features of damped oscillations (ii) around a linear growth (iii). Let us now study these properties in more detail in terms of κ (t, ∆t). We recall that this quantity can be thought of equivalently as the time derivative of p 2 (t, ∆t) or as the integral of the electric field correlator over the time difference up to ∆t, and that its ∆t → ∞ limit is the heavy quark diffusion coefficient.
The ∆t-dependence of κ (t, ∆t) at different times Qt = 1500, 3000, 5000 is shown in Fig. 3. The inset of Fig. 3 shows our result scaled by the constant hard scale Q only. The main plot shows the correlator divided by the heavyquark diffusion coefficient κ ∞ (t), i.e. the ∆t → ∞ limit. As we will show in Sec. V, the time dependence of κ ∞ (t) can be well described functional form of (Qt) −5/7 times a logarithm of time that we will motivate in Sec. IV (see the explicit expression given in Eq. (64)). The correlators in the main plot are plotted as a function of the time difference scaled by the plasmon frequency, ω pl (t)∆t. The plasmon frequency used for the rescaling is computed using the HTL formulas (4) and (5), which amounts to effectively rescaling the horizontal axis with a power law (Qt) −1/7 . As explained in Sec. II, the different values of t correspond to different ratios of the physical scales in the problem. The fact that the curves from different times t overlap as functions of rescaled time, clearly shows that the oscillations in κ (t, ∆t) happen at a scale determined by the plasmon frequency. The scaling of the amplitude of the oscillations in Fig. 3 shows that this same time dependence also describes the amplitudes of the oscillations as a function of ∆t, a sign of self similar evolution in κ (t, ∆t).
The oscillatory form in Fig. 3 can be fitted for ∆t 1/Λ separately for each t to a damped harmonic oscillator with a constant offset term From this fitting procedure, we extract the frequency ω fit (t) and show it in Fig. 4 together with the computed frequency ω HTL pl as functions of time. Since the frequencies are quantitatively close to each other, we conclude Also the heavy-quark diffusion coefficient κ ∞ (t) is extracted using the fit in (31). According to Eq. (26), the coefficient κ ∞ (t) is defined as the late relative time limit ∆t 1/γ pl of κ (t, ∆t). For finite ∆t, it corresponds to the offset κ fit ∞ (t) of the oscillations visible in Fig. 3 and incorporated into the fit function (31). Measuring κ ∞ (t) as κ fit ∞ (t) reduces the residual dependence on ∆t and will be used as our standard way to extract κ ∞ (t).
We have also studied lattice regularization effects on our extraction of κ ∞ (t), with results shown in Appendix A . We find that κ is insensitive to the IR cutoff. In the case of the UV cutoff, the results start to drift considerably for lattice spacings larger than Qa s > 0.6 and hence, we use smaller lattice spacings. The Debye scale stays well within the reach of our lattice at all times.

IV. UNDERSTANDING THE TIME DEPENDENCE OF THE CORRELATOR
To understand the observations in the previous section in terms of microscopic degrees of freedom in the system, we construct here two models for the ∆t-dependence of κ (t, ∆t) and compare them to our numerical results. A crucial role in our discussion here is played by the momentum space gluonic equal-time correlation function EE (t, t, p) in Coulomb gauge, which we interpret in terms of a single particle distribution of gluons. The input in the models that we want to construct is this single particle distribution that we extract from our simulation. This is an equal time correlator, meaning that it represents an integral over frequencies. It is also needed as a function of momentum, which means that in coordinate space it is not local. This single particle distribution is also not manifestly gauge invariant, but evaluated using field configurations in the Coulomb gauge. From this information we want to construct the heavy quark diffusion coefficient, which is local in coordinate space, i.e., involves an integral over gluonic momenta. The heavy quark diffusion coefficient is the zero frequency limit of an electric field correlator, i.e., requires correlators at unequal times.
We will here use two different approaches to go from equal time -unequal coordinate correlations to an unequal time -equal coordinate one. The basic idea of the first one, that we call here the spectral reconstruction (SR) method, is to assume that the spectral functions in the general momentum-frequency space are the ones given by (HTL) perturbation theory, and that they are related to the statistical function by a generalized fluctuation-dissipation relation. This connection enables us to relate the statistical functions in different parts of phase space to each other, via the intermediary of the spectral function. The underlying assumptions are backed up by our previous numerical results in Ref. [72]. The second one, of course related to the first one in the appropriate parametric regime, is to use known perturbative calculations of the heavy quark diffusion coefficient in kinetic theory, and simply substitute our measured single gluon distribution in such a calculation. We will first construct these two models in this section, and then compare them to the numerical result for the electric field correlator κ (t, ∆t) in Sec. V.

A. Equal time electric field correlator
We start here with the statistical correlation function. In general it is defined as the anticommutator of Heisenberg field operators with x ≡ (t, x). In the classical-statistical approximation it becomes the expectation value of the product On a periodic lattice the Fourier transform with respect to the relative coordinate x − x averaged over the whole lattice can conveniently be computed by Fourier transforming the electric fields as EE (t, t , p) = E a j (t, p), (E a k (t , p)) * /V . We refer to its Fourier transform to frequency space as EE jk (t, ω, p), neglecting the difference between fixed t and fixedt = (t + t )/2 as discussed in Sec. III.
It is necessary to distinguish transverse and longitudinal projections, which are defined as 2P T jk = δ jk −p j p k /p 2 and P L jk = p j p k /p 2 , respectively, where 2 counts the number of transverse polarizations. Correlators can then be decomposed into polarizations In our previous publication [72] we observed that the equal-time statistical correlation function EE T,L (t, t, p) is enhanced compared to HTL expectations at low momenta for both polarizations. We show the numerically extracted EE T,L (t, t, p) correlators in Fig. 5 in the selfsimilar regime at time Qt = 1500 as solid curves. The figure also shows in dashed lines a fit to these numerical results. To enable a more efficient evaluation of some integrals appearing below in our models for the timedependent correlator, we will in practice use these fits instead of our original numerical data. Moreover, using EE T (t, t, p), we can define the distribution function as in our previous publication [72] as This is the definition that we used to extract the values (9) in Sec. II, with an iterative procedure to simultaneously extract both f (t, p) and m using Eqs. (4) and (36).
In HTL theory at leading order, the electric field correlator would, for low momenta p m D , be expected to approach a constant Note that the second "≈" is in fact an equality at all p for the transverse polarization but only at p = 0 for the longitudinal one: for a discussion of the spectral function see Appendix B. The value of this constant, i.e., the effective temperature of the infrared modes 2 at low momenta, is conventionally denoted by T * . We construct a parametrization of this expected behavior by taking the parametrization of our data at high momenta and smoothly matching it to a constant value determined by the temperature T * calculated using Eq. (7). This "thermal IR" parametrization is shown by the dot-dash curves in Fig. 5. We emphasize that this second parametrization is meant to represent a scenario without the "infrared enhancement" seen in the correlator and discussed above, but keeping the large momentum degrees of freedom as close to the ones present in the lattice calculation as possible. We can then use both, the parametrization of our data including the infrared enhancement, and the one where it has been removed, to construct a model for the time-dependent correlators both with and without the infrared enhancement. Using this approach we will in fact argue below that the excess of gluons at low momenta is the main reason for the oscillations in ∆t observed in p 2 and κ.

B. Spectral Reconstruction (SR) method
The aim of this section is to develop a spectral reconstruction (SR) method, which is based on our previous measurements of the gluonic spectral functions, the EE T,L (t, t, p) equal-time correlators and expectations from perturbation theory, and that can be compared to our measurement of the unequal time electric field correlators, specifically to our numerical result for κ(t, ∆t). 2 One way to see that one expects a constant is to note that the equipartition of energy in thermal equilibrium at temperature T * for a classical noninteracting theory would correspond to an expectation value T * /2 for every quadratic term in the Hamiltonian, which in this case are (half) the squares of the components E a j (t, p).
For this, let us reformulate κ (t, ∆t) using the full correlation functionsρ(t, ω, p) and EE (t, ω, p), respectively. We discussed the statistical correlation function in the previous subsection. The spectral function (strictly speaking the time derivative of the spectral function, but we will callρ the spectral function in an abuse of language below) is defined as the commutatoṙ In the classical-statistical approximation, it can be computed using the retarded propagatorĠ R jk (x, x ) = θ(t − t )ρ jk (x, x ), which can be extracted from simulations with linear response theory [72]. We denote the Fourier transform ofρ with respect to relative time and spatial coordinates asρ jk (t, ω, p). Transverse and longitudinal polarizations can be distinguished as for the statistical correlator in Eq. (35). For thermal equilibrium the fluctuation-dissipation relation states that the ωdependence of these functions is the same, i.e., the ratio EE T,L (t, ω, p)/ρ T,L (t, ω, p) is a known function that only depends on the momentum p and the temperature. We can easily generalize this to a nonequilibrium situation by taking this function of p to be one determined by the equal-time statistical function (and the equal time spectral functionsρ T,L (t, t, p), whose expressions are written in Appendix B). Explicitly, we start by assuming that even out of equilibrium. We have indeed observed numerically in Ref. [72] that this relation holds in the selfsimilar regime.
With these considerations, we can write Eq. (23) as 3κ (t, ∆t) We now have to determine the spectral functions. Note that they appear as ratios, which are normalized to unity According to HTL calculations at LO, each spectral function can be decomposed in frequency space into parts that are associated with Landau damping and to quasiparticle excitations, resulting iṅ We have seen this same structure in our numerical calculations of the spectral functions in Ref. [72]. The quasiparticle contributions can be written aṡ in terms of the dispersion relations ω T,L (p) of transversely and longitudinally polarized quasi-particles, the residues Z T,L (p) and the functions h T,L (ω, p) that are normalized to unity and correspond to the quasiparticle peaks. Perturbatively, the damping rate of the quasiparticles is of the order g 2 T * . Thus, at leading order HTL the quasiparticle peaks correspond to delta functions h T,L (ω, p) → δ(ω). More realistically, the quasiparticle peaks exhibit a Lorentzian shape with finite damping rate γ T,L (p) which we have also verified numerically in [72]. The explicit leading order expressions for dispersion relations ω T,L (p), quasiparticle residues Z T,L (p) and Landau damping contributionsρ Landau T,L (t, ω, p)/ρ T,L (t, t, p) are written in the Appendix B.
Combining Eqs. (40) and (42), we can split the diffusion coefficient (40) into four parts, corresponding to transverse / longitudinal and to Landau / quasiparticle contributions. Note that for each of these contributions, the frequency integration simplifies. The Landau damping contributions only have support for |ω| < p, i.e., ρ Landau T,L (t, ω, p) ∝ θ(p 2 −ω 2 ). Thus the frequency integration becomes ∞ −∞ dω → p −p dω for these contributions. For the quasiparticle contributions the frequency integration can be done analytically as All the remaining integrals are in general performed numerically. We have observed deviations from the HTL expressions at LO for ω T,L (p), γ T,L (p) and EE T,L (t, t, p). Thus, in our numerical calculations of κ SR (t, ∆t) we will use the following forms, which we have extracted in our linear response framework [72,80]: • Since we observed for the dispersion relations ω T,L (p) some deviations from the respective HTL expressions at LO, we use fits to our data from [72].
• We use the damping rates γ T,L (p) extracted in Ref. [72]. However, as will be shown in Fig. 8, setting γ T,L (p) = 0 does not change the results considerably.
• For the statistical EE T,L (t, t, p) correlation function we observe significant enhancement over the HTL expectation in the infrared. To study this effect we use both, a parametrization of our numerical result for EE T,L (t, t, p) and one where this enhancement has been removed, as discussed above in Sec. IV A and shown in Fig. 5.
For the other ingredients needed in the calculation: for the functional form of the quasiparticle peak h T,L we use the expression (44), and for the Landau damping contributionρ Landau T,L (t, ω, p) and the quasiparticle residue Z T,L (p) we use the standard forms from the literature that can be found explicitly in Appendix B, with the value of the Debye mass obtained using the iterative procedure discussed above. Also note that in all of these expressions, a time dependence enters both due to EE T,L (t, t, p) and due to the time dependence of the Debye mass m D (t).
We will evaluate the full ∆t dependence numerically with the procedure described above. It is also important to compute κ ∞ (t), which emerges in the limit ∆t → ∞, i.e., the actual heavy quark diffusion coefficient. In this limit the sinc functions in the frequency integral (40) become delta functions lim ∆t→∞ 2 sin (ω∆t) ω → 2πδ(ω).
The transverse Landau damping contribution vanishes in the limit ω/p → 0 (see Eq. (B3) in Appendix B) and thus does not contribute. The longitudinal Landau cut, on the other hand, has a finite limit ω/p → 0 (see Eq. (B4) in Appendix B) and indeed gives a nonzero, in fact the only nonzero, contribution to κ ∞ (t).
One can see in two ways that the quasiparticle contributions vanish. One way is to simply use the delta function representation (46) and then note that the particle contribution toρ in Eq. (43) is proportional to ω. Alternatively, if one first integrates over ω as in Eq. (45), it is the quasiparticle damping term e −γ(p)∆t that vanishes in the ∆t → ∞ limit.
Thus, for the heavy quark diffusion coefficient κ SR ∞ (t) using the spectral reconstruction model method, we are left with only the Landau damping contribution, resulting in A straightforward evaluation yields κ SR ∞ (Qt = 1500) = 7.4 × 10 −5 Q 3 (48) κ SR ∞,th.IR (Qt = 1500) = 6.5 × 10 −5 Q 3 for our measured EE L (t, t, p) correlator and for the "thermal IR" parametrization in Fig. 5 where the infrared enhancement has been removed, respectively. To obtain a heavy-quark diffusion coefficient that is accurate to leading logarithmic order, we can use the simpler parametrization Then we would get where in the last line we have left only the leading logarithmic contribution. We leave this equation here as a reference since we will come back to it in the context of the kinetic theory expression discussed in the following.

C. Kinetic theory framework
We can also estimate κ ∞ (t) in the kinetic theory framework. Here the physical picture is quite intuitive: the heavy quark gains momentum from independent kicks by gluons in the medium. The scattering can be described by a perturbative Qg → Qg matrix element that, in the limit of large quark mass, is dominated by t-channel exchange of a gluon. To arrive at a scattering rate one additionally needs the gluon distribution for the incoming gluons, and the Bose enhancement factor for the outgoing ones (noting that for us f 1 so that we can neglect the unity in the Bose enhancement factor f + 1). We can easily obtain a quantitative expression for the diffusion coefficient following, e.g., the discussion in [18] as where M is the mass of the heavy quark. The incoming heavy quark momentum is (M, p) and the outgoing heavy quark momentum is (M, p ). The gluon momenta before and after a collision are given by (k, k) and (k , k ). The transferred momentum is given by q = p −p. Compared to the calculation in [18], which takes place in a thermal background, we have neglected scatterings from quarks and replaced the thermal gluon distributions with general ones.
The process is dominated by t-channel gluon exchange. We take the matrix element squared to be where C H = (N 2 c − 1)/(2N c ) is the color Casimir of the heavy quark and where we introduced an infrared regulator in terms of m 2 D . The angle between k and k can be expressed as Carrying out the integrals that are possible using the delta functions and an angular integral, we can transform this to This is now an integral that we can carry out numerically. Using the definition of the gluon distribution in Eq. (36), we get, again at Qt = 1500, the result κ KT ∞,th.IR (t = 1500) = 1.9 × 10 −5 Q 3 .
The first one of these is a little bit smaller than, but roughly in line with the results that we obtain from the SR method, see Eqs. (48) and (49). The second one is smaller by a factor of more than two. The effect of removing the infrared enhancement (i.e. the "thermal IR" approximation) is larger here than it is in the SR method. This is somewhat paradoxical, since conceptually the kinetic theory calculation is based on the integral over k in (56) being dominated by UV particle like degrees of freedom in f (k), whereas the expression (47) in the SR depends on the longitudinal correlator in the soft momentum region. The effect of the infrared enhancement is large enough that such parametric estimates start becoming unreliable quantitatively, even if they remain true at the leading logarithmic level, as we will discuss next. Let us now discuss how the KT and SR approaches are equivalent at the leading logarithmic order for large Λ/m D , again roughly following the discussion in [18]. We start from the first form in Eq. (55) that has an integral over both k and q. Now the momentum of the gluon k is typically of order Λ, while q is of the order of the Debye scale m D . In the limit Λ m D we can therefore assume that k q. This limit enables several modifications that decouple the k and q integrals. First, at leading order we can replace the upper limit in the q integral by the UV scale Λ. We can also drop terms O q 2 /k 2 in the integrand. In terms of the original kinetic theory this corresponds to a small-angle approximation where the scattering angle between k and k is cos θ kk ≈ 1. With these approximations we arrive at the same result as the leading logarithmic limit in the SR picture, Eq. (51) where we used the integral defining T * from Eq. (7) and the assumption that the longitudinal EE LL L (t, t, p) correlator has the form ∝ T * that we expect in HTL, Eq. (50).
To summarize, there is a common limit for our SR and KT models. From the KT side, it can be reached by systematically taking the limit of large Λ/m D , implying small angle scatterings and neglecting the exact kinematical limits for the scattering process. From the SR model reaching the common limit requires removing the infrared enhancement that we see in the equal time correlator and replacing it by a functional form that is proportional to ∼ T * which is determined by the hard modes (by the integral (7)). For the values of t that we have studied, i.e., the values of the scale separation Λ/m D , we see no indication of a disappearance of this infrared enhancement. However, even with the IR enhancement the leading logarithmic limit Eq. (51) should still be valid, the effect being on the value of the constant under the logarithm.
Incidentally, the equivalence of the common limit κ ∞,LL (t) in the KT and SR frameworks in Eq. (58) justifies a posteriori the form of IR regulation we used in Eq. (53). One arrives at the known leading logarithmic expression of κ KT ∞,LL (t) in thermal equilibrium with temperature T [18,27] by setting T * → T and Λ → T .
We can use the leading logarithmic limit to estimate the expected scaling behavior for the diffusion coefficient with t, which controls the magnitudes of the relevant scales Λ, T * and m D as discussed in Sec. II. Based on the scaling behavior of m D ∼ Q(Qt) −1/7 and g 2 T * ∼ Q(Qt) −3/7 in the self-similar regime, the scaling behavior of the diffusion coefficient is expected to be up to logarithmic corrections.
For a generalization to finite time κ KT (t, ∆t), we allow for a possibility of a nonzero frequency ω = k 0 −k 0 = k − k, corresponding to finite energy transfer between gluons and the heavy quark. Here one should in principle use the proper HTL expression. We, however, simplify the form by employing the same IR regulator as for κ KT ∞ (t) in Eq. (53). In this case the matrix element becomes and the frequency dependent expression can be written as We carry out this integral by a change of variables from q to k . The integrals over the angular variables can be carried out analytically. After Fourier transforming, our expression for the ∆t dependent diffusion coefficient is The remaining integrals are evaluated numerically. One recovers κ KT (t, ∆t → ∞) = κ KT ∞ (t). This procedure has enabled us to get a ∆t-dependent expression from the KT description.
We would only expect this result in Eq. (62) to describe the physics at large ∆t. In terms of the different physical ingredients included in the SR description, this finite-∆t calculation in kinetic theory still includes only the effect of the longitudinal Landau cut, and even that quite crudely. The gluon exchanged in the kinetic theory calculation is longitudinal and has a very spacelike fourmomentum, because the quark is infinitely heavy and the scattering gluon interacts only with its Coulomb field. As we will see explicitly in the following, the contribution of We also show the extracted values from the spectral reconstruction and kinetic theory frameworks in the infinite time limit. In addition, we show fits to the lattice data and kinetic theory whose form is given by (64).
the other parts of the spectral function: the quasiparticle poles and the transverse Landau cut, are essential to reproduce the full ∆t dependence. The effective Debye mass in the matrix element in Eq. (60) has been tuned such that it gives the leadingorder correct late-time limit but as it has a trivialized frequency structure it cannot capture the full leadingorder accuracy at finite ∆t. In order to arrive to full leading order description of the time evolution of κ one could further improve the estimate by replacing the simple mass regularization in Eq. (60) by the fully resummed HTL expression (as is done in, e.g., [17]).

V. RESULTS
Here we present our numerical results for the evolution of the heavy-quark diffusion coefficient κ ∞ (t) and the correlator κ (t, ∆t) and compare them to results using the SR and KT methods introduced in the previous section. With these tools, we are able to understand distinctive features of these quantities like the origin of the oscillations of κ (t, ∆t). The comparison to our full numerical data also helps us to assess the quality of the (mainly perturbative) SR and KT methods employed here.

A. Time dependence of κ∞(t)
We have extracted the values of κ ∞ (t) from the numerical calculation using the fitting procedure described above (see Eq. (31)). 3 The resulting κ ∞ (t) is shown in Fig. 6 as a function of t.
In Fig. 6 we also show the heavy quark diffusion coefficients computed using the SR and KT methods in the Eqs. (47) and (55), respectively, as discussed in the previous section. Both the SR and the KT methods have a time dependence that agrees rather well with our full numerical extraction. As for the normalization, the SR model slightly overestimates and the KT one underestimates it. Both methods agree with the lattice extraction with roughly 30 % accuracy. In the light of the discussion of the common leading logarithmic limit, Eq. (58), we could interpret the lower value for the KT model as being due to importance of the infrared enhancement: the KT model has an implicit assumption that the soft electric field modes have a thermal distribution with a temperature T * determined by the hard modes as in Eq. (7), while the SR model actually uses our measured correlator of electric fields at soft momenta, which is larger.
From the leading logarithmic limit we extracted an expectation for the scaling of the diffusion coefficient with time t in (59). The resulting power law, κ ∞ (t) ∼ t −5/7 , however, receives logarithmic corrections that can be large. Also at the earliest values of t, it is possible that κ ∞ (t) suffers from a larger transient effect from the system not yet being fully in the scaling regime. Using a simple t −5/7 times logarithmic fit to our data, the evolution of κ ∞ (t) can be effectively described as It can be instructive to compare our result for κ ∞ (t) to the corresponding thermal value κ therm ∞ . We emphasize that the classical-statistical framework here cannot be used to follow the system all the way to thermal equilibrium, since one has to remain in the regime of large occupation numbers. At the times where our system starts to approach a thermal equilibrium one, it will simultaneously fall out of the regime of validity of the classical statistical approximation that we are using here. To get a parametrical estimate, this will happen when occupation numbers are of order unity f (t therm , Λ) ∼ 1. Due to the scaling behavior g 2 f (t, Λ) ∼ (Qt) −4/7 this will occur at times that are parametrically Qt therm ∼ g −7/2 . Until this time, the non-equilibrium value is larger than the thermal one κ ∞ (t) κ therm ∞ by an inverse power of the coupling g 1.
Let us now try to make this comparison to a corresponding thermal system more quantitative. The choice to be made in comparing our result to thermal equilibrium is to define what is meant by a "corresponding" thermal system. One way to perform such a comparison is to choose a value for the coupling constant g and the energy density ε and keep them the same for the two systems. For a thermal system these are enough to determine the temperature T and other physical quantities. In our case, the combination g 2 ε of the coupling and the energy density is fixed by our initial conditions. But subsequently the system also has a time dependence that affects the magnitudes of the different physical scales in the system. Since we cannot follow our simulation to equilibrium we have to choose a meaningful value of t for a comparison to a thermal result. This choice should be made by looking at the value of some other physical quantity.
Since a defining feature of both the thermal system and our overoccupied cascade is the scale separation between a hard scale and a soft Debye scale, a natural way to compare to a thermal system is to compare the Debye scales. At fixed g 2 ε we can do this in two ways. One option is to compare our system to a thermal one at the same value of the Debye mass, or equivalently the ratio m 2 D / √ ε. For a thermal system the energy density and Debye mass are related to the temperature by ε = (N 2 c − 1)π 2 T 4 /15 and m 2 D = (N c /3)g 2 T 2 . We can quantify this Debye scale comparison by defining an "effective coupling" with the ratio m 2 D / √ ε as where we have used N c = 2 as in our numerical simulation for the explicit value. Another option is to target similar values of the scale separation between the hard and soft scales. This scale separation can be parametrized by the ratio m D /Λ. For a thermal system with N c colors and no dynamical fermions the energy density, coupling, Debye and hard scales are related by ε = (N 2 c − 1)π 2 T 4 /15 ≈ 1.97 T 4 ≈ 0.031 Λ 4 and m 2 D = (N c /3)g 2 T 2 . Note that here we define, as in Sec. II, the hard scale Λ as the momentum scale for which the integrand of the energy density is maximal. For a thermal system this condition (d/dp (p 3 f BE (p)) = 0) leads to (Λ/T ) = 3(1 − exp(−Λ/T )) and consequently Λ ≈ 2.82 T . We again parametrize the scale separation m D /Λ in terms of another effective coupling with again N c = 2 for the numerical value. Both of the effective couplings (66) and (67) are constructed so that for a thermal system at leading order in g they just give back the gauge coupling:g 2 ε =g 2 Λ = g 2 . The thermal calculation at leading order (e.g. [19]) gives (with N c colors and no fermions) a value κ therm = C F g 4 T 3 /(18π)N c (ln(2T /m D )−0.64718). In terms of the energy density and the scale separation ratio this gives us where in the second equality we have taken explicitly N c = 2. For a comparison of our overoccupied system results with this number, we now want to extrapolate our result to a value of t where our measured time dependent Debye scale has a value corresponding to ag that is the same as in a thermal system, i.e.,g → g. Since our results are usually expressed in terms of Q, we start by relating this quantity to the energy density by g 2 ε = 2(N 2 c − 1) p p g 2 f (t = 0, p) ≈ 0.0762 Q 4 . Thus, we have Q 3 ≈ 6.9 (g 2 ε) 3/4 and Q = 1.90 (g 2 ε) 1/4 . Using the scaling laws in Sec. II and the numerical extractions in Eq.
We can then express our result (64) for κ ∞ (t) in terms ofg ε as or in terms ofg Λ analogously as κ ∞ (t) ≈ 0.0047 ln 1 g Λ + 0.177 g The expressions (70) and (71) involve different powers of the effective couplings due to the fact that the definitions (66) and (67) use ε and Λ that are, in our overoccupied cascade, proportional to g −2 and g 0 . It is interesting to see that the coefficients in (70) and (71) are almost the same. This is a consequence of the values of the Debye and hard scales in Eq. (9) fortuitously being such that they match the thermal values when extrapolated to late t. In other words, it is a coincidence that the conditions g ε = g andg Λ = g are satisfied at not just parametrically, but also numerically at the same values of t. Now, we can make a comparison of the diffusion coefficient in the overoccupied gluonic cascade to that in a thermal system by comparing (68) to (70) and (71). We firstly re-emphasize that in the overoccupied phase of the evolution, where our classical description is justified,g ε g andg Λ g; thus the value of κ ∞ (t) is parametrically larger than for a thermal system with the same energy density. We can then extrapolate the time evolution of our system to occupation numbers f ∼ 1, i.e. the parametric region Qt ∼ g −7/2 where the occupation number becomes order one and the classical statistical approximation breaks down. This is done in practice by settingg ε → g,g Λ → g. With these assignments we observe two features from Eqs. (68), (70) and (71). Firstly, as a consistency check, the powers and logarithms of g are the same in all three cases. Secondly, the value of the coefficient in (70) and (71) is smaller by a factor of ∼ 3. This is mostly a result of the IR enhancement and the way we defined the effective couplings using m D .
To be more precise: the change in the functional form of f (p) with the IR enhancement increases the value of m D by 25% (from a comparison between our result and the "thermal IR" parametrization). Since in our classical simulations we have m D ∼ (Qt) −1 /7 , we have to wait for (Qt) −1 /7 to be 25% larger to reach the limitg ε =g Λ = g, compared to a case without the IR enhancement. During this time, κ(t) ∼ (Qt) −5 /7 becomes smaller by a factor of 1.25 5 ≈ 3. In other words, while the IR enhancement increases the value of κ(t), it increases the value of m D even more, relative to a thermal distribution. This has the effect that when we express the g in κ ∼ g 4 ε in terms of m D , the coefficient is smaller in our case than for a thermal distribution. As long as values are consistently expressed in terms of the appropriate physical scales of the problem, our overall result is thus quite consistent with the perturbative understanding of the microscopic picture of heavy quark diffusion, apart from these effects of the IR enhancement.
B. Understanding the ∆t dependence of κ (t, ∆t) The self-similar t dependence of the time derivative of p 2 (t, ∆t) , i.e., of κ (t, ∆t) has been discussed in Sec. III D. We also showed there that its dependence on ∆t is well described by a damped oscillating function with oscillation frequency ω pl shifted by κ ∞ (t), which was parametrized in Eq. (31). Here we study the origin of these oscillations by comparing our data to computa- Comparison of the different contributions to κ(t, ∆t) computed in the SR framework, where the EE with extracted IR enhancement is used. We observe that the source of the oscillations at the plasmon mass scale are the quasiparticle contributions, while neglecting the damping rates by setting γ = 0 has little effect on the oscillations. On the other hand, the major contribution to the heavy-quark diffusion coefficient κ∞ arises from the longitudinal Landau cut, around which the oscillations proceed.
tions within the SR and KT methods. The dependence of κ (t, ∆t) on ∆t is shown in Fig. 7 for Qt = 1500. The black solid curve shows our data while the other curves are obtained using the SR and KT models, Eqs. (40) and (62). The blue dashed curve corresponds to the SR model when using the parametrization of the equal time statistical correlation function which agrees with our extracted EE (t, t, p) correlator in the infrared. The green dotted curve corresponds to the SR calculation in which we use the "thermal IR" equal time statistical correlation function, where we have forced the EE (t, t, p) correlator to a constant ∼ T * in the infrared, see Fig. 5. We observe that the SR method quite well reproduces the oscillations seen in the data, where the shift to slightly larger values of κ can be attributed to the overestimation of κ ∞ (t) by the SR method, as observed above. However, if the "IR enhancement" in the equal time correlation function is removed, the oscillations practically disappear. Likewise, the oscillations seen in the data are not present in the KT model calculation, as could be anticipated based on the discussion in Sec. IV C. We take the comparison with the SR model with and without the infrared enhancement in the equaltime correlator as a confirmation of the existence of this infrared enhancement from a manifestly gauge invariant observable. This is one of the main conclusions of our paper.
We can also break down the SR model into different components contributing to κ SR (t, ∆t), as mentioned in Sec. IV B, to get an idea of their relative contributions. This decomposition is shown in Fig. 8, which includes the total κ SR (t, ∆t) together with the contributions from transverse and longitudinal quasiparticles and the transverse and longitudinal Landau cut. We can make two main observations from this figure.
Firstly, as discussed earlier, the only contribution to the ∆t → ∞ limit comes from the longitudinal Landau damping which, however, does not show any oscillations (the integration weight being mostly at very small ω in Eq. (40)). The transverse Landau cut only contributes at quite small ∆t, and does not exhibit oscillations.
Secondly, the oscillations are produced by both the transverse and longitudinal quasiparticle contributions, which, however, vanish at ∆t → ∞. Their amplitudes decrease due to the specific momentum integration in Eq. (40) and because of a finite damping rate γ T,L (p). The latter seems to have a smaller effect on the evolution. This can be seen in Fig. 8 by comparing κ SR (t, ∆t) to the curve κ SR γ=0 (t, ∆t), where the damping rates have been set to zero in its calculation. Therefore, even LO perturbative calculations, where quasiparticles are given by Delta functions in the spectral function, lead to damped oscillations with frequency ω pl if the IR enhancement is taken into account.

VI. CONCLUSIONS
In this paper we have measured the heavy quark diffusion coefficient κ ∞ (t) in a self-similar overoccupied cascade system, using numerical calculations in classical glu-odynamics. In addition to calculating the diffusion coefficient, we have observed strong oscillatory structures in ∆t in local gauge-invariant electric field correlators κ (t, ∆t) at the plasmon frequency scale.
The correlator κ (t, ∆t) is, in physical terms, the derivative of the momentum broadening p 2 (t, ∆t) of a heavy quark traversing the plasma with respect to the time ∆t. In general, we have observed that p 2 (t, ∆t) rises fast initially on a time scale of ∆t ∼ 1/Q, where Q is a typical hard scale of the plasma. This initial rise is followed by an approximately linear growth that can be interpreted as due to momentum diffusion. The linear growth is modified by additional damped oscillations with the plasmon frequency. While the initial rapid growth is a direct consequence of the decoherence of hard quasiparticles with momenta ∼ Q, in order to understand the observed oscillations, dynamics at lower momenta and frequencies have to be taken into account.
Prompted by these numerical observations, we have constructed two models to understand them, which start from a distribution of gluons in the system, obtained numerically from an equal-time correlator of electric fields. In the first method, that we call spectral reconstruction, we extrapolate our measured equal time correlators from (∆t = 0, p) to ω-space using a generalized fluctuationdissipation relation and assuming a HTL structure for the spectral function. This structure, which we have also observed numerically in [72], includes both damped quasiparticle peaks and a low-frequency Landau damping region. Both contributions are needed to explain the full ∆t dependence of the electric field correlators κ (t, ∆t). We found that the oscillations at the plasmon frequency are produced by the quasiparticle contributions, whereas the main contribution to the diffusion coefficient arises from the Landau damping of longitudinal modes in the plasma. In our second method we use a kinetic theory calculation where the diffusion coefficient is computed utilizing the fact that the process is dominated by t-channel gluon exchange.
Using these models we have argued that the oscillations at the plasmon scale can be explained by a larger occupation of infrared modes p m D in the system than expected from perturbation theory. This feature, that we call IR enhancement, had been observed earlier in the equal-time electric field correlator in Coulomb gauge. We have now demonstrated an independent confirmation for it in a gauge-invariant observable d p 2 (t, ∆t) /d∆t in a different part of phase space. These oscillations are also visible in the momentum broadening observable p 2 (t, ∆t) .
The time-dependence of physical quantities in the scaling system tends to follow simple, analytically derivable power laws as a function of time. We have seen that the time dependence of the heavy quark diffusion coefficient is consistent with the t −5/7 × log t-dependence expected from scaling arguments. Due to the IR enhancement, the contribution of small momentum quasiparticle-like modes below the plasmon frequency scale is larger than one would expect in perturbation theory. Recall that a crucial feature of our (SR) analysis was to assume an HTL form for the spectral function, which is not proportional to the number of quasiparticles in the system, together with a generalized fluctuation-dissipation relation. This means that the IR enhancement leads to an enhancement of the longitudinal statistical correlation function in the Landau damping region ω p over the perturbative expectation. Consequently the heavy quark diffusion coefficient is larger than one would perturbatively expect based on just the number of hard gluons in the system. In the case of heavy quark diffusion this enhancement is not a large correction to the the leading logarithmic behavior of the diffusion coefficient, but shows up as an oscillatory deviation from the diffusive behavior. In the kinetic theory framework the IR enhancement is a substantial correction. We take this as an indication that kinetic theory becomes less reliable due to the IR enhancement.
A similar argumentation can also be used to related transport coefficients like the jet quenching parameterq and the momentum broadening resulting from it. Therefore, we believe that the effects of overoccupied IR gluonic quasiparticle modes should be taken into account for calculations of other transport coefficients. Similar effects could be important also in different gluonic systems, including different initial states, anisotropies or in expanding geometry, if the underlying plasma involves infrared enhancement. It would therefore be interesting to study their effect, and thus the effects of such an infrared enhancement, on other phenomenological observables.
Our calculation has been done in an extremely weak coupling, classical field limit. To make a connection to phenomenologically relevant values of the parameters, one needs to scale this to realistic values of the coupling, keeping a meaningful set of physical quantities constant. Unsurprisingly, in the strongly overoccupied system the diffusion coefficient is enhanced by inverse powers of the coupling constant compared to a thermal system at the same energy density. We have done a more detailed comparison to a thermal system by extrapolating our overoccupied cascade system to a situation where the scale separation between the energy density and the Debye scale are similar. The results of this exercise, obtainable from a comparison of Eqs. (70) and (71) to (68), are intended to be usable by a brave phenomenologist to estimate the effect of an initial overoccupation of gluonic modes in the pre-equilibrium stage of a heavy ion collision on calculations of heavy quark diffusion. tent 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. J. P. acknowledges support for travel by the Jenny and Antti Wihuri Foundation. 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 ).

Appendix A: Lattice checks
The sensitivity to the lattice UV cutoff depends on the lattice spacing Qa s . The dependence of our extraction of the heavy quark diffusion coefficient at Qt = 1500 on the cutoff is shown in Fig. 9. The results remain relatively stable up to values Qa s = 0.6. For larger lattice spacings, we start to observe significant lattice UV cutoff effects, and for Qa s = 0.8 the observed value for κ has been nearly doubled compared to Qa s = 0.3. The main conclusion in this case is that one has to use lattice spacings Qa s ≤ 0.6. This also means that our standard choice Qa s = 0.5 is sufficiently small.
The infrared cutoff on the lattice is given by the finite size of the system L. The dependence of the extracted heavy quark diffusion coefficient at Qt = 1500 on this cutoff is shown in Fig. 10. The main result is that κ ∞ (t) is relatively independent of the IR cutoff, only for lattice sizes QL < 20 we see that our measurements start to break down. Typically our IR cutoff is QL = 130 (260 3 lattice with Qa s = 0.5), which is sufficiently large.
We want our lattice to be large enough that the Debye scale is resolved on the lattice. In our simulations we have m D ≈ 0.21 Q at Qt = 1500. The smallest momentum mode we have available is k min = 2π L ≈ 0.1 Q for a 128 3 lattice. Thus even for 128 3 lattices we can still ac- commodate the Debye scale. Because of the self-similar scaling, we have m D ∼ t −1 /7 and ultimately at some point the Debye scale will fall below our reach. Due to the slow power law evolution we expect this to happen roughly at Qt = (1500 × 2 7 ) ≈ 190000, which is much later than our simulation times.

Appendix B: HTL functions
In this appendix we will briefly recap the functional forms of the HTL functions which are necessary for the estimation of κ, in particular for the decomposition Eq. (42). For more details on these we refer the reader to textbooks on thermal field theory, see e.g., [81] (with a slightly different notation than here). Previously in [72] we have studied the spectral properties of classical Yang-Mills theory in the self-similar regime, and these functions and their characteristics are discussed there in more detail.
Equal time spectral functions are determined by the sum ruleṡ The transverse and longitudinal Landau cut contribu-tions are given by the functionṡ ρ Landau T (t, ω, p) ρ T (t, t, p) = π p m 2 where x = ω p . The residues of the quasiparticle peaks are given by the functions (B6) Next we will go through the explicit forms of the dispersion relation ω T,L (p). In temporal gauge, the retarded transverse and longitudinal propagators are given by where the transverse and longitudinal gluon self-energy tensors in the HTL framework at LO read Here Q 0 is the Legendre function of the second kind The dispersion relations of the transverse and longitudinal quasiparticle modes are given by the poles of the retarded propagator. In the general case one has to solve these numerically. It is, however, possible to find approximate solutions for small and large momenta. For the transverse dispersion relation we have ω T = ω 2 pl + 6 /5 p 2 , p m D (B10) while for the longitudinal dispersion relation the corresponding expressions are ω L = ω 2 pl + 3 /5 p 2 , p m D (B12) For the damping rates γ T,L (p), we parametrize the data we have previously extracted (see [72], Fig. 9). Since we could not observe any differences between transverse and longitudinal damping, we use the same parametrization for both polarizations. In perturbation theory the damping rate is a next to leading order effect. Perturbatively it has been computed at p = 0 [76] and also for momenta of the order of the hard scale [82].