Quasi-Normal Modes from Non-Commutative Matrix Dynamics

We explore the connection between the process of relaxation in the BMN matrix model and the physics of black holes in AdS/CFT. Focusing on Dyson-fluid solutions of the matrix model, we perform numerical simulations of the real time dynamics of the system. By quenching the equilibrium distribution we study the quasi-normal oscillations of scalar single trace observables, we isolate the lowest quasi-normal mode, and we determine its frequencies as function of the energy. Considering the BMN matrix model as a truncation of $\mathcal{N}=4$ SYM, we also compute the frequencies of the quasi-normal modes of the dual scalar fields in the AdS$_5$-Schwarzschild background. We compare the results, and we find a surprising similarity.


INTRODUCTION
One of the most fascinating phenomena in statistical systems is certainly the emergence of the macroscopic laws of physics. Worth mentioning are some results obtained from the study of random matrices: In the seminal paper [1] Dyson showed that the dynamics of the eigenvalues of a random matrix resembles that of a "Coulomb gas". More recently, Blaizot and Nowak pointed out [2] that the time evolution of the one-particle density of this Coulomb gas is related to the Burgers equation of fluid dynamics. In this work we consider a similar problem for the BMN matrix model [3], and we study the classical time evolution of an analogous "Dyson fluid" distribution, defined by the ensemble of initial conditions of each independent matrix element. We will show that qualitatively the Dyson fluid of the BMN matrix model vibrates as a black hole in AdS. The intuition behind the emergence of this unique type of collective behavior comes from the AdS/CFT correspondence, which we briefly review in the second part of this paper.
The dynamics of the classical BMN matrix model is non-commutative, and there are two features that have to be emphasized. In the first place, the dynamical system is chaotic [5], therefore during the time evolution any localized ensemble of initial conditions will spread over the allowed phase space. Secondly, the expectation values of the observables equilibrate at late times [6], meaning that all the elements in the ensemble populate the phase space according to a time-independent distribution. Such equilibrium distribution is not connected perturbatively to the vacuum configuration of the matrix model, and we generate it numerically by molecular dynamics. Once the system has equilibrated, we study fluctuations close-toequilibrium by activating a gauge invariant quench protocol which induces specific deformations of the equilibrium distribution. The quench protocol is novel, and allows us to perturb the system in a controlled manner. Ought to the expansive nature of the Hamiltonian flow, we then expect the ensemble to quickly approach a new equilibrium. The observables re-equilibrate via quasi-normal oscillations which are characterized by a complex frequency, i.e a damping rate and a ringing frequency. The lowest quasi-normal mode of twist two single trace scalar operators is isolated, and its frequency analyzed as a function of the parameter of the ensemble. Our results extend and consolidate into a unified statistical framework the interesting work of [7]. The BMN matrix model can be obtained as a classical consistent truncation of N = 4 SYM in 4d. Moreover, N = 4 SYM is famously known to admit an holographic description at strong 't Hooft coupling. Holography provides a unique realization of the idea that, for a class of gauge theories, the large-N limit of correlation functions is controlled by a master field configuration [8]. It is unique because the master field configuration is found to be a solution of a classical gravitational problem in Anti-de-Sitter space with prescribed boundary conditions. From independent field theory computations, there is nowadays a compelling evidence that in certain supersymmetric cases holography indeed provides the master field configuration at strong 't Hooft coupling [9][10][11]. However, under the generic assumptions of the AdS/CFT duality, any solution of the gravitation problem is in correspondence with a field configuration in the dual field theory. In the class of finite energy solutions, the most interesting ones are certainly black hole solutions undergoing non trivial real-time dynamics [12].
In the brane engineering of N = 4 SYM, the Dyson fluid represents a simple statistical model for a noncommutative ensemble of matrices at finite energy, in which the off-diagonal degrees of freedom of the branes are fluctuating. The phase space interpretation of the dynamics is complex, but we can look at simpler "geometrical" objects: these are the distribution of eigenvalues associated to gauge invariant operators, which we will show to be non trivial, and to have a finite large-N behavior. This is a very different notion of geometry compared to that offered by the AdS/CFT correspondence. Bearing in mind the impossibility of a direct comparison, we would like to understand, in a dynamical setting, how much these two notions of geometry are far from each other, and in particular we would like to ask how much the process of relaxation in the BMN Dyson fluid differs from that of a black hole in AdS. Since the late time relaxation of a black hole is determined by the spectrum of quasi-normal modes at the equilibrium, we shall answer our curiosity by considering the AdS 5 -Schwarzschild black hole as prototype. We will find that the parametric behavior of the quasi-normal frequencies in the two systems is surprisingly similar. We can interpret these similarities as an indication that a proper path integral formulation of the full N = 4 Dyson fluid might shed new light on the ensemble of black hole microstates. Even though a-priori our comparison has been rather disparate, another possibility, which we have not explored in this paper, would be to study the expectation value of operators of large R-charge, as the BMN operators, and look for the limiting behavior of their quasi-normal frequencies. We should emphasize that the computation we have performed in the matrix model in order to extract the quasi-normal frequencies is highly non trivial, and numerically challenging. The same would be true for the quasi-normal frequencies in the dual black hole background, since the latter is only known numerically [13].
The rest of the paper is organized as follows: In Section I, we fix our notation and we introduce the BMN matrix model. In Section II we describe the statistical framework which will be used to define our Dyson fluid. In particular, we discuss our "big-bang" initial condition, and we clarify some important aspects concerning the interplay between observables, gauge invariance, and noncommutativity. In Section III we study the dynamics of the system from the "big-bang" to the equilibration: We explain how the process of equilibration is related to chaos, and we prove that the joint probability distribution at the equilibrium does not factorize. In Section IV we describe our gauge invariant quench protocol, and we then study the quasi-normal oscillations of the scalar operators in the 20 of SO (6). In Section V, we take a short detour into the AdS/CFT, and we compute in the AdS 5 -Schwarzschild black hole background the quasi-normal modes of the scalar fields representing the 20 of SO (6). In Section VI we compare the behavior of the the frequencies of the quasi-normal modes in the Dyson fluid and in the black hole background. Finally, in Section VII we conclude with a summary and an outlook.

I. THE BMN MATRIX MODEL
We consider a set of N × N hermitian matrices X I (t) and the action: where and The notation is as follows.
No distinction is made between an upper and a lower index. We use capital letters I, J, to label the whole set of matrices {1, . . . , n tot } with n tot = 9, we use instead small letters of the type a and i to label two complementary subsets of {1, . . . , n tot }, respectively: 1 ≤ a ≤ r with r = 3 and r + 1 ≤ i ≤ n tot . The generalization to arbitrary n tot and r is straightforward. The values of n tot = 9 and r = 3 correspond to the bosonic truncation of the BMN supersymmetric matrix model [3]. The Lagrangian is known as the BFSS matrix model [14]. Mass terms in the BMN matrix model will not be important for the dynamics of the Dyson fluid, the main differences between BFSS and BMN come from the cubic interaction. The action (1) is invariant under local U (N ) gauge transformations, with U (t) a generic matrix in U (N ). The BFSS matrix model has maximal SO(9) global symmetry, whereas the global symmetry group of the BMN matrix model is reduced by the mass terms to SO(3) × SO (6).
Hence, the set of interacting degrees of freedom coincides with the set of hermitian and traceless matrices. It is then convenient to parametrize the X I in terms of the generators T m in the fundamental of su(N ), The r.h.s of the equations of motions can be written as The gauge field A t has no kinetic term, and its equation of motion becomes a constraint: The phase space P associated to our dynamical system is described by all the degrees of freedom in the variables X I=1...9 , and P I=1...9 , where A point X ∈ P represents a configuration of matrices and momenta denoted as X = {(X I , P I )} I=1...9 , or equiva- where the index m runs over the basis {T m } N 2 −1 m=1 of generators, and the vector notation stands for x m = (x 1 m , . . . , x 9 m ). In the Hamiltonian formalism, and the equations of motions (8)-(9) become equivalent to the 1st order system supplemented by the constraint (11): I [P I , X I ] = 0. The integral form of the equations of motion defines the Hamiltonian flow ϕ : I × P → P, where I ⊂ R is a interval of time, and ϕ is a map such that for each initial point X (0) in phase space, the path γ(t) := ϕ t (X (0) ) is the unique curve with initial condition γ(0) = X (0) . The flow commutes with the action of the global symmetries of the Hamiltonian.

Conserved Charges
In the BMN matrix model, the conserved quantities are the energy E and the Noether charges of the SO(3) × SO(6) symmetry group. The energy is obtained by evaluating the Hamiltonian H given in (13). The SO(3) and SO(6) charges, dubbed L c=1,2,3 and J q=1,...,15 , respectively, are where the matrices {A c } c=1,2,3 , and {Y q } q=1,...,15 generate rotations in R 3 , and R 6 , respectively. The canonical form of A k is that of an anti-symmetric matrix whose upper triangular part has 0 everywhere, except for +1 in the position (ā,c), corresponding to the plane (Xā, Xc) which is being rotated. Similarly for Y q . For example, the combination Tr X 5 P 9 − X 9 P 5 is the charge associated to rotations in the plane (X 5 , X 9 ). From the equations of motion (15)- (14), it can be explicitly checked that

II. STATISTICS AND DYNAMICS
In the previous section we have introduced the BMN matrix model as a dynamical system focusing on various aspects of the time evolution of a single configuration. In this section we describe a more general framework in which the degrees of freedom are interpreted as the microscopic elements of a statistical ensemble. The statistical framework that we are advocating is an example of a"Dyson fluid", and constitutes our starting point for the study of non-equilibrium dynamics in the BMN matrix model.
The definition of equilibrium that we shall use throughout the paper is the following. Given: a set of initial conditions, and an algebra of observables, a dynamical system is said to be at the equilibrium w.r.t. the given initial conditions, if the expectation value of any observable in the algebra is time independent. Let us point out that the when a statistical equilibrium is reached, the microscopic degrees of freedom are not necessarily steady, and in fact they can have an highly non trivial dynamics. As we are going to show, the process of equilibration in the BMN matrix model has precisely this feature.
Before presenting our numerical results we describe in great detail our choice of initial conditions, and we emphasize some important properties of the observables that we shall study. Several aspects of our analysis will be generic, therefore we expect our results to provide the key elements towards the understanding of the process of equilibration in non-commutative dynamical systems.

Big-Bang Initial Conditions
We consider the following class of out-of-equilibrium initial conditions: where TGU stands for "Traceless Gaussian Unitary ensemble". In practice, the momenta are parametrized as: and each coefficient p I m is extracted randomly from a single gaussian distribution, which we take to have standard deviation σ and zero mean. Because of this property, we expect [P I , P J ] = 0 for all of the randomly generated initial conditions, apart from a set of zero measure. The constraint (11) is satisfied at the initial time, and therefore A t will remain zero during the time evolution.
Intuitively, we designed the initial conditions (18) having in mind a big-bang in which the variables {X I } start at the origin with random momenta. As soon as the {X I } do not commute, interactions are turned on. This can be understood by looking at the force terms in (8)-(9), namely For a given initial condition we shall find [P I , P J ] = 0 at the initial time, and after an infinitesimal time step of δt, X I (δt) = X I (0) + δt P I (0) = δt P I (0) , thus F We generate a finite set of initial conditions of the form (18)- (19), which will be denoted by E(0), where E stands for ensemble. By definition E(t) = ϕ t (E(0)). Notice that E(0) belongs to a Lagrangian submanifold of the phase space, centered at X I = P I = 0, and extended only along the directions of the momenta. The numerical integration is achieved through an improved second order leap-frog algorithm [15]. The time step used in the simulation is δt = 0.05. The stability of the integration has been checked by increasing and decreasing δt of a factor of two. Once the ensemble E(t) is obtained, we compute the expectation value of any observable O through the estimator, where vol(E(t)) equals the number of configurations, and we determine the error on O by means of a jackknife analysis.
On E(t), the SO(3) × SO(6) charges vanish because the {X I = 0}, whereas the energy is non zero. Since the energy is purely kinematical at the initial time, we can calculate its expectation value and its standard deviation, analytically, It is convenient to define the parameter h and rescale σ in such a way that E does not depend on N , Then, E = h n tot . For concreteness we will always specify the dependence on the expectation value of the energy by referring to E = E h . The statistical framework developed so far has a general validity. For the BFSS matrix model, the study of the ensemble greatly simplifies because when m = 0 the equations of motion enjoy the scaling symmetry: t → λ −1 t with X I → λX I [7]. In particular, if two configuration in phase space are related as X 2 = λX 1 , we find that Therefore, given one ensemble of fixed h, any other ensemble of finite energy can be generated by means of the scaling symmetry. The BMN matrix model instead, depends on the value of m in a non trivial way. The strategy in this case is to keep fixed the value of m, and generate the one-parameter family E h by varying h. There is no loss of generality in doing so, because ensembles in which m has a different value can be obtained by the rescaling: t → λ −1 t, X I → λX I with m → λm. Let us mention that σ → λ 2 σ under the scaling, because σ generates momenta. In the numerical simulations, we will fix the value of the mass to m = 3.

Commuting Solutions
Non interacting solutions are described by the commuting ansatz [X I , X J ] = 0, ∀ I, J. In this case the equations of motion admit a simple set of time-periodic solutions. It is useful to identify such solutions, since comparing commuting versus non-commuting solutions will also clarify in which sense the latter are different.
Assuming [X I , X J ] = 0, a common basis of eigenvectors exists such that the X I=1,...9 are diagonal. The diagonal degrees of freedom are decoupled and become simple harmonic oscillators. Assuming X I = 0 at the initial time, t = 0, we find where p I m are the initial velocities, constrained by the requirement that X I is traceless. In the limit m → 0 we recover solutions of the BFSS matrix model with the given initial conditions, In the BFSS model, constant commuting matrices X I parametrize the flat directions of F (4) (and F (3) ). Such flat directions are lifted by the mass terms in the BMN potential, and therefore the solution acquires a non trivial time dependence. Let us mention that the BMN potential admits a set of zero energy vacua known as fuzzy spheres. These vacua are labelled by adjoint su (2) representations, and they are defined by the conditions The Hamiltonian flow of fuzzy spheres configurations has been discussed in [6].

Local Observables
The dynamics of E h (t) will be monitored by measuring the expectation value of gauge invariant observables as function of time. The simplest observables are the kinetic energy, and the singlets of SO(3) and SO (6), More interesting operators belong to non trivial irreducible representations (irrep) of the global symmetry group. In this work, we will consider mainly the symmetric traceless irrep of SO (6), of dimension 20. This irrep can be conveniently decomposed according to the SU (3) subgroup of SO (6), as 8 ⊕ 6 ⊕6. In particular, we will focus on, where M is polynomial which may depend both on the coordinates {X I } and the momenta {D t X I }. Gauge invariance is ensured by the transformation law In this formalism, the operators O s and N i=1,2 have a corresponding M of the form with f I a simple polynomial. When an observable O is defined as in (34), and furthermore the matrix M is hermitian, calculating O is equivalent to summing over the real eigenvalues of M. Two different observables, defined by M 1 and M 2 , respectively, will not in general have a common basis of eigenvectors, because the {X I } do not commute. However, even though the eigenvectors of M i=1,2 have no intrinsic meaning, the eigenvalues can be understood as particular functions of the dynamical degrees of freedom {x I m }, and as such, studied on their own. In this case, an interesting quantity to consider is the distribution of eigenvalues of a given M, integrated over the ensemble.
The relation between the eigenvalues of the {X I }, the observables, and the gauge invariance of the system, is more delicate. Some remarks are in order: 1. Given a configuration X , we may choose to diagonalize one (and only one) of the matrices {X I } by acting with a local gauge transformation U (t). The eigenvalues of this matrix are promoted to dynamical variables, and the price we pay is to introduce a non trivial gauge field, A t = −i(∂ t U )U † . The rotated configuration, X (U ) (t), carries the same information of X (t).
2. As we mentioned, observables whose matrix M is hermitian, are sensitive only to the eigenvalues of M. A stronger statements holds for observables whose matrix M has the form given in (36). In this case, the observables are sensible only to the eigenvalues of the {X I }, even though dynamically it is not possible to diagonalize simultaneously all the matrices. We could say that the eigenvalues of the matrices {X I } represents the first coarse-grained variables built out of the microscopic degrees of freedom {x I m }.

III. ASPECTS OF EQUILIBRATION
From the numerical results it is clear that the big bang is a far-from-equilibrium initial condition, but also that the system equilibrates at late times. The typical behavior of the observables is exemplified in Figure 1: Right after the big-bang these observables experience a highly with τ ≈ 700 in units of time. The dependence of τ as function of h can be inferred by dimensional analysis. The scaling we should keep in mind is In the case m = 0, we deduce immediately that τ h 1/4 = C(N ), where C is a dimensionless constant which only depends on N . When the mass is non vanishing, another dimensionless ratio can be considered. In this case the right parametrization is with f a function such that f → 0 as h 1/4 m. The same scaling argument can be used to infer the expectation value at the equilibrium of any observable O of dimension δ, as function of h and m: Both these formulae contain implicitly the assumption that the system equilibrates. We observed experimentally that the regimes of applicability are two: 1. the asymptotic regime h 1/4 m;

the intermediate regime h 1/4 ∼ m.
In the first case the dynamics tends to the BFSS matrix model because as h → ∞ the mass term in the BMN Lagrangian is negligible. In the second case instead, corrections due to the mass term become important. The case h 1/4 m is more complicated, since the mass term are not suppressed, and the time evolution could be comparable to that of harmonic oscillator (possibly for a time longer than our simulation time). In the next sections we will consider numerical simulations in which either 1) or 2) are valid.
The time evolution of the O . The trivial result N 1 = 0, also shown in Figure 1, was expected. Indeed, the time evolution commutes with the action of the global symmetries, and since the initial data at the big-bang are determined by a single gaussian, E h is symmetric under SO(3) × SO (6). It follows that for any configuration X ∈ E h , and group element g, the 'rotated' configuration g · X also belongs to E h , therefore, the expectation value of any charged observable vanishes. The results obtained so far extend in different directions those obtained in [6,7].

Equilibration and Chaos
Some initial conditions do not lead the system to a late time equilibrium. For example, the expectation value of any observable evaluated on the solutions (26) is obviously periodic in time. Thus, the process of equilibration must be triggered by the non-commutative interaction in the Lagrangian. However, this microscopic feature of the system cannot by itself explain the statistical equilibration of the expectation values of the observables, which is instead a collective phenomena. As we are going to see, chaos is the key towards our understanding of the process of equilibration.
In order to relate chaos and equilibration, we have to discuss an important point. On one hand, the BMN Hamiltonian is non dissipative, and any configuration in the ensemble will keep evolving under time evolution. On the other hand, the expectation value of the observables becomes time independent at late times. It is useful to resolve this apparent logical difficulty by first considering a simpler situation: We may decide to calculate where S is the submanifold in phase space allowed by the conserved charges, and dµ is the flow-preserving volume form on S. Acting with ϕ t on S we may calculate as well Even though the flow will move any single point in S, because ϕ t (S) = S, we conclude that O t = O 0 , and therefore O is at the equilibrium. The same conclusion would still hold true, if instead of S, we consider a region A dense in S, and we assume the flow to be such that for any t ≥ t eq the set ϕ t (A) is uniformly distributed in S. In the most general situation, equilibration takes place under the weaker condition that ϕ t (A) at late times defines a probability distribution which is time independent. This distribution, called P t hereafter, is obtained in the statistical sense by taking the continuum limit in the volume of E at fixed time, i.e.
It is therefore well defined both for small and large-N. The case of an uniform P t at the equilibrium is the case we mentioned in the example given above.
In our finite N simulations, we should also bear in mind another detail: The fluid, E h (t), is not restricted to a single submanifold S but covers the set ∪ γ S H where H is the energy of a single configuration in the ensemble. In fact, the energy of the ensemble is a gaussian variable with mean E = h n tot and standard deviation √ 2 E / n tot (N 2 − 1), as shown in (23). The property we have invoked about the flow is very much related to the definition of Lyapunov chaos [16]. An Hamiltonian flow ϕ is said to be chaotic in the sense of Lyapunov if two properties hold: 1. ϕ t is almost everywhere expansive, 2. under time evolution ϕ t (p) visits almost every point in S, i.e. ϕ t is topological transitivity 1 .
In particular, the well known Lyapunov exponent λ Lya quantifies how much ϕ t is expansive. We have measured λ Lya in our simulations, by using the strategy outlined in [5], for the BFSS matrix model. Within error we obtain the same result: λ Lya ≈ 0.29+O(1/N ). Even though the ensemble E(0) was gaussian and localized in phase space, as a consequence of the expansive nature of the flow, ϕ t (E(0)) populates the allowed phase at late times. This is an important feature of the dynamics of chaotic systems as opposed to that of integrable systems. In fact, by taking a constant distribution of initial conditions, it is not difficult to build a fine-tuned equilibration even for the harmonic oscillator. However, this type of equilibration is non generic and will fail upon specifying a different set of initial conditions. In this sense, equilibration for chaotic systems is generic for any given initial conditions.

Dynamics of Eigenvalues
Building on our previous remarks about the dynamics of observables of the type we can re-express the expectation value of these observables as the integral ∀ i, j, a, b, so that only two non trivial eigenvalue distributions exists, one for each symmetry group. Accordingly, the charged operators N i=1,2 have vanishing expectation values. In the BFSS model the symmetry group is enhanced to SO(9) and we would find in addition that ρ eig i = ρ eig a ∀ i, a. In the BMN model m = 0, and we do not expect this relation to hold. In fact, in the setup of Figure 1, we can explicitly verify that O The histogram of ρ eig (λ) in the SO(6) sector for N = 6 and E h=0.3 is shown in Figure 2. The solid curve is the TGU distribution, which can be computed analitically from the results of [17]. The numerical and the TGU distributions agree. For completeness, let us illustrate the analytic calculation: We quote from [17] the profile function The good agreement between ρ eig T GU (λ) and the distribution of eigenvalues of any X I at the equilibrium should not be confused with the fact that the ensemble at the equilibrium is TGU. In fact, as a consequence of the non trivial dynamics, the matrices are correlated, and the joint probability distribution does not factorize, i.e. P( x m , y n ) = P( x m )P( y n ). A simple convincing argument is the following: Let us assume that Z 1 and Z 2 are uncorrelated and TGU with variance σ. Then and we obtain the exact relation, .
The presence of dynamical correlations in E h (t) at the equilibrium can be detected by studying, for example, the time evolution of Y[X 3 , X 4 ]. In the top panel of Figure 3 we have plotted the expectation value of this observable for E h=0.3 and N = 10. At the equilibrium, Y[X 3 , X 4 ] deviates considerably from unity, and we conclude that E(t) is not TGU. Since the approach to the equilibrium is exponentially fast, the result is robust. From dimensional analysis we know that the value of Y[X 3 , X 4 ] at the equilibrium (hereafter Y eq ) is a function of N and m/h 1/4 . In particular, as h is increased, we expect to recover the result at m = 0 (which is h independent). From our simulations we are able to confirm the correctness of this argument. In fact, starting from the result of Figure 3, where Y eq ≈ 1.19 for h = 0.3, we have tested the converge of Y eq measuring Y eq ≈ 1.26 for h = 38.4, and Y eq ≈ 1.31 at m = 0. Differently from the X I , we can shown that the momenta belong to a TGU ensemble. This can be inferred from the bottom panel of Figure 3, where we have measured Y[P 3 , P 4 ]. At late times, Y[P 3 , P 4 ]−1 is zero within the error. Notice that at very early times instead, the time evolution of the coordinates is driven by that of the the momenta, and since the latter are initiated as TGU, the value of Y[X 3 , X 4 ] is still close to unity as long as the non-commutative interactions can be neglected.

IV. DYNAMICAL RELAXATION
The aim of this section is to study the close-to-the equilibrium relaxation of the observables. In order to do so, we devise a quench protocol that takes E h (t) at the equilibrium, and push it onto a phase space region whose distance from E h (t) can be tuned in a controlled way. We define the notion of distance between two ensembles by measuring the differences between the conserved charges. By assuming that the level sets of the conserved charges change in a smooth way, a given quench will generate a close-to-equilibrium configuration when the variation of the conserved charges before and after the quench is small. Instead, the quench will generate a far-from-equilibrium configuration when the conserved charges are changed abruptly. As we are going to show, our results will indeed confirm this phase space picture.

Quench Protocol
Let us begin by describing our three-steps protocol: 1. From the initial big-bang we wait until the system equilibrates, say at time t quench .
2. We perform a global gauge transformation with an unitary matrix U on (the history of) any configuration X , i.e.
and we choose U such that at t quench , and for a given index p, the matrix U X p U † = D is diagonal; 3. For any D, we order the eigenvalues of D = diag(λ 1 , . . . , λ N ) from least to greatest, and we perform the shift with quench parameters α=1,...N .
The constraint is not violated if at the same time we perform a deformation of the momentum Let us notice that in order for this deformation to actually take place as N → ∞, it may be necessary to consider, D ≡ D + diag(− 1 , .., − k , 0, .., 0, + k , .., + 1 ) , (58) while keeping k/N fixed. It is worth emphasizing at this point, that the quench is not symmetric under → − . This can be understood by considering the denominators of δ[P, D]. For example, the quantity −(D 11 − ) + D 22 depends on both the sign of , and the difference between two near-neighbour eigenvalues of D. For one sign , there exists a small value of − for which the denominator diverges, and as a result, we expect the quenched configuration to be a far-from-equilibrium configuration.
• Consider a quench in which D is inflated. This configuration is achieved by defining α = λ α , with a given . In matrix notation Notice that the traceless condition is automatically satisfied. The solution of δ[D, P ] also takes a very simple form, in particular the quench along the momentum reduces to P p . Again the traceless condition is automatically satisfied.
Finally, we remark that A t = 0 remains solution after the quench.

Features of the Quenched Ensemble
The outcome of the quench protocol is an initial condition for a new ensemble E , which depends explicitly on the quench parameters and implicitly on h. This new ensemble will evolve on a sub-manifold of the phase space which is not that of E h . As we mentioned, the notion of "distance" between ensembles is better quantified by looking at the variation of the conserved charges. For example, if we consider an edge-type quench we should calculate the expectation value of the following quantities, Since we are especially interested in generating close-toequilibrium initial conditions, we will tune the quench parameters in such a way that the conserved charges admit a perturbative expansion. Even in this regime, it is highly non trivial to compute their expectation value analytically. Nevertheless, simple arguments based on dimensional analysis and the scaling law (38) provide us with the right parametrization in terms of and h. Continuing with edge-type quenches, the variation of the energy ∆ E = 9h − H(E ) is then where the c i only depend on N . This is the most general expression compatible with both the scalings h → λ 4 h, → λ , and the limit → 0. For the BFSS matrix model, the conserved quantities will only depend on /h 1/4 , but in the generic case there are other contributions in h as long as h/m ∼ 1. Arguments based on dimensional analysis are also valid for inflation-type quenches, the only difference to bear in mind is that now plays the role of /h 1/4 . Some formulae simplify. For example, the order in the E is where m depends on whether p ∈ {a, i}, according to (3) or (4). G p represents the difference between kinetic and potential energy of the configurations (P p , X p ). As expected from considerations about the equipartition theorem, we checked that G p vanishes. Therefore, ∆ E is at least quadratic in for inflation-type of quenches. A more concrete way of understanding how the quench acts on the ensemble is to consider its consequences on the dynamics of the observables we are interested in. According to our previous discussion, from the knowledge of the distribution of eigenvalues ρ eig Deforming one of these distributions is precisely the task of the quench protocol. In fact, by rotating X p , we obtain N eigenvalues extracted from ρ eig p at the equilibrium, and by shifting D to D, in practice, we deform ρ eig p by changing its support. Therefore, a mismatch of the same order of magnitude of the quench parameters exists between the the shape of the distribution at t + quench and that of the equilibrium distribution at t − quench . This picture is particularly helpful if we want to visualize in which way E represents a close-to-equilibrium initial condition from the point of view of the observable.
In Figure 4 we compare the distribution just before and after the quench time for edge-type quenches. For any value of N , the equilibrium distribution in E is characterized by N distinct peaks. For small N the picture is simpler. At t + quench , we find that the position of the peaks at the edge of the distribution has moved from the inside out of order , whereas the bulk of the distribution has remained almost unchanged. Since the number of peaks increases with N , the same logic goes through. We may get a feeling about the large-N result by repeating the quench protocol in the case of a TGU matrix. As Figure 4 shows, the outcome of the quench produces a ripple in the equilibrium distribution.
For each event in E the SO(3) × SO(6) symmetry is broken to the subgroup of rotations that leave X p fixed. The breaking is explicit, and the full symmetry is not restored on the ensemble since the direction of the broken charges is the same for all configurations. Thus, if p ∈ {1, 2, 3} (or p ∈ {4, . . . , 9}) E preserves SO(2) × SO(6) (or SO(3) × SO(5)). After reaching the new equilibrium in E , we expect the set of ρ eig I with I = p to be equal, according to the preserved symmetries, whereas the distribution ρ eig p to be different in a way which is determined by the strength of the quench parameter. Given the asymmetry in the eigenvalues distributions, the charged operators N i=1,2 will acquire a non trivial expectation value. In the next section we will study the time evolution of these operators for closeto-equilibrium quenches. Some examples of edge-type quenches in which 1 are illustrated in Appendix B.

Quasi-Normal Modes
The time evolution of the observables after the quench has some notable features. In the top panel of Figure 5 we have plotted the expectation value of N 2 , from the quench to the new equilibration 2 . Two facts show up clearly: 1. N 2 relaxes fast, without experiencing any non linear transition typical of the big-bang.
2. The process of relaxation is driven by quasi-normal oscillations. These modes are collective excitations of the non-commutative dynamics.
Before discussing which ansatz describes the quasinormal ringing, it is interesting to ask how we should think about these oscillations in phase space. Let us consider that at the equilibrium E(t) is described by a distribution ρ t which does not change in time, so intuitively, any realization of E(t) covers the support of this distribution with the correct weight. As we pointed out, the equilibrium distribution after the quench has moved to a different support, and evidently E (t + quench ) does not cover enough of it. Under the assumption that E (t + quench ) is a close-to-equilibrium initial conditions, we expect the holes between E (t + quench ) and the support of the new equilibrium distribution, to populate quickly. During this process the dynamics of the observables is driven by quasi-normal oscillations. The ring-down of the observables is in one-to-one correspondence with the ring-down of the distribution of eigenvalues, since both are functions of the microscopic degrees of freedom.
The behavior of the O (t) after the quench is a superposition of quasi-normal modes. Close to the new equilibrium, the lowest quasi-normal mode is the dominant one, and the behavior of the observable is fitted into the following ansatz, where O stands for any of the observables we are considering. The details of the fitting procedure are reported in the Appendix A. Here, τ −1 and Ω are the damping and the ringing frequencies, respectively. In complex notation, it is convenient to define ω ≡ Ω − iτ −1 . Tuning the values of d and u, the ansatz accommodates both an over-and an under-damping behavior. It is important to emphasize that the quenched ensemble, by construction, has different properties compared to E h . In particular, we expect d, u and O ∞ to be directly proportionals to 2 We have used a volume of E of order 10 3 , and from each of these configurations at the equilibrium we have evolved approximately 10 3 quenches. the quench parameters because of the explicit symmetry breaking. Since we are interested in studying the fluctuations of the equilibrium distribution in E h , we will consider the limit → 0, for edge-type quenches, and → 0 for inflation-type quenches. The only quantities whose limiting values will be non trivial are the frequencies τ −1 and Ω. Once again, dimensional analysis provides us with the basic tools to understand the data. The complex frequency can be parametrized as, where c 0 is coefficients, and f is a function which vanish in the asymptotic regime. The limit h 1/4 m is very useful here, because as long as the mass parameter is irrelevant the BFSS and the BMN model have a similar dynamics. On one hand, we are free to study the BFSS model independently, setting m = 0, with the advantage that the model has a simpler dynamics. On the other hand, we can check that the values of ω, taken from the BMN model, asymptote those of BFSS model. Comparing the results of BFSS simulations at N = 10, 15, 20, within error we find the relation We conclude that in the asymptotic regime, the natural variable upon which ω and τ depend is h/N . We should mention that at smaller values of N , for example N = 5, 6, we have measured small deviations from the scaling regime (67). On the other hand, the dependence on h can be understood as follows: In the large-N limit we would find the relation h/N = ρ(λ H )λ H , where ρ is the distribution of eigenvalues of the BMN Hamiltonian H. Therefore, in order to have a well defined energy distribution we should keep h/N fixed. The scaling with N actually holds for the entire profile of the quasi-normal oscillation. This is shown in the bottom panel of Figure 5, where we have plotted the behavior of N 2 (t) for inflation-type quenches by keeping h and N fixed. The latter condition follows from the observation that Within error we cannot appreciate any difference on N 2 (t) among N = 10, 15, 20, and the deviations visible at N = 5 can be addressed to 1/N corrections. Outside the asymptotic regime, corrections due to the mass m cannot be neglected, and we expect them to organize in a series expansion of the form where p and c 1 need to be determined. We have found that p = 1/2 and c 1 is, within error, an N independent constant. In Figure 6 we have plotted ωh −1/4 as function of m 2 /h 1/2 . We can fit these two curves with the ansatz, with coefficients: It is worth mentioning that the numerical fit of (70) do not depend on the details of the quench protocol, in particular we do not find differences between edge-and inflation-type quenches when we look at the curves of ω h −1/4 . In the linear response regime, this is a consequence of the fact that the location of the poles of the Green functions are independent of the strength of the perturbation.
The BMN matrix model can be understood as a classical (supersymmetric) consistent truncation of N = 4 SYM on R × S 3 . This connection was established in [4], where the authors showed that the equations of motions of the N = 4 gauge multiplet (A µ , χ A=1,...,4 α , φ i=4,...,9 ), truncated to the lowest harmonics of S 3 , reduce to those of the BMN matrix model. The details of the truncation are as follows, where X i ⊕ V c+ ⊕ S β+ span the irrep (1, 1, 6) ⊕ (3, 1, 1) ⊕ (2, 1, 4) of the SU (2) + ⊗ SU (2) − ⊗ SU (4) R symmetry group. In (1) we considered only the scalar sector of this truncation, setting the fermions θ A β to zero, and redefining the radius R of the S 3 as m = 6/R. Recalling that the covariant derivatives acting on the gauge multiplet are D µ = ∇ µ − ig Y M A∧ and F = dA − ig Y M A ∧ A, and the quartic coupling is g 2 Y M [X I , X J ] 2 , we obtain the BMN action by considering the field redefinition X = X/g Y M . The final result is In the regime where g 2

Y M
1 and Xg Y M is kept fixed, the classical saddle point approximation is justified, and the study of the BMN matrix model carried out throughout sections I-IV can be reinterpreted as the study of N = 4 SYM in the BMN truncation in the classical limit. The only change we need to implement is to redefine the energy as: where H is the BMN Hamiltonian (13).
In the gauge theory, the Dyson fluid represents a finite energy ensemble of classical D3-branes in which the offdiagonal modes are fluctuating. Microscopically, the system is non-commutative and non-perturbative with respect to the dynamics of the diagonal degrees of freedom. Therefore, the weak-coupling "geometric" interpretation of the D3-branes has to be rediscovered through the looking glass of the gauge invariant operators of N = 4 SYM. In particular, a "geometric" picture should emerge from the description of the system offered by the eigenvalue distributions associated to each single trace operators. This picture is inspired by the AdS/CFT correspondence which, in the regime of strong 't Hooft coupling, relates the dynamics of field configurations to those of a gravitational problem in the Anti-de-Sitter space. Being aware that a direct comparison between AdS physics and our Dyson fluid would not be possible, we find important to study the two dynamics in a non trivial case, precisely with the aim of highlighting the major differences.
In the following we review basic aspects of the AdS/CFT correspondence, and we will compute quasinormal modes for the dual operators that we studied in the BMN truncation. Since at this point we are only interested in qualitative features, we will carry out our toy model computation in a simple black-hole geometry, and we will not go beyond the AdS 5 -Schwarzschild black hole.

AdS/CFT Correspondence and Black Holes
One of the greatest achievements of string theory has been the discovery of the AdS/CFT correspondence, or more generically, of the gauge/gravity dualities. For N = 4 SYM the dual space-time is an asymptotically AdS 5 × S 5 background. The two sides of the duality are related by the the D-brane construction of the field theory [19] which sets, where L AdS5 and L S 5 are the radii of AdS 5 × S 5 , g s is the dimensionless string coupling, and α is the string tension. The radii of AdS 5 and S 5 are equal and we shall refer to them simply as L. The reader unfamiliar with the AdS space might find useful to think about it as a "box" of constant negative curvature proportional to 1/L 2 . The quantity L 2 /α depends only on the 't Hooft coupling λ t ≡ g 2 Y M N . Keeping λ t fixed, and taking the large-N limit, the string theory becomes perturbative, i.e. g s 1. Then, two cases are well under control: 1. When λ t 1, the AdS 5 × S 5 background is classical and the string theory can be truncated to 10d classical gravity with Newton constant G 10 = 8π 6 g 2 s α 4 [20]. In this regime, the duality is very powerful and predicts that N = 4 SYM at strong coupling is dual to a gravitational theory in AdS 5 × S 5 . The holographic dictionary provides the concrete link between the two theories [21]. In particular, the 4d space-time, where the field theory lives, is identified with the boundary of AdS 5 , and chiral single-trace operators in the field theory are mapped to bulk fields of the 10d geometry. Intuitively, any field configuration at the boundary now acquires a bulk profile along the the fifth (extra) radial coordinate of AdS 5 . This bulk profile is obtained by solving classical gravitational equations of motion. The near-boundary behavior of bulk fields plays an important role in the holographic dictionary, and a more precise statement about it will be made in the next section.
2. When λ t 1 the field theory is perturbative. In this regime the curvature of the gravitational background is large, and the classical spacetime structure has to be modified by quantum gravity corrections.
It is important to realize that according to the basic principles of holography, asymptotically AdS solutions are specified uniquely by the set of conserved charges of the boundary field theory, together with the expectation values of charged operators. These are the same charges we used to initialize and characterize E, as well as more complicated ensembles like E . Therefore, in order to set up a comparison with E, we should at least look for a gravity solution with finite energy, zero momentum 3 , and zero SO(6) charges. Moreover, we should also look for a solution which supports quasinormal oscillations. Together, these two observations point to the simplest and most primitive candidate: the AdS 5 -Schwarzschild black hole.
Black holes have special properties: Since the seminal works of Bekenstein, Hawking and Gibbons, it has been recognized that thermodynamic variables can be assigned to stationary black holes [22,23], and via the AdS/CFT duality it has been understood that an asymptotically AdS black hole corresponds to a field configuration which is approximately thermal. Black holes do not only appear as static objects, but they are also characterized by important dynamical features, in particular by the spectrum of their quasinormal oscillations. In this respect, let us notice that AdS 5 would not have worked for our comparison, simply because at the linearized level it only allows for normal oscillations. By the AdS/CFT duality, the ringing frequency and the decay rate of the lowest quasinormal modes of the black hole control the scales of the process of "thermalization" in the field theory [25].
The AdS5-Schwarzschild black hole In our notation the AdS 5 -Schwarzschild black hole is described by the metric This background is a solution of the 5d Einstein action, where I GH is the standard Gibbons-Hawking surface term, and The 5d Newton constant G 5 is obtained as the ratio of G 10 over the volume of the S 5 , which is trivially fibered in the full 10d geometry of string theory. G 10 has a stringy expression in terms of α [20], thus upon substituting L for α, we obtain G 5 as written in (81). The radial coordinate r extends from the boundary at infinity, corresponding to R × S 3 , to the radius of the horizon, hereafter r h , which is defined as the greatest root of the equation The parameter M is called the non-extremality parameter. When M = 0 the geometry is that of pure AdS 5 in global coordinates. The radius of the boundary S 3 is measured in units of L, i.e. R = R/L, where R is the field theory quantity.

Boundary Energy and Thermodynamics
As we mentioned earlier, the mapping between gravitational solutions and field theory configurations goes through the correct identification of the conserved charges on both sides. In our case, the energy. From the bulk perspective, the quantity we need to evaluate is the time component of the stress-energy tensor integrated over the three-sphere at the boundary. This is properly defined as where u µ = (u t , 0, 0, 0, 0) is the time-like unit vector orthogonal to S 3 , and ξ ν is the Killing vector √ g tt u µ . The stress tensor is where h ab is defined from writing the metric as ds 2 5 = N 2 dr 2 + h ab (dσ a + V a dr)(dσ b + V b dr) (see [26] for our conventions). The expression of T reg µν includes counterterms which regulate well understood divergences in AdS. The result for T reg µν is Plugging this expression in (83) and using the relation between G 5 and L, we finally obtain where vol(S 3 ) = 2π 2 R 3 . The answer for E is the sum of a zero point (Casimir) energy proportional to 1/R 4 and a physical M -dependent energy which can then be rewritten in the form, The factors of g 2 Y M and vol(S 3 ) in this formula are the same as in the field theory definition of the energy, and the actual value of the energy is proportional to N times M . In the field theory E is computed thought the path integral with the insertion of the Hamiltonian, which is a single trace operator. Therefore, E will be proportional to N times the integrated distribution of the eigenvalues of the effective Hamiltonian at strong 't Hooft coupling, i.e. M . We recover in this way, the original meaning of the black hole parameter M , which was first derived in the seminal paper [24].
The energy E can also be interpreted as thermal by introducing the Hawking temperature, In this formulation, the partition function is given by the the on-shell value of the (regularized) gravitational action, the thermal energy E th , and the Bekenstein-Hawking entropy S BH can be obtained as follows, It is simple to check from (82) that E th = E. For a given temperature T there exist two solutions of the equation T H (r h ) = T , and the corresponding black holes are dubbed as "large" if R r h /L > 1/ √ 2, and "small" if R r h /L < 1/ √ 2. Large black holes have positive specific heat and are dual to thermal states in the field theory. Small black holes are always thermodynamically irrelevant since an Hawking-Page transition between large black holes and thermal AdS takes place at R r h = 1. In the field theory, this Hawking-Page transition has been interpreted as a second order transition [27].

Massive Quasinormal Modes of the AdS5-Schwarzschild black hole
In this section we study bulk quasinormal modes of the AdS 5 -Schwarzschild black hole which are dual to the scalar operators in the representation 20 of SO (6). Let us recall that under SU (3) ⊂ SO(6) the representation 20 decomposes into 8 ⊕ 6 ⊕6 and contains the operators N i=1,2 and C i=4,...9 that were analyzed earlier in the context of the BMN matrix model.
On the gravity side, the SO(6) R-symmetry is realized geometrically as isometries of the S 5 . Each scalar operator sitting in a representation of SO(6) is mapped to a specific harmonic deformation of the S 5 . The harmonics of the S 5 sitting in the representation 20 are charged under the isometries of the 5-sphere, and therefore must be coupled to bulk gauge fields. These bulk fields sit in the irrep 15 of SO(6), and their time-components are in one-to-one correspondence with the boundary charges J q=1,... 15 . It is perhaps useful to sketch how these fields are realized in concrete as deformations of the S 5 . Following the notation of [28], the 20 is parametrized by the matrix T ij = T ji in SL(6, R), and the metric of the S 5 is written as The coordinates µ i=1,...6 are subject to the constraint µ i µ i = 1, and the 1-forms Dµ i are the covariant derivatives Dµ i = dµ i + A ij µ j , where the matrix of bulk gauge fields A ij = −A ji represents the 15. The covariant derivative on T is DT = dT + [A, T ]. The round S 5 is recovered from the trivial configuration T = diag(1, 1, 1, 1, 1, 1) with all gauge fields turned off. Two are the 5d solutions in which T can be taken to be trivial: empty AdS 5 , and the AdS 5 -Schwarzschild black hole. The authors of [28] found out the fully non-linear consistent truncation of type IIB supergravity, which only retains gravity and the fields in the 20 ⊕ 15. Any charged black hole in this sector, bold or hairy, can in principle be constructed. For hairy black holes, the field theory configuration is further characterized by the expectation value of the corresponding operators in the 20.
Turning on expectation values for operators of the type C i1,...in T r(X i1 . . . X in ), with n > 2 and C i1,...in totally symmetric and traceless, cannot be done within consistent truncation ansätze, but requires non separable 10d gravitational backgrounds. Coulomb branch solutions are an example, [29].
In general, quasinormal modes materialize as linearized perturbations of black hole backgrounds. For the AdS 5 -Schwarzschild black hole, we are interested in perturbing the T ij sector. Setting T ij = δ ij + Φ ij , we truncate the equations of motion at order . The Φ ij decouple and each component satisfies the equation, where m 2 = −4/L 2 . According to the holographic dictionary we expect ∆(∆ − 4) = m 2 L 2 , where ∆ is the dimension of the dual operator 4 . For the case at hand ∆ = 2 We shall look for solutions of (91) with the following profile, Different harmonics of the spatial S 3 could also be considered, and the calculations would go through with minor modifications. We have taken the lowest harmonic Y (0,0) 0 so to reproduce (72) at the boundary. Changing variables to z = L 2 /r 2 simplifies the equation of motion to where b(z) = 1 + z/R 2 − M z 2 . In this new coordinate z, the boundary is placed at z = 0 and the horizon is The near-boundary behavior of φ is with A and B constants. By rescaling z → R 2z , it is simple to show that w ≡ ωRL, and x = M R 4 are the only parameters entering the equation of φ. By the holographic dictionary, A will be interpreted as being proportional to O , whereas the log z term will be interpreted as a source in the field theory [30]. Here O is any operator in the 20, since as we mentioned, at the linearized level the Φ ij decouple. We are interested in studying how the perturbation of the black hole relaxes when no boundary sources are turned on 5 . Then, we shall find solutions of (93) such that B(w, x) = 0. In order to properly define the Cauchy problem for Φ, we also have to specify a "regularity" condition at the horizon. Near the horizon, a generic solution would behave like where the +(−) sign corresponds to out-going (in-going) modes. Classically, matter can only fall into the horizon, therefore the only allowed solution has to be in-going at the horizon. With this second boundary condition the Cauchy problem is well defined. As shown in the pioneering work [31], it is a general fact that the condition B(w, x) = 0, for fixed x, admits a discrete set of complex solutions of the form w = w Im − iw Re with w Re > 0. These are the so called quasi-normal frequencies, and the lowest one characterizes how the perturbation equilibrates at late times. As it should be, the decay rate w Re comes out positive. We have solved the equation of motion for φ numerically, interpolating between a series solution at the boundary and a series solution at the horizon. Scanning through the complex frequency plane w, we have found the quasinormal frequencies as function of x, i.e we have found the curve w = w(x). It is interesting to recover from this result the relation between ω and the scales of the problem, M and R, according to dimensional analysis arguments. In fact, it is simple to show that in the limit of large R there is only one possibility, w(x) = c 0 x 1/4 , and therefore ωL = c 0 M 1/4 . The finite R dependence instead is non trivial, and follows from the actual numerical solutions. The result fit into the ansatz, where ωL is the proper frequency at the boundary. The numerical results for the lowest and the next-to-lowest quasi-normal modes are shown in Figure 7. For the lowest and O (6) s . However, these two operators are non-chiral and do not appear in the spectrum of supergravity.
and for next-to-lowest quasi-normal mode,

VI. DISCUSSION
Connecting the dots of what we have done so far, we have two parallel situations in N = 4 SYM, a statistical ensemble on one side, a black hole on the other side, and two dynamical quantities which we can now compare: their lowest quasinormal frequency. Let us here emphasize that there is a-priori no reason why we should expect some kind of resemblance, but according to our findings, their qualitative behavior is surprisingly similar. The parameters, h and m in the BMN matrix model, and M and R in the black hole, determine the frequencies very much in the same way. In particular, it is an independent outcome of the two calculations that the leading correction to the flat or zero mass limit, written as an expansion in powers of m 4 /h and 1/M R 4 , is fixed by the power p = 1/2. This result cannot be obtained by dimensional analysis. In the absence of a better analytical understanding, we can only appreciate the unexpected beauty of the precision fit in Figure 6 while considering its striking similarity with that of Figure 7.
Following the work of [34,35], it would be interesting to compute semiclassical corrections to our Dyson fluid distribution, and see how the quasi-normal frequencies are modified at small but finite 't Hooft coupling.

VII. CONCLUSIONS AND OUTLOOK
In this paper we have analyzed Dyson fluid solutions of the BMN matrix model living on Lagrangian submanifolds of the phase space. We have characterized the corresponding ensemble by symmetries, conserved charges, and at the equilibrium, by the expectation values of local operators. We have defined a novel gaugeinvariant quench protocol which allowed us to deform the equilibrium distribution in a controlled way, and we have shown that the expectation values of twist two scalar operators in the SO(6) sector re-equilibrate via quasinormal oscillations. Finally, we have determined the numerical dependence of the complex frequency ω of the lowest quasi-normal mode as function of the energy h and the dimensionless parameter m 4 /h.
The interesting features of the Dyson-fluid are triggered by the non-commutative nature of the microscopic dynamics. The study of eigenvalue distributions of observables offers an alternative geometric picture of the dynamics. In particular, the complexity of an observable provides the tool to detect different type of correlations in the ensemble. In the second part of the paper we have compared the emergent geometric structure of our classical Dyson-fluid and the holographic geometry of a black hole at strong 't Hooft coupling. Our main result has been to point out an unexpected similarity between the parametric behavior of the quasi-normal frequencies in the Dyson-fluid and in the black hole background.
There are similarities between the phase space picture of the quasi-normal oscillations, and the entropic principle of [36]: A rigorous way to understand the process of equilibration would be to consider a covering of the allowed phase space and verify that the population in each patch does not change in time. In this language, we are then led to conjecture the existence of an effective actions whose equations of motions determine the equilibrium distributions of the operators. As it happens for the case of the Brownian motion [2], it would be interesting to rewrite the time evolution of the distributions in terms of fluid-like equations. be actually needed to improve the stability of the results. For simplicity, in this work we have focused mainly on the lowest one.
The errors on the frequencies are estimated from the spread of distribution of the fit parameters obtained from the different jacknife samples used.