Eigenbranes in Jackiw-Teitelboim gravity

It was proven recently that JT gravity can be defined as an ensemble of L x L Hermitian matrices. We point out that the eigenvalues of the matrix correspond in JT gravity to FZZT-type boundaries on which spacetimes can end. We then investigate an ensemble of matrices with 1<<N<<L eigenvalues held fixed. This corresponds to a version of JT gravity which includes N FZZT type boundaries in the path integral contour and which is found to emulate a discrete quantum chaotic system. In particular this version of JT gravity can capture the behavior of finite-volume holographic correlators at late times, including erratic oscillations.


Introduction
In finite-volume holography, there is a deep tension between discreteness of the boundary theory, and quasi-normal decay in the holographic bulk. This tension is one version of the information paradox, due to Maldacena [1]. A bulk quantum gravity explanation for this behavior remains to some degree an open question, though important steps have been taken in [2,3,4]. More in general it is not yet completely clear how semiclassical gravitational arguments are to be augmented in a way that resolves the information paradox, for recent progress see though [5,6,7,8]. Two-dimensional Jackiw-Teitelboim (JT) gravity [9,10] has attracted a lot of attention in recent years. This is largely due to its exact solubility and its relevance as the low-energy universal sector of the SYK model [11,12,13,14,15]. JT gravity comes in different versions, as we will highlight in section 2. Most activity thus far has been in a version of JT gravity dual to Schwarzian quantum mechanics, which comes with a topologically trivial bulk. The Schwarzian has a continuous spectrum, so this version of JT gravity has fairly little to do with holographic discreteness. More recently there has been interest in a version of JT gravity that includes a sum over higher genus topologies in the bulk [4]. This model can be defined non-perturbatively as a double-scaled matrix integral. It is therefore an ensemble average over discrete systems. Due to the averaging, the spectrum remains continuous. Nonetheless the averaging does not eradicate all traces of discreteness, in particular late-time holographic correlators do not generically decay to zero in this version of JT gravity [4,16,17,18]. Motivated by these developments, we want to point out that there exists a further alternative version of JT gravity which resembles the behavior of a single discrete holographic system. In particular it captures the behavior of the spectral form factor of a discrete system for all times, including the late-time erratic oscillations. The new feature is to include of set of fixed energy boundaries in the gravitational path integral on which Riemann surfaces can end, to be distinguished from the asymptotic boundaries. Each energy labels corresponds to the energy of a state in the discrete system we aim to emulate. These boundaries, which we will refer to as eigenbranes, are related to unmarked FZZT boundaries and correspond to fixed eigenvalues in the matrix ensemble of [4].
This work is structured as follows. In section 2 we introduce a simple probe of discreteness. Using the freedom to choose the contour of the path integral over metrics in quantum gravity, we discuss three possible definitions of JT gravity. One has only disk topologies, one has all topologies and is completed as an ensemble of Hermitian random matrices. The final one includes eigenbranes and is dual to an ensemble of Hermitian random matrices with a large number of eigenvalues fixed. In section 3 we review and discuss calculational techniques for spectral densities in matrix integrals, or multi-boundary correlators in gravity. In section 4 we consider an ensemble of random matrices with certain eigenvalues kept fixed, and derive its perturbative interpretation as a gravity path integral with surfaces ending also on a number of fixed-energy boundaries. We proof the extent to which the resulting version of JT gravity resembles a discrete system. In particular, we show that the spectrum essentially reduces to a series of delta functions, and we show how different asymptotic boundaries essentially disconnect. In section 5 we briefly discuss a possible gravitational interpretation of the delta functions as due to boundary mergers [19] and touch on the generalization to JT supergravity [20]. The appendices contain some of the technical material.

Diagnosing discreteness
Consider a particular discrete maximally chaotic quantum mechanical system with an Ldimensional discrete Hilbert space: We choose the model in such a way that its coarse-grained spectrum matches the JT spectrum on the disk up to some large energy cutoff Λ 1 [11,14,2,21]: The system is taken to be quantum chaotic, because we aim for it to be dual to quantum black holes [2,22,23,24]. Specifying further to a system without time-reversal invariance, this implies its local level statistics should be those of a random matrix taken from an appropriately rescaled Gaussian unitary ensemble (GUE) [25]. At low energies E < Λ or late times, we might expect that this discrete quantum mechanical system has an effective pure JT gravity bulk description. 1 In the remainder of this work we collect evidence in favor of this.
We will be interested in thermal correlation functions, for example the two-point function: 2 At early times, the Fourier transform is insensitive to the fine structure of the spectrum, and coarse-grains. As a consequence, the thermal two-point function will be well approximated by a JT disk calculation. The late-time Fourier transform on the other hand is highly sensitive to the fine structure. Therefore late-time correlators are in general suitable probes of discreteness [1].
In particular we will be interested in the simplest such probe, a local version of the twopoint function where we integrate both E 1 and E 2 over only a narrow energy interval bin(E). We take 1/ρ 0 (E) |bin(E)| 1, this is small enough such that ρ 0 (E) is approximately constant, and large enough such that it contains a large number of eigenvalues. Because the system (1) is quantum chaotic, we known from the eigenvalue thermalization hypothesis (ETH) that | E 1 | O |E 2 | 2 are also slowly varying in this bin [26,27,17,2]: We are interested only in the t-dependence of this object, so we can remove the constant matrix element and take β = 0. We are left with a local version of the spectral form factor [3]: Labeling the eigenvalues of (1) within bin(E) as λ 1 . . . λ N , with 1 N L, we get: A log-log plot for a representative sample gives: 3 S E (t) = As compared to the usual spectral form factor, there are additional low-frequency oscillations, but the general shape is the same [2]. In particular, it has the same ramp-and plateau structure including erratic oscillations, with plateau height N . We will reproduce these erratic oscillations via a bulk JT gravity calculation in section 4.2.

Models of JT gravity
JT gravity is a model of 2d dilaton gravity with action [9,10]: The Euler character χ(M) comes from the Einstein-Hilbert term in 2d. 4 The quantity S 0 = Φ 0 /4G corresponds to the extremal entropy, and is a free parameter from the gravity point of view. Integrating out Φ localizes the metrics g on hyperbolic Riemann surfaces, or patches of the Poincaré disk, with asymptotically NAdS 2 boundary conditions. This boils down in JT gravity to fixing the total length of the asymptotic boundary to β/ , and the boundary value of the dilaton to 1/2 [13,14,15]. Defining the integration space of g boils down to specifying which surfaces to count. For the JT gravity partition function, which corresponds to the insertion of a single holographic boundary in the path integral, we have schematically: The spectral form factor is related in JT gravity to a correlation function Z JT (β 1 , β 2 ) with two asymptotic boundaries of respective lengths β 1 / and β 2 / : From this, one calculates the local spectral form factor (5) in gravity as: In the remainder of this section we highlight how different definitions of ??? can lead to structurally very different theories, using the spectral form factor (11) as probe. The different models we will discuss can be specified by the integration space in the path integral over metrics: In particular, we want to point out that the last definition which takes the energies λ 1 . . . λ n as input, resembles the discrete system with spectrum (1).

Version 1: Disks
The simplest definition is to restrict ??? to Riemann surfaces which are topologically disks. The spectrum and correlation functions of JT gravity on the disk have been extensively 4 For a 2d manifold of genus g with b boundaries we have studied in recent years, resulting in several complementary perspectives [28,29,30,31,32,33,34,35,36,37]. Its spectrum is continuous [11,14,2,21]: To probe the spectrum, we focus on the spectral form factor (5). For disk topologies, the gravitational spectral form factor factorizes ρ(E 1 , E 2 ) = ρ 0 (E 1 )ρ 0 (E 2 ), and we get: This dependence correctly captures the part of the curve (7) before the dotted red line, including the relatively slow oscillations. At later times, the Fourier transform is able to distinguish the coarse-grained disk spectrum (13) from the discrete spectrum (1).

Version 2: Genus expansion and random matrices
A second possible definition of JT gravity is to allow for Riemann surfaces of arbitrary genus to end on the asymptotic boundaries. This version was introduced and discussed in [4]. 5 The JT gravity partition function now has a genus expansion: The same is true for the spectral density, and all other observables: We will refer to (15) as the "perturbative" definition of JT gravity. 6 It is very feasible to calculate each term in this series, as we briefly review in section 3.1. The resulting perturbative series turns out to be asymptotic, and hence requires a non-perturbative definition. We can define JT gravity non-perturbatively as a double-scaled matrix integral with genus 5 Aspects of JT gravity on higher genus Riemann surfaces were also discussed in [38,16,39,40,20,17]. Let us note that an alternative mathematically consistent definition corresponds to integrating over Teichmüller space instead of the moduli space of Riemann surfaces [38]. This is more natural from a first-order BF point of view, but it gives divergences at higher genus and hence is not a workable version of JT gravity as sum over topologies. 6 It is perturbative in the string coupling e −S0 but non-perturbative in the Newton constant G ∼ 1/S 0 . zero spectral density (13) [4]. 7 In the weakly coupled regime e S 0 1, it turns out one can neglect all perturbative contributions, and that the leading correction is nonperturbative in the string coupling: 8 Such oscillatory non-perturbative contributions will in general not have a geometrical interpretation as counting Riemann surfaces. 9 The partition function of an ensemble of L × L Hermitian matrices with bare potential V (M ) is defined as: A more convenient way to write this is in terms of the eigenvalues λ i of the matrices: 10 Here ∆(λ 1 , . . . , λ L ) is the Vandermonde determinant, accounting for eigenvalue repulsion. An intuitive way to think about about such matrix integrals is as the steady state of the Brownian motion of L charged particles in an external potential V (x), a so-called Dyson gas [43,44,45,46,47,48]. The Vandermonde determinant then represents the electrostatic repulsion.
Typical observables in the matrix model are products of the spectral density ρ(E) or the macroscopic loop operator Z(β): Correlators are calculated as ensemble averages, for example: 11 7 see also [41,42]. 8 This is only true when E e −2S0/3 . We will operate under this assumption throughout the main text. Otherwise, we have to resort to an exact analysis of the Airy model, as we do in appendix A.4. 9 The nonperturbative contribution in the forbidden region E < 0 does seem to count Riemann surfaces that end on the appropriate boundaries, stretching between asymptotic boundaries and ZZ branes in the bulk [4]. These are to be distinguished from the FZZT branes discussed in the remainder of this work. 10 Any contour C represents a definition of a matrix model. In case of JT gravity, for stability reasons the contour cannot be chosen to extend along the real energy axis all the way up to −∞ [4]. The part of the contour that deviates from the negative real axis is not important for the content of this work though, so we will drop the subscript in what follows. See appendix A.4 for fixed eigenvalues in the forbidden region. 11 They are normalized such that 1 = 1.
The two-level spectral density is defined as ρ(E 1 , E 2 ) = ρ(E 1 )ρ(E 2 ). Ensemble averaging leads to correlation as connected contributions, for example for the two-loop operator defined as Z(β 1 , β 2 ) = Z(β 1 )Z(β 2 ): Comparing to the perturbative JT gravity definition of Z(β 1 , β 2 ) in (10), which counts all Riemann surfaces that end on the two asymptotic boundaries, one sees that connected correlators correspond to connected geometries: Again, it is not hard to calculate both perturbative and nonperturbative contributions to ρ(E 1 )ρ(E 2 ) ), see section 3, appendix A and [4]. The only significant perturbative contributions are due to the disconnected disks, and annulus connecting the two boundaries.
There are also significant nonperturbative contributions: 12 The second contribution here is of the same oscillatory type as (17). Unlike that contribution though, it isn't particularly small due to the multiplicative pole, and cannot be neglected in our analysis.
The spectral form factor is calculated as in (11). The ρ 0 (E 1 )ρ 0 (E 2 )-contribution in (23) gives the power-law decay (14). The Dirac-delta yields the constant plateau contribution N , and the other contributions add up to the sine-kernel (48) which gives a variant of the ramp at late times: 13 , 12 In addition to these contributions we find a generalization of the wiggles (17) in appendix A. Such wiggles are genuinely small corrections though, unlike (24) and are not relevant for our discussion. 13 In more detail, the contribution of the sine kernel is: This is a low-frequency filtered version of the usual ramp. Qualitatively, at t 2πρ(E) there will be significant smoothening of the onset of the usual ramp. In the regime of interest where N 1, the kernel acts as a Dirac-function and one obtains the linear ramp with plateau time t plateau = 2πρ(E) . This resulting function follows (7) before the dotted red line, and the blue curve at late times.

Version 3: Eigenbranes
The version of JT gravity discussed above is a double-scaled matrix integral. This means we take L → ∞ and simultaneously zoom in on a region E < Λ near the edge of the spectrum, keeping the total average number of eigenvalues 1 N L in this region fixed. We can visualize this as: The blue region represents JT gravity with spectral density (17). By definition, our discrete system (1) can be thought of as a single Hamiltonian M of such a matrix ensemble. Let us denote its lowest N eigenvalues by λ 1 . . . λ N . We expect that the IR behavior of our system (1) is accurately described by a modified matrix ensemble where N eigenvalues are kept fixed to λ 1 . . . λ N . This corresponds to a Dyson gas of charged particles equilibrating in an external potential around N static point charges. These fixed charges repel the charged gas, resulting in a void. We expect the spectral density of this new ensemble to essentially follow the spectrum of the discrete system (1) for E < Λ and that of the original ensemble (27) for E > Λ: In the remainder of this work, we will make this picture precise and pinpoint its JT gravity interpretation. In particular, we will see that each eigenvalue λ corresponds to a fixed energy boundary with label λ hovering in the Euclidean bulk. The contour (12) in the gravitational path integral is hence over all Riemann surfaces that end on the union of the asymptotic boundaries and on N fixed energy boundaries with labels λ 1 . . . λ N , as shown in (55). This version of JT gravity is able to capture the IR discreteness of (1). In particular, we will recover the spectral form factor (7) including erratic oscillations. In this picture, smooth geometry in the bulk is never in jeopardy: it is provided by our ignorance of the UV part of the system (1), which corresponds to the L N eigenvalues that remain in the continuum of the matrix integral.

Multi-boundary correlators
This section prepares for section 4, where we will encounter multi-spectral density correlators ρ(E 1 ) . . . ρ(E n ) . We discuss an efficient way to calculate all significant perturbative and nonperturbative contributions for e S 0 1 based on [49].

Genus expansion
In JT gravity, it is natural to consider fixed-length boundary conditions, as discussed around (9). These correspond to the insertion of macroscopic loop operators in the matrix integral, and are the Laplace transforms of the multi-spectral densities: This relation is an efficient tool to calculate the perturbative contributions to ρ(E 1 . . . E n ) [4]. For example, the genus g contribution to the n-loop correlation function is: is the volume of the moduli space of the n-holed sphere with g handles. This Weil-Petersson volume is a polynomial in b 2 1 , b 2 2 etc and is easily calculated recursively [50]. The single twisted Schwarzian partition function Z JT (β, b) is just a Gaussian: 14 The Gaussian b-integrals yield polynomials in β 1 , β 2 etc multiplied by (β 1 . . . β n ) 1/2 . Inverse Laplace transforming then gives us the spectral densities, which are polynomials in 1 The only exceptions to this polynomial behavior are the disk-and annulus topologies, for which the Weil-Petersson volumes V 0,1 (b 1 ) and V 0,2 (b 1 , b 2 ) are undefined. 15 The disk density of states is (13), and the annulus amplitude is: 14 See [21,4,36] for the evaluation and interpretation of the twisted Schwarzian partition function. 15 One could take Its contribution to the two-level spectral density is: where we approximated |E 1 − E 2 | 1. We see that with the exception of the disks and annuli, all the perturbative contributions are small corrections as long as we stay far enough from the spectral edge. 16 The annulus contribution itself can become large and comparable to the size of the disk contribution ρ 0 (E 1 )ρ 0 (E 2 ) for |E 1 − E 2 | ∼ 1/ρ 0 (E), the typical eigenvalue spacing in the ensemble. For much larger separations it is negligible. As ρ 0 (E) ∼ e S 0 and we are steering clear of the spectral edge, the only significant contributions from the annulus arise well within the range |E 1 − E 2 | 1. This validates using the second equality in (32) throughout.
In conclusion, when away from the spectral edge and in the regime e S 0 1, all perturbative contributions to the spectral densities are negligible except for those associated with the disk-and annuli topologies.
As the inverse Laplace transforms of fixed length correlators in JT gravity, the spectral densities ρ(E 1 . . . E n ) correspond to imposing certain fixed energy boundary condition at the n boundaries of the Riemann surfaces. These boundary conditions are closely related to the FZZT boundary conditions in Liouville theory [51]. 17

Exact answer
On top of the perturbative contributions discussed in the previous subsection, an exact matrix integral analysis reveals nonperturbative contributions to ρ(E 1 . . . E 2 ). We will now discuss an efficient way to calculate all significant contributions in the regime e S 0 1, perturbative and non-perturbative, with more details in appendix A.
Brane operators in our matrix ensemble are defined as: We can use this to extract and write the dependence on λ 1 of the Vandermonde determinant in (19) as: 16 We should take E e −2S0/3 , which ensures the topological suppression prevails over the poles in the spectral densities for small energies, the weakest and hence most important of which is (E 1 . . . E 2 ) −3/2 . This is confirmed from the exact analysis near the spectral edge in appendix A. 4 where E e −2S0/3 is the condition for higher loop contributions to the Airy function to become negligible [49]. 17 See for example [4,52]. or more generally: This essentially decomposes the measure of the matrix ensemble (19): For the branes that feature in this formula, the product in (33) is over the L − n remaining eigenvalues. 18 This basic formula allows us to extract exact formulas for spectral densities, which we can easily calculate exactly. Let us demonstrate this case-by-case.
• 1 eigenvalue. We can use the symmetries of the ensemble, and the property (36) to rewrite (21) as: From this, we read off the spectral density [4]: 19 Both sides of this equality can be calculated independently in JT gravity: in appendix A we calculate the double brane correlator using techniques of [49], and the spectral density was calculated using related techniques but via a different computation in [4]. We find: Comparison gives us a recursion relation for the matrix integral partition function at large L: 20 This will enable us to eliminate any dependence on Z L from the calculations that follow.
• 2 eigenvalues. The 2-boundary correlator decomposes as: Using the recursive formula (40), we end up with: These types of formulas are well-known in the random matrix literature. In fact they are referred to simply as the correlation functions [25]: 21,22 These are smooth functions. The operators in (43) represent the repulsive force exerted by a set of charges charges at E 1 . . . E 2 on the remainder of the Dyson gas.
• 3 eigenvalues. An equally easy calculation holds for for the 3-level spectral density. We find: This is readily generalized to any number of boundaries. 21 The constant in formula (6.1.1) of [25] is 1/Z L . The average in (6.1.2) generates a factor Z L−n , the recursion relation (40) removes the combinatorial prefactors. 22 Brane operators in the matrix integral are closely related to exponentiated spacetimes attached to a brane, see appendix A. In this sense, formulas of the type (42) are quite surprising, since they say that a brane-pair correlator actually corresponds to a single (fixed-energy) boundary in gravity.
The delta-functions that appear in these expressions of this kind are contact terms. Whereas a geometric interpretation of the correlation functions R(E 1 . . . E n ) is obvious from the discussion of section 3.1, the interpretation of these terms is somewhat obscure. We will come back to this in the concluding section 5.
It is convenient to extract from the correlation functions R(E 1 . . . E n ) the fully connected contribution T (E 1 , . . . E n ) known as the cluster function. The remaining disconnected pieces are then products of cluster functions at lower values of n. For example [25]: 23 Following the logic around (23), one immediately deduces that the clusters T (E 1 , . . . E n ) correspond to the nonperturbative completion of the gravitational genus expansion starting with the n-holed sphere. 24 The cluster functions have the property that they vanish when the spacing of two of its arguments is large compared to the average eigenvalue spacing. This means the only significant contributions of the cluster functions to the correlation functions for all energies in a cluster. The perturbative disk-and annuli contributions discussed in section 3.1 are part of the terms T (E i ) respectively T (E i , E j ) that contribute a generic correlator R(E 1 . . . E n ). As mentioned earlier, these are the only significant perturbative contributions to R(E 1 . . . E n ) away from the spectral edge.
An exact calculation of the correlators R(E 1 . . . E n ) in JT gravity reveals a set of nonperturbative contributions similar to those in (24). An efficient way to calculate these exactly in JT gravity is via formula (43). We do so in appendix A in detail for R(E 1 ) and R(E 1 , E 2 ) and discuss certain aspects of the calculation for R(E 1 , E 2 , E 3 ). The general trend is the appearance of significant non-perturbative contributions to R(E 1 . . . E n ) of the type: 23 The minus signs are convention [25]. 24 The precise version of this statements follows from the decomposition of the correlators ρ(E 1 ) . . . ρ(E n ) into cluster functions, including contact terms. The resulting cluster functions correspond precisely to the nonperturbative completion of the n-holed sphere genus expansion, which will generate the same contact terms. For example, the three-holed sphere with all corrections gives: This follows directly from (44), but also follows intuitively from the discussion on merging boundaries in 5.
It is convenient to extract from this the cluster functions, which as explained above can be evaluated for |E i − E j | 1. We find: The sine kernel S(E 1 , E 2 ) also appears in higher clusters, for example: This is very unsurprising. It is a widely held conjecture [25] for any Hermitian matrix ensemble that cluster functions are exactly equal to the universal GUE cluster functions when its arguments are close together |E i − E j | 1. The latter are known exactly [25] and feature only the sine kernel. In the brane calculations these arise due to the contributions of the type (47). The calculations of appendix A merely reassure us that this conjecture is true in JT gravity. We are then free to ship in the GUE clusters to calculate ρ(E 1 ) . . . ρ(E n ) in JT gravity.

Fixing eigenvalues or introducing boundaries
In this section we investigate a matrix ensemble with a series of eigenvalues fixed to consecutive ones of (1), and specify the integration space in the JT gravity path integral over metrics (12) associated to this ensemble. The specific contour follows from formula (43) combined with (36): each fixed eigenvalue corresponds to an additional fixed-energy boundary in the bulk on which Riemann surfaces can end. A matrix ensemble with n eigenvalues fixed to λ 1 . . . λ n (assumed all different) is obtained from the original ensemble (19) by including appropriate deltas in the measure: 25 The partition function replacing (19) is then: Here we used (36) and the generalization of (44) to n boundaries. Notice that the contact terms vanish because the eigenvalues of (1) are all different. Perturbatively, this partition function is counting Riemann surfaces of the type: where n eigenbrane boundaries are present, but no asymptotic boundary insertions.

Delta spikes and a void
The expectation value of the spectral density ρ We immediately obtain: This is a conditional probability. As announced, this corresponds to a version of JT gravity where each fixed eigenvalue of the matrix integral translates into the introduction of a fixedenergy boundary on which Riemann surfaces in the path integral are to end. As explained before, in the genus expansion disks and annuli dominate the regime of interest. From (54) we read off the type of geometries contributing significantly to the JT gravity path integral: There are also multi-annulus configurations where eigenvalue boundaries connect to other eigenvalue boundaries. Three-holed spheres and handle-body geometries contribute, but not significantly.
Using a generalization of formula (44) we can rewrite (54) as: At this point our discussion of the previous section comes into play: we can immediately write down the exact answer for a given n using the cluster functions (48) etc.
As a consistency check on the normalization, we can take the integral over E of (56) using the result of Appendix B: We see that the number of eigenvalues in the smooth continuum is exactly down by n as compared to ρ(E), and these eigenvalues are accounted for by the delta functions. As discussed in section 3, the contributions from the annuli connecting the asymptotic boundary to the eigenvalue boundaries is negligible when |E − λ i | 1/ρ(E), and the same holds for all nonperturbative contributions. Therefore, all effects due to the fixed eigenvalues are short-ranged and one has: 26 To fully understand the physics in the exact formula (56), let us do a small case-by-case study. 27 • 1 eigenvalue. We have from (56): Close to the fixed eigenvalue this looks like: This exhibits eigenvalue repulsion: the fixed charge repels the particles of the gas, as modeled here by the Vandermonde factor (E − λ) 2 .
• 2 eigenvalues. Using the GUE cluster functions (48) in (56), we find a less elegant answer for the case of two fixed eigenvalues. 28 A plot close to the fixed eigenvalues is much more intuitive: There is a relatively low probability for another eigenvalue to be found in between λ 1 and λ 2 , provided they are close enough. In general we can think of the initial coarse-grained density as a low-frequency approximation to the series of delta-functions. The maximal frequency here is the typical eigenvalue spacing 1/ρ(E). We see therefore manifestly that we are not changing any early-time t ρ(E) physics by fixing eigenvalues.
• A bin of eigenvalues. It is not hard to plot (56) exactly for an increasing number of consecutive eigenvalues of (1) in some region. For example, for n = 8 we find: We are starting to see the features claimed in formula (28). Firstly, fixing a large number of consecutive eigenvalues will create to good approximation a void in the continuum spectral density in the interval I where the eigenvalues are situated: Secondly, the effect is not felt far outside of I: the effect dies out over a range ∼ 1/ρ(E) 1.
As mentioned before, in the region closer to the spectral edge, where the spectral density ρ(E) changes rapidly, we can no longer trust the sine-kernel type GUE cluster functions (48). Fortunately, in that region, we have available the exact results of the Airy model. Using the method of [49] to calculate brane correlators, it is straightforward though slightly tedious to recover the known Airy cluster functions. We do so in appendix A.4 for the case T (E 1 , E 2 ) and recover the Airy kernel [25]. We then study the spectral density with one fixed eigenvalue close to the spectral edge. The behavior is very similar to that of (60). It would be straightforward to extend this to multiple fixed eigenvalues, but we will refrain from doing so.
All this points in the direction of the picture (28): by inserting the 1 N Λ fixed energy boundaries discussed in section 2.1.3 in JT gravity, we get a version of JT gravity with a spectral density that is essentially completely discretized in the region E < Λ, matching that of the abstract discrete system (1). We note that it is possible that the calculations of the brane correlators presented in appendix A are more subtle when n ∼ e S 0 . In particular, the limit e S 0 1 used in [49] to obtain the semiclassical brane correlators could be more subtle. It would be valuable to understand if this happens, and how the technical calculation is modified. There is no reason though to expect any qualitative deviations from the picture (28) and our conclusions. In particular we expect no sizeable modification of the correlation function R(E 1 . . . E n ) away from the GUE answer.

Erratic oscillations
To stack up the claim that introducing these eigenbranes in JT gravity allows one to capture the E < Λ features of the discrete system with spectrum (1), we would like to reproduce the local spectral form factor (7) from a JT gravity calculation. For this we will investigate ρ(E 1 )ρ(E 2 ) λ 1 ...λn , with emphasis on the terms that contribute to the plateau region t > 2πρ(E).
Using the ensemble with n fixed eigenvalues (51), one immediately writes down: Geometrically, we are calculating the correlator of two fixed-energy boundaries in a version of JT gravity that has n fixed-energy boundaries hovering in the bulk. The only significant perturbative contributions are due to the disk and annuli, for example: As in (55) the eigenvalue boundaries that don't connect to the asymptotic boundaries don't need to be capped off by disks, there can be annuli between them.
Using the exact formulas for the multi-spectral densities discussed in section 3, we obtain: Again, using the JT spectral density and the GUE cluster functions it is easy to calculate and plot this recursively for increasing n. Numerically investigating the continuous contributions to (67) it quickly becomes obvious that if we fix a large number of eigenvalues of (1), then in the region I where the eigenvalues are positioned, to good approximation: 29 If we take the region I large enough such that |bin(E)| |I| then we trivially recover the discrete version of the local spectral form factor (6) in JT gravity, including all erratic oscillations in (7). We would like to understand in a bit more detail the approach of the local spectral form factor to this erratic behavior though. Let us focus on the plateau region t > 2πρ(E) and take only a few fixed eigenvalues. 30 In the averaged version of JT gravity, the plateau behavior is only due to the first term in (23): In appendix C we point out that only the first and penultimate contributions to (67) are relevant for the spectral form factor at t > 2πρ(E): This formula nicely interpolates between the averaged variant (69) and the discretized variant (68). The first term contributes a constant plateau of height N . 31 The second term generates ever more erratic oscillations for increasing number of eigenvalues: For n = 0 this is the usual random matrix theory answer, for n = N we recover the discrete answer (6).

Disconnection
The connected part of the two level spectral density is defined as: From (68) and (64), we get to good approximation: This is trivial for a discrete system, but it entails a nontrivial equality in bulk gravity.
To appreciate this, consider the geometries that contribute to the two level spectral density ρ(λ 1 ) . . . ρ(λ 2 ) ρ(E 1 )ρ(E 2 ) λ 1 ...λn and compare this to the geometries that contribute to ρ(λ 1 ) . . . ρ(λ 2 ) ρ(E 1 ) λ 1 ...λn ρ(λ 1 ) . . . ρ(λ 2 ) ρ(E 2 ) λ 1 ...λn . If we strip off geometries that contribute to both, we end up in the former with connected geometries such as the annulus between the two asymptotic boundaries. In the latter we are left with configurations where the boundaries are indirectly connected via matching pairs of eigenbranes. The sum of what remains in either quantity is non-zero. We can calculate the exact answer for each quantity independently for increasing n using the techniques of section 3.2. It turns out that these quantities match for E 1 , E 2 ∈ I. This proves the following geometric property: This factorizing property is slightly surprising in the sense that the geometries on the right hand side are never counted in the original perturbative JT gravity path integral prescription for ρ(E 1 )ρ(E 2 ) λ 1 ...λn . One might have expected that connected contributions to the gravity analogue of a discrete system vanish. We find that instead they are non-zero, but their sum can be replaced by a sum over disconnected contributions. From the perspective of the matrix integral this disconnection is more intuitive. If we deplete a large energy region |I| 1 then locally the randomness of the initial Hamiltonian is lost. Indeed, a close relation between geometric disconnection and lack of randomness is generally expected [3,4]. We note that (74) looks somewhat like introducing a "complete set of baby universes" between E 1 and E 2 as hinted towards in [17], though the status of these eigenbranes as "states" in bulk JT gravity is at the moment unclear.

Concluding remarks
It would be interesting to understand what these eigenbranes mean for a Lorentzian observer probing the gravitational bulk. Can he somehow obtain information about the branes hovering deep in the bulk? One way to work towards this would be to investigate boundary correlators in the matrix ensemble, see for example [17,18]. It would be valuable to understand if we can construct bulk observables within JT gravity as a sum over these more complicated geometries, using geodesic localizing, in analogy to the construction of local bulk observables in the disk version of JT gravity [16,54].
We end this work with three remarks.

Boundary mergers
The Dirac delta's that appear in the exact answers for the spectral densities in section 3.2 have an a posteriori interpretation as eigenvalue boundaries merging with the asymptotic boundaries. 32 Considering for example the last equality in (41), the first term can be read as counting Riemann surfaces which end on the merger of the two original boundaries 32 See for example [19,55]. of lengths β 1 and β 2 , resulting in a boundary of total length β 1 + β 2 . Let us pretend here to take that interpretation seriously, and count Riemann surfaces that end on a merged boundary. It is convenient to introduce the JT gravity disk amplitude between a fixed length state |β and a fixed energy state |E : 33 The merger of an asymptotic disk with an eigenbrane results in the genus zero amplitude ρ 0 (λ) β|λ . 34 This merged boundary can connect to the other eigenvalue boundaries, and develop handles. In taking the sum, the overlap β|λ i is a spectator. We end up with a factor that cancels precisely the denominator in (54), and we are left only with β|λ i . As pointed out in section 4, all other contributions to the JT gravity partition function add up to zero. This suggests the net gravitational effect of fixing all eigenvalues (1) in JT gravity is the following replacement: It is tempting to imagine how this could extend to holographic correlation functions. 35 In JT gravity in its gauge-theoretic BF formulation, boundary correlators correspond to Wilson lines traversing the Riemann surfaces [57,33,38,37,58]. The Wilson line separates the Riemann surface into two disconnected pieces, each connected to a piece of boundary: 36 (78) 33 These are the states used in [35,17,56], with |β the Hartle-Hawking state of the JT gravity disk. We have: 34 This is the inverse Laplace transform of the boundary with length β 1 + β 2 with respect to β 2 . 35 See [17,18] for recent discussions. 36 A similar such configuration with a vacuum Wilson line does not contribute to the JT gravity partition function, because the eigenvalues are chosen not to be degenerate, and merging two fixed energy boundaries to a fixed length boundaries results in an amplitude proportional to a Dirac delta on those energies.
Even though the final expression (78) is the two point function of a discrete system for any operator O, it is nontrivial to find one that corresponds to a boundary-anchored Wilson line in the perturbative definition of JT gravity. One hint that (78) might be correct comes from the spectral form factor. This can be obtained from an analytic continuation of the Schwarzian boundary two-point function with β − it → β 1 1 and it → β 2 1. The Wilson line in (78) then effectively pinches off the disk. 37

Other ensembles
The analysis of this work is readily generalized to JT N = 1 supergravity [20]. The general Altland-Zirnbauer ensembles are defined by the integration measure [59]: in terms of two integers (α, β). It was shown in [20] that these are related to JT supergravity either on orientable or orientable plus nonorientable surfaces depending on the choices of α and β. Taking β = 2 and α = 1 is the orientable case (no time-reversal symmetry). For the nonorientable case (time-reversal symmetry), there is a possible divergence from the crosscap moduli space, with the exception of the case β = 2 and α = 0, 2. We focus hence only on these cases. In analogy to (39) and (42) one has: These can be calculated efficiently by generalizing the brane computation of Appendix A to including the crosscap Xcap(e) = − α−1 4 log(−E). Contributions to the spectral density then for example are: In this case: It is not difficult to considered fixed eigenvalues or conversely additional fixed energy boundaries in these models, one merely has to ship in the appropriate cluster functions for these (α, β) ensembles. Close to the spectral edge ρ 0 (E) ∼ 1/ √ E JT supergravity reduces to the exactly solvable Bessel model, which is the super-analogue of the Airy model [20].
A gravitational hint of ensemble averaging?
The statistical ensemble we found from the matrix integral, was interpreted gravitationally in terms of multiple boundaries. Here we illustrate that starting with gravity directly, one can get hints of this underlying ensemble, reversing the logic of this work. We start from a property of the n-boundary correlator in JT gravity (29): Let us take the length of one of the boundaries to zero. The single-trumpet partition function Z JT (β) (30) becomes a delta-function for β → 0. Taking β 1 → 0 therefore localizes on spacetimes where the neck length b 1 vanishes. Due to the twist factor b 1 and the polynomial behavior of the Weil-Petersson volumes, we see that every perturbative contribution vanishes except for the case when the β 1 -boundary is capped off by a disk. We end up with: Doing an n-fold inverse Laplace transform of this equation, we find: Recursively, one gets from this: This suggests to think of ρ(λ 1 . . . λ n ) Z JT (0) −n = w(λ 1 . . . λ n ) as the weight function of a statistical ensemble. This is strengthened by (85) and its generalization to multiple E i : correlators in JT gravity can be calculated as averages in this statistical ensemble. In particular the observable that calculates ρ(E) is extracted from (85) as: This corresponds to the quantity we considered in the main text.

Acknowledgements
We thank Phil Saad for useful discussion. AB and TM gratefully acknowledge financial support from FWO Vlaanderen.

A Some brane calculations
In this appendix, we calculate brane pair correlators of the type ψ 2 (E 1 ) . . . ψ 2 (E n ) in JT gravity, with a single brane defined as (33). As discussed in the main text (43), this is an efficient way to calculate objects such as R(E) and R(E 1 , E 2 ).
We can rewrite the brane operator (33) as: The operator in the exponential corresponds to the insertion of an unmarked fixed energy boundary in JT gravity [4]: This is the precise analogue of an unmarked FZZT boundary brane in Liouville theory [51,60]. Equation (88) is slightly misleading in combination with (89) though. The original brane correlator (33) is an analytic function of E, whereas the FZZT brane (89) has a discontinuity on the positive real axis. Consequently, to each energy E there correspond two different FZZT boundaries in gravity, depending on how we approach the real axis. This is equivalent to specifying the value √ −E for E > 0. Let us introduce a new variable e = i √ E, then √ −E = ±e for E > 0. Depending on this sign, exponentiating the FZZT boundary (89) gives two distinct gravitational brane correlators ψ(±e) . This raises the question which gravitational brane corresponds to inserting the brane operator (33) in the matrix integral. The answer was given in [49]. 38 The brane correlators have an exact expression for finite e S 0 as a Kontsevich matrix integral, or an appropriate JT gravity generalization thereof. 39 For e S 0 1 we can use the method of Laplace on this Kontsevich matrix integral. Depending on whether the energy parameters E 1 . . . E n are positive or negative, different saddles contribute due to Stokes phenomena. Each such saddle, and the loop corrections around it, correspond to a gravitational brane. It turns out that for all energies E 1 . . . E n positive, we need to sum over all possible corresponding gravitational branes with equal weight. For each such saddle, if we furthermore take E e −2S 0 /3 only 38 See also [4]. 39 See [19].
the exponentiation of disks and annuli contributes significantly. 40 We will be working throughout in the regime e S 0 1. In sections A.1 and A.2 we furthermore assume E e −2S 0 /3 . In section A.4 we calculate the correlators close to the spectral edge using the Airy model. The Airy calculations are exact for any e S 0 and by construction coincide with the JT gravity answers for E 1. For e S 0 1 these regions overlap, and we have an exact answer for all E.
we end up with: This matches the result of the resolvent-based brane dipole calculation of R(E) in [4].

A.2 Two pair
Next we calculate the two-brane-pair correlator ψ 2 (E 1 )ψ 2 (E 2 ) . For notational purposes consider ψ 2 (E)ψ 2 (K) with e = i √ E and k = i √ K. Summing the 16 saddles gives: The generic brane correlator is similar to (91): Secondly, there are 8 mixed terms: The remaining 4 terms have opposite signs within each brane pair. This terms requires a double use of L'Hpital's rule: One recognizes the first term as the product of two perturbative disks and the second term as the perturbative annulus (32). Adding these three contributions and multiplying by (E − K) 2 /4π 2 gives the exact pair density correlator for e S 0 1 away from the spectral edge. We can distill from the exact answer T (E, K) in JT gravity: This connected contribution contains the perturbative annulus as first term. The remainder is the non-perturbative contribution to the two-holed sphere (i.e. the annulus).
To uncover the GUE structure (48), we focus on |E − K| 1. 41 This simplifies things considerably: Furthermore: Collecting everything, we find for |E 1 − E 2 | 1: The first term is the well known GUE result. the second term is small and oscillatory. It is the analogue of the wiggles (17) in R(E). For the purpose of our story in the main text, these wiggles are negligible but for the fact that R wiggle (E 1 , E 1 ) = 0, such that R(E 1 , E 1 ) = 0 as demanded by eigenvalue repulsion in the ensemble.

A.3 More pairs
This procedure readily extends to generic n. The perturbative contribution is found by picking opposite signs within each brane pair, only then is there no oscillatory contribution. For example for n = 3 after a tedious three-fold application of L'Hpital's rule one recognizes the perturbative disks and annuli: It is reassuring to see these and only these perturbative contributions appear. Notice for example that there is no perturbative three holed sphere contribution, nor are there handlebody geometries. This is consistent, as discussed in section 3.1 those geometries don't contribute significantly. On the other hand, the full R(E, K, M ) does for example contain the nonperturbative corrections associated to the genus expansion seeded by the three-holed sphere, which are significant. 42

A.4 Fixing eigenvalues near the spectral edge
Close to the spectral edge |E| 1, JT gravity reduces to topological gravity or the Airy model with spectral density: 43 This theory is identical to the (2, 1) minimal string. The (p, 1) minimal strings are topological, and for these models the multi-brane correlators can be calculated exactly for any 42 We could find more perturbative contributions by including the exponentials of these surfaces in the brane correlators such as (91). 43 We have rescaled the energy, removing the e S0 dependence here.
value of the string coupling. 44 This is the content of formula (1.11) in [49]. In the case of the (2, 1) minimal string, we have: Multi-brane correlators are then calculated as formula (1.11) in [49]: 45 It is again straightforward, but slightly tedious to calculate the multi-brane-pair correlators that get the Airy cluster functions. We we'll show how this goes for R(E) and R(E 1 , E 2 ), and investigate the Airy spectral density with one fixed eigenvalue ρ(E) λ .
For the two-brane correlator, we have: Setting x 1 → x 2 , one finds: Inserting the solution (104), and using the Airy equation Ai (x) = xAi(x), this becomes: This is proportional to the Airy spectral density: 46 44 It is solvable because we can solve the two coupled differential equations that define the single brane correlator. This function is known as a Baker-Akhiezer function of the KP hierarchy. For more on that see for example [49] or the lecture notes [61]. 45 We have ψ(x, τ ) = Ai(x + τ ), therefore ∂ τ = ∂ x , which is one of the two differential equations that defines the Baker-Akhiezer function. The other one is the Airy equation. By rescaling the energy axis we can eliminate the τ -dependence. 46 The normalization of the wavefunction (104) is chosen different from that in (93), hence the different proportionality factor.
To calculate the two-brane-pair correlator, we are led to consider (105): The partial derivatives generate a total of 64 terms, of which some cancel, but 24 remain. For example the first term we would write down is: .
Each such term has 6 derivatives to distribute among the Airy functions, with a maximum of 3 per Airy. Taking x 1 → x 2 = x and x 3 → x 4 = y, one ends up with terms as: Each term has now 8 derivatives to distribute among the Airy functions, with a maximal of 4 per Airy. Repeatedly applying the Airy equation, one finds after what is very much a bookkeeping exercise: with K(x, y) the well-known Airy kernel: K(x, y) = Ai (x)Ai(y) − Ai(x)Ai (y) x − y .
This replaces the role of the sine kernel S(E i , E j ) for GUE away from the spectral edge, also in higher cluster functions [25]. Now that we have the appropriate clusters near the spectral edge, we can redo the analysis of section 4 and fix eigenvalues in this region, as formulas (56) and (67) are completely general. For example: with R(E) from (109). For an eigenvalue not too close to the spectral edge, we have: One recognizes the same features as in (60). It is interesting to see what happens when we insert an eigenvalue very close to the spectral edge or even in the forbidden region E < 0. Note that this is a very unlikely situation since the total spectral density in the forbidden region is much less than one. Hence, when inserting an eigenvalue in the forbidden region, we expect a depletion in the continuum of essentially the entire forbidden region and of the region closest to the spectral edge. Armed with our exact Airy formula (114) we find this is indeed the case. For example when fixing an eigenvalue at the origin, one finds:
This formula appeared in [49]. It means that we can add eigenvalues to an ensemble by introducing pairs of branes ψ 2 (λ) and integrating out λ.

C Details on late-time contributions
We check that certain contributions to the two level spectral density (67) are irrelevant for the behavior in the plateau region t > 2πρ(E). Concerning the terms on the second line of (67), taking the Fourier transform, we are led to: The leading factorizing piece decays as a power law: . (122) Using (48), we see that the E dependence of the connected piece is due to terms of the type: The integral over E gets a function f (τ ) which is finite and has compact support. Because of this we can get the late time behavior by Taylor expanding (t − τ ) −1 = 1/t + τ /t 2 + . . . . The leading term is of the imaginary part of: This exists and is finite. The result are contributions with the same type of t-dependence as (122). Concerning the final contribution of (67), we have a leading contribution which gets the averaged ramp (25). This is negligible beyond the plateau time. The other terms have E 1 and E 2 dependence either of either type: These result in corrections of the averaged ramp behavior, but does not contribute significantly beyond the plateau time.