An effective matrix model for dynamical end of the world branes in Jackiw-Teitelboim gravity

We study Jackiw-Teitelboim gravity with dynamical end of the world branes in asymptotically nearly AdS$_2$ spacetimes. We quantize this theory in Lorentz signature, and compute the Euclidean path integral summing over topologies including dynamical branes. The latter will be seen to exactly match with a modification of the SSS matrix model. The resolution of UV divergences in the gravitational instantons involving the branes will lead us to understand the matrix model interpretation of the Wilsonian effective theory perspective on the gravitational theory. We complete this modified SSS matrix model nonperturbatively by extending the integration contour of eigenvalues into the complex plane. Furthermore, we give a new interpretation of other phases in such matrix models. We derive an effective $W(\Phi)$ dilaton gravity, which exhibits similar physics semiclassically. In the limit of a large number of flavors of branes, the effective extremal entropy $S_{0,\text{eff}}$ has the form of counting the states of these branes.


Introduction
Euclidean wormholes have been at the center of a spate of recent progress in understanding quantum gravity [1][2][3][4][5][6][7]. These non-perturbative contributions were shown to reproduce late time behavior of the spectral form associated to level repulsion in chaotic systems (referred to as the ramp) [4], and the Page curve for entanglement entropy of an AdS black hole evaporating into an external bath [1,2]. In pure Jackiw-Teitelboim (JT) gravity [8,9] the non-perturbative Euclidean path integral can be computed exactly [10], and precisely matched with a particular double-scaled matrix model.
The latter fact highlights the curious aspect that the quantum mechanics dual is an ensemble average, rather than a particular Hamiltonian theory. This seems at odds with AdS/CFT in higher dimensions. In particular, although there are no further constraints of consistency on positive hermitian Hamiltonians in quantum mechanics, conformal field theories in higher dimensions are highly constrained by microlocality and come only in sporadic families. Therefore there cannot exist continuous ensembles of conformal field theories with a number of parameters that scales with the dimensionality of the Hilbert space, as happens in the matrix model for JT gravity, otherwise any potential matrix model dual of higher dimensional AdS quantum gravity would allow nonlocal heavy operators that does not obey crossing symmetry because their matrix elements could be tuned arbitrarily by the overwhelming number of parameters, and that contradict with AdS/CFT because bulk dual of heavy operators are black holes localized in a finite region.
For these reasons it is of great interest to know if JT gravity itself can be the long distance effective theory of some UV modified theory that behaves more conventionally, with an ordered dual possessing a discrete spectrum, as might be expected to arise in a string compactification to 1+1 dimensions. We will not answer that question here, but will explore consequences of a minimal ingredient that such a UV "completion" of JT gravity would require: dynamical branes that can end spacetime.
In 1+1 dimensions, the fundamental objects that could form the microstates of black holes are codimension 1 branes, in other words, end of the world (EOW) branes. Dynamical objects of that type are also required to solve the factorization problem [11] starting from canonical quantization of Lorentzian JT theory. At an even most basic level, there are simply no Lorentzian configurations in pure JT gravity with a single nearly AdS 2 boundary that could appear in the Hilbert space of a putative dual quantum mechanics.
EOW branes could have a variety of microscopic realizations in a compactification to JT gravity, including intrinsic end of the world branes like the Horava-Witten brane in M-theory and Kaluza-Klein solitons where part of an internal manifold smoothly shrinks. In the long distance effective theory these can be treated as boundary conditions in JT gravity. We will study the simplest possibility, which is characterized simply by the brane tension (equivalently the rest mass of this 0+1 dimensional object).
Branes of this type were used in [2] to model black hole microstates of pure JT gravity, however in that context they were not treated as dynamical objects added to the theory. Because of this, the loops of EOW branes were not considered in [2] even though the number of flavors K was as large as e S 0 . The Page curve of [2] indicates a phase transition when K is on the order of e S 0 . Here we will be interested in the effects of summing over EOW branes in the path integral, which will be important when K is of order e S 0 .
We will see that positive tension EOW branes are always cloaked behind horizons. Nevertheless, the EOW branes modify the spectral density because they act as gravitational instantons that correct the long distance effective action in the form of a more general W (Φ) gravity.
Similarly to [11], the classical phase space may be parametrized in a geodesic slice gauge in terms of the (renormalized) length L of the spatial slice and its conjugate momentum, leading to a Morse potential quantum mechanics. The spectrum of the Hamiltonian is a continuum, so it is unsurprising that the theory with EOW branes remains an ensemble.
The continuous spectrum of the Morse potential quantum mechanics, which is associated to the L → ∞ limit, leads to a logarithmic divergence in the contributions of gravitational instantons. The large L limit can equivalently be understood as a b → 0 limit of high temperature in the reference frame of the EOW brane, where b is the length of the EOW brane loop. As such the calculation is sensitive to UV physics deep in the throat; the large red shift allows it to contribute at a fixed energy as measured at the AdS boundary. Therefore an EFT understanding is crucial for the proper treatment of this divergence.
In the gravity calculation, we obtain finite results by using a regulated EOW brane such that the b → 0 pole is cancelled and the theory becomes pure JT in the far UV. The low energy spectral density ρ(E) becomes universal as the regulator is removed, up to a single parameter. One characterization of that UV sensitive parameter is the zero point energy, E 0 , at which the spectrum begins. In this way, it behaves like a relevant parameter in the Wilsonian sense.
The effect of the EOW branes in the dual quantum mechanics is to introduce K vectors in the SSS matrix model [10] (as in [2]). These are the states in the Hilbert space produced by a given EOW brane. The branes here are dynamical objects in the gravity theory, and one must integrate over the vector degrees of freedom in the matrix model.
For the purpose of computing the spectral density and its correlation functions, the vectors can be integrated out, leading to a modified matrix ensemble for the Hamiltonian, whose matrix potential differs from that of pure JT gravity by δV (E). The deformation of the matrix model potential by various types of EOW branes has been considered in a recent work [12].
The same δV can equally be computed from the trumpet partition function with fixed AdS boundary energy, which is the inverse Laplace transformation of Z trumpet (b, β). The full path integral is a sum over the number of closed loops of EOW brane boundaries, which is reproduced by the exponentiation of the additional potential, e −KTrδV (H) .
There is a trivial UV divergence in δV , which can be absorbed into the overall normalization of the ensemble measure. It is straightforward to compute the change in the tree-level spectral density δρ induced by the change in the potential δV . However, the large E behavior of δV (E) implies that the integral of δρ(E) over its support (which is noncompact) diverges. More severely, without the appropriate fine tuning of the UV sensitive parameter, E 0 would go to −∞. This means that an infinite number of eigenvalues that would have been pushed to infinity in the double-scaling limit of the SSS matrix model remain at finite energies.
Regulating the theory at short distance leads to a modified δV that decays at large energy, allowing a δρ that is normalizable. This is a special feature that is associated to this model's rapid approach to the pure JT spectrum at high energies. We exactly match the result with the gravitational calculation. More importantly, from the IR point of view, E 0 is a free parameter in both the matrix model and gravity sides, which is sensitive to the details of the bulk UV physics. Figure 1: The phase diagram of matrix model of JT gravity with K flavors of EOW branes. Green region is "Y" shaped phase, blue region is one-cut phase and red curve is critical line.
The contribution of near cusps, which arise in the limit of an infinitely long trumpet, corresponds to an effective W (Φ) ∝ e −2πΦ in the general dilaton gravity considered in [13,14]. The matrix model potential in the double scaling limit is in fact unchanged from SSS, and the only difference is via the IR parameter E 0 which determines how the double scaling limit is taken.
For large K e S 0 in the matrix model with heavy EOW branes, the effective extremal entropy S 0,eff defined as ρ(E) ≈ e S 0,eff √ E − E 0 for E − E 0 1 scales as log K. This is a form of induced gravity, in which K flavors of heavy EOW branes count the microscopic states of black hole.
As pointed out in [10], the matrix model dual to pure JT gravity is non-perturbatively unstable because the effective potential becomes arbitrarily negative away from the support of the spectrum. This non-perturbative instability causes negativity in the spectrum of JT gravity with too many flavors of deficit angles [15] and in JT gravity with EOW branes as well. To have a well-defined matrix model, one needs to extend the contour of each eigenvalue through the largest saddle point of the effective potential on the real axis and into the complex plane [10]. This completion promotes H to a complex matrix because an order e −e S 0 fraction of eigenvalues will be complex even in the pure JT model.
We will adopt this completion and study how the spectrum changes as we increase the number of flavors K of EOW branes. As K increases, the qualitative behavior depends on how the UV divergence of the EOW branes is renormalized, which comes with a renormalization parameter λ. There are three cases to consider, corresponding to the zero point energy E 0 being positive (λ < 0), zero (λ = 0) or negative (λ > 0) when K is slightly positive. See Fig. 1 as the phase diagram. For the first case, there is no critical K and no phase transition. For second case, a continuous phase transition occurs and E 0 moves to positive after K going beyond a critical value. For the last case, instead of a negative spectrum, our non-perturbative completion implies that the spectrum undergoes a continuous phase transition to a "Y" shaped on the complex energy plane when K is over a critical value (see Fig. 10).
We interpret these e S 0 order of complex eigenvalues as unstable black hole states that could decay to lower energy objects other than EOW branes. This matches the Lorentzian analysis in the effective W (Φ) dilaton gravity [13], for which the spectrum is unbounded from below. This is also related to the existence of a minimal temperature below which no stable black hole exists.
Because of this, we regard JT gravity with such EOW branes as incomplete, requiring other stable objects with lower energy (such as Dirichlet-Neumann branes [12]) to which these unstable black holes decay. Using the effective W (Φ) gravity description, we explicitly show how this occurs as a Hawking-Page phase transition.
The paper is organized as follows. In Section 2, we find the phase space of JT gravity with one dynamical EOW brane and quantize it canonically. We compute the partition function of this system and identify the measure of EOW branes in the Euclidean path integral. Using this measure to sum over arbitrary numbers of EOW branes, we compute the tree-level spectral density. For the E 0 = 0 case, we find the effective W (Φ) dilaton gravity by integrating out the EOW branes. In Section 3, we derive the change in the matrix model potential δV induced by the EOW branes and solve for the one-cut solution of the spectral density that matches with our gravitational computation. In Section 4, we study the "Y" shaped spectrum in the matrix model and the phase transition from the one-cut solution when the number of flavors, K, of EOW branes exceeds a critical value. In Section 5, we study the effective W (Φ) dilaton gravity for a gas of cusps that is related to K heavy EOW branes. By requiring smoothness of the Euclidean metric, we find complex saddles when K is beyond the critical value. These complex saddles exhibit similar "Y" shaped spectra. We interpret the complex energies as unstable black holes and study the Hawking-Page phase transition after including lower energy Dirichlet-Neumann branes. We conclude in Section 6 with a few discussions.
2 Ends of the world in 2d gravity 2

.1 Classical solution and phase space with a boundary brane
The JT gravity [8,9] action with boundaries is where S 0 = φ 0 4G is the extremal entropy and κ = (8πG) −1 . In this paper, we will set κ = 1 for notational simplicity. u and v are affine parameters along the AdS boundary and EOW brane respectively. As analyzed in Appendix A, the equations of motion are and the boundary conditions are For an AdS boundary, we fix the asymptotic metric h ab δg ab = 0 and the value of dilaton δΦ = 0; for an EOW brane, we allow the metric and dilaton to fluctuate but impose ∇ a n a = 0, n a ∂ a Φ = µ (2.4) The general solution to the equations of motion can be written as where σ ∈ [0, π], T ∈ R and α, β, γ are real numbers. Physically, Φ + φ 0 must be nonnegative as it would represent the area of codimension 2 surfaces in higher dimensions if JT gravity were obtained via dimensional reduction. As the solution is periodic in T , the physical region Φ + φ 0 ≥ 0 appears periodically along T axis. There is a SL(2) SO(2, 1) isometry for AdS 2 metric. To see this, we can write AdS 2 as a hyperplane in higher dimension Then the dilaton solution can be written in a simpler form Under SO(2, 1) transformation, α, β, γ will change according to the action on V µ . However, there is one invariant describing the solution, whose sign will separate the solutions into two types. Indeed, the existence of a saddle for Φ depends on γ 2 − α 2 . Differentiating with respect to T and σ, we get If |γ| ≤ |α|, the saddle exists and if |γ| > |α|, it does not. In this paper, we will only consider cases with a Φ saddle, which is the analog of the near extremal horizon "area" in JT gravity. The SO(2, 1) is a gauge symmetry and we fix it by choosing the dilaton solution to be Figure 2: The solution of JT gravity. The light green region is physical in which Φ + φ 0 ≥ 0. Red curves are two AdS boundaries with T ∈ [−π/2, π/2]. The dashed lines are horizon. Blue curve is the geodesic of EOW brane.
The physical region with AdS boundary is chosen to be T ∈ [−π/2, π/2] and Φ h is the saddle value of Φ (see Fig. 2). The AdS boundary condition is [16] For fixed Φ = A value, its solution is given by In large A limit, we would like to put AdS boundary near σ = π. In small expansion with fixed metric, we have The EOW brane boundary condition is given in (2.4). For any brane world line σ − f (T ) = 0, the normal vector is The normal derivative to Φ on the brane is On the other hand, ∇ a n a = 0 implies that it is a geodesic whose general solution is Note that we require the world line to be timelike, namely |f (T )| < 1. This restricts r ∈ [−1, 1]. Using coordinate Y µ , we can rewrite geodesic solution (2.19) as where U µ can be rescaled with any nonzero real number without change the solution. As SO(2, 1) is a gauge symmetry, we are free to choose V µ as above to partly fix it. In (2.12), we rotate V µ to (0, Φ h , 0) and clearly there is an unfixed SO(1, 1) between the V −1 and V 1 components. This SO(1, 1) can be further fixed by rotating the U −1 and U 1 components to set U −1 = 0 (note that |r sin θ| < 1 makes this always possible). This is equivalent to set sin θ = 0 (θ = 0). Hence, after fixing all gauge symmetries, the solution for the geodesic is Plugging this solution into n a ∂ a Φ = µ, we get As shown above, the EOW brane boundary condition completely determines its worldline from the parameters µ and Φ h . There is another way to characterize the brane geodesic from the AdS boundary point of view, which will be important for our phase space description. It is the length of a spacelike geodesic shooting from AdS boundary at u = 0 and ending normally on the brane. It is simplest to calculate this in the Y µ coordinate. Note that different geodesics can be transformed to each other via SO(2, 1) isometry. In particular, the brane geodesic (2.21) can be written as On the other hand, the geodesic U µ Y µ = 0 is very simple, namely σ = π/2. We can calculate the spacelike geodesic connecting it with boundary and then do a SO(2, 1) transformation to find the case related to (2.21). Note that orthogonality and geodesic length is preserved under isometry. It is clear that for σ = π/2, the orthogonal spacelike geodesic is T = T 0 for any T 0 . This spatial slice intersects the AdS boundary at (T 0 , σ 0 ). The geodesic length from the brane to (T 0 , σ 0 ) is To get the spatial slice orthogonal to (2.21), we simply do the transformation 25) Under this transformation, cos σ 0 transforms as The boundary location (σ 0 , T 0 ) is given by (2.14) with A = φ b / . Taking this into (2.24) gives divergent result and needs to be renormalized as [11] Using the expression for r in (2.22) and T 0 in (2.16), we have As the lower bound of Φ h is zero, we must require µ ≥ 0 for L to be well-defined. The AdS boundary stress tensor is where γ µν is the induced metric. In 2d case, there is only one component of stress tensor, namely Hamiltonian. Using (A.10), we get Using (2.14) (AdS boundary has A = φ b / ) and (2.18) (here we need to flip the sign for n a pointing outward), we can evaluate it on AdS boundary as This Hamiltonian is nonnegative as the ADM energy of a black hole should be. On the other hand, the brane Hamiltonian is vanishing because we choose the brane boundary condition such that T vv brane ∝ n a ∂ a Φ − µ = 0. This is consistent with gravity being dynamical on the EOW brane.
Similarly to the two-sided 2D JT gravity [11], the phase space is two dimensional, which is characterized by u 0 and H, where u 0 is the boundary time constant corresponding to T = 0 slice. As H is the Hamiltonian on AdS boundary, its canonical conjugate is time translation. This implies the symplectic form in phase space is The dynamics in phase space is given bẏ The symplectic form can be written as where we see that w is the conjugate coordinate for Φ h . Note that using parameters w and Φ h is not quite a good description of the phase space because Φ h is restricted to be a nonnegative number. In order to find nice phase space coordinates with range R 2 , we need to do a canonical transformation.
The canonical transformation to the L variable is easy to find. Solving for w in (2.28), we find where the sign of w depends on that of u 0 . This implies that Solving for Φ 2 h leads to the Hamiltonian in terms of L and P This is the Hamiltonian with Morse potential. It is obvious that the Hamiltonian is nonnegative for all L only when µ > 0. Indeed, for µ < 0, we can do a similar phase space analysis and end up with the same Hamiltonian. In that case, Hamiltonian is negative for L larger than a critical value but still lower bounded. In particular, there is a stable minimal energy point located at P = 0, L = − log |µ| with ground energy E = −µ 2 /φ b . Such Hamiltonian allows bound states, whose geometric meaning is a naked EOW brane, rather than a black hole. In this paper, we will mainly focus on µ > 0 case.

Quantization with a boundary brane
As the Hamiltonian is simply a particle in the Morse potential, its quantization is straightforward by replacing P → −i∂ L . It follows that the energy eigen-functions f E (L) obey Using a new variable z = 4e −L , we can rewrite the equation as For µ > 0, the spectrum of solutions is continuous. The general solution is given by Whittaker function Note that here we restrict k ≥ 0 because of identity W a,b (z) = W a,−b (z). The normalization is given by a flat measure integration over L ∈ R After some algebra, we can work out the normalized eigen-function The propagator for the Morse potential quantum mechanics is If we set L 1 = L 2 and integrate it over R, we will get the partition function of JT gravity with a single EOW brane, which is divergent. This is because the Morse potential quantum mechanics has a continuous spectrum. As we will now show, we may alternatively interpret this divergence as a UV divergence associated to an EOW brane loop that shrinks to zero size. To see this, use identity (6.647-1 in [17]) to write (2.43) as It follows that the partition function with a single EOW brane is where we used identity (6.621-3 in [17]). Redefining w = e −b/2 for positive b, we have where Z trumpet (β, b) is the Euclidean path integral of the trumpet bounded by one AdS boundary of regularized length β and one geodesic boundary of length b in [10], The integral in (2.47) represents a Euclidean spacetime that ends on a geodesic EOW brane in the bulk, and the length of the EOW brane b is integrated with a measure M(b) given by This integral clearly diverges due to its behavior as b → 0.
To obtain a well-defined path integral, we need to regulate the small b behavior of the measure M(b) such that (2.47) converges. In particular, we require that the regulated M(b) is bounded as b → 0. The UV divergence in (2.47) can be cancelled by a cusp-like counterterm, which physically corresponds to a closed geodesic whose length is zero in the limit that the regulator is removed. This is a modification of the UV physics that regulates the contributions of EOW branes to the path integral of the effective IR theory. The specific form of the regulator and the cusp-like counterterm is far from unique. One option is to regulate the UV divergence of (2.47) by cutting off the b integral at b = , adding a counterterm proportional to a trumpet smeared over b ∈ [0, ], and then taking → 0. The regulated path integral with one AdS boundary and one EOW brane becomes where λ represents the finite part of the cusp-like counterterm. We can thus make the following replacement in the definition of M(b), which is bounded for all b ≥ 0 as required. As usual, the renormalization procedure forces us to view λ as an additional parameter of the theory. The value of λ depends on UV physics that becomes relevant in the small b limit. Modifying this UV physics corresponds to modifying the large-eigenvalue behavior of the potential of the dual matrix model, as we will discuss in Section 3.3.3. We will find it convenient to express our results in terms of the inverse Laplace transform of M(b), which we denote by m(α). We have that The m(α) corresponding to the regulated measure (2.51) could be computed similarly, though its explicit form is not important for our calculations. We will replace the regulated measure with the unregulated measure whenever doing so leads to a finite result. The property (2.53) allows us to exchange the order of a limit and integral in (3.47).

Dynamical branes in path integral
In this section, we compute the Euclidean path integral with one AdS boundary and arbitrary numbers of EOW branes at genus zero. In general, given a fixed number of AdS boundaries, we are required to sum over all geometries with arbitrary numbers of handles and EOW branes (see Fig.  3). As shown in [10], path integrals can be computed by gluing trumpets to hyperbolic Riemann surfaces with closed geodesic boundaries. An EOW brane corresponds to a geodesic boundary whose length b is integrated over with measure M(b), given in (2.49). The b → 0 divergence can be regulated as shown in (2.50). In the Euclidean action, there is an Einstein-Hilbert term −S 0 χ, where χ is Euler characteristic χ = 2 − 2g − n of Riemann surface, where g is genus and n is number of boundaries. This topological term endows every Riemann surface with a weight e S 0 χ . For fixed genus, adding one more EOW brane means increasing n by one and leads to a factor of + + + + Z(β)= + Figure 3: The partition function Z(β) in which all genera and EOW brane loops are summed. Red curve is AdS boundary and blue curves are EOW branes.
e −S 0 . On the other hand, if we have K flavors of EOW brane, each loop will contribute a factor of K (modulo a possible permutation symmetry that we will specify later). Therefore, for an O(1) number of flavors, any effects of the EOW brane loops are suppressed by e −S 0 . To enhance the effects of the EOW branes, we will assume K ∼ O(e S 0 ). Furthermore, we will assume that S 0 is large so that we may restrict our attention to genus zero surfaces. For simplicity, we assume all flavors of EOW branes have the same value of µ. Our method for resumming EOW branes in the Euclidean path integral will closely follow the method used in [15]. The leading term is a disk where there is no EOW brane. As computed in [10], the disk partition function is The zero genus partition function is a series where Z n (β) is the Euclidean path integral over the surface with n EOW brane loops and one AdS boundary In this formula, V 0,n+1 (b, b 1 , · · · , b n ) is the WP volume for genus zero and n + 1 holes computed in [18], where J 0 (x) is Bessel function of first kind and u(x) is defined implicitly via where I 1 (x) is modified Bessel function of first kind. For a given x, this equation has infinitely many solutions. We take the largest u because we need u(0) = 0 in (2.58). This defines u(x) smoothly in a neighborhood of x = 0, and for all negative x. However, for negative u, √ uI 1 (2π √ u) oscillates (see Fig. 4) with increasing amplitude as |u| 1/4 , which leads to a piecewise continuous function u(x). As we will show later, for certain values of the parameters, this discontinuity will imply that the sum over EOW branes is divergent. This multi-valueness of u(x) also leads to non-perturbative instability as analyzed in [19]. In (2.57), C n is a symmetry factor that accounts for identical EOW branes. We assume EOW branes to be indistinguishable for the same flavor but distinguishable for different flavors. This leads to Let us denote the last infinite sum in (2.56) asẐ(β). Using (2.58), it can be written aŝ and defined Note that the boundedness of the regulated measure M(b) for small b also guarantees a finite f (u) whereas using unregulated measure (2.49) leads to divergent f (u). To sum this series to a closed form, we need to use the Lagrange inversion theorem (see Appendix B). It turns out that where g(x) is defined implicitly as and a must be a parameter such that u(a) → +∞. Clearly, we choose a = −∞. It turns out that the second and third terms in (2.64) are exactly first two terms in (2.56). For the second term, we can change the variable x to u to get From (2.59), we have and (2.66) becomes For the third term, we can recover the Bessel functions using (2.62) and (2.63), and change integration variable to u to get (2.69) Changing variable u → x 2 and using identity we can rewrite it as Redefining ξ(−x) ≡ u(x + g(x)) and using (2.59) and (2.65), we can write the partition function as (2.73) Equation (2.73) is the "string equation" that shows how the inclusion of EOW branes affects the partition function. Starting from the regulated and renormalized measure in (2.51) and taking the → 0 limit, we find that f (u) becomes f λ (u), which is defined by and H(x) is the Harmonic number. As the constant λ 0 has no physical meaning, we will shift λ → λ + λ 0 in the rest of paper for notational simplicity.
In the limit that the regulator is removed, any regulated measure leads to the same f λ (u) for some value of renormalization parameter λ. For example, we used the cusp-like counterterm such that f λ (0) = λ in (2.74). The derivative of f λ (u) agrees with the derivative of (2.63) with m(α) taken to be the unregulated measure (2.54). In the limit that the regulator is removed, the string equation becomes It is clear from the string equation that our theory contains three independent parameters: K, µ, and λ.

Spectral density
We perform an inverse Laplace transform to get the spectral density where the zero point energy E 0 is given by x = 0 in (2.76) The largest intersection point gives zero point energy E 0 . (a) yellow, green and red means increasing K and E 0 moves rightward as K increases; (b) yellow is K < K = cr , green is critical K = K = cr and red is K > K = cr . E 0 = 0 when K ≤ K = cr , and E 0 > 0 when K > K = cr ; (c) yellow is K > K > cr , green is critical K = K > cr , red is K < K > cr and purple is for K too large such that no intersection exists (this purple curve also has smaller µ than the other three). E 0 moves leftward when K increases and has a jump when K > K > cr . In the plot, we set 2φ b = 1 and e S 0 = 1.
As f λ (u) has singularities at −(1/2 + µ + n) 2 for n ∈ N, we will assume 2φE 0 > −(1/2 + µ) 2 throughout this paper. Using (2.63) and (2.51), we can also write the second term in (2.77) as In → 0 limit, the derivative f λ (u) = f (u) is finite, which means that the dependence of ρ(E) on the renormalization parameter λ is through the zero point energy E 0 only. We are interested in how the spectrum changes when K is increased from zero. For K = 0, the spectrum reduces to the pure JT gravity result For nonzero K, because √ uI 1 (2π √ u) passes through the origin as shown in Fig. 4, we can organize the problem into three cases based on the value of λ: f λ (0) < 0, f λ (0) = 0 and f λ (0) > 0.
Note that f λ (u) is a monotonically decreasing function. When f λ (0) < 0, λ < 0, the largest solution to (2.78) is positive, which implies E 0 > 0 (see Fig. 5a). As K increases, E 0 increases as well and the spectrum is pushed to right (see Fig. 6a).
At this critical point, the spectral density near the edge scales like ρ(E) ∼ E 3/2 rather than the generic case where ρ(E) ∼ (E − E 0 ) 1/2 as shown in Fig. 6b. To see this, let us define At the critical point, we have G(0) = G (0) = 0. Close to the critical point, expanding E around E 0 , we have where a 1 = 0 at the critical point. As we discussed before, E 0 > 0 for K > K = cr , and E 0 = 0 for K < K = cr . Taking this into (2.77), we find that the spectral edge scales as where the leading order scaling is (E − E 0 ) 1/2 unless a 1 = 0. When f λ (0) > 0, λ > 0, the largest solution to (2.78) is negative, which implies E 0 < 0 (see Fig. 5c). As K increases, E 0 decreases as well and the spectrum is pushed to left (see Fig.  6c). However, as √ xI 1 (2π √ x)/(2π) is an oscillatory function when x < 0, there will be another critical point at some negative x. Similar to (2.84), we can show that at this critical point the near edge spectrum scales as (E − E 0 ) 3/2 . Though K > cr does not have an analytic closed form, we can easily confirm that E 0 would have a jump when K > K > cr because of the oscillatory feature of 5c). However, such a jump of E 0 leads to an unphysical spectrum because there is a range of E where ρ(E) < 0 (see the yellow curve in Fig. 6c). Such a negative spectrum has been observed also in JT gravity with deficit angles [14,15] and it indicates a breakdown of the gravitational computation. To be precise, the gravitational sum over loops of EOW branes, as a perturbative series in K, will not be convergent for K > K > cr because u(x), defined in (2.59), has a branch cut beginning at a positive value of x. 1 Indeed, when K increases past K > cr , the largest two real solutions of (2.78) collide and then split into two complex conjugate values, which are separated by the branch cut. Therefore, the rule to take 2φ b E 0 as largest real solution to (2.78) is not analytic.
From the dual matrix model point of view, this is interpreted as a nonperturbative instability [10,19] that requires a nonperturbative completion of the divergent sum to yield a reasonable spectral density. In Sections 3 and 4, we will consider a nonperturbatively well-defined matrix model whose tree-level spectral density agrees with the results in this section for K < K > cr but remains well-defined for K > K > cr . We will find that beyond this critical point, the spectrum undergoes a special phase transition to include both real and complex energies.
Before we move on to next section, we will point out a worse problem for the spectrum formula given by (2.77) and (2.78) when K > K > cr . Our regulated f λ (x) in (2.74) has singularities at x = −(1/2 + µ + n) 2 for all nonnegative integer n. If we take K large enough, there could be no real solution to (2.78) at all because f λ (x) blows up near its first singularity too fast (see the purple curve in Fig. 5c). This feature does not exist in JT gravity with deficit angles [14,15], where real solution of E 0 always exists. This is another evidence that gravitational computation breaks down when we include too many EOW branes if λ > 0.

Integrating out EOW branes
In this subsection we discuss the effective action for the metric and dilaton that remains after integrating out K species of EOW branes with mass µ. Having calculated the change of the spectral density due to the EOW branes, we now want to find a dilaton gravity theory whose disk path integral reproduces the total spectral density. Because the effective action associated with cusps is already known [14,15], we will focus on the theory where E 0 = 0, or where λ = 0 in (2.74). For small enough K, before any phase transition occurs, the spectral density follows from (2.79), (2.85) In [13], Witten provides a formula for the spectral density ρ(E; U ) associated with the dilaton potential W (Φ) = 2Φ + U (Φ), for a certain class of potentials. The result is This result was derived for a class of potentials given by For certain potentials, (2.86) predicts that the spectral density can be negative for some energies, similar to Fig. 6c. In this case, we expect that (2.86) breaks down. The authors of [19] studied an example where a naive application of (2.86) predicts negativity in the spectral density, while the correct spectral density has E 0 = 0. For now, we assume that Ke −S 0 is sufficiently small such that (2.88) then (2.86) coincides with (2.85). We claim that (2.88) is the correction to the JT action from integrating out the EOW branes. However, (2.88) is not manifestly of the form of (2.87). In the spirit of effective field theory, we will expand (2.88) in a power series in 1 µ . The first few orders are given by Each term in the above expansion takes the form of (2.87) for a particular limit of the parameters. Thus, (2.86) can be used to compute the spectral density associated with (2.88) order by order in the 1 µ expansion. Alternatively, we can use the result that the string equation associated to a potential of the form (2.87) is [13,15] ξ(x) 2π It follows that the string equation that corresponds to (2.89) is This agrees with (2.74) for λ = 0 to the same order in 1 µ . When µ is large and Ke −S 0 is order one, the EOW branes have a small effect on the spectral density, and one cannot reach a phase transition. As shown in Figure 5b, there is a phase transition when K > K = cr , which is µ-dependent. To investigate the validity of (2.88) across this phase transition, we note that for sufficiently large K, the largest zero of W (Φ), which we denote φ * , discontinuously jumps from zero to a positive value. As K increases further, φ * increases monotonically. The zero temperature entropy of a black hole, evaluated from the on-shell action, is given by 2πφ * . We can compare this semiclassical entropy with our exact calculations by examining the low-energy behavior of the spectral density, (2.77). For We define the effective zero-temperature entropy (or effective extremal entropy), S 0,eff , by the coefficient of The effective zero-temperature entropy, which we simply read off from the normalization of the edge of the spectrum, represents the change in the zero-temperature entropy due to integrating out EOW branes. For order one values of Ke −S 0 , S 0,eff and 2πφ * disagree; indeed, for K < K = cr , S 0,eff is nonzero while 2πφ * is zero. 2 However, for large values of K both S 0,eff and 2πφ * diverge as log K. This provides a qualitative check of (2.88) for K > K = cr . The fact that S 0,eff increases as log K for large K suggests that in this limit, the zero-temperature entropy is counting the number of species of EOW branes. Thus, we can view JT gravity with EOW branes as a model of induced gravity, since part of the black hole entropy in the effective dilaton gravity theory is induced by dynamical EOW branes that have been integrated out.
Indeed, such induced gravity interpretation is a general feature for heavy EOW branes in large K. For λ < 0 and large µ limit, we can use (2.78) and (2.93) to see the same scaling S 0,eff ∼ log K when Ke −S 0 1 even though it does not undergo any phase transition. For λ > 0 where the gravitational computation breaks down when K > K > cr , we will see the same scaling in Section 4.3 and 5.1 for K K > cr .

Effective matrix model
The duality between JT gravity and a double-scaled matrix model [10] can be expressed as the following identity where H is a N × N hermitian matrix. In this equation, Z n (β 1 , · · · , β n ) is the Euclidean path integral over all Riemann surfaces with n AdS boundaries with inverse temperatures β 1 to β n , and Tre −β 1 H · · · Tre −βnH is the expectation value of n operator insertions Tre −β i H in the matrix model in the double scaling limit. Z is the matrix model partition function with no insertion of operators. The double scaling is the special limit in which we zoom in the lower edge of spectrum in the large N limit while fixing the total number of eigenvalues ∼ e S 0 in any finite energy range. On both sides of (3.1), we have a topological expansion, in which each order is related to a genus g Riemann surface weighted by e (2−2g−n)S 0 . The beautiful work in [10] derives (3.1) by showing that both topological expansions are equivalent order by order. In this section, we will work out how dynamical EOW branes modify the matrix model potential V (H) and the genus zero spectral density.

Potential deformation by EOW branes
Our strategy is to find the "inverse trumpet"Z(β) such that 3 UsingZ(β), we can compute a path integral with n loops of EOW branes by integrating the path integral of pure JT gravity with at least n AdS boundaries against n inverse trumpets. Invoking the duality between pure JT gravity and a matrix integral, the Euclidean path integral over all surfaces with m AdS boundaries and n EOW brane loops is given by Equation (3.5) includes unwanted terms that may be interpreted as disconnected vacuum bubble geometries with EOW branes; these may be cancelled by multiplying Z by an overall constant. Finding the inverse trumpetZ(β) is straightforward. The inverse Laplace transformation of trumpet is given by (2.48), which written in energy variable is Using the orthogonality of cosine function This integral is not well-defined near b ∼ 0. We can regularize it by, for example, replacing 1/b with b α and taking α → −1 + in the last step. Using the identity we find for (2.49) Note that the leading divergence and constant log e 2µγ E 2π can be absorbed into the normalization of the matrix model partition function, and thus does not affect any correlation functions. Therefore, we find that summing over dynamical EOW brane loops in the path integral amounts to shifting the matrix model potential by KδV , with (3.12) From this matrix model analysis, it appears that we do not need to introduce any counterterm to the matrix model potential because the b → 0 divergence of (3.9) does not seem to affect any normalized correlators. However, just as the computations in Section 2 required regularization, the analysis of the double-scaled matrix model with the potential deformed by (3.12) requires additional input to be well-defined. In Section 3.3.3, we will explain why and interpret the meaning of the cusp counterterm introduced in section 2.3 from the perspective of the matrix model. The derivative of δV (E) is a sum of poles: Throughout this paper, we will assume that the support of the spectrum of the matrix model does not contain any of the singularities in (3.13). The deformation of potential (3.12) could also be rewritten as an integral over complex vector fields. Using identity Tr log A = log det A, we have where each Q i is a N -dimensional complex vector for i = 1, · · · , K. Alternatively to a deformation of the potential, we could understand each flavor of EOW loop as being dual to a coupling vectors in the matrix model. Each of the vectors Q i has the interpretation of a state with a single EOW brane in the gravity theory.

Saddle equation for the matrix model with deformed potential
In this section, we discuss how the shift in the potential δV affects the tree-level spectral density of the matrix model. Instead of using the explicit form of δV given in (3.12), we will use the more general formula (3.9) that applies for general geodesic length measures. This is because it will turn out to be necessary to use a regulated M(b) in (3.9) when determining the location of the lower endpoint of the spectrum.
The total potential of the matrix model is where V JT is the potential of the matrix model dual to pure JT gravity and δV is given in (3.9). Next, we define the effective potential 16) and the action where D ρ is the support of ρ(λ). As in [20], we call the second term in (3.16) the Coulomb gas replusive potential between pairs of eigenvalues. In the large N limit the saddle equation for variation of each λ ∈ D ρ gives V eff (E) = 0 for E ∈ D ρ , namely where the integral is understood as its principal value. Here ρ(λ) is normalized as Dρ dλρ(λ) = N . This means that the effective potential on the support of spectrum is a constant. Moreover, this constant is the minimum value of the effective potential along the path over which all eigenvalues are integrated. Physically, this means that all eigenvalues tend to stay in the lowest energy configuration which is a balance between the external force −V and the internal Coulomb repulsion for each eigenvalue. Let us define ρ(λ) = ρ JT (λ) + δρ(λ) for λ ∈ D ρ where ρ JT (λ) is the spectral density in the original SSS matrix model. Note that ρ(λ) and ρ JT (λ) may have different support and thus δρ(λ) could be a continuous but only piecewise differentiable function. Define the resolvent as One can show that R(E) is a double cover map from E ∈ C ∪ {∞} to C with a branch cut along D ρ . It follows from the definition and (3.18) that As the double cover map has a sign ambiguity along the cut, for a sensible matrix model solution, we need to assign the phase around the cut such that the measure ρ(E)dE is nonnegative along . Then it is easy to see that Applying this to SSS matrix model, we can write (3.18) as . This is a nice formula in which the LHS is known except the range of D ρ , and RHS has the same form as (3.18). We can define the LHS in terms of a new potential NṼ (E) and write (3.23) as This is the saddle-point equation for a single-matrix model with potential NṼ (E). The normalization condition for δρ is Standard techniques of solving for the resolvent and spectral density can be applied to the new potentialṼ . For example, we can define the variation of the resolvent and it follows that One can also use the Tricomi relation [20] to write δR in terms of contour integral of the potential where in the second step we used (3.27). In the double scaling limit, the rightmost end a 2s → +∞, and the formula reduces to

One-cut solution of spectral density
To simplify the computation of the tree-level spectral density of JT gravity coupled to EOW branes, let us assume that the lower edge of the spectrum E 0 is negative. This simplifies (3.23) as the third term in LHS vanishes due to D * = ∅. By (3.29), δR is linear in NṼ and thus is the sum of two terms that correspond to the first two terms of NṼ in (3.23). We call the first term the potential related piece (with subscript "K") and the second term the universal piece (with subscript "U ") respectively. For the potential related piece, we will consider a general measure for the EOW branes (2.52) giving an extra force For the unregulated measure (2.54), we recover (3.13). Because it is a convergent sum, using any regulated measure leads to the same result in → 0 limit. In the rest of this paper unless specified, we will use unregulated measure (2.54) in all integrals over α whenever it is convergent (and thus it makes no difference which regulated one is used, in the limit the regulator is removed). The minus sign in (3.30) shows that the force is leftward for 2φ b E > −(1/2 + µ) 2 and one might have conclude E 0 < 0. However, although δV (E) is independent of regularization and λ in the → 0 limit, E 0 still depends on the regularization in a subtle way that we will explain in Section 3.3.3.
Let us still assume here that E 0 < 0 for computational simplicity. In the end, the result (3.36) and (3.37) has an analytic form that can be easily continued to E 0 > 0 case. On the other hand, solving (3.23) with E 0 > 0 requires dealing with the third term in (3.23) as D * = ∅, but this leads to the same result.

Universal piece
As ρ JT (E) in (2.80) has a branch cut along the negative real axis, the universal piece of NṼ (λ) is proportional to sin(2π 2φ b (−λ)) and only supported for λ < 0. This implies that the contour integral of (3.29) can be reduced to surrounding Expand the sine function in Taylor series, in which each term can be evaluated by moving the contour to infinity. It turns out that the integral becomes a sum over residue at E and ∞, The first term has a branch cut on the positive real axis. Thus, along the negative real axis, δρ U (E) = 0. For E > 0, the contribution to δρ U (E) = −ρ JT (E) which cancels out ρ JT (E). Therefore the spectral density comes from the second integral around infinity. We can do a coordinate transformation λ → 1/z, and the integral becomes where in the second step we used Taylor expansion and picked out the coefficient of z n−1 . Here " " means ignoring the analytic part of R U (E) that does not contribute to spectral density. Let us denote the sum over n 0 as S n (E 0 /E).
The sum over n is straightforward As I 0 (x) is an entire function, the integrand is analytic for all E. Therefore, the discontinuity of δR U (E) at E ± is purely determined by the factor , we can easily see that (3.36) matches with the first term of gravity calculation (2.77).

Potential related piece
Using (3.30) and (3.29), we have where in second line we deformed the contour to poles at λ = E and λ = −α 2 /(2φ b ) and computed it using the residue theorem. This is justified by our assumption below (3.13) that 2φ b E 0 > −(1/2 + µ) 2 . It is clear that only the second term of (3.37) contributes to ρ K which exactly matches with (2.79).

Determining E 0
In the gravity computation of section 2.3, the ambiguity in the finite part of the cusp-like counterterm leads to a free parameter λ. In the following, we will show that this ambiguity naturally corresponds to an ambiguity in the zero point energy E 0 of the tree-level spectrum of the matrix model, which is sensitive to additional UV deformations to potential (3.12). Such UV deformations are not unique but play an equivalent role to the cusp-like counterterm in gravitational computation. From the string equation (2.76), we see that λ can be adjusted to set E 0 to a desired value. At first glance, this seems to be at odds with our understanding of matrix models. Ordinarily, the endpoints of a single-cut spectrum in a matrix model are not free parameters. They are fixed by the condition that the spectral density is normalized to the appropriate value. From (3.25) and the  Figure 7: The blue curve is spectral density ρ Λ (E) that matches with ρ JT (E) for E < Λ; the yellow curve is full spectral densityρ Λ (E) with E 0 < 0; the red curve is the difference δρ Λ (E). Due to normalization (3.43), δρ Λ (E) has a large negative "bump" above Λ.
simplifying assumption that E 0 ≤ 0, the normalization condition is As detailed in the previous subsection, δρ is a sum of two terms, δρ U = ρ U − ρ JT and δρ K . From (3.36), (3.38), and (2.54), we will have a different large E behavior of δρ U and δρ K , respectively, as 4 Thus, (3.39) cannot be satisfied for any choice of E 0 because δρ U and δρ K have different asymptotic behavior. The upshot is that the normalization condition can no longer determine the tree-level spectral density of the double-scaled matrix model, where the spectrum is assumed to have noncompact support on [E 0 , ∞). On the other hand, in a non-double-scaled matrix model, where the spectral density has compact support, the potential and normalization condition uniquely determine the spectrum. Thus, to determine E 0 , we must also specify how the change in the potential δV approaches (3.12) in the double scaling limit. We claim that for any E 0 , it is possible to add a term to the potential of the matrix model dual to JT gravity that converges pointwise to (3.12) in the double scaling limit such that the resulting spectral density is supported on [E 0 , ∞) and given by ρ JT + δρ, where δρ = δρ U + δρ K . There is one exception: when the spectral density corresponding to a given value of E 0 has negativity, that value of E 0 is ruled out.
To see this, consider a family of matrix models labelled by Λ whose tree-level spectra ρ Λ have compact support and converge pointwise to ρ JT in the double scaling limit Λ → ∞. For concreteness, we may assume that ρ Λ (E) agrees with ρ JT (E) for E < Λ. Next, define δρ Λ (E) to agree with δρ for E 0 < E < Λ. 5 For E > Λ, let δρ Λ be a smooth function that interpolates between δρ(Λ) and 0 in a way such that the normalization condition is obeyed, whereρ Λ ≡ ρ Λ + δρ Λ is the full spectral density. See Figure 7 for an illustration of δρ Λ . Similar to (3.23), the saddle point equation for δρ Λ is where δF Λ is the extra force that must be applied to the eigenvalues such that the resulting spectrum isρ Λ . Given δρ Λ , one simply uses the above equation to calculate δF Λ . Let a Λ+ be the location of the upper end of Dρ Λ . Then (3.41) becomes From the large E behavior of δρ(E), we see that the first term above diverges as Λ 1/2 log Λ. This controls how large the "bump" in Figure 7 can be. This justifies that normalization does not hold within the double-scaled regime. At fixed E, the second term of (3.42) goes to zero in the large Λ limit. We are left with We have thus demonstrated that it is possible to pick a δF Λ that converges pointwise to the derivative of (3.12) such that E 0 takes on any desired value (modulo the possibility that the resulting spectrum might have negativity). Thus, E 0 should be viewed as a free parameter, just as λ is a free parameter in (2.76). The above argument mirrors in the dual matrix model our treatment of JT gravity with EOW branes as an effective field theory. The fact that δF Λ (E) can be tuned for E > Λ to achieve different values of E 0 reflects the concept that the UV physics of the model can be tuned to achieve a desired IR theory. Of course, the exact shape of the bump in Figure 7 is unimportant; different choices represent irrelevant UV modifications that belong to the same IR universality class. As explained above, we have a family of double-scaled matrix models, each of which is labelled by a free IR parameter E 0 . The spectral density of each matrix model agrees with the result of a gravitational calculation, (2.77). To match a particular matrix model to a choice of the cusp counterterm λ in the gravity theory, we may use (2.78) to match the spectral density on both sides. This IR perspective neglects the UV details that determine these parameters. 6 In section 4, we will consider the regime λ > 0 and K > K > cr , where the gravity computation breaks down but the matrix model computation does not (once we specify a nonperturbative completion of the model, which is not unique). Thus, it is also worthwhile to consider a UV perspective, where we explicitly consider how the IR parameter E 0 is determined by UV data. The previous paragraph implies that this UV data is equivalent to the details of how the double-scaling limit is taken, and E 0 is determined by the normalization condition (3.40). To simplify computations, we have identified two alternative methods to determine E 0 from UV data that are easier to use in practice. First, instead of deforming the potential by (3.12), one can use (3.30) to define a δV associated with a regulated measure (such as (2.51)), and then impose the condition that δρ is a normalizable function. The UV data is contained in the details of the regulated measure. The second method is to take δV to be (3.12), but impose a condition on the asymptotic behavior of δρ that fixes E 0 (we will specify this below). The benefit of these alternative conditions is that they allow us to bypass the use of a non-double-scaled matrix model to define what we mean by the UV data of the matrix model. Also, they explicitly relate the UV data of the matrix model to the choice of cusp counterterm λ in the gravity theory. However, these alternative methods can only be fully justified by checking that they reproduce (2.78), which determines E 0 in the gravity computation. We explain these two alternative methods below in more detail. Let us first expand δρ(E) in large E. Note that ρ JT (E) corresponds to E 0 = 0 and K = 0 We have δρ(E) in large E limit (3.46) If we use a regulated measure, we may take the large E limit inside the integral of α and have where the α integral is just f λ (2φ b E 0 ) by (2.63), which is finite and involves λ from the regulated measure. The E −1/2 order integrates to E 1/2 divergence for large E that violates (3.39) within the double scaling limit. Imposing that δρ(E) is normalizable in the double scaling limit requires that the order E −1/2 term in (3.47) vanishes.
, the same zero point equation (2.78) follows after taking → 0 in last step.
Note that the final result for δρ(E) after the double-scaling limit is taken is normalizable but does not obey dEδρ(E) = 0. In general, the normalization could be nonzero (but finite in the double-scaling limit). Given zero point equation (2.78) before removing the regulator, one can show that It is clear that normalization depends on f λ (∞). From (2.53) and (2.63), we know that f λ (ξ) ∼ ξ −1/2 in large ξ. Normalization holds exactly only for specially designed f λ (x) such that On the other hand, if we use unregulated measure (2.54) in (3.46), taking the large E limit in the integral of α is illegal because it is divergent. As we mentioned before, this term has an To match with the gravity result, we can specify the UV data of dual matrix model as that for E → ∞ where γ E is Euler constant. This requirement leads to which fixes E 0 exactly as in gravitational answer in (2.78).

A "Y" shaped phase
As shown in Fig. 6c, the one-cut solution with λ > 0 is unphysical for K > K > cr as it has negative spectral density. On the gravity side, there is no obvious way to obtain a sensible result beyond the critical point. Fortunately, the dual matrix model can be non-perturbatively well-defined, and the method to compute spectral density for a given potential is known. It turns out to has a "Y" shaped cut in the complex plane as we will see shortly.

Matrix model on complex contours
Before we discuss the phase transition, let us review some basic features of the non-perturbative definition of matrix models. The partition function of a matrix model of N × N hermitian matrices is defined as where ∆(λ) = i<j (λ i − λ j ) is Vandermonde determinant and domain of integration of all eigenvalues λ i is R N . We can generalize this definition to N contours in the complex plane Γ = ⊗ N i=1 Γ i ∈ C N and The partition function of this matrix model is It is clear that the hermitian matrix ensemble is just the special case Γ i = R for all i. Absolute convergence for all eigenvalues is required to have a well-defined integral. With this condition, we restrict each integration contour Γ i to end at poles of V with an appropriate angle such that V (x) → +∞ for x → ∂Γ i , where ∂Γ i means the two ends of Γ i . Here we assume V (x) tends to its singularities faster than a logarithm so that the Vandermonde determinant (scaling as e 2N log x for each eigenvalue x → ∞) is always subleading and does not affect the integral's convergence. Though it seems that there are infinitely many ways to choose contours, contours that differ by smooth deformations do not change the integral due to analyticity. It follows that the space of independent contours is isomorphic to the homology space H 1 (e −V (x) dx). We can pick a basis γ i for this homology space, and for each eigenvalue λ i , its integration contour can be chosen as an integer coefficient linear combination of this basis where c ij is the number of times λ i is integrated along γ i and the sign defines the direction of integration.
For example, if V (x) = x 4 , its pole is at infinity and there are four directions approaching infinity with V (x) → +∞ (i.e. arg x ∈ (nπ/2 − π/8, nπ/2 + π/8) and |x| → ∞ for n = 0, 1, 2, 3), which leads to a basis of three independent integration contours, say γ n,n+1 for n = 0, 1, 2, denoting  Figure 9: (a) The effective potential of SSS potential for E < 0 where the first saddle is at First a few basis γ i (that are infinite many) on complex E plane for the matrix model dual to JT gravity. Here we use different colors to distinguish independent basis. Along each γ i to infinity V eff → +∞ (except being constant along positive real axis D ρ JT ) and along the shaded regions between contours V eff → −∞. To have the spectrum ρ JT (E) starting at origin, we must choose the defining contour Γ i for every eigenvalue λ i on one of the rightmost two contours that are labeled by bright red and green respectively.
contours connecting e inπ/2 ∞ to e i(n+1)π/2 ∞. There are only three basis contours because γ 0,1 + γ 1,2 + γ 2,3 + γ 3,0 = 0 by analyticity. In general, for a potential with a finite total degree of its poles, the dimension of H 1 (e −V (x) dx) is finite, while for potentials with infinite degree of poles (infinitely many poles or an essential singularity), the dimension of H 1 (e −V (x) dx) is countably infinite. For the N dimensional integration in Γ, to properly count the number of basis contours, we also need to quotient by permutations of the eigenvalues. This leads to a basis of dimension dim(H 1 (e −V (x) dx) ⊗N /V sym ), where V sym is the volume of the permutation symmetry group of the eigenvalues. Given a matrix model potential and assignment of contours Γ from this basis of H 1 (e −V (x) dx), for each eigenvalue λ i ∈ Γ i , one can use the saddle equation V eff (λ i ) = 0 (and thus setting V eff (λ i ) = 0) for λ i ∈ D ρ and V eff (λ i ) ≥ 0 for λ i ∈ Γ i \D ρ to solve for the tree-level spectral density ρ. Here we implicitly used the fact that saddle contour Γ i must have nonzero overlap with D ρ . If the spectrum is multi-cut, we also need to use the information of how much fraction of eigenvalues lying on each basis γ i to determine the spectral density on each cut.
Conversely, given a tree-level spectral density ρ, we can derive an effective potential using analytic continuation of ρ (see (3.22)). As V eff (x) has the same non-logarithmic singularities as V (x) by definition (3.16), we can use it to define a basis γ i as follows. For any saddle point m such that V eff (m) = 0 (including the ends of spectrum), we define an ascending path γ(m) emanating from m such that V eff (x) = V eff (m) and V eff (x) is strictly increasing for x away from m. Such a path will end on a singularity of V (x) or a point p on D ρ . In the latter case, we further extend γ(m) along D ρ (on which V eff is constant) to reach another end m of D ρ and then continue to a singularity of V (x) along an ascending path γ(m ). Define a descent-ascent path [20] as the union of two different ascending paths starting at the same saddle point m. See Fig. 8 for illustration. Our γ i is chosen as a complete and independent set from all descent-ascent paths. Each such contour γ i has a unique minimum value of V eff (x) and extends to singularities of V (x) monotonically.
Let us apply above construction of the basis γ i to the matrix model dual to JT gravity. The effective potential away from D ρ JT = R + is [10] which has infinitely many saddles at E = −n 2 /(8φ b ) for n ∈ Z (see Fig. 9a). Using the construction of descent-ascent paths, we can choose the basis γ i as plotted in Fig. 9b, where four nearby basis contours are separated around the local maximum of V JT eff (E) at E = −(2n − 1) 2 /(8φ b ) for n ∈ N. As there are infinite many n, the dimension of the basis is infinite.
Note that the interval around the saddles at −k 2 /(2φ b ) for k ∈ Z has lower effective potential than that on D ρ JT (see Fig. 9a). This means that the integral contour Γ for the spectrum of JT gravity cannot be along the real axis, otherwise there would be infinitely many eigenvalues filling these energy wells on R − . This is a known non-perturbative instability [10,19,21]. Indeed, in order to have ρ JT (E) supported only on the positive real axis, we must choose the defining contour Γ i for every eigenvalue λ i in (4.3) to be one of the rightmost two contours in Fig. 9b that are labeled by bright red and green respectively. Each of these two contours consists of two pieces (one along the real axis and the other extending into upper/lower half complex E plane) joined at the first saddle of V JT eff (E) at E = E saddle = −1/(8φ b ). Other choices of contours will lead to nonzero support of spectrum in some regions of R − . Given these two contours, there are still many ways to choose how many eigenvalues are assigned to each contour. Each choice defines a nonperturbative completion of the model. The most natural choice of Γ is to assign half of the eigenvalues to the upper contour and the other half to the lower contour, as this leads to a real partition function. This is the prescription we will use to study the phase transition.
Adding the extra potential from EOW branes changes the spectrum as seen in the gravitational computation. Increasing K in (3.15) from zero will also smoothly deform Γ as the basis defined as descent-ascent roads also smoothly moves. Let us analyze the effective potential near the critical point K > cr when λ > 0. Using the near edge spectral density (2.84), the effective potential for Before reaching the critical point, V eff (E) is positive in a small neighborhood of E 0 but it becomes negative at the critical point when a 1 = 0. For a 1 > 0 but close to zero, we can solve for the rightmost saddle of V eff (E) in the above approximation as At critical point, a 1 vanishes and E saddle coincides with E 0 . Note that increasing K does not change the direction of left-pointing extra force from δV and the two integration contours in Γ each has two pieces joined at E saddle . Given the fact that the end points of the tree-level spectrum of matrix Figure 10: The support of "Y" shaped spectrum. The two complex conjugate ends are E ± and the joint point is E 0 . γ ± are two arcs connecting E 0 to E ± respectively. models must deform continuously (not necessarily smoothly at phase transition point), 7 this means that past the critical point, the spectrum cannot extend along the real axis into the region with negative effective potential but rather it goes into the upper or lower complex E plane along Γ. Therefore, the support of spectrum after phase transition must be a "Y" shape with two complex conjugate pieces joined with a piece along real axis (see Fig. 10). 8 Such transition to "Y" shaped phases have been studied in other complex matrix models (for example in [23]). For readers not familiar with this topic, we examine a simple matrix model in Appendix C with cubic potential V (x) = x 3 /3 − tx 2 /2 and an even distribution of eigenvalues on complex conjugate contours C + and C − , where C ± is a contour connecting ±e 2iπ/3 ∞ to +∞. In Appendix C, we show that there are two critical values ±t cr = ±2 2/3 · 3 1/2 ≈ ±2.749 such that for t > t cr the spectrum is one-cut; for −t cr < t < t cr the spectrum is "Y" shaped; for t < −t cr the spectrum becomes to one-cut again.

"Y" shaped solution of the spectral density
Let us denote the two complex ends of the spectrum as E ± , the intersection on real axis as E 0 and the two arcs connecting E 0 to E ± as γ ± (see Fig. 10). As we assume eigenvalues are evenly distributed on the upper and lower contour, E ± and γ ± are complex conjugates. 9 Though eigenvalues can be complex, the spectral density ρ(E)dE must be real and nonnegative along γ ± . This is a strong condition that determines the shape of γ ± . Assuming E 0 < 0, it is easy to see that same equation (3.29) and the contour integral techniques of the one-cut case in Section 3.3 are both applicable to solve for δR(E). As usual, δR(E) is the sum of a universal term and a potential-dependent term.
The computation is similar to the one-cut case, and we relegate it to Appendix D. The result is as follows. We need to place the branch cut of δR(E) along γ For E ∈ γ ± , ρ(E) has the same formula as (4.8) and (4.9) but with nontrivial branch-cut choice of which is determined by the locus of γ ± by requiring ρ(E)dE to be real and nonnegative along γ ± .
We need three equations to determine three real parameters of the spectrum E 0 , E + and E + . By a similar argument to Section 3.3.3, there are two ways to impose appropriate UV data to determine these three IR parameters. One can either use a regulated measure and impose normalizability of δρ(E), or use the unregualted measure (2.54) and specify the same large E asymptotic behavior in (3.50). They are equivalent in the limit that the regulator is removed just as in the one-cut case. Let us take the large E behavior defined by (3.50) as an example. Indeed, it is very easy to see that (4.9) only has three non-normalizable pieces at large E, namely E 1/2 , E −1/2 log E and E −1/2 . Including the contribution from δρ U (E) at large E, the first piece should vanish and the remaining two pieces should obey (3.50). As computed in Appendix D, matching these three orders gives two equations It is a little surprising that three orders give us only two equations. This means that the E −1/2 log E piece is universal and independent of the phase transition in the IR. The last equation comes from requiring ρ(E)dE to be a real non-negative measure on γ ± . From the general expression for ρ(E) given in (4.8) and (4.9), we have where r and i subscripts represent real and imaginary parts respectively. Even though we do not know where the cut is a priori, the above decomposition holds up to an overall ± sign because ρ(E) is a two fold map. Being a real measure gives an ordinary differential equation which has a unique solution given an initial point E = E ± . On the other hand, we require γ ± to join E ± and E 0 , which gives another initial condition for this differential equation, leading to the final condition for determining the three parameters of the spectrum. On the other hand, since ρ(E) is analytic except at the branching point, the integral from E 0 to E ± is invariant under path deformation. Therefore, we can pick a simple one, for example the straight line from E 0 to E + , and require the integral to be real This condition is independent on the choice of branch cut as ρ(E) is a two fold map [20] which only leads to a sign ambiguity that is irrelevant in (4.16). The same applies to the path from E 0 to E − . This leads to the third equation required to determine the three parameters. In Appendix D, we rewrite it explicitly as

Phase transition to the "Y" shape
Before we solve for the "Y" shaped configuration numerically, let us first show that this "Y" shaped solution continuously connects to one-cut solution at critical point E 0 = E + = E − . The last condition for solving the ODE is trivial and we only need to check the first two conditions (4.11) and (4.12). Taking E 0 = E + = E − in (4.11), we have  Figure 11: (a) The numerical solution of E 0 (blue) and E r (yellow) as a function of E i . As E i increases, E 0 increases monotonically and E r decreases in the beginning then increases. At around E i ≈ 1.2, E r becomes larger than E 0 . (b) The numerical solution of ζ as a function of E i . As E i increases, ζ increases exponentially. In both pictures, we set 2φ b = 1.
To do numerics, we will take K ∼ O(e S 0 ), λ ∼ O(1) and µ 2 Ke −S 0 to simplify computation. This corresponds to a large number of heavy EOW branes. Let us first look at one-cut case. In large µ limit, we have Ke −S 0 f λ (x) ∼ Ke −S 0 µ −2 → 0 in (2.77). The only contribution to ρ(E) is from the universal piece. For the zero point equation (2.78), using (2.74) we see f λ (x) = λ is just a constant. In this limit, we can find the critical point for λ > 0 in Section 2.4 easily. Define ζ ≡ Ke −S 0 λ > 0. As we increase ζ from zero to the critical value, the solution of E 0 from (2.78) moves from zero to some negative E cr that is the largest solution to The simplification is similar in the "Y" shaped case. In (4.9) we have ρ K ∼ Kµ −4 ρ U ∼ O(e S 0 ) and thus can be neglected. The three conditions (4.11), (4.12) and (4.17) simply to 2+x−y ) 1/2 = 0 (4.23) One may naively think that each of (4.21) and (4.22) contain two real equations as E ± are complex variables. However, one can show that (4.21) is always real given the symmetry u → 1 − u and the integral in (4.22) is always imaginary using (4.21) (still under transformation u → 1−u). Therefore, they only give two real equations.
Our numerical strategy is to solve (4.21) and (4.23) for a given E i and use E i to compute ζ using (4.22). This reduces the numerical work to solving only two integral equations. The result is plotted in Fig. 11a and 11b, where we set 2φ b = 1 for simplicity. At E i = 0, both E 0 and E r start at E cr . As E i increases, E 0 increases monotonically and E r first decreases, then increases. Before reaching E i ≈ 1.2, E r is less than E 0 , and it becomes larger for E i 1.2. While E 0,r seems to increase in roughly the same order as E i in the numerical range, ζ increases exponentially with E i . It is interesting that as E i increases (corresponding to increasing ζ), both E 0 and E r become positive. As we will see later, this implies that the canonical ensemble partition function for fixed temperature will decrease for growing ζ.
Given the solution for E 0,r,i , we can find the support of the density γ ± by solving ODE (4.16). We illustrate the result for two parameter choices in Fig. 12a and 12b. In these two pictures, we plot the vector field (ρ r (E), −ρ i (E)) on the upper half complex E plane. The red curve connecting two black dots is the support of the eigenvalues γ + on which ρ(E)dE is nonnegative and real. γ − is the reflected curve with E → − E. In Fig. 12a, E i is small and the whole γ + has negative real part, for which particularly E r < E cr and E 0 > E cr ; in Fig. 12b, E i is large and the whole γ + has positive real part. This is consistent with the behavior in Fig. 11a. As E 0 is a triple branch point, near E 0 we have ρ(E)dE ∼ d(E − E 0 ) 3/2 , which implies that the angle between any two branches is 2π/3. To have an intuitive impression of ρ(λ)dλ on γ ± , using dλ ∝ ρ r (λ) − iρ i (λ) for λ ∈ γ ± , we have where λ(λ i ) is the locus of γ ± with affine parameter λ i . Note that we need to choose the right sign of ± such that ρ(λ)dλ is nonnegative. In Fig. 13a, we plot |ρ(λ(λ i ))| 2 /|ρ i (λ(λ i ))| along γ + for E i = 1. The total number of eigenvalues on γ + (which is the same as that on γ − ) is computed in (D.21) by integration of ρ(E)dE along γ + (or equivalently the straight line) from E 0 to E + . As a function of ζ, we plot it in Fig. 13b. From the plot, it seems that the number of eigenvalues on γ + quickly becomes a power law growth in ζ (note that Fig. 13b is log-log plot). We now turn to the interpretation of the complex energy spectrum resulting from the matrix model analysis in terms of black hole physics. In next section, we will use the effective W (φ) dilaton JT gravity to illustrate a similar "Y" shaped phase transition using semiclassical analysis of Euclidean black holes. Based on that observation, we will interpret the complex energy states in spectrum as unstable black holes that decay to eigenbranes, which are not included in current matrix model. See section 5.3 for more details. The bottom line is that we can treat the real energy part of spectrum as stable black hole states. For this piece of the spectrum, we would like to define the lowest energy black hole entropy (similar to zero-temperature entropy in Section 2.5) as This S lowest could also be understood as the effective extremal entropy S 0,eff of stable black holes due to integrating out EOW branes along the lines of the discussion in Section 2.5, although the full complex spectrum does not start at E 0 . In the large µ limit, the leading contribution to ρ(E) is the universal piece ρ U (E). Using (4.8), we have We plot exp(S lowest ) as a function of ζ in Fig. 14. From this plot, it is clear that exp(S lowest ) scales linearly in ζ for large ζ, which is relevant for K e S 0 . In other words, numerics suggest the lowest energy black hole entropy scales as As mentioned in Section 2.5, it is very natural to interpret this as a state counting result, where K flavors of EOW branes behind the horizon are microstates of the lowest energy black holes.

Gas of cusps for heavy EOW branes
In Section 2.5, we derived an effective W (Φ) dilaton gravity for the λ = 0 case by comparing spectral densities. In this section, we will study the effective W (Φ) related to heavy EOW branes with λ = 0.  The large µ limit basically receives only contributions from small b of order 1/µ in the EOW brane measure M(b). On the other hand, the b → 0 limit shrinks the EOW loops into cusps as discussed in Section 2.2. Therefore, adding heavy EOW branes to JT gravity can be understood as deformation by cusps, the 2π deficit angle, which is described by a new potential for the dilaton Φ [14,15] W (Φ) = 2Φ + 2χe −2πΦ (5.1) with total Euclidean action On the other hand, large µ leads to f λ (x) = λ. Comparing (2.78) with (8.4) in [14], we can determine Given a generic dilaton potential W (Φ) with asymptotic behavior W (Φ) ∼ 2Φ for Φ → +∞, the classical solutions of the equations of motion were studied in [24] and recently in [13]. Choosing the diffeomorphism gauge Φ = φ h + r, the solution for the 2d metric is where φ h is the horizon value of the dilaton and r ∈ [0, +∞). Matching with the AdS boundary condition (2.13) at r = r ∞ , we have The ADM energy of such a black hole is Taking r ∞ → ∞ limit, and using (5.1) and (5.3), we have where ζ = Ke −S 0 λ as defined in Section 4.3. The black hole temperature is determined by requiring smoothness of the metric at r = φ h . Setting y = √ r, and A(r) = W (φ h )r + o(r), we have for y 1 For t ∈ R, this metric has a conical singularity at y = 0 unless In the following, we will mainly consider λ > 0 (namely ζ > 0). For small ζ, W (Φ) = 0 has two real roots (see Fig. 15a). Based on [13], we should take the larger root as the lower bound for φ h because we require A(r) > 0 for all r > 0. The red part of Fig. 15a shows the allowed range of φ h , which by (5.7) determines the spectrum. In particular, the black hole temperature can take any nonnegative value. There is a critical value of ζ cr where these two roots coincide at φ cr (see Fig. 15b) This gives a critical zero point energy If we compare with the exact result (4.20), we see they are off to some extent but not too much.
As this is only a semiclassical analysis, we should not expect an exact match with the full quantum computation. They match qualitatively, especially in that for some positive ζ we get negative critical energy E cr . If we increase ζ past the critical value, we will have all φ h ∈ R allowed in the spectrum by the condition A(r) > 0 for r > 0. In this case, as W (φ) is lower bounded by a positive value, the minimum temperature of black holes is nonzero and given by where φ h = φ * = (2π) −1 log 2πζ with W (φ * ) = 0. For a given temperature T > T * , there are two solutions for black holes with different φ h = φ h,± . It has been shown in [13] that the free energy The support γ * of energy spectrum on complex E plane. It has two branches γ * ,± conjugate to each other.
difference between these two solutions is which implies that the larger black hole with φ h,+ is thermodynamically stable. Equivalently, one can also show that the specific heat dE/dT is negative for the smaller black hole with φ h,− . Note that we have set κ = (8πG) −1 equal to 1 throughout this paper. Restoring the G N dependence, the contribution to the partition function from smaller black holes is relatively O(e −1/G N ) suppressed in the small G N limit. Therefore, if we only look at the spectrum of canonically stable black holes, the energy will be lower bounded by This physically corresponds to the lower bound E 0 of the real energy spectrum in the "Y" shaped phase of the matrix model computation.
The entropy of this lowest energy stable black hole is S * = 2πφ * = log 2πζ (5.15) which scales logarithmically with K. This also matches with the numerics for S lowest in the large K limit as shown in Fig. 14.
Similarly, one can easily show that for λ < 0, the zero point energy is given by W (φ * ) = 0, which leads to same logarithmic scaling S * ∼ log K for Ke −S 0 1. This, at the semiclassical level, justifies the induced gravity interpretation mentioned at the end of Section 2.5.
One can continue the condition W (φ h ) = 0 that defined the endpoint of the spectrum for ζ < ζ cr past the critical value, where φ h becomes complex. This corresponds to the complex part of spectrum in the matrix model. We could consider complex saddles of Euclidean action simply  Figure 17: (E * ,+ ), E * and log 10 ζ as a function of (E * ,+ ). We set 2φ b = 1 in both plots.
by analytic continuation of φ h to complex numbers in (5.4). As φ h is a free parameter, it is obvious that such a complex metric and dilaton solve the equations of motion (see Fig. 16a). For this complex solution to be a saddle, we still need to require the complex metric is smooth around φ h . Expanding in small r, (5.8) still holds. For t ∈ R, the metric has conical singularity at y = 0 unless W (φ h ) ∈ R and t have the same periodicity as (5.9). The condition W (φ h ) ∈ R gives the allowed complex φ h , which defines the support of the complex part of the spectrum γ * via E(φ h ) (see Fig.  16b). The endpoint of γ * corresponds to W (φ h ) = 0 as it is the singularity of metric (5.8).
Let us check what follows from this complex saddle. For a given ζ, the solution of W (φ h ) = w is given by where W(z) is product logrithmic function (or Lambert W function) defined implicitly as We W = z. This is a multivalued function but we will only choose branches that connect to real φ h for ζ < ζ cr . This leads two solutions φ h,± (w) conjugate to each other and thus splits γ * into two pieces γ * ,± joint at E * , similar to γ ± in matrix model (see Fig. 16b). Setting w = 0 and using in (5.7), we can solve the two ends of γ * ,± To have a comparison with matrix model result, we plot (E * ,+ ), E * and ζ as a function of (E * ,+ ) in Fig. 17. Qualitatively, this matches with the exact result in Fig. 11.

Dirichlet-Neumann EOW branes
We found that for ζ > ζ cr , no real solutions for W (φ h ) = 0 exist. This means that there is no zero temperature black hole. It is natural to ask what happens if we decrease the temperature below T * . This phenomenon is not unique to JT gravity. There is also a minimum temperature T min for d > 3 dimensional AdS-Schwarzschild black holes. For any fixed temperature above T min , two black holes solutions exist and only the bigger one is thermodynamically stable.
Hawking and Page studied this property in their well-known paper [25] and found a first order phase transition from black holes to a thermal gas at some higher temperature T HP > T min . At T = T HP , if we decrease the total energy, black holes will be in equilibrium with a thermal gas in the microcanonical ensemble. For T < T HP , one transitions completely to the thermal gas phase.
By analogy with higher dimensional black holes, we consider how a "thermal gas" could exist in W (Φ) dilaton gravity for temperatures below a Hawking-Page temperature T HP > T * . The simplest example is to consider another type of EOW brane that leads to spacetimes different from black holes. All types of EOW branes in pure JT gravity are classified in [12] based on their boundary conditions. We can perform a similar classification for W (Φ) dilaton gravity. The Dirichlet-Neumann (DN) EOW brane (in the language of [12]), corresponding to a location where the dilaton reaches −Φ 0 , is the dimensional reduction of the smooth origin of polar coordinates in the center of empty AdS. Thus it is the natural candidate for states that behave like the thermal gas.
Let us consider the Euclidean action where d is a constant defined on the DN brane world line ∂. The variation of action gives equations of motion as well as boundary conditions (see Appendix A), which on DN brane are We fix Φ = φ DN (Dirichlet) and also n a ∂ a Φ + d = 0 (Neumann). The bulk solution has the same form as (5.4) but with a different "horizon" parameter φh. Note that this spacetime does not have any real horizon at φh because we truncate the spacetime at a new boundary (DN brane) before it, which restricts φ DN > φh.
Using solution (5.4), we have The Hamiltonian defined in (2.30) only depends on the asymptotic behavior of Φ and thus is given by the same formula as (5.6) with φ h → φh. 10 As d is a fixed constant, (5.20) determines φh and completely fixes the Hamiltonian to be This corresponds to the eigenbrane first introduced in [26] that has fixed energy. In the dual matrix model, it gives delta function in spectrum [12].
To discuss dominance for a given temperature, we will evalute the Euclidean on-shell action for DN brane and black hole respectively. In both cases, Euclidean time is periodic 10 No Hamiltonian arises on the DN brane due to the boundary conditions.
For the DN brane, extrinsic curvature is The action on the DN brane is Similarly, the extrinsic curvature on the AdS boundary is The action on the AdS boundary is The bulk action is Altogether we find that the Euclidean on-shell action for DN branes is which is as expected because the DN brane is an energy eigen state and thus has zero entropy S = βE − I = 0. Note that there is no topological term linear in S 0 in I DN because the topology is a cylinder. On the other hand, the Euclidean action for black holes only has two terms, I AdS and I bulk , because the geometry is smooth at the horizon φ = φ h . For black holes, φ h is not a free parameter and is determined by β using (5.9). We have As the topology is a disk, we do have the topological term −S 0 = −2πφ 0 . This leads to where we used (5.9) to write I bh in the form of βE − S.  Figure 18: (a) Plot of I bh + 2πφ 0 as a function of T for different ζ, where blue is for ζ < ζ cr , yellow is for ζ = ζ cr , and green is for ζ > ζ cr ; (b) and (c) are both plots of I bh and I DN as a function of T , where blue curve is for I bh , yellow curve is for I DN with E DN < E DN,max and green is for I DN with E DN > E DN,max . We set 2φ b = 1 and S 0 = 10 in these plots.

Hawking-Page phase transition
The Hawking-Page phase transition is determined by dominance of e −I DN versus e −I bh . A quick observation is that for any temperature, for large enough d, the DN brane is always dominant. This corresponds to very low energy eigenbranes.
There is no principle that determines d a priori, thus the Hawking-Page phase transition is model dependent. However, the upper bound of eigenbrane energies is for d = 0 and we can compute a lower bound for the phase transition temperature corresponding to this case.
For W (Φ) potential (5.1), we have where φ h (β) is found by solving (5.9) The plot of I bh + 2πφ 0 as a function of T is shown in Fig. 18a for a few different ζ. For ζ ≤ ζ cr , I bh + 2πφ 0 has a maximum at finite temperature and tends to negative infinity for both zero and infinite temperature. For ζ > ζ cr , I bh + 2πφ 0 still tends to negative infinity for infinite temperature but ends at a finite value at T * . The temperature dependence of I DN is very simple, namely inversely proportional to T . For large T , it is clear that I DN > I bh and black holes dominate at high temperature. For small T , in order to have I DN < I bh , we must require φ 2 DN − ζ π e −2πφ DN < φ h (β) 2 − ζ π e −2πφ h (β) for β = 1/T * (ζ > ζ cr ) or β → +∞ (ζ ≤ ζ cr ). Taking the T derivative of I DN − I bh , we have because φ h (β) is a monotonically increasing function of T in the allowed range of temperature. This means that there is only one Hawking-Page phase transition point if it exists. In other words, a necessary condition for a Hawking-Page phase transition is that the energy of DN brane must be below the stable black hole spectrum. For ζ ≤ ζ cr , the lowest energy of a stable black hole is E * ,± in (5.17) (it takes real values in this case). For any E DN below it, the Hawking-Page phase transition always exists because is an increasing but bounded function of β. In other words, the upper bound of E DN is E DN,max = E * ,± for ζ ≤ ζ cr . For ζ > ζ cr , the upper bound E DN,max depends on the value of φ 0 because Hawking-Page phase transition must occur above T * . To find it, we need to solve I DN = I bh at β = 1/T * . It turns out that We plot with a few examples in Fig. 18 to show that Hawking-Page phase transition occurs only when E DN < E DN,max . The lesson of above discussion is that including eigenbranes with low enough energy in the black hole system is necessary to have a canonical ensemble description at all temperatures. If we isolate the black holes, zero temperature (defined by the periodicity of Euclidean time) can only be achieved by complex geometries.Therefore, it is natural to interpret these complex black holes as metastable states that decay to eigenbranes in finite time.
On the matrix model side, adding eigenbranes corresponds to fixing some eigenvalues λ i by inserting i δ(λ i − λ DN,i ) with a set of fixed numbers λ DN,i in the matrix integral [26]. If we take λ DN,i to be real numbers to the left of the support of spectrum, they will exert a rightpointing repulsive force on all other eigenvalues, which can balance the left-point force exerted by the potential of EOW branes in (5.3). It would be interesting to see how this would change the solution of the matrix model. In particular, we expect to see no complex part of spectrum if we assign eigenbranes properly. This question is beyond the scope of this paper and we will leave it for future works.

Discussion and conclusion
Matrix model dual to JT gravity with deficit angles The matrix model dual to JT gravity with deficit angles has been studied in [14,15]. For a single type of deficit angle 2π − θ (0 ≤ θ < π) with weight ε, the genus zero one-cut spectral density is where the zero-point energy E 0 is determined by the largest real solution to If θ = 0, the second term in brackets of (6.1) scales like ∼ e θ √ ξ ξ −3/4 in the large ξ limit. This term leads to a spectral density that differs by an exponential amount at large E from that of pure JT gravity, ρ JT (E), in contrast to the behavior of our matrix model dual to EOW branes. In matrix model language, this corresponds to a huge left-pointing extra force δF Λ that pushes a large number of eigenvalues into the double-scaled IR regime. In the following, we will study this matrix model by deriving the corresponding potential deformation δV . For JT gravity with deficit angles, the measure is formally an analytic continuation of a delta function, δ(b − iθ). The way to define this analytic continuation in the matrix model is as follows. First, we assume the measure M(b) to be δ(b − b 0 ) for b 0 > 0, which leads to where b 0 > 0 and infinitesimal δ > 0. 11 For any rational function F (α) that is integrated against m(α), adding a lower half plane infinite arc for k does not change the result because b 0 > 0. Therefore, we can define the analytic continuation of b 0 based on a slightly different m(α) with this extra lower half plane infinite arc where C − ∞ is the clockwise contour consisting of the real axis and the lower half plane infinite arc. Using this definition, the integral of F (α) against m(α) becomes contour integrals of k around the singularities of F (δ − ik) and thus b 0 can be safely continued to complex numbers. Taking this definition in (3.30), we have For generic complex b 0 , we will have a complex potential along the real axis. Putting (6.4) into the first line of (3.37), where we need to assign the contour of λ at infinity inside of the contour of −α 2 /(2φ b ), 12 (3.38) follows and using the first line of (2.79), we have (assuming ξ > 0) For ξ < 0, we simply deform the contour of k in the first line appropriately such that the second 11 Here the complex delta function is a formal device that just means α is taken to a specific value when integrated over D. A mathematically rigorous formulation of δ(α − α0) could be line holds. It is clear that (6.6) matches with the second term in (6.1) after continuing b 0 → iθ. From (6.5), it is also clear that the potential exponentially grows for large E for b 0 → iθ, and this huge potential deformation explains the deviation of ρ from ρ JT as mentioned below (6.2).
As discussed in Section 3.3.3, to find the equation determining the zero point energy E 0 in the matrix model, we need to input UV data of either a regulated measure or asymptotic large E behavior of δρ(E). While it is unclear how to define a regulated measure for a deficit angle, figuring out required asymptotic large E behavior is straightforward. Indeed, we can rewrite (6.6) as Comparing with the first term in (3.46), it obvious that we need to require that for E → +∞ Continuing b 0 → iθ, it matches with (6.2). Note that for generic complex b 0 , the zero point energy E 0 can take multiple complex values because (6.9) has multiple solutions. We will pick the one smoothly connecting to the E 0 of b = iθ case, which is the largest real solution to (6.2). Given a complex E 0 , the support of the spectrum D ρ is defined by requiring ρ(E)dE to be a real nonnegative measure along D ρ asymptotically extending to infinity. This ODE completely fixes the one-cut spectrum of the matrix models. Indeed, the above construction defines a class of matrix models labelled by the complex number b 0 , which contains two special cases with geometric duals in deformed JT gravity: with deficit angles (b 0 = iθ) and with fixed length EOW branes (b 0 > 0). At first glance, specifying asymptotic behavior (6.8) looks like a fine-tuning starting with the SSS matrix model because a deviation at any E p (p > −1/2) order will strongly violate the zero point equation although it seems negligible compared to the exponentially large behavior at large E for b 0 = 0. However, the SSS matrix model is equally fine-tuned in this sense, since its spectral density also grows exponentially at large E in double scaling limit.
The method used in [14,15] to solve JT gravity with deficit angles by summing over topologies only applies to one-cut case. Given the dual matrix model defined above, it would be straightforward to study the spectrum beyond the critical point in this model requiring the same asymptotic behavior (6.8) of the spectral density. With this definition, we expect to see a phase transition similar to our "Y" shape because heavy EOW branes with λ > 0 can be regarded as the special case θ = 0 for deficit angles as discussed in Section 5.1. Our method should thus offer a reasonable solution to the negative spectrum (and non-perturbative instability) problem in [14,15,19].

Comparison with previous work
Recently, a different non-pertubative completion for JT gravity and JT supergravity has been discussed [21,[27][28][29] that in some cases avoids non-pertubative instability. In their method, the Hermitian matrices are promoted to squares of complex matrices and the spectrum is nonperturbatively bounded from below. Alternatively, one could use the original Hermitian matrix but introduce a hard wall in the potential. This completion was further generalized to JT gravity with deficit angles [19].
There, all eigenvalues are still integrated along the real axis but the matrix model potential is modified (in particular, by fixing a specific asymptotic behavior for the function u(x) that appears in the string equation as x → ±∞). The computational method uses an auxiliary Schrödinger equation that is associated with the orthogonal polynomial formalism of double-scaled matrix models with even potentials [30]. The standard orthogonal polynomial technique only applies to the case where all eigenvalues are on the same integration contour. Generalizing to the case of multiple distinct contours may require multi-matrix techniques, in which the coupling among eigenvalues on different contours from the Vandermonde could be treated as interactions among different matrices.
Our method to nonperturbatively define a double-scaled matrix model starts with the observation that in a non-double-scaled matrix model with analytic potential, the tree-level spectral density plus a choice of contours for the eigenvalues completely specify the model. The doublescaled matrix models in this paper are defined such that this property still holds. 13 We promote the Hermitian matrices to normal matrices along some contour for each eigenvalue, and the effective potential is defined everywhere by analytically continuing the spectral curve. The nonperturbative definition is thus supplied by the spectral curve and the choice of eigenvalue contours. This is along the same lines of [10].
More importantly, one of our results is this completion naturally applies to the multi-cut case and avoids the non-perturbative instability of [19]. Our completion allows order e S 0 eigenvalues to be complex in the "Y" shaped phase, and it seems impossible to interpret each matrix in the ensemble as an instantiation of a Hermitian Hamiltonian. As argued in Section 5, we can interpret the complex eigenvalues as metastable black hole states that decay in finite time. Indeed, this implies that JT gravity with EOW branes is not a complete theory and one must include other objects with lower energy, e.g. DN branes, to have an unitary theory with Hermitian Hamiltonian.
On the other hand, the formalism of [21,[27][28][29] can easily compute higher genus contributions. We only considered the tree-level (genus zero) spectral density that solves the saddle equation for the matrix model potential. It would be interesting to study higher genus contribution to our matrix model in future works. Phase space of W (Φ) dilaton gravity Given that heavy EOW branes with λ > 0 lead to an effective W (Φ) dilaton gravity, it is also interesting to study the phase space for this theory along the lines of Section 2 and [11]. Let us consider the case of one-sided black holes in W (Φ) dilaton JT gravity with one AdS boundary and one µ = 0 EOW brane boundary, for simplicity. This is equivalent to half of a single sided black hole. The Hamiltonian is given by (5.6). Using the solution from [24], one can show that the regularized geodesic length connecting the AdS boundary at u = 0 to the EOW brane normally is given by where φ 0 is the dilaton value at the end of geodesic on the EOW brane, which depends on the location of the other end at the AdS boundary in a complicated way. After some algebra, one can show that the Hamiltonian can be written in a canonical form in terms of L and its conjugate momentum P as where φ 0 (L) is defined implicitly by (6.10). One can check that for W (r) = 2r, it reduces to the pure JT result (2.37) with µ = 0. For W (Φ) given by (5.1), we plot the potential in Fig. 19a for different ζ. It turns out that for 0 < ζ < ζ cr , the potential is similar to Liouville potential (ζ = 0) but with a negative shift of the ground energy. For ζ > ζ cr , the potential is double valued for some L and thus ill-defined. Also, the potential is unbounded from below and the allowed range of L is finite. This matches with the expectation from Euclidean computations (5.7) that the spectrum is unbounded from below for ζ > ζ cr because φ h can take any real value.
In the Lorentzian picture, this ill-defined potential reflects the fact that L is no long an appropriate phase space variable because of the ambiguity of two geodesics emanating from same AdS boundary point but ending on two different points on the µ = 0 EOW brane normally (see Fig.  19b). In other words, using L as a radial coordinate is no longer a good gauge for the metric. It would be interesting to find an appropriate gauge in this case and see how the quantization leads to results consistent with the semiclassical analysis of Section 5.

Modified inner product from geodesic separable objects
As formulated in Section 2, the Hilbert space for gravity in spacetimes with one AdS boundary and one EOW brane can be found by quantizing in the regularized geodesic length |L basis. Canonical quantization requires orthogonality of this basis L 2 |L 1 = δ(L 1 − L 2 ). Such Lorentzian analysis is based on a fixed spacetime topology whereas full quantum gravity allows for contributions from all topologies. To take this into account, we can define a modified inner product between |L states using Euclidean path integrals including handles and loops of EOW branes between the two geodesic slices as shown in Fig. 20a.
Similarly to the gluing of partition function in Fig. 3, the idea is to find the geometric building block Ψ 3 in Fig. 20e that contains a geodesic loop with length b between two geodesics L 1 and L 2 that connect the AdS boundary with the EOW brane. The contribution of other topologies to modified inner product L 2 |L 1 is obtained from Ψ 3 by gluing with all possible geometric objects (handles and EOW brane loops) along the geodesic loop of length b. In other words, Ψ 3 plays a role similar to the trumpet in the computation of partition functions.
To evaluate Ψ 3 , we need some gluing and pinch-off surgery of the Euclidean propagator G µ,β (L 2 , L 1 ) in (2.43). Diagrammatically in Fig. 20c, G µ,β (L 2 , L 1 ) is the same as the Euclidean path integral over the region bounded by an AdS boundary arc with length β (black curve), two geodesics L 1 and L 2 (orange curves) and an EOW brane (blue curve). We will first pinch off the region (Ψ 1 in Fig. 20b) bounded by the AdS boundary arc and the geodesic connecting its two ends (green curve). Indeed, Ψ 1 is the Hartle-Hawking state in the two-sided AdS system in JT gravity [31] where we defined ρ 0 (k) ≡ 2k sinh πk π , y ≡ 4e − /2 (6.13) Using normalization of Bessel function one can show that gluing two Hartle-Hawking states along geodesic leads to disk partition function Pinching off Ψ 1 from G µ,β leads to a piece Ψ µ ( ; L 1 , L 2 ) bounded by four geodesics as labelled in Fig. 20c. It can be written as where z ≡ 4e −L and ϕ(k) is given in (2.42), because one can check easily that Then, we will glue two Ψ µ along the geodesic to get Ψ 2 in Fig. 20d Lastly, we glue L 3 and L 4 of Ψ 2 , which forms an EOW brane loop where Ψ 3 is depicted in Fig. 20e, and pinch off the EOW brane measure M(b) on that loop to get where Ψ k (z) is the normalized eigen function from (2.42), and we used the technique from (2.43) to (2.47) to separate M(b) from RHS of (6.19). It follows that the modified inner product in Fig. 20a is where X(b) represents all topologies of handles and EOW brane loops that are separated by the geodesic loop of length b. Indeed, (6.21) can be applied to any geometric objects that are separable by a geodesic loop as they are all accountable by some X(b). We can rewrite cos kb as the inverse Laplace transformation of the trumpet using (2.48). Note that the partition function can be written in a similar form as It follows that where in the second line we assumed that the b integral is interchangeable with the other two, and in the last line ρ(E) is the full spectral density for E = k 2 /2φ b and inverse Laplace transformed from the partition function Z(β). If ρ = ρ JT , we get back to L 2 |L 1 = δ(L 1 − L 2 ) as expected. This result shows that the modified inner product among |L basis can still be diagonalized in terms of the positive energy basis of the Morse potential quantum mechanics even after taking all nontrivial topologies into account. In particular, the contribution from each energy is proportional to its spectral density.
On the other hand, if the full spectral density allows negative energy, like the tree-level spectrum with λ > 0 in this paper, the states |L are blind to it. Indeed, if we naively extend Ψ k (z) to imaginary k, this wave function becomes non-normalizable.
Importantly, if the full spectral density is discrete [26], |L becomes overcomplete and there are infinitely many null states, namely the energy eigenstates that are not on the support of the spectrum. It would be interesting to understand what geometric objects are required to have such discrete spectra and how they could lead to a unique α state of baby universes in such a modified JT gravity, in the language similar to pure topological 2D gravity [32]. The variation of the boundary action becomes δg n a − (∇ a n a h bc − ∇ b n c + n a ∇ a n c n b )δg bc Φˆ Using the fact that where∇ a is the covariant derivative constructed from the induced metric on ∂M , we can write the first term in (A.7) as The third term in δI ∂M is nontrivial in general dimensions but vanishes in 2d. This can be easily seen using h ab = −t a t b . Putting the bulk and boundary variations together, we have (∇ a n a − C)δΦˆ (A.10) The boundary condition for a well defined phase space must be one of the following two types    n c ∂ c Φ = CΦ + D or h ab δg ab = 0 ∇ a n a = C or δΦ = 0 (A.11)

B Lagrange inversion theorem
Consider two smooth functions g(x) and f (x) that satisfy f (0) = g(0) = 0 and f (g(x)) = x. The Lagrange inversion theorem is a formula for the Taylor series coefficients of g(x) in terms of the Taylor series coefficients of its inverse, f (x). The formula is given by (B.1) This formula can be slightly generalized to compute h(g(x)) for any smooth function h defined in a neighborhood of 0. The result is We could choose f (x) = x/φ(x) for φ(0) = 0, which leads to Let f λ (x), g λ (x) be a one-parameter family of pairs of functions obeying the same properties as f (x) and g(x) above. Let f λ (x) = x/φ(x + λ). We have We may define g λ implicitly in terms of φ as follows: We now write where g λ (y) = Φ(g λ (y) + λ, y), Φ(λ, y) = yφ(λ).
Consider a family of variables y i and functions φ i (x), indexed by i. We can generalize the above as follows: We have not discussed the radius of convergence of the above series. In practice, we need to sum the series within its radius of convergence and then analytically continue y as needed.

C Phase transition in matrix model with cubic potential
Let us consider the matrix model with cubic potential V (x) = x 3 /3 − tx 2 /2. For x → ∞, there are three directions where V (x) → +∞, namely x → +∞ and x → ±e 2iπ/3 ∞. As discussed in Section 4.1, there are three contour choices for each eigenvalue. We will assign half of the eigenvalues to the contour C + and the other half to the contour C − , where C ± is defined to be the contour that connects ±e 2iπ/3 ∞ to +∞ (see Fig. 21). In this appendix, we use a slightly different convention and divide both ρ(E) and R(E) in (3.21) and (3.22) by N such that the density of states is normalized to unity. For large |t|, we see that the potential has a deep well along the real axis (see Fig. 22), and we expect the eigenvalues to stabilize in the well and form a one-cut solution along the real axis. Because we have evenly assigned the eigenvalues to the two contours C ± , the spectrum must be invariant under a reflection across the real axis. Aside from a single cut along the real axis, another possibility that respects this symmetry is a single cut along a complex curve that is invariant under complex conjugation (see Fig. 21 (a) and (b)). However, this alternate possibility is forbidden by our contour choice because the real part of the effective potential (3.16) must be minimized on D ρ ∩C ± , where D ρ is the support of the spectrum [20]. However, a general fact of matrix models [20] is that V eff decreases along the direction normal to D ρ (that is, V eff decreases for some amount along the blue and green dashed lines on the real axis in Fig. 21 (b)).
Let the support of the one-cut spectrum be [a, b]. For a one-cut solution of a polynomial potential V with degree d + 1, we will closely follow the techniques in Section 3.3.3 of [20]. As R(x) is a two-fold covering of the complex plane with a cut along [a, b], the topology of R(x) is a Riemann sphere. A biholomorphism between the complex plane (with coordinate z) and the -1 Figure 22: The potential V (E) (blue), (rescaled) spectral density 2πρ(E) (shaded yellow), and effective potential V eff (E) (dashed red) in the one-cut case for different t.
Riemann sphere (with coordinate x) is given as follows: This biholomorphism is called the Joukowsky map. The upper (lower) sheet is mapped to the exterior (interior) of the unit disc, and the cut is mapped to the unit circle |z| = 1. Sending z → 1/z amounts to switching the sheet for fixed x. The inverse map is where the plus (minus) sign is for the upper (lower) sheet. Using (C.1), we can write where we choose the sign for the square root such that near x → ∞, σ(x) → x ∼ δz on the upper sheet, and σ(x) → −x ∼ −δ/z in the lower sheet. It follows that is a rational function of z with poles at z = 0 and z = ∞ because both V (x) and M (x) are polynomials in x. Let us denoteR(z) = R(x(z)). From the definition (3.19), R(x) → 1/x+O(1/x 2 ) for x → ∞ on upper sheet, which implies thatR(z) ∈ C[1/z], or where we used the fact that both V (x) 2 and M (x) 2 σ(x) have degree 2d. As V (x) = R(x + i ) + R(x − i ) for x ∈ [a, b], we have V ( a + b 2 + δ(z + 1/z)) =R(z) +R(1/z) = d k=0 v k (z k + z −k ) (C.6) Using (C.5) and (C.6), one can determine a and b and thus the spectrum. Taking our potential into (C.6), we have v 0 = 1 8 (3a 2 + 2ab + 3b 2 − 4(a + b)t) (C.7) Solving with the condition (C.5) leads to a = a(t) and b = b(t), though an explicit analytic expression is not available. It follows that For x ∈ [a, b], the spectral density is given by For x > b or x < a, the derivative of the effective potential is given by an analytic continuation of ρ(x), V eff (x) = M (x) σ(x), (C.12) from which we can compute the effective potential. The critical point occurs when x + a+b−2t 2 = 0 for x = a, which means that near a, the spectrum scales like (x − a) 3/2 . Solving this condition with a(t) and b(t), we find two cases: where t cr = 2 2/3 √ 3 ≈ 2.749. For t ∈ (−t cr , t cr ), one can show that no real solution exists for a(t) and b(t). We plot the one-cut spectrum for all cases with |t| ≥ t cr in Fig. 22, where we see that for |t| > t cr , the effective potential V eff (x) (dashed red lines) in the region x < a increases and then decreases. For |t| = t cr , V eff (x) monotonically decreases for x < a and has zero derivative at x = a. This is similar to our matrix model of JT gravity with EOW branes in (4.6).
Because a one-cut solution does not exist for t ∈ (−t cr , t cr ), we need to consider two-cut solutions. There are three possibilities that respect the symmetry of C ± : two cuts symmetric under complex conjugation (Fig. 21 (c)), one cut along the real axis plus one cut along a curve symmetric under complex conjugation (Fig. 21 (d)), and a special "Y" shaped case of the previous type obtained by joining the two cuts at a junction (Fig. 21 (e)). 14 The first case cannot be stable because the Coulomb repulsive force between the two cuts cannot be balanced by V (x). The second case is not allowed due to the requirement that V eff is minimized on D ρ ∩ C ± . The third case is the only choice. Physically, for |t| small, the potential well along the real axis is not deep enough to hold all of the eigenvalues; this naturally leads to complex eigenvalues along C ± . Let c ± label the two complex ends of the "Y" shaped solution, let b label the real end, and let a label the junction (see Fig. 23a). Let Y ± denote the two complex branches. As the degree of M 2 σ is d = 2, M must be a constant, θ. It follows that )). (C.14) Let c ± = c r ± ic i . The asymptotic condition for R(x) when x → ∞ requires that R(x) → 1/x, which implies that θ = 1, (C.15) To solve for c r , we still need one more equation. This is given by requiring that ρ(x)dx is a real and nonnegative measure along Y ± . This condition fixes the locus of Y ± and also implies that Note that the above integral is invariant under deformations of the integration path. We choose a straight line for computational simplicity. We plot c r , c i , a, b as functions of t ∈ (−t cr , t cr ) in Fig.  23. From the plot, we see that at the critical point t = ±t cr , c i = 0 and c r = a, which means that the "Y" shaped phase continuously connects to the one-cut phase. Given a t ∈ (−t cr , t cr ), the locations of the branches Y ± are determined by the differential equation (ρ(x)dx) = 0 ⇐⇒ dx r /dx i = −ρ r (x)/ρ i (x), (C. 19) where x = x r + ix i and ρ(x) = ρ r (x) + iρ i (x), and the boundary conditions x r (c i ) = c r (or x r (−c i ) = c r ) and x r (0) = a. As an example, we plot D ρ of the "Y" shape for t = 3/2 in Fig. 23a and the spectral density in Fig. 23b. Note that on Y ± we can project the density measure to the real axis as follows: ρ(x)dx = |ρ(x)| 2 /ρ r (x)dx r (C.20) where we have used (C.19).

D "Y" shaped spectrum
Universal piece For the universal piece, we have where C is the contour circling the branch cut γ + ∪ γ − ∪ [E 0 , 0] anti-clockwise. We expand the sine function in a Taylor series and evaluate each term by moving the contour to infinity. It turns out that the integral becomes the sum over the residues at E and ∞, (D. 2) The first term contributes to δρ U (E) = −ρ JT (E) for E > 0 and cancels out ρ JT (E) as in the one-cut case. Therefore the spectral density only comes from the second integral around infinity. Under the coordinate transformation λ → 1/z, it becomes where in the second step we used a Taylor expansion and picked out the coefficient of z n−1 . Let us denote the sum over all n i as S n (E i /E) (i = 0, ±). We will rewrite it in a form similar to the one-cut case. First, we expand each (E i /E) n i = ((E i /E − 1) + 1) n i in powers of (E i /E − 1) where in the second line we defined s i = n i − k i and changed the order of sum. In the third line we used the identity k i ≤m (a 1 ) k 1 · · · (a n ) kn k 1 ! · · · k n ! = Γ(1 + m + n i=1 a i ) Γ(1 + n i=1 a i )Γ(1 + m) , (D. 5) and in the fourth line we used the identity Γ(k 1 ) · · · Γ(k n ) Γ(1 + k 1 + · · · + k n ) = Vn t k 1 −1 1 · · · t kn−1 n n i=1 dt i , V n ≡ {(t 0 , · · · , t n−1 )|t i ≥ 0, n−1 i=0 t i ≤ 1} (D.6) and in the last line we computed the sum over k i which gives a simple power function of order n − 1. The sum over n is straightforward, (D.7) As I 1 (x) is an entire function, the integrand is analytic for all E. Therefore, the discontinuity of δR U (E) is purely determined by the factor (E − E + )(E − E − )(E 0 − E). We choose the branch cut to be along γ + ∪ γ − ∪ [E 0 , +∞). For E ∈ [E 0 , +∞), we have (D.8) For E ∈ γ ± , ρ(E) differs from the above expression by a phase that depends on the locus of γ ± , which is determined by requiring that ρ(E)dE is real and nonnegative along γ ± .

Potential related piece
Using a similar trick, we have where D ρ ≡ γ + ∪ γ − ∪ [E 0 , +∞) and C circles around D ρ anti-clockwise. 15 Similarly, we can move the contour C and this leads to the residues at E and −α 2 /2φ b . The residue at E obviously does not contribute to δρ K (E). The residue at −α 2 /2φ b gives 10) where (λ − E + )(λ − E − ) in the denominator contributes a minus sign by our branch cut prescription. For E ∈ [E 0 , +∞), this leads to and for E ∈ γ ± the expression depends on the locus of γ ± .

Large E limit
Taking the large E limit of ρ(E) − ρ JT (E) only requires (D.8) and (D.11) and not the specific solution of γ ± . One way to see the large E behavior is to go back to (D.3). Using the identity The large E limit of ρ K (E) is straightforward We can take the large E limit inside the dα integral on the second line only when using a regulated measure.
Integration of ρ(E) from E 0 to E + along the straight line Define E ± = E r ± iE i , E r,i ∈ R and E = uE + + (1 − u)E 0 . The variable in the Bessel function in x. Solving t ± in terms of x and y, the integral of the universal piece ρ U (E) from E 0 to E + along the straight line is dt 0 F(E 0 − E 0r y + iE i x)(x + y) 2 F 1 (− 1 2 , 3 2 , 2, where we defined q = x+y 2−2t 0 +x−y in the last line. The integral of the potential related piece is where F 1 is Appell function.