Understanding the approach to thermalization in a non-Abelian gauge theory

Sayak Guina,∗, Harshit Pandeya , Sayantan Sharmaa

a The Institute of Mathematical Sciences, a CI of Homi Bhabha National Institute, Chennai, 60113, India

Abstract
We measure the maximal Lyapunov exponent λL of physical states in a SU(2) gauge theory both in and out-of-thermal equilibrium
conditions using ab-initio lattice techniques. A non-equilibrium state is realized starting from an over-occupied initial condition for
low momentum, soft gluons whereas the thermal state comprises of strongly interacting soft gluons at temperatures where these
are well separated from the hard momentum modes. A spectra of positive Lyapunov exponents is observed in both these states,
similar to a chaotic dynamical system. From the Kolmogorov-Sinai entropy rate measured in terms of this spectrum, we estimate
arXiv:2501.04397v1 [hep-lat] 8 Jan 2025

the typical time-scale of ∼ 0.50(3) fm/c to achieve thermalization at T ∼ 600 MeV starting from the non-thermal state. We also
measure, for the first time, the λL for long wavelength critical modes of SU(2) using the out-of-time-ordered correlation functions
of a classical Z2 scalar field theory, which shares the same universal behavior near the deconfinement phase transition. The λL is
observed to maximize at the transition temperature.
Keywords: Non-equilibrium quantum field theory, QCD, non-Abelian gauge theory, chaos, out-of-time-ordered correlation

1. Motivation states near to the color-deconfinement phase transition and at

high temperatures, deep in the deconfined phase. Measuring
Understanding the mechanism of thermalization in strongly
the KS entropy rate in a typical thermal state at ∼ 600 MeV
interacting systems is a challenging research problem in theo-
from the Lyapunov exponents, we for the first time provide an
retical physics [1, 2] as it would be able to explain physics at
estimate of the typical thermalization time ∼ 0.50(3) fm/c for
different length scales, e.g. in heavy-ion collisions, cold atom
the system to relax into it, starting from the classical attractor
traps or in the early universe. Starting from a non-equilibrium
regime. This estimate is consistent with the recently reported
state, a necessary condition for eventually reaching a thermal
result [15] for the thermalization time measured in terms of the
state is a mixing between different regions of the phase space. A
evolution of the scale that characterizes the infrared magnetic
mixing in the phase space can occur if there is a mechanical in-
modes of QCD starting from the same non-thermal state. Apart
stability in the system [3], which is characterized by at least one
from reinforcing the stochastic features of non-Abelian gauge
positive Lyapunov exponent. Furthermore in the thermal state
theories [16] already strongly developed in its classical modes,
the entropy or the information content [4] is maximized [5].
our study has a practical relevance as well. It provides an esti-
Except for a few simple systems [6, 7] it is not very well un-
mate of how long it takes to achieve local thermal equilibration
derstood how the entropy evolves as a function of time, starting
i.e. to evolve to a quark gluon plasma phase with initial tem-
from a chaotic non-equilibrium state and eventually maximiz-
peratures 300-600 MeV starting from a highly occupied gluon
ing itself in a thermal state. One important observable is the rate
(glasma) phase believed to be formed in the early stages of a
of change of entropy, known as the Kolmogorov-Sinai (KS) en-
typical heavy-ion collision event [17].
tropy rate [8], which equals the sum of all positive Lyapunov
exponent for Hamiltonian systems [9]. This is known to be true
for dissipative systems that has an attractor solution. Hence cal- 2. Theoretical background
culating the Lyapunov exponents in chaotic dynamical systems
is crucial to understand if, and in what time-scales these would We begin with a discussion about how to measure the (maxi-
eventually thermalize. mal) Lyapunov exponent in a non-Abelian SU(2) gauge theory
In this Letter, we revisit the question whether non-Abelian which is a non-linear quantum field theory. We recall that in
gauge theories are chaotic dynamical systems [10, 11, 12], how- a classical dynamical system, starting from two closely spaced
ever with a goal to understand the phase space mixing and initial coordinates and performing a Hamiltonian evolution if
the process of entropy maximization in an interacting quantum the separation between the two trajectories increases exponen-
field theory. For this purpose, we measure the maximal Lya- tially with time, then the dynamics is defined as chaotic. The
punov exponent of 3 dimensional SU(2) gauge theory using ab- rate at which the trajectories diverge from each other is given
initio lattice techniques in three distinct regimes; for a particu- by the Lyapunov exponent, λL = limt→∞ 1t ln d(0)
d(t)
where d(t) de-
lar non-thermal classical attractor state [13, 14] and for thermal notes the separation between two classical trajectories at time
t. This idea can be easily generalized in a quantum field the-
∗ Corresponding author, email: sayakg@imsc.res.in ory like SU(2), in the semi-classical regime, where the gauge
Preprint submitted to Physics Letters B January 9, 2025
links Uµ (x) are analogous to the space coordinates in a classical 3. The algorithm
configuration space. Using non-perturbative lattice techniques
one can generate ensembles of gauge links using Monte-Carlo 3.1. Studying highly occupied SU(2) gauge theory under non-
algorithms. In thermal equilibrium, two different gauge config- equilibrium conditions
urations which are separated by infinitesimal Monte-Carlo time For a specific realization of a non-equilibrium state of SU(2)
steps can be evolved in real-time using the SU(2) Hamiltonian gauge theory we start from an initial condition where the mo-
in order to estimate the Lyapunov exponent. mentum p modes of the gauge fields are occupied according to
−|p|2
Alternately the Lyapunov exponent can also be motivated by
a phase-space distribution given by, f˜(p) = g2 f (p) = n0 Q|p|s e 2Q2s
studying the following out-of-time ordered correlation (OTOC)
where Q s is the gluon-saturation scale which is typically be-
defined in one-dimensional classical dynamics as,
tween 1-2 GeV and g is the gauge coupling. We consider
!2 Q s = 1.5 GeV in this work. Here and in the subsequent sec-
δx(t) tions we would be considering the theory to be discretized on
{x(t), p(0)} =
2
∼ e2λL t , (1)
δx(0) a 3 dimensional spatial box with N sites along each spatial di-
rection and a is the lattice spacing. The quantity n0 sets the
The OTOCs can be easily generalized for quantum field theo- initial occupation number of the infrared modes of the gauge
ries. For real-valued scalar fields, the expression in the LHS fields which we have chosen such that it is close to the typi-
of Eq. 1 in terms of position and momenta can be replaced by cal equilibrium energy densities at high temperatures T > 600
a commutator of the field operator ϕ̂(x, t) and their conjugate MeV. Since the occupation numbers are large (owing to g being
momentum fields π̂(x, t). However for scalar fields the OTOC small), the gauge links behave classically. For generating these
can be also equivalently defined as [ϕ̂(x, t), ϕ̂(0, 0)]2 [18]. In a classical configurations we solve for the Hamiltonian’s equa-
quantum system, OTOC grows exponentially upto some time tion of motion for the gauge links Ux (t) and its conjugate fields
tE , called the Ehrenfest time [19] beyond which its growth sat- Exib (t) along directions i = 1, 2, 3 and with color components
urates. This is the timescale beyond which the wavefunction b = 1, 2 on the lattice,
spreads over the entire volume of the system and there is a tran-
N3
∆t ∆t ∆t X h b i j
! !
sition from a wave-like behavior to a particle-like behavior of i
the quantum system [20]. However in the semi-classical limit, Exib t + = Exib t − +2 tr iΓ (Vx + c.c) ,
2 2 a j,i
the OTOC will not saturate. Such a scenario is realized when ∆t ∆t
x (t+ 2 )Γ
Eib b
scalar fields have large occupation numbers, where one can re- Uxi (t + ∆t) = ei a Uix (t) , (5)
place the commutator with the following expression,
where Γb are generators of SU(2) and Vxi j is the spatial staple
!2 defined in the i- j spatial plane as the derivative of the gauge
δϕ(x, t)
C(x, t) = ⟨{ϕ(x, t), ϕ(0, 0)} ⟩ = ⟨
2
⟩. (2) plaquette U P defined using path ordering P in color space as,
δπ(0, 0)
δU P  
Vxi j = , U ij
= P U i
U j
U U j†
x x+î x+ĵ x .
i†
(6)
We henceforth denote the classical fields ϕ(x, t) as ϕx (t) and δUxi P

similarly for its conjugate fields. If indeed the system is chaotic,

the (maximal) Lyapunov exponent λL can be extracted from the We use a symplectic leap-frog integrator with a time-step ∆t =
OTOC using the relation, 0.0125Q s for the time evolution in order to respect time reversal
invariance. We study different lattice volumes with typical ex-
1 tent along the spatial directions N = 128, with spacings aQ s =
λL = lim ln C(x = 0, t) . (3) 1.115, 0.9, 0.8 such that the spatial length of theh box in physicali
t→∞ 2t
units is ∼ 16 fm. The Gauss law ∆ta i triΓa Uxi (t) − Ux+ i
P
î
(t)
An OTOC typically spreads ballistically within the light cone was implemented with a precession of 10 at t = 0 and it was
−15
along rays of constant speed v = x/t along the x-t plane with a checked to remain so at later times. With such a choice of ini-
scaling exponent ν [21] as, tial condition the system eventually enters a self-similar scaling
" !ν # regime at late times where the gluon distribution function takes
v the form f˜(p, t) = (Q s t)−4/7 f s ((Q s t)−1/7 p), where f s is a homo-
C(x, t) ∼ e2λ(v)t , λ(v) = λL 1 − . (4)
vB geneous function describing a static distribution [22]. Within
such a regime the hard (ultraviolet), electric and the magnetic
From the spectrum of Lyapunov exponents λ(v) in the veloc- momentum scales, the later two corresponding to scales set by
ity space of propagating OTOCs, one can define the butterfly the Debye mass and the square root of the spatial string tension
velocity vB from the relation λ(vB ) = 0. The quantities λL , vB respectively are well separated similar to a thermal plasma [23].
and exponent ν are all, in general, temperature dependent. In In order to measure the degree of chaoticity in such a sys-
Section 4 we will extract λL , vB for SU(2) gauge theory under tem in a gauge-invariant manner, we follow the procedure as
thermal equilibrium for a wide range of temperatures as well as outlined in Ref. [15]. First we implement an infinitesimal
in a particular non-equilibrium state. The ab-initio lattice algo- shift in the initial condition to n0 + δn0 where δn0 = 10−5
rithms used for this purpose are detailed in the next section. and perform a Hamiltonian evolution of the gauge links Ux′ (t)
2
as a function of time. We next measure the separation be- 3.3. Chaotic dynamics near deconfinement phase transition
tween two gauge trajectories at a time t starting with these Near the deconfinement phase transition temperature T c the
two infinitesimally close initial conditions through the gauge- hard, electric and the magnetic scales coincide, hence the effec-
invariant distance measure defined in terms of the spatial pla- tive thermal field theory described in the preceding section can
quettes d(t) = 1/NP P |trU P (t) − trU P′ (t)|/2, where NP are the
P
no longer describe the gauge field dynamics. Hence we change
number of independent plaquettes. We particularly choose the our strategy to estimate the Lyapunov exponent. Since SU(2)
time interval to be within the self-similar scaling regime, in or- gauge theory shares the same universality class as the 3D Ising
der to extract the Lyapunov exponent. model, hence its long-wavelength or low energy critical modes
modes will undergo the same dynamics as a Z2 scalar field. We
3.2. Studying thermal SU(2) plasma at high temperatures
thus simulate a relativistic Z2 scalar field theory to understand
At finite temperatures, apart from the hard scale πT present in
the chaotic dynamics of the critical modes near the deconfine-
the problem, there are two new scales that are generated dynam-
ment phase transition. We start with the Hamiltonian which on
ically; the electric scale gT and the magnetic scale g2 T/π [24].
the lattice is defined as,
When these three scales are well-separated, one can integrate
out the hard and electric gluons, with momenta ≳ gT to obtain
 
X  π2 1 X m 2 2
a
!
λ s

an effective Hamiltonian [25] which describes the dynamics of H= x
 − ϕx ϕy + + 3 ϕ2x + ϕ4x  . (8)
2 2 y∼x 2 4!
the magnetic (soft) gluons where the electric fields Ẽxi satisfy x
the equation of motion P
Here y∼x implies the sum over all nearest neighbors y around
−∂t Ẽxi + [D j , F̃ (x)] =
ji
σẼxi + ζ̃xi (t) , (7) the lattice site x. In our work, we set the values of m2 a2 =
where σ is the color conductivity and ζ̃xi (t) are the −1 , λ s = 1. We first generate thermal ensembles using the
Langevin algorithm, in which we consider the time evolution of
stochastic color-force fields which satisfy ⟨ζ̃xic1 (t1 )ζ̃xjb2 (t2 )⟩ =
the scalar and its conjugate momentum fields through equations
2T σδi j δcb δ3 (x1 − x2 )δ(t1 − t2 ) in accordance to the fluctuation-
of motion,
dissipation relation. Discretizing the Eq. 7, the electric fields,
noise fields and color-conductivity can be written in dimen- ∂H
∂t ϕx (t) = = πx (t) , (9)
sionless units on the lattice in terms of Exi = gẼxi a2 , ζxi = ∂πx
a3 gζ̃xi (t) and σa = σT .T a respectively. We set σ/T to its per- ∂H
∂t πx (t) = − γπx (t) + ηx (t) 2γT .
p
turbative values [26] where the coupling g at each tempera- − (10)
∂ϕx
ture T is determined from the 2-loop β function. We numer-
ically implement Eq. 7 in dimensionless units which is simi- where γ is the damping term and ηx are Gaussian noise fields
lar to the electric field evolution in Eq. 5 but with additional defined on each lattice site with a zero mean value which sat-
contributions due to the noise and isfies ⟨ηx (t)ηy (t′ )⟩ = δxy δ(t − t′ ). The relation between the
i damping terms given by
damping factor γ and the coefficient of the noise term due
h
−∆t/a aσExib (t − ∆t) − ζxib (t − ∆t) . The evolution in Eq. 7 is
similar to a Langevin equation of magnetic modes of the gauge to thermal fluctuation can be obtained using the fluctuation-
fields in presence of a heat bath where the noise term represents dissipation theorem. To initialize the algorithm we select ini-
the effects of the hard modes which forms the heat-bath and tial values of ϕx (0) and πx (0) on each lattice site x. We then
the damping dynamically brings the system to a thermal equi- numerically evolve the above differential equations using the
librium. This classical system consists of 6N 3 oscillators of symplectic leap-frog algorithm with a time step ∆t = 0.01. Af-
energy aT and has an energy density on the lattice [12] given ter some Langevin time tL , we expect the system to thermalize
since the Langevin algorithm is ergodic and follows the condi-
by N63 k |k| aT |k| , where the sum is over all allowed lattice mo-
P
menta k, which we have also verified in our calculations. On tion of detailed balance. Once we achieve thermalization using
the other hand, the energy density of quantum SU(2) gauge the- the Langevin algorithm, we use the thermal configuration thus
ory at temperatures ∼ 2 times than the deconfinement temper- obtained, switch off the noise and damping terms and perform
ature is close to its Stefan-Boltzmann limit π2 T 4 /5 [27]. The a real-time classical evolution of the scalar fields. We also ad-
lattice spacing in the effective theory is set in physical units ditionally select infinitesimal perturbation about the initial ther-
from the criterion that the measured energy density matches mal configurations denoted by δϕx (0) = 0 and δπx (0) = rδx0 ,
with the quantum theory at each T . This results in a condi- where r is a Gaussian random number between 0-0.01. Since
tion T.a = (30/π2 )1/3 , which we use to set the lattice spacing in δπx (t) = ∂t δϕx (t) one can re-express the initial momentum pro-
physical units. Once the gauge links are thermalized at a tem- file in the configuration space at x ≡ (0, 0, 0), t = 0 and then
perature T , there is an onset of a plateau in the values of the evolve the fluctuations of the scalar field through
plaquette term as a function of Langevin time. We then con- λs 2
sider two thermal configurations infinitesimally separated as a ∂2t δϕx (t) = ∇2 δϕx (t) − m2 δϕx (t) − ϕ (t) δϕx (t) . (11)
2 x
function of Langevin time, switch off the noise and the damp-
ing terms in the Hamiltonian and evolve them subsequently as In the evolution equation Eq. 11 one requires the values of the
a function of real-time. We then calculate the gauge invariant scalar fields ϕx (t) which are obtained through a classical Hamil-
distance measure introduced in the earlier section to calculate tonian evolution of initial thermalized values of ϕx (0) at each
the (maximal) Lyapunov exponent. temperature T , generated earlier using the Langevin algorithm.
3
 0.009 λL/Tc (T/Tc)-2.8(3)

0.007
T/Tc=0.907
25 =0.96
6.8(7) =1
(T/Tc) ln(C(x=0,t))
=1.02
15

0.005 5
-5 t.Tc

0 2500 5000

0.88 0.93 0.98 1.03 1.08

T/Tc

Figure 1: The Lyapunov exponent as a function of T/T c for a Z2 scalar field

Figure 2: The temperature dependence of the diffusion coefficient D near T c
theory which from universality represent the critical long-distance modes of a
extracted using Eq. 4. The inset shows the ballistic spread of the OTOC function
S U(2) thermal state near the deconfinement temperature T c The inset shows
of Z2 scalar fields in x-t plane at T c . The color profile from green to red denotes
the growth of OTOC functions with time t.T c .
increasing values of OTOC in the x-t plane.

We further evolve the initial momentum field fluctuations δπx (t) fly velocity through the relation D = v2B /λL [29]. As evident
through ∂t δϕx (t) = δπx (t). We thus have all ingredients to calcu- from the figure, the diffusion coefficient is constant below T c
late the OTOC function in semiclassical limit using Eq. 2, from but decreases with increasing temperatures as (T/T c )−3.7(2) for
which we extract a positive Lyapunov exponent, discussed in T > Tc.
the next section. At higher temperatures, T ≫ T c , the gauge invariant dis-
tance measure d(t) = 2N1 p P |trU P (t) − trU P′ (t)| were measured
P

4. Results on lattice of N = 32 spatial sites and for different temperatures

results for which are shown in the inset of Fig 3. The initial
We first measure the OTOC functions for Z2 scalar field the- d(t) ∼ 10−15 measures the typical difference between the pla-
ory at different temperatures close to the phase transition, which quettes in two thermal configurations which, when evolved in
is shown in the inset of Fig 1. After an initial dip, The OTOCs time grows exponentially eventually saturating at late times to
exhibit an exponential rise as a function of time. By perform- ∼ 10−2 , since the trace of plaquette is bounded to unity and
ing a fit to the OTOCs at late times t.T c > 1000 with an expo- the gauge group is compact. The lattice size are chosen large
nential ansatz, we extract the (maximal) Lyapunov exponent at enough for this purpose such that the values are ∼ 15(6) fm at
different temperatures, results of which are shown in Figure 1. temperatures 0.6(1.5) GeV respectively. Extracting the (maxi-
Here a.T c = 9.37 in dimensionless unit and agrees very well mal) Lyapunov exponent λL from the exponentially rising part
with the values in the literature [28]. The Lyapunov exponent of d(t) with time, we study its dependence on temperature, re-
thus extracted varies as λL ∼ T 6.8(7) below T c , whereas in the sults of which are shown in Fig. 3. The λL varies linearly with
range (1, 1.05)T c its temperature variation can as characterized temperature between 0.6-3 GeV with a slope λL /T ∼ 0.52
as T −2.8(3) , which was determined after performing a fit to our which satisfies the conjectured bound λL ≤ 2πT in quantum
data. The deconfinement phase transition temperature in SU(2) many-body chaotic systems [30] and also satisfies λL ∼ g2 E/6,
is T c ∼ 300 MeV [27]. where E is the average energy per plaquette [11]. The diffusion
Further in the inset of Fig 2 we show the heat-map of the coefficient extracted from λL thus vary inversely with temper-
OTOC at T = T c in a two dimensional plane denoted by the ature as (T/T c )−0.95(3) shown in Fig 4. At these temperatures
time and one of the spatial coordinates, performing a sum over the distance measure spreads from an initial tiny fluctuation at
the other spatial coordinates. The OTOC spreads ballistically t.T = 0 ballistically over the space-time and grows exponen-
within the light cone from an initial small perturbation at the tially within the light-cone, evident from the inset of the same
origin as evident from the color gradient of the heat-map. In figure, the color gradient from green to red indicates increasing
general the λL is a function of the velocity and it can be mea- values of the distance measure d(t). Across these high temper-
sured from the OTOC function using Eq. 4. As evident in Fig 1, atures, the butterfly velocity is vB = 0.8c and does not vary
the Lyapunov exponent is maximum along v = 0 i.e. along the significantly. Compared to its value at T c , the butterfly velocity
time direction. We next calculate the butterfly velocity vB by reduces by only ∼ 10% which implies that once thermalized,
fitting the OTOC C(x, t) to the ansatz given in Eq. 4. The vB most of the configuration (color) space gets occupied.
varies across the phase transition and is found to be maximum We next extract the Lyapunov exponent within the self-
at T = T c . In Fig 2, we have shown the diffusion coefficient similar scaling regime of the non-equilibrium SU(2) plasma as
as a function of temperature, which is related to the butter- outlined in Sec. 3.1. In order to compare with a typical ther-
4
 0.7
5.5
60
λL/Tc

t.T
4.5 30
0.5
0.999(3)
3.5 (T/Tc) 0
-20 -10 0 10 20
-0.95(3)
2.5 -5 D.Tc α (T/Tc) x.T
ln(d(t))
0.3
1.5 -15 T/Tc=4
=5
=6
0.5 40 Qs.t 80
0.1
0 2 4 6 8 10 12 2 4 6 8 10
T/Tc T/Tc

Figure 3: The Lyapunov exponent as a function of T/T c for a SU(2) gauge Figure 4: The temperature dependence of the diffusion coefficient D in SU(2)
theory at high temperatures in thermal equilibrium shown as triangles. These gauge theory at high temperatures T ≫ T c . The inset shows the ballistic spread
are compared with the Lyapunov exponents in a non-thermal state of SU(2) of the Lyapunov exponents in a typical thermal state of SU(2) at T = 1.2 GeV.
with comparable energy densities (circles). The inset shows the growth of the The color gradient from green to red denotes increasing values of OTOC in the
gauge invariant distance measure with time t.T c for the thermal SU(2) state for x-t plane.
3 different temperatures.

is not particularly related to the details of microscopic interac-

mal state at temperature T ∼ 2 GeV, we set the initial density
tions between the degrees of freedom but to the symmetries of
of gluons to be n0 = 16, and consider lattice box of spatial ex-
the system.
tent ∼ 16 fm such that the energy densities are similar in the
thermal as well in the non-thermal scaling regime. The (maxi- The diffusion coefficient D which measures how fast the cor-
mal) Lyapunov exponent λL = 0.66 Q s is similar in magnitude relation between the fields increases with time, has many inter-
to that in a thermal state at T ≫ T c . We have also shown the esting features as a function of temperature. The λL at T c shown
values of the maximal Lyapunov exponents in Fig 3 in this non- in Fig. 1 correspond to the long-wavelength critical modes of
thermal state, as a function of an effective temperature which is SU(2) gauge theory which are a subset of the infrared mag-
the fourth root of the energy density, as blue data points. The netic modes, whose momenta are bounded by |p| ≤ g2 T/π.
data points lie on the same curve determining the temperature The fact that D is independent of temperature below and at T c
dependence of the thermal Lyapunov exponents. For lower en- and decreases sharply above T c implies that these critical gluon
√ modes at the phase transition are maximally spread in the con-
ergy densities ε(∼ n0 ≲ 1) it was earlier conjectured [31] and
later observed [15] from classical-statistical lattice calculations figuration space. This spread is an inherent chaotic phenomena,
in the self-similar non-thermal scaling regime of SU(3), that which also manifests itself in the stochastic properties [37] of
the (maximal) Lyapunov exponent λL varies as ε1/4 . Our cal- the eigenvalue density of Wilson loop operator in the limit of
culations were performed in SU(2) non-thermal scaling state large number of colors. In this limit, the quantum fluctuations
√ are sub-dominant, and an order-disorder transition [38, 39] in
at larger gluon energy densities (∼ n0 ≳ 4) where we ob-
serve a linear scaling of λL with ε. The linear scaling of λL terms of the area of Wilson loops has similar features as the
at higher energies similar to a thermal plasma can be under- (de)confinement transition.
stood in terms of linear temperature dependence of the plasmon The temperature dependence of D at T ≫ T c in Fig. 4 can be
damping rate [32]. interpreted from the fact that D represents a diffusion in the con-
figuration space of color gauge links and hence can be related to
the inverse of color conductivity σ. The σ calculated perturba-
5. Physical implications of our results tively in a thermal non-Abelian plasma, increases linearly with
temperature [26]. In the effective dynamics of the soft (mag-
Our results demonstrate that the infrared (soft) modes of the netic) modes given in Eq. 7, the σ represent how these modes
SU(2) have chaotic dynamics both in thermal equilibrium at rearrange in response to the ultra-relativistic hard modes. The
a wide range of temperatures 0.9-10 T c as well in a particu- soft modes receive random kicks from the thermal bath of hard
lar non-thermal fixed point as evident from the positive (max- modes and a finite σ effectively dampens them, allowing for an
imal) Lyapunov exponent λL . Interestingly the λL decreases efficient thermalization.
near critical temperature T c with temperatures T > T c , maxi- Next we address the implications of our study for understand-
mizing at the deconfinement phase transition. Such a property ing the mechanism of thermalization in gauge theories. When
of λL for the critical modes has been observed earlier in spin we discuss about thermalization, to be more precise, we imply
systems [33, 34, 35, 36] as well as in scalar field theories in how fast the magnetic (soft) modes equilibrate. Starting from
2 dimensions [18]. The temperature dependence of λL thus, a non-equilibrium state the entropy increases until it reaches its
5
maximum in a thermal state. The KS entropy density rate ṡKS modes during the entire evolution of the system. Interestingly
can thus be used to estimate the time required in this process. our estimate of the thermalization time is consistent with the
From the sum of all the positive Lyapunov exponents λi [9], one early switch-on time of hydrodynamics to describe flow observ-
can estimate this quantity on the lattice as ables in heavy-ion collisions [43] but is smaller than estimates
(≳ 2.5 fm/c) based on perturbative scattering process [44].
λi a
P
ṡKS a4 = − i 3 ∼ −4λL .a , (12) In our approach we can only predict the time taken by the
N soft modes to thermalize if the initial non-thermal and the final
for the range of temperatures 2-10 T c , considering that about thermal states are known. However understanding the micro-
1/3 of 18N 3 number of the Lyapunov exponents are posi- scopic mechanism that drives the system away from the self-
tive [12] and their values can be obtained from the relation similar non-thermal scaling regime would require a more de-
λ(v) ∼ λL [1 − (v/vB )2 ] evident in our data, where we have tailed understanding of how hard particles are created in this
taken the values of v as randomly distributed between 0 and vB . regime and how their interactions with the soft modes can be
We have verified that our estimate does not vary significantly if formulated within an effective theory. Possible mechanisms in-
the values of v are chosen to be uniformly distributed between clude the onset of plasma instabilities [45, 46, 47, 48, 49, 50,
[0, vB ]. Starting from a thermal state ṡKS measures the rate at 13, 14, 51], and the process of thermalization is an inherently
which information about the initial state of the system is lost. non-perturbative phenomenon [52]. We would like to address
We use the thermodynamic relation sth T = ε + p to obtain the these topics in a future work.
entropy density of a thermal state of SU(2) to be sth a3 = 43 ε.a
4
T.a
at temperatures T ≫ T c where the energy density and pressure 7. Acknowledgements
are related by ε ≃ 3p. On the other hand the entropy density
for the non-thermal state inR the self-similar We are grateful to Sumilan Banerjee, Rajiv Gavai and So-
iregime described in
eren Schlichting for discussions related to this work. The au-
h
section 3.1, snon-th a = − f (|p|) ln f (|p|) d3 (ap) ≃ 4.8 can
3 ˜ ˜
be calculated by performing a numerical integration using the thors acknowledge the computing time provided by the High
non-thermal distribution f obtained in our computations. Us- Performance Computing Center at the Institute of Mathemati-
ing the entropy densities of the thermal and non-thermal states cal Sciences.
one can derive the thermalization time tth = 0.8/λL (T ), by in-
tegrating over Eq. 12. Thus starting from a highly occupied References
gluon state in a non-thermal scaling regime, a thermal state at
T ≃ 600 MeV can be reached within a time, tth ∼ 0.50(3) fm/c, [1] J. Berges, Introduction to nonequilibrium quantum field theory, AIP Conf.
Proc. 739 (1) (2004) 3–62. arXiv:hep-ph/0409233.
whereas it would take tth ∼ 0.70(5) fm/c to reach a thermal state [2] J. Berges, M. P. Heller, A. Mazeliauskas, R. Venugopalan, QCD ther-
at T ≃ 450 MeV. malization: Ab initio approaches and interdisciplinary connections, Rev.
Mod. Phys. 93 (3) (2021) 035003. arXiv:2005.12299.
[3] N. Krylov, J. Migdal, Works on the Foundations of Statistical Physics,
6. Outlook Princeton Series in Physics, Princeton University Press, 2014.
[4] C. E. Shannon, A mathematical theory of communication, The Bell Sys-
tem Technical Journal 27 (3) (1948) 379–423.
In this letter, we provide a detailed understanding of the dy- [5] E. T. Jaynes, Information theory and statistical mechanics, Phys. Rev. 106
namical properties of soft (magnetic) modes of non-Abelian (1957) 620–630.
SU(2) gauge theory as a function of temperature as well as [6] V. Latora, M. Baranger, Kolmogorov-sinai entropy rate versus physical
in a particular non-thermal attractor state. These magnetic entropy, Phys. Rev. Lett. 82 (1999) 520–523.
[7] M. Dzugutov, E. Aurell, A. Vulpiani, Universal relation between the
modes interact non-perturbatively even at asymptotically high kolmogorov-sinai entropy and the thermodynamical entropy in simple liq-
temperatures [40] and influence dynamical properties like the uids, Phys. Rev. Lett. 81 (1998) 1762–1765.
sphaleron transition rate [25, 41], and other transport properties [8] A. Kolmogorov, Viscosity, black holes, and quantum field theory, Dokl.
of the non-Abelian plasma [42]. We show here that these modes Akad. Nauk SSSR 119 (1958) 861.
[9] J. B. Pesin, Ljapunov characteristic exponents and ergodic properties of
play an important role in the process of thermalization in gauge smooth dynamical systems with an invariant measure, in: Hamiltonian
theories. In a high temperature plasma, where a clean separa- Dynamical Systems, CRC Press, 2020, pp. 512–515.
tion between the hard (ultraviolet momentum) modes and these [10] G. Savvidy, The yang-mills classical mechanics as a kolmogorov k-
soft modes is possible, resulting in an effective stochastic de- system, Physics Letters B 130 (5) (1983) 303–307.
[11] B. Müller, A. Trayanov, Deterministic chaos in non-abelian lattice gauge
scription of the latter [25], we demonstrate the chaotic behavior theory, Physical review letters 68 (23) (1992) 3387.
inherent in the system by extracting the spectra of positive Lya- [12] C. Gong, Lyapunov spectra in su (2) lattice gauge theory, Physical Review
punov exponents. This allows us to measure the Kolmogorov- D 49 (5) (1994) 2642.
Sinai entropy rate which determines how fast the entropy flows [13] J. Berges, K. Boguslavski, S. Schlichting, R. Venugopalan, Turbulent
thermalization process in heavy-ion collisions at ultrarelativistic energies,
into the system eventually attaining a maximum. We have used Phys. Rev. D 89 (7) (2014) 074011. arXiv:1303.5650.
this rate to measure a thermalization time ∼ 0.7-0.5 fm/c to [14] J. Berges, K. Boguslavski, S. Schlichting, R. Venugopalan, Universal at-
reach to a thermal state at 450-600 MeV starting from a partic- tractor in a highly occupied non-Abelian plasma, Phys. Rev. D 89 (11)
(2014) 114007. arXiv:1311.3005.
ular non-thermal state of soft gluons which shows a particular [15] H. Pandey, R. Shanker, S. Sharma, Understanding the approach to ther-
self-similar behavior [14]. The only assumption that goes into malization from the eigenspectrum of non-Abelian gauge theories (7
our calculation is the conservation of energy density of the soft 2024). arXiv:2407.09253.

6
[16] L. Ebner, B. Müller, A. Schäfer, L. Schmotzer, C. Seidl, X. Yao, Entangle- [44] R. Baier, A. H. Mueller, D. Schiff, D. T. Son, Does parton saturation at
ment properties of su (2) gauge theory, arXiv preprint arXiv:2411.04550 high density explain hadron multiplicities at RHIC?, Phys. Lett. B 539
(2024). (2002) 46–52. arXiv:hep-ph/0204211.
[17] A. Andronic, P. Braun-Munzinger, K. Redlich, J. Stachel, Decoding the [45] R. Baier, A. H. Mueller, D. Schiff, D. T. Son, ’Bottom up’ thermalization
phase structure of QCD via particle production at high energy, Nature in heavy ion collisions, Phys. Lett. B 502 (2001) 51–58. arXiv:hep-ph/
561 (7723) (2018) 321–330. arXiv:1710.09425. 0009237.
[18] A. Schuckert, M. Knap, Many-body chaos near a thermal phase transition, [46] P. B. Arnold, J. Lenaghan, G. D. Moore, QCD plasma instabilities
SciPost Phys. 7 (2) (2019) 022. arXiv:1905.00904. and bottom up thermalization, JHEP 08 (2003) 002. arXiv:hep-ph/
[19] J. Rammensee, J. D. Urbina, K. Richter, Many-body quantum interfer- 0307325.
ence and the saturation of out-of-time-order correlators, Physical Review [47] P. Romatschke, M. Strickland, Collective modes of an anisotropic quark
Letters 121 (12) (2018) 124101. gluon plasma, Phys. Rev. D 68 (2003) 036004. arXiv:hep-ph/
[20] K. Hashimoto, K. Murata, R. Yoshii, Out-of-time-order correlators in 0304092.
quantum mechanics, Journal of High Energy Physics 2017 (10) (2017) [48] A. Rebhan, P. Romatschke, M. Strickland, Hard-loop dynamics of non-
1–31. Abelian plasma instabilities, Phys. Rev. Lett. 94 (2005) 102303. arXiv:
[21] V. Khemani, D. A. Huse, A. Nahum, Velocity-dependent lyapunov ex- hep-ph/0412016.
ponents in many-body quantum, semiclassical, and classical chaos, Phys. [49] P. B. Arnold, J. Lenaghan, G. D. Moore, L. G. Yaffe, Apparent thermal-
Rev. B 98 (2018) 144304. ization due to plasma instabilities in quark-gluon plasma, Phys. Rev. Lett.
[22] J. Berges, K. Boguslavski, S. Schlichting, R. Venugopalan, Universal 94 (2005) 072302. arXiv:nucl-th/0409068.
attractor in a highly occupied non-abelian plasma, Physical Review D [50] A. Kurkela, G. D. Moore, Bjorken Flow, Plasma Instabilities, and Ther-
89 (11) (2014) 114007. malization, JHEP 11 (2011) 120. arXiv:1108.4684.
[23] J. Berges, K. Boguslavski, L. de Bruin, T. Butler, J. M. Pawlowski, Order [51] S. Schlichting, S. Sharma, Chiral Instabilities and the Fate of Chirality
parameters for gauge invariant condensation far from equilibrium, Phys. Imbalance in Non-Abelian Plasmas, Phys. Rev. Lett. 131 (10) (2023)
Rev. D 109 (11) (2024) 114011. arXiv:2307.13669. 102303. arXiv:2211.11365.
[24] M. Laine, A. Vuorinen, Basics of Thermal Field Theory, Vol. 925, [52] Y. V. Kovchegov, Can thermalization in heavy ion collisions be de-
Springer, 2016. arXiv:1701.01554. scribed by QCD diagrams?, Nucl. Phys. A 762 (2005) 298–325. arXiv:
[25] D. Bodeker, On the effective dynamics of soft nonAbelian gauge fields at hep-ph/0503038.
finite temperature, Phys. Lett. B 426 (1998) 351–360. arXiv:hep-ph/
9801430.
[26] P. B. Arnold, L. G. Yaffe, High temperature color conductivity at next-
to-leading log order, Phys. Rev. D 62 (2000) 125014. arXiv:hep-ph/
9912306.
[27] J. Fingberg, U. M. Heller, F. Karsch, Scaling and asymptotic scaling in
the SU(2) gauge theory, Nucl. Phys. B 392 (1993) 493–517. arXiv:
hep-lat/9208012.
[28] D. Schweitzer, S. Schlichting, L. von Smekal, Spectral functions and dy-
namic critical behavior of relativistic z2 theories, Nuclear Physics B 960
(2020) 115165.
[29] M. Blake, Universal charge diffusion and the butterfly effect in holo-
graphic theories, Physical review letters 117 (9) (2016) 091601.
[30] J. Maldacena, S. H. Shenker, D. Stanford, A bound on chaos, Journal of
High Energy Physics 2016 (8) (2016) 1–17.
[31] B. Chirikov, D. Shepelyanski, Stochastic oscillations of classical yang-
mills, JETP Lett 34 (4) (1981).
[32] T. Biró, C. Gong, B. Müller, Lyapunov exponent and plasmon damping
rate in non-abelian gauge theories, Physical Review D 52 (2) (1995) 1260.
[33] T. Bilitewski, S. Bhattacharjee, R. Moessner, Temperature dependence
of the butterfly effect in a classical many-body system, Physical review
letters 121 (25) (2018) 250602.
[34] T. Bilitewski, S. Bhattacharjee, R. Moessner, Classical many-body chaos
with and without quasiparticles, Physical Review B 103 (17) (2021)
174302.
[35] S. Ruidas, S. Banerjee, Many-body chaos and anomalous diffusion across
thermal phase transitions in two dimensions, SciPost Physics 11 (5)
(2021) 087.
[36] P. Butera, G. Caravati, Phase transitions and lyapunov characteristic ex-
ponents, Physical Review A 36 (2) (1987) 962.
[37] J.-P. Blaizot, M. A. Nowak, Large N(c) confinement and turbulence, Phys.
Rev. Lett. 101 (2008) 102001. arXiv:0801.1859.
[38] B. Durhuus, P. Olesen, The Spectral Density for Two-dimensional Con-
tinuum QCD, Nucl. Phys. B 184 (1981) 461–475.
[39] R. Narayanan, H. Neuberger, Universality of large N phase transitions in
Wilson loop operators in two and three dimensions, JHEP 12 (2007) 066.
arXiv:0711.4551.
[40] S. Tah, et. al., The spatial string tension and its effects on screening cor-
relators in a thermal QCD plasma (2025). arXiv:tobepublished.
[41] G. D. Moore, The Sphaleron rate: Bodeker’s leading log, Nucl. Phys. B
568 (2000) 367–404. arXiv:hep-ph/9810313.
[42] H. B. Meyer, Transport Properties of the Quark-Gluon Plasma: A Lattice
QCD Perspective, Eur. Phys. J. A 47 (2011) 86. arXiv:1104.3708.
[43] U. W. Heinz, Thermalization at RHIC, AIP Conf. Proc. 739 (1) (2004)
163–180. arXiv:nucl-th/0407067.