Holographic Kolmogorov-Sinai entropy and the quantum Lyapunov spectrum

In classical chaotic systems the entropy, averaged over initial phase space distributions, follows a universal behavior. While approaching thermal equilibrium it passes through a stage where it grows linearly, while the growth rate, the Kolmogorov-Sinai entropy (rate), is given by the sum over all positive Lyapunov exponents. A natural question is whether a similar relation is valid for quantum systems. We argue that the Maldacena-Shenker-Stanford bound on quantum Lyapunov exponents implies that the upper bound on the growth rate of the entropy, averaged over states in Hilbert space that evolve towards a thermal state with temperature T, should be given by πT times the thermal state’s von Neumann entropy. Strongly coupled, large N theories with black hole duals should saturate the bound. To test this we study a large number of isotropization processes of random, spatially homogeneous, far from equilibrium initial states in large N, N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 Super Yang Mills theory at strong coupling and compute the ensemble averaged growth rate of the dual black hole’s apparent horizon area. We find both an analogous behavior as in classical chaotic systems and numerical evidence that the conjectured bound on averaged entropy growth is saturated granted that the Lyapunov exponents are degenerate and given by λi = ±2πT. This fits to the behavior of classical systems with plus/minus symmetric Lyapunov spectra, a symmetry which implies the validity of Liouville’s theorem.


Introduction
A quantum mechanical description of chaotic many body systems is of interest for a multitude of research areas in physics, especially in the context of condensed matter physics, heavy ion physics, thermalization and quantum information theory. In classical physics the question of "how chaotic" a system is, can be quantified by examining the rate with which phase space trajectories X i (t), with initial separation δX j (0), diverge from one another. In chaotic systems the distance between paths in phase space δX j (t) grows exponentially and the singular values of the matrix δX i (t)/δX j (0) grow or decrease as e λt . In the late time limit t → ∞ the exponents represented by λ are referred to as Lyapunov exponents. In quantum theories this behavior is encoded in out-of-time-order correlators (OTOCs) [8]. The quantum Lyapunov exponents can be extracted from the exponential growth rate of OTOCs ∼ e λ OTOC t at late times 1 T t, where T is the asymptotic temperature of the system, defined by its average energy density.
In recent years OTOCs and these exponents, encoding the speed with which quantum systems scramble information, received much attention [9][10][11][12][13][14][15], especially after it was shown that there exists an upper bound [1] for λ OTOC given by λ OTOC ≤ 2πT . This implies an upper bound on the speed of the development of quantum chaos. Systems that are holographic duals to Einstein gravity have been found to saturate this bound [1,5,16,17]. However, to the best of our knowledge, there exists no formal proof as to exactly which quantum systems show this behavior for which operators. Therefore, we focus an the best established holographic system, namely AdS 5 /CFT 4 , where this bound is known to be saturated.
For classical chaotic systems the Kolmogorov-Sinai entropy provides information about the Lyapunov spectrum at time scales that are relevant for thermalization or dissipation t ∼ 1 T . In order to avoid confusion we will henceforth refer to the Kolmogorov-Sinai JHEP01(2022)165 entropy as the Kolmogorov-Sinai entropy rate, since (contradicting its name often found in the literature) it is not an entropy, but rather an entropy growth rate. The heuristic idea behind the relation between entropy growth rate and the Lyapunov spectrum is the following: starting from some initial ensemble of configurations in phase space, whose evolution is described by a chaotic, Hamiltonian system, the volume of this ensemble spreads throughout phase space while branching out and evolving towards a fractal shape. More and more phase space cells (used to evaluate the coarse grained entropy S) are required to cover its shape, while Liouville's theorem ensures that its total volume stays unchanged. Those directions in phase space for which the Lyapunov exponents are positive, will contribute to the growth of the number of needed cells. Thus they contribute to (coarse grained) entropy growth, 1 such that naively we have In the time interval during which entropy grows linearly, dS dt = S KS is referred to as the Kolmogorov-Sinai entropy rate. Empirically we know that this relation between entropy growth and Lyapunov exponents of chaotic systems is only correct with further specifications. For a thermalizing system this statement can clearly only be true for some intermediate time period before the system reaches thermal equilibrium. Moreover, the right hand side of the above equation is independent of initial conditions, while the time derivative of the growing entropy will depend on the initial state. Thus, in general only the ensemble averaged entropy, where we average over a large ensemble of initial phase space configurations that are far from equilibrium, will allow us to determine the sum over all positive Lyapunov exponents, as demonstrated in [21].
In [3] the authors speculated that the Maldacena-Shenker-Stanford (MSS) bound implies an upper bound on entropy growth. This bound should be saturated for conformal theories with Einstein gravity duals in the bulk, which implies that black holes are also the fastest entropy generators with dS/dt = i,λ i >0 2πT . In this contribution we propose a slightly modified version of this conjecture, guided by observations made in classical statistical mechanics [21]. We are going to argue that strongly coupled, large N theories, holographically dual to Einstein gravity, fulfill where · denotes the Hilbert space average over states that initially are far from equilibrium 2 and evolve towards the same thermal state with temperature T . Relation (1.2) is supposed to hold until thermal equilibrium is almost reached, at which point dS dt gradually decreases to 0, in analogy to classical chaotic systems. We examine this numerically by studying, via holography, far-from-equilibrium istoropization in N = 4 Super Yang Mills JHEP01(2022)165 theory (SYM). We determine the number of Lyapunov exponents from the number of degrees of freedom of the equilibrated black hole, which is taken to be (albeit in all generality not proven to be) the same as its Bekenstein-Hawking entropy S eq . For a (semi-)classical 3 Yang-Mills theory the Lyapunov spectrum computed on the lattice can be split into three parts of which each belongs to one third of the degrees of freedom. 4 One third of the Lyapunov exponents λ + i is positive and their sum is the Kolmogorov-Sinai entropy rate, one third is negative with λ − i = −λ + i , and the remaining third corresponds to the unphysical degrees of freedom (e.g. longitudinal polarizations), which have zero Lyapunov exponent, see section 3. Thus, the classical phase space volume is constant as long as we don't smear or coarse-grain. As argued in [25] any measurement provides such a coarse graining due to the quantum mechanical uncertainty relation and thus leads to net entropy growth. For the field theory part of the AdS 5 /CFT dual it thus depends crucially on the precise question which is asked and how entropy is defined, whether the latter grows or not, i.e. whether information gets lost or not. In this contribution we focus on just one specific detail of this highly complex topic which can be clarified numerically.

Quantum chaos, Lyapunov exponents and the Kolmogorov-Sinai entropy rate
Chaos, information scrambling and operator growth in quantum theories can all be studied with the help of OTOCs of general hermitian operators V and W separated by time t: where · T is the thermal expectation value at temperature T . By studying the growth rate of this quantity at times before the Ehrenfest time but well after the dissipation time, we can quantify "how chaotic" a quantum system is via the exponent λ OTOC in In [1] it has been famously shown that λ OTOC , given in natural units, is bounded from above by 2π times the temperature T . In addition, large N conformal field theories (CFTs) dual to Einstein gravity are known to saturate this bound In general, the exponent λ OTOC as defined in eq. (2.2) corresponds to the largest Lyapunov exponent of an entire Lyapunov spectrum (which may be extracted by considering OTOCs of suitable operators) with which the operators V and W overlap. For classical physics Liouville's theorem forces the Lyapunov spectrum to be symmetric, implying that for every

JHEP01(2022)165
λ i there is a λ j with λ j = −λ i . It is worth noting at this point that, besides the fact that it would appear natural, there is no proof that this symmetry holds for a quantum version of Lyapunov exponents. For chaotic quantum systems the authors of [1] determined the exponent λ OT OC defined in eq. (2.3) via an auxiliary function F (t), which decreases with the same exponential rate as C(t) increases. For large N conformal field theories with gravity duals the function F at times 1 T t can be written as [1,5,16,17] where f 0 and f 1 are positive order O(1) constants depending on the choices for V and W . While the operators V (0) and W (0) are hermitian operators, which can be described as a sum of products, which contain only O(1) degrees of freedom. The thermal one point functions of V and W should vanish. The relation (2.4) suggests that in this case every positive Lyapunov exponent is maximal λ = 2πT . Moreover, at first glance it also appears that eq. (2.4) demands that all Lyapunov exponents are positive and that the symmetric structure of the classical Lyapunov spectrum is lost. However, this point is quite subtle. The question is not whether the commutator of generic operators W , V has overlap with an exponentially growing mode, but whether there exist very specific operators, the commutator of which has zero overlap with any such mode. It could have happened that the modest assumptions for V and W that led to eq. (2.4) implied overlap with at least one mode with positive Lyapunov exponent. The idea behind this paper is that, while it might be impossible to extract negative Lyapunov exponents from eq. (2.4), in practice it is still possible to decide whether they exist by determining the fraction of modes which have positive Lyapunov exponents. In the N → ∞ and strong coupling limit equation eq. (2.4) suggests that all positive Lyapunov exponents are expected to be equal, such that, granted eq. (1.2) holds, the ensemble averaged total entropy growth rate is uniquely determined by this fraction. The total entropy growth rate, however, can be calculated as growth rate of the apparent horizon 5 on the gravity side of the duality. Ordering these thoughts, the following statement can be made: if one is completely agnostic as to what to expect one can imagine three possible scenarios, namely • The Lyapunov spectrum keeps the plus/minus symmetry of classical, chaotic systems. It is given by ±2πT , and the relation between Lyapunov exponents and ensemble averaged entropy growth naturally generalizes from statistical mechanics to quantum systems as described by eq. (1.2). Then, the averaged entropy growth rate dS/dt is πT times the number of physical degrees of freedom N DOF (i.e. the total number of Lyapunov exponents).
• Eq. (1.2) is correct, however there are no modes with negative Lyapunov exponent. Then dS/dt is equal to 2πT · N DOF . 5 The reason why we focus on the apparent horizon instead of the event horizon is twofold. On the one hand as argued in 6.6 of [40] the holographic dual of the coarse grained entropy S(t) is proportional to the area of the apparent horizon at time t in infalling coordinates, not the event horizon. On the other hand, despite the highly symmetric setup we consider, it is numerically still easier to determine the apparent horizon, allowing us to compute large ensemble averages. • Either not all Lyapunov exponents associated with physical degrees of freedom have the maximal absolute value 2πT , or the fraction of negative Lyapunov exponents is some number other than 0 or 1/2. Then dS/dt /(πT N DOF ) could be any number between 0 and 2. (It could also happen that eq. (1.2) is not the correct relation between quantum Lyapunov exponents and the quantum Kolmogorov-Sinai entropy rate (i.e. state-averaged entropy growth for the quantum system). In this case we could, of course, not make any statement about dS/dt /(πT N DOF ).) We cannot hope to obtain an indisputable answer, since dS/dt fixes only the sum of all Lyapunov exponents and not their size distribution. However, given what we know about the symmetric shape of the classical Lyapunov spectrum (see next section) and its relation to the classical Kolmogorov-Sinai entropy rate, finding dS/dt ≈ πT · N DOF would be a strong hint that the quantum Lyapunov spectrum is also plus/minus symmetric.

The Lyapunov spectrum of classical SU(N ) Yang Mills theory
As explained in the last section, the central questions to be answered are whether the symmetry between positive and negative Lyapunov exponents persists, i.e. whether Liouville's theorem stays valid, and whether all positive Lypunov exponents approach λ max = 2πT , assuming that the proposed relation (1.2) is correct. For these questions some intuition can be gained from studying classical Yang-Mills theories. This is such a natural thing to do and in fact was already done such a long time ago that we do not feel competent to decide who investigated this question first. Instead, we cite the review [18]. Earlier work JHEP01(2022)165 can be found there. The rational motivating the study of classical Yang-Mills theory is that many examples demonstrate that if a classical theory is chaotic the quantized theory is so too, and that many fundamental properties are related (A typical examples are scars for quantum billiards). For classical Yang Mills theory it was shown by numerical studies that to high accuracy SU(2) even fullfills the criteria for a globally hyperbolic (Anosov) system [19]. These criteria concern the dependence of the uncertainty of the numerically obtained Kolmogorov-Sinai entropy rate on system size and sampling time, a topic we will address below when discussing the precision of our results. These numerical simulations were made by solving the classical Hamilton equations on a finite three dimensional grid, for which one has only a finite number of degrees of freedom, such that it is possible to determine all Lyapunov exponents. Typical results are show in figures 9, 10, and 11 of [18]. One third of the Lyapunov exponents is positive, one third negative, with the same distribution of absolute values, and one third is zero. The latter is due to the fact that a spin 1 field has three degrees of freedom, but for a massless gauge field only two of these are physical while the third is a gauge degree of freedom. These figures show also that even for very small systems, reaching numerically the asymptotic limit, in which all Lyapunov exponents of gauge degrees of freedom are really zero requires very long simulation time.
For systems with finite energy density which equilibrate in finite time, such long fitting windows cannot be realized, and the "intermediate" Lyapunov exponents of the gauge degrees of freedom are numerically still non-zero, see figures 1 and 3 in [20]. This shows that the Kolmogorov-Sinai entropy rate cannot be calculated exactly for finite energy density, which was explored in detail in [19]. In this contribution we will numerically determine the Kolmogorov-Sinai entropy rate from the holographic dual, by analyzing the time dependence of the apparent horizon. However, we expect that also in the dual picture the length of time until saturation is reached is relevant and that the precision with which the slope of the growing apparent horizon area can be determined depends on this time. Let us comment on a feature of the results obtained in [20], which might otherwise be confusing for a careful reader of that paper: in these numerical calculations space was discretized. The lattice spacing, δ, was tuned to a finite value to obtain the energy density of the quantum theory also for the classical theory (rather than infinity), i.e. the continuum limit δ → 0 was not taken. This finite discretisation implied that spatial derivatives were substituted by quotients of differences, i.e. spatial derivatives became non-local, leading in turn to a violation of local gauge symmetry. This artefacts resulted in the Lyapunov exponents of the gauge degrees of freedom becoming non-zero. This effect is barely visible in figure 3 of [20]. To illustrate it more clearly we also performed such an analysis of classical SU(2) and SU(4) theory, see figure 2. (Time derivatives were not affected which lead in addition to different effects for the F 0j (electrical) and F ij (magnetic) components of the field strength tensor.) The symmetry between positive and negative Lyapunov exponents, which is the only point relevant for our discussion, is obviously not affected by this artefact. The bottom line of the discussion in this section is that (semi-) classical 6 YM calculations imply that JHEP01(2022)165 it would at the very least be an unexpected feature of quantum Lyapunov spectra, if the ± degeneracy of the spectrum of Lyapunov exponents wouldn't be observed there, too.

Thermalization and holographic isotropization
We can simulate an isotropizing, initially far from equilibrium, strongly coupled SYM plasma, by using the dual gravitational description. Following the pioneering work of [27,28], there is a large amount of literature on the topic of studying out of equilibrium SYM plasmas via numerical holography (see for example [29][30][31][32][33][34][35]), which we cannot do full justice here. Reference [29] contains a detailed, pedagogical description of the type of calculation we perform in the following. In this section we are going to briefly review the most important points before continuing to perform ensemble averages over a multitude of isotropization processes.
The authors of [28] studied the numerical evolution of an anisotropic initial state in the CFT, produced by a time dependent shear deformation of the metric coupling to the stress energy tensor of the CFT. The holographic duality relates states produced by a time dependent, four dimensional metric h αβ in the large N , large 't Hooft coupling CFT to solutions of five dimensional, classical AdS gravity in the bulk with the time dependent 4D metric as asymptotic boundary. Moreover, the AdS dictionary allows us to determine the expectation value of the CFT stress energy tensor from the bulk metric where g µν is given in Eddington-Finkelstein coordinates (4.2) and g (2) µν represents the second order coefficient of the bulk metric's expansion around the boundary ρ = 0. Assuming JHEP01(2022)165 spatial homogeneity, this metric ansatz reads with det(ĝ ij ) = 1. Near the boundary ρ → 0 one has A(ρ, t) ∼ 1 2ρ 2 , Σ(ρ, t) ∼ 1 ρ ,ĝ ij (ρ, t) ∼ δ ij . For simplicity we follow [29], where the action on the state by the time dependent, arbitrary shear deformation of the boundary metric is replaced by an arbitrary choice of the anisotropy function B(ρ, t) on the initial Cauchy surface t = 0, where B(ρ, t) is given viâ With the ansatz (4.2) the Einstein equations can be written as a nested system of differential equations on null slices [29]. In the case of spatial homogeneity, the following data on time slice t T 00 (t) ,ĝ(ρ, t) (4.4) is sufficient to uniquely solve the system of ordinary differential equations with vanishing spatial gradients, which is most conveniently done using spectral methods [24]. The equations of motion of the boundary stress energy tensor together with knowledge ofĝ(ρ, t), ∂ tĝ (ρ, t) − ρ 2 A(ρ, t)∂ ρĝ (ρ, t) and A(ρ, t) (the latter two of which are functions we solved the nested system of equations for) allow us to compute the data (4.4) on the next time slice. We apply the fourth order Runge Kutta method to obtain a solution to the Einstein equations in the bulk. In the setting we consider, one has 7 T 00 (t) = T 00 (0) = 3 8 N 2 π 2 T 4 and a flat boundary metric. The time dependent pressure components of the stress energy tensor are anisotropic in the beginning and relax towards the equilibrium value N 2 π 2 T 4 /8 on time scales t 1/T . A constant energy density and a non trivial, arbitrary radial dependence of the bulk functions B(ρ, 0) (arbitrary up to an appropriate near boundary behavior B(ρ, 0) ∼ ρ 4 for ρ → 0) on the initial time slice, can be seen as the result of an appropriate, time dependent, spatially homogeneous deformation of the boundary metric with compact support restricted to t < 0, such that for 0 < t the 7 This restricts the region in Hilbert space from which we draw the configurations over which we ensemble average. However, this point could also be made, if we performed a similar calculation as in [29]. The most general case of an arbitrary time and space dependent boundary metric deformation is numerically not practical, as we are interested in large ensemble averages. The entropy growth rate for a single configuration at early times can be split into the Kolmogorov-Sinai entropy rate plus some initial state dependent term dS/dt = S KS + dSi/dt, where S KS is constant until thermal equilibrium is approached, while dSi/dt will in general not be constant and cancels after ensemble averaging. Our ensemble average does yield a constant growth rate for the entropy density s, while individual samples do not exhibit this property of the entropy density (see figure 4). Nonetheless, in future works it might be interesting to further test our numerical results, by computing (if necessary smaller) ensemble averages of more complicated cases with a similar, albeit non-stochastic version of the calculations and numerics described in [26].

JHEP01(2022)165
asymptotic boundary is Minkowski space time without deformations and T 00 is constant. As in [29] we use the radial shift invariance of the metric ansatz (4.2) to keep the apparent horizon at a constant ρ h πT = 1, where T refers to the equilibrium temperature. We use the scaling symmetry ρ → αρ and x i → αx i to set ρ h = 1 throughout our numerical simulations. Physical quantities are given in units of (4.6)

Results
Our aim is to compute the Kolmogorov-Sinai entropy for N = 4 SYM in the large N limit using holography. We consider an ensemble of out of equilibrium, spatially homogeneous but anisotropic states, that evolve towards the same thermal state, in strongly coupled N = 4 SYM. Different initial anisotropic states correspond to different choices for the anisotropy function B introduced in section 4. Using the holographic principle, we can determine the number of physical degrees of freedom (DOF) per unit volume of the SYM plasma. 8 Once the plasma has reached thermal equilibrium the dual description is a Schwarzschild black hole geometry and its number of DOF N DOF is usually taken to be its area measured in Planck length cubed. 9 Thus N DOF can be computed from the Bekenstein-Hawking entropy of the black hole. The number of DOF per unit volume is This implies that the Kolmogorov-Sinai entropy rate density s KS = S KS /V should be given by Where c in (5.2) is equal to 2, if every Lyapunov exponent is positive and maximal, and equal to 1, if every positive Lyapunov exponent is maximal and the Lyapunov spectrum stays plus/minus symmetric. This is assuming that the averaged entropy growth rate of black holes actually saturates the upper bound, derived from the upper bound on quantum Lyapunov exponents. We can compute s KS via holography by ensemble averaging over the growth rates of the apparent horizons' volume element. The ensemble consists of numerical simulations of the isotropizing SYM plasma, described in section 4, with different initial anisotropy functions. Let g(ρ h , t) be the determinant of the metric induced on the apparent horizon on timeslice t, then the CFT entropy density grows as   where · is the state average described in the introduction, section 1, and · {φ} is the average over the ensemble {φ}. 11 We generate multiple ensembles by choosing φ ∈ C ∞ ([0, 1]) randomly with B(ρ, 0) = ρ 4 φ(ρ).

Ensemble elements per ensemble
The set of different ensembles we consider can be classified into two main categories: on the one hand we choose φ(ρ) from finite dimensional subspaces of the function space C ∞ ([0, 1]) generated by polynomials and Gaussians with random coefficients. On the other hand we generate random points for each element of an equidistant grid on [0, 1], interpolate and 10 With the constraint that we are only interested in states that evolve to the same thermal state. 11 Rigorously showing that an ensemble is 'good' in the sense described above appears close to impossible.
Necessary requirements on ensembles include that known Hilbert space averages should be matched by the approximation · {φ} to good accuracy. Therefore, we checked that for the ensembles we studied Tij {φ} matches the equilibrium value N 2 π 2 T 4 /8 up to a small error. For instance in the case of ensemble (IV) we find a maximal deviation from the equilibrium value of 1.6%. 3 3 as a function of time µt, which corresponds to the ensemble (II) averaged CFT entropy density in units of N 2 2π µ 3 with cut-off parameter a = 3. For the ensemble average (red curve) we obtain a constant slope in the interval µt ∈ [0.05, 0.5]. The grey dashed curves show the corresponding plots for a selection of single simulations in our ensemble. Our results only depend negligibly on the cut-off a (see figure 6) as long as a is chosen within a ∈ [0, 5], such that contributions with very large initial slopes, which skip stage (2), are suppressed. smooth the resulting function. We smooth via filtering out large radial derivatives in order to improve numerical stability. In total we collected data from 5 different ensembles ranging in size from several hundred thousand to 3 million simulations. Ensemble (I) was generated by choosing N p random real numbers r i between −15 and 15, while N p itself is a random integer between 5 and 199. The resulting list of points {i/N p , r i } is filtered with a Gauss-filter of width 3/N p . The random function φ is then found by fitting the filtered list with a polynomial of order 15. Ensemble (II) was generated from initial anisotropy functions of the form

JHEP01(2022)165
with β 1,2 ∈ [−10, 10], w 1,2 ∈ [−5, 5], ρ 1,2 ∈ [−0.5, 0.5], a 0,1 ∈ [−4, 4], a 2 ∈ [−1, 1] drawn from uniform distributions. Ensemble (III) is generated analogously to ensemble (II), but with 3 instead of 2 Gaussians, where the corresponding parameters ρ 3 , β 3 and w 3 have the same range as ρ 1,2 , β 1,2 and w 1,2 . For ensemble (IV) we chose N p random points between −5 and 5, where N p is again a random integer between 5 and 99. We then apply a low pass filter onto the list of random points and transform directly from the equidistant grid to a Chebyshev grid via spectral methods. Finally ensemble (V) is generated analogously to (I), but now also the order of the interpolating polynomial is randomly chosen between 5 and 25 and the Gauss-filter has width 4/N p . In each case we represent φ(ρ) as a vector of 26 values on a Chebyshev grid and solve the system of differential equations using spectral methods. In figure 5 we show the density plots for the evolution of the rescaled apparent horizon volume element for ensemble (II) and (V). Figure 5. The density plots for the evolution of the rescaled apparent horizon volume element Σ 3 (ρ h , t) 5 of ensemble (II) and (V) are shown above. We chose the cutting parameter a to be 5 in both cases (Note that the result (5.7) is given for a cutting parameter of a = 0). Besides the density of curves (thin, blue curves), their average (bold, orange curve) and median (purple crosses) we display the densities' maximum on each time slice (black dots). Interestingly for ensembles generated from basis functions on ρ ∈ [0, 1] with random coefficients (ensemble (II) and (III)) median, maxima and average are close and appear to have the same slope. For Ensembles (I), (IV) and (IV), which are generated from interpolated random points on [0, 1], maxima and average differ noticeably. For both density plots there is a small percentage of curves present, which at some time 0 < µt have a slope that is larger than the cutoff a = 5 at µt = 0. The contribution of these curves is small (they make up ∼ 2%). In general we observe a larger variance of the horizon area for ensembles with initial conditions generated via smoothed random points compared with ensembles, whose initial conditions are given by sums of Gaussians and polynomials with random coefficients.

JHEP01(2022)165
For classical, chaotic systems we know what one might expect for the behavior of the ensemble averaged entropy. There, the entropy S(t) follows a general pattern (see [21]). In the first, short stage (1) the behavior of S(t) is dominated by the initial distributions and no general statement can be made. 12 In the second stage (2) S(t) grows approximately linearly JHEP01(2022)165 Figure 6. Here we show the constant slope during the linear growth phase of the ensemble averaged volume element of the apparent horizon Σ(ρ h , t) 3 a , the dual of the CFT entropy, in units of N 2 2π µ 4 as a function of the cutting threshold a. The subscript a in · a indicates that the ensemble average is taken over all histories, for which d dµt Σ(ρ h , t) 3 t=0 does not exceed a. The results displayed above are computed for averages over simulations from our 5 different ensembles described in the text, which all start far from equilibrium, i.e. Σ(ρ h , t = 0) 3 ≈ 0.1 (filled symbols) or Σ(ρ h , t = 0) 3 ≈ 0.01 (empty symbols). For each ensemble and each value of a we determine a linear fit for the averaged entropy growth. The results of the slopes of those fits are displayed above. The average of our results is very close to 1 (central dashed line) with a variance of ≈ 0.16 at a = 0. For very large initial growth rates a the entropy jumps close to equilibrium within a very short initial time span t 1/(aµ). Stage (2) is skipped in this case and linear growth cannot be observed. Thus, for a > 5 some ensembles averages (mainly those starting at Σ(ρ h , t = 0) 3 ≈ 0.1) get 'spoiled' by contributions to the average similar to those depicted in figure 7. and the growth rate corresponds to the Kolmogorov-Sinai entropy rate, i.e. the sum over all positive Lyapunov exponents. Then in the third stage (3) the entropy tends asymptotically towards its equilibrium value. 13 However, for individual runs with inconvenient choices of initial conditions (i.e. initial configurations that already start close to thermal equilibrium, or for which the entropy grows so fast in the first stage that the system is brought close to equilibrium already there) stage (1) and stage (3) might merge, skipping stage (2) in which we are interested. To avoid these pathological contributions to our ensemble average, we consider initial conditions which are far from equilibrium (we both consider averages of runs for which Σ 3 (ρ h , t = 0) = 0.1 ± 0.01 and Σ 3 (ρ h , t = 0) = 0.01 ± 0.01) and focus on initial configurations for which the entropy density at the starting time t = 0 doesn't grow faster than the threshold N 2 π 3 T 4 a/2. We display our results as a function of the threshold or cut-off parameter a in figure 6. We indicate averages over runs for which the slope of the entropy density ∂ t s| t=0 at initial time t = 0 does not exceed 14 N 2 π 3 T 4 a/2 with a subscript a by · a . We find that the growth rate of the ensemble averaged entropy density during 13 In our case thermal equilibrium is synonymous with Σ(ρ h , t) = 1.
14 Or, put differently, for which d dµt Σ 3 (ρ h , t) t=0 < a. Here we exclusively average over those simulations with very large initial growth rates of the entropy density, specifically for which 10 ≤ d dµt Σ(ρ h , t)| t=0 . The initial slope in this example is so large, that the system reaches near-equilibrium before the linear growth phase can start. The results shown above correspond to ensemble (II). the time period in which it grows linearly does not depend on the cut-off parameter a for a wide variety of different cut-off choices a ∈ [0, 5]. For large cut-off values a > 5 pathological contributions to the ensemble average, which skip stage (2), have non-negligible influence on some ensemble averages. See caption of figure 6 for more details. We find

JHEP01(2022)165
The large error is due to our ignorance of which type of ensembles gives the best approximation to the average over all states in Hilbert space that evolve towards the same thermal state with temperature T . Nonetheless it is intriguing that we find a result of approximately 1 in (5.7), which neatly fits to the physical intuition, that the ensemble averaged entropy growth rate for AdS/CFT saturates the theoretical maximal value (1.2), that the Lyapunov spectrum of N = 4 SYM inherits the ±-symmetric structure of the classical YM theory (see section 3) and that in the holographic limit all positive Lyapunov exponents are maximal. The heuristic explanation for the symmetry of the Lyapunov spectrum of classical YM theory is its time reversal invariance. This is equivalent to a constant microscopically resolved, fine grained entropy, i.e. a constant phase space volume, which implies that for every direction in phase space, in which the phase space volume grows with rate e λt there has to be another direction in which it contracts with rate e −λt . Thus, one could make a point that any reasonable 15 generalization of the classical Lyapunov spectrum to a quantum Lyapunov spectrum should inherit this symmetry as long as the quantum theory is unitary.

JHEP01(2022)165 6 Conclusion
By ensemble averaging over a multitude of isotropization processes we found that the Kolmogorov-Sinai entropy density rate of large N , N = 4 SYM at strong coupling is given by We argued that the two most plausible shapes of the Lyapunov spectrum of strongly coupled, large N CFTs with Einstein gravity duals (such that the MSS-bound is saturated) are either all Lyapunov exponents being positive and maximal λ = 2πT , or a degenerate spectrum λ = ±2πT , such that λ = 0, with the Lyapunov spectrum keeping the ± symmetry that we are used to in the case of classical YM theory, or classical physics in general. In the case of a degenerate spectrum, the result (6.1) implies that the intuition, that large N , N = 4 SYM at strong coupling has the largest possible Kolmogorov-Sinai entropy rate fulfilling (during the linear growth phase of the ensemble averaged entropy) 2πT, (6.2) appears to be correct. 16 One interesting statement derived from AdS/CFT is that the quark gluon plasma, produced during heavy ion collisions, thermalizes very quickly [29] on time scales that are just a fraction of 1 fm/c. Even when finite 't Hooft coupling corrections are taken into account [36,37] or non-trivial transverse fluctuations of the energy density are considered [38,39], both of which roughly doubling the thermalization time, one still ends up with a result below 1 fm/c, which does not contradict experimental observations, but rather estimates from weakly coupled, N = 3 YM-calculations [20]. Granted that QCD at high temperatures strongly resembles large N , N = 4 SYM, which according to (6.1) is likely to actually saturate the possible upper bound on the (ensemble averaged) entropy production rate, this mismatch between weak coupling results on the one side and phenomenology and holography on the other side is not surprising.
In future works we will further test the results obtained in this paper by considering simulations with non-homogeneous initial conditions and arbitrary boundary metrics. Other interesting questions to explore in this context are on the one hand, how the Kolmogorov-Sinai entropy rate behaves at finite coupling. This can be done by either including Gauss-Bonnet coupling corrections as in [35] and [42], or tackle the more attractive, but challenging case of α 3 corrections. On the other hand, it would be interesting to clarify how/whether the situation changes, when we replace the entropy computed via the apparent horizon area by the entanglement entropy of some boundary region, weighted by the measure of this boundary area (for spatially homogeneous, anisotropic geometries, entanglement entropy has already been computed in [41]). Again we expect to see a linear growth rate for the (ensemble averaged) entanglement entropy that is given by the sum JHEP01(2022)165 over all positive Lyapunov exponents. This expectation is both motivated by this work and by [43], where the authors found a proof for the relation between the entanglement entropy growth and (classical) Lyapunov exponents for unstable quadratic Hamiltonian describing a bosonic system.