Solving puzzles in deformed JT gravity: phase transitions and non-perturbative effects

Recent work has shown that certain deformations of the scalar potential in Jackiw-Teitelboim gravity can be written as double-scaled matrix models. However, some of the deformations exhibit an apparent breakdown of unitarity in the form of a negative spectral density at disc order. We show here that the source of the problem is the presence of a multi-valued solution of the leading order matrix model string equation. While for a class of deformations we fix the problem by identifying a first order phase transition, for others we show that the theory is both perturbatively and non-perturbatively inconsistent. Aspects of the phase structure of the deformations are mapped out, using methods known to supply a non-perturbative definition of undeformed JT gravity. Some features are in qualitative agreement with a semi-classical analysis of the phase structure of two-dimensional black holes in these deformed theories.


Introduction
Low dimensional gravity theories have proven to be extremely fruitful playgrounds for exploring important questions in quantum gravity that are substantially more difficult to address in higher dimensions. One of the simplest examples is Jackiw-Teitelboim (JT) gravity [1,2], a two-dimensional theory including a scalar (or "dilaton") with the following Euclidean action: where χ is the Euler characteristic of the spacetime and S 0 a constant. This simple theory captures the low energy sector of both near-extremal black holes and branes [3][4][5][6][7][8] and the SYK model [9][10][11][12] (see ref. [13] for a review). Most importantly, the full quantum theory,

JHEP04(2021)030
where the spacetime is allowed to change topology, has been recently solved to all orders in perturbation theory in the topological expansion [14]. This was accomplished by using a variety of powerful mathematical results [15][16][17] to show that the full topological expansion of the Euclidean partition function: is reproduced by a double-scaled [18][19][20] random Hermitian matrix model. Given this success, it is reasonable to explore more complicated two-dimensional dilaton theories, obtained by deforming the JT gravity action (1.1). Studies of such deformations were recently initiated in refs. [21,22], by considering the following class of theories: (1.3) For this to be a deformation of JT gravity, U (φ) is required to vanish at φ → +∞ and so α i ≤ 1. Further constraining α i ∈ (0, 1/2), the Euclidean partition function of this theory was computed in refs. [21,22], where the exponential terms in U (φ) were interpreted as inserting sharp defects in the ordinary JT path integral. In a similar fashion to the work of ref. [14], the topological expansion of the partition function was shown to agree with that of a double scaled Hermitian matrix model. 1 The matrix model expansion is defined from the leading genus spectral density ρ 0 (E), obtained from the gravitational disc partition function after an inverse Laplace transform. All higher genus contributions are then obtained from ρ 0 (E) via recursion relations derived from the matrix model.
While for generic deformations U (φ) this procedure works nicely, it was noted in refs. [21,22] that in certain cases the leading genus spectral density ρ 0 (E) provided by the disc partition function is not positive definite. Given that ρ 0 (E) is the seed of the recursion relation, starting from a sick spectral density results in an ill-defined perturbative expansion. This puzzle was left unresolved in refs. [21,22], but it was speculated that either non-perturbative effects or phase transitions might resolve the issue. This paper was motivated by the desire to better understand these deformations and determine whether the problems can be resolved using recently developed techniques that provide an alternative formulation of the matrix model description of JT gravity in terms of certain combinations of minimal string models [24,25]. As emphasized in refs. [25][26][27], this formulation is particularly well-suited to enable the extraction of non-perturbative physics. Fully specifying the matrix model is equivalent to determining a function u(x) for x ∈ R, which is the potential of an associated quantum mechanics problem with = e −S 0 . While the function u(x) satisfies a complicated non-linear ordinary differential equation (called the "string equation"), its leading piece as → 0, u 0 (x), satisfies a simple algebraic constraint.
A central result of our work is the observation that the non-positive spectral density ρ 0 (E) arising in some JT gravity deformations is only a symptom of the sickness of JHEP04(2021)030 the models, and the actual virus is a multi-valued potential u 0 (x). Sometimes the problem can be eradicated by identifying a phase transition, and at other times it cannot. In such cases, we give evidence that the non-perturbative physics (at least using the definitions employed here) does not exist. Furthermore, we explore the physics of the deformations in the semi-classical gravity approximation, the leading behavior of the matrix model, and the full non-perturbative definition, contrasting the three approaches and identifying several phenomena. We now summarize our core results in more detail, briefly explaining the salient features.

Summary of results
We start in section 2 by introducing the matrix model technology used throughout this work. An important issue we address is the non-perturbative formulation. The Hermitian matrix model definition [14] of JT gravity is afflicted by non-perturbative instabilities long known [28,29] to be present in certain double-scaled Hermitian matrix models. In ref. [25] an alternative matrix model definition of JT gravity was provided which is free from such instabilities and reproduces the perturbative physics of ref. [14]. Here we present an extension of ref. [25]'s framework that is able to naturally incorporate the deformed JT gravity physics and serve as a non-perturbative definition. It reduces to the nonperturbatively stable definition when the deformations are removed and will therefore act as a sharp tool for determining if some deformations are non-perturbatively ill-defined (at least within a framework that works for ordinary JT gravity). The matrix model origins of the string equation are discussed.
Section 3 discusses a particular example of a deformation of JT gravity given by: which satisfies U (0) = 0. To start, we compute the Euclidean partition function in the semi-classical approximation and uncover an interesting variety of phase transitions of the two-dimensional black holes of the theory (see figures 2 and 3). We continue by studying the full Euclidean partition function to leading order in genus (disc topology), computed in refs. [21,22]. The leading genus spectral density obtained from the results in ref. [22] is given by: (1.5) where ≡ e −S 0 and the threshold energy is E 0 = 0. Expanding this expression for small energies, one finds ρ 0 (E) ∝ √ E + O(E 3/2 ), where the proportionality factor in the leading term is negative for λ > λ c , see (3.15).
To understand the origin of this issue from the matrix model perspective, the model is reinterpreted as a particular combination [24] of minimal models, which enables a formulation to all orders in and beyond [25]. The central quantity in this formalism is the string equation, that determines a potential u(x). The leading genus string equation associated to this model is given by:

JHEP04(2021)030
(where I n is the nth modified Bessel function), which provides an implicit definition of the leading genus potential u 0 (x) ≡ lim →0 u(x). We first show that ρ 0 (E) in (1.5) goes negative whenever the corresponding potential u 0 (x) defined through equation (1.6) becomes multi-valued, see figure 4. From this perspective, the issue is easily solved by picking the (unique) solution to the implicit equation (1.6) that yields a single valued function u 0 (x) in the region x < 0. This modifies the leading spectral density from (1.5) to: where E 0 ≥ 0 is determined from the largest solution to R 0 [E 0 , 0] = 0. While for deformations for which E 0 = 0 the integral (1.7) reduces to (1.5), when E 0 > 0 the correct positive definite answer is instead obtained from the integral expression. The transition of E 0 between the two regimes is not analytic (see figure 5) and corresponds to a phase transition which qualitatively matches with results obtained from the semi-classical analysis. 2 At this point we must comment on the resolution of this first puzzle and its relation to the calculations in ref. [21]. In that work the disc partition function was shown to be given by (1.7) irrespective of the value of U (0). As a result, the gravitational computations in ref. [21] already have the correct spectral density when U (0) = 0, given by equation (1.7) instead of equation (1.5). However, the issue is relevant from the perspective of ref. [22], where the gravitational calculation with U (0) = 0 unequivocally leads to the poorly defined result in (1.5). It would therefore be interesting to understand how the derivation of (1.5) in ref. [22] must be modified in order to yield the correct result.
Section 3 concludes with a computation of the full spectral density ρ(E) including both higher genus and non-perturbative contributions (see figures 6 and 7). Non-perturbative effects generate corrections to ρ 0 (E) in (1.7) analogous to those previously computed for ordinary JT gravity in refs. [25,27], with the difference that in this case E 0 can be non-zero. Section 4 repeats the analysis for an example of a different class of deformation of JT gravity. Here, U (0) = 0, and: After computing the Euclidean partition function in the semi-classical approximation, the full partition function to leading order in genus [21,22] is studied. The spectral density in this case is given by 3 with E 0 again determined from the largest solution to R 0 [E 0 , 0] = 0. While the integral expression is the same answer as given in (1.7) for U (0) = 0, when U (0) = 0 the result of 2 We should note that the semi-classical and matrix model computations apply to different αi regimes, αi ∼ 1 and αi ∈ (0, 1/2) respectively. While certain quantitative aspects do not match exactly (like the exact location and order of the transition), the qualitative agreement between the two approaches is quite remarkable. 3 We are using the same conventions as in ref. [21]. This is obtained from ref. [22] by changing the integration variable in eq. (8.12) to u0 = (1 − s 2 )E + s 2 E0.

JHEP04(2021)030
the integral is not necessarily positive (see figure 9 for some examples). As in section 3, the issue originates in a multi-valued solution u 0 (x) to the leading genus string equation (1.9), but in this case in the region x < 0 (see figure 9). This makes the issue substantially different and ultimately results in a breakdown of the perturbative expansion of the matrix model. We also show that certain models which might appear to be well defined (since they have ρ 0 (E) positive definite) still suffer from this perturbative breakdown.
To determine whether the model can be made sense of even when the perturbative expansion breaks down, the system is studied non-perturbatively. We find that the nonperturbative definition also breaks down precisely when u 0 (x) becomes multi-valued in the region x < 0. Hence, we argue that deformations of JT gravity for which u 0 (x) presents such multi-valuedness issues are both perturbatively and non-perturbatively unstable.
The paper concludes in section 5 with a brief discussion and some thoughts about future directions. Several appendices contain details and results used in the main text.

Matrix models
This section introduces some of the random matrix model tools used throughout this work. The introduction will be rather brief (for longer treatments see e.g., refs. [30][31][32]) but will help establish some notation, and contextualize key aspects of the perturbative and non-perturbative formulations.

Double-scaling limit
Both random Hermitian matrix models and random complex matrix models will be relevant here, but it is traditional to begin with the former. The model is defined through a probability measure dH e −N Tr V (H) with H an N ×N Hermitian matrix and V (H) a polynomial potential generalizing Wigner's prototype Gaussian case. Expectation values of matrix observables O are computed as: One of the central observables of the theory is the spectral density where α i ∈ R are the eigenvalues of H. 4 In the N → ∞ limit ρ(E) becomes a smooth density function, beginning at some classical threshold energy E 0 . The central idea is to use random matrix models to compute universal properties of sums over random two-dimensional surfaces. Those surfaces can be seen [33,34] to be present in the definition above by expanding the integral Z in terms of Feynman diagrams (the quadratic term of V (H) gives a propagator, higher order terms give vertices), which are graph duals of tesselations of the surfaces into polygons. The power N χ is carried JHEP04(2021)030 by a graph/polygonization with Euler characteristic χ = 2 − 2g − b where g and b are the number of handles and boundaries, respectively. Universal physics from sum over surfaces (i.e., properties that do not depend on model-dependent details such as whether triangles or pentagons were used for the discretization) are picked out by taking the double-scaling limit [18][19][20]. This involves taking N → ∞, while simultaneously tuning the polynomial coefficients in the potential V (H) to a critical point dominated by smooth surfaces that are very large compared to the number of constituent polygons. The 1/N corrections yield the topological expansion. The critical point is characterized by the rate of vanishing of the leading spectral density ρ(E) at one of its endpoints (see e.g., ref. [35]). The generic behavior is Wigner's semi-circle law: ρ ∼ (E − E 0 ) 1 2 . An additional (k − 1) zeros coincident at the endpoint E 0 define the critical behavior corresponding to a distinct model [36], the kth model: ρ ∼ (E − E 0 ) k− 1 2 . 5 Non-perturbative physics that lies beyond the topological 1/N expansion is accessible too. A starting point for extracting it [37] is a family of polynomials P n (α) = α n + lower powers (where n = 1 · · · N ), that are orthogonal with respect to the measure dµ and for which a recursion relation can be written (e.g. for even V (H)) as: αP n (α) = P n+1 (α) + R n P n−1 (α). All matrix model observables can be re-written in terms of the recursion coefficients R n and identities satisfied by the defining matrix integral can be used to derive recursion relations for them. The P n (α) actually define a finite dimensional Hilbert space description of the system, built from |n ∝ P n (α), with n| m = δ nm , within which observables, as discussed in equation (2.1), are associated with operators.
In the double scaling limit, the Hilbert space supplied by the P n (α) becomes infinite dimensional; |n → |x with x ∈ R, and an effective quantum mechanical system emerges from which matrix model observables can be computed [18,28]. It is governed by a simple Schrödinger Hamiltonian: withx |x = x |x , and the scaling part of 1/N in the limit. The potential u(x) arises as the scaling part of the R n in the limit, and its functional form is determined by an ordinary differential equation called the "string equation" (to be discussed below) which is the scaling part of the aforementioned recursion relations for the R n . An important observable is the expectation value of a "macroscopic loop" (a single boundary) of length β, which in this formalism is: It is important to notice the integration range for x. The double scaling limit captures the universal physics to be found in the infinitesimal neighbourhood of the spectral density's JHEP04(2021)030 endpoint (where the critical behaviour was specified earlier). Recall that the index n labeled the ordered N energies, and so at large N they are labeled by a continuous coordinate X = n/N ∈ [0, 1]. The coordinate x ∈ (−∞, +∞) describes the local physics of the X = 0 endpoint, and in the present conventions they are related as X = 0 − xδ p (where p is some positive power and δ → 0 as N → ∞). Crucially, the x < 0 regime (integrated over in equation (2.5)) is the place where the surviving part of the spectral density classically has support, and x = 0 is the location of the endpoint. Very importantly, note that the full quantum mechanics is defined on the whole of the range of x. The physics from the x > 0 regime controls non-perturbative aspects of the system. The part of the spectral density that survives in the double scaling limit, ρ(E), can be written in this language too. Starting with equation (2.5) and inserting a complete set of eigenstates of the Hamiltonian, H |ψ E = E |ψ E : where E 0 is the lowest or "threshold" energy in the spectrum. The macroscopic loop expectation value is denoted Z(β) in the above (as done in equation (1.2)) because it is of the same structural form as JT gravity, i.e. a sum over all surfaces with a single boundary of fixed length β. Evidently, the central object in this formalism is the potential u(x) that determines the Hamiltonian (2.4). For the appropriate choice of u(x), equation (2.6) will be a definition of the JT gravity matrix model that is alternative to the recursion methods of ref. [14]. The core idea is that the full non-perturbative behavior can be accessed from u(x) obtained as the full solution to a differential equation, the string equation.
The all-orders perturbative expansion for u(x) is also nicely packaged in the string equation. However, it is important to note here that the string equation arising from Hermitian matrix models is non-perturbatively unstable (this will be discussed in section 2.2). Therefore, the Hermitian matrix model string equation that we now write, should only be used perturbatively in [18][19][20]: (2.7) The objects R k [u] are kth order polynomials in u(x) and its derivatives called the Gel'fand-Dikii polynomials [38], normalized so that the coefficient of u k is unity. They satisfy a simple recursion relation [38]: where R 0 [u] = 1 and primes are derivatives with respect to x. The coefficients t k in equation (2.7) indicate how much the kth model contributes to the form of the potential u(x) and essentially define the double scaled model. To construct the matrix model equivalent to JT gravity, a natural strategy is to determine the combination of t k that yields the leading perturbative ρ 0 (E), known from disc level computations in JT gravity. This determines the leading contribution to u(x), denoted u 0 (x). This is described next.

Leading perturbative analysis
Writing a power series expansion for the potential u(x) = ∞ n=0 u n (x) n and inserting it into the string equation (2.7), a perturbative analysis of the physics can be performed. The leading behavior (the classical limit in the quantum mechanics of equation (2.4)) of the potential u 0 (x) satisfies the following simple algebraic equation: Subtleties regarding the implicit definition of u 0 = u 0 (x) from this equation will play an important role in this paper. A simple and useful formula can be written for the loop expectation value (the putative JT gravity partition function) Z(β) defined in equation (2.5) at leading (disc) order in genus Z 0 (β) = lim →0 Z(β) in terms of u 0 (x). Using that the momentump = −i ∂ x and position operatorsx commute to leading order in and inserting a complete set of momentum eigenstates |p we can solve the Gaussian integral in p and find: where we have used p|x = e ixp/ / √ 2π . Changing the integration variable to u 0 and using (2.9), we can obtain the following formula for the spectral density ρ 0 (E) after applying an inverse Laplace transform: This formula has played a central role in recent JT gravity explorations [24,25]. In fact, in the mathematical literature this integral transform (when E 0 = 0) is known as the Abel transform, and in such cases, it can be explicitly inverted. This is described in appendix A, the final result is: This very useful simple formula allows us to write the leading genus string equation from the spectral density ρ 0 (E) and ultimately identify the coefficients t k that define the model. For example, in the case of undeformed JT gravity the result is [24,25]: and explicitly in terms of the t k this is: (2.14) As will become clear later, deformations of JT gravity will yield more general ρ 0 (E). This will imply different formulae for the t k , and hence different matrix model definitions (i.e., new combinations of the underlying minimal models). The resulting u 0 (x) can then be used as the seed for higher order corrections, as part of the boundary conditions of the string equation that supplies non-perturbative physics. The latter will be described next.

Non-perturbative completion
A complete definition of JT gravity should include a non-perturbative sector, and in this approach this means seeking contributions to u(x) that cannot be captured in an expansion in . As already mentioned, the Hermitian matrix model string equation (2.7) captures perturbation theory, but fails to include consistent non-perturbative physics. This issue has a long history, and can be traced to the fact that the even k minimal models defined by the Hermitian matrix model are individually non-perturbatively unstable [28,29]. Since the Hermitian matrix model proposed by ref. [14] includes (as reviewed in the previous section) the contribution of the k even models, they inherit the instability. 6 Overall, this means that the non-perturbative definition of the double scaled Hermitian matrix model is unstable and must be modified. The function u(x) must satisfy a different differential equation. In this work, the differential equation used can be obtained from two different perspectives: from an ensemble of complex (instead of Hermitian) random matrices or from a scaling argument. Our proposal is a slight generalization of the nonperturbative definition of JT gravity recently given in ref. [25], that incorporates cases for which E 0 = 0.
A different family of non-perturbately stable minimal models can be obtained from an ensemble of complex matrices M with probability measure dM dM † e − N γ Tr V (M M † ) . The double scaling limit of this system was obtained in [39][40][41], where much of the formalism described in section 2.1 for Hermitian matrix models also applies (for a thorough exposition, see ref. [40]). The crucial difference is given by the fact that the potential depends on the combination M M † ≥ 0, and so its eigenvalues α i are non-negative. The spectral density is therefore non-zero only on the positive real line. It is instructive to think of this as structurally the same problem provided by the Hermitian matrix model, but with a "wall" at zero that prevents the energies from going below zero (presently we will extend the discussion to E 0 = 0). The procedure for tuning to a critical point is analogous to what was done before. In terms of the endpoint of the distribution ρ(E), for the kth model there are again (k − 1) additional zeros at the endpoint, but also the wall is present, which produces a key new feature. In terms of the variable x, the regime x < 0 contains the physics of what happens all the way up to the wall, while the x > 0 sector yields the non-perturbative (stabilising) physics supplied by the wall. It is the presence of this latter regime that makes complex matrix models radically differ from Hermitian matrix models, which contain no wall. 7 6 A suggestion for how to avoid this issue was given in ref. [27], obtained by thinking of the JT definition as a limit of a (stable) k odd model, with the even models as perturbations of it. However, it is not clear if this will work. Moreover, the overall picture developed in this paper suggests that the needed solution to the differential equation simply does not exist. See the remark at the end of section 4.3. 7 In fact, it is easy to see that if one insists on using Hermitian matrix models, one can simply place a wall into such a model by hand, and reproduce the physics under discussion, as discussed in e.g., ref. [42]. Alternatively since Hermitian matrix models are not really any better motivated for studying sums over two-dimensional surfaces than complex matrix models, one might as well start with complex matrix models which supply the wall "naturally". The string equation that results from taking the double-scaling limit of these models is [40,42,43]:

JHEP04(2021)030
where R is defined as in (2.7). Note that the Hermitian matrix model string equation R = 0 is a solution of this equation, which results in the potential u(x) agreeing to all orders in perturbation theory when expanding (2.15) in the x < 0 regime, as is appropriate for JT gravity. The proposal of this paper is that for the deformations of JT gravity (as described in section 1) the equation (2.15) gives the appropriate non-perturbative definition.
For E 0 = 0 it generalizes the previous proposal of ref. [25], used in ref. [27] to explicitly compute the ordinary JT gravity non-perturbative potential u(x) and spectral density ρ(E) shown in figure 1. 8 Let us comment on the E 0 contribution in (2.15), that is the novel piece with respect to the proposal in ref. [25]. The leading order behavior u 0 (x) obtained from (2.15) as → 0 is given by If we want to have a single-valued and non-trivial function u 0 (x) that is defined in the whole real line x ∈ R, none of the regimes is enough by itself, meaning there must be a transition between them. Since leading genus observables are determined by u 0 (x) in the region x < 0 (2.10), we should take u 0 (x) defined from R 0 [u 0 , x] = 0 when x < 0. In this way, leading genus observables computed from (2.15) agree with those obtained from the original Hermitian matrix model string equation. The full solution for u 0 (x) is then

JHEP04(2021)030
given by In order to have a continuous solution, we must fix E 0 as R 0 [u 0 (0) = E 0 , x = 0] = 0, which gives the following condition The threshold value E 0 is the energy at the endpoint of the spectral density, where the wall is, at x = 0. The origin of the shifted equation (2.15), which first appeared in ref. [43], can be understood as a simple modification of the complex matrix model by shifting the matrix combination M M † → (M M † + E 0 ), so that the eigenvalues lie above the threshold E 0 . The region x > 0 is beyond the wall, and so classically u 0 (x) should stay frozen at E 0 , the value it reached at x = 0. The boundary conditions used to solve the full string equation (2.15) are obtained from the asymptotic expansions for x → ±∞ starting from each of the two regimes in (2.17).
There is a different way of motivating the full differential equation (2.15), following the work in refs. [40,42]. Instead of going through the double scaling limit of complex matrices, one seeks an equation that has the same perturbative physics as given by the Hermitian matrix model (already known to agree with JT gravity perturbatively) but then assumes that it preserves the following two established properties of u(x) not just in perturbation theory but non-perturbatively: where the first line refers to the explicit dependence of the potential u(x) on the coefficients {t k }. Its evolution as a function of t k is governed by the KdV flow equation [28,45]. The second line is a simple scaling relation that specifies how u(x, {t k }, E 0 ) changes with a rescaling of the parameters by s. 9 Appendix B reviews how to use these conditions, with (2.17), to arrive at the string equation (2.15). All things considered, the non-perturbative definition given in this section is a well motivated proposal that allows for explicit computations of physical quantities in deformations of JT gravity.

Deformed JT gravity: model A
In this section we study a particular deformation of JT gravity in (1.3), where the dilaton potential is deformed according to

JHEP04(2021)030
where U (0) = 0 and for simplicity we assume α 1 < α 2 . In the following sections we analyze the Euclidean partition function using the semi-classical approximation as well as the matrix model description. Several interesting characteristics of this theory will be uncovered. Overall, it will emerge as both perturbatively and non-perturbatively well defined.

Semi-classical approximation
A standard semi-classical analysis of the deformed JT gravity theory involves computing the Euclidean partition function through a saddle point approximation: µν and φ (i) are on-shell solutions to the equation of motion with renormalized boundary length β. This procedure yields a good approximation when α i ∼ 1 (see appendix D of ref. [21]).
The first step in computing (3.2) is solving the equations of motion with a dilaton potential given by W (φ) = 2φ + U (φ) and the appropriate boundary conditions. Luckily, this has been analysed a long time ago and the general solution can be written as [46][47][48]: where r ≥ φ h and t ∼ t + β. 10 Black hole solutions are labeled by φ h , which must satisfy the following constraint in order to have a well defined signature: Evaluating the on-shell action (3.2), a standard computation yields the following thermodynamic quantities [48]: From this we can compute the specific heat and find the following stability condition for black hole solutions Let us now specialize to the dilaton potential in equation (3.1). The black hole solutions and their behavior can read off just by plotting the potential W (φ) = 2φ + U (φ) for fixed (α 1 , α 2 ) ∼ 1 and λ. Figure 2 shows sample characteristic behavior of the potential depending on the value of λ. Segments of the curve in red, blue and green indicate whether the black hole solution with φ h = φ is nonexistent (violates (3.4)), unstable (negative specific heat) or stable (positive specific heat). The behavior can be quite different depending on the sign of λ, so each case is analysed separately. 10 The extra factor of two in the time coordinate is such that the boundary conditions are satisfied with β matching with the conventions in ref. [21]. Positive λ: since in this case the potential has W (±∞) = +∞ there are two branches of black hole solutions, one stable and the other unstable. While at low enough temperatures there are only stable black holes, as we increase the temperature the unstable branch comes into play. In this regime, the favored black hole that dominates the partition function (3.5) is the one with lower free energy. A simple analysis using (3.5) (see also section 3 of ref. [48]) shows the stable black hole in the green branch always dominates the partition function over the unstable one in blue. As a result, the overall system is stable in the semi-classical approximation.

JHEP04(2021)030
Note that there is an interesting behavior for the zero temperature black hole solution φ h = φ 0 , determined by the largest root of W (φ 0 ) = 0. While for small enough λ we always have φ 0 = 0, there is a transition at some finite value λ c in which φ 0 becomes non-zero and positive. The critical value λ c can be easily computed as .
Solving for the zero temperature black hole solution φ 0 = φ 0 (λ) in a series expansion around λ c we find the following non-analytic behavior: The full function φ 0 (λ) can be plotted numerically by solving W (φ 0 ) = 0, so that we get the left diagram in figure 3. This corresponds to a phase transition of the zero temperature black hole solution. The order of the transition is obtained from the free energy F = E−T S, which at T = 0 agrees with the energy (3.5):   .7), where we plot φ h = φ 0 (λ) obtained by numerically solving W (φ 0 ) = 0. In the right diagram is plotted the black hole energy as a function of the temperature for λ = −40, taking into account the phase transition (left diagram in figure 2). We observe a gap in the spectrum generated by the phase transition at λ below λ c in (3.11).
find that: 11 , the semi-classical transition of the zero temperature black hole solution is third order.
Negative λ: in this case an interesting structure arises for negative enough values of λ. The potential develops a local maximum and minimum, see leftmost diagram in figure 2. This happens for λ < λ c given by: It is straightforward to check that when λ = λ c the potential has an inflection point meaning that for λ < λ c a local maximum and minimum develops. The behavior of the dilation potential W (φ) in this regime induces a phase transition of a different kind from the one discussed for positive λ. In the canonical ensemble (at fixed temperature T = W (φ h )/2π), we identify a shaded region in the left diagram in figure 2 where there is a finite range of temperatures with three black hole solutions. The dominant black hole in this regime is the one with the lowest free energy F = E − T S. From the equations in (3.5), a simple analysis shows that seeking the dominant black hole produces JHEP04(2021)030 a discontinuous jump from φ 1 to φ 2 , as indicated in the leftmost component of figure 2. In the right diagram of figure 3 is plotted the energy as a function of the temperature (see eqs. (3.5)), taking into account the leap in φ h . At the critical temperature we observe a discontinuous jump in the energy. This generates a gap E gap for which there are no dominant black hole solutions in a finite range of energies. This kind of first order phase transition (where three black holes compete) was first observed for charged AdS black holes in ref. [49]. Along with other transitions, they have been recently explored in the current JT gravity context in ref. [48].

Perturbative genus expansion
While the semi-classical analysis of the previous section showed a rich and interesting structure, it was obtained from the crude saddle point approximation to the Euclidean path integral (3.2). To compute Z(β) more carefully we use the results in refs. [21,22], where it was shown the full partition function can be computed from a double scaled matrix model. In this subsection we focus on the leading genus contribution, Z 0 (β), given by the disc topology: where we are neglecting additional higher genus and non-perturbative contributions that will be studied in the next subsection. The leading genus disc spectral density (1.5) derived in ref. [22], that we rewrite here for convenience, is given by: (3.14) This formula has a problem, since ρ 0 (E) is not always positive definite. This can seen by expanding for small energies, where we find: 12 . (3.15) Here λ c > 0 since we are assuming α 1 < α 2 . The spectral density goes negative for λ > λ c and small enough energies. 13 This certainly does not fit within the matrix model analysis and framework outlined in section 2.2. In particular, the definition of ρ(E) in equation (2.2) is manifestly positive, and large N and double-scaling limits cannot change that, as is clear from equation (2.6). To understand the origin of the issue from the matrix model perspective, consider the behavior of u 0 (x), defined from the leading genus string equation. Using equation (2.12) we can easily compute R 0 and find (3.16) 12 We have suggestively defined λc, which is different from the semi-classical definition in (3.7). We comment on the relation between the two at the end of this subsection. 13 As already explained in the introduction, this issue does not arise from the perspective of ref. [21]. To have a well defined single valued function u 0 (x) on the real line x ∈ R, the implicit function theorem imposes the following constraint

JHEP04(2021)030
From (3.16) we can easily check this condition is never satisfied, not even for undeformed JT gravity (i.e., λ = 0). This might seem strange, since ordinary JT gravity has a well defined disc spectral density (2.13). However it is important to note that the leading genus partition function Z 0 (β) (2.10) is built from the regime x < 0. Therefore we only need R 0 [u 0 , x] = 0 to define a unique single-valued function u 0 (x) in the range x < 0.
To understand the situation more clearly it is useful to plot the curve x = x(u 0 ) defined through equation (3.16). The undeformed JT gravity case is the curve shown in the λ = 0 diagram in figure 4. While for u 0 < 0 there are an infinite number of oscillations around x = 0 (marked in red), for u 0 > 0 there is a single curve marked in green going from u 0 = 0 to u 0 = +∞. This green curve defines the unique single valued function u 0 (x) in x < 0. As can be seen from the other plots in figure 4, the behavior for λ < λ c is quite similar: the upper (green) branch can be used to define a unique single-valued function defined in the range x < 0. Using that in all these cases u 0 (0) = E 0 = 0, the integral (2.11) can be solved explicitly, recovering the starting result (3.14), which is positive for such λ.
Consider now λ > λ c . From the rightmost diagram in figure 4 it is clear that there is an interesting feature: the value of u 0 (0) is no longer zero. The shift is caused by the fact that exactly at λ = λ c defined in (3.15), the derivative of the leading genus string equation vanishes at u 0 = 0 (∂ u 0 R 0 ) (u 0 ,λ)=(0,λc) = 0 . that results in a single valued function u 0 (x) for x < 0, i.e. the green portion of the diagram, the ρ 0 (E) obtained from integral (2.11) is explicitly positive (since ∂ u 0 R 0 ≥ 0 in this case). The primary new feature is the value of u 0 (x) at the origin, the threshold energy E 0 = u 0 (0). For λ > λ c , instead of having E 0 = 0 it is a non-zero value obtained from the largest solution to the leading genus equation at x = 0 This equation can be easily solved numerically to yield figure 5. The behavior of E 0 as a function of λ is non-analytic at λ c and corresponds to a phase transition. We can determine the order of the transition using (3.19) to expand E 0 (λ) on either side of λ c , so that: (3.20) this is a first order phase transition. This phase transition solves the negativity of the leading genus spectral density, which for λ > λ c is given by (2.11) with u(0) = E 0 = 0 instead of the problematic formula (3.14).
It is interesting to compare these leading order matrix model results with those obtained in the previous subsection from the semi-classical gravity analysis. The most important similarity is given by the fact that in both cases we observe a phase transition of the threshold energy E 0 for some positive value of λ. By comparing figures 2 and 4 we also note some resemblance in the mechanisms triggering the transition. There are however important quantitative differences, since the exact location of the transition λ c in (3.7) and (3.15) does not agree. Moreover, while the matrix model transition is first order, the semi-classical analysis yields a third order transition. When doing this comparison, keep in mind there is no overlap between the α i regimes of validity of each computation. Moreover, the saddle-point approximation of the path integral (3.2) does not take into account quantum corrections that might be important at low energies. All things considered, it is quite pleasing (and perhaps surprising) that the semi-classical analysis captures the transition of the vacuum solution at all.

JHEP04(2021)030
On the other hand, the semi-classical gravity analysis yielded a phase transition among black hole solutions that is present at negative λ (below λ c in (3.11)). It suggested the existence of a gap E gap opening up in the available energy of black holes. There is no sign of such a gap in the spectrum in the matrix model analysis at disc level, but what about higher orders in in the perturbative expansion?
A simple way of showing perturbative corrections are not able to account for E gap is by computing the first genus correction to the disc spectral density ρ 0 (E) in (3.14). The partition function at this order is given in equation (9.19) of ref. [22]. Carefully applying an inverse Laplace transform (dealing with divergences in the same way as in footnote 9 of ref. [21]) we find the corrections to spectral density are of the form c 0 /E 3/2 + c 1 /E 5/2 for some constants c 0 and c 1 . Contributions of this form are not able to explain the energy gap E gap , which should arise from a low energy behavior given by: Neither the leading genus spectral density ρ 0 (E) in (3.14) nor its first correction contribute with the term ∼ 1/ √ E that is necessary to resum the series and obtain the gap. Overall, the energy gap E gap identified in the semi-classical analysis is not observed when computing the Euclidean partition function including contributions of arbitrary genus.

Non-perturbative effects
So far, the Euclidean partition function of the theory has been studied in one or other approximation scheme: either the semi-classical gravity limit (3.2) or the perturbative expansion of the matrix model (3.13). Using our double-scaled matrix model definition presented in section 2.3, we now compute Z(β) exactly, including all higher order genus and non-perturbative effects.
The first step is to identify the coefficients t k (λ) in (2.7) that define the matrix model. These can be obtained by expanding the leading genus string equation (3.16) in a power series in u 0 , giving: Using this, the combination R[u, x] in (2.7) can be constructed, and hence the non-linear ordinary differential equation (2.15) that determines the potential u(x). The boundary condition is the full = 0 solution for u 0 (x), which is given by the u 0 (x) discussed in the previous section for x < 0, and for x ≥ 0 by u 0 (x) = E 0 . Since all the t k (λ) are non-zero, the differential equation is of infinite order, but (as demonstrated in the case of λ = 0 in ref. [27]) there is a meaningful sense in which it can be truncated to finite order, by including all t k contributions up to a maximum, t kmax :

JHEP04(2021)030
This is a controlled truncation, as can be seen from the fact that the t k rapidly decrease in relative size as k increases. The truncation gives a good estimate for E 0 , λ c , and moreover gives an energy spectrum that matches well to the exact one. Adding higher k improves the truncation accuracy. It was found that k max = 7 gave good accuracy for the range of energies where non-perturbative effects are important. This required a solution of a highly non-linear 14th order differential equation, which turns out to be possible. 14 After finding an accurate enough numerical solution for u(x), the spectrum of the Schrödinger Hamiltonian H[u] of equation (2.5) can be computed (again, using numerical techniques), and ultimately the full spectral density ρ(E) computed using (2.6). 15 Picking the values (α 1 , α 2 ) = (1/4, 1/3), the next step is to analyse each sign of λ separately.
Positive λ: the potential u(x) was found, and the spectrum computed to enable the construction of ρ(E) for positive λ in the k max = 7 truncation. In the left diagram of figure 6 is plotted the potential u(x) for several values of λ, where the phase transition (3.20) in E 0 is manifest in the shift of the boundary condition u(+∞) = E 0 . Comparing with the ordinary JT gravity solution with λ = 0 (blue curve), it can be appreciated how the potential u(x) is deformed.
On the right diagram of figure 6 is plotted the spectral density ρ(E) computed from each of the potentials. For λ = 0 it is known that the threshold density value ρ(0) is non-zero due to non-perturbative effects [25,27]. While for λ ∈ (0, λ c ) we also observe ρ(0) non-zero (the slight numerical jitter near the threshold energy can be safely ignored), it takes a smaller value as compared to that of ordinary JT gravity (blue and red curves). As λ increases past the transition value λ c , the threshold energy E 0 shifts to the positive values obtained from equation (3.19). The value of ρ(E 0 ) is still non-zero and increasing with λ (compare green and black curves).
The behavior of ρ(E 0 ) is related to the value of the non-perturbative potential u(x) at x = 0. The spectral density ρ(E) is computed from (2.6) by integrating over the (modulussquared) tail of the E = E 0 wavefunction that penetrates into the x < 0 region. How much of that tail penetrates depends upon the value of u(x) − E 0 at x = 0. Here, note that u(x) is the full non-perturbative potential -its leading behavior u 0 (x) − E 0 of course vanishes at x = 0 by construction. Comparing the diagrams in figure 6 we observe that ρ(E 0 ) increases whenever u(x) − E 0 decreases, since |ψ(x, E 0 )| 2 (integrated to find ρ(E 0 ) in (2.6)) has a larger support in the region x < 0.
Negative λ: let us now consider negative λ, where there is no phase transition in the threshold energy E 0 = 0. On the diagrams in figure 7 we plot the non-perturbative 14 In practice, a first derivative of the string equation (and then dividing by R) results in an equation that is structurally less complicated, so a 15th order equation was solved for kmax = 7. The boundary conditions, obtained by fixing the appropriate number of derivatives ∂ n x u(x) ≡ u (n) (x) at x = ±xmax, are fixed to u (n) (±xmax) = u (n) 0 (±xmax). See the appendix of ref. [27] for detailed tips on how to numerically solve the equations, and to solve the spectral problem for H[u] once u(x) is found. Throughout, the choice = 1 was made, although explorations were done with < 1, with qualitatively similar results. 15 The same numerical techniques as the ones developed in refs. [25][26][27] were used, and so they will not be described here. See an appendix of ref. [27] for suggestions about solving the equations and extracting the spectrum of H. numerical solutions to u(x) and ρ(E) for several values of λ. The quantity ρ(0), the spectral density at the threshold energy E 0 = 0, increases for increasingly larger negative λ. From the discussion immediately above, this is in agreement with the value of the potential u(x) at the origin x = 0.

JHEP04(2021)030
Overall, we conclude that the deformation of JT gravity according to (3.1) is both perturbatively and non-perturbatively well defined (using the non-perturbative scheme adopted in this paper). For λ positive the phase structure obtained from the semi-classical analysis of the two-dimensional black holes solutions is in qualitative agreement with the full partition function as captured by the matrix model. For λ negative (we explored out to λ = −100 to ensure that the patterns of figure 7 persist) this does not appear to be the case, as the energy gap E gap observed in the semi-classical approximation is invisible in the matrix model description, even after accounting for non-perturbative effects. 16

Deformed JT gravity: model B
This section will repeat the previous the analysis of Euclidean partition function but for a different deformation of JT gravity, given by: where now U (0) = 0. We show that for certain values of λ this deformation is both perturbatively and non-perturbatively inconsistent, using our definition in section 2.3. 16 The first version of this manuscript reported effects that were interpreted as the appearance of a nonperturbative Egap. They were due to a subtle error in the code used to numerically solve the full string equation. This error has been corrected.

Semi-classical approximation
As done in section 3.1, the semi-classical approximation can be used to compute the Euclidean partition function (3.2) of the theory with potential W (φ) = 2φ + U (φ), but now with U (φ) given in equation (4.1). The black hole solutions are written in equation (3.3) and depend on the single parameter φ h . The conditions derived on W (φ) have the same interpretation as before, and the thermodynamic quantities are computed from equation (3.5). 17 In figure 8 the potential is plotted for α = 0.9 ∼ 1 and several values of λ in order to illustrate its key qualitative features. As before, segments of the curve in red, blue and green indicate whether the black hole solution with φ h = φ is non-existent, unstable or stable. For λ negative there is a single branch of stable black hole solutions, with the zero temperature solution determined by φ h = φ 0 (λ) > 0. The value of φ 0 (λ) changes smoothly with λ and there are no finite temperature phase transitions. For positive λ there are two separate branches with stable and unstable black holes. At λ = λ c > 0 there is a transition to a phase for which there is no zero temperature black hole solution. The critical value associated to this transition can be easily computed and is given by , where e is Euler's number. Note that this is a different kind of transition than the one analysed around equation (3.7). While in that case the transition was between two different zero temperature regimes, in this case there is a transition between a system with a T = 0 regime and one without.

Perturbative genus expansion
Consider now the matrix model computation of the partition function at leading (disc) order, as in equation (3.13). The spectral density computed in refs. [21,22] The top row of figure 9 displays the spectral density ρ 0 (E) (extracted numerically from (1.9) with (4.3)) for α = 1/4 and several values of λ. Strikingly, ρ 0 (E) develops negative regions for a finite window of λ approximately given by λ ∈ (0.035, 0.122). As discussed before, this clearly signals a problem, and again its origins lie in the behavior of u 0 (x) as a solution to equation (4.3). The place to look at is the implicit definition of u 0 (x) from the string equation (4.3). In the second row of figure 9 is plotted the curve x = x(u 0 ) for the same values of α and λ as before. From the first two diagrams starting from the left, we note that after λ 0.034, the threshold energy E 0 jumps discontinuously as a function of λ. This zeroth order phase transition generates a multi-valued potential u 0 (x) in the crucial region x < 0 as can be seen at e.g., λ = 0.045. This ultimately results (for the reasons discussed after equation (3.16)) in a negative spectral density. Amusingly, while a phase transition in E 0 solved the issues of ρ 0 (E) in the previous section, in this case it is the origin of the problem.
Since the potential u 0 (x) is multivalued in the region x < 0, the issue is more severe than the one discussed in the previous section. This is because leading order observables in the matrix model are computed from the behavior of u 0 (x) in x < 0. Recall for instance, the disc partition function Z 0 (β), given in terms of u 0 (x) in equation (2.10). That expression makes no sense if the function u 0 (x) is multivalued in the region x < 0, meaning the → 0 limit that defines Z 0 (β) does not exist. Since any single trace observable (like the spectral density) can be computed from Z 0 (β), it means none of these observables have a well defined leading genus behavior. This results in the breakdown of the perturbative expansion of the matrix model. 18 This type of failure is not extremely unusual, for instance observables JHEP04(2021)030 in a multi-cut Hermitian matrix model do not have the ordinary power series expansion in 1/N [50].
It must be emphasized that the multi-valuedness of u 0 (x) for x < 0 is the root cause of the physical problems here. Negativity of ρ 0 (E) is only a symptom. A more careful comparison of the diagrams in figure 9 illustrates this. Since for λ = 0.15 the leading genus spectral density is positive, one might think the matrix model perturbative expansion in this case is sensible. However, a closer inspection shows that u 0 (x) for λ = 0.15 is still multi-valued in the region x < 0, meaning the leading genus spectral density obtained from equation (2.10) is not defined at all. This shows the negativity of the spectral density is not a robust indicator when assessing the health of the matrix model.

Non-perturbative effects
The breakdown of the perturbative expansion calls for an appeal to non-perturbative physics. Here again, the definition of section 2.3 is employed. To do so, the coefficients t k appearing in R[u, x] in equation (2.7) must be obtained. These can be identified by expanding the disc level equation (4.3), giving: Using these coefficients we can write R[u, x] and the non-perturbative differential equation (2.15) for u(x). As discussed in the case of section 3.3, explicit solution can be obtained via a truncation to some highest k = k max as in equation ( to pick a truncation that captures the phase transition of E 0 shown in figure 9. Only odd truncations k max are useful in this regard, with k max = 5 the minimum truncation with the necessary structure, given the properties of the t k (λ). To obtain a more accurate spectral density for a larger range of energies, the truncation k max = 7 is used. The curve x = x(u 0 ) obtained from R (7) 0 [u 0 , x] = 0 is analogous to the ones shown in figure 9, with the transition occurring at λ c 0.03480. For k max = 7, the differential equation is 14th order but as before (see footnote 14) some persistence can yield highly accurate numerical solutions for u(x) (errors smaller than 10 −4 ).
Before the transition: first is the case of λ < λ c , where λ c is the critical value triggering the transition in E 0 . Similarly to the previous section solutions can be found for the nonperturbative potential u(x) and its associated spectral density ρ(E) computed. The results are shown in figure 10 for several values of λ < λ c . There is no phase transition for E 0 , and its value changes continuously with λ, as can be seen from the behavior at u(+∞) = E 0 . The spectral density ρ(E) does not show any particularly novel behavior for λ < λ c (in contrast to some of the phenomena seen in the previous section).
Across and beyond the transition: next is a study of the non-perturbative behavior of the model as λ crosses the transition value λ c where E 0 has a finite jump. In the k max = 7 truncation the range of λ for which the perturbative expansion breaks down is given by λ ∈ (0.03480, 0.36733), with the left edge corresponding to the zeroth order phase transition in E 0 .
As soon as λ crosses the transition the string equation (2.15) for u(x) becomes extremely unstable. While for λ arbitrarily close to λ c before the transition, well behaved solutions can be constructed, the numerics becomes unpredictable as soon as λ goes be-JHEP04(2021)030 yond λ c . 19 Before jumping to the conclusion that the deformed JT gravity theory is non-perturbatively unstable beyond the transition, it is prudent to carefully consider other alternatives that might explain the issue. The first obvious possibility is a lack of numerical precision: the discontinuous jump in E 0 might be making the differential equation harder to solve, even though a well behaved solution might in principle still exist. Deforming away from a neighboring solution already found (a common trick for numerically exploring difficult equations) does not work here since E 0 jumps. To explore this, a toy model was constructed that has the same features as the deformed JT gravity theory with the difference that the jump in E 0 can be controlled and made small. The model has the following leading genus string equation: where b ∈ R and a > 0. While for b = 0 this model has u 0 (0) = E 0 = 0, a zeroth order phase transition is triggered when b > 0, where the threshold energy jumps to E 0 a and u 0 (x) becomes multi-valued in the region x < 0. The plot of u 0 (x) for this toy model is analogous to the first two plots in the second row of figure 9. By tuning the parameter a, the discontinuity in E 0 across the phase transition can be made as small as needed. The full non-perturbative string equation (2.15) associated to this toy model is written after identifying the coefficients t k (a, b) from (4.5). Attempting to numerically solve the resulting equation for b > 0, we find the instabilities are still present, no matter how small the parameter a is taken to be. This shows the instability issue is not addressed by reducing the magnitude of the jump in E 0 .
There is an analytical argument that has been previously used in refs. [44,51] to argue in favor of the existence of solutions to the differential equation (2.15) (and variants of it),

JHEP04(2021)030
when the right hand side has the constant 2 Γ 2 instead of 0. It involves taking a t'Hooft limit as → 0 such that Γ becomes large while the combination q = Γ remains finite. In this limit, the string equation simplifies to the following algebraic constraint: can be easily solved numerically and visualized, see figure 11 for a plot for a fixed q and λ > λ c . The solutionû(x) automatically satisfies the appropriate boundary conditions, and for large q a smooth single-valued solution exists. However for small enough q the solution u(x) is multi-valued within a certain range in x (which interestingly occurs for x > 0). The emergence of the multi-valuedness (for small enough q) suggests that there is no path connecting possible smooth stable solutions of the full differential equation that may exist at finite Γ to the ones needed at Γ = 0. 20 More generally, it is difficult to see how a multivalued solution of an algebraic equation, our u 0 (x), (or even one smoothed out by the q-deformation above) can be completed into a solution of a boundary value problem for an ordinary differential equation once is turned on. The multi-valuedness overconstrains the information supplied by the boundary conditions. As a last remark, it is with noting that a multivalued classical potential u 0 (x) also makes no sense from the point of view of the quantum mechanical problem on the real line x ∈ R. At a given point x there would be an ambiguity as to the value of the force F x = −∂ x u 0 (x).
All the evidence therefore points towards a non-perturbative instability or nonexistence of the deformed JT gravity theory for λ across the transition, when u 0 (x) becomes multi-valued in the x < 0 regime. Since u 0 (x) is multi-valued in the finite window λ ∈ (0.03480, 0.36733), this means we should be able to construct numerical solutions for λ on the other side of this window, i.e. λ > 0.36733. Figure 12 shows that this is indeed the case, displaying regular solutions to u(x) for λ beyond the multi-valued window of u 0 (x). The associated spectral densities are also well-behaved and shown in that figure. Overall, there is enough evidence to reach the following conclusion about the leading u 0 (x) (defined, possibly piecewise, as a continuous function over the whole domain −∞ ≤ x ≤ +∞): For λ with u 0 (x) multivalued, the system is non-perturbatively unstable This is the central result of this section.
It is worth remarking that the above statement/result has a more general validity. In the case of λ = 0, ordinary JT gravity, the non-perturbative string equation from a Hermitian matrix model definition is just R = 0. There, the full u 0 (x) is multivalued in 20 Such a path can be found when there is no multi-valuedness, by for example using a solution-generating transformation [52,53] that acts on solutions of the string equation to generate new functions that are solutions of the string equation with Γ changed by an integer. It is in fact a kind of "Bäcklund" transformation, known in the integrable systems literature for changing the soliton number of certain kinds of solutions to integrable systems, here the relevant integrable system is KdV. The transformation can be made explicit, and an extension of it for cases with non-zero E0 described and derived in appendix C.  9). Here λ = 4 (purple), λ = 4.5 (red), λ = 5 (green) and λ = 5.5 (black). On the right the associated spectral densities obtained from (2.6) and plotted using the same colors. the x > 0 regime, and the statement about multi-valuedness above is consistent with the fact that the string equation does not seem to have a stable solution. The non-perturbative definition of ref. [25] (and its extension to this paper) is such that the piece of u 0 (x) in the x > 0 regime is replaced by u 0 = E 0 , removing the multi-valuedness, and the string equation (2.15), has a nice (unique) solution. It is interesting to speculate if there are alternative non-perturbative definitions that somehow repair the multi-valuedness in the x < 0 regime when it occurs, maintain perturbation theory, and yield sensible solutions to some new differential equation.

Final remarks
In this work we have investigated and resolved some puzzles regarding the double scaled matrix models shown in refs. [21,22] to describe certain deformations of JT gravity. In doing so, we have supplied a non-perturbative completion of the double scaled matrix model physics that extends the proposal for ordinary JT given in ref. [27]. Studying the deformations in detail, we have uncovered an interesting phase structure, in some special cases in qualitative agreement with the semi-classical analysis of the two-dimensional dilaton gravity. Let us mention a few open questions.
Beyond sharp defects: the matching between deformations of JT gravity and double scaled matrix models in refs. [21,22] was performed for potentials U (φ) in (1.3) with α i ∈ (0, 1/2). This corresponds to inserting sharp defects to the JT path integral. The issue with extending this range further is due to the fact that the decomposition of the hyperbolic surfaces in the topological expansion as done in ref. [14] breaks down for α i ∈ [1/2, 1). This technical difficulty has been very recently overcome in ref. [54]. 21 Building on earlier studies JHEP04(2021)030 of deformations of minimal models [55,56], the authors were able to compute the leading genus string equation for α i ∈ (0, 1). For a single defect (i.e. r = 1 in (1.3)) it can be written as [54] (5.1) While for α ∈ (0, 1/2] there is a single term in the summation and we recover (1.9) as given in refs. [21,22], this equation provides a very interesting and non-trivial generalization in the range α ∈ (1/2, 1).
This string equation supplies the necessary ingredients to extend our matrix model analysis from section 4 to α ∈ (0, 1). In particular, it raises the following question: are the non-positive spectral density ρ 0 (E) and multi-valued potential u 0 (x) observed in figure 9 still present in the regime α ∈ (1/2, 1)? A simple analysis shows the answer is affirmative. For instance, taking α = 3/4 the threshold energy E 0 obtained from the largest solution to R 0 [E 0 , 0] = 0 above, exhibits a zeroth order phase transition around λ ∼ 0.17. In a completely analogous way as in figure 9, the potential u 0 (x) becomes multi-valued in the region x < 0 and ρ 0 (E) non-positive. Overall, this suggests that while the disc string equation (5.1) becomes quite complicated in the extended range of α, our analysis in section 4 does not change that much. That being said, studying the equation (5.1) in further detail is an interesting question we would like to explore in future work.
Normalization of the dilaton potential: as explained in ref. [22], the computation of the Euclidean partition function from the classical definition of the action in (1.3) entails some ambiguities. These have to do with different choices of operator orderings and renormalization procedures in the quantum theory. This results in an ambiguity regarding the normalization of the potential U (φ), already seen by comparing equations (1.6) and (D.2) in refs. [22] and [21] respectively. In this work we have used the normalization of ref. [22]. The choice of normalization turns out being important, specially for cases in which U (φ) in (1.3) contains several terms, as it affects the behavior of U (φ = 0). In particular, it determines whether the phase transition at T = 0 discussed around (3.7) is present or not in the semi-classical analysis. It would be interesting to better understand the correct normalization of the dilaton potential (see section 5.1 in ref. [54] for progress in this direction).

A Inverting the Abel transform
In this appendix we invert the Abel transform and derive the simple formula for R 0 [u 0 , x] given in (2.12). Our starting point is (2.11) where in the second equality we have defined g(y) = 1/ √ y, written the integral as a convolution and the prime is a derivative with respect to u 0 . Applying the Laplace transform, defined in the usual way and using the convolution theorem we find where in the second equality we used L{f }(β) = β L{f }(β) − f (0) and solved the Laplace transform of g. From the leading genus string equation (2.9) we can evaluate the boundary term R 0 u 0 =0 = x, given that t 0 must vanish so that u 0 (0) = E 0 = 0. Solving for the Laplace transform of R 0 we find where we have again used the convolution theorem and used L{1}(β) = 1/β. Applying the inverse Laplace transform we obtain the explicit formula for R 0 in terms of ρ 0 (E) given in (2.12).

B Further details of the non-perturbative definition
Based on the work in refs. [40,42], in this appendix we motivate in a simple way the non-perturbative differential equation for u(x) in (2.15) using the assumptions in (2.19). Differentiating the scaling relation (2.19) with respect to s, setting s = 1 and using the KdV flow equation we find where a prime is a derivative with respect to x. In the second step we have used the recursion relation satisfied by R k [u] (2.8) and we have written everything in terms of R defined in (2.7).
Our additional condition is that whatever equation we obtain must perturbatively agree with the Hermitian matrix model, i.e., it must have the solution R = 0 in x < 0 perturbation theory. Given that E 0 scales the same way as u, one of the simplest options JHEP04(2021)030 is that ∂u/∂E 0 is proportional to R . Powers of R and other powers of R are in principle possible, which would need to be combined with powers of x to achieve the correct scaling. However, requiring that the boundary condition for u(x) at large positive x must be fixed as u(x → +∞) = E 0 (the perturbative continuity discussed in equation (2.17)), means that ∂u/∂E 0 must go to unity as x → +∞. Indeed, from the boundary condition u(x → +∞) = E 0 we have lim x→+∞ R = u (x) ∞ k=0 t k ku(x) k−1 + 1 = 1, where we used ∂ n x u(x) vanishes for n > 0 and large positive x. No other combinations of R and R (combined with powers of x, for scaling) can satisfy this additional constraint. Therefore, 22 Using this we can multiply both sides in (B.1) by R, giving a total derivative, then integrate and obtain: where we have fixed the integration constant to zero, so that R = 0 (the Hermitian model string equation) is a solution. We have arrived at the non-perturbative string equation (2.15) for u(x). Some interesting features of the physics are apparent from this derivation (some of which have been studied in ref. [57]). For example, equation (B.2) shows that as E 0 changes (as it will under some deformations) by some amount δE 0 , the corresponding change to the function in the asymptotic region (because R → 1 there) is δu = δE 0 , preserving the boundary condition as it should. More generally, equation (B.2) is a new non-perturbative flow equation (joining the KdV ones) for u as a function of E 0 . Recall also that the parameter E 0 (the end of the matrix model spectrum) defines "the wall". It arises from the potential that defines the probability measure, now as a function of (M M † + E 0 ) ≥ E 0 . In the limit of E 0 → −∞ the model is increasingly equivalent to an Hermitian matrix model (arbitrarily negative eigenvalues are now allowed), and indeed the dominant part of the equation (B.3) in that limit is simply E 0 R 2 = 0, which forces u(x) to solve the Hermitian matrix model equation R = 0.

C Deriving a Bäcklund transformation
In this appendix we derive a Bäcklund transformation that relates solutions to the differential equation (2.15) when the right-hand side has the constant 2 Γ 2 instead of zero. This is a slight generalization of the derivation in refs. [52,53], where this transformation was derived for E 0 = 0.
Let us start by assuming we have u(x), a solution to the differential equation (2.15). We can use this solution to define two functions v ± (x) according to While there are other more complicated scale invariant terms that could be written down that vanish to all orders in perturbation theory (such as e − 2 /R ) , it is not clear if any of them can be produced from a matrix model, as is the case for the one given in equation (B.2).

JHEP04(2021)030
Using X ± = 0 together with the fact u(x) satisfies the differential equation (2.15) (with 2 Γ 2 on the right-hand side) we have the following trivial identity that we can use to solve for u(x) in terms of v ± (x) This is the Miura transformation, mapping between the KdV and mKdV hierarchies, generalized to non-zero E 0 . The differential equations satisfied by the functions v ± (x) can be obtained combining (C.3) and (C.1), so that we find where we have defined This is the string equation in the mKdV hierarchy (arising from unitary [58] and Hermitian multi-cut [59] matrix models) generalized to non-zero E 0 . Let us introduce an additional notation, by adding a subscript to u(x) and v ± (x) that indicates the parameter that appears on the right hand side of their respective differential equations. The solutions to (2.15) and (C.4) in this notation are given by While this property was noted in refs. [52,58] for E 0 = 0 it is non-trivially extended to non-zero E 0 . For any arbitrary value of k > 1 it is straightforward to check it continues to hold. This turns out to be very useful, as it gives us another way of relating solutions to the differential equation in (C.4). Given any solution v c (x) (with constant c on the right hand side), we can use the reflection property satisfied by S k to generate a solution with opposite constant, i.e. v −c (x) = −v c (x). We now have all the ingredients necessary to derive the Bäcklund transformation for u Γ (x), shown below in (C.11). Let us prove it for (Γ − 1), the other case being completely analogous. The starting point is (C.3), which in the notation in (C.6) becomes The Bäcklund transformation is obtained after rewriting the first term in this expression using (C.1): This is a useful transformation in a number of situations. For example, sometimes it is easier to solve the string equation at another value of Γ (because for example, it has a shallower potential well), and then afterwards apply the transformation to change Γ. This was sometimes done in ref. [27] as a convenience (e.g., solving for Γ = 1 and then reducing to Γ = 0 afterwards), although the current paper did not need to employ this method since the numerical methods used have since become more powerful for solving directly at Γ = 0. It is also possible to imagine using this transformation successively to reduce Γ from some high integer value down to zero, having started with it large, or even formally infinite.
Finally, note that in the 't Hooft limit of section 4.3 where the solution is algebraic, it is easy to see that there is no path to smaller Γ using these transformations. Since here u Γ (x) isû(x), a solution to equation (4.6) valid in the assumed limit, the same limit should be used in the transformation. This means that the transformation becomes: where in the last step equation (4.6) was used to write q 2 /R 2 0 [û, x] =û − E 0 . In other words, the 't Hooft limit solutionû(x) is invariant under the Bäcklund transformation.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.