Liouville theory and Matrix models: A Wheeler DeWitt perspective

We analyse the connections between the Wheeler DeWitt approach for two dimensional quantum gravity and holography, focusing mainly in the case of Liouville theory coupled to $c=1$ matter. Our motivation is to understand whether some form of averaging is essential for the boundary theory, if we wish to describe the bulk quantum gravity path integral of this two dimensional example. The analysis hence, is in a spirit similar to the recent studies of Jackiw-Teitelboim (JT)-gravity. Macroscopic loop operators define the asymptotic region on which the holographic boundary dual resides. Matrix quantum mechanics (MQM) and the associated double scaled fermionic field theory on the contrary, is providing an explicit"unitary in superspace"description of the complete dynamics of such two dimensional universes with matter, including the effects of topology change. If we try to associate a Hilbert space to a single boundary dual, it seems that it will contain less information compared to the non-perturbative bulk quantum gravity path integral and MQM.


Introduction
The studies of SYK and its low energy (hydrodynamic) limit described by the one dimensional Schwarzian theory [1,2,3] revealed a holographic connection with a bulk two dimensional Jackiw-Teitelboim (JT) gravity theory. In fact recent work [8] elucidated that this connection continues to hold for bulk topologies other than the disk, and that the complete bulk genus expansion can be resummed using a particular limit of the (double scaled) Hermitean matrix model The argument supporting this connection, is that the usual double scaling limit of a single Hermitean matrix model can describe the (2, p) minimal models coupled to gravity, and the physics of JT gravity can be reached as a p → ∞ limit of these models 1 . Actually it is quite reasonable to expect such a limiting connection between the JT gravity and Liouville theory. As an example the Liouville equation appears naturally after employing two steps, first the identification of the boundary mode Schwarzian action with the Kirillov coadjoint orbit action on M = Dif f /SL(2, R) that is then identified with a 2d bulk non-local Polyakov action together with appropriate boundary terms [4]. The classical solutions of this latter action are then in correspondence with those arising from Liouville theory, albeit in this case the conformal mode of the metric is a non-dynamical non-normalisable mode fixed by imposing certain appropriate Virasoro constraints. This is an indirect way to say that the only dynamical degrees of freedom left in this problem are those of the fluctuating boundary arising from large diffeomorphisms 2 . This then indicates that the various existing models of Liouville quantum gravity coupled to matter [31,32,22], are in fact richer examples of two dimensional bulk quantum gravity theories. Since a lot is known for these models both at a perturbative and non-perturbative level, it is both conceptually interesting and feasible to elucidate the properties of their holographic boundary duals. That said, we should clarify that from this point of view the various matrix models do not play the role of their boundary duals, but should be instead thought of as providing directly a link 3 to a "third quantised description" of the bulk universes splitting and joining in a third quantised Hilbert space [13]. This interpretation is even more transparent in the c = 1 case for which there is a natural notion of "time" in superspace in which universes can evolve. Simply put, the target space of the c = 1 minimal string plays the role of superspace in which these two dimensional geometries are embedded. The matrix models provide a quite powerful description, since it is possible to use them in order to obtain the partition function or other observables of the boundary duals -from the matrix model point of view one needs to introduce appropriate loop operators that create macroscopic boundaries on the bulk geometry. Let us briefly discuss the case of the partition function. In this case for the precise identification, one should actually use a loop/marked boundary of fixed size that is related to the temperature β of the holographic dual theory [8]. The Laplace transform of this quantity then gives the expression for the density of states (dos) of the boundary dual. In the concrete example corresponding to the (2, p) models, this was shown to reduce to the Schwarzian density of states in the limit p → ∞ [8]. Since the Schwarzian theory captures only the IR hydrodynamic excitations of the complete SYK model, it is then natural to ponder whether and how one could connect various integrable deformations of the (2, p) matrix models with corresponding corrections to the "hydrodynamic" Schwarzian action.
In particular a matrix model with a general potential of the form V (H) = k t k H k is still an integrable system, and it is known that its partition function corresponds to a τ -function of the KP-Hierarchy [31,32]. Similar things can be said about two-matrix models (2MM), with which one can describe the more general (q, p) minimal models [41,43,42,22]. The partition function of such two matrix models takes the general form τ N = Z(N ) = dM dM e −N (tr(MM)+ k>0 (t k tr M k +t k trM k )) , (1.2) and is a τ -function of the Toda integrable Hierarchy with t k 's ,t k 's playing the role of Toda "times". Very interesting past work on the integrable dynamics of interfaces (Hele-Shaw flow) has revealed a deep connection between the dynamics of curves on the plane and this matrix model [51]. In fact the Schwarzian universally appears in the dispersionless limit of the Toda hierarchy when = 1/N → 0, and can be related with a τ -function for analytic curves which in turn is related to (1.2). We will briefly review some of these facts in appendix A.1, since they are related tangentially to this work. The main focus of the present paper will be the case of c = 1 Liouville theory having a dual description in terms of Matrix quantum mechanics of N -ZZ D0 branes [36]. We emphasize again that even though in this case there is a natural interpretation of the theory as a string theory embedded in a two dimensional target space, the Liouville theory being a worldsheet CFT, in the present paper we will take the 2-d Quantum Gravity point of view [34,32] where the worldsheet of the string will be interpreted as the bulk spacetime. This is in accordance with our previous discussion and interpretation of Jackiw-Teitelboim gravity and the minimal (q, p) models. In short we will henceforth interpret the combination of c = 1 matter with Liouville theory as a bone fide quantum gravity theory for the bulk spacetime. For the c = 1 case at hand, at the semiclassical level the dynamical degrees of freedom are then the conformal mode of the two dimensional bulk metric -Liouville field φ(z,z) -together with that of a c = 1 matter boson which we denote by X(z,z). The Euclidean bulk space coordinates will then be denoted by z,z.
The possibility for connecting this bulk quantum gravity theory with holography is corroborated by the fact that once we introduce macroscopic boundaries for the Liouville CFT (corresponding to insertions of macroscopic loops of size on a worldsheet in the usual point of view), near such boundaries the bulk metric can take asymptotically the form of a nearly-AdS 2 space. The adjective nearly here corresponds to the fact that we can allow for fluctuations of the loop's shape keeping its overall size fixed. This is quite important, since it allows for a holographic relation between the bulk quantum gravity theory with a quantum mechanical system on the loop boundary akin to the usual AdS/CFT correspondence 4 . It could also pave the way to understand the appropriate extension and interpretation of the correspondence in the case of geometries having multiple asymptotic boundaries (Euclidean wormholes). Even though several proposals already exist in the literature [18,19,8,20], it is fair to say that no ultimate consensus on the appropriate holographic interpretation of such geometries has been reached (and if it is unique). We provide more details on what we have learned about this intricate problem in the conclusions.
Another natural question from the present point of view, is the role of the original matrix quantum mechanics (MQM) of N-D0 branes (ZZ branes) described by N ×N Hermitean matrices M ij (x). As we describe in the main part of the paper, the dual variable to the loop length that measures the size of macroscopic boundaries of the bulk of space, is a collective variable of the matrix eigenvalues λ i (x) of M ij (x), while the coordinate x is directly related to the matter boson in the bulk. Hence if we interpret the compact as the inverse temperature β of the boundary theory, we are then forced to think of the collective matrix eigenvalue density ρ(λ) as describing the energy spectrum of the dual theory. This is in accord with the duality between JT gravity and the Hermitean one matrix model of (1.1). From this point of view the D0 branes now live in "superspace" and the second quantised fermionic field theory is actually a "third quantised" description of the dynamics of bulk universes. This means that MQM and the associated fermionic field theory provide us with a specific non-perturbative completion of the c = 1 bulk quantum gravity path integral 5 , as the one and two matrix models do in the simpler cases of JT-gravity and (p, q) minimal models.
This discussion raises new interesting possibilities as well as questions. To start with, one can now try to understand at a full quantum mechanical level various asymptotically AdS 2 bulk geometries such as black holes (together with the presence of matter excitations) directly on what was previously interpreted as the worldsheet of strings 6 . In fact one can go even further, using the free non-relativistic fermionic field theory. Based on the analysis of [13] it is a very interesting and unexpected fact that this field theory is non-interacting but can still describe the processes of bulk topology change. This is made possible due to the fact that the field theory coordinate λ is related with the conformal mode of the bulk geometry φ via a complicated non-local transform [34] 7 . Therefore it seems that we have managed to find a (quite simple) non-disordered quantum mechanical system that performs the full quantum gravity path integral and automatically sums over bulk topologies. Remarkably this system is dynamical and defined on the superspace instead of being localised on a single boundary of the bulk space. More precisely, it is an integrable system on superspace where "time" is related to the c = 1 matter field. This point of view then inevitably leads to the following questions: Does the bulk theory actually contain states that can be identified with black holes? Can we then describe complicated processes such as those of forming black holes -what about unitarity if we can create baby universes? Is there any notion of chaos for the bulk theory even though the superspace field theory is an integrable system? Is there a quantum mechanical action defined directly on the boundary of the bulk space (or at multiple boundaries) that can encapsulate the same physics? Would this have to be an intrinsically disorder averaged system such as SYK or could it be a usual unitary quantum mechanical system with an associated Hilbert space? Our motivation hence is to try to understand and answer as many of these questions as possible.
Structure of the paper and results -Let us now summarise the skeleton of 5 A more recent understanding of the non-perturbative completion/s of the model and the role of ZZ-instantons is given in [45]. 6 Similar considerations have previously been put forward by [48,49]. 7 It is an interesting problem whether some similar transform could encode compactly the process of topology change in higher dimensional examples, at least at a minisuperspace level. our paper as well as our findings. In section 2 we review some facts about Liouville theory with boundaries, such as the various solutions to equations of motion and the minisuperspace wavefunctions that are related to such geometries. In appendix A we briefly describe the matrix models dual to the (p, q) minimal models and their integrable deformations and describe the connections with conformal maps of curves on the plane and the Schwarzian. The main focus of our analysis is the c = 1 case. We first review how the fermionic field theory can be used as a tool to extract various correlators in section 3, and then move into computing the main interesting observables: First the boundary dual thermal partition function Z dual (β) both at genus zero and at the non-perturbative level in section 4.1 and then the dual density of states ρ dual (E) in 4.2.
At genus zero the partition function corresponds to the Liouville minisuperspace WdW wavefunction, while the non-perturbative result cannot be given such a straightforward interpretation, but is instead an integral of Whittaker functions that solve a corrected WdW equation encoding topology changing terms. We then observe that there is an exponential increase ∼ e 2 √ µE in the dos at low energies (µ plays the role of an infrared cutoff to the energy spectrum such that the dual theory is gapped), that transitions to the Wigner semicircle law ∼ E 2 − 4µ 2 at higher energies, combined with a persistent fast oscillatory behaviour of small amplitude. These oscillations might be an indication that whatever the boundary dual theory is, it has a chance of being a non-disorder averaged system. Similar non-perturbative effects where also found in [8] and argued to be related to a doubly-exponential nonperturbative contribution to various observables. We will comment on a possible interpretation of such non-perturbative effects from the bulk point of view in the conclusions section 7.
We then continue with an analysis of the density of states two-point function and its fourier transform, the spectral form factor SF F (t) = Z(β + it)Z(β − it) in section 5. The geometries contributing to these quantities are both connected and disconnected. The correlator of energy eigenvalues exhibits the universal short distance repulsion of matrix model ensembles such as the GUE, and resembles very much the sine-kernel. Nevertheless, its exact behaviour deviates from the sine-kernel slightly, indicative of non-universal physics. For the SFF, the disconnected part quickly decays to zero and at late times it is the connected geometries that play the most important role. These are complex continuations of Euclidean wormholes corresponding to loop-loop correlators Ŵ ( 1 , q)Ŵ( 2 , −q) = M 2 ( 1 , 2 ) from the point of view of the matrix model 8 . The relevant physics is analysed in section 5.4 where it is found that at genus zero (and q = 0), they lead to a constant piece in the SFF. The non-perturbative answer can only be expressed as a double integral with a highly oscillatory integrand. A numerical analysis of this integral exhibits the expected increasing behaviour that saturates in a plateau, but on top of this there exist persistent oscillations for which ∆SF F c /SF F c → O(1) at late times. On the other hand for non-zero q the behaviour is qualitatively different. The nonperturbative connected correlator exhibits an initially decreasing slope behaviour, that transitions into a smooth increasing one relaxing to a plateau at very late times. The main qualitative difference with the q = 0 case is that the behaviour of the correlator is much smoother.
In section 6, we analyse the cosmological interpretation of the wavefunctions, as WdW wavefunctions of two-dimensional universes. In order to do so we follow the analytic continuation procedure in the field space proposed in [15,21]. In our description this corresponds to using loops of imaginary parameter z = i . The wavefunctions at genus zero are Hankel functions ∼ 1 z H (1) iq (z). Nevertheless, it is known from the work of [26], that all the various types of Bessel functions can appear, by imposing different boundary conditions and physical restrictions on the mini-superspace WdW wavefunction. We provide a review of the two most commonly employed ones (no-boundary and tunneling), in appendix B. The non-perturbative description seems to encode all these various possibilities in the form of different large parameter limits of the non-perturbative Whittaker wavefunctions (4.3). This corresponds to complexifying the bulk geometries and choosing different contours in field space. Obstructions and Stokes phenomena are then naturally expected to arise, for the genus zero asymptotic answers.
We conclude with various comments and suggestions for future research.
A summary of relations -We summarize in a table the various relations that we understand between the various physical quantities from the bulk quantum gravity (Liouville) and matrix model point of view. More details can be found in the corresponding chapters. Empty slots correspond either to the fact that there is no corresponding quantity, or that we do not yet understand the appropriate relation. Notice that the bulk quantum gravity theory can also be interpreted as a string theory on a target space. The acronyms used are DOS : for the density of states and SF F : for the spectral form factor.

Liouville theory
We begin by briefly reviewing some general facts about Liouville theory, focusing mainly in the c = 1 case. In this latter case, the Liouville action is to be completed with the action of one extra free bosonic matter field which we will label X(z,z). In course we will also delineate the points of departure of the usual interpretation of the theory as a string theory embedded in the linear dilaton background. The Liouville action on a manifold with boundaries is [23,24] (2.1) with K the extrinsic curvature and the parameters µ, µ B the bulk-boundary cosmological constants. This interpretation for these parameters stems from the fact that the simplest bulk operator is the area A = M d 2 z √ ge 2bφ and the simplest boundary operator is the length of the boundary = du g 1/4 e bφ(u) , with u parametrising the boundary coordinate of the surface. Notice that since we keep the conformal mode dynamical we are free to choose a gauge in which g = 1 where the scale factor is absorbed into the Liouville mode. One finds the following relations between the parameters (µ KP Z is the so called KPZ scaling parameter appearing in correlation functions) .
If this action is completed together with a c = 1 boson which we will denote by X(z,z), then c matter = 1 ⇒ b = 1, Q = 2 and one finds a renormalization of µ KP Z , µ B such that µ B = 2 √ µ cosh(πσ) becomes the correct relation between the bulk and boundary cosmological constants in terms of a dimensionless parameter σ.
To get into contact with the holographic picture, it is first important to discuss the various boundary conditions for the metric and matter fields φ and X and then the properties of the relevant boundary states. The matter field being a free boson, can satisfy either Dirichlet or Neumann conditions in the usual fashion. It is easy to see from (2.1), that the analogous possibilities for the Liouville mode are (n is the unit normal vector) 3) The Dirichlet boundary conditions φ| bdy = φ b are conformally invariant only asymptotically for φ = ±∞. In the limit φ → −∞ we have the weakly coupled region where the metric acquires an infinitesimal size and thus this is the regime in which we describe local disturbances of the bulk space. On the other hand for φ → ∞ (strongly coupled region of Liouville) distances blow up and we probe large scales of the bulk metric.
In addition, for a Holographic interpretation that is in line with the AdS/CF T correspondence, there ought to be a possibility for the bulk geometry to asymptote to AdS 2 . In fact this is precicely a solution of the Liouville theory equations of motion for which the asymptotic boundary is at φ → ∞ e 2bφ (z,z) = Q πµb This solution is that of the constant negative curvature metric on the Poincare disk. The metric is invariant under the Moebius transformations of the group P SL(2, R).
Let us now discuss more general solutions. On a quotient of hyperbolic space H 2 /Γ with Γ a discrete Fuchsian group, the general metric that solves the Liouville equations of motion is defined in terms of two arbitrary functions There exist three types of monodromy properties of A, B near non-trivial cycles of the manifold (three SL(2, R) conjugacy classes). For the hyperbolic one the curve surrounds a handle and this class can be used to describe higher topologies. We also have the elliptic monodromy class, which corresponds to surfaces on which the curve surrounds local punctures and the parabolic class for which the curve surrounds macroscopic boundaries. The solution (2.4) is of the parabolic class. Solutions of the elliptic class can be useful to describe the presence of singularities in the bulk of space 9 . The boundary of (2.5) is now at the locus A(z) = B(z), instead of the previous |z| = 1. Of course these more general metrics capture also the sub-case of Nearly-AdS 2 geometries [4] that correspond to slightly deforming the shape of the boundary. We notice that the important coordinate independent condition on the possible deformations that we allow (this holds even away from the strict Schwarzian limit), is that they keep the SL(2, R) conjugacy class near the boundary to be of the parabolic type so that one still finds a macroscopic boundary on which the holographic dual can reside. Another point of view for understanding such types of geometries, is to relate them to solutions to the bulk minisuperspace WdW equation [32] 10 In the expression above q is the momentum conjugate to the zero mode x of the matter field X(z,z) and thus a real number. In case the surface has a boundary of size , this can be expressed in terms of the zero mode φ 0 as = e φ 0 which is kept fixed. Notice that by fixing only the overall boundary size we are still allowing the possibility of other non-trivial non-zero mode deformations that change its shape.
The wavefunctions corresponding to loops with macroscopic sizes (at genus-zero) are given by These are (delta-function) normalizable wavefunctions with the norm which are exponentially damped for → ∞ and oscillate an infinite number of times for → 0. Different topologies of the bulk geometries are not described by these wavefunctions, but we will later see that the matrix model is able to resum the topologies automatically and express the wavefunctions in terms of integrals of Whittaker functions. On the other hand, the microscopic states correspond to wavefunctions that diverge as → 0 and vanish as → ∞, given by an analytic continuation of the previous solutions where ω is a real number. These are non-normalisable and correspond to local punctures (short distance bulk singularities). Notice that we can assign two different interpretations for such wavefunctions. One is from the point of view of the bulk surface and in this case they correspond to solutions with elliptic monodromy around a puncture. The second is to consider an analytic continuation for the matter field x → it with a subsequent interpretation of ω as a superspace frequency. In particular from the point of view of third quantisation, these are good asymptotic wavefunctions to use, if we wish to describe the scattering of universes in superspace, and in particular they probe the region near what we would call a "Big-Bang" or "Big-Crunch" singularity. Let us finally briefly discuss the boundary conditions for the matter field X(z,z) and the various types of D-branes that appear in Liouville theory [36,37,38]. As we noted these can be either Neumann or Dirichlet, which corresponds to the two possible choices of quantising matter fields on Euclidean AdS 2 . Neumann boundary conditions for X(z,z) correspond to ZZ boundary conditions and are relevant for describing D0 branes. Such branes are localised in the large φ → ∞ region. The FZZT-branes have Neumann conditions for φ (and either condition for X), stretch from φ → −∞ and dissolve at the region of φ ∼ − log µ B , so they can even penetrate the strongly coupled (large scale of the geometry) region depending on the value of µ B . This is for the fixed µ B ensemble (that corresponds to unmarked boundary conditions according to [8]). The relevant wavefunction (for c = 1) is If we instead perform a Laplace transform in order to keep the length of the boundary = e φ 0 fixed, we are then describing surfaces with fluctuating boundaries of fixed size. Using the relation µ B = 2 √ µ cosh(πσ), the fixed wavefunction is [27] This is again a solution to the minisuperspace WdW equation (2.6) upon identifying q = 2ν. This was given as an argument [23] for the minisuperspace description being exact (up to overall normalisations of the wavefunctions). Of course this WdW equation is exact only for the genus zero result, whilst the non-perturbative result given in section 4.1 can be related to a corrected version of the WdW that incorporates a term related to changes in topology. We conclude this section briefly mentioning the last type of D-branes, the Dinstantons [44,45] that correspond to Dirichlet conditions both in X(z,z) and φ(z,z). Their wavefunctions are labelled by two integers m, n (as for the D0-branes). These instantons are important for properly defining the vacuum state of the theory (here this is the third quantised vacuum) and were also argued to be related to fragmented AdS 2 spaces but we will not analyse them further here.

Matrix quantum mechanics and fermionic field theory
Let us now pass to the double scaled free fermionic field theory description of the model [33,31,32] that allows for an exact computation of various observables. After diagonalising the variables of matrix quantum mechanics and passing to the double scaling limit, the dynamics can be equivalently described in terms of the second quantised non-relativistic fermionic field action The double scaled fermi fields are defined using the normalised even/odd parabolic cylinder functions C ψ s (ω, λ), s = ± as (summation over the signs s, is implicit) where the fermi-sea vacuum |µ (µ is a chemical potential), is defined bŷ We should mention at this point that there exist various choices of defining the vacuum. In the old works the two common vacua are the one related to the bosonic string that has only one side of the potential filled, and the 0B vacuum that has both sides of the potential filled [39,40]. In the recent work [45] a new definition of the bosonic vacuum appeared in which there is no rightgoing flux from the left side of the fermi sea. These choices affect the non-perturbative physics of the model and as we will see, a fermi sea having support only on one side gives a non-perturbative WdW wavefunction that has an easier interpretation as a partition function of a Euclidean theory on the boundary of AdS 2 . On the other hand a fermi sea with both "worlds" interconnected via non perturbative effects, seems to be better suited to describe after an analytic continuation, wavefunctions of geometries of a cosmological type.
So far we described the theory in real time/energy t ↔ ω. We can also pass to the Euclidean description via t → ix, ω → −iq, but notice that this analytic continuation has a priori nothing to do with the bulk space/time notion of time. The natural interpretation of the matrix model time in our point of view is that of a time variable in superspace in which the universe is embedded. This means that in agreement with the discussion in the introduction, the natural interpretation of (3.1), is that of a third quantised action describing the evolution of two dimensional universes in superspace described by the coordinates (t, λ). Notice also that this action can capture the process of topology change for the bulk geometry (since the observables that one can compute using it are known to incorporate a resummed bulk genus expansion). At the same time this action is simply quadratic in the superspace field ψ(t, λ). This is quite interesting and unexpected since the third quantised superspace analysis in [13] indicated that one needs a non-linear modification of the WdW equation in order to describe topology changing processes. The reason of why this is possible is that the third quantised action is in terms of λ which is related via a complicated non-local transform to the metric conformal Liouville mode φ. From this point of view the condition that the bulk theory is a CFT, fixes the possible geometries of the superspace in which the bulk geometry can be embedded into.
Working now in Euclidean minisuperspace signature, we first define Matrix operators 1 N tr f (M (x)) and their fourier transform A simple example is the macroscopic loop operator (L is a discrete lattice variable) There is another description of the form [28] W This description keeps fixed a chemical potential z = µ B dual to the loop size. This is the boundary cosmological constant in the Liouville side as described in section 2.
In terms of the fermions, the most basic operator is the density operator In the double scaling continuum limit, we can employ the second quantised fermionic formalism to rewrite the Matrix operators in terms of the basic density operator. Since λ is the coordinate parametrising the matrix eigenvalues, the general relation is as followsÔ In particular the macroscopic loop operators (3.5) with length , have as a function f (λ, ) = e − λ . We can also describe local operators in terms of these, by shrinking the loop operators → 0 and dividing by appropriate powers of . A technical complication is that the support of the density ρ, depending on the non-perturbative completion of the model, can be on either side of the quadratic maximum of the inverted oscillator potential (−∞, −2 √ µ]∪[2 √ µ, ∞) so typically it is more convenient to consider the operatorsŴ Figure 1: The perturbative expansion of the WdW wavefunction in powers of g st ∼ 1/µ. The dashed lines indicate that we keep fixed only the overall size of the loop . and then wick rotate z = ±i for the various pieces of the correlator that have support in either side of the cut so that the corresponding integrals are convergent 11 . We will denote the expectation value of this loop operator as and so forth for the higher point correlators M n (z i , x i ). The details of the computation of such correlation functions are reviewed in appendices C and D. The results for the correlation functions of a compact boson X can be obtained from those of the corresponding non-compact case, via the use of the formulae of appendix E. We will now directly proceed to analyse the one and two point functions of loop operators.

One macroscopic loop 4.1 The WdW equation and partition function
We first start analysing the result for one macroscopic loop, or the non-perturbative WdW (Hartle-Hawking 12 ) wavefunction We first notice that this expression does not depend on q, the superspace momentum dual to the matter field zero mode x, in contrast with all the higher point correlation functions. This means that this wavefunction will obey a more general equation compared to (2.6), but for q = 0. To find this, we can compute the derivative of this integral exactly in terms of Whittaker functions Secondly, the functions that appear in this expression obey the WdW equation which is a generalisation of the minisuperspace Liouville result (2.6). The last term in particular was argued [33] to come from wormhole-like effects that involve the square of the cosmological constant operator ∼ ( e 2φ ) 2 . This is also consistent with the fact that the genus zero WdW equation (2.7), is missing precicely this term, while the exact result resummed the various topologies as shown in fig. 1. To get the genus zero result from the exact (4.1), one should perform a 1/µ expansion of the integrand and keep the first term to obtain the genus zero or disk partition function While the genus zero wavefunction has an exponential decaying behaviour for large 13 , the exact wavefunction has an initially fast decaying behaviour that transitions to a slowly decaying envelope with superimposed oscillations, see fig. 2. This can be also studied analytically by employing a steepest descent approximation (see appendix F) to the integral (4.1). In the first quadrant of the complex ξ plane the integrand vanishes exponentially fast at infinity and hence we can rotate the contour at will in the region ξ > 0. The steepest descent contour goes actually along the direction of the positive imaginary axis, so we should be careful treating any poles or saddle points. The poles at ξ = i2πn are combined with a very fast oscillatory behaviour from the exponent due to the factor coth ξ/2 14 . There are two types of leading contributions to the integral as → ∞. The first is from the region near ξ = 0. The integral around this region is approximated by the genus zero result (4.4) and decays exponentially. In order to study possible saddle points, it is best to exponentiate the denominator and find the saddle points of the function There is a number of saddle points at u * = Arc [sin( 2 )] + 2nπ or u * = −Arc [sin( 2 )] + (2n + 1)π. The integral along the steepest descent axis is now expressed as To get the leading contribution we just need to make sure that our contour passes through the first saddle point u * ( ). From that point we can choose any path we like in the first quadrant, since this only affects subleading contributions. An obvious issue is that this is a movable saddle point since it depends on , which we wish to send to ∞. To remedy this according to the discussion of appendix F, we define u = u * u and perform the saddle point integral to find the leading contribution This expression is in an excellent agreement with the numerical result plotted in fig. 2. The amplitude decays as 1/ √ log 2 and the phase oscillates as e i −2iµ log asymptotically for large .
We can also compute exactly the wavefunction with the insertion of a local operator V q , the result employing the exact wavefunctions for general η (4.3) This reduces to the genus zero result in accordance with the Liouville answer for the wavefunctions describing microscopic states (punctures of the surface). Multiple insertions can be computed with the use of the more general formula for correlation functions (D.6) by shrinking the size of all exept one of the loops and picking the appropriate terms. We close this subsection with some remarks on the wavefunctions and the consequences of the identification Ψ W dW ( ) = Z dual (β = ). The non-perturbative wavefunction is of the Hartle-Hawking type (see appendix B and especially eqn. (B.12)). It is real by construction and exhibits a decaying behaviour at small indicative of being in the forbidden (quantum) region of mini-superspace 16 . The genus zero large decaying behaviour is indicative of describing a space that reaches the Euclidean vacuum with no excitations when it expands to infinite size. This holds at any fixed genus truncation of the non-perturbative result. On the contrary, for large the non-perturbative fast oscillatory behaviour is indicative of a semi-classical space interpretation (see again the appendix B for some discussion on that). This is reasonable since large geometries are indeed expected to have a semi-classical description, while small geometries are highly quantum mechanical. Quite remarkably the same behaviour also appears in the case of positive cosmological constant upon continuing = −iz and we will comment on this unexpected relation in section 6. This is an indication that the non-perturbative wavefunction does not correspond to the Euclidean vacuum when the geometries are large, but has a more natural interpretation in a cosmological setting.
On the other hand, if we are to interpret the wavefunction as the finite temperature partition function of a dual system, we run into the difficulty of having an oscillatory partition function as a function of the dual temperature β = for small temperatures T ∼ 1/ . This seems to be related to the well known problem of the difficulty of assigning a probabilistic interpretation to the WdW wavefunction. Related to this, as we will see in the next section 4.2, the dual density of states ρ dual (E) is manifestly positive definite, but the spectral weight has support both on E > 0 and E < 0 and is an even function of E. The non-perturbative effects are precisely the ones that make these "two worlds" communicate. We can remedy this by fiat, demanding the spectral weight to have support only at positive energies E > 0. We will describe the consequences of this choice for the partition function in subsection 4.2.2. The infinite temperature limit → 0 in (4.1) is singular, (needs a UV-cutoff Λ corresponding to putting the inverted oscillator in a box). Nevertheless such a limit still makes sense physically, from a third quantised point of view. The reason is that the geometries then smoothly close and reproduce a closed string theory partition function of compact manifolds [31].

The density of states of the holographic dual
We will now compute the density of states of the dual quantum mechanical theory, assuming that the renormalised wavefunction of loop length corresponds to the dual thermal partition function according to our interpretation. This also dovetails with the fact that the loop length is related to the zero-mode of the conformal mode of the metric = e φ 0 and hence its dual variable is the trace of the boundary theory stress-tensor (energy in the case of a purely quantum mechanical system). Since the length = β of the loop is related to the inverse temperature of the dual field theory, we can define the density of states through the inverse Laplace transform 17 It is more clarifying instead of directly performing the inverse Laplace transform (or fourier transform in terms of z) in the final expression for the wavefunction (4.1), to go back to the original definition of the Loop operator (3.9) and Laplace transform it. This then results into the following expression for the density of states in terms of parabolic cylinder wavefunctions (see also appendix C.4) (4.10) We therefore observe the remarkable fact that the energy of the dual field theory E corresponds to the fermionic field theory coordinate λ, since it is a conjugate variable to the loop length . In addition the natural interpretation of the parameter µ from this point of view is that of an IR mass gap to the energy spectrum (since the inverted oscillator potential provides an effective cutoff [2 √ µ, ∞) to the allowed range of the eigenvalues). In more detail, by expanding the parabolic cylinder functions in the region E √ µ, the eigenvalue distribution is found to follow a smooth Wigner-Dyson envelope on top of which many rapid oscillations of small amplitude are superimposed This behaviour is similar to the large-N limit of a random Hermitian matrix Hamiltonian, but the small rapid oscillations is an indication that we have also incorporated some additional non-universal effects. Moreover, the oscillations become more pronounced in the limit µ → ∞ and diminish as µ → 0 when the sum over topologies breaks down. In addition, there is a tail of eigenvalues that can penetrate the forbidden region outside the cut [2µ 1 2 , ∞). In this region E 2 µ, there is an exponentially growing density of states This is a Hagedorn growth of states, but only for a small window of energies. In other words it is important that there is a transition from an exponential to an algebraic growth, that makes the density of states Laplace transformable and the WdW wavefunction well defined. In fig. 3 the complete behaviour of the density of states is depicted. The exact integral (4.10) can be further manipulated to give an integral expression (c is an infinitesimal regulating parameter) This then matches the fourier transform of (4.1) as expected. The complete spectral weight ρ dual (E) is manifestly an even function of E (as well as positive definite). It would be very interesting to try to give a Hilbert space interpretation for this density of states but we will not attempt that here. The only comment we can make is that since we also have the presence of negative energy states, the dual boundary theory needs to have a fermionic nature so that one can define a Fermi/Dirac sea and a consistent ground state.

Comparison with minimal models and JT gravity
Let us now compare this result with the density of states found in the Schwarzian limit of the SYK model dual to JT gravity, as well as with the result for the minimal models. The first point to make is that the growth of states near the edge of the support of the spectrum is in fact faster than that of JT gravity, since where now γ provides an energy scale and in the exponent one finds merely a square root growth with the energy.
In fact we can also make a further comparison, using our perturbative expansion in 1/µ of the exact density of states (4.10), with the ones related to minimal strings. This is because the inverse Laplace transform of the genus zero result (4.4) and in fact of every term in the perturbative expansion of the exact partition function (4.1) can be performed via the use of the identity ( This gives a spectral curve for each genus that is very similar to the ones related to the minimal strings, discussed in [7]. Nevertheless none of them can capture the non-perturbative ∼ e −µ effects that give rise to the oscillations both in the exact partition function and density of states plotted in 3. These are also effects that make the two worlds E ≶ 0 communicate through eigenvalue tunneling processes in the matrix quantum mechanics model. Similar non-perturbative effects were also discussed in the context of JT gravity and SYK [8]. In that case it has been argued that they are doubly non-perturbative in exp(c e N SY K ), while in the present example this could arise only if µ was a parameter having a more microscopic description (as happens in the model of [50] 18 ). We will come back to a discussion of a possible interpretation of such a double layered asymptotic expansion from the point of view of the bulk theory in the conclusions section 7.

The one sided Laplace transform
As we mentioned previously, we can define the dual density of states to have support only for E > 0. This is equivalent to demanding that no eigenvalues can penetrate to the other side of the potential. In such a case the dual partition function is given by the Laplace transform of (4.10), with support at E > 0 so that This is a positive real function exhibiting a decaying behaviour with no large amplitude oscillations 19 . By restricting the spectral weight to positive energies ρ (+) dual (E), the WdW wavefunction does acquire a well defined probabilistic interpretation. This modification though, does not seem to make sense for the analytic continuation z = i to the cosmological regime of section 6, since the wavefunction still remains non oscillatory and decays for large . We believe that this is an indication that nonperturbative microscopic models of AdS 2 and dS 2 geometries, can be very different.

The case of two asymptotic regions
We will now analyse observables in the presence of two asymptotic boundaries. One should include both disconnected and connected geometries. We focus mainly in the density two-point function and in the spectral form factor (SFF) related to its fourier transform.

Density two-point function
We first analyse the density two-point function µ|ρ dual (E 1 )ρ dual (E 2 )|µ . A quite thorough discussion of similar correlation functions from the point of view of quantum gravity can be found in [6]. This correlation function is defined in appendix C.5. The disconnected part is given by the product of (4.10) with itself.
The genus zero result can be computed analytically and is missing any oscillatory behaviour. A plot is given in fig. 4. The non-perturbative result is plotted in figs. 4 for the short spacing behaviour δE = E 1 − E 2 of energy eigenvalues, and 5 for larger spacings. It has the behaviour of the sine-kernel indicative of short range eigenvalue repulsion and chaotic random matrix statistics for the eigenvalues. For an easy comparison we have also plotted the sine kernel on the right hand side of fig. 5. We observe a qualitative similarity, but there do exist differences that will become more pronounced in the SFF leading to a slightly erratic oscillatory behaviour. These results combined with the ones of section 4.2, indicate that if we would like to endow the boundary dual with a Hilbert space and a Hamiltonian, its spectrum is expected to be quite complicated and resemble those found in quantum chaotic systems.

Spectral form factor due to disconnected geometries
Another interesting quantity we can compute is the spectral form factor (SFF), first due to disconnected bulk geometries. This corresponds to the expression In the genus zero case we can use the analytic expression (4.4) to compute it. This results in a power law ∼ 1/t 3 decaying behaviour at late times, with a finite value at t = 0, due to the non-zero temperature β.
The non perturbative result can be computed at three different limits using a steepest descent analysis. The first is for µ t, β that is equivalent to the genus zero result. The early time limit is for β t, µ that gives a result that is again very similar to the genus zero answer. The last is the late time result for t β, µ which is plotted on the right hand side of fig. 6. The decay in this case, has the scaling behaviour ∼ 1/t log 4 t as t → ∞.

Euclidean wormholes and the loop correlator
We now pass to the case of the connected loop-loop correlator Ŵ ( 1 , q)Ŵ( 1 , q) = M 2 ( 1 , q, 2 , −q). From the bulk quantum gravity point of view according to [7,8], this observable corresponds to a correlator of partition functions that does not factorise. For the SYK model this happens due to the intrinsic disorder averaging procedure when computing this observable. In non-disordered theories with complicated spectra it was argued that this could arise from Berry's diagonal approximation [10] that effectively correlates the two partition functions, even though the exact result does factorise. On the other hand in [19] it was proposed that such multi-boundary geometries could in fact correspond to a single partition function of a system of coupled QFT's. A previous work by [18] also considered such geometries and discussed various possibilities for their interpretation such as that the dual partition function is not well defined, or that at the full non-perturbative level of the quantum gravity path integral, such observables indeed factorise 20 . Since there is no ultimate resolution given to this question yet, it is of great importance to analyse such correlators in the present simple example where we can compute them (in principle) at the full non-perturbative level.
In particular for two loops we find the following expression for the derivative of the correlator ∂M 2 (z 1 , q, z 2 , −q)/∂µ [33] The prescription for analytic continuation is z i → i i together with → i/2. One can also express the second integral over s as an infinite sum of Bessel functions giving 21 4i n n n 2 − q 2 J n (x) sinh(nξ/2) , (5.3) with x = z 1 z 2 / sinh(ξ/2). This expression can be used to extract the genus zeroresult. In addition the integral (5.2) does have a nice behaviour for large ξ, and all the corresponding integrands vanish for large ξ exponentially. A similar property holds for large-s for each term of the s-integral at the corresponding quadrant of the complex-s plane. One can also pass to the position basis q ↔ ∆x via the replacement e −|q|s ↔ s π(s 2 + (∆x) 2 ) (5.4) These two bases reflect the two different choices for the matter field X, either Dirichlet x =fixed, or Neumann q = fixed at the boundary and hence to the two basic types of correlation functions. At genus zero the expression for the correlator simplifies drastically. In particular it can be written in the following equivalent forms The expressions above elucidate different aspects of this correlation function. The first expression has an interpretation in terms of a propagation of states between macroscopic boundary wavefunctions Ψ macro p ( ) (2.7). The result can also be expanded in an infinite sum of microscopic states as the second line indicates. The final expression is a superposition of single wavefunctions corresponding to singular geometries (microscopic states). It also shows that there might be an interpretation for which the complete result corresponds to a single partition function of a coupled system. It would be interesting to see whether the exact result (5.2) can be manipulated and written in a similar form, for example using the exact wavefunctions (4.3) or (4.7). We now turn to the study of the spectral form factor arising from such connected geometries.

Spectral form factor due to connected geometries
In this subsection we analyse the part of the SFF corresponding to a sum over connected bulk geometries. According to the discussion in appendix C.5, for the unrefined SFF it is enough to study the limit q → 0 of the more general expression for the momentum dependent loop correlator, eqn. (5.2). Before doing so, we first define the parameters of the spectral form factor through z 1,2 = i(β ± it) We can then distinguish the three basic timescales: t β, t ≈ β and β t. We will denote them as early, median and late time-scale respectively. The spectral form factor can then be expressed as a double integral using (5.2) as ∞ 0 ds ∞ 0 dξ ξ sinh(ξ/2) e iµξ−i(β 2 −t 2 ) coth(ξ/2)−i(β 2 +t 2 ) coth(ξ/2) cosh s sin (β 2 + t 2 ) sinh s .  The genus zero result is shown in fig. 7 and is found to be a time independent function. The graph can also be obtained by directly integrating (5.5) for q = 0. One can notice that the g = 0 SFF captures the plateau behaviour, unlike the case in [8] where the g = 0 SFF captures the ramp behaviour, for t >> β. The plateau behaviour arises due to the repulsion among neighboring energy eigenvalues. On the other hand, the ramp behaviour arises due to the repulsion among eigenvalues that are far apart 22 [6]. The difference between our case and the one in [8] is not surprising, since in our case the genus zero part is obtained in the limit µ → ∞ where effects involving eigenvalues that are far apart are suppressed. The double integral describing the exact result contains a highly oscillatory integrand, which needs further manipulation so that it can be computed numerically with good accuracy using a Levin-rule routine. It is useful to perform the coordinate transformation u = coth ξ/2, that effectively "stretches out" the oscillatory behaviour near ξ = 0. The numerical result is then contrasted with the genus zero answer in fig. 7 and results in a ramp -plateau like behaviour with erratic oscillations near the onset of the plateau, that become more regular at late times as seen in fig. 8. At late times the oscillations are of the same order as the function itself: ∆SF F c (t)/SF F c (t) → O(1), t → ∞. This is an indication that the boundary dual is a theory with no disorder averaging.
We have also studied a refined SSF, or better said the q = 0 correlator. Unexpectedly, its behaviour is qualitatively different and much smoother from the q = 0 case, even for small values such as q = 0.1. The genus zero result is shown in fig. 9. It exhibits a slope at early times, that very slowly reaches a plateau, see the left side of fig. 9.
On the other hand the non-perturbative connected correlator after the initially decreasing slope behaviour, transitions into a smooth increasing one that at late times again relaxes to a plateau.

Comments on the cosmological wavefunctions
In this section we analyse the possibility of giving a dS 2 or more general cosmological interpretation for the WdW wavefunctions, after discussing the various possibilities for analytically continuing the AdS 2 results 23 .
One such analytic continuation in our notation and language, is going back to the parameter z = i in section 3, and using the fourier transformed operators to compute the partition function and correlators. In [21], the authors explained why this analytic continuation describes "negative trumpet" geometries by analysing how the geometries change in the complex field space. Let us first define the dS 2 global metric ds 2 = −dτ 2 + cosh 2 τ dφ 2 , (6.1) and consider then the case of complex τ , so that we can describe Euclidean geometries as well. The usual Hartle-Hawking [14] contour for dS 2 involves gluing a half-sphere to dS 2 . This is indicated by the blue and black lines in fig. 11. On the other hand the geometries obtained by the continuation z = i are "negative trumpet" geometries that again smoothly cap-off much similarly to what happens in the no-boundary geometries of Hartle and Hawking. These are indicated via the red line in fig. 11. Even though these are not asymptotically dS 2 geometries, nevertheless they can be used to define an appropriate no-boundary wavefunction Ψ W dW (z = i ) = tr e izĤ , whereĤ is now the generator of space translations at the boundary. In addition according to fig . 11, one can reach the same point in field space (describing a large dS 2 universe), either through the usual Hartle-Hawking prescription for the geometries, or according to a different contour that passes through the "negative trumpet" geometries which then continues along the imaginary axis so that it connects them to the dS 2 geometry 24 . One can also imagine the presence of obstructions, in the sense that the two paths in field space might not commute and therefore give different results for the wavefunction. In fact this is precisely what happens in the present example, for the genus-zero wavefunctions. The mathematical counterpart for this, are the properties and transitions between the various Bessel functions as we analytically continue their parameters.
In order to clarify this further, we should also mention a slightly different approach of analysing bulk dS 2 geometries proposed in [26] and [27]. In the latter case the author performed an analytic continuation of the Liouville theory: b → ib, φ → iφ, resulting in the supercritical case for which c ≥ 25. The appropriate states describing the dS 2 geometries (the conformal boundary is still at φ → ∞), are given by the analytic continuation of those in (2.7) and result into the dS 2 Hankel wavefunctions Ψ iq (z). These wavefunctions are disk one-point functions that describe an asymptotically large dS 2 universe that starts at a Big-Bang singularity (whose properties are determined via the vertex operator inserted at the disk). If there is no operator insertion one simply obtains the corresponding Hartle-Hawking wavefunction.
In the present c = 1-model at the level of the genus zero wavefunctions, the approaches of [21] and [27] do coincide, and give the same type of Hankel functions as the appropriate dS 2 WdW wavefunctions. This can be seen using the formula that holds in particular for z ∈ R. If we apply this formula to the genus zero ... Figure 12: The summation of all the possible smooth Euclidean geometries that asymptote to a "trumpet" geometry.
Its plot is shown in the left panel of fig. 13. At the non-perturbative level, we can analytically continue the wavefunction of eqn. (4.1). This result is resumming the geometries plotted in fig. 12. A plot of this non-perturbative wavefunction is given on the right panel of fig. 13. The behaviour is again qualitatively similar to the genus zero Hankel functions but with more rapid oscillations. This is in contrast with the Euclidean AdS 2 case, leading to the conclusion that the most reasonable non-perturbative choice describing the dual of the AdS 2 geometries is the one having only one side of the potential filled as described in subsection 4.2.2, whilst the twosided fermi-sea is better suited for describing cosmological types of geometries. This is also in line with the fact that in Ψ W dW (z = i ) = tr e izĤ ,Ĥ is performing space translations and there is no constraint on the positivity of this trace.
So far we discussed a very particular analytic continuation. Nevertheless as we mentioned above, the story is more complicated and there exist various other physical choices in choosing an appropriate wavefunction. In particular at genus zero, all the rest of the Bessel functions appear and play the role of various other choices for states depending on the boundary conditions imposed [26]. These possibilities are summarised here: • Ψ n.b. ∼ J ν (z), ν ∈ R is the Hartle-Hawking no-boundary wavefunction. It is real, goes to zero at small z and oscillates at large z (see appendix B for more details on the no-boundary wavefunction).
ν (z) corresponds to the tunneling proposal. In particular it is complex and increasing at small z (see appendix B for more details on the tunneling wavefunction).
• Finally there is the option of demanding only expanding universes at small scales with the resulting wavefunction Ψ exp. ∼ J iν (z), ν ∈ R. This has a more fitting minisuperspace interpretation as giving the "birth of a superspace quantum" [26], that can be used in a "scattering process in superspace".
A natural question regards the possibility of realising all these wavefunctions as various limits of a complete non-perturbative description. This is what seems to be happening, since the more general Whittaker WdW wavefunctions appearing in (4.3), (4.7), can be related to all the various types of Bessel functions as we send µ → ∞, along different regions of the complex µ, z, q planes. We then expect the presence of interesting Stokes phenomena and transitions between the various asymptotic genus zero wavefunctions. We plan to revisit this problem in the future.

Conclusions
We will now comment on our findings and discuss some interesting future directions. We have discussed MQM from the point of view of selecting a non perturbative completion of the c = 1 Liouville bulk quantum gravity path integral (which is not unique). The bulk quantum gravity path integral seems not to factorise even at the full non-perturbative level and hence Euclidean wormholes and the physics of multiple boundaries are inevitable consequences of this "third quantised" point of view.
The Schwarzian alone is known not to be a consistent quantum mechanical model, since the path integral on the circle cannot be interpreted as tr e −βĤ for anyĤ [46,47]. On the other hand the boundary duals of Liouville theory coupled to matter, have a possibility of being consistent quantum mechanical theories of their own. In the present c = 1 case, this is a problem that relies on a Hilbert space interpretation of the exact function for the density of states given by eqn. (4.10) or equivalently of the dual resolvent of appendix C.4 (the usual resolvent can be interpreted in terms of the inverted oscillator Hamiltonian). Such a system could have both a continuous and a discrete spectrum, the small persistent oscillations being indicative of a microscopic discrete part. Since the density of states resembles that of large-N random matrices, any dual system is expected to have approximate chaotic properties. This is also evident from an analysis of the exact connected SFF of eqn. (5.7), displaying a ramp plateau behaviour with the presence of erratic oscillations.
An important property of the spectral form factor is that the ramp and plateau are not self-averaging. It has been argued that for a fixed Hamiltonian system, even though we are summing over many energy levels, the result should be a function with O(1) fluctuations. The smooth ramp and plateau are results of a time or disorder average, but the exact function is expected to be erratic [7]. In our example we have remnants of such an erratic behaviour even in the double scaling limit. While they are not as pronounced as the ones expected in higher dimensional examples they are still present and affect both the DOS and the SFF pointing to a non-disorder averaged dual theory.
It is also interesting to notice that while there exists a complete description of the bulk path integral in superspace in the form of MQM, a single boundary dual theory alone does not seem to be able to capture all its intricacies. For example, while the interpretation of Euclidean wormholes is rather straightforward from a superspace perspective, it is not obvious how and if this information is encoded in a single boundary dual theory. Perhaps the most reasonable attitude on this, is that one needs to first choose the types of boundary conditions in superspace that are allowed and depending on this choice there could be a single boundary dual description. The reason is that from a third quantised point of view the dual partition function and density of states correspond to operators, and one can form various expectation values out of them. The complete Hilbert space is hence naturally larger than naively expected. Nevertheless this is not an argument against holography. Holography simply seems to hold slightly differently than expected, involving matrix degrees of freedom and an enlarged third quantised Hilbert space. By imposing certain restrictions on that complete Hilbert space we can then use a more usual holographic correspondence.
A preliminary analysis of the cosmological regime described after the analytic continuation z = i in section 6, reveals that while there is a single WdW wavefunction that takes the form of an integral of a Whittaker function, there exist various asymptotic expansions/limits of the parameters of the wavefunction leading to interesting Stokes phenomena. This is the reason behind the appearance of several types of semi-classical wavefunctions in superspace and boundary conditions, such as the Hartle Hawking or the no-boundary proposal etc., while all of them seem to arise from a unique progenitor in this simple two dimensional example. This is clearly a point that deserves further study.
Let us also mention that we could also have considered other types of "target space" minisuperspace geometries, that would correspond to other types of matter content from the two dimensional bulk point of view. For example one could try to analyse WZW-models in the presence of boundaries. It would be interesting to see what are the differences with the present example.
Worldsheets inside worldsheets and double layered expansions -Another question that we briefly alluded to in the introduction, is the possibility of having a double layered expansion, or what one could call "doubly non-perturbative effects". One natural quantum gravitational setting in which such an expansion could arise is the idea of having "worldsheets inside worldsheets" or "geometries inside geometries". For example, it would be interesting to understand whether the two dimensional Liouville theory worldsheets (quantum gravity target spaces from our WdW point of view) can emerge themselves from an underlying string theory in a fashion similar to the proposal in [54] 25 . This would be based on the sequence of mappings σ,σ → z,z → X, φ and result in a third quantised theory of universes inside of which strings propagate. That said, the resulting theory should a priori have a two parameter expansion: the genus expansion of the strings g s as well as the topology expansion of the resulting target space, which in our interpretation was the 1/µ expansion. Then for a fixed topology of the target space, one would have to sum over all the string worldsheet genera. The doubly non-perturbative effects then take the form exp(−µe −1/gs ), since one can introduce both boundary branes for the space as well as worldsheet boundaries.
It is much harder to concoct a dual model that can realise such a setup. If there is any matrix/tensor model realisation of this idea, it should involve two free parameters in the appropriate continuum scaling limit. While the tensor models naturally exhibit similar rich scaling limits [58], they are much harder to study. Something analogous is also realised in rectangular N × M matrix models [56,57], where one can tune the parameters M, N → ∞ together with a coupling constant g independently, to reach new critical points. There is an obvious problem in this idea related to the fact that when the shape of the matrices becomes very narrow the surfaces degenerate into branched polymers. To overcome this obstacle, we can use several large-N limits in conjuction with the idea of breaking the original Hermitean matrix M ij into blocks, for example In this splitting there is a breaking of the U (N ) symmetry into U (n)×U (N −n), when the off-diagonal elements are zero, that results in two distinct geometrical objects (surfaces) if we take the double scaling limit in both of them. On the other hand when N ≈ n the off-diagonal elements become very narrow (leading to a surface interacting with a particle through branched polymers in the scaling limit). This then means that we can try to break the original matrix into n-blocks N = N 1 + N 2 + ...N n and introduce a chemical potential µ N for N and µ n for n. While µ N governs a genus expansion in superspace as before, the parameter µ n playing the role of a genus expansion for any fixed value of µ N . More hierarchical matrix embeddings could lead into analogous hierarchical surface embeddings. Such embeddings also always satisfy a version of the stringy exclusion principle [60]: submatrices are always smaller than the matrix they are embedded into.
Higher dimensions -There is a crucial difference of the present models with higher dimensional examples. All the theories in two dimensions could be related to some limit of string theories with some particular form of matter content. For example it is not clear whether we can consistently define higher dimensional models of geometries propagating in some form of superspace. This is a first crucial point to understand if one wishes to extrapolate the present results into higher dimensions, and makes even more pertinent the analysis of higher dimensional examples even at a truncated mini-superspace level. The most probable conclusion consistent with the present results is that the bulk quantum gravity theory is in fact richer than that of a single quantum field theory (in the present case boundary quantum mechanics). Imposing certain restrictions and boundary conditions at the superspace level, one can indeed reduce the path integral in a subsector that is dual to a single QFT as in the usual implementations of holography. In the most general case though the non-perturbative bulk path integral can compute more complicated objects than that of a single QFT partition function, that might or not have a single QFT interpretation. This is not in contradiction though, either with the holographic interpretation at fixed genus given in [19], or with the possibility that there could exist a complicated matrix/tensor model, such as MQM in the present case, that is able to capture all the non-perturbative quantum gravity information. Such a model would not be the boundary dual of the bulk quantum gravity theory in the usual holographic sense, but a model describing directly the dynamics of geometries in superspace. In fact variants of the BFSS matrix model [61] could very well have the potential to serve as a dual of this kind. acknowledge the hospitality provided by APC and ENS Paris during the initial stages of this work. This work is supported in part by the Advanced ERC grant SM-GRAV, No 669288.

A. Matrix models for minimal models A.1 Conformal maps and Integrable hierarchies
There exists a well studied relation of contour dynamics in two dimensions with the dispersionless limit of the Toda Hierarchy [51]. We now briefly review the results of these works and then discuss their relation to the two dimensional quantum gravity path integral.
An equation for a curve on the complex plane F C (z,z) = 0 can be resolved locally with the help of the Schwarz function asz = S(z). We assume this curve to bound a simply connected domain D + and we label the exterior domain with D − . The problem of multiple domains is described in [52] (in terms of the Schottky double). The Schwarz function obeys a unitarity conditionS(S(z)) = z and can be decomposed into two functions S(z) = S + (z) + S − (z) that are holomorphic in the interior/exterior of the domain. Let us also define a conformal map z(w) that maps the exterior of the unit disk to the exterior domain D − .
We then define a function Ω(z) via S(z) = ∂ z Ω(z). It plays the role of the generating function of the canonical transformation from the unit disk to the region bounded by the curve. Its differential defines a multi-time Hamiltonian system through In this formula and t k are "time variables" corresponding to the moments of the region and t is the zero time dual to the area, for more details see [51]. This indicates a dual way of describing the curve via the moments C(t, t k ,t k ) and a dual "prepotential" F (t, t k ,t k ). The relevant set of equations that governs this system corresponds to the dispersionless limit of the Toda hierarchy → 0. In this dispersionless limit, the τ function is then defined as τ = exp(F/ 2 ) and the Baker-Akhiezer wavefunction as Ψ = exp(Ω/ ). A particular representation of a τ -function of the Toda hierarchy is in terms of a two-matrix model The interface dynamics is then described by the dispersionless limit of this matrix model i.e. = 1/N → 0, which is the large-N limit. Its free energy is thus the prepotential F (t k ,t k ). An equivalent definition of the τ -function is encoded in the Schwarzian derivative of the conformal map w(z) through the following relation A.2 The (p, q) minimal models The matrix model (A.3) describes also the dynamics of the (p, q)-minimal models coupled to gravity [41,42,43]. The (p, q) minimal models arise when the square of the Liouville parameter becomes a rational b 2 = p/q. In this series of works it was understood that all the possible (p, q) models can be described in terms of a M (p,q) Riemann surface at the perturbative level. Let us first define the dual cosmological which describe a symmetry of the physical observables under b → 1/b. If we then use the parameters with Z F ZZT the disk partition function of the FZZT-brane, the duality means that x = y andỹ = x. The set of equations that determine the partition functions is the second integral corresponds to an integral through the pinched cycles of M (p,q) . This is where the ZZ-branes reside. An equivalent way of rewriting all these equations is through This means that the function F describes a curve corresponding to M (p,q) . This surface has (p − 1)(q − 1)/2 singularities when that correspond to pinched cycles of the surface where the ZZ-branes reside. The presence of a background of such branes has the effect of opening up the pinched cycles. At the non-perturbative level in the string coupling, Stokes phenomena change the picture drastically and the surface M (p,q) is replaced by the simple complex plane C [42]. For example the exact FZZT partition function is described by an Airy function Ai(x + 1) which is an entire function of x. In order to understand precicely how this happens, the matrix model comes at rescue, since one can compute the appropriate loop operator expectation values. The loop operator is defined through In particular at genus zero W (x) = Z F ZZT corresponds to the disk amplitude and y = ∂ x Z F ZZT to the resolvent of the matrix model. The full non-perturbative FZZT brane corresponds to an exponential of the loop operator It is known that the expectation value of the determinant operator, in the double scaling limit that corresponds to forming continuous surfaces, corresponds to a Baker-Akhiezer function that is an entire function of x and hence the complete non-perturbative moduli space of the FZZT branes is the complex plane C.

A.3 Deformations
We can also turn on an arbitrary number of closed string couplings t k ,t k to deform the closed string background as seen from the matrix model (A.3). The notation here is due to the fact that the deformed quantities such as the partition functions are generically related to τ -functions of the KP (a single set of "times") or Toda hierarchies with t n 's playing the role of "times".
The new differential to be integrated is then defined via a deformation of ydx: where H m (x) are mutually commuting "Hamiltonians" dual to the time variables t k and so forth for the bar quantities. In particular one can now replace the previous differential ydx by dΦ in all the formulae (A.7). Notice the equivalence to the Hamiltonian system A.1 that allows to transition between the physics of the minimal models and that of interface dynamics. We therefore conclude that the dispersionless limit captures a universal sector of the interface dynamics as well as of (p, q)-minimal models coupled to gravity. We can further consider a small sector of the "Goldstone hydrodynamic modes" that describe small ripples of the boundary curve (interface) geometry. Nevertheless there is a huge class of integrable-deformations governed by the time parameters t k ,t k , that capture finite deformations of the geometry.

B. Properties of the WdW equation
In this appendix we will describe some basic properties of the Wheeler DeWitt equation as well as the most common boundary conditions employed in the literature. Some reviews can be found in [13]. The complete superspace WdW equation is not really well defined and suffers from various problems such as infinite configuration space, operator ordering ambiguities, action unbounded from below etc. We will hence restrict to a minisuperspace formulation of the problem which reduces the degrees of freedom to a finite number, fortunately the two dimensional example of Liouville theory is very well described just by the minisuperspace approximation due to the small number of physical degrees of freedom in two dimensions.
The general minisuperspace action takes the form where q a (r) are the finite number of variables, N (r) is a non-dynamical Lagrange multiplier and G ab (q) is a reduced form of the metric in superspace. The momentum constraints are satisfied trivially in the minisuperspace ansatze, so what is left is the Hamiltonian constraint coming from the Lagrange multiplier N (r) Upon quantising one should replace the momenta with but generically there is a non-trivial operator ordering problem because of the metric G ab (q) on minisuperspace. Assuming a freedom of field redefinitions, the appropriately ordered operator takes the form In this expression the covariant Laplacian is computed using the minisuperspace metric. This minisuperspace metric has generically indefinite signature and hence one can find both exponential and oscillatory solutions. Let us also note that there is a path integral representation of the wavefunction. First there is a reparametrisation symmetry due to N (r). It can be shown that only the zero mode plays a role and the path integral takes the form (in the Lorentzian case) Ψ(q a rmax ) = dN Dp a Dq a e iS L (pa,q a ,N ) = dN Ψ(q a rmax , N ) (B.5) where the boundary conditions are q a (r max ) = q a rmax and the other variables are free. There are also extra boundary conditions to be specified at any other boundary of the r variable. The wavefunction Ψ(q a rmax , N ) satisfies the time dependent Schroendinger equation with N playing the role of the time variable. In particular so that we need to take either an infinite contour in the N -space or a closed contour. Some usual physical requirements on Ψ(q a rmax ) employed in the literature are that it is peaked near the classical configurations, and the interference between two configurations is small (decoherence).
To pass over to the classical limit it is convenient to search for a semi-classical WKB ansatz for the wavefunction of the Hamilton-Jacobi form Ψ = e −I R (q)/ +iS(q)/ , with I R (q), S(q) real. Plugging this inside (B.4) one then finds the set of equations The second equation is equivalent to the usual condition of steepest descent analysis: the steepest contours of the real part are orthogonal to those of the imaginary part.
As for the first equation, in the regions where U (q) > 0 we find The subscript c.a stands for classically allowed region of the minisuperspace and c.f for classically forbidden or tunneling. The classically forbidden case exhibits a tunneling wavefunction, for which quantum effects are important. On the other hand the classically allowed oscillatory wavefunctions are strongly peaked on a set of classical solutions of the Hamilton-Jacobi equation (B.7). In addition if |∇ G I R | |∇ G S| the amplitude is changing very slowly and the classical approximation becomes extremely accurate. Let us now discuss some caveats that one should be careful about. A first complication comes from the fact that the metric in minisuperspace is not positive definite (due to the conformal mode of the metric) and generically the boundary value problem is of the hyperbolic type. Another issue related to this is the choice of contour in mini-superspace which is quite subtle. The most common ones are: The no-boundary proposal that passes from the Euclidean regime and posits a smooth Euclidean completion of the boundary geometry, it is CPT invariant and consists out of two WKB modes, so that it corresponds to a real wavefunction.
The tunneling wavefunction that posits an outward probability flux near the singular boundary of superspace that gives only one WKB mode 26 . The probability current cannot be defined globally but for the WKB approximation it takes the form J ∼ −∇S.
Let us now give a few more details about these two most common options.
The tunneling wavefunction can be described is the classically forbidden/allowed regions via which in the semi-classical approximation becomes The no-boundary proposal of Hartle-Hawking, consists of the following solutions in the classically forbidden/allowed regions 12) or in terms of the semi-classical variables Ψ HH (q) ∼ e −I R (q)/ cos(S(q)/ ) .

(B.13)
This is a manifestly real wavefunction.
So far we had been careful not to use Lorentzian vs. Euclidean geometries, what is important is that the tunneling behaviour incorporates both classically forbidden modes, whilst the no-boundary has both classically allowed oscillatory modes. Moving in the field space it is also possible to cross Stokes-lines, so extra care needs to be taken in the choice and interpretation of the contour. In addition a complex saddle can have different interpretations: if we fix two points in superspace we can choose various contours to connect them. This is used in the main text in section 6.

C. Parabolic cylinder functions
These are the eigenfunctions of the inverted harmonic oscillator time independent Schroendinger equation. The specific differential equation is

C.1 Even -odd basis
A useful basis of solutions for real λ are the delta function normalised even/odd parabolic cylinder functions [33] which we will denote by ψ ± (ω, λ)
Another useful representation is in terms of Whittaker functions Φ(ω, λ) = 2 −iω/2 e πω/4+iφ 2 /2 e i5π/8 The complex solutions obey the following relations so they form an appropriate orthonormal set as well (up to the factor of −i).

C.3 Resolvent and density of states
In this subsection we define the usual resolvent and density of states corresponding to the inverted oscillator equation (C.1) that defines the HamiltonianĤ. The resolvent operator is defined asR The fixed energy amplitude between two position states is then The expression for the propagator is given by the Mehler formula for parabolic cylinder functions [33]: (C.10) which is computing the real-time T inverted H.O. propagator. This holds for −π < T < 0 or T = 0 with T = 0. To prove it one can use the general expression (7.694) in [53]. Notice that the same expression is also related to the Euclidean propagator for the normal oscillator upon double analytic continuation.
Another useful formula is (Sokhotski-Plemelj) The density of states is then given by (c is an infinitesimal number)

C.4 Dual resolvent and density of states
The expressions of the previous section were useful for performing computations using the usual interpretation of the theory as a string theory in a linear dilaton background. In our interpretation it will be instead useful to define a dual version of the resolvent that is related with the notion of time and energy for the boundary theory "living" on the macroscopic loop. This "duality" is in a sense exchanging the two indices of the Parabolic cylinder functions, since the new notion of energy is conjugate to the loop length and hence is directly identified with the coordinate label of the parabolic cylinder function. If we wish to explicitly compute the integrals below it turns out that instead of using the even/odd parabolic cylinder function basis, it is more convenient to use the complex basis defined by (C.5).
We first define the dual resolvent operator aŝ whereĤ dual is the position operatorλ. Once should also understand the role of the chemical potential µ from this point of view, though. What we find is that it provides an IR mass-gap for the eigenvalues, at the semi-classical level 27 . We can define again a fixed "energy" amplitude now as a transition amplitude between the previous oscillator energy eigenstates (C.14) with the new dual version of the Mehler formula being This is most easily computed in the complex basis (C.5) where it takes the form By expressing the complex wavefunctions in terms of the Whittaker function W µ,ν (z), one can then use (6.643) of [53] to express them in terms of Bessel functions. By changing the order of the integrals one finds the following expression for (C.16) (C.17) The first thing to notice is that in the limit z → 0 we recover δ(ω 1 − ω 2 ) due to the orthonormality of the complex wavefunctions. Unfortunately we are not aware of a more compact way of expressing the dual Mehler kernel (C.16).
The new dual density of states (E ≡ λ) is This is exactly equivalent to the formula (4.10) of the main text. We also observe the importance of the mass gap µ and the definition of the fermionic vacuum state |µ for obtaining a non-trivial result. Our final expression, is the definition of the density of states as an operator. This is obtained from the general time dependent fermionic bi-linear operator (3.7) upon Euclidean time averaging i.e.
This then means that we can compute any correlator of the dos operator upon sending q i → 0 in eqn. (D.4) of appendix D.

C.5 Density correlation functions
Using the previous formulae it is also easy to compute the density correlation functions or resolvent correlation functions. It is useful to use (C.3) and (C.4) to resolve the identity operator accordingly. For example the two point function of the dual resolvent is As another example the two point function of the dual density of states is given by Nevertheless in this case it is most convenient to use the general formula eqn.(D.4) sending q i → 0 and use the inverted oscillator propagator to get sinh(s 1 +s 2 ) sinh s 1 sinh s 2 The connected piece to the spectral form factor is given by the fourier transform of this expression that is found to match the expression in the main text 5.7.
In this appendix we review the computation of various correlation functions through the use of the fermionic field theory described in section 3. It was first described in detail in [33], which we closely follow almost verbatim.
Using this, the generic formula for the n-point density correlator then takes the form  with ξ = n i=1 s i and Q σ k = q σ(1) + ... + q σ(k) . From this expression it is possible to describe the general correlation function, expressed through the formula (q i , z i are the correlator parameters) ∂ Ô (q 1 , z 1 )...Ô(q n , z n ) ∂µ = i n+1 δ( 2 + 1 2 ωλ 2 so that in the end one analytically continues ω to get the result for the inverted oscillator. One can also analytically continue x to get the real time result. One then finds a compact expression for the derivative of the n-point loop correlator in terms of an n-fold integral  So far in this expression we have used both sides of the inverted oscillator potential, so that ξ ∈ (−∞, ∞). If we wish to describe the bosonic theory, we can focus on the one side of the potential ξ ∈ [0, ∞). In the integral (D.6), we can analytically continue z i = i i in the positive ξ region and z i = −i i in the negative ξ region to obtain convergent answers. This is the expression we use in the main text in the specific cases of a single or two loops.

E. Correlation functions for a compact boson
The analysis so far was performed in the case of a non-compact boson X(z,z). If we take this to be compact with a radius R, it is still possible to compute nonperturbatively all the correlation functions as before [35]. In order to do so one needs to use the formalism of free fermions at finite temperature. The thermal vaccum satisfies (notice that β = 2πR is not related with , the temperature of the holographic boundary theory, but is related with a "temperature in superspace") µ|b † s (ω)b s (ω)|µ R = δ ss 1 e β(µ−ω) + 1 , (E. 1) which means that in all the formulas of appendix D one just needs to replace the strict occupation θ(ω −µ) with the fermi distribution f (ω) = 1 e β(µ−ω) +1 . It also results in having discrete frequencies ω n = (n + 1 2 )/R, which certifies the anti-periodicity of the correlator around the compact Euclidean time direction with ξ = n i=1 s i and Q σ k = q σ(1) + ... + q σ(k) . This simply means that one can just replace δ( q i ) → Rδ q i , adding an extra factor ξ/2R sinh(ξ/2R) in the formulae of the previous appendix D. For example one finds the simple interesting relation between the infinite/finite radius observables The genus expansion is achieved using

F. Steepest descent
We will briefly review here the method of steepest descent that can be used to obtain asymptotic expansions of integrals as one parameter goes to infinity. Let us define I(t) = C h(s)e tρ(s) ds , (F.1) an integral supported in a contour C of the complex plane. We are interested in the t → ∞ asymptotics of this integral. Depending on the analyticity properties of the functions h(s), ρ(s) we can deform the contour to another contour C . If ρ(s) = φ(s) + iψ(s) with φ(s), ψ(s) real functions, it is useful to deform the contour such that we fix ρ(s) = const.. This is because the ψ = const. contours are parallel to the steepest contours in φ. Since the full asymptotic expansion of the integral is governed by the neighbourhood of s on which φ(s) acquires its maximum value, it is then useful to use precicely these steepest contours to approximate the integral. This idea can fail in case ρ (s * ) = 0. In such a case distinct steepest contours can intersect at the point s * and one needs to make a careful choice depending on the precise way one takes the t → ∞ limit on the complex t plane. This can lead to the Stokes phenomenon. This phenomenon is simply the fact that the analytic continuation of an asymptotic expansion of a function does not agree with the asymptotic expansion of the exact function's analytic continuation. A special case that will be of interest here is the case of a movable maximum. Such a case arises for example when at the saddle point φ (s) = 0 the function h(s) goes to zero exponentially fast. In such cases one is instructed to find the maxima of the total exponent log(h(s)) + t φ(s). These maxima are then dependent on t, s * = s * (t). In such a case one needs to rescale s appropriately, so that the maxima no-longer depend on t. Another special case is when at the maximum φ(s * ) blows up. This case can be treated in the same fashion with the movable maxima.

F.1 Stationary phase approximation
In the main text we encountered various integrals with highly oscillatory behaviour as we send t → ∞. Here for concreteness we collect some useful results of an asymptotic analysis of such integrals. This covers a subcase of integrals for which we can apply the more general steepest descent method. the integral whose asymptotic properties as t → ∞ we wish to analyse. Asumming that there is no stationary point g (s) = 0 in the region s ∈ [a, b], then we can use the general Riemann-Lebesgue lemma and its corollaries • If f (s) is integrable at the region of support and if g(s) is continuously differentiable and non-constant in any sub-interval then I(t) → 0 as t → ∞.
• In particular integration by parts gives the leading asymptotic behaviour of (F.4) provided that f (s)/g (s) is smooth for s ∈ [a, b], that is • If there is a stationary point of g(s) denoted by c, then if g (c) = g (c) = ... = g (p−1) (c) = 0 and g (p) (c) = 0 the leading asymptotic behaviour is given by where the sign of the phase is the same as the sign of g (p) (c).