Highly occupied gauge theories in 2 + 1 dimensions: A self-similar attractor

K. Boguslavski , A. Kurkela , T. Lappi , and J. Peuron 6 Institute for Theoretical Physics, Technische Universität Wien, 1040 Vienna, Austria Department of Physics, University of Jyväskylä, P.O. Box 35, 40014 University of Jyväskylä, Finland Theoretical Physics Department, CERN, Geneva, Switzerland Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway Helsinki Institute of Physics, P.O. Box 64, 00014 University of Helsinki, Finland European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*) and Fondazione Bruno Kessler, Strada delle Tabarelle 286, I-38123 Villazzano (TN), Italy


I. INTRODUCTION
A characteristic feature of many highly occupied systems is that they often approach universal self-similar attractors, also referred to as nonthermal fixed points (NTFP) [1,2]. Examples have been found with classical field methods in various theories in three spatial dimensions (3D) including non-Abelian gauge theories, relativistic and nonrelativistic scalar field theories [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16], and in two-dimensional scalar systems [17][18][19][20]. Nonthermal fixed points have recently also been found experimentally in ultracold atom experiments [21,22]. A kinetic theory description of the underlying theory is often a natural way to explain the existence and properties of such fixed points [1,14,[23][24][25][26][27]. These NTFPs appear because the interaction rate of the initial conditions is faster than that of the final equilibrium state. Therefore, the system loses memory of its initial conditions faster than it reaches thermal equilibrium and, hence, stays in a state that is not thermal yet but does not remember details of its initial conditions.
Much less is known about two-dimensional (2D) gauge theories. 1 Differently from the three-dimensional case where an effective kinetic theory has been formulated to leading order accuracy [23], IR effects play a stronger role in 2D due to the lower dimensionality. As we will discuss, the hard (thermal) loop (HL) treatment used to regulate the Coulomb divergence of elastic scatterings in 3D is insufficient in the two-dimensional case. Thus, it is a priori not obvious whether or to what extent quasiparticle descriptions are applicable and whether the system can exhibit self-similar behavior.
Apart from these theoretical questions, this uncertainty has also conceptual consequences for our understanding of the thermalization (hydrodynamization) process in ultrarelativistic heavy-ion collisions. In this context, nonlinear interactions of gluons produced at central rapidities have been argued to lead to a transverse momentum scale Q s ≫ Λ QCD up to which gluonic fields are of order A ∼ 1=g [29], where g is the gauge coupling. If this saturation scale is sufficiently hard, the system is weakly coupled α s ðQ s Þ ≡ g 2 =ð4πÞ ≪ 1 and consists of highly occupied "Glasma" color fields [30]. These are initially approximately boost-invariant along the beam axis and can be described by a 2D classical Yang-Mills field theory. Therefore, it is interesting to ask to what extent hard-loop theory and quasi-particle approximations are applicable also to extremely anisotropic media and to 2D theories. Adjacent to this is the question of what is the earliest time during the thermalization process when kinetic theory can be used to describe the dynamics.
In this paper, we will study whether highly occupied two-dimensional gauge theories approach a universal selfsimilar attractor and whether their scaling properties can be understood with simple kinetic theory arguments. We will consider two related SU(2) gauge theory systems using a classical lattice formulation. 2 The first one of these systems is a 2 þ 1-dimensional gauge theory (2D). The second theory also includes a scalar field in the adjoint representation of the gauge group and will be denoted as "2D þ sc." The latter corresponds more closely to the situation in the initial effectively two-dimensional stage of a heavy-ion collision, since it is the theory one obtains by dimensional reduction starting from a 3 þ 1-dimensional pure gauge theory and restricting it to field configurations (and gauge transformations) that do not depend on the longitudinal spatial coordinate.
Our main result is that we indeed observe a self-similar scaling behavior for the hard modes of both theories that can be explained using parametric considerations in kinetic theory. Some evidence for such scaling behavior was seen in [34], where the focus was more on the determination of the plasmon frequency, and in the present work we establish with different methods the existence of the NTFP. While we focus here on the dynamics of hard modes, questions concerning HL and quasiparticle descriptions of soft momentum modes p ∼ m D will be studied with unequal-time correlation functions in classical-statistical simulations in a forthcoming work. This paper is structured as follows. In Sec. II, we briefly discuss the two theories we are studying, the initial conditions used, and the numerical algorithm. Then in Sec. III we present the results from the numerical calculations. In Sec. IV, we derive the observed scaling exponents from a kinetic description. We conclude and outline some potential future work in Sec. V. The appendices cover details of our approach and of our analysis.

A. Theories and initial conditions
We consider non-Abelian SUðN c Þ gauge theories with N c ¼ 2 in d ¼ 2 spatial dimensions. The starting point is the classical gauge field action c − 1 and Lorentz indices μ; ν ¼ 0; …; d imply summation over them. Using the generators Γ a of the suðN c Þ algebra, the gauge field in fundamental representation reads A μ ¼ A a μ Γ a . We study the following two theories: 2D gauge theory in d ¼ 2 spatial dimensions, with the Yang-Mills action (1). 2D þ sc the same gauge theory supplemented with an additional scalar field in the adjoint representation of the gauge group. This corresponds to a classical action with an adjoint scalar field ϕ a and Here the summation is over j ¼ 1, 2 and the covariant derivative is D ab j ¼ δ ab ∂ j − gf abc A c j . This theory can be obtained from Yang-Mills theory in three spatial dimensions by dimensional reduction, assuming that the field configurations do not depend on the third coordinate x 3 . To maintain this symmetry, also gauge transformations are not allowed to depend on x 3 , turning the third component of the gauge field into a scalar A a 3 ≡ ϕ a . Note that in 2D the dimensionalities of the fields and coupling constants are different from the 3D case. The action must be dimensionless ½S YM ¼ ½S 2D YM ¼ ½S 2D ϕ ¼ 0, from which one easily deduces that ½g ¼ 1=2 and ½A ¼ ½ϕ ¼ 1=2. The dimensionality of the interaction term of the covariant derivative has to be that of a derivative ½gA ¼ 1, as in three spatial dimensions.
The theory 2D þ sc is the nonexpanding space-time analogy of the boost-invariant Glasma, while the 2D theory also drops the adjoint scalar contribution. Therefore, both theories can be considered as simplifications of the Glasma state. Note that there is only d pol ¼ 1 transverse polarization in 2D, while the 2D þ sc theory, originating from a 3D system, has d pol ¼ 2 transverse polarizations: one from gauge d.o.f. and one from adjoint scalars.
The systems are initialized, using the method described in Sec. II B, with a field configuration that has a chosen single-particle distribution function fðt; pÞ at the initial time Qt ¼ 0. Here Q is a conserved momentum scale characterizing the system and will be defined in (5). We consider weakly coupled g 2 =Q ≪ 1 but highly occupied f ≫ 1 initial conditions of the form 2 We expect the results to carry over to SU(3) theories as well. The qualitative agreement of weakly coupled SUðN c Þ theories far from equilibrium for N c ¼ 2, 3 has been observed for different phenomena [31][32][33].
for gauge and scalar fields. Unless stated otherwise, we will use n 0 ¼ 0.1, for which p 0 ¼ Q for our chosen definition as detailed below. As we will show in Sec. III, the exact form of the initial conditions and the values of p 0 and n 0 separately are not relevant after a transient time, since the systems will approach an attractor solution that only depends on Q.
To define the characteristic momentum scale Q, note that the energy density is a conserved quantity in the systems studied here and can be computed in a gauge-invariant way in classical-statistical field theory (e.g., in [5,9]). The combination that we have access to in the classical field formulation is the energy density scaled with the coupling g 2 ε, which has the momentum dimension ½g 2 ε ¼ 4. This allows us to define a constant momentum scale in a gaugeinvariant way as with the number of transverse polarizations d pol and the dimension d A ¼ N 2 c − 1 of the adjoint representation. The constant C is taken as C ¼ 20 ffiffiffiffiffi ffi 2π p ≈ 50, a value that has no physical meaning and has been chosen for convenience such that for n 0 ¼ 0.1 one indeed has p 0 ¼ Q. Recall that the coupling constant g is now dimensionful: if one keeps the dimensionless combination g 2 =Q constant, it is easy to see that (5) leads to the proportionality Q ∝ ffiffi ffi ε 3 p that is natural for a scale derived from a two-dimensional energy density. The scale Q will be used to measure all dimensionful quantities.

B. Semiclassical simulations
At high occupation numbers, we can use the classicalstatistical approximation to study the dynamical evolution of systems far from equilibrium [35,36]. A description of this standard technique can be found, for instance, in Refs. [9,37]. Then all fields are classical and discretized on a cubic lattice of size N 2 s with lattice spacing a s . The real-time dynamics results from solving a gauge-covariant formulation of the classical Hamiltonian equations of motion in temporal gauge A 0 ¼ 0 in a leapfrog scheme for the gauge-covariant link fields U j ðt; xÞ ≈ expðia s gA j ðt; xÞÞ and chromoelectric fields E j a ¼ ∂ t A a j . For the theory 2D þ sc, we use j ¼ 1, 2, 3 in the same scheme.
The fields can be initialized in momentum space by requiring that the transversely polarized fields 3 at Qt ¼ 0 follow with π a ≡ E 3 a , while other combinations for two-point functions vanish, as well as hAi ¼ hEi ¼ 0. 4 Since such initial conditions violate the Gauss law constraint, the latter is restored with the algorithm from [38] before starting the dynamical evolution. Alternatively, we also started with initial conditions with E j ¼ 0 but twice the amplitude n 0 , where the Gauss constraint is automatically satisfied and the energy density approximately the same. We found that both lead to the same dynamics after a short transient time.
We will be especially interested in observables in momentum space. For that, we fix the gauge to Coulomb-like gauge ∂ j A j ¼ 0 at the measurement time (see also Appendix A) and Fourier transform the fields according to Aðt; pÞ ¼ R d d xe −ip·x Aðt; xÞ. A central quantity of interest is the single-particle distribution function fðt; pÞ, which provides the occupation number density of the system at a given time and momentum. One can define the distribution function using different correlators, such as in Eqs. (6)- (9). Unless stated otherwise, our standard definition will be where we will set the dispersion ωðpÞ ¼ p neglecting soft scale effects since we are primarily interested in the dynamics at high momenta. We also used an abbreviated notation and similarly for other correlators. 3 The fields at each momentum p ¼ ðp 1 ; p 2 Þ can be split into a transverse and longitudinal part A a j ðpÞ ¼ A a T;j ðpÞ þ A a L;j ðpÞ, such that A a T ðpÞ · p ¼ 0 while jA a L ðpÞ · pj ¼ jA a L ðpÞjp. 4 In practice, this is achieved by setting A a j ðt ¼ 0; pÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi fðt ¼ 0; pÞ=p p c a ðpÞv T;j ðpÞ, and similarly for the other fields, with complex-valued Gaussian random numbers c a ðpÞ and the transverse polarization vector v T ðpÞ.
We note that the classical-statistical framework with the initial correlators [Eqs. (6)-(9)] corresponds to an accurate mapping of the corresponding quantum field theory onto a classical-statistical field theory in the limit of weak coupling g 2 → 0 for g 2 fðt; pÞ held fixed [35,36,39]. In this limit, the vacuum 1=2 contribution to the distribution function is suppressed by g 2 and is negligible for the observables studied here. Therefore, contributions at the order of the lattice cutoff p ∼ 1=a s ≫ Λ are nonphysical in our case and we use different discretizations to verify that our results are not sensitive to a s (or to the lattice volume).

A. Universality and self-similarity
In this section, we demonstrate numerically that starting from initial conditions with high occupation numbers, both theories 2D and 2D þ sc approach a common nonthermal fixed point at high momenta where the distribution function follows a self-similar evolution fðt; pÞ ¼ ðQtÞ α f s ððQtÞ β pÞ: ð13Þ In order to constitute a universal nonthermal fixed point, the scaling exponents α, β and the scaling function f s ðpÞ in Eq. (13) must be the same for different initial conditions. The time evolution at the fixed point only depends on a single conserved quantity, which is the energy density ε in the case of an energy cascade to higher momenta [1,25,39], as is the case here. Hence, the distribution function becomes insensitive to details of the initial conditions after a transient evolution when rescaled with the only dimensionful scale Q determined by ε. This attractor property is illustrated in Fig. 1 where we show the gauge distribution f E as a function of momentum in a double-logarithmic plot for both theories. Dashed lines correspond to different initial conditions at Qt ¼ 0, where the fields are constructed to reproduce the chosen momentum distribution according to Eqs. (6)- (9). Each initial condition was used in both theories and the figure shows their f E also at the later time Qt ¼ 4000, where dashed-dotted lines correspond to the 2D þ sc theory and full lines to the 2D theory. Although resulting from different initial conditions and theories, all six distributions at Qt ¼ 4000 lie indistinguishably on a single curve. This demonstrates that after some transient time that depends on details of the initial conditions, systems in both theories get close to the same attractor. There they follow a universal evolution, insensitive to their original initial conditions.
In this universal regime, the distribution function becomes self-similar, following Eq. (13). This is demonstrated in Fig. 2 for the 2D theory. The upper panel depicts the distribution function in the universal regime at several vastly different times 5 Qt ¼ 75-16000. The lower panel shows this same data in rescaled coordinates: the rescaled gauge distribution ðt=t r Þ −α g 2 f E =Q is shown as a function of rescaled momentum ðt=t r Þ β p=Q. The scaling indices α and β have been numerically extracted to produce the best overlap of the distribution functions at the different times employing a least-square fit procedure [9,14] as detailed in Appendix B, leading to best-fit values The first combination results from energy conservation and that its best-fit value is consistent with zero is a consistency check of the procedure. For Fig. 2 (as well as for all the following figures), we use the values that will be derived in Sec. IV, FIG. 1. The gauge distributions f E at the initial time Qt ¼ 0 for three different initial conditions are shown by dashed lines; note that these same initial conditions are used for both 2D and 2D þ sc theories. At a later time Qt ¼ 4000, full lines show the distributions from these initial conditions in the 2D theory and dashed-dotted lines in the 2D þ sc theory. These six curves overlap so well that they are indistinguishable in this plot, demonstrating the attractor property of the common nonthermal fixed point for both theories. 5 For Qt ≤ 4000 a 768 2 lattice with lattice spacing Qa s ¼ 1=12 has been used, for the later time we used a 256 2 lattice with Qa s ¼ 1=16, where the first two points of the latter were omitted due to volume artefacts. These artefacts occur when the lattice is too small to contain the screening mass m D . In this situation, the smallest momentum modes are artificially enhanced. We checked that simulations of both discretizations coincide otherwise for Qt ≤ 4000.
which are consistent with the extracted ones in Eqs. (14) and (15). The good overlap of the different curves obtained at different times demonstrates scaling behavior.
A similar conclusion can be drawn for the 2D þ sc theory. In Fig. 3, the rescaled gauge and scalar distributions f E and f π are shown in the upper and lower panels, respectively. 6 For comparison, the original curves are depicted in the insets. For the rescaling of amplitudes and momenta, the same exponents α and β have been used as in Fig. 2 for the 2D theory. One indeed observes that at high momenta the rescaled gauge distributions as well as the scalar curves lie on top of each other within error bars.
This form also agrees with the 2D theory, which can be seen by comparison to the dark-blue curve.
The scaling function f s ðpÞ of the gauge distribution consists of a power law ∝ ðp=QÞ −σ at lower momenta and a steep drop at high momenta. This closely resembles the nonthermal fixed point in 3D theory [3,5,8], which also consists of a power law at low momenta and a steep decrease at high momenta. The power law at small momenta is consistent with σ ¼ 1, which can be seen in the lower panel of Fig. 2, where a power law with σ ¼ 1 is also displayed. Small deviations from this power law occur at very low momenta below the Debye mass that is indicated by the blue arrow. This value σ ¼ 1 corresponds to a classical thermal distribution T eff =p at low momenta with a time-dependent effective temperature T eff [24]. Analytical considerations in effective kinetic theory suggest that the form of the distribution function in the infrared should take this form also out of equilibrium [25]. However, numerical classical Yang-Mills simulations in 3D theory have found large corrections to this quantity at early times leading to extractions of σ ≈ 1.3 [3,5]. Such corrections to the classical thermal distribution seem to be absent in these 2D theories.
The main difference between gauge and scalar distributions is that the scalar curves f π are enhanced for low momenta p ≲ m D roughly below the Debye mass that will be discussed in Sec. III C. This IR enhancement can be seen in the lower panels of Fig. 3 and, further below, of Fig. 6. It may seem similar to the case of OðNÞ-symmetric scalar field theories where an IR region has been observed [2,17] or may even be connected to nontrivial topological structures [40]. However, since this enhancement is not part of the self-similar region at high momenta, a detailed study is beyond the scope of this work and is left for further study elsewhere.

B. Gauge-invariant hard scales
So far, we have observed that the gauge-fixed distribution functions lose their memory on details of their initial conditions and approach a self-similar attractor. To confirm this behavior with gauge-invariant observables, we also compute the time evolution of manifestly gauge-invariant measures of the hard scale [5,8] hðD ab k D bc k ϕ c ðt; xÞÞðD ad l D de l ϕ e ðt; xÞÞi: These provide typical hard momentum scales that are expected to grow as Λ 2 ðtÞ ∼ t −2β in the self-similar regime. This can be seen from their perturbative expressions where all higher powers in the field amplitude were neglected and the Coulomb gauge condition was used. Note that because of Q 4 ∝ g 2 ε, the hard scales can be interpreted as ratios ∝ R d 2 pp 2 ωðpÞf= R d 2 pωðpÞf, characterizing the momentum scale that dominates the energy density.
The gauge-invariant hard scales Λ 2 ðtÞ, rescaled with t 2β , are shown in Fig. 4 in a linear-logarithmic plot as dashed lines for the gauge and scalar sectors of 2D and 2D þ sc theories for p 0 ¼ Q initial conditions. The data points of matching color indicate the respective perturbative expressions Λ 2 pert ðtÞ that are obtained by integrating the gauge-fixed distribution functions according to Eq. (19). The good agreement between points and lines of the same color confirms our interpretation of the hard scales as E=π ≈ Λ 2 pert;E=π . Moreover, hard scales from different sectors and theories agree well Λ 2 E ≈ Λ 2 π . One observes that for Qt ≳ 75, the rescaled hard scales become approximately constant, indicating Λ 2 E=π ðtÞ=Q 2 ∝ ðQtÞ −2β . The onset time of self-similar scaling and the value for β employed in Fig. 4 are the same as used for the self-similar evolution in Fig. 2. A similar power law evolution of the hard scale in 2D þ sc theory has been observed in Ref. [34]. This consistency between gauge-invariant and gauge-fixed observables confirms the emergence of a self-similar attractor.
In general, the approach to the attractor depends on the initial conditions and on the observables. This is illustrated for the 2D theory in Fig. 5, which shows the evolution of the hard scale in the 2D theory for initial conditions with different values of p 0 . The flattening of this observable indicates the onset of a power-law evolution. One observes that for the hard scales, the onset of scaling becomes slower with larger p 0 or, equivalently, with lower amplitude n 0 (cf., Fig. 1).

C. Longitudinal polarization and Debye mass
While distribution functions are useful for a comparison to kinetic theory, one can extract further information about a system by studying more general correlation functions. Due to our definition of distribution functions in Eq. (10), we have already discussed the evolution and properties of the transversely polarized equaltime correlator hE T E Ã T iðt; pÞ ≡ pf E ðt; pÞ and its scalar analogy hππ Ã iðt; pÞ ≡ pf π ðt; pÞ. Let us now focus on the longitudinally polarized equal-time correlation function hE L E Ã L iðt; pÞ. It is shown in Fig. 6 at fixed time Qt ¼ 2000, together with the transverse polarization and the scalar correlator for 2D in the upper and for 2D þ sc theory in the lower panel. Every correlation function is shown for two different sets of discretization parameters, written in terms of the lattice length L ¼ a s N s and the lattice spacing a s . The good agreement among curves of different volumes and lattice spacings indicates the absence of numerical lattice artefacts.
The correlators hE T E Ã T i in Fig. 6 are flat up to a high momentum p ∼ Λ beyond which they decrease rapidly, which is, of course, equivalent to our previous observation that fðt; pÞ ∼ 1=p up to a hard scale. Similarly, we have observed the IR enhancement of the hππ Ã i correlator and its agreement with hE T E Ã T i at higher momenta already at the example of f π . The longitudinally polarized correlator hE L E Ã L i approaches hE T E Ã T i at the lowest momenta, while strongly differing for momenta above the Debye mass m D . Indeed, as known in thermal equilibrium [5] and also observed far from equilibrium in the 3D theory [37], the longitudinal correlation function follows the form: for momenta p ≲ Λ. In the late-time limit and in thermal equilibrium, one then expects δ ¼ 0. We have fitted this form to hE L E Ã L i and included the fits in Fig. 6 as black dashed lines. They are seen to accurately describe the correlator.
Fitting this form to the longitudinal correlator at different times, we have extracted the evolution of the fitting parameters A, δ and m D . As expected from our previous discussions, the amplitude quickly approaches an A ∼ t α−β power law and the deviation δ monotonously decreases from δ ≈ 0.2 to 0.3 at early times to δ ≈ 0.08-0.12 at time Qt ¼ 2000 for both theories.
Most interestingly, the fitting procedure allows us to extract an estimate for the Debye mass m D from the p-dependence of the correlator. Its time evolution is shown in Fig. 7. In the main figure, the normalized m 2 D =d pol is plotted as a function of time on a double-logarithmic panel. One observes that the curves stemming from the different theories almost coincide while they lie far apart in the inset where m 2 D is depicted. This indicates that m 2 D scales with the number of d.o.f. d pol , which are 1 for 2D and 2 for 2D þ sc theory. Moreover, m 2 D is observed to approach a t 2β power law evolution that is represented by a black dashed line. Its power law evolution sets in roughly at the same time scale as for the hard scale in the upper panel of Fig. 4. These observations suggest a relation where we used energy conservation in the last line and where f is the amplitude at hard momenta p ∼ Λ. Since The leading order HL expression for the Debye mass is [41]  where f π ≡ 0 for the 2D theory and where fðt; pÞ is an average distribution. Because of f ∼ 1=p at low momenta, or even steeper for the scalar distribution, the integral diverges in the IR in two spatial dimensions and needs to be regularized by some cutoff at the scale m D . This leads to bringing a logarithmic correction to the estimate (21). The HL expression (22) was used in Ref. [34] as one of three methods to extract the mass scale. In all the methods employed, the mass followed an approximate power law evolution with m D =Q ∼ ðQtÞ β with values for β that are roughly consistent with our results. We will return to the discussion of the physical interpretation of this logarithm in Sec. IV.

IV. SCALING BEHAVIOR IN A KINETIC THEORY PICTURE
The nonequilibrium evolution of gauge theories can also be studied using an effective kinetic theory setup. In [23], an effective kinetic theory has been formulated for d ¼ 3 spatial dimensions which is defined by a Boltzmann transport equation where f is the distribution function of gluons and where the effective collision kernel C½f is the sum over the relevant elastic and inelastic scattering processes between the particles. Many of the features of the over-occupied UVcascading system have been well understood in terms of such a kinetic theory description [4,5,25]. This effective kinetic theory describes the evolution of modes at momentum scales well above the screening scale p ≫ m D . With the assumption that scattering against modes that carry soft momenta is subdominant compared to the scattering with the hard particles, this effective description may be used to follow the time evolution of the hard scale at late times when a scale separation between the soft and hard scales has developed. The soft scale makes its entrance to the kinetic theory because of the Coulomb-divergent t-and u-channel elastic scattering amplitudes 7 jMj 2 vacuum ∼ g 4 =ðq 2 ⊥ Þ 2 appearing in the elastic part of the collision kernel pÞ and (d þ 1)-momenta P. In medium, the Coulomb divergences are regulated by the physics of screening, taking place at the momentum transfer scale q ⊥ ∼ m D . In d ¼ 3 dimensions, the particles at the hard scale screen the most, and the scale separation between the soft and the hard scales allows one to describe the screening in simple effective theory, the hard-loop effective theory. Because of this simplification, the effective screened matrix elements may be solved analytically to complete the effective kinetic theory.
Similarly, in d ¼ 2 dimensions, the soft and hard scales are parametrically separated at late times allowing for a quasiparticle description of the hard modes. In contrast to three dimensions, however, the equation for the Debye mass (22) in two dimensions gets equal contributions from all momentum scales m D < p < Λ, such that soft modes contribute equally to screening. This implies that the modes at the soft scale m D interact among each other in a nonperturbative way, reminiscent of the ultrasoft magnetic  Fig. 6 and shown for 2D and 2D þ sc theories. The black dashed line corresponds to a power law t 2β with β ¼ −1=5, leading to Eq. (21). 7 Here jMj 2 is expressed in a nonrelativistic normalization related to the usual relativistic normalization by jMj 2 ¼ jMj 2 = ð16pkp 0 k 0 Þ. scale in three spatial dimensions. The practical implication of this is that the dynamics of the soft modes is no longer governed by a simple effective hard-loop theory. Because of this complication, it is not straightforward to analytically find the forms of the required effective matrix elements and formulate a leading-order accurate kinetic theory.
Nevertheless, while we do not have access to a formulation of the kinetic theory that would be accurate at a numerical level, we may still include parametric considerations in this picture following Ref. [24]. In particular, even while we do not know analytically how the effective elastic scattering matrix element looks like, the elastic scattering amplitude for t-channel is regulated by the soft scale and to parametric accuracy is where q ⊥ is the transverse momentum transfer and the regulator m D is at the correct scale but its form is simplified and not correct at a numerical level. The hard particles moving in the medium experience then soft scatterings with other particles in the medium with a (soft) elastic scattering rate of This is a two-dimensional version of the usual kinetic theory relation expressing the rate in terms of the cross section and number density of scattering targets. In two dimensions, there is only one transverse direction, and thus the squared amplitude gives a rate differential in a onedimensional transverse momentum dq ⊥ . The factor R d 2 pf accounts for the number of density particles in the medium against whom the collision occurs and ð1 þ fÞ ∼ f is the final state Bose enhancement factor. For a thermal-like infrared spectrum f ∼ T Ã =p with an effective temperature T Ã , the integral over dp gets equal contributions from all scales, such that collisions with soft particles are equally frequent as those with hard ones. This is in contrast to 3D, where because of the higher dimensionality and larger phase space at high p most of the scatterings take place against hard particles. Because the kinetic theory describes the evolution of hard modes only, the kinetic theory framework does not numerically describe the collisions against the soft particles. We will, however, neglect this complication here and proceed with our analysis assuming that p ∼ Λ in Eq. (27), with the understanding that our accuracy is purely parametric.
The repeated collisions with the medium particles lead to an integrated momentum transfer of Δp 2 ∼qt, with the momentum diffusion coefficient ofq In two dimensions, the integral over the momentum transfer q ⊥ is dominated by the softest collisions with q ⊥ ∼ m D . Using further that m 2 D ∼ g 2 fΛ from Eq. (21), and reminding the reader that the coupling is dimensionful in 2D with ½g 2 ¼ 1, the momentum broadening coefficient is parametrically of order Energy conservation dictates a relation between the amplitude and the hard scale corresponding to α ¼ 3β. With this, the momentum transport coefficient at a given value of the hard scale readŝ During the nonequilibrium cascade, elastic scatterings push hard scales to harder momenta and the highest momenta reached at time t are given by We expect that, as in 3D [24], the inelastic scattering rate at the hard scale plays an important role as well. An important qualitative feature of inelastic scatterings is that they do not conserve particle number. As a consequence, there is no strong buildup of particle number in the infrared.
While the UV-cascade shares commonalities with the one in 3D, we emphasize one important difference between 2D and 3D. In 3D, elastic scatterings of large (q ⊥ ∼ Λ) and of small (q ⊥ ∼ m D ) momentum transfers both lead to the same scaling exponent β ¼ −1=7. This would also be the case in 2D if hard scatterings with q ⊥ ∼ Λ were dominant. This can be seen from Eq. (28), where the hard scales q ⊥ ∼ Λ give a contribution ∼Λðg 2 fÞ 2 to the estimate (30). If this contribution from the hard scales were the dominant one forq, the resultingq ∼ Λðg 2 fÞ 2 together with Eqs. (31) and (33) would lead to Λ ∼ QðQtÞ 1=7 as in 3D. Such a hard scattering dominance is clearly disfavored both by our numerical results and the above analytical arguments. Instead, collisions with small momentum transfer dominate the momentum diffusion, resulting in a different value for β.
Since the 2D þ sc theory results from a dimensionally reduced three-dimensional theory, and the scalar field is identified with the originally third component of the gauge field ϕ ≡ A 3 , we expect similar arguments to apply for the scalar distribution at high momenta.

V. CONCLUSION
We have studied a far-from-equilibrium attractor in twodimensional highly occupied gauge theories using real-time classical field simulations. We found that these systems exhibit self-similar evolution of the distribution function that corresponds to an energy cascade toward higher momenta and that its scaling properties can be understood using parametric estimates in kinetic theory.
In particular, we have studied the time dependence of equal-time field correlators and that of the hard scale Λ and of the screening scale m D . We have used different observables, including manifestly gauge-invariant measures, to extract the scaling exponents. We found that both the 2D and 2D þ sc theories exhibit self-similar cascades that bring energy toward the UV and are insensitive to the initial conditions. The cascades are characterized by the evolution of the hard scale Λ ∼ t −β , whose time evolution is described by a scaling exponent β ¼ −1=5. Moreover, the Debye scale is extracted from a longitudinally polarized correlator of chromoelectric fields and is shown to decrease with time as m D ∼ t β toward low momenta. While in 2D, we do not have access to a leading-order accurate kinetic theory description, these scaling exponents can be understood in terms of parametric consideration within a kinetic theory setup. A crucial difference to three dimensions is that soft scattering is enhanced compared to the hard scattering. Therefore, unlike in 3D, in order to derive the correct scaling exponents, screening effects have to be taken into account.
These findings are consistent with a description of hard momentum modes in terms of quasiparticle d.o.f. even if a hard-loop theory is insufficient to describe the dynamics of soft modes ∼m D . To learn more about soft dynamics, further numerical studies are required. These include also unequal-time correlators that can be studied numerically with methods that have been developed recently [42] and successfully used for three-dimensional systems [37,43]. We will report on results from such studies in a subsequent work. In addition, it would be interesting to better understand the origin of the observed IR enhanced region of the scalar field correlator.
With regard to heavy-ion collision phenomenology, this attractor might emerge at times that are too late to be reached in a collision, since other phenomena like plasma instabilities can set in earlier [28,[44][45][46][47][48][49]. However, our observation that the evolution of the self-similar attractor can be understood in terms of kinetic estimates is relevant nonetheless for the understanding of the dynamics at early times in heavy-ion collisions: a kinetic description, and thus a description in terms of quasiparticles, can be used to describe two-dimensional plasmas despite the break-down of hard-loop resummations. Hence, we show that kinetic descriptions can be valid already at the early times of the evolution of the two-dimensional Glasma. The exact time when such descriptions become valid can be estimated in further studies.

ACKNOWLEDGMENTS
We are grateful to P. Arnold In Ref. [34], large values for the amplitude n 0 ≳ 1 were used and it was observed that the initial state changes considerably after the gauge fixing procedure was used. This problem is primarily caused by the nonlinear mapping between suðN c Þ algebra elements A k and SUðN c Þ group elements U k . In this work, we circumvent this problem in several ways simultaneously. First of all, we employ smaller initial amplitudes n 0 < 1. Second, we construct the initial link field U k ðxÞ in such a way that its Fourier transformed anti-Hermitian traceless part ½U k ah ðpÞ 8 (and not the logarithm as in Ref. [34]) of the link is constructed to reproduce the desired momentum distribution (4). Thus, here −g jk ∂ j ½U k ah ðxÞ ¼ 0 ðA1Þ is correct to machine precision initially without the need of an additional gauge fixing procedure. This avoids the issue of the 8 This is defined by ½V ah ≡ −i 2 ðV − V † − 1 N c TrðV − V † ÞÞ for an SUðN c Þ matrix V. exponentiation of the algebra element spoiling the transversality of the field that was problematic in the study [34].
Also note that since no fields depend on the x 3 coordinate, the Coulomb gauge condition only applies to k ¼ 1, 2, while the adjoint scalars of theory 2D þ sc do not enter the condition.

APPENDIX B: SCALING EXPONENTS FROM LIKELIHOOD ANALYSIS
To extract the scaling exponents α and β, we employ the self-similarity analysis of Ref. [9] and its modification [14] for the self-similar evolution of the 2D theory displayed in Fig. 2. We define a rescaled distribution f test ðt; pÞ ¼ ðt=t r Þ −α fðt; ðt=t r Þ −β pÞ: By construction, this rescaled distribution is f test ðt r ; pÞ ≡ fðt r ; pÞ for the reference time Qt r ¼ 500. In case of self-similarity, f test ðt; pÞ is time independent. Hence, its difference to the distribution at t r is a good measure of deviation from a self-similar evolution. We can quantify the deviation by computing where χ 2 ðα 0 ; β 0 Þ ≡ χ 2 min takes its minimal value. The normalization N is chosen to satisfy R dαdβWðα; βÞ ¼ 1. We integrate over one of the exponents to obtain an estimate for the distribution of the other exponent, e.g., WðβÞ ¼ R dαWðα; βÞ. We extract an estimate for the uncertainty σ β for every m by fitting the resulting distributions to Gaussian functions ∝ exp½−ðβ − β 0 Þ 2 =ð2σ 2 β Þ. The statistical error σ χ β of the χ 2 fit is estimated as the maximal value of σ β among the different m, giving σ χ β ¼ 0.012, σ χ α ¼ 0.019. We can also extract a systematical error by varying m and requiring that all β 0 values for different values of m lie in the interval ½β 0 − σ sys β ;β 0 þ σ sys β and similarly forα. This leads to the error estimates σ sys β ¼ 0.004 and σ sys α ¼ 0.0035. The statistical χ 2 errors are clearly the larger of these. The mean values and error estimates quoted in Eqs. (14) and (15) are obtained by combining and rounding the mean values and error estimates.