The black hole entropy has been observed to generically turn negative at exponentially low temperatures Tsim eS0 in the extremal BekensteinHawking entropy S0 a seeming pathology often attributed to missing nonperturbative effects. In fact we show that this negativity must happen for any effective theory of quantum gravity with an ensemble description. To do so we identify the usual gravitational entropy as an annealed entropy Sa and prove that this quantity gives S0 at extremality if and only if the groundstate energy is protected by supersymmetry and diverges negatively otherwise. The actual thermodynamicallybehaved quantity is the average or quenched entropy Sq whose calculation is poorly understood in gravity it involves replica wormholes in a regime where the topological expansion breaks down. Using matrix integrals we find new instanton saddles that dominate gravitational correlators at Tsim eS0 and are dual to semiclassical wormholes involving dynamical branes. These brane solutions give the leading contribution to any black hole very near extremality and a duality with matrix ensembles would not make sense without them. In the nonBPS case they are required to make Sq nonnegative and also enhance the negativity of Sa both effects consistent with matrix integrals evaluated exactly. Our instanton results are tested against the onshell action of D3branes dual to multiply wrapped Wilson loops in mathcalN4 superYM and a precise match is found. Our analysis of lowenergy random matrix spectra also explains the origin of spectral gaps in supersymmetric theories not only when there are BPS states at zero energy but also for purely nonBPS supermultiplets. In the former our prediction for the gap in terms of the degeneracy of BPS states agrees with the Rcharge scaling in gapped multiplets of mathcalN2 superJT gravity.
['High Energy Physics - Theory', 'General Relativity and Quantum Cosmology', 'Quantum Physics']
Entropy and Spectrum of NearExtremal Black Holes semiclassical brane solutions to nonperturbative problems
{"Sergio Hern'andez-Cuenca": 'Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \nE-mail:  \nsergiohc@mit.edu', 'Abstract:': 'The black hole entropy has been observed to generically turn negative at exponentially low temperatures T ∼ e -S 0 in the extremal Bekenstein-Hawking entropy S 0 , a seeming pathology often attributed to missing non-perturbative effects. In fact, we show that this negativity must happen for any effective theory of quantum gravity with an ensemble description. To do so, we identify the usual gravitational entropy as an annealed entropy S a , and prove that this quantity gives S 0 at extremality if and only if the ground-state energy is protected by supersymmetry, and diverges negatively otherwise. The actual thermodynamically-behaved quantity is the average or quenched entropy S q , whose calculation is poorly understood in gravity: it involves replica wormholes in a regime where the topological expansion breaks down. Using matrix integrals we find new instanton saddles that dominate gravitational correlators at T ∼ e -S 0 and are dual to semiclassical wormholes involving dynamical branes. These brane solutions give the leading contribution to any black hole very near extremality, and a duality with matrix ensembles would not make sense without them. In the non-BPS case, they are required to make S q non-negative and also enhance the negativity of S a , both effects consistent with matrix integrals evaluated exactly. Our instanton results are tested against the on-shell action of D3-branes dual to multiply wrapped Wilson loops in N = 4 super-YM, and a precise match is found. Our analysis of low-energy random matrix spectra also explains the origin of spectral gaps in supersymmetric theories, not only when there are BPS states at zero energy, but also for purely non-BPS supermultiplets. In the former, our quantitative prediction for the gap in terms of the degeneracy of BPS states agrees with the R-charge scaling in gapped multiplets of N = 2 super-JT gravity.', '1 Introduction': 'The Euclidean gravitational path integral is generally understood to compute a canonical partition function, an identification that provides a powerful window into black hole thermodynamics [1]. Given a black hole at inverse temperature β , the gravitational partition function Z ( β ) can be used to obtain an entropy via the standard thermodynamic identity  \nS ( β ) = (1 -β∂ β ) log Z ( β ) . (1.1)  \nThe insights that this relation and generalizations thereof have yielded into the quantum gravitational physics of black holes cannot be overstated [2-10]. However, that eq. (1.1) actually computes the thermal entropy of a standard quantum system is a statement that can only be established reliably via Euclidean methods in field theory, not in gravity. Hence that S always is and should behave as such an entropy is an assumption which may lead to drawing potentially incorrect physical conclusions. For instance, one may expect S ( β ) to be non-negative for all β , and argue that some contributions may be missing from Z ( β ) if S ( β ) happens to turn negative for some β . In fact, as first realized by [11], quite generally S is not and should not be expected to behave like a thermal entropy in effective quantum gravity.  \nIn this paper, we refer to the quantity S that eq. (1.1) computes as the gravitational entropy. As will be explained, generically S ( β ) must behave drastically different from a thermal entropy at large β , a regime which gravitationally corresponds to near-extremal black holes. Crucially, this means that to reliably study the statistical physics of black holes near extremality one must revisit eq. (1.1) in pursuit of two main goals:  \n- 1. Identify precisely what kind of entropy S is and what its properties are near extremality.\n- 2. Understand how to compute thermal entropies in gravity and, in particular, establish whether semiclassical gravity alone can succeed in describing black holes very near extremality.  \nA brief preview of our findings in addressing these goals is as follows:  \n- 1. This S is an annealed entropy, which provably must diverge negatively in the extremal limit unless supersymmetry protects the ground-state energy. This phenomenon is unrelated to the degeneracy of BPS states, which we quantitatively relate to the size of spectral gaps.\n- 2. Thermal entropies which behave non-negatively are computed by a replica trick, which is successfully implemented by novel near-extremal saddle points. Gravitationally, these are identified as semiclassical solutions for dynamical branes which dominate on-shell near extremality and provide a non-perturbative resummation of off-shell topological expansions.  \nThe introductory sections that follow contextualize these issues, explain the basic notions needed to tackle them, and summarize our main results. For the sake of conciseness, we keep the presentation here minimal, and relegate the discussion of additional details to appendices.', '1.1 Effective Gravity Ensembles': 'It is expected that a quantum gravitational system, such as a black hole, should be describable from the outside as a standard, unitary quantum system with a discrete spectrum [12-18]. This is well motivated in a complete theory of quantum gravity; however, an effective theory need not adhere to such an expectation. Indeed, it has been observed that effective quantum gravity does not behave like any single quantum system, but like an ensemble thereof. This correspondence has been realized very explicitly in low-dimensional models of gravity [19-29], and argued for not only in effective gravity [30-37], but also in general effective theory [38].  \nA simple class of observables which illustrates this phenomenon consists of gravitational correlation functions involving multiple asymptotic boundaries. Denote the gravitational path integral by P and consider imposing thermal boundary conditions Z ( β ) at different β . Then for an m -point function one finds  \nP ( Z ( β 1 ) Z ( β 2 ) · · · Z ( β m )) = ⟨ Z ( β 1 ) Z ( β 2 ) · · · Z ( β m ) ⟩ , (1.2)  \nwhere on the right-hand side ⟨ · ⟩ is an average over an ensemble of non-gravitational quantum theories with canonical partition functions Z ( β ) with the corresponding β values. 1 For instance, the object Z appearing in eq. (1.1) corresponds to  \nZ ( β ) ≡ P ( Z ( β )) = ⟨ Z ( β ) ⟩ . (1.3)  \nThe emergence of an ensemble description may be attributed to the breaking of factorization on the gravity side. This occurs because P allows for contributions from all topologies with the requisite boundary conditions, including multi-boundary wormholes connecting different boundaries.  \nGiven that the gravitational path integral provides only an effective description of quantum gravity, there is also a simple heuristic reason why such an ensemble description ought to arise: if an effective theory admits multiple microscopic completions, its predictions should be consistent with but maximally ignorant of the microscopics of any one of them [37-39]. In other words, as soon as one begins to ignore details about a complete theory, all that an effective theory may have access to is expectation values capturing the statistics across the ensemble of all possible completions of the theory. 2 Hereon, we thus take the viewpoint that effective quantum gravity behaves like an ensemble of quantum theories in the sense of eq. (1.3).  \nSince a canonical partition function is given in terms of the Hamiltonian H of a theory by  \nZ ( β ) ≡ Tr e -βH , (1.4)  \none can always think of an ensemble of theories as specified by an ensemble of Hamiltonians. In turn, since Z ( β ) actually only depends on the energy spectrum of H , one can generally focus on  \nstudying ensembles of spectra and think of ⟨ · ⟩ as computing expectation values in them. 3 All in all, correlations across theories are seen to induce correlations across their energy levels, the statistics of which is what spectral ensembles capture. In these ensembles, the discrete energy levels of each theory combine into probabilistic distributions which are generically continuous. The drastic consequences this phenomenon has on entropic quantities cannot be overlooked, and one of the goals of this paper is to characterize these effects.', '1.2 Annealing and Quenching': "It is important to carefully identify which quantities are actually of interest when faced with an ensemble of theories. Consider an observable which for the partition function Z of a single theory is computed by F ( Z ). In an ensemble, there are two natural classes of quantities one may compute:  \n- · Annealed : These involve averaging Z over theories, then evaluating F on just ⟨ Z ⟩ :  \nF a ( Z ) ≡ F ( ⟨ Z ⟩ ) . (1.5)  \n- · Quenched : These involve evaluating F on every Z , then averaging F ( Z ) over theories:  \nF q ( Z ) ≡ ⟨ F ( Z ) ⟩ . (1.6)  \nBy definition, F q ( Z ) will clearly preserve every property of F ( Z ) that the linearity of the expectation value respects. In contrast, nothing guarantees that F a ( Z ) should preserve any such property of F ( Z ), and thus in general it will not. For instance, if F ( Z ) ≥ 0 always then clearly F q ( Z ) ≥ 0 as well, but one should not be conflicted if F a ( Z ) happens to take negative values.  \nWe wish to apply the distinction between annealing and quenching to the thermal entropy. In a single quantum theory with canonical partition function Z ( β ), this is given by 4  \nS ( β ) ≡ (1 -β∂ β ) log Z ( β ) . (1.7)  \nHence, in an ensemble of theories, eq. (1.5) defines the annealed entropy as  \nS a ( β ) ≡ (1 -β∂ β ) log ⟨ Z ( β ) ⟩ , (1.8)  \nwhereas eq. (1.6) defines the quenched entropy as  \nS q ( β ) ≡ (1 -β∂ β ) ⟨ log Z ( β ) ⟩ . (1.9)  \nComparing eq. (1.1) with eq. (1.8), it immediately follows by eq. (1.3) that the gravitational entropy is an annealed entropy [11]. We emphasize that the identification of S with S a also holds for a single theory, and thus makes no assumption about effective gravity generally having a non-trivial  \nensemble description (cf. footnote 1). In other words, the gravitational entropy can always be said to be an annealed entropy. This realization is the starting point for addressing goal (1) above, and will allow us to make concrete predictions about how S ( β ) must behave as a function of β in effective quantum gravity. For instance, we prove that it generically diverges negatively as the temperature goes to zero, implying that the gravitational entropy of generic black holes must become negative near extremality. Regarding goal (2), as the reader may suspect, the entropic quantity with the desired thermodynamic behavior is the quenched entropy.  \nThe discrepancy between S a and S q clearly boils down to the difference between annealing and quenching the logarithm function. It will thus at times be useful to consider the free energy,  \nF ( β ) ≡ -1 β log Z ( β ) , (1.10)  \na physical quantity that is linear in log Z ( β ) and in terms of which eq. (1.7) simply reads  \nS ( β ) = β 2 F ' ( β ) . (1.11)  \nClearly S is non-negative if and only if F is monotonically non-decreasing, which means that negativity of the annealed entropy is equivalent to non-monotonicity of the annealed free energy. In a single theory, the condition F ' ( β ) ≥ 0 is equivalent to log Z ( β ) ≥ β∂ β log Z ( β ), which gives the physical requirement F ( β ) ≤ ⟨ E ⟩ β , where ⟨ E ⟩ β is the statistical average energy of the thermal system. Our results thus imply that this property is violated when annealing.", '1.3 The No-Replica Trick': "Annealed quantities are usually easier to calculate, since their linearity in the partition function means that standard path integral techniques are applicable. In contrast, quenching is generally non-trivial and has to be addressed case by case depending on the pertinent observable. This is why annealing is sometimes used as a proxy for quenching. However, the success of such an approximation relies on whether the observable at hand happens to be self-averaging, a property which is 'a priori hard to ascertain. For instance, we will see that annealed entropies reliably approximate quenched ones at high temperatures, but radically disagree at low temperatures.  \nGiven the possibility of discrepancy, in the ever surprising context of quantum gravity it would seem strictly more desirable to always consider quenched quantities. 5 Unfortunately though, except for the recent attempts of [11, 40], the tools for quenching quantities such as the entropy directly within quantum gravity remain completely mysterious. The difficulties one faces when attempting such a computation are conceptual, technical, and one may worry potentially fundamental. A central corollary of this paper in regard to goal (2) is that, in fact, no fundamental obstruction seems to exist: effective gravity does have the requisite ingredients for calculating quenched entropies.  \nTo understand the difficulties, note that given an ensemble of theories with some measure, the calculation in eq. (1.9) is conceptually clear even if technically involved. In gravity, however, an  \nobject like P (log Z ( β )) cannot even be directly defined. This is because there is no separate notion of an ensemble of theories and its measure; the gravitational path integral defines them all at once. This hurdle can be overcome by means of a well-known replica trick for the logarithm function,  \n⟨ log Z ( β ) ⟩ = lim m → 0 ⟨ Z ( β ) m ⟩ -1 m = lim m → 0 d dm ⟨ Z ( β ) m ⟩ , (1.12)  \nthe gravitational implementation of which was pioneered by [11]. The no-replica limit m → 0 in eq. (1.12) is very different from the limit m → 1 in the standard replica trick for the von Neumann entropy, which is already well-understood in gravity [4]. The two identities in eq. (1.12) can be used interchangeably, but consistency clearly requires both limits to exist and agree. The intermediate expression emphasizes that ⟨ · ⟩ must be a unit-normalized expectation value.  \nIn gravity, the computation of the moments ⟨ Z ( β ) m ⟩ in eq. (1.12) can be performed by the gravitational path integral P ( Z ( β ) m ) with identical boundary conditions Z ( β ) on m disjoint boundaries. For integer m ≥ 1 these are well-defined m -point gravitational correlators like those in eq. (1.2). The replica trick in eq. (1.12), however, requires an analytic continuation to real m> 0. Ambiguities in this continuation posed technical difficulties in [11], and a semiclassical resolution was proposed in [40]. Nonetheless, the need to resort to numerics still rendered the gravitational implementation of this method inconclusive. There thus remained the troubling possibility that effective gravity would not even be able to compute quenched entropies at all.  \nIn this work, we address this question by making use of the duality between theories of JackiwTeitelboim (JT) gravity [41, 42] and matrix integrals [19, 20, 43, 44], in which the ensembles of theories that emerge correspond to a random matrix Hamiltonian. In particular, this duality is established at the level of a topological expansion of the gravitational correlators on the left-hand side of eq. (1.2) and a corresponding perturbative expansion of a double-scaled matrix integral for the right-hand side. While this correspondence is a precise match to all orders in the topological expansion, the random matrix formulation in general has access to non-perturbative details that may not be accessible from the gravity side. 6 Our goal in this paper is to study the replica trick in eq. (1.12) using the concrete framework of random matrix theory, but making sure to only rely on ingredients which are accessible to the gravitational path integral. 7 As a result, we are able to pinpoint the minimal ingredients necessary to succeed in this task, singling out the key role of a specific class of branes which effective gravity can indeed accommodate.", '1.4 Gravitational Entropy': 'While the subject of quenched entropies in gravity is very recent and underexplored, extensive work has been devoted to calculating the annealed entropy S that eq. (1.1) gives. It is of particular interest to study S ( β ) for black holes near extremality, corresponding to the parametric hierarchy β ≫ S 0 ≫ 1. This regime can be explored by looking at the largeβ behavior of S ( β ) on top  \nof a largeS 0 perturbative expansion. The computation of this near-extremal entropy including quantum gravitational effects is a formidable endeavor which has been pursued by [51-66]. We note that all of their results can be captured by the following general expression: 8  \nS ( β ) = S 0 -s log β +#∆ e -β ∆ +# β -1 e -β ∆ + O ( β -2 e -β ∆ ) , (1.13)  \nwhere we are using # to denote O (1) numerical factors whose precise value will not play a role in our discussion. Here, S 0 is the quantum-corrected Bekenstein-Hawking black hole entropy at extremality, 9 s ∈ Q is an important parameter determined by supersymmetry, and ∆ ≥ 0 is the size of a possible gap in the spectrum above the ground state. In particular, one observes that s = 0 when the extremal black hole is BPS [58] and strictly s > 0 otherwise [57]; e.g. theories which are purely bosonic yield s = 3 / 2. As for ∆, in general one finds ∆ > 0 if and only if the extremal limit is BPS and has an extensive ground-state degeneracy N 0 ∼ O ( e S 0 ); otherwise, ∆ = 0 when N 0 ∼ O (1). Summarizing, eq. (1.13) is found to involve:  \n         Non-near-BPS: s > 0 , ∆ = 0 , Near-BPS: s = 0 ,    N 0 ∼ O (1) : ∆ = 0 , N 0 ∼ O ( e S 0 ) : ∆ > 0 , (1.14)  \nwhere we are using the term near-BPS to refer to near-extremal black holes which preserve some supersymmetry in the extremal limit. Note that at this point eq. (1.14) is a heuristic classification obtained by inspecting the form of eq. (1.13) for various gravitational results. In this paper, we show that eq. (1.14) holds in full generality, and explain the quantitative origin of both s and ∆.  \nFocus for now on the differences between the near-BPS ( s = 0) and non-near-BPS ( s > 0) cases. The behavior of S ( β ) in eq. (1.13) is drastically different at large β depending on supersymmetry due to the -s log β term. For near-BPS black holes, the β →∞ limit simply yields a finite extremal entropy given by S 0 . Generically, this is a large number interpreted as quantifying the ground-state degeneracy by e S 0 , a satisfying result since one does expect a large number of degenerate BPS states based on symmetry arguments. That this extremal entropy indeed has a statistical interpretation has been verified by direct enumeration of BPS states in string theory [12, 17]. On the other hand, for black holes which are not near-BPS, we observe S becoming smaller with β and eventually turning negative for β ≳ O ( e S 0 ). The monotonic decrease of S towards low temperatures has been interpreted positively, considering that in the absence of a symmetry arguments non-BPS black holes should not be highly degenerate [57, 61]. As for the eventual negativity of S ( β ), it has been argued that for such large β ≳ O ( e S 0 ) some other quantum gravitational object must dominate Z ( β ), giving non-perturbative corrections to eq. (1.13) which would take over and keep S ( β ) nonnegative [64, 65]. Statistically, this would effectively mean that non-near-BPS black holes at such low temperatures do not exist. Our findings in this paper suggest otherwise.  \nThe qualitative features of eq. (1.13) that follow from eq. (1.14) will serve all along this paper as both an example and a consistency check for many results. Let us demonstrate this with a brief  \npreview of the general picture that arises from our spectral analysis. What turns out to characterize near-BPS black holes is the fact that the ground-state energy is protected by supersymmetry across the ensemble, whereas for black holes which are not near-BPS this energy is allowed to vary. As we will show, this simple property of a spectral ensemble alone already implies that the annealed entropy should capture the degeneracy of the ground-state energy when protected, and diverge negatively when not. Hence the behavior eq. (1.13) exhibits is perfectly consistent with S being an annealed entropy, including the fact that it becomes negatively divergent in the β →∞ limit. In other words, any potential non-perturbative correction to eq. (1.13) should not make S nonnegative, and nonetheless nothing compromises the existence of such non-near-BPS black holes.  \nAnother important fact we will prove is that how degenerate the ground state is actually plays no role on the negativity issue from the spectral point of view; all that matters is whether its energy is protected. This is already suggested by eq. (1.14), where we see that the sufficient condition for non-negativity s = 0 holds in the near-BPS case regardless of whether N 0 is extensive in e S 0 or not. In particular, we will show that the annealed entropy asymptotes to log N 0 at extremality so long as the ground-state energy is protected, no matter how small or large N 0 is. On the other hand, we prove that this entropy becomes negative at low temperatures if the ground-state energy is not protected, regardless of how degenerate it is. 10 Black holes in supersymmetric theories which break all supersymmetry will nicely illustrate this fact. Hence what makes the behavior of the gravitational entropy of near-BPS black holes so different is the fact that the ground-state energy is protected, not that it is highly-degenerate. The characteristic feature that a ground state with high degeneracy does give rise to, as will be explicitly seen, is the appearance of a gap in the spectrum above its energy. In other words, the ∆ gap in eq. (1.13) does arise precisely because the ground state is highly degenerate, and would be zero otherwise (cf. eq. (1.14)).', '1.5 Summary of Results': 'As explained in section 1.2, the gravitational entropy S computed in effective quantum gravity through eq. (1.1) is an annealed entropy. Using this fact, our results for goal (1) are as follows:  \n- 1. The standard gravitational path integral gives a black hole entropy S which is non-negative at low temperatures if and only if the extremal limit is supersymmetric, and diverges negatively otherwise. This is shown to be an exact, non-perturbative statement for any effective theory of gravity with an ensemble description. This result turns the generic negativity of S from a pathological behavior that ought to be fixed by non-perturbative corrections into a required feature which should be robust against non-perturbative corrections.  \n⟨ Z ( β ) ⟩ = ( e βδ -1) e S 0 βδ + √ πe S 0 2 β 3 / 2 .  \nOne easily verifies that the strict δ → 0 case, corresponding to a ground-state energy protected to E = 0 across the ensemble, yields S a ( β ) → S 0 as β → ∞ . Note that this fact is independent of how large or small S 0 is, and that S a ( β ) ≥ 0 in this case even though there is no gap in the spectrum. In contrast, for any δ > 0 one finds S a ( β ) →-∞ as β →∞ , no matter how large e S 0 is; in this toy model, the negativity of the annealed entropy occurs for β ≳ e S 0 /δ .  \n- 2. If the extremal limit is a BPS black hole, the usual gravitational entropy stays non-negative and asymptotes to the Bekenstein-Hawking entropy. Nonetheless, away from the strict extremal limit, this quantity continues to behave unlike the thermal entropy of any quantum system. We identify the supersymmetry protection of the ground-state energy as the sole mechanism behind the qualitative negativity difference between the BPS and non-BPS limits; whether or not the ground state is highly degenerate also has consequences which are interesting but unrelated.  \nThese results follow primarily from a general analysis of spectral ensembles in section 2. This provides the framework needed to address the differences between annealed and quenched entropies given an ensemble of theories, and make very general predictions about the behavior of these entropies in effective quantum gravity. These predictions can be compared with explicit gravitational results for near-extremal black holes that we review in section 2.4 and appendix A. The consequences of supersymmetry depending on whether the ground-state energy is protected across Hamiltonians in the ensemble are stated in theorem 1 and corollary 1. The discrepancy between annealed and quenched entropies even with supersymmetry, as well as the actual role of large BPS degeneracy, are more explicitly explored later in sections 3.6.2, 4.5 and 5.2.  \nAs for the quenched entropy, its calculation in gravity requires computing gravitational correlators of Z ( β ) including wormhole contributions to the gravitational path integral. However, the physics of interest occur at β ∼ O ( e S 0 ), a non-perturbative regime often regarded as inaccessible to effective gravity. Our findings using general ensembles of random matrix theory dual to effective gravity strongly suggest otherwise; in fact, already semiclassical gravity suffices to capture these effects. Using matrix integrals to study ⟨ Z ( β ) m ⟩ , our results for goal (2) are as follows:  \n- 1. For a completely general random matrix ensemble, we find that at β ∼ O ( e S 0 ) a new saddle arises which exists and dominates at large e S 0 . This saddle is a single-eigenvalue instanton with a clear holographic dual: dynamical brane solutions. The saddle-point treatment of the matrix integral is thus understood to be dual to a semiclassical treatment of effective gravity where brane actions are evaluated on-shell. As a non-trivial test of this claim, we show that our instanton results for β ∼ O ( e S 0 ) identically match the on-shell action of gravitational D3-branes holographically dual to 1 / 2 BPS Wilson loops in N = 4 supersymmetric Yang-Mills with a random matrix description.\n- 2. These brane solutions provide a controlled resummation that includes contributions of arbitrary genus to the gravitational path integral in a β ∼ O ( e S 0 ) regime where the usual topological expansion breaks down. In addition, our instantons continuously interpolate between small β ≪ O ( e S 0 ) and large β ∼ O ( e S 0 ), suggesting that semiclassical gravity alone should also be able to capture the transition to a brane regime. For near-extremal black holes, this transition suggests that the near-horizon throat gets cut off by a brane at a finite proper distance as soon as β ∼ O ( e S 0 ) already at the semiclassical level. 11\n- 3. Our instantons are found to be necessary and sufficient to obtain a non-negative quenched entropy via the no-replica trick, motivating its gravitational implementation including brane replica  \nwormholes. In addition, they capture non-perturbative effects which enhance the generic negative divergence of the annealed entropy at large β . Our analysis of the replica trick resolves a subtlety that has hindered previous attempts even within matrix integrals: 12 the largeβ and m → 0 replica limits directly compete and thus do not commute. This dooms any approach that involves approximating integer m ≥ 1 moments at large β , and then analytically continuing to m → 0.  \nThe general framework of random matrix ensembles and its constructs relevant for the analysis of random Hamiltonians with and without supersymmetry is the subject of section 3. In particular, as explained in sections 3.2 and 3.3 (see also appendix C), our matrix integral approach applies to all matrix ensembles relevant to quantum gravity: Wigner-Dyson (no supersymmetry), AltlandZirnbauer (broken supersymmetry), and their Wishart generalizations (unbroken supersymmetry). Matrix scalings to spectral edges and double scaling limits are also accounted for (see section 3.6).  \nThese tools allow us to carry out a general exploration of the near-extremal limit in section 4. Once the relevant matrix dynamics are identified to be eigenvalue instantons, an appropriate continuum formulation of the matrix integral for ⟨ Z ( β ) m ⟩ is performed in sections 4.2 to 4.4. The consequences of these instantons for entropic quantities are analyzed in section 4.5. A crucial aspect of our treatment throughout section 4 is the restriction to employing solely random matrix theory tools with a clear semiclassical holographic dual. 13 This is what allows us to elucidate the key gravitational ingredients at play in section 4.6. Furthermore, by working at a saddle-point level, the lessons we arrive at make sense also in higher-dimensional gravity, where an off-shell gravitational path integral is not even available. Explicit examples of our instanton results for various matrix ensembles are presented in section 5, including the cases relevant to 1 / 2 BPS Wilson loops in section 5.1 and JT gravity in section 5.4.  \nIt is worth emphasizing that many of our conclusions benefit from the robustness of universality results in random matrix theory, and thus are independent at a semiclassical level from any potentially non-unique choice of non-perturbative completion of the theory. 14 More specifically, we find the dynamics of the relevant near-extremal instantons to be governed by edge statistics which are universal to leading order at large e S 0 . By analyzing universal aspects of the low energy spectrum of general matrix ensembles, we are able to make the following general observations and predictions about the near-extremal black hole spectrum in effective quantum gravity:  \n- 1. The spectral statistics of Wigner-Dyson (non-supersymmetric) and Altland-Zirnbauer (supersymmetric) ensembles are often associated to soft Airy edges and hard Bessel edges, respectively. In fact, we find that in the more general Wishart construction for supersymmetric ensembles, the lower edge of the spectrum becomes soft and governed by Airy statistics. This is a consequence  \nof the natural scaling of Wishart ensembles, which gives rise to a spectral gap above zero energy. This gap pushes the spectrum off the hard Bessel edge at zero, causing positive eigenvalues to have a soft Airy onset above the gap.  \n- 2. By studying random Hamiltonians with at least N = 2 supersymmetry, we identify two distinct physical mechanisms explaining the origin of gaps in the spectrum. For supermultiplets with an extensive degeneracy of BPS states at zero energy N 0 ∼ ν e S 0 , the spectrum exhibits a gap ∆ to the first non-BPS state which behaves as ∆ ∼ ν 2 for small gaps and as ∆ ∼ ν for large gaps. This gap arises due to the Vandermonde repulsion exerted by the BPS states at zero energy. For purely non-BPS sectors, we observe that what causes the shift of the spectrum to higher energies for higher supermultiplets is a monotonic growth extensive in e S 0 in the number of non-BPS states.  \nThe universal gap behavior we observe is consistent with gravitational results. In particular, the quadratic scaling of ∆ with the ratio N 0 /e S 0 matches the behavior of the gap in supermultiplets of N = 2 super-JT, where ∆ ∼ q 2 when N 0 ∼ q e S 0 (see sections 2.4 and 3.6.2). When N 0 > 0 but not extensive in e S 0 , the spectrum becomes gapless and develops a hard edge again even though a discrete ground state protected by supersymmetry still remains. This is again consistent with gravity results for deformed models of N = 1 super-JT with a non-extensive number of BPS states, and a particularly anomalous sector of the theory with N = 2 supersymmetry. The argument for spectral gaps in purely non-BPS supermultiplets is first made in section 3.3, and a random matrix construction is given in appendix C.', '2 Spectral Ensembles': 'We introduce here the general structure in an ensemble of theories that is relevant to the study of annealed and quenched entropies. This allows us to make completely general statements about how these quantities behave and differ from each other. Known gravitational entropy results are then presented and reviewed through the lens of spectral ensembles. The tools we employ here will also emerge and be explicitly realized in section 3 when studying random matrix ensembles.', '2.1 Basics': 'Given a theory with Hamiltonian H , the canonical partition function Z ( β ) is given by eq. (1.4). Let H have a discrete spectrum of eigenvalues E ≡ { E k ∈ R } . The spectral density of H is  \nρ E ( E ) ≡ ∑ k δ ( E -E k ) , (2.1)  \nand Z ( β ) can more explicitly be written  \nZ ( β ) = ∑ k e -βE k = ˆ R dE ρ E ( E ) e -βE . (2.2)  \nAs anticipated below eq. (1.4), given that Z ( β ) only depends on the spectrum of the Hamiltonian, we can generally reduce our theory ensembles to spectral ensembles. Specifically, an ensemble of  \ntheories given by a probability density function (PDF) over some space of Hamiltonians H yields a unique joint PDF over their eigenvalues E upon diagonalization. This joint PDF defines the measure p ( E ) of the spectral ensemble. Equipped with it, the expectation value in the spectral ensemble of any function f of E is defined by  \n⟨ f ⟩ ≡ ˆ d E p ( E ) f ( E ) . (2.3)  \nFor N eigenvalues, one can generally write p ( E ) in terms of a symmetric function p N as  \np ( E ) = p N ( E 1 , . . . , E N ) . (2.4)  \nThe joint PDF induced on any number 1 ≤ m ≤ N of eigenvalues can be obtained by simply integrating p N over the other set of N -m eigenvalues,  \np n ( E 1 , . . . , E n ) ≡ ˆ R N -m dE m +1 · · · dE N p N ( E 1 , . . . , E N ) . (2.5)  \nThese can be understood as m -point eigenvalue correlation functions, and simply related to correlators of the discrete spectral density in eq. (2.1) by 15  \n⟨ ρ ( E 1 ) · · · ρ ( E n ) ⟩ = N m p n ( E 1 , . . . , E n ) . (2.6)  \nFor m = 1, this defines the average spectral density in terms of the 1-point function as  \n⟨ ρ ( E ) ⟩ = Np 1 ( E ) , (2.7)  \nWhen clear from context, the average ⟨ ρ ( E ) ⟩ will be referred to simply as the spectral density.  \nThe definitions of annealed and quenched entropies in eqs. (1.8) and (1.9) can now be made more explicit. For the logarithm function, annealing leads to  \nlog ⟨ Z ( β ) ⟩ = log   N ˆ R dE p 1 ( E ) e -βE   , (2.8)  \nwhere we have simply used eqs. (2.2) and (2.7). Hence the annealed logarithm solely depends on the 1-point function of the spectral ensemble, which is why it is simply expressible in terms of the spectral density. In contrast, when quenching we are stuck with  \n⟨ log Z ( β ) ⟩ = ˆ R N dE 1 · · · dE N p N ( E 1 , . . . , E N ) log ( N ∑ k =1 e -βE k ) , (2.9)  \nwhich makes manifest that the quenched logarithm depends on the full spectral N -point function. Clearly, the annealed entropy in eq. (1.8) is a much coarser quantity determined by only 1-point statistics, whereas the quenched entropy in eq. (1.9) is sensitive to all higher-point statistics. Needless to say, the quenched entropy is not computable from just knowledge of the average spectral density, unless the m -point functions are all trivial. 16 More explicitly, eqs. (2.8) and (2.9) are equal only in the trivial case in which all m -point correlators in eq. (2.6) factorize as  \n⟨ ρ ( E 1 ) · · · ρ ( E n ) ⟩ = ⟨ ρ ( E 1 ) ⟩ · · · ⟨ ρ ( E n ) ⟩ , (2.10)  \ni.e., an uninteresting ensemble of independently and identically distributed random energies.', '2.2 Quenched Expectations': 'The logarithm expressions in eqs. (2.8) and (2.9) imply that annealed and quenched entropies are generally very different quantities in non-trivial ensembles. The only room for agreement between them is if in some regime in β the quenched logarithm becomes dominated by 1-point statistics. Expanding log Z ( β ) about β = 0, clearly a term of O ( β k ) can only possibly depend on up to k distinct eigenvalues. Hence at small β ,  \n⟨ log Z ( β ) ⟩ = log N -β ⟨ E ⟩ + O ( β 2 ) , (2.11)  \nwhich clearly matches log ⟨ Z ( β ) ⟩ up to O ( β ). Contributions to eq. (2.11) from 2-point functions begin at O ( β 2 ), and thus the annealed logarithm begins missing the effect of their connected parts already at this order. Similarly, k -point functions begin contributing to ⟨ log Z ( β ) ⟩ at O ( β k ), and their connected parts cannot possibly be captured by log ⟨ Z ( β ) ⟩ . For any non-trivial spectral ensemble, the annealed logarithm thus unavoidably begins to differ from the quenched one already at O ( β 2 ). As a result, we expect the annealed entropy to approach the quenched entropy at high temperatures, but become a worse approximation to it at lower temperatures if the connected parts of k -point functions contribute non-negligibly.  \nSince higher-point functions become important at higher orders in β , we learn that ⟨ log Z ( β ) ⟩ becomes sensitive to the full statistics of the ensemble as β increases. At large β one obtains  \n⟨ log Z ( β ) ⟩ = ⟨ log N 0 ⟩ -β ⟨ E 0 ⟩ + 〈 O ( e -β ∆ ) 〉 , (2.12)  \nwhere E 0 ≡ min { E 1 , . . . , E N } should be understood as the ground-state energy, N 0 ≥ 1 is its degeneracy, and the exponentially suppressed corrections depend on the statistics of the separation ∆between E 0 and the next higher eigenvalue. Although N 0 may sometimes be a fixed constant, the expectation value of E 0 does generically depend non-trivially on the full spectral N -point function. In other words, the largeβ behavior of ⟨ log Z ( β ) ⟩ can by no means be reproduced by log ⟨ Z ( β ) ⟩ .  \nThere is a simple exception to this which turns out to be of much interest. In particular, it is possible to make annealed and quenched logarithms agree to leading order at largeβ as follows. Suppose some symmetry of the original ensemble of Hamiltonians makes their ground-state energy a protected quantity. 17 In this case the expectation value of the ground-state energy in eq. (2.12) just yields ⟨ E 0 ⟩ = E 0 . On the other hand, one easily sees that the greatest lower bound on the support of the average spectral density is inf supp ⟨ ρ ⟩ = E 0 , and that  \n⟨ ρ ( E ) ⟩ ⊃ ⟨ N 0 ⟩ δ ( E -E 0 ) . (2.13)  \nConsequently, in this non-generic case the annealed logarithm in eq. (2.8) would behave just like the quenched one to leading order at large β , up to a potential discrepancy between ⟨ log N 0 ⟩ and log ⟨ N 0 ⟩ . In the concrete realizations of this setting that we will study using matrix ensembles, N 0 is in fact fixed and thus this constants agree. Even in this case though, discrepancies between annealing and quenching do still arise at large β at orders higher than just discussed.  \nIt is illustrative to be more explicit about how ⟨ ρ ⟩ preserves a discrete ground state as in eq. (2.13) when its energy is protected. The average spectral density can be thought of as the  \nresult of smearing the discrete energy levels of each Hamiltonian in the ensemble as p ( E ) dictates. Hence consider letting θ k be the map on E which yields 1 if say E 1 is the k th energy level, and zero otherwise. The sum over all these indicator functions θ k adds up to 1, and thus  \nρ k ( E ) ≡ N ˆ d E p ( E ) θ k ( E ) ρ E ( E ) , (2.14)  \ngives the probability distribution of values that the k th energy level takes across the ensemble. The average spectral density is then suggestively a superposition of these for every energy level, 18  \n⟨ ρ ( E ) ⟩ = N ∑ k =1 ρ k ( E ) . (2.15)  \nThis expression grants the interpretation of ⟨ ρ ( E ) ⟩ as capturing the statistical uncertainty in the energy levels of the random Hamiltonian. The non-generic case described above where the groundstate energy E 0 is protected corresponds to having no uncertainty in the lowest eigenvalue, i.e.,  \nρ 1 ( E ) = ⟨ N 0 ⟩ δ ( E -E 0 ) . (2.16)  \nBy construction, in this case no ρ k with k ≥ 2 can possibly have a tail with support below E = E 0 , and thus we indeed verify that inf supp ⟨ ρ ⟩ = E 0 .  \nThe observations here about log Z ( β ) lead us to expect large discrepancies between annealing and quenching when β is large, and also point at a non-generic situation in which they may behave similarly. In trying to understand how the gravitational entropy compares to a thermal entropy, we are thus drawn to study how annealed entropies do actually behave at large β in these situations. As we will show in section 2.3, eq. (2.16) turns out to be necessary and sufficient for the annealed entropy to stay non-negative at low temperatures. Before we turn to this though, let us conclude with two brief remarks.  \nAs emphasized above, while log ⟨ Z ( β ) ⟩ only captures 1-point statistics, ⟨ log Z ( β ) ⟩ depends on all-point statistics and is thus sensitive to much finer details of the spectral ensemble. The investigation of quantities probing higher-point functions has proven remarkably insightful in quantum gravity, as most prominently demonstrated by calculations of the spectral form factor [19, 67-72]. However, this object basically involves the second moment ⟨ Z ( β ) 2 ⟩ , which only depends on up to 2-point functions. The regime in which the spectral form factor becomes interesting corresponds precisely to when ⟨ log Z ( β ) ⟩ begins to differ from log ⟨ Z ( β ) ⟩ . Being sensitive to even higher-point statistics, one may thus expect quenched entropies to capture even more interesting physics.  \nFinally, note that our analysis above says that ⟨ log Z ( β ) ⟩ becomes sensitive to m -point functions with increasingly higher m as β becomes large. On the other hand, the replica trick in eq. (1.12)  \n̸  \nsays that to compute this quantity we should consider m -point functions in the m → 0 limit. This is our first hint that β → ∞ and m → 0 are competing limits whose order matters. Since the latter is strict in the replica trick, this already tells us that doing largeβ approximations at finite m before taking the m → 0 replica limit may be a terrible idea.', '2.3 Negative Entropies': 'Since this section focuses solely on annealed quantities, the notation ⟨ · ⟩ will be temporarily dropped for the sake of clarity. Annealed entropies just depend on the average spectral density, and thus can be studied in terms of the properties of a general such object. For the sake of generality, let this spectral density be given by an arbitrary map  \n˜ ρ : R → [0 , ∞ ] , (2.17)  \nwhere we allow for distributions by including infinity in the codomain. We demand ˜ ρ to be locally integrable, and assume its support is bounded from below to define a ground-state energy by  \nE 0 ≡ inf supp ˜ ρ. (2.18)  \nThe canonical partition function is then given by 19  \nZ ( β ) = ˆ R dE ˜ ρ ( E ) e -βE = ˆ ∞ 0 -dE ρ ( E ) e -βE , (2.19)  \nwhere in the second equality we have introduced a canonically renormalized spectral density,  \nρ ( E ) ≡ e -βE 0 ˜ ρ ( E 0 + E ) . (2.20)  \nThis is useful because now Z takes the form of a Laplace transform. In what follows, we refer to the canonically shifted ρ whose ground-state energy is zero as the spectral density.  \nAllowing for ρ to be discontinuous and distributional, the general analysis of its Laplace transform can get cumbersome. Motivated by section 2.2, we will assume that ρ admits a decomposition into distributional and continuous parts of the general form  \nρ ( E ) = ∆( E ) + g ( E ) , ∆( E ) ≡ N 0 δ ( E ) + ∑ k ∈ I N k δ ( E -E k ) , (2.21)  \nfor some set I ⊆ N (possibly empty or infinite), N 0 ≥ 0, N k > 0, all E k > 0 distinct, and where g is now an actual real function, assumed to be non-negative and piecewise continuous. By the canonical shift in eq. (2.20), it follows that either inf supp g = 0, or N 0 > 0, or both. Generally g may contain discontinuities and singularities, so when studying it we will generally be interested in keeping track of one-sided limits, and will not require two-sided limits to exist. Our convention will be to allow for limits to exist in the extended non-negative reals [0 , ∞ ].  \nThe Laplace transform for ρ of the general form in eq. (2.21) reads  \nZ ( β ) = ˆ Z ( β ) + Z g ( β ) ,        ˆ Z ( β ) ≡ N 0 + ∑ k ∈ I N k e -βE k Z g ( β ) ≡ ˆ ∞ 0 + dE g ( E ) e -βE (2.22)  \nWe will only be interested in spectral densities ρ whose Laplace transform Z exists, which requires ˆ Z and Z g to independently exist. For ˆ Z , if its form in eq. (2.22) involves an infinite series, we assume it converges for any β > c for some finite c > 0. As for Z g , we assume g is locally integrable and of exponential type, 20 which guarantees that lim β →∞ e -βE g ( E ) = 0 for any β > c for finite c > 0. Note that the local integrability condition allows for integrable singularities where g diverges. In particular g ( E ) may diverge as E → 0 + and, so long as it does so slower than E -1 , Z g will still exist.  \nWith a well-defined Z , one can then proceed to study the annealed entropy S a from eq. (1.8). Note that this object is invariant under shifts of ρ like eq. (2.20) or, at the level of Z , under Z ( β ) → e -βE 0 Z ( β ) for any E 0 ∈ R . This demonstrates that our canonical shift to always fix inf supp ρ = 0 was without loss of generality. The main result of this section may now be stated:  \nTheorem 1. Let S a ( β ) be the annealed entropy defined in eq. (1.8) , with ⟨ Z ( β ) ⟩ given by the Laplace transform of an average spectral density ρ ( E ) of the general form in eq. (2.21) . Then:  \nS a ( β ) ≥ 0 ⇐⇒ ρ ( E ) ⊃ N 0 δ ( E ) with N 0 ≥ 1 .  \nIf instead N 0 < 1 , then S a ( β ) < 0 , ∀ β > β ∗ for some β ∗ > 0 , and lim β →∞ S ( β ) = -∞ for N 0 = 0 . 21  \nThe proof of this statement is given in appendix B. In words, theorem 1 says that the necessary and sufficient condition for the annealed entropy to be non-negative is that the average spectral density preserve at least one discrete ground state. This both confirms and explains the general behavior that the gravitational entropy is observed to exhibit in eq. (1.13). In particular, by theorem 1 the negativity of S ( β ) at large β when no supersymmetry protects the energy at extremality is not an artifact of the perturbative treatment, but a genuine property that the gravitational entropy must exhibit as an annealed entropy. This negativity can thus be read as a remarkable consistency check of the ensemble description dual to the gravitational path integral. The physical restatement of theorem 1 may be phrased as follows:  \nCorollary 1. Let S ( β ) be the gravitational entropy defined in eq. (1.1) for a partition function Z ( β ) given by the Euclidean gravitational path integral. Then S ( β ) is non-negative for all β > 0 if and only if the ground-state energy is protected. Otherwise, S ( β ) →-∞ continuously as β →∞ .', '2.4 Black Hole Spectra': 'The spectral densities that realize the near-extremal gravitational entropies in eq. (1.13) nicely illustrate these results. The β dependence relevant to our discussion comes from strongly coupled gravitational dynamics that arise in the throat of black holes near extremality. These are universally governed by a quantum mechanical model known as the Schwarzian theory and its supersymmetric generalizations [56-59, 61-63, 65, 66, 73-75] (see appendix A for more details). The path integral of these Schwarzian theories has been shown to be 1-loop exact and fully solved, thus providing a robust framework for the study of near-extremal black hole spectra [58, 67, 76-81].  \nFigure 1 . Black hole spectra (top) and gravitational entropies (bottom) from each N = 0 , 1 , 2 , 4 superSchwarzian theory (left to right). Spectra are respectively given by eqs. (2.24), (2.25), (2.27) and (2.29), and the corresponding entropies are obtained by Laplace transform and evaluation of eq. (1.7). For the N = 2 theory we have chosen ˆ q = 3 and δ = 0. The plots for N = 2 , 4 show curves corresponding to different supermultiplets separately, and also for the total spectral density of the theory (black). Discrete ground-state contributions are symbolized by an arrow at zero of 10 times the degeneracy value for ease of visualization. These BPS states are accounted for when computing the entropy for the total spectral density, which clearly is not just a sum over supermultiplet entropies. In particular, note that the entropies of supermultiplets with no BPS states ( q = 7 / 2 for N = 2 and all J > 0 for N = 4) all diverge negatively at large β , just like the entropies for the N = 0 , 1 theories. Note the logarithmic scale of the horizontal axis.  \n<!-- image -->  \nHere, we review the gravitational spectral densities that arise from these Schwarzian theories. All cases, near-BPS or not, and gapped or not, can be concisely captured starting from an average spectral density which takes the simple form 22 [58-66]  \nρ c ( E ) = e S 0 δ ( E ) δ p ≤ 0 + e S 0 √ 2 π ( √ cE 2 πc ) p -1 I p -1 ( 2 π √ cE ) , (2.23)  \nwhere I n is the modified Bessel function of the first kind. The parameter p ∈ Q matches s = p for the non-near-BPS case of eq. (1.13), and takes values p ≤ 0 for near-BPS black holes. In consistency with theorem 1, the annealed entropy stays non-negative precisely only when the discreteness of the ground-state energy is protected by supersymmetry.  \nFor non-near-BPS black holes, one may choose units such that c = 1 and eq. (2.23) captures the complete spectral density of states. For near-BPS black holes, one can obtain the full spectral density by adding up the contribution of all possible saddles that supersymmetric localization gives [58, 77]. Each such saddle gives a spectral density of the form of eq. (2.23) with a corresponding value of c , which are to be summed over and appropriately weighted. This series can be reorganized into a sum over supermultiplets with different R-charges. Supermultiplets are statistically independent, and thus their spectral densities are the ones of physical interest [43].  \nWe now evaluate eq. (2.23) for the specific values of p found in gravitational calculations, and illustrate the results in fig. 1. In purely bosonic theories of gravity with no supersymmetry ( N = 0),  \nas mentioned below eq. (1.13), one finds p = 3 / 2 and eq. (2.23) yields [20, 57, 61, 62, 65, 66, 77, 79]  \nρ N =0 ( E ) = e S 0 sinh ( 2 π √ E ) 2 π 2 , Z N =0 ( β ) = e S 0 + π 2 /β 2 √ πβ 3 / 2 . (2.24)  \nNote that ρ N =0 ( E ) ∼ √ E near E = 0 + . Certain supersymmetric theories contain black holes which break all supersymmetry at extremality and whose near-extremal physics is governed by the N = 1 super-Schwarzian. In this case one finds p = 1 / 2 and eq. (2.23) yields [20, 37, 77, 79]  \nρ N =1 ( E ) = e S 0 cosh ( 2 π √ E ) π √ E , Z N =1 ( β ) = e S 0 + π 2 /β √ πβ 1 / 2 . (2.25)  \nNote that ρ N =1 ( E ) ∼ 1 / √ E near E = 0 + . The near-BPS case requires the emergence of a theory with more supersymmetry near extremality. For the N = 2 super-Schwarzian one finds p = 0, and the saddle-point contributions labelled by c in eq. (2.23) read [43, 58-60, 77, 79, 82]  \nρ N =2 c ( E ) = e S 0 δ ( E ) + e S 0 √ 2 πc I 1 ( 2 π √ cE ) √ E . (2.26)  \nThese contributions can be reorganized into spectral densities for supermultiplets of the theory. Here we quote the results and defer a more detailed discussion to section 3.3. The general N = 2 theory can be specified by an odd integer ˆ q > 0 and a constant δ ∈ [0 , 1), and its supermultiplets labelled by a parameter q ∈ Z + δ -1 / 2. 23 The spectral density of the q supermultiplet is [43]  \nρ N =2 q ( E ) = e S 0 sin( πq/ ˆ q ) 2 δ ( E ) δ | q | < ˆ q + e S 0 sinh ( 2 π √ E -∆ q ) 4 πE θ ( E -∆ q ) , ∆ q ≡ q 2 4ˆ q 2 , (2.27)  \n̸  \nwhere θ is the Heaviside step function. Only supermultiplets with 0 < | q | < ˆ q contain BPS states, which give rise to a discrete ground state with degeneracy proportional to e S 0 . The continuous part of the spectrum is only supported on energies E = ∆ q ≥ 0 above the ground state. The gap size ∆ q depends on the supermultiplet and is generically nonzero. We will show in section 3.6.2 that the emergence of such a positive gap can be explained from the perspective of matrix integrals not only when there is a high degeneracy of BPS states at zero energy, but also for purely nonBPS supermultiplets (see also section 3.3 and appendix C). In addition, that the growth of ∆ q is quadratic in q for supermultiplets with BPS states will be shown to be universally predicted by random matrix theory. The smallest such ∆ q determines the onset of the continuous density and thus the overall gap ∆ ≡ min q ∆ q for the full spectrum. Generically ∆ > 0, except for the specific anomalous case with δ = 1 / 2, which is the only one that allows for a gapless q = 0 supermultiplet [58, 77, 83]. With this anomaly, the ground state can become completely depopulated if ˆ q = 1; otherwise, there are always supermultiplets with O ( e S 0 ) BPS states. Note that ρ N =2 q ( E -∆ q ) ∼ √ E -∆ q near E = ∆ + q for the generic q = 0 case (cf. eq. (2.24)), whereas for the anomalous q = 0 case ρ N =2 0 ( E ) ∼ 1 / √ E near E = 0 + (cf. eq. (2.25)).  \nFinally, for the N = 4 super-Schwarzian one finds p = -1 and the saddle-point contributions from eq. (2.23) read [58, 60, 84]  \nρ N =4 c ( E ) = e S 0 δ ( E ) + e S 0 (2 π ) 3 / 2 c I 2 ( 2 π √ cE ) E . (2.28)  \nThe N = 4 theory is unique, and its supermultiplets labelled by a half-integer spin J have [43]  \nρ N =4 J ( E ) = e S 0 δ J, 0 + e S 0 J sinh ( 2 π √ E -∆ J ) 2 π 2 E 2 θ ( E -∆ J ) , ∆ J ≡ J 2 . (2.29)  \nIn this case BPS states are all spin zero, and non-BPS states have J ≥ 1 / 2. In contrast with N = 2, the spectrum for N = 4 does always exhibit a positive gap ∆ = 1 / 4, and note that ρ N =4 J ( E -∆ + J ) ∼ √ E -∆ J near E = ∆ + J always in this case [43]. 24 The gapped spectral densities for both N = 2 , 4 give rise to gravitational entropies of the form in eq. (1.13). When gapless, this reduces to the behavior in eq. (1.13) with s = 0.  \nThe intermediate case in eq. (2.25) of a supersymmetric theory that arises for non-near-BPS black holes is physically interesting. In particular, this spectral density diverges as 1 / √ E in the E → 0 + limit, but the gravitational entropy nonetheless becomes negative at large β in consistency with theorem 1. The same happens in the anomalous N = 2 theory that allows for a gapless supermultiplet. This divergence at the ground-state energy may be understood as a supersymmetric transition from a non-degenerate to a highly-degenerate ground state. Intuitively, it captures how the degeneracy of the near-BPS cases gets lifted as soon as supersymmetry is broken, leaving ground-state energies highly populated but unprotected to values above zero. In addition, N = 1 also exhibits a dip between its divergent ground state and its exponential growth at large E , anticipating the formation of a spectral gap with supersymmetry.  \nBesides the obvious fact that eqs. (2.24) to (2.29) have continuous parts, we have thus far said nothing about how these gravitational spectra are associated to spectral ensembles. The realization of these as average spectral densities is more explicitly established by dimensional reduction to a universal AdS 2 sector of the geometry that emerges in the throat of near-extremal black holes. This reduction allows for the definition of 2-dimensional theories of JT gravity where the corresponding Schwarzian governs large diffeomorphisms on the AdS 2 geometry. Formulating JT as a gravitational path integral with a sum over topologies leads to a non-perturbative topological expansion in e -S 0 for observables like those in eq. (1.2). Through topological recursion relations, these have been shown to precisely match double-scaled matrix integrals with expansion parameter e -S 0 to all orders [19, 20, 43, 44]. This result is what leads to the identification of eqs. (2.24) to (2.29) as capturing the statistics of spectral ensembles in the sense of section 2.1.  \nIn the context of JT gravity and random matrices, the above results are all exact perturbatively at large S 0 , but only hold to leading non-perturbative order at large e S 0 . The non-perturbative effects that the genuinely quantum gravitational sum over topologies gives rise to appear suppressed by factors of e -S 0 . These non-perturbative corrections are what makes the duality between JT gravity and matrix ensembles very non-trivial. In this paper, however, we mostly work to leading  \norder at large e S 0 , where the emergence of an ensemble description is not as obvious. Namely, oftentimes one interprets higher-order wormhole contributions to gravitational correlators and their consequent lack of factorization as the root cause for such an ensemble description to arise. It is thus very remarkable that to leading order at large e S 0 , Schwarzian theories alone without wormhole contributions give the spectral densities in eqs. (2.24) to (2.29) all in perfect consistency with universal expectations from random matrix ensembles. Let us briefly comment on these universal features in anticipation of the results in section 3.  \nWith regard to corollary 1, it is worth emphasizing that in random matrix integrals the generic negativity of the annealed entropy is an exact statement, i.e., exact to all perturbative and nonperturbative orders in S 0 . As understood in the previous sections, this negativity is caused by the ensemble generically having a non-trivial distribution of ground-state energies, unlike the nongeneric case in eq. (2.16). Hence the fact that the Schwarzian theories with N = 0 , 1 lead to negative gravitational entropies at low temperatures while the N = 2 , 4 super-Schwarzians do not may be understood as a remarkable consistency check of the ensemble description of effective gravity already to leading non-perturbative order. In particular, in this sense the Schwarzian already knows about the ensemble nature of the gravitational path integral, even though to leading order at large e S 0 we are not even including wormhole contributions. This point of view should be contrasted with the discussion in appendix A, where we review the traditional interpretation of these negativities as a pathology caused by the neglect of non-perturbative corrections that arise at β ∼ O ( e S 0 ). While we agree that such corrections indeed become important at non-perturbatively low temperatures, the random matrix framework says that the qualitative negativity of the gravitational entropy must remain at non-perturbative orders. 25  \nThe consistency between eqs. (2.24) to (2.29) and other well-known features of the spectral statistics of matrix ensembles is also noteworthy. In all cases, the continuous part of the spectrum behaves in one of two characteristic ways near its lower edge. The spectral densities for N = 0, N = 2 with ∆ q > 0, and N = 4, all go to zero as ρ ∼ √ E , whereas those for N = 1, and N = 2 with ∆ q = 0, diverge as ρ ∼ 1 / √ E . These are precisely the two universal behaviors that arise in the analysis of edge statistics in random matrix theory, as we review in section 3. Another general fact about matrix integrals we will address is that a gap in the spectrum arises if and only if the degeneracy of the zero eigenvalue scales with the size of the matrix, and that the size of the gap is quadratic in their ratio. The perfect consistency between this phenomenon and what is observed in eqs. (2.27) and (2.29) is rather remarkable. In both cases, supermultiplets with a positive gap ∆ q > 0 always involve a discrete ground state with a degeneracy that is extensive in e S 0 , while those N = 2 supermultiplets with ∆ q = 0 do not have such a ground state [43, 58]. 26', '3 Random Matrices': 'Our main motivation for studying random matrix ensembles is to understand whether the replica trick in eq. (1.12) can be implemented in effective gravity and, if so, to what perturbative or nonperturbative order in S 0 one has to work to succeed. 27 To pursue this goal in broad generality, we consider integrals over arbitrary matrix ensembles, thus making our predictions apply to any quantum gravity theory with such a dual description. This section reviews the basics of the framework necessary for our discussion, including how Hamiltonians are constructed out of the random matrix in different ensembles, the manifestation of supersymmetry, the matrix scalings relevant to gravity, and the relation between spectral gaps and zero-eigenvalue degeneracy. 28 The application of this technology to the replica trick from eq. (1.12) is the subject of section 4.', '3.1 Probability Densities': 'Let M ⊆ C N × N be some space of square matrices with real eigenvalues. An ensemble of matrices drawn from M can be defined by a PDF on this space, P : M → R ≥ 0 , normalized via  \nˆ M dM P ( M ) = 1 , (3.1)  \nwhere dM is a Lebesgue measure on the algebraically independent entries of M . A standard way of expressing P , up to a normalizing factor that ensures eq. (3.1), is as follows:  \nP ( M ) ∝ e -N Tr V ( M ) , (3.2)  \nwhere V is known as the matrix potential, 29 and the prefactor of N is a conventional scaling (see section 3.6). Writing the diagonal matrix of eigenvalues of M as Λ ≡ diag { λ i } N i =1 , we have  \nTr V ( M ) = Tr V (Λ) = N ∑ i =1 V ( λ i ) . (3.3)  \nOne may understand P ( M ) as a joint PDF for all entries of M under the measure dM . However, we are only interested in the spectral statistics M . To obtain the induced distribution on the space of eigenvalues, one simply has to integrate out all other degrees of freedom. This amounts to computing a Jacobian for whichever transformation diagonalizes M . The upshot is p (Λ), a joint PDF for the N eigenvalues of M . Denoting the Jacobian matrix of the diagonalization of M by J (Λ), this joint PDF against the Lebesgue measure on the eigenvalues may be written  \np (Λ) ≡ 1 Z e -N Tr V (Λ) | det J (Λ) | , Z ≡ ˆ R N d Λ e -N Tr V (Λ) | det J (Λ) | , (3.4)  \n28  \nSee e.g. [85-87] for more details on random matrix theory and its applications to physics relevant to this paper.  \nwhere d Λ ≡ ∏ N i =1 dλ i , and Z is the partition function of the matrix ensemble, which ensures p is normalized to unity. In general, p (Λ) is a full permutation-symmetric function of all eigenvalues.  \nThis is how random matrices realize the general structure of spectral ensembles in section 2.1. In terms of p (Λ), the definition of expectation values ⟨ · ⟩ , n -point functions p n , and in particular the average spectral density ⟨ ρ ⟩ , just parallel eqs. (2.3), (2.5) and (2.7), respectively. For now, however, p (Λ) here captures the spectral statistics of an arbitrary random matrix M , not necessarily the Hamiltonian. To construct physically meaningful Hamiltonians out of M , as we do in section 3.3, we first review in section 3.2 the spaces of matrices M of interest to physics.', '3.2 Ensemble Measures': "Some important ensembles in random matrix theory are those where the space M is diagonalizable under the action of some symmetry group, which also leaves the distribution P invariant. These are referred to as invariant ensembles, and those exhausting the ten discrete symmetry classes of topological invariants are of particular physical relevance to us [88-90]. Three of them are the wellknown Wigner-Dyson (WD) ensembles [91-93], and the other seven will be referred to here as the Altland-Zirnbauer (AZ) ensembles [88, 89, 94]. In contrast with more general matrix constructs, for invariant ensembles symmetry allows one to work out Jacobian measures explicitly.  \nThe three WD ensembles [91-93] are characterized by invariance under the adjoint action of one of the classical Lie groups. They are known as the orthogonal ensemble (OE), unitary ensemble (UE), and symplectic ensemble (SE). Respectively, these describe real symmetric matrices with a distribution invariant under O ( N ), complex Hermitian matrices with a distribution invariant under U ( N ), and quaternionic anti-Hermitian matrices with a distribution invariant under Sp ( N ). 30 For these three WD ensembles, the Jacobian determinant gives the eigenvalue measure  \n| det J (Λ) | WD = | ∆(Λ) | β , (3.5)  \nwhere ∆(Λ) is the Vandermonde determinant  \n∆(Λ) ≡ det { λ i -1 j } N i,j =1 = ∏ i<j ( λ j -λ i ) , (3.6)  \nand the parameter β = 1 , 2 , 4 respectively for the OE, UE, and SE. The other seven AZ ensembles [88, 89, 94] are characterized by additional discrete symmetries: charge conjugation C , spatial reflections R , and time reversal T . 31 The eigenvalue measure these give takes the form  \n| det J (Λ) | AZ = ∏ j<i ∣ ∣ λ 2 i -λ 2 j ∣ ∣ β ∏ k | λ k | α . (3.7)  \nThe specific values of the pair ( α , β ) for which eq. (3.7) descends from a random matrix integral with specific discrete symmetries can be found in [96]. 32 Since eq. (3.7) is invariant under sign changes, the canonical integral over eigenvalues can be taken to be non-negative and eq. (3.7) turned into a measure over squared eigenvalues ω i ≡ λ 2 i for Ω ≡ Λ 2 ,  \n| det J (Ω) | AZ = 2 -N | ∆(Ω) | β ∏ k | ω k | α -1 2 . (3.8)  \nAlthough for AZ ensembles with α = 1 this reduces to eq. (3.5), one should not be misguided into thinking that such AZ ensembles will give equivalent results to WD ensembles. In the latter we are interested in the eigenvalues of the random matrix, which lie unconstrained on R , whereas for the former eq. (3.8) only arises for squared eigenvalues, which are non-negative (cf. appendix D.5).  \nSome of these AZ ensembles admit invariant generalizations of much interest. In four of them, M is a rank-2 tensor of some type, with the continuous symmetry G ( N ) being one of the classical Lie groups. The other three AZ ensembles describe a random N × N matrix M which transforms as a bifundamental of the direct product G ( N ) × G ( N ), i.e., each index of M is separately acted upon by one copy of G ( N ) in the fundamental representation.  \nThe latter three can very naturally be generalized to ensembles of rectangular matrices by considering products of symmetry groups of different rank. In particular, letting the symmetry group be G ( N ' ) × G ( N ) allows one to describe a rectangular random matrix M of size N ' × N . As will be seen shortly in section 3.3, this construction precisely reproduces the requisite structure to describe the supermultiplet spectra from eqs. (2.27) and (2.29), with the number of BPS states given by ¯ ν ≡ | N ' -N | . 33 For this reason, the generalization of standard AZ ensembles to rectangular matrices with ¯ ν > 0 will be referred to here as BPS ensembles.  \nThe spectrum of a rectangular matrix M is given by its singular values, which by definition are real and non-negative for any complex matrix. For the construction above, singular value decomposition gives a measure of the AZ form in eq. (3.7), but with the replacement [43]  \nα → α + β ¯ ν = β (¯ ν +1) -1 , (3.9)  \nand β specified as usual by the corresponding simple group. On the right, we used that α = β -1 for the AZ ensembles with bifundamental symmetry that generalize to BPS ensembles.  \nWithout loss of generality, let the rectangular matrix M in our BPS ensembles be ( N +¯ ν ) × N with ¯ ν > 0. As a random matrix, M has rank N with unit probability, so consider the generic case in which the N singular values of M are all positive. Using M , one can construct positive semidefinite square matrices MM † and M † M of sizes N +¯ ν and N , respectively. Here M † denotes the Hermitian conjugate, transpose, or quaternionic conjugate of M , depending respectively on whether one is interested in unitary, orthogonal, or symplectic ensembles. In either case, the rank of both of the resulting square matrices is just N as inherited from M , implying that only M † M is  \ngenerically full-rank, whereas MM † is always singular and deterministically has ¯ ν zero eigenvalues. Despite this qualitative difference, the positive eigenvalues of both MM † and M † M are simply the squares of the singular values of M , and thus identically distributed. In particular, their measure is the same and given by eq. (3.8) with the appropriate value of β and α as specified by eq. (3.9). When addressing black hole spectra, this fact will help explain why the onset of the continuous part of the spectrum in eqs. (2.27) and (2.29) is separated from zero by a positive gap, not only for supermultiplets with BPS states, but also for purely non-BPS ones.  \nBecause the measures of all matrix ensembles of interest end up looking like eqs. (3.5) and (3.8), we will henceforth take the Jacobian determinant to have the following general form: 34  \n| det J (Λ) | = | ∆(Λ) | β ∏ k | λ k | ν . (3.10)  \nThis reproduces all invariant ensemble measures of interest as follows:  \n       WD: ν = 0 , AZ: ν = α -1 2 , BPS: ν = β 2 (¯ ν +1) -1 , (3.11)  \nwhere ( α , β ) are the standard parameters determined respectively by the discrete and continuous symmetry groups of the corresponding WD or AZ ensemble, and ¯ ν is the parameter specifying the difference in rank between groups in the bifundamental of the corresponding BPS ensemble. Despite all WD, AZ, and BPS ensembles being accounted for by eq. (3.10), their spectral statistics end up being qualitatively very different and rich in structure.  \nLet us conclude this section by noting that the construction of square matrices above in terms of a random rectangular matrix M is a well-studied one in random matrix theory. In particular, the full-rank M † M defines what is traditionally known as a Wishart matrix [99, 100], and MM † is often referred to as a singular Wishart matrix [101]. For this reason, BPS ensembles may be understood as generalizing certain AZ ensembles by a Wishart construction. Note, however, that the literature on Wishart ensembles usually defines M by making its entries follow Gaussian distributions. In contrast, our BPS ensembles allow for arbitrary non-Gaussian matrix potentials.", '3.3 Random Hamiltonians': 'The physical motivation for considering matrix ensembles is to study the spectrum of random Hamiltonians. Non-supersymmetric theories can be described by WD ensembles where the random matrix M itself is identified as the Hamiltonian,  \nH ≡ M. (3.12)  \nCorrespondingly, the eigenvalue measure for such Hamiltonians is given by eq. (3.10) with ν = 0.  \nSupersymmetric theories require the richer structure of AZ and BPS ensembles. To make the discussion of supersymmetric Hamiltonians as general as possible, it is thus convenient to allow for  \nthe random matrix M to be ( N +¯ ν ) × N for any integer ¯ ν ≥ 0. This way, both the square matrices of AZ ensembles (¯ ν = 0) and the rectangular matrices of BPS ensembles (¯ ν > 0) are accounted for. The algebraic structure of supersymmetry makes supercharges the natural objects to be actually identified with M . As usual in supersymmetric theories, the Hamiltonian is then constructed from anti-commutators of conjugate supercharges which here involve M and M † matrices. The resulting random Hamiltonian naturally consists of the Wishart combinations MM † and M † M described in the previous section. More specifically, the Hamiltonian generally decomposes into statistically independent supermultiplets, each described by its own ensemble. Hence the overall spectrum of the theory can be obtained by combining Hamiltonians acting on supermultiplets of the form  \nH = ( MM † 0 0 M † M ) . (3.13)  \nA detailed discussion of how this structure arises in supersymmetric random matrix theory is provided in appendix C. For both of the positive semidefinite sub-matrices in eq. (3.13), the spectrum of positive eigenvalues is governed by the measure in eq. (3.10) with the corresponding AZ and BPS parameters given by eq. (3.11). The difference between the two types of ensembles is that only the latter are able to actually make MM † singular with a zero eigenvalue of arbitrarily high degeneracy ¯ ν . Interestingly, this degenerate zero eigenvalue does not always capture zeroenergy BPS states. As we show in appendix C and explain below, the parameter ¯ ν acquires two possible meanings in supersymmetric random matrix theory.  \nWith N = 1 supersymmetry, eq. (3.13) captures the full Hamiltonian of the theory and thus determine its spectrum. As a result, in this case the ¯ ν zero eigenvalues that MM † has genuinely correspond to BPS ground states. This is the random matrix mechanism behind the non-generic spectral ensembles alluded to in section 2.2 and giving rise to a discrete ground state in eq. (2.16). Indeed, this rectangular matrix construction results in a singularly distributed ground-state energy which in the spectral density becomes ¯ ν δ ( E ) (cf. the near-BPS spectra in section 2.4). The other N positive eigenvalues of MM † are equal to those of M † M , thus giving a 2-fold degenerate spectrum of non-BPS states. This degeneracy is the result of pairs of non-BPS states of the theory forming irreducible doublets. In other words, the 2 N positive eigenvalues correspond to a (reducible) supermultiplet of N doublets consisting of pairs of non-BPS states with the same energy (see appendix C for more details). One may think of the full Hilbert space of non-BPS states of the N = 1 as forming a single supermultiplet.  \nWith N = 2 supersymmetry and higher, the situation is different. In particular, the supersymmetry algebra now has a non-trivial R-symmetry which allows one to organize the Hilbert space into different representations. Let us for simplicity here assume that states have R-charge k ∈ Z and the supercharges Q and Q † have ± 1 R-charge. 35 Then the Hilbert space reads H = ⊕ k H k , where H k is the subspace of states with R-charge k . The restriction of the supercharges to these subspaces gives Q k : H k → H k +1 and Q † k : H k +1 → H k . This way, the supercharges decompose into Q = ∑ k Q k and Q † = ∑ k Q † k , and the Hamiltonian H = { Q,Q † } becomes H = ∑ k H k with Q k Q † k + Q † k Q k . As it turns out, every Q k describes the spectrum of a statistically independent supermultiplet involving states with R-charge values ( k, k +1). What eq. (3.13) captures in this  \nsetting is the spectrum of a single supermultiplet, not of the full Hamiltonian. In other words, every H k term in H is described by a different ensemble realization of eq. (3.13) with the supercharge operator Q k identified as Q k ∼ M in terms of the random matrix M .  \n̸  \nThe important consequence this has on the meaning of ¯ ν can be easily described. Consider some state ψ k +1 ∈ H k +1 . 36 If ψ k +1 is a BPS state, then it is annihilated by both supercharges and, in particular, Q † k ψ k +1 = 0. In terms of the random matrix M ∼ Q k , this means that ψ k is in the kernel of MM † and thus contributes to a positive value ¯ ν > 0. This is the familiar case we already encountered where ¯ ν is naturally associated to the number of BPS states. However, suppose now that ψ k +1 is a non-BPS state, and recall that supercharges are exact operators satisfying Q 2 = Q † 2 = 0. This supersymmetry condition implies that ψ k +1 must be annihilated by at least one supercharge, 37 while the assumption that it is non-BPS means that ψ k must be annihilated by only one supercharge. Out of the two cases, consider first Qψ k +1 = 0, which in particular means Q † k ψ k +1 = 0. For M ∼ Q k this implies ψ k +1 is not in the kernel of MM † , and thus we learn nothing about ¯ ν and its meaning. Hence consider now the case in which Q † ψ k +1 = 0, such that we have Q † k ψ k +1 = 0. Then M ∼ Q k actually does make the kernel of MM † non-trivial and ¯ ν > 0. In other words, here we see the non-BPS state ψ k +1 contributing to ¯ ν simply because it is annihilated by Q k , even though it is not BPS. 38 Indeed, ψ k +1 gives a zero eigenvalue for the supermultiplet Hamiltonian H k , but it would give a non-zero eigenvalue for H k +1 . The matrix ensemble description of H k hence feels a zero eigenvalue from ψ k +1 , even though there is no zeroenergy BPS states associated to it in the theory.  \nIn summary, we see that in supersymmetric random matrix theory eq. (3.13) generally captures the spectrum of single supermultiplets, and that its ¯ ν zero eigenvalues may or may not be associated to BPS states. In the thorougher discussion from appendix C, we find that ¯ ν is always associated to either one or the other. In particular, when there are BPS states associated to the supermultiplet R-charges, ¯ ν is precisely the number of such BPS states. Otherwise, ¯ ν can generally be associated to the varying number of states in non-BPS supermultiplets. More explicitly, if N + k is the number of non-BPS states in ( k, k + 1) doublets and say N + k -1 < N + k < N + k +1 , then the random matrix description for Q k ∼ M turns out to involve ¯ ν = N + k -1 .  \nRegarding the black hole spectra from section 2.4, the relation to random matrices is already rather manifest. The random Hamiltonians of WD ensembles describe non-supersymmetric theories like the bosonic Schwarzian theory. The seven standard AZ ensembles can describe supersymmetric Hamiltonians of theories with no unbroken supersymmetry, and thus are particularly suited for the N = 1 super-Schwarzian theory, or the specially anomalous q = 0 gapless supermultiplet of the N = 2 super-Schwarzian theory with δ = 1 / 2 and ˆ q = 1. Finally, the richer Wishart matrix structure of BPS ensembles allows to describe the more general cases of the N = 2 , 4 super-  \n̸  \n̸  \nSchwarzian theories, both for supermultiplets with and without BPS states. In particular, the ¯ ν parameter of the latter will turn out to be crucial to explain not only the spectral gap when there are BPS states, but also for purely non-BPS sectors where the onset of the continuous spectrum continues to rise to higher energies for higher supermultiplets (see section 3.6.2).', '3.4 Equilibrium Conditions': "It is often convenient to raise the Jacobian determinant in eq. (3.4) into the exponential so as to write the matrix partition function as an action integral,  \nZ = ˆ R N d Λ e -NI (Λ) , I (Λ) ≡ Tr ( V (Λ) -1 N log | J (Λ) | ) . (3.14)  \nwhere we have used the matrix identity log det M = Trlog M to rewrite the J (Λ) contribution. Correspondingly, the joint eigenvalue PDF in terms of this action reads  \np (Λ) = e -NI (Λ) Z . (3.15)  \nFor the general Jacobian in eq. (3.10), the matrix action can be written out to be  \nI (Λ) = N ∑ i =1 V ν ( λ i ) -β N N ∑ i,j =1 j>i log | λ j -λ i | , (3.16)  \nwhere we have absorbed the ν -dependent term into a redefinition of the matrix potential to  \nV ν ( λ ) ≡ V ( λ ) -ν N log | λ | . (3.17)  \nMinima of I (Λ) are known as Fekete points and obey the equilibrium conditions  \n∂I (Λ) ∂λ k ∣ ∣ ∣ ∣ Λ=Λ ∗ = 0 , ∀ k = 1 , . . . , N. (3.18)  \nThese points, which dominate the matrix integral at large N , motivate the notion of an effective potential characterizing the interaction of each eigenvalue with all others. Such an effective potential ̂ V for the eigenvalue λ k can be easily constructed by simply keeping those terms in I (Λ) which depend on λ k . Applying this to eq. (3.16),  \n̂ V Λ ( λ k ) ≡ V ν ( λ k ) -β N N ∑ k = i =1 log | λ k -λ i | , (3.19)  \n̸  \nwhere the Λ subscript emphasizes the dependence of ̂ V Λ on all eigenvalues. Upon differentiation,  \n̂ V ' Λ ( λ k ) = V ' ν ( λ k ) -β N N ∑ k = i =1 1 λ k -λ i . (3.20)  \n̸  \nBy construction, this gives the same equilibrium conditions as eq. (3.18) and thus  \n̂ V ' Λ ( λ k ) ∣ ∣ ∣ Λ=Λ ∗ = 0 , ∀ k = 1 , . . . , N, (3.21)  \nat Fekete points. A crucial property of the effective potential is that it exists even when the action does not, as happens in section 3.6.4 for the double scaling limit relevant to gravity. 39  \nOf course, the system of N saddle-point equations that eq. (3.21) defines for a general matrix ensemble is highly non-trivial. A standard procedure in random matrix theory to evaluate such largeN matrix integrals, whether it is at a saddle-point level or not, is to first take a continuum limit. In the continuum, the Fekete conditions in eq. (3.21) for effective potentials like eq. (3.19) become well-studied differential equations which can be solved exactly.", '3.5 Continuum Formalism': "The continuum limit can be easily obtained by first turning summations into integrals via  \nN ∑ k =1 ( · ) = N ˆ R dx ˆ ρ Λ ( x ) ( · ) , ˆ ρ Λ ( x ) ≡ 1 N N ∑ k =1 δ ( x -λ k ) . (3.22)  \nIn the continuum, the discrete measure ˆ ρ Λ is replaced by a continuous measure ˆ ρ , in terms of which the discrete λ k eigenvalues are now captured by a continuous x variable. Doing so, ˆ ρ takes on the dynamical role of the variables Λ, thereby turning the matrix integral from an eigenvalue integral to a functional integral over unit-normalized measures,  \nˆ R N d Λ -→ N ˆ D ˆ ρ. (3.23)  \nApplying this to eq. (3.14) leads to the continuum limit  \nZ = N ˆ D ˆ ρe -N 2 ̂ I [ˆ ρ ] , (3.24)  \nwhere the action I (Λ) in eq. (3.16) has been replaced by a functional N ̂ I [ ˆ ρ ] given by  \n̂ I [ ˆ ρ ] ≡ ˆ R dx ˆ ρ ( x ) V ν ( x ) -β 2 ˆ R dxdy ˆ ρ ( x ) ˆ ρ ( y ) log | x -y | . (3.25)  \nSimilarly, in the continuum the effective potential in eq. (3.19) becomes  \n̂ V [ ˆ ρ ; x ] = V ν ( x ) -β ˆ R dy ˆ ρ ( y ) log | x -y | , (3.26)  \nwhere the notation emphasizes that the effective potential ̂ V itself depends not only on the eigenvalue x , but also on all the other eigenvalues through the ˆ ρ measure (cf. ̂ V Λ in the discrete case). Differentiating eq. (3.26) we also obtain the continuum version of eq. (3.20),  \n∂ x ̂ V [ ˆ ρ ; x ] = V ' ν ( x ) -β R dy ˆ ρ ( y ) x -y , (3.27)  \n̸  \nwhere ffl denotes a principal value integral. More abstractly, in general ̂ V is defined in the continuum by functional variation as  \n̂ V [ ˆ ρ ; x ] ≡ δ ̂ I [ ˆ ρ ] δ ˆ ρ ( x ) . (3.28)  \nThe extremization of the action for a saddle-point approximation of the integral in eq. (3.24) thus clearly involves the effective potential. Since such an extremization is over measures satisfying the normalization condition ' R dx ˆ ρ ( x ) = 1, the variational problem is better solved by introducing a Lagrange multiplier for this constraint. The upshot is that extremality of ̂ I [ ˆ ρ ] against functional variations requires ̂ V to be constant along the support of ˆ ρ or, equivalently,  \n∂ x ̂ V [ ˆ ρ ∗ ; x ] = 0 , ∀ x ∈ supp ˆ ρ ∗ , (3.29)  \nfor ˆ ρ = ˆ ρ ∗ a unit-normalized saddle point of ̂ I [ ˆ ρ ]. This is the continuum form of the Fekete conditions in eq. (3.21). The equilibrium measure ˆ ρ ∗ that eq. (3.29) defines is recognized as the leading average spectral density of the matrix ensemble at large N ,  \nˆ ρ ∗ ( x ) ≡ lim N →∞ ⟨ ρ ( x ) ⟩ N . (3.30)  \nThe equilibrium condition in eq. (3.29) can be made more explicit using eq. (3.27),  \nR dy ˆ ρ ∗ ( y ) x -y = V ' ν ( x ) β , ∀ x ∈ supp ˆ ρ ∗ . (3.31)  \nObtaining the leading spectral density thus amounts to finding a function ˆ ρ ∗ that satisfies eq. (3.31). This is a type of singular integral equation whose solutions remarkably admit a closed-form expression in terms of a general matrix potential. The derivation of such a general solution for ˆ ρ ∗ is reproduced in appendix D.3. The final result is  \nˆ ρ ∗ ( x ) = + √ η ( x ) β π 2 Σ dy y -x V ' ν ( y ) + √ η ( y ) , x ∈ Σ , (3.32)  \nwhere Σ = supp ˆ ρ ∗ is specified by eq. (D.17) and the conditions in eqs. (D.29) and (D.30), and the function η is given by eqs. (D.18) and (D.19) for WD and AZ ensembles, respectively. The notation + √ · instructs one to approach the square-root branch cut from above as in eq. (D.26).  \nHaving ˆ ρ ∗ , it is useful to define the leading effective potential that appears in eq. (3.29) by  \n̂ V ∗ ( x ) ≡ ̂ V [ ˆ ρ ∗ ; x ] . (3.33)  \nBesides vanishing along supp ˆ ρ ∗ by eq. (3.29), the gradient of ̂ V ∗ also defines an important object upon analytic continuation to C off supp ˆ ρ ∗ known as the spectral curve,  \ny ( z ) = -̂ V ' ∗ ( z ) β , z ∈ C ∖ supp ˆ ρ ∗ . (3.34)  \nAs we show in appendix D.2, y has important analytic properties and is intimately related to ˆ ρ ∗ itself. In particular, denoting the analytic continuation of ˆ ρ ∗ off its support to C by ˆ ρ c ∗ ,  \ny ( z ) = -iπ ˆ ρ c ∗ ( z ) , z ∈ C ∖ supp ˆ ρ ∗ . (3.35)  \nFurthermore, approaching the branch cut of y along supp ˆ ρ ∗ from opposite sides leads to  \nˆ ρ ∗ ( x ) = ± i π lim ϵ → 0 + y ( x ± iϵ ) , x ∈ supp ˆ ρ ∗ . (3.36)  \nIn words, eq. (3.35) says that y is fully determined by ˆ ρ ∗ , while eq. (3.36) says that ˆ ρ ∗ is fully determined by y . The fact that both objects contain equivalent information about the spectral ensemble will turn out to be very useful for us.", '3.6 Matrix Scalings': 'Matrix scalings refer to largeN limits which amplify certain features of the spectral statistics of a matrix ensemble. These scalings may be single N →∞ limits that zoom into the bulk [102] or the edge [103] of the spectrum, or double limits in which the parameters of the matrix integral are simultaneously tuned towards critical values as N →∞ [104-106].  \nAs part of these procedures, one generally introduces a zooming variable ¯ x that replaces the continuum matrix eigenvalue variable x . This is done by setting x ≡ f N (¯ x ) and keeping ¯ x ∼ O (1) as N →∞ , where the choice of scaling function f N determines which spectral features the zooming variable captures in the limit. The resulting construct describes a spectral ensemble in its own right, capturing the statistics of the effective eigenvalue variable ¯ x in the original matrix ensemble at large N . Certain matrix ensembles are actually more naturally defined in their scaled version and suitable parameters therein, with no allusion to finite matrices.  \nWhen constructing Hamiltonians out of random matrices as in section 3.3, we are faced with a choice to either interpret x or ¯ x as the random variable associated to eigenvalues of the Hamiltonian. If the Hamiltonian is naturally identified with the random matrix before any scaling, then x is the desired variable. For instance, if the Hamiltonian is the random matrix itself, a Boltzmann factor would read e -βx . However, if we are interested in the statistics of a specific zooming variable, then applying the scaling after would lead to e -βf N (¯ x ) . This would defeat the purpose, since such a Boltzmann factor would remain sensitive to the statistics of x rather than of ¯ x in the limit. Instead, we will be addressing situations where the Hamiltonian spectrum is naturally associated to whichever zooming variables are relevant to the problem, and thus particularly sensitive to the spectral statistics that these highlight. In other words, we take scaled matrix ensembles as the starting point when constructing spectral ensembles, and take ¯ x as the only eigenvalue variable that matters. This way, if e.g. the Hamiltonian is the random matrix itself, a Boltzmann factor would read e -β ¯ x and thus successfully capture the statistics of ¯ x .  \nThe above implies that in general there will be no simple mapping between quantities involving Hamiltonians for matrix ensembles before and after scaling, or between different scalings. Namely, Hamiltonian spectra defined in terms of the inequivalent zooming variables x and ¯ x will give rise to different largeN physics which we will not be interested in relating. Since matrix integrals in different scaling limits will generally be treated as different constructs, it will be simpler to think of matrix scalings as replacements of the form x → f N ( x ) before taking the N →∞ limit. This way it will not be necessary to deal with the notation ¯ x for zooming variables, and we can unambiguously use x to refer to whichever eigenvalue variable remains O (1) in the largeN limit. Let us now describe in more detail the matrix scalings relevant to this paper.', '3.6.1 Single Limits': "We begin by addressing single limits in which all parameters of the matrix integral are kept fixed as N →∞ . Already the factor of N in eq. (3.2) was introduced with hindsight as a form of matrix scaling. Namely, this factor guarantees that for potentials with a quadratic term the typical smallest and largest eigenvalues are of O (1) size in x . In other words, had we not introduced this factor of N , the size of these eigenvalues would be extensive and O ( N 1 / 2 ). Keeping x ∼ O (1) at large N guarantees that the leading spectral density given by eq. (3.30) converges to a bounded support of finite size in x in the limit. For instance, this largeN limit for Gaussian WD matrices leads to the well-known Wigner semicircle law,  \nˆ ρ ∗ ( x ) GWD = 1 π β √ 2 β -x 2 , -2 β < x < 2 β , (3.37)  \nwhich is illustrated in fig. 2. Given that there are N eigenvalues in total, a consequence of such a scaling is that in the x variable the typical spacing between eigenvalues is O ( N -1 ). Since this goes to zero at large N , the variable x is clearly not suitable for capturing eigenvalue fluctuations. A better zooming variable would be obtained by replacing x → x/N , in which case the typical separation between eigenvalues would be expected to remain O (1) at large N . This heuristic argument gives the provably correct scaling one should perform in order to study the statistics of eigenvalues in the bulk of the spectrum.  \nWhen studying random Hamiltonians, such a bulk scaling would focus on the statistics of the energies of generic excited states. However, by the observations in section 2, we are actually more interested in understanding the spectral statistics of energies near the ground state. This corresponds to studying the tail of eigenvalues near the lower edge inf supp ˆ ρ ∗ of the leading spectrum. Interestingly, as it turns out the eigenvalue fluctuations near the edge are often larger than O ( N -1 ). For e.g. the WD ensembles the typical fluctuations around the x = -√ 2 β edge are O ( N -2 / 3 ), and would thus blow up in the largeN limit that focuses on bulk statistics. As a result, studying edge statistics generally requires both a precise shift to the desired edge and a subtler scaling in N (see appendix D.4 for explicit examples).  \nEdge scalings generally result in non-normalizable spectral functions in the N → ∞ limit. When zooming on spectral edges it is thus useful to introduce an eigenvalue density parameter e S 0 that effectively replaces the diverging total eigenvalue number N . This can be done by simply using Ne -S 0 instead of just N in the appropriate zooming variable. For the GUE this procedure defines the Airy matrix model, whose exact spectral density in e S 0 is [107] 40  \nˆ ρ Airy ( x ) = e -S 0 / 3 ( Ai ' ( -e 2 S 0 / 3 x ) 2 + e 2 S 0 / 3 x Ai ( -e 2 S 0 / 3 x ) 2 ) , x ∈ R . (3.38)  \nAt large e S 0 one obtains the leading spectral density of the Airy model,  \nˆ ρ Airy ∗ ( x ) ≡ lim S 0 →∞ ˆ ρ Airy ( x ) = √ x π , x > 0 , (3.39)  \nFigure 2 . Illustration of spectral densities for the Gaussian WD ensemble with β = 2, i.e., the GUE. Left: scaling to a compact distribution, with the N = 10 result (red) obtained via the orthogonal polynomials techniques from appendix D.4, and the N → ∞ limit (black) corresponding to the Wigner semicircle law from eq. (3.37). Right: scaling to the lower edge, with the exact N → ∞ limit (red) given by the Airy model from eq. (3.38), and its leading soft edge (black) given by eq. (3.39).  \n<!-- image -->  \nwhich is illustrated against eq. (3.38) in fig. 2. This precisely matches the behavior of the leading spectral density of the GUE, given by eq. (3.37) for β = 2, near the x = -2 edge of the spectrum.  \nThese scalings for Gaussian potentials are presented in more detail in appendix D.4, for all WD, AZ, and BPS ensembles. For the AZ case, the counterpart of eq. (3.37) is the singular Marˇcenko-Pastur law [108],  \nˆ ρ ∗ ( x ) GAZ = 1 2 π β √ 4 β -x x , 0 < x < 4 β , (3.40)  \nwhich diverges as x → 0 + (see fig. 3). This is in stark contrast with the regularity near the lower edge of the spectrum in the WD case. Spectral edges where the spectral density goes to zero (or some finite value) like in eq. (3.37) are often referred to as soft edges, and those where it diverges like in eq. (3.40) as hard edges. As we explain in section 3.6.3, this turns out to give a universal characterization of the edge statistics of spectral ensembles.  \nNotice how eq. (3.40) in this single largeN limit has no dependence on the ν parameter, even though it is typically nonzero for AZ ensembles. In the standard AZ ensembles ν takes fixed O (1) values and the leading result in eq. (3.40) is unavoidable; in the BPS ensembles, however, an interestingly different story arises. The BPS case is explained in section 3.6.2.  \nThe scaling that zooms onto the x = 0 edge of the spectrum in eq. (3.40) for β = 2 defines a family of Bessel matrix models, whose exact spectral density in e S 0 is [109]  \nˆ ρ Bessel ν ( x ) = e S 0 2 ( J ν ( e S 0 √ x ) 2 -J ν +1 ( e S 0 √ x ) J ν -1 ( e S 0 √ x )) , x > 0 . (3.41)  \nAt large e S 0 , this reproduces the hard edge of the leading spectral density in eq. (3.40),  \nˆ ρ Bessel ∗ ( x ) = lim S 0 →∞ ˆ ρ Bessel ν ( x ) = 1 π √ x , x > 0 . (3.42)  \nNote, however, that this singular behavior only arises in the limit of eq. (3.42); in fact, eq. (3.41) gives a perfectly regular spectral density as x → 0 + for any ν ≥ 0. It is only for ν = -1 / 2, a  \nFigure 3 . Illustration of spectral densities for the Gaussian AZ ensemble with ( α , β ) = (2 , 2). Left: scaling to a compact distribution, with the N = 10 result (red) obtained via the orthogonal polynomials techniques from appendix D.4, and the N →∞ limit (black) corresponding to the singular Marˇcenko-Pastur law from eq. (3.40). Right: scaling to the lower edge, with the exact N →∞ limit (red) given by the Bessel model from eq. (3.41), and its leading hard edge (black) given by eq. (3.42).  \n<!-- image -->  \nspecific case allowed by eq. (3.11) for AZ ensembles with α = 0, that eq. (3.41) actually is singular. Explicitly, near x = 0 one finds ˆ ρ Bessel ν ( x ) ∼ x ν , which vanishes in the x → 0 + limit for ν > 0, gives the finite value e S 0 / 2 for ν = 0, and diverges as 1 / √ x for ν = -1 / 2. In other words, the exact spectral density of the Bessel model generically does not actually diverge as x → 0 + , despite the hard edge that the leading result in eq. (3.42) exhibits. The behavior ˆ ρ Bessel ν ( x ) ∼ x ν actually suggests that the growth of these models off x = 0 can be made arbitrarily slow by increasing the ν parameter, which is certainly possible in BPS ensembles. Indeed, as we describe next, ν turns out to soften the hard edge of these ensembles by pushing their leading spectrum off x = 0 and creating a spectral gap ∆ > 0.", '3.6.2 Spectral Gaps': 'Recall that by eq. (3.11) in BPS ensembles ν is linearly related to the free parameter ¯ ν ∈ N which determines the size of the rectangular ( N +¯ ν ) × N random matrix M . A natural largeN scaling of rectangular matrices would preserve their aspect ratio ( N +¯ ν ) /N , which requires ¯ ν to grow linear in N . For the Wishart matrices constructed out of M as described in section 3.2, this is indeed the standard scaling in the random matrix literature [100, 110]. Letting ¯ ν ∼ O ( N ) makes ν and ¯ ν coincide up to factors of 2 from β to leading order at large N . It is thus of interest for BPS ensembles to generally consider a scaling where ν = νN with ν ∼ O (1). Given that ν ∼ ¯ ν at large N , for these ensembles we may at times refer to ν itself as the degeneracy of the zero eigenvalue.  \nWith this scaling, the potential V ν from eq. (3.17) does not trivialize to just V at large N , but preserves a logarithmic pole at zero in the continuum, 41  \nV ν ( x ) = V ( x ) -ν log | x | . (3.43)  \nRemarkably, this turns eq. (3.40) into a regular Marˇcenko-Pastur law,  \nˆ ρ ∗ ( x ) GBPS = √ ( a + -x )( x -a -) 2 π β x , a -< x < a + , (3.44)  \nρ  \nFigure 4 . Illustration of spectral densities for the Gaussian BPS ensemble with β = 2 and ν = 1 2 + √ 2 such that eq. (3.46) gives gap size ∆ = 1. Left: scaling to a compact distribution, with the N = 10 result (red) obtained via the orthogonal polynomials techniques from appendix D.4, and the N →∞ limit (black) corresponding to the regular Marˇcenko-Pastur law from eq. (3.44). Right: scaling to the lower edge, with the exact N = 50 result (red) approaching the Airy model from eq. (3.38) as N →∞ , and its leading soft edge (black) approaching eq. (3.39) as N →∞ (cf. the strict Airy limit in fig. 2).  \n<!-- image -->  \nwhere the endpoints of the support interval are given by  \na ± = β ( 1 ± √ 1 + 2 ν β ) 2 . (3.45)  \nThis is illustrated in fig. 4. Since eq. (3.40) is recovered in the limit ν → 0, we may generally refer to eq. (3.44) also when addressing standard AZ ensembles. In contrast with eq. (3.40), we see that eq. (3.44) is now non-singular and has no support on an open interval (0 , a -) above zero for ν > 0. Since zero is the smallest eigenvalue possible for the positive semidefinite random Hamiltonians in these ensembles, a nonempty interval of this kind constitutes a spectral gap. This gap in the spectrum separates the discrete ground state with zero energy from all other higher-energy states. Its size is  \n∆ ≡ a -=    ν 2 β + O ( ν 3 ) , ν ≪ 1 2 ν + O (1) , ν ≫ 1 (3.46)  \nwhere we have written out the scaling of ∆ at small and large ν = ν /N . In particular, note that the gap is never extensive in N , but determined by the ratio between the degeneracy of the zero eigenvalue and the number of excited states. The gap can be numerically seen to be populated by O (1) eigenvalues in total as N → ∞ , thus becoming relatively depleted in the limit. When the relative number of zero eigenvalues ν → 0 at large N , corresponding to a non-extensive degeneracy, the spectrum becomes gapless.  \nIn supersymmetric terms, this means that there arises a gap between the zero-energy ground state and the energy of the lightest non-BPS states if and only if ¯ ν is extensive in N . In supermultiplets where ¯ ν is the number of BPS states, this explains the origin of the spectral gap and predicts its size in terms of the BPS ground-state degeneracy. 42 In consistency with our observations here, we see that the N = 2 , 4 super-Schwarzians that arise in near-BPS black holes indeed exhibit an extensive degeneracy of BPS states given by e S 0 when their spectrum is gapped. In particular, our  \nresults precisely reproduce the behavior of N = 2 supermultiplets in eq. (2.27): at small q ≪ ˆ q , the BPS degeneracy gives a ratio ¯ ν /e S 0 ≈ q/ ˆ q , and correspondingly the gap is seen to be quadratic in q/ ˆ q (the identification N ∼ e S 0 will be explained shortly).  \nRemarkably though, ¯ ν clearly must also be playing an important role even in purely nonBPS supermultiplets. Namely, the onset of the continuum in eqs. (2.27) and (2.29) keeps growing for higher supermultiplets, even though these are unrelated to BPS states. This suggests that we should match the quadratic growth of ∆ q (or ∆ J ) in the R-charge q (or J ) to the large-gap limit of eq. (3.46). This gives the identification ν ∼ q 2 asymptotically, and thus ¯ ν ∼ q 2 e S 0 . The interpretation of this extensive value of ¯ ν requires in this non-BPS case following the discussion in section 3.3 and results of appendix C. In particular, the quadratic growth of ¯ ν in q must be related to a similar growth in the number of non-BPS states in higherq supermultiplets. Letting L q e S 0 be the number of non-BPS states in a q -multiplet, explicitly we would expect L q ∼ q 2 (see comments before eq. (C.9)).  \nAnother important consequence of the ν ∼ O ( N ) scaling is that now the leading spectral density in eq. (3.44) behaves regularly near both edges of the spectrum. In particular, eq. (3.44) has a soft edge at x = ∆ where ˆ ρ ∗ just goes to zero. The hard edge of the AZ case in eq. (3.40) can be recovered in the ν → 0 limit. More generally, in terms of the original ν parameter, the hard edge arises if ν ∼ O (1) or potentially if its growth with N is slower than linear in the largeN limit. 43 Summarizing, the BPS ensembles give rise to two qualitatively distinct scenarios: a gapless spectrum with a hard edge if ν ∼ O (1), and gapped spectrum with a soft edge if ν ∼ O ( N ),  \nThe scaling limit that zooms onto the hard edge in the gapless case again gives rise to the Bessel model in eq. (3.41). The gapped case requires first a shift of the zooming variable to center x = ∆, and an inequivalent scaling in N to capture the statistics of the soft edge. Remarkably though, the appropriate scaling here turns out to be of the same form as in the WD ensembles, and yields the same Airy model from eq. (3.38) in the largeN limit [113], 44 as illustrated in fig. 4. As a result, for BPS ensembles with a positive gap, the spectral statistics of the lower edge are in fact not governed by the Bessel models in eq. (3.41), but by the Airy model in eq. (3.38). The surprising return of the BPS ensemble to Airy statistics is not a coincidence, but explainable by universality results in random matrix theory.', '3.6.3 Edge Universality': 'Scaling limits were described in previous sections referencing eqs. (3.37), (3.40) and (3.44), which are leading spectral densities for our invariant ensembles whose precise form assumes a Gaussian potential. In fact, universal properties of random matrix statistics imply that upon edge scalings, the results in eqs. (3.38) and (3.41) actually arise much more broadly for general matrix ensembles and potentials. 45 Indeed, universality results hold beyond just invariant ensembles, and take the  \nform of local spectral statistics which in specific scaling limits are fully determined by the eigenvalue measure and independent of details of the matrix potential [117].  \nUniversality applies to bulk statistics [102, 118] and to edge statistics [103, 116, 119, 120]. Given our interest in understanding low-energy fluctuations, we are in particular interested in edge universality. Our spectral ensembles give rise to edge statistics that fall into two universality classes: Airy and Bessel [121-123], shown in figs. 2 and 3. We emphasize however that these two classes are not in one to one correspondence with the distinction between supersymmetric and non-supersymmetric ensembles. Rather, the universality class is determined by where the relevant edge of the spectrum is soft or hard, which as explained in section 3.6.2 depends on the number of BPS states in the supersymmetric case.  \nFor WD ensembles, edge statistics fall under the universality class of Gaussian WD matrices. The eigenvalues of such matrices famously asymptote to eq. (3.37), and exhibit the soft edge in eq. (3.39) by which the spectral density goes to zero as √ x . The spectral statistics of ensembles in this universality class are governed by Hermite polynomials, which upon the largeN scaling that focuses on this soft edge give rise to Airy processes like in fig. 2 [103, 107, 124] (cf. the nearextremal black holes spectrum in eq. (2.24) from the N = 0 Schwarzian). For the standard AZ ensembles, the eigenvalues in the Gaussian case asymptote at large N to eq. (3.40), and exhibit the hard edge in eq. (3.42) by which the spectral density diverges as 1 / √ x . The spectral statistics of ensembles in this universality class are governed by Laguerre polynomials, which upon the largeN scaling that focuses on this hard edge give rise to Bessel processes like in fig. 3 [109, 124-127] (cf. the near-non-BPS black holes spectrum in eq. (2.25) from the N = 1 super-Schwarzian).  \nBy our understanding from section 3.6.2 about spectral gaps in the BPS ensembles, we learn that the Wishart generalization of the AZ ensembles is actually able to accommodate both Bessel and Airy statistics [112, 121-123]. The gapless BPS spectrum falls under Bessel universality, with a characteristic 1 / √ x hard edge just like in the standard AZ ensembles (cf. the near-BPS black hole spectrum in eq. (2.27) from the N = 2 super-Schwarzian for gapless supermultiplets). In contrast, the gapped BPS spectrum falls back under Airy universality, recovering the characteristic √ x soft edge of WD ensembles as illustrated in fig. 4 (cf. the near-BPS black hole spectrum in eqs. (2.27) and (2.29) from the N = 2 , 4 super-Schwarzians for gapped supermultiplets).', '3.6.4 Double Scaling Limits': "The largeN scalings with a finite e S 0 parameter that gave rise to the Airy and Bessel models above are somewhat trivial examples of a more sophisticated type of matrix scalings that go under the name of double scaling limits (DSL) [85, 104-106]. Such limits consist of zooming on the edge of the spectrum while simultaneously tuning the coefficients in the matrix potential to certain critical values as N becomes large. The motivation for doing so is the emergence of the structure of topological surfaces which connect matrix integrals to gravity in the continuum N →∞ limit.  \nA perturbative largeN expansion of the matrix partition function Z can be represented by double line matrix diagrams which come with factors of N χ g , where χ g = 2 -2 g is the Euler character of the corresponding graph, and g its genus. Organizing diagrams this way gives an  \nexpansion for Z itself,  \nZ = ∞ ∑ g =0 N χ g Z g , (3.47)  \nwhere Z g contains all contributions from genusg diagrams. The single largeN limits considered above thus suppress all terms but Z 0 , meaning that only planar diagrams survive. When the couplings γ in the matrix potential are allowed to vary, however, contributions from higher genus diagrams can be preserved at large N . As it turns out, the Z g functions in the genus expansion all exhibit divergences at the same critical values γ = γ c of the couplings, near which they behave as  \nZ g ≈ Z c g ( γ -γ c ) hχ g , h > 0 , (3.48)  \nfor some finite constant Z c g independent of γ . Near such critical points, the expansion in eq. (3.47) can thus be rewritten as  \nZ ≈ ∞ ∑ g =0 e χ g S 0 Z c g , e S 0 ≡ N ( γ -γ c ) h . (3.49)  \nThis expression realizes the possibility of a taking N →∞ strictly while preserving higher genus contributions by simply demanding that the new effective expansion parameter e S 0 stays finite in the limit. As defined, this clearly requires that γ → γ c as N → ∞ in the precise way specified by eq. (3.49). By coordinating this coupling limit with the appropriate edge zooming variable as N →∞ , the DSL succeeds in capturing the edge statistics exactly in e S 0 of largeN matrix models where diagrams of arbitrary topology contribute. Furthermore, there is a precise sense in which the diagrammatic contributions may be interpreted as actual topological surfaces in the DSL. Note that the number of vertices n in a given diagram appears as the power of the coupling γ , so the expectation value of vertex number at fixed g is computed by  \n⟨ n ⟩ g ≡ γ ∂ ∂γ log Z g ≈ γ hχ g γ -γ c . (3.50)  \nSince ⟨ n ⟩ g →∞ in the DSL, we see that diagrams may be understood as approaching a continuum limit of surfaces given by graphs with infinitely many vertices. In this sense, the partition function of a matrix ensemble in the DSL takes the form of a sum over surfaces of arbitrary topology of very much the same flavor as quantum gravity. In particular, the tuning of matrix couplings to different critical values in the DSL nicely corresponds to introducing different choices of matter content coupled to gravity. Of particular interest for us are the DSL matrix ensembles dual to the models of JT quantum gravity describing the physics of near-extremal black holes [19, 20, 43, 44].  \nCrucially, it can be shown that for matrix integrals in the DSL, the expansion of the spectral density to all orders in e -S 0 is fully determined by its leading from at large e S 0 , together with a discrete choice out of the ten invariant matrix ensembles. 46 This is accomplished through topological recursion relations which can be extracted from the loop equations [20, 128, 129] (see [130] for an  \noverview on topological recursion). As mentioned in section 2.4, in JT the leading spectral density at large e S 0 is computed on a trivial topology, where the dynamics reduce to a Schwarzian theory governing the fluctuations of the boundary of a hyperbolic disk. In other words, the near-extremal black hole spectral densities in eqs. (2.24) to (2.29) coming from Schwarzian theories precisely correspond to the leading spectral densities of the dual DSL matrix ensembles. More explicitly, with our choice of 'normalization' in section 3.6.1 (cf. eq. (3.30)),  \nˆ ρ N ∗ ( E ) ≡ e -S 0 ρ N ( E ) . (3.51)  \nwith the quotation marks emphasizing that actually ˆ ρ ∗ is not even normalizable, and all we are doing is stripping off the extensive factor of e S 0 from ρ N . With this convention, one easily verifies that near the lower edge of the continuous parts of their spectra, ˆ ρ N ∗ exhibits exactly the soft edge in eq. (3.39) for N = 0 and gapped N = 2 , 4, and the hard edge in eq. (3.42) for N = 1 and gapless N = 2. This is how edge universality manifests itself in the DSL, i.e., in the near-edge behavior of the DSL spectral density, which itself is a near-edge spectral density. In particular, we expect the exact DSL spectral densities to be governed by Airy statistics near soft edges, and by Bessel statistics near hard edges.", '4 Near-Extremal Saddles': 'The goal of this section is to investigate the replica trick that quenches the thermal entropy in a way that elucidates its workings in quantum gravity. To do so we pursue the calculation in the nearextremal regime employing only those tools from random matrix theory with a clear gravitational dual description. Our results show that a treatment strictly to leading order suffices, the crucial realization being that a new saddle point arises at large β ≳ O ( e S 0 ) and becomes dominant. For the replica trick in eq. (1.12), the behavior of this saddle point is shown to be determined by spectral edge statistics, implying that to leading order the quenched entropy behaves universally and is completely independent of non-perturbative details. On the gravity side, this saddle point is seen to be realized by a particular type of branes associated to energy eigenvalues. Because these branes are needed at leading order, our findings imply that their inclusion on the gravity side is non-optional for a genuine duality between JT gravity and matrix integrals. 47  \nTo implement eq. (1.12) we need to study the moments of Z ( β ) in a matrix ensemble,  \n⟨ Z ( β ) m ⟩ = ˆ R N d Λ p (Λ) ( Tr e -β Λ ) m , (4.1)  \nwhere p is given in eq. (3.15) and Λ denotes the diagonal matrix of N eigenvalues of the random Hamiltonian. As in eq. (3.14), it is again convenient to express this as an action integral,  \n⟨ Z ( β ) m ⟩ = 1 Z ˆ R N d Λ e -NI m (Λ) , I m (Λ) ≡ I (Λ) + m N log Tr e -β Λ , (4.2)  \nwith the matrix action I given by eq. (3.16). This integral admits a saddle-point approximation at large N . One could follow the same strategy as in section 3.5 to study this integral in the  \ncontinuum limit to leading order. If done na¨ıvely, however, the resulting functional integral and variational problem would in fact be ill-posed for any m> 0. This is because the Z ( β ) continuum one would get by replacing a discrete measure ˆ ρ Λ by a continuous one ˆ ρ (cf. eqs. (3.22) and (3.23)),  \nˆ R dx ˆ ρ Λ ( x ) e -βx -→ ˆ R dx ˆ ρ ( x ) e -βx , (4.3)  \nwould generally cause the action to be unbounded from below. This fundamental issue would be caused by shallow tails of the continuum measure ˆ ρ extending to arbitrarily negative values. 48 But such tails would clearly not make sense from the discretum: there must be a lowest eigenvalue which, in a saddle-point approximation, would lie at some finite value. The effect of Z ( β ) in eq. (4.2) thus has to be better understood before going to the largeN continuum.  \nValuable intuition can be gained from a finiteN evaluation of the object whose largeN limit gives the saddle point we are after. Namely, thinking of eq. (4.2) altogether as just another matrix integral, we may define its joint eigenvalue PDF by (cf. eqs. (3.14) and (3.15))  \np m β (Λ) = e -NI m (Λ) Z ⟨ Z ( β ) m ⟩ , (4.4)  \nand correspondingly by eqs. (2.5) and (2.7) obtain an average spectral density via  \n〈 ρ m β ( λ ) 〉 ≡ N ˆ R N -1 dλ 2 · · · dλ N p m β ( λ, λ 2 , . . . , λ N ) . (4.5)  \nThis is precisely the object we need to understand in the largeN limit. In particular, by eq. (3.30),  \nˆ ρ m β ∗ ( x ) ≡ lim N →∞ 〈 ρ m β ( x ) 〉 N , (4.6)  \nis the equilibrium measure that a saddle-point evaluation of eq. (4.2) must give if the continuum limit is taken correctly. As an example, we integrate eq. (4.5) exactly for the GUE at finite N and show the result of varying m and β in fig. 5. The observed behavior of 〈 ρ m β ( x ) 〉 suggests that in a largeN limit where β ∼ O ( N ) the lowest eigenvalue naturally behaves in a discrete manner, explaining why eq. (4.2) does not admit a na¨ıve continuum limit. With the intuition that fig. 5 provides, we now derive a general expression for eq. (4.2) that actually admits such a limit.', '4.1 Low Temperatures': "The form of log Z ( β ) at small and large β can be read off from eqs. (2.11) and (2.12), respectively. At small β , this goes as log Z ( β ) ∼ log N , whereas at large β it behaves as log Z ( β ) ∼ -β minΛ.  \n̂ I m [ˆ ρ ] ! = 1 2 ˆ R dx ˆ ρ ( x ) x 2 -ˆ R dxdy ˆ ρ ( x )ˆ ρ ( y ) log | x -y | -m N 2 log   N ˆ R dx ˆ ρ ( x ) e -βx   .  \nτ  \nConsider letting ˆ ρ ( x ) decay like e -1 τ ( | x | a ) as x →-∞ for 1 < τ < 2. The lower bound guarantees that the integral in the β -dependent term converges for any β , while the upper bound implies that this term will grow negative faster in a than the potential term grows positive. Then one finds that ̂ I m behaves as -τ -1 τ m ( aβ ) τ τ -1 at large a , which implies the action is unbounded from below since ̂ I m can be made arbitrarily negative by taking a →∞ .  \nFigure 5 . Illustration of ⟨ ρ m β ( x ) ⟩ by exact evaluation of the integral in eq. (4.5) for the GUE at finite N . The red curve corresponds to N = 5 and is plotted for the values of mβ shown; here m = 2, but the relevant qualitative features are seen to only depend on the ratio mβ/N . The black curve corresponds to N = 4 with m = 0 fixed, and serves as a reference for the behavior of the red curve. As mβ/N increases, the peak associated to the lowest eigenvalue of the N = 5 curve separates from the other N -1 eigenvalues, which converge down to the N = 4 curve. In other words, only one eigenvalue depends significantly on β , while the others are mostly insensitive to the Z ( β ) insertion. The single dynamical eigenvalue can be analytically shown to asymptote with mβ/N to a Gaussian of variance 1 /N around -mβ/N . Importantly, this peak shrinks as N →∞ , giving no overlap with the bulk of eigenvalues at large N for any constant mβ/N > 0.  \n<!-- image -->  \nHence at large N , if β is kept fixed, the contribution from the β -dependent term in eq. (4.2) gets suppressed by an extra factor of N relative to I (Λ). In a largeN treatment of this expectation value integral, the saddle points would thus be solely determined by the original matrix integral and completely unaffected by the insertion of Z ( β ) m .  \nThe more interesting regime arises at large β , when log Z ( β ) ∼ -β minΛ. In this case, by making β ∼ O ( N ), the two action terms in eq. (4.2) become of the same order in the largeN limit. This must lead to a saddle point for the expectation value integral strictly different from that of the original matrix integral due to the effect of the Z ( β ) m insertion. The new saddle point clearly appears for any m> 0 and gives a distinct role to the lowest eigenvalue min Λ in the integral, as observed also in fig. 5. In particular, we expect this lowest eigenvalue to be particularly sensitive to β , and thus capture the dependence of ⟨ Z ( β ) m ⟩ on it in the largeβ regime we are interested in. Since β ∼ O ( N ) translates into β ∼ O ( e S 0 ) in the DSL, note that this is precisely the nearextremal regime where the gravitational entropy generically begins to turn negative, as discussed in section 1.4. In other words, the appearance of this new saddle point is perfectly consistent with the expectation that annealed and quenched entropies must depart from each other near extremality.  \nTo track the behavior of the lowest eigenvalue, consider splitting Λ into 49  \n̸  \nˆ λ ≡ minΛ , ̂ Λ ≡ diag { λ k = ˆ λ } N k =1 . (4.7)  \nSeparating eigenvalues this way, the canonical partition function may equivalently be written  \nlog Z ( β ) = -β ˆ λ +log   1 + ∑ λ ∈ ̂ Λ e -β ( λ -ˆ λ )   , (4.8)  \nan expression which is natural for large β but holds for any β . Any global eigenvalue degeneracies coming from the symmetries of the ensemble can be factored out from this treatment. As for dynamical degeneracies, recall that the Vandermonde determinant exerts a repulsion across all eigenvalues which forbids them. In particular, this means that λ > ˆ λ can be assumed strictly for all λ ∈ ̂ Λ. In general though, every eigenvalue λ ∈ ̂ Λ can get arbitrarily close to ˆ λ , so it is not 'a priori clear whether the sum in eq. (4.8) can be assumed to be exponentially suppressed in β . The key observation is that the typical separation between eigenvalues that the matrix integral induces is no smaller than O ( N -1 ), and often larger than that near spectral edges (cf. section 3.6.1). Since we will be interested in scaling β linearly in N , this statistical separation between eigenvalues will always be counteracted by β in the exponent. Hence we expect the logarithmic term in eq. (4.8) to statistically be at most O (log N ), and thus subleading relative to the β ∼ O ( N ) term. As a result, the largeβ approximation  \nlog Z ( β ) ≈ -β ˆ λ, (4.9)  \nholds generically across the ensemble. In fact, notice that this -β ˆ λ term will always be favoring strictly smaller values for ˆ λ in the matrix integral than for any λ ∈ ̂ Λ. Therefore we actually expect the typical difference λ -ˆ λ ≳ O ( N -1 ), thereby enhancing the approximation in eq. (4.9) to be accurate up to at most O (1) corrections, as our results turn out to confirm (cf. fig. 5). The strategy in what follows is to use the largeβ behavior of Z ( β ) in eq. (4.9) to break down ⟨ Z ( β ) m ⟩ into matrix integrals which are computable using a saddle-point approximation.", '4.2 Eigenvalue Splitting': 'Given the special status the lowest eigenvalues acquires, it will be useful to identify it inside the matrix integral. As expressed throughout, spectral integrals run over all of R for every eigenvalue in Λ, so any one of them can be the smallest in some integration domain. However, all integrals are fully symmetric under permutations of the eigenvalues. This allows one to reorder them so as to always keep track of which eigenvalue is the lowest. Using the notation in eq. (4.7), one easily arrives at the following identity  \nˆ R N d Λ( · ) = N ˆ R N -1 d ̂ Λ min ̂ Λ ˆ -∞ d ˆ λ ( · ) . (4.10)  \nThis is useful because under this reorganization of the eigenvalue integral we can identify ˆ λ as a specific integration variable corresponding to min Λ, and thus study what happens to the lowest  \nFigure 6 . Illustration of F m β (red) against ⟨ ρ m β ⟩ (black) by direct integration of eqs. (4.5) and (4.13), respectively, for the GUE at N = 3 (cf. fig. 5). Here β is being varied at fixed m = 6, but the qualitative features are seen to only be sensitive to mβ/N . As desired, the eigenvalue splitting from eq. (4.11) allows for F m β to isolate the PDF of just the statistically lowest eigenvalue within the total spectral density ⟨ ρ m β ⟩ . At large β , F m β ( ˆ λ ) tends to simply F ( ˆ λ ) e -mβ ˆ λ , while the other eigenvalues converge to the N = 2 spectrum.  \n<!-- image -->  \neigenvalue inside the full matrix integral. The representation in eq. (4.10) can be read as upper bounding the ˆ λ eigenvalue by all other ̂ Λ eigenvalues. However, one could also think of all eigenvalues ̂ Λ as being upper bounded by ˆ λ , which gives the alternative representation  \nˆ R N d Λ( · ) = N ˆ R d ˆ λ ˆ R N -1 ≥ ˆ λ d ̂ Λ( · ) . (4.11)  \nBoth eqs. (4.10) and (4.11) become useful and interchangeable at the saddle-point level if one trades between an upper bound constraint for ˆ λ and a lower bound constraint for ̂ Λ. In particular, eq. (4.11) allows one to define an induced PDF on the lowest eigenvalue. Namely, applying eq. (4.11) to the joint PDF of eigenvalues p (Λ) = p ( ˆ λ, ̂ Λ) gives the PDF for the lowest eigenvalue, 50  \nF ( ˆ λ ) ≡ N ˆ R N -1 ≥ ˆ λ d ̂ Λ p ( ˆ λ, ̂ Λ) . (4.12)  \nApplying this splitting to eq. (4.1) gives the PDF induced on ˆ λ under the insertion of Z ( β ) m ,  \nF m β ( ˆ λ ) ≡ N ⟨ Z ( β ) m ⟩ ˆ R N -1 ≥ ˆ λ d ̂ Λ p ( ˆ λ, ̂ Λ) ( e -β ˆ λ +Tr e -β ̂ Λ ) m . (4.13)  \nAn illustration of this PDF for finite N is shown in fig. 6. In consistency with eq. (4.9), we observe that F ( ˆ λ ) e -mβ ˆ λ approximates F m β ( ˆ λ ) very well as β → ∞ . However, making such a largeβ approximation from the onset would be too drastic and not very illuminating. Instead, our strategy will be to take a continuum limit, use saddle-point methods to simplify the integral, and only make large N and β at the controlled level of the equations of motion. 51  \nWith this strategy in mind, consider breaking up the integral in eq. (4.2) in terms of the decomposition of Λ into ˆ λ and ̂ Λ. Firstly, writing out Z ( β ), the total action is  \nI m (Λ) = I (Λ) -m N log ( e -β ˆ λ +Tr e -β ̂ Λ ) . (4.14)  \nThe original action I splits into  \nI (Λ) = I ( ̂ Λ) + V ν ( ˆ λ ) -β N log | det( ̂ Λ -ˆ λ ) | , (4.15)  \nwhere the interaction is the Vandermonde determinant. Recalling eq. (3.19), we recognize the ˆ λ -dependent terms in eq. (4.15) as precisely giving the effective potential for the ˆ λ eigenvalue,  \n̂ V ( ˆ λ ) = V ν ( ˆ λ ) -β N log | det( ̂ Λ -ˆ λ ) | . (4.16)  \nPutting this altogether in eq. (4.2) and using the integral representation from eq. (4.11) gives  \n⟨ Z ( β ) m ⟩ = N Z ˆ R d ˆ λe -NV ν ( ˆ λ ) ˆ R N -1 ≥ ˆ λ d ̂ Λ e -NI ( ̂ Λ) | det( ̂ Λ -ˆ λ ) | β ( e -β ˆ λ +Tr e -β ̂ Λ ) m . (4.17)  \nThis matrix integral is now in a form that is amenable to a continuum treatment. Crucially, in contrast with eq. (4.2), the distinct behavior of the ˆ λ eigenvalue has now been isolated such that a continuum limit can be taken for the ̂ Λ integral alone (cf. footnote 48). The reader may recognize the determinant in eq. (4.17) as the first hint that branes will be important in the gravitational avatar of this calculation.', '4.3 Single-Eigenvalue Instanton': "As discussed so far, inserting Z ( β ) m inside the matrix integral and making β ∼ O ( N ) causes the lowest eigenvalue ˆ λ to behave differently from the rest. This single eigenvalue behaves in such a way that it dynamically separates from the others at large β , a phenomenon which we already illustrated for the GUE in fig. 5. In more general matrix integrals this may be recognized as the characteristic behavior of a single-eigenvalue instanton (see e.g. [71]). The purpose of this section is to work out the largeN instanton that characterizes ⟨ Z ( β ) m ⟩ at large β . 52  \nTo tackle the integral in eq. (4.17) by saddle-point approximation, it becomes relevant to take a continuum largeN limit first. Although at large N the distinction between N and N -1 becomes immaterial to leading order, it will be useful to distinguish such terms in order to more reliably connect results back to finite N later on. Following section 3.5, in the continuum the integral over N -1 eigenvalues ̂ Λ bounded from below by ˆ λ becomes a functional integral over spectral densities ( N -1) ' D ˆ σ where ˆ σ is unit-normalized and restricted to have supp ˆ σ ⊆ [ ˆ λ, ∞ ). Correspondingly, the action I ( ̂ Λ) becomes a functional ( N -1) ̂ I [ˆ σ ] where the summations over N -1 eigenvalues are replaced by integrals ( N -1) ' R dx ˆ σ ( x ). The upshot is 53  \n̂ I [ˆ σ ] = ˆ R dx ˆ σ ( x ) V ν ( x ) -N -1 N β 2 ˆ R dxdy ˆ σ ( x )ˆ σ ( y ) log | x -y | , (4.18)  \nwhere we could have equivalently written the integrals over just [ ˆ λ, ∞ ) given the constrained support of ˆ σ . By eq. (3.28), functionally varying ̂ I in ˆ σ gives the effective potential  \n̂ V [ˆ σ ; x ] = V ν ( x ) -N -1 N β ˆ R dy ˆ σ ( y ) log | x -y | . (4.19)  \nAs for the determinant in eq. (4.17), using the identity log det M = Trlog M , in the continuum  \n| det( ̂ Λ -ˆ λ ) | β -→ Ψ[ˆ σ ; ˆ λ ] ≡ exp   ( N -1) β ˆ R dx ˆ σ ( x ) log | x -ˆ λ |   . (4.20)  \nwhich can also be written in terms of the effective potential in eq. (4.19) as  \nΨ[ˆ σ ; x ] e -NV ν ( x ) = e -N ̂ V [ˆ σ ; x ] . (4.21)  \nThe choice of notation here is suggestive of the brane operators of [19], which correspond to the full left-hand side of eq. (4.21). Finally, in the continuum Z ( β ) as written in eq. (4.17) becomes  \nZ ( β ) -→ Z [ˆ σ ; ˆ λ ] ≡ e -β ˆ λ +( N -1) ˆ R dx ˆ σ ( x ) e -βx , (4.22)  \nwhere the β -dependence is left implicit. This continuum limit of Z should be contrasted with the na¨ıve one in eq. (4.3). The complete continuum form of eq. (4.17) thus reads  \n⟨ Z ( β ) m ⟩ = N ( N -1) Z ˆ R d ˆ λe -NV ν ( ˆ λ ) ˆ D ˆ σ e -N ( N -1) ̂ I [ˆ σ ] Ψ[ˆ σ ; ˆ λ ] Z [ˆ σ ; ˆ λ ] m . (4.23)  \nIn the largeN limit, note how the action term above scales as N 2 , whereas the Ψ insertion only goes as N , as can be seen from eq. (4.20). This means that, to leading order, the saddle point ˆ σ = ˆ σ ∗ that dominates the functional integral in eq. (4.23) is insensitive to Ψ. Intuitively, the determinant insertion of a single eigenvalue that Ψ encodes has a subleading effect on the dynamics of the O ( N ) eigenvalues that participate in ̂ I ; the effect of the latter on the former will however not be negligible. Additionally, from the discussion below eq. (4.8) we know that Z is dominated at large β by the lowest eigenvalue ˆ λ alone with ˆ σ only contributing exponentially suppressed corrections, so ˆ σ ∗ will also be insensitive to β at leading order. Altogether, the extremization problem that defines ˆ σ ∗ does not depend on ˆ λ at all to leading order at large N and β , in perfect consistency with the examples in figs. 5 and 6.  \nExtremizing eq. (4.23) thus amounts to applying the equilibrium condition, eq. (3.29), to the effective potential in eq. (4.19). As a result, the saddle point ˆ σ ∗ corresponds to the unit-normalized measure that solves the singular integral equation  \nR dy ˆ σ ∗ ( y ) x -y = N N -1 V ' ν ( x ) β , ∀ x ∈ supp ˆ σ ∗ . (4.24)  \nJust like eq. (3.31), this equation can be solved exactly for ˆ σ ∗ , with the general solution taking a form analogous to eq. (3.32). When approximating the ˆ σ integral in eq. (4.23) by its saddle point ˆ σ ∗ , consistency with the unit-normalization of the expectation value at m = 0 requires to also evaluate Z on the same saddle. The saddle-point approximation Z ∗ for Z reads  \nZ ∗ ≡ N ( N -1) e -N ( N -1) ̂ I [ˆ σ ∗ ] ̂ Z ∗ , ̂ Z ∗ ≡ ˆ λ 0 ˆ -∞ d ˆ λ Ψ[ˆ σ ∗ ; ˆ λ ] e -NV ν ( ˆ λ ) = ˆ λ 0 ˆ -∞ d ˆ λe -N ̂ V ∗ ( ˆ λ ) , (4.25)  \nFigure 7 . Illustration of F m β ∗ (red, solid) given by eq. (4.28) evaluated on the continuous ˆ σ ∗ saddle that solves eq. (4.24) for the GUE at N = 5 (dashed). Together, F m β ∗ ( x ) + ( N -1)ˆ σ ∗ ( x ) is an approximation to ⟨ ρ m β ( x ) ⟩ in which the lowest eigenvalue is only sensitive to the saddle form of the higher eigenvalues. The black curve corresponds to the exact result for ⟨ ρ m β ⟩ using eq. (4.13) (cf. fig. 5), and which also captures F m β ( x ) from eq. (4.13) exactly (cf. fig. 6). Here m = 1, but once again, qualitative features mostly depend on mβ/N . Already at such small N = 5, we see that F m β ∗ rapidly converges to F m β at large β , even though ˆ σ ∗ is a very coarse approximation to the N = 4 spectral density that the higher eigenvalues asymptote to.  \n<!-- image -->  \nwhere we have used eq. (4.21), the notation V ∗ refers to the leading effective potential with ˆ σ = ˆ σ ∗ (cf. eq. (3.33)), and we have defined the lower edge of the leading spectrum  \nˆ λ 0 ≡ inf supp ˆ σ ∗ . (4.26)  \nThe new constant ̂ Z ∗ is convenient because after evaluating the saddle-point approximation for eq. (4.23), all other constants cancel out and one is left with  \n⟨ Z ( β ) m ⟩ ∗ = 1 ̂ Z ∗ ˆ λ 0 ˆ -∞ d ˆ λe -N ̂ V ∗ ( ˆ λ ) Z [ˆ σ ∗ ; ˆ λ ] m . (4.27)  \nNote that ⟨ Z ( β ) m ⟩ ∗ is indeed consistently normalized to unity for m = 0. The integrand above provides a largeN approximation to the lowest eigenvalue PDF we encountered in eq. (4.13),  \nF m β ∗ ( ˆ λ ) = 1 ̂ Z ∗ ⟨ Z ( β ) m ⟩ ∗ e -N ̂ V ∗ ( ˆ λ ) Z [ˆ σ ∗ ; ˆ λ ] m , ˆ λ ≤ ˆ λ 0 . (4.28)  \nThis PDF is illustrated for the GUE in fig. 7. Physically, eq. (4.28) captures the statistics of the ground-state energy ˆ λ in the presence of a classical thermal gas of O ( N ) higher energies that localize in the ˆ σ ∗ saddle. Importantly, ˆ λ must by construction lie below the lower edge of ˆ σ ∗ , sometimes referred to as the classically forbidden region. However, we emphasize that there is nothing non-classical or forbidden about this instanton: it is a leading saddle (i.e., classical), and its location is perfectly allowed due to the largeβ back-reaction that Z ( β ) sources.  \nIf the result for the integrand in eq. (4.27) is sufficiently simple, one could now attempt to directly evaluate the integral over ˆ λ to obtain a final answer for ⟨ Z ( β ) m ⟩ ∗ . However, since eq. (4.27) itself already is the result of a largeN saddle-point approximation, it is in fact more natural to also approximate the ˆ λ integral by its saddle-point value. Indeed, this is a very reasonable approximation at large N , since the associated PDF for the lowest eigenvalue in eq. (4.28) clearly has a characteristic width which shrinks as 1 /N at large N (cf. fig. 5). By reducing this PDF to the dynamics of where it peaks, we will be treating the lowest eigenvalue as a point particle and understanding its motion as a function of β . Recalling that we are in a regime in which β and N  \nare simultaneously large, the extremization problem we are after will be sensitive to the presence of the Z factor. Fortunately, from eq. (4.8) we see that Z at large β will only depend on ˆ λ up to exponentially suppressed corrections from higher eigenvalues. Hence using the approximation in eq. (4.9) and extremizing, the saddle point ˆ λ = ˆ λ ∗ we are after is the solution to 54  \nd d ˆ λ ∣ ∣ ∣ ∣ ˆ λ = ˆ λ ∗ e -N ̂ V ∗ ( ˆ λ ) -mβ ˆ λ = 0 , ˆ λ ∗ ≤ ˆ λ 0 , (4.29)  \nThis reduces to a simple equation for the location of ˆ λ ∗ in terms of the effective potential it feels,  \n̂ V ' ∗ ( ˆ λ ∗ ) = -mβ N , ˆ λ ∗ ≤ ˆ λ 0 , (4.30)  \nwhere recall that, on the support of ˆ σ ∗ , the equilibrium condition from eq. (3.29) says  \n̂ V ' ∗ ( ˆ λ ) = 0 , ∀ ˆ λ ∈ supp ˆ σ ∗ . (4.31)  \nIntuitively, this means that the lowest eigenvalue feels a force proportional to mβ/N pushing it to smaller values away from the rest, which stay at equilibrium on supp ˆ σ ∗ feeling no force. Written out, eq. (4.30) gives the equation of motion for the instanton that we were after:  \nV ' ν ( ˆ λ ∗ ) + N -1 N β ˆ R dx ˆ σ ∗ ( x ) x -ˆ λ ∗ = -mβ N , ˆ λ ∗ ≤ ˆ λ 0 . (4.32)  \nIt is not necessary to write the integral as a principal value now because ˆ λ ∗ is required to lie off the support of ˆ σ ∗ anyway. This general analysis explains the behavior exemplified in figs. 5 to 7, and realizes eq. (4.6) as  \nˆ ρ m β ∗ ( x ) = δ ( x -ˆ λ ∗ ) + ( N -1) ˆ σ ∗ ( x ) . (4.33)  \nBut recall that off supp ˆ σ ∗ , the left-hand side of eq. (4.30) is precisely what defines the spectral curve in eq. (3.34). Hence in terms of y , eq. (4.30) reduces to simply 55  \ny ( ˆ λ ∗ ) = mβ β N , ˆ λ ∗ ≤ ˆ λ 0 . (4.34)  \nThis final expression for the location of ˆ λ ∗ is particularly powerful for several reasons. We were after a saddle-point evaluation of eq. (4.23), which required us to find saddle points ˆ λ ∗ and ˆ σ ∗ for each integral. The equation of motion we found for ˆ σ ∗ is eq. (4.24), which at large N describes nothing but the usual leading spectral density of the ensemble and is explicitly solved by eq. (3.32). As for ˆ λ ∗ , what we found is the algebraic equation in eq. (4.34), where y follows directly from analytically continuing ˆ σ ∗ as in eq. (3.35). All in all, this means that our single-eigenvalue instanton for ⟨ Z ( β ) m ⟩ is fully characterized just from knowledge of the leading spectral density of the ensemble.  \nVery importantly, our final results do not depend on objects like the matrix action and potential, which are not defined in scaling limits. Instead, all we need is access to the leading spectral  \ndensity of the matrix ensemble, which is well-defined under matrix scaling. This obviously includes matrix ensembles in the DSL, as relevant for gravity. For these, not only are the requisite ingredients well-defined on the matrix side, but also accessible from the gravity side. Indeed, as we elaborate on in section 4.6, all one needs is to account for a particular type of brane effects that the theory requires.", '4.4 Instanton Action': 'We would like to finally obtain an expression for eq. (4.2), the matrix integral that computes moments of Z ( β ). At the saddle-point level, this amounts to evaluating the action of the singleeigenvalue instanton found in section 4.3. This instanton corresponds to the full saddle-point approximation of both ˆ σ and ˆ λ applied to eq. (4.17). Recalling the saddle-point result for ˆ σ ∗ from eq. (4.27) and further extremizing in ˆ λ gives  \n⟨ Z ( β ) m ⟩ ∗ = e -N ̂ V ∗ ( ˆ λ ∗ ) e -N ̂ V ∗ ( ˆ λ 0 ) Z [ˆ σ ∗ ; ˆ λ ∗ ] m , (4.35)  \nwhere in evaluating ̂ Z ∗ we have to approximate eq. (4.25) by a saddle point ˆ λ 0 for ˆ λ corresponding to m = 0. This is the solution ˆ λ ∗ = ˆ λ 0 to eq. (4.34) for m = 0, which is given by eq. (4.26) as can be easily seen by comparing eqs. (4.30) and (4.31). For m = 0 we thus obtain unity, in consistency with the proper normalization of the saddle-point approximation. Using now eq. (4.21), we reduce eq. (4.35) to the final expression  \n⟨ Z ( β ) m ⟩ ∗ = e -NI ∗ ( m ) , (4.36)  \nwhere we have defined the instanton action,  \nI ∗ ( m ) ≡ ̂ V ∗ ( ˆ λ ∗ ) -̂ V ∗ ( ˆ λ 0 ) -m N log Z [ˆ σ ∗ ; ˆ λ ∗ ] . (4.37)  \nHere, recall that ˆ σ ∗ and ˆ λ ∗ are the solutions to eqs. (4.24) and (4.30), respectively, and that ˆ λ ∗ implicitly depends on m , with ˆ λ ∗ = ˆ λ 0 for m = 0. This is our general result for the saddle-point approximation to ⟨ Z ( β ) m ⟩ at large β and N for any value of m . Written out in full glory,  \n⟨ Z ( β ) m ⟩ ∗ = e -N ( ̂ V ∗ ( ˆ λ ∗ ) -̂ V ∗ ( ˆ λ 0 ))   e -β ˆ λ ∗ +( N -1) ˆ R dx ˆ σ ∗ ( x ) e -βx   m . (4.38)  \nBecause the instanton eigenvalue ˆ λ ∗ depends non-trivially on m , this result clearly does not factorize for m ≥ 1 moments. This is consistent with the expectation that factorization in m should only occur at small β . Although we are primarily interested in large β , note that our results also capture this asymptotic factorization in the smallβ regime. Namely, at small β the insertion of Z ( β ) does not modify the integral saddle, and indeed we see from eqs. (4.30) and (4.31) that to leading order the solution for ˆ λ ∗ recedes to the spectral edge at ˆ λ = ˆ λ 0 . Hence eq. (4.38) reduces to  \n⟨ Z ( β ) m ⟩ ∗ ≈   e -β ˆ λ 0 +( N -1) ˆ R dx ˆ σ ∗ ( x ) e -βx   m ≈   N ˆ R dx ˆ σ ∗ ( x ) e -βx   m , β ≪ N (4.39)  \nwhere in the last approximation we used large N . Not only do we see factorization at small β , but also consistency with the finiteβ result that one can easily compute for m = 1. This is the simple case in which the matrix integral only depends on the average spectral density (cf. eq. (2.8)). A saddle-point approximation for this case would just give  \n⟨ Z ( β ) ⟩ ∗ = N ˆ R dx ˆ ρ ∗ ( x ) e -βx , (4.40)  \nwith ˆ ρ ∗ the solution to eq. (3.31). Clearly eq. (4.39) agrees with eq. (4.40) for m = 1, up to the separate treatment of the lowest eigenvalue in eq. (4.39), or in the last approximation up to the mismatch between N and N -1 in the equations for ˆ ρ ∗ and ˆ σ ∗ , negligible at large N . Namely, note that eq. (4.32) to leading order at large N with β kept fixed reduces to eq. (3.31). If β is large, however, the spectral density term in eq. (4.38) in fact becomes small against the eigenvalue instanton contributions. In this regime ˆ λ ∗ also does not necessarily lie close to ˆ λ 0 , so the largeβ behavior of eq. (4.38) turns out to actually be dominated by  \n⟨ Z ( β ) m ⟩ ∗ ≈ e -mβ ˆ λ ∗ -N ( ̂ V ∗ ( ˆ λ ∗ ) -̂ V ∗ ( ˆ λ 0 )) . (4.41)  \nClearly, this behavior would be completely missed were the eigenvalue instanton to be ignored. In particular, this has implications for any moment m , and thus not only for the quenched logarithm, but also for the annealed one, as we discuss in section 4.5.  \nWe emphasize that the largeN result for ⟨ Z ( β ) m ⟩ ∗ obtained here applies in full generality to any matrix integral with an arbitrary matrix potential and eigenvalue measure of the general form posed in eq. (3.10). Furthermore, notice that eq. (4.37) only depends on quantities like the effective potential, which are well-defined in the DSL ensembles relevant for gravity.  \nAs explained in section 3.6.4, for DSL integrals the typical starting point is the leading spectral density ˆ ρ ∗ itself. We already pointed out at the end of section 4.3 that ˆ ρ ∗ is actually all one needs to define the equation of motion for ˆ λ ∗ , given that ˆ ρ ∗ uniquely determines the spectral curve. Furthermore, for the same reason, we see that ˆ ρ ∗ itself also defines the effective potential through eq. (3.34) up to an integration constant. This constant is irrelevant, as it is cancelled out by the difference of effective potentials in eq. (4.37). It thus follows that everything that is needed to define and evaluate the instanton action is ˆ ρ ∗ , which can be taken for granted in DSL ensembles.  \nAlthough ˆ ρ ∗ in the DSL is generally non-normalizable due to the scaling, the correct identification with ˆ σ ∗ is unambiguous. Namely, the relevant DSL spectral density is to be obtained by scaling a normalized spectral density, thus giving a ˆ ρ ∗ that is non-extensive in the DSL parameter e S 0 (cf. the left-hand side of eq. (3.51)). As a result, eq. (4.38) in the DSL becomes simply  \n⟨ Z ( β ) m ⟩ ∗ = e -e S 0 ̂ V ∗ ( ˆ λ ∗ )   e -β ˆ λ ∗ + e S 0 ˆ R dx ˆ ρ ∗ ( x ) e -βx   m , (4.42)  \nwhere we have used the convention that ̂ V ∗ (0) = 0 at the leading spectral edge.', '4.5 Quenched Universality': 'The crucial question that remains is whether our results for ⟨ Z ( β ) m ⟩ ∗ succeed in implementing the replica trick in eq. (1.12), which takes the subtle m → 0 limit. In other words, we still have to  \nunderstand whether working with just a saddle-point approximation for simultaneously large N and β suffices to restore the non-negativity of the quenched entropy, eq. (1.9), against the arbitrary negativity of the annealed one, eq. (1.8).  \nAll of our results follow from a largeN treatment while accounting for a β ∼ O ( N ) scaling, but no assumption was ever made about m . In this paper m started life as a non-negative integer determining the moment m of Z ( β ) under consideration. But in fact, notice that beginning already in eq. (4.1), our equations and limits all make sense for any real m ≥ 0. Our analytic continuation of m from integers to reals is thus unambiguously unique and well-defined from the onset. 56 A controlled way of understanding what happens to our results as m → 0 is to study the equations that govern our saddle points in this limit (cf. the approach of [40] in semiclassical gravity).  \nThe pertinent equations of motion are eqs. (4.24) and (4.34), where clearly only the latter depends on m . Actually, both m and β appear just in eq. (4.34), and only through the effective parameter mβ . It would thus seem that the largeβ regime we are interested in is in direct conflict with the m → 0 replica limit (cf. the discussion at the end of section 2.2). In fact, every largeβ approximation we made was perfectly consistent with m being arbitrarily small. In addition, the replica trick crucially takes m → 0 at fixed β , so there is no ambiguity in the order of limits. Namely, m becomes arbitrarily small in the replica trick, so the effective parameter mβ must be treated as small, even at large β . 57 This observation has far-reaching consequences which we elaborate on in section 4.6.  \nAs m → 0, we expect the saddle ˆ λ ∗ to return to the lower edge of the spectrum. Indeed, the solution to eq. (4.34) in this limit is ˆ λ ∗ = ˆ λ 0 as can be easily seen by comparing eqs. (4.30) and (4.31). By continuity of the spectral curve away from supp ˆ σ ∗ , for small m ≳ 0 we thus expect ˆ λ ∗ to be in the neighborhood of ˆ λ 0 . The universality discussion in section 3.6.3 implies that the near-edge region of the spectral curve of any matrix integral is solely determined by the ensemble measure, and independent of any choice of matrix potential. This already leads us to expect that the behavior of our largeβ saddle in the replica trick will be universal, and thus yield a universal result for the quenched logarithm. Importantly, this in particular means that the kind of nonperturbative details in matrix integrals that gravity does not have access to will be completely irrelevant in our leading-order results.  \nApplying the replica trick in eq. (1.12) to eq. (4.36),  \n⟨ log Z ( β ) ⟩ ∗ = lim m → 0 d dm ⟨ Z ( β ) m ⟩ ∗ = -N lim m → 0 dI ∗ ( m ) dm , (4.43)  \nwhere we have used that I ∗ (0) = 0. The evaluation of the total derivative in m is particularly simple if one recalls that I ∗ is precisely a stationary point of the matrix action in ˆ λ . Hence any implicit dependence on m coming from the stationary value ˆ λ = ˆ λ ∗ will give zero, implying that  \ndI ∗ ( m ) dm = ∂I ∗ ( m ) ∂m . (4.44)  \nSince the only explicit dependence of I ∗ ( m ) on m is through the prefactor of the logarithm, the result for the quenched logarithm in eq. (4.43) is simply  \n⟨ log Z ( β ) ⟩ ∗ = lim m → 0 log Z [ˆ σ ∗ ; ˆ λ ∗ ] = log Z [ˆ σ ∗ ; ˆ λ 0 ] , (4.45)  \nExplicitly, using eq. (4.22) our final expression is  \n⟨ log Z ( β ) ⟩ ∗ = log   e -β ˆ λ 0 + N ˆ R dx ˆ σ ∗ ( x ) e -βx   , (4.46)  \nwhere recall ˆ λ 0 = inf supp ˆ σ ∗ and ˆ σ ∗ is the leading spectral density of the matrix ensemble. If normalizable, ˆ σ ∗ is unit-normalized, but if not (as in the DSL), ˆ σ ∗ is a spectral density with no scaling in the scale parameter e S 0 (cf. the discussion around eq. (4.42)). In extreme regimes of scaling β against N , this quenched result behaves as 58  \n⟨ log Z ( β ) ⟩ ∗ ≈            log N +log ˆ R dx ˆ σ ∗ ( x ) e -βx , β ≪ N, -β ˆ λ 0 + N ˆ R dx ˆ σ ∗ ( x ) e -β ( x -ˆ λ 0 ) , β ≫ N. (4.47)  \nNote that the last integral is a Laplace transform for a canonically shifted spectral density, and thus tends to zero in the largeβ limit by the results of section 2.3. These quenched result should be contrasted with the annealed one, which corresponds to taking the logarithm of the m = 1 moment. This case corresponds to just eq. (4.38) with m = 1. At small or finite β , we already argued around eq. (4.39) that this approximately reduces to  \nlog ⟨ Z ( β ) ⟩ ∗ ≈ log N +log ˆ R dx ˆ σ ∗ ( x ) e -βx , β ≪ N, (4.48)  \nthereby matching the smallβ behavior of the quenched logarithm in eq. (4.47). If β is large, however, this term becomes small relative to the eigenvalue instanton, leading to  \nlog ⟨ Z ( β ) ⟩ ∗ ≈ -β ˆ λ ∗ -N ( ̂ V ∗ ( ˆ λ ∗ ) -̂ V ∗ ( ˆ λ 0 )) , β ≫ N, (4.49)  \nwhich is drastically different from the quenched result in eq. (4.47), as will more explicitly be seen in concrete examples below. This largeβ behavior of ⟨ Z ( β ) ⟩ also differs significantly from what one would obtain by solving the largeN matrix integral without scaling β , which would just give eq. (4.48) as well. So we observe that our largeβ eigenvalue instanton is not only necessary to obtain the quenched logarithm, but it also qualitatively modifies the behavior of ⟨ Z ( β ) ⟩ from eq. (4.48) to eq. (4.49) at large β . Intuitively, this is because our instanton captures statistics of the tail of eigenvalues in the classically forbidden region, which always matter at large β .  \nIt only remains to apply our results for the logarithm of Z ( β ) to the calculation of entropies. The extreme β ≫ N is of particular interest, given the expectation that the quenched entropy  \nmust remain non-negative and that the annealed entropy must become arbitrarily negative. For the quenched entropy, applying eq. (1.9) to eq. (4.46),  \nS q ( β ) ≈ N ˆ R dx ˆ σ ∗ ( x ) e -β ( x -ˆ λ 0 ) , β ≫ N. (4.50)  \nAs mentioned below eq. (4.47), this term goes to zero as β → ∞ , giving a vanishing quenched entropy in the extremal limit. This is precisely the expected result for the non-degenerate eigenvalues that the leading spectral density captures, and demonstrates the success of our saddle-point analysis in implementing the replica trick for the quenched entropy. For the annealed entropy, applying eq. (1.8) to eq. (4.49),  \nS a ( β ) ≈ -N ( ̂ V ∗ ( ˆ λ ∗ ) -̂ V ∗ ( ˆ λ 0 )) , β ≫ N. (4.51)  \nFor well-defined matrix integrals, ̂ V ∗ generally grows unbounded off ˆ σ ∗ , and so does the magnitude of its gradient. Hence it follows from eq. (4.30) that S a ( β ) →-∞ as β →∞ , as expected. Explicit examples of matrix integral calculations of quenched and annealed entropies using our saddle-point analysis are given in section 5.', '4.6 Lessons for Gravity': "In connection with JT gravity, [19] interpreted single-eigenvalue instantons as the analogue of Liouville theory ZZ branes [132], and determinants like that in eq. (4.20) as FZZT branes [133, 134]. Our single-eigenvalue instantons require a reconciliation of these two notions of brane effects. The key point is that the determinant in eq. (4.17) is not inserted here as a probe with eigenvalue ˆ λ fixed by hand (unlike FZZT branes), but arises naturally inside the integral with ˆ λ actually being dynamical (like ZZ branes). The gravitational dual of our single-eigenvalue instantons should thus correspond to branes with the energy eigenbrane boundary conditions of [135, 136], but treated dynamically like the end-of-the-world (EoW) branes of [137].  \nAs understood in section 4.2, what gives rise to these branes is the largeβ back-reaction of Z ( β ) itself on the integral (with no additional probes), and their effect alone dominates ⟨ Z ( β ) m ⟩ at sufficiently large β to leading order. In other words, the brane is not a non-perturbative object we add by hand to complete the theory and probe ⟨ Z ( β ) m ⟩ ; the brane is already part of the theory and ⟨ Z ( β ) m ⟩ does not even make sense at large β without it. This narrative translates directly to the gravity side, and says that EoW eigenbranes are necessarily part of any quantum gravity theory with a matrix dual. Once again, this is not a statement about optional non-perturbative corrections; it is a leading-order effect without which the theory does not make sense. 59 The inclusion of such branes would not supplement the theory with a fixed discrete spectrum, much like there is no unique discrete spectrum in a matrix integral. Rather, their inclusion is needed for consistency with a random discrete spectrum and an ensemble description.  \nA gravity calculation of ⟨ Z ( β ) m ⟩ would correspond to the usual gravitational path integral P ( Z ( β ) m ) where one inserts m boundaries all with Z ( β ) boundary conditions. The inclusion of branes in the theory gives the option for geometries to end not only on Z ( β ) boundaries, but also on EoW eigenbranes. It is well-understood that the contribution of fully-connected m -boundary wormholes to P ( Z ( β ) m ) becomes important at large β , but that such geometries are off-shell and do not admit a semiclassical treatment [11, 19]. In fact, not even the off-shell contributions make sense at sufficiently large β , since the topological expansion breaks down when β becomes comparable to e S 0 . 60 Our findings suggest that a new semiclassical regime arises for β ∼ O ( e S 0 ), effectively capturing a resummation of the full topological expansion to all genus. In particular, the new largeβ saddle point to leading order at large e S 0 should be matched on the gravity side by a semiclassical geometry that dominates the gravitational path integral. Of course, such a geometry must necessarily involve eigenbrane boundaries, the preferred energy of which has to be dynamically determined by extremizing the action. 61 For the m = 1 case that more directly descends from the near-horizon geometry of a near-extremal black hole, the dominance of a brane configuration would correspond to the black hole throat ending at a finite distance, rather than becoming infinitely deep as β →∞ .  \nThe emergence of an effective gravitational description in terms of branes when some other perturbative topological expansion breaks down is certainly not unheard of in string theory [140]. An example which beautifully illustrates this in a duality between gravity and matrix integrals is that of [141-143]. This setting involves the computation of Wilson loops in N = 4 supersymmetric Yang-Mills using AdS/CFT at large N . For 1 / 2 BPS Wilson loops the boundary calculation can be performed by a GUE matrix integral. Holographically, these can be computed from the world-sheet of a fundamental string anchored to the loop asymptotically. When the Wilson loop is multiply wrapped with winding number k ≥ 1, the strings may coincide and interact among themselves. For k ∼ O (1), corrections are technically hard but in principle computable by including higher genus contributions. For k ∼ O ( N ), the topological expansion does not make sense anymore, and the many coincident strings are better described in terms of the dynamics of the D3-brane they end on. The classical action of the D3-brane solution is then seen to match precisely the result of the matrix dual. This is not merely analogous to what we are arguing for: a precise match between the winding number k and our inverse temperature β can be made and, as we show in section 5.1, our single-eigenvalue instantons exactly match the D-brane results of [141-143]. This provides a highly non-trivial consistency check of our saddle-point analysis of matrix integrals at large β , and the expectation of a dual description in terms of semiclassical geometries.  \nIn trying to reproduce our largeβ results with a purely gravitational calculation, some previous explorations of JT gravity with branes and its matrix dual are useful. The work of [135, 144] studied how N eigenvalues react to probe branes fixing O (1) of them. Instead, we are interested in how a dynamical brane eigenvalue reacts to the other N -1, which to leading order at large β are fixed. The half-wormholes or correlators between FZZT branes and Z ( β ) operators of [135, 144] would seem to provide the right tool if the eigenbrane energies are treated dynamically rather than  \nfixed by hand. A purely gravitational exploration of our largeβ instantons in terms of branes is currently underway [145].  \nThe discussion above applies to ⟨ Z ( β ) m ⟩ for any integer m ≥ 1, but the m → 0 replica trick for quenching requires extra care. That branes are indispensable to capture the largeβ behavior of the quenched entropy is already clear. 62 How to make sense of non-integer m in gravitational correlators, however, is not so clear. A natural approach one may be tempted to follow consists of computing moments for integer m ≥ 1, find their limiting behavior at large β and, if their m -dependence happens to be sufficiently simple, analytically continue them to m → 0. As we explain next, even if the analytic continuation of the largeβ moments is unique (which typically is not the case), this strategy is fundamentally flawed and bound to give wrong results.  \nThe basic reason the above does not work was already anticipated at the end of section 2.2 and in footnote 57: β →∞ and m → 0 are opposing limits, and the order in which they are taken does not commute. This is obvious from the saddle-point equations of motion, where m and β only appear in eq. (4.34) combined into mβ . However, it is much obscurer from any explicit expression for the moments of Z ( β ) like eq. (4.38), where m and β appear separately and also implicitly in ˆ λ ∗ . In other words, at the level of ⟨ Z ( β ) m ⟩ it is unclear how large β and small m compete with each other, and what approximations are allowed at large β given that the strict m → 0 limit must be taken first. This hidden competition in ⟨ Z ( β ) m ⟩ explains the failure of previous attempts in the literature at performing this replica trick even on the matrix side. Indeed, as [45] observed, the pursuit of this strategy leads to a pathological m → 0 limit even when the matrix integral is treated fully non-perturbatively. In particular, one finds that the continuation in m of the largeβ approximation for ⟨ Z ( β ) m ⟩ gives a power-law divergence as m → 0, in perfect inconsistency with the requirement that ⟨ Z ( β ) m ⟩ → 1 in the no-replica limit. A 'replica-scaling' on top of the DSL was proposed by [45] in order to obtain a sensible replica limit; unfortunately, by our general arguments above, we do not expect this to work.  \nThe question thus remains of how to make sense of ⟨ Z ( β ) m ⟩ for non-integer m in gravity. A promising direction is that of [40], where the semiclassical quotient construction of [4] was used to make the analytic continuation in m well-defined. By studying the replica trick at the level of the equations of motion, this gravitational prescription allows for a controlled treatment of the competition between m and β just like on the matrix side. As understood here though, it would be essential to include branes in the calculation of [40] in order to capture the largeβ regime.  \nFinally, let us comment on an additional subtlety that may arise in gravity in relation with replica symmetry. In gravity, the moments ⟨ Z ( β ) m ⟩ correspond to gravitational correlators involving all possible wormholes consistent with boundary conditions, connected and disconnected. On the matrix side it is simple to study these moments directly, without even having to think about breaking them down into more fundamental cumulants ⟨ Z ( β ) m ⟩ c . On the gravity side, however, the quantities one naturally computes are connected correlators involving connected wormholes, which correspond precisely to ⟨ Z ( β ) m ⟩ c . Hence to implement the replica trick in gravity one has to work harder, and figure out how to combine connected correlators in a way that can be analyt-  \nFigure 8 . Illustration of the relative contribution to ⟨ Z ( β ) m ⟩ from different connected correlator terms involving 〈 Z ( β ) k 〉 c for 1 ≤ k ≤ m . The plot considers m = 4, for which 〈 Z ( β ) 4 〉 is equal to the sum of all terms appearing on the legend. The figure shows the ratio of each of these terms relative to 〈 Z ( β ) 4 〉 as a function of β for the GUE with e S 0 ∼ N = 10 6 . The small and large β limits can be seen to follow eq. (4.52). The transition between 〈 Z ( β ) 4 〉 being dominated by the replica-symmetric contributions 〈 Z ( β ) 4 〉 c (fully connected) and ⟨ Z ( β ) ⟩ 4 (fully disconnected) occurs at β ∼ e 2 S 0 / 3 . This regime receives contributions of comparable order by partially-connected correlators which necessarily break replica symmetry.  \n<!-- image -->  \nd in m . This raises the question of which 〈 Z ( β ) k 〉 c contributions with 1 ≤ k ≤ m dominate ⟨ Z ( β ) m ⟩ , and how to even make sense of this question when m < 1. This problem was addressed by [40], where a gravitational ansatz to account for replica symmetry breaking (RSB) was also proposed.  \nWhether or not RSB plays a role in the gravitational implementation of the no-replica trick has remained a puzzle since the work of [11]. Our random matrix results seem to answer this in the affirmative: the m → 0 limit is fundamentally dominated by RSB configurations. To reach this conclusion, we may follow the strategy in [146] of decomposing ⟨ Z ( β ) m ⟩ into connected contributions 〈 Z ( β ) k 〉 c for each 1 ≤ k ≤ m as a function of m and β . The results of this exercise are illustrated in fig. 8. 63 In general, for any integer m ≥ 1, one easily observes the following behavior:  \nlim β →∞ ⟨ Z ( β ) m ⟩ = ⟨ Z ( β ) m ⟩ c , lim β → 0 ⟨ Z ( β ) m ⟩ = ⟨ Z ( β ) ⟩ m . (4.52)  \nIn gravitational terms, and without branes, this would mean that ⟨ Z ( β ) m ⟩ is dominated by a fully connected wormhole at low temperatures, but by a fully disconnected topology at high temperatures. These two topologies are the only ones with only Z ( β ) boundaries which can possibly accommodate full replica symmetry, which means that the transition between the two β limits must involve RSB with partially connected wormholes being important (see fig. 8). In the largeβ regime, eq. (4.52) would suggest approximating ⟨ Z ( β ) m ⟩ by ⟨ Z ( β ) m ⟩ c and analytically continuing to m → 0. But as we already explained, the β →∞ and m → 0 limits do not commute, and this approximation would indeed lead to pathological replica limit [46]. On the other hand, were we to  \n〈 e tZ ( β ) 〉 = exp ( ∞ ∑ k =1 t k k ! 〈 Z ( β ) k 〉 c ) .  \napproximate ⟨ Z ( β ) m ⟩ by ⟨ Z ( β ) ⟩ m , we would obtain an annealed answer. Since these are the only two replica-symmetric options without branes, an immediate corollary is that RSB is unavoidable in the m → 0 limit. With branes, however, it is plausible that replica symmetry could be restored. We leave an exploration of RSB and branes in the gravitational implementation of this replica trick to future work [145].", '5 Explicit Examples': 'We first do the exercise for Gaussian WD and BPS matrix ensembles without scaling, then for the soft edge of the universal Airy limit, and finally for the JT gravity matrix model in the DSL.', '5.1 Gaussian WD': "For a WD ensemble, the eigenvalue measure for the random matrix is the general one in eq. (3.10) with ν = 0, and according to eq. (3.12) the Hamiltonian is taken to be the random matrix itself. With a Gaussian potential V ( x ) = 1 2 x 2 , the exact unit-normalized solution to eq. (4.24) for the saddle ˆ σ ∗ is just a Wigner semicircle,  \nˆ σ ∗ ( x ) = N N -1 1 π β √ 2 β ( 1 -1 N ) -x 2 = 1 π β √ 2 β -x 2 + O ( N -1 ) , (5.1)  \nwhen the radicand is non-negative, and zero otherwise. The corresponding spectral curve off the support of ˆ σ ∗ is given by  \ny ( x ) = √ x 2 -2 β ( 1 -1 N ) = 1 β √ x 2 -2 β + O ( N -1 ) (5.2)  \nand the solution to the eigenvalue instanton condition, eq. (4.34), is  \nˆ λ ∗ = -√ 2 β ( 1 -1 N ) + ( mβ N ) 2 = -√ 2 β + ( mβ N ) 2 + O ( N -1 ) . (5.3)  \nAs it should, note that ˆ λ ∗ lies below the support of ˆ σ ∗ , and satisfies ˆ λ ∗ → inf supp ˆ σ ∗ as m → 0. Together, eqs. (5.1) and (5.3) completely characterize our single-eigenvalue instanton. Dropping pesky factors that are O ( N -1 ) and thus negligible at large N anyway, the leading effective potential for the spectral density in eq. (5.1) can be written  \n̂ V ∗ ( x ) = β 2 ( 1 -xy ( x ) + 2 log ( -x β -y ( x ) )) . (5.4)  \nOn the instanton eigenvalue from eq. (5.3) this evaluates to  \n̂ V ∗ ( ˆ λ ∗ ) = β ( 1 -log β / 2 2 + κ √ 1 + κ 2 -sinh -1 κ ) , κ ≡ mβ √ 2 β N , (5.5)  \nwhere κ is generally O ( m ) given that both N and β scale together. It only remains to obtain an expression for eq. (4.22) on the instanton solution,  \nZ [ˆ σ ∗ ; ˆ λ ∗ ] = e β √ 2 β √ 1+ κ 2 +2 N I 1 ( √ 2 β β ) √ 2 β β , (5.6)  \nwhere I 1 is the modified Bessel function of the first kind. Evaluating eq. (4.37), upon a convenient rewriting we arrive at the following instanton action  \nI ∗ ( m ) = -β ( κ √ 1 + κ 2 +sinh -1 κ ) -m N log ( 1 + 2 Ne -β √ 2 β √ 1+ κ 2 I 1 ( √ 2 β β ) √ 2 β β ) . (5.7)  \nThe final result for eq. (4.38) at the saddle-point level is thus  \n⟨ Z ( β ) m ⟩ ∗ = ( 1 + 2 Ne -β √ 2 β √ 1+ κ 2 I 1 ( √ 2 β β ) √ 2 β β ) m e N β ( κ √ 1+ κ 2 +sinh -1 κ ) . (5.8)  \nAt large β one may further reduce eq. (5.8) to  \n⟨ Z ( β ) ⟩ ∗ = e NF ( κ ) , F ( κ ) ≡ β κ √ 1 + κ 2 + β sinh -1 κ. (5.9)  \nFor the GUE case ( β = 2) and m = 1, this provides a remarkable consistency check of our formalism: eq. (5.9) identically matches the results of [141-143] for 1 / 2 BPS Wilson loops of N = 4 supersymmetric Yang-Mills at large N , 't Hooft coupling λ , and winding number k ∼ O ( N ), upon the identification β = k √ λ/ 2. Here F was obtained from the action of a single-eigenvalue instanton that arises at β ∼ O ( N ). In [141-143], F is the on-shell action of the D3-brane holographically dual to multiply wrapped Wilson loops at winding number k ∼ O ( N ). This is strongly suggestive that our largeβ instantons may generally be describable gravitationally in terms of dynamical branes. Furthermore, since eq. (5.8) is also valid for β ∼ O (1), our results provide an interpolation between the effective brane description that arises at k ∼ O ( N ), and the worldsheet description in terms of weakly interacting fundamental strings that is natural for k ∼ O (1).  \nWe would like to use eq. (5.8) to compare quenched and annealed entropies. For the former, using the replica trick from eq. (1.12) leads to a result of the universal form in eq. (4.46), which for this Gaussian WD matrix integral reads  \n⟨ log Z ( β ) ⟩ ∗ = log [ e √ 2 β β + N 2 I 1 ( √ 2 β β ) √ 2 β β ] ≈    log N + √ 2 β N β + β 4 β 2 , β ≪ N, √ 2 β 2 N β + 2 N √ 2 π (2 β ) 3 / 4 β -3 / 2 , β ≫ N. (5.10)  \nUsing eq. (1.9), the extreme β limits of the quenched entropy for this matrix integral are  \nS q ( β ) ≈    log N -β 4 β 2 , β ≪ N, 5 N √ 2 π (2 β ) 3 / 4 β -3 / 2 , β ≫ N. (5.11)  \nClearly, at small β we recover the expected maximal entropy for a system with N degrees of freedom, whereas at large β we obtain a consistently non-negative result that goes to zero in consistency with having a non-degenerate lowest eigenvalue. See fig. 9 for an illustration.  \nThese quenched results shall be compared to annealed ones, obtained by simply taking the logarithm of the m = 1 expectation value in eq. (5.8),  \nlog ⟨ Z ( β ) ⟩ ∗ = -I ∗ (1) ≈    log N + √ 2 β N β + β 4 β 2 , β ≪ N, β 2 2 N + β N log β N , β ≫ N. (5.12)  \nFigure 9 . Comparison between annealed (black) and quenched (red) entropies for the GUE using the instanton result from eq. (5.8). These entropies begin to differ at β ∼ O ( N 2 / 3 ), which determines the transition to a low-temperature regime where the largeβ instanton dominates both quantities and cannot be neglected (cf. fig. 10). Here N = 10 6 ; in the thermodynamic limit N →∞ , with an appropriate scaling S q ( β ) becomes non-differentiable at β ∼ O ( N 2 / 3 ), signaling a second-order phase transition. The extreme regimes β ≫ N and β ≪ N reproduce eqs. (5.11) and (5.13), and are respectively consistent with theorem 1 and the statistical expectation that e S asymptote to the total number of degrees of freedom N as β → 0.  \n<!-- image -->  \nwhich for the annealed entropy give the limiting behaviors  \nS a ( β ) ≈    log N -β 4 β 2 , β ≪ N, -β 2 2 N + β N log β N , β ≫ N. (5.13)  \nIn consistency with the general arguments of previous sections, we see that quenched and annealed results agree at small β , but strongly disagree at large β . In particular, eq. (5.13) gives a negatively divergent annealed entropy as β →∞ , exemplifying theorem 1. This is shown in fig. 9.  \nIt is also illustrative to compare these results with the na¨ıve answer one would obtain by neglecting the largeβ instanton. In other words, if one simply evaluates Z ( β ) on the WD semicircle saddle from eq. (3.37) (cf. using eq. (4.40)), one obtains  \nlog ⟨ Z ( β ) ⟩ na¨ıve ≈    log N + β 4 β 2 , β ≪ N, √ 2 β β -3 2 log β, β ≫ N, (5.14)  \nwhich for the entropy gives  \nS na¨ıve ( β ) ≈    log N -β 4 β 2 , β ≪ N, -3 2 log β, β ≫ N, (5.15)  \nAs previously pointed out, this gives the correct behavior at small β , but fails to capture the leading largeβ behavior of even the annealed entropy. In other words, our largeβ instantons are also important for annealed quantities, as we illustrate in fig. 10.  \nAs a matter of fact, the na¨ıve result in eq. (5.15) is actually the one that should be compared to the gravitational entropy in eq. (1.13). Indeed, eqs. (1.13) and (5.15) exhibit the exact same logarithmic divergence at large β for the bosonic s = 3 / 2 case, whereas the leading divergence of  \nFigure 10 . Free energies (left) and entropies (right) for the GUE with N = 100. Na¨ıve annealed results (blue) correspond to using the standard Wigner saddle from eq. (3.37) which does not take into account the largeβ instanton (cf. eqs. (5.14) and (5.15)). Exact annealed results (yellow) correspond to evaluating ⟨ Z ( β ) ⟩ exactly for N = 100 using the orthogonal polynomial techniques from appendix D.4. Annealed results (green) correspond to using the largeβ instanton form for ⟨ Z ( β ) ⟩ ∗ from eq. (5.8) at m = 1 (cf. eqs. (5.12) and (5.13). Quenched results (red) correspond to using the largeβ instanton form for ⟨ Z ( β ) m ⟩ ∗ from eq. (5.8) in the m → 0 replica trick (cf. eqs. (5.10) and (5.11)). Only quenched quantities preserve physical properties such as the monotonicity of the free energy and the corresponding non-negativity of the entropy. Importantly, our instanton is seen to also be necessary for capturing the correct β dependence of the annealed entropy at large β . In particular, although the na¨ıve finiteβ saddle gives the desired qualitative negativity, only the largeβ instanton gives a quantitative approximation to the exact result at large β .  \n<!-- image -->  \nour annealed result from eq. (5.13) is actually a quadratic power law. The reason the gravitational entropy behaves like eq. (5.15) is simply that eq. (1.13) is obtained to leading order at large e S 0 at finite β , which generally corresponds to a lower-bounded leading spectrum with a classically forbidden region below it. By working at large β ∼ e S 0 , our single-eigenvalue instanton is effectively accounting for corrections at higher O ( e -S 0 ), which capture the statistics of the tail of eigenvalues in the classically forbidden region. A gravitational calculation including O ( e -S 0 ) corrections would make the gravitational entropy reproduce eq. (5.13).", '5.2 Gaussian BPS': 'For BPS ensembles, the random matrix M is ( N +¯ ν ) × N and the Hamiltonian H can be constructed out of M as in eq. (3.13). 64 The parameter ¯ ν determines the number of zero eigenvalues, which corresponds to the degeneracy of the ground state. The measure for the positive eigenvalues of H is the general one in eq. (3.10). The natural largeN scaling for rectangular matrices corresponds to fixing their aspect ratio. Hence we let ν ≡ Nν and consider taking the largeN limit keeping ( N + ν ) /N = 1 + ν fixed. This amounts to keeping the proportion of zero eigenvalues against nonzero eigenvalues constant. The case in which the ground-state degeneracy ν is kept finite in the limit can be simply recovered by taking ν → 0.  \nWith a Gaussian potential for the singular values of M , the squared eigenvalues of H feel a linear potential, so in particular V ( x ) = 1 2 x in eq. (3.43). The solution ˆ σ ∗ to eq. (4.24) for a  \nsupport interval supp ˆ σ ∗ = [ a -, a + ] takes the Marˇcenko-Pastur form of eq. (3.44)  \nˆ σ ∗ ( x ) = N N -1 √ ( a + -x )( x -a -) 2 π β x . (5.16)  \nUnit normalization of the measure requires  \n√ a + = √ a -+2 √ β √ N -1 N , (5.17)  \nwhile direct evaluation of the principal value integral in eq. (4.24) leads to  \n[ a -,a + ] dy ˆ σ ∗ ( y ) x -y = N N -1 1 2 β ( 1 -√ a -a + x ) . (5.18)  \nMatching to the potential gives a simple relation between the endpoints and the ν parameter,  \n√ a -a + = 2 ν. (5.19)  \nThe explicit solution for the endpoints can be compactly written (cf. eq. (3.45))  \na ± = β N -1 N ( 1 ± √ 1 + N N -1 2 ν β ) 2 = β ( 1 ± √ 1 + 2 ν β ) 2 + O ( N -1 ) . (5.20)  \nAs usual, notice that for ν > 0 the lower end of the [ a -, a + ] interval is strictly positive, corresponding to a positive ∆ ≡ a -> 0 gap in the spectrum between the zero eigenvalues and the non-trivial eigenvalues that ˆ σ ∗ describes. The ν → 0 limit recovers the gapless ∆ → 0 spectrum associated to a ground-state degeneracy ¯ ν that does not scale with N . In this case eq. (5.16) degenerates into  \nˆ σ ∗ ( x ) ν = 0 = 1 2 π β √ 4 β -x x + O ( N -1 ) , (5.21)  \nwhich exhibits no gap and develops a pole at x = 0. As explained in appendix D.5, this is just the distribution that a Gaussian WD ensemble induces on squared eigenvalues (cf. eq. (D.58)) Here we study eq. (5.16) for general ν ≥ 0, from which eq. (5.21) is just the limiting case ν → 0. The spectral curve corresponding to eq. (5.16) for x > 0 off the support of ˆ σ ∗ is  \ny ( x ) = √ ( a + -x )( a --x ) 2 β x . (5.22)  \nThe corresponding eigenvalue instanton that solves eq. (4.34) is  \nˆ λ ∗ = a -+ a + -√ ( a + -a -) 2 +4 a -a + ξ 2 2 (1 -ξ 2 ) = 2 ( β + ν -√ β 2 + ν (2 β + νξ 2 ) ) 1 -ξ 2 + O ( N -1 ) , ξ ≡ 2 mβ N . (5.23)  \nAs required, note that ˆ λ ∗ lies below the support of ˆ σ ∗ but still on the allowed region of non-negative eigenvalues, i.e., 0 ≤ ˆ λ ∗ ≤ a -. In the limit ν → 0 one has ˆ λ ∗ → 0, consistent with the tightening  \nFigure 11 . Comparison between annealed (black) and quenched (red) entropies for the non-zero eigenvalues of the Gaussian BPS ensemble with β = 2 and ν = 1. The instanton solution from eqs. (5.16) and (5.23) is in this case evaluated numerically for the moments ⟨ Z ( β ) m ⟩ ∗ setting N = 10 6 . As in the GUE case from fig. 9, the soft edge that the gapped BPS ensemble exhibits gives a transition to the largeβ instanton regime also at β ∼ O ( N 2 / 3 ). The extreme regimes β ≫ N and β ≪ N reproduce eqs. (5.27) and (5.28).  \n<!-- image -->  \nof these bounds due to the spectral gap shrinking to zero size. Integrating eq. (3.34) gives the effective potential up to an unimportant constant as  \n̂ V ∗ ( x ) = -β xy ( x ) -1 2 ( a -+ a + ) tanh -1 √ a --x a + -x + √ a -a + tanh -1 √ a + ( a --x ) a -( a + -x ) . (5.24)  \nThe evaluation of the effective potential for the instanton in eq. (5.23) does not simplify in any particular way at finite ν and ξ even if N →∞ , so it is not worth quoting explicitly. Furthermore, the integral in eq. (4.22) does not take a simple form either for the spectral density in eq. (5.16) for general finite ν > 0. We nonetheless evaluate these results numerically in fig. 11.  \nThe BPS ensembles are interesting both for ν > 0 and ν = 0, corresponding to a positive gap and no gap in the spectrum, respectively. These are qualitatively distinct regimes and, while ν = 0 can be studied as the ν → 0 limit of the former, the converse is not true. Namely, it is not possible to study finite-gap ν > 0 results as an expansion about the gapless ν = 0 model. The reason is simply that the spectral density in eq. (5.21) has support all the way to zero, and so will capture any quantity computed by expanding eq. (5.16) about ν = 0. For instance, the allowed range for the eigenvalue instanton ˆ λ ∗ will remain to be just zero at any order in ν , which is inconsistent with the existence of a finite gap in the spectrum for any finite ν > 0. Another way of seeing this failure is to simply note that the terms in an expansion of the spectral density in eq. (5.16) about ν = 0 are not even integrable on the region of real support beyond second order.  \nSince a small-gap expansion makes no sense, we will only be interested in keeping finite ν > 0 or setting ν = 0 exactly. The latter is also better understood as a limit of the former, for otherwise at ν = 0 the instanton dynamics become trivial as its only allowed value is zero. In order to obtain explicit results, we will consider the largeβ limit which, after all, is the regime of interest for us. As was the case in the WD example, our quantities of interest in this limit are also here governed by the largeβ instanton. For the moments of Z ( β ), using the largeβ result from eq. (4.41) and  \nfurther expanding the exponent at large β we get  \n⟨ Z ( β ) m ⟩ ∗ = e NF ( ξ ) , F ( ξ ) ≡ a -+ a + 2 ( tanh -1 √ a -a + -1 2 ξ ) -ν log ( 8 νξ a + -a -) . (5.25)  \nThis result should be compared to the one in eq. (5.9) for the WD case. In particular, it would be interesting to understand if F in eq. (5.25) similarly matches the on-shell action of any known gravitational brane solutions.  \nAs emphasized previously, while eq. (5.25) makes sense for integer m ≥ 1, this largeβ expansion is incompatible with taking the m → 0 replica limit. Instead, directly applying eq. (4.47),  \n⟨ log Z ( β ) ⟩ ∗ = -β ∆+ O ( Ne -β ∆ ) , β ≫ N, (5.26)  \nwith ∆ = a -, and the exponentially small errors coming from the spectral density integral. It is worth keeping track of these as they emphasize that already the next-to-leading order corrections at largeβ are exponentially suppressed. This contrasts with the largeβ approximation for the annealed result that follows from eq. (5.25) for m = 1, where corrections appear at all orders in 1 /β . The entropy that results from quenching behaves as  \nS q ( β ) = O ( β e -β ∆ ) , β ≫ N, (5.27)  \nwhereas the annealed entropy yields  \nS a ( β ) = -νN log β N + O (1) , β ≫ N. (5.28)  \nOnce again, this verifies the general expectation from section 2.3 that the annealed entropy negatively diverges as β →∞ , whereas the quenched entropy remains non-negative.  \nThese are the results one obtains for the continuous spectrum of positive eigenvalues in BPS ensembles. Remember though that the spectrum of such ensembles also includes a discrete zero eigenvalue with degeneracy ν = νN . In the above we have analyzed the dynamical eigenvalues, and the entropic properties of their spectrum excluding the highly degenerate ground state from the partition function. However, if we want to study the full spectrum of the Hamiltonian altogether, we may want to also take into account the non-dynamical zero eigenvalues. Denoting the partition function with zero eigenvalues included by Z 0 , its moments are simply related to those of the partition functions above by  \n〈 Z 0 ( β ) m 〉 ∗ = ⟨ ( ν + Z ( β )) m ⟩ ∗ . (5.29)  \nAs expected from section 2.3, the presence of a nonzero number of discrete ground states, extensive or not in N , drastically changes the largeβ behavior of the annealed entropy from negatively divergent to non-negative and finite. In particular, for the extensive case above we obtain  \nS 0 a ( β ) = log ν + O ( β -ν ) , β ≫ N (5.30)  \nwhereas the quenched entropy becomes 65  \nS 0 q ( β ) = log ν + O ( β e -β ∆ ) , β ≫ N (5.31)  \nFigure 12 . Entropies for Gaussian BPS ensembles with β = 2 including the ground state with degeneracy ν = νN in the spectrum. On the left ν = 1, i.e., a non-extensive degeneracy giving a gapless BPS spectrum as N → ∞ ; on the right ν = 1, corresponding to an extensive degeneracy ν = N giving a gapped BPS spectrum as N →∞ . In these plots N = 10 6 . In contrast with the situation in fig. 11, in this case both the annealed (black) and quenched (red) entropies give log ν in the strict extremal limit. However, the approach to this asymptotic value is power-law for S a and exponential for S q , as quoted in eqs. (5.30) and (5.31), respectively. We demonstrate this discrepancy here by showing that as β increases, ( β/N ) ν +1 ( S ( β ) -log ν ) diverges monotonically for S a but goes to zero for S q . The scale of the vertical axis is irrelevant; on the horizontal axis, however, the entropies are generally observed to depart from each other at β ∼ O ( N/ν ).  \n<!-- image -->  \nHence we see that both entropies now asymptote to the same value capturing the degeneracy ν of the ground state of the Hamiltonian. This is precisely what happens when supersymmetry protects the ground-state energy, and why gravity despite computing an annealed quantity succeeds in computing the extremal entropy of BPS black holes. Nonetheless, as is clear from eqs. (5.30) and (5.31), the approach to the strict β → ∞ value of annealed and quenched entropies is very distinct qualitatively. In particular, the quenched entropy stays exponentially close to its extremal value up to β -1 scales of the order of the gap size ∆. In contrast, the annealed entropy departs from its extremal value polynomially in β -1 at order ν . Hence the behavior of the entropies and their derivatives at gap scales will differ significantly whether one is quenching or annealing, as illustrated in fig. 12.  \nThe attentive reader might wonder why the gravitational entropy in eq. (1.13), despite being an annealed entropy, exhibits the exponential suppression of a gapped quenched entropy. This lucky accident occurs for the same reason explained at the end of section 5.1. Namely, eq. (1.13) is a leading-order result at large e S 0 and finite β , the underlying spectrum of which has a classically forbidden region E < ∆. The statistics of eigenvalues in the gap region require O ( e -S 0 ) effects which would make the gravitational entropy in eq. (1.13) actually behave like eq. (5.30), and not like eq. (5.31).  \nTo address the gapless case it is preferable to reintroduce ν = Nν , and take the largeN limit keeping ν finite. The gap size is identically zero for ν = 0 at any order in N , but if ν > 0 then the gap size ∆ ∼ O ( N -2 ) at large N , i.e., the gapless spectrum only occurs in the strict N →∞ limit. Without a gap, the spectral density is given by eq. (5.21), for which the integral in eq. (4.38) is  \neasily computable and reads  \nˆ R dx ˆ σ ∗ ( x ) e -βx ν = 0 = ( I 0 (2 β β ) + I 1 (2 β β )) e -2 β β . (5.32)  \nBoth annealed and quenched quantities behave equivalently at leading order at small β and are governed by eq. (5.32). The smallβ results can be easily worked out to be  \n⟨ log Z ( β ) ⟩ ∗ ≈ log ⟨ Z ( β ) ⟩ ∗ ≈ log N -β β + 1 2 β 2 β 2 , β ≪ N, (5.33)  \nfor the logarithms, and  \nS q ( β ) ≈ S a ( β ) ≈ log N -1 2 β 2 β 2 , β ≪ N, (5.34)  \nfor the entropies. At large β the eigenvalue instanton once again dominates and the contributions from eq. (5.32) become subleading. The results above still apply so long as ν > 0, a parameter that can just be restored by resetting ν = ν /N . In other words, for ν > 0 one can just take the limit N → ∞ at finite ν in the expressions above to obtain gapless results. For quenched quantities, the m → 0 limit sends the instanton to the spectral edge ˆ λ ∗ → ∆, which does not depend on β , and then the largeN limit at finite ν sends ∆ → 0. The upshot is a quenched logarithm of the form of eq. (4.46) with ˆ λ 0 = 0 and a continuous part given by involving eq. (5.32). The resulting quenched entropy is clearly non-negative, as usual. For annealed quantities, however, m = 1 and the instanton saddle ˆ λ ∗ remains non-trivial and β -dependent. Scaling β with N , this gives rise to a regime in which the annealed logarithm becomes arbitrarily negative, which again results in an arbitrarily negative annealed entropy. This regime, however, only makes sense at nonzero ν > 0, whether or not ν scales with N .  \nThe case ν = 0 is strictly qualitatively different. According to eq. (3.11), this case corresponds to AZ ensembles with α = 1 or BPS ensembles with ¯ ν = 1 and β = 1 (a BPS ensemble with ¯ ν = 0 is just an AZ ensemble). In either case, the corresponding eigenvalue measure can be understood as coming from a WD ensembles upon squaring the random matrix. As such, the location of zero in the spectrum is not characterized by any particular eigenvalue, but rather by whichever eigenvalue is closer to zero in the original WD matrix. In other words, the statistics that govern the lowest eigenvalue in this case are not edge statistics, but rather bulk statistics of the eigenvalues that were squared. For completeness, we note that ν = -1 / 2 is also allowed for AZ ensembles with α = 0. Our main motivation for studying BPS ensembles in this paper was to explore the gapped case, so we will not discuss these cases any further. A thorough exploration of the spectral statistics of ensembles which include ν = 0 can be found in [37].', '5.3 Airy Edge': 'As explained in section 3.6.3, the statistics of soft edges are universally captured by the Airy model in eq. (3.38). This includes WD and gapped BPS ensembles upon edge scaling, not only for the Gaussian cases studied in sections 5.1 and 5.2, but for any matrix potential. Recalling eq. (3.39), the leading spectral density of the Airy model is supported on x > 0 and reads  \nˆ ρ ∗ ( x ) = √ x π . (5.35)  \nThe corresponding spectral curve along x < 0 is simply  \ny ( x ) = √ -x. (5.36)  \nBy eq. (4.34) we obtain the instanton solution  \nˆ λ ∗ = -κ 2 , κ ≡ mβ 2 e S 0 , (5.37)  \nwhere since β can be absorbed in the zooming variable to leading order, we have simply chosen to fix β = 2. Once again, the edge of the spectrum ˆ λ ∗ → 0 is approached as m → 0. The effective potential along x < 0 can be obtained by integrating eq. (3.34), which gives  \n̂ V ∗ ( x ) = 4 3 ( -x ) 3 / 2 . (5.38)  \nFor the instanton in eq. (5.37) this evaluates to  \n̂ V ∗ ( ˆ λ ∗ ) = 4 3 κ 3 . (5.39)  \nThe canonical partition function from eq. (4.22) yields  \nZ [ˆ σ ∗ ; ˆ λ ∗ ] = e βκ 2 + e S 0 2 √ π β 3 / 2 . (5.40)  \nWith a convenient rearrangement of terms, the resulting instanton action is  \nI ∗ ( m ) = -2 3 κ 3 -me -S 0 log ( 1 + e S 0 -βκ 2 2 √ π β 3 / 2 ) . (5.41)  \nThe final result for eq. (4.42) in the Airy model is  \n⟨ Z ( β ) m ⟩ ∗ = ( 1 + e S 0 -βκ 2 2 √ π β 3 / 2 ) m e 2 3 e S 0 κ 3 . (5.42)  \nUsing the replica trick from eq. (1.12) on eq. (5.42) gives the extreme limits  \n⟨ log Z ( β ) ⟩ ∗ = log ( 1 + e S 0 2 √ π β 3 / 2 ) ≈    S 0 -3 2 log β +2 √ πe -S 0 β 3 / 2 , β ≪ e S 0 , e S 0 2 √ π β -3 / 2 , β ≫ e S 0 . (5.43)  \nFor the quenched entropy, these give  \nS q ( β ) ≈    S 0 -3 2 log β, β ≪ e S 0 5 e S 0 4 √ π β -3 / 2 , β ≫ e S 0 (5.44)  \nIn consistency with quenched universality, this largeβ result matches precisely that of the GUE, i.e., the Gaussian WD result for β = 2 in eq. (5.11). In this case though, the entropy grows unbounded as β → 0, as expected from the non-normalizability of the spectral density of scaled matrix integrals (see fig. 13). For the annealed entropy,  \nS a ( β ) ≈    S 0 -3 2 log β, β ≪ e S 0 -e -2 S 0 6 β 3 , β ≫ e S 0 (5.45)  \nFigure 13 . Comparison between annealed (black) and quenched (red) entropies for the Airy model using the instanton result from eq. (5.42), and setting e S 0 = 10 6 . The transition to the instanton-dominated regime occurs at β ∼ e 2 S 0 / 3 . Under the identification of the scaling parameter N ∼ e S 0 , this is consistent with the examples in figs. 9 and 11 given that both fall under the same soft edge universality class of the Airy model. The extreme regimes β ≫ N and β ≪ N reproduce the expectations from eqs. (5.44) and (5.45). In particular, the entropies grow unbounded as β → 0, as expected from the N →∞ limit of edge scalings.  \n<!-- image -->  \nOnce again, as expected, smallβ behaviors match, while for large β the annealed entropy diverges negatively. Finally, a na¨ıve application of eq. (4.40) to eq. (5.35) yields  \nS na¨ıve ( β ) ≈    S 0 -3 2 log β, β ≪ e S 0 -3 2 log β, β ≫ e S 0 (5.46)  \nAs was also observed in section 5.1, this na¨ıve result reproduces the largee S 0 gravitational entropy from eq. (1.13) for the bosonic s = 3 / 2 case, but fails to capture the O ( e -S 0 ) corrections that our largeβ instanton accounts for in eq. (5.45).', '5.4 JT Gravity': "The leading spectral density for pure JT gravity is given by eq. (2.24). With the desired scaling from eq. (3.51) for DSL matrix integrals,  \nˆ ρ ∗ ( x ) = 1 2 π 2 sinh ( 2 π √ x ) , (5.47)  \nwhere x > 0. Near x = 0, eq. (5.47) reduces to the Airy soft edge from eq. (5.35), as expected from universality. The corresponding spectral curve along x < 0 is  \ny ( x ) = 1 2 π sin ( 2 π √ -x ) . (5.48)  \nThis spectral curve is oscillatory and bounded to 2 πy ( x ) ∈ [ -1 , 1], so eq. (4.34) would seem to admit either infinitely many or no instanton solutions. This is a manifestation of the fact that eq. (5.47) alone does not actually completely define a matrix integral, but requires a nonperturbative completion [19, 147]. However, regardless of the non-perturbative completion, we know the only correct instanton solution for m = 0 (or β = 0) is eq. (4.26), i.e., ˆ λ ∗ must lie at the x = 0 edge of the spectrum. Hence independently of how one completes the model nonperturbatively, we can safely just focus on instanton solutions which are continuously connected  \nto the edge of the spectrum as m → 0. In other words, when solving eq. (4.34) we want to restrict our solutions to a domain connected to x = 0 where the spectral curve y is injective. Since y attains its first maximum at x b = -1 / 16, our desired domain is x ∈ [ x b , 0]. The maximum here is y ( x b ) = 1 / 2 π , and thus the only possible solutions without a non-perturbative completion are  \nˆ λ ∗ = -( θ 2 π ) 2 , sin θ ≡ π mβ e S 0 ≤ 1 , 0 ≤ θ ≤ π 2 , (5.49)  \nwhere we have used β = 2 for ordinary JT (cf. footnote 46). The parameter constraint above implies that these real instantons can only capture finite m ≥ 1 moments if β is not too large. 66 Fortunately, this still allows one to probe β ∼ O ( e S 0 ), even if there is a finite upper bound on β/e S 0 . Furthermore, this restriction does not affect the replica trick at all: as the limit m → 0 is taken, the real instanton exists for arbitrarily large β . In other words, the results we obtain for quenched quantities are valid at any β , including the limit β →∞ .  \nBy eq. (3.34), integrating eq. (5.48) we find the leading effective potential on x < 0 to be 67  \n̂ V ∗ ( x ) = 1 2 π 3 ( sin ( 2 π √ -x ) -2 π √ -x cos ( 2 π √ -x )) . (5.50)  \nOn the instanton from eq. (5.49), this evaluates to  \n̂ V ∗ ( ˆ λ ∗ ) = 1 2 π 3 (sin θ -θ cos θ ) . (5.51)  \nThe canonical partition function from eq. (4.22) gives  \nZ [ˆ σ ∗ ; ˆ λ ∗ ] = e β ( θ 2 π ) 2 + e S 0 + π 2 /β 2 √ π β 3 / 2 . (5.52)  \nThe second term is the familiar 1-loop exact canonical partition function from eq. (2.24) that JT gravity yields on the disk. The first term is a novel contribution from the largeβ instanton. When θ ≪ 1, corresponding to mβe -S 0 ≪ 1, this contribution behaves as  \ne β ( θ 2 π ) 2 = e γ 2 + O ( θ 3 ) , γ ≡ mβ 3 / 2 2 e S 0 (5.53)  \nThis approximation clearly holds at large S 0 for any β < O ( e 2 S 0 / 3 ), but also for arbitrarily large β in the m → 0 replica limit. Gravitational effects that are doubly non-perturbative in S 0 are typically expected to give small corrections of O ( e -e S 0 ); eq. (5.53) has a different flavor. From its form we see that in eq. (5.52) the instanton contribution grows exponentially in γ 2 , while the disk contribution goes to zero as 1 /γ . The two terms exchange dominance when γ ∼ O (1) or β ∼ O ( e 2 S 0 / 3 ) at finite m > 0. In particular, in this case the instanton contribution becomes dominant for β ≳ O ( e 2 S 0 / 3 ). Hence there is a parametrically large regime of values of β for which the instanton in eq. (5.49) exists and dominates eq. (5.52), namely,  \nO ( e 2 S 0 / 3 ) ≲ β ≲ O ( e S 0 /m ) . (5.54)  \nFigure 14 . Comparison between annealed (black) and quenched (red) entropies for JT gravity using the instanton result from eq. (5.56), and setting e S 0 = 10 6 . Both entropies match precisely the corresponding ones in the Airy model from fig. 13 at large β , and only differ from those at small β ∼ O (1) when S ∼ O ( S 0 ). At this point S a and S q are already indistinguishably uninteresting and dominated by bulk statistics.  \n<!-- image -->  \nOnce again, the upper bound above matters for finite m ≥ 1 moments, but disappears as m → 0. The instanton action that results from eqs. (5.51) and (5.52) can be written  \nI ∗ ( m ) = -1 2 π 3 [ θ cos θ -( 1 -θ 2 2 ) sin θ ] -me -S 0 log   1 + e S 0 + π 2 /β -β ( θ 2 π ) 2 2 √ π β 3 / 2   , (5.55)  \nwhich leads to the following final result for the moments in JT gravity:  \n⟨ Z ( β ) m ⟩ ∗ =   1 + e S 0 + π 2 /β -β ( θ 2 π ) 2 2 √ π β 3 / 2   m e e S 0 2 π 3 [ θ cos θ -( 1 -θ 2 2 ) sin θ ] . (5.56)  \nThe term in the exponential above agrees with the leading m = 1 contribution that [148, 149] finds in what they call the 't Hooft expansion of log ⟨ Z ( β ) ⟩ . In the m → 0 replica limit, this JT instanton result precisely matches the Airy form of eq. (5.43) at large β , as expected from section 4.5. As for the quenched entropy, the results for JT gravity exactly match eq. (5.44) in both extreme β regimes. For integer m ≥ 1 moments, given the largeβ bracket from eq. (5.54), it is of interest to specifically consider the β ∼ O ( e 2 / 3 S 0 ) regime. The leading form of eq. (5.56) this gives is  \n⟨ Z ( β ) m ⟩ ∗ ≈ ( e γ 2 / 3 + me -2 γ 2 / 3 4 √ πγ ) m , β ∼ O ( e 2 / 3 S 0 ) . (5.57)  \nwhere γ was defined in eq. (5.53). This controlled largeβ regime still allows for a successful implementation of the replica trick, i.e., eq. (5.57) still yields the Airy result from eq. (5.44). Intuitively, this is because β ∼ O ( e 2 / 3 S 0 ) precisely corresponds to the transition between the disk and instanton dominating; had we trivialized the disk contribution, the replica trick would fail as discussed in section 4.6. As for the annealed entropy, both eqs. (5.56) and (5.57) also yield consistent results which precisely match the Airy expressions we obtained in eq. (5.45). The only differences with the Airy model occur at small β , a regime no longer dominated by the universality of edge statistics. This is illustrated in fig. 14.  \nAs it turns out, β ∼ O ( e 2 / 3 S 0 ) is also a controlled regime in the topological expansion of JT gravity. The first attempt at implementing the quenched replica trick in gravity made use of precisely this largeβ limit [11]. A difficulty [11] faced is that the growth in m of these moments is O ( e m 3 ), i.e., much faster than exponential. As a result, an analytic continuation to m → 0 would be on general grounds highly non-unique. In our case, what makes the analytic continuation unique is the fact that eq. (5.57) is a saddle-point result obtained by solving equations of motion which hold for all m (cf. the strategy of [40]). As a result, eq. (5.57) actually applies automatically for any real m ≥ 0. Because eq. (5.57) captures a transition to dominance of branes where the topological expansion should still make sense, we expect this result to provide a unique analytic continuation for the analysis of [11] which should succeed in implementing the replica trick.  \nTo relate our results here to those of [135, 144] one must treat branes dynamically rather than as probes. A plausible expectation is that the θ -dependent terms in eq. (5.56) may be sourced by halfwormhole contributions to the gravitational path integral. Such effects are generally studied offshell, which makes them practically inaccessible beyond 2-dimensional models of quantum gravity like JT. However, our treatment here is strictly to leading order at large e S 0 in terms of new saddle points that arise at large β ∼ O ( e S 0 ). This strongly suggests that a semiclassical description of these half-wormholes may actually arise in the appropriate regime. In other words, this kind of non-perturbative effects would seem to admit an on-shell analysis which would easily generalize to higher-dimensional theories of gravity. A purely semiclassical calculation of these effects in quantum gravity will be pursued in [145].", 'Acknowledgements': "It is a pleasure to thank Netta Engelhardt, Daniel Harlow, Tom Hartman, Gary Horowitz, Luca Iliesiu, Finn Larsen, Samuel Leutheusser, Henry Maxfield, Marija Tomaˇsevi'c, Misha Usatyuk, Nico Valdes-Meller, Evita Verheijden, and Wayne Weng for stimulating discussions. This research was supported by the Templeton Foundation via the Black Hole Initiative (award number 62286) and the MIT Center for Theoretical Physics. The opinions expressed in this publication are those of the author and do not necessarily reflect the views of the John Templeton Foundation.", 'A Near-Extremal Black Holes': 'This section gives a general overview of the current understanding of the gravitational entropy of near-extremal black holes, commenting on the results and expectations in semiclassical, perturbative, and non-perturbative temperature regimes relevant to this paper.', 'A.1 Semiclassical Entropy': 'The approach to black hole thermodynamics using the Euclidean gravitational path integral was pioneered at a semiclassical level by Gibbons-Hawking [150], who using eq. (1.1) successfully re-  \nproduced the famous Bekenstein-Hawking entropy formula [1, 151],  \nS BH ( β ) ≡ A ( β ) 4 G N , (A.1)  \nwhere A ( β ) is the horizon area of the black hole at a given temperature. Generically, this area remains non-zero even in the extremal limit β → ∞ , and grows linearly to leading order at low temperature. As a result, the general form of eq. (A.1) at large β is 68  \nS BH ( β ) = S 0 + # β + O ( β -2 ) , (A.2)  \nwhere S 0 is often referred to as the extremal entropy. Addressing the precise microscopic meaning of this entropy and relation or lack thereof to counting extremal states is one of the goals of this paper. A standard statistical interpretation of S 0 as counting microstates would imply a degeneracy of the ground state given by e S 0 . Since the extremal black hole area is typically macroscopically large in Planck units, e S 0 would in particular be a very large number. Indeed, the result in eq. (A.2) itself should be understood as an approximation S ( β ) ≈ S BH ( β ) capturing a largeβ expansion of just the classical, leading largeS 0 contribution to Z in eq. (1.1).  \nThe Strominger-Vafa enumeration of BPS states provided a celebrated confirmation of this result for a class of extremal black holes with supersymmetry in string theory [12, 17]. On symmetry grounds, such a large degeneracy is not surprising for black holes which preserve supersymmetry. However, that extremal black holes with no underlying symmetries could generally have such a high ground-state degeneracy of O ( e S 0 ) has historically been a puzzling possibility [152-159]. 69  \nThere is actually a simple argument for why one should be wary of the semiclassical analysis near extremality. Semiclassically, as the temperature of a black hole goes to zero, so does its heat capacity. Hence at low temperatures the emission of a typical Hawking quantum would drastically change the black hole temperature. But this would contradict the adiabaticity assumption of the semiclassical derivation of Hawking radiation [57, 152, 157]. On dimensional grounds, this argument implies that the semiclassical result in eq. (A.1) should only be trusted for β ≪ S 0 .', 'A.2 Perturbative Corrections': "The near-extremal regime we are interested in corresponds to β ≫ S 0 , where the new large parameter β implies that quantum effects that are subleading at large S 0 may no longer be negligible. Because eq. (1.1) gives the thermal entropy linearly in log Z , leading quantum corrections in a largeS 0 expansion are logarithmic in 1-loop determinants. The computation of these has been tackled through a variety of methods for a plethora of black hole spacetimes, with and without supersymmetry, and at and near extremality [51-56, 58, 61-63, 65, 66, 162-169]. In general this turns out to be an extremely technical and laborious endeavor, with the tools at one's disposal strongly depending on the specific setting and on which contributions among those from matter  \nfields, gravitons, or even strings, are being accounted for. Fortunately though, the functional form of these logarithmic corrections happens to be quite universal and straightforward to characterize.  \nWith a single exception to be described shortly, all 1-loop contributions near extremality yield logarithmic corrections to the extremal entropy S 0 of the form [51-55, 162-169]  \nS log ≡ c log log S 0 . (A.3)  \nwhere c log is an O (1) factor that only depends on the number of massless fields (massive fields contribute at subleading orders). This is the type of logarithmic correction generated by all nonzero modes (both bosonic and fermionic, including gravitons), and also every zero mode but one kind [51-53, 169]. 70 The only exception is a single family of relevant metric tensor deformations which become strongly coupled at low temperatures. These are exact zero modes of the extremal metric that have to be included, but the integral over which is formally divergent. Regulating their contribution requires flowing off the strict β → ∞ limit, which gives as a result a temperaturedependent logarithmic correction of the form [56, 58, 61-63, 65, 66, 170] 71  \nS Sch ( β ) ≡ -s log β, (A.4)  \nwhere s ≥ 0 is a rational number which is found to only depend on supersymmetry considerations. In particular, s = 0 for black holes which are BPS at extremality [58], but s > 0 for any black hole which is not [57]. For instance, in purely bosonic theories one finds s = 3 / 2.  \nAltogether, combining the near-extremal results from eqs. (A.2) to (A.4), we obtain  \nS ext ( β ) = S BH ( β ) ︸ ︷︷ ︸ classical + S log + S Sch ( β ) ︸ ︷︷ ︸ 1-loop = S 0 + c log log S 0 -s log β + # β , (A.5)  \nwhich gives the approximation S ( β ) ≈ S ext ( β ) at large β of eq. (1.1) with Z capturing up to 1-loop order in a largeS 0 expansion. In fact, as we will review in appendix A.3, the result that gives rise to eq. (A.4) is 1-loop exact, meaning that the β dependence of this approximation should be reliable to all orders in perturbation theory at large S 0 . Non-perturbative effects of O ( e -S 0 ), about which the result in eq. (A.5) has nothing to say, will be addressed shortly.  \nClearly, the corrections to the extremal entropy coming from eq. (A.3) cannot possibly resolve the degeneracy puzzle raised above. Since these are temperature-independent, we will hereon absorb them in S 0 itself by a redefinition of its Bekenstein-Hawking value in eq. (A.2) to S 0 + S log . Hence, in the absence of symmetry principles, the only quantum correction which can possibly counter the large contribution from S 0 to the entropy near extremality must come from eq. (A.4). In beautiful agreement with expectation, S Sch leaves unchanged the independently confirmed extremal entropy of BPS black holes ( s = 0), but provides the desired negative contribution to compete against S 0 for those which are not BPS at extremality ( s > 0). We refer to near-extremal black holes which are BPS at extremality as near-BPS black holes.  \nFor near-extremal black holes which are not near-BPS, remarkably not only does S Sch reduce the near-extremal entropy, but it also would make it negative at very low temperatures starting at β ∼ O ( e S 0 ), 72 a seemingly nonsensical result. Different potential resolutions to this pathology have been suggested in the literature. Since for non-near-BPS black holes Z ( β ) → 0 as β → ∞ , one option is that when β ∼ O ( e S 0 ), some other object must dominate the partition function of the theory [64, 65]. 73 From the viewpoint of statistical physics, this is tantamount to saying that very near-extremal black holes which are not near-BPS do not exist. Alternatively, at such large β ∼ O ( e S 0 ) one enters a regime where corrections which are non-perturbative in S 0 become important. These could potentially be able to tame the negativity induced by eq. (A.4) so as to restore a non-negative result for eq. (1.1). Put differently, despite perturbative exactness, nonperturbative effects would constrain the range of validity of eq. (A.5) to e S 0 ≫ β ≫ S 0 [57, 61, 65], hence leaving room for non-perturbative effects of O ( e -S 0 ) to prevent a pathology in eq. (1.1). 74  \nHowever plausible, these ideas assume that the quantity S computed by eq. (1.1) in an effective theory of gravity ought to behave like a standard thermodynamic entropy and thus be non-negative. A central goal of this paper is to clarify that this assumption is incorrect: in effective quantum gravity S is in fact not the thermal entropy of any quantum system. Moreover, we show that S ( β ) must diverge negatively in the β →∞ limit for non-near-BPS black holes whenever the Euclidean gravitational path integral is used to compute Z ( β ). To identify what kind of quantity S is, we now turn to a brief review of the origin of the quantum effects that give rise to S Sch in eq. (A.4).", 'A.3 Non-perturbative Corrections': 'As mentioned above eq. (A.4), the contribution S Sch to the near-extremal entropy comes from metric fluctuations of the black hole geometry which are nearly-zero modes near extremality and become exact zero modes in the extremal limit. These correspond to large diffeomorphisms supported on the throat of near-extremal black holes, a near-horizon region which is found to universally develop an AdS 2 factor with an enhanced SL (2 , R ) symmetry in the extremal limit [171-173]. These strongly coupled gravitational dynamics in the throat are governed by a Schwarzian theory of reparameterization of the AdS 2 asymptotics or its supersymmetric generalizations [73-75], which can be argued to arise on general grounds from the emergent conformal symmetry at extremality and its breaking for β < ∞ [56-59, 61, 170]. These Schwarzian theories turn out to be solvable quantum mechanical models with path integrals which are 1-loop exact, and the computations of which have been carried in various contexts [67, 76-81].  \nMost often, a Schwarzian theory is arrived at by performing a dimensional reduction in the throat to the emergent AdS 2 factor. 75 The 2-dimensional theory one obtains is JT gravity [41, 42] or a generalization thereof, where a dilaton field on the AdS 2 geometry captures fluctuations  \nof the volume of the transverse dimensions. The non-trivial dynamics of this JT dilaton are well understood to reduce to an asymptotic boundary mode governed by a Schwarzian theory [178, 179]. As a result, JT gravity theories themselves are frequently taken as the starting point for the study of the low-energy dynamics of near-extremal black holes [44, 80, 81, 172, 173, 176-185].  \nIt is natural to treat the JT theory one obtains this way as a quantum gravitational theory in its own right. Formulated as a Euclidean gravitational path integral, such a theory involves integrating over all possible 2-dimensional Riemannian manifolds with prescribed boundary conditions on AdS 2 boundaries. 76 The JT dilaton path integral imposes a constant negative curvature everywhere, thus making the gravity path integral run only over rigid hyperbolic spaces. These are Riemann surfaces which can be topologically classified, and whose moduli space within each topology class is finitedimensional. Hence the geometric part of the JT path integral is well-defined, and boils down to a sum over all topologies consistent with boundary conditions, weighted by the finite Weil-Petersson volumes of corresponding moduli spaces. In addition, each topology comes with a topological suppression in O ( e -S 0 ) by the Euler characteristic of the corresponding Riemann surface. In the general settings of interest in JT gravity, one allows for spaces to have at least one asymptotic boundary. This way, on top of every geometric contribution, there also appear Schwarzian modes along every boundary governing the asymptotic dynamics of the JT dilaton. All such contributions are, once again, known exactly in S 0 . As worked out in [19], the upshot is a JT gravitational path integral which is explicit at large S 0 , perturbatively exact as a quantum expansion in 1 /S 0 , and non-perturbatively asymptotic as a topological expansion in powers of e -S 0 .  \nThis remarkably powerful arena captures precisely the O ( e -S 0 ) effects that one needs to potentially probe non-perturbatively small temperatures with β ∼ O ( e S 0 ). 77 Since, the expansion at small e -S 0 is an asymptotic series, the competition between the topological suppression and large β ∼ O ( e S 0 ) has to be handled carefully. In particular, the asymptotic expansion parameter that becomes effective at large β turns out to be the combination β s e -S 0 , which ceases to be small when β ∼ O ( e S 0 /s ). According to eq. (A.5), this is actually also precisely the non-perturbatively scaling of β for the onset of the negativity of S ( β ). A reliable analysis would thus seem impossible given that we begin to lose non-perturbative control as soon as we reach the regime of interest. In fact, the partial sums of the asymptotic series still provide tight enough bounds to show that the pathological negative behavior of S ( β ) persists non-perturbatively, as first shown by [11]. In other words, one reliably concludes that the quantity in eq. (1.1) with Z ( β ) computed by the JT gravitational path integral genuinely becomes negative at sufficiently low temperatures. A central goal of this paper is to establish that this negativity occurs universally for non-near-BPS black holes and holds non-perturbatively in effective quantum gravity. A key insight for this realization follows from the non-perturbative correspondence between JT gravity theories and random matrix ensembles [19], which we make use of extensively.', 'B Proofs for Negative Entropies': "In this section we prove theorem 1, as well as three additional requisite results. The notation relevant to the results and proofs below is defined in section 2.3.  \nWe first establish the general region of convergence of Z g ( z ) → 0 as | z | → 0 for a piecewisecontinuous function g with a well-defined Laplace transform Z g given by eq. (2.22). This result is needed for understanding how a canonical partition function behaves at large β .  \nLemma 1. If g is locally integrable and of exponential type such that Z g in eq. (2.22) exists, then lim β →∞ Z g ( βe iα ) = 0 for any α ∈ ( -π 2 , π 2 ) .  \nProof. Assuming Z g ( βe iα ) exists, it can be written  \nZ g ( βe iα ) = lim ϵ → 0 + ˆ ∞ ϵ dE g ( E ) e -βe iα E .  \nBy assumption of exponential type, | g ( E ) | ≤ O ( e cE ) for finite c > 0 and sufficiently large E . Hence one can pick a finite β ∗ > 0 satisfying β ∗ cos α ≥ c such that Z g ( βe iα ) exists for all β cos α with β > β ∗ . Additionally, by integrability of g , the function f ( E ) ≡ | g ( E ) | e -β ∗ cos αE is also integrable and dominates the integrand above, i.e., | g ( E ) e -Ee iα | ≤ f ( E ) for all β ≥ β ∗ and E ≥ ϵ . The dominated convergence theorem can thus be used to take the β →∞ limit inside the integral,  \nlim β →∞ lim ϵ → 0 + ˆ ∞ ϵ dE g ( E ) e -βe iα E ≤ ˆ ∞ 0 + dE g ( E ) lim β →∞ e -βe iα E = 0 . (B.1)  \nNote that this result allows for g ( E ) to diverge as E → 0 + . For instance, suppose g ( E ) ⊃ E p . The Laplace transform of this term exists if and only if p > -1, and gives p ! /β p +1 , which indeed tends to zero as β → ∞ for all allowed p > -1. As another example, if g ( E ) ⊃ -log E , the Laplace transform of this term gives ( γ +log β ) /β , which again goes to zero in the limit.  \nThe result that follows is needed to specifically address the behavior of the entropy at large β .  \nLemma 2. Let f ( e x ) = ∂ x log( -log Z g ( e x )) . If lim x →∞ f ( e x ) = λ , then for any small ϵ > 0 , there exists c ± , β ∗ > 0 such that for all β > β ∗ ,  \ne -c -β λ -ϵ > Z g ( β ) > e -c + β λ + ϵ .  \nProof. If lim x →∞ f ( e x ) = λ , then for any ϵ > 0 there exists a large 0 < x ∗ < ∞ such that  \nλ -ϵ < f ( e x ) < λ + ϵ, ∀ x > x ∗ . (B.2)  \nUnder a semidefinite integral from x ∗ to x > x ∗ , this gives  \n( λ -ϵ )( x -x ∗ ) < log ( log Z g ( e x ) log Z g ( e x ∗ ) ) < ( λ + ϵ )( x -x ∗ ) , (B.3)  \nExponentiating twice and letting β = e x , this leads to  \nexp ( log Z g ( β ∗ ) ( β/β ∗ ) λ -ϵ ) > Z g ( β ) > exp ( log Z g ( β ∗ ) ( β/β ∗ ) λ + ϵ ) , (B.4)  \nwhere we used that log Z g ( β ) < 0 at large β = β ∗ , since lim β →∞ Z g ( β ) = 0 by lemma 1. Letting c ± ≡ -log Z g ( β ∗ ) β ∗ ( λ ± ϵ ) and noting that c ± > 0, we obtain the desired result.  \nWe now prove a mild generalization of the initial value theorem for Laplace transforms:  \nLemma 3. If g is of exponential type and has Laplace transform Z g ( β ) , then there is a constant 0 < ν ≤ 1 such that lim E → 0 + E 1 -ν g ( E ) = ˆ g (0 + ) < ∞ exists and  \nlim β →∞ β ν Z g ( β ) = Γ( ν ) ˆ g (0 + ) .  \nFor ν = 1 , this is known as the initial value theorem.  \nProof. For Z g to exist, g must be locally integrable, which requires lim E → 0 + g ( E ) E = 0. Hence there must exist 0 < ν ≤ 1 such that  \ng ( E ) = ˆ g ( E ) E 1 -ν and lim E → 0 + ˆ g ( E ) = g (0 + ) < ∞ . (B.5)  \nUsing this and changing variables to ξ ≡ βE , the integral form of Z g ( β ) becomes  \nZ g ( β ) = β -ν ˆ ∞ 0 dξ e -ξ ˆ g ( ξ/β ) ξ 1 -ν , (B.6)  \nHence for g of exponential type, using dominated convergence as in lemma 1,  \nlim β →∞ β ν Z ( β ) = lim β →∞ ˆ ∞ 0 dξ e -ξ ˆ g ( ξ/β ) ξ 1 -ν = ˆ g (0 + ) ˆ ∞ 0 dξ e -ξ ξ 1 -ν = Γ( ν ) ˆ g (0 + ) . (B.7)  \nFinally, the main theorem and its proof go as follows:  \nTheorem 1. Let S a ( β ) be the annealed entropy defined in eq. (1.8) , with ⟨ Z ( β ) ⟩ given by the Laplace transform of an average spectral density ρ ( E ) of the general form in eq. (2.21) . Then:  \nS a ( β ) ≥ 0 ⇐⇒ ρ ( E ) ⊃ N 0 δ ( E ) with N 0 ≥ 1 .  \nIf instead N 0 < 1 , then S a ( β ) < 0 , ∀ β > β ∗ for some β ∗ > 0 , and lim β →∞ S ( β ) = -∞ for N 0 = 0 . 78  \nProof. Applying lemma 1 to ⟨ Z ⟩ = ˆ Z + Z g from eq. (2.22) and recalling E k > 0 gives  \nlim β →∞ ˆ Z ( β ) = N 0 , lim β →∞ Z g ( β ) = 0 , (B.8)  \nand therefore  \nlim β →∞ Z ( β ) = N 0 . (B.9)  \nSuppose first that N 0 > 0. Then starting from eq. (1.8) one easily obtains  \nlim β →∞ S a ( β ) = log N 0 -1 N 0 lim β →∞ βZ ' g ( β ) . (B.10)  \nNote that Z ' g ( β ) is equal to the Laplace transform of -g ( E ) E . Local integrability of g ( E ) including E = 0 requires lim E → 0 + g ( E ) E = 0. Hence by lemma 3, the Laplace transform of -g ( E ) E satisfies lim β →∞ βZ ' g ( β ) = 0. As a result, for N 0 > 0,  \nlim β →∞ S a ( β ) = log N 0 . (B.11)  \nThis already shows that, for N 0 > 0 arbitrarily small, S a ( β ) can be made arbitrarily negative as β → ∞ . It thus seems reasonable to expect that S a ( β ) will diverge negatively as β → ∞ for N 0 = 0, but taking the limit N 0 → 0 above would be unjustified.  \nSuppose now that N 0 = 0, and rewrite eq. (1.8) as  \nS a ( β ) = log Z ( β ) (1 -f ( β )) , f ( e x ) = ∂ x log( -log Z ( e x )) . (B.12)  \nThen, if we define the limit  \nλ ≡ lim β →∞ f ( β ) , (B.13)  \na sufficient condition for S ( β ) to diverge negatively as β →∞ is λ < 1. We now show that indeed only λ < 1 is possible. Using lemma 2, the limiting behavior of f in eq. (B.13) would imply there exist finite c ± , β ∗ > 0 such that  \ne -c -β λ -ϵ > Z g ( β ) > e -c + β λ + ϵ , (B.14)  \nfor any β > β ∗ and arbitrarily small ϵ > 0. For a function g , lemma 1 implies  \nlim β →∞ Z g ( βe iα ) = 0 , (B.15)  \nfor any α ∈ ( -π 2 , π 2 ) , which on the asymptotic form for Z above this imposes λ ≤ 1. 79 Hence it only remains to show that λ = 1 cannot happen.  \nBecause N 0 = 0, the distribution ∆ in eq. (2.21) has no support near E = 0, so the canonical shift that defines ρ in eq. (2.20) implies inf supp g = 0. Recall lim E → 0 + g ( E ) ≡ g (0 + ) ∈ [0 , ∞ ], with divergence allowed so long as g remains integrable. Suppose g (0 + ) > 0. By piecewise continuity, there is a sufficiently small ϵ > 0 such that (0 , ϵ ) ⊂ supp g and g 0 ≡ min g (0 , ϵ ) > 0. Hence,  \nZ g ( β ) ≥ ˆ ϵ 0 + dE g ( E ) e -βE ≥ g 0 β ( 1 -e -βϵ ) . (B.16)  \nThis lower bound on Z g ( β ) is only consistent with the upper bound in eq. (B.14) for λ = 0. Finally, for the sake of contradiction, suppose λ = 1. By lemma 2, for this to hold in the β → ∞ limit, it must be the case that lim β →∞ e cβ Z g ( β ) is finite for some c > 0. However, this would require g to have no support on (0 , c ), for otherwise the limit would diverge. This would contradict the construction, which fixes inf supp g = 0.  \nAltogether, we have shown that lim β →∞ S a ( β ) = log N 0 < 0 for N 0 > 0, and that for N 0 = 0 the limit diverges negatively. If N 0 ≥ 0, then S a ( β ) ≥ 0 for all β by continuity and monotonicity of the Laplace transform. That for N 0 < 1 including N 0 = 0 there exists some finite β ∗ > 0 such that S a ( β ) < 0 for all β > β ∗ follows for the same reason.  \nFor the N 0 = 0 case, note that one could consider writing 80  \nS a ( β ) = log ⟨ Z ( β ) ⟩ -β ⟨ Z ' ( β ) ⟩ ⟨ Z ( β ) ⟩ . (B.17)  \nAs β → ∞ , the first term clearly diverges negatively. However, Z ' ( β ) < 0 for any finite β , so in principle the second term could cancel the divergence of the first. For this second term we know lim β →∞ Z ( β ) = 0 and also lim β →∞ βZ ' ( β ) = 0. In particular, the latter follows from ⟨ Z ⟩ = ˆ Z + Z g by noting that all terms in ˆ Z ' ( β ) decrease exponentially in β , and using lemma 3 for Z ' g ( β ) as the Laplace transform of -Eg ( E ). Then one would hope to apply L'Hˆopital, but this turns out to be useless: both numerator and denominator continue to give zeroes at arbitrary derivative order. Not only is the behavior of this term in the limit not obvious, but it can also be inconclusive. Namely, while in many cases one can show this ratio is finite in the limit, there are examples for which it is divergent and could indeed compete with the first term. These correspond to cases in which the parameter in eq. (B.13) satisfies 0 < λ < 1, which occurs when ρ and all of its derivatives vanish at zero (cf. 1 /ρ ( E ) having an essential singularity at E = 0). For instance, ρ ( E ) ∼ e -1 /E s gives λ = 1 / (1 + s -1 ) for s ≥ 0; functions with non-trivial derivatives at zero give λ = 0.", 'C Supersymmetric Random Matrix Theory': 'This section discusses the relevant details for the construction of supersymmetric Hamiltonians that justifies our general use of eq. (3.13) throughout this paper. Since the implementation of higher supersymmetry in random matrices becomes increasingly intricate, let us begin with the smallest amount of supersymmetry.', 'C.1 N = 1': 'Random matrix theory with N = 1 supersymmetry is discussed in e.g. [20, 186-189]. Its simple structure involves a single self-adjoint supercharge Q and a Hamiltonian given by H = Q 2 . The canonical construction of Q out of the random matrix M takes the form 81  \nQ = ( 0 M M † 0 ) , (C.1)  \nin terms of which the N = 1 random Hamiltonian H = Q 2 is identically given by eq. (3.13). The symmetry group G ( N + ¯ ν ) × G ( N ) naturally realizes the structure of a Z 2 -graded algebra, where  \n˜ Q = ( 0 M 0 0 ) , ˜ Q † = ( 0 0 M † 0 ) ,  \nwhich are nilpotent, ˜ Q 2 = ˜ Q † 2 = 0, and give the Hamiltonian in the anti-commutator form H = 1 2 { ˜ Q, ˜ Q † } .  \nthe sectors transforming under G ( N +¯ ν ) and G ( N ) may respectively be identified as bosonic and fermionic. The grading is implemented in terms of the usual fermion number operator F by  \n( -1) F ≡ ( 1 N +¯ ν 0 0 -1 ¯ ν ) , (C.2)  \nwhich clearly satisfies the desired relations [ H, ( -1) F ] = 0 and { Q, ( -1) F } = 0. Because there are N +¯ ν bosonic and N fermionic states, the supersymmetric index is ¯ ν . In addition, because ( -1) F commutes with H , we have the Witten index  \nTr ( -1) F e -βH = ¯ ν , (C.3)  \nat any inverse temperature β . Hence ¯ ν = 0 is associated to the breaking of all supersymmetry, whereas ¯ ν > 0 corresponds to having ¯ ν states that preserve supersymmetry. That these BPS states precisely capture the degenerate ground state of H is easily seen as follows.  \nIf M is a square AZ matrix of size N , then H is a square matrix of size 2 N whose only N distinct eigenvalues are those of M squared. A global 2-fold degeneracy is in this case expected given that M † M and MM † have identical eigenvalues. Besides these, the resulting Hamiltonian H will generically exhibit no other degeneracies and will, in particular, be non-singular. If instead M is an ( N + ¯ ν ) × N rectangular BPS matrix, then H still is a square matrix, but now of size 2 N +¯ ν . The Hamiltonian H is then singular with at least ¯ ν > 0 zero eigenvalues. This means that generically every H drawn from this ensemble deterministically has ¯ ν degenerate ground states all with the same zero energy, corresponding to the BPS states mentioned above.', 'C.2 N = 2': 'Although N = 1 random matrices can clearly accommodate BPS states, recall that the N = 1 super-Schwarzian spectrum from eq. (2.25) has none. It is thus of interest to also understand the N = 2 matrix construction for the supermultiplet spectra of the N = 2 super-Schwarzian theories in eq. (2.27), which do exhibit BPS states. The presentation here will be concise, and we refer the reader to [43] for more details on the relevant structure of N = 2 supersymmetry.  \nThe N = 2 superalgebra has a Hermitian conjugate pair of nilpotent generators, the supercharges Q and Q † , and a Hamiltonian H = { Q,Q † } they commute with (cf. footnote 81 for N = 1). The algebra also has an R-symmetry given by a U (1) outer automorphism group. There is no unique N = 2 theory, but a two-parameter family thereof. A general N = 2 theory can be specified by an odd integer ˆ q > 0 and, to account for a possible anomaly, by a constant δ ∈ [0 , 1). The R-charges of Q and Q † are respectively ± ˆ q + δ , 82 and states have R-charge k ∈ Z + δ . The Hilbert space is a direct sum H = ⊕ k H k , where H k consists of all states ψ k with R-charge k . Correspondingly, the supercharge Q = ∑ k Q k with the restriction Q k : H k → H k +ˆ q , and similarly for Q † with Q † k : H k +ˆ q → H k . The supersymmetry condition Q 2 = Q † 2 = 0 then implies Q k +ˆ q Q k = Q † k Q † k +ˆ q = 0 for all k .  \nThe algebra decomposes into two types of irreducible supermultiplets. Any state ψ k satisfying Qψ k = Q † ψ k = 0 forms a singlet by itself; these are BPS states. Any pair of states ( ψ k , ψ k +ˆ q )  \n̸  \nobeying Qψ k = λψ k +ˆ q and Q † ψ k +ˆ q = λ ∗ ψ k with λ, λ ∗ = 0 forms a doublet; these are non-BPS states. All BPS states have zero energy since Hψ k = 0, whereas non-BPS states come in degenerate pairs with energy Hψ k = Hψ k +ˆ q = λ ∗ λ > 0. Because doublets always involve consecutive R-charge values ( k, k + ˆ q ), a convenient parameter turns out to be their average R-charge q k ≡ k + ˆ q 2 . In terms of q k , a ( k, k + ˆ q ) doublet involves R-charges q k ∓ ˆ q 2 . All irreducible supermultiplets contain disjoint sets of states, so we can decompose H k into  \nH k = H -k ⊕H 0 k ⊕H + k , (C.4)  \nwhere H 0 k consists of all BPS states in H k , H -k of all non-BPS states in H k forming ( k -ˆ q, k ) doublets with states in H k -ˆ q , and H + k of all non-BPS states in H k forming ( k, k +ˆ q ) doublets with states in H k +ˆ q . By definition, a state is in H + k if and only if it has a doublet pair in H -k +ˆ q , so clearly N ± k ≡ dim H ± k satisfy N + k = N -k +ˆ q . This leads to the following reorganization of the Hilbert space:  \nH = ⊕ k H 0 k ⊕H ↔ q k , H ↔ q k ≡ H + k ⊕H -k +ˆ q . (C.5)  \nHere we have combined all non-BPS states forming ( k, k +ˆ q ) doublets into H ↔ q k , which itself forms a (reducible) supermultiplet that we will refer to as a q k -multiplet. Since H ↔ q k contains an even number of states, we let 2 N ↔ q ≡ dim H ↔ q and note that N ↔ q k = N + k = N -k +ˆ q . Dropping the k from q k we may sometimes use q as a standalone parameter. The definition of q k makes q ∈ Z + δ -1 2 , since ˆ q is always odd. Each of the Hilbert subspaces H 0 k and H ↔ q defines a supermultiplet which turns out to be statistically independent of the rest. In particular, there is a fixed number of zeroenergy BPS states in each H 0 k , and there is an independent probability distribution governing the spectrum of energies of non-BPS states in each H ↔ q . Interestingly though, the spectral statistics of a q -multiplet do happen to be sensitive to the total number of BPS and certain non-BPS states in H 0 q ∓ ˆ q/ 2 . 83 Because these supermultiplets have statistically independent spectra, they are the right target for a matrix ensemble description. Indeed, what eq. (2.27) describes is the spectral density of a q -multiplet, including a potential contribution from BPS states in H 0 q ± ˆ q/ 2 . 84  \nTo understand these spectra consider first expanding the Hamiltonian as  \nH = ∑ k Q k Q † k + Q † k Q k . (C.6)  \nWe want to isolate the operators that act on and thus determine the energies of each q k -multiplet. Since states in H q k only have R-charge k or k + ˆ q , let ψ k and ψ k +ˆ q respectively denote any two of them. Those ψ k ∈ H k can be acted upon by Q k and Q † k -ˆ q , whereas those ψ k +ˆ q ∈ H k +ˆ q can be acted upon by Q k +ˆ q and Q † k . The BPS states are annihilated by any supercharge operator that acts on them, whereas the non-BPS states all come in doublets satisfying Q k ψ k = λψ k +ˆ q and Q † k ψ k +ˆ q = λ ∗ ψ k . But by the supersymmetry conditions Q k +ˆ q Q k = Q † k -ˆ q Q † k = 0 the latter are also annihilated as Q k +ˆ q ψ k +ˆ q = Q † k -ˆ q ψ k = 0. Hence in fact the only supercharge operators that act  \nnon-trivially on H q k are Q k and Q † k . This means one can single out from eq. (C.6) the term  \nH q k ≡ Q k Q † k + Q † k Q k = ( ( QQ † ) q +ˆ q/ 2 0 0 ( Q † Q ) q -ˆ q/ 2 ) . (C.7)  \nas it solely determines the spectrum of the q k -multiplet. In the second expression we have simply made use of the fact that H q acts just on the direct sum H q -ˆ q/ 2 ⊕ H q +ˆ q/ 2 to make eq. (C.7) suggestively similar to eq. (3.13).  \nSince H q k only depends on Q k and its Hermitian conjugate, one has to understand the structure of Q k operators. Recall each of these is a linear map Q k : H k →H k +ˆ q , so letting N k ≡ dim H k one can represent Q k as an N k +ˆ q × N k matrix. Comparing eq. (C.7) to the N = 1 matrix Hamiltonian from eq. (3.13), one would like to identify Q k with an N k +ˆ q × N k random matrix M in an AZ or BPS ensemble. Unfortunately this is not immediately possible given the supersymmetry constraints Q k Q k -ˆ q = 0 that these must satisfy. Nonetheless, by analyzing the interplay between the Jacobian measure of a matrix integral over Q k and a functional constraint imposing Q k Q k -ˆ q = 0, [43] showed that Q k continues to be describable by an independent AZ or BPS random matrix on the orthogonal complement of the image of Q k -ˆ q . However, their focus on short sequences of R-charge values suggested that the effective random sub-matrix of Q k would only be rectangular when H k ⊕H k +ˆ q contained BPS states, and would be a square matrix otherwise. As we discuss in section 3.6.2, this would explain the formation of spectral gaps on q -multiplets when H q ± ˆ q/ 2 contains BPS states, but would not explain why the onset of the spectrum keeps shifting to higher energies when there no longer are any BPS states nearby. In what follows we depart from the analysis of [43] to address this issue.  \nThe strategy will be to analyze an ordered sequence of R-charge values, providing a random matrix representation for each Q k that consistently implements the supersymmetry conditions Q 2 = 0. As we will see, the random matrix description that arises for the spectrum of an arbitrary q -multiplet becomes manifestly intrinsic to the structure of the Hilbert spaces H q ± ˆ q/ 2 relevant to the supermultiplet, irrespective of the order followed in the construction.  \n̸  \nSuppose a sequence of R-charge values begins at k , i.e., H k = ∅ but Q k -ˆ q = 0. In general H k may contain BPS states in H 0 k and non-BPS states in H + k , but clearly H -k = ∅ , so N -k = 0. The first non-trivial supercharge operator Q k : H k →H k +ˆ q may be represented by an N k +ˆ q × N k matrix. Since Q k -ˆ q is trivial, there are actually no constraints on Q k from Q k Q k -ˆ q = 0. Hence Q k should be identified with a random AZ or BPS matrix of size N k +ˆ q × N k . Interestingly, the cases N k +ˆ q > N k and N k +ˆ q ≤ N k turn out to be qualitatively distinct.  \nSuppose first that N k +ˆ q ≤ N k . Then generically rank Q k = N k +ˆ q and ker Q k = N k -N k +ˆ q . Because Q † k -ˆ q = 0, any state in H k annihilated by Q k is automatically a BPS state, giving  \nN 0 k = N k -N k +ˆ q , (C.8)  \nwhich could also be just zero. Since N -k = 0, the remaining states in H k must be forming doublets with states in H k +ˆ q , so N + k = N k +ˆ q . But by the doublet relation between H + k and H -k +ˆ q we know N + k = N -k +ˆ q = N ↔ q k , which in particular implies N k +ˆ q = N ↔ q k . In other words, H k +ˆ q = H ↔ q k solely consists of non-BPS states in ( k, k + ˆ q ) doublets. As a result, H k +ˆ q contains neither BPS states  \nnor non-BPS states that could possibly continue the R-charge sequence beyond k + ˆ q . Having terminated with H = H k ⊕H k +ˆ q , we simply have supercharge Q = Q k and Hamiltonian H = H q k as in eq. (C.7). This simple case reproduces the N = 1 setting with a random AZ or BPS matrix M of dimensions ( N + ¯ ν ) × N having N = N ↔ q k , ¯ ν = N 0 k , and random supercharge Q = M † .  \nSuppose now that N k +ˆ q > N k . Then generically rank Q k = N k and the kernel of Q k in H k is trivial. Hence there are no BPS states in H k and since both N -k = N 0 k = 0 we have N k = N ↔ q k , i.e., every state in H k is in a non-BPS doublet with some state in H k +ˆ q . While the content of H k +ˆ q remains to be determined, let us emphasize that at this point Q k already is rectangular even though there are no BPS states in H k . The random ( N + ¯ ν ) × N matrix M here has N = N ↔ q k , ¯ ν = N k +ˆ q -N ↔ q k , and is identified with Q k = M . The corresponding q k -multiplet spectrum is described by the Hamiltonian H q k given by the right-most expression in eq. (C.7) with Q = M . Crucially, even though the matrix MM † is singular with ¯ ν degenerate zero eigenvalues, these now have nothing to do with BPS states. In fact, this ¯ ν -dimensional kernel of ( MM † ) k +ˆ q simply corresponds to all those states in H k +ˆ q which are not in H ↔ q k . Namely, this term corresponds to Q k Q † k , which annihilates any state in H k +ˆ q that is not in a ( k, k + ˆ q ) doublet, whether or not it is BPS. We thus learn that the statistics of this purely non-BPS q k -multiplet do actually require a rectangular random matrix with ¯ ν > 0 from a BPS ensemble for their description. In addition, the spectrum of non-BPS states will behave just like that of a supermultiplet with ¯ ν BPS states with zero energy, even though in this case there actually are no such states at zero.  \nContinuing the sequence to H k +ˆ q , we know dim H -k +ˆ q = N ↔ q k , and that the remaining N k +ˆ q -N k states in H k +ˆ q must be either BPS or non-BPS in doublets with states in H k +2ˆ q . To address these, we must now study Q k +ˆ q : H k +ˆ q →H k +2ˆ q . Although Q k +ˆ q can be represented by an N k +2ˆ q × N k +ˆ q matrix as usual, we now have to take into account the supersymmetry constraint Q k +ˆ q Q k = 0. Because rank Q k = N ↔ q k , this projects out a subspace of codimension N ↔ q k in H k +ˆ q from the domain of Q k +ˆ q . In other words, Q k +ˆ q effectively becomes representable by a sub-matrix of dimensions N k +2ˆ q × ( N k +ˆ q -N ↔ q k ). This reduced matrix is actually now unconstrained and thus identifiable with a random AZ or BPS matrix. The situation here essentially reduces to our starting point at the beginning of the sequence on H k . More explicitly, note that the subspace of H k +ˆ q that Q k +ˆ q Q k = 0 removes corresponds precisely to the non-BPS states in the q k -multiplet H ↔ q k . Hence the reduced random matrix for Q k +ˆ q is representing a map acting only on H 0 k +ˆ q ⊕H + k +ˆ q ⊂ H k +ˆ q . This is just like the situation we encountered when starting the sequence at R-charge k , where obviously H -k = ∅ and we had H k = H 0 k ⊕ H + k . Hence at this point one can simply iterate the logic we followed thus far, but starting from H k +ˆ q and an operator Q k +ˆ q restricted to an effective domain of codimension N ↔ q k .  \nAbove we chose to analyze a sequence of R-charges starting from some smallest value k , iteratively imposing supersymmetry constraints of the form Q k +ˆ q Q k = 0. However, this order was arbitrary and we could have equally started from some largest value k +ˆ q , and iteratively imposed supersymmetry constraints of the form Q † k -ˆ q Q † k = 0. In fact, the pattern that arises for an arbitrary q k -multiplet is clear, and can be easily described solely in terms of the number of BPS states and structure of non-BPS doublets in H k ⊕H k +ˆ q . Explicitly, the spectrum of the q k -multiplet is captured by Q k : H k → H k +ˆ q , which before imposing any supersymmetry conditions Q 2 = 0 is represented by an N k +ˆ q × N k matrix with N k = N -k + N 0 k + N + k and N k +ˆ q = N -k +ˆ q + N 0 k +ˆ q + N + k +ˆ q ,  \nwhere N + k = N -k +ˆ q ≡ N ↔ q k is the number of non-BPS forming ( k, k + ˆ q ) doublets, and at most one of N 0 k or N 0 k +ˆ q is non-zero and counts the number of BPS states of the corresponding R-charge. The supersymmetry conditions Q k +ˆ q Q k = Q k Q k -ˆ q = 0 restrict the rank of Q k to N ↔ q k , and thus its resulting AZ or BPS random matrix representation has N = N ↔ q k . As for the ¯ ν parameter, if there are BPS states in H k ⊕H k +ˆ q then ¯ ν = max { N 0 k , N 0 k +ˆ q } , whereas if there are no BPS states N 0 k = N k +ˆ q = 0 and then ¯ ν = min { N -k , N + k +ˆ q } = min { N ↔ q k -ˆ q , N ↔ q k +ˆ q } . 85 In the BPS case, if N 0 k > 0 the R-charge sequence terminates at the smallest R-charge k , and if N 0 k +ˆ q > 0 it terminates at the largest R-charge k + ˆ q . In the non-BPS case the sequence need not terminate, unless it turns out N k -ˆ q = 0 or N k +2ˆ q = 0. In either case, we see that the q k -multiplet spectrum is indeed described by eq. (3.13) with the random matrix M drawn from a BPS ensemble with ¯ ν > 0, or possibly an AZ ensemble if ¯ ν = 0.  \nLet us note in passing that if the number of states in q -multiplets 2 N ↔ q e.g. monotonically increases with the R-charge, then we generally find that ¯ ν q = N ↔ q -ˆ q for the q -multiplet. Hence if N ↔ q grows say quadratically in q , we would similarly expect a monotonic growth of ¯ ν q ∝ q 2 . Even if there are no BPS states involved, this intuitively means that the spectrum of the non-BPS states in the q -multiplet with R-charge q -ˆ q/ 2 feels the absence of those states with the same R-charge which actually fell into the ( q -ˆ q )-multiplet.', 'C.3 N = 4': 'For N = 4 supersymmetry, the algebra contains four Hermitian supercharges { Q I } 4 I =1 which transform as vectors of SO (4) and give the Hamiltonian as  \n{ Q I , Q J } = 2 Hδ IJ . (C.9)  \nIn this case the R-symmetry is SU (2) and the Hilbert space can be organized into supermultiplets labelled by a half-integer spin J , which is what the spectral density in eq. (2.29) describes. Once again, multiplets at different values of J are expected to be statistically independent and describable by separate ensembles. Although the specific random matrix construction would now be technically more involved (see [43]), the conclusion seems to be similar to that for the N = 1 , 2 cases: the spectrum of different supermultiplets would be captured by eq. (3.13) with M a random matrix. In this case, however, M should always be drawn from a BPS ensemble where we expect ¯ ν > 0 to be related to the varying number of non-BPS states associated to supermultiplets of different spin.', 'D Random Matrix Tools': 'This section reviews some standard tools for the analysis of random matrix integrals.', 'D.1 The Resolvent': "The resolvent formalism is particularly suited to the study of spectral properties of operators, including but not limited to matrices. Given an operator U , the resolvent operator is  \nR ( z ) ≡ ( z 1 -U ) -1 . (D.1)  \nLet ρ denote the spectral density of U . The support supp ρ defines the spectrum of U . The resolvent R is only well-defined on C ∖ supp ρ , an open set called the resolvent set. Taking the trace of R leads to 86  \nR ( z ) ≡ Tr R ( z ) = ˆ R du ρ ( u ) z -u , z ∈ C ∖ supp ρ. (D.2)  \nHenceforth, R will be referred to simply as the resolvent. Assume the spectrum of U is real, and define ρ on all R by letting ρ ( u ) = 0 for any u ∈ R ∖ supp ρ . The integral representation in eq. (D.2) may be recognized as the Stieltjes transform of a non-negative measure of density ρ . 87 Importantly, as defined away from supp ρ , note that R is an analytic function. For a general point z in the resolvent set, eq. (D.2) can be evaluated as a contour integral closing around the pole at z . 88 Near the real axis, letting z = x + iϵ , separating real and imaginary parts of R , and taking the limit ϵ → 0, one obtains a representation of the resolvent in terms of the Sokhotski-Plemelj theorem 89  \nlim ϵ → 0 + R ( x ± iϵ ) = R du ρ ( u ) x -u ∓ iπρ ( x ) , (D.3)  \nwhere ffl denotes a principal value integral. The inverse to the Stieltjes transform immediately follows from this expression. Explicitly, the spectral density ρ can be extracted from the imaginary part of eq. (D.3) or, equivalently, from the discontinuity of R across the real axis,  \nρ ( x ) = ± 1 π lim ϵ → 0 + Im R ( x ∓ iϵ ) = i 2 π lim ϵ → 0 + ( R ( x + iϵ ) -R ( x -iϵ )) . (D.4)  \nIt is in this way that the resolvent of an operator precisely captures its spectral density.  \nIn the context of random matrix theory, the operators of interest are functions of the random matrix. Oftentimes though, one is just interested in the spectral properties of the matrices in the ensemble, and thus studies the resolvent for the random matrix itself. However, more generally, R  \nR ( u ) = ∑ i 1 u -u i .  \nNote that if the spectrum is a finite set, ρ integrates to the total number of eigenvalues.  \n87 In solving matrix integrals, oftentimes resolvent methods involving R are not needed, and it suffices to make use of the analyticity properties of the Stieltjes transform of the spectral density (see e.g. the derivation of the loop equations [20, 128, 190]). Nonetheless, following conventions, we will continue to refer to R as the resolvent.  \n89 Even when the definition of R is only formal (see footnote 88), its imaginary part can still be evaluated if the limit ϵ → 0 is taken before sending the contour to infinity. The real part of the resolvent may diverge either way. Its imaginary part diverges for any finite ϵ , but is well-defined if the ϵ → 0 limit is taken first.  \nstands for the resolvent of any matrix operator U , and ρ is the spectral density of that operator. Since U may be any non-trivial functional of the random matrix, ρ may generally be very different from the spectral density of the random matrix. For example, the construction of Hamiltonians for AZ ensembles in section 3.3 makes the matrix operator of interest quadratic in the actual random matrices, and the leading spectrum of these turns out to exhibit a hard edge which would not be present otherwise (see section 3.6.1 and appendices D.4 and D.5).  \nGiven that the resolvent of an operator encodes its spectral density and can be inverted to obtain it, both can be used to extract statistical properties of the ensemble. In particular, correlation functions of spectral densities and resolvents are related by  \n⟨ R ( z 1 ) · · · R ( z n ) ⟩ = ˆ R [ n ∏ i =1 du i z i -u i ] ⟨ ρ ( u 1 ) · · · ρ ( u n ) ⟩ , ∀ z i ∈ C ∖ supp ρ. (D.5)  \nThese resolvent correlation functions are particularly useful for their role in the loop equations that allow for solving matrix integrals perturbatively in a 1 /N expansion [20]. 90 While we will not be using this machinery here, it will still be very important for us to study leading order results at large N . Expanding both sides of eq. (D.5) and taking the largeN limit for n = 1 gives  \nR ∗ ( z ) = ˆ R du ˆ ρ ∗ ( u ) z -u , z ∈ C ∖ supp ˆ ρ ∗ , (D.6)  \nwhere eq. (3.30) is the leading average spectral density, unit-normalized by absorbing factors of N on both sides. These largeN objects correspond precisely to the continuum saddle-point results studied in section 3.5. In particular, by eq. (3.31), the spectral density ˆ ρ ∗ for the Hamiltonians of interest for WD (energies linear in x ) and AZ/BPS (energies quadratic in x ) ensembles obeys the equilibrium condition  \nR dy ˆ ρ ∗ ( y ) x -y -V ' ν ( x ) β = 0 , ∀ x ∈ supp ˆ ρ ∗ . (D.7)  \nUsing this, for the leading spectral density the Sokhotski-Plemelj representation of the resolvent in eq. (D.3) can be written  \nlim ϵ → 0 + R ∗ ( x ± iϵ ) = V ' ν ( x ) β ∓ iπ ˆ ρ ∗ ( x ) , x ∈ supp ˆ ρ ∗ . (D.8)  \nIn turn, this provides a discontinuity relation for the resolvent in terms of the matrix potential,  \nV ' ν ( x ) = β 2 lim ϵ → 0 + ( R ∗ ( x + iϵ ) + R ∗ ( x -iϵ )) x ∈ supp ˆ ρ ∗ , (D.9)  \nThese important relations allow one to formulate the study of the spectrum of matrix integrals as a Riemann-Hilbert problem, which we turn to in appendix D.3.", 'D.2 The Spectral Curve': "The analytic continuation of the left-hand side of eq. (D.7) off the support of the spectral density ˆ ρ ∗ defines a particularly important object known as the spectral curve,  \ny ( z ) ≡ R ∗ ( z ) -V ' ν ( z ) β , z ∈ C ∖ supp ˆ ρ ∗ . (D.10)  \nRecalling that eq. (D.7) comes from eq. (3.29), y is related to the leading effective potential by  \n̂ V ' ∗ ( z ) = -β y ( z ) , z ∈ C ∖ supp ˆ ρ ∗ . (D.11)  \nIntegrating gives the leading effective potential off the support of ˆ ρ ∗ as  \n̂ V ∗ ( z ) = -β ˆ z dwy ( w ) , z ∈ C ∖ supp ˆ ρ ∗ , (D.12)  \nup to a constant which may be fixed by making ̂ V ∗ ( z ) vanish as z → supp ˆ ρ ∗ (cf. eq. (3.29)).  \nThe spectral curve is a useful object due to its analyticity properties. Since the matrix potential V is typically analytic, it does not contribute any poles to y . In WD ensembles ν = 0 so V ν is also analytic, and y just inherits the branch cut of R ∗ along the real support of ˆ ρ ∗ . For the AZ/BPS ensembles we are actually interested in the spectrum of squared eigenvalues, which because the integral measure dx = d ( x 2 ) 2 x induces a 1 /x pole on y ( x ) at the origin (cf. eq. (D.57)). In addition, if ν = νN ∼ O ( N ) as in the BPS ensembles, a logarithmic divergence ν log x survives the largeN limit in V ν , giving a total 1 /x 2 pole to y ( x ). This does not happen for AZ ensembles since ν ∼ O (1) and this pole is ignored at infinite N .  \nGenerally V ν does not have any branch cuts, so applying eq. (D.8) to eq. (D.10) we have  \nlim ϵ → 0 ± y ( x + iϵ ) = ∓ iπ ˆ ρ ∗ ( x ) , x ∈ supp ˆ ρ ∗ , (D.13)  \ni.e., the discontinuity of y across the support of ˆ ρ ∗ precisely captures ˆ ρ ∗ (cf. eq. (D.4)). Because eq. (D.13) corresponds to a square-root branch cut, 91 it follows from the above that for WD ensembles y ( x ) 2 defines an analytic function everywhere. For AZ ensembles we can still obtain an analytic function by considering the combination y ( x ) 2 x . As for BPS ensembles, there actually remains a 1 /x pole at x = 0 for y ( x ) 2 x . One could obviously consider instead y ( x ) 2 x 2 , but this turns out to not be necessary. The logarithmic divergence of V ν as x → 0 + means the support of ˆ ρ ∗ cannot possibly extend all the way to x = 0. Since this BPS pole is off the the support of ˆ ρ ∗ , for BPS ensembles we will still be interested in the analyticity properties of y ( x ) 2 x .  \nIn general, we see that y ( x ) 2 is an analytic function except for possible poles at x = 0, and thus the spectral curve itself, y ( z ), defines a double cover of the complex plane branched over the support of ˆ ρ ∗ . An important corollary of this is that the value y attains off the support of ˆ ρ ∗ is given by the analytic continuation of ˆ ρ ∗ itself (cf. eq. (3.35)). Correspondingly, by eq. (D.12), the effective potential for ˆ ρ ∗ is given by an integral of the analytic continuation of ˆ ρ ∗ itself.  \nLet us end this section by quoting the spectral curve for each of our ensembles of interest with a Gaussian potential. For WD ensembles, eq. (3.37) gives  \ny ( z ) GWD = 1 β √ z 2 -2 β , z ∈ C ∖ [ -√ 2 β , √ 2 β ] , (D.14)  \nwhich indeed is analytic everywhere it is defined. For AZ ensembles, eq. (3.40) gives  \ny ( z ) GAZ = 1 2 β √ z -4 β z , z ∈ C ∖ [0 , 4 β ] . (D.15)  \nSince z is now referring to squared eigenvalues, this reproduces the expected 1 /x pole at x = 0. Finally, for BPS ensembles, eq. (3.44) gives  \ny ( z ) GBPS = √ ( a + -z )( a --z ) 2 β z , z ∈ C ∖ [ a -, a + ] , (D.16)  \nwith the cut endpoints given by eq. (3.45). Once again, since z refers to squared eigenvalues, this objects exhibits the expected 1 /x 2 pole at x = 0. However, as pointed out above, this pole indeed lies off the [ a -, a + ] support of ˆ ρ ∗ .", 'D.3 General Equilibrium Measure': "This section shows how to solve eq. (3.31) to obtain an explicit expression for the equilibrium measure ˆ ρ ∗ that minimizes the action ̂ I in eq. (3.25) for an arbitrary matrix potential V .  \nAccording to [191, 192], eq. (D.7) is a singular integral equation which has a unique solution under fairly general assumptions which have been implicit throughout. 92 Furthermore, the unique equilibrium measure ˆ ρ ∗ that attains the infimum of ̂ I can be shown to be absolutely continuous and supported on a finite union of intervals,  \nΣ ≡ supp ˆ ρ ∗ = q ⋃ j =1 [ a j , b j ] , (D.17)  \nwhere a j , b j ∈ R are ordered such that a j < b j < a j +1 < b j +1 for all j = 1 , . . . q -1 93 . By eq. (D.13), we know that across each of these intervals, the spectral curve y has a sign-change discontinuity, which must take the form of a square-root branch cut. For the WD ensembles, this same branch cut structure can be reproduced by the square root of the following function:  \nη ( z ) WD ≡ q ∏ j =1 ( z -a j )( z -b j ) . (D.18)  \nFor the AZ ensembles, we know y ( x ) 2 additionally has a 1 /x pole at the origin (cf. the discussion below eq. (D.13)), so for these we instead define  \nη ( z ) AZ ≡ 1 z 2 q ∏ j =1 ( z -a j )( z -b j ) . (D.19)  \nIn writing this expression, we are using the fact that these ensembles typically have a 1 = 0. Hence one needs a factor of 1 /z 2 in order for η to develop the right 1 /z pole. For BPS ensembles the  \npotential diverges at x = 0, thus effectively repelling eigenvalues away from the origin towards positive values. This effect forces the endpoint closest to x = 0 to obey a 1 > 0 strictly. As a result, to cancel the non-analyticities of y one can still just use eq. (D.19): a nonzero a 1 here implies the desired 1 /z 2 pole for η . As a result for BPS ensembles eq. (D.19) actually exhibits a 1 /z 2 pole, precisely matching the behavior of y ( z ) as well. Hence eq. (D.19) successfully reproduces the pole structure of the spectral curve of both AZ and BPS ensembles.  \nThese η functions can thus be used to cancel out any non-analyticity of y . Defining the square root of η with its branch cuts along Σ, the combination y ( z ) / √ η ( z ) no longer has branch cuts or poles anywhere, and is in particular analytic everywhere. The upshot is that for a contour C around Σ, the integral  \n1 2 πi ˛ C dw w -z y ( w ) √ η ( w ) = 0 , (D.20)  \nfor any point z ∈ C not enclosed by C . By analyticity of the integrand, the contour C can be deformed to surround Σ arbitrarily tightly, which means eq. (D.20) in fact holds for any z ∈ C ∖ Σ. This result can be applied to eq. (D.10) in order to eliminate the dependence on the spectral curve, and obtain a direct relation between the resolvent and the matrix potential. To do so, consider using Cauchy's integral formula to write  \nR ∗ ( z ) √ η ( z ) = 1 2 πi ˛ γ z dw w -z R ∗ ( w ) √ η ( w ) , z ∈ C ∖ Σ , (D.21)  \nfor γ z a counterclockwise contour around a small neighborhood of z . Since R ∗ ( z ) clearly goes as 1 /z at infinity according to eq. (D.6), the integrand above is guaranteed to go to zero faster than 1 /z at infinity. Together with the analyticity of this expression away from Σ, this means we can deform γ z past infinity without picking up any contributions. Retracting the contour back from infinity we obtain a contour that goes around Σ with clockwise orientation. This can be matched onto the contour C used in eq. (D.20), leading to  \nR ∗ ( z ) √ η ( z ) = -1 2 πi ˛ C dw w -z R ∗ ( w ) √ η ( w ) , z ∈ C ∖ Σ , (D.22)  \nwhere we are sticking to C running counterclockwise. Substituting eq. (D.10) for R ∗ on the righthand side above and using eq. (D.20), we arrive at 94  \nR ∗ ( z ) = -1 β √ η ( z ) h ( z ) , h ( z ) ≡ 1 2 πi ˛ C dw w -z V ' ν ( w ) √ η ( w ) , z ∈ C ∖ Σ . (D.23)  \nUsing eq. (D.10) allows to express the spectral curve explicitly in terms of the matrix potential,  \ny ( z ) = -V ' ν ( z ) β -1 β √ η ( z ) h ( z ) , x ∈ C ∖ Σ . (D.24)  \nBy eq. (D.13), to obtain ˆ ρ ∗ it just remains to understand the limit in which z approaches Σ. The non-trivial aspect is the evaluation of the contour integral for h in eq. (D.23). Firstly, flattening the contour around Σ, one easily arrives at  \nh ( z ) = 1 π ˆ Σ dy y -z V ' ν ( y ) + √ η ( y ) , z ∈ C ∖ Σ , (D.25)  \nwhere it has been convenient to introduce the following branch cut prescription:  \n+ √ η ( x ) ≡ lim ϵ → 0 + Im √ η ( x + iϵ ) . (D.26)  \nTaking z towards Σ, we obtain a Sokhotski-Plemelj representation of h (cf. eq. (D.3),  \nlim ϵ → 0 + h ( x ± iϵ ) = 1 π Σ dy 1 y -x V ' ν ( y ) + √ η ( y ) ± i V ' ν ( x ) + √ η ( x ) , x ∈ Σ . (D.27)  \nEquipped with this result, we can finally obtain an explicit expression for the equilibrium measure. Taking the limit of z towards Σ in eq. (D.24) and using eq. (D.13),  \nˆ ρ ∗ ( x ) = + √ η ( x ) β π 2 Σ dy y -x V ' ν ( y ) + √ η ( y ) , x ∈ Σ , (D.28)  \nwhere recall the support Σ = supp ˆ ρ ∗ takes the form of eq. (D.17). 95 For instance, one can easily verify the expressions for WD, AZ, and BPS ensembles with Gaussian potentials respectively in eqs. (3.37), (3.40) and (3.44) by applying eq. (D.28) to the single-cut case q = 1 using the appropriate function from eqs. (D.18) and (D.19). The determination of the location of the a j , b j endpoints of the support of ˆ ρ ∗ from eq. (D.17), as well as the correct number of cuts q , however, has to be worked out separately.  \nEssentially, what eq. (D.28) provides is a solution to the Euler-Lagrange equation for the matrix integral action, which is guaranteed to give a local minimum. The solution to the extremization problem that is unique is the infimum of this action, which requires finding the appropriate values of q and the a j , b j endpoints for a global minimum [193]. For a general matrix potential V , one may wonder how to even determine the number of intervals q the measure ˆ ρ ∗ will be supported on, and their specific endpoints. In general q is not known, and one has to try different values. 96 Given a choice of q , one attempts to find the endpoints of the intervals as follows. Expanding eq. (D.23) at large | z | , and using the fact that the resolvent R ( z ) goes as 1 /z at large | z | , one obtains the following relations,  \nΣ dx x l V ' ( x ) + √ η ( x ) = β π δ l ˜ q , l = 0 , . . . , ˜ q, (D.29)  \nwhere ˜ q = q for the WD ensembles and ˜ q = q -1 for the AZ/BPS ones. This difference between ensembles comes from the fact that at large | z | eq. (D.18) gives η ( z ) = z q + O ( | z | q -1 ) for the former, whereas eq. (D.19) gives η ( z ) = z q -1 + O ( | z | q -2 ) for the latter. Additionally, for the AZ ensembles we know that a 1 = 0 from non-negativity, and for the BPS ensembles a 1 > 0 implies  \nΣ dy + √ η ( y ) 1 y -x = 0 , x ∈ Σ ,  \nwhich can be used to turn eq. (D.28) into [100]  \nˆ ρ ∗ ( x ) = 1 β π + √ η ( x ) Q ( x ) , Q ( x ) ≡ 1 π Σ dy + √ η ( y ) V ' ( x ) -V ' ( y ) x -y , x ∈ Σ .  \neq. (D.29) also vanishes for l = -1 due to the enhanced pole. Altogether, this means that in all cases we have a system of q +1 equations for the 2 q endpoints of the intervals that make up Σ.  \nAnother q -1 equations are needed to fully fix these. Consider the effective potential from eq. (D.12) along the real axis. By the equations of motion we know ̂ V ∗ must vanish at the a j , b j endpoints of the support intervals. In its indefinite integral expression from eq. (D.11) this means ̂ V ∗ must attain the same constant value on all of them. This gives us the desired additional q -1 equations as  \n0 = ̂ V ∗ ( b j ) -̂ V ∗ ( a j +1 ) = β ˆ a j +1 b j dxy ( x ) , j = 1 , . . . , q -1 . (D.30)  \nAltogether we obtained 2 q conditions for the 2 q variables determining the endpoints of supp ˆ ρ ∗ . These can be shown to be linearly independent (cf. remark 11.2.6 from [100]), thus completing the characterization of the unique equilibrium measure ˆ ρ ∗ for a general matrix integral.", 'D.4 Orthogonal Polynomials': "A standard treatment in random matrix theory proceeds by factoring the matrix potential term into the eigenvalue determinant and then performing elementary row operations so as to build a suitably useful sequence of orthogonal polynomials. The upshot of such manipulations is a repacking of the joint PDF of eigenvalues into determinants,  \np n ( x 1 , . . . , x n ) = 1 N ! det { K N ( x i , x j ) } N i,j =1 . (D.31)  \nHere K N is a symmetric integral kernel given by  \nK N ( x, y ) = N -1 ∑ k =0 φ k ( x ) φ k ( y ) , (D.32)  \nwith { φ k } N -1 k =0 a family of polynomials orthonormalizing the matrix integral,  \nˆ dxφ j ( x ) φ k ( x ) = δ jk . (D.33)  \nThis kernel has the following important properties:  \nˆ duK N ( x, u ) K N ( u, y ) = K N ( x, y ) , K N ( x, x ) = ρ ( x ) . (D.34)  \nThe Christoffel-Darboux formula reduces the series in eq. (D.32) to  \nK N ( x, y ) = c N φ N ( x ) φ N -1 ( y ) -φ N ( y ) φ N -1 ( x ) x -y , (D.35)  \nwhere the N -dependent constant can be extracted from  \nc N ≡ lim x →∞ xφ N -1 ( x ) φ N ( x ) . (D.36)  \nFor the GUE, the matrix integral involves the measure in eq. (3.10) with β = 2 and ν = 0, and the matrix potential is V ( x ) = 1 2 x 2 . The { φ k } polynomials in this case are built so as to orthonormalize the sequence { x k e -1 2 x 2 } over R . This is accomplished by  \nφ k ( x ) = 1 (2 π ) 1 / 4 √ k ! h k ( x ) e -x 2 4 , h k ( x ) ≡ ( -1) k e x 2 2 d k dx k e -x 2 2 , (D.37)  \nwhere h k are standard Hermite polynomials. For these the constant in eq. (D.36) is c N = √ N , and the spectral density can be readily computed using eq. (D.34). The Wigner semicircle follows from the following scaling limit:  \nˆ ρ ∗ ( x ) ≡ lim N →∞ 1 √ N ρ ( √ Nx ) = 1 2 π √ 4 -x 2 , -2 < x < 2 . (D.38)  \nIn this paper, this is actually the canonical scaling that our conventional factor of N in front of the matrix potential in eq. (3.2) implements. 97 The Airy spectral density is obtained by the following scaling limit which zooms on the lower edge of the spectrum [107]:  \nˆ ρ ∞ ( x ) ≡ lim N →∞ 1 N 1 / 6 ρ ( -2 √ N + x N 1 / 6 ) = Ai ' ( -x ) 2 + x Ai( -x ) 2 . (D.39)  \nThe exact spectral density of the Airy model is defined to be  \nˆ ρ Airy ( x ) ≡ e -S 0 / 3 ˆ ρ ∞ ( e 2 S 0 / 3 x ) , (D.40)  \nwhere e S 0 captures the eigenvalue density. Its leading spectral density at large S 0 is  \nˆ ρ Airy ∗ ( x ) ≡ lim S 0 →∞ ˆ ρ Airy ( x ) = √ x π , x > 0 , (D.41)  \nthereby reducing to the behavior of eq. (D.38) near the x = -2 edge of the spectrum. Because ˆ ρ Airy ∗ goes to zero at its edge, the Airy spectrum is said to exhibit a soft edge.  \nA similar exercise can be done for the AZ/BPS ensembles. In particular, consider the measure in eq. (3.10) applied to squared eigenvalues for β = 2 and any ν . With a Gaussian potential the squared eigenvalues experience a linear potential of the form V ( x ) = 1 2 x . For this integral, the { φ k } polynomials are built so as to orthonormalize the sequence { x k x ν e -1 2 x } over R + . This is accomplished by  \nφ k ( x ) = √ k ! 2 ν +1 ( k + ν )! L ν k ( x/ 2) x ν / 2 e -x 4 , L ν k ( x ) ≡ x -ν e x k ! d k dx k ( x k + ν e -x ) , (D.42)  \nwhere L ν k are generalized Laguerre polynomials. For this reason, these ensembles are often referred to as Laguerre unitary ensembles. These polynomials give c N = -√ 2 N ( N + ν ) for the constant in eq. (D.36). A singular Marˇcenko-Pastur law is obtained in the following scaling limit:  \nˆ ρ ∗ ( x ) = lim N →∞ ρ ( Nx ) = 1 4 π √ 8 -x x , 0 < x < 8 . (D.43)  \nThe divergent behavior of this distribution as x → 0 is commonly referred as a hard edge, in contrast with the soft edge of eq. (D.38). Zooming onto this hard edge gives rise to Bessel spectral densities [109],  \nˆ ρ ∞ ν ( x ) ≡ lim N →∞ 1 N ρ ( x 2 N ) = 1 2 ( J ν ( √ x ) 2 -J ν +1 ( √ x ) J ν -1 ( √ x )) . (D.44)  \nThe exact spectral density of the Bessel model is defined to be  \nˆ ρ Bessel ν ( x ) ≡ e S 0 ˆ ρ ∞ ν ( e 2 S 0 x ) , (D.45)  \nwhose behavior at large S 0 yields the leading spectral density  \nˆ ρ Bessel ∗ ( x ) = lim S 0 →∞ ˆ ρ Bessel ν ( x ) = 1 π √ x , x > 0 . (D.46)  \nAn alternative standard scaling limit of the Laguerre ensembles involves making the parameter ν extensive in N by letting ν ≡ νN and keeping ν ∼ O (1). Doing so, the analogous limit to eq. (D.43) leads in this case to a non-singular Marˇcenko-Pastur law,  \nˆ ρ ν ∗ ( x ) = lim N →∞ ρ ν ( Nx ) = √ ( a + -x )( x -a -) 4 πx , a -< x < a + , (D.47)  \nwhere the endpoints of the leading spectrum are given by  \na ± ≡ 2 ( 1 -√ 1 + ν ) 2 . (D.48)  \nThe previous case in eq. (D.43) is recovered in the limit ν → 0 in which a -→ 0 and the pole at x = 0 becomes again accessible. However, for finite ν > 0, the distribution ˆ ρ ν ∗ has now a characteristic gap (0 , a -) with no support. Additionally, x = a -now becomes a soft edge of the same √ x form as in the Airy model.", 'D.5 WD/AZ Ensemble Relation': "In section 3.2 we replaced the natural AZ eigenvalues λ i by their squares ω i ≡ λ 2 i at the level of the measure. This allowed us to obtain a general measure in eq. (3.10) whose form captures all WD and AZ ensembles together. However, we had to keep in mind that for the WD ensembles the eigenvalues λ i could extend over all R , whereas for AZ ensembles the new ω i variables had to be restricted to being non-negative. Additionally, at the level of the potential this meant that e.g. a quadratic term in λ i would just give a linear term in ω i . Here we make this process more transparent by working out the AZ matrix integral in its original eigenvalues λ i which extend over all R and are subject to the original AZ measure in eq. (3.7).  \nFirstly, let us quote that for the WD ensembles the equilibrium condition is simply  \nR dy ˆ ρ ( y ) x -y -V ' ( x ) β WD = 0 , x ∈ supp ˆ ρ ∗ , (D.49)  \nwhere we just set ν = 0 in eqs. (3.27) and (3.29), which followed from extremizing the effective potential in eq. (3.26). On the other hand, using eq. (3.7), the continuum effective potential for the AZ ensemble reads  \n̂ V [ ˆ ρ ; x ] AZ = V α ( x ) -ˆ R dy ˆ ρ ( y ) log | x 2 -y 2 | β , (D.50)  \nwhere V α uses the notation from eq. (3.17), and thus eq. (3.29) gives the equilibrium condition  \nR dy ˆ ρ ∗ ( y ) x 2 -y 2 -V ' α ( x ) 2 x β AZ = 0 , ∀ x ∈ supp ˆ ρ ∗ . (D.51)  \nThis is more closely related to eq. (D.49) than it may seem at first glance. Using the identity  \n1 x 2 -y 2 = 1 2 x ( 1 x -y + 1 x + y ) (D.52)  \none easily rearranges eq. (D.51) into  \nR dy ( ˆ ρ ∗ ( y ) + ˆ ρ ∗ ( -y ) x -y ) -V ' α ( x ) β AZ = 0 , ∀ x ∈ supp ˆ ρ ∗ . (D.53)  \nComparing this to eq. (D.49), the following result is immediate: if ˆ ρ WD ( β ) ∗ is the leading spectral density for a WD ensemble with parameter β and matrix potential V α , then the leading spectral density for the AZ ensemble with the same potential and parameter β must obey 98  \nˆ ρ AZ ( β ) ∗ ( x ) + ˆ ρ AZ ( β ) ∗ ( -x ) = 2ˆ ρ WD (2 β ) ∗ ( x ) . (D.54)  \nFor instance, if α ∼ O (1) in the largeN limit, then V α → V and e.g. for a Gaussian potential V ( x ) = 1 2 x 2 the solution to eq. (D.51) is  \nˆ ρ ∗ ( x ) GAZ = 1 2 π β √ 4 β -x 2 , -2 √ β < x < 2 √ β , (D.55)  \nwhich recalling eq. (3.37) verifies our finding in eq. (D.54). In this sense, we thus see that for general matrix potentials, WD and AZ/BPS ensembles are in fact rather similar at leading order, particularly if α /N → 0 in the largeN limit. What actually makes them quite different in the context of random Hamiltonians, whether or not α is extensive in N , is the fact that for AZ ensembles we are interested in the spectrum of squared eigenvalues (see section 3.3).  \nDenote the leading spectral density of the Hamiltonian H for an AZ ensemble by ˆ µ ∗ , and that of the random AZ matrix by ˆ ρ ∗ , as above. By the construction of H in section 3.3, the eigenvalues of the former χ are related to those of the latter x simply by χ = x 2 . Clearly ˆ µ ∗ ( χ ) = 0 for any χ < 0. To get an expression for µ on the non-negative real axis, we demand that for any even function f , the expectation values computed with respect to either density agree. In particular, this must be true for a constant function,  \nˆ ∞ -∞ dx ˆ ρ ∗ ( x ) = ˆ ∞ 0 d ( x 2 ) 2 x (ˆ ρ ∗ ( x ) + ˆ ρ ∗ ( -x )) , (D.56)  \nfrom which one reads off 99  \nˆ µ ∗ ( x 2 ) = ˆ ρ ∗ ( x ) + ˆ ρ ∗ ( -x ) 2 x , x > 0 . (D.57)  \nHence, on general grounds, the spectral density ˆ µ ∗ ( χ ) of the square of an operator will exhibit a 1 / √ χ singularity near χ = 0 + , and a discontinuity at the origin since ˆ µ ∗ ( χ ) = 0 for χ < 0. For e.g. the Gaussian case in eq. (D.55), the result in eq. (D.57) yields (cf. eq. (3.40))  \nˆ µ ∗ ( χ ) GAZ = 1 2 π β √ 4 β -χ χ , 0 < χ < 4 β . (D.58)  \nThe procedure above shows rather explicitly how the leading spectral densities of WD and AZ/BPS ensembles are directly related in terms of their original eigenvalue variables. This is particularly relevant for all the classical AZ ensembles where α is a fixed constant and therefore the potential V α appearing in eq. (D.51) at large N reduces to the same potential V that would appear in the WDequilibrium condition. Additionally, we see that for general AZ ensembles, the hard edge that arises when α /N → 0 at large N is essentially a consequence of squaring eigenvalues, rather than an intrinsic feature of the AZ ensembles. In other words, the original eigenvalues of an AZ matrix integral exhibit no hard edges, just like a WD one.", 'References': '- [1] S. W. Hawking, Particle Creation by Black Holes , Commun. Math. Phys. 43 (1975) 199.\n- [2] S. Ryu and T. Takayanagi, Holographic derivation of entanglement entropy from AdS/CFT , Phys. Rev. Lett. 96 (2006) 181602 [ hep-th/0603001 ].\n- [3] V. E. Hubeny, M. Rangamani and T. Takayanagi, A Covariant holographic entanglement entropy proposal , JHEP 07 (2007) 062 [ 0705.0016 ].\n- [4] A. Lewkowycz and J. Maldacena, Generalized gravitational entropy , JHEP 08 (2013) 090 [ 1304.4926 ].\n- [5] T. Faulkner, A. Lewkowycz and J. Maldacena, Quantum corrections to holographic entanglement entropy , JHEP 11 (2013) 074 [ 1307.2892 ].\n- [6] N. Engelhardt and A. C. Wall, Quantum Extremal Surfaces: Holographic Entanglement Entropy beyond the Classical Regime , JHEP 01 (2015) 073 [ 1408.3203 ].\n- [7] G. Penington, Entanglement Wedge Reconstruction and the Information Paradox , JHEP 09 (2020) 002 [ 1905.08255 ].\n- [8] A. Almheiri, N. Engelhardt, D. Marolf and H. Maxfield, The entropy of bulk quantum fields and the entanglement wedge of an evaporating black hole , JHEP 12 (2019) 063 [ 1905.08762 ].\n- [9] G. Penington, S. H. Shenker, D. Stanford and Z. Yang, Replica wormholes and the black hole interior , JHEP 03 (2022) 205 [ 1911.11977 ].\n- [10] A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, Replica Wormholes and the Entropy of Hawking Radiation , JHEP 05 (2020) 013 [ 1911.12333 ].\n- [11] N. Engelhardt, S. Fischetti and A. Maloney, Free energy from replica wormholes , Phys. Rev. D 103 (2021) 046021 [ 2007.07444 ].\n- [12] A. Strominger and C. Vafa, Microscopic origin of the Bekenstein-Hawking entropy , Phys. Lett. B 379 (1996) 99 [ hep-th/9601029 ].\n- [13] T. Banks, W. Fischler, S. H. Shenker and L. Susskind, M theory as a matrix model: A conjecture , Phys. Rev. D 55 (1997) 5112 [ hep-th/9610043 ].\n- [14] J. M. Maldacena, The Large N limit of superconformal field theories and supergravity , Adv. Theor. Math. Phys. 2 (1998) 231 [ hep-th/9711200 ].\n- [15] S. S. Gubser, I. R. Klebanov and A. M. Polyakov, Gauge theory correlators from noncritical string theory , Phys. Lett. B 428 (1998) 105 [ hep-th/9802109 ].  \n- [16] E. Witten, Anti-de Sitter space and holography , Adv. Theor. Math. Phys. 2 (1998) 253 [ hep-th/9802150 ].\n- [17] J. R. David, G. Mandal and S. R. Wadia, Microscopic formulation of black holes in string theory , Phys. Rept. 369 (2002) 549 [ hep-th/0203048 ].\n- [18] D. N. Page, Hawking radiation and black hole thermodynamics , New J. Phys. 7 (2005) 203 [ hep-th/0409024 ].\n- [19] P. Saad, S. H. Shenker and D. Stanford, JT gravity as a matrix integral , 1903.11115 .\n- [20] D. Stanford and E. Witten, JT gravity and the ensembles of random matrix theory , Adv. Theor. Math. Phys. 24 (2020) 1475 [ 1907.03363 ].\n- [21] A. Belin and J. de Boer, Random statistics of OPE coefficients and Euclidean wormholes , Class. Quant. Grav. 38 (2021) 164001 [ 2006.05499 ].\n- [22] N. Afkhami-Jeddi, H. Cohn, T. Hartman and A. Tajdini, Free partition functions and an averaged holographic duality , JHEP 01 (2021) 130 [ 2006.04839 ].\n- [23] A. Maloney and E. Witten, Averaging over Narain moduli space , JHEP 10 (2020) 187 [ 2006.04855 ].\n- [24] J. Cotler and K. Jensen, AdS 3 gravity and random CFT , JHEP 04 (2021) 033 [ 2006.08648 ].\n- [25] N. Benjamin, C. A. Keller, H. Ooguri and I. G. Zadeh, Narain to Narnia , Commun. Math. Phys. 390 (2022) 425 [ 2103.15826 ].\n- [26] S. Collier and A. Maloney, Wormholes and spectral statistics in the Narain ensemble , JHEP 03 (2022) 004 [ 2106.12760 ].\n- [27] J. Chandra, S. Collier, T. Hartman and A. Maloney, Semiclassical 3D gravity as an average of large-c CFTs , JHEP 12 (2022) 069 [ 2203.06511 ].\n- [28] S. Collier, L. Eberhardt and M. Zhang, Solving 3d gravity with Virasoro TQFT , SciPost Phys. 15 (2023) 151 [ 2304.13650 ].\n- [29] S. Collier, L. Eberhardt and M. Zhang, 3d gravity from Virasoro TQFT: Holography, wormholes and knots , 2401.13900 .\n- [30] S. R. Coleman, Black holes as red herrings: Topological fluctuations and the loss of quantum coherence , Nucl. Phys. B 307 (1988) 867.\n- [31] S. B. Giddings and A. Strominger, Loss of incoherence and determination of coupling constants in quantum gravity , Nucl. Phys. B 307 (1988) 854.\n- [32] S. B. Giddings and A. Strominger, Baby Universes, Third Quantization and the Cosmological Constant , Nucl. Phys. B 321 (1989) 481.\n- [33] S. W. Hawking, The Alpha parameters of wormholes , Phys. Scripta T 36 (1991) 222.\n- [34] R. Bousso and M. Tomaˇsevi\'c, Unitarity From a Smooth Horizon? , Phys. Rev. D 102 (2020) 106019 [ 1911.06305 ].\n- [35] D. Marolf and H. Maxfield, Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information , JHEP 08 (2020) 044 [ 2002.08950 ].\n- [36] R. Bousso and E. Wildenhain, Gravity/ensemble duality , Phys. Rev. D 102 (2020) 066005 [ 2006.16289 ].\n- [37] C. V. Johnson, The Microstate Physics of JT Gravity and Supergravity , 2201.11942 .  \n- [38] S. Hern\'andez-Cuenca, Wormholes and Factorization in Exact Effective Theory , 2404.10035 .\n- [39] J. de Boer, D. Liska, B. Post and M. Sasieta, A principle of maximum ignorance for semiclassical gravity , JHEP 02 (2024) 003 [ 2311.08132 ].\n- [40] V. Chandrasekaran, N. Engelhardt, S. Fischetti and S. Hern\'andez-Cuenca, A tale of two saddles , JHEP 11 (2022) 110 [ 2207.09472 ].\n- [41] R. Jackiw, Lower Dimensional Gravity , Nucl. Phys. B 252 (1985) 343.\n- [42] C. Teitelboim, Gravitation and Hamiltonian Structure in Two Space-Time Dimensions , Phys. Lett. B 126 (1983) 41.\n- [43] G. J. Turiaci and E. Witten, N = 2 JT supergravity and matrix models , JHEP 12 (2023) 003 [ 2305.19438 ].\n- [44] T. G. Mertens and G. J. Turiaci, Solvable models of quantum black holes: a review on Jackiw-Teitelboim gravity , Living Rev. Rel. 26 (2023) 4 [ 2210.10846 ].\n- [45] C. V. Johnson, Low Energy Thermodynamics of JT Gravity and Supergravity , 2008.13120 .\n- [46] K. Okuyama, Quenched free energy in random matrix model , JHEP 12 (2020) 080 [ 2009.02840 ].\n- [47] K. Okuyama, Quenched free energy from spacetime D-branes , JHEP 03 (2021) 073 [ 2101.05990 ].\n- [48] O. Janssen and M. Mirbabayi, Low-temperature entropy in JT gravity , JHEP 06 (2021) 074 [ 2103.03896 ].\n- [49] C. V. Johnson, On the Quenched Free Energy of JT Gravity and Supergravity , 2104.02733 .\n- [50] C. V. Johnson, Quantum Gravity Microstates from Fredholm Determinants , Phys. Rev. Lett. 127 (2021) 181602 [ 2106.09048 ].\n- [51] A. Sen, Logarithmic Corrections to N=2 Black Hole Entropy: An Infrared Window into the Microstates , Gen. Rel. Grav. 44 (2012) 1207 [ 1108.3842 ].\n- [52] A. Sen, Logarithmic Corrections to Rotating Extremal Black Hole Entropy in Four and Five Dimensions , Gen. Rel. Grav. 44 (2012) 1947 [ 1109.3706 ].\n- [53] A. Sen, Logarithmic Corrections to Schwarzschild and Other Non-extremal Black Hole Entropy in Different Dimensions , JHEP 04 (2013) 156 [ 1205.0971 ].\n- [54] S. Bhattacharyya, A. Grassi, M. Marino and A. Sen, A One-Loop Test of Quantum Supergravity , Class. Quant. Grav. 31 (2014) 015012 [ 1210.6057 ].\n- [55] A. Castro, V. Godet, F. Larsen and Y. Zeng, Logarithmic Corrections to Black Hole Entropy: the Non-BPS Branch , JHEP 05 (2018) 079 [ 1801.01926 ].\n- [56] A. Ghosh, H. Maxfield and G. J. Turiaci, A universal Schwarzian sector in two-dimensional conformal field theories , JHEP 05 (2020) 104 [ 1912.07654 ].\n- [57] L. V. Iliesiu and G. J. Turiaci, The statistical mechanics of near-extremal black holes , JHEP 05 (2021) 145 [ 2003.02860 ].\n- [58] M. Heydeman, L. V. Iliesiu, G. J. Turiaci and W. Zhao, The statistical mechanics of near-BPS black holes , J. Phys. A 55 (2022) 014004 [ 2011.01953 ].\n- [59] J. Boruch, M. T. Heydeman, L. V. Iliesiu and G. J. Turiaci, BPS and near-BPS black holes in AdS 5 and their spectrum in N = 4 SYM , 2203.01331 .  \n- [60] L. V. Iliesiu, S. Murthy and G. J. Turiaci, Black hole microstate counting from the gravitational path integral , 2209.13602 .\n- [61] L. V. Iliesiu, S. Murthy and G. J. Turiaci, Revisiting the Logarithmic Corrections to the Black Hole Entropy , 2209.13608 .\n- [62] N. Banerjee and M. Saha, Revisiting leading quantum corrections to near extremal black hole thermodynamics , JHEP 07 (2023) 010 [ 2303.12415 ].\n- [63] A. A. H., P. V. Athira, C. Chowdhury and A. Sen, Logarithmic correction to BPS black hole entropy from supersymmetric index at finite temperature , JHEP 03 (2024) 095 [ 2306.07322 ].\n- [64] G. J. Turiaci, New insights on near-extremal black holes , 2307.10423 .\n- [65] I. Rakic, M. Rangamani and G. J. Turiaci, Thermodynamics of the near-extremal Kerr spacetime , JHEP 06 (2024) 011 [ 2310.04532 ].\n- [66] D. Kapec, A. Sheta, A. Strominger and C. Toldo, Logarithmic Corrections to Kerr Thermodynamics , Phys. Rev. Lett. 133 (2024) 021601 [ 2310.00848 ].\n- [67] J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker et al., Black Holes and Random Matrices , JHEP 05 (2017) 118 [ 1611.04650 ].\n- [68] J. Liu, Spectral form factors and late time quantum chaos , Phys. Rev. D 98 (2018) 086026 [ 1806.05316 ].\n- [69] D. Stanford, More quantum noise from wormholes , 2008.08570 .\n- [70] P. Saad, S. H. Shenker and D. Stanford, A semiclassical ramp in SYK and in gravity , 1806.06840 .\n- [71] K. Okuyama, Eigenvalue instantons in the spectral form factor of random matrix model , JHEP 03 (2019) 147 [ 1812.09469 ].\n- [72] P. Saad, Late Time Correlation Functions, Baby Universes, and ETH in JT Gravity , 1910.10311 .\n- [73] K. Jensen, Chaos in AdS 2 Holography , Phys. Rev. Lett. 117 (2016) 111601 [ 1605.06098 ].\n- [74] J. Maldacena, D. Stanford and Z. Yang, Conformal symmetry and its breaking in two dimensional Nearly Anti-de-Sitter space , PTEP 2016 (2016) 12C104 [ 1606.01857 ].\n- [75] J. Engelsoy, T. G. Mertens and H. Verlinde, An investigation of AdS 2 backreaction and holography , JHEP 07 (2016) 139 [ 1606.03438 ].\n- [76] D. Bagrets, A. Altland and A. Kamenev, Sachdev-Ye-Kitaev model as Liouville quantum mechanics , Nucl. Phys. B 911 (2016) 191 [ 1607.00694 ].\n- [77] D. Stanford and E. Witten, Fermionic Localization of the Schwarzian Theory , JHEP 10 (2017) 008 [ 1703.04612 ].\n- [78] V. V. Belokurov and E. T. Shavgulidze, Exact solution of the Schwarzian theory , Phys. Rev. D 96 (2017) 101701 [ 1705.02405 ].\n- [79] T. G. Mertens, G. J. Turiaci and H. L. Verlinde, Solving the Schwarzian via the Conformal Bootstrap , JHEP 08 (2017) 136 [ 1705.08408 ].\n- [80] A. Kitaev and S. J. Suh, Statistical mechanics of a two-dimensional black hole , JHEP 05 (2019) 198 [ 1808.07032 ].\n- [81] Z. Yang, The Quantum Gravity Dynamics of Near Extremal Black Holes , JHEP 05 (2019) 205 [ 1809.08647 ].  \n- [82] C. V. Johnson, A Non-Perturbative Random Matrix Model of N = 2 JT Supergravity , 2306.10139 .\n- [83] H. W. Lin, J. Maldacena, L. Rozenberg and J. Shan, Looking at supersymmetric black holes for a very long time , SciPost Phys. 14 (2023) 128 [ 2207.00408 ].\n- [84] L. V. Iliesiu, M. Kologlu and G. J. Turiaci, Supersymmetric indices factorize , JHEP 05 (2023) 032 [ 2107.09062 ].\n- [85] P. H. Ginsparg and G. W. Moore, Lectures on 2-D gravity and 2-D string theory , in Theoretical Advanced Study Institute (TASI 92): From Black Holes and Strings to Particles , pp. 277-469, 10, 1993, hep-th/9304011 .\n- [86] B. Eynard, T. Kimura and S. Ribault, Random matrices , 1510.04430 .\n- [87] D. Anninos and B. Muhlmann, Notes on matrix models (matrix musings) , J. Stat. Mech. 2008 (2020) 083109 [ 2004.01171 ].\n- [88] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures , Phys. Rev. B 55 (1997) 1142 [ cond-mat/9602137 ].\n- [89] M. R. Zirnbauer, Riemannian symmetric superspaces and their origin in random-matrix theory , J. Math. Phys. 37 (1996) 4986 [ math-ph/9808012 ].\n- [90] A. Kitaev, Periodic table for topological insulators and superconductors , AIP Conf. Proc. 1134 (2009) 22 [ 0901.2686 ].\n- [91] E. P. Wigner, Characteristics vectors of bordered matrices with infinite dimensions ii , Annals of Mathematics 65 (1957) 203.\n- [92] F. J. Dyson, Statistical theory of the energy levels of complex systems. I , J. Math. Phys. 3 (1962) 140.\n- [93] F. J. Dyson, The Threefold Way. Algebraic Structure of Symmetry Groups and Ensembles in Quantum Mechanics , J. Math. Phys. 3 (1962) 1199.\n- [94] M. R. Zirnbauer, Symmetry classes , 1001.0722 .\n- [95] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder and S. Ryu, Classification of topological quantum matter with symmetries , Reviews of Modern Physics 88 (2016) .\n- [96] D. A. Ivanov, Random-matrix ensembles in p-wave vortices , arXiv: Superconductivity (2001) 253.\n- [97] P. Desrosiers and P. J. Forrester, Hermite and laguerre β -ensembles: Asymptotic corrections to the eigenvalue density , Nuclear Physics B 743 (2006) 307-332.\n- [98] I. Dumitriu and A. Edelman, Matrix models for beta ensembles , Journal of Mathematical Physics 43 (2002) 5830-5847.\n- [99] J. Wishart, The generalised product moment distribution in samples from a normal multivariate population , Biometrika (1928) 32.\n- [100] L. Pastur and M. Shcherbina, Eigenvalue Distribution of Large Random Matrices , Mathematical surveys and monographs. American Mathematical Society, 2011.\n- [101] H. Uhlig, On Singular Wishart and Singular Multivariate Beta Distributions , The Annals of Statistics 22 (1994) 395 .\n- [102] L. Erdos, S. Peche, J. A. Ramirez, B. Schlein and H.-T. Yau, Bulk universality for wigner matrices , 0905.4176 .  \n- [103] M. J. Bowick and E. Br\'ezin, Universal scaling of the tail of the density of eigenvalues in random matrix models , Physics Letters B 268 (1991) 21.\n- [104] D. J. Gross and A. A. Migdal, Nonperturbative Two-Dimensional Quantum Gravity , Phys. Rev. Lett. 64 (1990) 127.\n- [105] E. Brezin and V. A. Kazakov, Exactly Solvable Field Theories of Closed Strings , Phys. Lett. B 236 (1990) 144.\n- [106] M. R. Douglas and S. H. Shenker, Strings in Less Than One-Dimension , Nucl. Phys. B 335 (1990) 635.\n- [107] C. A. Tracy and H. Widom, Level spacing distributions and the Airy kernel , Phys. Lett. B 305 (1993) 115 [ hep-th/9210074 ].\n- [108] L. A. Pastur and V. A. Marˇcenko, Distribution of Eigenvalues for Some Sets of Random Matrices , Math. USSR Sb. 1 (1967) 457.\n- [109] C. A. Tracy and H. Widom, Level spacing distributions and the Bessel kernel , Commun. Math. Phys. 161 (1994) 289 [ hep-th/9304063 ].\n- [110] A. Perret and G. Schehr, Finite n corrections to the limiting distribution of the smallest eigenvalue of wishart complex matrices , Random Matrices: Theory and Applications 05 (2016) 1650001.\n- [111] C. V. Johnson and M. Usatyuk, God of the Gaps: Random matrix models and the black hole spectral gap , 2407.17583 .\n- [112] J. Ram\'ırez, B. Rider and B. Vir\'ag, Beta ensembles, stochastic airy spectrum, and a diffusion , Journal of the American Mathematical Society 24 (2011) 919-944.\n- [113] T. H. Baker, P. J. Forrester and P. A. Pearce, Random matrix ensembles with an effective extensive external charge , Journal of Physics A: Mathematical and General 31 (1998) 6087-6101.\n- [114] P. J. Forrester and A. K. Trinh, Functional form for the leading correction to the distribution of the largest eigenvalue in the gue and lue , Journal of Mathematical Physics 59 (2018) .\n- [115] P. J. Forrester and A. K. Trinh, Optimal soft edge scaling variables for the gaussian and laguerre even β ensembles , Nuclear Physics B 938 (2019) 621-639.\n- [116] L. Erd"os, Universality of wigner random matrices: a survey of recent results , Russian Mathematical Surveys 66 (2011) 507-626.\n- [117] T. Tao and V. Vu, Random matrices: The universality phenomenon for wigner ensembles , 1202.0068 .\n- [118] L. Pastur and M. Shcherbina, Universality of the local eigenvalue statistics for a class of unitary invariant random matrix ensembles , Journal of Statistical Physics 86 (1997) 109.\n- [119] P. Deift and D. Gioev, Universality at the edge of the spectrum for unitary, orthogonal, and symplectic ensembles of random matrices , Communications on Pure and Applied Mathematics 60 (2006) 867-910.\n- [120] P. Bourgade, L. Erdos and H.-T. Yau, Edge universality of beta ensembles , Communications in Mathematical Physics 332 (2014) 261-353.\n- [121] I. Dumitriu and A. Edelman, Global spectrum fluctuations for the β -hermite and β -laguerre ensembles via matrix models , Journal of Mathematical Physics 47 (2006) .  \n- [122] P. J. Forrester and N. S. Witte, Asymptotic forms for hard and soft edge general β conditional gap probabilities , Nuclear Physics B 859 (2012) 321-340.\n- [123] M. Bertola and M. Cafasso, Universality of the matrix airy and bessel functions at spectral edges of unitary ensembles , Random Matrices: Theory and Applications 06 (2017) 1750010.\n- [124] P. J. Forrester, The spectrum edge of random matrix ensembles , Nucl. Phys. B 402 (1993) 709.\n- [125] B. V. Bronk, Exponential ensemble for random matrices , Journal of Mathematical Physics 6 (1965) 228.\n- [126] T. Nagao and M. Wadati, Correlation Functions of Random Matrix Ensembles Related to Classical Orthogonal Polynomials , Journal of the Physical Society of Japan 60 (1991) 3298.\n- [127] T. Nagao and K. Slevin, Nonuniversal correlations for random matrix ensembles , Journal of Mathematical Physics 34 (1993) 2075 [ https://pubs.aip.org/aip/jmp/article-pdf/34/5/2075/19245150/2075 1 online.pdf ].\n- [128] B. Eynard, Topological expansion for the 1-Hermitian matrix model correlation functions , JHEP 11 (2004) 031 [ hep-th/0407261 ].\n- [129] B. Eynard and N. Orantin, Invariants of algebraic curves and topological expansion , Commun. Num. Theor. Phys. 1 (2007) 347 [ math-ph/0702045 ].\n- [130] B. Eynard, A short overview of the \'Topological recursion\' , 1412.3286 .\n- [131] O. Janssen, M. Mirbabayi and P. Zograf, Gravity as an ensemble and the moment problem , JHEP 06 (2021) 184 [ 2103.12078 ].\n- [132] A. B. Zamolodchikov and A. B. Zamolodchikov, Liouville field theory on a pseudosphere , hep-th/0101152 .\n- [133] V. Fateev, A. B. Zamolodchikov and A. B. Zamolodchikov, Boundary Liouville field theory. 1. Boundary state and boundary two point function , hep-th/0001012 .\n- [134] J. Teschner, Remarks on Liouville theory with boundary , PoS tmr2000 (2000) 041 [ hep-th/0009138 ].\n- [135] A. Blommaert, T. G. Mertens and H. Verschelde, Eigenbranes in Jackiw-Teitelboim gravity , JHEP 02 (2021) 168 [ 1911.11603 ].\n- [136] A. Goel, L. V. Iliesiu, J. Kruthoff and Z. Yang, Classifying boundary conditions in JT gravity: from energy-branes to α -branes , JHEP 04 (2021) 069 [ 2010.12592 ].\n- [137] P. Gao, D. L. Jafferis and D. K. Kolchmeyer, An effective matrix model for dynamical end of the world branes in Jackiw-Teitelboim gravity , JHEP 01 (2022) 038 [ 2104.01184 ].\n- [138] B. Post, J. van der Heijden and E. Verlinde, A universe field theory for JT gravity , JHEP 05 (2022) 118 [ 2201.08859 ].\n- [139] A. Altland, B. Post, J. Sonner, J. van der Heijden and E. P. Verlinde, Quantum chaos in 2D gravity , SciPost Phys. 15 (2023) 064 [ 2204.07583 ].\n- [140] J. Polchinski, Dirichlet Branes and Ramond-Ramond charges , Phys. Rev. Lett. 75 (1995) 4724 [ hep-th/9510017 ].\n- [141] N. Drukker and D. J. Gross, An Exact prediction of N=4 SUSYM theory for string theory , J. Math. Phys. 42 (2001) 2896 [ hep-th/0010274 ].  \n- [142] N. Drukker and B. Fiol, All-genus calculation of Wilson loops using D-branes , JHEP 02 (2005) 010 [ hep-th/0501109 ].\n- [143] K. Okuyama, Connected correlator of 1/2 BPS Wilson loops in N = 4 SYM , JHEP 10 (2018) 037 [ 1808.10161 ].\n- [144] K. Okuyama and K. Sakai, FZZT branes in JT gravity and topological gravity , JHEP 09 (2021) 191 [ 2108.03876 ].\n- [145] S. Hern\'andez-Cuenca and W. W. Weng, \'Work in progress.\'\n- [146] K. Okuyama, Replica symmetry breaking in random matrix model: a toy model of wormhole networks , Phys. Lett. B 803 (2020) 135280 [ 1903.11776 ].\n- [147] C. V. Johnson, Nonperturbative Jackiw-Teitelboim gravity , Phys. Rev. D 101 (2020) 106023 [ 1912.03637 ].\n- [148] K. Okuyama and K. Sakai, JT gravity, KdV equations and macroscopic loop operators , JHEP 01 (2020) 156 [ 1911.01659 ].\n- [149] K. Okuyama and K. Sakai, Multi-boundary correlators in JT gravity , JHEP 08 (2020) 126 [ 2004.07555 ].\n- [150] G. W. Gibbons and S. W. Hawking, Action Integrals and Partition Functions in Quantum Gravity , Phys. Rev. D 15 (1977) 2752.\n- [151] J. D. Bekenstein, Black holes and entropy , Phys. Rev. D 7 (1973) 2333.\n- [152] J. Preskill, P. Schwarz, A. D. Shapere, S. Trivedi and F. Wilczek, Limitations on the statistical description of black holes , Mod. Phys. Lett. A 6 (1991) 2353.\n- [153] C. Teitelboim, Action and entropy of extreme and nonextreme black holes , Phys. Rev. D 51 (1995) 4315 [ hep-th/9410103 ].\n- [154] G. W. Gibbons and R. E. Kallosh, Topology, entropy and Witten index of dilaton black holes , Phys. Rev. D 51 (1995) 2839 [ hep-th/9407118 ].\n- [155] S. W. Hawking, G. T. Horowitz and S. F. Ross, Entropy, Area, and black hole pairs , Phys. Rev. D 51 (1995) 4302 [ gr-qc/9409013 ].\n- [156] S. W. Hawking and G. T. Horowitz, The Gravitational Hamiltonian, action, entropy and surface terms , Class. Quant. Grav. 13 (1996) 1487 [ gr-qc/9501014 ].\n- [157] J. M. Maldacena, J. Michelson and A. Strominger, Anti-de Sitter fragmentation , JHEP 02 (1999) 011 [ hep-th/9812073 ].\n- [158] D. N. Page, Thermodynamics of near extreme black holes , hep-th/0012020 .\n- [159] A. Dabholkar, A. Sen and S. P. Trivedi, Black hole microstates and attractor without supersymmetry , JHEP 01 (2007) 096 [ hep-th/0611143 ].\n- [160] N. Engelhardt and A. C. Wall, Decoding the Apparent Horizon: Coarse-Grained Holographic Entropy , Phys. Rev. Lett. 121 (2018) 211301 [ 1706.02038 ].\n- [161] N. Engelhardt and A. C. Wall, Coarse Graining Holographic Black Holes , JHEP 05 (2019) 160 [ 1806.01281 ].\n- [162] S. N. Solodukhin, The Conical singularity and quantum corrections to entropy of black hole , Phys. Rev. D 51 (1995) 609 [ hep-th/9407001 ].  \n- [163] S. N. Solodukhin, On \'Nongeometric\' contribution to the entropy of black hole due to quantum corrections , Phys. Rev. D 51 (1995) 618 [ hep-th/9408068 ].\n- [164] A. Sen, Extremal black holes and elementary string states , Mod. Phys. Lett. A 10 (1995) 2081 [ hep-th/9504147 ].\n- [165] R. B. Mann and S. N. Solodukhin, Universality of quantum entropy for extreme black holes , Nucl. Phys. B 523 (1998) 293 [ hep-th/9709064 ].\n- [166] R. K. Kaul and P. Majumdar, Logarithmic correction to the Bekenstein-Hawking entropy , Phys. Rev. Lett. 84 (2000) 5255 [ gr-qc/0002040 ].\n- [167] S. Das, P. Majumdar and R. K. Bhaduri, General logarithmic corrections to black hole entropy , Class. Quant. Grav. 19 (2002) 2355 [ hep-th/0111001 ].\n- [168] A. Sen, Black Hole Entropy Function, Attractors and Precision Counting of Microstates , Gen. Rel. Grav. 40 (2008) 2249 [ 0708.1270 ].\n- [169] S. Banerjee, R. K. Gupta and A. Sen, Logarithmic Corrections to Extremal Black Hole Entropy from Quantum Entropy Function , JHEP 03 (2011) 147 [ 1005.3044 ].\n- [170] S. Pal and J. Qiao, Lightcone Modular Bootstrap and Tauberian Theory: A Cardy-like Formula for Near-extremal Black Holes , 2307.02587 .\n- [171] H. K. Kunduri, J. Lucietti and H. S. Reall, Near-horizon symmetries of extremal black holes , Class. Quant. Grav. 24 (2007) 4169 [ 0705.4214 ].\n- [172] P. Nayak, A. Shukla, R. M. Soni, S. P. Trivedi and V. Vishal, On the Dynamics of Near-Extremal Black Holes , JHEP 09 (2018) 048 [ 1802.09547 ].\n- [173] U. Moitra, S. P. Trivedi and V. Vishal, Extremal and near-extremal black holes and near-CFT 1 , JHEP 07 (2019) 055 [ 1808.08239 ].\n- [174] J. M. Maldacena and A. Strominger, Universal low-energy dynamics for rotating black holes , Phys. Rev. D 56 (1997) 4975 [ hep-th/9702015 ].\n- [175] A. Castro and F. Larsen, Near Extremal Kerr Entropy from AdS(2) Quantum Gravity , JHEP 12 (2009) 037 [ 0908.1121 ].\n- [176] A. Castro, F. Larsen and I. Papadimitriou, 5D rotating black holes and the nAdS 2 /nCFT 1 correspondence , JHEP 10 (2018) 042 [ 1807.06988 ].\n- [177] U. Moitra, S. K. Sake, S. P. Trivedi and V. Vishal, Jackiw-Teitelboim Gravity and Rotating Black Holes , JHEP 11 (2019) 047 [ 1905.10378 ].\n- [178] A. Almheiri and J. Polchinski, Models of AdS 2 backreaction and holography , JHEP 11 (2015) 014 [ 1402.6334 ].\n- [179] A. Almheiri and B. Kang, Conformal Symmetry Breaking and Thermodynamics of Near-Extremal Black Holes , JHEP 10 (2016) 052 [ 1606.04108 ].\n- [180] G. S\'arosi, AdS 2 holography and the SYK model , PoS Modave2017 (2018) 001 [ 1711.08482 ].\n- [181] F. Larsen and Y. Zeng, Black hole spectroscopy and AdS 2 holography , JHEP 04 (2019) 164 [ 1811.01288 ].\n- [182] S. Sachdev, Universal low temperature theory of charged black holes with AdS 2 horizons , J. Math. Phys. 60 (2019) 052303 [ 1902.04078 ].  \n- [183] J. Hong, F. Larsen and J. T. Liu, The scales of black holes with nAdS 2 geometry , JHEP 10 (2019) 260 [ 1907.08862 ].\n- [184] A. Castro and V. Godet, Breaking away from the near horizon of extreme Kerr , SciPost Phys. 8 (2020) 089 [ 1906.09083 ].\n- [185] A. M. Charles and F. Larsen, A one-loop test of the near-AdS 2 /near-CFT 1 correspondence , JHEP 07 (2020) 186 [ 1908.03575 ].\n- [186] W. Fu, D. Gaiotto, J. Maldacena and S. Sachdev, Supersymmetric Sachdev-Ye-Kitaev models , Phys. Rev. D 95 (2017) 026009 [ 1610.08917 ].\n- [187] T. Li, J. Liu, Y. Xin and Y. Zhou, Supersymmetric SYK model and random matrix theory , JHEP 06 (2017) 111 [ 1702.01738 ].\n- [188] T. Kanazawa and T. Wettig, Complete random matrix classification of SYK models with N = 0 , 1 and 2 supersymmetry , JHEP 09 (2017) 050 [ 1706.03044 ].\n- [189] F. Sun and J. Ye, Periodic Table of the Ordinary and Supersymmetric Sachdev-Ye-Kitaev Models , Phys. Rev. Lett. 124 (2020) 244101 [ 1905.07694 ].\n- [190] A. A. Migdal, Loop Equations and 1/N Expansion , Phys. Rept. 102 (1983) 199.\n- [191] N. I. Muskhelishvili, Singular Integral Equations . Courier Corporation, 1977.\n- [192] L. A. Pastur, Spectral and Probabilistic Aspects of Matrix Models , pp. 207-242. Springer Netherlands, Dordrecht, 1996. doi.org/10.1007/978-94-017-0693-3 10.\n- [193] J. Harnad, Random matrices, random processes and integrable systems . Springer, 1, 2011, 10.1007/978-1-4419-9514-8.\n- [194] P. Deift, T. Kriecherbauer and K.-R. McLaughlin, New results on the equilibrium measure for logarithmic potentials in the presence of an external field , Journal of Approximation Theory 95 (1998) 388.'}