Gas of baby universes in JT gravity and matrix models

It has been shown recently by Saad, Shenker and Stanford that the genus expansion of a certain matrix integral generates partition functions of Jackiw-Teitelboim (JT) quantum gravity on Riemann surfaces of arbitrary genus with any fixed number of boundaries. We use an extension of this integral for studying gas of baby universes or wormholes in JT gravity. To investigate the gas nonperturbatively we explore the generating functional of baby universes in the matrix model. The simple particular case when the matrix integral includes the exponential potential is discussed in some detail. We argue that there is a phase transition in the gas of baby universes


Introduction
It has been shown by Saad, Shenker and Stanford [1] that the genus expansion of a certain matrix integral generates the partition functions of Jackiw-Teitelboim (JT), [2,3], quantum gravity on Riemann surfaces of arbitrary genus with an arbitrary fixed number of boundaries. It is shown in [1] that an important part of JT quantum gravity is reduced to computation of the Weil-Petersson volumes of the moduli space of hyperbolic Riemann surfaces with various genus and number of boundaries for which Mirzakhani [4] established recursion relations. Eynard and Orantin [5,6] proved that Mirzakhani's relations are a special case of random matrix recursion relations with the spectral curve y = sin(2πz) 4π. This is a natural extension of results on topological gravity [7][8][9][10]. Relation of random matrices and gravity, including black hole description, has a long history, see [11,12] and refs therein.
The results of [1] provide a nonperturbative approach to JT quantum gravity on Riemann surfaces of various genus and perturbative description of boundaries. We use an extension of this result for nonperturbative studying of gas of baby universes in JT gravity. To investigate the boundaries nonperturbatively we explore the generating functional of boundaries in the matrix model and in JT gravity. One interprets the generating functional as the partition function of gas of baby universes in grand canonical ensemble in JT multiverse with the source function describing the distribution of boundaries being treated as the chemical potential. The interaction is presented by splitting and joining of baby universes 1 .
In this note we consider the generating functional for the gravitational correlation functions Z grav n (β 1 , ..., β n ; γ) Z grav (J; γ) =  where J(β) is a source function. An appropriate generating functional in matrix theory 1 One can compare this picture with string interactions and using this analogy closed strings describe baby universes without boundaries, meanwhile the baby universes with boundary correspond to open strings. An analogue of matrix theory is given by string field theory [13][14][15]. As has been noted in [1] there is an essential difference in coupling constant in SFT and JT.
has the form Here Z(β) = Tre −βM where M is a random N × N Hermitian matrix. This amounts to shifting the potential in the matrix model We define the generating functional for connected correlation functions take the double scaling limit introducing the parameter γ and obtain the relation between JT gravity and the matrix model in terms of the generating functionals: The "≃" symbol indicates the equality in the sense of formal series.
The paper is organized as follows. In Sect.2 the generating functional in matrix theory Z matrix (J) is discussed. Here J is the source function. In Sect.3 the generating functional of boundaries in JT gravity Z grav (J) is considered. In Sect.4 we investigate the double scaling limit in matrix models with a particular choice of the source J(β) which leads to the change of the potential V (x) → V (x) − Je ωx . In Sect.5 the matrix model with the exponential potential is investigated. In Sect.6 the matrix model with the spectral curve y = sin(2πz) 4π and the source is discussed and phase transition is observed. In Sect.7 the discussion of obtained results is presented.

Generating functional in matrix models
Generating functional. We consider ensemble of N × N Hermitian matrices [16][17][18][19] with potential V (M ). Let where λ i are eigenvalues of the matrix M . The n-point correlation function of Z(β) in the matrix model is given by Its generating functional can be presented as This amounts to shift the potential V (x) → V (x) −J(x).
One expands G(J) = log Z matrix (J) to get the connected correlation functions Particular case. We consider a special case In this case the consideration of Z matrix (J) is equivalent to dealing with the matrix model with a deformed potential In this case the singular integral equation defining the eigenvalues distribution has the form 2 10) and the spectral density is Double scaling limit. 3 All the correlation functions Z matrix n (β 1 , ..., β n ) in principle could be derived if the potential V (x) or the spectral density/spectral curve ρ(µ) is known. To get a connection of the matrix model with JT gravity one has to go to the double scaling limit, see [1]. Consider a matrix model with a non-normalized spectral density where S 0 is a constant. Now, shifting E → E − a and sending a → ∞ we get the spectral density of the double-scaled matrix model The correlation functions Z matrix, d.s. g,n (β 1 , ..., β n ) and the constant S 0 will be used in the next section to describe the connection of the matrix model with JT gravity. The double scaling limit will be discussed also in Sect.4.

Resolvents.
Similarly one has generating functional for correlation functions of resolvents where f (z) is a test function and Expanding log R (matrix) (f ) in the series on f one gets connected correlation functions < R(z 1 )...R(z n ) > c . We perform the double scaling limit which admits an expansion One defines the correlation functions which satisfy the loop equations [28][29][30] and will be used in the next section, and the generating functional

Generating functional in JT gravity
The Euclidean action of JT gravity [2,3,31] has the form Here g µν is a metric on a two dimensional manifold M, φ is a scalar field (dilaton) and the constant S 0 was mentioned in the previous section. The path integral for Riemann surface of genus g ≥ 2 with n boundaries with lengths β 1 , ..., β n reads Z grav g,n (β 1 , ..., β n ) = e −S 0 χ where χ is the Euler characteristic χ = 2 − 2g − n andÎ JT is the JT action with the first S 0 term left out.
The following relation between the matrix model and JT gravity holds [1]: (e −S 0 ) 2g+n−2 Z grav g,n (β 1 , ..., β n ) (3.5) It was found in [1] that the partition function has the form (g ≥ 2) is the Weil-Petersson volume of the moduli space of a genus g Riemann surface with n geodesic boundaries of lengths b 1 , ..., b n and From this we get and the generating functional Finally, one obtains the relation Similarly, the correlation functions W matrix g,n are related with volumes of the moduli spaces as The generating functional is (3.14) Baby universes. In cosmology [33][34][35][38][39][40], one usually deals with baby universes that branch off from, or join onto, the parent(s) Universe(s). In matrix theories one parent is a connected Riemann surface with arbitrary number of handles and at least one boundary. We assume that the lengths of boundaries of baby universities are small as compare with the boundary length of the parent, see Fig.1. Baby universes are attached to the parent by necks that have restricted lengths of geodesics at which the neck is attached to the parent, We assume that the lengths of the boundaries of baby universes are small compared to the length of the boundary of the parent, see Fig.  1. Baby universes are attached to the parent with the help of thin necks. Thickness of the neck is defined as the geodesic length of the loop located at the thinness point of the neck, and this length is assumed to be essentially smaller than the length of theboundary of the parent, see Fig.1.b and Fig.1.c. There are also restrictions on the area of the surface of baby universes, see [41,42] for more precise definitions. Cosmological baby universes in the parent-baby universe approximation interact only via coupling to the parent universes, that themselves interact via wormholes. In matrix models the baby universes always interact via their parents too and parents interact via wormholes, Fig.2. One can expect that at large number of baby universes interaction between different parts of the system increases and this leads to phase transition (an analog of the the nucleation of a baby universe in [35]). We interpret the matrix partition function Z, defined by equation (2.16) as a partition function of the gas of baby universes. Figure 2. Two parents connected by the wormhole 4 Double scaling limit 4.1 Double scaling limit for the GUE There are various notions of the double scaling limit in matrix theory [22][23][24][25][26][27][43][44][45][46][47][48][49]. A special double scaling limit was considered in [1] at the level of spectral density. Here we discuss it at the level of the potentials. We will see that the linear term in the potential plays a special role.
Let us start with the Wigner distribution for the Gaussian Unitary Ensemble (GUE) [16][17][18][19]. The ordinary Wigner distribution 2 √ a 2 − λ 2 πa 2 is supported on the interval [−a, a] and is obtained from a matrix model with the potential V (x) = m 2 x 2 2, where m 2 = 4 a 2 . We want to make a shift and get a distribution on the interval [0, Λ], Λ > 0. To this end we consider the gaussian model with an external source 4 and we make parameters m and j depending on Λ to put the measure support on [0, Λ]. The singular integral equation defining the density takes the form Here ffl means the Cauchy principal value of the integral. The solution of (4.2) is given by the shifted Wigner distribution Note that the constant j from the linear term in the potential does not enter into the expression for the spectral density. This is valid for any potential. The constant appears only throughout the normalization and consistency conditions that in this case read We get that the eigenvalue density (4.3) is supported by the potential in the sense that where η = 3 2. So, we get an expected result in two steps. First we send N → ∞ and then Λ → ∞.
One can use also another procedure. Set Λ = t N and in this case one has x (4.10)

Double scaling limit for the cubic interaction
Let us consider the matrix model with the potential and we will make parameters m, j, g depending on Λ to put the support of the spectral density on [0, Λ]. Again first we take the limit N → ∞ and then the limit Λ → ∞.
In the large N limit one gets the singular integral equation that defines the spectral density We write solution in the form given by (4.3), and we get (4.14) The normalization condition is achieved by a suitable choice of m 2 , (4.15) and the consistency condition by the suitable choice of j The solution of these equations is evidently which should be substituted into the potential. Now we take the large Λ limit and get More reach picture appears for the quartic interaction with negative mass square. In this case for suitable choice of parameters there is a double cut solution [51,52], universality of the double scaling limit in this case has been proved in [45], see also [53] for more general multi-cut solutions.
5 Deformation by an exponential potential 5

.1 Exponential potential
Here we consider a deformation of the Wigner distribution by the insertion of the exponential potential 5 There are 4 different choices of signs of J and ω, see Fig.3. We see that only J < 0 may produce some non trivial effects due to an appearance of potential instability, Fig.3.c and Fig.3.d. The choice of sign of ω is irrelevant.   We see that the minimum of the potential disappears when m decreases, Fig.4.a, or J increases for negative J, Fig.4.b.  The shift of λ → λ + δλ in the potential (5.1) produces the linear term in the potential and multiplies the current J on the constant J → Je δλ . We parametrize our potential as and fix the parameters in (5.2) in an agreement with the measure localization on the segment [0, Λ]. To this purpose we first find the non-normalized measure as a sum of non-normalized ones corresponding to the shifted Wigner and exponential potential distributions ρ W (λ, Λ) and ρ e wλ (λ, Λ, ω).
The forms of non-normalized measures ρ e wλ (λ) for positive and negative w presented are presented in Fig.5. We can compare the contribution to the non-normalized density from the Wigner semi-circle for mass equal to 1 and the exponential potential taken with arbitrary current J. We see that this sum always defines the positive density for J > 0 and becomes negative for J < J cr (Λ, ω) < 0. Then for J > J cr (Λ, ω) we find relation between m 2 and Λ, ω and J from normalization condition.
In Fig.6 the appearance of the phase transition at negative J is presented. In Fig.  7.a relations between m 2 and J for fixed Λ and ω are shown for different values of Λ and in Fig.7.b. We can present the eigenvalues distribution corresponding to the potential (5.2) as a sum 6 and fix constant C and m 2 from the consistency condition and normalization, respectively, The forms of non-normalized measures ρ e wλ (λ) for positive and negative w presented are presented in Fig.5. We see that these two segment distributions are nonsymmetric under the centre of the segment. The distribution of eigenvalues for the case (e) (f) Figure 6. The plot of non-normalized density for the quadratic potential deformed by the exponential potential for different values of the regularization parameter Λ and ω = 1.
of negative ω is pressed to the left boundary of the segment, and the for the case of positive ω it is pressed to the right one. Applying this deformation with a positive J to the GUE we "activate" the left or right part of of eigenvalues. Applying the same with negative J we can destroy the constructed solution. We compare the contribution to the non-normalized density from the Wigner semi-circle for mass equal to 1 and the exponential potential taken with arbitrary current J. We see the this sum always define the positive density for J > 0 and becomes negative for J < J cr (Λ, ω) < 0.
Then J > J cr (Λ, ω) we find relation between m 2 and Λ, ω, and J from normalization condition.
In Fig.6 is presented the appearance of the phase transition at negative J for different Λ. We see that for chosen parameters, m 2 = 1, ω = 1 the critical J cr decreases with increasing Λ. To find the real mass that supports the normalized solution, we find n(λ) from the normalization condition, so and assume that the mass is given my In Fig.7.a the dependence of m 2 on J for fixed Λ and ω is shown for different values of Λ. Here J > J cr (Λ). We see that mass (in our parametrization of the potential of the model) decreases with increasing J. The slow of the mass is more fast for larger Λ. In Fig.7.b the dependence of the critical current J cr on mass is shown. We see that J cr goes to zero when m 2 goes to zero, that corresponds to increasing Λ.

Fine tuning
It is obvious that fixing from the beginning the location of the eigenvalue one immediately gives restrictions on parameters of the potential of the matrix model. If we want to shift the location of the eigenvalues, S [a,b] → S [0,b−a] we have to make a shift in the potential, V (x) → V (x + a). For the quadratic potential this shift produces the linear term ∆V (x) = jx and j can be determined from the location of the left point of the cut. For higher polynomial interaction the shift produces the linear term in the LHS of singular equation, as well as change of coupling constants. The shift in the exponential potential produces just a multiplication on positive constant.
We have also seen that if we want to deform a given distribution by a new potential, that has the same locations of eigenvalues, we can just take the sum of the given two distributions and multiply all coupling constants of two initial model on the same parameter to fix the normalization condition for the distributions that is the some of given two distributions. As to the consistency condition it follows from consistency conditions of individual distributions. More precisely, if we know that Note that in both integrals the segment is the same. Taking we can claim that ρ T solves equations The consistency condition is automatically satisfied.
One can put a coupling constant in front of the second potential, say J. In our previous example this were the coupling constant g in the perturbation of Gaussian model by qubic term, or J in the case of the exponential potential. For positive J we keep the positivity condition for the sum of two distribution, meanwhile we can lost it for the case of big negative current. This loss of positivity leads to the destruction of the large N expansion of the model and can be interpreted as a phase transition.
As to double scaling limit, we can consider it in three steps. First we bring the support of the eigenvalue distribution on the interval [0, Λ] by fine tuning the linear term in the potential. Then one goes to the limit N → ∞ and after that Λ → ∞.
6 Matrix model for JT gravity 6

.1 Potentials for non-normalized density distribution
For fixed Λ we consider the eigenvalue distribution normalized to 1 This form of the eigenvalue distribution in the SYK model has been obtained in [55,56] and is nothing but the Bethe formula for the nuclear level density [57]. For large Λ one has n(Λ) ≈ 8π 3 e −2π √ Λ √ Λ. The distribution normalized to N has the form compare with (2.13). It is evident that e S 0 = N n(Λ) To recover the potential that supports the distribution ρ norm,1 (E), we write This gives Here Ei is the exponential integral that for real non zero values of x is defined as Ei has an expansion and faster converging Ramanujan's series has the form Ei can be bounded by elementary functions as follows The potential up to a constant is and it is presented in Fig.8.

Effective energy
The effective energy E ef f , [18], is evaluated on the normalized density and it is given by the following formula In Fig.9 the effective action E ef f as function of Λ is shown. We see that it can be approximated by a log 10.60 + 0.959 log x (the next approximation is given by the double-log, 10.76 + 1.050 log x − 0.393 log log x).

Phase transition
To study the phase transition we consider the deformation of ρ sinh It is interesting to compare this density with the density (6.14) In Fig.10 we plot the density (6.12) and (6.13) for ω = −1 and different J and Λ. We see that for negative J > J cr (Λ) domains on the segment 0, Λ, where ρ nn and ρ nn,s become negative, appear. We interpret this as a destruction of the solution, or in other words, as an appearance of a phase transition at J cr . It is interesting that the critical value J cr are the same for both ρ nn and ρ nn,s . In Fig.11 the locations of critical points depending on Λ are shown.
In Fig.12 the dependence of n(Λ) on J for fixed Λ and ω for different values of Λ is shown. For (a),(b) and (c) these plots are for Λ = 5, 7.5, 10. Here J > J cr (Λ). We see normalization factor increases with increasing J and this dependence is linear. This dependence differs from the dependence for the Gaussian perturbed ensemble (see Fig.7.a), where it is given by a decreasing nonlinear function. The plot in Fig.12.b) shows the dependence of critical current J cr on n(Λ).  Relations between n(Λ), normalization factor, and J for fixed Λ for ρ(λ) given by (6.12). The lines start at the points presented at Fig.11.a.

Discussion and conclusion
The generating functional for the correlation functions of the boundaries of the Riemann serfaces is considered in JT gravity and in matrix theory. The matrix integral [1] provides a nonperturbative completion of the genus expansion in JT gravity on Riemann surfaces with fixed number of boundaries for all genus. The generating functional Z(J) considered in this paper, gives completion also for infinitely many number of boundaries.
By using this formulation with several matrix models, including double scaling limit for the Gaussian model, relevant to topological gravity, cubic model, and JT gravity have been investigated. In all these cases (here we have presented only results for topological and JT gravities), in the corresponding gas of baby universes the phase transition is observed. By analogy with states of matter one could expect that this phase change is condensation from the gaseous state of JT multiverse to liquid state.
To study this phase transition, we consider generating functional that we associate with baby universes. The generating functional for a matrix model with potential V (x) is obtained from the partition function of the model just by shift V (x) → V (x) −J(x) whereJ(x) is the Laplace transform of J(β). Baby universes should look as point-like objects with very small length of the boundary. Therefore the distribution J(β) should have a peak somewhere near zero. In our model with J(β) = Jδ(β − ω) baby universes correspond to small ω. The shifted potential in this case is V (x)−Je −ωx . Note that one could expect < Z(β) n >≈< Z(β) > n for large n similarly to the infinite replica number limit considered in [67].
We have obtained that in all cases there are critical negative J cr , such that for J < J cr the corresponding solution is destroyed. The mass at the critical points decreases with increasing of Λ. In this paper we have considered the deformation of the density function by the exponential potential, but is also possible to study all momenta deformations, as well as det(E − M ), compare with [58].
Note that the deformation of potential by the linear term or nonlinear ones as a tool for investigation of phase transitions and spontaneous symmetry breaking in SYK-like models was used in [64,65,72]. In particular, a nonlocal interacting of two-replica by a nonlocal term proportional to an external current permits to reveal nonperturbative effects in the SYK model. This nonlocality in some sense is analogous to an external nonlocal source J(x) (cosmological daemon) in cosmology [69], there it specified boundary conditions.
There are numerous investigations of wormholes and baby universes in cosmology and particle physics, including the Giddings-Strominger wormhole solution [35] and Coleman's approach to the cosmological constant problem [38]. There are many open questions in theory of wormholes and baby universes. In particular, it would be interesting to see whether in JT gravity or its generalizations there is a mechanism of suppression [70] the probability of creation of giant wormholes and big baby universes. By using the wormhole/baby universe approach it was found that the probability for the universe to undergo a spontaneous compactification down to a four-dimensional spacetime is greater than to remain in the original homogeneous multidimensional state [39]. It is interesting to find an analog of this in the context of JT gravity.
Using the recent result of [79] according which JT gravity in dS 2 is an analytic continuation of JT gravity in Euclidean AdS 2 it would be interesting to understand the meaning of the phase transition considered here in the dS case.
It would be also interesting to to compare the Hartle-Hawking construction of string baby universes related with AdS 2 × S 2 geometry, using free fermionic formulation [80] with the fermionic interpretation of determinant in matrix models, and also find a slot for baby universes in an quasi-classical [81,82] or exact quantization of JT proposed in [83].