JT gravity, KdV equations and macroscopic loop operators

We study the thermal partition function of Jackiw-Teitelboim (JT) gravity in asymptotically Euclidean AdS2 background using the matrix model description recently found by Saad, Shenker and Stanford [arXiv:1903.11115]. We show that the partition function of JT gravity is written as the expectation value of a macroscopic loop operator in the old matrix model of 2d gravity in the background where infinitely many couplings are turned on in a specific way. Based on this expression we develop a very efficient method of computing the partition function in the genus expansion as well as in the low temperature expansion by making use of the Korteweg-de Vries constraints obeyed by the partition function. We have computed both these expansions up to very high orders using this method. It turns out that we can take a low temperature limit with the ratio of the temperature and the genus counting parameter held fixed. We find the first few orders of the expansion of the free energy in a closed form in this scaling limit. We also study numerically the behavior of the eigenvalue density and the Baker-Akhiezer function using the results in the scaling limit.


Introduction
The Sachdev-Ye-Kitaev (SYK) model [1][2][3] and its holographic dual Jackiw-Teitelboim (JT) gravity [4][5][6][7][8][9] are useful testing ground to study various issues in quantum gravity and holography. This duality is based on the fact that the 1d Schwarzian theory, which arises from the Nambu-Goldstone mode of the spontaneously broken time-reparametrization symmetry of the SYK model, also appears as the boundary mode dynamics of JT gravity on JHEP01(2020)156 asymptotic AdS 2 . This duality tells us that the random average of the thermal partition function Z(β) = Tr e −βH SYK of the SYK model reduces, at large number N SYK of fermions and at low energy, to the partition function of JT gravity on Euclidean AdS 2 which is topologically a disk with renormalized boundary length β.
Recently, Saad, Shenker and Stanford [10] found that one can go beyond the strict large N SYK limit and actually compute the partition function of JT gravity including the contribution of various topologies adding handles (or Euclidean wormholes) to the disk. They proposed that the partition function Z JT (β) of JT gravity on asymptotically AdS 2 space is defined by a certain double-scaled random matrix integral Tr e −βH , where the Hamiltonian of the SYK model H SYK is replaced by a random hermitian matrix H. Then the sum over topologies is reproduced from the 1/N expansion of the matrix integral with N ∼ e N SYK .
Their proposal comes from the following facts [10]: the path integral of JT gravity reduces to the contribution of the Schwarzian mode describing the boundary wiggles, together with the Weil-Petersson volume V g,1 (b) of the moduli space of Riemann surfaces with g handles and one geodesic boundary of length b. The crucial point is that the recursion relation obeyed by the Weil-Petersson volume found by Mirzakhani [11] is equivalent to the topological recursion of a double-scaled matrix model with the spectral curve y = 1 2 sin(2z) [12]. Moreover, the genus-zero eigenvalue density ρ 0 (E) corresponding to this spectral curve is exactly equal to the eigenvalue density computed from the Schwarzian theory [13]. Since the topological recursion of Eynard and Orantin [14] is essentially determined by the data of spectral curve (or ρ 0 (E)) only, the above observations imply that the boundary Schwarzian theory "knows" how to perform the sum over topologies on the bulk JT gravity side. It is further argued in [10] that this relation between JT gravity and the matrix model is generalized to arbitrary number of boundaries.
In this paper we will study the proposal in [10] more closely, focusing on the single boundary case. We find that the matrix model of JT gravity in [10] is nothing but a special case of the old matrix model of 2d gravity coupled to c ≤ 1 matter [15][16][17][18] (see also [19] for a review). The important difference from the old story is that in the JT gravity case infinitely many closed string couplings t n are turned on in a specific way: t 0 = t 1 = 0, t k = (−1) k (k−1)! (k ≥ 2). 1 We then introduce a natural two-parameter generalization of Z JT (β) by releasing t 0 and t 1 from the above constraint. Using this we find that the partition function of JT gravity Z JT (β) is written as the expectation value of the macroscopic loop operator Tr(e βQ Π) [22], where Q = ∂ 2 x + u(x) is the Lax operator and Π is the projection to the states below the Fermi level. We will show that this expression of Tr(e βQ Π) naturally includes both of the Schwarzian contribution and the Weil-Petersson volume.
This rewriting of the partition function using the Lax operator Q is not just a formal expression, but is very useful in practice for the actual computation of the genus expansion. We will develop a systematic method of computing the higher genus corrections to Z JT using the Korteweg-de Vries (KdV) equations obeyed by u(x) and ∂ x Z JT , generalizing the JHEP01(2020)156 approach of Zograf [23]. Using this method we have computed the genus expansion of Z JT up to g = 46. It turns out that this genus expansion of Z JT is valid in the high temperature regime (β 1). In the low temperature regime, on the other hand, we can compute Z JT as a series expansion in T = β −1 using the same KdV equations as above. We find that this low temperature expansion can be rearranged by taking a scaling limit, which we will call the 't Hooft limit: and we find the first few terms of F n (λ) in a closed form.
Another interesting quantity to consider is the Baker-Akhiezer (BA) function ψ(E), which is a solution of the Schrödinger equation −Qψ(E) = Eψ(E) and interpreted as the wavefunction of FZZT brane [24,25]. We find that the Laplace transform ψ(λ) of the BA function has a natural expansion in the 't Hooft limit. We also study the behavior of the eigenvalue density ρ(E) and the BA function ψ(E) by numerically evaluating the inverse Laplace transform of e F (λ) and ψ(λ). We confirm the oscillating behavior of ρ(E) and ψ(E) in the classically allowed region E > 0 discussed in [10] which is non-perturbative in the coupling . This paper is organized as follows. In section 2, we develop a technique of the genus and the low temperature expansions of Z JT based on the KdV equations generalizing the approach of [23]. Along the way, we show that Z JT is written as the expectation value of the macroscopic loop operator Tr(e βQ Π). In section 3, we consider the low energy expansion of ρ(E) and ψ(E), as well as the corresponding low temperature expansion of Z JT and ψ in the 't Hooft limit. In section 4, we study the behavior of ρ(E) and ψ(E) numerically. In section 5, we comment on the connected correlator Z(β 1 )Z(β 2 ) conn and its analytic continuation known as the spectral form factor. Finally we conclude in section 6 with some discussions for the interesting future directions. In appendix A we summarize the known facts in the Airy case described by the spectral curve y = z. In appendix B we consider a partial resummation of the genus expansion. In appendix C we consider the so-called string equation for the JT gravity case. In appendix D we summarize useful properties of the resolvent and the wave functions for the Schrödinger equation.

General properties of partition function
In this section we will show that JT gravity is realized as the conventional 2d topological gravity in the background where infinitely many couplings are turned on in a specific way. We will consider the partition function of JT gravity on Riemann surfaces with one boundary and introduce its two-parameter generalization. The generalized partition JHEP01(2020)156 function is closely related to the tau-function for the KdV hierarchy. Using this relation we will derive a simple differential equation which uniquely determines the partition function both in the genus and the low temperature expansions.

JT gravity as 2d gravity in specific coupling background
Before discussing the partition function of JT gravity, let us first recall some useful properties of the partition function of the general 2d topological gravity which we will use shortly.
(See e.g. [26] for a recent review). Let Σ be a closed Riemann surface of genus g with n marked points p 1 , . . . , p n and let M g,n be the moduli space of Σ. We are interested in the intersection numbers ψ 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 . M g,n is the Deligne-Mumford compactification of the moduli space M g,n . Note that (2.1) vanishes unless m + d 1 + · · · + d n = 3g − 3 + n.
For the above correlation functions one can introduce the formal generating function It is proved in [27] (see also [26]) that the intersection numbers involving both κ and ψ's can be obtained from those involving ψ's only. More specifically, let F be the formal generating function that involves ψ's only G is then given by By using this property we will see that JT gravity is nothing but the special case of the topological gravity with t k = γ k .

JHEP01(2020)156
Let us now consider the partition function of JT gravity on two-dimensional surfaces of arbitrary genus with one boundary. In [10] this partition function is evaluated as the one-point correlation function where Z(β) = Tr e −βH is the thermal partition function of a certain Hermitian matrix model. The genus-zero part is to be evaluated separately. Let us first begin with the g ≥ 1 part. In [10] it is evaluated as 2 (2.8) Here, V g,1 (b) is the Weil-Petersson volume of the moduli space of a genus g surface with one geodesic boundary of length b and Z trumpet Sch (β, b) comes from the path integral of the Schwarzian mode on the "trumpet" geometry. V g,1 (b) is given by (2.9) As mentioned below (2.2), the correlation function κ k ψ l 1 g,1 (k, l ∈ Z ≥0 ) has the following property κ k ψ l 1 g,1 = 0 unless k + l = 3g − 2. (2.10) One can thus expand V g,1 as (2.11) By plugging this expression into (2.8) and evaluating the integral, one obtains where we have identified the genus counting parameter as On the other hand, the genus-zero part comes from the path integral of the Schwarzian mode on the disk, which is expressed as [10] (2.14) 2 Throughout this paper we fix the normalization of the Weil-Petersson form by setting α = 1 (see [10]).

JHEP01(2020)156
We see from (2.12) and (2.14) that it is convenient to absorb γ into the normalization of β, β 2π 2 γ → β, (2.15) or equivalently, one can simply set γ = 1/2π 2 . Doing this, we find (2.16) The last expression is obtained from (2.12) with the help of the property (2.10). Putting these expressions together we obtain Plugging this into (2.17) we finally obtain We have thus shown that the partition function of JT gravity on surfaces with one boundary is expressed entirely in terms of the general topological gravity in a specific background t k = γ k .

Generalized partition function and KdV constraints
The relation (2.19) of JT gravity with the general topological gravity provides us with a better understanding of Z(β) as well as an efficient algorithm of computing it. It is well known that the partition function of the topological gravity obeys the KdV constraints [28][29][30]. In fact, Zograf proposed an efficient algorithm of computing the Weil-Petersson volume by making use of the KdV equation [23]. In what follows we will generalize his idea and present a more direct application of the KdV constraints to JT gravity. Let us first recall how the KdV constraints occur in the general topological gravity. It was conjectured by Witten [28] and proved by Kontsevich [29] (see also [30]) that e F with F defined in (2.4) is a tau function for the KdV hierarchy. This means that

JHEP01(2020)156
satisfies the (generalized) KdV equations where R k are the Gelfand-Dikii differential polynomials of u , · · · . (2.22) Here we have introduced the notations (2.23) For k = 1, (2.21) gives the traditional KdV equation (2.24) R k are determined by the recursion relation with the initial condition R 0 = 1. Integrating (2.21) once in t 0 we have In this paper we call the above relations obeyed by F the KdV constraints.
We would like to make use of the KdV constraints to study the JT gravity partition function (2.19). To do this, it is better not to fix the value of t i completely as in (2.19) but rather leave t 0 and t 1 as parameters. In what follows we regard F as a function in t 0 , t 1 (and also in g s ) (2.27) As we will see, at least locally around (t 0 , t 1 ) = (0, 0) one can introduce such a twoparameter deformation. One should also keep in mind that there is no guarantee that F (t 0 , t 1 ) is well-defined for arbitrary values of (t 0 , t 1 ). For our purposes it is convenient to introduce the rescaled parameters := 1 √ 2 g s , x := −1 t 0 , τ := −1 t 1 (2.28) and the notation We then introduce a two-parameter deformation of the partition function (2.19) as We have added the term βt 0 in the definition of Z JT (t 0 , t 1 ) in (2.30) so that we obtain a simple relation has a beautiful interpretation. It is expanded as with coefficients being again the Gelfand-Dikii polynomials (2.35) In this notation R k are written as (2.36) With change of notation u → −u these R k are identified precisely with the original polynomials appeared in the paper of Gelfand and Dikii [31]. This means that their generating function R is the resolvent [31] for the Schrödinger equation

JHEP01(2020)156
Here, |x is the coordinate eigenstate. Note that Q is nothing but the Lax operator L for the KdV equation, which we will discuss later. By taking the inverse Laplace transform of (2.37) we obtain the formal expression From the relation ∂ x Z JT = W in (2.32), we find Introducing the projector Π by Topological gravity and other models of 2d gravity coupled to matter are described by a double-scaling limit of the general matrix model, in which Tr(e βQ Π) is known as (the expectation value of) the macroscopic loop operator [22]. 3 We have thus shown that the partition function of JT gravity on surfaces with one boundary is identified with a single macroscopic loop operator of the matrix model. In this sense JT gravity is merely an example of the old 2d gravity (see [19] for a review). What is special about JT gravity, when compared with the previously known examples, is that infinitely many couplings t n are turned on with a specific value t n = γ n in (2.6).

Lax formalism and master differential equation
As is well known, the KdV equation admits the Lax formalism. This enables us to derive a simple differential equation for W , which can be used to compute Z JT . A crucial fact about the resolvent R is that it is written as [32] (see also appendix D) where ψ ± are certain two independent solutions to the auxiliary linear problem Here Our definition of the sign of x is opposite from that in [22]; in [22] the projector is given by Π = ∞ x dx |x x | while in our definition Π is given by (2.42 (2.49) The first equation is equivalent to the recursion relation (2.25), which is written for R k as From the second equation it immediately follows thaṫ We have thus derived a simple, linear differential equation for W = ∂ x Z JT . Explicitly in terms of Z JT it is expressed as (2.52)

Genus expansion
We can use the differential equation (2.51) together with the KdV equation (2.47) to compute Z JT as a power series expansion in g s . For this purpose, it is convenient to rewrite these equations in such a way that the g s -dependence is manifest As we will see below, the genus expansion of u and W are completely determined by these equations. Prior to the practical computation it is useful to recall the following fact [30]:

JHEP01(2020)156
with u 0 := ∂ 2 0 F 0 . In our present case with t k = γ k (k ≥ 2), where γ k is given in (2.6), it is convenient to introduce the new variables y := u 0 , t := 1 − I 1 (2.59) and the functions Here, J n (z) is the Bessel function of the first kind. One then finds that I n (n ≥ 2) are identified as Therefore, in our case F g (g ≥ 2) is a polynomial in B n (y) (n ≥ 1) and t −1 . Note that B n satisfies and also The old variables (t 0 , t 1 ) and the new ones (y, t) are related as The first equation of (2.64) simply follows from the definition of I 1 in (2.58), while the second equation of (2.64) comes from the classical, → 0 limit of the string equation (C.10) [30] u This relation (2.65) can be interpreted as the stationarity condition ∂F 0 ∂u 0 = 0 of the genuszero free energy [30] In terms of (y, t) the differentials ∂ 0,1 are written as In terms of Λ, I k is written as [33,34] It is interesting to observe that going from t k to I k amounts to shifting Λ 2 → Λ 2 − 2u0. Note that γ k in (2.6) is written as a contour integral on the spectral curve y = 1 2 sin(2

JHEP01(2020)156
This is identical to the change of variables introduced by Zograf (see [23]). One can show that the "on-shell" value (t 0 , t 1 ) = (0, 0) corresponds to (y, t) = (0, 1). It is important to note that the map (2.64) from (y, t) to (t 0 , t 1 ) is neither one-to-one nor onto. This gives rise to a somewhat delicate issue: Z JT and related quantities we are going to solve become multivalued and not globally well-defined (at least as real functions) when viewed as functions in (t 0 , t 1 ). But as far as local expansions around the "on-shell" value are concerned, one can ignore such intricacies. With these preparations let us first consider the genus expansion of u. By construction u has the expansion of the form One can easily show that the differential equation (2.53) is written as the recursion relation Then u g are obtained by recursively solving this equation with the initial condition u 0 = y. First two of them are Precisely speaking, the equation (2.69) by itself does not determine the "integration constant," i.e. there is freedom to add a term linear in t −1 at each step of the recursion. Such a term is, however, forbidden by the relation u g = ∂ 2 0 F g . For instance, the above results of u g can also be obtained from the results of F g (g = 1, 2) in [30] and one can verify that u g (g = 1, 2) do not contain any terms linear in t −1 . More generally, since F g (g ≥ 2) are polynomials in t −1 and the action of ∂ 0 increases the degree of t −1 at least by one, u g (g ≥ 2) cannot have any terms linear in t −1 . Therefore, one can in fact unambiguously determine u g by recursively solving (2.69). Note that the polynomial structure of F g also allows us to compute it unambiguously from u g .
Let us next consider the genus expansion of W . In contrast to u, W depends not only on y, t, g s but also on β, though β does not appear explicitly in the differential equation (2.51). The genus-zero part of W is obtained from (2.40) by ignoring the commutator of ∂ x and u(x) and by performing the integral with respect to the momentum

JHEP01(2020)156
Alternatively, this form can also be obtained from (2.32) by noticing the following property of the Gelfand-Dikii polynomials Let us then expand W as We have chosen the overall factor so that W 0 = 1. By plugging the genus expansions (2.74) and (2.68) into the differential equation (2.51) we obtain the recursion relation where∂ 0 is given by∂ The equation (2.75) by itself again does not determine the t-independent part of W g . However, if we express W in terms of u using (2.32) and (2.22) and consider the genus expansion, we see that the only possible source of t-independent term is u 0 with no derivatives acting on it. As we saw above, the contribution of such u 0 is entirely captured by the overall factor in (2.74) and consequently W g (g ≥ 1) does not contain any t-independent term. We can therefore unambiguously compute W g by recursively solving (2.75), starting from the initial condition W 0 = 1. For instance, we find It is also easy to prove that W g is a polynomial of weight (3g, −2g) in the generators B n (n ≥ 1), β and t −1 , to which weights (n, 1), (1, 0) and (0, −1) are assigned respectively. Finally let us consider the genus expansion of Z JT From the relation ∂ x Z JT = W , one can show that Z g and W g are related by where∂ 0 is defined in (2.76). Note that the extra factor β −1 comes from the difference of the powers of β in the prefactor of W (2.74) and Z JT (2.78). One can prove that when written as a polynomial in t −1 .
It then follows from (2.79) that Z g (g ≥ 1) has the structure Z g = 5g−3 g . Then using (2.79) we can easily compute Z g (g ≥ 1) from the result of W g by iteratively JHEP01(2020)156 in descending order with respect to k. For example, from the result of W 1 in (2.77) we find We have checked that the on-shell value of Z g reproduces the known intersection numbers One can see that Z g becomes small when β 1, i.e. at high temperature T = β −1 1. In this sense the genus expansion of Z JT in (2.78) can be thought of as a high temperature expansion. As we will see in section 2.6 we can consider the opposite low temperature limit T 1. To study the low temperature regime it is useful to define In the last equality we have used the selection rule (2.10). The first three terms are (2.83) Using the above algorithm we have computed Z g up to g = 46. 5 These data provide us with valuable information of the large genus behavior of the genus expansion. In particular, we find the all-genus result of the intersection number κ ψ 3g−2− 1 g,1 with fixed where P (g) is a degree-2 polynomial of g. The = 0 case is computed in [30] with the famous result P 0 (g) = 1. The = 1 case has appeared in [36] with the result

JHEP01(2020)156
From our data of Z g we find P (g) for ≥ 2, which are not known in the literature: (2.86) One can see from (2.84) that κ ψ 3g−2− 1 g,1 does not exhibit the usual (2g)! growth with fixed . The (2g)! growth comes from the opposite end κ 3g−2−d ψ d 1 g,1 with fixed small d [23,37]. One can also see that the sum over genus of (2.84) is convergent, which we will study in detail in section 2.6.
From our data of Z g up to g = 46, we have extracted numerically the large genus asymptotics of κ 3g−2−d ψ d 1 g,1 using the technique of Richardson transformation. We find (2.87) For d = 0 this agrees with the result in [23,37]. Plugging (2.87) into the definition of V g,1 (b) in (2.11), we find the large genus asymptotics of V g,1 (b) (2.89) f 0 (b) agrees with the result in [10]. Note that the above f n (b) vanishes at b = 2πi which is consistent with the property V g,1 (2πi) = 0 [37]. The large genus asymptotics in (2.88) implies that there is a non-perturbative correction of the form This is interpreted in [10] as the effect of ZZ brane [38] sitting at E = − π 2 4 and the instanton action agrees with the value of the effective potential V eff (− π 2 4 ) = π 2 . (See (3.18) for the explicit form of V eff (E)).

Genus-zero part of Z JT
Let us revisit the genus-zero part of Z JT using our expression of the macroscopic loop operator (2.41). At genus-zero, (2.41) is reduced to

JHEP01(2020)156
As discussed in [30], the x-dependence of u 0 is determined from the classical string equation (2.65). Recalling our definition t 0 = x in (2.28), (2.65) is rewritten as For the on-shell value of the coupling t n = γ n (n ≥ 1), this becomes Then we can change the integration variable in (2.91) from x to u 0 via the relation (2.93) (2.94) Here I 0 (2 √ v) denotes the modified Bessel function of the first kind, which should not be confused with I 0 (u 0 , {t k }). Finally, using the relation (2.96) When u 0 = 0 this reduces to the familiar form of the eigenvalue density of the Schwarzian theory [13] ρ 0 (E) = E 0 dv 2π To summarize, our expression of Z JT as the macroscopic loop operator Tr(e βQ Π) in (2.43) automatically includes the contribution of disk topology (g = 0) as well as the higher genus (g ≥ 1) corrections. In other words, we do not have to treat the disk and other contributions separately as in (2.16). Our expression Z JT = Tr(e βQ Π) captures all contributions in one shot.

JHEP01(2020)156
In the large N limit of the matrix model, the genus-zero resolvent obeys an algebraic equation which defines the so-called spectral curve. In our normalization of E and in (2.97), the spectral curve of the matrix model is written as where E, ξ and z are related by As utilized heavily in [10], the genus expansion of matrix model is essentially determined by the spectral curve via the topological recursion [14]. However, to perform the actual computation of the genus expansion with a fixed number of boundaries, the topological recursion turns out to be a very slow algorithm since to compute V g,1 (b) we need to know all the data of V g ,n with g + n ≤ g + 1 (g ≥ 0, n ≥ 1). As emphasized in [23], the method of KdV equation in section 2.4 provides us with a very fast algorithm for the computation of the genus expansion at a fixed number of boundaries.

Low temperature expansion
As we saw in section 2.4, the intersection number κ ψ 3g−2− 1 g with small can be computed for all genus. We observed that the values for fixed are governed by the polynomial P (g). Based on this observation we expect that one can write the expansion of Z(β) = Z JT (0, 0) as (g s β 3/2 ) 2g P (g) 24 g g! . (2.101) When going from the first line to the second line of (2.101) we have removed the restriction 3g − 2 ≥ since P (g) = (24) g g! κ ψ 3g−2− 1 g,1 vanishes when 3g − 2 < . In the last equality of (2.101) we used the property P (0) = 1 to extend the summation to g = 0.
It is natural to expect that Z JT (t 0 , t 1 ) also admits a similar low temperature expansion, namely a power series expansion in T with g s β 3/2 being fixed. In what follows we will see that such an expansion indeed exists and can be computed with the help of the differential equation (2.51) and the results of the genus expansion of u.
To begin with, let us consider the all-genus resummation of the = 0 term in (2.101)

JHEP01(2020)156
We notice that this is essentially the partition function in the Airy case and is related to g s by (2.28). (See appendix A for the summary of the Airy case). In what follows we will use as the genus counting parameter instead of g s . One can generalize the exponential factor in (2.103) to the off-shell value, as follows. When computing Z JT from W = ∂ x Z JT , the above = 0 term at each genus is originated from the O(β 3g ) term in W g . We observe that the O(β 3g ) term, which is of highest order in β, appears in W g as (2.105) This gives rise to the factor Thus it is natural to make an ansatz where we have factored out the genus-zero part (2.72) and the exponential in (2.106) as the prefactor, so that we have w 0 = 1. In the rest of this section we regard T and h as the independent parameters. Plugging this ansatz into (2.54), we find where D is given by and u = u − y is given by Here u g is determined by the recursion relation (2.69). By plugging the expansion (2.108) into the differential equation (2.109) and integrating the O(T ) part with respect to t, one can recursively compute w starting with the initial condition w 0 = 1. To determine w JHEP01(2020)156 uniquely one needs to require that w ( ≥ 1) does not contain any O(t 0 ) term. This can be shown as follows: since W g (g ≥ 1) is a polynomial in t −1 without any O(t 0 ) term, W L − 1 = exp(− h 2 12 t −2 ) ∞ g=0 g 2g s W g − 1 is a formal power series in t −1 without any O(t 0 ) term and so does w ( ≥ 1). Consequently, (2.109) unambiguously determines w . The first two terms of w are given by (2.112) It turns out that w is actually a polynomial of weight ( , 0) in the generators B n (n ≥ 1), t −1 and h, to which weights (n, 1), (0, −1) and (0, 1) are assigned respectively.
Finally, let us expand Z JT as The relation ∂ x Z JT = W is rewritten as Comparing the coefficient of T on both sides of (2.114) we find the recursion relation for z ( ≥ 0) where we formally set z −1 = 0. For instance, using w 0 = 1 and the result of w 1,2 in (2.112) we find z 0 = t, (2.116) For the on-shell value (y, t) = (0, 1), (2.113) becomes (2.119) We have computedz (h) up to = 50. 6 The method above serves as another very efficient algorithm for computing the intersection number κ ψ 3g−2− 1 g,1 , in particular for large g with fixed .
Note thatz (h) is obtained from P (g) as one can easily calculate P (g) fromz (h) by iteratively determining the coefficients of g k (0 ≤ k ≤ 2 ) in descending order. We have thus obtained P (g) also up to = 50. As we explained below (2.101), P (g) vanishes for {g ∈ Z >0 | 3g − 2 < }. We have verified that our results indeed satisfy this property.

Various limits in the low temperature regime
In the low temperature regime we can take various limits of Z JT and the BA function ψ(E). In this section we will consider the low energy limit of ρ(E) and ψ(E), and also the 't Hooft limit of Z JT and the Laplace transform ψ of the BA function. In the rest of this paper we will turn off the deformation parameter t 0 = t 1 = 0 and consider the partition function Z JT and related quantities at the on-shell value of t n = γ n (n ≥ 0).

Low energy expansion of ρ(E)
In this subsection we will consider the low energy expansion of the eigenvalue density ρ(E) in the limit In this limit the genus-zero part ρ 0 (E) in (2.97) is expanded as

JHEP01(2020)156
Let us consider the first term of this expansion It is well-known that this term corresponds to the Airy case as reviewed in appendix A. In this case the genus-zero eigenvalue density in (3.3) is promoted to the full eigenvalue density ρ Airy (E) in (2.104). It turns out that each term of this expansion (3.2) has its own all-genus completion where (η) is defined in such a way that it reduces to the -th term in the expansion of ρ 0 (E) in (3.2) in the classically allowed region E > 0 up to an oscillatory correction. As we discussed above, 0 (η) is given by ρ Airy (E) in (2.104) up to a normalization factor We can determine the higher order terms (η) by matching the low temperature expansion of Z JT in (2.117) From the definition of eigenvalue density

JHEP01(2020)156
For instance, from the result ofz 1 (h) in (2.119), 1 (η) is given by Here the negative power of ∂ η should be understood as the integration with respect to η. One might think that there is an ambiguity in the integration constant, but (η) is actually determined unambiguously by requiring (3.10). Using the data ofz (h) in (2.119) we find (3.13) This procedure enables us to find the all-genus completion of the eigenvalue density order by order in the small E expansion (3.2). In general, (η) is written as a combination of the Airy function Ai(η) and its derivatives. This implies that (η) is exponentially small in the classically forbidden region η > 0, which is indeed necessary for the convergence of the integral (3.10). In appendix B, we consider a partial resummation of this expansion of ρ(E).

't Hooft expansion of Z JT
In the low temperature regime we can take the 't Hooft limit (1.1). As we will see shortly, the relation between Z JT and the spectral curve becomes manifest in this limit. We can rearrange the low temperature expansion in terms of the parameters λ and in (1.1). Plugging the relation into the low temperature expansion of Z JT in (3.8), we find that the free energy is expanded as (1.2). From the data ofz (h) obtained in the previous section, we can compute F n (λ) in (1.2) as a power series expansion in λ. By matching the first few orders of this series expansion, we find the closed form of F 0 (λ) One can show that this is written as

JHEP01(2020)156
Namely, F 0 (λ) is given by the integral of one-form ξdy on the spectral curve. Recall that the effective potential V eff (E) for the eigenvalue is given by the integral of another one-form ydξ As we will discuss in section 3.4, the appearance of the "dual" one-form ξdy in (3.16) can be understood from the Laplace transformation. From the data of series expansion, we also find the closed form of F 1,2 (λ) In section 3.4, we will see that F 2 (λ) can be obtained analytically from the result of topological recursion. Apparently, the above form of F n (λ) becomes singular at λ = 1, and (3.19) can be trusted only in the region λ < 1. If we analytically continue F n (λ) to complex λ, there is a cut running from λ = 1 to λ = +∞ along the real axis of complex λ-plane. It is interesting to understand the physical origin of the singularity at λ = 1.
Before closing this section, we comment on the genus expansion of free energy F = log Z JT . In the original parameters (g s , β) without taking any particular limit, the free energy admits the ordinary genus expansion This expansion is valid in the high temperature regime β 1. On the other hand, in the low temperature regime in the 't Hooft limit (1.1), the free energy is expanded as (1.2). One can recognize that the high temperature expansion (3.20) is "closed string" like, while the low temperature expansion (1.2) is "open string" like. g s and can be thought of as the closed string coupling and the open string coupling, respectively.

Low energy expansion of ψ(E)
In this section we will consider the low energy expansion of BA function ψ(E) in the limit (3.1). This expansion is easily obtained from the expansion of W = x|e βQ |x by using the relation  .

(3.23)
As in the case of (η) in (3.11), ψ(E) 2 can be formally written as From this expansion we can easily find the expansion of ψ(E) in the low energy limit (3.1) The first few terms of the differential operators Ψ read Ψ 0 = 1, More explicitly, ψ(λ) is written as In (3.27) the integration contour C is chosen as the so-called Airy contour running from e − πi 3 ∞ to e

JHEP01(2020)156
From the data of Ψ in (3.26), we can compute G n (λ) as a series expansion in λ. By matching the first few orders of the series expansion, we find the closed form of G n (λ) where F 0 (λ) is given by (3.15). Again, in the next section we will see that G 2 (λ) can be obtained analytically from the topological recursion.

WKB expansion of ψ(E) and ρ(E) from topological recursion
In this section we will systematically compute the semi-classical -expansion (WKB expansion) of ψ(E) and ρ(E) from the topological recursion.

WKB expansion of ψ(E)
Let us first consider the WKB expansion of ψ(E). Once we know the WKB expansion of ψ(E), the expansion of ψ(λ) can be obtained from the saddle point approximation of the integral where E and ξ are related by (2.100). The BA function has the following WKB expansion It is well-known that the leading term S 0 (ξ) is given by the integral of one-form ydξ on the spectral curve (2.99) (see appendix D for a review) where V eff is given by (3.18). In the limit → 0, we can evaluate the integral (3.32) by the saddle point approximation. The saddle point ξ * of (3.32) is given by Then the leading term G 0 (λ) in the -expansion of ψ(λ) in (3.30) becomes As advertised, the integral of dual one-form ξdy naturally arises from the saddle point approximation of (3.32).

JHEP01(2020)156
Let us proceed to the higher order corrections. Using the fact that the BA function is the expectation value of the determinant operator S n (ξ) can be computed from the connected correlators of the operator X = Tr log(E − H) [14] S n (ξ) = As demonstrated in [10], these correlators can be computed systematically by the topological recursion. 7 For instance S 1 (ξ) comes from the cylinder amplitude Then the order O( 0 ) term of the integral (3.32) is obtained by evaluating the Gaussian integral around the saddle point One can check that this reproduces the result in (3.31).
One can easily generalize this calculation to higher order corrections. To do this we set and perform the integral of φ perturbatively by the Wick contraction with respect to the Gaussian measure around the saddle point ξ * where φ 2m is given by Let us compute G 2 (λ) using this formalism. At this order we need S 2 (ξ), which is easily obtained from the topological recursion as , (3.46) 7 Our normalization of the spectral curve y = 1 2 sin(2z) = z + O(z 3 ) is the same as the Airy curve y = z near z = 0. Thus the first few orders of resolvent g 2g−2+n Wg,n(z1, · · · , zn) have the same coefficients as the Airy case W0,1(z) = 2zy(z), W0,2(z1, z2) = 1 (z1 − z2) 2 , W0,3(z1, z2, z3) = where z is the uniformization coordinate defined in (2.100). From our general formula (3.44), G 2 (λ) is given by (3.47) One can check that this reproduces the result in (3.31). We can in principle compute G n (λ) up to any desired order using this formalism.

WKB expansion of ρ(E)
We can repeat the same analysis in the previous subsection for the eigenvalue density ρ(E). It turns out that the 't Hooft expansion of Z JT is related to the WKB expansion of ρ(E) in the forbidden region E < 0. Let us consider the WKB expansion of ρ(E) S 0 (z) and S 1 (z) are given by [10] As discussed in [10], S n≥2 (z) is written as some combination of the connected correlator of the operator Y where Y is given by Here the sign of z in E(±z) distinguishes the two sheets of the spectral curve. In other words, Y is defined by integrating the resolvent from −z to +z. Again, one can compute S n (z) systematically from the topological recursion. For instance S 2 (z) is given by .
One can check that the saddle point approximation of the integral correctly reproduces the free energy F n (λ) in the 't Hooft limit in (3.19). This computation is completely parallel to that in the previous subsection 3.4.1, so we will not repeat it here.

JHEP01(2020)156 4 Numerical analysis of ρ(E) and ψ(E)
In this section we will numerically study the behavior of ρ(E) and ψ(E) as a function of E. Let us consider the integral representation of ρ(E) given by the inverse Laplace transform of Z JT (β) Here we take the contour C to be homotopic to the Airy contour. We will approximate this integral by keeping the free energy F n (λ) up to n = 2 obtained in (3.16) and (3.19) This truncation might be justified when the coupling is small 1. To avoid the cut of F n (λ) running from λ = 1 to λ = +∞ along the real axis, we choose the contour C to cross the real axis in the region 0 < Re(λ) < 1. In practice, in order to evaluate the integral numerically we choose C as a union of three straight segments and use NIntegrate in Mathematica to evaluate the integral. In figure 1 we show the numerical plot of the integral (4.2). As expected, ρ(E) approaches the genus-zero value ρ 0 (E) in the allowed region E > 0. It turns out that the genus-zero part comes from the integral around the origin λ = 0. Although our integration contour (4.3) does not encircle the origin, we can deform the contour to pick up the contribution around λ = 0. However, we should emphasize that the contour (4.3) is completely fixed in the actual numerical computation of the integral (4.2). Near λ = β = 0, we can go back to the original expression (4.1) using β as the integration variable. In the limit → 0 with fixed β, only F 1 (λ) and the first term Note that this is the same as the first two terms in the high temperature expansion (3.20). It is interesting that the genus-zero term 1/β in the original expansion (3.20) becomes a part of F 2 (λ) after taking the 't Hooft limit. Put differently, in order to reproduce the genus-zero part ρ 0 (E) numerically we have to include F 2 (λ) in the approximation (4.2). Then the contribution around β = 0 is evaluated as Using the explicit form of the modified Bessel function Next consider the difference between ρ(E) and ρ 0 (E) We can estimate this difference by the saddle point approximation of (4.2). When E > 0, we can pick up the contribution of two saddle points on the imaginary axis of complex λ-plane λ ± = ±i sinh(2 by deforming the contour C within the homotopy class of Airy contour. Adding the contributions of two saddle points (4.8) we find , (E > 0), (4.9) where the prefactor comes from the Gaussian integral around the saddle points. This agrees with the result of [10] obtained from a different method. In figure 2 we show the plot of ρ np (E). One can see that the numerical value of ρ np (E) fits nicely with the analytic expression in (4.9).
In a similar manner we can numerically compute the BA function ψ(E) in the approximation of keeping G 0 (λ) and G 1 (λ) in (3.31) in the 't Hooft expansion of ψ(λ) (3.30) (4.10) In this case we do not have to include G 2 (λ) for the purpose of numerical analysis since G 2 (λ) = O(λ) in the small λ expansion and hence there is no non-trivial contribution from λ = 0. Again, in the allowed region E > 0 there are two saddle points λ ± in (4.8). Adding the contributions of these saddle points we find , (E > 0). In figure 3 we show the plot of ψ(E). One can see that the numerical value of ψ(E) agrees well with the saddle point result (4.11) in the allowed region E > 0. Let us consider the behavior of ψ(E) in the forbidden region E < 0. Naively, when E < 0 there is a saddle point λ * = sin(2 √ −E) (4.12) and it contributes to ψ(E) as where the effective potential V eff (E) is given by (3.18). It is argued in [10] that this model is non-perturbatively unstable since V eff (E) is not positive definite and ψ(E) blows up as E → −∞. However, we do not see this pathological behavior in the numerical plot of ψ(E) in figure 3. As we can see from figure 4, Re[F 0 (λ)] is negative in the region Re(λ) > 0. Thus, the real part of the leading term Eλ + F 0 (λ) in the WKB expansion (4.10) is negative for E < 0 with an appropriate choice of contour C. This suggests that the JHEP01(2020)156 integral representation of ψ(E) in (3.28) and its approximation (4.10) are convergent and well-defined in the region E < 0 as long as the contour C lies on the right half Re(λ) > 0 of the complex λ-plane, under the condition that C crosses the real axis in the region 0 < λ < 1 to avoid the cut of F 0 (λ). Now suppose that we decrease the value of E from E = 0 toward the negative E direction. At the begging E ∼ 0 the saddle point λ * in (4.12) lies on the positive half plane λ * > 0, but λ * turns negative at E = − π 2 4 and it ceases to contribute to the integral below E = − π 2 4 . This suggests that V eff (E) in (3.18) cannot be trusted for E < − π 2 4 . It is tempting to speculate that this model is actually non-perturbatively stable. It would be very interesting to understand the non-perturbative instability discussed in [10] better.

Comment on the spectral form factor
One can easily generalize our expression of the macroscopic loop operator Z JT = Tr(e βQ Π) to the case of arbitrary numbers of boundaries by applying the general formula in [22] to the JT gravity case t n = γ n . Of particular interest is the connected correlator of two macroscopic loops and its analytic continuation known as the spectral form factor. The spectral form factor is extensively studied in the literature as a useful diagnostics of the quantum chaos of the SYK model and its bulk gravity dual [39][40][41][42].

JHEP01(2020)156
As a function of t, g(β, t) exhibits characteristic features called ramp and plateau. These features naturally correspond to the following decomposition of (5.1) 3) The first term is independent of t and it sets the value of plateau. On the other hand, the second term is a non-trivial function of t and it is expected that this term gives rise to the linearly growing ramp. It is interesting to show this explicitly for JT gravity. (See appendix A for the computation of the spectral form factor in the Airy case). It would also be interesting to consider the bulk gravity picture of plateau. The first term of (5.3) might be interpreted on the bulk gravity side as a geometry where the two boundary circles are merged into a single boundary. Such a geometry was considered before in the context of 2d gravity (see figure 20 in [19]), but its status in the bulk geometry is not clear as mentioned in [19].

Conclusions and outlook
In this paper we have seen that the partition function of JT gravity Z JT (β) = Z(β) is written as the expectation value of the macroscopic loop operator Tr(e βQ Π) in the matrix model of 2d gravity in the closed sting background t n = γ n (2.6). By deforming this background by the two parameters (t 0 , t 1 ), one can utilize the KdV equation to compute the genus expansion of Z JT in a very efficient way. We have also shown that the low temperature expansion of Z JT as well as its 't Hooft limit (1.1) can be obtained systematically. By evaluating the inverse Laplace transformation numerically, we have confirmed the oscillating behavior of ρ(E) and ψ(E) in the region E > 0 as discussed in [10]. Interestingly, the oscillating cosine term arises by adding the contributions of two saddle points (4.8). On the other hand, we do not see any evidence of the pathological behavior of ρ(E) and ψ(E) in the region E < 0 within our approximation. It would be very interesting to understand the non-perturbative instability discussed in [10]. It is desirable to perform more detailed numerical analysis of ψ(E) along the lines of [43]. 8 There are many open questions and interesting future directions. First, it is interesting to understand the physical meaning of the background t n = γ n corresponding to JT gravity. Naively one can imagine that the asymptotic AdS 2 is "built" by this background. To see this more quantitatively, it would be useful to study the Kontsevich's matrix Airy integral [29] corresponding to the background t n = γ n . In the modern interpretation [43][44][45], the Kontsevich's model and Witten's topological gravity [28] are related by the open/closed duality; the Kontsevich's model arises as the open string theory on the FZZT branes while the closed string background t n is obtained by replacing the insertion of FZZT branes with the deformation of matrix model potential. It would be interesting to understand the configuration of background FZZT branes corresponding to t n = γ n (see also footnote 4).

JHEP01(2020)156
It is very important to understand the analytic properties of the genus expansion of Z JT and its non-perturbative completion. Apparently, the 't Hooft expansion of the free energy becomes singular at λ = 1, and the analytic form of F n (λ) in (3.19) can be trusted only in the region λ < 1 when λ is real. It is interesting to understand what happens at λ = 1 (or β = −1 ). Also, it is very important to see if JT gravity is non-perturbatively well-defined. One possible avenue is to study the string equation for u(x) in the background t n = γ n , which we will discuss briefly in appendix C.
In section 3.1 we have constructed the full eigenvalue density ρ(E) as a low energy expansion in the limit (3.1) starting from the Airy case ρ Airy (E). It would be very significant if we can find the exact eigenvalue density ρ(E). It is argued in [46,47] that the eigenvalue density of the SYK model is closely related to the q-Hermite polynomials. It would be interesting to see if the double scaling limit of the q-Hermite polynomials has some connection to the exact eigenvalue density ρ(E) of the JT gravity case.
In section 5 we have briefly commented on the spectral form factor. Using the result of [22] it is straightforward to write down the connected correlator of two macroscopic loops (5.1). It would be interesting to compute it at least in the genus expansion. To this end, we need to know not only the diagonal matrix element W = x|e βQ |x but also the nondiagonal part x|e βQ |y . Fortunately, it is known [48] that the non-diagonal matrix element x|e βQ |y is written as some combination of the derivatives of tau-function, hence it is possible to generalize the method of KdV equation in our paper to the computation of multiboundary correlators. We will report on the computation of multi-boundary correlators elsewhere [49]. The result of the spectral form factor in the Airy case (A. 16) indicates that in the double-scaled matrix model the time scale of the transition to plateau diverges as β → 0, which deserves further investigation.
Finally, it is interesting to extend our approach to more general settings, including JT supergravity [50], adding gauge fields to the bulk theory [51], and a possible analytic continuation to the 2d de Sitter space [52,53], to name a few.

JHEP01(2020)156
and the corresponding classical eigenvalue density is This is realized by a double scaling limit of the Gaussian matrix model by zooming in on the edge of the Wigner semi-circle (see [43] and references therein). In this case u(x) = x and Q is given by In this appendix we will use the normalization t 0 = x, which differs from (2.28) by a factor of . The BA function obeying (Q + E)ψ(E) = 0 is given by the Airy function One can show that ψ(E) in (A.4) is normalized as Now let us consider the one-point function of macroscopic loop operator where the eigenvalue density ρ Airy (E) is given by This can be thought of as the generating function for the intersection numbers ψ 3g−2 1 g,1 . Next consider the connected correlator of two macroscopic loops where Π is the projector The general n-loop amplitude n i=1 Z(β i ) conn has been computed in [54] and the result for the two-loop correlator reads where Erf(z) denotes the error function It is interesting to compare the connected part and the disconnected part of the two-loop correlator as a function of β Here we have set β 1 = β 2 = β for simplicity. In figure 5 we show the plot of the twoloop correlator in (A.13) for = 1/10. One can see that at high temperature (small β) the disconnected part (blue curve) is dominant, while at low temperature (large β) the connected part (orange curve) becomes dominant. As we lower the temperature there occurs the exchange of dominance between the disconnected and the connected part at some critical value β = β crit . A similar phenomenon was observed in the coupled SYK model [55] and it was interpreted as the Hawking-Page transition on the bulk gravity side. Since both Z(β) 2 dis and Z(β) 2 conn depend on β only through the combination 2 β 3 , the critical temperature scales as 9 β crit ∼ − 2 3 . (A.14) We can also study the spectral form factor in the Airy case by analytically continuing the result in (A.11) It turns out that the time derivative of g(β, t) has a simple form (A.16) 9 A similar scaling behavior of βcrit also appeared in the Gaussian matrix model before taking the doublescaling limit [56].

JHEP01(2020)156
This ∂ t g(β, t) decays exponentially at large t and the spectral form factor approaches a constant plateau at late times lim t→∞ g(β, t) = e 2 12 (2β) 3 2 √ π (2β) 3/2 = Z(2β) . (A.17) From (A.16) one can read off the time scale for the transition from ramp to plateau as We notice that t plateau depends on β. This is in contrast to the situation in the Gaussian matrix model before taking the double scaling limit where t plateau is independent of β [57].
It is interesting to observe that the β → 0 limit of g(β, t) is singular due to the one-loop factor (β 1 + β 2 ) −3/2 in (A.11). In [52,53] the analytically continued two-loop correlator Tr e −i H Tr e i H in the JT gravity case was interpreted as the inner product of Wheeler de Witt wave functions of the 2d de Sitter space. This inner product naively corresponds to g(0, ), which is divergent. Interpretation of the singularity of g(β, t) at β = 0 is unclear at present.

B Partial resummation of the eigenvalue density
In this appendix we consider a partial resummation of the low energy expansion of ρ(E). We observe from (2.119) that lim By the inverse Laplace transformation we find a simple closed form expression of the eigenvalue density ρ partial (E) for Z partial (β) This is expected since Z partial (β) in (B.2) reduces in the limit → 0 to the genus-zero part of Z JT in (2.16).
To see that this is indeed the case, in figure 6 we show the plot of ρ partial (E) in (B.3) for = 1/30. One can see that ρ partial (E) agrees with ρ 0 (E) in (2.97) in the allowed region E > 0 up to an oscillatory correction. This implies that the genus-zero part ρ 0 (E) (2.97) is completely accounted for by ρ partial (E) and the difference from the true density ρ(E) has only oscillatory contribution in the region E > 0 ρ(E) − ρ partial (E) = (oscillatory), (E > 0). (B.5)

C String equation for JT gravity
In this appendix we consider the so-called string equation for u(x) (see [32] for a review). It is known that the genus-zero relation (2.92) can be promoted to the all-genus string equation [58] [P, Q] = , (C.1) which arises from the compatibility condition for the following set of equations obeyed by the BA function Qψ = ξψ, P ψ = ∂ ξ ψ.

(C.2)
Here Q is defined in (2.39). To find P , we start with the relation where t k is defined by t k = t k − δ k,1 . Here the subscript + of Q k+ 1 2 + indicates that we truncate the pseudo-differential operator Q k+ 1 2 to its differential part. In this notation M in (2.45) is written as M = M 1 = 2 3 Q 3/2 + . From (C.3) and (C.5) we find that P is given by The compatibility of the flow equation and Qψ = ξψ leads to the following relation where we used u = ∂ 2 0 F and (2.26). Then the string equation (C.1) becomes which can be integrated as (C.10) Using R 0 = 1 this is more compactly written as ∞ k=0 t k R k = 0. (C.11) From the behavior of R k in (2.73), one can see that (C.10) reduces to (2.65) in the classical limit → 0. Note that the shift of t 1 in (C.4) is important to recover the classical equation (2.65). This equation (C.10) determines the x-dependence of u(x). For instance, the string equation for the pure gravity t 0 = −R 2 is known as the Painlevé I equation. The so-called minimal string theory (2d gravity coupled to a minimal model CFT) [20,21] is obtained by turning on a finite number of couplings t k , in which case the string equation can be solved at least numerically [43,59,60].
For the JT gravity case t n = γ n (2.6), P in (C.7) becomes P = 1 2 sin 2Q (C.12) In the classical limit, this reduces to the spectral curve in (2.99) by the replacement P → y, Q → ξ. Eq. (C.12) can be thought of as the "quantum spectral curve" for the JT gravity. It would be interesting to study the property of (C.12) along the lines of [61,62]. Let us consider the string equation (C.10) for the JT gravity case. We set t n = γ n for n ≥ 1 and leave t 0 as a free parameter. Then (C.10) becomes

JHEP01(2020)156
Due to the fact that infinitely many couplings t n are turned on, (C.13) is no longer a differential equation for u; it is rather thought to be a certain non-linear difference equation for u. Based on this expectation, we would like to write down the string equation (C. 13) in the form (C.14) The operator D n can be found from the recursion relation of R k (2.25). The first two terms are ) . indicates that D n should be regarded as a difference operator rather than a differential operator. It would be interesting to find the general structure of D n .

D Resolvent and wave functions
In this appendix we summarize useful properties of the resolvent and the wave functions for the Schrödinger equation.
As discussed in [36,63]  Let us consider the classical limit of S. On general ground, we expect that S cl is written as where y is given by the classical limit of P in (C.7) Evaluating the integral in (D.10) we find On the other hand, we can take the classical limit of (D.9) directly. At the classical level = 0, one can see from (D.4) that R has a square-root branch cut Plugging (D.13) into (D.9) we find S cl = 1 dt 0 ξ − u 0 . (D.14) Using the classical string equation we can rewrite the t 0 -integral to u 0 -integral t k 2 k ξ k+1/2 (2k + 1)!! .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.