Broad excitations in a 2+1D overoccupied gluon plasma

Motivated by the initial stages of high-energy heavy-ion collisions, we study excitations of far-from-equilibrium 2+1 dimensional gauge theories using classical-statistical lattice simulations. We evolve field perturbations over a strongly overoccupied background undergoing self-similar evolution. While in 3+1D the excitations are described by hard-thermal loop theory, their structure in 2+1D is nontrivial and nonperturbative. These nonperturbative interactions lead to broad excitation peaks in spectral and statistical correlation functions. Their width is comparable to the frequency of soft excitations, demonstrating the absence of soft quasiparticles in these theories. Our results also suggest that excitations at higher momenta are sufficiently long-lived, such that an effective kinetic theory description for 2+1 dimensional Glasma-like systems may exist, but its collision kernel must be nonperturbatively determined.


Introduction
In the weak coupling framework, the initial stage of a relativistic heavy-ion collision in the high-energy limit is dominated by nonperturbatively strong boost invariant gluon fields [1]. At the characteristic momentum scale Q s , these strong fields have nonperturbatively large occupation numbers f ∼ 1/g 2 and are thus effectively classical [2,3]. The success of a hydrodynamic description of the later stages of the evolution indicates that the matter evolves towards sufficient isotropy quickly. One of the major open questions in the field is to understand in a controlled theoretical framework the evolution from a classical fielddominated system with approximately boost-invariant fields at τ 1/Q s to a state close to local thermal equilibrium [3].
Several studies in the past have addressed the evolution of the initial boost-invariant field during the timescale τ ∼ 1/Q s , which can be understood in terms of 2+1-dimensional lattice simulations [4][5][6][7] or analytical calculations [2,[8][9][10]. These approaches have been combined with hydrodynamic evolution in the later stages, producing a phenomenologically successful descriptions of heavy-ion collisions [11][12][13]. However, these calculations have not been based on a detailed microscopical description of the evolution from boost invariant classical fields to hydrodynamics.
In practice, the initial boost invariance is broken by different mechanisms. At finite collision energy it is broken by the finite thickness of the nucleus [14][15][16]. Another source of violation are rapidity dependent unstable quantum fluctuations [17][18][19][20][21]. After the boost invariance has been broken and the system has become more dilute with an occupation number 1 f 1/g 2 , the classical fields can also be described using kinetic theory [22][23][24][25]. This stage has been studied extensively using real time lattice simulations [25][26][27][28][29][30][31]. More recent studies, which match kinetic theory to hydrodynamics [32][33][34], seem to favor a bottom-up type isotropization scenario [35]. However, it is unclear if the kinetic theory approach can be used all the way back to the purely 2+1-dimensional (2+1D) initial stage. This is one of the questions we aim to address in this paper, which is a continuation of our earlier study of equal-time correlation functions in the same systems in Ref. [36].
Formulating a kinetic theory in 2+1D faces certain problematic issues [36]. The kinetic theory relies on a separation of modes, into soft modes described by classical field theory and hard momentum modes described as particles that undergo scattering. In order to find leading order expressions of the medium-modified scattering cross-sections of the hard particles, the dynamics of the soft sector need to be solved. In three dimensions the dominant medium modification to the soft fields arises from the interaction with the hard modes, and this interaction can be analytically solved giving rise to the Hard-Loop (HTL) theory [37,38] and the Debye mass scale m D . However, in two spatial dimensions, the lower dimensionality puts more weight on infrared modes in momentum integrals, and consequently, the soft and hard scales contribute equally to the Debye screening mass. As a consequence of their self-interactions, already modes at the soft scale m D become nonperturbative, bearing resemblance to the magnetic scale in three dimensions. This means that it is not obvious whether the soft modes can be described using HTL theory. This makes it difficult to find leading order accurate matrix elements for kinetic theory in 2+1D analytically.
Nevertheless, our recent studies in 2+1D classical lattice gauge theory [36] indicated that not only isotropic 3+1D but also 2+1D systems exhibit a self-similar attractor behavior, for which we extracted the scaling exponents. We also found that the scaling exponents can be understood using a simple kinetic theory analysis. This suggests that 2+1D systems involve quasiparticle excitations, at least at high momenta.
Inspired by these results, our main motivation in this paper is to study the excitation spectrum of a gluonic plasma in extreme anisotropy, such that it is effectively in 2+1 dimensions. For this reason, we work in a theoretically clean environment, which is not complicated by additional phenomenological ingredients, and we will study theories with and without a scalar field in the adjoint representation. Especially the former case resembles the boost invariant initial state of heavy-ion collisions at high energy, as the scalar field arises by dimensional reduction from 3+1 dimensional pure gauge theory. In the future, it will be interesting to study excitation spectra of expanding systems and the influence of different initial conditions that are used to model heavy-ion collisions, while here we are looking at universal properties of the spectrum that are independent of the details of the initial conditions.
In practice, the excitation spectrum is studied using the linear response framework built in [39] and first applied in [40] to an isotropic 3+1D gluonic system. A similar classicalstatistical linear response framework has been recently used in scalar theories at self-similar attractors [41,42] and extends classical-statistical simulations in thermal equilibrium that use a fluctuation-dissipation relation explicitly [43,44]. Our work is a natural extension of our previous study of the excitation spectrum of isotropic 3+1D gluodynamics at a classical self-similar attractor [40]. There we have observed narrow quasiparticle excitation peaks in the spectrum for all momenta and a generalized fluctuation-dissipation relation between spectral and statistical correlation functions.
While we also find a generalized fluctuation-dissipation relation, our main result is that the gluonic correlation functions in 2+1D systems exhibit a broad excitation peak of non-Lorentzian shape, in contrast to 3+1D. In fact, we show that the peak width is of the order of and scales with the mass. We separately confirm in Appendix B that this happens also in classical thermal equilibrium. This indicates that gluonic excitations for momenta below the mass are too short-lived to form quasiparticles. Moreover, expressions obtained from HTL theory mostly do not provide a good description of our simulation results. This is in line with the breakdown of the HTL theory at the soft scale. The interaction between soft excitations is genuinely nonperturbative in 2+1D, and the dynamical description has to take this into account. Our results also suggest that gluonic and scalar excitations at higher momenta are sufficiently long-lived, such that we may have an effective kinetic description also in the strongly anisotropic systems which arise in the early stages of ultrarelativistic heavy ion collisions, albeit with nonperturbatively detemined collision kernels.
We start in Sec. 2 by introducing the considered theories and relevant correlation functions, after which we describe our numerical method and introduce the main HTL expressions. Our numerical results for non-equilibrium 2+1D systems are shown in Sec. 3. We conclude in Sec. 4. The analytical forms for the HTL functions in 2+1D are derived in Appendix A, and a numerical study in classical thermal equilibrium is discussed in Appendix B.

3+1D, 2+1D and Glasma-like 2+1D theories
We consider non-Abelian SU(N c ) gauge theories with N c = 2 in d spatial dimensions. Their classical action reads where repeated color indices a = 1, ..., N 2 c − 1 and Lorentz indices µ, ν = 0, ..., d imply summation over them. Using the generators Γ a of the su(N c ) algebra, the gauge field can be written as a fundamental representation matrix A µ = A a µ Γ a . We consider the following theories: • '3D' or '3+1D': in d = 3 spatial dimensions, such that • '2D+sc' or 'Glasma-like 2+1D': originally in d = 3 spatial dimensions where no field depends on the coordinate x 3 . This results in a 2+1D theory with the classical action with an adjoint scalar field φ a ≡ A a 3 having the action with summation over µ = 0, 1, 2 and with the covariant derivative D ab µ = δ ab ∂ µ − gf abc A c µ . The length in the x 3 direction L 3 drops out of the classical dynamics. We have also factored out a constant momentum scale Q and rescaled the gauge coupling and all fields as g Q 1/2 → g, A Q −1/2 → A, . . . , such that the momentum dimensions of the action, the coupling constant and the fields match the ones for 2+1D. This The momentum scale Q used in the rescaling is in principle arbitary, but we will use a value related to the conserved energy density, as discussed in more detail below.
We will refer to this theory as Glasma-like 2+1D, since the Glasma occurring at the initial stages at ultrarelativistic heavy ion collisions is also invariant in one spatial dimension (rapidity) [2,6,8]. Unlike the expanding Glasma usually employed in models of the initial collision dynamics, the 2+1D theory we consider here is defined in a nonexpanding coordinate system, and our initial condition corresponds to a positive (but small) longitudinal pressure P L > 0. Also in the expanding Glasma the expansion becomes less important at later times, as reflected e.g. in the fact that P L becomes positive at τ 1/Q s . Thus our Glasma-like system could be interpreted as a nonexpanding model of this later Glasma state.

Spectral and statistical correlation functions
Here we will give a brief overview of the spectral and statistical correlation functions measured in this work. For a more comprehensive introduction and description of the methods used we refer the reader to [40].
The statistical correlation function is defined as the anticommutator of two field operators with spatial components j, k = 1, . . . , d, the dimension of the adjoint representation d A = N 2 c − 1 and spacetime coordinates x ≡ (t, x). Because of E j (x) = ∂ t A j (x), the correlators EE and AA are related via time derivatives. In analogy, the statistical correlation function for scalar components of the Glasma-like theory are defined as φφ (x, The scalar fields are related to the fields of the original theory via φ ≡ A 3 and π ≡ E 3 . In the classical limit the statistical correlation functions are easy to measure, since the anticommutator of Heisenberg field operators reduces to a product of two classical fields.
The spectral function, on the other hand, is defined as the commutator of two field operators and analogously for the scalar components of the Glasma-like theory. We will mostly study its time derivativeρ, which we will refer to as the "dotted spectral function". In the classical limit the commutator corresponds to the Dirac bracket, which generalizes the Poisson bracket for systems with constraints (such as the Gauss' law). Therefore, a direct measurement of the spectral function is a quite involved task. However, the spectral function is intimately related to the retarded propagator The latter can be extracted numerically by using our linear response framework [39,40]. There we introduce a small instantaneous source j for different momentum modes on top of the classical fields A(x), E(x). The source generates small response fields a(x), e(x). These follow linearized equations of motion constructed to satisfy the Gauss law constraint exactly. The spectral function ρ and the dotted spectral functionρ are computed as correlation functions of j with a and e, respectively. In (2.8) we implied a spatial Fourier transform with respect to x − x . Apart from working in momentum space, it is beneficial to compute the correlation functions in frequency space. For that, a Fourier transform with respect to the relative time ∆t = t − t for fixed central timet = (t + t )/2 has to be performed, respecting the even parity of EE (t, ∆t, p) andρ(t, ∆t, p) with respect to ∆t. In practice, as discussed in Refs. [40,45], we approximate this by a numerically more efficient transform where the lower limit t is kept fixed, and one Fourier transforms with respect to an upper limit t + ∆t: and analogously forρ. We justify this approximation by the observation that our systems at sufficiently late times depend on t only weakly within a relative time window of the order of the inverse plasmon mass ∼ 1/ω pl , which is the relevant timescale for ∆t. If not stated otherwise, the resulting curves are made smoother by employing standard signal processing techniques, as in Ref. [42]. In this work, we will not Fourier transform the correlation functions AA and ρ for 2+1D theories directly, because these functions turn out to oscillate around a non-zero value in the time domain ∆t. In stead, we obtain the frequency space spectral function by Fourier-transforming the time derivativeρ(t, ∆t, p), and subsequently dividing by ω. Note that ρ(t, ∆t, p) and ρ(t, ω, p) are odd functions in ∆t and ω respectively, and are related by a Fourier sine transform in stead of the cosine transform (2.9). Thus we have explicitly ρ(t, ω=0, p) = 0. However, as we will see in Sec. 3.5, the oscillations around a non-zero value at large ∆t lead to the actual spectral function not approaching this limit smoothly when ω → 0.
We will be interested in these correlation functions for transverse and longitudinal polarizations of the gluon field. In 2 spatial dimensions, they can be computed in momentum space by projections on transverse or longitudinal vectors, respectively, as for The (single-particle) distribution function is an important quantity that we can use to characterize the dynamical state of the system. We will use the same definitions as in our previous publications [36,40,45] The employed temporal gauge leaves room for gauge transformations that only depend on the spatial coordinates. Since all the correlation functions discussed in this subsection are not manifestly gauge-invariant observables, we remove the residual gauge freedom for equal time measurements by fixing to Coulomb-type gauge ∂ j A j = 0 t at the time t of the measurement, as often employed in classical-statistical gauge simulations [25][26][27][28][29][30]. In this physical gauge, gauge fields are always transversely polarized while electric fields may include longitudinal terms. For unequal-time correlators we also impose the same condition at the initial time t for the measurement. During the time evolution, the system then gradually shifts away from the gauge condition so that it is not exactly satisfied at the second measurement time t + ∆t. However, as argued in [40], this effect is expected to be small for sufficiently short relative times ∆t t, which we also assumed in order to perform the Fourier transform (2.9) in a finite ∆t interval.

Self-similarity and nonthermal fixed points
Many highly occupied systems exhibit non-thermal fixed points (also referred to as universal classical attractors). These have been seen in classical non-Abelian gauge theories [25,27,28,36,46], scalar theories [47][48][49][50][51], longitudinally expanding systems [30,52,53] and recently also in ultra-cold atom experiments [54][55][56]. A characteristic feature of these attractors is that the system forgets the details of its initial conditions and can usually be understood using a very simple set of scaling functions and power laws, which tremendously simplifies the theoretical description.
In our previous work [36], we established that pure gauge theories also in 2+1D exhibit self-similarity, as defined by and extracted the universal scaling exponents We observed there that the scaling exponents in the 2+1D systems can be understood using relatively simple kinetic theory arguments, even if a kinetic theory description is not expected to work quantitatively. This result was one of the main motivations for the present paper, since it hints that a quasiparticle description of 2+1D theories at high momenta might be possible.
To understand the physical meaning of Eqs. (2.12) and (2.13), we define the timedependent hard scale Λ(t) as the momentum scale that contributes the most to the perturbative estimate of the energy density ε ∼ d 2 p p f (t, p). Then the scaling relation (2.12) and the scaling exponents imply where we included the expected scaling of the soft scale represented by the Debye mass from Eq. (A.3). Thus, the hard scale grows and the soft scale decreases with time, while energy density is conserved. Throughout this work all the measurements are performed at sufficiently large values of time t to be in the self-similar regime.

Initial conditions
For the 2+1D systems, we use the same initial condition as in [36]. We consider weakly coupled g 2 /Q 1 but highly occupied f 1 systems. The initial single-particle distribution function f (t, p) at the initial time t = 0 both for gauge and scalar excitations (i.e., f E = f π = f ) is given by Here Q is a gauge invariant momentum scale, which is more precisely defined by .
We consider systems in a fixed size box, and thus the energy density is a conserved quantity.
Since on the classical level it is possible to carry out field simulations without making explicit reference to the precise value of the coupling g, the quantities we have access to are the energy density scaled with the coupling g 2 ε, which has the momentum dimension [g 2 ε] = 4, and g 2 f . Unless stated otherwise, we will use the initial occupation number n 0 = 0.1, for which p 0 = Q for our chosen definition as detailed below. We have previously shown [36] that the form of the initial conditions is unimportant since the systems will approach an attractor solution that only depends on Q in the sense of Eq. (2.12). The number of non-longitudinal polarizations in (2.15) is d pol , with d pol = 1 for the 2+1D and d pol = 2 for the Glasma-like 2+1D theories. The constant C is taken as C = 20 √ 2π ≈ 50, which we merely choose for convenience such that for n 0 = 0.1 one has p 0 = Q.
In 2+1D the coupling constant g is dimensionful: if one keeps the dimensionless combination g 2 /Q constant, one observes that (2.16) leads to the proportionality Q ∝ 3 √ ε, which is natural for a scale derived from a 2-dimensional energy density. Combining the definition (2.16) with a perturbative estimate ε ≈ d pol (N 2 c − 1) d 2 p/(2π) 2 p f (t = 0, p) for the energy density, one obtains For the 3+1D theory we employ the same isotropic initial condition as in [40] with the distribution function initially given by Here the coupling is dimensionless. In this work we use n 0 = 0.2 for the spatially threedimensional case and define the characteristic momentum scale as Q = 4 √ 5n 0 p 0 ∝ 4 g 2 ε. In our figures, all dimensionful quantities are rescaled by appropriate powers of Q of the respective theory to make them dimensionless, unless a different rescaling prescription is stated explicitly.

HTL expressions
Diagrammatically, the HTL approximation corresponds to an all-order resummation with the kinematic approximation that the external lines are considered to be soft compared to the hard momenta flowing in the internal lines of the diagrams. This corresponds to the non-Abelian generalization of the Vlasov equations, the Wong equations, where the soft modes are considered to be classical fields and the hard modes correspond to classical particles [38]. While we expect the interaction between the soft and hard modes to be described by the HTL theory as is the case in 3+1D, the parametric counting for 2+1D in [36] suggests that in addition to the interaction between soft and hard modes, there is also a leading-order interaction between the soft modes among themselves, which is not captured by the HTL approximation. Therefore, we expect deviations from HTL expressions arising from the nonperturbative soft-soft interactions. We compare our numerical results to the HTL calculations to quantify the deviations.
Here we summarize the main HTL results that are relevant for the comparison with our data. Details are written in Appendix A. The HTL spectral function can be decomposed into a transversely polarized, a longitudinally polarized and, for the Glasma-like theory, a scalar contribution, denoted below by the index α = T , L, or φ, respectively. Each contribution can be further split into a Landau damping part for |ω| < p and a quasiparticle part The Landau damping expressions in 2+1D read with x = ω/p. The Debye mass m D entering the expressions is determined within HTL at leading order in Eq. (A.3) and depends on the distribution function f (t, p). It is connected with the asymptotic mass via m 2 D = m 2 HTL in 2+1D theories and m 2 D = 2m 2 HTL for the 3+1D system. In this work, we compute the asymptotic mass like in [40] in d spatial dimensions as which is a self-consistent generalization of Eq. (A.3) and reduces the dependence on the definition of the distribution function [40]. The extracted and normalized transverse statistical correlation function EE T (t, ω, p)/ EE T (t, ∆t=0, p) is shown for a) the 2 + 1 dimensional theory, b) for the Glasmalike theory and c) for an isotropic 3 + 1 dimensional system. Additionally, we show the scalar correlator ππ (t, ω, p)/ ππ (t, ∆t=0, p) of the Glasma-like system in d). Note that all amplitudes have been cut off at the same fixed value 22/Q corresponding to the red region. For comparison, we added the HTL dispersion relations ω HTL T /φ (p) as black dashed lines and a relativistic dispersion ω rel (p) = ω 2 pl + p 2 as a continuous gray line with the extracted values for ω pl as in the text.
The dispersion relations in 2+1D are where we definedp = p/m HTL . The expressions for the quasiparticle residues Z α (p) are written in Appendix A.

Numerical results
In this section, we present our numerical results for the excitation spectrum of nonequilibrium 2+1D theories by extracting spectral and statistical correlation functions. We briefly discuss some results in (classical) thermal equilibrium in Appendix B.
For our simulations, we use the lattice size 1024 2 and lattice spacing Qa s = 1/8 for the 2+1D system and Qa s = 1/4 on a 512 2 lattice for the Glasma-like 2+1D system. We showed in [36] for equal-time correlation functions that for the considered initial conditions, these discretization parameters are sufficient to avoid lattice artifacts. Our resulting correlation functions will be compared with those from isotropic 3+1D systems computed in [40] with Qa s = 0.7 on a 256 3 lattice, which has been demonstrated not to show any significant lattice artifacts for this choice of parameters.
While respecting the Gauss law constraint to almost machine precision, the linearized fluctuations on top of the inhomogenous, time-dependent, background field do not have conservation laws of quantities like energy, momentum and angular momentum. Thus, there is no guarantee that they would not grow strongly with time, leading to increasing error bands. Indeed, we observe such a growth, more strongly in two than in three spatial dimensions. This limits in practice our ability to extract the spectral function in the time domain. In practice we do not extend our extraction of the spectral function for 2+1D theories to more than ∆t ≤ ∆t max = 80/Q. No such restrictions apply for statistical correlation functions, which allows us to simulate to much later times. While it suffices to use ∆t max = 80/Q also for transverse and longitudinal statistical correlators in 2+1D theories, unless stated otherwise, we employ ∆t max = 200/Q for the scalar statistical correlator of the Glasma-like theory to correctly capture the long-lived scalar excitations at low momenta.
All the results in the figures are shown at Qt = 500 for the 2+1D systems and Qt = 1500 for the 3+1D theory, unless stated otherwise. From our earlier papers [36,40] we see that the systems at these times are well in the universal scaling regime discussed in Sec. 2.3. We will discuss how the correlators depend on time t in Sec. 3.6.
Correlation functions as functions of frequency or relative time for fixed momenta are shown with error bars. These result from statistical averaging of the correlation functions over 320 -400 configurations in frequency space or in the relative time domain, respectively. Figure 1 shows the extracted normalized statistical correlation functions EE α (t, ω, p)/ EE α (t, ∆t=0, p) as functions of frequency and momentum for transverse excitations in 2+1D theory (a), Glasma-like 2+1D (b), and isotropic 3+1D (c), and for scalar excitation in the Glasma-like 2+1D theory (d).

Transverse, longitudinal and scalar correlations
All frequencies and momenta of the correlators are normalized using the mass ω pl , which is the plasmon frequency of gluonic polarizations. While there are multiple ways of determining ω pl [57,58], we extract it from the dotted spectral functionρ T (t, ω, p=0) at vanishing momentum using a (Gaussian) fit function, as explained in Sec. 3.3. At the considered times, the values of the plasmon frequency read ω 2D pl = 0.12 Q for the 2+1D theory, ω 2D+sc pl = 0.2 Q for the Glasma-like 2+1D theory and ω 3D pl = 0.13 Q for the 3+1D theory. The continuous gray lines in Fig. 1 correspond to a relativistic dispersion ω rel = ω 2 pl + p 2 with the numerically extracted ω pl . We show the transverse and scalar HTL dispersion relations ω HTL T /φ (p) from Eq. (2.24) and Eq. (2.26) as black dashed lines. Note that the relativistic dispersion is an ad hoc ansatz with the parameter ω pl extracted by fitting the data. The HTL dispersion relation, in contrast, has a functional form determined by the leading order HTL calculation, and a mass scale determined by the distribution function f (t, p) using Eq. (2.23). For the HTL mass, the resulting values are m 2D HTL = 0.14 Q, m 2D+sc HTL = 0.21 Q and m 3D HTL = 0.15 Q for the considered parameters and times of the different theories.
When considering the correlation functions in Fig. 1 for fixed momentum, one observes an excitation peak as a function of frequency for each of the different non-Abelian gauge theories and their transverse gluonic or scalar contributions. The additionally shown relativistic dispersion ω rel (p) and HTL dispersion relations lie well within the width of the peak. Interestingly, the scalar excitation shows some deviations from a relativistic dispersion, which actually is the leading order HTL prediction (2.26). These deviations are at their strongest at p ∼ ω pl . However, at larger momenta p ω pl the dispersion relation coincides with its HTL prediction, which also agrees with the transverse HTL dispersion for large momenta. As we will see in Sec. 3.4, also the shape of the spectral function agrees between transverse and scalar polarizations at higher momenta.
An important difference to the expected form of the excitations within the hard-loop framework lies in the width of the peaks (considered as functions of ω), which in hard loop theory is a next-to-leading order effect and is thus supposed to be small. We also observe that for gluonic excitations the peaks in 2+1D theories are much wider than in 3+1D or for low-momentum scalar excitations, whose width seems to be comparable to the gluonic width in 3+1D. 1 Transverse gluonic excitations of 2+1D theories at low momenta p ω pl also show sizeable contributions all the way to ω = 0. We will discuss this observation in more detail below in Sec. 3.5. Figure 2 shows the normalized longitudinal statistical correlator EE L (t, ω, p)/ EE L (t, ∆t=0, p) of 2+1D (a), Glasma-like 2+1D (b), and isotropic 3+1D theory (c) for comparison. The gray dash-dotted lines correspond to the free dispersion ω = p, the gray continuous lines to ω rel (p) and the black dashed curves to the HTL dispersion relations ω HTL L (p) from Eq. (2.25). One observes that the peaks at low momenta p ω pl have a similar width and magnitude as the corresponding transversely polarized gluonic peaks of Fig. 1, and agree with the longitudinal HTL dispersion reasonably well. At all momenta, one finds that the longitudinal correlations involve a finite valued continuum for ω p, which corresponds to the Landau cut contribution. Different from transverse dotted spectral functions, the quasiparticle peaks of longitudinal correlators are strongly suppressed at higher momenta p ω pl and Landau damping becomes the dominant contribution, as visible in the figure.
These properties will be studied in more detail by comparing with HTL calculated curves in the following subsection. Note that while we have shown the statistical correlation functions EE α (t, ω, p) here, they appear to be related to the respective dotted spectral functionsρ α (t, ω, p) by the generalized fluctuation-dissipation relation (3.1) even in our farfrom-equilibrium situation, as will also be shown in the following. Therefore, the discussions of this subsection are also valid for the dotted spectral functions in the considered theories.

Comparison of spectral and statistical correlation functions
In the isotropic 3+1D theory, we observed in [40] that the spectral and statistical correlation functions in our considered far-from-equilibrium system are related by the generalized fluctuation-dissipation relation . fluctuation-dissipation relation, which for low momenta p T below the temperature reads EE α (ω, p) = Tρ α (ω, p). In contrast, out of equilibrium their connection may be much more complicated or even absent, as observed for scalar theories [41,42,59]. The HTL framework predicts a generalized fluctuation-dissipation relation [60] with a timedependent effective temperature T eff (t) = EE α (t, ∆t=0, p)/ρ α (t, ∆t=0, p) to hold for sufficiently small momenta p Λ(t). Here Λ(t) is a time-dependent hard scale that dominates the energy density and plays the role of a temperature scale for momenta in the non-equilibrium setting. However, we emphasize that we observe the relation (3.1) to be satisfied in general. This also includes higher momenta p Λ(t), where EE α (t, ∆t=0, p)/ρ α (t, ∆t=0, p) are non-constant functions of t and p.
The validity of the generalized fluctuation-dissipation relation (3.1) in 2+1D non-Abelian plasmas is one of the main results of our work. It is demonstrated in Figs. 3  and 4, which show the normalized statistical and spectral correlation functions for three different momenta for the 2+1D (left panels) and for the Glasma-like theory (right panels) for transverse and longitudinal polarizations, respectively. The normalized statistical cor-  We also observe in these figures that the width and shapes of the correlation functions look very similar in both the 2+1D and Glasma-like 2+1D theories for transverse and longitudinal excitations, when plotted as functions of ω/Q. Note that the scale Q is related to the energy density per polarization state. Thus, adding the scalar degrees of freedom to the system without changing the distribution f (p) of the gauge particles does not seem to modify the shape or the width. In contrast to the width of the peak, the mass in the Glasma-like system is always larger, as seen particularly well by the shift of the peak at lower momenta. Thus, the squared plasmon mass is proportional to the number of polarization states, as expected from HTL in Eq. (2.23). Taking together these two observations seems to indicate that, to a first approximation, the scalar sector contributes to the mass but does not contribute to the damping of gluonic excitations.
If we plotted the correlation functions at the same p/ω pl as functions of ω/ω pl , i.e., normalized by the plasmon mass of the corresponding theory instead of an energy-related QΔt max = 150 Gauss fit Sech fit BW fit Figure 5. The normalized transverse statistical correlator EE T (t, ω, p)/ EE T (t, ∆t=0, p) of Glasma-like 2+1D system for momentum p = 0.3 Q. In the left panel, different ways to compute the Fourier transform are compared, which are seen to agree well. In the right panel, we compare different fit functions for the peak shape, as explained in the main text. momentum scale, the peak positions would overlap by construction. However, the peaks in the 2+1D Glasma-like system would appear more narrow than in the 2+1D theory. This effect, which we have already seen in Figs. 1 and 2, means that in the Glasma-like theory the quasiparticle excitations are slightly longer lived than in pure 2+1D. We now compare the correlators to the transverse and longitudinal spectral functions calculated in HTL (2.19), which consist of a Landau cut contribution and a quasiparticle delta peak. They are shown as green dashed lines in Figs. 3 and 4. Importantly, different from 3+1D gluonic plasmas, the extracted dotted spectral functions are so broad that an accurate distinction between Landau cut region and an excitation peak becomes difficult for both polarizations. Therefore, the HTL expressions mostly do not provide a good description of the nonperturbative simulation results of the 2+1D systems.
One exception is the large-momentum region p m HTL of the longitudinally polarized dotted spectral function, which concerns the central and lower panels of Fig. 4. According to HTL expectations that are discussed in Appendix A, the Landau damping term is expected to dominate the spectral function for these momenta. We indeed observe that the HTL Landau damping curves for low frequencies ω p agree well with the numerical data. This observation, which validates the HTL formulas for this specific case, is nontrivial, since the HTL framework is not expected to work well for the 2+1D theories due to the importance of the missing soft-soft interactions, as discussed in Sec. 2.5.
Note that the sharp peaks of the HTL Landau contributions around ω ∼ p are smoother in the numerically extracted correlators. We have reported of a similar smoothening of the Landau cut region for 3+1D plasmas in [40], which may result from including more mode interactions than considered in HTL at leading order.

Shape of the peaks
Let us now discuss the shape of the excitations in more detail. Computations within the HTL formalism (Sec. 2.5) suggest that the normalized correlation functions can be written as a sumρ with an excitation peak h(ω, p) that originates from a pole at frequency ±ω α (p) with residue Z α (p) and becomes a Delta function at leading order in HTL. A finite peak width should thus be a subleading-order effect. Since the peak results from poles in the retarded propagator (Eqs. (A.10)-(A.12)), its shape could be expected to be of Breit-Wigner form, which becomes a Lorentzian distribution for narrow peaks. A Lorentzian shape for the excitation peaks has been found in 3+1D gluonic systems [40] and would usually correspond to quasiparticle excitations. The peak form is investigated in Fig. 5. In the right panel, the normalized transverse statistical correlator for the larger time window Q∆t max = 150 is shown in red for the Glasma-like theory at momentum p = 0.3 Q. It is compared to fits h(ω, p) that are modeled as a Gaussian, hyperbolic secant and Breit-Wigner distribution, respectively, given by with sech(x) = 1/cosh(x). Here [ω R → −ω R ] denotes the previous term with the indicated substitution. 2 We neglect the Landau cut contribution in the fitting procedure because it is subleading and does not describe our data well at frequencies ω ∼ p, which lie within the width of the excitation peak. Comparing the different fits in the plot, one observes that the Breit-Wigner distribution fails to describe the peak faithfully, which is particularly well visible at the right tail. In contrast, the shape of the peak is well described by the non-Lorentzian distributions h Gauss and h Sech . Indeed, we have included corresponding fits to our data of a Gaussian distribution h Gauss for different momenta and polarizations already in the Figs. 3 and 4 as black lines. The fits are seen to agree well with the shape of the dotted spectral functions for all momenta, where such a comparison is sensible. Note that for longitudinal spectral functions, quasiparticle excitations are expected to be suppressed at higher momenta p m D . Therefore, we only included fits for small momenta p m D in Fig. 4.
We have checked that the non-Lorentzian shape is not an artefact of a finite time window in the Fourier transform or of signal processing techniques. In the left panel of Fig. 5 we show the same correlation function as in the right panel for the time windows Q∆t max = 150 and Q∆t max = 80 with Hann windowing as well as for Q∆t max = 80 without windowing, which implies h(∆t) = 1 in Eq. (2.9). One observes that all curves accurately agree within uncertainties.
This non-Lorentzian peak, thus, appears to be a distinct feature of 2+1D theories in contrast to 3+1D plasmas, where Lorentzian excitations emerge instead [40]. A similar non-Lorentzian shape has been encountered in single-component non-relativistic and relativistic scalar models at low momenta [41,42]. It was also used to distinguish the corresponding excitations from the usual Lorentzian peaks that dominate in O(N )-symmetric scalar models for a large number of components N 1 [42]. For the 2+1D gauge theories considered here, the origin for the observed non-Lorentzian form is currently unknown.

Scalar excitation
We can also compare correlation functions in the (relative) time domain. Consistently with Eq. (3.1), the generalized fluctuation-dissipation relation is also satisfied in ∆t, where α denotes the polarizations T, L, φ. This is demonstrated in the upper panels of Fig. 6 for the 2+1D Glasma-like theory, where the normalized statistical EE α (t, ∆t, p)/ EE α (t, ∆t=0, p) and (dotted) spectral correlation functionsρ α are seen to nicely coincide for scalar and transverse polarizations separately. Interestingly, both polarizations even agree with each other at high momenta p m D (right panels), while they are seen to differ at lower momenta p m D (left panels).
This observation implies that the distinction between transverse and scalar excitations is only relevant for low momenta. We confirm this interpretation in the frequency domain in the lower panels of Fig. 6. In the right panel we show the correlation functions for the same higher momentum and also include a (Gaussian) fit to the scalar correlator given by (3.4). All correlators lie on top of each other, revealing the same shape, dispersion and width among transverse and scalar polarizations. In contrast, at the lowest momentum p = 0, a clear difference is visible between the gluonic EE T and the scalar correlation ππ . The scalar excitation is much more narrow, leading to longer-lived quasiparticles than for transverse gluonic excitations. We have observed this property also in Fig. 1d) for low momenta p m D , which is similar to the gluonic excitations in a spatially threedimensional system. Furthermore, the peak position is shifted towards a slightly higher frequency. Interestingly, this agrees qualitatively with the HTL dispersion relations (Appendix A.3), which predict that the peak positions of the gluonic and scalar excitations at p = 0 are ω pl and m D , respectively, with ω pl < m D .

Low frequency behavior
We now turn to the actual spectral functions ρ α , only having discussed their time derivativė ρ α so far. For the scalar polarization, we see in Fig. 6 thatρ φ as a function of ∆t always oscillates around zero and has a narrow excitation peak in the frequency domain with a vanishing value for ω → 0. As visible in Fig. 7, the scalar spectral function ρ φ also oscillates around zero. This confirms that its Fourier transform ρ φ (t, ω, p) smoothly approaches zero as ω → 0.
In contrast, the dotted spectral function of gluonic excitationsρ T (t, ω, p) is seen to behave differently than the scalar one. This is visible in frequency domain in Figs. 3 and 6, where the spectral function approaches non-vanishing values for ω → 0 at low momenta p m D . As a consequence, the resulting spectral function shown in Fig. 8 (computed as ρ T =ρ T /ω as discussed earlier) behaves as ∼ 1/ω for ω → 0. Thus, although ρ T (t, ω=0) = 0 due to the odd symmetry, the spectral function actually behaves as ρ T (t, ω) ∼ 1/ω for ω → 0, for both positive and negative ω. This is consistent with the behavior of ρ T as a function of ∆t for low momenta. This can be seen in the left panel of Fig. 7 where we observe that ρ T does not oscillate around zero but possibly approaches a nonzero value for ∆t → ∞. This behavior is surprising since, based on the analytical structure of the HTL self-energies, we would expectρ T and ρ T to smoothly vanish as ω → 0 and, consistently, ρ T (t, ∆t, p) to oscillate around zero.  Figure 8. Transverse spectral function ρ T computed asρ T /ω for small momenta p = 0 Q and p = 0.1 Q as a function of ω for 2+1D (left) and Glasma-like 2+1D (right) theories. Due to its symmetry, the spectral function satisfies ρ T (t, ω=0) = 0 identically while it follows ρ T ∼ 1/ω for ω → 0 due to the finite values ofρ T .
Such a behavior could, for instance, follow from having an additional excitation at ω ≈ 0 or a conserved quantity; however, we do not expect either here. Another possible explanation would be if it was related to time reversal non-invariance, since our system is evolving in time. To check this we have performed a similar calculation in the time reversal invariant case of classical thermal equilibrium in Appendix B. However, we also see a finite value forρ T (t, ω=0, p) there for low momenta. Thus, the small violation of time reversal invariance in our system does not seem to be the explanation either.
Another possibility is that this observed feature is associated with our gauge fixing procedure that is discussed in the end of Sec. 2.2. Assuming ∆t t, we fix to a Coulomb-type gauge at the time t, but the system does not exactly remain in Coulomb gauge at t + ∆t. To check whether this approximation may lead to such an effect, one could perform the calculation fully in Coulomb gauge, which would avoid the mentioned approximation but require reintroducing a temporal component for the gauge potential. Moreover, it would be interesting to look at the heavy quark diffusion coefficient (as in Ref. [45]), which is a gauge invariant observable sensitive to the same physical scales. We leave these studies to a further paper. To summarize, we do not have a clear interpretation of the finite value ofρ at ω = 0.

Time dependence of dispersion relations & damping rates
Now we will study how the dispersion relations ω α (t, p) and excitation widths γ α (t, p) of the 2+1D systems depend on the time t. We recall that in the scaling solution (see Sec. 2.3) this can be used to obtain information about the dependence on the soft and hard scales, which follow separate power laws in t. Figure 9 shows our numerical results for the dispersion relations and damping rates at different times Qt for the transverse gluon excitation in the 2+1D theory and for the transverse and scalar excitations of the Glasma-like 2+1D theory. They were extracted by fitting h Gauss (ω, p) to the normalized transverse gluon and scalar excitations, as explained in Sec. 3.3. The main observation here is that all curves lie on top of each other when plotted in terms of the time dependent mass ω pl (t). For the dispersion relation this is relatively trivial, and just shows that the functional form of the dispersion relation does not change in time and only depends on the current value of the mass. We have added HTL and relativistic dispersion relations as black dashed and continuous gray lines, respectively. We observe that both agree sufficiently well with the fitted values to the numerical data. Small deviations lie within the width of the excitation peaks, as discussed in Sec. 3.1 and are visible in Fig. 1.
The Qt depedence of the width, shown in the right panels of Fig. 9, is much more interesting. Although the scaling is not quite as good as for the dispersion relation, the plots show that both the p-dependence and the magnitude of the damping rate follow the time dependence of the plasmon mass. This is different from the 3+1D case, where γ(t) was found to decrease faster than ω pl (t) with time [40]. Thus, while for the 3+1D theory the quasiparticle peaks get narrower as the scale separation between the dynamical hard and soft scales grows with increasing time, this does not happen in the 2+1D theories: the width of the excitation peak relative to the mass remains constant.
As a consequence of this scaling, correlators rescaled with appropriate powers of ω pl (t) remain time-independent as functions of frequency ω/ω pl (t) and momentum p/ω pl (t). This is demonstrated in Fig. 10 at the example ofρ T at p = 0 for both 2+1D systems. Over a wide range of times Qt = 200, 500, 2000, the appropriately rescaled correlation function is indeed stationary. Note that we have not used self-similarity explicitly but merely rescaled all dimensionful quantities with appropriate powers of ω pl (t). In fact, also in 2D classical thermal equilibrium the width of the excitations is of the order of the mass, as we show in Appendix B. This is suggestive of a general relation γ α ∼ ω pl for 2+1D systems.
The behavior γ α ∼ ω pl should be contrasted with the expectation from HTL, where the width is expected to be proportional to the "effective temperature of the soft modes" p / Q γ ϕ / Q Figure 9. The peak position (dispersion relation) ω α (t, p) and the peak width (damping rate) γ α (t, p) as functions of momentum p at different times for (top:) the transversely polarized gluons of the usual 2 + 1-dimensional theory and for (center:) the transversely polarized gluons and (bottom:) the scalar correlator of the Glasma-like 2 + 1-dimensional theory. Dispersion relations, widths and momenta are normalized by the time dependent mass scale ω pl (t) in the main plots and by the time independent energy scale Q in the insets. Additionally, we show, in the left panels, the same HTL dispersion relations ω HTL T /φ (p) and relativistic dispersions ω rel (p) as in Fig. 1 as black dashed and continuous gray lines, respectively. g 2 T * (t). Indeed, in the 3+1D case we have observed [40] that the time dependence of γ α is consistent with this expectation. For the 2+1D theories a simple parametric estimate would lead to the scaling g 2 T * ∼ Q(Qt) −2/5 [36], which is a much stronger time dependence than what we observe since ω pl ∼ Q(Qt) −1/5 . Thus, the behavior of γ α (t) in the 2+1D theories indicates that the phenomenon of plasmon decay happens in a qualitatively different way than in 3+1D. As a consequence, excitations at low momenta p m D (t) for different times remain too short-lived to be considered as quasiparticles. Therefore, while an effective kinetic theory description may exist for larger momenta, its collision kernel, depending on the soft modes, must be determined nonperturbatively.

Conclusion
In this paper we have studied the excitation spectrum of 2+1D classical gluodynamics in the self-similar regime. Our aim was to understand whether the boost invariant system created at the initial stages of ultrarelativistic heavy-ion collisions can be understood in a quasiparticle picture. We have studied two theories. The first one is a genuine 2+1D theory where the longitudinal direction is completely absent. The second theory is Glasma-like, in the sense that it is obtained from a 3+1D theory that is invariant in the longitudinal direction. As a consequence, the gauge field in this direction becomes a scalar field.
One of our main observations is that excitations are broader in 2+1D than in isotropic 3+1D systems. This is seen in the damping rates, which correspond to the widths of the excitation peaks in statistical and spectral correlation functions. In the 3+1D theory, the damping rate is a subleading effect, and an increasing scale separation between the hard and soft scales Λ(t) ω pl (t) (with increasing time t in our case) leads to an increasing separation between plasmon mass and width, ω pl (t) γ(t). In the 2+1D cases, in contrast, the damping rate is of the order of the plasmon mass γ(t) ∼ ω pl (t) at all times. This indicates that the damping rate is not determined by the hard degrees of freedom, but results from nonperturbative interactions between the soft gauge degrees of freedom. This is a qualitatively different dynamical picture than in three spatial dimensions and has the important consequence that gluonic transverse and longitudinal excitations below the mass scale are too short-lived to form quasiparticles.
Unlike the gauge fields, the scalar fields do have narrow excitations at low momenta. At larger momenta, the damping rates saturate and the transverse and scalar excitation peaks become identical.
Similarly to gluonic plasmas in 3+1D, we observe that the statistical and spectral correlation functions obey a generalized fluctuation-dissipation relation. However, in contrast to 3+1D, the excitation peaks have a non-Lorentzian shape, that is reasonably well described by Gaussian or hyperbolic secant distributions.
We have also performed simulations in classical thermal equilibrium in 2+1D and observed qualitatively similar behaviour as along the self-similar attractor. We interpret this as a sign that the qualitative features may be generic to 2+1D gauge theories.
Our results indicate that an effective formulation of kinetic theory in 2+1 dimensions should indeed be possible for large momenta p ω pl . However, this is complicated by the fact that in lower dimensions infrared effects become more important. In particular, in 2+1D the dynamically generated screening scale gets equal contribution from all momenta, whereas in 3+1D the dominant contribution comes from hard particles. As a consequence, the HTL approximation breaks down already at the Debye mass scale (instead of the magnetic scale in 3+1D). This picture is supported by our simulations. Therefore obtaining effective matrix elements for kinetic theory is a nonperturbative problem where the dynamics have to be deduced in a self-consistent way taking into account contributions from the plasmon mass scale.
This work also contributes to the study of thermalization in heavy-ion collisions concerning the role of plasma instabilities, which is still under debate in the literature. In particular, there is a persistent discrepancy between classical Yang Mills simulations [29,30] and corresponding HTL-theory calculations (including simulations and analytical calculations) [17,18,61,62]. In the case of full 3+1D Yang-Mills theory, instabilities seem to play a much smaller role in the isotropization than in the case of HTL. This observation has been used as an argument to justify kinetic-theory descriptions without instabilities in phenomenological applications [32-34, 63, 64]. Our results show that, as expected, gluonic excitation modes at low momenta are not correctly described by HTL perturbation theory in 2+1D and 3+1D plasmas at extreme anisotropy. The presence of these nonperturbative corrections in the purely 2+1D case suggests that there may be significant nonperturbative corrections present also in systems with finite but large anisotropy, which are relevant to phenomenological applications and may alter the role of instabilities. For instance, this may cure some issues encountered in HTL perturbation theory where observables cannot be computed due to emerging instabilities in anisotropic systems. To approach these questions in more detail, we plan to extend our simulations to expanding and anisotropic plasmas.
CoG-681707, by the EU Horizon 2020 research and innovation programme, STRONG-2020 project (grant agreement No 824093) and by the Academy of Finland, project 321840. This work was funded in part by the Knut and Alice Wallenberg foundation, contract number 2017.0036. The authors wish to acknowledge CSC -IT Center for Science, Finland, for computational resources. The computational results presented have been also achieved in part using the Vienna Scientific Cluster (VSC). The content of this article does not reflect the official opinion of the European Union and responsibility for the information and views expressed therein lies entirely with the authors.

A.1 Polarization tensor
In the HTL formalism, a crucial quantity is the polarization tensor Π µν . For a non-Abelian SU(N ) gauge theory in d spatial dimensions, its leading order (LO) expression reads [37,38] where ∂/∂ p 0 = 0, with Minkowski metric g νβ = diag(1, −1), → 0, V = (1, p/ω p ), K = (ω, k) and with the approximation ω p ≈ p for the dispersion relation, which is valid for sufficiently large momenta. In general, the frequency ω and momentum k are independent. The polarization tensor is gauge invariant, symmetric Π µν = Π νµ and transverse with K µ Π µν (K) = 0 for all ν. For the latter to be true, the distribution function has to vanish at large momentum lim p i →∞ f p = 0 for directions i, which is fulfilled in practice.
We have summed different bosonic contributions into the distribution functionf = N c λ f (λ) =: N c d pol f (t, p), where each f (λ) corresponds to the distribution of the fields with non-longitudinal polarization λ. For the 2+1D theory, one has d pol = 1 and the integration reads p = d 2 p/(2π) 2 . When considering the Glasma-like theory, the distribution function has the form f 3D (t, p) = 2πδ(p z ) f (t, p ⊥ ), such that integration becomes identical to the 2+1D theory. Since the scalar contributions correspond to a polarization in z direction, the distribution function f (t, p) becomes an average over gauge and scalar distributions, i.e., d pol f (t, p) = f G + f φ with d pol = 2 in the Glasma-like theory. Note that the resulting expressions are the same for both 2+1D and Glasma-like theories with only a different d pol .
Using the isotropy of f (t, p) and performing an integration by parts, the expression (A.1) can be cast into the form with i, j = 1, 2 and The transverse and longitudinal polarizations of Π ij (K) can be computed as with the transverse vector q = (−k 2 , k 1 ) such that q · k = 0. Evaluating the residual integral [65], one arrives at For the Glasma-like theory, we also need the scalar component Π φ = Π zz . To obtain it, we start from (A.1) and perform an integration by parts resulting from p z = k z = 0.

A.2 Retarded propagator
Next we consider the retarded propagator in temporal gauge G HTL ij (ω, p), as employed in this work, and ask how these components of the polarization tensor emerge. Using a tensor decomposition as suggested in Refs. [17,66] with where for the Glasma-like theory the normal vector is simply the unit vector in z direction n = e z , the propagator can be written as with transverse, longitudinal and scalar polarizations (A.12) One immediately observes the same screening properties in the static limit as in the 3 + 1dimensional case, with additional to the static screening of scalar fields with lim ω→0 −ω 2 + p 2 + m 2 D = p 2 + m 2 D . This gives the interpretation of m 2 D as the Debye mass, as anticipated by its notation.

A.3 Spectral function
Each polarization of the spectral function can now be computed as the imaginary part of the respective retarded propagator 3 ρ HTL α (ω, p) = 2 Im G HTL α (ω + i , p) , (A.14) with α = T, L, φ denoting the polarization. They satisfy the sum ruleṡ which are the same as in 3+1D for T and L.
The HTL spectral functions can be further split into a Landau damping part for |ω| < p and a quasiparticle part with the dispersion relations ω HTL α (p) and residues Z α (p) of quasiparticle excitations as discussed below.
The where the dispersion relation for scalar particles is simply the relativistic dispersion of free particles. It is interesting to study their behavior at low and high momenta. These are  The leading asymptotic behavior is similar to the d = 3 case. Like there, one has ω HTL T (p) > ω HTL L (p) for p > 0. The large-momentum behavior for transversely polarized quasiparticles agrees with the scalar dispersion, and hence, with a relativistic dispersion relation, which shows that m D can also be interpreted as the asymptotic mass. On the other hand, the longitudinally polarized quasiparticles do not have an asymptotic mass, just as in the 3+1dimensional case. An important difference to the latter is, however, that their dispersion relation approaches the ultra-relativistic limit ω p considerably slower than in the d = 3 case where the approach is exponential.
The residues at the quasiparticle peaks are defined as where we dropped the HTL label and the dependence on momentum of the dispersion relations to shorten the expressions. Thus, also in the 2+1D case one expects contributions from longitudinal quasiparticles to decrease at high momenta. However, it is a power law decrease, as opposed to the exponential decrease of 3+1D gauge theory. Note that the leading contribution of Z T for p m D and p m D and Z L for p m D are the same in 2+1D and 3+1D. This implies that the sum rules (A.15)-(A.16) for the spectral functions are dominated by quasiparticle contributions for p m D for both polarizations and for p m D for transverse polarizations, while the longitudinal sum rule for p m D is dominated by the Landau damping contribution.  Figure 11. The correlatorρ T as a function of frequency for momenta p = 0 T (left) and p = 0.6 T (right) for a 2+1D theory in classical thermal equilibrium. Time translation invariance is shown by comparing the correlations at two different times. One finds the same properties as for the non-equilibrium case considered in the main part of the paper. Note that in the right panel the curve for T t = 100 is barely visible because it accurately coincides with the curve for T t = 400.

B Cross check: classical thermal equilibrium
Here we compute the correlators of the 2+1D theory in frequency space in classical thermal equilibrium, and compare the results to the non-equilibrium systems discussed above. To obtain the classical thermal state, we initialize our system with f (t = 0, p) = T 0 /p on a 256 2 lattice with lattice spacing T 0 a = 0.5 and coupling g 2 = 0.1 T 0 . After the initialization we restore the Gauss law constraint and let the system evolve. We find that the equaltime correlators EE T /L (t, ∆t=0, p) become stationary for times T 0 t 20, marking the onset of classical thermal equilibrium. Due to the small Debye mass m 2 D ∼ g 2 T T 2 , the resulting thermal state has almost the same temperature as the initial temperature parameter T ≈ T 0 .
Next, we extract correlation functions at different times T t = 100 and 400 to demonstrate time-translation invariance explicitly, as should be the case in thermal equilibrium. Our results are shown for the transverse dotted spectral functionρ T (t, ω, p) in Fig. 11 for momenta p = 0 (left) and p = 0.6 T (right). One observes that the curves corresponding to different extraction times lie indeed on top of each other within statistical uncertainties. For higher momenta, where error bars are small, the curves can be barely distinguished by eye (right panel). This confirms that we have reached classical thermal equilibrium.
Our numerical results qualitatively agree with our findings in the non-equilibrium systems above. We can summarize our most important findings as follows: Sec. 3.5 is not related to the violation of time reversal invariance, and should have some other physical origin. 4. The functional form is again non-Lorentzian and HTL does not provide a good description of the data (except for the longitudinal polarization at high momenta).
All of this indicates that the observed features in non-equilibrium 2+1D theories are not special to the considered self-similar attractor but also arise in the time-translation invariant setting of classical thermal equilibrium, and, presumably, for other non-equilibrium states.