Gas of Baby Universes in JT Gravity and Matrix Models

Abstract

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.

1. 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\left(2\pi z\right)/4\pi$. 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. 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 (SFT) [13,14,15]. As has been noted in [1] there is an essential difference in coupling constant in SFT and JT.

Let $Z_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n)$ be the JT gravity path integral for Riemann surface of genus $g\ge 2$ with n boundaries with lengths ${\beta }_1,\cdots ,{\beta }_n$. Consider a generating function for these functions

$${\mathcal{Z}}_n^{grav}({\beta }_1,\cdots ,{\beta }_n;\gamma )\simeq \sum _{g=0}^{\infty }{\gamma }^{2g+n-2}{\mathcal{Z}}_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n)$$

(1)

where $\gamma$ is a constant which in notations of [1] is $\gamma =e^{-S_0}$.

The following remarkable relation between correlation functions in matrix model and JT gravity holds [1]:

$${\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n;\gamma )\simeq Z_n^{grav}({\beta }_1,\cdots ,{\beta }_n;\gamma )$$

(2)

Here ${\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n;\gamma )$ is the double scaling (d.s.) limit of the correlation function in a matrix model with the spectral curve mentioned above. This form of the curve was obtained in [1] by computing the JT path integral for the disc.

In this note we consider the generating functional for the gravitational correlation functions $Z_n^{grav}({\beta }_1,\cdots ,{\beta }_n;\gamma )$

$${\mathcal{Z}}^{grav}(J;\gamma )=\sum _{n=0}^{\infty }\frac{1}{n!}{\int }_0^{\infty }d{\beta }_1...{\int }_0^{\infty }d{\beta }_n{Z}_n^{grav}({\beta }_1,\cdots ,{\beta }_n;\gamma )J\left({\beta }_1\right)...J\left({\beta }_n\right),$$

(3)

where $J\left(\beta \right)$ is a source function. An appropriate generating functional in matrix theory has the form

$$\begin{array}{c}\hfill {\mathfrak{Z}}^{matrix}\left(J\right)=<{\mathrm{e}}^{N\int Z\left(\beta \right)J\left(\beta \right)d\beta }>\\ \hfill =\sum _{n=0}^{\infty }\frac{N^n}{n!}\int d{\beta }_1...\int d{\beta }_n{\mathcal{Z}}_n^{matrix}({\beta }_1,\cdots ,{\beta }_n)J\left({\beta }_1\right)...J\left({\beta }_n\right).\end{array}$$

(4)

Here $Z\left(\beta \right)=\mathrm{Tr}{e}^{-\beta M}$ and M is a random $N\times N$ Hermitian matrix. This amounts to a shift of the potential in the matrix model $V\left(x\right)\to V\left(x\right)-\tilde{J}\left(x\right)$, where $\tilde{J}\left(x\right)$ is the Laplace transform of $J\left(\beta \right)$, see Section 2.

We define the generating functional for connected correlation functions

$${\mathfrak{G}}^{matrix}\left(J\right)=-\frac{1}{N^2}log{\mathfrak{Z}}^{matrix}\left(J\right),$$

(5)

take the double scaling limit introducing the parameter $\gamma$ and obtain the relation between JT gravity and the matrix model in terms of the generating functionals:

$$d.s.lim{\mathfrak{G}}^{matrix}\left(J\right)\simeq {\mathfrak{Z}}^{\mathrm{grav}}(J;\gamma ).$$

(6)

The “≃” symbol indicates the equality in the sense of formal series.

The paper is organized as follows. In Section 2 the generating functional in matrix theory ${\mathfrak{Z}}^{matrix}\left(J\right)$ is discussed. In Section 3 the generating functional of boundaries in JT gravity ${\mathfrak{Z}}^{grav}\left(J\right)$ is considered. In Section 4 we investigate the double scaling limit in matrix models with a particular choice of the source $J\left(\beta \right)$ which leads to the change of the potential $V\left(x\right)\to V\left(x\right)-J{e}^{\omega x}$. In Section 5 the matrix model with the exponential potential is investigated. In Section 6 the matrix model with the spectral curve $y=sin\left(2\pi z\right)/4\pi$ and the source is discussed and phase transition is observed. In Section 7 the discussion of obtained results is presented.

2. Generating Functional in Matrix Models

Generating functional. We consider ensemble of $N\times N$ Hermitian matrices [16,17,18,19] with potential $V\left(M\right)$. Let

$$Z\left(\beta \right)=\mathrm{Tr}{e}^{-\beta M}=\sum _{i=1}^{N}exp(-\beta {\lambda }_i),\beta >0,$$

(7)

where ${\lambda }_i$ are eigenvalues of the matrix M. The n-point correlation function of $Z\left(\beta \right)$ in the matrix model is given by

$$\begin{array}{c}\hfill {\mathcal{Z}}_n^{matrix}({\beta }_1,\cdots ,{\beta }_n)\equiv <Z\left({\beta }_1\right)...Z\left({\beta }_n\right)>=\frac{\int Z\left({\beta }_1\right)...Z\left({\beta }_n\right){\mathrm{e}}^{-N\mathrm{Tr}V\left(M\right)}dM}{\int {\mathrm{e}}^{-N\mathrm{Tr}V\left(M\right)}dM}.\end{array}$$

(8)

Its generating functional can be presented as [20]

$$\begin{array}{c}\hfill {\mathfrak{Z}}^{matrix}\left(J\right)=<{\mathrm{e}}^{N\int Z\left(\beta \right)J\left(\beta \right)d\beta }>=\sum _{n=0}^{\infty }\frac{N^n}{n!}\int {\mathcal{Z}}_n^{matrix}({\beta }_1,\cdots ,{\beta }_n)\prod _{i=1}^{n}J\left({\beta }_i\right)d{\beta }_i,\end{array}$$

(9)

or

$${\mathfrak{Z}}^{matrix}\left(J\right)=\int {\mathrm{e}}^{N{\sum }_{j=1}^N\tilde{J}\left({\lambda }_j\right)}d{\mu }_N({\lambda }_1,\cdots ,{\lambda }_N),$$

(10)

where

$$\begin{array}{ccc}\hfill d{\mu }_N({\lambda }_1,\cdots ,{\lambda }_N)& =& \frac{1}{Z_N}\prod _{j>k}{({\lambda }_j-{\lambda }_k)}^2\prod _{j=1}^N{\mathrm{e}}^{-NV\left({\lambda }_j\right)},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill \tilde{J}\left(\lambda \right)& =& \int d\beta J\left(\beta \right){\mathrm{e}}^{-\beta \lambda }d\beta .\hfill \end{array}$$

(12)

This amounts to shift the potential $V\left(x\right)\to V\left(x\right)-\tilde{J}\left(x\right)$.

One expands $\mathfrak{G}\left(J\right)=log{\mathfrak{Z}}^{matrix}\left(J\right)$ to get the connected correlation functions

$$\begin{array}{c}\hfill {\mathfrak{G}}^{matrix}\left(J\right)=\sum _{n=0}^{\infty }\frac{1}{n!}\int {\mathcal{Z}}_{n,conn}^{matrix}({\beta }_1,\cdots ,{\beta }_n)\prod _{i=1}^{n}J\left({\beta }_i\right)d{\beta }_i.\end{array}$$

(13)

Particular case. We consider a special case

$$J\left(\beta \right)=-J\delta (\beta +\omega ).$$

(14)

In this case the consideration of ${\mathfrak{Z}}^{matrix}\left(J\right)$ is equivalent to dealing with the matrix model with a deformed potential

$$U\left(x\right)=V\left(x\right)+J{\mathrm{e}}^{\omega x}.$$

(15)

In this case the singular integral equation (See [21,22] for the theory of singular integral equations.) defining the eigenvalues distribution has the form

$$U^{\prime }\left(\lambda \right)=2{⨏}_{S}d\mu \frac{\rho \left(\mu \right)}{\lambda -\mu }$$

(16)

and the spectral density is

$$\rho \left(\mu \right)=\underset{N\to \infty }{lim}\frac{1}{N}<\sum _{i=1}^N\delta (\mu -{\lambda }_i)exp\{-NJ\sum _{j=1}^N{\mathrm{e}}^{\omega {\lambda }_j}\}>.$$

(17)

Double scaling limit. This limit in matrix models has been introduced in [23], see [24,25,26] for review and refs therein. All the correlation functions ${\mathcal{Z}}_n^{matrix}({\beta }_1,\cdots ,{\beta }_n)$ in principle could be derived if the potential $V\left(x\right)$ or the spectral density/spectral curve $\rho \left(\mu \right)$ 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

$${\rho }^{\mathrm{nnorm}}\left(E\right)=\frac{{\mathrm{e}}^{S_0}}{{\left(2\pi \right)}^2}sinh\left(2\pi \sqrt{\frac{a^2-E^2}{2a}}\right),-a<E<a,$$

(18)

where $S_0$ is a constant. Now, shifting $E\to E-a$ and sending $a\to \infty$ we get the spectral density of the double-scaled matrix model

$${\rho }_{d.s.}^{\mathrm{nnorm}}\left(E\right)=\frac{{\mathrm{e}}^{S_0}}{{\left(2\pi \right)}^2}sinh\left(2\pi \sqrt{E}\right),E>0,$$

(19)

and the correlation functions in the double scaled limit

$$d.s.lim{\mathcal{Z}}_{n,conn}^{matrix}({\beta }_1,\cdots ,{\beta }_n)={\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n).$$

(20)

The limiting correlation functions ${\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n)$ have an expansion of the form

$${\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n)\simeq \sum _{g=0}^{\infty }{\left(e^{-S_0}\right)}^{2g+n-2}{\mathcal{Z}}_{g,n}^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n).$$

(21)

The double scaling limit of the generating functional is

$$\begin{array}{c}\hfill d.s.limlog{\mathfrak{Z}}^{matrix}\left(J\right)=\sum _{n=0}^{\infty }\frac{1}{n!}\int {\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n)\prod _{i=1}^{n}J\left({\beta }_i\right)d{\beta }_i.\end{array}$$

(22)

The correlation functions ${\mathcal{Z}}_{g,n}^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_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 Section 4.

Resolvents. Similarly one has generating functional for correlation functions of resolvents

$${\mathfrak{R}}^{\left(matrix\right)}\left(f\right)=<{\mathrm{e}}^{N\int R\left(z\right)f\left(z\right)dz}>,$$

(23)

where $f\left(z\right)$ is a test function and

$$R\left(z\right)=\mathrm{Tr}{(z-M)}^{-1}.$$

(24)

Expanding $log{\mathfrak{R}}^{\left(matrix\right)}\left(f\right)$ in the series on f one gets connected correlation functions $<R\left(z_1\right)\cdots R\left(z_n\right){>}_c$. We perform the double scaling limit

$$d.s.lim<R\left(z_1\right)...R\left(z_n\right){>}_c=<R\left(z_1\right)...R\left(z_n\right){>}_c^{d.s.},$$

(25)

which admits an expansion

$$<R\left(z_1\right)...R\left(z_n\right){>}_c^{d.s.}\simeq \sum _{g=0}^{\infty }\frac{1}{{\left({\mathrm{e}}^{S_0}\right)}^{2g+n-2}}{R}_{g,n}(z_1,\cdots ,z_n).$$

(26)

One defines the correlation functions

$$W_{g,n}^{matrix}(z_1,\cdots ,z_n)={(-1)}^n{2}^n{z}_1...z_n{R}_{g,n}(-z_1^2,\cdots ,-z_n^2),$$

(27)

which satisfy the loop equations [20,27,28] and will be used in the next section, and the generating functional

$${\mathcal{W}}_g\left(f\right)=\sum _{n=0}^{\infty }\frac{1}{n!}\int d{z}_1...\int d{z}_n{W}_{g,n}^{matrix}(z_1,\cdots ,z_n)f\left(z_1\right)...f\left(z_n\right).$$

(28)

3. Generating Functional in JT Gravity

The Euclidean action of JT gravity [2,3,29] has the form

$$I_{JT}=-\frac{S_0}{2\pi }\left[\frac{1}{2}{\int }_{\mathcal{M}}\sqrt{g}R+{\int }_{\partial \mathcal{M}}\sqrt{h}K\right]-\left[\frac{1}{2}{\int }_{\mathcal{M}}\sqrt{g}\varphi (R+2)+{\int }_{\partial \mathcal{M}}\sqrt{h}\varphi (K-1)\right].$$

(29)

Here $g_{\mu \nu }$ is a metric on a two dimensional manifold $\mathcal{M}$, $\varphi$ 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\ge 2$ with n boundaries with lengths ${\beta }_1,\cdots ,{\beta }_n$ reads

$$Z_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n)=e^{-S_0\chi }\int \frac{\mathcal{D}{g}_{\mu \nu }\mathcal{D}\varphi }{\mathrm{Vol}\left(\mathrm{diff}\right)}{e}^{-{\widehat{I}}_{JT}[g_{\mu \nu },\varphi ]},$$

(30)

where $\chi$ is the Euler characteristic $\chi =2-2g-n$ and ${\widehat{I}}_{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]:

$${\mathcal{Z}}_{g,n}^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n)=Z_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n),$$

(31)

or

$${\mathcal{Z}}_n^{matrix,d.s.}({\beta }_1,\cdots ,{\beta }_n)\simeq Z_n^{grav}({\beta }_1,\cdots ,{\beta }_n),$$

(32)

where

$${\mathcal{Z}}_n^{grav}({\beta }_1,\cdots ,{\beta }_n)\simeq \sum _{g=0}^{\infty }{\left(e^{-S_0}\right)}^{2g+n-2}{\mathcal{Z}}_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n).$$

(33)

It was found in [1] that the partition function has the form $(g\ge 2)$

$$Z_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n)={\int }_0^{\infty }\cdots {\int }_0^{\infty }{V}_{g,n}(b_1,\cdots ,b_n)\prod _{i=1}^n{b}_i{Z}_{\mathrm{Sch}}^{\mathrm{trumpet}}({\beta }_i,b_i)d{b}_i,$$

(34)

where $V_{g,n}(b_1,\cdots ,b_n)$ is the Weil-Petersson volume of the moduli space of a genus g Riemann surface with n geodesic boundaries of lengths $b_1,\cdots ,b_n$ and

$$\begin{array}{c}\hfill Z_{\mathrm{Sch}}^{\mathrm{trumpet}}(\beta ,b)=\frac{1}{2{\pi }^{1/2}{\beta }^{1/2}}{e}^{-\frac{b^2}{4\beta }}.\end{array}$$

(35)

From this we get

$$\begin{array}{c}\hfill {\mathcal{Z}}_g^{grav}\left(J\right)=\sum _{n=0}^{\infty }\frac{1}{n!}{\left(e^{-S_0}\right)}^{2g+n-2}{\int }_0^{\infty }d{\beta }_1...{\int }_0^{\infty }d{\beta }_n{Z}_{g,n}^{grav}({\beta }_1,\cdots ,{\beta }_n)J\left({\beta }_1\right)...J\left({\beta }_n\right)\\ \hfill =\sum _{n=0}^{\infty }\frac{1}{n!}{\left(e^{-S_0}\right)}^{2g+n-2}{\int }_0^{\infty }{b}_{1}db...{\int }_0^{\infty }{b}_{n}d{b}_n{V}_{g,n}(b_1,\cdots ,b_n)\widehat{J}\left(b_1\right)\cdots \widehat{J}\left(b_n\right),\end{array}$$

(36)

where

$$\widehat{J}\left(b\right)={\int }_0^{\infty }d\beta \frac{1}{2{\pi }^{1/2}{\beta }^{1/2}}{e}^{-\frac{b^2}{4\beta }}J\left(\beta \right),$$

(37)

and the generating functional

$${\mathfrak{Z}}^{grav}\left(J\right)\simeq \sum _{g=0}^{\infty }{\mathcal{Z}}_g^{grav}\left(J\right).$$

(38)

Finally, one obtains the relation

$$d.s.limlog{\mathfrak{Z}}^{matrix}\left(J\right)\simeq {\mathfrak{Z}}^{\mathrm{grav}}\left(J\right).$$

(39)

Similarly, the correlation functions $W_{g,n}^{matrix}$ are related with volumes of the moduli spaces as

$$W_{g,n}^{matrix}(z_1,\cdots ,z_n)={\int }_0^{\infty }\cdots {\int }_0^{\infty }{V}_{g,n}(b_1,\cdots ,b_n)\prod b_i{e}^{-b_i{z}_i}d{b}_i.$$

(40)

The generating functional is

$$\begin{array}{cc}& {\mathcal{Z}}_{R,g}\left(f\right)=\sum _{n=0}^{\infty }\frac{1}{n!}\int \cdots \int d{z}_n{W}_{g,n}^{matrix}(z_1,\cdots ,z_n)\prod _{i=1}^{n}f\left(z_i\right)d{z}_i\hfill \\ =& \sum _{n=0}^{\infty }\frac{1}{n!}{\int }_0^{\infty }\cdots {\int }_0^{\infty }{V}_{g,n}(b_1,\cdots ,b_n)\prod _{i=1}^n{b}_{i}d{b}_i\tilde{f}\left(b_i\right),\tilde{f}\left(b\right)=\int dz{e}^{-bz}f\left(z\right).\hfill \end{array}$$

Baby universes. In cosmology [30,31,32,33,34,35,36,37,38], 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 Figure 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 Figure 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 Figure 1b,c. There are also restrictions on the area of the surface of baby universes, see [39,40] 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, Figure 1b,c). 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 [32]). We interpret the matrix partition function $\mathfrak{Z}$, defined by equation (22) as a partition function of the gas of baby universes.

4. Double Scaling Limit for the GUE

There are various notions of the double scaling limit in matrix theory [23,24,25,26,41,42,43,44,45,46,47,48,49,50,51]. 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\sqrt{a^2-{\lambda }^2}/\pi a^2$ is supported on the interval $[-a,a]$ and is obtained from a matrix model with the potential $V\left(x\right)=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,\mathsf{\Lambda }]$, $\mathsf{\Lambda }>0.$ To this end we consider the gaussian model with an external source. Note that the Gaussian matrix model with an arbitrary matrix source has been studied in [52]. Here we consider the case corresponding in notations [52] to $A=jI$, I is the $N\times N$ unit matrix.

$$V\left(x\right)=m^2\left(\frac{1}{2}{x}^2+jx\right)$$

(41)

and we make parameters m and j depending on $\mathsf{\Lambda }$ to put the measure support on $[0,\mathsf{\Lambda }]$. The singular integral equation defining the density takes the form

$$m^2(\lambda +j)=2{⨏}_{[0,\mathsf{\Lambda }]}d{\lambda }^{\prime }\frac{{\rho }_W\left({\lambda }^{\prime }\right)}{\lambda -{\lambda }^{\prime }}.$$

(42)

Here $⨏$ means the Cauchy principal value of the integral. The solution of (42) is given by the shifted Wigner distribution

$$\begin{array}{ccc}\hfill {\rho }_W\left(\lambda \right)& =& \frac{8}{\pi {\mathsf{\Lambda }}^2}{S}_{\mathsf{\Lambda }}\left(\lambda \right),S_{\mathsf{\Lambda }}\left(\lambda \right)=\sqrt{\lambda (\mathsf{\Lambda }-\lambda )}.\hfill \end{array}$$

(43)

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: $m^2{\mathsf{\Lambda }}^2=16,\mathsf{\Lambda }=-2j.$ We get that the eigenvalue density (43) is supported by the potential $V_{\mathsf{\Lambda }}\left(x\right)=(2x^2-4\mathsf{\Lambda }x)/{\mathsf{\Lambda }}^2$ in the sense that

$${\rho }_W\left(\lambda \right)=\underset{N\to \infty }{lim}{\rho }_N^{V_{\mathsf{\Lambda }}}\left(\lambda \right),\text{where}{\rho }_N^{V_{\mathsf{\Lambda }}}\left(x\right)=\frac{1}{N}<\mathrm{Tr}\delta (x-M){>}_{{}_{V_{\mathsf{\Lambda }}}}.$$

(44)

$<f{>}_{{}_V}$ in (44) means averaging with potential $<f{>}_{{}_V}=\int f{e}^{-N\mathrm{Tr}V}dM/\int e^{-N\mathrm{Tr}V}dM.$

Now we take the limit $\mathsf{\Lambda }\to \infty$ and write

$$\underset{\mathsf{\Lambda }\to \infty }{lim}{\mathsf{\Lambda }}^{\eta }{\rho }_W\left(\lambda \right)=\frac{8}{\pi }\sqrt{\lambda },\lambda \in [0,\mathsf{\Lambda }],$$

(45)

where $\eta =3/2$. So, we get an expected result in two steps. First we send $N\to \infty$ and then $\mathsf{\Lambda }\to \infty$.

One can use also another procedure. Set $\mathsf{\Lambda }=tN$ and in this case one has

$$\underset{N\to \infty }{lim}{\left(N\right)}^{3/2}{\rho }_N^{V_{tN}}\left(\lambda \right)=\frac{8}{t^{3/2}\pi }\sqrt{\lambda },\text{where}{V}_{tN}\left(x\right)=\frac{2}{t^2{N}^2}{x}^2-\frac{4}{tN}x.$$

(46)

5. Deformation by an Exponential Potential

Here we consider a deformation of the Wigner distribution by the insertion of the exponential potential. It is interesting to compare this model with the model [53] that represents planar graphs with dynamical holes of arbitrary sizes. In this model there is spontaneous tearing of the world sheet, associated with the planar graph, which gives a singularity at zero coupling constant of string interaction.

$$U\left(x\right)=V_0\left(x\right)+V_1\left(x\right)=\frac{m^2{x}^2}{2}+J{e}^{\omega x}.$$

(47)

There are 4 different choices of signs of J and $\omega$, see Figure 2. There is no instability for $J>0$, see Figure 2a,b. Only $J<0$ may produce some non trivial effects due to an appearance of potential instability, Figure 2c,d. The choice of sign of $\omega$ is irrelevant. Instability regions also appeared in Higgs-type models [54,55,56].

The shift of $\lambda \to \lambda +\delta \lambda$ in the potential (47) produces the linear term in the potential and multiplies the current J on the constant $J\to J{\mathrm{e}}^{\delta \lambda }$. We parametrize our potential as

$$V\left(x\right)=m^2\left(\frac{1}{2}{x}^2+Cx+J{\mathrm{e}}^{\omega x}\right)$$

(48)

and fix the parameters in (48) in an agreement with the measure localization on the segment $[0,\mathsf{\Lambda }]$. 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 ${\rho }_W(\lambda ,\mathsf{\Lambda })$ and ${\rho }_{e^{w\lambda }}(\lambda ,\mathsf{\Lambda },\omega )$.

The forms of non-normalized measures ${\rho }_{e^{w\lambda }}\left(\lambda \right)$ for positive and negative w are presented in Figure 3.

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}(\mathsf{\Lambda },\omega )<0$. Then for $J>J_{cr}(\mathsf{\Lambda },\omega )$ we find relation between $m^2$ and $\mathsf{\Lambda },\omega$ and J from normalization condition.

In Figure 4 the appearance of the phase transition at negative J is presented. In Figure 5a relations between $m^2$ and J for fixed $\mathsf{\Lambda }$ and $\omega$ are shown for $\mathsf{\Lambda }=6$ and in Figure 5b. More plots for different $\mathsf{\Lambda }$ one can find in the extended version of this paper [57].

We can present the eigenvalues distribution corresponding to the potential (48) as a sum. Note that the potential ${\mathrm{e}}^{\omega x}$ does not have a solution to the singular equation and does not itself defines the eigenvalues distribition, but $V=jx+J{\mathrm{e}}^{\omega x}$ does in the case $jJ\omega <0$. By ${\rho }_{e^{w\lambda }}$ we mean this distribution with corresponding choice of the linear term.

$$\begin{array}{ccc}\hfill {\rho }_{nn}\left(\lambda \right)& =& {\rho }_W(\lambda ,\mathsf{\Lambda })+J{\rho }_{e^{w\lambda }}(\lambda ,\mathsf{\Lambda },\omega )\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill {\rho }_W\left(\lambda \right)& =& \frac{1}{2{\pi }^2}\sqrt{\lambda (\mathsf{\Lambda }-\lambda )}{\int }_0^{\mathsf{\Lambda }}\frac{d\mu }{\sqrt{(b-\mu )(\mu -a)}}\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill {\rho }_{e^{w\lambda }}\left(\lambda \right)& =& \frac{\omega }{2{\pi }^2}\sqrt{\lambda (\mathsf{\Lambda }-\lambda )}{\int }_0^{\mathsf{\Lambda }}\frac{(e^{\omega \lambda }-e^{\omega \mu })d\mu }{(\lambda -\mu )\sqrt{(b-\mu )(\mu -a)}}\hfill \end{array}$$

(51)

and fix constant C and $m^2$ from the consistency condition and normalization, respectively,

$$\begin{array}{c}\hfill 0={\int }_0^{\mathsf{\Lambda }}\frac{\lambda +C+J\omega e^{\omega \lambda }}{S_{\mathsf{\Lambda }}}d\lambda ,\frac{1}{m^2}={\int }_0^{\mathsf{\Lambda }}\left({\rho }_W(\lambda ,\mathsf{\Lambda })+J{\rho }_{e^{w\lambda }}(\lambda ,\mathsf{\Lambda },\omega )\right)d\lambda \end{array}$$

(52)

The forms of non-normalized measures ${\rho }_{e^{w\lambda }}\left(\lambda \right)$ for positive and negative w presented are presented in Figure 3. We see that these two segment distributions are nonsymmetric under the centre of the segment. The distribution of eigenvalues for the case of negative $\omega$ is pressed to the left boundary of the segment, and the for the case of positive $\omega$ 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}(\mathsf{\Lambda },\omega )<0$. Then $J>J_{cr}(\mathsf{\Lambda },\omega )$ we find relation between $m^2$ and $\mathsf{\Lambda },\omega$, and J from normalization condition.

In Figure 4 is presented the appearance of the phase transition at negative J for different $\mathsf{\Lambda }$. We see that for chosen parameters, $m^2=1,\omega =1$ the critical $J_{cr}$ decreases with increasing $\mathsf{\Lambda }$. To find the real mass that supports the normalized solution,

$$\rho \left(\lambda \right)=n\left(\lambda \right){\rho }_{nn}\left(\lambda \right)$$

(53)

we find $n\left(\lambda \right)$ from the normalization condition, so that $n\left(\mathsf{\Lambda }\right){\int }_0^{\mathsf{\Lambda }}{\rho }_{nn}\left(\lambda \right)d\lambda =1$ and assume that the mass is given by $m^2=n\left(\mathsf{\Lambda }\right)$.

In Figure 5a the dependence of $m^2$ on J for fixed $\mathsf{\Lambda }$ and $\omega$ is shown for different values of $\mathsf{\Lambda }$. Here $J>J_{cr}\left(\mathsf{\Lambda }\right)$. We see that mass (in our parametrization of the potential of the model) decreases with increasing J. The slope of the mass is more fast for larger $\mathsf{\Lambda }$.

In Figure 5b 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 $\mathsf{\Lambda }$.

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]}\to S_{[0,b-a]}$ we have to make a shift in the potential, $V\left(x\right)\to V(x+a)$. For the quadratic potential this shift produces the linear term $\mathsf{\Delta }V\left(x\right)=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

$$\begin{array}{ccc}\hfill V_1^{\prime }\left(\lambda \right)& =& 2{⨏}_S\frac{{\rho }_{V_1}\left({\lambda }^{\prime }\right)}{\lambda -{\lambda }^{\prime }}d{\lambda }^{\prime },{\int }_S{\rho }_{V_1}\left(\lambda \right)d\lambda =1\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill V_2^{\prime }\left(\lambda \right)& =& 2{⨏}_S\frac{{\rho }_{V_2}\left({\lambda }^{\prime }\right)}{\lambda -{\lambda }^{\prime }}d{\lambda }^{\prime },{\int }_S{\rho }_{V_2}\left(\lambda \right)d\lambda =1,\hfill \end{array}$$

(55)

Note that in both integrals the segment is the same. Taking

$${\rho }_T=\frac{1}{2}({\rho }_{V_1}+{\rho }_{V_2})$$

(56)

we can claim that ${\rho }_T$ solves equations

$$\begin{array}{ccc}\hfill V_T^{\prime }\left(\lambda \right)\equiv \frac{1}{2}(V_1^{\prime }\left(\lambda \right)+V_2^{\prime }\left(\lambda \right))& =& 2{⨏}_S\frac{{\rho }_T\left({\lambda }^{\prime }\right)}{\lambda -{\lambda }^{\prime }}d{\lambda }^{\prime },{\int }_S{\rho }_T\left(\lambda \right)d\lambda =1\hfill \end{array}$$

(57)

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,\mathsf{\Lambda }]$ by fine tuning the linear term in the potential. Then one goes to the limit $N\to \infty$ and after that $\mathsf{\Lambda }\to \infty$.

6. Matrix Model for JT Gravity

6.1. Potentials for Non-Normalized Density Distribution

For fixed $\mathsf{\Lambda }$ we consider the eigenvalue distribution normalized to 1

$${\rho }_{norm,1}\left(E\right)=n\left(\mathsf{\Lambda }\right){\rho }_{nn}\left(E\right),{\rho }_{nn}\left(E\right)=\frac{1}{{\left(2\pi \right)}^2}sinh\left(2\pi \sqrt{E}\right),E>0.$$

(58)

This form of the eigenvalue distribution in the Sachdev-Ye-Kitaev (SYK) model has been obtained in [58,59] and is nothing but the Bethe formula for the nuclear level density [60]. For large $\mathsf{\Lambda }$ one has $n\left(\mathsf{\Lambda }\right)\approx 8{\pi }^3{e}^{-2\pi \sqrt{\mathsf{\Lambda }}}/\sqrt{\mathsf{\Lambda }}$. The distribution normalized to N has the form

$${\rho }_{d.s.}^{\mathrm{nnorm}}\left(E\right)\equiv {\rho }_{norm,N}\left(E\right)=\frac{{\mathrm{e}}^{S_0}}{{\left(2\pi \right)}^2}sinh\left(2\pi \sqrt{E}\right),E>0.$$

(59)

compare with (19). It is evident that ${\mathrm{e}}^{S_0}=Nn\left(\mathsf{\Lambda }\right)$

To recover the potential that supports the distribution ${\rho }_{norm,1}\left(E\right)$, we write

$$\begin{array}{ccc}\hfill V\left(\mu \right)& =& n\left(\mathsf{\Lambda }\right)V_{nn}\left(\mu \right),\hfill \end{array}$$

(60)

where $V_{nn}^{\prime }\left(\mu \right)$ is defined by

$$\begin{array}{ccc}\hfill V_{nn}^{\prime }\left(\mu \right)& =& \frac{1}{{\left(2\pi \right)}^2}{\int }_0^{\mathsf{\Lambda }}\frac{sinh2\pi \sqrt{\lambda }}{\mu -\lambda }d\lambda .\hfill \end{array}$$

(61)

This gives

$$\begin{array}{ll}{V}_{nn}^{\prime }\left(\mu \right)& =\frac{1}{{\left(2\pi \right)}^2}[e^{-2\pi \sqrt{\mu }}\left(e^{4\pi \sqrt{\mu }}\mathrm{Ei}\left(2\pi \left(\sqrt{\mathsf{\Lambda }}-\sqrt{\mu }\right)\right)+\mathrm{Ei}\left(2\pi \left(\sqrt{\mathsf{\Lambda }}+\sqrt{\mu }\right)\right)\right)\\ & - e^{-2\pi \sqrt{\mu }}\left(\mathrm{Ei}\left(2\pi \left(\sqrt{\mu }-\sqrt{\mathsf{\Lambda }}\right)\right)+e^{4\pi \sqrt{\mu }}\mathrm{Ei}\left(-2\pi \left(\sqrt{\mathsf{\Lambda }}+\sqrt{\mu }\right)\right)\right)]\end{array}$$

(62)

Here Ei is the exponential integral that for real non zero values of x is defined as

$$Ei\left(x\right)=-{⨏}_{-x}^{\infty }\frac{e^{-t}}{t}dt.$$

(63)

The potential up to a constant is

$$V\left(\mu \right)={\int }^{\mu }{V}^{\prime }\left(\lambda \right)d\lambda$$

(64)

and it is presented in Figure 6.

6.2. Effective Energy

The effective energy $E_{eff}$, [18], is evaluated on the normalized density and it is given by the following formula

$$\begin{array}{ccc}\hfill E_{eff}& =& n^2\left(\mathsf{\Lambda }\right){\mathcal{E}}_{eff}=\frac{{\mathrm{e}}^{2S_0}}{N^2}{\mathcal{E}}_{eff},\hfill \\ \hfill {\mathcal{E}}_{eff}& =& \int d\lambda {\rho }_{nn}\left(\lambda \right)V_{nn}\left(\lambda \right)-\int \int d\lambda d{\lambda }^{\prime }{\rho }_{nn}\left(\lambda \right){\rho }_{nn}\left({\lambda }^{\prime }\right)ln{(\lambda -{\lambda }^{\prime })}^2\hfill \end{array}$$

(65)

In Figure 7 the effective action $E_{eff}$ as function of $\mathsf{\Lambda }$ is shown. We see that it can be approximated by a log $10.60+0.959logx$. The next approximation is given by the double-log, $10.76+1.050logx-0.393loglogx$.

6.3. Phase Transition

To study the phase transition we consider the deformation of ${\rho }_{sinh\sqrt{\lambda }}$ by ${\rho }_{e^{w\lambda }}$

$$\begin{array}{c}\hfill {\rho }_{nn}\left(\lambda \right)={\rho }_{sinh\sqrt{\lambda }}(\lambda ,\mathsf{\Lambda })+J{\rho }_{e^{w\lambda }}(\lambda ,\mathsf{\Lambda },\omega )\end{array}$$

(66)

It is interesting to compare this density with the density

$$\begin{array}{c}\hfill {\rho }_{nn,s}\left(\lambda \right)={\rho }_{sinh\sqrt{\lambda (1-\lambda /\mathsf{\Lambda })}}(\lambda ,\mathsf{\Lambda })+J{\rho }_{e^{w\lambda }}(\lambda ,\mathsf{\Lambda },\omega ),\end{array}$$

(67)

where

$${\rho }_{sinh\sqrt{\lambda (1-\lambda /\mathsf{\Lambda })}}(\lambda ,\mathsf{\Lambda })=\frac{n_s\left(\mathsf{\Lambda }\right)}{{\left(2\pi \right)}^2}sinh2\pi \sqrt{\lambda (1-\lambda /\mathsf{\Lambda })}.$$

(68)

In Figure 8 we plot the density (66) and (67) for $\omega =-1$ and different J and $\mathsf{\Lambda }$. We see that for negative $J>J_{cr}\left(\mathsf{\Lambda }\right)$ on the segment $(0,\mathsf{\Lambda })$ domains, where ${\rho }_{nn}$ and ${\rho }_{nn,s}$ become negative, appear. Note that the critical value $J_{cr}$ are the same for ${\rho }_{nn}$ and ${\rho }_{nn,s}$. In Figure 9 the locations of critical points depending on $\mathsf{\Lambda }$ are shown. This broken of positivity of the density means an instability of the solution (phase) under a perturbation by the exponential potential. We interpret this instability as an appearance of a phase transition at $J=J_{cr}$.

The value of $J_{cr}$ depends on $\omega$ and $\mathsf{\Lambda }$. In particular, the dependence of $J_{cr}$ on $\mathsf{\Lambda }$ at $\omega =-1$ is presented at Figure 9a). The values of $J_{cr}$ for different $m^2$ and for the same set of $\mathsf{\Lambda }$ as at Figure 9a) are presented at Figure 9b).

7. Discussion and Conclusions

The generating functional for the correlation functions of the boundaries of the Riemann surfaces 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 $\mathfrak{Z}\left(J\right)$ considered in this paper, gives completion also for infinitely many number of boundaries.

By using this formulation several matrix models, including double scaling limit for the Gaussian model relevant to topological gravity, cubic model, and JT gravity are investigated in this paper. In all these cases, indications to phase transitions are found. In particular, in JT gravity the phase transition in the gas of baby universes has been establised. By analogy with phase transitions between different states of matter, it can be expected that the phase transition discussed above describes the transition from the gaseous state of the JT multiverse to the 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\left(x\right)$ is obtained from the partition function of the model just by shift $V\left(x\right)\to V\left(x\right)-\tilde{J}\left(x\right)$ where $\tilde{J}\left(x\right)$ is the Laplace transform of $J\left(\beta \right)$. Baby universes should look as point-like objects with very small length of the boundary. Therefore the distribution $J\left(\beta \right)$ should have a peak somewhere near zero. In our model with $J\left(\beta \right)=J\delta (\beta -\omega )$ baby universes correspond to small $\omega$. The shifted potential in this case is $V\left(x\right)-J{e}^{-\omega x}$. Note that one could expect $<Z{\left(\beta \right)}^n>\approx <Z\left(\beta \right){>}^n$ for large n similarly to the infinite replica number limit considered in [61].

We have obtained that in all considered 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 $\mathsf{\Lambda }$. 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 [62,63], as well matrix models with other potentials considered in [64].

The holographic boundary dual of the JT gravity in spacetime with M boundaries is the M-replica SYK model in the low energy limit [65,66,67,68]. The studies of [61,69,70,71,72,73] show that the nonperturbative completion of SYK involves nontrivial replica-nondiagonal saddle points. The replica-nondiagonal structures in SYK with replica interaction [73,74,75,76] demonstrate nontrivial phase structures and symmetry breaking patterns.

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 [69,70,73]. 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\left(x\right)$ (cosmological daemon) in cosmology [77], 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 [32] and Coleman’s approach to the cosmological constant problem [36]. 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 [78] 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 [37]. It is interesting to find an analog of this in the context of JT gravity.

In the study of SYK from gravity perspective, a crutial role plays consideration of wormholes in JT gravity [74,76,79,80,81,82,83]. A interplay between baby universes and wormholes could lead to nontrivial effects and this subject requires further investigations.

Using the recent result of [84] according which JT gravity in $d{S}_2$ is an analytic continuation of JT gravity in Euclidean $Ad{S}_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 compare the Hartle-Hawking construction of string baby universes related with $Ad{S}_2\times S^2$ geometry, using free fermionic formulation [85] with the fermionic interpretation of determinant in matrix models, and also find a slot for baby universes in an quasi-classical [86,87] or exact quantization of JT proposed in [88].