Spectral function for overoccupied gluodynamics from real-time lattice simulations

We study the spectral properties of a highly occupied non-Abelian non-equilibrium plasma appearing ubiquitously in weak coupling descriptions of QCD matter. The spectral function of this far-from-equilibrium plasma is measured by employing linear response theory in classical-statistical real-time lattice Yang-Mills simulations. We establish the existence of transversely and longitudinally polarized quasiparticles and obtain their dispersion relations, effective mass, plasmon frequency, damping rate and further structures in the spectral and statistical functions. Our new method can be interpreted as a non-perturbative generalization of hard thermal loop (HTL) effective theory. We see indications that our results approach leading order HTL in the appropriate limit. The method can also be employed beyond the range of validity of HTL.

We study the spectral properties of a highly occupied non-Abelian non-equilibrium plasma appearing ubiquitously in weak coupling descriptions of QCD matter. The spectral function of this far-from-equilibrium plasma is measured by employing linear response theory in classical-statistical real-time lattice Yang-Mills simulations. We establish the existence of transversely and longitudinally polarized quasiparticles and obtain their dispersion relations, effective mass, plasmon frequency, damping rate and further structures in the spectral and statistical functions. Our results are consistent with hard thermal loop (HTL) effective theory but also indicate effects surpassing its leading order. The method can be employed beyond the range of validity of HTL.

I. INTRODUCTION
Strong Yang-Mills fields are a ubiquitous feature of weak coupling descriptions of Quark-Gluon Plasma and heavy ion collisions. They are created through different mechanisms and affect a broad range of phenomena. For instance, in thermal equilibrium the infrared tail of thermal Bose distribution f ∼ T /ω makes low frequency modes highly occupied, corresponding to strong, long-wavelength fields. The interaction of these fields with hard modes at the scale T becomes nonperturbative at the soft or asymptotic mass scale m ∼ (α s d 3 p f /p) 1/2 ∼ gT . The physics of these soft modes affects the weak coupling expansion of the equation of state starting from the O(g 3 ) level [1,2]. These same modes contribute to the thermal photon and low-mass dilepton production rates already at leading order (LO) because of the infrared sensitivity of the t-channel exchange [3]. Similarly the effective kinetic theory description used to describe near-equilibrium transport in weakly coupled QCD is sensitive to the physics of the soft modes already at LO [4].
Strong classical Yang-Mills fields also appear in the weak coupling description of the initial stages of heavyion collisions. The initial condition of the post-collision debris in the midrapidity region is dominated by nonperturbatively strong gluon fields at the characteristic momentum scale Q s with mode-occupancies f ∼ 1/g 2 [5,6]. Once the dominant modes of these fields have diluted because of the combined effect of expansion and nonperturbative self-interactions, the subsequent approach towards local thermal equilibrium is described by the effective kinetic theory that is, again, sensitive to strong infrared fields [7].
In most of the above mentioned cases, there is a scale separation between the highly occupied soft (p ∼ m) modes, and the hard (p ∼ Λ) modes with which they interact nonperturbatively. In this case, the real-time dynamics of the soft modes can be followed in the Hard (Thermal) Loop (HTL) effective theory [8,9] in which the soft classical fields are coupled to hard ballistically propa-gating point-like color charges that interact with the soft classical fields through colored Vlasov-Wong equations [9]. The Hard Loop dynamics can then be solved perturbatively in the self-interactions of the soft modes (while still treating the interaction with the hard particles nonperturbatively, i.e., performing the HTL resummation). The small parameter in this expansion is given by the scale separation m/Λ, which in thermal equilibrium is gT /T . This expansion parameter can be significantly larger than that of the hard sector (α s /4π), and the extent to which the weak coupling expansions of the various quantities are reliable is dependent on how well the soft sector is described by the LO HTL expressions.
Computing complicated observables and higher order corrections in HTL theory is technically involved and is a major limiting factor in the advancement of the program of thermal weak coupling calculations [10,11]. For space-like correlation functions powerful Euclidean techniques exist [12,13]. These techniques have made high order calculations of the equation of state at finite temperature possible [14], and even allowed for a numerical evaluation of the full contributions arising from the soft sector through the simulation of an effective theory [15]. They have also enabled several next-to-leading order (NLO) computations, including the photon and smallmass dilepton production rates [16,17] among other results.
Lattice formulations of HTL in real time (see e.g. [18][19][20][21][22][23][24]) rely on an explicitly separate description of the soft modes as gauge fields and of the hard modes as classical particles. There are two possible approaches used in the literature. One is to initialize the fields with a (classical) thermal distribution, where most of the field energy resides in the modes close to the lattice UV cutoff 1/a s . In this case the inability to renormalize this system prevents following the time evolution of the soft fields in a controlled way and prohibits a numerical evaluation of time-like correlation functions (see however [22]). The other option used e.g. in studies of plasma instabilities [25][26][27] is to initialize the classical field with a sufficiently UV safe distribution. In this latter case, in order to correctly resolve the soft physics free of cutoff artefacts, one needs m 1/a s . On the other hand one needs to simultaneously have Λ 1/a s in order for the hard modes to be sufficiently localized to justify their description as classical particles. These requirements make the simulations rather expensive. More importantly they make it impossible to study systems with a smaller scale separation, i.e. larger values of m/Λ, and thus inherently difficult to extend the HTL theory beyond strict leading order. We would like to argue here that a formulation in terms of classical fields and linearized fluctuations provides a method to describe a similar physical system in a way which can also be used in a regime where the scale separation m/Λ need not be small. This enables one to go beyond the strict leading order in m/Λ also in nonperturbative lattice calculations for time-like correlation functions.
The aim of this paper is twofold. First, we study a specific isotropic and over-occupied far-from-equilibrium system that, starting from an overoccupied initial condition, is undergoing a cascade of energy towards the UV in a self-similar regime [28][29][30][31][32][33][34]. This system exhibits an increasing separation of scales between a hard scale Λ that dominates the energy density, and the softer scale ∼ m where the interaction with the hard scale becomes non-perturbative. This system is similar to thermal equilibrium at weak coupling in the sense that there is a scale separation m Λ (or gT T ). On the other hand, it differs from a thermal system in the sense that this scale separation is not given by the coupling, but rather increases with time as m/Λ ∼ (Qt) −2/7 , where Q is a momentum scale constant in time. In this system, also the hard sector is over-occupied and hence, the time evolution can be followed numerically within a classical Yang-Mills simulation. We compare the expectations of LO HTL theory to the full numerical time evolution. Doing so, we confirm and quantify to what extent HTL at LO is a good approximation to describe the physics of soft modes of this far-from-equilibrium system. Similar comparisons have been done for scalar theory [35][36][37], while for gauge theory previous studies of the quasiparticle dispersion relation [31,[38][39][40] have only used the behavior of the background field (i.e. the statistical function) at equal times.
Secondly, we turn the argument around and interpret the classical simulation as a non-perturbative simulation of HTL. As HTL theory is insensitive to the detailed form of the hard sector, all isotropic equilibrium and nonequilibrium systems exhibiting the large scale separation are related within HTL theory. Therefore the results obtained from the controlled classical simulation of classical fields can be directly applied to thermal equilibrium where no other controlled non-perturbative methods exist. As an application we determine the plasmon damping rate generalizing the classic result from Braaten and Pisarski [41] to finite momentum.
We will start in Sec. II by introducing the considered theory, relevant correlation functions, predictions from the HTL theory and our numerical method. Our results are shown in Sec. III. We will conclude in Sec. IV.

A. Classical Non-Abelian gauge theory
We consider an SU(N c ) gauge theory 1 with N c = 2 in temporal A 0 = 0 gauge. The classical equations of motion in continuum are given by with the gauge field A i and the (chromo-)electric field E i , (spatial) vector components i, j = 1, 2, 3 and the adjoint group indices a, b = 1, . . . , N 2 c − 1. The covariant derivative and the field strength tensor are given by D ab Here g is the coupling constant and the structure constants f abc for SU (2) are given by the totally anti-symmetric Levi-Civita tensor f abc = abc .
The fields are discretized on a cubic lattice with N s lattice sites with lattice spacing a s in each of the 3 dimensions. In order to preserve gauge invariance, one uses link variables U j (t, x), that are related to the gauge fields by where Γ a are the generators of the su(N c ) algebra normalized in the standard way, i.e. Tr Γ a Γ b = 1 /2 δ ab . The time evolution then follows from the discretized equations of motion where the plaquette is given by , withî and unit vectors in the i, j directions. The antihermitian traceless part of a matrix V is defined as The discretized equations of motion preserve the Gauss law condition at every time step where the discretized covariant (backward) derivative is (suppressing the time variable for brevity) given by )/a s . The parallel transporting link field in adjoint representation reads U ab . To work in Fourier space, we use the common definition of the spatial Fourier transform of a field w as w(x) = 1 V p e ix·p w(p), with the volume V = a 3 s N 3 s . The discrete momenta are given by p j = 2πj/(a s N s ). On the other hand, the backward derivative is defined We define the (complex-valued) backward momentum on the lattice as its eigenvalue The forward derivative leads to the momentum p F j = (p B j ) * and the discretized second derivative −∂ 2 s sin 2 (p j a s /2). The magnitude of a lattice momentum will be denoted by p = |p B | = |p F | ≡ j p 2 j in the following. This allows us to define a basis of vectors v (λ) (p) that will be referred to as polarization vectors v (λ) (p) with polarization λ. The longitudinal polarization λ = 3 is v (3) (p) = p F /p, while the remaining vectors with λ = 1, 2 are orthonormal to it and to each other and are commonly referred to as transversely polarized. One can now define transverse and longitudinal projection operators 2 as while one has 2P T ij (p) + P L ij (p) = δ ij . The gauge and electric fields are initialized in momentum space with only transverse polarizations as with complex Gaussian random numbers that satisfy and similarly forc, while c * c cl = 0, where · cl denotes a classical average over the distribution of these random numbers. We choose an isotropic initial single-particle distribution The fields constructed in this way are not guaranteed to satisfy the Gauss law (4), which thus has to be separately imposed by projecting back to the constraint surface using the algorithm in [45]. The form (9) of the initial conditions, with the initial amplitude n 0 and momentum scale p 0 , is expected to be close to the attractor distribution encountered in highly occupied plasmas [34]. Moreover, for weak coupling g 1 this initial condition guarantees high occupation numbers and we always compute g 2 f , whose evolution is independent of the coupling constant in classical-statistical simulations since g drops out of the classical equations of motion for these initial conditions. For convenience, we also define a characteristic energy scale where is the energy density. Unless stated otherwise, all dimensionful quantities will be expressed in terms of Q.
At output times t > 0, we define the distribution function with transversely polarized fields as and employ f (t, p) = f EE (t, p) as our standard definition of the distribution function. In all definitions the distribution function is averaged over transverse polarizations and over adjoint gauge components, where d A = N 2 c − 1 is the dimension of the adjoint representation. Since the considered systems are isotropic, we additionally average over momentum modes with the same magnitude p when computing observables, which improves the statistical accuracy of our results. The mass entering these expressions is the (asymptotic) mass m, computed iteratively as where we refer to it as m HTL to indicate its connection to the HTL formalism. While in the corresponding HTL expression at LO the mass does not occur on the right hand side, Eq. (12) takes some NLO corrections into account. The different distribution functions defined in (11) can lead to slightly different results in HTL computations. This is demonstrated at the example of the mass m HTL that is shown in Fig. 1 as a function of the iteration step for different definitions of the distribution function. In the first iteration, we use m HTL = 0, which corresponds to the LO formula for the mass. Since m 2 HTL + p 2 drops out of the calculation of the mass when f AA is used, the corresponding value does not change with the iterations, while m HTL computed with other definitions jumps to a lower value. One observes that the value for m HTL ranges between 0.145 Q and 0.156 Q, with f EE leading to a value in-between m HTL (Qt = 1500) = 0.149 Q.
We use m HTL computed with the distribution function f EE for all plots where m HTL is used in Sec. III because it typically provides values for the mass located between the ones from the other definitions. However, we emphasize that this value is not precise and other definitions provide values that differ by typically 5 %. Similarly, we use the spread in other HTL computed observables to estimate its error.

B. Spectral and statistical correlation functions
We are mainly interested in properties of two-point correlation functions, especially at unequal time, which encode information about the quasiparticle character, the excitation spectrum and the distribution of quasiparticles.
For this purpose, one can define the statistical correlation function and its double time derivative, which are symmetric and real-valued functions, as anticommutators of two field operators with x = (t, x) and where we used that ∂ tÂ b i (x) =Ê i b (x) (in the continuum). Since we consider spatially homogeneous systems, F does not depend on the central spatial coordinates (x + x )/2 but only on the relative coordinates ∆x = x − x and one can perform a spatial Fourier transform with respect to ∆x, arriving at F (t, t , p) and F (t, t , p). Further employing the system's isotropy, one arrives at F (t, t , p),F (t, t , p).
In classical-statistical simulations, the anticommutator of Heisenberg fields reduces to a product of classical fields 1 /2 {Â,Â } → A A cl . Hence, the definitions (14) become where we employed the operators P T /L jk from (6) to average over transverse or longitudinal polarizations but to simplify notation, we omitted the subscript T or L. Both F andF also include averaging over the d A = N 2 c − 1 adjoint gauge components. From Eq. (15) it becomes apparent that the statistical correlation function is closely related to the distribution function. Indeed, our definitions of the latter in Eq. (11) can be expressed as f EE (t, p) =F T (t, t, p)/ m 2 HTL + p 2 and similarly for the others.
Another important correlation function is the spectral function (and its time derivative) that can be defined as the commutator of two field operators Following the same steps as provided below Eq. (14), one arrives at The spectral function is anti-symmetric and real valued.
While the first relation is also valid for the longitudinal component lim t→t ρ L (t, t , p) = 0, the equivalent of the second relation will be discussed in Sec. II C. These transverse and longitudinal equal-time relations for the spectral function fully determine the initial conditions for its evolution. Note that this is different from the statistical correlation function whose initial conditions can be chosen by specifying an initial distribution f (t = 0, p).
Since in a classical-statistical framework commutators are mapped to Poisson-brackets (or Dirac-brackets in our case), a direct computation of the spectral function as in the case of the statistical correlator is more involved. Instead, we will use linear response theory for its computation, which will be discussed in Sec. II D. There it is described how one can compute the retarded propagator G R (t, t , p), which is closely related to the spectral function via Thus, having computed G R and exploiting the fact that for positive time differences it coincides with ρ, we simultaneously get the spectral function.
To discuss quasiparticle excitations, it is useful to transform the correlation functions to frequency space. This is done by first transforming the time variables tō t = (t+t )/2 and ∆t = t−t and Fourier transforming the correlation functions with respect to ∆t while keepingt fixed. Ideally, we would compute but in the interest of practicality we approximate this by (21) and always replacet by t. This approximation is justified when taking a sufficiently large time t such that ∆t max t andt ≈ t. Moreover, values at small ∆t typically give the dominant contributions to the integral because of the damping of oscillations in the correlation functions. Note that we have defined a real-valued spectral function in frequency space in Eqs. (20) and (21) by dropping a factor of i.
Because of the separation of (time) scalest −1 ω, where ω is a typical frequency we are interested in, the correlation functions change much faster as functions of relative time ∆t than of central timet. Thus, ωρ is approximatelyρ and similarly, the statistical correlation function ω 2 F becomesF . 3 Having introduced the relevant correlation functions, we provide a short summary of expressions that are derived within the HTL formalism at LO [8,46]. Although developed primarily for thermal equilibrium, the HTL formalism can also be applied to systems out of equilibrium [4,47].
An important quantity is the polarization tensor whose transverse and longitudinal components are functions of the ratio x = ω/p. For isotropic non-equilibrium systems with a scale separation between m and the hard scale Λ, the polarization tensors are given by with the Legendre function of the second kind At strict leading order, the mass m can be computed by the LO version of the HTL expression (12). While the polarization tensor is gauge invariant at LO, the exact form of the retarded propagator G HTL R depends on the chosen gauge. For the temporal gauge A 0 = 0 as here employed, its transverse and longitudinal components read The components of the spectral function are obtained by which is consistent with our formerly stated relations in (19) and (21). We note that in the literature the longitudinally polarized spectral function is usually defined asρ to compensate for the gauge dependent prefactor p 2 /ω 2 in (24). However, we will use ρ HTL L instead ofρ HTL L in comparisons with our data since ρ L is what we obtain in both time and frequency domains numerically.
Quasiparticles emerge as poles of the retarded propagator. Therefore, by finding the roots of the denominator in (24) numerically, one obtains the transverse and longitudinal dispersion relations ω HTL T /L (p). Their low-and high-momentum expansions read Here the plasmon frequency in the HTL framework is given by and is approached for both dispersion relations in the limit of low momenta. At high momenta, the transverse dispersion relation corresponds to a relativistic dispersion with asymptotic mass m while for longitudinal momenta, one essentially has an ultrarelativistic dispersion ω HTL L ≈ p.
The quasipartice peak enters the spectral function as a Delta function δ(ω − ω HTL T /L ) with a prefactor. Due to the imaginary part of the Legendre function (23), the polarization tensor obtains an imaginary part for low frequencies ω 2 ≤ p 2 , which also enters the spectral function and is referred to as the Landau cut. Hence, we can write the spectral function as a sum of the Landau cut region β T,L (ω, p) and the quasiparticle peak The Landau cut region with x = ω/p is given by for the transverse case and by for the longitudinal polarization. The spectral function satisfies the sum ruleṡ The transverse sum rule is equivalent to (18) while the longitudinal sum rule results from ). The HTL framework also predicts that the statistical correlation function F is not independent of the spectral function ρ. Instead, they are connected at soft momenta and frequencies ω, p Λ via with T * given by and the integrals where we generalized g 2 J = m 2 HTL /(2N c ) to include the self-consistent resummation of the frequency denominator f /p → f / m 2 HTL + p 2 , as we did for the mass. Moreover, for systems with large occupation numbers f (t, Λ) 1 one can make the approximation f (f + 1) ≈ f 2 . On the other hand, if f is given by the thermal Bose-Einstein distribution (e ω/T − 1) −1 with ω ≈ p, one obtains T = T * and Eq. (34) becomes the fluctuationdissipation relation in thermal equilibrium for ω Λ. The relation (34) can be turned to an equal-time relation by multiplying the equation by ω 2 and integrating over the frequency, which leads tö F HTL (t, ∆t = 0, p) = T * (t)ρ HTL (t, ∆t = 0, p) . (37) If the distribution function follows a 1/p power law at momenta p Λ, then T * can be understood as the effective temperature of low momenta where f (t, p) ≈ T * /p.
While in a thermal system m and T are constants in time, in a nonequilibrium system the distribution function depends on timet, which for instance implies a slow time dependence of the mass (12). Important assumptions for this are an existing scale separation between the hard scale Λ which dominates the energy density and the asymptotic mass m as well as that relevant loop diagrams are dominated by modes of the order of Λ. Although in general this may pose constraints for the form of the distribution function, for the case of a highly occupied non-Abelian plasma close to its self-similar scaling solution that will be studied in Sec. III these conditions are expected to be satisfied up to higher order effects.

D. Linear response theory
To compute the retarded propagator G R numerically, we study the response of the non-Abelian plasma to a small perturbation with a source j k b (x) = j k b (t, x). Then the plasma field can be split into two partŝ which are written as field operators in the Heisenberg picture. If no source is applied, the response is zero â = ê = 0. Otherwise, it is given by 4 [46] â with the retarded propagator Since we consider spatially homogeneous systems, G R does not depend on the central spatial coordinates (x + x )/2 but only on the relative coordinates ∆x = x − x . In Fourier space, Eq. (39) then reads Using for the source an instant perturbation of a mode p at time t pert one arrives at From this, we would like to compute the retarded propagator G R . Because of the summation over indices, we cannot simply divide this expression by the source term. Moreover, since we consider a linear response to a perturbation, we can set multiple momentum modes simultaneously for the source j 0 . To deduce the retarded propagator, we choose the perturbation to satisfy 5 where the number of polarizations is d λ = 2 for the transverse and d λ = 1 for the longitudinal case and where · j is a classical average over a set of sources. Indeed, similar to the initial conditions in our simulations in Eq. (7), the relations (44) can be achieved by choosing 4 Note that we use a phase convention where G R (x, x ) is realvalued. 5 We are grateful to A. Piñeiro Orioli for sharing with us this method in a private communication. This method has been developed for non-relativistic scalar field theories in Ref. [37].
with a random phase c (8) with the replacement . cl → . j . The summation over polarizations λ in the initialization of the source is chosen to involve only transverse modes (λ = 1, 2) when discussing the transversely projected spectral and retarded correlation functions and only the longitudinal polarization (λ = 3) in the longitudinal case. As usual, we will occasionally omit the subscripts T, L to simplify notation.
The retarded propagator, averaged over adjoint components and polarizations, can then be computed as Similarly, we define the time derivative of the retarded propagator aṡ such that one hasĠ R (t, t , p) = θ(t − t )ρ(t, t , p). Therefore, and due to the relation (19), the retarded propagator is fixed at t → t pert in the same way as the spectral function, which serves as the initial condition for its evolution.

E. Linearized fluctuations
To compute the linear response in the discretized classical-statistical framework, the gauge and chromoelectric fields are split according to where A b k and E k b will from now on be referred to as background fields and a b k and e k b are the linearized fluctuation fields. They transform under a gauge transformation V (x) in the same way as the electric background field, which is E k a (x) → V ab (x)E k b (x) in the adjoint representation. These fluctuations are the expectation values of the fluctuation field operators of Sec. II D, where in the classical approximation the expectation value is calculated as an average of the classical distribution and the formulas for G R andĠ R have to be changed accordingly.
The equations of motion for the fluctuations are derived in Ref. [48] by linearizing the lattice equations of motion (2) and additionally demanding that the Gauss law condition (4) is satisfied in the linearized framework. For an SU(2) theory, the a-field update is given by where we defined the matrix U 0j (t, x) = exp(i dt a s gE j a (t, x)Γ a ). All electric fields on the right hand side are evaluated at (t, x) while the linearized gauge field appears there at (t − dt/2, x). Moreover, we split the linearized fields in parts parallel and transverse to the electric background field in color where electric fields and the source j on the right hand side are at time t while the gauge and link fields are taken at t + dt/2. For parallel transported fields we used the notation . We discuss now the initial conditions of the linearized fields. Before the time of the perturbation t pert , the linearized fields a and e are zero. Using a source term according to (42) the electric field at time t pert becomes while the gauge field stays zero a b k (t pert , p) = 0. Therefore, initializing just the linearized e field is equivalent to perturbing the system with a source j at time t pert . Moreover, it can be easily checked that these initial conditions also satisfy the initial relations for the transverse retarded propagator and thus, also of the spectral function (18). 6 To measure gauge-dependent observables such as momentum space correlations, one needs to fix the gauge. The interpretation in terms of physical degrees of freedom is the clearest in a Coulomb gauge 6 The initial condition forĠ T is satisfied exactly before the Gauss law restoration algorithm is applied. After its application, however, this is not the case any more, especially for longitudinally polarized modes, and the initial conditions have to be set by hand by rescaling the amplitude accordingly for each momentum mode. This is possible for well separated excited momentum modes because the equations are linear for the fluctuations.
This implies that the gauge field is always transversely polarized while the electric field may have longitudinal contributions. For equal-time correlation functions the gauge is fixed always at the time of measurement, as for the distribution function f (t, p). Here, however, we are measuring unequal-time correlators of the fields and the fluctuations. Fixing the gauge separately at each measurement time would mean that the fields at the different times would be taken from different gauge trajectories. As we show in Appendix A, this leads to gauge artefacts in the spectral function outside of the quasiparticle peak. In order to avoid this effect, we only fix the Coulomb gauge condition at the time t pert when the fluctuation is introduced, but not at later times. Thus the system gradually shifts away from the gauge condition (54). However, the timescale of this deviation (which is of ordert) is expected to be long compared to the relevant time separations ∆t and the effect should therefore not be large. A more thorough investigation of the gauge dependence of our results is left for a future study. The gauge fixing is efficiently done with a Fourier-accelerated algorithm [49].
Because of (53), a sensible source j has to also satisfy the Gauss law, whose violation in the linearized framework of Ref. [48] is defined by and is preserved by the linearized equations of motion. Therefore, to restore the Gauss law after the perturbation at t pert , we employ the Gauss law restoration algorithm of Ref. [45] for the Gauss violation c a only at time t pert . The algorithm proceeds by iterating the transformation e k a (x) → e k a (x) + γD F,ab k c b (x), with the corresponding forward derivativeD F,ab k and the same parameter γ as in [45]. Our Gauss law violation is typically reduced to 10 −9 − 10 −12 in suitable units of Q, which is close to machine precision. Our results are observed to be insensitive to the chosen Gauss law precision.

F. Numerical setup
In the following sections we show our numerical results. Our standard choice of parameters is n 0 = 0.2 and Qt pert = 1500, while our time extent for the Fourier transform is typically Q∆t max = 200, see also Sec. II D.
If not stated otherwise, we employ a 256 3 lattice with lattice spacing Qa s = 0.7, averaged over 5 simulations. To show that our results are insensitive to changes of volume and lattice spacing, in some figures we compare to a smaller 192 3 lattice with the finer spacing Qa s = 0.47, averaged over 10 simulations. For the time step we use the ratio dt/a s = 0.05, while we have checked that smaller ratios down to dt/a s = 0.01 do not change the results. Unless stated otherwise, we initialize only transverse modes.
Using the method described in Sec. II D, we initialize the source at multiple momenta p. To gain more statistics, we additionally bin our correlation functions linearly in momentum p = |p B |, with the bin size being 1/8 of the smallest momentum on the lattice p min = 2/a s sin(π/N s ). We also checked that smaller bin fractions do not change the results.
We initialize all modes within a momentum bin, which amounts in initializing a thin spherical shell with radius p in momentum space. Moreover, we do not initialize all momentum modes of a linearized fluctuation field. The reason for this is that, especially for low momentum modes, the evolution of G R becomes more noisy when neighboring momentum bins are also excited, since the initialized modes typically disperse and interfere with noise from dispersing neighboring modes if the latter were also excited. Achieving accurate results would require more statistics in that case. Therefore, we initialize single momentum bins separated by at least 0.15 Q, thus reducing noise from neighboring points considerably. As noted above, our simulation results are averaged over 5 − 10 realizations, where each realization has an independent background and source. The error bars correspond to the standard error of the mean. We will explicitly note where data without further averages is taken.
Finally, to further reduce computational costs, we initialize multiple linearized fluctuations independently, exciting different modes, but with the same background field. In this way, we can extract the entire spectrum from a single simulation by initializing sufficiently many linearized fluctuation fields that most momentum bins become excited.

A. Self-similar attractor in non-Abelian plasmas
We start the discussion of our numerical results by presenting some properties of the well-studied attractor in non-Abelian systems far from equilibrium. Starting the (background) non-Abelian plasma at Qt = 0 from the initial conditions of Eq. (9), the system quickly approaches a nonthermal fixed point, where the distribution function follows a self-similar evolution with exponents α = −4/7 and β = −1/7. The selfsimilar behavior is shown in Fig. 2. The occupation numbers rescaled by t −α f are plotted as functions of the rescaled momentum t β p. Since these curves for different times fall on top of each other, this shows that the system follows a self-similar evolution, and the stationary curve corresponds to the scaling function f S (p). For the simulated time in classical simulations, one observes f S (p) ∼ (p/Q) −κ with an exponent κ ≈ 1.3 [28,30] while a comparison with kinetic theory predicts that it should eventually approach κ = 1 at late times [34].  The physical interpretation of the self-similar evolution is a direct energy cascade, where energy density ∼ d 3 p ω(p)f (t, p) ∼ Λ 4 f (Λ) = const is conserved and transported to higher momenta as Λ/Q ∼ (Qt) 1/7 , with Λ being the hard momentum scale that dominates the energy density. Indeed, the evolution in this time regime only depends on the energy density [34], while details of the initial conditions are washed away by the transient evolution to the nonthermal fixed point [30,50]. Since our definition of Q relies on the energy density, all quantities should become independent of the initial conditions at late times Qt.
This is demonstrated at the example of the HTL mass of (12). Because of the parametric relation m 2 HTL ∼ d 3 p f (t, p)/p, the self-similar evolution of f (t, p) should lead to m HTL /Q ∼ (Qt) −1/7 . We show its evolution in Fig. 3 for different initial amplitudes n 0 . One observes that the mass indeed becomes proportional to the expected power law behavior for all amplitudes at sufficiently late times. When expressed in units of Q, the masses for different n 0 are seen to even collapse to a single curve at late times, which signals that the dynamics becomes insensitive to the initial conditions. This is already the case at Qt pert = 1500 for n 0 = 0.2, which is the set of parameters we will be using in the following. For comparison, we show m 2 HTL /p 2 0 as a function of p 0 t in the inset, where the curves stay clearly apart and do not coincide.
The different power laws of the mass and hard scale lead to an increasing scale separation with time m HTL /Λ ∼ (Qt) −2/7 . Since a scale separation is one of the main assumptions in the HTL formalism that was discussed in Sec. II C, one expects this formalism to be applicable to the considered highly occupied system. In the following, we will show our results and compare them to the respective HTL predictions at LO in m HTL /Λ. In principle, we could also tune m HTL /Λ by extending simulations to late times to measure subleading effects. However, computational costs restrict such an analysis since the required lattices would also need to increase to prevent lattice artifacts.

B. Comparing spectral and statistical correlation functions
We start our discussion of spectral and statistical correlation functions by studying their equal-time correla-tionsF (t, ∆t = 0, p) andρ(t, ∆t = 0, p), the latter of which is fixed by (33). In Fig. 4 we show the transverse and longitudinal correlation functionsF T /L (t, ∆t = 0, p) at times Qt = 250 and 1500. Gray areas show the HTL prediction in (37) with error bands. Instead of following a straight line, as predicted by HTL, the transverse correlation depends on momentum, following an approximate power law T * (p/Q) −0.3 between m HTL and Λ and surpassing the gray band denoting T * at low momenta. This power law is connected to the f ∼ (p/Q) −κ behavior with κ > 0 discussed in Sec. III A while κ = 1 is predicted in HTL for a perfect scale separation.
Similarly, the longitudinal correlation shows deviations from the expected T * ρL (t, ∆t = 0, p) evolution. At momenta below the mass p m HTL , bothF T andF L are enhanced and approximately coincide. At high momenta p Λ they decrease exponentially, which is beyond the HTL prediction since the latter is only expected to hold for p Λ. As time proceeds, the longitudinal correla-tionF L is observed to approach the lower bound of the HTL predicted gray band while the transverse correlation functionF T typically comes closer to the upper bound of the HTL prediction.
The HTL relation between the equal-time correlation functions (37) that we were testing is based on the more general fluctuation-dissipation relation (34) that predicts the same connection T * betweenF andρ as functions of relative time ∆t and of ω. Focusing here on the transverse polarization, we will therefore study the relation betweenF T andρ T in the considered system. The normalization of the spectral function is fixed by the sum rule (33). Similarly,F T (t, ∆t, p)/F T (t, ∆t = 0, p) satisfies the same normalization relation. Their evolution is shown in Fig. 5 as functions of relative time ∆t at time Qt pert = 1500. One observes that they lie on top of each other to good approximation and show damped oscillations. Their corresponding Fourier transforms are presented in Fig. 6. They are seen to also coincide in frequency space, establishing the relation and equivalently as functions of ∆t. This relation corresponds to a generalized fluctuation-dissipation relation for our far-from-equilibrium situation by constraining the connection between the correlation functions not to depend on frequency or on the relative time. Compared to the HTL relation (34) at LO, this translates to allowing T * (t) →F T (t, ∆t = 0, p) to depend on time and momentum, which is consistent with what we observed in Fig. 4. These observations show that the correlation function F T (t, ω, p) in the considered highly occupied system is similar to its LO HTL expectation by satisfying a generalized fluctuation-dissipation relation but shows deviations for its amplitudeF T (t, ∆t = 0, p) that becomes momentum dependent instead of being equal to T * . One of the reasons for this is that m HTL /Λ is not negligibly small, which affects the form of f (t, p) and thus, the form of F T (t, ∆t = 0, p) and the values of m HTL and T * . We also observe deviations for momenta below m HTL that cannot be explained by the HTL expressions at hand even qualitatively and may be influenced by processes at the magnetic scale. Moreover, the spectrum in Fig. 6 includes a narrow peak for each momentum, which follows a Lorentzian form 7 (58) Here A is a normalization constant and ω T and γ T are the transverse dispersion relation and damping rate, respectively, all of which can in general depend on momentum and time. This shape is demonstrated in Fig. 7 at the example of p = 0.7 Q. To the correlation functions shown in Fig. 6 we added ω 2 F T /F T (t, ∆t = 0, p) and ωρ T , which confirm that time derivatives become frequency factors in frequency space despite the residual central time dependence. Most importantly, their shape indeed matches with (58), which establishes the existence of quasiparticles in the considered system. 8 Their dispersion relation and damping rate are discussed in the following subsection. 7 The small deviations from this form at low momenta like p = 0.09 Q can be attributed to the finite time window employed and are expected to vanish with a larger time window for the Fourier transform. 8 Note that because of ω T γ T , the spectral function ρ T also follows a Lorentzian form in frequency space with approximately the same ω T and γ T . We have checked this explicitly.

C. Dispersion relation, damping rate and Landau cut
Proceeding with our discussion of correlation functions at unequal time, we study the transverse dispersion relation ω T of the observed quasiparticle peak and compare it to the expected HTL curve at LO ω HTL T discussed in Sec. II C. While the latter involves the plasmon frequency ω HTL pl and asymptotic mass m HTL given by (29) and (12), in general, the plasmon frequency is defined as the frequency of the zero mode, where ω L (p) is the longitudinal dispersion relation that will be discussed and studied in Sec. III D. Since at p = 0 there is no distinction between transverse and longitudinal polarizations, both dispersion relations have to coincide. Similarly, the asymptotic mass is defined as the mass in the relativistic dispersion relation m 2 + p 2 that is expected to be seen for high momenta p m. One way of expressing this relation is Both definitions are consistent with the Taylor expanded LO HTL dispersion relation in (27). However, they also enable us to measure ω pl and m independently and to compare their values to the HTL predictions. More generally, we can now compare our results to the expected LO curves from HTL. We start with the transverse dispersion relation ω T (p) in the upper panel of Fig. 8. It is deduced fromρ T (t pert , ω, p) by finding the frequency of the maximum of the quasiparticle peak for each momentum, but we have also checked that the dispersion resulting from the maxima of ρ T (t pert , ω, p) is consistent within uncertainties. One observes that the extracted dispersion relation does not show sensitivity to the different discretization parameters employed, which indicates that this observable is close to its continuum limit. This can even be observed at low momenta. For instance, the smallest momentum of the Qa s = 0.47 system is located at p ≈ 0.48 m HTL but it is barely visible because the corresponding point of the Qa s = 0.7 data lies almost on top of it.
For comparison, we included a relativistic dispersion relation with a fitted mass value m rel = 0.132 Q into the upper panel of Fig. 8 . We also show there the numerically computed transverse HTL dispersion relation ω HTL T , where we use the mass parameter m HTL = 0.149 Q of Eq. (13), which comes from a computation within the HTL formalism with the distribution function f EE while other definitions lead to deviating values. Interestingly, both functional forms provide good overall descriptions of the Upper: Transverse dispersion relation ωT (p) deduced fromρT (tpert, ω, p) by finding the location of the maximum of the quasiparticle peak for each momentum in a logarithmic plot. The dashed black curve depicts a fit to a relativistic dispersion ω rel T while the red solid curve shows the expected dispersion relation from HTL calculations at LO with HTL mass mHTL = 0.149 Q. For comparison, we also show the ultra-relativistic dispersion relation ωT = p as a gray dashed curve. Lower: The expression ω 2 T − p 2 is shown for the same data in a linear plot. We also included the zero mode frequency ω(0) from the Qas = 0.7 and Qas = 0.47 discretizations of the longitudinally polarized systems in Fig. 12 by an orange triangle and purple circle, respectively. data. For the HTL curve, we get the same value for m HTL when fitting ω HTL T to our data. This justifies a posteriori our choice of the mass value for m HTL in our previous and following figures where data is compared to HTL predicted curves.
On the other hand, there are discrepancies at low and high momenta. They are better visible in the lower panel of Fig. 8 where we show ω 2 T − p 2 as a function of momentum for the curves of the upper panel. While the data points from this combination are quite noisy, one observes some systematic behavior. If the transverse dispersion relation followed a simple relativistic form ω rel T , the data points would be constant and equal to m rel . Instead, one observes that apart from the zero mode this combination steadily grows overshooting this mass value at high momenta and being smaller at low momenta. This is qualitatively similar to the behavior of the HTL curve but the data points are systematically lower at high momenta and higher at low momenta. This systematics can be also observed in the upper panel of Fig. 8 when compared to the HTL curve, i.e., frequencies are shifted to slightly larger values at low momenta. Moreover, we will see in Sec. III D that the longitudinal dispersion has the same systematic behavior when compared to the corresponding HTL curve. After this qualitative discussion of the functional form of the extracted transverse dispersion relation, we can deduce the values for the plasmon frequency and the asymptotic mass from our data. We measure the plasmon frequency directly by finding the maximum in frequency space of the spectral function at p = 0 ω fit pl (Qt = 1500)/Q = 0.132 ± 0.002.
The corresponding points at p = 0 are depicted in the lower panel of Fig. 8. They correspond to different discretizations, which lead to the same value within the given uncertainty. This value is larger than the predicted value ω HTL pl = 0.122 Q but agrees with the value from the relativistic dispersion fit ω rel pl = m rel . The asymptotic mass m is deduced in our simulations by fitting ω HTL T (p) to the observed dispersion relation ω T (p) for high momenta p min ≤ p ≤ 1 Q = 6.7 m HTL . In practice, we first compute the HTL dispersion relation at various momenta p, interpolate the solution and fit the interpolation function to the data. Varying the lower momentum between p min = 2 m HTL − 4 m HTL and considering the maxima of both ρ T andρ T provides an estimate for the systematic error while the (fit) value is the typical error of the fit m fit /Q = 0.138 ± 0.002 (sys) ± 0.0015 (fit).
The fitted value for the asymptotic mass is smaller than the value used for m HTL , which is consistent with the lower panel of Fig. 8 because the ω 2 T − p 2 data points lie below the corresponding HTL curve. According to (27), ω rel T approximates the HTL curve ω HTL T at large momenta. However, the value from the relativistic dispersion fit m rel is lower than the extracted value m fit , which is also consistent with the lower panel of Fig. 8.
With this, we can compute the relation between the plasmon frequency and the mass. At LO in the HTL framework one has ω HTL pl /m HTL = 2/3 ≈ 0.8165, while for the relativistic dispersion one has ω rel pl /m rel = 1. We find ω fit pl m fit (Qt = 1500) = 0.957 ± 0.028 , which is even closer to 1 than to the HTL expected 2/3. However, we emphasize that the extracted transverse dispersion agrees well with both functional forms and shows Transverse damping rate γT (p) deduced froṁ ρT (tpert, ∆t, p) by fitting to the form of a damped oscillator (65) within a time range of Q∆t ≤ 70. The error bars are calculated as a sum of the error of the mean and the average fitting error. The expected HTL value γHTL(p = 0) is also included. For the data points at p = 0 we fitted a damped oscillator to already averaged data within Q∆t ≤ 100, therefore not extracting error bars.
deviations from both forms as well. On the other hand, as we will see, we will observe the HTL predictions of a Landau cut region and a different dispersion relation for longitudinally polarized modes. This, together with the agreement with the transverse dispersion relation, shows that the HTL formalism provides a good overall description of our data, while systematic deviations can be taken as an indication of effects beyond this formalism at LO.
Let us now proceed with the discussion of the damping rate of transverse quasiparticles. Within the HTL framework, poles of the transverse retarded propagator at leading order (24) correspond to Delta functions δ(ω − ω HTL T (p)) in ρ T . However, realistic quasiparticle peaks involve a finite width, which corresponds to a non-vanishing damping rate γ T (p), as was observed in Sec. III B. This is an effect of subleading order in the HTL framework and may contain non-perturbative contributions. So far, the damping rate within HTL has only been calculated at p = 0 [41]. Our employed numerical framework provides a unique opportunity to access the momentum dependence of this quantity out of equilibrium.
Instead of working in Fourier space, one can measure γ T (p), as well as ω T (p), also directly in the time domain ∆t. A Lorentzian peak Fourier transforms to a damped oscillator, which reads at the example ofρ T . Using this fitting curve directly in the time domain, we have checked that the dispersion relation ω T (p) extracted in this way coincides with the data shown in Fig. 8. However, it turns out that this procedure provides more accurate values for the damping rate than the corresponding measurement in frequency space, where in practice, the finite time window for the Fourier transform leads to deviations from Lorentzian peaks especially at low momenta, introducing additional uncertainties on the damping rates.
The damping rate extracted using (65) is shown in Fig. 9 for different discretizations that lie on top of each other within uncertainties. The data points are obtained by averaging over the damping rate of different simulations while the error bars are a sum of the error of the mean and the average fitting error. For low momenta p 0.15 Q (i.e., p m HTL ), the damping rate is observed to increase before it eventually flattens at higher momenta. We have also included the expected value at p = 0 within the HTL formalism that has been computed in Ref. [41] γ HTL (0) = 6.63538 where we have replaced T → T * in the original formula and estimated the error bars by employing different definitions of the distribution function in (11). Our data is roughly consistent with the HTL predicted value. The extracted damping rate γ T (p) in Fig. 9 is one of the main results of this work. In our considerations, we have only talked about the dominant quasiparticle peak so far. However, the transverse spectral function has a much richer structure, also involving the Landau cut as discussed in Sec. II C. Our results for the spectral function ρ T are shown in Fig. 10 for different momenta, where we also included the expected HTL curves. While the additional ω-factor suppressed the Landau cut region in Fig. 6, the corresponding structure at low frequencies for ω ≤ p is clearly visible in Fig. 10. Moreover, it agrees well with the HTL curves at LO. Note that the only parameter in these curves is the mass m, for which we employ m HTL to be consistent with Fig. 8. Therefore, no free parameter is left and the agreement with the curves indicates that the HTL framework provides a valid description for the Landau cut region even far from equilibrium.
Small discrepancies at the lowest frequencies ω 0.02 Q can be attributed to the finite resolution in frequency space due to limitations of the time window for the Fourier transform. On the other hand, the quasiparticle peak of low momenta like p = 0.05 Q is located at a slightly larger frequency than expected at LO in HTL. This behavior is consistent with our observation of the dispersion relation in Fig. 8.

D. Longitudinal correlations
Apart from the transversely polarized correlation functions, we can also study longitudinal correlation functions at unequal time. Our data is shown for the spectral functionρ L in Fig. 11 in relative time and frequency domains for different momentum modes, where we choose it to satisfy the sum rule (33). In the time domain, one observes damped oscillations similar to what was observed for the transverse case. For larger momenta p 0.15 Q the oscillations decrease very quickly and oscillations are barely visible at later times.
Their Fourier transforms are shown in the right panel of Fig. 11 together with the longitudinal spectral function computed within the HTL formalism at LO that was discussed in Sec. II C. For small momenta p 0.15 Q, one sees clear quasiparticle peaks and a small Landau cut region, which is consistent with the HTL curve. At larger momenta p 0.15 Q, the HTL Landau cut region becomes more pronounced and starts to dominate the sum rule (33) while the residue of the quasiparticle peak gets exponentially suppressed as ∼ exp(−p 2 /m 2 ). At the same time, the longitudinal dispersion relation of the quasiparticle peak gets exponentially close to the light cone ω L ≈ p, where the Landau cut region starts. Therefore, it becomes difficult to distinguish between the Landau cut region and the quasiparticle peak numerically. Moreover, the Landau cut region gets smeared a little around the light cone, such that one sees practically only the Landau cut for p 0.29 Q.
We also show the normalized longitudinal statistical correlation functionF L /F L (t, ∆t = 0, p)ρ HTL L (t, ∆t = 0, p) as orange curves in both relative time and frequency domains. This quantity is noisy but one observes an overall good agreement with the spectral function. Although it seems as for p = 0.12 Q it does not capture the Landau cut region correctly, at larger momenta where the Landau cut region starts to dominate the spectrum, it agrees quite well with the spectral function. This confirms the existence of an approximate generalized fluctuationdissipation relation similar to the transverse case in (57) also for the longitudinal correlations. The connection be-tweenF L andρ L is provided byF L (t, ∆t = 0, p) that was depicted in Fig. 4.
The longitudinal dispersion relation ω L (p) can be extracted from the spectral function by, for instance, finding its maximum after subtracting the Landau cut. Since numerically, the Landau cut region appears as smeared around the light cone, while the analytic function has a steep increase there, the subtraction is not precise in the region where the Landau cut starts dominating the spectrum. Therefore, the longitudinal dispersion extracted   Fig. 9 for the latter), extracted by fitting to a damped oscillator (65). Error bars consist of a sum of statistical and fitting error. For the data points at p = 0 we fitted a damped oscillator to already averaged data within Q∆t ≤ 100, thus not obtaining error bars.
this way is not more than an estimate for the true longitudinal dispersion relation for momenta p m. Our results for this are shown in the upper panel of Fig. 12 for two different discretizations. The transverse dispersion relation ω T (p) from Fig. 8 as well as the longitudinal and transverse dispersion relations from the HTL formalism at LO are added to the figure. One observes that the longitudinal dispersion estimate is quite distinct from the transverse one and exhibits very similar momentum dependence as the corresponding HTL curve. As expected, there are differences to the latter around p ∼ m HTL , which may have resulted from the subtraction of the analytical Landau curve instead of the nu- merical one. On the other hand, at low momenta one sees a similar deviation from the HTL curve at LO as for the transverse distribution in Fig. 8, where the extracted frequencies are systematically larger. This behavior is seen for both discretizations. Moreover, both dispersion relations come close, smoothly approaching ω pl with the value (62), which is larger than expected from LO HTL, as discussed in Sec. III C. Hence, ω L (p) confirms the observations in that section.
For momenta where the Landau cut does not dominate the spectral function, one can extract the longitudinal damping rate γ L (p) by fitting to a damped oscillator (65). We show the extracted data points in the lower panel of Fig. 12, where we also include the transverse damping rate from Fig. 9 for comparison. For larger momenta p m HTL , the Landau cut becomes important, eventually dominating the ∆t evolution ofρ, and the fitting procedure provides information about the Landau cut rather than providing γ L (p). Therefore, the lower panel of Fig. 12 describes only low momenta p m HTL correctly. For these, one observes that the damping rate does not depend on the polarization within uncertainties. A similar observation was made for the equal-time correlation functionF (t, ∆t = 0, p) in Fig. 4.

E. Variation of initial parameters
We finally check how our results depend on a variation of the initial amplitude n 0 and time t pert . We will discuss our findings at the example of the transverse dispersion relation and damping rate.
Our results on the dispersion relation are shown in Fig. 13 for simulations with n 0 = 3.2 at three different times Qt pert = 400, 750 and 1500 and for our n 0 = 0.2 curve at time Qt pert = 1500 from Fig. 8. All simulations have been performed on 256 3 lattices with Qa s = 0.7 and the n 0 = 3.2 simulations have not been ensemble averaged. The curves are plotted in units of m HTL that is computed as in Sec. II A with f EE (t, p). While the mass is time-dependent as shown in Fig. 3, the functional form of the dispersion relation stays the same since the curves are seen to lie on top of each other. This is predicted by the HTL formalism, where the mass is the only parameter of the dispersion relation, and hence, expressing frequencies and momenta in terms of the mass should lead to the same form. We note, however, that the agreement between the curves implies that all curves share similar deviations from the HTL prediction at LO as discussed in Sec. III C. Whether these deviations are even the same and, thus, time-independent cannot be answered accurately with the present data and is left for future work. The transverse damping rate γ T (p) is shown in Fig. 14 for different times, rescaled by (t/t ref ) 3/7 with Qt ref = 1500 as a function of rescaled momentum (t/t ref ) 1/7 p. The exponent 1/7 corresponds to the evolution of the mass 1/m HTL , while the exponent 3/7 in the damping rate reflects the evolution of the effective temperature which follows from the self-similar evolution of the distribution function. We checked explicitly that g 2 T * (t) follows this power law in our simulations. Since the rescaled data coincide while the original one does not as seen in the inset, this indicates that the damping rate is proportional to γ T (p) ∝ g 2 T * for all momenta and not only for the zero mode as given by Eq. (66). However, note that these are single-run simulations and therefore, more statistics is needed to confirm this observation beyond doubt. A precision measurement of the damping rate may also show possible influence of the magnetic scale, which is of the same subleading order g 2 Q.

IV. CONCLUSION
In this work, we have studied the spectral and statistical correlation functions ρ and F at unequal times of highly occupied classical Yang-Mills theory at the selfsimilar regime to get further insight into the dynamics. To compute the spectral function, we have used linear response theory with our recently developed formalism to simulate linearized fluctuations on top of a classical background. With these techniques, we were able to access information about the highly populated system that had not been accessible previously by only considering equaltime correlation functions. Our results were compared to the HTL formalism at LO.
We showed that the functional forms of the transverse and longitudinal spectral functions ρ T and ρ L agree well with the HTL-predicted curves, including the Landau cut region of low frequencies ω 2 ≤ p 2 . For larger frequencies, both ρ T and ρ L involve a quasiparticle peak with a dispersion relation and damping rate. The dispersion relations agree well with the corresponding HTL curves but also with a relativistic dispersion in the transverse case while the damping rate is observed to be proportional to an effective temperature T * and to roughly agree with the predicted value in [41] at p = 0. While transverse and longitudinal dispersion relations follow distinct functional forms, the corresponding damping rates are observed to agree within error bars for momenta below the mass.
For the dispersion relations, we found some deviations at low and high momenta. Extracting the plasmon frequency ω pl and asymptotic mass m from our data, we were able to quantify some of these deviations. Instead of the HTL prediction ω pl /m = 2/3, we found a value close to 1 for the ratio. Moreover, we observed that the spectral and statistical correlation functions are not connected by the effective temperature T * as predicted by the HTL formalism but by a momentum dependent function that is of the same order of magnitude for momenta m p Λ.
Many general properties of the HTL correlation functions at LO in m/Λ have been observed in our simulations. The deviations from the HTL expressions may have resulted from higher order corrections since for the studied times, we did not have a perfect scale separation between the mass and the hard scale. Moreover, we were able to measure the transverse and longitudinal damping rates for the first time, extending the result from [41] to finite momenta. In this sense, our simulations can also be regarded as a non-perturbative way of measuring observables that are hard to obtain analytically within the HTL formalism. We plan to use this numerical method also for similar cases, where HTL quantities are hard to access.
Moreover, in the future we aim to extend our studies to the case of a background field with a very anisotropic momentum distribution. Here the fluctuations are known to be unstable, i.e., one has to map out both the real and imaginary parts of the dispersion relation, and with different angles of orientation between the momentum vector and the anisotropy. The advantage of the formulation in terms of linearized fluctuations is that this method does not rely on a scale separation between the hard and soft modes, although we have here operated in a regime where such a separation exists. Thus, for the isotropic case, the calculation can be extended to small values of Qt before the onset of the self-similar regime. More importantly for the phenomenological context, the calculation can be performed in a longitudinally expanding coordinate system and extended all the way to τ = 0. This will enable one to follow the growth of the quantum fluctuations in a realistic geometry corresponding to a relativistic heavy ion collision. We intend to pursue these avenues in future work. therefore follows a Lorentzian form and matches with our usual computation of the spectral function. This implies that both the dispersion relation as well as the damping rate of the peaks inF are close to those measured above. The amplitude modulations and thus, the side peaks of the spectral function are quite sensitive to the lattice spacing and the volume, showing clearer side peaks and stronger modulations with increasing volume.