Multi-boundary correlators in JT gravity

We continue the systematic study of the thermal partition function of Jackiw-Teitelboim (JT) gravity started in [arXiv:1911.01659]. We generalize our analysis to the case of multi-boundary correlators with the help of the boundary creation operator. We clarify how the Korteweg-de Vries constraints arise in the presence of multiple boundaries, deriving differential equations obeyed by the correlators. The differential equations allow us to compute the genus expansion of the correlators up to any order without ambiguity. We also formulate a systematic method of calculating the WKB expansion of the Baker-Akhiezer function and the 't Hooft expansion of the multi-boundary correlators. This new formalism is much more efficient than our previous method based on the topological recursion. We further investigate the low temperature expansion of the two-boundary correlator. We formulate a method of computing it up to any order and also find a universal form of the two-boundary correlator in terms of the error function. Using this result we are able to write down the analytic form of the spectral form factor in JT gravity and show how the ramp and plateau behavior comes about. We also study the Hartle-Hawking state in the free boson/fermion representation of the tau-function and discuss how it should be related to the multi-boundary correlators.


Introduction
Jackiw-Teitelboim (JT) gravity [1,2] is a very useful toy model to study various issues in quantum gravity and holography. As discussed in [3][4][5][6], JT gravity is holographically dual to the low energy Schwarzian sector of the Sachdev-Ye-Kitaev (SYK) model [7,8]. In a recent paper [9] Saad, Shenker and Stanford showed that the partition function of JT gravity on asymptotically Euclidean AdS spacetime is equal to the partition function of a certain double-scaled random matrix model and the contributions of higher genus spacetimes originated from the splitting and joining of baby universes is captured by the 1/N expansion of the matrix model. See also [10][11][12][13][14][15] for related works in this direction. This opens up an interesting avenue to study the effect of topology change in holography using the powerful techniques of the random matrix theory. This connection between JT gravity and the random matrix model comes from the fact that the density of states in Schwarzian theory is exactly equal to the planar (genus-zero) eigenvalue density of the random matrix model which arises in the topological recursion of the Weil-Petersson volume [16]. This connection is very interesting from the viewpoint of holography. It clearly shows that JT gravity is dual to an ensemble of boundary theories and the partition function on asymptotic AdS spacetime with renormalized boundary length β is interpreted as the ensemble average Z(β) = Tr e −βH over the random Hamiltonian H.
In our previous paper [12], we showed that JT gravity is nothing but a special case of the Witten-Kontsevich topological gravity [17,18] and we studied the partition function Z(β) of JT gravity on spacetime with a single asymptotic boundary in detail. In particular, we found that Z(β) is written as the expectation value of the macroscopic loop operator in 2d gravity [19]. The important difference of JT gravity from the known example of 2d gravity is that infinitely many couplings t k are turned on with a specific value t k = γ k with By generalizing the approach of Zograf [20], we found that the contributions of the higher genus topologies can be systematically computed by making use of the KdV constraint obeyed by the partition function. As emphasized in [20], this method serves as a very fast algorithm for the higher genus computation compared to the Mirzakhani's recursion relation for the Weil-Petersson volume [21]. We also found that in the low temperature regime the genus expansion can be reorganized in the following scaling limit, which we call the 't Hooft limit → 0, β → ∞ with λ = β fixed, (1.2) where ∼ e −S 0 is the genus-counting parameter. In this limit the free energy log Z(β) admits the 't Hooft expansion log Z(β) = ∞ k=0 k−1 F k (λ) (1.3) and we found the analytic form of the first few terms of F k (λ). It turned out that leading term F 0 (λ) is given by the Legendre transform of the effective potential for a probe eigenvalue in the matrix model. In the present paper, we will study the partition function of JT gravity on spacetimes with multiple boundaries by generalizing the method of KdV equation in [12]. We find that the KdV constraints for the connected part of the multi-boundary correlator i Z(β i ) c is obtained by acting the "boundary creation operators" B(β i ) to the original KdV equation for the potential u(x). Our boundary creation operator B(β) is the same as the one discussed in the old 2d gravity literature [22,23] which is based on the idea that the macroscopic loop operator is expanded in terms of the microscopic loop operators in the limit β → 0, up to the so-called non-universal terms which scale with negative powers of β. We can systematically compute the genus expansion of the correlator i Z(β i ) c by solving this KdV constraints recursively. Most of the computation can be done away from the "onshell" value of the couplings (1.1). In particular, we define the off-shell generalization of the effective potential and its Legendre transform, the off-shell free energy. We find that the multi-boundary correlators can be written in terms of a certain combination of the off-shell free energy in the 't Hooft limit. We also study the WKB expansion of the Baker-Akhiezer (BA) functions. 1 In this paper we will focus on the two-point function and compute its genus expansion using the above formalism. We also compute its low temperature expansion and study its behavior in the 't Hooft limit as well. It turns out that the two-point function in JT gravity is expressed in terms of the error function, which is a natural generalization of the known result of pure topological gravity [24]. From the bulk gravity viewpoint, the connected part of the two-point function Z(β 1 )Z(β 2 ) c corresponds to a Euclidean wormhole (also known as the "double trumpet" [25]) connecting the two asymptotically AdS boundaries with renormalized lengths β 1 , β 2 . The analytic continuation of the two-point function Z(β + it)Z(β − it) c , known as the spectral form factor (SFF), is of particular interest in the context of quantum chaos and the SFF is widely studied in the SYK model and JT gravity [25][26][27][28]. We find the analytic form of the SFF in the 't Hooft limit and show that the SFF in JT gravity exhibits the characteristic feature of the so-called ramp and plateau, as expected for a chaotic system with random matrix statistics of eigenvalues.
In a recent interesting paper [29], Marolf and Maxfield considered the boundary creation operators in the context of the AdS/CFT correspondence and made some interesting argument on the baby universe Hilbert space building upon the earlier works by Coleman [30] and by Giddings and Strominger [31,32]. The argument in [29] is mostly based on the intuition coming from a simple toy model, which is not a full-fledged JT gravity. It is interesting to ask how our boundary creation operator B(β) fits into the story in [29], but we do not have a clear understanding of it. We make some preliminary remarks on this problem in section 6 and leave the details for a future work. This paper is organized as follows. In section 2, we compute the genus expansion of multi-boundary correlators using the KdV constraint obeyed by these correlators. In section 3, we consider the WKB expansion of the Baker-Akhiezer functions and the 't Hooft expansion of the multi-boundary correlators. Along the way, we define the off-shell extension of the effective potential and the free energy. In section 4, we compute the low temperature expansion of the two-boundary correlator. In section 5, we study the spectral form factor in JT gravity and show that it exhibits the ramp and plateau behavior as expected for chaotic system. In section 6, we consider the free boson/fermion representation of the τ -function and discuss the boundary creation operator in this formalism. Finally we conclude in section 7. In appendix A, we consider the wavefunction of microscopic loop operators in the 't Hooft limit.

Basics and conventions
In this paper we will generalize our method [12] developed for one-boundary partition function to the case of multi-boundary correlators. To begin with, let us summarize basics, notations and conventions.
As we showed in [12], JT gravity can be regarded as a special case of the general Witten-Kontsevich topological gravity [17,18]. In this model the intersection numbers play the role of correlation functions. They are defined on a closed Riemann surface Σ of genus g with n marked points p 1 , . . . , p n . We let M g,n denote the moduli space of Σ and M g,n the Deligne-Mumford compactification of M g,n . κ (often denoted as κ 1 in the literature) is the first Miller-Morita-Mumford class, which is proportional to the Weil-Petersson symplectic form ω = 2π 2 κ. ψ i is the first Chern class of the complex line bundle whose fiber is the cotangent space to p i and τ d i = ψ d i i . The intersection number in (2.1) obeys the selection rule which we will use frequently. For the above correlation functions one can introduce the generating functions 2 G is actually expressed in terms of F as [33,34] G(s, where γ k is defined in (1.1). Using this property we showed [12] that JT gravity is nothing but the special case of the general Witten-Kontsevich gravity with t k = γ k . Conversely, we can define a natural deformation of JT gravity by (partially) releasing t k from the constraint t k = γ k and regard them as deformation parameters. This is one of our main ideas in [12] and enables us to investigate JT gravity using the techniques of the traditional 2d gravity. In what follows we will apply this prescription to multi-boundary correlators.
In this paper we study the n-boundary connected correlator of JT gravity We consider two kinds of its generalizations, Z n (β 1 , . . . , β n ; t 0 , t 1 ) and Z n (β 1 , . . . , β n ; {t k }). They are obtained respectively by releasing only t 0 , t 1 or all {t k } from the constraint t k = γ k . We often call them "off-shell" correlators. They are related to the "on-shell" correlators (2.5) as As in [12] we introduce the notations and As discussed in [12], g s is the genus-counting parameter in the high temperature regime while is the natural coupling constant in the low temperature 't Hooft limit (1.2). It is also convenient to introduce We will work mostly on the partially constrained case t k = γ k (k ≥ 2), leaving t 0 and t 1 as parameters. In this case we further introduce the new variables y := u 0 , t := 1 − I 1 (2.11) and the functions B n (y) := J n (2 √ y) Here, J n (z) is the Bessel function of the first kind. B n are related to I n as and satisfy (2.14) The old variables (t 0 , t 1 ) and the new ones (y, t) are related as The differentials ∂ 0,1 are then written in the new variables as 3 The on-shell value (t 0 , t 1 ) = (0, 0) corresponds to (y, t) = (0, 1). At this value B n becomes B n (0) = 1 n! , n ≥ 0. (2.17)

Multi-boundary correlators of general topological gravity
In our previous paper [12] we formulated how to compute the genus expansion of the partition function Z 1 = Z JT of JT gravity on a surface with one boundary. In this section we generalize the method to the case of multi-boundary correlators of general topological gravity. Let us start with the fact that the n-boundary correlator of general topological gravity is given by [22] Here F is defined in (2.3) and the operator B(β) is given by As mentioned in section 1, (2.19) is based on the idea that the macroscopic loop operator Z(β) is expanded in terms of the microscopic loop operator τ d in the limit β → 0 2.20) and the insertion of τ d is represented by the derivative ∂ d when acting on the free energy F . B(β) can be thought of as the "boundary creation operator." We put the symbol " " in (2.18), meaning that the equality holds up to an additional non-universal part [22] when 3g − 3 + n < 0. Such a deviation appears, however, only in the genus-zero part of n = 1, 2-boundary correlators, which we will discuss separately in section 2.3. Note that the complex dimension of the moduli space M g,n is 3g − 3 + n which becomes negative for (g, n) = (0, 1) and (0, 2). As in the case of single boundary the genus expansion can be computed by solving a differential equation which follows from the KdV equation. To see this, let us first introduce Recall that u satisfies the KdV equation [17,18] u = uu + 1 6 u . (2.22) Integrating this equation once in x = −1 t 0 we havė Since B(β i ) commutes with˙= ∂ τ and = ∂ x , we immediately obtain a differential equation for W n by simply acting B(β 1 ) · · · B(β n ) on both sides of the above equation. For instance, by acting B(β 1 ) on both sides of (2.23) we obtaiṅ for W 1 (β 1 ). This is nothing but the differential equation for W in [12] with the identification W 1 = W . By further acting B(β 2 ) on both sides of (2.24) we obtaiṅ (2.25) In general the differential equation for W n may be written aṡ It is important to note that the non-universal parts are entirely absent in (2.26) because all the elements other than W 0 = u 2 appearing in (2.26) are equal to or higher than the third derivative of F 0 (see the discussion in the next subsection).
Finally Z n is obtained by merely integrating W n once in x. This can be done order by order in the genus expansion. As a demonstration we will study in detail the case of two-boundary correlator of JT gravity in subsection 2.4.

Genus zero part
In this section let us consider the genus zero part of the multi-boundary correlator Z g=0 n and calculate W g=0 n = ∂ x Z g=0 n . In fact, Z g=0 n has been already calculated in the literature [22,23,36]. In what follows we will reproduce the results in our notation.
Restricting (2.18) to the genus zero part, we have Recall that F 0 is expressed as [35] F 0 = 1 2 where I 0 and u 0 are defined in (2.9)-(2.10). Note that for v = u 0 we have Using these relations we obtain (2.30) As mentioned above, the last expression is only reliable up to the non-universal part. The non-universal part arises because the correlator is not fully constrained by the intersection number of quantum gravity which is defined only for 3g − 3 + n ≥ 0. However, by taking derivative with respect to t k we can insert the microscopic loop operator τ k into the bracket of the intersection number and increase n by one. By repeating this procedure we can map the computation of the partition function precisely to the integral over the well-defined moduli space of punctured Riemann surfaces. For the one-point function we can remove the non-universal part by differentiating twice in t 0 is obtained by integrating the above relation twice in t 0 . We impose the boundary condition that Z g=0 This is naturally understood from our viewpoint [12] that In other words the true genus-zero part of the one-point function (2.33) including the nonuniversal term is obtained from (2.30) by replacing the region of integration from [0, u 0 ] to (−∞, u 0 ]. Note that if we set t k = γ k (k ≥ 2) the above expression gives the result for JT gravity By further setting t 1 = 0 this reproduces our previous result obtained in [12].
In a similar manner, we can compute the genus zero part of the two-boundary correlator (2.35) We can remove the non-universal part by differentiating once in t 0 By imposing the boundary condition Again the true two-point function is obtained from (2.35) by extending the integration region to (−∞, u 0 ]. Given this expression we can easily determine the genus-zero part of the n-point function by induction in n where we have used the genus-zero version of the KdV flow 4 Finally, the genus zero part W g=0 n = ∂ x Z g=0 n is obtained from (2.31), (2.32) and (2.39) as (2.41)

Two-boundary correlator of JT gravity
In this section we focus on the two-boundary correlator of JT gravity and demonstrate in detail how to compute the genus expansion by the method described in section 2.2. Before explaining our method, let us first briefly recall how to compute the correlator by using the method of [9]. The correlator is evaluated by the path-integral of JT gravity on two-dimensional surfaces of arbitrary genus with two boundaries. As shown in [9], the n-boundary correlator is written as a combination of simple building blocks: the partition function of Schwarzian mode on the "trumpet geometry" Z trumpet (β i , b i ) and the Weil-Petersson volume V g,n (b 1 , · · · , b n ) of the moduli space of Riemann surface with geodesic boundaries with lengths b i (i = 1, · · · , n) where γ is the asymptotic value of the dilaton field at the boundary of spacetime. Then the genus sum of two-boundary correlator is written as where g = 0 and g ≥ 1 parts are evaluated respectively as (2.44) Here we have set as in [12] and we have used the selection rule (2.2). From the above expressions one obtains .
This expansion can be computed up to arbitrary genus in principle given the data of Let us now move on to explaining our method described in section 2.2. Using this method the genus expansion (2.46) can be computed very efficiently, as in the case of oneboundary partition function [12]. Regarding the genus zero result (2.41) we first expand u, W 1 , W 2 as 5 (2.47) The genus zero coefficients are given respectively as By plugging the expansions (2.47) into the differential equations (2.22), (2.24) and (2.25) we obtain the recursion relations where we have introduced the notation ∂ 0,β := e −βy ∂ 0 e βy = ∂ 0 + βt −1 . (2.50) The higher genus coefficients u g , W g,1 , W g,2 are computed by solving these recursion relations with the initial data (2.48). We next expand Z 2 as The coefficient Z g,2 is then obtained from the relation As in [12], the integration in t 0 can be done unambiguously assuming that Z g,2 (g ≥ 1) is a polynomial in t −1 without t-independent term. We find Setting (y, t) = (0, 1) with the on-shell values (2.17) of B n one can check that (2.51) with (2.53) reproduces the expansion (2.46).

On multi-boundary correlator of JT gravity
Using the method of [9] the n-boundary connected correlator for n ≥ 3 is obtained by combining the contribution of n trumpets and the Weil-Petersson volume in (2.42) where we have set γ and g s as in (2.45) and have used the selection rule (2.2). This expression is reproduced from (2.18) as follows. For n ≥ 3 we have (2.55) Note that the non-universal part is absent for n ≥ 3 since the dimension of the moduli space 3g − 3 + n is non-negative in this case. Using the relation (2.4) between F and G we have [12].

WKB and 't Hooft expansions
In this section we study the WKB expansion of the Baker-Akhiezer wave function and the 't Hooft expansion of the multi-boundary correlators. Our new formalism is much more efficient than our previous method based on the topological recursion [12].

Baker-Akhiezer function
As we saw in [12] the Baker-Akhiezer functions ψ ± (ξ; {t k }) are certain two independent solutions of the Schrödinger equation and u is defined in (2.21). More specifically, in terms of the resolvent Let us introduce (3.5) From the differential equation [37,38] 2RR we see that v ± are solutions to the equation Using this equation we can compute the WKB expansion of v ± as follows. Let us assume that v admits the expansion By plugging this form as well as the genus expansion of u at the leading and the next to the leading orders. From these we find where we have introduced the notation At the order of n (n ≥ 2) (3.7) is written as the recursion relation Solving this recursion relation with the initial condition v 0 = +z one can easily calculate the WKB expansion of v + . (Recall that the genus expansion of u is calculated by solving (2.49).) In the same way one can calculate v − starting from the initial condition v 0 = −z, but instead v − is obtained from v + by merely replacing z with −z. The same arguments hold for A ± and ψ ± . Therefore, in what follows we omit the subscript "+" and write with the understanding that As a side remark, note that (3.7) is viewed as a Miura transformation. From this viewpoint v can be viewed as a solution to the modified KdV equatioñ It is also possible to compute the WKB expansion of v by directly solving this equation with the initial condition (3.11).
Using the above method we obtain Next, let us consider the WKB expansion of A = log ψ. We expand A as so that we have By solving this we find Here A 1 is immediately obtained from (3.11) up to the integration constant − 1 2 log(4π). This constant is universal in the sense that it does not depend on the background. Thus it can be determined by the asymptotic expansion of the Airy function which is the BA function for the pure topological gravity corresponding to the trivial background t n = 0 (n ≥ 1). 7 A n (n ≥ 2) is also easily obtained assuming that A n is a polynomial in t −1 without t-independent term. Getting A 0 is less trivial, but one can explicitly check that A 0 given in (3.21) satisfies (3.20). One can also check that is regarded as the off-shell generalization of the effective potential V eff (ξ) discussed in [9,12]: In this way one can in principle compute the WKB expansion of ψ = exp ∞ n=0 n−1 A n up to any order. We note in passing that the on-shell BA function and its derivative are expanded as where η = − 2 3 z 2 . By matching the WKB expansion of ψ, ∂ x ψ and the asymptotic expansion of the Airy function, we find the first few terms of Ψ n , Ψ n (3.25) 7 As reviewed in appendix A of [12], the BA function for the pure topological gravity is given by ψ = Above Ψ n agrees with the result in our previous paper [12] obtained by a different method.

Trace formula
In [12] we showed that the one-boundary partition function is expressed as where we have introduced the projector As shown in [39] the general connected correlator is given by . (3.28) For instance, two-and three-boundary correlators are written explicitly as In general, Z n is a sum of the multi-boundary correlator Tr(e β 1 Q Π · · · e βnQ Π) =: exp K (n) (β 1 , . . . , β n ) . (3.31) In [12] we saw that K (1) = F admits the 't Hooft expansion in what follows we explicitly show that K (n) admits the 't Hooft expansion (3.32)

Darboux-Christoffel kernel
Let |ξ be the energy eigenstate corresponding to ψ(ξ) in section 3.1, namely |ξ is normalized so that By inserting n copies of (3.34) with variables ξ i (i = 1, . . . , n), the multi-boundary correlator is expressed as e K (n) = Tr(e β 1 Q Π · · · e βnQ Π) is the Darboux-Christoffel kernel. Since ψ(ξ) satisfies the Schrödinger equation (3.1), we see that (3.37) The Darboux-Christoffel kernel then becomes Plugging this expression into (3.35) and using the genus expansion of A(ξ) calculated in section 3.1, one can compute the 't Hooft expansion of K (n) , as we see below.

Saddle point calculation
In [12] we considered the 't Hooft expansion of F = log Z JT and calculated the first three coefficients F k (k = 0, 1, 2) at the on-shell value (y, t) = (0, 1).
In what follows let us generalize the calculation to the off-shell as well as the multi-boundary cases. Let us first consider the off-shell generalization of the above free energy. Using the technique developed in the previous sections we have (3.40) The above integral is evaluated by the saddle point approximation. The saddle point ξ * is given by the condition This is equivalent to where z * := √ ξ * − y and V eff (ξ) is the off-shell effective potential defined in (3.22). Inverting this relation we obtain (3.43) As in [12] let us introduce a new variable φ as The integral (3.40) is then written as By expanding the integrand in , the integral in φ can be performed as Gaussian integrals. One can in principle calculate F k up to any order. Evaluating the integral up to O( ) for instance, we obtain (3.46) -18 -In particular, F 0 (λ) is given by the Legendre transform of the effective potential V eff (ξ) = −2A 0 (ξ). By using (3.42) and (3.43) we see that they are expanded as and the results (3.46) reproduce those obtained in [12]: In the same manner as above one can calculate the 't Hooft expansion of the twoboundary correlator. We start from It is clear that the saddle point ξ i * (i = 1, 2) is given by This is the same relation as in the one-boundary case (3.41) and thus λ i and ξ i are related as in (3.42)-(3.43). Evaluating the integral by the saddle point approximation we find where F 0 is given in (3.46). At the on-shell value (y, t) = (0, 1) the above results reduce to (3.53) In the same way one can calculate K (n) k for n ≥ 3. We find that K (n) 0,1 (λ 1 , . . . , λ n ) for n ∈ Z >0 take the universal form where the subscript i should be identified mod n.
We also find that K 2 at the on-shell value is given by where F 2 at the on-shell value is given in (3.49).

Low temperature expansion of two-boundary correlator
In this section let us consider the low temperature expansion of the two-boundary correlator. More specifically, we consider the situation where is small and calculate the expansion of Z 2 (β 1 , β 2 ) in T . To begin with, one can observe that the leading order term in each coefficient Z g,2 (2.53) of the genus expansion (2.51) is independent of y and has the form O(β 3g−1 t −2g ). We find that they can be summed over genus √ β 1 β 2 2π e (β 1 +β 2 )y 1

(4.2)
Here Erf(z) is the error function and we have introduced the notation At the on-shell value (y, t) = (0, 1) (4.2) precisely reproduces the result of the Airy case [24] (discussed also in our previous paper [12]). More generally, including the subleading corrections the two-point function at the on-shell value (2.43)-(2.44) is written as Regarding the above result, as in the case of Z 1 [12] it is natural to make an ansatz The subleading parts z (2) ( ≥ 1) can also be estimated from the data of the genus expansion (2.51), (2.53). We find that z (2) has the structure Interestingly, the above z coincides with the coefficient z of the low temperature expansion of Z 1 studied in [12]. In [12] we saw that z can be calculated by solving a set of recursion relations following from the KdV constraint. Similarly, one can derive recursion relations for g , as we will see below. We should emphasize that from (4.5) z and g contain the all-genus information of the intersection numbers at the fixed power of κ, i.e. κ . Let us first express the small T expansion of Z 2 as (4.11) Using the properties it is not difficult to see that the small T expansion of W 2 = ∂ x Z 2 takes the form where w (h) is the expansion coefficient for W 1 (β) introduced in [12] and b is some polynomial in h, r, t −1 , B n (n ≥ 1). The small T expansions of W 1 (β 1 ) and W 1 (β 2 ) are explicitly written as (4.14) As we saw in [12], w (h) is determined by the KdV constraint for the one-point function Similarly, b can be computed from the KdV constraint for the two-point function by equating the terms proportional to B. Note that the last term of (4.16) is proportional to B. If we formally set B = 0, we obtain the homogeneous equation for W 2 which is equivalent to the KdV constraint for the one-point function (4.15). This justifies the ansatz (4.13) for W 2 (and thus our original conjectures (4.6) and (4.8) for Z 2 ) where the expansion coefficient of the first term is given by that of the one-point function w . Plugging (4.13)-(4.14) into (4.16) and using the relations (4.12) one finds that the recursion relation for b is given by (4.17) Starting from b 0 = 0 one can compute b . For instance, the first term is From ∂ x Z 2 = W 2 , one can show that Starting from g 0 = 0 one can calculate g ( ≥ 1) up to arbitrary high order. We have verified that this indeed reproduces our conjectured results (4.9) estimated from the genus expansion.
A few remarks are in order. First, it is worth noting that the two-point function admits a low-temperature expansion of the form The structure of Z 2 (β 1 , β 2 ) in (4.10) is naturally understood from (4.20) by expanding the error function in T . Note that (3.26) and (4.20) imply that the two terms in (3.29) correspond to Tr(e (β 1 +β 2 )Q Π) = Z 1 (β), One can rewrite (4.20) as  From these expressions one immediately obtains D = D 0 c . Second, as discussed in [12], given the result of the low temperature expansion it is straightforward to take the 't Hooft limit (1.2) and one can rearrange the low temperature expansion as the 't Hooft expansion. From the above results one can compute the 't Hooft where D n is obtained as a double series expansion in (λ 1 , λ 2 ) = ( β 1 , β 2 ). Alternatively, from the relation Erfc( and the result of K (2) in (3.52), one can calculate D n as exact functions. Here, the complementary error function can be expanded in with the help of the asymptotic formula For instance, the leading term is given by The higher order corrections D n≥1 can also be easily obtained from the result of K (2) n≥1 in (3.52). We verified at the on-shell value (y, t) = (0, 1) that the series expansions of D n obtained from (4.27) are in perfect agreement with the exact expressions of D n (n = 1, 2) obtained through (4.29).
Third, one might think that g would be interpreted as the expansion coefficients for Tr(e β 1 Q Πe β 2 Q Π). This intuition, however, is not precise. Rather, by using (3.29), (3.26) and (4.10) Tr(e β 1 Q Πe β 2 Q Π) is rewritten as (4.32) This clearly shows that not only g but also z are involved in the low temperature expansion of Tr(e β 1 Q Πe β 2 Q Π). By rearranging the low temperature expansion as the 't Hooft expansion and using the asymptotic expansion formula (4.30), we explicitly verified at the on-shell value (y, t) = (0, 1) that (4.32) is indeed in agreement with e K (2) = e ∞ k=0 k (k = 0, 1, 2) given in (3.53).

Spectral form factor in JT gravity
In this section we will study the spectral form factor (SFF) in JT gravity using our result of two-point function. The SFF is extensively studied in the SYK model as a useful diagnostic of the late-time chaos [25][26][27][28]. The SFF of chaotic system exhibits a characteristic behavior called the ramp and the plateau. From the bulk gravity perspective, the ramp comes from the Euclidean wormhole connecting the two boundaries. The plateau behavior, on the other hand, is a doubly non-perturbative effect with respect to the Newton's constant whose bulk gravity interpretation is still missing. From the random matrix model picture, the origin of the plateau can be traced back to the universal eigenvalue correlation given by the so-called sine-kernel formula. However, this argument is based on the matrix model before taking the double-scaling limit and the analytic form of the SFF in the JT gravity case has not been obtained yet as far as we know. Using our result in the previous section, we can explicitly write down the analytic form of SFF in JT gravity and see how the ramp and the plateau come about. The SFF is defined by analytically continuing the two-boundary correlator Z(β 1 )Z(β 2 ) c to a complex value of the boundary length β 1,2 = β ± it. It is convenient to define the normalized SFF by Using our result in section 4, this is given by the error function (4.20) We are interested in the late-time behavior of the SFF at the timescale of order t ∼ −1 .
To study this regime, it is natural to take the 't Hooft limit 8 As we have seen in section 4, D is expanded as (4.28) in this 't Hooft limit. To see the behavior of the ramp and the plateau, it is sufficient to take the first term of the 't Hooft expansion where F 0 (λ) is given by (3.49). In Fig. 1, we show the plot of SFF in the approximation (5.4) for = 1/30 with several different values of λ. One can see that the SFF exhibits the characteristic feature of the ramp and the plateau. We observe that the timescale of the transition from ramp to plateau depends on λ as in the pure topological gravity case [12].
In [9] it is argued that the ramp is reproduced from the genus-zero part of the connected correlator In Fig. 2 we show the plot of the genus-zero part (orange dashed curve) and the full result (blue solid curve) for the SFF with = λ = 1/30 as an example. One can see that the genus-zero part captures the growing ramp behavior of SFF at early times. This agreement at early times can be shown analytically using the Taylor expansion of the error function and the small λ behavior of F 0 (λ) ∼ λ 3 12 . The appearance of the plateau behavior is almost guaranteed by the functional form of the error function. However, one can pin down the origin of plateau by looking closely at the late-time behavior of the second term Tr(e β 1 Q Πe β 2 Q Π) in the connected correlator Z 2 (β 1 , β 2 ) in (3.29). As we have seen in section 3.4, this term can be evaluated by the saddle point approximation. For β 1,2 = β ± it, the saddle points are given by 6) and the saddle point value is given by This contribution decays exponentially at late-times and the SFF approaches the plateau value given by the first term Tr(e (β 1 +β 2 )Q Π) in (3.29). These saddle points (5.6) can be thought of as the eigenvalue instantons sitting at the complex conjugate pair of points E = −ξ * 1,2 and the transition from ramp to plateau is induced by the pair creation of eigenvalue instantons as advocated in [40].
Another interesting phenomenon is that the connected and the disconnected contributions exchange dominance as we lower the temperature. This transition is observed in a coupled SYK model [41] and it is expected to occur in JT gravity as well. To see this, let us compare the disconnected part Z(β) 2 and the connected part Z(β) 2 c and study their behavior as a function of β. Here we set β 1 = β 2 = β for simplicity. Since JT gravity becomes a good approximation of the SYK model in the low energy limit, it is useful to study the behavior of two-boundary correlator in the 't Hooft limit. At the leading order in the 't Hooft expansion we find In Fig. 3, we show the plot of (5.8) for = 1/30. One can see that at high temperature the disconnected part is dominant, but as we lower the temperature the connected part becomes dominant below some critical temperature. Thus we succeeded to reproduce the transition observed in [41] directly from the JT gravity computation. In the bulk gravity picture, this is an analogue of the Hawking-Page transition between two different topologies of spacetime. At high temperature the two disconnected Euclidean black holes are dominant while at low temperature the Euclidean wormhole connecting the two boundaries becomes dominant.

Boundary creation operator and Hartle-Hawking state
As we have seen in section 2, we can write the connected n-point amplitude as where F denotes the free energy (2.3) and the operator B(β) is given by (2.19). B(β) can be thought of as the "boundary creation operator." The same operator has been considered in the context of 2d gravity in [22]. (6.1) should be understood as the equality up to the non-universal terms at genus-zero for the one-and two-point functions, which should be treated separately. In a recent paper by Marolf and Maxfield [29], the idea of boundary creation operators is also discussed. The important property of the boundary creation operators is that they all commute and hence can be diagonalized simultaneously. It is argued in [29] that the simultaneous eigenstate of the boundary creation operators, the so-called α-state, can be thought of as a member of an ensemble and the correlator i Z(β i ) is interpreted as the ensemble average. Moreover, by reinterpreting the earlier discussion of baby universes [30][31][32] from the viewpoint of AdS/CFT duality, it is argued that one can define the baby universe Hilbert space from the data of correlators i Z(β i ) and this Hilbert space includes many null states due to the bulk diffeomorphism invariance. 9 To demonstrate these properties, a simple toy model is studied in [29] where the action S of the model has only the topological term given by the Euler characteristic χ of the 2d spacetime.
In this section we will consider whether the proposal in [29] can be generalized to the JT gravity case. Firstly, the boundary creation operator B(β) defined in (2.19)  and hence one can try to diagonalize B(β)'s simultaneously. One immediate problem is that B(β) in (2.19) does not look like a hermitian operator, thus its eigenvalue is not necessarily a real number. According to the proposal in [29], this problem might be resolved on the physical Hilbert space, which is obtained by taking the quotient of the original Hilbert space by the space of null states. We do not have a clear understanding of how this happens. In the rest of this section, we will examine how the proposal of [29] is generalized or modified in the case of JT gravity.
To study the proposal of [29] in JT gravity, it is convenient to use the free boson or free fermion representation of the Witten-Kontsevich τ -function (see e.g. [38,[43][44][45] and references therein) τ = e F = t|V , (6.3) where the state t| is given by the coherent state of free boson t| = 0| exp ∞ k=0 t k α 2k+1 g s (2k + 1)!! (6.4) with α n obeying the usual commutation relation of the free boson [α n , α m ] = nδ n+m,0 , 0|α n = 0 (n < 0). (6.5) Note that only the odd-modes α 2n+1 appear in the state t| in (6.4) since the KdV hierarchy associated with the Witten-Kontsevich τ -function is a mod 2 reduction of the KP hierarchy. Note that the derivative ∂ k with respect to the coupling t k is mapped to the operator α 2k+1 when acting on the state t| ∂ k t| = t| α 2k+1 g s (2k + 1)!! . (6.7) In other words, the microscopic loop operator τ k corresponds to α 2k+1 up to a normalization constant.
To write down the state |V in (6.3), it is useful to introduce the free fermions ψ r , ψ * r (r ∈ Z + 1 2 ) obeying the anti-commutation relation {ψ r , ψ * s } = δ r+s,0 , ψ r |0 = ψ * r |0 = 0 (r > 0). (6.8) They are related to α n by the bosonization Then |V is written as In general any choice of A m,n defines a τ -function of KdV hierarchy, but the Witten-Kontsevich τ -function for the topological gravity corresponds to a specific choice of A m,n . The generating function of A m,n for the Witten-Kontsevich τ -function is obtained in [46][47][48]: where a(z) and b(z) are given by The series a(z) and b(z) appear in the asymptotic expansion of the Airy function and its first derivative, respectively. The important property of the state |V is that it satisfies the Virasoro constraint [49,50] L n |V = 0 (n ≥ −1), (6.13) where the Virasoro generator L n is given by : α 2k+1 α 2n−2k−1 : + 1 16 δ n,0 . (6.14) One can show that [L n , L m ] = (n − m)L n+m for n, m ≥ −1. The constant term 1 16 in L 0 can be thought of as the conformal weight of the Z 2 twist field. Note that L n is written as In other words, the linear term − 1 2gs α 2n+3 in (6.14) arises from the shift of t 1 = t 1 − 1 [35,51].
One can show that |V in (6.10) satisfies the Virasoro constraint order by order in the g s -expansion. To see this, let us expand the state |V as From the result of A m,n in (6.11), the first two terms are given by |0 . (6.17) As an example, let us check the constraint L n |V = 0 at the order O(g 0 s ). At this order the constraint L n |V = 0 reads One can easily show that this is indeed satisfied for |V 0,1 in (6.17). The higher order constraint can be shown in a similar manner. In [52,53], the Virasoro constraint of matrix model is interpreted as the gauge symmetry of closed string field theory in a minimal model background. This suggests that the Virasoro constraint is the analogue of the bulk diffeomorphism invariance discussed in [29]. Now let us consider the Hartle-Hawking state |HH [54]. As discussed in [55], it is natural to identify the Hartle-Hawking state |HH as "the most symmetric state." In the present case, |V is such a state since |V is invariant under the Virasoro generators (6.14). |V can be thought of as the SL(2, C) invariant vacuum corresponding to the identity operator and it is a natural candidate for the no-boundary state. Thus we propose to identify the Hartle-Hawking state |HH with the state |V in (6.10) In particular, this state satisfies the equation L 0 |HH = 0 which corresponds to the Wheeler-DeWitt equation.
Next we consider the interpretation of the correlator i Z(β i ) in JT gravity. The correlator here refers to the full correlator including both the connected and the disconnected parts. One can generalize (6.1) to the full correlator by acting B(β)'s on the τ -function instead of the free energy Z(β 1 ) · · · Z(β n ) B(β 1 ) · · · B(β n ) t|V t|V =: t| B(β 1 ) · · · B(β n )|V t|V . (6.20) By using the dictionary (6.7) the operator B(β) is written as It turns out that the non-universal terms at genus-zero are correctly incorporated by extending the summation to all n ∈ Z. Namely we define the operator Z(β) by Then the full correlator is given by where we used our identification |HH = |V . To see that this is the correct prescription, let us consider the genus-zero part of the one-point function Here we have used t k = (−1) k (k−1)! (k ≥ 1) and Similarly, the genus-zero part of the two-point function becomes .
(6.26) (6.24) and (6.26) agree with the known result of the genus-zero part in JT gravity. Using the relation (6.25) one can show that Z(β)'s commute at least formally Our proposal (6.23) is consistent with the identification of the one-point function Z(β) as the wavefunction of the Hartle-Hawking state, which is usually assumed in 2d gravity literature (see e.g. [23] and references therein) where β| is given by More generally, the multi-point correlator is written as Z(β 1 ) · · · Z(β n ) = β 1 , · · · , β n |HH , β 1 , · · · , β n | = t| Z(β 1 ) · · · Z(β n ) t|HH . (6.30) Our expression (6.23) is different from the proposal in [29] Z(β 1 ) · · · Z(β n ) = HH| Z(β 1 ) · · · Z(β n )|HH HH|HH . (6.31) This difference comes from the fact that the bra and the ket are treated asymmetrically in the free boson/fermion representation of the τ -function (6.3). In other words, our expression (6.23) corresponds to a special (Euclidean) time-slicing of the spacetime where the initial state has no boundary and all the boundaries are on the final state. At present, it is not clear to us how to reconcile our (6.23) and the proposal (6.31) in [29].

Conclusions and outlook
We have studied the multi-boundary correlators in JT gravity using the KdV constraints obeyed by these correlators. Along the way, we have defined the off-shell generalization of the effective potential and have studied the WKB expansion of the Baker-Akhiezer functions as well. In particular, we have computed the genus expansion of the connected two-boundary correlator Z(β 1 )Z(β 2 ) c as well as its low temperature expansion. We have found that the two-point function is written in terms of the error function and the ramp and plateau behavior of the SFF in JT gravity is explained by the functional form of this error function. We have also confirmed the picture put forward in [40] that the transition from ramp to plateau is induced by the pair creation of eigenvalue instantons. There are many interesting open questions. In section 6 we briefly discussed a possible connection to the recent work by Marolf and Maxfield [29] which clearly deserves further investigation. It would be interesting to construct the α-state which simultaneously diagonalizes the operator Z(β) in (6.22) and see how the argument in [29] is generalized to the JT gravity case. In particular, it is interesting to see what the non-factorized contribution Z(β 1 )Z(β 2 ) c coming from the Euclidean wormhole [56,57] looks like in the α-state. The pure topological gravity would be a good starting point to study this problem since the explicit form of the n-point correlator is known in the literature [24,58,59].
It is emphasized in [29] that non-perturbative effects are important to realize the massive truncation of the Hilbert space by the diffeomorphism invariance. The free fermion representation of the state |V in (6.10) is defined by the asymptotic expansion in g s and hence it only makes sense as a perturbative expansion (see (6.12)). However, it is possible to include the effect of D-instanton corrections systematically within this framework [60][61][62]. It would be interesting to study such D-instanton effects in JT gravity and see how they affect the argument of diffeomorphism invariance in JT gravity.
In [63,64] it is argued that the Page curve for the black hole evaporation is correctly reproduced if we include the contribution of replica wormholes in the computation of entropy of Hawking radiation using the replica method in the gravity path integral. One can immediately apply our formalism to compute the contribution of the replica wormholes in pure JT gravity sector. To model the black hole microstates one can add the end of the world (EOW) branes to JT gravity [29,63]. It would be interesting to construct a generalization of the JT gravity matrix model which incorporates the degrees of freedom of the EOW branes.
As discussed in [65,66], the matrix model description of JT gravity can be generalized to the 2d de Sitter space by analytically continuing the boundary length β to imaginary value β → ±i . In [67] the boundary creation/annihilation operators are considered in this de Sitter setting. It would be interesting to see how they are related to our discussion in section 6.
Finally, it would be interesting to generalize our computation in this paper to JT supergravity [10]. In particular, the genus expansion of JT supergravity on orientable surfaces without time-reversal symmetry can be computed from the Brezin-Gross-Witten τ -function [68]. We will report on the computation of JT supergravity case elsewhere.
The derivative of F 0 (λ) with respect to the coupling t n can be computed by using the fact that V eff (ξ) and F 0 (λ) are related by the Legendre transformation. Thus we find ∂F 0 (λ) ∂t n λ fixed = ∂ ∂t n λ fixed λξ * − V eff (ξ * ) = λ ∂ξ * ∂t n − V eff (ξ * ) ∂ξ * ∂t n − ∂V eff (ξ * ) ∂t n ξ * fixed One can go beyond the leading order and compute the higher order correction to the wavefunction of microscopic loop operators by using the off-shell generalization of the free energy F in (3.47). After some algebra, we find the first order correction to the 't Hooft expansion (A.7) where n i > 0 (i = 1, · · · , m).