$T\overline{T}$-deformed free energy of the Airy model

Sharpening the holographic correspondence of Jackiw-Teitelboim (JT) gravity and its dual matrix model description at a finite radial cutoff $\lambda$ through the $T\overline{T}$ deformation is of interest. To proceed, we simplify the problem by considering the Airy model and deform Airy correlators in the same way as in $T\overline{T}$-deformed JT gravity. We use those correlators to compute the annealed and quenched free energies for both $\lambda>0$ and $\lambda<0$ from an integral representation of the replica trick. At the leading order in $\lambda$ and low temperatures, we confirm that the genus-zero quenched free energy monotonically decreases as a function of temperature when perturbation theory is valid. We then study the all-genus quenched free energy at low temperatures, where we discover and discuss subtleties due to non-perturbative effects in the Airy model, as well as the contributions from the non-perturbative branch under the $T\overline{T}$ deformation.


Introduction
An interesting question to ask is: does the AdS/CFT correspondence hold in a finite patch in spacetime? This question has thoroughly been investigated through the lens of the irrelevant T T deformation. Recently, the T T deformation [1][2][3] has attracted much attention due to being a well-defined operator with complete solvability along the RG flow. Historically, unlike relevant and marginal deformations, irrelevant deformations modify the UV behavior of a theory and are mathematically challenging to determine their UV behavior along the RG flow as generically infinitely many operators are included.
Given a seed action, for example Euclidean action S E (0), the T T deformation is defined by point-splitting up to total derivatives of local operators 1 where the right-hand-side depends on the determinant of the two-dimensional stress tensor T ij (λ) of the deformed seed theory, and T ij (λ) is related to its trace via An application of the T T deformation (1.1) relevant to this paper is probing the low temperature limit of JT gravity and its Airy model description. 2 We will now go over the salient features of JT gravity in the low temperature limit, and then discuss more of the T T deformation in holography.
JT gravity [6,7] is a special case of two-dimensional dilaton gravity, and on AdS 2 which is of concern in this paper 3 , it is described by the action in Euclidean signature with boundary terms where the dual description is a Hermitian random matrix model shown by Saad, Shenker and Stanford [12]. Here R fixed to be −2 is the Ricci scalar of the metric g µν , Φ is the dilaton, and S ∂M is the one-dimensional Gibbons-Hawking-York boundary action.
To further motivate the Airy model, we list a few reasons why it is interesting and useful to study. Firstly, the genus expansion can be summed via [13] allowing one to make definite statements on nonperturbative corrections. Secondly, JT gravity at low temperatures, or more precisely the 't Hooft limit: is dual to the Airy model at all genus known from [14]. Here the authors of [14] showed all the non-trivial information of the spectral curve of JT gravity is still preserved under 1 For a modern review of the T T deformation regarding applications to holographic systems as well as several other scenarios, see [4,5] and references therein. 2 Hereafter, we will refer to the Airy limit of Gaussian matrix models as the "Airy model". 3 On (nearly) dS2, "+2" in (1.3) will become "−2". The T T deformation in dS3 has been studied recently by [8][9][10][11] where the trace flow equation (1.2) is modified by an additive term ∝ 1 λ .
the 't Hooft limit (1.4). Additionally, there is an important caveat as [12] finds: the exact eigenvalue density has an exponential leakage when E < 0 which is denoted as the "classically forbidden" region making the system unstable. This non-perturbative instability is not special for the Airy model as the same phenomena occurs for the matrix model dual of JT gravity. However, a recent proposal by [15] improves the non-perturbative behavior of JT gravity by removing the non-perturbative instabilities. Therefore, one should only expect this relation between JT gravity and the Airy model confidently holds perturbatively. Finally, there are subtleties for the Airy model's quenched free energy not monotonically decreasing as function of temperature using directly the replica trick as done by [16]. Fortunately, Okuyama [17] showed this failure of monotonicity arises due to analytical continuation issues in the correlators Z n to Z n=0 and proposed an alternative formulation of the replica trick to correctly give a monotonically decreasing quenched free energy in the low temperature limit after summing over all genus. With these motivations of the Airy model, we wish to investigate how some of these features change under the T T deformation with a holographic picture in mind to understand JT gravity and its matrix model dual.
The first study of probing holographic systems under the T T deformation was initiated by [18] in AdS 3 gravity and inspired several follow-up investigations which computed partition functions and correlators of the gravitational and boundary field theory [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35], improving the T T -deformed holographic correspondence. More specifically, through the lens of the T T deformation, it is interesting to probe observables between JT gravity and its double-scaled matrix model dual description at a finite cutoff governed by λ, as done in [36][37][38][39][40][41]. The purpose of this work is to further sharpen the correspondence between JT gravity at low temperatures and its Airy model description by computing the correlators and quenched free energy. We now comment on the deformed energy spectrum in JT gravity.
Gross et al. [36,37] confirmed that the T T -deformed energy spectrum of JT gravity and its dual one-dimensional Schwarzian quantum mechanics match. 4 Here the 2D definition in (1.1) and (1.2) is dimensionally reduced via T τ φ = 0 to yield the following differential equation for the one-dimensional stress scalar T τ τ = f ± λ (E): 5 with two solutions where λ ∈ R and E is the undeformed energy. 4 JT gravity can be written in terms of a BF gauge theory. The deformed BF formalism and related supersymmetric quantum mechanics was found in [42,43]. 5 An alternative perspective on deriving this flow equation may be found through the Wheeler-de Witt equation derived in [23] for higher-dimensional large-N CFTs. For two-dimensions, results in [36] match (1.6) from the Wheeler-de Witt equation perspective. For family of dilaton gravity theories which are more general than JT gravity, e.g., [44], the authors of [45] derived the deformed energy spectrum at a finite cutoff.
As can be seen in the deformed energy spectrum (1.6), the sign of λ > 0 violates unitarity when the undeformed energy is E > 1/8λ. 6 The authors of [38] carefully dealt with this violation of unitarity and restored it in their rigorous non-perturbative treatment for the deformed partition function from a Wheeler-de Witt wavefunctional perspective. In short, one writes a linear combination of the wavefunctional of the two branches (1.6) such that the density of states stays real for all energies.
Other treatments of the deformation parameter λ's sign are addressed by [40] from a double-scaled matrix model perspective of JT gravity with an attempt to define the dual deformed matrix model description at a finite cutoff. Unfortunately, the analysis of [40] was unable to match the T T -deformed correlators between JT gravity and the dual matrix model, but did provide several alternative methods on how one could properly make this correspondence well-defined. We will discuss the importance of [38] and [40] more throughout various places of the paper. 7 Having laid out the motivations for studying the T T deformation in JT gravity and its dual matrix model, we will now summarize results found in this paper.

Summary of results and outline
It is helpful to point out in this paper, we will think about JT gravity and its T T deformation from the point of view of the matrix model description. This is because, as shown in [12], higher topology contributions in JT gravity are captured entirely by a double-scaled matrix model. The partition function of JT gravity at any genus and with any number of boundaries can be computed from the correlators in its dual matrix model. From the point of view of boundary theory, we only know that JT gravity on a disk is dual to Schwarzian theory on the boundary. Also, the result of [12] showed the connected n-point function in JT gravity -namely the partition function on a connected surface with n boundaries -does not factorize. The boundary theory was shown to be an ensemble of theories, not a single theory. Since the replica trick is essential to compute quenched free energy and requires knowledge on n-point correlators, we will think mostly from the point of view of matrix model and its correlators rather than directly from the boundary theory.
In §2.1, we first explain how the T T -deformed partition function of various topologies in JT gravity are found through an integral transformation of the undeformed partition function. Additionally in §2.1, we then extend this integral transformation to find a relation between the deformed and undeformed correlators via: 6 Throughout this paper, we frequently use the terminology "bad sign" for λ > 0 as the deformed energy spectrum is complex-valued when E > 1/8λ. The "good sign" corresponds to λ < 0 and the deformed energy spectrum is always real for all energy values. 7 A recent proposal in [46] successfully matched the T T -deformed partition functions of N = 1 type 0A and type 0B JT supergravity with the associated matrix models. Additional evidence of the duality from [46] was calculating the deformed T T -deformed matrix model correlators via topological recursion relations.
where ρ(E 1 , · · · , E n ) is defined by ensuring Only when the integration range is E i ∈ 1 8λ , ∞ and λ < 0, we can use the integration to be reviewed in §2.1 to conveniently compute the deformed correlators: This is the case when we look at any particular order in the genus expansion. However, when we look at the exact solution of the Airy model, non-perturbative effects there extend the integration range to E i ∈ (−∞, ∞), reflecting the translational invariance of the underlying effective Hamiltonian [15]. As we will later encounter the non-perturbative effect of the T T deformation, hereafter we refer to this as the non-perturbative instability to distinguish the two.
With the density of states for the Airy model, the undeformed one-point function is given by where, from [47,48], the density of states is and Ai(·) is the Airy function of the first kind [49], defined as In this Airy model, we should use the following instead of (1.9), which is invalid because the integral transformation (1.10) of e −βE diverges when E < 1 8λ . From determining the deformed correlators in §2.1, an immediate application is to calculate the annealed and quenched free energies in JT gravity at low temperatures where certain approximations are feasible. In §2.2, we review a new way of the computing the quenched free energy from an integral representation of replica trick due to Okuyama [17]: where e −Zx is fully determined by connected correlators Z n c . We elaborate later in §2.2 why we use this integral representation (1.15) instead of directly using the replica trick directly in the low temperature regime and we comment more on the non-perturbative contributions from the T T deformation in §2.3.
However, it is generally difficult to evaluate the correlators of JT gravity under the T T deformation let alone sum over all the connected correlators in order to use Okuyama's formula (1.15). Since JT gravity is known to have a matrix model dual, we will simplify the problem by studying a simpler matrix model, the Airy model, to make progress. Here, it is important to notice that although the T T deformation of random matrix models have been studied in [40], the correlators there do not match those of T T -deformed JT gravity studied in [38]. Since ultimately our goal is to understand JT gravity and its dual matrix model at a finite cutoff, we will not follow the T T deformation of the matrix model defined in [40]. Instead, we simply require the correlators of the Airy model to transform in the same way as the correlators in JT gravity do under the T T deformation, and we will take this as a working definition for our version of T T deformation applied to the Airy model. We will not attempt to completely explore this version of deformation for matrix models and hope to revisit this problem in future work.
It is also important to point out one caveat of our simplification. Naively, the doublescaled matrix model dual to JT gravity shares the same non-perturbative instability as the Airy model we considered; for instance, see [12,15]. However, there have been studies on how to improve the non-perturbative behavior of JT gravity and remove this undesired feature [15]. Therefore, it would be interesting to extend our work to JT gravity at general temperatures with the improved non-perturbative behavior.
In §3, we use (1.15) to numerically evaluate the quenched free energy in the Airy model for both λ > 0 and λ < 0. As a warm-up, we first compute the quenched free energy F q,λ (T ) at genus-zero at the leading order of λ in perturbation theory. We confirm that F q,λ (T ) is a monotonic function in T at low temperature when a leading order approximation is valid. Additionally, we find that the good sign λ < 0 deformation decreases the quenched free energy while the bad sign λ > 0 increases it. Intuitively, this sign of the deformation parameter λ > 0 corresponds to JT gravity in a finite box with Dirichlet boundary conditions at r c = πλ 4G . The increase of the quenched free energy for the λ > 0 theory is related to the fact that the T T deformation cuts off the spectrum. On the gravity side, one can think of this phenomena as a gravitational redshift. For an object in a gravitational potential, the energy measured at the conformal boundary (which is the undeformed case) is reduced compared to the energy measured at the particle location. For example, in the extreme case where a particle is sent towards a black hole's horizon, the energy at infinity vanishes and is negative when the particle is inside the event horizon. Here, in the cutoff gravitational theory, we are measuring this local energy at the finite cutoff boundary, which is closer to the particle compared to the conformal boundary. Thus the amount the redshift decreases and the energy we measure increases. 8 Next, we use the low temperature approximation Z(β) n c Z(nβ) of the Airy 8 We thank Per Kraus for explaining this gravitational interpretation.
model to compute the quenched free energy for finite λ. Notice that this computation include contributions from all genus as well as non-perturbative effects in the Airy model, which makes ρ Airy (E) = 0 even for E < 0.
For the bad sign λ > 0, as studied in [38,41], there are contributions from the nonperturbative branch. We compute the quenched free energy with this branch included and excluded. In both cases, we find the quenched free energy to be divergent and a careful analytical examination can be found in the Appendix A.
For the good sign of λ < 0, the issue of the complex-valued energy always arises regardless of the value of λ since ρ Airy (E) has support on the entire real axis. Although the issue of complex energy of the good sign λ < 0 has been noticed before in [1] for 2d CFT since the ground state energy is − c 12 , we emphasize this is different from our case, as the energy spectrum is still bounded below by − c 12 . One can avoid the complex energy by simply choosing |λ| ≤ 3 2c in this case. However, for the Airy model, the issue of complex energy will always arise regardless of the value of λ. This has not been noticed before for the good sign of the T T deformation. A simple solution will be to impose a hard cut-off in E, namely we simply remove these states with complex-valued energy. However, as we will see in §3.4, this option would lead to a violation of the T T flow equation by a boundary term in the integral. Another option would be to include those states with complex-valued energy, however, to demand the partition function to be real-valued, one must include the other branch as well. By a careful choice of coefficient, one can make sure that the boundary terms cancel each other properly so the flow equation is satisfied. We then compute the quenched free energy for both cases: if we exclude the non-perturbative contribution, we find the quenched free energy to be finite and monotonic at low temperatures, and it decreases under the deformation. If we include the non-perturbative contribution, we find the quenched free energy diverges, and a careful analytical dissection is included in Appendix B.

T T deformed correlation functions in JT gravity
Here we review how one computes deformed partition functions in JT gravity via an integral transformation on the partition function in the undeformed theory. As already alluded to in the introduction, the deformed energy spectrum is given by (1.6) and given the density of states, one can immediately compute the deformed partition function systematically as For the moment, we will only consider the contribution from the perturbative branch f − λ (E). We will return to the potential contribution from the non-perturbative branch f + λ (E) later in §2.3.
The partition function Z λ (β) satisfies a flow equation derived by [38]: which is closely related to the usual inviscid Burgers' equation from the 2D T T -deformed energy spectra (1.5) (e.g., see [4]). The deformed partition function (2.1) may be written in terms of an integral transformation involving a kernel and the undeformed partition function [36] when ρ E < 1 8λ = 0: where K(β, β ) is a kernel determined from the deformed energy spectrum (1.6): namely the inverse Laplace transform of the Boltzmann factor after deformation for λ < 0.
With the integral transform (2.3) and its kernel (2.4) at hand, one can proceed to compute the partition function of the deformed JT gravity on disk, trumpet and other topologies for λ < 0. For example, the map between the undeformed and deformed disk and trumpet partition functions respectively are where we have used an identity for the modified Bessel functions of the second kind when λ < 0 and 8aλ + β 2 > 0. We comment on one of the convergence conditions 8aλ + β 2 > 0 in (2.6) implies that the disk partition function (where m = 3/2, a = π 2 ) is well-defined up to the Hagedorn temperature T H = (−8π 2 λ) −1/2 as shown by [36]. Fortunately, for the trumpet partition function, there is no Hagedorn temperature (which takes place when m = 1/2, a = −b 2 /4). As pointed out in [38] for the bad sign λ > 0, the trumpet partition function diverges when the total proper length of the boundary equals the geodesic length. This divergence is removed by taking into account the non-perturbative branch thus making the trumpet's density of states real and converts the modified Bessel function of the second kind K 1 (·) into the first kind I 1 (·) as done by [38]. Likewise from [38], the same logic schematically holds for the disk partition function thus converting the modified Bessel for the second kind K 2 (·) into I 2 (·). We refer the reader to equations (4.7)-(4.8) and Appendix B in [38] for more elaborate discussions.
As in [38], knowing Z λ (β) D and Z λ (b, β) T allows us to build general correlation functions in the deformed JT gravity. The connected correlators on a hyperbolic Riemann surface with n boundary components and genus g are where S 0 is the two-dimensional Einstein-Hilbert action, and Z g,n;λ (β) are defined from [41] as follows: with the Weil-Petersson volume V g,n (b 1 , . . . , b n ) of a Riemann surface Σ g,n (i.e., with genus g and n distinct marked points p i ) defined in [50] as 9 Here M g,n is the Deligne-Mumford compactification of the moduli space M g,n of Σ g,n of complex dimension (3g − 3 + n), ψ i ≡ c 1 (L i ) is the first Chern class 10 of the tautological line bundle L i over M g,n whose fiber at the point (C, x 1 , . . . , x n ) ∈ M g,n is the cotangent line to the curve C at x i , ω is the Weil-Petersson symplectic form on M g,n , and b i is the length of the ith geodesic boundary component of Σ g,n .
From the definition in (2.8), one might wonder if the Weil-Petersson volumes should flow under the T T deformation? One possible way to see why Weil-Petersson volumes do not flow under the deformation 11 is that the flow equation (2.2) should be satisfied on each asymptotic boundary component with proper legnth β i . But by definition, the flow equation only contains derivatives with respect to λ and β, not b i , the length of geodesic boundary component to be glued together. This fact is also adopted in [40,41].
Then for generic Z g,n;λ (β), we can write the deformed partition functions as A few of Z g,n;λ (β) has already been computed in [41]. Also, see [12,51] for a review of Weil-Petersson volumes in 2D topological gravity and matrix models. 10 ψi are also called "ψ-classes", "Witten classes" or "gravitational descendants". 11 However, in the context of topological recursion, both the resolvent R g,n;λ and function W g,n;λ is deformed [41], while their relation to each other W g,n;λ (z1, · · · , zn) ≡ (−2) n z1 · · · znR g,n;λ (−z 2 1 , · · · , −z 2 n ) is intact. So in the deformed theory, Vg,n is no longer the Laplace transform of W g,n;λ . The topological recursion formula in terms of W g,n;λ [52] is covariant under the deformation, and retains the same form.
To conclude this section, we derive a differential operator presentation of the T T deformation similar to [37], which is sufficient for computing perturbative expansions in λ and will be used in §3.1, by rewriting the exponential where we have used . (2.14) It is then straightforward to compute partition functions of generic topologies using this differential operator. However, for multiple boundary components, one starts from the undeformed partition function with different inverse temperatures β i on each boundary component, and then applies the differential operator D y i ;λ | y i =β i to each boundary component separately We will use this differential operator presentation to perform perturbation calculation in the leading order of λ in §3.1.

Quenched free energy
In this subsection, we quickly review a recent novel way of performing the replica trick 12 following [17]. The replica trick (2.17) has been shown in [17] to be written as a rather convenient integral representation such that the analytical continuation from Z n to Z n=0 remains unambiguous. From (2.18), the first term is the annealed free energy while the second term encodes the contribution from Euclidean replica wormholes with the interpretation that the operator e −Zx creates spacetime boundary components, so-called "spacetime D-brane" or "SD-brane" introduced in the context of baby universes [53]. Additionally, the term e −Zx can be rewritten as e −Z(x) in terms of the following generating function of connected correlators This will turn out to be an important quantity when we compute the T T -deformed annealed and quenched free energies of the Airy model in the upcoming section and the appendices. A motivation to study a simple observable, such as the free energy, is to see how do Euclidean replica wormholes contribute to the Euclidean gravitational path integral (2.20) with spacetime boundary B, metric measure [dg] and JT gravity action S[g]. An easy way to determine the presence of Euclidean wormholes is to see whether correlation functions of the partition function cease to factorize among n boundary components. It turns out not to be the case due to the factorization failure as shown in (2.18), and this fact can be confirmed by directly computing the annealed and quenched free energies 13 as done in [16]: at inverse temperature β. They are shown not to be the same indeed, clearly indicating the factorization failure already hinted at by (2.18). Unfortunately, the authors of [16] computed F q (β) with direct usage of the replica trick (2.17) and their analysis at low temperature found that F q (β) is not monotonically decreasing as a function of temperature. This is fundamentally due to the non-uniqueness of analytically continuing Z n to Z n=0 . Given this conundrum at low temperature, the correct analytical continuation was performed by the author of [17] without directly using the replica trick (2.17), and there F q (β) was shown to be a monotonically decreasing function of temperature. We now review how this is done following the proof in [17]. The correlator Z n can be expanded in terms of connected correlators Z k c as [59]: so now the analytic continuation of Z n is unambiguous due to its being rewritten as a polynomial in n up to an overall Z n . To further find the integral representation (2.18), we generalize the above expansion to where i ≥ 2 and j k 's constitute an integer partition n k=1 kj k = n . (2.25) Now, using the standard prescription for analytical continuation and then the identity the quenched free energy F q (β) can now be written as the integral representation (2.18), whose integration range inherits that of (2.27).

Comments on non-perturbative contributions
In this subsection, we clarify more of the non-perturbative features arising from the T T deformation appearing in this paper. 14 The perturbative branch is denoted by the negative branch in the energy spectrum (1.6) since lim λ→0 f − λ (E) = E. In contrast, the λ → 0 limit of the positive branch for (1.6) diverges so, as expected, most papers omit this branch in their perturbative analysis and only consider the negative branch. We plot the entire deformed energy spectrum in Figure 1 to show both perturbative and non-perturbative branches. Unfortunately, when λ > 0, the spectrum along the flow becomes complexvalued for large enough energies. To resolve this issue, as explicitly shown in [38], one is forced to include the non-perturbative contribution such that the partition function is real with appropriate constraints on the density of states ρ ± (E) for JT gravity and is a solution to the flow equation (2.2). The appropriate constraints on ρ ± (E) are explained more in-depth by [38].
14 Another non-perturbative effects we will encounter is in the Airy model, whose density of states ρAiry(E < 0) = 0. This is already present in the undeformed theory and is unrelated to [41]. Yet another non-perturbative effect takes place upon summing over all genus. We believe the non-perturbative effects in summing over genus results in the non-perturbative instability. We will often refer this effect as the non-perturbative instability in order to distinguish the two. But based on context, the readers should have not trouble in distinguishing the two. Figure 1: We plot the entire deformed energy spectrum as a function of the undeformed energy E when |λ| = 1 and |λ| = 0.01. Here f + λ>0 (E) (blue curve) and f + λ<0 (E) (green curve) together describe the nonperturbative branch. While f − λ>0 (E) (orange curve) and f − λ<0 (E) (red curve) describe the perturbative branch. For either λ > 0 or λ < 0, we find a smoothly connected curve. We also see that when |λ| decreases, the non-perturbative branch becomes further away from the origin.
An alternative approach to naturally incorporate non-perturbative effects is through a resurgent analysis. In [41], the disk and trumpet deformed JT partition functions written as power series in λ are Borel resummed to obtain non-perturbative results, which are used to further study how the partition functions summed over topologies (i.e., topological recursion) and spectral form factor are modified under the T T deformation. Alas, despite these technical advances of the T T deformed JT gravity, the resurgent analysis still contains a spectral density not positive-definite.

Quenched free energy for deformed Airy model
In this section, we compute the T T -deformed annealed and quenched free energies in the Airy model. It is important to point out that our deformation of the double-scaled matrix model is different from the ones considered in [40], whose T T deformation of the doublescaled matrix model dual to JT gravity does not completely match the correlators with the T T -deformed JT gravity.
Our ultimate goal is to understand the quenched free energy in T T -deformed JT gravity. However, this is a rather difficult problem since even for the one-point function Z(β) JT of the undeformed theory, there is no analytical expression which includes all-genus contributions, let alone non-perturbative effects. In contrast, there has been much progress with numerical calculations [15,[60][61][62][63][64]. Therefore, we want to make a simplification and study the Airy model instead, which is known to be the low energy approximation of JT gravity. Since it is a double-scaled matrix model, one could apply the T T deformation defined in [40]. However, we already know that the T T deformation defined and studied there does not provide an exact match between JT gravity and its matrix model dual in terms of general correlators. Hence, this is not a good choice if our ultimate goal is to understand the quenched free energy in T T -deformed JT gravity.
In order to investigate the Airy model with the hope of retaining some essential features of the T T -deformed JT gravity, we will have to work with a different deformation from the one considered in [40]. Since one can construct T T -deformed JT correlators with any number of boundary components and genera from basic observables like the deformed disk and trumpet partition functions by gluing together pants decomposition [38], we can recast the T T deformation of the correlators of JT gravity in various ways (e.g., using the differential operator D y;λ | y=β ) as reviewed in §2.1. Now in order to match the T T -deformed matrix model dual to JT gravity, correlators on both sides must be deformed in the same fashion. We can then apply the same recipe of deforming the matrix model correlators to the Airy model, instead of using the one defined in [40]. This is the deformation we will adopt and study in this paper.
To be more specific, we will take the deformed n-point functions in the Airy model as It is important to notice that if we are interested in computing the deformation at any given genus, there will not be non-perturbative effects (from the Airy model itself before T T deforming) so the integration range of E i is [0, ∞), and we can safely apply the integral transformation (2.3) on (3.2) to obtain (3.1) for convenience. However, if we are interested in the exact result, for instance, Z(β) Airy,λ , then we cannot use the integral transformation, no matter if λ is positive or negative, because even for λ < 0, ρ Airy has support over the entire real axis of E.
There are a few other important caveats we need to mention. First, carrying over the T T deformation of JT gravity correlators to matrix model correlators only determines the perturbative branch of the T T deformation. There will be non-perturbative contributions of the T T deformation that need to be addressed if we are interested in results with finite, not infinitesimal, λ. We will analyze those non-perturbative effects case-by-case as we encounter them by matching the deformed Airy correlators with either the genus-zero result or the flow equation. A possible systematical treatment would be adapting the resurgent analysis performed in [41]. However, in our study, eventually we will have to sum over not only all orders in λ, but all genus g as well. As hinted in the introduction of this paper, it is wellknown that there is non-perturbative instability in the exact spectral density in the Airy model related to the genus expansion. Hence, we expect this "double" resurgent analysis in λ and g (both variables on the same footing) to be rather complicated and will leave it for future work.
Second, related to the non-perturbative instability of the Airy model, the same nonperturbative instability appears in the naïve matrix model dual of the JT gravity as well [12,15]. There have been works on how to improve the non-perturbative feature of the JT gravity and remove this undesired feature [12,15,65]. Hence, there are possibilities that the non-perturbative instability can qualitatively affect the behavior of the quenched free energy in the deformed theory and would be interesting as well as important to understand if this is the case. A possible way to figure this out is to extend our study to the JT gravity at general temperatures with the non-perturbative instability taken care of.
Third, it should be emphasized that we do not provide a complete description of our version of the T T deformation for the matrix model. We simply require that in the deformed theory, every correlator has to transform in a way in order to match the gravity side. One can take this as the working definition of our version of the T T deformation for matrix models. 15 To further illustrate the difference between our version and the one in [40], let us start with a Hermitian matrix integral: In [40], it is assumed that the deformed matrix model still takes the form of (3.3), and the only physical quantity that changes is the matrix model potential which shifts from As a somewhat unfortunate consequence of this assumption, the n-point functions in the deformed matrix model still do not match the n-point functions in the deformed JT gravity. While in our case, we start by requiring that the n-point functions in the deformed matrix model match the results in deformed JT gravity leading to a different deformation on the matrix model side. In particular, naively the integration measure [dM ] of matrices will receive corrections from the T T deformation as well if we want to keep the potential Tr V (M ) as a single trace operator as implicitly assumed in [40]. To see this, consider the correlators (3.5) in the undeformed matrix model. To match the gravity result, its deformation should assume the following form: To find the deformed matrix integral, we consider a change of variables for To further illustrate this fact, we consider the change of variables y i = f − λ (x i ) as well as neglect for the moment issues with the branch cut from the square root in f − λ (x i ) and potential non-perturbative subtleties related to the change of variable for M . 16 Then  Either way, the deformation presented here violates the implicit assumptions in [40] and would be the starting point on the matrix model side. However, it is unclear whether one should first T T deform and take the double-scaling limit or vice versa. Further analysis is beyond the scope of this paper, and we will leave the details of the T T -deformed matrix model to a future study. 17

Genus-zero quenched free energy at leading order of λ
We begin our analysis on the deformed quenched free energy F q,λ (T ) at genus-zero perturbatively. To determine F q,λ (T ), we first compute the deformed genus-zero multi-boundary connected correlators Z(β 1 ) · · · Z(β n ) g=0 c,λ . The undeformed connected correlators at genuszero have been computed in [14] using the genus-zero Korteweg-De Vries (KdV) flow: where g s ≡ √ 2 is the genus-counting parameter. 18 The deformed correlators can be computed directly using the integral transformation (2.3) and, by construction, solve the flow equation (2.2). For instance, for the good sign λ < 0, the deformed one-point, two-point, and n-point correlators are given by (3.12) Despite the n = 2 integral below can be easily integrated numerically, to our best knowledge, we do not know how to evaluate it in closed form; the β 2 -integral contains a piece not contained in standard tables of integrals, such as [68]. A possible way to evaluate the n = 2 integral is by writing the Taylor series expansion with the appropriate analytic coninuation such that the right-hand-side of (3.14) is welldefined. Therefore, we see a similar pattern with modified Bessel functions of the second kind (3.15) Generally, how one evaluates the above sums of products of modified Bessel functions 18 Here in fact gs = √ 2 and we will set to unity. Adopting the conventions by [17], is the genuscounting parameter in the Airy limit of matrix models while gs is the natural genus-counting parameter in 2D topological gravity. Also refer to earlier discussions in [14,67].
in (2.19) as in [17]. Therefore, instead, we consider working with perturbation theory and just keep the leading order in λ. For this purpose, it is convenient (in comparison with our "working definition" introduced in the beginning of §3) to express the T T deformation in terms of differential operators and only keep the leading order term in λ so that the deformation acts as where c 1 = 2 from (2.14). Then, we find On the other hand, can be evaluated using the similar trick as in [17]. Specifically, focusing on the leading order piece in λ, we have the following sum: Let z = x β 3 π , and the sum can be decomposed into three separate pieces: satisfying where W (z) is the Lambert function defined by the following Taylor series expansion: (3.23) The above differential equations (3.21) can be solved by making the ansatz 24) and using the property of the Lambert function . (3.25) We then find A(W (z)) − 10B(W (z)) + 12C(W (z)) .
where we have changed the integration variable from z to W (z) and used Notice that in the last line of (3.28), we expanded in λ for the logarithm and exponential since our result is only valid in the leading order of λ. This integral (3.28) can be evaluated numerically and below we plot the quenched free energy F q (β) at genus-zero with T T deformation coupling λ = −1/20 against temperature T . We also plot the genus-zero quenched free energy for various different signs of λ, again using differential operators as in (3.17), see Figure 4. Notice that in the perturbative expansion of λ, λ also always appears as λT . Hence, for the leading order approximation to hold, λT must be small. In the numerical calculation, we find F q,λ (T ) monotonically decreases as a function of T for the good sign λ > 0. However, monotonicity can break down when λT is too large for the bad sign λ > 0 and this is likely due to λT exceeding the range of validity of leading order approximation in λ.

All genus quenched free energy in low temperature limit
Now using the low-temperature approximation as in [17], we compute the quenched free energy starting from the relation in the undeformed theory [17,59]: Figure 4: We plot the quenched free energy F q,λ (T ) for the deformed Airy model as a function of temperature T . We notice that a deformation with the good sign λ < 0 lowers the quenched free energy while the one with the bad sign λ > 0 increases F q (T ). For both signs of λ, F q,λ (T ) monotonically decreases as T increases for T < 1; however, for the good sign of λ, this breaks down for λ = 1/20 and at around T = 1. Notice that the perturbative expansion of λ always appears as λT . In order for the first order approximation to hold, λT should be small. Therefore, this breakdown of monotonicity when T becomes large is likely due to the fact that we go beyond the validity of perturbation theory.
Under this low temperature approximation, the deformation of Z(β) n c Z(nβ) is easily computable. This is seen from expressing the undeformed correlator for n boundary components the term dependent on β essentially factorizes. Thus, the deformed n-point function in low temperatures is However, we must show that the change from the T T deformation for Z(nβ) c,λ − Z(nβ) c,0 is of lower order compared to the correction Z(β) n c − Z(nβ) to the approximation (3.30). This is the case because the correction to the approximation Z(β) n Z(nβ) is exponentially suppressed as e −c 0 β 3 in the low temperature limit β → ∞, where c 0 is some positive constant. One can explicitly check this approximation for n = 2, 3, where the exact expression of the partition functions can be conveniently found in [13]. 19 For instance, with n = 2, we have Using the asymptotic expansion of the error function: it is clear that the correction is suppressed by e −β 3 /2 in the low temperature limit for (3.33). Similarly, for n = 3: From the definition of Owen's T function: one can see in the large h limit (i.e., large β limit), it is indeed suppressed as e − 1 2 h 2 = e − 3 2 β 3 . For instance, we can manipulate the integral (3.36) as (3.37) One can also confirm this from numerical integration. In the following plot, we can see clearly, the correction is suppressed by e − 3 2 β 3 in the low temperature limit. The connected correlator has a well-known integral representation 20 given by [13] and is known to be a closed form only for n = 1, 2, 3. We can easily check that Z(β) n c Z(nβ) when n = 1, 2, 3 for small temperatures, but proving this for n > 3 is numerically difficult. We will content ourselves and assume Z(β) n c Z(nβ) is true for all n. Meanwhile, even at the leading order of λ, the T T deformation will give polynomial corrections in β for low temperatures. Hence, the corrections in the approximation Z(β) n 0 Z(nβ) 0 can still be neglected even when we consider the T T deformation. 20 Essentially a cyclic linear combination of n! integrals, each of which is an n-dimensional Laplace transform of the product of n so-called Airy kernels [70][71][72]  There is one caveat, however, the all-genus density of states does not vanish for E < 0. Instead, it is exponentially suppressed when E < 0: (3.39) This means for both good and bad signs of λ, there will be states with complex-valued energy. We discuss the two cases separately and in each case, there will be a plausible non-perturbative contribution from the f + λ (E) branch defined in (1.5). The deformed quenched free energy is defined and computed as

Bad sign: λ > 0
In this subsection, we study the quenched free energy of the deformed Airy model with λ > 0. In the deformed theory's spectrum, as explained in the introduction, the deformed energy f − λ E > 1 8λ is complex-valued. Therefore, a cut-off in E at 1 8λ is needed to remove the complex-valued energy states. Furthermore, there could be contributions arising from non-perturbative states with energy given by the other branch f + λ (E). These kind of nonperturbative effects have been rigorously studied in [38,41]. Their analyses lead us to conjecture what the non-perturbative contribution could be in the deformed Airy theory for the bad sign λ > 0. We will study the quenched free energy with and without the nonperturbative contributions. In both cases, we find the quenched free energy F q,λ (T ) diverges for every λ > 0 and T in our low temperature approximation. It is unclear to us whether the low temperature approximation which causes F q,λ (T ) to diverge or this is the feature of the deformation itself. Additionally, we do not know if repeating the same calculation in JT gravity will lead to the same problem. We will leave this for future investigations.
For the bad sign λ > 0, the issue with complex-valued energy is fixed by imposing a cutoff in the energy E such that the deformed energy spectrum is real-valued. Furthermore, the resurgent analysis in §4.2 of [41] indicates that there could be non-perturbative contributions to the partition function coming from the other branch f + λ (E). More specifically, the deformed genus-zero partition function of JT gravity is given by is just the usual JT gravity density of states for the disk. After a change of variables for the two terms in (3.41) respectively, they combine into a single integral (see for example [41]) where the deformed density of states is (3.45) Motivated by this, and following our working definition of T T deformation in the beginning of this section, it is natural to expect that for the Airy model deformed by λ > 0, we should have the one-point function upon an identical change of variables as before 47) and the second term in the first line signals non-perturbative effects. Note that here, instead of using the integration transformation for λ < 0, we must resort to our working definition introduced at the beginning of this section as the new prescription for the T T deformation which agrees with (2.28).
One can use this result together with the approximation Z(β) n Z(nβ) to compute the quenched free energy. However, in this case, we find that the quenched free energy actually diverges. To see this, notice that the difference between the quenched free energy F q (β) and the annealed free energy F a (β) is given by .

(3.49)
Hence, in order for the x-integral to converge, the difference D(x) ≡ e −Z λ (x) − e −x Z λ must at least go to zero as x → ∞ but at least faster than 1/ ln x. However, as we will show in the Appendix A, this is not the case for any λ > 0 and T . Here, we numerically plot D(x) for λ = 1/15 and T = 1/12. As one can see, the function D(x) is monotonically decreasing with x at the beginning, but turns around and starts monotonically increasing at a very large value x ∼ 2 × 10 19 . One might wonder if the non-perturbative branch with the energy f + λ (E) causes the integral to diverge. We can certainly only include perturbative branch when computing the quenched free energy. However, as we show further justification in Appendix A, the integral still diverges. One can see this divergence numerically from D(x) asymptotic to some finite number in the x → ∞ limit. Though not as dramatic as Figure 6, as one can see D(x) still does not go to zero as x → ∞.
It is unclear to us if such divergence is intrinsically physical, or due to any of our approximations. There are several possibilities. For instance, Okuyama's formula (2.18) may fail for the deformed theory in general. The derivation of (2.18) in [17] requires one to exchange the integral with an infinite sum which is not absolutely convergent. This may lead to the failure of (2.18) in the deformed theory. Another possibility could be that the divergence is due to the non-perturbative instability of the Airy model. A reliable way to rule out some of these possibilities is to extend our work to JT gravity with the proper improvement of its non-perturbative behavior. We wish to investigate this in JT gravity per se in the future work to see if this divergence still persists.

Good sign: λ < 0
In this subsection, we compute the quenched free energy for the good sign λ < 0. Since non-perturbatively the density of states of Airy model extends to E = −∞, the deformed theory will have unitarity issues caused by complex-valued energy which seems not to be discussed before in the previous literature of deformed JT gravity. All the densities of states in the literature have lower bounds, i.e., ρ(E ≤ E 0 ) = 0. Thus, by considering ρ(E) = ρ(E − E 0 ), we can also have a well-defined spectrum for all λ < 0. This is why λ < 0 is referred as the good sign in the literature 21 [73]. Here, the non-perturbative effect makes ρ(E) = 0 for all E ∈ R. Therefore, the deformed energy spectrum will be complexvalued for E < 1 8λ . There are two possible treatments. The first option is to impose a cut-off in the deformed energy spectrum up to when it becomes complex-valued. This cutoff resolves the unitarity issues, however, this leads to a violation of the flow equation (2.2) of Z(β) λ as a boundary term arises at the cut-off E = 1 8λ . Alternatively, we may include these states with complex-valued energy, but we must then include their corresponding states from the non-perturbative sectors as well as to ensure that the partition function is real. We will refer this part as non-perturbative. By carefully choosing coefficients, the contributions in the boundary piece of this term cancels the boundary term (3.55) in the first option. Thus, the flow equation of the one-point function Z(β) λ will be satisfied in this case. We will study the quenched energy with and without the non-perturbative contribution. Excluding the non-perturbative contribution, we numerically confirm the quenched free energy is monotonically decreasing with temperature T at a given λ < 0. We also find the quenched free energy monotonically decreases as we increase the absolute value of λ. Including the non-perturbative branch, unfortunately, we find that the quenched free energy computed from Okuyama's formula (2.18) diverges in general and we will illustrate this subtlety numerically in this subsection and include an analytical analysis in Appendix B.
We start with discussing in detail how to treat the complex-valued energy states in the deformed spectrum. One may expect that the correct answer is given by the exact recipe for the bad sign λ > 0, i.e., we cut off the spectrum below E < 1 8λ where the deformed energy becomes complex-valued and include the other branch for the remaining spectrum. However, there are two objections. The first objection is that if we consider the genus expansion, then the spectrum at genus-zero does not extend to E < 0. As a result, the deformed spectrum remains unchanged and the good sign of the T T deformation is welldefined making no additional branch required. If we included the other branch, through the genus expansion, the genus-zero partition function will receive corrections from the other branch as well which lead to inconsistencies. The second objection is that the Boltzmann weight e −βf + λ (E) diverges as e β √ E/2|λ| when E → ∞. Hence, the contribution from the other branch is divergent. These two reasons suggest we should not include the contribution from the other branch for the real-valued energy region. Thus, one might conclude the deformed partition function for the good sign λ < 0 is simply given by truncating the spectrum with complex-valued energy: where However, there is one caveat. The deformed partition function should satisfy the differential equation (2.2) derived in [38]: For convenience, we introduce the differential operator F ≡ 4λ∂ λ ∂ β + 2β∂ 2 β − 4λ β − 1 ∂ λ . Consider a change of variables E =Ẽ + 1 8λ so that the bound of the integral does not depend on λ: As one can show, F acting on the integrand leads to a total derivative (3.54) Therefore, (3.55) Now we see the problem: our guess Z(β) λ,guess violates the flow equation (3.52) due to the appearance of the boundary term (3.55). This is just another manifestation of the nonperturbative effect of the Airy model. If ρ Airy (E) had been supported on [E 0 , ∞), we can shift the ground state energy such thatρ Airy (E) = ρ Airy (E − E 0 ) to remove the complexvalued energy and makeρ Airy ( 1 8λ ) = 0 such that the flow equation is satisfied. However, this is not possible due to the non-perturbative effects as ρ Airy (E) has support on the entire real axis.
Thus, to make sure that Z(β) λ,guess satisfies the flow equation (3.52) while keeping it finite, we must include the complex-valued energy region where E ∈ (−∞, 1/8λ) to cancel the unwanted boundary term. To make sure the deformed partition function is real, we must also add its complex conjugate, i.e., the contribution from the other branch. Therefore, the partition function Z(β) λ for λ < 0 should be given by (3.56) where the sum of exponentials in the second integrand can be rewritten in terms of cosine as and is a highly oscillatory integral when λ or T is small, but can still fairly easily be numerically evaluated to be finite. As one can check, the boundary terms indeed cancel with each other and the flow equation (3.52) is satisfied.
One can then use Okuyama's formula (2.18) to numerically compute the quenched free energy. In this case, we can express Z λ (x) using (3.56) as the following: .
(3.58) Since the second integral can be rewritten as we will refer to it as Z cos λ (x). 22 , and the part I(E, x, λ, T ) will be the key ingredient in Appendix B.
We can turn off the non-perturbative effect by including only the first term Z pert,< as well in the one-point function Z(β) Airy,λ . Next, we will study the quenched free energy first with and then without the non-perturbative contribution.

Without the non-perturbative contribution
We first study the quenched free energy without the contribution from the non-perturbative part. In this case, the numerical calculation is straightforward without subtleties. We are able to numerically confirm that the deformed quenched free energy F q,λ (T ) monotonically decreases as T increases. Furthermore, we find the quenched free energy F q,λ (T ) monotonically decreases as the absolute value of λ increases and present our numerical results below in Figures 8 and 9.

With the non-perturbative contribution
Here, similarly to what has been done in §3.3 (Figures 6 and 7 in particular), we numerically demonstrate the integral diverges by showing the difference D(x) ≡ e −Z λ (x) − e −x Z λ does not vanish in the large x limit. Notice that e −x Z λ → 0 in the large x limit and we can show that D(x) diverges as x → ∞ by proving Z λ (x) is oscillating with its amplitude growing rapidly with x. All the analytical detail will be provided in Appendix B. Here, instead, we will numerically plot Figure 10 to illustrate this fact by showing the depth of the three valleys 23 deepening as x increases which shows the divergence of D(x). Figure 9: We plot the quenched free energy F q,λ (T ) without the contribution from the non-perturbative branch as a function λ for fixed T = 0.1, 0.2, 0.4, 0.6. As shown, F q,λ (T ) is monotonically decreasing as |λ| increases. One can see that the valleys' depths deepening towards the right makes e −Z λ (x) [and therefore D(x) and F q,λ (T )] unbounded in the large x limit. Yet, it is not obvious how to rigorously derive the locations of these valleys in terms of x and this could be an interesting exercise in (numerical) analysis.
Incidentally, we can also resort to a commonly used trick for numerical improper integrals [74]. We already know that the divergence comes from the x ∈ [1, ∞) part of the x-integral (3.40). On the other hand, because (3.40) is of the form ∞ a f (x) x dx, we can convert the integral to 1/a 0 f (1/t) t dt by a change of variables t = 1/x. One can apply this to (3.40) with a = 1 to observe its divergence over a compact interval t = 1/x ∈ [0, 1].

Conclusion
In this paper, we studied the T T -deformed correlators for JT gravity and its dual matrix model. Additionally, we computed the quenched free energy of the Airy model under the T T deformation for both signs of λ and their non-perturbative features.
We briefly summarize our numerical results. At genus-zero and at the leading order of perturbation theory in λ, we confirmed the quenched free energy F q,λ (T ) is a monotonic function in T for a given λ within the validity domain of the leading order approximation. We also find the good sign λ < 0 deformation decreases F q,λ (T ) while the bad sign λ > 0 increases F q,λ (T ).
For all genus and in the low temperature approximation, we computed the quenched free energy F q,λ (T ) using Okuyama's formula (2.18) which diverges regardless whether we include the non-perturbative contribution of the T T deformation or not when λ > 0. For the good sign λ < 0, we are able to numerically compute F q,λ (T ) without including the potential non-perturbative contributions and confirm the monotonicity of F q,λ (T ) at low temperature. Additionally, we find F q,λ (T ) decreases by the deformation and matches the result in perturbation theory at genus-zero. However, when including the possible contributions from the non-perturbative branch, we find F q,λ (T ) computed from Okuyama's formula (2.18) diverges again.
We conclude with a few open questions and future directions relevant to extending this work.
• The first question would be to understand the source of the divergence in the deformed quenched free energy F q,λ (T ). While it is possible that Okuyama's formula (2.18) may fail or require modification for the deformed theory. We suspect that the analytical continuation or exchanging order of the infinite sum and integral in the derivation may not hold for the T T -deformed theory. Additionally, it is also likely that the divergence is caused by the non-perturbative instability suffered from the Airy model. There have been works on how to improve the non-perturbative feature of JT gravity to remove this undesired feature [12,15,[60][61][62][63][64]75]. An interesting extension of our work is to study the deformed JT gravity at general temperatures and see if the same divergence we found in the Airy model appears as well as check if the non-perturbative instability would affect the quenched free energy qualitatively. For the good sign λ < 0 without the non-perturbative contributions from the T T deformation, we expect the finite F q,λ (T ) to decrease hold in JT gravity at general temperatures.
• As we have pointed out, there are non-perturbative effects in the Airy model which endow the density of states ρ Airy (E) with support on the entire real axis making the integral transformation (2.3) no longer hold. However, if we consider the genus expansion then it is straightforward to compute the deformed one-point correlator. For instance, when λ < 0, we have the deformed one-point correlator relevant to the genus expansion Now, a natural question to ask is if one can sum over all Z(β) g Airy,λ to produce the deformed one-point correlator Z(β) Airy,λ ? Unfortunately, as one can easily check, Airy,λ diverges. We naively expect a resurgent analysis to possibly help understand this.
• As we have explained, our version of the T T deformation for the matrix model seems to require the measure in the matrix integral to be deformed as well. It would be important to further explore and figure out how to treat the branch cut in the half change of variables (3.6)-(3.7).

Acknowledgments
We particularly thank Alexander Frenkel and Jiabao Yang for numerous discussions on JT gravity and matrix models. A Full rigorous details on the quenched free energy with bad sign λ > 0 Here we will show that the quenched free energy defined in (3.40) always diverges when λ > 0.

A.1 Perturbative and nonperturbative branches combined
As already shown numerically near Figures 6 and 7, the difference D(x) = e −Z λ (x) −e −x Z λ between the two numerators in the integrand (3.40) does not asymptote to zero as x → ∞. This renders the x-integral to be divergent. In this appendix, we analytically examine the large-x behavior of the integrand in (3.40) and show that the integral diverges when λ > 0.
We will see that this conclusion does not change regarding the fine details of the undeformed density of states ρ Airy (E).
The second term e −x Z λ in (3.40) surely vanishes as x → ∞, and we will show that the exponent −Z λ (x) in the first term vanishes as x → ∞. We define the portion of the E-integrand in (3.49) as where β has been replaced by the inverse temperature 1/T , and we will work with T in both Appendices A and B. First, we notice that, because x, λ, T > 0 then Next, it is straightforward to see that the function h is monotonically decreasing with E: respectively, so as temperature T decreases, h becomes steeper and steeper at E = 0 and the other way around (although very inconspicuously) when λ decreases. Note that according to (A.4), the slope ∂h/∂E has a limit − xe −x T as λ → 0 and because of the exponentially suppressing factor e − 1 λT −xe − 1 2λT in (A.7), when x = 1, the graph of h becoming steeper near E = 0 as λ decreases is not obvious at all as shown in Figure 11. In contrast, the slope is not bounded as T decreases as shown in Figure 12.
Overall, since we are restricting to small T and λ, we can determine a rough picture of the integral in (3.49) as follows: • Starting from x = 1, and the value h(0, 1, λ, T ) at E = 0 is roughly 1 − 1/e (with small T and λ); • Then if we decrease T or increase λ, we see clearly that the slope becomes steeper; 24 One warning is that when λT is sufficiently small, around 10 −3 , numerical calculations (by Mathematica) are no longer reliable. To improve precision, we use MinRecursion→ 9 or WorkingPrecision→ 10 for NIntegrate.
-��� • Most importantly, if we keep increasing x beyond a sufficiently large value (which parametrically depends on T and λ), h will eventually keep shifting to the left; • Now (3.49) is an E-integral over the product of h(E, x, λ, T ) and ρ Airy (E), which is essentially a "convolution" between h and ρ Airy (E) [despite the changing shape of h, and the integration range here being E ∈ (−∞, 1/8λ) instead of E ∈ (−∞, +∞)].
Since ρ Airy (E) has a fixed shape, this "convolution" only depends on the shape and position of h(E, x, λ, T ). We can see that as x → ∞, h(E, x, λ, T ) moves to the far left, and consequently the overlap between h(E, x, β, λ) and ρ Airy (E) approaches zero making the E-integral in (3.49) vanish. Figure 13: The graph of h against E monotonically shifts to the left when x is sufficiently large. We set λ = 0.1 and T = 1.

A.2 Perturbative branch alone
Now let us take a step back and examine if we only include the perturbative branch, will the quenched free energy F q,λ (T ) still diverge? Below we run an analysis in parallel with the previous subsection's, and we will see that the answer is a resounding yes.
This time we only include the first term from the integrand in (3.49): and we define We notice that Next, it is straightforward to see that the E-integrand is monotonically decreasing with E: Now let us examine how the value of h pert (0, x, λ, T ) = 1 − e −x changes with x: so as x increases, the graph of h against E keeps moving to the right (different from the previous subsection), but it eventually reaches a limit at E = 1/8λ as shown in Figure 14. Since this limiting shape, a reflected and shifted Heaviside step function 1 − θ E − 1 8λ , does not depend on either λ or T there are no need for the plots similar to Figures 11 Figure 14: We show the graphs of h pert against E monotonically shifts to the right when x is sufficiently large. We set λ = T = 0.1.

B Proof of the divergence of the good-sign quenched free energy with non-perturbative effects
Since inside the x-integral of (3.40), e −x Z λ → 0 as x → ∞ due to the finiteness of Z λ even with its nonperturbative part in (3.56) included, to discuss the convergence of this x-integral, again we only need to examine e −Z λ (x) . The first term in its exponent Z λ (x) as in (3.58) is benign, so we focus our attention on the second term (3.59), namely a highly oscillatory E-integral because of I(E, x, λ, T ). In Figure 15, the plot of I(E, x, λ, T ) against E shows that it has infinitely many "mountains" of the same height, and within each "mountain", there are many peaks. 25 Figure 15: We plot I E < 0, x = 1 2 , λ = − 1 6 , T = 1 3 and show the first three "mountains", each of which contain several sharp peaks. All "mountains" have the same height e xe − 1 4λT , and they are all symmetric with respect to their vertical axes. Here we choose small x merely to make plotting easy, but yet the peaks are already sharp and localized.
We then notice that I(E, x, λ, T ) is a product of a "waveform" cos xe − 1 4λT sin √ 8λE−1 4λT and an "envelope" e −xe − 1 4λT cos √ 8λE−1 4λT , and the number and height of peaks increase with x. For large x, the "mountains" (i.e. maxima of the "envelope") are well localized around the local minima of cos On the other hand, the asymptotics of ρ Airy (E) in the E-integral (3.59) can be determined from the asymptotics of the Airy function This means Ai(−z) and Ai (−z) behave as e 3 E 3 as E → −∞. Therefore, for a given x, it quickly wins over the amplitude 25 The state-of-the-art methods for infinite highly oscillatory integrals are Longman's method [76] and double exponential quadrature method by Ooura and Mori [77]. Although we use neither of them, our treatment from now on is inspired by the former.
of I(E, x, λ, T ) in (3.59) as E decreases. Hence, to work out the asymptotics of Z λ (x) as x → ∞, we only need to focus on the first or rightmost mountain centered at E * where n = 0 in (B.1). We define the horizontal range of the mountain to be between two adjacent maxima of cos √ 8λE−1 4λT in the exponent of the "envelope", namely from E L = (8λT π) 2 +1 8λ to E R = 1 8λ . Notice that all these treatment of that "mountain" is independent of the "waveform".
We furthermore write E as (4λT π) 2 +1 8λ + with a very small , then in the "envelope" and the "waveform", Due to large x, this interval is much narrower than the range (between E L and E R ) of the first "mountain" just defined above. We will adopt this new truncation hereafter. Furthermore, because the first "mountain" becomes infinitely narrow as x → ∞, the density of states ρ Airy (E) in (3.59) can be approximated by ρ Airy (E * ). As a result, up to a multiplicative factor, the E-integrand there becomes .

(B.5)
Empirically, can we ignore the 2 term in cosine and yet approximate well? From Figure  (16) below, we see that the first-order is a bit off so we include up to the second order in 2 .
After setting xe − 1 4λT = a and 4πλT 2 = b, we only need to integrate 26 the underbraced part in (B.5) over the "mountain" range [ , − ] is:  Figure 16: We show the comparison between the plots of the original waveform (green) and its first-order (blue) and second order (green) approximations in , where λ = −1/10, T = 1/2, x = 1 (again, to make his figure visually discernible, a small x is chosen). We see that the first-order one clearly deviates from the other two plots.
(B.17) Finally, we plug in θ from (B.12) to take the real part of the exponential and arrive at π 2 a 1 + π 2 − 8 + i We conclude that a multiplicative part Z cos λ (x) of Z λ (x) in the x-integral (3.40) is oscillating with an exponentially diverging amplitude as x → ∞. Since we know that the integral (3.40) 27 We remark that other famous estimates, which are resoundingly accurate numerically, such as Identity where f k (x, y) = 2x(1 − cos 2xy cosh ky) + k sin 2xy sinh ky , g k (x, y) = 2x sin 2xy cosh ky + k cos 2xy sinh ky , (B. 16) and | (x, y)| ≤ 10 −16 |erf(x + iy)|, is not useful here. The first line alone is terribly far from being sufficient because for quite a few small k > 0, the absolute values of the real and imaginary parts in each summand of the second line grow as e −x 2 2πx e ky−k 2 /4 when x, y → ∞ makes it more dominant than the second term in the first line of (B.15).
converges if we restrict to Z pert,< λ (x) in Z λ (x), including Z cos λ (x) will make the new integral divergent, rendering no hope of numerically evaluating the quenched free energy with the non-perturbative effect.
We finish by remarking that in principle, one can find a better estimate on ρ Airy (E) over the narrow range [E * + , E * − ] by modelling it as a polynomial in instead of a constant ρ Airy (E * ) and discover more accurate asymptotics than (B.20). To this end, below, we collect some indefinite integrals 28 of the products of monomials in up to degree 3, Gaussian function and cosine functions: (π + i) 7/2 π 3/2 √ 2a(3i − π(a − 3))b 2 F √ a(πb + i(i + π) ) b √ i + π + √ π + i π(π(a − 2) − 2i)b 2 + π(1 − iπ)ab − (π + i) 2 a 2 , (B.23) where F (·) denotes the Dawson function F (x) ≡ e −x 2 x 0 e t 2 dt. They are absolutely not necessary to our proof of divergence of (3.40) with Z cos λ (x) included. 28 Compact results in terms of the Dawson function can be obtained for monomials up to 7 .