Black holes as random particles: entanglement dynamics in infinite range and matrix models

We first propose and study a quantum toy model of black hole dynamics. The model is unitary, displays quantum thermalization, and the Hamiltonian couples every oscillator with every other, a feature intended to emulate the color sector physics of large-$\mathcal{N}$ matrix models. Considering out of equilibrium initial states, we analytically compute the time evolution of every correlator of the theory and of the entanglement entropies, allowing a proper discussion of global thermalization/scrambling of information through the entire system. Microscopic non-locality causes factorization of reduced density matrices, and entanglement just depends on the time evolution of occupation densities. In the second part of the article, we show how the gained intuition extends to large-$\mathcal{N}$ matrix models, where we provide a gauge invariant entanglement entropy for `generalized free fields', again depending solely on the quasinormal frequencies. The results challenge the fast scrambling conjecture and point to a natural scenario for the emergence of the so-called brick wall or stretched horizon. Finally, peculiarities of these models in regards to the thermodynamic limit and the information paradox are highlighted.


Introduction and description of results
The purpose of this article is twofold. We first want to propose and study a quantum model that captures several important features of black hole dynamics and it is still tractable from a unitary microscopic perspective. The desire is to collect enough and robust intuition through a concrete solvable model. The second objective is to extend this intuition to exact descriptions of black holes, such as large-N Matrix Models 1 [1,2,3]. The essential features we want our toy model to possess are the following: • The model must display quantum thermalization [4].
• The model must have the strongest form of non-locality: the Hamiltonian must couple every oscillator with every other oscillator through couplings of the same size [5]. Throughout the article we will call these type of models democratic, to distinguish between other possible forms of non-locality 2 .
The first bullet point is an obvious necessary feature. Since the groundbreaking contributions of Bekenstein [8] and Hawking [9], it is widely accepted that black holes are thermal objects. With the advent of the AdS/CFT correspondence [3], the previous statement has a firm theoretical ground, since AdS/CFT maps black hole physics into a many-body quantum system at finite temperature [10]. The second bullet point is less obvious. Within conjectured non-perturbative definitions of string theory [2,3], black holes are mapped to high energy states of large-N matrix models. These matrix models have interaction terms of the following form: comprising interactions between fundamental and adjoint (matrix valued) degrees of freedom, see [11] for an explicit example applied to quantum black holes. In the high energy sector of the model, the matrix-valued field A ij has all its entries thermally excited. From the point of view of the vector field π i , the previous term connects every oscillator π i with every other π j democratically. Another argument supporting the use of democratic models comes from the recently proposed microscopic model of certain black hole dynamics [12], in which a set of fermions interact through random quartic interactions.
Given the previous set of arguments and requirements, the toy model of black hole dynamics we propose in this article is defined by the following Hamiltonian: where α and η are parameters with dimensions of energy, c † i and c i are creation and annihilation operators of spinless free fermions (with the usual anticommutator relations), and V ij are independent random gaussian real numbers with zero mean and unit variance 3 . Therefore, the matrix (ηV ) ij ≡ ηV ij is a random matrix taken from the Gaussian Orthogonal Ensemble (GOE) with deviation σ ηV = η. For a beautiful and modern treatment of random matrices see [13]. Notice that considering bosons instead of fermions is straightforward. In second-quantized formalism all the computations are similar, and one obtains the same conclusions.
The eigenstate correlations and entanglement patterns of the family of Hamiltonians (1.2) were studied in [4], where it was analytically proved that the model satisfies the so-called Eigenstate Thermalization Hypothesis (ETH) [14] 4 . This hypothesis concerns the origin of statistical mechanics in closed quantum systems. As remarked in [4], the family of Hamiltonians (1.2) shows that ETH is 'typical' within the space of gaussian systems. Non-linear interactions are not needed in this regard, contrary to common belief.
Below we will study two different types of non-equilibrium unitary processes which are particularly interesting from the point of view of black hole physics. We first consider the typical black hole thought experiment, in which we 'throw in' some small system to a black hole. This amounts to set an initial state in which a small subsystem is factored out from the thermal ensemble (2.1). We then compute the time evolution of two-point correlation functions and entanglement entropies (2.21), (2.22), (2.23), and (2.28). The fact that our Hamiltonian is gaussian implies that the computation of all the two-point correlation functions provides every n-point correlation function as well. This will allow us to discuss properly any subtle question regarding information spreading. The second type of initial states we consider are completely factorized states (2.29). This non-equilibrium process intends to emulate 'collapse' scenarios in black hole physics, in which a set of particles with definite quantum numbers collapse beyond their Schwarzschild's radius. This situation turns out to be severely more complicated. Although we are not able to solve for the evolution of occupation densities exactly, we will show interesting relations that the evolution of entanglement entropies satisfy, see (2.35), (2.36), (2.39) and (2.40). These special relations were argued to hold generically for democratic systems in [17].
In section (3) we apply the results to black hole physics. We first analyze the toy model. Without taking into account Poincaré recurrences, the model displays three types of thermalization time scales. The first one is the local relaxation time scale, the analog of quasinormal decay. It is given by t local = 1/R, where R 2 ≡ 4Nη 2 ∼ O(α) can be taken to be independent of the system size. The second one is the scrambling time scale [18,5], defined equivalently as the time by which entanglement entropies become thermal or the time by which information spreads through the whole system, as it was defined in the original reference [5]. It can also be defined as the time at which the mean field approximation ceases to apply. It turns out to be given by t scrambling = t local = 1/R, a result that challenges the fast scrambling conjecture, as we describe properly below. Finally, we notice the existence of another time scale. It is the time by which deviations from thermality are uniform over the entire system. This 'randomization' time scale turns out to be t random = N 1/3 /R for our toy model.
In the second part of section (3), we extend the results and intuition to exact descriptions of black hole dynamics, i.e to large-N matrix models or the SYK model [12]. The extension is possible because the results of the toy model are rooted in a single feature: large-N factorization during time evolution in high energy states. This feature causes extensive entanglement dynamics. For the SYK the extension is trivial. For the matrix model it works in Fourier space, providing a gauge invariant notion of entanglement entropy associated to a given set of 'generalized free field' modes 5 . If we have a set A of generalized free fields modes O i ω,k with associated number operator n O i ω,k , the entanglement entropy is given by: As for the toy model, the previous analysis shows the existence of three time scales. The first is the well-known local relaxation time scale, associated to quasinormal ringing. It is given by the inverse of the imaginary part of the quasinormal frequency, typically proportional to the inverse of the temperature. The second is the scrambling time scale, which again turns out to be independent of the system size N. It is controlled by the lowest quasinormal frequency. The third is the time by which deviations from equilibrium are uniform throughout the system. This time scale turns out to be: is the imaginary part of the lowest quasinormal frequency. As we expand below, this time scale coincides with the time it takes a freely falling observer to cross the near horizon region until it hits the so-called brick wall [21] or stretched horizon [22], as pointed out in [23,24].
In the last section, we comment on the information paradox [25] in the light of our results. We show that the thermodynamic limit of democratic systems is essentially different from the thermodynamic limit of common local systems, being effectively described by non-unitary dynamics. The large-N limit acts as a coarse-graining by neglecting information originating in subleading correlation functions. It generates entropy without the need of tracing out part of the system. This observation naturally reconciles the old solution described in [25] with unitary microscopic descriptions of black holes in which the interaction structure is democratic.

Non-equilibrium unitary processes with random free fermions
In this section, we study non-equilibrium dynamical processes using the family of Hamiltonians (1.2). Within the context of many-body quantum mechanics, there are many out of equilibrium initial conditions one can consider. Relaxation times, global or local (to be explicitly defined below), depend on the choice of the initial state. In this article, we consider two important classes of initial states. The first initial state we consider intends to emulate the typical black hole thought experiment which consists of 'throwing in' a small subsystem A to a black hole. Subsystem A qualifies as an 'unentangled perturbation' of the black hole state. The combined system is left to relax through unitary evolution. Any information contained in A is shared with the large number of internal degrees of freedom of the black hole and radiation [18,5], and becomes only accessible through fine-grained measurements of the global quantum state. Notice that the previous 'natural' and 'physical' assertion, concerning at the same time both information mixing and information conservation, has proven to be very difficult to observe in a specific model. There exist an obvious tension between the need of complicated dynamics to chaotically mix the initial information and the ability of actually solving the proposed model. As we show below, the family of Hamiltonians (1.2) overcome this tension.
The second initial state corresponds to a 'collapse scenario', in which there is no black hole whatsoever to start with. In this context, an extreme situation arises by choosing as an initial state N/2 unentangled particles, a factorized state belonging to the highest entropic sector of the theory. In this case, the evolution of entanglement entropy is significant for any subsystem size.

"Throwing in" an unentangled particle
For the first out-equilibrium scenario, we will consider the initially unentangled subsystem to be the first oscillator. The initial quantum state is: In the previous expressions, ρ 1 corresponds to the reduced density matrix of the first degree of freedom. It is itself a thermal density matrix with occupation number n. The rest of the degrees of freedom form the 'black hole' and are set in the thermal ensemble ρ β = e −βH /Z, associated to the Hamiltonian (1.2). This thermal ensemble is fully characterized by its mean occupation number f . The computation of the average occupation number f at temperature β for the Hamiltonian (1.2) can be found in [4]. Notice that gaussianity ensures that the two-point correlation function characterizes the full evolution of the state (2.1), including entanglement entropies [26]. Moreover, due to Wick's decomposition principle, stabilization of two-point correlation functions to the stationary value, whatever it is, implies stabilization of all higher point correlation functions as well. All time scales associated to the thermalization process are then encoded in the evolution of the time evolved two-point correlation matrix: The unitary evolution operator is defined as usual by To compute C kl (t), notice that the 'free' nature of the Hamiltonian (1.2) allows an exact solution via diagonalization of the matrix V . If ψ a , for a = 1, · · · , N, are the eigenvectors of V with eigenvalues ǫ a : and we define new creation and annihilation operators d † a and d a by: then the Hamiltonian can be equivalently written as: In the basis d a , the Hamiltonian is a set of decoupled fermionic oscillators. Their time evolution is: Since V is a random matrix, the creation and annhilitation operators create 'random particles'. Their properties can be unravelled by using the theory of random matrices, which deals with the statistical properties of eigenvectors and eigenvalues of matrices such as ηV . For an extensive treatment of random matrices see [13] and [27]. In relation to the eigenvalues, we will only need the widely known Wigner's semicircle law, accounting for the probability distribution of having an eigenvalue equal to λ. It reads: where R 2 ≡ 4Nη 2 , and where we remind that it concerns the eigenvalues of η V , a matrix of size N with deviation equal to η, see the Hamiltonian (1.2). On the other hand, the statistical properties of eigenvectors are less widely known. The main assertion is that the orthogonal matrix of eigenvectors (ψ 1 , · · · , ψ N ) is distributed according to the Haar measure on the orthogonal group O(N). This means that the eigenvectors have independent and random gaussian entries, up to normalization: wherein what follows we will use [p] to denote the average of the random variable p over the random matrix ensemble. With this previous statistical information about eigenvalues and eigenvectors, we can compute the evolution of the correlation matrix C kl (t) for a 'typical' Hamiltonian belonging to the family (1.2). Notice that: the previous equation can be simply written as: The intuition coming from the previous equation is the following. At the initial time S ij (t in ) = δ i1 δ j1 , we recover the unentangled initial state (2.1), as we should. As time evolves S ij (t) → 0, and the final state is the global thermal distribution ρ β = e −βH /Z, characterized by its mean occupation number f . The next step is to compute the statistical properties of C ij (t). In terms of the statistical properties of S ij , they are given by: and (2.14) To compute the mean value [S ij (t)] and the deviations from the mean (σ S ij (t)) 2 = [S 2 ij (t)] − [S ij (t)] 2 we use (2.8) and (2.9). We want to remark that if the deviations from the mean are small (as they will), we are indeed computing the typical properties of single Hamiltonian realizations taken from the family (1.2). This is due to self averaging, a known feature in random matrix theory.
The computation divides into different computations, as shown by the structure of the following matrix: where we mean that C 1j (t) = C 1k (t), C jj (t) = C kk (t) and C ij (t) = C kl (t) for j, k, i, l = 1. Let's start by S 11 (t). The mean is given by: where in the second line we have used the fact that eigenvenctors and eigenvalues are not correlated in the large N limit. The leading term comes from the following contraction , since in this case we do not kill any sum. Noticing that the a = b part of the sum is time independent: The averages over the spectrum are given by: where J 1 (x) is the Bessel function of the first kind. The final result reads: We remind that the Bessel function J 1 (Rt) is an oscillatory decaying function. For large argument R t ≫ 1, the previous expression behaves as: ) . (2.20) By an analogous procedure, we can compute all the other mean values and deviations. The computations are somewhat tedious, especially for [S 11 (t) 2 ], so we will just quote the results in matrix form. The mean is given by: while the deviation reads: Joining all results together, we can write the evolution of the correlation matrix C ij (t) schematically as: (2.23) From these formulas we conclude that the initial out of equilibrium state (2.1) thermalizes and evolves towards a state with homogeneous average occupation given by f . The final state is a 'global black hole'.
Besides, since the system is gaussian, the entanglement entropy of any subsystem A can be computed using the framework developed in [26]. This framework writes the entanglement entropy in terms of the eigenvalues of the correlation matrix of subsystem A. If C ij , for i, j ∈ A, is such a correlation matrix and λ m are its eigenvalues, the entanglement entropy can be written as: The key point to observe now is that off-diagonal terms in the correlation matrix (2.23) vanish in the thermodynamic limit. Consider a subsystem A composed of M oscillators. We can express the correlator matrix as: where n ij (t) is the occupation density of the oscillator i, and ∆C is a random matrix taken from the GOE ensemble with zero mean and deviation σ C ∼ O(1/N 3/2 ). To compute the entropy we need to compute: This can be approximately computed by using the techinques of Appendix A in [20]. The trace is a sum over the eigenvalues of C. These can be expressed as We conclude that the Von Neumann entropy in the thermodynamic limit is extensive at all times, even for big subsystems. For a susbystem B not including the first degree of freedom it is just the thermal entropy of B, independent of time. For a subsystem B ≡ 1 ∪ A including the first degree of freedom, it is just the thermal entropy of A plus S 1 (t). Therefore, in this scenario, the most interesting entanglement entropy is the one associated to the first degree of freedom. It is just given by: Using (2.23) we observe that S 1 (t) evolves towards the thermal entropy of the Fermi distribution, saturating at a distance of O(1/N) from it. Notice that stationarity, the time in which the initial factorization of the density matrix (2.1) is broken, is reached in a time scale t r = 1/R independent of the system size. It is clear from this result that mean field approximations of democratic systems in time dependent scenarios have a very limited range of applicability. We will expand more on this issue and on its applications to black hole physics below. The previous equations (2.21), (2.22), (2.23), and (2.28) constitute the main results of this section. In section (3) we will use them to make several remarks about information spreading in democratic systems, and comment on their applications to black hole physics.

Entanglement dynamics from factorized initial states
The second class of initial states we consider is the class of completely factorized states. More concretely, we can have an initial state in which the firstN oscillators are excited, while the N −N left over are not: This situation emulates a 'collapse' scenario, in whichN initially unentangled particles interact with each other, forming a black hole at late times. To analyze the time evolution we need to write the state in terms of the decoupled oscillator basis d † a : The time evolved correlator matrix is then: The computation is severely more complicated. There is an exponentially growing number of terms, with respect to the number of initial particles. This is due to the vacuum (2N + 2)-point correlation function, after applying Wick's decomposition principle.
Although we were not able to compute the full statistical properties of the previous quantity, it is possible to observe that it is diagonal on average: On average, the only surviving entries are the diagonal ones, corresponding to the expectation values of the number operators. This could have been expected from generic considerations, since the process must be similar to the previous case, in regards to information spreading. The intuition is the following. Given the randomness of V in the Hamiltonian (1.2), from the point of view of a single degree of freedom, the rest of the system behaves as a thermal bath at all times. Relation (2.32) is an example of large-N factorization in fully time-dependent scenarios.
Since there is an effective permutation symmetry between the degrees of freedom in the Hamiltonian (1.2), the decay of the number operator of the initially excited particles is the same for all of them [ c † ↑ c ↑ ](t) = n ↑ (t). The same can be said about the decays of the number operator associated to the oscillators that were initially non-excited [ c † ↓ c ↓ (t)] = n ↓ (t). The average value of the correlator matrix has then the following diagonal form: where there areN entries with n ↑ (t) and N −N entries with n ↓ (t). We will assume in what follows that the off-diagonal deviations σ 2 2 have the same structure as in the previous case, and are of O(1/N p ), for some positive p 6 . This is a natural assumption given the results of the previous section and the intuition coming from large-N factorization in matrix models. At any rate, we stress that the following statements rely on such assumption.
If such assumption is correct, we can repeat the computation of the previous section, formulas (2.25) and (2.27). Considering a subsystem A composed of M oscillators we obtain: From the previous relations we conclude that the deviation from entanglement extensivity is subleading for any M N/2 if p ≥ 2, which is a natural expectation as commented before. Notice that for M = aN, the leading term is proportional to a while the subleading term is proportional to a 2 .
Therefore, up to subleading corrections, the average entanglement entropy of a subset of M = M ↑ + M ↓ oscillators satisfy the following scaling law: is the entanglement entropy associated to one oscillator that was initially excited, and where a similar relation holds for S n ↓ (t). The previous scaling relation was found by generic arguments concerning information spreading in democratic systems in [17]. It seems to be a distinctive feature of democratic systems, when compared to local ones. Relation (2.35) shows that the global properties of information spreading, such as the growth of entanglement entropies of large subsystems are fully controlled by the growth of entanglement entropies of each oscillator separately. But from (2.35) we gain a further insight. Since the single particle entanglement entropies are themselves controlled by the decay of the number operators n(t), the global properties of information spreading are encoded in the local relaxation properties as well.
Given the previous results, it is of obvious interest to find n(t) as a function of theN , in particular its characteristic time scale. This can be computed without too much trouble for the case of one excitationN = 1, in which: 37) decaying to zero at large times. This is expected since the final state has negligible energy. The entanglement entropy associated to the first oscillator grows at intial times, attains its maximum possible value for [ c † ↑ c ↑ ] = 1/2 and then decays to zero for times t ≫ 1/R. Notice again that this simple calculation shows that entanglement is created in times of O(1/R). Mean field methods fail to explain the dynamics at these early time scales.
Also, notice that the simple computation ofN = 1 could have been predicted with the solution of the previous section. Namely, if we take (2.23) and set f = 0 and n = 1 we obtain (2.37) as the leading term in the thermodynamic expansion. This suggests that each degree of freedom sees the rest of the system as a thermal bath already by times t ≪ 1/R much shorter than the relaxation time scale. Looking at the evolved state (2.30), we might expect that the state becomes similar to the random state in theN -particle sector 7 when the phases randomize. This gives a time scale of the order of the inverse of the energy sum in (2.30). In turn, due to the central limit theorem, the energy sum is a gaussian random variable with squared deviation given by σ 2 E =N R 2 . The time scale for dephasing is 7 For a characterization of random states in sectors with a fixed number of particles see [4].
Rt ∼ O(1/ √N ), dying in the thermodynamic limit. It seems natural to expect that at very early times t ∼ R/ √N , each degree of freedom sees the rest of the system as a thermal bath, and, therefore, decays according to (2.23) 8 Finally, notice that the results of this section are not expected to depend on the initial factorized state we consider. For a generic factorized state of N fermions: , and the evolution of entanglement entropy of a set A of degrees of freedom is expected to be generically given by: ) .  .40) show is that these relations are not only valid at stationary states. Due to the democratic structure of interactions, these relations are valid at all times in the thermodynamic limit. 8 The previous assertion needs to be rigorously proven. Indeed, the previous argument might be supporting the claim that the global time scale for relaxation is Rt ∼ O(1/ √N ) in these collapse scenarios. We leave this issue for future work.

Black hole physics, random particles and matrix models
In this section we apply the results of the previous section, concerning the dynamical behavior of random particles, to the physics of black holes and large-N Matrix models. These large-N Matrix models hold key clues towards understanding the emergence of black hole geometric backgrounds [2,3,19], and we would like to understand in more detail several aspects of their color dynamics.
One crucial aspect to analyze in this regard is how perturbations (and the information associated with them) spread through the system as time evolves [5]. A related aspect of primary interest in the community is to study these thermalization processes through entanglement entropy [28,29] since entanglement entropy has a direct geometrical meaning [30]. Finally, it is interesting to see what are the peculiarities of democratic models in regards to the information paradox [25].
In this section we analyze aspects of those three questions. We first do it for the random particle toy model, the main objective being to gain the necessary intuition in a simplified solvable scenario. Besides, as we will see, the results are interesting in their own right. Later on, we show how the intuition naturally extends to large-N matrix models, which furnish exact conjectured descriptions of quantum black holes.

Information transport and scrambling with random particles
In the context of black hole physics, information transport in thermal processes was coined 'information scrambling' in [18,5]. A lower bound was conjectured for such information scrambling in physical systems: where β is the temperature of the system or other significant time scale associated with the Hamiltonian or initial state considered, and N the number of black hole degrees of freedom.
Let's begin with defining in various different ways information scrambling. In [5] it was broadly defined as the time scale by which a given perturbation spreads or contaminates the whole system. In local theories the intuition leads naturally to consider diffusion processes. It was then more rigorously (and abstractly) defined by the 'Page's test' [31], as the time by which subsystems are maximally entangled with the environment, see [5,32,24]. The logic behind this last definition is well known in the context of quantum thermalization, see the recent review [16]. If we choose a generic subsystem A smaller than half of the system, with reduced density matrix ρ A , the non-equilibrium unitary process will drive this reduced state towards the thermal distribution at temperature β. This implies that the entanglement entropy associated to A runs towards the thermal entropy as well. Asking for such an entanglement relaxation to all possible subsystems is a global notion of thermalization. The problem with this abstract notion of thermalization is that its connection to actual 'information transport' is to some extent opaque. This gap was filled in [33,17], in which information transport was carefully defined by the evolution of the structure of correlations between subsystems in the global state. In this setup, the scrambling time is naturally defined by the time scale in which correlations uniformize in size so that information is transported to the whole system democratically. This democratic structure of correlations is a characteristic feature of random pure states [33,20].
We can analyze most types of definitions with the global solution (2.23) 9 . The first thing we observe is that correlations spread instantaneously. From (2.23) it is clear that any property of global relaxation can be characterized by analyzing the number operator associated to the first degree of freedom: In C 11 (t) appears the essential time scale characterizing information spreading, which is given by t relax = 1/R. Indeed, in the thermodynamic limit: For t ≫ 1/R the correlation matrix is simply the global thermal density matrix: up to computed subleading corrections. The information lost by the first degree of freedom is instantaneously spread through the whole system, since correlations are uniform through the system C 1j (t) ∼ C 1k (t). This is of course due to the microscopic non-locality. The fact that correlations are at all times uniform over the system does not immediately imply scrambling. A clear and extreme example of this situation arises when all correlations are zero, as happens for the initial state considered before: where there is a clear factorization between the first degree of freedom and the thermal ensemble. Scrambling of the information associated with the first degree of freedom requires correlations to be uniform over the system plus the first degree of freedom to be maximally entangled with the environment. Since any created correlations are homogeneously spread through the system, the time of entanglement production, or analogously the time in which the initial factorization (3.5) breaks down is the scrambling time itself. The entanglement entropy was computed in the previous section to be: which takes the Von Newmann entropy at the initial time: to the Von Newmann entropy associated to the global thermal state: The characteristic time scale for stationarity of entanglement entropy is thus t local = 1/R. Also, as shown in the previous section, notice that any other entanglement entropy associated to other subsystems reach the entanglement plateaux at this time scale as well.
At this point, we want to remark that this result, if correct, furnish a counterexample of the mean field bound presented in [32]. In Ref. [32] it was claimed that there is a lower bound for the time scale of entanglement production in gaussian systems with fully non-local interactions. This lower bound was claimed to be proportional to log N. More specifically, it is claimed that for models such as (1.2) the time evolution of the N oscillator system can be approximated by: for some suitably defined 'mean field local unitary evolution' U MF i (t), acting in each oscillator separately. The claim is that this approximation is valid until times proportional to log N. If we would apply such approximation to our initial state (2.1), it would imply that until such large times entanglement could not be created. We believe that our results, formulas (2.23), (2.28) and (2.37), which follow from relatively simple computations, show that production and stationarity of entanglement entropy are attained in a time scale independent of the system size. If correct, they furnish a counterexample to the arguments presented in [32], since this early time entanglement cannot be explained with a mean field approximation, such as (3.9). It is obvious that this issue deserves further consideration. Finally, in view of our results, notice that there is indeed some mean field approximation applying here. In the large-N limit, and before times of O(β log N), the evolution of the global state can be approximated by: where ρ(t) is not obtained from ρ(0) by unitary evolution but by some super-evolution operator which change the on-site Von Neumann entropy of the reduced on-site state. This is another path to understand (2.39) and (2.40).
Coming back to scrambling, and to make the statement even clearer, we can also use the approach to information scrambling proposed in [33,17]. Due to the effective permutation symmetry of the random free fermions Hamiltonian (1.2), for any pair of subsystems A and B with the same size, the Mutual Information between the first degree of freedom and those subsystems is equal: For a subsystem A of M degrees of freedom, using (2.27) we obtain: We now join all features together. The information lost by the first degree of freedom, quantified by ∆S 1 (t), satisfies two properties. First, it is non-locally spread through the whole system, the Mutual Information being of the same size for any given pair of subsystems with the same size (3.11). This is due to the democratic structure of interactions in the microscopic Hamiltonian (1.2). Second, this information can only be recovered by looking to big enough subsystems, since for small enough subsystems the Mutual Information vanishes in the thermodynamic limit (3.12). It is of obvious interest to quantify how big enough the subsystem must be to recover the information about the first degree of freedom. However, we want to stress that subsystems with a size not scaling with N in the thermodynamic limit will not achieve this goal. From this precise point of view, the information lost by the first degree of freedom is fully scrambled over the environment. Given (3.6), by t local = 1/R the first degree of freedom has shared an O(1) amount of its information content with the environment. Therefore, after t local = 1/R has elapsed, an O(1) amount of information has been shared non-locally with an O(N) amount of degrees of freedom. We can repeat all the arguments for the factorized initial state (2.29). The only thing we need to assume to arrive at the same set of conclusions is again that the mean squared deviation of the correlation matrix is of O(1/N p ), for p ≥ 1. As explained in the previous section this is a natural assumption given the large-N matrix structure of the model and given the exact results obtained with the first initial state. If such assumption is correct, entanglement entropies evolve extensively, even for subsystems with extensive size M N/2. This implies that saturation at the plateaux of the subsystem occurs when one degree of freedom saturates by itself, i.e the time in which the on-site entanglement is created. This is equivalent to the time by which the mean field approximation breaks down, which is just given by the local relaxation time scale.
Summarizing, from the behavior of entanglement entropies, Mutual Information and correlation functions, we conclude that the time scale associated to: • Information lost by on-site degrees of freedom, quantified by their associated entanglement entropy, • Information spreading through the black hole degrees of freedom, characterized by correlations and Mutual Information being uniform through the system, • Relaxation of all entanglement entropies to thermal values at leading order, • Breakdown of the mean field approximation is size independent and given by: We conclude that the family of Hamiltonians (1.2) studied in this article furnish an specific example of the generic statement described in [17]: • For democratic systems, information spreads instantaneously. The characteristic time scales appear already at the level of local relaxation of perturbations.
On the other hand, these results challenge the fast scrambling conjecture [5]. The fast scrambling conjecture states that the time scale for a global spreading of information is bounded from below by t scrambling β log N, while here we find a time scale independent of the system size. We want to stress again that the lower bounds presented in Ref. [32] seem not to be valid generically. The breakdown of the mean field approximation, or analogously the breakdown of the factorization of the initial state (3.5), occurs on a time scale given by t = 1/R, independent of the system size, see (3.6). In such a situation, given that the spreading of information is instantaneous, the model (1.2) seem to furnish a counterexample of the claimed fast scrambling lower bound.

Scrambling vs randomization time scales
In the previous section we have shown that an O(1) amount of information about the initial state is spread in a time of O(1/R) over the entire system. But has the thermalization process reach equilibrium in all possible senses? Looking to the exact form of the solution (2.21), (2.22),(2.23), we see that after the relaxation time scale O(1/R) has elapsed, there is still a long way towards 'global equilibration of the deviations from thermality'. To be precise, the set of number operators corresponding to the black hole internal degrees of freedom (the oscillators with i = 1), are shifted from the thermal expectation value by: (3.14) already at times of O(1/R). On the other hand, the deviation from thermality of the number operator corresponding to the first degree of freedom is: It is now natural to define a 'randomization' time scale t random , by asking the deviations from thermality to be uniform over the system. Notice that this is not equivalent to ask for correlations to be uniform over the system. From a more physical perspective, it is defined as the time scale at which the leading term in (3.15), the term driving the decay, ceases to be valid due to O(1/N) corrections. For random free fermions we obtain: Summarizing, from the present perspective there are three natural time scales one can consider: local relaxation time scales, information transport or scrambling, and near equilibrium randomization. For usual local systems, scrambling and randomization are equal, since the time that information takes to spread through the entire system is sufficiently big due to causality constraints. For non-local systems, such as random free fermions, the situation is different. Since information transport is instantaneous, local relaxation and scrambling are governed by the same time scale, while randomization takes a longer time.
In the next section we extend all this intuition to large-N matrix models.

Information transport and scrambling with matrix models
Let us recapitulate the important results obtained so far. In the large-N limit of the random free fermion model, the correlator matrix is diagonal, see (2.23) and (3.3). The deviation from this diagonal behavior can be computed. This large-N factorization happens even in time-dependent scenarios. It is rooted in the democratic structure of the microscopic Hamiltonian (1.2). Since the diagonal entries are given by the local expectation values of number operators c † i c i (t) = n i (t), the entanglement entropy of a subset A of M N/2 degrees of freedom is simply given by: For random particles, the entanglement entropy is an extensive quantity even in timedependent scenarios. The deviation from extensivity can be computed with the deviation of the correlator matrix from its diagonal mean. This computation shows that the deviation is subleading in the thermodynamic limit, and allows the study of entanglement in every subsystem, no matter its size, by the entanglement of on-site subsystems. As commented in the introduction, the generalization of the Hamiltonian (1.2) to the bosonic case is straightforward. The only change is the connection between entanglement entropy and occupation densities, the formula being slightly different. For random free bosons it would read: We claim that the previous formulas extend to large-N matrix models and the SYK model [12]. The reason is that this formula is only based on the effective gaussianity of reduced subsystems in democratic models. In other words, formulas (3.17) and (3.18) are information theoretic versions of large-N factorization and therefore apply to any theory satisfying such property (such as large-N matrix models and the SYK model). The only difference between the toy model and real microscopic descriptions of black holes lies in the specific time dependence of occupation densities. For random particles, we have a polynomial decay of number operators (3.15) while for black holes we expect quasinormal ringing [34]. Therefore, for the SYK model the extension is trivial. Since the Hilbert spaces are isomorphic, the nature of subsystems is the same. Entanglement evolution is just given by (3.17).
For large-N field theories one has to take care due to mainly two reasons. The first is that the democratic interactions appear only in color space. The second is that subsystems are not defined by partitions of the Hilbert space, since these partitions are not gauge invariant. Subsystems are better defined by operator algebras themselves. To avoid these two issues, the easiest path is to consider the theory in Fourier space. Our subsets will be sets of generalized free field modes 10 O ω,k . Generalized free fields are traces of products of the fields of the theory, in which we are not allowed to consider products with O(N) fields. Generalized free fields are gauge invariant operators by construction, due to the trace, a feature that avoids the problem of defining entanglement in a gauge invariant manner. They are called 'generalized free fields' because they generate a Fock space of free particles in the large-N limit, but they accomplish so without satisfying any free wave equation. This last feature is possible due to large-N factorization [1]. At finite temperature, for any set of generalized free fields O, large-N factorization is the following statement 11 Tr( ρ β O(π 1 )O(π 2 )) · · · Tr( ρ β O(π 2n−1 )O(π 2n ))+O(1/N) .
(3.19) As shown in [19], relation (3.19) allows us to define creation and annihilitation operators O ω,k and O † ω,k satisfying the usual algebra of thermal free oscillators. Assuming for concretenes a bosonic generalized free field we have: together with where connected higher point correlation functions vanish. This implies that the two point correlation matrix C ωk,ω ′ k ′ is diagonal up to subleading corrections, sclaing as O(1/N). This whole discussion does not change if the quantum state ρ is not the thermal density matrix, but an evolving pure/mixed ρ(t) 12 : So just by assuming (3.19), and repeating the argument of the previous section to bound the errors when computing entanglement entropy, the time dependent and gauge invariant entanglement entropy of any set A of M N/2 generalized free fields is just given by: Notice that due to large-N factorization we are able to define entanglement in a gauge invariant manner. The clue is to consider entanglement entropy to be associated to observables or operator algebras 13 , such as generalized free fields O. Given an operator algebra, there are many observations one can do (all possible correlation functions). Finding entanglement entropy from all correlator functions in the operator algebra can be quite complicated in general. Crucially, for large-N matrix models the situation simplifies due to large-N factorization, and entanglement entropy can just be found by looking at the covariance matrix. This covariance matrix is diagonal in Fourier space and it is controlled only by occupation densities (one point functions), see (3.24). At it should, formula (3.24) is coherent with the expectations for random pure states, where it correctly provides the thermal entropy at the appropriate temperature: Besides, formula (3.24) implies that the conclusion obtained for random particles extends to matrix models: • Entanglement entropies are fully controlled by the decay of occupation densities. Global scrambling of information is controlled by local relaxation of perturbations.
We want to stress again at this point that the only assumption we used to derive this result is large-N factorization in thermal ensembles (3.19). This is well supported by the toy model of random particles and by the black hole geometric description of large-N matrix models. Given this feature, reduced subsystems are products of uncorrelated gaussian density matrices. We can then use the covariance matrix approach for gaussian systems [26], since this is the only variable surviving the large-N limit. Finally, bounds on subleading terms can be found as in the previous section, from the size of the corrections to the two-point functions. Formula (3.24) is a remarkably simple result with many potential consequences. In what follows, we will just focus on its implications to information transport, scrambling, and randomization. To this end we first need to argue for the actual functional form of n O i ω,k (t) in large-N matrix models. Notice that on generic grounds, quasinormal ringing implies that the field mode decays at late times as: where Ω I ω,k and Ω R ω,k are the real and imaginary parts of the quasinormal frequency Ω ω,k . For a harmonic oscillator satisfying (3.22), the number operator is basically the square of the field, so that we expect: Although the argument leading to (3.27) is heuristic, a similar law has indeed been found recently in [36] for a related occupation density described first in the context of holography in [37]. More importantly, we have recently verified that this type of law describes the evolution of Bekenstein-Hawking entropy [38]. Plugging the previous expression (3.27) into (3.24), we obtain the evolution of entanglement entropy of a subset A of generalized free fields. The formula is extensive, but each mode decays with a different Ω ω,k . To analyze its properties we just need to analyze the evolution of single perturbations. Defining f ω,k (t) ≡ A ( cos(2 Ω R ω,k t +θ)+1 ), the evolution of entanglement entropy associated to a single mode is given by: We cannot trust the previous solution to small enough times since in that regime we should include higher quasinormal frequencies into (3.27). But notice that for 2Ω I t ≪ 1, the difference between the entanglement entropy at time t and the one at time t = 0 grows linearly with time: On the other hand, for large times 2 Ω I t ≫ 1, where formula (3.27) can be trusted, we enter the so-called entanglement plateaux, which is predicted to be: Relations (3.28) and (3.30) are analogous to relations (2.39) and (3.6) for the case of random free fermions, and are part of the main results of the article. The difference between the toy model and the real one lies only in the functional form of the decay of occupation densities. A similar law was recently found to describe Bekesntein-Hawking entropy evolution, see Ref. [38]. Finally, just for completeness, notice that the Mutual Information between gauge invariant generalized free field modes is again zero at all times: due to the extensivity of entanglement entropies through unitary evolution. This is an information theoretic expression of large-N factorization. For the Mutual Information between bigger subsets of generalized free field modes one can use (3.24). Joining all results together, the conclusions regarding scrambling and local relaxation are universal to all democratic models, including SYK and large-N field theories. The time scale associated to: • Information lost by the generalized free field mode O ω,k , quantified by its associated entanglement entropy, • Information spreading through the matrix model, characterized by the vanishing of connected correlators and Mutual Information between small subsystems at all times, • Relaxation of all entanglement entropies to thermal values at leading order, • Breaking of the mean field approximation, characterized by the growth of entanglement of one generalized free field mode.
is operator dependent, size independent and given by: a result which put under question marks the fast scrambling conjecture [5]. We want to stress again that although some aspects seem complicated, everything is based on large-N factorization, which implies factorization of reduced subsystems, and gaussianity. Information theoretically, large-N factorization translates into equations (3.18), (3.17) and (3.24).

Scrambling vs randomization time scales
We have shown that random particles, SYK and large-N matrix models behave very similarly in regards to information spreading. Entanglement entropy evolution is extensive for all of them, and it is fully controlled by the decay of occupation densities. This is rooted in large-N factorization. What differs in each model is the time dependence of the occupation densities themselves. Quite generically, we expect the occupation densities to decay exponentially fast, while for random particles we found only a polynomial decay 14 . This difference can be seen through the randomization time. This was defined as the time scale by which deviations from thermality are uniform throughout the system. It can be found by looking at the properties of entanglement entropies at the entanglement plateaux, or equivalently by looking at the occupation densities. Physically, it is the time scale by which we cannot trust the time decay, since O(1/N) corrections kick in. For large-N matrix models, it is just given by: Taking into account the geometric observations done in [23,24] concerning black hole geometric backgrounds, this randomization time scale is of the same order of magnitude as the flight crossing time through the near horizon region until the so-called brick wall [21]. Our results suggest that dynamical processes in the near horizon region correspond to near equilibrium evolution, fully controlled by quasinormal frequencies, and cut dynamically at the time scale (3.34). The remarkable point is that the saturation of the process at such randomization time scale is not a phenomenological assumption, unlike in the gravity side. In the present perspective, the saturation at (3.34) is a purely dynamical feature, due to finite N effects coming from unitarity of the microscopic model. Notice also that we expect a different randomization time for each perturbation. The emergence of the brick wall might be not universal after all. In this framework, it definitely depends on the specific observable we use to probe the black hole.

Comments on unitarity breaking in the thermodynamic limit
In this section we comment on another interesting and distinguishing feature of democratic systems. Going to (2.23), and taking the large-N limit we obtain: Since the resulting correlation matrix is diagonal, the Von Neumann entropy associated to it is given by: a time-dependent function. In the large-N limit, the entropy associated with the global state changes with time. Its evolution cannot be driven by a unitary operator, and should be represented by some super-evolution operator $, in the lines of [25]. In this precise sense, unitarity is globally broken in the thermodynamic limit. We stress that this is essentially different from local theories, where superoperators and non-unitary behavior appear only locally when tracing out some part of the system. Taking the thermodynamic limit N → ∞ does not affect the entropy globally as time evolves in local systems. We conclude that there is in principle no contradiction whatsoever between unitarity of quantum gravity and non-unitarity of the Einstein-Hilbert action (its semiclassical limit), whenever the microscopic theory is democratic, such as Matrix Models [2], AdS/CFT [3], or the recently introduced infinite range fermionic model [12].

Conclusions
In this article we first have proposed and studied a quantum toy model of black hole physics. While capturing essential features of quantum black holes, the model is still exactly solvable. The two properties we wanted our system to possess are quantum thermalization and democracy of interactions. As described in the introduction, 'democracy of interactions' means the strongest form of non-locality, which amounts to have every oscillator interacting with every other oscillator through couplings of the same size. The first requirement is an obvious necessary feature since there are very strong reasons to consider black holes as many body systems at finite temperature [9,10]. The second property intends to emulate realistic black hole models, such as large-N matrix models [2,3,5,12]. The color sector physics of large-N matrix models furnish an example of such democratic systems. The proposed toy model was defined by the Hamiltonian (1.2). It is a gaussian system of coupled fermionic oscillators in which the mass matrix belongs to the Gaussian Orthogonal Ensemble of random matrices 15 . This system was described in [4], as an analytical example of the Eigenstate Thermalization Hypothesis [14,15,16]. The advantages of choosing such a model are the following: • The model is exactly solvable. Corrections of O(1/N) can be properly taken into account.
• The model has enough complexity to display quantum thermalization. It allows a precise study of information spreading in democratic systems, at the global and local levels.
• The extension to bosons is straightforward.
Within the toy model, we have chosen to analyze two generic physical processes. The first scenario is the usual black hole thought experiment, in which one throws in a particle initially unentangled with the black hole (2.1). Subsequent evolution entangles both particle and hole, and information about the initial state n is mixed/spread through the system (but conserved due to unitarity). This unitary time evolution can be fully accounted by the two-point correlation functions of the theory, equations (2.21), (2.22),(2.23), and (2.28). The second scenario is akin to a 'collapse' scenario, in which the initial state is a factorized many-particle state (2.29). For this case, the computations are harder, but we could show that the typical correlation matrix is again diagonal (2.33).
As described in the second part of the article, this type of covariance matrix appears also in realistic models of black holes, such as SYK or large-N matrix models. The property that underlies this covariance matrix is large-N factorization, a property characteristic of both SYK and large-N matrix models that seems to be a defining feature of democratic models.
Given large-N factorization, extensivity of entanglement evolution follows naturally, with subleading computable corrections to it. The formulas for all the models are equivalent, and given by (3.17) for random free fermions or the SYK model, by (3.18) in the random bosonic cases and by formula (3.24) for matrix models, where it is more convenient to work with generalized free field modes in Fourier space. An intuitive explanation of the extensivity of entanglement evolution is the following. Due to the non-locality of the theory, correlations between a given degree of freedom and any other are all of the same size. Since there are O(N) degrees of freedom interacting at the same time and given monogamy of entanglement, the shared correlation between any pair of them must be of O(1/N). Large-N factorization seems thus to be related with monogamy of entanglement in democratic systems.
In every case, the main conclusions are: • In the thermodynamic limit, the behavior of entanglement entropies is extensive at all times, see (3.18) and (3.17). The evolution of entanglement entropies is fully controlled by local relaxation of occupation numbers.
• Information spreading is instantaneous. Any lost information is instantaneously scrambled since the Mutual Information between small enough subsystems is zero at all times.
• The time scale by which the entanglement entropy of a subset A reach the entanglement plateaux is independent of the system size and it is controlled by on-site quantities. The size A can scale with the size of the total system, as long as it is not bigger than half of it.
• The mean field approximation, defined as the time of on-site entanglement production is broken at a time independent of the system size, although another kind of mean field approximation is valid for longer times (3.10).
The previous results challenge the fast scrambling conjecture [5]. As described in section (3), this conjecture states that there is a minimum lower bound for the time scale associated to 'information scrambling'. The conjectured lower bound is of order O(β log N). Information scrambling was defined by global equilibration of entanglement entropies, together with information spreading through the entire system. The previous results suggest that both requirements are fulfilled in a time t scrambling = 1 2 Ω I , where Ω I is the lowest quasinormal frequency associated with the subset A that has been brought out of equilibrium. Although there is still room for discussion of the present results, and probably subtle aspects concerning information theory might appear given that we are dealing with very unconventional Hamiltonians, the present results seem to be making a strong case against the fast scrambling conjecture and they invalidate the mean field bounds presented in [32].
We conclude that the only difference between the different models lies in the specific time dependence of number operators. Whereas for random particles we found polynomial decays, for matrix models we have exponential decays. Therefore, the time scale by which O(1/N) corrections need to be taken into account is different in each case. This time scale was coined 'randomization time'. It is the time scale by which the deviations from thermality are uniform over the system. For random particles, we found a randomization time scaling polynomially with N. For matrix models, the result turns out to be: Interestingly, this time scale sets the causality bounds for crossing the near horizon region until the so-called brick wall or stretched horizon, see [23,24]. It is tempting to conclude that the emergence of the near horizon geometry is related to near equilibrium evolution, characterized by the properties of the entanglement plateaux, or analogously by the properties of the occupation densities. These two quantities are controlled by the famous quasinormal modes, so it is tempting to conclude that the emergence of near horizon regions is fully encoded in the quasinormal ringing of the dual field theory, i.e in the poles of retarded correlation functions [34,40]. This speculative proposal provides a dynamical reason for the appearance of the brick wall, as the time by which the evolution of occupation densities and entanglement entropies is affected by O(1/N) corrections. In the last section, we made a remark about the information paradox [25] in the light of the toy model. We noticed that unitarity of the toy model, characterized for example by a time independent Von Neumann entropy of the global state, is lost in the large-N limit at all times. Taking the large-N limit affects the entropy globally. This is a characteristic feature of democratic systems, whose semiclassical approximation at high energies is effectively described by non-unitary dynamics. The old and famous proposal of [25] seems to naturally emerge in theories of quantum gravity which have a democratic structure of interactions.