Spacing Statistics of Energy Spectra: Random Matrices, Black Hole Thermalization, and Echoes

Recent advances in AdS/CFT holography have suggested that the near-horizon dynamics of black holes can be described by random matrix systems. We study how the energy spectrum of a system with a generic random Hamiltonian matrix affects its early and late time thermalization behaviour using the spectral form factor (which captures the time-dependence of two-point correlation functions). We introduce a simple statistical framework for generating random spectra in terms of the nearest neighbor spacing statistics of energy eigenvalues, enabling us to compute the averaged spectral form factor in a closed form. This helps to easily illustrate how the spectral form factor changes with different choices of nearest neighbor statistics ranging from the Poisson to Wigner surmise statistics. We suggest that it is possible to have late time oscillations in random matrix models involving $\beta$-ensembles (generalizing classical Gaussian ensembles). We also study the form factor of randomly coupled oscillator systems and show that at weak coupling, such systems exhibit regular decaying oscillations in the spectral form factor making them interesting toy models for gravitational wave echoes. We speculate on the holographic interpretation of a system of coupled oscillators, and suggest that they describe the thermalization behaviour of a black hole geometry with a membrane that cuts off the geometry at the stretched horizon.


Background and Context
Black holes serve as an important testing ground and challenge for our understanding of the quantum nature of gravity. The application of quantum field theory on black hole backgrounds leads to the conclusion that pure states can evolve into mixed states through the mechanism of Hawking radiation [1,2]. On the other hand, theoretical developments in string theory and AdS/CFT suggest that quantum gravity should be a unitary theory. The conflict with the expectation that quantum gravity should be a unitary theory and the apparent violation of unitarity by Hawking's calculation has given rise to the black hole information problem [3][4][5]. Explorations into the information problem have led to the idea of quantum gravity effects modifying the near horizon description of black holes, most notably in the form of some kind of microstructure near the horizon (e.g. fuzzballs, firewalls, non-locality etc.) [3,[5][6][7][8][9][10][11].
The idea of quantum effects modifying the near horizon physics of a black hole have recently also been considered from an experimental perspective. If quantum effects near the horizon give rise to the possibility of partial reflection rather than complete absorption of a perturbation then it could give rise to "echoes" in the ringdown behaviour of a black hole after being perturbed. In particular, such echoes may be detected by LIGO in gravitational waves emitted during black hole binary merger events [12][13][14][15][16][17][18][19][20][21][22]. In most studies, these modifications are modelled by cutting off the black hole geometry near the horizon. At the cutoff, semi-reflective boundary conditions are enforced and perturbations on the cutoff background are studied 1 [23][24][25][26][27][28]. Studies that do this make an implicit assumption that the quantum effects that encode the unitary nature of black holes are "localized" near the horizon and that they effectively manifest themselves as a cutoff near the horizon with non-standard boundary conditions. Although such models are completely classical, they highlight the interesting possibility of finding observable imprints of horizon scale quantum gravity effects in the ringdown behaviour of black holes.
In this paper, we will explore the idea of echoes in the context of the thermalization behaviour of quantum chaotic systems. The motivation behind this study is the idea that the near horizon dynamics of a black hole is generally captured by quantum chaotic degrees of freedom. When a perturbation begins to probe and excite the microstructure near the horizon, it could lead to observable deviations from the standard ringdown as predicted by quasi-normal mode decay. If the microstructure effectively behaves as a semi-reflective cutoff near the horizon, then one would get echoes in the unitary description of the black hole. The goal of this paper is to explore the properties a quantum chaotic system, if it were to exhibit echoes in its thermalization behaviour. To aid in this exploration, we will make use of recent advances and tools in AdS/CFT that have placed interesting constraints 1 In these studies the time scale over which the echoes manifest after the black hole is perturbed is called the "echo time scale." It depends on how close the cutoff is placed relative to the horizon. For a cutoff placed a proper radial Planck length from a Schwarzschild black hole it can be shown the echo time takes the form of the scrambling time β ln(SBH ) [23] (e.g. for a 10 solar mass Schwarzschild black hole the scrambling or echo time is of the order of 0.2 secs.). on the energy spectrum of black holes. In the next subsection, we will give a brief overview of these developments and tools.

Black Holes as Quantum Chaotic Systems
The AdS/CFT correspondence, which was first formulated by Maldacena [29] suggests that (quantum) gravitational systems in AdS have a non-perturbative formulation in terms of a strongly coupled CFT [30][31][32]. This in turn provides a way to explore various aspects of the quantum nature of black holes from the perspective of a thermal CFT [3,5,7,33,34].
One of the earliest explorations into the quantum spectrum of a black hole from the perspective of AdS/CFT involved the thermofield double (TFD) state [33]. The TFD is a state which is conjectured to describe a two sided eternal black hole in AdS and is written as: |T F D = 1 Z(β) n e −βEn/2 |n L ⊗ |n R , (1.1) the state above describes two identical systems (one living on the left boundary and the other on the right) which are entangled. The states |n L,R are eigenstates of a Hamiltonian with eigenvalues E n , and the quantity Z(β) is the partition function. Using this, one should expect the decay of quasi-normal modes of the dual black hole to be consistent with the decay of two-point thermal correlators. It was noted in [33], that the spectrum of a CFT which is dual to a black hole (with a spherical horizon) lives on a compact space and as a result its spectrum will be discrete. The discreteness of the spectrum implies that correlation functions cannot decay to zero which is an indication of unitary evolution in the bulk [4,33,35,36]. The inability for correlators to decay to zero comes as a surprise from the bulk perspective because quasinormal modes for an AdS black hole decay exponentially to zero [34]. It was argued in [33] that correlators will decay in a manner consistent with the semi-classical calculation of quasi-normal modes of the black hole until the correlator roughly becomes of the order e −S , where S is the entropy. Afterwards, the discreteness of the spectrum will become important and the semi-classical calculation will no longer agree with the boundary calculation of the two-point function. In a semi-classical context, it was shown that to prevent the continual decay in the bulk, one can include other saddle point geometries in the Euclidean path integral approach [33]. It was also demonstrated that one can avoid continual decay of the two-point function and make the spectrum discrete by cutting off the spacetime near the horizon [35][36][37].
Quantum chaos is another important feature of thermal CFT systems that are dual to black holes. In particular, CFT systems that are dual to black holes are conjectured to be maximally chaotic 2 [42]. Quantum chaos manifests itself in the spectrum of a quantum system through nearest neighbor eigenvalue "repulsion." In particular, it is conjectured that many statistical properties of the spectrum of a quantum chaotic system follow the spectrum statistics of random matrices drawn from a suitably chosen ensemble [43][44][45]. The connection between gravity and random matrix theory has been further explored in the context of 2D Jackiw-Teitelboim (JT) gravity by showing the partition function of JT gravity can be expressed as a random matrix integral [46][47][48][49]. These recent results corroborate the conjecture that random matrix theories are capable of describing the statistical properties of the spectrum of quantum chaotic systems such as black holes.
A particularly useful observable that has been computed in the context of quantum chaos, is the normalized spectral form factor which is defined by: , (1.2) where Z(β) is the partition function in the canonical ensemble with inverse temperature β. It is a useful quantity since it can be used to diagnose the discreteness of a spectrum.
In the context of the TFD state, if one considers the two point function of a Hermitian operator of the form I R ⊗ A (where I R is the identity on the right boundary and A is a Hermitian operator living on the left boundary) we have: we can see the time dependence (up to matrix elements which will depend on the energy) of the form factor and the two point function are governed by the details in the energy differences in the spectrum. If the matrix elements are smooth functions of energy that vary slowly, we can use the form factor as a proxy for understanding how a perturbation thermalizes for a TFD state. This is interesting because by identifying a Hamiltonian and its spectrum of states we can construct the associated TFD state. If we conjecture that the TFD is dual to a two-sided black hole geometry, then we can roughly interpret the early time behaviour of the spectral form factor as a description of the ringdown of the black hole after being perturbed (i.e. the quasi-normal modes). From this perspective, different choices of Hamiltonian give different types of ringdowns of the conjectured dual black hole. Usually, at early times, the discrete spectrum can be replaced by a coarse grained smooth density -and one initially sees decay in the form factor (controlled by the decay of quasinormal modes of the black hole). At very late times, the behaviour of the spectral form factor is governed primarily by the small energy differences between nearby eigenvaluesand in general it produces very erratic oscillations which never decay to zero.
In the context of a random matrix ensemble (RME), one is interested in the normalized averaged spectral form factor (throughout this paper we shorten the name and call it the "form factor" unless otherwise stated) given by: where · RM E is an average over a random matrix ensemble 3 . One can compute this numerically, by taking an average over many samples of the form factor which is constructed using the eigenvalues of sample matrices drawn from the ensemble. In the context of a random spectrum, the early time behaviour of the form factor has the property of being "self-averaging" which means that the form factor of a single sample is close to the average over many samples. However, at late times the form factor of a single sample will deviate significantly from the average and will no longer be self-averaging. In the case of classical Gaussian ensembles whose averaged form factor was studied in [50], it was shown that the averaged late time behaviour after the initial decay took the form of a ramp followed by a plateau(e.g., see Figure 1). The existence of the ramp can be attributed to the near- We plot the infinite temperature form factor for a spectrum generated by 100 × 100 matrices pulled from the Gaussian unitary ensemble (GUE). The blue line is the form factor for a single sample matrix from the GUE and the yellow line is the 500 sample average. One can see that at late times, the averaged form factor exhibits a linear ramp followed by a plateau.
est neighbor spacing statistics of the eigenvalues of the matrices drawn from the classical Gaussian ensembles. Furthermore, the spectrum of the SYK model (which is a toy model for near extremal black hole physics) had similar features to the spectrum of random matrices belonging to classical Gaussian ensembles 4 . These findings led to the idea that large AdS black holes are described by theories that have a spectrum resembling the spectrum 3 Note that we are taking the average over the numerator and denominator separately. This particular way of averaging is called "annelled" disorder averaging and this is in contrast to "quench" disorder average which is obtained by doing an average of the entire normalized expression for the form factor. At infinite temperature annealed and quenched averages are the same. 4 In Section 3.3 we will discuss how the SYK model form factor is differs in certain details with the form factor associated with matrices pulled from classical Gaussian ensembles.
of random matrix theories. In particular, the results reviewed in this subsection suggest that black holes are dual to quantum systems that have a discrete and chaotic energy spectrum 5 . In light of this, it is natural to ask if it is possible for a quantum chaotic system to exhibit echoes in its thermalization behaviour. If the answer is yes, then one might view such a system as potential candidate for describing black holes with microstructure near the horizon 6 . In our exploration we will primarily be interested in how various aspects of the spectrum of a system manifest in the early and late time behaviour of the form factor and use this as a way to understand how perturbations to such systems thermalize. In the context of black hole echoes, we will be interested in finding a system with a quantum chaotic spectrum that gives rise to regular but decaying oscillations in the form factor which might be interpreted as echoes being generated by some microstructure near the horizon of a black hole in the bulk.

Overview of the Paper
In Section 2, we will introduce a model for a random spectrum and derive closed form expressions for the averaged spectral form factor in terms of the nearest neighbor spacing statistics. In particular, Section 2.1, will go over the details of the model we adopt to describe a random spectrum and also go over the basics of computing various averaged quantities in such a model. In Section 2.2, we use our model to derive a closed form expression for the averaged spectral form factor in terms of integrals which involve the nearest neighbor spacing (NNS) distribution.
In Section 3, we consider various common NNS distributions and compute the associated averaged form factor using the formulas we derived in Section 2. In Section 3.1, we analyze the form factor of a non-random evenly spaced spectrum and show the form factor is periodic with a period that is inversely proportional to the spacing between the energy levels. In Section 3.2, we analyze the form factor associated with having a Poisson NNS distribution (which would serve as a model for the spectrum of a generic quantum integrable system). We show that the averaged form factor monotonically decreases and saturates to a non-zero value at late times. In Section 3.3, we compute the averaged form factor of a system which has NNS statistics following the Wigner surmise (which serve as a model for the spectrum of a quantum chaotic system that follows the statistics of classical Gaussian ensembles). We find that the form factor has an initial dip followed by a ramp and plateau, which is consistent with expectations from the studies of classical Gaussian ensembles, as shown in Figure 1. In Section 3.4, we introduce the gamma distribution as a simpler alternative to the Wigner surmise NNS distribution. We show that, in the appropriate regimes, the form factor associated with the gamma distribution will still have a ramp and plateau. 5 It should also be noted that there is an additional constraint that the system should saturate the Lyapunov bound discussed in [42]. We will address this point in the discussion of echoes in quantum chaotic systems at the end of the paper. 6 Here we are specifically referring to microstructure that gives rise to echoes. It is possible that a black hole could have microstructure near its horizon but not exhibit echoes when perturbed. We will discuss this possibility in more detail at the end of the paper.
We also show that there are other regimes in which the NNS distribution contains a large "gap." We show the existence of the gap results in regular decaying oscillations in the averaged form factor after the initial dip. We argue that similar oscillations in the form factor will occur for large values of the Dyson index in the Wigner surmise distribution (Equation 3.12). These larger values of the Dyson index can be obtained from generalizations of classical Gaussian ensembles and are commonly refereed to as "β-ensembles," which are know to arise from tri-diagonal random matrices [51].
In Section 4, we study the spectrum of a many-body system composed of identical oscillators coupled to each other by a matrix belonging to a classical Gaussian ensemble. In the weak coupling regime, we study the role the chaotic interactions play in the splitting of degenerate energy states. We show that within each degenerate sector, the splitting of the energy levels will generally give rise to chaotic spectrum statistics within each sector. We verify this by numerically studying the form factor and spectral density in the weak and strongly coupled regime. We show that the averaged spectral form factor of such systems at weak coupling exhibit regular decaying oscillations at early times followed by a ramp and plateau at late times -which is consistent with random matrix theory models that describe black holes. We then propose such systems of weakly coupled oscillators as toy models for describing a membrane which gives rise to decaying echoes. The echoes repeat at approximately the fundamental frequency of the oscillators. We identify the parameters in the oscillator system with those of the bulk geometry+membrane system. Through this identification, we show that the question of whether the oscillator system is strongly or weakly coupled depends on how far the membrane is placed from the horizon.
In Section 5, we conclude by summarizing the major findings and discuss future research directions.

A Simple Model for Random Spectra
As stated in the introduction, quantum chaotic systems are generally conjectured to have spectrums that follow the statistics of random matrix theories (most notably Gaussian ensembles). One of the most important pieces of statistical data in the spectrum of a random matrix ensemble is the NNS statistics between eigenvalues. The NNS statistics encode non-trivial correlations between eigenvalues that give rise to quantum chaotic or integrable dynamics. In the following subsections, we will introduce a simple model for a random spectrum which allows us to specify the NNS statistics of the system we are trying to model. Such an approach is unorthodox from the perspective of random matrix theory. Usually, one specifies an ensemble of random matrices and then derives the spectral statistics of the ensemble. Here, we take a different approach and generate a random spectrum by specifying the NNS statistics. The choice of NNS statistics is then assumed to correspond to some choice of random matrix ensemble (for example, by using the Wigner surmise for the NNS distribution, our model should approximately describe the statistical behaviour of eigenvalues pulled from a Gaussian ensemble). In Section 3 we will address the issue of how well these models reproduce various aspects of the spectral form factor of certain matrix ensembles.

Random Spectrum with Fixed NNS ("i.i.d Model")
We begin the construction of our random spectrum model by expressing its energy levels as follows: where E gs is the ground state energy, which we assume is fixed to a constant E gs = E 0 , and δE k are independent-identically-distributed (i.i.d) random variables which follow a probability distribution P. By assuming that P is a distribution with non-zero support for δE k ∈ R + we get an ordered set of energy levels in the spectrum Throughout the rest of this paper we will refer to this model as the i.i.d model.
Since we choose δE k to be i.i.d random variables we can express the joint probability density of the energy levels {E k } N k=1 of the spectrum as 7 : Using this, we can define the average of some function, f (E 1 , .., E N ), associated with the random spectrum: For example, the spacing distribution between two eigenvalues E i and E j in terms of the joint probability distribution function is given by the following integral expression: Using the formula for the spectrum defined by Eq. (2.1), we can see that the spacing distribution between the adjacent energy levels E m and E m+1 is 8 : This proves that the NNS distribution between any nearest neighbor pair is given by P in the i.i.d model. 7 Note that under the change of variables δE k = E k − E k−1 we have the equivalence of probability measures dE1 · · · dEN = dδE1 · · · dδEN . This is due to the Jacobian being a lower(or upper) triangular matrix filled with 1's. 8 The result that the NNS distribution is equal to P relies on the assumption that P(x) = Θ(x)P(x) i.e. the probability distribution can only have non-zero support for non-negative x values, this allows us to know the ordering of the variables E0 ≤ E1 ≤ · · · ≤ EN .
We can also define the average spectral density by doing the following integral: (2.6) The averaged double spectral density is also another useful quantity defined by: (2.7) Using this we can also define the connected double spectral density given by: (2.8)

Averaged Spectral Form Factor in the i.i.d Model
In this subsection, we will derive a closed form expression for the annealed averaged normalized spectral form factor associated with the type of spectrum we discussed in subsection 2.1. The normalized annealed averaged spectral form factor is given by the following expression: (2.9) Following the derivation provided in Appendix A, we will find the following expression for the averaged spectral form factor 9 : (2.10) We can simplify the complicated expression under the assumption that |a| < 1, |b| < 1, and N 1 (thermodynamic regime). In this case, at leading order, we ignore terms of the order O(a N ) and O(b N ) to obtain the following approximation for the normalized form factor: (2.11) Note that the approximation given in Eq. (2.11) is only expected to be useful at finite temperatures for sufficiently large values of N . It is expected to fail when β = 0. When β = 0 we formally take the limit of Eq. (2.10) above as a → 1 to get the following normalized average form factor at β = 0: (2.12) 9 Note, the expression for the form factor does not rely on the assumption of P begin zero for negative arguments. So in principle one could actually have P with support at negative values (you will just extend the lower limit of integration to −∞ for the expressions describing the form factor) but then the interpretation of P being the NNS distribution is not as transparent. Although, we suspect that as long as P is an even function then it can still roughly be interpreted as the NNS distribution between eigenvalues in this model. In cases where b and b * go to zero as t → ∞ it is easy to see that 10 : At infinite temperature we identify the entropy as S = ln(N + 1) and conclude that after a sufficiently long time, the infinite temperature averaged form factor will plateau toward a value e −S . In a similar manner, we can approximate the expression for the large time value of the form factor (again we assume b and b * vanish in the t → ∞ limit) at finite temperature for large values of N using the approximation given in Eq. (2.11): (2.14)

Spectral Form Factor of Common NNS Distributions
Thus far, we have proposed the i.i.d model for a random spectrum which involves specifying the NNS statistics of the system we are interested in modelling. In the following sections we will use the i.i.d model to compute the averaged spectral form factor associated with various choices of NNS statistics. In particular, we will investigate how well the i.i.d model (with a Wigner surmise NNS distribution) reproduces the features of the spectral form factor for Gaussian ensembles in Section 3.3. Aside from the more common NNS statistics that are found in integrable and quantum chaotic systems, we will also study generalizations of these statistics in Section 3.4 and explore how the behaviour of the form factor changes.

Delta Function (Non-Random Evenly Spaced Spectrum)
The most trivial spectrum that one can generate using the i.i.d model is an evenly spaced spectrum (with no degeneracy) where the NNS distribution is given by the Dirac delta distribution below: where ∆E ≥ 0, in this case the integrals in Eq. (2.10) are simply given by 11 : The normalized averaged spectral form factor equals: 10 In many of the examples we consider the assumption that b and b * vanish in the large t limit will be valid. However, there are examples where the assumption fails. In Sec. 3.1 we have P(x) equal to a delta function, which results in an oscillating form factor that does not settle any particular value in the large t limit. 11 It is clear from this that b and b * do not go to zero in the large t limit this is an example in which the form factor does not plateau towards a fixed value.
In this case we can see from the expression above that the form factor is periodic with a period:

Poisson Distribution (Quantum Integrable Systems)
For sufficiently large diagonal matrices whose diagonal elements are i.i.d random variables, it is known that the NNS distribution of eigenvalues follows a Poisson distribution [45,52]. Furthermore, the Poisson distribution also appears when studying the energy spacing statistics of a wide variety of integrable systems [53][54][55]. This gives us reason to consider the spectral form factor of a spectrum generated by the i.i.d model whose NNS distribution is given by the Poisson distribution 12 : The average spacing between adjacent eigenvalues is given by δE = σ. We compute the integral expressions for a, b, b * , they are given by: .
Notice that, in this case b and b * do go to zero in the limit as t → ∞, this means that the normalized average spectral form factor at infinite temperature goes to (N + 1) −1 . At infinite temperature the form factor is given as: For finite temperature, the expression for the form factor is not as simple -so we will not explicitly write it here. However, because |b| < 1 and |a| < 1 we can use the approximated form factor in the large N regime given by Eq. (2.11) to give a simple approximation for the form factor: We compare the approximated expression for the form factor given in Eq. We can see that for sufficiently large N the approximated expression converges toward the exact expression for the form factor. At lower temperatures, convergence will occur at smaller values of N . We also consider the averaged spectral form factor at various temperatures for a fixed value of N = 100 in Figure 3 13 . The expression for the asymptotic value of the form factor is given by the following expression: We can see that the form factor decays from its initial normalized value towards a non-zero value given by Eq. (3.9). Another interesting point to make based on the plots we made is that the plateau phase occurs at a time scale t σ −1 = δE −1 , where δE is the average spacing between eigenvalues. At infinite temperature, we see oscillations, these are similar to the oscillations in the form factor discussed in [50] they arise at high temperatures due to the spectral edges of the spectral density 14 .
It is possible to explicitly calculate the averaged spectral density in the i.i.d model with Poisson NNS statistics (details are given in Appendix B) it is given by the following expression: We can rewrite the sum in terms of the incomplete gamma function as follows: (3.11) By taking the derivatives of incomplete gamma function we can show that the function is monotonically decreasing and has an inflection point at E − E 0 = (N − 1)σ. We expect the function to very slowly decrease up until one gets close to the inflection point. Then after the inflection point we expect the spectral density to be very small. In Figure 4 we plot the spectral density and verify this. Ignoring the delta function at E = 0, we can see the spectral density is approximately constant then quickly goes to zero near the inflection point. We can also see that as N increases the nearly constant phase is extended. The nearly constant phase can be regarded as a natural consequence of the i.i.d spacing assumption in our model. We will generally expect that spectral densities generated by the i.i.d model will have an approximately constant portion far from the edges of the spectrum. We can also see that width of the density scales with N .

Wigner Surmise (Quantum Chaotic Systems)
In this section we will compute the spectral form factor in the i.i.d model when the NNS distribution is given by the Wigner surmise. The Wigner surmise is known to be a good approximation to the NNS distribution for classical Gaussian ensembles [45,54,56]. The Wigner surmise NNS distribution is defined below 15 : The value of c (sometimes called the Dyson index) depends on the type of Gaussian ensemble one wants to consider. In particular, the three main classical Gaussian ensembles are 15 The value of σ fixes the average spacing between eigenvalues and is usually chosen so that the spacing is unity. Just as in the Poisson spacing case we will not make assumptions on the value of σ and simply keep σ in our expressions.
orthogonal (c = 1), unitary (c = 2), and symplectic (c = 4). In Figure 5 we plot the distributions for the three ensembles for σ = 1. For the sake of simplicity we will mainly focus on the Gaussian unitary ensemble, c = 2 16 . In this case we will have the following NNS distribution: The average spacing between adjacent eigenvalues is given by δE = 2σ √ π ∼ 1.1σ. The expressions for a, b, and b * in this case are given by: (3.14) One can easily check that for β > 0 we have |a| < 1 and |b| < 1. This means that the thermodynamic approximation given by Eq. (2.11) will be valid for sufficiently large N at finite temperatures. Also, we can check that lim t→∞ b = 0, which tells us that the asymptotic value of the form factor at infinite temperature is (N + 1) −1 . We can plot the form factor to get a sense of the general features that appear for various choices of parameters.
In Figure 6, we compare the thermodynamic approximation given by Eq. (2.11) with the exact result given by Eq. (2.10). As expected, for sufficiently large N the thermodynamic approximation converges to the exact result, and one can further verify that for lower temperatures the convergence will occur more quickly at lower values of N .
In Figure 7, we plot the exact form factor for fixed N at different temperatures. We can see that the form factor initially decreases and eventually reaches a minimum value 17 . After reaching the minimum, in a time scale roughly given by t ∼ δE −1 , the form factor increases along a "ramp" and finally plateaus toward an asymptotic value that is given by taking the t → ∞ limit.
The ramp and plateau also appears in numerical calculations of the averaged form factor of the SYK model discussed in [50]. This occurs since the spectrum spacing statistics of the SYK model can be understood in terms of the spacing statistics of random matrices pulled from Gaussian ensembles. Although it is true that the form factors of the SYK and the classical Gaussian ensembles have many features in common, it is important to note that there are differences. One important distinction to make is the time scale over which the "ramp" begins to manifest in the form factor. In the case of the form factor generated by random matrices pulled from the GUE (and also our i.i.d. model see Appendix C) the beginning of the ramp occurs on time scales comparable to the Heisenberg time, which is 17 At high temperatures we see oscillations similar to what we saw for Poisson spacing, this is again due to the spectral edges being probed at high temperatures (not to be interpreted as "echoes"). O( δE −1 ). This is in stark contrast to the SYK model where the ramp phase begins on time scales which are much shorter than the Heisenberg time 18 .
In Appendix C, we do a numerical study which compares the form factor generated by eigenvalues from a Hermitian random matrix to the form factor of a random spectrum generated by the i.i.d model with Wigner surmise NNS distribution. We show that if one focuses on the eigenvalues near the centre of the spectrum then our model is in reasonable agreement with actual numerical computations of the form factor associated with the truncated spectrum near the centre of the spectral density.
It is not possible to get a closed form expression for the spectral density when we have Wigner surmise spacings. However, we can numerically compute the averaged spectral density by generating a large number of eigenvalues from the i.i.d model for Wigner surmise spacing. As an example, in Figure 8, we generate a histogram from the eigenvalues generated by the i.i.d model for 10 4 samples with each sample having 100 eigenvalues (for simplicity we set σ = 1 and E 0 = 0). Unsurprisingly, we find that the density away from the edges is approximately constant and similar to the Poisson case. However, it will differ from the Poisson density near the ground state where it will be zero due to the repulsion between the ground state and the rest of the excited states (Poisson case exhibits "attraction toward the ground state"). Another important point to make is that the averaged spectral density 18 [57] estimates the ramp in the SYK form factor to start at times scales of the order O(N 1/2 ln(N )). This is much shorter than the Heisenberg time for the SYK model which is of the order O(e N ) [58]. For more physical quantum chaotic theories, the time scale after which the universal ramp and plateau manifest , the so-called "Thouless time", is not universal (i.e. it changes depending on the model and the observable. We thank Julian Sonner for pointing this out). of the spectrum generated by the i.i.d model is not the semi-circle as it is for classical Guassian ensembles (e.g., see Figure 17). This is not surprising since our model only captures correlations between nearest neighbors and not the longer range correlations that conspire to give the semi-circle. Nonetheless, the i.i.d model still captures the general features of the ramp and plateau in the averaged form factor.
In the next subsection, we will introduce a new NNS distribution which retains the interesting features of the Wigner surmise (i.e. repulsion in NNS statistics) but will be easier to analytically handle in our model when computing the averaged spectral density in the i.i.d model. We will also discuss generalizations of classical Gaussian ensembles and how the late time behaviour of the form factor of such ensembles differ from the usual ramp and plateau behaviour of the form factor of classical Gaussian ensembles.

Form Factor of Gamma NNS Distribution and β-Ensembles
In the previous subsection, we looked at the Wigner surmise as a canonical example of what the NNS distribution of eigenvalues of a chaotic system looks like. Although we are able to compute the form factor, many other quantities of interest such as the spectral density are difficult to compute in closed form. This is primarily due to the e −x 2 in the NNS density makes the integrals in Eq. (2.6) difficult to compute. To facilitate more explicit calculations -while still retaining the essential repulsion between eigenvalues -we will consider a slightly different NNS distribution. We define the gamma distribution for NNS: The value of q (plays the role of the Dyson index in the Wigner surmise) fixes the degree of repulsion between eigenvalues. For a fixed value of q, σ fixes the average spacing between nearest neighbor pairs to δE = (1 + q)σ. At leading order, it is clear that such a distribution will contain the same repulsion behaviour near zero spacing as the Wigner surmise with c in Eq. (3.12) identified with q in Eq. (3.15). The major difference being the tail; the Wigner surmise has a Gaussian tail whereas the gamma distribution has an exponentially decaying tail. The advantage to using this is that it still contains the important repulsion of a chaotic system for any q > 0. Furthermore, due to the exponential tail of the gamma function the spectral density can be computed exactly.
We begin by computing the necessary integrals that define the averaged form factor in the i.i.d model: We can clearly see that |a|, |b| < 1 at finite temperature for q ≥ 0, therefore in the large N regime we can use the thermodynamic expression of the form factor given in Eq. (2.11). We find: The plateau at sufficiently large N and finite temperature regime has an approximate height equal to: We plot the averaged form factor in the i.i.d model using Eq. (2.10) for various choices of parameters. The energy scale, σ, is related to the average energy spacing between nearest neighbor eigenvalues, δE , through the following relation δE = (q + 1)σ.
In Figure 9, we plot the averaged form factor at infinite temperature with N = 100 and vary q. At lower values of q > 1, we still see that ramp and plateau (like the Wigner surmise), however as we increase the value of q we start to see oscillations after the initial dip before saturation to the plateau. In Figure 10, we verify that these late time oscillations, at large q persist at lower temperatures with the period roughly being the same as the temperature varies.
The existence of these oscillations for larger values of q can be attributed to the NNS distribution localizing far from the origin and forming a substantial "gap" between the origin and the main distribution. One can check that the form factor of any system that has NNS density with a large gap will exhibit oscillations (although depending on the details of the NNS distribution the shape of the oscillations as well has how long they persist will vary). In particular, if one allows the Dyson index, c, in Eq. (3.12) to be sufficiently large one can also check that similar oscillations will arise for the Wigner surmise.
It is also useful to compare the averaged form factor with the form factor of a single sample in our model. In Figure 11, we can see that the regular late time oscillations are not self averaging. This is similar to the nature of the ramp and plateau behaviour studied in form factors of Gaussian ensembles [50]. The ramp and plateau manifest most clearly after averaging over many samples but is more difficult to ascertain when one analyzes only a single sample.
Due to the simple expression of the form factor in the large N and finite temperature regime it is possible to get a sense of where the oscillations occur by finding local extrema. By taking the first derivative of the expression in Eq. (3.17) and setting it to zero we find the following condition: In general, the solution to this equation cannot be found exactly (in Appendix D we discuss the special cases when q = 1, 2, 4 where exact solutions can be found). In the q 1 regime it can be argued (see Appendix D for the argument) that the period of the visible oscillations in the form factor will roughly be given by: Above, we used that the average spacing between adjacent energy levels in our model is δE = (q + 1)σ ≈ qσ. Now that we have discussed the general features of the form factor we consider the averaged spectral density -which is given by (see Appendix E for details of the calculation): Ignoring the delta function at the ground state, we can see the left most edge of the spectral density near E = E 0 is dominated by the first term in the sum over m and vanishes at zero. To get an estimate for how wide the spectral density is, we analyze the m = N term in the sum. This will give us a rough sense of what the tail of the spectral density looks like: This function will give a bell shaped curve. We will be interested in the right most inflection point of the bell curve to give us an idea of where the edge of the spectral density is. This involves solving d 2 dE 2 ρ tail (E) = 0 which yields a simple quadratic equation. The larger root gives us an estimate for where the right edge of the spectral density is: So, for N 1 and q ≥ 1, the width of the spectral density roughly scales as N (q + 1). In Figure 12 we can verify these findings by a simple plot of the spectral density with σ = 1 and N = 20 for various values of q. We can see that for larger values of q the low energy part of the spectrum exhibits visible oscillations. This is due to the fact that there is little overlap between the lower m terms that make up the spectral density. A straightforward calculation of the local maximum of each term reveals that the distance between the visible local maxima is roughly q + 1.
It is interesting to ask what type of random matrix ensembles give rise to NNS statistics with a more general Dyson index. One answer to this question can be found in the study of "β-ensembles" which are discussed in [51]. These ensembles have been shown to arise from N × N tri-diagonal random matrices of the form 19 : where N (0, 2) denotes random real variables pulled from a normal distribution with variance of 2 and mean of 0. The χ m denotes positive real random variables pulled from the Chidistribution with m degrees of freedom 20 . It can be proved that the joint eigenvalue density (up to a normalization constant) of these matrices is given by: The joint eigenvalue density for the classical Gaussian ensembles are reproduced when β = 1, 2, 4 -but more general values of β can also occur 21 . Since β controls the degree of repulsion between eigenvalues, it seems reasonable to identify it with the Dyson index, c, in the Wigner surmise for NNS distribution. This conjecture was explored numerically in [59], where it was found that a (generalized) gamma distribution provided reasonable approximation to the NNS distribution of β-ensembles. Based on the results we obtained using the i.i.d model in this section, we suggest that with sufficiently large values of β (here β is not inverse temperature it is the Dyson index which labels the β-ensemble) the form factor of β-ensembles will have late time oscillations in the ensemble averaged form factor. In the next subsection, we briefly explore a class of simple oscillator systems that demonstrate sharp gaps in the NNS statistics. In these systems one can also see regular decaying oscillations in the late time behaviour similar to the ones we studied here.

Chaotic Perturbation of a Harmonic Oscillator
One concrete example of a class of systems that have large gaps in the NNS distribution are interacting oscillator systems whose Hamiltonian may be written in the form: where H 0 is a Hamiltonian with a regularly spaced spectrum. H chaos is a chaotic Hamiltonian in the sense that the NNS statistics of the eigenvalues of H chaos obey a Wigner surmise type distribution. The parameter λ controls the relative strength of the two terms.
One familiar example of a system that has an evenly spaced spectrum is the one dimensional Harmonic oscillator vibrating at a frequency ω 0 . The introduction of the a chaotic term in such a scenario would be analogous to adding an additional potential term in the Hamiltonian which generates chaotic dynamics. Concretely, we write: where ω 0 is some characteristic energy scale of the unperturbed oscillator and is a characteristic energy scale for the chaotic Hamiltonian. In terms of matrices we let H 0 be an N × N diagonal matrix with entries (H 0 ) ij = N −1 (i − 1)δ ij (we choose to normalize by N so that the width of the unperturbed spectrum does not change). We model H chaos using a random matrix. Roughly speaking, the NNS distribution between adjacent eigenvalues will be delta functions when = 0 once the chaotic term is turned on we will expect the delta functions to broaden and look similar to a kind of translated Wigner surmise. This results in a large gap in the NNS distribution which would lead to oscillations in the form 21 It was shown in [51] that the matrices in the classical Gaussian ensembles can be written in the tridiagonal form. For example, a matrix pulled from the GOE ensemble can be tridiagonalized and the matrix will belong to a "β = 1-ensemble." factor. To get approximation for what the NNS distribution might look like, it is useful to compute the NNS distribution for 2 × 2 matrices. This is a straightforward exercise which is done explicitly in Appendix F, the final result is given by: (3.28) Using the i.i.d model with the NNS distribution given by Eq. (3.28) we can compute the form factor. Figure 13 plots the form factor at infinite temperature for various choices of /ω 0 . We can see that oscillations appear for smaller values of /ω 0 and begin to go away as the chaotic term becomes larger. The decay of the oscillations occur due to the random chaotic term in the Hamiltonian. In this scenario the parameter ω 0 / characterizes the sharpness and size of the gap in the NNS distribution. Using the results we found in the previous subsection for gamma NNS statistics, we expect the period of the oscillations in the form factor in the /ω 0 1 regime to be approximated by: where we used the fact that δE ≈ ω 0 /2 when /ω 0 1.

Systems with Self-Averaging Oscillations in the Form Factor (Echoes)
In the previous section, we studied features of the averaged form factor and spectral densities associated with different choices of NNS distributions using the i.i.d model. For the i.i.d.
model of random spectra we generally found that the spectral form factor decays until a time scale roughly given by t ∼ δE −1 , where δE is the average nearest neighbor spacing between eigenvalues. At time scales t δE −1 the spectral density can be approximated as smooth and gives rise to self averaging behaviour at early times. After this time scale, the discreteness of the system manifests and the form factor is no longer self averaging. The averaged late time behaviour over many samples depends on the details of the NNS distribution for the eigenvalues of the system. From the examples we studied, we had three classes of behaviour after the initial dip: 1. Continues to decrease and plateaus to a value (e.g., Poisson NNS where eigenvalues exhibit "attraction").
2. Increases for a duration of time along a "ramp" followed by a plateau (e.g., Wigner surmise NNS where eigenvalues exhibit "repulsion").
3. Increases and exhibits damped oscillations toward a plateau value (e.g., large q gamma NNS where the eigenvalues exhibit "enhanced" repulsion resulting in large/sharp gaps between adjacent eigenvalues).
Case 1 describes the NNS statistics of generic quantum integrable systems and has no "ramp" in the form factor 22 . Case 2 describes NNS statistics of quantum chaotic systems whose spectral form factors contain a "ramp." In both these cases, we can see that there are no regular oscillations. In Case 3, there are oscillations in the averaged spectral form factor. These oscillations, however, occur in the late time behaviour of the form factor when it is no longer self-averaging (e.g., Figure 11). Furthermore, the time scale over which these oscillations manifest are of the order of the Heisenberg time, δE −1 . In the context of black holes one usually regards the energy spacing between black hole microstates to be of the order e −S . This gives a Heisenberg time of the order e S . So if we want to interpret the late time oscillations as echoes then we would have to wait a time scale of the order of the Heisenberg time which is extremely long for a macroscopic black hole and for practical purposes unobservable 23 .
In this section, we shall revisit one of the key assumption of Sections 2 and 3, i.e. that the NNS distribution are independent identically distributed (i.i.d). As an example, the i.i.d assumption can be significantly violated if there are approximate degeneracies in the spectrum -causing many states to cluster around the same value while also being widely 22 Systems with constant energy spacings, such as harmonic oscillator are exceptions. Such systems have form factors that oscillate forever and never saturate to a particular value. 23 There might be a "caveat" to these findings. separated from other clusters of states. An example of this phenomenon is the hydrogen atom, where the splitting of degenerate states is controlled by the fine structure constant, α. We can see that the assumption of "independent-identically-distributed" is not valid anymore, as the energy spacing is different by a factor of α, depending on whether states belong to the same, or a different degenerate sub-sector. We shall see that this structure can also lead to self-averaging oscillations in the spectral form factor. What one needs is to have an averaged spectral density that contains regularly spaced "spikes" (or neardegenerate sub-sectors). The introduction of random interactions (regardless of how small) will generally lead to chaotic level splittings (a broadening of the "spikes") within each degenerate sub-sector. We will discuss one simple but explicit example in the following subsection which illustrates this.

Randomly Coupled Identical Fermionic Oscillators
In this subsection, we will analyze a system of N/2 Fermionic oscillators (N is even) which are coupled by a random matrix pulled from the GUE 24 . The Hamiltonian describing N/2 uncoupled Fermionic oscillators, which we denote as H F HO , can be written as: where b k and b † k are operators which satisfy the anti-commutation relations, and ω k is the frequency of oscillation of the k-th oscillator (set = 1). Before analyzing the effect of the GUE coupling we review the eigenstates of H F HO .
We denote the ground state of the FHO Hamiltonian as |0 , it satisfies the following identity: b k |0 = 0. The operators b k and b † k play the role of annihilation and creation operators respectively. We can check that the following identities are true: We can create any state in the Fock space by starting with the ground state and acting with the creation and annihilation operators. There are exactly 2 N/2 distinct states in the Hilbert space since we can only act once on the ground state using b † k for a particular k. These states (by convention) can be written as:

4)
24 Although we choose to work with fermionic oscillators, similar techniques can also be used to model a system of bosonic oscillators. We primarily choose fermionic degrees of freedom since they can be discussed in the context of the SYK model (a toy model describing the near horizon dynamics of extremal black holes) where the eigenstates of fermionic oscillator are the pure states studied in [60]. Furthermore, eigenstates can also be mapped to binary strings making this an interesting "qubit" model.
where n 1 , .., n N/2 = 0 or 1 (we define b † k 0 = I, the identity operator). When we act with a creation operator we find that: = (−1) σ k δ 0,n k |n 1 , n 2 , .., n k + 1, .., n N/2 , Similarly acting with an annihilation operator gives:  We define the k-th number operator N k as: This means the Fock states are eigenstates of the number operator and since the Hamiltonian can also be written in terms of the number operator the Fock states diagonalize the Hamiltonian. In particular, we have: (4.8) To simplify our considerations we will set ω k = ω 0 for all k -sometimes referred to as a mass term. In this case, we will have a regularly spaced spectrum which is easy to analyze. We define p = 0, 1, 2, .., N/2 as the number of "occupied" sites for a particular state (i.e. the number of 1s in the ket). Then the energy of such a state is given by: The number of microstates with the energy E(p) is given by: (4.10) Using these results we can find the canonical partition function for the system given by: p=0 Ω(p)e −βE(p) = e The spectral form factor for this system will oscillate with a period given by: Now that we have discussed the details of the FHO, we address what happens when we add a random chaotic coupling. In particular, we will consider: We know that the unperturbed energies are given by Eq. (4.9) and we know that these energies are degenerate. We will label each degenerate sector by p, which is equal to the number of occupied sites. We define the degenerate eigenstates in the p-th degenerate sector as {| n p } Ω(p) p=1 . Within the p-th degenerate sector we would go about computing the matrix elements: m p |H chaos | n p . (4.14) Diagonalizing a matrix with the matrix elements given above would tell us the information about how the degenerate energies within the p-th sector split. At leading order in degenerate perturbation theory, the splitting of the energy levels only depend on the details of the sub-block in the degenerate sector. This means that first order theory will contain no information about correlations between different degenerate energy sectors. In other words, at leading order we can treat our random matrix as block diagonal with each block having the same size as the particular degenerate sector. More specifically, we can write: where I(p) is an identity matrix of dimension Ω(p) × Ω(p) and H chaos (p) is a sub-block of the chaotic interaction Hamiltonian of dimension Ω(p) × Ω(p).
To simplify our considerations, we will assume that the H chaos is given by a random matrix pulled from the GUE 25 . If H chaos is a random matrix from the GUE, then we know the splitting of degenerate states within a sector will occur so that the NNS statistics are consistent with the GUE (i.e. you will have Wigner surmise spacing). Furthermore, we expect the spectral density to spread out from a delta function at the degenerate energies to "semi-circles" centered around the degenerate energies. The width of the spread for the p-th degenerate sector will roughly go as: (4.16) This allows us to define a rough criterion in which we expect perturbation theory to hold: The elements of the GUE are Hermitian matrices. The off-diagonal elements will be of the form x + iy where x and y are random variables pulled from a Gaussian distribution with a mean of zero and variance of 1/2. The diagonal elements are real random variables pulled from a Gaussian with mean of zero and variance of 1. In this case it is well known that the width of the semi-circle spectral density goes as 4 √ N [45]. This width of the semi-circle can be scaled by multiplying the matrices in the ensemble by a factor c then the new width is 4c √ N . This essentially states we want to work in a regime where the spread has a width much smaller than the spacing between energy levels of the unperturbed oscillator. Figure 14 depicts numerically generated spectral density at N = 22 for different choices of the parameter /ω 0 . One can check that in the specific case we plotted the condition given in Eq. (4.17) for being in the perturbative regime translates to /ω 0 10 −2 . We can see that the average spectral densities on the top row of Figure 14 satisfy the condition and the spread within the degenerate sectors is much smaller than the spacing between different degenerate sectors. In contrast, the plots on the bottom row do not satisfy the perturbative condition and we can see that at /ω 0 = 10 −2 , the spread within the degenerate sectors start to become comparable to the spacing between them. The spread will increase as we increase /ω 0 and there will be substantial overlap between sectors, eventually leading to a point where the overall density simply looks like that of the GUE at /ω 0 = 10 −1 . In Appendix G, we do a numerical analysis of the exact width of the spread within the degenerate sectors and compare to what we expect from first order degenerate perturbation theory, finding reasonable agreement for /ω 0 ≤ 10 −2 . Now that we have discussed the spectral density, we can turn our attention to the spectral form factor. Within the perturbative regime, there are two distinct time scales.
One is given by the oscillator frequency ω −1 0 and the other scale is given by −1 which gives the strength of the chaotic perturbation (roughly related to the average spacing of eigenvalues within a degenerate sector). At time scales t −1 , the form factor will be dominated by the coarse grained spectrum which is a sharply spiked spectrum with regular spacing. Due to this we expect to see oscillations in the form factor with a period that is well described by the unperturbed form factor, τ = 2πω −1 0 . We will also expect these oscillations to gradually decay due to the repulsion of adjacent eigenvalues within a degenerate sector. For time scales t −1 , the form factor will become sensitive to fine structure of the degenerate sectors. In this regime, the form factor should contain a ramp and a plateau consistent with the late time behaviour of a form factor in the GUE, as we have seen before. In Figures 15 and 16, we numerically plot the form factor (averaged over 100 samples) for different values of /ω 0 at βω 0 = 0. For the bottom right plot we compare the numeric early time behaviour to the behaviour predicted by Eq. (4.19), which is the form factor of a semi-circle spectral density.
The the perturbative regimes (weakly coupled) occurs when /ω 0 = 10 −4 , 10 −3 and the non-perturbative regimes (strongly coupled) occur at /ω 0 = 10 −2 , 10 −1 . Figure 15 depicts the early time behaviour of the form factor and it shows that the form factor oscillates at early times with a frequency 2πω −1 0 (the red dots indicate the value of the form factor at regular time steps given by the period τ = 2πω −1 0 ). The oscillations have a decaying amplitude bounded by an envelope which can be numerically fit by: where n 0 and n are fit parameters. We can see from the legend in Figure 15 that the amplitude of the oscillations in the weakly coupled regime decays through a power law with a power n ∼ 3, furthermore we can see that n 0 ∼ Cω 0 / where C is a proportionality constant that will depend on N 26 . In the non-perturbative regime (the bottom two plots of Figure 15), we can see that the regular oscillations start to disappear and give way to the more familiar features present in the early time behaviour of the form factor for the GUE. The early time behaviour of the form factor in the case when /ω 0 = 10 −1 is well approximated by the form factor associated with a semi-circle density (SCD) given by: 19) where E 0 is the half the width of the semi-circle density (in our case we used E 0 /ω 0 = √ 2 13 /10 ≈ 9.05).
In Figure 16, we depict the averaged form factor on a Log-Log scale at longer time scales (averaged over 100 samples). The early time behaviour contains oscillations in the weakly coupled regime which we already described in our discussion of Figure 15. The averaged late time behaviour in all coupling regimes contains the ramp and plateau. This late time behaviour occurs due to the manner in which each degenerate sector of the uncoupled oscillator theory splits in the presence of the GUE matrix coupling.

Randomly Coupled Oscillators as a Toy Model for the Membrane
Let us now see how our system of chaotically coupled oscillators, can act as a quantum mechanical toy model to describe gravitational waves echoes (e.g., [21]).
To do this, we recall that gravitational wave echoes are thought to arise due to quantum gravity effects near the horizon which result in the partial reflection of perturbations. These reflections are usually modelled by cutting off the semi-classical black hole geometry close to the horizon and then placing semi-reflective boundary conditions at the cutoff. In the context of large AdS black holes, one can also place a cutoff near the horizon and echoes will arise due to the reflection of perturbations near the horizon and also at the conformal boundary [23]. In this paper, we adopt the view that the cutoff with modified boundary conditions describe some kind of dynamical membrane deep in the bulk near the horizon. This membrane has its own internal spectrum which should resemble the chaotic and discrete spectrum of a black hole.
If one begins with a perfectly reflective membrane, then perturbations to the geometry will not decay and the geometry will have normal modes (rather than quasi-normal modes). The dynamics of the perturbations on the spacetime may be roughly modelled by some kind of oscillator system similar to the one discussed in the previous section. The characteristic frequency of the oscillators in the previous subsection would be identified with the spacing between normal modes in the cutoff spacetime. For example, for a spherical Schwarzschild AdS black hole (with r H /L 1), if the cutoff is a proper radial length, l prop , from the horizon, then we expect massless scalar perturbations to make a round trip (from the cutoff to the conformal boundary and back) in a time scale roughly given by [23]: (4.20) where r H and L is the horizon and AdS radii, respectively. We identify ω 0 = 2πτ −1 with the characteristic oscillator frequencies 27 . This will give rise to a periodic form factor with oscillations that have a period equal to τ . To introduce dissipation into the system, we can further add chaotic interactions between the oscillators. In the weak coupling regime, these interactions are responsible for the decay of the regular oscillations at early times in the form factor and also will be responsible for the ramp and plateau at late times making late time thermalization similar to SYK theory.
To interpret the energy scale in Eq. (4.13), we turn off the oscillator term by sending the cutoff to the horizon (i.e. l prop = 0 in Eq. (4.20)) so that ω 0 = 0, and we are left with the chaotic term. This gives rise to the usual semi-circle spectral density and gives a standard SYK/JT gravity type model for a black hole -where no oscillations are present at early times. In this case, the energy scale, , controls the average spacing between adjacent energy levels which is of the order δE ∼ e −S/2 where S is the entropy of the black hole 28 . 27 Notice that the frequencies have a linear temperature dependence, this could arise as a thermally induced mass from a renormalization procedure in QFT at finite temperature [61,62]. 28 We obtain e −S/2 by recalling that the width of the semi-circle scales with the size of the matrix √ 2 N/2 . Then we take the width of the semi-circle and divide by the total number of states which is 2 N/2 , to get an average spacing δE ∼ 2 −N/4 . Finally we identify S ∼ ln 2 N/2 to obtain δE ∼ e −S/2 . At the same time, if we consider the first law of black hole thermodynamics it suggests that the average spacing between adjacent levels is of the order δE min = δS min /β. We will be agnostic about what δS min should be. We know the smallest possible value of δS min is of the order e −S but there might also be other choices (for example, δS min ∼ 1 might also be sensible when discussing entropy changes by emission of a single Hawking quanta).
Identifying the δE expression from the random matrix model to the δE min from the first law gives ∼ β −1 δS min e S/2 . When we do this and consider the dimensionless ratio between and ω 0 we have: If δS min ∼ e −S (the smallest possible change in entropy) and l prop ∼ p , where p is the Planck length, then our oscillator system is in the weakly coupled regime and will exhibit echoes that manifest after a scrambling time scale. Such echoes would be observable within a reasonable window of time for astrophysical black holes, making this model interesting in the discussion of gravitational wave echoes.

Discussion and Conclusion
The AdS/CFT correspondence formulates quantum gravitational systems in AdS in terms of CFT systems on the boundary. In such a scenario, the statement that a black hole is a quantum chaotic system is understood in terms of the high energy eignestates of the CFT.
In the context of spacing statistics of the CFT spectrum, one could imagine diagonalizing the Hamiltonian H to find a complete orthogonal basis of eigenstates |n and eigenvalues E n . Any state can then be written in terms of the eigenstates of the Hamiltonian including the black hole. As one goes to "very high" energies, one should expect the spectrum of the CFT to exhibit chaotic spacing statistics (i.e. eigenvalue repulsion). If one assumes that the a black hole is a thermal ensemble of such microstates, then how the black hole relaxes after being perturbed should be described by thermal correlation functions. The question of how quantum aspects of a black hole manifest in the ringdown and whether they are potentially observable could be addressed using this picture.
In this paper, we explored the concept of echoes due the black hole microstructure from the perspective of the thermalization behaviour of quantum chaotic systems. The form factor served as a convenient proxy to study this behaviour for systems with varying spectral statistics. We employed a simple model for a random spectrum that involved summing a set of independent-identically-distributed (i.i.d) random variables which represent the spacing between adjacent eigenvalues in the spectrum. By defining a random spectrum in this manner, we were able to specify the resulting spectrum's Nearest Neighbor Spacing (NNS) statistics. Despite the simplicity of the i.i.d model, we demonstrated that it could reproduce various important late time features seen in the spectral form factor of strongly coupled chaotic systems. Most notably, by using Wigner surmise NNS statistics we were able to obtain the ramp and plateau behaviour. We also used the i.i.d model to study more general NNS distributions where the "Dyson index" is allowed to be any real positive number, rather than the canonical values of 1, 2, and 4 (corresponding to the GOE, GUE, and GSE ensembles respectively).
We briefly discussed how more general choices of the Dyson index could correspond to more general random matrix ensembles called "β-ensembles," which are constructed from certain types of tri-diagonal matrices. We also argued that for a large Dyson index, the averaged form factor of a β-ensemble would have damped oscillations at late times before approaching a plateau. The origin of these oscillations are the large regularly spaced "gaps" in NNS distribution which are expected to arise when the Dyson index becomes large. We further showed that such gaps could be realized in systems with an evenly spaced spectrum -like a harmonic oscillator with an additional chaotic potential. Although such oscillations might be thought of as "echoes" in the black hole context, such echoes would manifest on extremely long time scales of the order e S , disqualifying them as being observable in the context of astrophysical black hole microstructure. Although, the late time oscillations found in Section 3.4 may not be interesting in the astrophysical context, they may be connected to the discussion of form factors in more exotic theories of gravity. In particular, in the work [63], the form factor of Narain's family of free boson CFTs (which is dual to "U(1) gravity" in AdS 3 ) has been shown to have a form factor that contains oscillations at late time similar to the ones we found in Section 3.4. Our results suggest that such oscillations occur due to enhanced "repulsion" between eigenvalues in the spectrum. It would be interesting to investigate if the late time oscillations showing up in the form factor of Narain ensemble CFTs can also be traced back to enhanced repulsion statistics in the spectrum.
Studies using the i.i.d. model in Sections 2 and 3 suggested that the quantum aspects of the thermalization behaviour of black holes usually manifest at very late times which are unobservable in the astrophysical context. A common aspect of the spectra generated by the i.i.d. model was that there was a single large cluster of microstates with very little coarse-grained structure. To get non-trivial behaviour on observable time scales, we thus proposed that there must be additional structure in the coarse grained density of states. In particular, for echoes to manifest, we showed that the coarse-grained density had to have large clusters of states widely separated from other clusters. In Section 4, we considered a many body system of identical fermionic oscillators coupled by an interaction term modelled by a random matrix in the GUE. In the weakly coupled regime, we used degenerate perturbation theory to argue that degenerate sectors of H F HO would split in a way consistent to the random matrix statistics of the GUE. This naturally led to regularly spaced decaying oscillations at early times followed by a ramp and plateau in the averaged form factor. We then discussed the possibility of using our model of coupled oscillators as a toy model that describe near horizon modifications that give rise to gravitational wave echoes.
Although we were able to show that one can construct quantum chaotic systems with echoes, there are still a number of issues that need to be addressed before we can view them as reasonable models for describing black holes with microstructure near the horizon. Aside from the condition that a black hole has a quantum chaotic energy spectrum, it must also satisfy additional constraints which are usually discussed in the context of out-of-time-order correlators (OTOCs) [42,64,65]. In particular, quantum chaotic systems describing black holes are expected to also saturate the Lyapunov bound discussed in [42]. We did not demonstrate that the coupled oscillator systems introduced in Section 4 satisfy this constraint. A good place to start exploring the issue of OTOCs (and also other types of correlators) is with the following concrete oscillator model 29 : where χ's are Majorana fermions: {χ α , χ β } = δ αβ , and j αβµν is a completely anti-symmetric coupling tensor where each component is an independent random Gaussian variable with a mean of zero. The quadratic term in Eq. (5.1) describes N/2 free identical fermionic oscillators expressed in the Majorana basis 30 and the quartic interaction is the standard SYK model Hamiltonian. As before, the quantity /ω 0 , would control the strength of the interactions between the oscillators. One would then have a family of theories labeled by /ω 0 . In the limit when /ω 0 → ∞ we should recover the SYK model with no echoes which saturates the Lyapunov bound (i.e. the Lyapunov exponent λ L = 2π β ). In the opposite limit, when /ω 0 → 0, we would have free fermionic oscillators (non-decaying periodic behaviour in form factor) and the Lyapunov exponent will be zero. In the intermediate regime, we should have regular decaying oscillations; similar to the ones discussed in Section 4.1. We expect the Lyapunov exponent in intermediate regimes to interpolate between the values of 0 and 2π/β (perhaps even in a discontinuous manner). The central question is what the Lyapunov exponent is when echoes manifest in the form factor. If we are able to show that in the regime where echoes exist the Lyapunov bound is nearly saturated 31 then it might be plausible that the quantum chaotic system described by a Hamiltonian given in Eq. (5.1) describes a black hole with microstructure near the horizon which exhibits echoes when perturbed. If we find that echoes only exist in coupled oscillator models when the Lyapunov exponent is nearly zero then the modified black hole interpretation is not plausible, we will likely have to consider more complicated models where we can tune the value of the Lyapunov exponent. We leave this exploration to future work.
An important point to make in the discussion of modifications to the near horizon description of black holes and their imprints on the ringdown is that they may not manifest as dramatically as described by the models discussed in Section 4. In particular, the model of identical oscillators coupled by a random interaction can give sharp echoes like the ones shown in Figure 15 at weak coupling. These models have evenly spaced clusters of 29 It is interesting to point out that a coupled oscillator system is similar to the kinds of models studied in [40], although the quadratic term is random in those models and of a more general form. 30 The transformation between the creation and annihilation operators defined in Eq. (4.1) and the The reason we say the Lyapunov bound is nearly saturated is due to results of a work that computes the Lyapunov exponent in a fuzzball geometry by analyzing geodesic motion in the vicinity of the fuzzball [66]. The work suggests that fuzzballs tend to have a slightly smaller Lyapunov exponent when compared to their black hole counter parts. This may also be true more generally for highly compact objects that resemble black holes far away but deviate near the horizon.
states. If the clusters were not evenly spaced one might get less dramatic echoes or no discernible echoes at all. Presumably, a completely random placement of these clusters would completely "wash out" a clear echo signal in the form factor. It would be interesting to explore where on the spectrum between evenly spaced clusters and randomly spaced clusters echoes persist. This might give a sense of how "resilient" the "echo" phenomenon may be to deviations from the ideal evenly-spaced clusters scenario.
Although we expect the high energy sector of a holographic CFT to exhibit chaos, the exact structure of the spectrum is unknown. The results in this paper suggest a wide variety of deviations can occur at early times before the universal ergodic behaviour manifests (i.e. ramp and plateau). The origin of possible deviations at intermediate and early times may be additional coarse-grained structure in the spectrum of states that give rise to clustering effects on energy separations β −1 e −S . In our randomly coupled oscillator models, degenerate states were responsible for this clustering. This is interesting because degeneracies arise due to "symmetries" in a system. The symmetries thus act as "seeds" for additional structure at intermediate energy spacings. If there are symmetries in a holographic CFT that have the same effect, then this may also give rise to additional structure in the spectrum which would imprint itself in the thermalization behaviour of a black hole and potentially be detected as echoes (or other deviations from GR predictions) that could be searched for in future gravitational wave experiments. The i.i.d. model might serve as a useful phenomenological tool to model the wave-forms for unitary black hole ringdowns.
To conclude our discussion, it is worth clarifying a few points. In this paper, we explored the possibility of echoes from black hole microstructure in the context of unitary quantum chaotic systems. Our primary conclusion is that echoes may appear, either due to enhanced repulsion between individual microstates, or the occurrence of regular spacing between clusters of microstates. Appearance of echoes in the unitary quantum description on time scales comparable to the scrambling time suggests the existence of Planck-scale microstructure near the horizon. The reason we asserted this was because the echo time scale found in the classical models with a cutoff depends on how close the cutoff is placed from the horizon. In particular, cutoffs placed a proper-radial Planck length from the horizon (i.e. r cutoff ∼ r h + 2 p /β) give echoes around a scrambling time scale β ln(S), 32 [23]. The randomly coupled oscillator example in Section 4 served as a simple toy model to illustrate how self-averaging echoes can arise due to non-trivial clustering of microstates. The oscillations found in Section 3.4, however, are of a different type. They are not self-averaging and arose from enhanced repulsion statistics between individual adjacent microstates. Furthermore, in the context of the i.i.d. model, these non-self averaging oscillations occur on time scales of the order of the Heisenberg time O(e S ), which is typically much longer than the scrambling time 33 . The physical interpretation of non-self averaging echoes that occur on the Heisenberg time scale is more subtle and might require a more careful understanding in terms various saddles (or even off-shell effects) that occur in the Euclidean path integral. While the time scale of e S is what one may expect from quantum tunneling from structure 32 Note that the S ∼ L/ p ∼ N for large AdS black holes. 33 Assuming the average spacing between adjacent microstates was of the order e −S .
behind the horizon, this is more likely an upper limit, given that black hole evaporation happens on polynomial times. We hope to explore these ideas in more detail in future work.

A Deriving the Spectral Form Factor of the i.i.d Model
In this section we derive Eq. (2.10). We begin with the expression for the partition function for a spectrum given by Eq. (2.1). It is given by (in this derivation we will treat β as a general complex parameter): Now we define the following recursive indexed quantity: where m = 1, 2, ..., N − 1, N . Using these indexed quantities, we may express the partition function and its complex conjugate as follows: where β * is the complex conjugate of β. We can then write: To compute Z 1 (β) , we will find a general formula for Z m (β) . We use the recursion relation given by Eq. (A.2) to find the following result for Z N (β) : (A.5) For Z N −1 (β) , we have: Continuing inductively we can show that: Using these results we can conclude that: where b * is the complex conjugate of b. We can immediately conclude that: In a similar manner, we consider the quantity B m = Z m (β)Z m (β * ) : (A.10) Using the fact that Z N +1 (β) = 0, we obtain the following result for B N : Using the fact that Z m+1 is only a function of {δE i } m+1 i=N , we have the following recursion relation: (A.12) Using the inductive relation above it is straightforward to show that: One can then explicitly evaluate the geometric sums involved and arrive at the following expression for B m : (A.14) We can check that this expression satisfies the recursion relation given in Eq. (A.12). By setting m = 1 we can obtain B 1 : (A.15) Now that we have calculated all the necessary averages we can plug them into the Eq. (A.4) to find the following result: Making the replacements β → β + it and β * → β − it gives us the expression for spectral form factor given in Eq. (2.10).

B Average Spectral Density For Poisson Spacing
In this section we go over the details of the computing the integrals involved in computing expression for the average spectral density for a Poisson spacing distribution. To begin we can write the JPDF for the energy levels given as: Using this we can write the average spectral density as: With some work we can derive the following identities: Using these integral identities we can show that: We therefore conclude that the average spectral density is given by: This gives the averaged spectral density given in Eq. (3.10). By integrating ρ(E) over E we correctly get N + 1 which is the total number of states.

C Form Factor of GUE vs i.i.d Model with Wigner Spacing
In this section we will do a brief numerical study of how well the Wigner Surmise distribution given by Eq. (3.13) fits the NNS distribution of eigenvalues of a 100×100 random Hermitian matrix (i.e. matrices in GUE). We will numerically compute the averaged form factor of the Gaussian unitary ensemble and compare to the form factor expression given by Eq.
(2.10) with a, b, b * given by Eq. (3.14). By doing this, we will get a sense of how well our naive model can capture certain aspects of the true form factor associated with matrices in GUE.
The numerical calculation is done by defining 10 5 , 100×100 random Hermitian matrices. The diagonal entries M ii are real and are pulled from the following Gaussian distribution: Histogram of Eigenvalues for 100 X 100 Matrix from GUE Figure 17. Histogram for eigenvalues of a 100 × 100 random matrix pulled from the GUE. The histogram is generated from the eigenvalues of 10 5 random samples.
The off diagonal entries are complex entries of the form x + iy and the real and imaginary parts of the entries are pulled from the following Gaussian distribution: We can diagonalize all 10 5 matrices to get 10 7 eigenvalues and we will obtain the histogram (i.e. averaged spectral density) for the eigenvalues depicted in Figure 17. From our sample of eigenvalues we can determine the NNS statistics. An interesting quantity to analyze is the sample averaged spacing of nearest neighbors, and how this changes between which nearest neighbor pairs we choose. We order the eigenvalues so that E 1 ≤ E 2 ≤ · · · ≤ E 100 , then we define the average spacing between the i-th pair as: where i = 1, 2, 3, .., 99. This represents the average (over 10 5 samples) of the spacing between nearest neighbor pairs throughout the spectrum. We can do a discrete plot of this as a function of i to obtain Figure 18. We can see that for eigenvalue pairs located near the center, the average spacing between nearest neighbor pairs is roughly constant but as we approach the edges the spacing changes more quickly between one pair and the next.   Figure 18. We plot the average (over 10 5 samples) spacing, ∆ n , between the n-th nearest neighbor pairs of eigenvalues in the GUE (100 × 100 matrices). We can see that near the edge of the GUE spectrum, the average spacing between nearest neighbor pairs changes quickly. Near the centre of the spectral density of the GUE the average spacing changes slowly.
average spacing changes slowly). Therefore, we expect the model introduced in Section 2.1 to give an accurate description of the spectrum of Gaussian ensembles near the centre of the spectrum and far from the edges 34 . Although, we do not expect our model to describe the edges of Gaussian ensembles to a high degree of accuracy it is still interesting to explore how the form factor of our model spectrum compares to the form factor of an actual Gaussian ensemble as we include eigenvalues near the edge of the spectrum.
To do this we define the following numerically calculated quantity: This is computing the form factor in the GUE which contains the a-th through b-th eigenvalues anneal averaged over S samples. This is essentially describing a "microcanonical" anneal averaged form factor which focuses on eigenvalues within a window [a, b]. We want to compare this numerical result to our model which involves the i.i.d assumption on the NNS distribution. We must specify the following two parameters c and σ for the NNS distribution given in Eq. (3.13). Since we are working with the GUE we should set c = 2 which leaves us one free parameter σ. Recall that σ controls the average spacing between eigenvalues in our model. In particular, for our model the with c = 2 the relation between the average spacing, ∆, and the parameter σ is simply: We will fix σ by analyzing the ∆ i 's in the numerical simulation of the GUE. As we already showed in Figure 18 the average spacing between eigenvalue pairs is roughly constant near the centre of the spectrum. Motivated by this observation we define following quantity which involves data of the eigenvalues E a , E a+1 , ..., E b−1 , E b : In the above equation, a and b label the index of the a-th and b-th eigenvalues (which . Then for a particular set of eigenvalues in the GUE with indices in the closed interval [a, b] we identify the σ =σ [a,b] . This gives the following NNS distribution for our model spectrum: We use this NNS distribution along with N = b − a and plug all this into Eq. (2.10) and plot it along-side the "microcanonical" GUE form factor defined by Eq. (C.4). At infinite temperature we get the plots given in Figure 19. The top left plot in Figure 19 depicts the numerically averaged form factor of the GUE (blue line) for the entire spectrum (i.e. it includes E 1 to E 100 ). We compare with our i.i.d model form factor (the yellow line) with N = 99 and NNS given by Eq. (C.7) with σ [1,100] ≈ 0.345. We can see that the plateau of the model matches the plateau of the model but at earlier times it deviates away from the GUE.
The top right plot depicts the GUE form factor with the edges of the spectrum cut off -only retaining the eigenvalues E 20 to E 80 . We compare with the i.i.d model form factor with parameters N = 60 and σ [20,80] ≈ 0.292. We see that the plateaus still match but now the early time behaviour is matching more closely than before.
The bottom left plot depicts the GUE form factor with the edges of the spectrum cut off even further than, only retaining the eigenvalues E 30 to E 70 . We compare with the i.i.d model form factor with parameters N = 40 and σ [30,70] ≈ 0.285.
Finally, the bottom right plot depicts the GUE, retaining the eigenvalues E 40 to E 60 . We compare with the i.i.d model form factor with parameters N = 20 and σ [40,60] ≈ 0.281.
Together the four plots in Figure 19 plots of the form factor illustrates our claim that our model captures the behaviour spectrum of Gaussian ensembles near the centre of the semi-circle far from the edge. We can see that as we shrink our window of eigenvalues closer to the centre of the spectrum the form factor of the truncated GUE (blue lines) matches more closely with the form factor of our model (yellow line).

D Local Extrema of Form Factor for Gamma Distribution NNS
In this section, we will analyze Eq. (3.19) which is given by: By inspection, it is clear that t = 0 is a valid solution and corresponds to an extremal point. Furthermore, in the limit that t → ∞, the equation is also satisfied, this represents the saturation of the form factor to a constant value (the plateau phase). Aside from these trivial solutions we should expect at least one other solution in the case when q > 0. This is because for q > 0 there is repulsion and therefore a ramp after the initial dip downward. This will give a local minimum in the form factor. Furthermore for very large values of q 1 we also expect oscillations to exist which will give further extremal points in the form factor.
To understand where these non-trivial extremal points occur it is useful to write the complex variable z in the polar coordinate representation: In this representation, Eq. (3.19) can be written in the following form: This is difficult to solve in general so we will consider the special cases when q = 1, 2, 4 where we can solve the equation exactly. We will also make some estimates in the q 1 regime. q = 1 Case: In this case one will find five distinct solutions. One solution is the trivial solution when t = 0, the other four are non-trivial but only one is real and positive for βσ > 0, it is given by: (D.4) q = 2 Case: In this case one will find seven distinct solutions. One solution is the trivial solution when t = 0, the other six are non-trivial but only one is real and positive for βσ > 0, it is given by: In both cases we can see that up to a pre-factor the local minimum before the ramp at low temperatures scales with β this is contrast to the high temperature case where it scales as β 1/4 . q = 4 Case: In this case we will have two non-trivial real solutions which have closed form expressions. The exact expressions are complicated so we will not write them here, instead we will write the solutions as series expansion for small and large βσ. We start with the smaller solution which represents the initial dip in the form factor before the ramp given by: The next real solution corresponds to the small "kink" in the form factor that connects the ramp to the plateau which is also seen in GSE. The expression for the time when the kink occurs is well approximated in all temperature regimes by the following expression: Overall, we see that for q = 1, 2, 4 the inital dip before the plateau has the same temperature dependence up to an order one pre-factor. Furthermore, we are able to show that for the q = 4 case the form factor has an additional local maximum which represents the "kink" which connects the plateau to the ramp seen in Figure 9. q 1 Case: We know that for very large values of q, the form factor exhibits an initial dip followed by oscillations before saturation to a plateau. In this case, we should expect to find more than one non-trivial extremal point which describe the peaks and troughs of the oscillations. In the case when q 1, Eq. (D.3) at leading order will read: A q − A −q sin(qθ) = 2 sin(θ). (D.8) If we define x = tσ/(1 + βσ) we can explicitly write: (1 + βσ) q (1 + x 2 ) q/2 − (1 + βσ) −q (1 + x 2 ) −q/2 sin [q arctan(x)] = 2x √ 1 + x 2 . (D.9) We can clearly see that the left hand side of Eq. (D.9) will diverge like x q sin(qπ/2) as x → ∞. The left hand side saturates to a value of 2 in the x → ∞ limit. This means that there are a finite number of solutions for a given q 1 -which is in agreement with the plots we have made. Clearly x = 0 is a solution, but we expect there to be more.
To understand how many solutions there are, we start by analyzing the number of roots (excluding x = 0) to the expression on the left hand side of Eq. (D.9). This is simply given by the roots of sin[q arctan(x)] which occur when q arctan(x) = mπ where m ∈ N ∪ {0}. Since 0 ≤ arctan(x) ≤ π/2, the number of zeros is controlled by q. In particular, we denote N root as the number of zeros (excluding the zero at x = 0) we have: (D.10) Now consider the right hand side of Eq. (D.9). Since the right-hand side is a monotonically increasing function which takes on values in the interval (0, 2). It is not difficult to see that for a sufficiently large q there will be at least N root − 1 non-trivial solutions to Eq. (D.9) 35 . Whether there is an additional solution depends on the limits as x → ∞. In particular, one can argue that the number of solutions (excluding the x = 0 solution) is given by: .

(D.11)
Now that we have quantified the number of solutions we estimate that the solutions for sufficiently large q occur roughly at: x n = t n σ 1 + βσ = nπ q , n = 0, 1, 2, 3, .., N sol. (D.12) To illustrate how accurately this approximation captures the location of the extremal points of the form factor described by Eq. (3.17). We make Figure 20, which depicts the averaged form factor with q = 50 and β = 0.01 with the red dots representing the value of the form factor at t n σ = (1 + βσ)nπ/q = (1.01)πn/50. We can see that our approximation for the location of local extrema representing the peaks and troughs of the regular decaying oscillation is not perfect but it does give a reasonable estimate. Similar plots can also be made at other temperatures as well.
Based on these results we conclude that the period of the oscillations in the form factor are roughly given by: (D.13) 35 By sufficiently large q, we mean that for a given βσ we must have q such that the left hand side of Eq.

E Average Spectral Density For Gamma Distribution Spacing
In this section we compute the average spectral density for the Gamma Distribution spacing distribution given in Eq. (3.15). To begin, we can write the JPDF for the energy levels given as: (E.1) Using this we can write the average spectral density as: With some work we can derive the following identities: Using the identities above with some additional work one will find: This gives us the average spectral density given by Eq. (3.21). Furthermore, setting q = 0 we will recover the spectral density for the Poisson NNS model discussed in Section 3.2.

F NNS Distribution for Oscillator with Chaotic Interactions 2 × 2 Case
We go over the calculation of the NNS distribution for the following matrix: where the real entries are drawn from the following distributions: Where x 1 , x 2 , x 3 , and x 4 are real random variables which follow the following distributions: by considering what happens when = 0. In this case, there is no perturbation and the Hamiltonian is diagonal. If we took r samples of the exact same diagonal matrix and ordered all the eigenvalues in a list we would find there are exactly |A p | = rΩ(p) eigenvalues with the value E(p). Adding a small perturbation would lead to a small spreading in the eigenvalues within each degenerate sector so the set A p is the set of eigenvalues that have "split" in the p − th degenerate sector. We define the numerical width of the p − th sector as: min(A p ) = E k min (p) . (G. 3) The numerical average energy within the p − th sector is given by: Applying this numerical procedure we obtain the following plots which compare various numerical and estimated quantities when r = 100 and N = 22. ϵ/ω0 = 10 -1 Figure 21. Above are plots of the average energy within degenerate sectors labelled by p at various coupling regimes ranging from the weakly coupled regime with /ω 0 = 10 −4 to the strongly coupled regime with /ω 0 = 10 −1 . The solid circles are numerical computations of the average and the "+" is the energy of the degeneracy sector in the free oscillator case.
We see that in the weakly coupled regime which is given by the top two plots of Figure  21 the numerical average of the energies within the sector p − th sector is close to the degenerate energy of the free oscillators which is consistent with predictions of first order degenerate perturbation theory. When /ω 0 = 10 −2 we can start to see small deviations in the average energy from the unperturbed energy. By the time we get to the strongly coupled regime where degenerate sectors strongly mix there are large deviations. In the weakly coupled regime which is given by the top two plots of Figure 22, the numerical width of the spread in the energies within the p−th sector is close to the estimate given by Eq. (4.16) which is obtained by applying first order degenerate perturbation theory. When /ω 0 = 10 −2 , we can start to see small deviations from the estimate. By the time we get to the strongly coupled regime where degenerate sectors strongly mix there are large deviations.