Microstructure in matrix elements

We investigate the simple model of Pennington, Shenker, Stanford and Yang for modeling the density matrix of Hawking radiation, but further include dynamics for EOW branes behind the horizon. This allows interactions that scatter one interior state to another, and also allows EOW loops. At strong coupling, we find that EOW states are no longer random; the ensemble has collapsed, and coupling constants encode the microscopic matrix elements of Hawking radiation. This suggests strong interior dynamics are important for understanding evaporating black holes, without any ensemble average. In this concrete model the density matrix of the radiation deviates from the thermal state, small off-diagonal fluctuations encode equivalences between naively orthogonal states, and bound the entropy from above. For almost evaporated black holes the off-diagonal terms become as large as the diagonal ones, eventually giving a pure state. We also find the unique analytic formula for all Renyi entropies.

: Normalized version of the density matrix (1.2) for k " 20 and e S " 50 , 5 , 1 (left, middle, right). Orange is positive and blue negative. Intensity of colors reflects the magnitude of individual matrix elements. Off-diagonal elements are less/not suppressed for black holes that have almost evaporated (middle/right).
We believe that more realistic models of quantum gravity, like those typically imagined in AdS/CFT, require no ensemble averaging. The price for unitarity, is a less simple bulk description. See [32][33][34] for a concrete recent example that averaging over string theories is not required.
In those realistic models, the Page curve is unitary, and the density matrix of the radiation follows the standard rules of quantum mechanics; therefore the state does not remain mixed. We would like to understand the bulk phenomena that explain the deviations of the density matrix from being maximally mixed. To accomplish this, we study the gravity description of one member of the ensemble dual to JT gravity with non-dynamical EOW branes [6], focusing on the description of density matrix elements.
The partition function for the ensemble dual to JT gravity with non-dynamical EOW branes, is [6]  One member of this ensemble is described by a LˆL Hamiltonian H 0 , describing the bulk gravitational degrees of freedom; and furthermore by a Lˆk matrix C 0 , describing the interior states [6]. In section 3 we explain that in one member of the ensemble the density matrix of the radiation is essentially Here we took the microcanonical ensemble with e S black hole states. We want to understand the gravity dual to a theory with fixed C 0 , which produces this density matrix; this is described in section 2. The numerical structure of these matrix elements is discussed in section 3. The state is plotted for fixed C 0 in Fig. 1. 2 The gravity interpretation of fixing the random Hamiltonians to H 0 is a decoupled problem, aspects of which have been understood in [12,13,35,36]. This is not the focus of the present paper; nevertheless we include a short discussion on the associated extra ingredients in section 4.

Summary, structure and main lessons
In section 2 we are invited by gravity considerations to investigate deformations of the matrix integral 3) The actual model for which we construct the gravitational dual in section 2 is slightly more complicated; here we simplify for presentation purposes. The gravitational interpretation of this deformation is to introduce scattering interactions from one EOW state into another, with coupling constants g ij that depend on γ and C 0 . See section 2. (1.6) The stronger the interactions the less random the matrix C, and the more realistic the quantum gravity model under consideration. This is one key lesson of this work; in these two dimensional models, realistic gravity systems involve strong interior dynamics. See section 2.4.
Furthermore, the microscopic data of the theory, here represented by the non-random matrix C 0 , are encoded in the specific coupling constants g ij for the EOW brane interactions. See section 2.3.
2 Similar plots were made for SYK with fixed couplings [11], using a different representation for the density matrix (1.2). 3 This is technically rather similar to deformations considered recently in [35]. We briefly suppress the H matrix integral. 4 Overall constants are irrelevant. xψ j |ψ i y " δ ij``k ÿ k 1 "1 k 1`k ÿ k 1 ,k 2 "1 k 2 k 1`.
. . For weak coupling γ ! 1, the matrix elements acquire small off-diagonal components (1.7) The first term, δ ij , corresponds to the usual rule for summing over EOW branes connected to asymptotic boundaries. If the brane flavor is the same on both ends of the boundary, the EOW particle can freely propagate. The second term accounts for scattering interaction between branes of different flavor, this leads to nonzero off-diagonal matrix elements [6,11,37,38]. See section 2.2 and Fig. 2.
In the strong coupling regime γ " 1, one recovers the non-random density matrix (1.2). Although it may not be obvious, the off-diagonal terms are small relative to the diagonal ones in (1.2) when e S " 1.
In realistic gravity models, without averaging, the density matrix (1.2) is not maximally mixed. The matrix C 0 has dimensions e Sˆk , which means that the rank of the density matrix is upper bounded by both k and e S [10]; this suffices to understand the Page curve. See section 3.1 and section 3.2.
When we are ignorant about the microstructure of our system, we believe the states |ψ i y are linearly independent; but in reality there are equivalence relations between them, there are null states. This is captured by numerical plots in Fig. 1. The leading order approximation to the density matrix retains only the dominant diagonal in Fig. 1 (left) and suggests a maximally mixed state; however the smaller off-diagonal matrix elements become competitive when k 9 e S and encode equivalence relations. When the black hole has almost evaporated e S 9 1, the off-diagonal matrix elements are actually no longer suppressed, and the density matrix is far from being maximally mixed. Ultimately, for e S " 1 the state becomes pure again. See Fig. 1 (right) and section 3.2.
The Page transition is caused by the collective behavior of many small off-diagonal matrix elements of order Ope´S {2 q [6,11,37,38], but the fact that the density matrix becomes pure again actually relies on the off-diagonal matrix elements ultimately becoming large at the end of evaporation.
In section 3.3 we obtain the unique analytic continuation for the Renyi entropies within the planar approximation, as an aside. In section 4 we summarize our main findings, and discuss generalizations.
In section 2.1 we review the simple model of [6] both from the gravity and matrix model perspective. In section 2.2 brane interactions are introduced on the gravitational side, and it is shown that previously orthogonal inner products receive small overlaps due to these new interactions. In section 2.3 the matrix model deformation that gives the brane interactions is introduced and analyzed, and in 2.4 the strong coupling limit is considered where the matrix C is fixed to C 0 . 5

Replica wormholes
Following [6], let |ψ i y be the state of a black hole with a non-dynamical EOW brane with flavor i behind the horizon. We would like to think of this state as modeling an interior mode of the radiation.
Furthermore, let |iy be a basis of an auxiliary system, modeling the early outgoing Hawking modes.
Evaporating black holes can be emulated by considering the entangled, unnormalized state [6] k ÿ i"1 Naively |ψ i y form a k dimensional basis; and concordantly this state is maximally mixed. Black hole evaporation is modeled by increasing k. Naively, increasing k can make the entanglement between the black hole and radiation arbitrarily large, even exceeding the Bekenstein-Hawking entropy S -this is the maximal entanglement the black hole can have with the exterior. This contradiction is a version of the information paradox, similar to the Page curve.
The resolution is that the states |ψ i y are not a k dimensional basis; they have small nonzero overlaps that conspire to put an upper bound e S on the dimension of their span [6,11,39,40]. The way this gets diagnosed in [6], is by computing the Renyi entropies associated with the reduced density matrix of the The state |ψ i y is prepared by shooting in an EOW particle with flavor i and mass µ at the boundary and then evolving for a thermal time β{2, giving the black hole some finite temperature. The bra works similarly, and because the EOW branes are non-dynamical there is no way for them to change flavors along their trajectory. To calculate inner products we apply the rules of the gravitational path integral.
The inner product (2.2) introduces boundary conditions consisting of an asymptotic boundary of length β connected to two EOW branes of flavors i and j. We sum over all geometries consistent with these boundary conditions. This results in the following gravitational amplitude, to leading order in e S 5 The Hamiltonians H remains random throughout this section, see section 4. 6 Throughout this section the density matrix is unnormalized.
The black circle denotes an asymptotic JT gravity boundary of length β, and the orange line denotes the EOW particle with mass µ. 7 These amplitudes significantly simplify in the microcanonical ensemble, the black circle then denotes a microcanonical boundary condition [41]. Details about gravitational amplitudes are gathered in appendix A. 1. 8 This looks like a maximally mixed state with the states |ψ i y spanning a k dimensional basis due to the δ ij . However, this conclusion changes when considering higher moments of matrix elements ρ ij ρ kl " xψ j |ψ i y xψ l |ψ k y . (2.4) The particle with flavor i can end up where one detects the outgoing particle with flavor j, however it could also end up where one detects the flavor l. Therefore we have two contributing Feynman diagrams This is not simply the square of (2.3), because JT gravity with EOW branes is not one quantum system, but an ensemble of unitary quantum systems (1.1).
The expectation value for off-diagonal elements xψ j |ψ i y vanishes, but their standard deviation does not, as captured by the second geometry in (2.5). These smaller off-diagonal terms become important and severely constrain the span of |ψ i y when k exceeds e S .
Consider for example the purity, which in the leading order approximation becomes [7] 9,10 The second term stems from the replica wormhole in (2.5) and dominates when k " e S , this effectively places an upper bound e S on the dimension of the span of interior states |ψ i y. The ensemble averaged theory knows about the finite dimensionality of the interior state space.
The ensemble in question consists of a random kˆL matrix C, describing the EOW brane degrees of 7 We exclude brane labels on most gravitational diagrams for presentation purposes. 8 We suppress subleading higher genus wormholes for reader comfort throughout. 9 There are corrections in the denominators from contracting with the normalization of the density matrix [42]. 10 We define Rn " Trpρ n q throughout, with ρ normalized. These are not quite the Renyi entropies but play the same role.
freedom; and some random LˆL Hamiltonian H describing the bulk gravitational degrees of freedom (1.1). In this ensemble language, the random interior states and corresponding matrix elements are [6] 11 The vectors |ay represent a fixed rigid basis, analogous to the spin basis in SYK [43], the Hamiltonian H in this basis is a matrix of random numbers. The brane states |ψ i y are some random linear combination of the fixed basis states, specified by the matrix C. Wick contractions of C are EOW branes in gravity, see section 2.3.

Interior dynamics
In a unitary gravity theory the density matrix elements are just numbers without any variance. This raises the question of what the gravity interpretation of these numbers is. We partially address this, by describing a gravity model where the ensemble over random matrices C collapses to a fixed matrix C 0 .
We now enrich this model by allowing EOW interactions. EOW branes are boundaries on which two dimensional spacetimes end. This severely limits the set of EOW brane dynamics that we can introduce: 1. There can be interaction vertices for 1 Ñ 1 EOW particle scattering where a particle of flavor i scatters to a particle of flavor j, potentially accompanied by the emission or absorption of particles into the bulk spacetime. We restrict ourselves to one type of interaction, weighted with coupling Our specific choice for interaction vertices, which is of the first kind above, is detailed below and in appendix A.2.
2. An EOW particle can propagate between two distinct points, these are either interaction vertices or points where the EOW particle ends on a bra or ket. We have the liberty to include an extra 11 The˘means we multiply both signed gamma functions. We implicitly use the double scaled Hamiltonian throughout.
overall factor G weighting every EOW propagator G , (2.10) which corresponds to adding a constant to the EOW brane action (2.8). This is somewhat ad-hoc, but entirely similar to the introduction of S B in [10]. One can imagine more complicated quantum systems living on EOW branes which give these extra factors. We forsake the details since we are only interested in constructing a proxy for more general dynamical interiors.
3. When 1 Ñ 1 interactions are included, there exist clearly also 0 Ñ 2 creation events, and 2 Ñ 0 annihilation events with the same coupling constants g ij . Furthermore there is nothing preventing dynamical EOW branes from forming closed vacuum loops, making holes in the spacetimes. This exhausts all options for EOW brane interactions. 12 Let us summarize the ingredients of our theory. There are k flavors of EOW particles with the same mass µ, and we have black hole states |ψ i y for each particle flavor. There are interaction vertices where a particle of flavor i turns into a particle with flavor j. Because this theory is dynamical, any number of interactions is allowed, and we must sum over all possible interactions when calculating amplitudes.
Also, there are closed loops of EOW particles, with and without [47] interaction vertices on them.
The main new ingredient are interaction vertices where branes can change flavor; we must choose a specific way to model these interactions and deduce the corresponding JT gravity boundary conditions.
Then each flavor of EOW particles corresponds with a D-brane, and interaction vertices correspond naturally with insertions of boundary operators T n ij ; these are open string Tachyons stretching between D-branes, with Chan-Paton indices i and j. Only the simplest of the chiral vertex operators T 1 ij , known as marking operators, have a known JT gravity interpretation [51,[54][55][56] when stretching between Dbranes with FZZT boundary conditions [8,[57][58][59]. 13 Marking operators then correspond with the β " 0 limit of a fixed length boundary. This extends to mass µ boundary particle segments, relevant for EOW particles, since these are linear combinations of FZZT segments. See appendix A.2.
We have obtained sensible JT gravity boundary conditions for the EOW interaction vertices. This allows us to calculate any desired amplitude. We consider two examples to clarify the rules. 12 For example 1 Ñ 2 brane scattering interactions clearly cannot represent the boundary of some nonsingular spacetime. 13 Recently an educated guess was made for the interpretation of the other boundary operators [54].
We first consider something we call the D-brane partition function, consisting of all interacting EOW brane loops, and all spacetimes ending on them. These are the vacuum fluctuations of our system, they are modded out in all interesting calculations of matrix elements. Going through this first, lightens the presentation of matrix elements later.
The D-brane partition function Z has the following expansion 14 This is a sum over the number n of scattering interactions, within each closed EOW particle loop. The 1{n symmetry factor is because cyclic permutations of the flavors describe the same Feynman diagram, which should be counted only once. The log reflects the fact that in the D-brane partition function Z, we can have any number m n of those closed EOW particle loops with n interactions and an identically ordered set of flavors; which are therefore indistinguishable.
To compute the full D-brane partition function we must exponentiate the disks (2.11) and then fill in bulk geometries; this includes cylinders connecting disks and more general wormhole topologies.
It is no accident that we use the same notation Z as for the matrix integral partition function, these are the same modulo disconnected spacetimes, like the sphere, which can be ignored; see section 2.3.
To compute these diagrams we must first isolate the EOW brane Feynman rules, meaning the factors G and g ij , from the basic JT gravity amplitudes; then simply compute the latter.
We find it convenient to denote gravitational boundary conditions by their matrix integral counterparts. The JT gravity calculation is insensitive to the flavors of the EOW particles, and only depends on the number of interactions n. In that case there are n segments of mass µ boundary particles, separated by marked points; as explained in appendix A.2 this corresponds with the operator insertion Tr´Γ´µ´1{2˘iH 1{2¯n¯. (2.12) The double scaled matrix integral dual of the loop without marking operators was deduced in [47] and Tr log´Γ´µ`1{2˘iH 1{2¯¯. (2.13) When translating the partition function (2.11) to gravity calculations, the EOW particle Feynman weights come out as prefactors for the gravity amplitudes. Because of the summation over flavor indices, the couplings g ij combine nicely into traces, and we obtain the gravitational "boundary conditions" log Z " k G Tr log´Γ´µ`1{2˘iH 1{2¯¯`8 ÿ n"1 G n n Tr´g n¯T r´Γ´µ´1{2˘iH 1{2¯n¯. (2.14) There is an exponential of JT gravity boundaries in Z, represented by the traces, we should sum over all spacetimes ending on them. The way to proceed with the gravity calculations is to use the general identity for expectation values of observables in any theory, crucial to understand D-branes [8,60] log xexppxqy " To good approximation one can then only include the exponential of disk shaped topologies, and annulus shaped topologies connecting to EOW particle boundaries [8]. 16 Because we will not need any detailed answers for the point we are trying to make in this paper, we omit the resulting expression.

Matrix elements
Next consider matrix elements ρ ij " xψ j |ψ i y, which give the gravitational boundary conditions discussed in [6] xψ j |ψ i y " . (2.16) As always we should sum over all possible Feynman diagrams ending on the boundary conditions. This includes EOW particle dynamics, and all gravitational spacetimes ending on the resulting diagrams.
Let us first consider the leading order amplitudes in small g ij perturbation theory, ignoring vacuum loops of EOW particles. For diagonal matrix elements one obtains up to order g ij where in the second equality we applied the EOW particle Feynman rules, and rewrote the gravity amplitudes by the corresponding observables in random matrix theory. For off-diagonal matrix elements, 16 The others topologies have negative Euler character and therefore contribute negligibly assuming that e S 0 " 1.

(2.18)
So off-diagonal elements in the interacting theory are nonzero, unlike in the non-interacting model (2.3).
This is the first sign that EOW particle interactions are important for understanding matrix elements in any non-random gravitational theory.
Higher orders in g ij are obvious; ignoring the vacuum loops this is an expansion in the number of scattering interactions on the EOW brane Note that unlike in (2.11) there is no 1{n for the diagram with n interactions, because the bra and the ket break the cyclic permutation symmetry. We translate this to pure gravity amplitudes by extracting the EOW particle Feynman weights. The JT gravity amplitudes are insensitive to the flavors, the sum over intermediate flavor indices combines the couplings into traces

2.20)
We should also include the effects of the EOW particle loops Z. This means we will have contributions to (2.19) with extra holes in the spacetimes; these holes are the EOW particle loops that the spacetime wormholes are connecting to, for example 17 There can be any number such holes in each portion of spacetime; if some of those have the same labels they are again indistinguishable. Notably one should not include the possibility of EOW particle loops that connect via spacetimes to each other, but not to any probe boundaries; those are normalized away with the 1{ xZy.
It is straightforward to compute the matrix elements order by order in small g ij perturbation theory; these are JT gravity amplitudes, which are known exactly. One all-encompassing example is discussed The generalization to products of matrix elements like ρ ij ρ kl " xψ j |ψ i y xψ l |ψ k y, relevant for Renyi entropies, is clear. There are scattering interactions on all the EOW branes, and a nonzero answers for all values of i, j, k and l. Furthermore, there can be holes due to EOW loops in all pieces of spacetime.
The Lorentzian interpretation is that there are interactions in the interior, where these EOW branes reside [6,43,47]; schematically the associated Lorentzian spacetimes are

Dual matrix integral
We next discuss the matrix integral dual of this model of JT gravity with interaction EOW particles.
Using this, we can consider large g ij . This strong coupling limit selects one member C 0 of the ensemble.
Consider first the matrix dual to JT gravity with non-dynamical EOW branes with partition function (1.1), and with matrix elements (2.7) Z " ż dC dC : expˆ´Tr´C : C¯˙ż dH expˆ´L Tr´V pHq¯˙(2.23) The Gaussian integral over the complex matrix C reduces to standard Wick contractions For one matrix elements ρ ij " xψ j |ψ i y, the ensemble average over C therefore gives For two copies of the matrix element ρ ij ρ kl " xψ j |ψ i y xψ l |ψ k y, summing over Wick contractions gives 18 As explained in appendix A.1, these operator insertions in the H matrix integral correspond in gravity with the EOW brane geometries shown in formula (2.3) and (2.5) respectively. Each Γ`µ´1{2˘iH 1{2ȓ epresents a geodesic boundary segment with a mass µ EOW particle, and each factor e´β H corresponds with a fixed length β segment. Segments inside each trace form a closed loop. We see that every Wick contraction of matrix elements of C becomes an EOW particle propagator in gravity. Therefore, (2.23) and (2.24) corresponds indeed with JT gravity with non-dynamical EOW branes, see appendix D in [6].
We claim that the model of section 2.2 corresponds with the deformed matrix integral where the coupling constants for EOW particle scattering are determined by g and C 0 For now the matrix C 0 are just complex numbers parameterizing the coupling constants, but eventually it will represent the matrix that the C ensemble collapses to, see section 2.4. We next prove that (2.29) is equivalent to the gravity model of section 2.2, by computing the same quantities and matching them.

Brane partition function
We start with the brane partition function Z. The integral over C in (2.29) remains a simple Gaussian, this can immediately be computed by completing the square 19 19 We discard some irrelevant overall normalization constant which depends only on G.
C 0 , C : 0 . Concordantly, the gravitational interpretation of this insertion is not immediately clear; the dictionary between operator insertions in random matrix theory and gravitational boundary conditions concerns only U pLq invariant operators. To transform the above into U pLq invariant operator insertions we can diagonalize the random matrix H with random unitaries U [61] 20 with dU the Haar measure on U pLq; and then explicitly compute the integral over the random unitaries.
We are led in (2.31) to calculate the following integral over Haar random unitaries This Harish-Chandra integral can be computed exactly [62,63]. However, to link with the expansion of section 2.2, it is actually more practical to instead apply (2.15) to correlators in the Haar ensemble [35] log The unitary group integrals on the right evaluate to a double sum over permutations 21 This sum over permutations is weighted with Weingarten functions Wgpσ, Lq, which are known explicitly [64,65]. We are interested in continuum JT gravity, where L " 8. Using the leading large L behavior of Weingarten functions; one checks that these correlators of Haar unitaries reduce, for large L, to the Wick contractions σ " τ generated by a Gaussian complex matrix U with variance L [13]. Explicitly In this Gaussian approximation one therefore finds The following manipulations follow those previously used in [35], where more details can be found. 21 Trσ´A n¯" ź σ i The combinatorial prefactor counts the cycles of length n in S n , only these contribute to the connected correlator. Combining this with (2.34) and using the dictionary between the couplings (2.30), the brane partition function (2.31) becomes This matches the D-brane partition function of our JT gravity theory with interacting EOW particles (2.14), modulo the first term in (2.14), which represents EOW loops without interactions. Those have nothing to do with the C matrix integral, and are therefore of little interest here. We can include them by deforming the potential V pHq by the first term in (2.14), before doing the C or U integrals in (2.29).
In summary, we have shown that the deformed matrix integral partition function (2.29) is equivalent for weak coupling γ, with the D-brane partition function of JT gravity with interacting EOW branes discussed in section 2.2. At strong coupling there are modification to this picture, as the approximation (2.38) breaks down; we discuss in the discussion section 4 how this affect the gravitational description.
We define the strongly coupled version of the gravity theory in section 2.2 as the gravity dual to (2.29).

Matrix elements
We now check that the matrix model (2.29) also reproduces the matrix elements (2.20) that we computed in the gravitational theory, hence establishing the full fledged duality of the matrix and gravity models.
These matrix elements are computed in the deformed matrix integral (2.29) via the dictionary (2.24) We first calculate the Gaussian integral over the matrix C. The two matrices C and C : in the operator insertion can either contract with one another, or with the C and C : in the deformation term in (2.29).
The first Wick contraction gives where Z represents (2.31). The Wick contraction with the exponential deformation gives Neither of these integrands is U pLq invariant, to give this integral a gravitational interpretation we again diagonalize the Hamiltonian as in (2.32), and compute the resulting integral over random unitaries U via the Gaussian large L approximation as explained around (2.38). For the first contribution (2.41) this is the same calculation as (2.33) because the operator insertion on the second line of (2.41) is U pLq invariant, therefore we find with Z representing (2.39). The second contribution (2.42) requires the large L Gaussian approximation for the integral We need to account for Wick-contractions between the operator insertion and the exponential, in other words this calculation can be rewritten as (2.46) Combining everything one finds that the second contribution (2.42) becomes Combining with (2.43), we see that the matrix elements are indeed computed via the gravitational rules discussed in section 2.2; since we precisely recover (2.20). The exponential within this H integral, and the corresponding normalization with 1{Z reflects the fact that we also include geometries where EOW particle loops are connected with the probe boundaries; as discussed below (2.20). These calculations generalize in an obvious way to products of matrix elements like xψ j |ψ i y xψ l |ψ k y.
In summary, we have shown that all observables in the deformed matrix integral (2.29) are equivalent for weak coupling γ, with those of JT gravity with interacting EOW branes; establishing their duality. and a gravity model with the matrix C fixed to one member of the ensemble at strong coupling γ " 8 Z "

Strong coupling and non-random states
ż dC dC : expˆ´Tr´C : C¯˙δ´C´F pHq 1{2 C 0¯δ´C :´C: The stronger the interactions, the less random the matrix C, and the more realistic the quantum gravity model under consideration. This is one key lesson of this work, in these two dimensional models, realistic gravity systems appear to feature strong interactions/have strongly interacting interiors g ij " 1.
Equally important, the microscopic details of the theory, here represented by the non-random matrix C 0 , are encoded in the specific coupling constants g ij for the interior mode interactions, via (2.30). These 22 In this section we are only interested in the C matrix integral and suppress the H ensemble for presentation purposes.
coupling constants take typical values respecting the fact that C 0 is a typical draw of the ensemble.
Concordantly, one can also interpret the completely random theory (2.49) as an ensemble of gravity theories, with the ensemble average over the specific coupling constants of the theory. The interpretation of (2.49) as bulk models with random couplings constants agrees with ideas of Coleman and company [10,24,[66][67][68]. Here too, the theory with random couplings is "simpler" than that with fixed couplings.
In the most realistic models, where γ " 8, the matrix elements of Hawking radiation (2.24) become xψ j |ψ i y "´C : This is obvious from (2.50), but can also be seen explicitly in (2.41)  This extends to products of matrix elements, EOW branes without any interaction vertices on them are suppressed because G " 0; for the product of two matrix elements one obtains for example xψ j |ψ i y xψ l |ψ k y "´C : where H remains random but C 0 is non-random. The non-randomness of C 0 does not imply that (2.52) is numerically close to factorizing, or to the matrix elements in a completely non-random gravity theory.
Since H remains random there remains large correlation between two density matrix elements. This is largely because of the random unitaries U ; not so much the more common eigenvalue correlation. 23 The bulk gravity description of the matrix elements of Hawking radiation in realistic incarnations of these two dimensional quantum gravity models, described by non-random matrices H 0 and C 0 , involves at minimum these interior mode interactions.
On top of that, these involve extra ingredients that capture the non-random matrix H 0 . What these ingredients are is an orthogonal question that goes beyond our current scope; progress has been made in this regard in [12,13,35,36], see section 4. Regardless of the specific details of those extra ingredients, the result is some gravity theory whose matrix elements are good-old-fashioned non-random numbers xψ j |ψ i y "´C : The structure in these numbers for a typical draw of the ensemble deserves some attention. 23 These random unitaries did not play any role for the non-dynamical theory (2.29), because its action is U pLq invariant. However this U pLq invariance is broken by the deformation (2.48) because C0 is some fixed matrix not proportional to the identity.
In this section we examine the consequences of collapsing the theory to a single member of the ensemble, by studying numerical features of the matrix elements (2.53), and how they reproduce the Page curve.
Throughout this section, we work in a microcanonical ensemble centered around E, with e S states.
Since F pH 0 q is approximately constant for eigenvalues of H 0 within the microcanonical window we have The kernel F pEq drops out of the normalized density matrix and all derivative quantities like the Renyi entropies and the von Neumann entropy, this means we can immediately drop it and continue with Crucially, the matrix C 0 must be interpreted as typical representative of the undeformed C ensemble (2.49). This means the entries of C 0 are complex numbers with typical norm-squared one. 24 Formula (3.2) is a toy model for the matrix elements of evaporating black holes in quantum gravity.
The ensemble is an incredibly powerful tool that simplifies calculations and presents a simpler effective picture, but is has also caused confusion. Mistaking the ensemble for proper quantum mechanics, one apparently finds that the density matrix is maximally mixed, as in Hawking's calculation [1].
In realistic theories, represented by the toy model (3.2), the density matrix is not maximally mixed.
The complex matrix C 0 has dimensions e Sˆk , and concordantly the rank of the matrix is upper bound by both k and e S [10]; this suffices to understand the Page curve. The point is that when we are ignorant about the microstructure of our system, that we believe the states |ψ i y are linearly independent; where in reality there are equivalence relations between them, there are null states. The remaining question is then how microstructure gets realized in the bulk, this work has been a step in that direction. 25 We plotted the density matrix in Fig. 4, the purity in Fig. 5, and the von Neumann entropy in Fig.   6 for one single matrix C 0 . Comparison with the averaged answers confirms the latter are self-averaging.
In the remainder of this section we explain analytically why (3.2) produces these plots. 24 Fixing to non-typical members is physically non-sensible, since those systems are not even remotely accurately described by the ensemble to begin with. 25 Factorization is another interesting question, this remains geometrically nontrivial even in microscopic realizations; we have not investigated that here [6, 8, 11-13, 30, 36, 69, 70].
Repρ ij q

Density matrix
The unnormalized density matrix of radiation is plotted in Fig. 1 and Fig. 4, and reads explicitly To understand this figure, let us estimate the typical size of the matrix elements. The entries of C 0 are complex numbers with typical norm squared one, which means that the diagonal matrix elements are real numbers with typical size e S , because they are the sum of e S real numbers of typical size one.
The off-diagonal matrix elements are complex numbers with a "random" phase, and typical amplitude of size e S{2 . This is because they are the sum of e S complex numbers with random phases and typical amplitudes of size 1, combined they are to be interpreted as a "random" walk in the complex plane with unit step length and e S steps 26 .
The typical radial distance traveled by this random walk is the square root of the number of steps e S{2 .
The off-diagonal terms correspond to sums over random complex numbers that do not constructively 26 Random is quoted because there is nothing uncertain about the number in (3.5), "typical" might be more appropriate.
interfere, the phases in (3.4) align but those in (3.5) do not, hence they give small contributions relative to the diagonal terms as long as e S " 1.
Alternatively one could directly estimate the norm squared for the off-diagonal matrix elements The most sizable contribution comes from the terms with α 1 " α 2 . These are e S real numbers of typical size one, adding up to give some real contribution of typical size e S . There are other terms that pair up to form real numbers, combining terms where we exchange α 1 and α 2 . The typical size of the sum of such two terms is ? 2; but there is typically the same number of terms with overall positive sign, and overall negative sign. Therefore these are subleading. Crucially, the leading contributions always come from terms where the phases precisely cancel.
Typical values of matrix elements are most efficiently estimated by using the ensemble average, this is what typicality means. In the ensemble averaged description, one indeed computes xψ i |ψ j y aver " δ ij e S , |xψ j |ψ i y| 2 aver " δ ij e 2S`eS , (3.7) reproducing the above typical estimates. These formulas also highlight the discussion below (3.2). The leading order approximation for the density matrix is the maximally mixed diagonal Hawking state; but there are subleading non-self-averaging contributions of order e S{2 in the unnormalized density matrix ρ " e S k ÿ i"1 |iy xi|`Ope S{2 q . (3.8) We cannot emphasize enough that these subleading corrections can, and do save unitary [6,11,37,38].
We will show in the following section that these subleading corrections are important for observables.
Their effect is to cause linear relations between the different vectors |ψ i y The off-diagonal matrix elements do not need to be leading order for such linear relations to exist. The density matrix (3.3) shown in Fig. 4 and Fig. 1, is an extremely concrete example of how this happens.
The leading order approximation is diagonal, but by construction the dimension of the span of the |ψ i y, which equals the rank of C 0 ; is upper bound by both k and e S . The off-diagonal matrix elements are therefore responsible for the Page transition at late times.
Circling back to the gravitational picture of section 2.2, notice in (2.30) that the typical off-diagonal couplings g ij are subleading compared to the diagonal couplings g ii . This intuitively explains why the diagonal matrix elements remain bigger than the off-diagonal ones, even in the interacting theory.

Higher moments
Next we examine how calculations of various entropies in a typical member of the ensemble are consistent with the answers given by the replica wormhole calculations of [6]. We will find that off-diagonal matrix elements constructively add up and give large contributions to observables, corresponding to the effects of replica wormhole geometries.
The simplest entropy observable is the purity, which is plotted in Fig. 5 and reads explicitly here Trpρq is computed using (3.3). As explained below (3.6) the leading contributions come from terms in the sums where two phases align, meaning that complex matrix elements of C 0 pair up as In the purity summation (3.10), this happens if i " j and or α 1 " α 2 . Using the leading order estimate Trpρq " k e S one then estimates the typical purity 27 The first contribution in (3.12) comes from terms with i " j, corresponding with diagonal matrix 27 These orange "Wick contractions" represent the ways that complex matrix elements can pair up in real numbers (3.11).
elements or the disk geometries in (2.6). By itself these diagonal terms claim the purity of the radiation decreases forever as 1{k. However, there is another important contribution from the terms with α 1 " α 2 .
This term counts the norm squared of the non-self-averaging fluctuations in all of the matrix elements ρ ij , represented by the final contribution in (3.7). The noise in the off-diagonal elements constructively interferes when computing the purity, this corresponds with the replica wormhole in (2.6) [6,11].
There are k 2 matrix elements with fluctuations, whereas there are only k diagonal matrix elements whose self-averaging behavior gives the Hawking answer. For late enough times k the non-self-averaging fluctuations win over the self-averaging terms and cause the Page transition [6,11,37,38]. Here this is represented by the second contribution in (3.12) winning over the first one [6], see also It is hence not true that small corrections save unitarity [37,38], large off-diagonal matrix elements do.
The last term in (3.12) estimates the size of non-self-averaging fluctuations, it comes from terms in the sum where a 1 ‰ α 2 and i ‰ j. Each term in this sum is a "random" complex number with typical norm squared one. There are roughly k 2 e 2S such terms, therefore the sum represents a "random" walk that travels a typical radial distance k e S . For large k and e S these fluctuations are small, meaning that the purity is self-averaging. Ensemble averages accurately compute entropies, but not matrix elements.
This generalizes to the other moments R n " Trpρ n q{ Trpρq n  For n ą 2 there are more ways the matrix elements C 0 can pair with their complex conjugates (3.11), corresponding with the sum over Wick contractions in the ensemble averaged calculation, and with the different replica wormhole geometries in the model of [6]. For example for n " 3, 4, one estimates (3.14) Here we are no longer tracking the non-self-averaging fluctuations and we dropped terms of subleading order in either k or e S . This corresponds with using the planar approximation in replica wormholes [6].

Entropy and planar resummation
Using S "´Trpρ log ρq . (3.15) This entropy is plotted in Fig. 6. We now reproduce this figure using the planar approximation for the moments discussed above, by explicitly applying the replica trick S "´B n R n n"1 .

(3.16)
This is nontrivial as it requires finding the unique analytic expression for the Renyis as function of n.
The first difficulty is computing R n for arbitrary positive integer n by summing all planar diagrams.
This is in principle a difficult counting problem, the key to solving this implicitly was discussed in [6], following [71,72]. To obtain the Renyi entropies as function of n one should instead compute a generating function, for example the resolvent of the density matrix Rpλq " Trˆ1 λ´ρ˙" k λ`8 ÿ n"1 1 λ n`1 Trpρ n q . (3.17) Its Taylor expansion around infinite λ encodes the moments R n , and hence also the Renyi entropies.
This expansion, along with the structure of the planar geometries that contribute, makes it possible to write down a Schwinger-Dyson equation for Rpλq; in our microcanonical setup this becomes [6] Rpλq 2`ˆe S´k λ´k e S˙R pλq`k 2 e S λ " 0. (3.18) According to (3.17) the solution must behave as k{λ near λ " 8; furthermore recognizing the generating functional of Gegenbauer polynomials one then obtains the unique solution Using the relation between Gegenbauer polynomials and Jacobi polynomials, 28 we can rewrite this into From the expansion coefficients one therefore finds the moments 29 where changing coordinates as p " n`1´s highlights the symmetry under exchange of k and e S . This reproduces the correct answer for any integer n. One easily checks the simplest cases (3.12) and (3.14).
The second, we believe generally less appreciated difficulty, is finding a unique analytic continuation of this formula away from integer n. Usually, one hopes that one obvious analytic continuation presents itself; however here there are two obvious and inequivalent options. As 1{Γpn`1´sq " 0 for s ą n`1 one can extend the range of the first sum from s " 0 to s " 8. Via the same argument one may extend the second sum from p " 0 to p " 8, for positive integers n. 30 Both procedures give a hypergeometric function, 31 these agree on the positive integers but disagree elsewhere; resulting in two possible analytic continuations (as function of n) R n " 1 k pn´1q 2 F 1 p´n, 1´n, 2, k{e S q or R n " 1 e pn´1qS 2 F 1 p´n, 1´n; 2; e S {kq , (3.24) which are swapped when exchanging k and e S . Then which of these is the correct analytic continuation?
We need a theorem that specifies uniqueness of analytic continuation, given data at the positive integers. 29 This formula was also recently derived in [73,74]. 30 This corresponds with extending the first sum from s "´8 to s " n`1, doing both does not yield a convergent sum. 31 Hypergeometric functions are defined as semi-infinite sums, tread carefully with Mathematica here.
The only such theorem that we know of is due to Carlson, see [75]. If there is a function f pzq that is analytic for Repzq ě 0, that takes assigned values f n on the positive integers, that grows exponentially slower than sinpπzq for imaginary z, and no faster than exponential elsewhere; then this function f pzq is the unique one with these properties.
Conversely, the data f n is insufficient to uniquely specify an analytic function f pzq without further constraints, and Carlson's theorem gives the necessary constraints that uniquely select one function. It is not a priori obvious why the moments Rpzq should satisfy these constraints, but Carlson's theorem proves that if they do not, then there is no unique function Rpzq; and therefore no unique von-Neumann entropy. That is clearly nonphysical, therefore we conclude that Rpzq must satisfy Carlson's theorem. 32 The first function in (3.24) satisfies this theorem as function of n when k ă e S , however it grows too quickly on the negative imaginary axis when k ą e S ; the second function in (3.24) satisfies the theorem for the same reason only when k ą e S . Therefore the unique analytic continuation in n is 33 R n " 1 k n´1 2 F 1 p´n, 1´n, 2, k{e S q θpk ă e S q`1 e pn´1qS 2 F 1 p´n, 1´n; 2; e S {kqθpk ą e S q . S " plogpkq´k{2e S q θpk ă e S q`pS´e S {2kq θpk ą e S q . This agrees excellently with a direct calculation of the entropy (3.15) using our density matrix (3.3), see Fig. 6. It is not surprising that the Page curve is self-averaging, the point is that we have produced it using the density matrices in Fig. 1, confirming the claims made about off-diagonal matrix elements in section 3.2. More importantly, we have given a gravitational interpretation for these matrix elements in section 2.2.

Concluding remarks
In this work we made progress towards understanding the bulk gravity dual to one quantum system. We investigated how the density matrix elements of evaporating black holes are computed in non-random gravity theories, and in particular what explain small off-diagonal density matrix components.
For this we investigated an enrichment of the model of Pennington, Shenker, Stanford and Yang, by including dynamics for EOW branes; namely brane flavor changing interaction vertices and loops of 32 Otherwise we can just add the function c sinpπzq to the moments with arbitrary c, which contributes πc to the entropy. 33 The Heavisides also conveniently save us from evaluating the hypergeometric on the branchcut it has in the last variable. 34 Analytically, using the representation of the hypergeometric as infinite sums (3.23) from s " 0 to s " 8, one takes the derivatives of the Gamma functions, then uses the poles and residues of the Gamma and Digamma functions to prove that only one term in the sum contributes. It is crucial that in the analytic continuation the sums are semi-infinite, otherwise the derivative is not well defined. EOW branes; and discovered a dual description as a deformed matrix integral Z " ż dC dC : expˆ´Tr´C : C¯´γ Tr´´C´F pHq 1{2 C 0¯´C :´C: where the coupling constants for the interaction vertices are related to the matrix deformation as For increasing values of the coupling constant γ, and hence g ij , the random matrix C gets gradually fixed to a non-random matrix C 0 . This means that the interior states |ψ i y are becoming less random.
Our main conclusions are the following.
1. Strong interactions behind the horizon are essential for understanding the microstructure of matrix elements of evaporating black holes, within this simple model.
2. The microscopic details of the theory, here represented by the non-random matrix C 0 , are encoded in gravity as coupling constants for interior interactions.
3. Large amounts of tiny off-diagonal matrix elements eventually overtake the bigger diagonal matrix elements, these correspond with replica wormholes and cause the Page curve transition.
4. For nearly evaporated black holes off-diagonal matrix elements are large, and the state approaches some pure state as required for unitarity.
In a gravity model where the density matrix of Hawking radiation is described by (2.51) one immediately sees that there are nontrivial off-diagonal matrix elements, without having to compute their variance. The raison d'être for the simplified ensemble averaged gravity theories is they are simple to compute with, this is the whole philosophy behind random matrix theory [61,76]. They were never meant to describe the microstructure of individual systems, we should not forget this. Random matrices are sufficiently smart to understand that the off-diagonal matrix elements are nonzero. But they are not a microscopic description of the theory, where we can actually understand why they are nonzero. The real universe is clearly not an ensemble average; no one would claim Navier-Stokes is the fundamental description of fluids, neither are random matrices the fundamental description of the bulk.
Simple effective description like JT gravity with non-dynamical EOW branes, Brownian motion, and pure Einstein-Hilbert gravity in higher dimensions [18,69] are ensembles. However, real fundamental description like deformed JT gravity [35] with dynamical EOW particles, atoms, and full-fledged string theory [33] are factorizing and unitary quantum systems without ensembles.
xψ j |ψ i y " δ ijF igure 7: Picture for (two-sided) matrix elements that generalizes to higher dimensions. One could imagine creating orthogonal states by preparing states with particles of different flavors. In microscopic models the particles could interact heavily in the interior, perhaps close to the singularity, resulting in off-diagonal matrix elements. This conclusion is similar to the effective half-wormholes for eigenvalue correlation in [30,36] and also to the picture of [12,13,35]. Strong interior interactions encode microstructure. This also applies in [5].
That being said, ensemble averaged descriptions are clearly extremely useful, precisely because they corresponds with simple gravitational duals; those simple duals suffice for many calculations.
We end this work with several comments, first and foremost about higher dimensional implications.

General lessons
We believe our findings are evidence that strong interactions in the interior will generically be important to capture the microstructure of higher dimensional black holes. These interactions factorize (replica) wormholes, because they collapse the ensemble, and because without the ensemble everything factorizes.
Knowing how to calculate the microscopic out-state of the radiation is equivalent to understanding how to factorize (replica) wormholes.
In our model factorization is less geometrically obvious than is the case with eigenvalue correlation [12,13], where there is some exclusion rule and concordantly a diagonal " cylinder identity [30]. In this setup, when we calculate in gravity the product of two matrix elements, there is the replica wormhole; but also other connected components, from both matrix elements connecting to the D-brane partition function (or EOW loops). They can both be connected to EOW loops via wormholes, or via the nonlocal interactions discussed below. Since the replica wormholes are not related to eigenvalue correlation, we believe that the nonlocal interactions might be the key. Somehow the replica wormholes should then be canceled by nonlocal interactions between different copies, restoring factorization. This must happen, because the ensemble is collapsed, nevertheless it would be interesting to make this more precise.
This picture we obtain here is, perhaps surprisingly, morally related to the one advocated in [30,36], where they discuss an effective description for (eigenvalue) microstructure. The information about that microstructure is located in the interior, perhaps even near the singularity. We have a different setup here, and are describing different aspects of black hole microstructure, namely the out-state of radiation.
Nevertheless the overall lesson is similar: strong interior interactions encode microstructure. This is also one possible interpretation of [12,13,35], which gives a precise description of eigenvalue microstructure.
See Fig. 7 The EOW branes studied in this paper are behind the horizon [43,47]. But one could consider an alternative version of this same model, with negative energy dynamical branes modeling random states, but these are outside of the horizon [47]. Their coupling constants would still encode microstructure, but that information would now be outside. Perhaps this is a sensible toy model for fuzzballs; it would be interesting to connect that literature better with the current developments concerning wormholes and ensembles [77][78][79].
Concerning the dependence of the couplings on C 0 we believe the generalized picture is the following.
Consider a UV complete theory of quantum gravity, like string field theory. These theories are probably rather unique and special, concordantly the couplings between the matter fields in the spectrum would take rather specific values. One could imagine integrating out most fields in the spectrum, leaving some "simplified" model with fewer fields; like dilaton gravity with EOW branes.
However integrating out the matter fields would leave its imprint, it would not leave something nice and simple like ordinary JT gravity, with non-dynamical EOW branes. Rather, one would obtain some highly deformed JT gravity with complicated dilaton gravity interactions, and dynamical EOW branes.
The details about the UV microstructure would get imprinted in all these interactions, see also [35].
Concerning the off-diagonal matrix elements there should obviously be a generalization to arbitrary black holes. States describing the interior of black holes should acquire nontrivial overlaps as the black holes grow old, to ultimately restore unitarity.
The mechanism by which these states become equivalent is an important open question, our results suggest that strong interior interactions are relevant for understanding this phenomenon. 35

Strong coupling
At strong coupling when g ij become big, the gravitational picture of section 2.2 is modified, since in the approximation from (2.36) to (2.38) and (2.45), we assume that terms with n of order L are suppressed.
When the coupling become big, that assumption is no longer valid; and so neither is the approximation. This Gaussian approximation fails because for n 9 L, the contributions from subleading Weingarten functions are not obviously suppressed [35]. For example, we can no longer trust the scaling formula which validated the Gaussian approximation. The combinatoric prefactors may also enhance naively subleading contributions at high order in the coupling constant.
This means that the multi-trace contribution in (2.36) might become relevant at strong coupling.
One would obtain multi-trace terms in the brane partition function (2.39). These are clearly interpreted in gravity as corresponding with multi-local interaction vertices, we could then represent these by also allowing dotted lines connecting dotted vertices and multiple local interaction vertices; which makes the gravitational expansion more involved. Feynman rules for those dotted diagrams contain further information on C 0 since these rules depend on multi-trace combinations of C 0 .
It would be interesting to obtain analytic control over these multi-trace deformations by scaling the couplings in certain specific ways.

Fixed Hamiltonians
Finally, we briefly mention the gravitational interpretation of fixing the random Hamiltonians H to one single Hamiltonian matrix H 0 . This was investigated in [35] using a deformed matrix integral similar to (2.48), but where H is coupled to an external matrix H 0 with coupling constant 1{σ 2 .
Whilst not our focus, the gravitational interpretation of non-random matrix elements (2.53) involves understanding how H 0 gets encoded in gravity in addition to C 0 . Therefore we briefly summarize the results of [35]. The gravity interpretation for the eigenvalues of H 0 only affects the bulk JT gravity spacetime description, not the behavior of EOW branes discussed throughout this work. This is why our discussion in earlier sections decoupled from fixing H 0 .
For weak coupling 1{σ 2 one finds a deformation of the JT gravity action which can be interpreted as inserting many local defects [51,52,[80][81][82]. These are the analogue of the interaction vertices discussed in this work. The associated coupling constants depend on H 0 , in line with point 2 of the main conclusion.
When the coupling increases, nonlocal bulk spacetime interactions become important, for precisely the same Weingarten reasons, giving a nonlocal dilaton gravity action. The analogue of terms with n 9 L becoming important, is that macroscopic operator insertions appear. These tear up the smooth spacetime with large holes [83].
For strong coupling we approach the eigenbrane picture [12,13] with many extra macroscopic holes in spacetime. The boundary conditions on these extra holes [41] encode the eigenvalues of H 0 . However, the theory with infinitely many eigenbranes is not under good control, and something far more drastic probably happens. Signs were found [35,84] of some branched polymer phase of gravity, where smooth spacetime is completely broken. Then the question is what replaces smooth spacetime; what is the true microscopic description of gravity? These works build towards deriving a concrete microscopic picture.
What remains is the more illusive gravitational interpretation of the random unitaries U . These are irrelevant for observables like partition functions, or the spectral form factor; but crucial for correlation functions [9,11,13,85,86] and density matrix elements. This is an important open problem [42].

A Gravitational amplitudes
Here we gather some gravitational details relevant in the main text.

A.1 Pinwheels
First consider the pinwheel geometry Z n pβq of [6], which is a disk geometry with n pieces of asymptotic boundary of length β, separated by n boundary segments that describe the geodesic trajectory of some particle of mass µ; these represent the EOW branes. For example This amplitude was computed using the techniques of [9,87] and gives We note in passing that this formula is also easily derivable in the BF formalism of [88,89], where the mass µ boundary particles are represented by Wilson lines and one recognizes F pEq as the 3j symbols with one trivial representation; because there is nothing on the other side of the particle.
In the matrix integral this pinwheel corresponds with the observable The leading order expectation value of Tr δpE´Hq equals the disk amplitude with fixed energy boundary conditions [8,12] xTr δpE´Hqy " e S 0 4π 2 sinh´2πE 1{2¯" ρpEq , (A.4) which indeed reproduces (A.2). Including handles on the pinwheel replaces the genus zero disk answer where xρpEqy is the exact spectral density in the matrix integral. This can be calculate order per order in the genus expansion using Weil-Peterson volumes, and nonperturbatively using D-branes [8,12].
When there are two pinwheels, we must include spacetime wormholes that connect them. Summing over all genus gives rise to the full spectral correlation xρpE 1 qρpE 2 qy of random matrix theory [8,12,61] xZ n 1 pβqZ n 2 pβqy " ż`8 which, for all intents and purposes, can be approximated as The generalization to multiple pinwheels is obvious [12].
Though these expressions are very explicit, it is more practical to work in a microcanonical ensemble, where things simplify even further. In some microcanonical ensemble centered around E, one computes for example 36 xZ n 1 pEqZ n 2 pEqy " where the total number of eigenvalues in this microcanonical bin computes the microcanonical entropy The second equality in (A.8) follows from the definition of the microcanonical ensemble, the width of the energy bin is much smaller than 1 but much bigger than the typical level spacing 1{ρpEq. The function F pEq varies on energy scales of order 1 and can therefore be approximated as constant within the bin.
Furthermore, on energy scales bigger than 1{ρpEq the sine kernel in (A.7) is essentially indistinguishable from the Dirac delta term; and these two contributions therefore cancel out, and to good approximation Using these results one computes, with the rules explained in section 2.1 ρ ij " δ ij F pEq e S , ρ ij ρ kl " δ ij δ kl F pEq 2 e 2S`δ il δ kj F pEq 2 e S . (A.11) These are the results mentioned in (2.3) and (2.5). Notice that the details of the EOW brane boundary conditions, captured by the kernel F pEq, are just overall normalization constants in these amplitudes; all this dependence drops out when we consider normalized density matrices, and compute normalizes quantities like Renyi entropies. Life in the microcanonical ensemble is simple.

A.2 Modeling interactions
We gather formulas about disk amplitudes with marking operators in minimal strings and in JT gravity, more details are contained in [51,[54][55][56].
First consider a circular FZZT boundary [8,[57][58][59], without marked points, which corresponds in random matrix theory with TrplogpE´Hqq " . (A.12) The random matrix observable includes a sum over genus in gravity, which we suppress for presentation purposes. Now consider the minimal string boundary three point function of marking operators 37 As written the marking operators intertwine between segments with FZZT boundary states respectively E 1 , E 2 and E 3 . 38 This corresponds in random matrix theory with the following observable [51,[54][55][56] Trˆ1 Now consider a disk with thermal boundary length α 1`β1`α2`β2`α3`β3 , which corresponds in random matrix language with Trpexpp´pα 1`β1`α2`β2`α3`β3 qHqq. Laplace transforming some segment of thermal boundary gives a segment of FZZT boundary, this is obviously true from the matrix 37 Boundary chiral vertex operators have three labels [90], the two extra labels denote the boundary states between which they intertwine. In string language these are the Chan-Patton indices for the two D-branes between which the open string operator stretches, generalized to non-coincident D-branes. 38 FZZT boundary conditions are technically a double cover of the energy axis and should be labeled by z with E "´z 2 , there is a unique Liouville primary corresponding with each z; this is not relevant here so we suppress it for reader comfort.
integral formulas. Therefore we have the correspondence The gravity amplitude mimics the pinwheel amplitude considered in appendix A.1 and in [6], but with FZZT boundary conditions instead of mass µ boundary conditions. FZZT and mass µ boundary states are linear combination of each other; one checks that their JT boundary wave functions form complete sets for certain complex ranges of E respectively µ [6,9], hence there is a basis transform between them.
Now we see that taking the thermal length of one of the segments in the pinwheel to zero reproduces amplitudes of the type (A.14). The random matrix dual clarifies that this is indeed the correct limit to take, if you send β 1 , β 2 and β 3 to zero in (A.15), you recover (A.14). This proves that the dilaton gravity interpretation of a marking operator T 1 E 1 E 2 corresponds with a piece of thermal boundary sandwiched between FZZT boundary segments with boundary conditions E 1 and E 2 , where the thermal length of the sandwiched segments is taken to zero.
The generalization to mass µ boundaries is straightforward, since these are just linear combinations of FZZT boundaries. We can consider for example the minimal string boundary three point function 39 and therefore in gravity with the β " 0 limit of the pinwheel diagrams studied in appendix A.1.
Notice that taking the β " 0 limit of the pinwheels, gives a finite answer for the second diagram in (2.11); representing a disk ending on a mass µ particle with a single marking operator inserted. Naively one might have thought that amplitude would vanish, since the boundary is a geodesic, but it does not.
In summary, if we model interactions by insertions of marking operators, we know the corresponding observables in random matrix theory and the corresponding boundary conditions in gravity, and we can compute all amplitudes we want. One could consider modeling interactions by more general minimal open string Tachyons T n µµ , however the random matrix dual of the corresponding boundary correlators 39 In stringy language the k flavors of interior modes are ordinary Chan-Paton indices, because all D-branes coincide at µ.
is not actually known, 40 and concordantly neither are the precise boundary conditions in dilaton gravity.
We expect the conclusions of this work to hold when working with these other models for interactions.

A.3 Computing amplitudes
Here we go through the JT gravity calculation for one amplitude that contributes to the matrix element.
The example which we choose is sufficiently complex so that the generalization to all amplitudes should be straightforward. We consider (2.21) where we already extracted the EOW particle Feynman rules, such that the diagram reflects a pure JT gravity calculation. The way to proceed is to first treat each boundary loop as analogous to a standard fixed length boundary; chopping up the surface by cutting off "trumpets" ending on each boundary [8], and on the unique geodesic inside the Riemann surface homologous to the boundary in question.
The remaining amputated amplitude with geodesic boundaries computes the Weil-Petersson volume, which can be calculated by further chopping up this Riemann surface into three holed spheres. This is explained in great detail in [8,9,13,25,85,91,92], and will not be repeated here. The newer ingredient is computing the trumpet ending on a boundary circle that includes interacting EOW branes.
Let us work through this in the above example. Cutting the Riemann surface on the blue geodesic of Each of the remaining pieces is a trumpet with one geodesic boundary and one boundary that involves segments of EOW particles.
Including the geodesic boundary is easy within the BF or first order formulation, where it is interpreted as introducing a hyperbolic defect; amplitude wise this simply replaces the sinh factor with a cosine [80] " ż`8 0 dE expp´βEq F pEq 4 1 2π The trumpet with one geodesic boundary and one triangle boundary therefore becomes " 1 pE 1´E2 q 2 " xρpE 1 qρpE 2 qy 0 , (A. 23) which is the genus zero contribution to the spectral correlation [8]. Combining the elements, we obtain Including any number of handles and nonperturbative effects in that genus expansion simply replaces the genus zero connected spectral correlator with the full correlator (A.7) of random matrix theory [9,13,85].
The generalization to arbitrary amplitudes should now be obvious.
Another new application of the BF formulation is the calculation of [47], which considers a trumpet with mass µ particle on the geodesic boundary. In the BF formulation, massive particles become Wilson lines, and the mass µ labels discrete series irreducible representations of SLp2,Rq. The geodesic length b, over which is integrated, labels hyperbolic conjugacy class elements of SLp2,Rq, and Wilson lines in the conjugacy class element basis contribute characters to BF amplitudes [88,89,93,95]. One then finds and the discrete series characters evaluated on hyperbolic conjugacy class elements are in this convention [96] χ µ pbq " expp´µbq 2 sinhpb{2q .