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

TL;DR

This paper investigates the non-equilibrium evolution of highly occupied SU(2) gauge theories in 2+1 dimensions, revealing a universal attractor and specific scaling behaviors relevant to heavy-ion collision initial states.

Contribution

It demonstrates the existence of a universal non-thermal fixed point in 2+1D gauge theories with scalar fields, with detailed scaling exponents and distribution functions derived from simulations.

Findings

Identification of a universal attractor in the evolution of the system.

Determination of scaling exponents $eta = -1/5$ and $eta = 3eta$ for the energy cascade.

Observation of IR enhancement of scalar correlations below the Debye mass.

Abstract

Motivated by the boost-invariant Glasma state in the initial stages in heavy-ion collisions, we perform classical-statistical simulations of SU(2) gauge theory in 2+1 dimensional space-time both with and without a scalar field in the adjoint representation. We show that irrespective of the details of the initial condition, the far-from-equilibrium evolution of these highly occupied systems approaches a unique universal attractor at high momenta that is the same for the gauge and scalar sectors. We extract the scaling exponents and the form of the distribution function close to this non-thermal fixed point. We find that the dynamics are governed by an energy cascade to higher momenta with scaling exponents $\alpha = 3\beta$ and $\beta = -1/5$. We argue that these values can be obtained from parametric estimates within kinetic theory indicating the dominance of small momentum transfer in the scattering processes. We also extract the Debye mass non-perturbatively from a longitudinally polarized correlator and observe an IR enhancement of the scalar correlation function for low momenta below the Debye mass.

Figures

Equations

$$S_{\mathrm{Y}M}[A]=-\frac{1}{4}\int d^{d+1}x{F}_a^{\mu \nu }{F}_{\mu \nu }^a,$$

$$S_{\mathrm{Y}M}^{2D+sc}[A]=S_{\mathrm{Y}M}^{2D}[A]+S_{\varphi }^{2D}[\varphi ]$$

$$S_{\varphi }^{2D}[\varphi ]=-\frac{1}{2}\int d^{2+1}x(D_j^{ab}{\varphi }^b)(D_{ac}^j{\varphi }^c).$$

$$f(t=0,p)=\frac{Q}{g^2}{n}_0{e}^{-\frac{p^2}{2p_0^2}},$$

$$Q\equiv \sqrt[4]{\frac{C{g}^2\epsilon }{d_{\mathrm{pol}}{d}_{\mathrm{A}}}},$$

$$\frac{1}{d_{\mathrm{A}}V}⟨A_T^a(0,\mathbf{p})(A_T^a(0,\mathbf{p}){)}^{\ast }⟩$$

$$\frac{1}{d_{\mathrm{A}}V}⟨E_T^a(0,\mathbf{p})(E_T^a(0,\mathbf{p}){)}^{\ast }⟩$$

$$\frac{1}{d_{\mathrm{A}}V}⟨{\varphi }^a(0,\mathbf{p})({\varphi }^a(0,\mathbf{p}){)}^{\ast }⟩$$

$$\frac{1}{d_{\mathrm{A}}V}⟨{\pi }^a(0,\mathbf{p})({\pi }^a(0,\mathbf{p}){)}^{\ast }⟩$$

$$f_{\mathrm{E}}(t,p):=\frac{⟨E_T{E}_T^{\ast }⟩(t,p)}{\omega (p)}$$

$$f_{\pi }(t,p):=\frac{⟨\pi {\pi }^{\ast }⟩(t,p)}{\omega (p)},$$

$$⟨E_T{E}_T^{\ast }⟩(t,\mathbf{p})=\frac{1}{d_{\mathrm{A}}V}⟨E_T^a(t,\mathbf{p})(E_T^a(t,\mathbf{p}){)}^{\ast }⟩,$$

$$f(t,p)=(Qt{)}^{\alpha }{f}_s\left((Qt{)}^{\beta }p\right).$$

$${\alpha }_{\mathrm{fit}}-3{\beta }_{\mathrm{fit}}$$

$${\beta }_{\mathrm{fit}}$$

$$\alpha =3\beta ,\beta =-\frac{1}{5},$$

$${\Lambda }_E^2(t)$$

$${\Lambda }_{\pi }^2(t)$$

$${\Lambda }_{\mathrm{pert},E/\pi }^2(t)=\int \frac{{\mathrm{d}}^{2}p}{(2\pi {)}^2}\frac{p^3}{Q^3}\frac{g^2{f}_{E/\pi }}{Q},$$

$$⟨E_L{E}_L^{\ast }⟩\approx \frac{A}{1+(p^2/m_D^2{)}^{1+\delta }},$$

$$m_D^2$$

$$\sim d_{\mathrm{pol}}{g}^{2}f\Lambda ,$$

$$m_{D,\mathrm{HL}}^2$$

$$=d_{\mathrm{pol}}{N}_{\mathrm{c}}\int \frac{{\mathrm{d}}^{2}p}{(2\pi {)}^2}\frac{g^{2}f(t,p)}{p},$$

$$m_{D,\mathrm{HL}}^2\sim g^{2}f\Lambda d_{\mathrm{pol}}{N}_{\mathrm{c}}\ln (\Lambda /m_{D,\mathrm{HL}}),$$

$$\frac{\partial f(t,p)}{\partial t}=-C[f](t,p),$$

$$C^{2\leftrightarrow 2}[f_{\mathbf{p}}]$$

$$=\frac{1}{2}{\int }_{\mathbf{k},{\mathbf{p}}^{\prime },{\mathbf{k}}^{\prime }}{\left|M\right|}^2(2\pi {)}^{d+1}{\delta }^{d+1}\left(P+K-P^{\prime }-K^{\prime }\right)$$

$$f_{\mathbf{p}}{f}_{\mathbf{k}}\left(1+f_{{\mathbf{p}}^{\prime }}\right)\left(1+f_{{\mathbf{k}}^{\prime }}\right)-\left(1+f_{\mathbf{p}}\right)\left(1+f_{\mathbf{k}}\right)f_{{\mathbf{p}}^{\prime }}{f}_{{\mathbf{k}}^{\prime }},$$

$$|M{|}^2\sim \frac{g^4}{(q_{\perp }^2+m_D^2{)}^2},$$

$$\frac{\mathrm{d}\Gamma }{\mathrm{d}{q}_{\perp }}\sim \frac{g^4}{(q_{\perp }^2+m_D^2{)}^2}\int {\mathrm{d}}^{2}pf(1+f).$$

$$\hat{q}$$

$$\sim g^4\int \mathrm{d}{q}_{\perp }\frac{q_{\perp }^2}{(q_{\perp }^2+m_D^2{)}^2}\int {\mathrm{d}}^{2}pf(1+f).$$

$$\hat{q}\sim \frac{{\Lambda }^2(g^{2}f{)}^2}{m_D}\sim {\Lambda }^{3/2}(g^{2}f{)}^{3/2}.$$

$$g^{2}f\sim \frac{Q^4}{{\Lambda }^3},$$

$$\hat{q}\sim \frac{Q^6}{{\Lambda }^3}.$$

Full text

††preprint: CERN-TH-2019-115

I Introduction

A characteristic feature of many highly-occupied systems is that they often approach universal self-similar attractors, also referred to as non-thermal fixed points (NTFP) Micha:2004bv ; Berges:2008wm . Examples have been found with classical field methods in various theories in three spatial dimensions (3D) including non-Abelian gauge theories, relativistic and non-relativistic scalar field theories Berges:2008mr ; Berges:2012ev ; Kurkela:2012hp ; Schlichting:2012es ; Lappi:2016ato ; Berges:2013eia ; Berges:2013fga ; Micha:2004bv ; Berges:2008wm ; Nowak:2010tm ; Nowak:2011sk ; Berges:2014bba ; Berges:2015ixa ; Orioli:2015dxa ; Schmied:2018upn ; Mazeliauskas:2018yef , and in two-dimensional scalar systems Gasenzer:2011by ; Schole:2012kt ; Karl:2016wko ; Deng:2018xsk . Non-thermal fixed points have recently also been found experimentally in ultra-cold atom experiments Prufer:2018hto ; Erne:2018gmz . A kinetic theory description of the underlying theory is often a natural way to explain the existence and properties of such fixed points Micha:2004bv ; Arnold:2002zm ; Kurkela:2011ti ; York:2014wja ; Orioli:2015dxa ; Walz:2017ffj ; Chantesana:2018qsb . 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 theories111See, however, Ref. Arnold:2007tr where the soft degrees of freedom are treated as 2D, but the hard ones as 3D.. Differently from the three-dimensional case where an effective kinetic theory has been formulated to leading order accuracy Arnold:2002zm , infrared (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 quasi-particle 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 ultra-relativistic heavy-ion collisions. In this context, non-linear interactions of gluons produced at central rapidities have been argued to lead to a transverse momentum scale $Q_{\mathrm{s}}\gg\Lambda_{\mathrm{QCD}}$ up to which gluonic fields are of order $\mathcal{A}\sim 1/g$ Gelis:2010nm , where $g$ is the gauge coupling. If this saturation scale is sufficiently hard, the system is weakly coupled $\alpha_{s}(Q_{s})\equiv g^{2}/(4\pi)\ll 1$ and consists of highly occupied “Glasma” color fields  Lappi:2006fp . 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 self-similar 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.222We expect the results to carry over to SU(3) theories as well. The qualitative agreement of weakly-coupled SU($N_{\mathrm{c}}$) theories far from equilibrium for $N_{\mathrm{c}}=2,3$ has been observed for different phenomena Ipp:2010uy ; Berges:2017igc ; Berges:2008zt . 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 Lappi:2017ckt , 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 quasi-particle descriptions of soft momentum modes $p\sim 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 will 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.

II Theoretical background

II.1 Theories and initial conditions

We consider non-Abelian $\mathrm{SU}(N_{\mathrm{c}})$ gauge theories with $N_{\mathrm{c}}=2$ in $d=2$ spatial dimensions. The starting point is the classical gauge field action

[TABLE]

with $F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}+g\,f^{abc}\,A^{b}_{\mu}A^{c}_{\nu}$, where repeated color indices $a=1,...,N_{\mathrm{c}}^{2}-1$ and Lorentz indices $\mu,\nu=0,...,d$ imply summation over them. Using the generators $\Gamma^{a}$ of the $\mathrm{su}(N_{\mathrm{c}})$ algebra, the gauge field in fundamental representation reads $A_{\mu}=A^{a}_{\mu}\Gamma^{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

[TABLE]

with an adjoint scalar field $\phi^{a}$ and

[TABLE]

Here the summation is over $j=1,2$ and the covariant derivative is $D^{ab}_{j}=\delta^{ab}\partial_{j}-gf^{abc}A_{j}^{c}$. This theory can be obtained from Yang-Mills theory in 3 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_{3}^{a}\equiv\phi^{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_{\mathrm{Y}M}]=[S_{\mathrm{Y}M}^{2D}]=[S_{\phi}^{2D}]=0$, from which one easily deduces that $[g]=1/2$ and $[A]=[\phi]=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_{\mathrm{pol}}=1$ transverse polarization in 2D, while the 2D+sc theory, originating from a 3D system, has $d_{\mathrm{pol}}=2$ transverse polarizations: one from gauge degrees of freedom and one from adjoint scalars.

The systems are initialized, using the method described in Sec. II.2, with a field configuration that has a chosen single-particle distribution function $f(t,\mathbf{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\ll 1$ but highly occupied $f\gg 1$ initial conditions of the form

[TABLE]

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 Kurkela:2012hp ; Berges:2013fga ). The combination that we have access to in the classical field formulation is the energy density scaled with the coupling $g^{2}\varepsilon$, which has the momentum dimension $[g^{2}\varepsilon]=4$. This allows us to define a constant momentum scale in a gauge-invariant way as

[TABLE]

with the number of transverse polarizations $d_{\mathrm{pol}}$, and the dimension $d_{\mathrm{A}}=N_{\mathrm{c}}^{2}-1$ of the adjoint representation. The constant $C$ is taken as $C=20\sqrt{2\pi}\approx 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\propto\sqrt[3]{\varepsilon}$ that is natural for a scale derived from a 2-dimensional energy density. The scale $Q$ will be used to measure all dimensionful quantities.

II.2 Semi-classical simulations

At high occupation numbers, we can use the classical-statistical approximation to study the dynamical evolution of systems far from equilibrium Aarts:2001yn ; Smit:2002yg . A description of this standard technique can be found, for instance, in Refs. Berges:2013fga ; Boguslavski:2018beu . Then all fields are classical and discretized on a cubic lattice of size $N_{s}^{2}$ 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,\mathbf{x})\approx\exp(ia_{s}gA_{j}(t,\mathbf{x}))$ and chromo-electric fields $E^{j}_{a}=\partial_{t}A_{j}^{a}$. 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 fields333The fields at each momentum $\mathbf{p}=(p_{1},p_{2})$ can be split into a transverse and longitudinal part $A_{j}^{a}(\mathbf{p})=A_{T,j}^{a}(\mathbf{p})+A_{L,j}^{a}(\mathbf{p})$, such that $A_{T}^{a}(\mathbf{p})\cdot\mathbf{p}=0$ while $|A_{L}^{a}(\mathbf{p})\cdot\mathbf{p}|=|A_{L}^{a}(\mathbf{p})|\,p$. at $Qt=0$ follow

[TABLE]

with $\pi^{a}\equiv E^{3}_{a}$, while other combinations for two-point functions vanish, as well as $\langle A\rangle=\langle E\rangle=0$.444In practice, this is achieved by setting $A_{j}^{a}(t=0,\mathbf{p})=\sqrt{f(t=0,p)/p}\;c_{a}(\mathbf{p})\,v_{T,j}(\mathbf{p})$, and similarly for the other fields, with complex-valued Gaussian random numbers $c_{a}(\mathbf{p})$ and the transverse polarization vector $v_{T}(\mathbf{p})$.

Since such initial conditions violate the Gauss law constraint, the latter is restored with the algorithm from Moore:1996qs 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 $\partial_{j}A_{j}=0$ at the measurement time (see also App. A) and Fourier transform the fields according to $A(t,\mathbf{p})=\int d^{d}x\,e^{-i\mathbf{p}\cdot\mathbf{x}}A(t,\mathbf{x})$. A central quantity of interest is the single-particle distribution function $f(t,\mathbf{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 (6) – (9). Unless stated otherwise, our standard definition will be

[TABLE]

where we will set the dispersion $\omega(p)=p$ neglecting soft scale effects since we are primarily interested in the dynamics at high momenta. We also used an abbreviated notation

[TABLE]

and similarly for other correlators.

We note that the classical-statistical framework with the initial correlators (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}\rightarrow 0$ for $g^{2}f(t,p)$ held fixed Aarts:2001yn ; Smit:2002yg ; Berges:2013lsa . 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\sim 1/a_{s}\gg\Lambda$ are non-physical in our case and we use different discretizations to verify that our results are not sensitive to $a_{s}$ (or to the lattice volume).

III Simulation results

III.1 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 non-thermal fixed point at high momenta where the distribution function follows a self-similar evolution

[TABLE]

In order to constitute a universal non-thermal fixed point, the scaling exponents $\alpha$, $\beta$ 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 $\varepsilon$ in the case of an energy cascade to higher momenta Micha:2004bv ; Berges:2013lsa ; York:2014wja , 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 $\varepsilon$.

This attractor property is illustrated in Fig. 1 where we show the gauge distribution $f_{\mathrm{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_{\mathrm{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 times555For $Qt\leq 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\leq 4000$. $Qt=75-16000$. The lower panel shows this same data in rescaled coordinates: the rescaled gauge distribution $(t/t_{r})^{-\alpha}g^{2}f_{\mathrm{E}}/Q$ is shown as a function of rescaled momentum $(t/t_{r})^{\beta}p/Q$. The scaling indices $\alpha$ and $\beta$ have been numerically extracted to produce the best overlap of the distribution functions at the different times employing a least-square fit procedure Berges:2013fga ; Orioli:2015dxa as detailed in App. B, leading to best-fit values

[TABLE]

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:

[TABLE]

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 behaviour.

A similar conclusion can be drawn for the 2D+sc theory. In Fig. 3, the rescaled gauge and scalar distributions $f_{\mathrm{E}}$ and $f_{\pi}$ are shown in the upper and lower panels, respectively.666For $Qt\leq 4000$ a $512^{2}$ lattice with spacing $Qa_{s}=1/8$ has been used, for the later time we used a $256^{2}$ lattice with $Qa_{s}=1/16$. We checked that both discretizations coincide for $Qt\leq 4000$. For comparison, the original curves are depicted in the insets. For the rescaling of amplitudes and momenta, the same exponents $\alpha$ and $\beta$ 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 $\propto(p/Q)^{-\sigma}$ at lower momenta and a steep drop at high momenta. This closely resembles the non-thermal fixed point in 3D theory Berges:2008mr ; Kurkela:2012hp ; Berges:2013eia , 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 $\sigma=1$, which can be seen in the lower panel of Fig. 2, where a power law with $\sigma=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 $\sigma=1$, corresponds to a classical thermal distribution $T_{\mathrm{eff}}/p$ at low momenta with a time-dependent effective temperature $T_{\mathrm{eff}}$ Kurkela:2011ti . 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 York:2014wja . However, numerical classical Yang-Mills simulations in 3D theory have found large corrections to this quantity at early times leading to extractions of $\sigma\approx 1.3$ Berges:2008mr ; Kurkela:2012hp . 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_{\pi}$ are enhanced for low momenta $p\lesssim m_{D}$ roughly below the Debye mass that will be discussed in Sec. III.3. 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 $\mathcal{O}(N)$-symmetric scalar field theories where an IR region has been observed Berges:2008wm ; Gasenzer:2011by or may even be connected to non-trivial topological structures Gasenzer:2013era . 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.

III.2 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 Kurkela:2012hp ; Berges:2013eia

[TABLE]

These provide typical hard momentum scales that are expected to grow as $\Lambda^{2}(t)\sim t^{-2\beta}$ in the self-similar regime. This can be seen from their perturbative expressions

[TABLE]

where all higher powers in the field amplitude were neglected and the Coulomb gauge condition was used. Note that because of $Q^{4}\propto g^{2}\varepsilon$, the hard scales can be interpreted as ratios $\propto\int\mathrm{d}^{2}p\,p^{2}\,\omega(p)f/\int\mathrm{d}^{2}p\;\omega(p)f$, characterizing the momentum scale that dominates the energy density.

The gauge-invariant hard scales $\Lambda^{2}(t)$, rescaled with $t^{2\beta}$, 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 $\Lambda^{2}_{\mathrm{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 $\Lambda^{2}_{E/\pi}\approx\Lambda^{2}_{\mathrm{pert},E/\pi}$. Moreover, hard scales from different sectors and theories agree well $\Lambda^{2}_{E}\approx\Lambda^{2}_{\pi}$. One observes that for $Qt\gtrsim 75$, the rescaled hard scales become approximately constant, indicating $\Lambda^{2}_{E/\pi}(t)/Q^{2}\propto(Qt)^{-2\beta}$. The onset time of self-similar scaling and the value for $\beta$ 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. Lappi:2017ckt . 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).

III.3 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 equal-time correlator $\langle E_{T}E_{T}^{*}\rangle(t,p)\equiv pf_{\mathrm{E}}(t,p)$ and its scalar analogy $\langle\pi\pi^{*}\rangle(t,p)\equiv pf_{\pi}(t,p)$. Let us now focus on the longitudinally polarized equal-time correlation function $\langle E_{L}E_{L}^{*}\rangle(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 $\langle E_{T}E_{T}^{*}\rangle$ in Fig. 6 are flat up to a high momentum $p\sim\Lambda$ beyond which they decrease rapidly, which is, of course, equivalent to our previous observation that $f(t,p)\sim 1/p$ up to a hard scale. Similarly, we have observed the IR enhancement of the $\langle\pi\pi^{*}\rangle$ correlator and its agreement with $\langle E_{T}E_{T}^{*}\rangle$ at higher momenta already at the example of $f_{\pi}$. The longitudinally polarized correlator $\langle E_{L}E_{L}^{*}\rangle$ approaches $\langle E_{T}E_{T}^{*}\rangle$ at the lowest momenta, while strongly differing for momenta above the Debye mass $m_{D}$. Indeed, as known in thermal equilibrium Kurkela:2012hp and also observed far from equilibrium in the 3D theory Boguslavski:2018beu , the longitudinal correlation function follows the form

[TABLE]

for momenta $p\lesssim\Lambda$. In the late-time limit and in thermal equilibrium one then expects $\delta=0$. We have fitted this form to $\langle E_{L}E_{L}^{*}\rangle$ 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$, $\delta$ and $m_{D}$. As expected from our previous discussions, the amplitude quickly approaches an $A\sim t^{\alpha-\beta}$ power law and the deviation $\delta$ monotonously decreases from $\delta\approx 0.2-0.3$ at early times to $\delta\approx 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_{D}^{2}/d_{\mathrm{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_{D}^{2}$ is depicted. This indicates that $m_{D}^{2}$ scales with the number of degrees of freedom $d_{\mathrm{pol}}$, which are $1$ for 2D and $2$ for 2D+sc theory. Moreover, $m_{D}^{2}$ is observed to approach a $t^{2\beta}$ 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

[TABLE]

where we used energy conservation in the last line and where $f$ is the amplitude at hard momenta $p\sim\Lambda$. Since Eq. (III.3) leads to $m_{D}/\Lambda\sim(Qt)^{2\beta}\ll 1$, the scale separation between the soft scale $m_{D}$ and the hard scale $\Lambda$ increases with time. This should in general allow for a perturbative (HL) expansion with the ratio $m_{D}/\Lambda$ as the expansion parameter.

The leading order HL expression for the Debye mass is Laine:2016hma

[TABLE]

where $f_{\pi}\equiv 0$ for the 2D theory and where $f(t,p)$ is an average distribution. Because of $f\sim 1/p$ at low momenta, or even steeper for the scalar distribution, the integral diverges in the IR in 2 spatial dimensions and needs to be regularized by some cutoff at the scale $m_{D}$. This leads to

[TABLE]

bringing a logarithmic correction to the estimate (III.3). The HL expression (III.3) was used in Ref. Lappi:2017ckt 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\sim(Qt)^{\beta}$ with values for $\beta$ 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 non-equilibrium evolution of gauge theories can also be studied using an effective kinetic theory setup. In Arnold:2002zm , an effective kinetic theory has been formulated for $d=3$ spatial dimensions which is defined by a Boltzmann transport equation

[TABLE]

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 UV-cascading system have been well understood in terms of such a kinetic theory description Kurkela:2012hp ; Berges:2012ev ; York:2014wja . This effective kinetic theory describes the evolution of modes at momentum scales well above the screening scale $p\gg 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 amplitudes777Here $|M|^{2}$ is expressed in a non-relativistic normalization related to the usual relativistic normalization by $|M|^{2}=|\mathcal{M}|^{2}/(16\,p\,k\,p^{\prime}\,k^{\prime})$.

$|M|_{\mathrm{vacuum}}^{2}\sim g^{4}/(q_{\perp}^{2})^{2}$ appearing in the elastic part of the collision kernel

[TABLE]

with $\int_{\mathbf{k}}\equiv\int\mathrm{d}^{d}k/(2\pi)^{d}$, $f_{\mathbf{p}}\equiv f(t,\mathbf{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_{\perp}\sim 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 (III.3) in two dimensions gets equal contributions from all momentum scales $m_{D}<p<\Lambda$, 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 non-perturbative way, reminiscent of the ultra-soft magnetic 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. Kurkela:2011ti . 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

[TABLE]

where $q_{\perp}$ 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

[TABLE]

This is a 2-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 1-dimensional transverse momentum $\mathrm{d}q_{\perp}$. The factor $\int\mathrm{d}^{2}pf$ accounts for the number density particles in the medium against whom the collision occurs and $(1+f)\sim f$ is the final state Bose enhancement factor. For a thermal-like infrared spectrum $f\sim T_{*}/p$ with an effective temperature $T_{*}$, the integral over $\mathrm{d}p$ 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\sim\Lambda$ 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 $\Delta p^{2}\sim\hat{q}t$, with the momentum diffusion coefficient of

[TABLE]

In two dimensions the integral over the momentum transfer $q_{\perp}$ is dominated by the softest collisions with $q_{\perp}\sim m_{D}$. Using further that $m_{D}^{2}\sim g^{2}f\,\Lambda$ from Eq. (III.3), and reminding the reader that the coupling is dimensionful in 2D with $[g^{2}]=1$, the momentum broadening coefficient is parametrically of order

[TABLE]

Energy conservation dictates a relation between the amplitude and the hard scale

[TABLE]

corresponding to $\alpha=3\beta$. With this, the momentum transport coefficient at a given value of the hard scale reads

[TABLE]

During the non-equilibrium cascade, elastic scatterings push hard scales to harder momenta and the highest momenta reached at time $t$ are given by

[TABLE]

This equation governs the evolution of the hard scale. Combining it with Eq. (32) gives the time evolution of the hard scale

[TABLE]

corresponding to the value $\beta=-1/5$. These values for the scaling exponents are summarized in Eq. (16) and were used in all the plots of Sec. III. They are consistent with the numerically determined values in Eq. (15).

We expect that, as in 3D Kurkela:2011ti , 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 build-up 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_{\perp}\sim\Lambda$) and of small ($q_{\perp}\sim m_{D}$) momentum transfers both lead to the same scaling exponent $\beta=-1/7$. This would also be the case in 2D if hard scatterings with $q_{\perp}\sim\Lambda$ were dominant. This can be seen from Eq. (28), where the hard scales $q_{\perp}\sim\Lambda$ give a contribution $\sim\Lambda\,(g^{2}f)^{2}$ to the estimate (30). If this contribution from the hard scales were the dominant one for $\hat{q}$, the resulting $\hat{q}\sim\Lambda\,(g^{2}f)^{2}$ together with Eqs. (31) and (33) would lead to $\Lambda\sim 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 $\beta$.

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 $\phi\equiv 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 2-dimensional 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 towards 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 $\Lambda$ 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 towards the UV and are insensitive to the initial conditions. The cascades are characterized by the evolution of the hard scale $\Lambda\sim t^{-\beta}$, whose time evolution is described by a scaling exponent $\beta=-1/5$. Moreover, the Debye scale is extracted from a longitudinally polarized correlator of chromo-electric fields and is shown to decrease with time as $m_{D}\sim t^{\beta}$ towards 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 quasi-particle degrees of freedom even if a Hard-Loop theory is insufficient to describe the dynamics of soft modes $\sim 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 Kurkela:2016mhu and successfully used for 3-dimensional systems Boguslavski:2018beu ; PineiroOrioli:2018hst . 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 Mrowczynski:1994xv ; Arnold:2003rq ; Rebhan:2004ur ; Romatschke:2005pm ; Arnold:2007tr ; Berges:2012iw ; Berges:2012cj . 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 quasi-particles, 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.

Acknowledgements.

We are grateful to P. Arnold, J. Berges, A. Mazeliauskas, A. Rebhan, P. Romatschke, S. Schlichting, M. Strickland and R. Venugopalan for valuable discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No ERC-2015-CoG-681707). 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. K. B. and J. P. would like to thank the CERN Theory group for hospitality during part of this work. The authors wish to acknowledge CSC - IT Center for Science, Finland, for computational resources.

Appendix A Note on initial conditions in Coulomb gauge

In Ref. Lappi:2017ckt , large values for the amplitude $n_{0}\gtrsim 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 $\mathrm{su}(N_{\mathrm{c}})$ algebra elements $A_{k}$ and $\mathrm{SU}(N_{\mathrm{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$. Secondly, we construct the initial link field $U_{k}(\mathbf{x})$ in such a way that its Fourier transformed anti-Hermitian traceless part $[U_{k}]_{\mathrm{ah}}(\mathbf{p})$ 888This is defined by $\left[V\right]_{\mathrm{ah}}\equiv\frac{-i}{2}\left(V-V^{\dagger}-\frac{1}{N_{\mathrm{c}}}\mathrm{Tr}\left(V-V^{\dagger}\right)\right)$ for an $\mathrm{SU}(N_{\mathrm{c}})$ matrix $V$. (and not the logarithm as in Ref. Lappi:2017ckt ) of the link is constructed to reproduce the desired momentum distribution (4). Thus here

[TABLE]

is correct to machine precision initially without the need of an additional gauge fixing procedure. This avoids the issue of the exponentiation of the algebra element spoiling the transversality of the field that was problematic in the study Lappi:2017ckt .

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 $\alpha$ and $\beta$, we employ the self-similarity analysis of Ref. Berges:2013fga and its modification Orioli:2015dxa for the self-similar evolution of the 2D theory displayed in Fig. 2. We define a rescaled distribution

[TABLE]

By construction, this rescaled distribution is $f_{\mathrm{test}}(t_{r},p)\equiv f(t_{r},p)$ for the reference time $Qt_{r}=500$. In case of self-similarity, $f_{\mathrm{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

[TABLE]

with $\Delta f(t_{i},p)=f_{\mathrm{test}}(t_{i},p)-f(t_{r},p)$ and with the exponent of the energy density $\tilde{\alpha}\equiv\alpha-3\beta$. Momentum integrals are performed in the interval $0.2\leq p/Q\leq 5$. The deviations $\chi^{2}_{m}$ are averaged over the test times $Qt_{i}=75,200,1500,4000,16000$ used in Fig. 2 for different moments with $m=2,\dots,5$. For brevity, we will omit the index $m$. We can now define a likelihood function

[TABLE]

where $\chi^{2}(\tilde{\alpha}_{0},\beta_{0})\equiv\chi^{2}_{\mathrm{min}}$ takes its minimal value. The normalization $\mathcal{N}$ is chosen to satisfy $\int\mathrm{d}\tilde{\alpha}\,\mathrm{d}\beta\;W(\tilde{\alpha},\beta)=1$. We integrate over one of the exponents to obtain an estimate for the distribution of the other exponent, e.g., $W(\beta)=\int\mathrm{d}\tilde{\alpha}\,W(\tilde{\alpha},\beta)$. We extract an estimate for the uncertainty $\sigma_{\beta}$ for every $m$ by fitting the resulting distributions to Gaussian functions $\propto\exp[-(\beta-\beta_{0})^{2}/(2\sigma_{\beta}^{2})]$.

The statistical error $\sigma_{\beta}^{\chi}$ of the $\chi^{2}$ fit is estimated as the maximal value of $\sigma_{\beta}$ among the different $m$, giving $\sigma_{\beta}^{\chi}=0.012,\ \sigma_{\tilde{\alpha}}^{\chi}=0.019$. We can also extract a systematical error by varying $m$ and requiring that all $\beta_{0}$ values for different values of $m$ lie in the interval $[\bar{\beta}_{0}-\sigma_{\beta}^{\mathrm{sys}},\bar{\beta}_{0}+\sigma_{\beta}^{\mathrm{sys}}]$ and similarly for $\tilde{\alpha}$. This leads to the error estimates $\sigma_{\beta}^{\mathrm{sys}}=0.004$ and $\sigma_{\tilde{\alpha}}^{\mathrm{sys}}=0.0035$. The statistical $\chi^{2}$ errors are clearly the larger of these. The mean values and error estimates quoted in Eqs. (14), (15) are obtained by combining and rounding the mean values and error estimates.