# A Flow in the Forest

Alexander Gorsky, Vladimir Kazakov, Fedor Levkovich-Maslyuk, Victor
Mishnyakov

► To cite this version:

Alexander Gorsky, Vladimir Kazakov, Fedor Levkovich-Maslyuk, Victor Mishnyakov. A Flow in the Forest. Journal of High Energy Physics, 2023, 2023 (03), pp.67. 10.1007/JHEP03(2023)067. hal-03837693

HAL Id: hal-03837693

https://hal.science/hal-03837693v1

Submitted on 12 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A flow in the forest

Alexander Gorsky,ª Vladimir Kazakov,ᵇ Fedor Levkovich-Maslyukᶜ,<fnref n="1" />
and Victor Mishnyakovᵈ,ᵉ,*f,g

E-mail: agorsky3@gmail.com, kazakov@lpt.ens.fr,
fedor.levkovich@gmail.com, mishnyakovvv@gmail.com

**ABSTRACT:** Using the matrix-forest theorem and the Parisi-Sourlas trick we formulate and solve a one-matrix model with non-polynomial potential which provides perturbation theory for massive spinless fermions on dynamical planar graphs. This is a lattice version of 2d quantum gravity coupled to massive spinless fermions. Our model equivalently describes the ensemble of spanning forests on the same graphs. The solution is formulated in terms of an elliptic curve. We then focus on a near-critical scaling limit when both the graphs and the trees in the forests are macroscopically large. In this limit we obtain one-point scaling functions (condensates), parameterized in terms of the Lambert function. Our results provide a rare example where one can explore the flow between two gravity models — in this case, the theories of conformal matter coupled to 2d gravity with $c = -2$ (large trees regime) and $c = 0$ (small trees regime). We also compute the disc partition functions with Dirichlet and Neumann boundary conditions in the same critical limit.

**KEYWORDS:** 2D Gravity, Matrix Models

ARXIV ePRINT: 2210.07176

<footnote>ªInstitute for Information Transmission Problems,
Moscow 127994, Russia</footnote>

<footnote>ᵇLaboratoire de Physique de l'École Normale Supérieure,
CNRS, Université PSL, Sorbonne Universités,
24 rue Lhomond, 75005 Paris, France</footnote>

<footnote>ᶜUniversité Paris Saclay, CNRS, CEA, Institut de physique théorique,
91191 Gif-sur-Yvette, France</footnote>

<footnote>ᵈLebedev Physics Institute,
Moscow 119991, Russia</footnote>

<footnote>ᵉNRC "Kurchatov Institute",
Moscow 123182, Russia</footnote>

<footnote>f MIPT,
Dolgoprudny 141701, Russia</footnote>

<footnote n="9">⁹Institute for Theoretical and Mathematical Physics, Lomonosov Moscow State University,
Moscow 119991, Russia</footnote>

<footnote n="1">¹Also at Institute for Information Transmission Problems, Moscow 127994, Russia.</footnote>
# Contents

1 Introduction 1

2 Definition of the model 4

3 Derivation of the matrix model 6

3.1 Derivation via Parisi-Sourlas representation 6

3.2 Combinatorial derivation of the matrix model 8

3.3 Boundary conditions for disc partition functions 9

4 Solution of the matrix model in the planar limit: saddle point 11

4.1 Fixing parameters 12

4.2 Exact results for correlators 14

5 The critical line 15

5.1 Merging of solutions and the physical region 16

5.2 Looking for further singularities 16

5.3 Decoupling of the heavy matter: $M \to \infty$ limit 17

5.4 Critical line at $M \to 0$ 18

6 New double scaling limit for 1-pt functions 19

6.1 The $\frac{1}{N}\langle\text{tr } \phi\rangle$ one-point function 21

6.1.1 Limiting regimes 22

6.1.2 Rational representation and rescaled form 23

6.1.3 “Asymptotic freedom” on the second sheet 24

6.2 The $\frac{1}{N}\langle\text{tr } \phi^3\rangle$ one-point function 24

7 Disc partition function 26

7.1 The resolvent $G(x)$ for tr $X^k$ correlators 27

7.1.1 Limiting cases 27

7.2 The resolvent $H(r)$ for tr $\phi^k$ correlators 28

7.2.1 Limiting cases 29

8 Conclusion 30

9 Further directions 31

9.1 Double scaling limit and sum over topologies 31

9.2 Generalization to other critical flows 32

9.3 Fragmentation of the RRG into finite number of trees and many-body localization 32

9.4 ZZ-brane instantons for one-cut solution 33

9.5 Kesten-McKay distribution and criticality without planarity 34

9.6 Analogy with QCD matrix model 35
<table><tr><td>A</td><td>Combinatoric explanation of the Kirchoff-Tutte matrix-tree theorem</td><td>36</td></tr><tr><td>B</td><td>Some properties of the spectral determinant</td><td>37</td></tr><tr><td>C</td><td>Elliptic functions</td><td>37</td></tr><tr><td></td><td>C.1 Some integrals</td><td>38</td></tr><tr><td></td><td>C.2 Relations used for integration of elliptic functions</td><td>38</td></tr><tr><td>D</td><td>Expansions at small M and p</td><td>38</td></tr><tr><td>E</td><td>Lambert function and brane insertion</td><td>39</td></tr><tr><td></td><td>E.1 Two examples with Lambert function</td><td>39</td></tr><tr><td></td><td>E.2 Lambert in the forest</td><td>42</td></tr></table>

# 1 Introduction

Matrix model approach to 2d quantum gravity, first proposed in [1–3], has undergone since then many important and diverse developments (see [4–6] for early results, and a more recent review [7]). They range from the study of critical exponents in the presence of conformal matter with a central charge $c < 1$ [8–12] to the construction of quantum field theory for non-critical strings via the double scaling limit [13–16]. As was first proposed in [4], and then generalized to the double-scaling limit [17–20], a special, limiting $c = 1$ case of non-critical string theory can be studied via the matrix quantum mechanics (MQM) [21]. An interesting peculiarity of $c = 1$ case are the integer critical exponents, resulting in the logarithmic scaling, as well as in the presence of Berezinsky-Kosterlitz-Thouless vortex excitations on the world-sheet of the $c = 1$ string [22–24], whose condensation can lead to the black hole formation [25]. The $c = 1$ MQM model has served as a basis of formulation and study of type-0 2d superstring models [26]. The other interesting cases with integer critical exponents and a logarithmic scaling are various models with matter central charges $c = -2$ and $c \to 0$, in particular, the $Q \to 0, 1$ limits of the Q-state Potts model [10, 27, 28]<fnref n="1" /> and the $n \to 0, -2$ limits of $O(n)$ spin model [11, 30–34],<fnref n="2" /> both on planar graphs. More recently the matrix models for the pure JT gravity [37] and JT gravity coupled with the massive scalar [38] have been suggested.

Beyond the scaling properties, certain important physical quantities have been studied and computed for various central charges of the matter $c$, such as the disc and cylinder partition functions and the multi-loop correlators for the minimal $(p, q)$ series, with $c = 1 - 6(p-q)^2/pq$ [31]. However, the questions of universal flows between critical points with different central charges of the matter hasn’t yet been carefully studied. A particular, well-studied case is the non-unitary $(2, 2n - 1)$ series [39, 40], where the precise map from the

<footnote n="1">¹The solution of the saddle point problem for the matrix model of Potts spins on random graphs was found by I.Kostov and V.K., published in [29].</footnote>

<footnote n="2">²See also [35, 36].</footnote>
matrix model to the minimal Liouville theory of 2d gravity with conformal matter [41–44], completing and old proposal [45], was formulated. The study on the matrix side is simplified here by formulation in terms of the one-matrix model with a polynomial potential [12, 46]. Certain flows have been observed between various p, q points within the $O(n)$ model of Kostov [31–33]. Apparently, these are only particular examples, and a thorough study of the whole variety of such flows, especially around logarithmic criticalities mentioned above, is largely missing.

The main subject of this paper is the study of such flows, and of the associated physical quantities (one point functions, disc partition functions) for an interesting case of 2d QG coupled to $c_M = -2$ matter, with a logarithmic scaling. This system can be represented as the 2d QG in the presence of complex Grassmannian scalar matter field $\theta(\sigma)$. Such model has been studied in both the matrix model approach, where the sum over geometries is represented by the sum over planar graphs [3, 30, 47], as well as in continuous Liouville field theory approach [48–50].

Massive spinless fermions on the fixed non-fluctuating torus have been explored in [51].<fnref n="3" /> If the four-fermion interaction for massive fermions is added at a particular value of four-fermion coupling this model on a fixed graph is identical to $q \to 0$ Potts model and $O(-2)$ sigma model [53]. The $q \to 0$ Potts model coupled to 2d gravity has been considered in [10, 54, 55] via matrix model representation. At a particular value of the four-fermion coupling the model gets reduced to the statistical ensemble of unrooted spanning forests. The matrix potential is non-analytic and a complicated pattern of phase transitions has been discovered [55].

In this work, we will study a generalization of this $c = -2$ QG matrix model where the matter field $\psi(\sigma)$ is massive. The disc partition function of this 2d QG theory in continuum is given by the functional integral

$$Z(\Lambda_0, \eta, m_0) = \int Dg(\sigma) \int D\psi(\sigma) D\bar{\psi}(\sigma) \exp \left[ - \int_D d^2\sigma \sqrt{g} \left( g^{\alpha\beta} \partial_\alpha \bar{\psi} \partial_\beta \psi + m_0^2 \bar{\psi} \psi + \Lambda_0 + \kappa R \right) - \eta_0 \int_{\partial D} ds [g(\sigma(s))]^{1/4} \right] \quad (1.1)$$

where $g$ is the 2d metric, $\Lambda_0$ is the bulk cosmological constant, $\eta_0$ is the boundary cosmological constant, $R$ is the Gaussian curvature and $\kappa = 2 - 2g$ is the Euler characteristics which in our paper will correspond to the topologies of the sphere $\kappa = 2$ or of the disc $\kappa = 1$. At zero mass $m = 0$, the matter field has the conformal central charge $c_M = -2$.

Using the matrix-forest theorem [56] and Parisi-Sourlas trick [57] we will derive the corresponding matrix model, with a specific non-polynomial potential, with two-parameters — the bare mass of fermions and the bare cosmological constant. In our matrix formulation, the functional integral over metrics will be represented by the sum over $\phi^3$ (i.e. trivalent) planar graphs, and the matter fields $\psi_j$ will be placed at the sites $j = 1, 2, \dots, V$ of each graph of size $V$, see eq. (2.1).

If one turns on the mass, at large distances (small cosmological constant) the matter is screened out and one “flows” in the infrared to the pure gravity $c_M = 0$ point, with a

<footnote n="3">³See also [52] for a related model of fermions coupled to gravity.</footnote>
renormalized cosmological constant. The study of this flow by the matrix model tools is the principal task of this paper. Solving the corresponding one-matrix model with a special, non-polynomial potential, we will establish the universal scaling functions describing the one-point function and the disc partition functions (with various boundary conditions).

We should stress here that, whereas a lot is known about the universal critical behavior of various physical quantities (correlation functions, disc and annulus partition functions, etc) at a given central charge of conformal matter [11, 28, 31, 34, 58–61], the universal flows between the fixed points with various central charges are not that well studied. In this paper, we present our results for the one-point function and the disc partition functions for such a flow, between the $c = -2$ and $c = 0$ fixed points.

As an example, we present in the introduction the t-parametric representation of the one point function of the type $\Phi_s = \langle \bar{\psi}\psi \rangle : (\text{see section 6 for details}^4)$

$$ \Phi_s \sim 2(\mathcal{J} + \text{const})(t^2 - 2t) + 3t^2 - 2t^3 \quad (1.2) $$

where (in certain units)

$$ \mathcal{J} = t - \log(m_*^2 t), \qquad (\mathcal{J} = \frac{\Lambda_*}{m_*^2}) \quad (1.3) $$

and $\Lambda_*, m_*^2$ are the lattice analogues of $\Lambda_0, m_0^2$ of (1.1), and const in (1.2) indicates a numerical constant (given in section 6). This equation is closely related upon exponentiation to the Lambert function which has appeared previously in related contexts. It was found in [62] that the Lambert function corresponds to the peculiar brane which provides the generating function for the multiple boundaries in the Airy limit of topological gravity at genus zero. Since the brane insertion shifts the closed moduli [63] it defines these shifts as well. Another example of the Lambert function providing the moduli shift has been identified in the topological string on $CP_1$ coupled to topological gravity [64]. In our study the Lambert function will play a similar role. Lambert function also emerges in various formulas related to asymptotics of Hurwitz numbers [65] and as the spectral curve of the type B topological string on $C^3$ in the limit of one-leg infinite framing.

One can easily see that this scaling function interpolates between pure 2d QG $c = 0$ regime $\Phi_s \sim (\Lambda_c - \Lambda)^{3/2}$ when $\Lambda \ll m^2$ and the 2d QG with $c = -2$ matter regime $\Phi_s \sim \Lambda^2 \log \Lambda$ when $\Lambda \gg m^2$, described by spanning trees on large planar graphs.

In addition to (1.2), we also computed another 1-pt function of a similar type, which turns out to differ only by the value of the additive shift of $\mathcal{J}$ in (1.2). This suggests a universality-like property that would be interesting to explore further.

We will also derive two universal disc partition functions describing the flow between the $c = -2$ and $c = 0$ fixed points, for massive spinless worldsheet fermions with Neumann and Dirichlet boundary conditions. The results are given by equations (7.6) and (7.16), respectively. In the limits mentioned above for the one-point function $\Phi_s$, they reproduce the known asymptotic behavior at the $c = -2$ and $c = 0$ fixed points [32].$^5$

<footnote>$^4$In the notation used there this correlator corresponds to tr $\phi^3$.</footnote>

<footnote>$^5$In the cited paper, the author presents another example of a flow between the same critical points, which he derives from the $O(-2)$ model on planar graphs. Apparently, that flow is different, it leads to a different universal scaling function.</footnote>
If we drop off the planarity condition, the sum over the graphs is nothing but the random regular graph (RRG) ensemble which enjoys some exact results in the limit of large number of nodes. In particular we can utilize the famous Kesten-McKay(KM) distribution for the resolvent and spectral density of RRG ensemble [66, 67]. Recently the RRG has attracted a lot of attention from the very different perspective which has nothing to do with the quantum gravity. It is considered as the toy model for the Hilbert space of the interacting many-body system and according to the conjecture [68] the one-particle Anderson localization of spinless fermion on RRG with diagonal disorder is equivalent to the many-body localization in the physical space-time (see [69] for the recent review and references therein). It was argued that the fragmentation of the Hilbert space graph into some number of “trees” is one of the key mechanisms of transition to the MBL phase (see [70] for review and references therein). It is not clear whether the problem we are solving in this work — the massive spineless fermions on *planar* graphs — is directly related to the abovementioned RGG ensembles of *generic* graphs. Possibly, the similar phenomena can occur in our model in the appropriate double scaling limit, summing up the large graphs (near criticality) of all topologies [13–16].

The paper is organized as follows. In section 2 we define our model of massive spinless fermions on the random planar graphs. It can be equally considered as the theory of anticommuting bosons. In section 3 we derive the one-matrix model for our theory using two different strategies. First we generalize the Parisi-Sourlas trick for the massive case and obtain the generalization of the matrix model for $c = -2$ theory. Secondly we shall utilize the matrix-forest theorem [56] for the massive determinant of the graph Laplacian and develop a combinatorial derivation of the matrix model. This derivation is somewhat similar to the derivation of the matrix model in [54, 55] for a slightly different model of spinless fermions coupled to 2d gravity. In section 4 we find the one-cut solution of this one-matrix model in the planar limit via the standard tools, in terms of the elliptic curve. In section 5 we derive the critical curve which bounds the regime of validity for the one-cut solution in the two-dimensional parameter space of the model (related to the $(\Lambda, m^2)$ space of (1.1) in the critical regime). In section 6 we consider in more detail the limit of small fermion mass and define a new double scaling regime in the space of two parameters. We will compute two 1-pt functions and find they are given by closely related scaling functions that interpolate between $c = -2$ and $c = 0$ critical regimes and are related to the Lambert function. In section 7 we derive the universal disc partition functions of the model for the Dirichlet and Neumann boundary conditions (involving the parameter analogous to the boundary cosmological constant $\eta_0$ from (1.1)). In section 8 we summarize the main results of the paper and lastly section 9 discusses several possible topics for future research. The appendices contain various technical details, while in appendix E we discuss the role of the Lambert function in our model and related contexts. We also attach a supplementary Mathematica notebook with some of the lengthy explicit results.

## 2 Definition of the model

In the spirit of the discrete, random lattice approach to 2d QG, we will construct a matrix model for which the perturbative expansion for the free energy, combined with the $1/N$
expansion, is given as the following sum over planar graphs G:

$$ \log \zeta = \sum_{G} N^{2-2g} \lambda^{|G|} \int \prod_{i \in G} d^2 \psi_i \prod_{<ij> \in G} e^{-(\bar{\psi}_i - \bar{\psi}_j)(\psi_i - \psi_j)} \prod_{i \in G} e^{-m^2 \bar{\psi}_i \psi_i} \quad (2.1) $$

where in the first exponent the propagators $e^{-(\bar{\psi}_i - \bar{\psi}_j)(\psi_i - \psi_j)}$ mimic the kinetic term of fermions in the action (1.1), the last product introduces the mass of fermions, $\log \lambda$ plays the role of bare bulk cosmological constant and $g$ represents the genus of the discretized worldsheet. In this paper we will focus on the strictly planar case. We see that this partition function looks as a lattice analogue of the partition function (1.1).<fnref n="6" />

Integrating in (2.1) over $\psi$'s we rewrite this partition function in the form

$$ Z \equiv \log \zeta = \sum_{G} N^{2-2g} \lambda^{|G|} \det[m^2 + \Delta(G)] \quad (2.2) $$

where $\Delta = -Q + A$ is the graph Laplacian with $Q = \text{diag}\{q_1, q_2, \dots\}$ where $q_j$ are the valencies of vertices and $A$ is the adjacency matrix. In what follows we will use, for definiteness, the 3-valent graphs, i.e. $q_i = 3$. For $m^2 = 0$ the model was solved in [3, 47] using the spanning trees representation and interpreted there as 2d QG with $c = -2$ matter. We will generalize it to the $m^2 \neq 0$ case and use for it the Parisi-Sourlas approach proposed in [30, 71].

Notice that the partition function of the $c = -2$ model coupled to gravity is given by the sum over graphs of the determinant of the graph Laplacian with the zero mode removed, which we denote by $\det'\Delta(G)$. Nicely, it is related in a simple way with the massive determinant<fnref n="7" /> since

$$ \det'\Delta(G) = \frac{d}{dm^2} \det(\Delta(G) - m^2)_{m^2=0}. \quad (2.3) $$

Hence we have an exact relation between the $c = -2$ and massive partition functions

$$ Z_{c=-2}(\lambda) = \frac{d}{dm^2} Z(\lambda, m^2 = 0) \quad (2.4) $$

The partition function (2.2) can be usefully rearranged via the generalization of the Kirchhoff matrix-tree theorem to the matrix-forest theorem for characteristic polynomial for the graph Laplacian which was obtained in [56] (see also [51])

$$ \det[m^2 + \Delta(G)] = \sum_{F=(F_1 \dots F_l) \in G} \prod_{i=1}^{l} m^2 V(F_i) \quad (2.5) $$

where $V(F_i)$ is the number of nodes in the tree $F_i$ and $m^2$ is the generating parameter for the number of trees in the forest. That is, our model can be considered as the partition function of spanning rooted forests interacting with 2d gravity.

<footnote n="6">⁶The boundary term for the disc partition function will be introduced later within the matrix model formalism.</footnote>

<footnote n="7">⁷The properties of this spectral determinant are discussed in appendix B.</footnote>
Notice that the work [54] considered a related model which involves massive spinless fermions supplemented with the four-fermion term coupled to 2d gravity. In that case the generalization of the matrix-forest theorem exists [53] but the model, contrary to our case, reduces to the spanning *unrooted* forests on fluctuating surface. Instead we have *rooted* forests and no four-fermion interaction.

Since the degrees of all nodes are equal, our model can be considered as a particular case of the random regular graph (RRG) ensemble with general exponential weight:

$$Z = \sum_{\text{RRG}} \exp(-\text{Tr}V[A(G)]) \qquad (2.6)$$

where the sum goes w.r.t. graphs G with fixed number of nodes with the same degree. This partition function describes the microcanonical ensemble with fixed area of the surface, contrary to the canonical ensemble with the cosmological constant. Usually in consideration of RRG ensemble no planarity condition for the graphs is implied. Since Tr here is taken w.r.t. indices of adjacency matrix of a graph, Tr$A^n$ gives the number of cycles of length $n$ on the graph. In our case, the potential will have a specific, determinant form $V(A) = -\log((m^2 - 3)\mathbb{I} - A)$, and furthermore we restrict to only planar graphs among the whole 3-valent RRG ensemble.

## 3 Derivation of the matrix model

In this section we will derive the Hermitian one-matrix model which gives our partition (2.1) and provides the interpolation between pure gravity and $c = -2$ conformal matter coupled with 2d quantum gravity. We will use two different strategies which generalize two approaches used for the $c = -2$ case. First, extending the approach of [30, 51] we will apply a massive version of the Parisi-Sourlas trick. Second, extending the derivation of the matrix model for the unrooted spanning forests in [55]$^8$ we will derive the matrix model for our case of spanning forests. We will show that the matrix model potentials we find in the two derivations coincide.

### 3.1 Derivation via Parisi-Sourlas representation

Using the Parisi-Sourlas approach [30, 49, 57, 73] we can represent our partition function as a super-matrix integral of the form

$$\zeta = \int D^{N^2} \Phi(\psi) e^{-NS(\Phi)} = \int d^{N^2} \phi d^{2N^2} \theta d^{2N^2} \epsilon e^{-NS(\Phi)} \qquad (3.1)$$

where the functional integration is performed w.r.t. the matrix superfield

$$\Phi(\psi) = \phi + \bar{\psi}\theta + \psi\bar{\theta} + \bar{\psi}\psi\epsilon, \qquad (3.2)$$

where $\bar{\theta}, \theta$ are the grassmann valued matrix fields, $\bar{\psi}, \psi$ are anticommuting parameters and $\epsilon$ and $\phi$ are commuting matrix fields. Note that the matrix superfield $\Phi$ is taken to be hermitian, hence $\phi$ is an hermitian matrix while $\epsilon$ is anti-hermitian. For our case, the matrix

<footnote>$^8$See also [72].</footnote>
model action producing the Feynman diagram expansion (2.1) with cubic ($q=3$) vertices is given by

$$S(\Phi) = \text{tr} \int d^2\psi \left( -\frac{1}{2}\Phi^2(\psi) + \frac{1}{2}\partial_\psi\Phi(\psi)\partial_{\bar{\psi}}\Phi(\psi) + \frac{\lambda}{3}e^{-m^2\bar{\psi}\psi}[\Phi(\psi)]^3 \right) \quad (3.3)$$

where $d^2\psi = d\psi d\bar{\psi}$. The superfield propagator is chosen to reproduce the kinetic term of the worldsheet fermion action, while the coefficient of the cubic vertex is responsible for the mass term. To reproduce the worldsheet action, notice that the kinetic term of the action (3.3) can be easily inverted, using:

$$(-1 - \partial_{\psi} \partial_{\bar{\psi}}) \exp\{-(\bar{\psi} - \bar{\chi})(\psi - \chi)\} = (\bar{\psi} - \bar{\chi})(\psi - \chi) = \delta^{(2)}(\psi - \chi). \quad (3.4)$$

This means that the propagator is indeed the desired function of the grassmannian distance:

$$\langle \Phi_{ij}(\psi) \Phi_{kl}(\chi) \rangle_{\lambda=0} = \delta_{il} \delta_{jk} \exp\{-(\bar{\psi} - \bar{\chi})(\psi - \chi)\}.$$

Here the expectation value is with respect to the Gaussian part of the action (3.3). Now we would like to integrate out the grassmannian variables to obtain a bosonic matrix model that would generate the same expansion (2.1). Taking the integrals we obtain:

$$S(\Phi) = \text{tr} \left[ -\frac{1}{2}\epsilon^2 - (\phi - \lambda\phi^2)\epsilon - \frac{1}{3}\lambda m^2 \phi^3 - \bar{\theta}\theta + 2\lambda\bar{\theta}\theta\phi \right]. \quad (3.5)$$

After integration over auxiliary fields $\epsilon$ (the positive sign of $\epsilon^2$ in $-S(\Phi)$ is due to anti-hermiticity of the integration variable) and $\theta$ we get:

$$\zeta = \int d^{N^2} \phi \det(1 - 2\lambda\phi) e^{N\text{tr}[-\frac{1}{2}(\phi - \lambda\phi^2)^2 + \frac{1}{3}\lambda m^2 \phi^3]} . \quad (3.6)$$

We have dropped an irrelevant overall numerical prefactor in front of the partition function. Finally doing a change of variables $X = \phi - \lambda\phi^2$ we come to the one-matrix model partition function:

$$\zeta = \int d^{N^2} X \exp \left( N \text{tr} \left[ -\frac{1}{2} X^2 + \frac{\lambda m^2}{3} (\phi(X))^3 \right] \right) \quad (3.7)$$

with $\phi$ related to $X$ by

$$\phi(X) = \frac{1}{2\lambda} \left( 1 - \sqrt{1 - 4\lambda X} \right) \quad (3.8)$$

where we have selected the proper root of the quadratic equation. Thus we see that the result is a 1-matrix model with a particular non-polynomial potential. This opens the way to solve the model using standard techniques. In the next subsection we will provide yet another derivation leading to this matrix model.

Instead of $(\phi(X))^3$ interaction in the potential in (3.7) we could take any polynomial potential. Then, instead of trivalent graphs, we will study the same problem of spinless massive fermions (or forests) on the corresponding collection of planar graphs. For the generic couplings of this potential the critical behavior we study below should not change (due to universality). However, we can have multicritical points of the kind known from [12, 32].
## 3.2 Combinatorial derivation of the matrix model

Here we will present another, combinatoric, derivation of the matrix model (3.7) describing our partition function.

Let us consider the graphical interpretation for our action (3.7) and demonstrate that it indeed reproduces the expansion (2.1). First, recall that the function (3.8) appearing in our matrix model potential

$$ \phi(X) = \frac{(1 - \sqrt{1 - 4\lambda X})}{2\lambda} = X + \lambda X^2 + 2\lambda^2 X^3 + 5\lambda^3 X^4 + \dots = \sum_{n=1}^{\infty} \frac{1}{n+1} \binom{2n}{n} \lambda^n X^{n+1} \quad (3.9) $$

is the generating function of rooted trees with node degree $q = 3$. Here the expansion coefficients $C_n = \frac{1}{n+1}\binom{2n}{n}$ are nothing but the Catalan numbers. Each vertex is weighted with the coupling $\lambda$, and the exterior branches are decorated with the matrix $X$ (‘Christmas Tree’). In our potential we have $\phi^3(X)$ which generates three such trees growing from the same vertex (call it ‘root vertex’).

Next, we expand the exponent in (3.7) under the matrix integral w.r.t. the 2nd term in the action. This expansion generates forests of such Christmas trees with marked vertex. Each marked vertex (corresponding to $\phi^3(X)$) is weighted with $\lambda m^2$, while the other vertices are weighted with $\lambda$. Hence each such forest has weight $\lambda^{|G|}m^{2l}$, where $|G|$ and $l$ denote the total number of vertices and the number of trees respectively.

As a next step, we see that the Gaussian integral over $X$ connects, via rainbow diagrams, these Christmas trees into graphs, so that we sum up over all “forests” of such trees on each graph. All graphs are planar, due to the matrix structure. We notice that this statistical-mechanical system, on each graph, is the same as given by the equation (13) of [53], where it is proven to be equal to the determinant $\det[m^2 + \Delta(G)]$ on each graph $G$. This is another confirmation of our derivation of the matrix model (3.7) in section 3.1 and its relation to the partition function (2.2).

To prove that the combinatorial factors, such as the Catalan numbers and symmetry factors of Feynman graphs in (3.7), do indeed match up with the Kirchoff theorem (the proof of its “massless” version is given in appendix A) and the expansion (2.2) we use the result of the paper [54] which considers a model that generates forests of “unrooted trees”, i.e trees without a marked point. Instead of (3.7) it corresponds to a similar potential:

$$ S_X = \mathrm{tr} \left[ -\frac{1}{2} X^2 + m^2 X^2 \tilde{V}(\lambda X) \right] \quad (3.10) $$

where

$$ \begin{aligned} \tilde{V}(z) &= \frac{1}{12z^2} (-6z^2 + (1 - 4z)^{3/2} + 6z - 1) \\ &= \sum_{n=1}^{\infty} \tilde{c}_n z^{n+2}, \quad \text{where } \tilde{c}_n = \frac{(2n)!}{n!(n+2)!}. \end{aligned} \quad (3.11) $$

The coefficients of the expansion are equal to the Catalan numbers $C_n = \frac{(2n)!}{n!(n+1)!}$ divided by $n+2$ giving the number of unrooted trees with $n$ vertices with cyclic symmetry factored
out. Notice now the following identity relating our potential with $\tilde{V}(z)$:

$$ \frac{\lambda}{3} \left( \frac{1 - \sqrt{1 - 4\lambda X}}{2\lambda} \right)^3 = \lambda \partial_{\lambda} (X^2 \tilde{V}(\lambda X)). \qquad (3.12) $$

In terms of coefficients of expansion for the l.h.s. this relation looks as follows:

$$ c_n = n \frac{C_n}{(n+2)}. \qquad (3.13) $$

In other words, the number of trees with $n$ marked vertices is equal to the number of rooted trees divided by the cyclic permutation order $n + 2$ of external legs (called “leafs” in [54]), which makes these trees unrooted, and multiplied by the $n$ — the number of ways to mark one vertex. This is another proof that our model (3.7) describes indeed the sum over forests of rooted trees over planar graphs, which is equal, according to the above-mentioned theorem from [53] (see equation (13) there), to the partition function of massive spinless fermions (our $\psi_i$'s in (2.1)) on planar graphs. In our notation, for each individual planar graph this theorem states

$$ \det[m^2 + \Delta(G)] = \sum_{k=1}^{\infty} m^{2k} F_k(G) \qquad (3.14) $$

where $F_k(G)$ is the number of forests consisting of $k$ unrooted trees, each with one marked vertex, on this graph $G$.

### 3.3 Boundary conditions for disc partition functions

As announced in the Introduction, we are going to be interested in computing disc partition functions. There are two types of disc partition functions we can compute in the matrix model (3.7). One corresponds to a boundary with the $X$-matrix, and another to the matrix $\phi(X)$ that generates trees. They correspond to two resolvents we will study: the resolvent $G(x)$ which generates $\langle \text{tr} X^k \rangle$ correlators,

$$ G(x) = \sum_{k=0}^{\infty} \frac{1}{x^{k+1}} \frac{1}{N} \langle \text{tr} X^k \rangle \qquad (3.15) $$

and the resolvent for $\langle \text{tr} \phi^k \rangle$ correlators,

$$ H(r) = \sum_{k=0}^{\infty} \frac{1}{r^{k+1}} \frac{1}{N} \langle \text{tr} \phi^k \rangle \qquad (3.16) $$

with $k$ in the correlator being the disc boundary length, while $\log x$ and $\log r$ play the role of the bare boundary cosmological constant.

One can expect that these two resolvents correspond to two different types of boundary conditions. We would like to formulate these boundary conditions in terms of the worldsheet action (1.1) for the anti-commuting bosons. To do this we first recall how the trees appear from the microscopic model on graphs (2.1).
Let us fix a graph G and consider the action for this graph. To produce trees we expand the kinetic terms from (2.1). Keeping in mind that our variables are anti-commuting we obtain:

$$ \prod_{<i,j>} e^{|\psi_{ij}|^2} = \prod_{<i,j>} (1 + |\psi_{ij}|^2). \qquad (3.17) $$

Expanding the product over edges for each edge we can either choose to include it into the tree, which corresponds to choosing the $|\psi_{ij}|^2$ term, or not, in which case we take the identity. Loops are forbidden because of the vanishing of the cyclic product $\prod_{<i,j> \in L} |\psi_{ij}|^2$.

Hence we obtain the expression for trees T on the graph:

$$ \prod_{<i,j>} e^{|\psi_{ij}|^2} = \sum_{T \in G} V(T) \prod_{<i,j> \in T} |\psi_{ij}|^2 \qquad (3.18) $$

where $V(T)$ is the number of vertices in a tree. However in order to obtain a non-zero expression we should also introduce boundaries for the trees. The first type of boundary corresponds to insertions of Tr $\phi^n$. In graph terms this means that we have a graph with trees that start from the boundary and some trees in the bulk. Trees in the bulk are generated by the mass terms, while boundary trees should be enforced by fixing the field $\phi$ on the boundary, as the tree-generating formula (3.8) suggests. This is done by inserting a specific boundary term under the integral:

$$ \left( \prod_{<i,j>} e^{|\psi_{ij}|^2} \prod_{<b,j>} e^{|\psi_{bj}|^2} \right) \psi_b \bar{\psi}_b \qquad (3.19) $$

where b denotes the boundary. In other words we fix the value of $\psi$ on the boundary.

The other type of boundary corresponds to operators Tr $X^n$. Now we have a boundary with edges going into the bulk, while trees appear already in the bulk and do not touch the boundary. One can describe such behaviour on the worldsheet by an insertion

$$ \prod_{<i,j>} e^{|\psi_{ij}|^2} \prod_{<b,j>} e^{|\psi_{bj}|^2} \prod_{<b,j>} |\psi_{bj}|^2 \qquad (3.20) $$

effectively fixing the derivative of the field on the boundary.

Having described the discrete picture it is natural to conjecture that the proper continuous boundary conditions are Dirichlet conditions

$$ \psi(\partial D) = 0 \qquad (3.21) $$

for the $\phi$-type boundary disc partition function $H(r)$ and Neumann conditions

$$ \partial_{\perp} \psi|_{\partial D} = 0 \qquad (3.22) $$

for the X-type boundary appearing in $G(x)$.

The critical behaviour of both disc partition functions, related to the continuous limit of 2d QG interacting with massive fermions, appears to be different. We will compute it in section 7.
<figure><img src="image_2.png" /><figcaption>Figure 1. The matrix model potential (4.1) at $\lambda = 6/100$, $M = 1/15$.</figcaption></figure>

## 4 Solution of the matrix model in the planar limit: saddle point

As discussed above, the problem we study reduces to a 1-matrix model with the potential

$$V(x) = \frac{1}{2}x^2 - \frac{3M}{16\lambda^2}(1 - \sqrt{1 - 4\lambda x})^3 \quad (4.1)$$

where we make a redefinition of the mass parameter for convenience:

$$M = \frac{2}{9}m^2 \quad (4.2)$$

In this section we describe how to solve this matrix model at large $N$ and present the 1-cut solution in detail.

We give a plot of the potential on figure 1. One can check that the shape remains the same regardless of the values of $M$ and $\lambda$ (with $M, \lambda > 0$) — namely, as we go from negative $x$ the potential has a local minimum at $x = 0$ and then a local maximum before falling off to some finite value at $x = 1/(4\lambda)$ which is the boundary of the allowed region due to the presence of the square root in the potential.

Since the potential has a local minimum at $x = 0$, it can be populated by eigenvalues and thus it is natural to look for a 1-cut solution which can be found by standard methods and which we will focus on in this paper. We parameterize the endpoints of the cut by $b$ and $a$ with $b < a$ and also denote the branch point of the square root by<fnref n="9" />

$$c = \frac{1}{4\lambda} \quad (4.3)$$

so that $b < a < c$ are the three branch points in our problem (with one more at infinity coming from the square root) meaning it is resolvable in elliptic functions. Introducing the density of eigenvalues $\rho(x)$ normalised to 1 on the interval $[b, a]$ and the resolvent

$$G(x) = \int_b^a \frac{dy \rho(y)}{x - y} \quad (4.4)$$

<footnote n="9">⁹We hope this notation will not create confusion with the notation $c$ for the central charge of the model.</footnote>
we can write the saddle point equation on the $[b, a]$ cut as

$$G(x + i0) + G(x - i0) = x + 9cM + \frac{9}{2} M \sqrt{c} \frac{x - 2c}{\sqrt{c - x}}. \quad (4.5)$$

The solution can be written easily as

$$G(x) = \sqrt{x-a}\sqrt{x-b} \int_{b}^{a} \frac{dy}{2\pi(x-y)} \frac{1}{\sqrt{(a-y)(y-b)}} \left( y + 9cM + \frac{9}{2} M \sqrt{c} \frac{y-2c}{\sqrt{c-y}} \right). \quad (4.6)$$

The integrals here can be taken explicitly (see appendix C) and we find the resolvent,

$$G(x) = \frac{9}{4\pi} \sqrt{cM} \sqrt{x-a} \sqrt{x-b} \left( \frac{2(x-2c)\Pi\left(\frac{a-b}{x-b}\right|1-p)}{\sqrt{c-b}(x-b)} - \frac{2K}{\sqrt{c-b}} \right) + \frac{1}{2} \left( -\sqrt{x-a}\sqrt{x-b} + 9cM + x \right) \quad (4.7)$$

where we denoted

$$p = \frac{c-a}{c-b} \quad (4.8)$$

and $K \equiv K'(p) \equiv K(1-p)$ is the elliptic integral with modulus$^{10}$ $1-p$.

Then we can find the density as the discontinuity of the resolvent on the $[b, a]$ cut,

$$G(x \pm i\epsilon) - G(x - i\epsilon) = -2i\pi\rho(x) \quad (4.9)$$

which gives the density in terms of elliptic functions as well,$^{11}$

$$\rho(x) = \sqrt{(a-x)(x-b)} \left( \frac{9\sqrt{cMK}}{2\pi^2\sqrt{c-b}} + \frac{1}{2\pi} \right) - \frac{9\sqrt{cM}(x-2c)}{2\pi^2\sqrt{c-b}} \sqrt{\frac{a-x}{x-b}} \left( K - \Pi \left( \frac{x-b}{c-b}, 1-p \right) \right). \quad (4.10)$$

## 4.1 Fixing parameters

We fix $a, b$ from the condition that for $x \to \infty$ we must have $G(x) \simeq 0 \cdot x^0 + \frac{1}{x} + O(x^{-2})$ which follows from the definition of the resolvent and the normalisation of the density. Introducing the variable

$$y = \frac{\sqrt{c-b}}{\sqrt{c}} \quad (4.11)$$

we get two relations fixing $y$ and $p$ (or equivalently the branch points $a$ and $b$) in terms of our original couplings $\lambda$ and $M$, namely the $x^0$ term in $G(x)$ gives

$$-\pi(p+1)y^3 - 18MEy^2 + 2\pi(9M+1)y - 18MK = 0 \quad (4.12)$$

<footnote>$^{10}$We denote the modulus by $1-p$ for convenience, so that $p \to 0$ corresponds to the modulus approaching 1 which is the limit we will study later. In Mathematica the functions we use correspond to EllipticK[1-p], EllipticPi[n,1-p], etc.</footnote>

<footnote>$^{11}$In terms of the expression in (4.10), the density we should integrate to get $G$ is to be taken as $\rho(x + i\epsilon)$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The density $\rho(x)$ given by (4.10) on the $[b, a]$ cut for $\lambda = 5/106 \approx 0.047$, $M = 1156/25 \approx 46$. We chose the small value of $\lambda$ and large $M$ so that the density has a nontrivial profile substantially different from a semicircle.</figcaption></figure>

relating $p$ and $y$, while from the $1/x$ term we have

$$ -\lambda^2 + \frac{(1-p)^2}{256} y^4 - \frac{3M}{64\pi}(-(p+1)E+2pK)y^3 + \frac{9M}{64\pi}(-(p+1)K+2E)y = 0. \quad (4.13) $$

Here we denoted the elliptic integral of the 2nd kind as $E \equiv E'(p) \equiv E(1-p)$. Notice that $\lambda$ only enters the second equation and is explicitly expressed as a function of $y, p$.

As these are two polynomial equations in $y$, we can exclude $y$ by taking their resultant<fnref n="12" /> which we denote as $P(\lambda, M, p)$, so we have

$$ P(\lambda, M, p) = 0. \quad (4.14) $$

This gives a lengthy, though explicit, equation linking $\lambda, M$ and $p$. The function $P$ is a polynomial of 6th order in $M$ and of 3rd order in $\lambda^2$, thus one can in principle write $\lambda(M, p)$ explicitly.

These equations are a complete system that allows one to obtain the 1-cut solution at a given value of our parameters $\lambda$ and $M$. In order to e.g. solve the system numerically, we can first solve (4.14) for $p$ and then plug the result into (4.13) which is then solved for $y$. As an example we give a plot of the density in figure 2. Note that in general these equations have multiple solutions and we should be careful to select the ‘physical’ one, for which the density is real and positive. We discuss this requirement in more detail in section 5.

It is instructive to plot the effective potential $V_{\text{eff}}$ felt by the eigenvalues, which as usual is given by

$$ V_{\text{eff}}(x) = V(x) - 2 \int_{0}^{x} dy (G(y - i0) + G(y + i0)) \quad (4.15) $$

From this definition it follows that it has a flat section on the $[b, a]$ cut where it vanishes. We show a plot of it on figure 3. Apart from the flat region, its shape is similar to the original matrix model potential.

<footnote n="12">¹²We recall that for two polynomials with roots $a_i$ and $b_j$ the resultant is defined as $\prod_{i,j} (a_i - b_j)$. As a symmetric function in both sets of roots, it is a polynomial in the coefficients of the original two polynomials. It is also implemented in Mathematica.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. The effective potential $V_{\text{eff}}(x)$ at $\lambda = 1/18$, $M = 16/100$.</figcaption></figure>

## 4.2 Exact results for correlators

Let us discuss how to compute the observables in this model. First, the correlators $\langle \frac{1}{N} \text{tr} X^k \rangle$ are encoded in the large $x$ expansion of the resolvent $G(x)$ and can be obtained directly by expanding (4.7). Another important set of observables are built from the matrix $\phi$ related to $x$ by (3.8). They can be written as

$$ \langle \frac{1}{N} \text{tr} \phi^k \rangle = \int_b^a dx \, \rho(x) \left( 2c \left( 1 - \sqrt{1 - \frac{x}{c}} \right) \right)^k, \quad k = 1, 2, \dots \quad (4.16) $$

These integrals are highly nontrivial since the density involves the elliptic $\Pi$ function. However, remarkably, we found a way to compute them in closed form for any given $k$. The tricky part is the integration of the square root with the elliptic integral of third kind in (4.10). To explain how this integral can be evaluated explicitly, first denote

$$ n = \frac{b-x}{b-c}. \quad (4.17) $$

The integral is then over $n$ between 0 and 1. The relevant piece is:

$$ \int_0^1 dn (n-1)^{k/2} \sqrt{\frac{1-p-n}{n}} (bn - b - cn + 2c) \Pi(n|1-p) \quad (4.18) $$

The key trick is to get rid of the integration of elliptic functions, by using the identities presented in appendix A (see section C.2 for details). The integration of $\Pi(n|1-p)$ is reduced to an integral of a derivative, while the remaining part contains only algebraic dependence on $n$ and can be easily integrated.

As an example, important for the further computations of this one point function in the critical regime, we give below the explicit result for $k=1$

$$ \langle \frac{1}{N} \text{tr} \phi \rangle = \frac{c^3 y}{2\pi^2} \left[ \frac{4}{15}\pi p(p+1)y^4 K - \frac{8}{15}\pi (p^2 - p + 1)y^4 E + \frac{1}{4}\pi^2(p-1)^2 y^3 \right. \\ \left. + M \left( K \left( 12y \left( py^2 + p + 1 \right) E - 3\pi \left( p \left( 2y^2 + 3 \right) + 3 \right) \right) \right. \\ \left. - 6pyK^2 - 6y \left( (p+1)y^2 + 3 \right) E^2 + 3\pi \left( (p+1)y^2 + 6 \right) E \right) \right]. $$
<figure><img src="image_5.png" /><figcaption>Figure 4. The critical curve $\lambda_c(M)$. The shaded area below the curve shows the allowed physical region on the $(M, \lambda)$ plane where the density is real and positive.</figcaption></figure>

We also give the result for $\langle \frac{1}{N} \text{tr} \phi^3 \rangle$, which is more lengthy, in a supplementary Mathematica notebook accompanying this paper. This observable is particularly important as it gives the derivative of the partition function w.r.t. the mass $m^2$, which in view of (2.1) means it describes the fermionic condensate of the type $\langle \bar{\psi} \psi \rangle$.

It would be also interesting to obtain an explicit result similar to (4.7) for the generating function of these correlators, which seems to be quite challenging but might be possible to do by using the methods of [33].

## 5 The critical line

The critical regime in our model, corresponding to large graphs, is obtained by going close to a singularity of the partition function. A (standard) shortcut to deriving the condition for criticality is requiring that the density has zero slope at the endpoint, $\rho'(a) = 0$. The reason for this is that generically we have $\rho(x) \sim \sqrt{a - \bar{x}}$ when $x \to a$, and since the density should be positive the coefficient in front of the square root must be nonnegative as well. The case when the coefficient becomes zero thus belongs to the boundary of the allowed parameter space where the model becomes singular.

From our explicit result for the density (4.10) we find that imposing $\rho'(a) = 0$ gives

$$ -\pi(1-p)py^3 - 9Mp(K-E)y^2 + 9M(-pK + E) = 0. \qquad (5.1) $$

Combining this with (4.12), (4.13) have three conditions linking the four variables $M, \lambda, p, y$, and thus we get a line on the $(M, \lambda)$ plane that we will call the critical line. We denote it by $\lambda_c(M)$ and we give a plot of it on figure 4.
In order to study the critical regime analytically we can take the resultant in $y$ of (4.12) and (5.1) which gives a 4th order polynomial in $M$ equation involving only $M$ and $p$, whose solution<fnref n="13" /> is $M_c(p)$ which can be written explicitly but is rather lengthy. Furthermore, we can solve (4.12) as a linear equation for $M$, then plugging the result into (4.13) and (5.1) and taking their resultant in $y$ gives a 4th order polynomial equation whose solution is $\lambda_c(p)$. Thus we have a parametric representation of the critical line $\lambda_c(p)$ in terms of two (explicit but very lengthy) functions $M_c(p)$ and $\lambda_c(p)$.

## 5.1 Merging of solutions and the physical region

Numerically we observe that when $\lambda < \lambda_c(M)$ our two constraint equations (4.12), (4.13) have two real solutions for $p, y$, with one of them being the actual physical solution while the other one should be discarded as it corresponds to a non-positive density. As we move close to the critical line, we find that these solutions get closer and finally merge at the value $\lambda = \lambda_c(M)$. Beyond that point, i.e. for $\lambda > \lambda_c(M)$, there is no solution with real positive density.<fnref n="14" /> Note that naturally the origin $(\lambda, M) = (0, 0)$ lies in the allowed region since $\lambda$ and $M$ are weights in our partition function which is well defined for small enough values of them.

As a technical consistency check, let us show how to derive $M_c(p)$ independently starting from our two original constraints (4.12) and (4.13). For that we consider the equation (4.14) which follows from them and reads $P(\lambda, M, p) = 0$. Considering its l.h.s. as a function of $p$ at fixed $\lambda, M$, we find that generically it has three roots $0 < p_1 < p_2 < p_3 < 1$ of which $p_2$ is the physical one. It merges with $p_1$ at some $\lambda = \lambda(M)$ which will be the critical value and is characterized by the condition that $\partial_p P(\lambda, M, p) = 0$. Since this condition and the original equation $P(\lambda, M, p) = 0$ are both polynomials in $\lambda$, we can take their resultant and find a polynomial equation now for $M(p)$, which is solved by the same function $M_c(p)$ we found using the shortcut $\rho'(a) = 0$, i.e (5.1). This shows that indeed (5.1) corresponds to the boundary of the allowed parameter region, as expected. Notice also that the condition $\partial_p P(\lambda, M, p) = 0$ we just discussed means that on the critical line

$$ \partial_p \lambda(M, p) = 0 \qquad (5.2) $$

which is<fnref n="15" /> another equivalent formulation of the criticality condition.

## 5.2 Looking for further singularities

It could happen that on the critical curve itself there are special points at which additional phase transitions take place and the structure of the solution changes. So far we have not found any features of this type. For example, since $M_c(p)$ is given by a solution to a polynomial equation, one source of singularities could be the crossing of its roots at particular values of $p$. Indeed we can identify several such values by looking at the discriminant of this equation and numerically they are $p \simeq 0.05, p \simeq 0.29$. We find numerically that they

<footnote n="13">¹³Note that we also have to pick the correct branch.</footnote>

<footnote n="14">¹⁴There could be a two-cut solution in that region, whose exploration we postpone to the future.</footnote>

<footnote n="15">¹⁵Notice that here we first differentiate $\lambda$ as a function of $M$ and $p$, and only then set $M = M_c(p)$.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5. The derivative $\partial_{pp}\lambda(M, p)$ evaluated on the critical line, shown as a function of $p$. We see no zeros or singularities for $0 < p < 1$.</figcaption></figure>

do not seem to correspond to a singularity of the physical solution but rather the crossing of singularities of unphysical roots. Similarly, for $\lambda_c(p)$ the discriminant of the equation that it solves vanishes at $p \simeq 0.09$ but again this does not seem to give a singularity in the physical solution.

To further make sure that we are not missing any singularities, we can plot the 2nd derivative $\partial_{pp}\lambda(p, M)$ at fixed $M$, evaluated on the critical line $M = M_c(p)$ (recall that the first derivative $\partial_p\lambda(p, M)$ vanishes on the critical line). We give its plot on figure 5. We find that it does not have any zeros or singularities and thus we would not expect to find any special points on the critical curve. A “physical” picture of eigenvalues behaving as Coulomb charges confined in a potential well also does not suggest the existence of any additional phase transition for finite $\lambda$ and $M$ along the critical line: the appropriate changes in these parameters only smoothly change the density of eigenvalues. The critical line corresponds to the values of $M_c(\lambda)$ when the eigenvalues start spilling over the top of the potential (see again figure 1). We leave a more careful analysis for the future.

## 5.3 Decoupling of the heavy matter: $M \to \infty$ limit

Let us discuss the limit of large mass when we expect the fermions to decouple so that we are left with the pure gravity limit. We see from the matrix model potential that when $\lambda \to 0$ the expansion of the square root yields the cubic term in the potential required for pure gravity. To get finite coupling in pure gravity we have to take a scaling limit $\lambda \to 0$, $M \to \infty$ with $M\lambda = \lambda_{\text{eff}}$ kept finite:

$$V(z) = \frac{1}{2}z^2 - \frac{3}{2}\lambda_{\text{eff}}z^3 + O(\lambda^2). \qquad (5.3)$$

Let us identify this limit in our critical curve. The regime $M \to \infty$, $\lambda \to 0$ corresponds to $p \to 1$ and from the numerical solution for the critical line we see that the scaling is
$M \sim \frac{1}{1-p}$, $\lambda \sim 1-p$, $y \sim 1$. Then from equations (4.12), (4.13), (5.1) we find

$$
\begin{align}
\lambda_c &= \frac{1-p}{16\sqrt[4]{3}} + \frac{(31\sqrt{3}-18)(1-p)^2}{768\ 3^{3/4}} + \dots, \tag{5.4} \\
M_c &= \frac{16}{9\sqrt{3}} \frac{1}{1-p} + \left(-\frac{4}{9} - \frac{28}{27\sqrt{3}}\right) + \dots. \nonumber
\end{align}
$$

Therefore the critical value for the effective coupling gives:

$$
\left( \frac{3\lambda_{\text{eff,crit}}}{2} \right)^2 = \frac{1}{108\sqrt{3}} \qquad (5.5)
$$

which is in perfect agreement with the well known results from [3] and [21]. This is a
nontrivial test of our calculation.

The emergence of the proper scaling parameter can be seen in the initial combinatorial
partition function as well

$$
Z(\lambda, m) = \sum_G \lambda^{|G|} \det[m^2 + \Delta(G)] = \sum_G \lambda^{|G|} \sum_{F=(F_1 \dots F_l) \in G} \prod_{i=1}^l m^2 V(F_i). \quad (5.6)
$$

Consider the $m \to \infty$ limit (we recall that $M = 2/9m^2$) and assume that the contribution from the term when all trees become the isolated nodes dominates. In this case $V(F_i) = 1$ and we immediately get the $(\lambda m^2)^{|G|}$ weight factor which means the logarithmic renormalization of the bare cosmological constant $\Lambda = -\log \lambda$ obtained from the matrix model.

5.4 Critical line at $M \to 0$

Having discussed the large *M* limit in the previous subsection, in the remaining part of the
paper we will mostly focus on the opposite limit of small *M*. Since *M* controls the number
of trees, for *M* = 0 we have only 1 tree in the forest and thus for each graph we are counting
the number of spanning trees on it with *λ* weighing the number of graph’s vertices. This
is a well studied problem, describing in the critical regime the 2d QG in the presence of
*c* = −2 matter, discussed in [3, 30, 47–50, 71] and the critical value is known to be *λ*<sub>*c*</sub> = 1/8.
Indeed we verified this is what we find from our critical curve (as can be seen already on
the plot in figure 4). This is another successful test of our calculation. In subsequent parts
of this paper we will make further contact with known results in this regime.

Let us also present here the expansion of the critical line near this point, which will be
important for our further calculations. For the critical line $M \to 0$ corresponds to $p \to 0$
and it is convenient to use $p$ as an expansion parameter. We will also use instead of $\lambda$ a
redefined coupling

$$
g = 256\pi\lambda^2. \tag{5.7}
$$

Then it is straightforward to obtain the expansion of the critical values $M_c(p)$ and $g_c(p)$.
For $M_c(p)$ we have

$$
M_c(p) = Ap + [B + C \log(p)] p^2 + [D + E \log(p) + F \log^2(p)] p^3 + \mathcal{O}(p^4 \log^3(p)) \quad (5.8)
$$
where the coefficients read

$$A = \frac{2\sqrt{2\pi}}{9}, \quad B = \frac{1}{18}\pi(24\pi - \sqrt{2}(41 + \log(16))), \quad C = \frac{\pi}{9\sqrt{2}}, \quad (5.9)$$

$$D = \frac{1}{288}\pi(64\pi(21\pi\sqrt{2} - 155 + \log(16)) + \sqrt{2}(8385 + 8\log(2)(173 - 168\log(2)))},$$

$$E = -\frac{1}{144}\pi(32\pi + \sqrt{2}(173 - 336\log(2))), \quad F = -\frac{7\pi}{12\sqrt{2}}.$$

We see that each power of $p$ is accompanied by an expansion of a growing number of powers of $\log p$. For the $g$ coupling we find

$$g_c(p) = A + [B + C \log(p)] p + [D + E \log(p) + F \log^2(p)] p^2 + \mathcal{O}(p^3 \log^3(p)) \quad (5.10)$$

with

$$A = 4\pi, \quad B = \frac{16}{3}(3\pi^2\sqrt{2} - \pi - 12\pi \log(2)), \quad C = 16\pi, \quad (5.11)$$

$$D = \frac{4}{3}\pi(96\pi^2 - 171\sqrt{2}\pi + 38 + 96\log^2(2) - 204\sqrt{2}\pi \log(2) + 568\log(2)),$$

$$E = \frac{4}{3}\pi(51\pi\sqrt{2} - 142 - 48\log(2)), \quad F = 8\pi.$$

Notice that the expansion starts with $g = 4\pi + \dots$, corresponding of course to the critical value $\lambda_c = 1/8$ discussed above.

We can also invert these expansions and find $g_c(M)$ which reads to leading order

$$g_c = 4\pi + (36\sqrt{2}L_M + 72\pi - 12\sqrt{2})M + O(M^2) \quad (5.12)$$

where

$$L_M = \log \frac{9M}{32\sqrt{2}\pi}. \quad (5.13)$$

This gives the shape of the critical curve near its $M = 0$ endpoint.

## 6 New double scaling limit for 1-pt functions

At a generic point on the critical curve we expect the continuum theory to be pure gravity. However, when $M$ is strictly zero we have the very different $c = -2$ theory. This suggests to explore a double scaling limit when we get close to the critical line, in the vicinity of its $M = 0$ endpoint. In this section we will define and study this near-critical regime in which we expect to see a nontrivial interpolation (flow) between the $c = 0$ and $c = -2$ theories.

In order to design an interesting scaling limit, let us look at the 1-pt function $\langle \frac{1}{N} \text{tr} \phi \rangle$ that we computed in closed form in (4.16). Technically it is convenient to use $p$ instead of $\lambda$ as an expansion parameter, as otherwise we would need to deal with iterated log log corrections. The endpoint of the critical line corresponds to $(M, p) = (0, 0)$, and we will consider an expansion near it with both $M$ and $p$ being small. As a first example, the
expansion of $\lambda$ or equivalently $g$ to the first few orders has the form

$$g = 4\pi - 16\pi p + 32\pi p^2 + O(p^3) \quad (6.1)$$

$$+ M \left[ \frac{(576 (\sqrt{2} - \pi) \pi - 54\sqrt{2}\pi L + O(L^2)) p}{2\pi} + \frac{48\pi\sqrt{2} + 144\pi^2 + 72\pi\sqrt{2}L + O(L^2)}{2\pi} + O(p^2) \right] + O(M^2) \quad (6.2)$$

with

$$L = \log \frac{p}{16}. \quad (6.3)$$

Thus at each order in $M$ we have a series in positive powers of $p$ and $\log p$. As a technical intermediate result, we also give the expansions of $a, b, c$ to higher order in appendix D.

Next, expanding (4.16) we find for the 1-pt function $\frac{1}{N}\langle \text{tr } \phi \rangle$

$$\frac{1}{N}\langle \text{tr } \phi \rangle = 4 - \frac{128\sqrt{2}}{15\pi} + \left(8 - \frac{96\sqrt{2}}{5\pi}\right)p \quad (6.4)$$

$$+ M \left[ \frac{12(124 + 11\sqrt{2}\pi - 15\pi^2)}{5\pi^2} + \frac{(432 - 90\sqrt{2}\pi)L}{5\pi^2} + p \left( -\frac{24(3\pi(9\sqrt{2} + 5\pi) - 464)}{5\pi^2} + \frac{(6192 - 945\sqrt{2}\pi)L}{10\pi^2} \right) \right] + \dots$$

where like in (6.1) we dropped various higher order terms.

We see in both expansions (6.1) and (6.4) a similar structure of a double series in $p$ and $M$ (with additional $\log p$ terms). This suggests to consider the limit when $p$ and $M$ both go to zero with their ratio fixed. Thus we define

$$z = \frac{p}{16M} = \text{finite}, \quad M \sim p \to 0 \quad (6.5)$$

and eliminate $p$ in favor of $z$. This gives a series now only in $M$, in which any given order receives contributions only from a finite number of terms in the original expansion and is thus straightforward to compute. Geometrically in this limit the branch points at $a$ and $c$ collide as they both approach 2 with their difference being of order $p$, while the $[b, a]$ cut remains of a finite size as $b \simeq -2$.

Below we will compute in this limit the 1-pt functions $\frac{1}{N}\langle \text{tr } \phi \rangle$ and $\frac{1}{N}\langle \text{tr } \phi^3 \rangle$. We will find that the results for them are closely related.
## 6.1 The $\frac{1}{N}\langle \text{tr } \phi \rangle$ one-point function

Let us first discuss the behavior of the 1pt function $\frac{1}{N}\langle \text{tr } \phi \rangle$ in this scaling limit. We will need the result for the 1pt function to quadratic order only, and it reads

$$ \frac{1}{N} \langle \text{tr } \phi \rangle = - \frac{4(32\sqrt{2}-15\pi)}{15\pi} \qquad (6.6) \\ + \frac{2M}{5\pi^2} \left( (-45\sqrt{2}\pi + 216)\log(Mz) + 320\pi^2 z - 768\sqrt{2}\pi z - 90\pi^2 + 66\pi\sqrt{2} + 744 \right) \\ + \frac{M^2}{20\pi^3} \left[ \begin{aligned} & -20480\sqrt{2}\pi^2 z^2 \log(Mz) - 13608\sqrt{2}\log^2(Mz) + 4455\pi \log^2(Mz) \\ & - 30240\sqrt{2}\pi^2 z \log(Mz) + 198144\pi z \log(Mz) - 73008\sqrt{2}\log(Mz) \\ & - 23004\pi \log(Mz) + 8100\pi^2\sqrt{2}\log(Mz) + 40960\pi^3 z^2 - 182272\sqrt{2}\pi^2 z^2 \\ & - 23040\pi^3 z - 41472\sqrt{2}\pi^2 z + 712704\pi z \\ & + 6480\pi^3 + 4536\pi^2\sqrt{2} - 74880\sqrt{2} - 108432\pi \end{aligned} \right]. $$

Let us also define instead of $g$ the rescaled and shifted cosmological constant

$$ J = 2\pi \frac{4\pi - g}{M} - (648\sqrt{2}\pi - 648\pi^2) M + 144\pi^2 + 48\sqrt{2}\pi \qquad (6.7) $$

which in the limit $M \to 0$ reads

$$ J = 512\pi^2 z - 72\sqrt{2}\pi \log(Mz) \\ + M \left[ \begin{aligned} & -81 \log^2(Mz) + 864\pi\sqrt{2}z \log(Mz) - 324(\sqrt{2}\pi - 1) \log(Mz) \\ & - 16384\pi^2 z^2 + 9216\pi(\pi - \sqrt{2})z \end{aligned} \right]. \qquad (6.8) $$

The first term $2\pi \frac{4\pi-g}{M}$ in the r.h.s. of (6.7) is a finite nontrivial quantity in the limit we consider (as follows from (6.1)) and can be viewed as a renormalized coupling in our regime. The other terms in that equation serve to remove from (6.8) the trivial nonsingular pieces that are polynomial in $M$ and do not contain $\log(M)$ or $z$, and thus will not affect the part of the 1-pt function we are interested in.

Similarly, subtracting from (6.6) the irrelevant regular terms we define its singular part as

$$ \frac{1}{N} \langle \text{tr } \phi \rangle_{\text{sing}} = \left\langle \frac{1}{N} \text{tr } \phi \right\rangle + \frac{4(32\sqrt{2}-15\pi)}{15\pi} - \frac{(5\sqrt{2}\pi-24)MJ}{20\sqrt{2}\pi^3} \qquad (6.9) \\ - \frac{(1134\pi^2\sqrt{2}-18720\sqrt{2}+1620\pi^3-27108\pi)M^2}{5\pi^3} \\ - \frac{12(11\pi\sqrt{2}-15\pi^2+124)M}{5\pi^2} $$
where the coefficient of the $MJ$ term is chosen so as to cancel the $M \log(Mz)$ term in (6.6). The final result reads

$$ \langle \frac{1}{N} \mathrm{tr} \phi \rangle_{\mathrm{sing}} = \frac{1}{\pi^3} M^2 \left[ \begin{aligned} & (486\pi^2\sqrt{2} - 3456\sqrt{2} - 1620\pi) \log(Mz) - 243(3\sqrt{2} - \pi)\log^2(Mz) \\ & + z \left( (10944\pi - 1728\sqrt{2}\pi^2)\log(Mz) + (5760\pi^2\sqrt{2} - 3456\pi^3 + 24576\pi) \right) \\ & + z^2 \left( (6144\pi^3 - 18944\sqrt{2}\pi^2) - 1024\sqrt{2}\pi^2\log(Mz) \right) \end{aligned} \right]. \quad (6.10) $$

Notice also that we can now retain in the definition (6.8) of $J$ only the leading terms since it is already given by 2nd order in $M$, and we denote the resulting quantity by little $j$,

$$ j = 512\pi^2 z - 72\sqrt{2\pi} \log(Mz). \quad (6.11) $$

To summarise, the singular part of the 1-pt function is given by (6.10) where $z$ is implicitly a function of the coupling $g$ determined via the intermediate variable $j$ in (6.11) (or equivalently $J$), related to the coupling via (6.7). These equations are one of our main results and give the novel 1pt function representing the flow between $c = -2$ and $c = 0$ regimes.

Notice that according to (6.7) the variable $j$ is essentially the coupling $g$, up to constant factors, a shift and a rescaling by $M$. To write the 1-pt function in terms of $M$ and $g$ one would need to invert (6.11) to obtain $z(j)$. Curiously, this inverse function is (up to numerical constants that can be scaled away) in fact the well known Lambert function which has a variety of interpretations from combinatorics to quantum field theory and which appears here in our new result for the 1-pt function. We discuss its role and origins in more detail in section E.

### 6.1.1 Limiting regimes

Let us now discuss the two limiting cases in more detail. In the $c = -2$ regime describing spanning trees we have to send $z \to \infty$ as we take $M$ to be small and $4\pi - g \gg M$. This gives $z \simeq 2\pi \frac{4\pi - g}{512\pi^2 M}$. Plugging it into (6.10) and retaining there the leading term (last line) we get the correct $c = -2$ scaling

$$ \langle \frac{1}{N} \mathrm{tr} \phi \rangle_{\mathrm{sing}} \approx -\frac{1}{32\sqrt{2\pi^3}} (4\pi - g)^2 \log(4\pi - g). \quad (6.12) $$

The dependence on $4\pi - g$, which is the parameter that measures deviation from criticality, is in complete agreement with predictions from [3, 47].

In the latter case $c = 0$ pure gravity regime we have to solve the equation (6.11) up to the 2nd order expansion around the critical$^{16}$ point $j'(z_c) = 0$. We find

$$ z_c = \frac{9}{32\sqrt{2\pi}} \quad (6.13) $$

which can be equivalently read off from (5.8), (5.9). Taking $z$ to be near this critical value so that

$$ z = \frac{9}{32\sqrt{2\pi}} + \sqrt{\epsilon}, \quad \text{where } \epsilon = \text{const} \times (g_c(M) - g) \quad (6.14) $$

<footnote>$^{16}$We recall that the critical line corresponds to $\partial g(M, z)/\partial z = 0$.</footnote>
we find the scaling for the one-point function

$$ \langle \frac{1}{N} \text{tr} \phi \rangle_{\text{sing}} \sim \epsilon^{3/2}. \qquad (6.15) $$

Remarkably, the term $\epsilon^{1/2}$ cancels, as it should be! This behavior thus perfectly matches the prediction from [1, 2].

### 6.1.2 Rational representation and rescaled form

Let us present some other useful representations of the singular part of the 1-pt function. One natural way to rewrite it is to exclude $\log(zM)$ from (6.10) using (6.11). Then we obtain

$$ \Phi_s \equiv 1152\pi^5 \frac{1}{M^2} \frac{1}{N} \langle \text{tr} \phi \rangle_{\text{sing}} = 16384\pi^3 j z^2 - 4608\pi^2 \sqrt{2} j z - 8388608\pi^5 z^3 + 1769472\pi^4 \sqrt{2} z^2 \qquad (6.16) $$

where we dropped the (polynomial in $j$) terms that do not contain $z$ and thus do not affect the critical behavior. We see that in terms of $z$ and $j$ the result is purely polynomial. This form can be used just as the one above for analysing the asymptotic regimes. For example, to find the $c = -2$ behavior we solve (6.11) iteratively and find

$$ z \approx \frac{j + 72\sqrt{2\pi} \log(M \frac{j}{512\pi^2})}{512\pi^2}. \qquad (6.17) $$

Then we plug this into (6.16) and take the leading term, which is precisely the result (6.12) we had before. The $c = 0$ limit (6.15) is also straightforward to take.

Another useful rewriting can be done by absorbing various coefficients into new rescaled variables. Thus we redefine the variables in (6.11) and (6.16) as

$$ z = \frac{t}{\frac{32\sqrt{2\pi}}{9}}, \qquad M = \frac{32}{9}\sqrt{2\pi\mu}, \qquad (6.18) $$

$$ g = 4\pi - 256\pi\Delta + 24\left(\sqrt{2} + 3\pi\right)M - 324\left(\sqrt{2} - \pi\right)M^2 \qquad (6.19) $$

and lastly

$$ \mathcal{J} = \frac{\Delta}{\mu} = \frac{j}{72\sqrt{2\pi}}. \qquad (6.20) $$

Then the coupling parameterization (6.11) becomes simply

$$ \mathcal{J} = t - \log(\mu t). \qquad (6.21) $$

In these variables the critical value of $z = z_c$ corresponds to $t_c = 1$. The physical range of $\mathcal{J}$ is then from $1 - \log \mu$ to infinity (corresponding to $z_c < z < \infty$). For the singular part of our one-point function from (6.16) we get

$$ \Phi_s = C[2\mathcal{J}(t^2 - 2t) + 3t^2 - 2t^3], \qquad (6.22) $$

where $C = 23328\pi^2\sqrt{2}$. Equations (6.22) and (6.21) represent the canonical parameterization of the (singular part of the) novel one point function $\langle \frac{1}{N} \text{tr} \phi \rangle$. This is one of the main results of this paper. We give a plot of the singular part of the 1-point function on figure 6.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Plot of singular part of $\langle \frac{1}{N} \text{tr} \phi \rangle$ given by (6.22) (without the overall constant $C$) as a function of normalized coupling, $\mathcal{J} = \Delta/\mu$ for fixed $\mu = 0.1$. At $\mathcal{J} \simeq 3.30259$ we see the $c = 0$ singularity $\langle \frac{1}{N} \text{tr} \phi \rangle \sim (g_c^{(c=0)} - g)^{3/2}$. The $c = -2$ singularity $\langle \frac{1}{N} \text{tr} \phi \rangle \sim (g_c^{(c=-2)} - g)^2 \log(g_c^{(c=-2)} - g)$ occurs at $\mu \to 0$ at fixed $g$, on the 1st sheet of the function (upper branch in the picture). The exponential singularity is on the 2nd sheet (w.r.t. the $c = 0$ branchpoint, lower branch on the picture).</figcaption></figure>

### 6.1.3 “Asymptotic freedom” on the second sheet

Let us consider the case $j \to 0$. This is of course already beyond the $c = 0$ criticality, which means that these are the “unphysical” values of parameters. However the underlying physical quantity may have singularities on the second sheet which correspond to subleading exponential corrections. Then we have from (6.10)

$$z \approx M^{-1} e^{-\frac{j}{72\sqrt{2\pi}}}. \qquad (6.23)$$

Plugging it into (6.16) and picking the leading term we find a behavior reminding asymptotic freedom

$$\Phi_s \approx -4608\pi^2\sqrt{2}(4\pi - g)e^{-\frac{4\pi - g}{72\sqrt{2\pi}M}}, \qquad (M \ll 4\pi - g). \qquad (6.24)$$

What are the excitations leading to these exponential effects? It cannot be the analogs of ZZ branes known for 2d gravity, since these ones have the $e^{-\text{const} \times N}$ behavior. They seem to be corrections to the regime of “almost spanning trees”, already for the leading order of planar graphs. We leave a more detailed exploration and interpretation of this regime for the future.

## 6.2 The $\frac{1}{N}\langle \text{tr} \phi^3 \rangle$ one-point function

Having studied above the 1-pt function $\frac{1}{N}\langle \text{tr} \phi \rangle$, here we will discuss another one, namely $\frac{1}{N}\langle \text{tr} \phi^3 \rangle$. This 1-pt function is particularly important as it is related to the derivative of the partition function in the mass $M$ and consequently to the fermion condensate of the type $\langle \bar{\psi} \psi \rangle$. As discussed above in section 4.2, it can be explicitly computed at generic values of the parameters and the result is given in the supplementary Mathematica file accompanying this paper. As an illustration, figure 7 shows a 3d plot of this observable as a function of $M$ and $\lambda$ in the physical region. Here we will study its expansion in the scaling limit we just discussed above.
<figure><img src="image_8.png" /><figcaption>Figure 7. Plot of the 1-pt function $\frac{1}{N}\langle \text{tr } \phi^3 \rangle$ as a function of $M$ and $\lambda$, in the physical region bounded by the critical line. We show the critical line in red, in the planes $\frac{1}{N}\langle \text{tr } \phi^3 \rangle = 1$ and $\frac{1}{N}\langle \text{tr } \phi^3 \rangle = 0$.</figcaption></figure>

For this correlator we find expansions similar to what we had before, starting with the small $M$ expansion

$$ \frac{1}{N}\langle \text{tr } \phi^3 \rangle = -\frac{59392\sqrt{2}}{105\pi} + 256 \qquad (6.25) \\ + \frac{64M(-90(21\sqrt{2}\pi - 94)\log(Mz) + 105\pi^2(128z - 27) - 8\sqrt{2}\pi(3760z - 333) + 16700)}{35\pi^2} \\ - \frac{32M^2}{7\pi^3}(3(7\sqrt{2}\pi^2(512z^2 + 5616z - 1215) - 6\pi(30112z - 2997) + 47562\sqrt{2})\log(Mz) \\ + 243(159\sqrt{2} - 70\pi)\log^2(Mz) + 2(-63\pi^3(2048z^2 - 864z + 135) \\ + \sqrt{2}\pi^2(315520z^2 + 24624z - 7047) + \pi(96624 - 692800z) + 57624\sqrt{2})) \qquad (6.26) $$

Then like before we define the singular part as

$$ \frac{1}{N}\langle \text{tr } \phi^3 \rangle_{\text{sing}} = \left[ \frac{1}{N}\langle \text{tr } \phi^3 \rangle - \frac{16(21\pi - 47\sqrt{2})}{7\pi^3}MJ \right] \Big|_{M^2} \qquad (6.26) $$

where we indicate that we take the term of order $M^2$ (dropping the terms of order $M^0$ and
$M^1$ which like for $\frac{1}{N}\langle \text{tr } \phi \rangle$ are regular). This gives

$$ \frac{1}{N} \langle \text{tr } \phi^3 \rangle_{\text{sing}} = -\frac{16M^2}{7\pi^3} \left( \begin{aligned} & 42(\sqrt{2}\pi^2(512z^2+6048z-1377)-96\pi(289z-36)+6432\sqrt{2})\log(Mz) \\ & +567(143\sqrt{2}-63\pi)\log^2(Mz)+4(-21\pi^3(10240z^2-4896z+405)+9\sqrt{2}\pi^2(56448z^2-14672z-783)+\pi(96624-476224z)+57624\sqrt{2})) \end{aligned} \right) \quad (6.27) $$

Finally rewriting this in terms of $\mathcal{J}$ we find, up to an overall factor, the result analogous to (6.22) for this correlator,<fnref n="17" />

$$ \Phi_{3,s} = 2 (\mathcal{J} - 6\sqrt{2}\pi + 8) (t^2 - 2t) + 3t^2 - 2t^3. \quad (6.28) $$

Notice that this function up to an overall multiplier coincides with the one for tr $\phi$ defined in (6.22) up to an overall multiplier and redefinition $\mathcal{J} \to \mathcal{J} + \text{const}$. This is a remarkable property which indicates a kind of universality for these 1-pt functions that would be important to elucidate further. As consequence, this function just like $\langle \text{tr } \phi \rangle$ also interpolates perfectly between the $c = -2$ scaling (6.12) and $c = 0$ scaling (6.15), with the limiting result having exactly the same form and the only difference being the overall constant factor. Furthermore, we see that both 1-pt functions are of course written in terms of the Lambert function which technically originates from the relation between $\mathcal{J}$ and the couplings in (6.21) which is the same for all observables of this type.

In principle one should be able to compute further 1-pt functions of the type $\langle \text{tr } \phi^n \rangle$ at least for fixed $n$ explicitly, and we leave this to future work.

# 7 Disc partition function

While in section 6 we have studied extensively the 1-pt function, here we will discuss a more complicated observable — the disc partition function. We will explore its critical behavior and the interpolation (flow) between $c = -2$ and $c = 0$ regimes.

We will study two resolvents: the resolvent $G(x)$ which generates $\langle \text{tr } X^k \rangle$ correlators,

$$ G(x) = \sum_{k=0}^{\infty} \frac{1}{x^{k+1}} \frac{1}{N} \langle \text{tr } X^k \rangle = \int_b^a dy \frac{\rho(y)}{x-y} \quad (7.1) $$

and another resolvent

$$ H(r) = \sum_{k=0}^{\infty} \frac{1}{r^{k+1}} \frac{1}{N} \langle \text{tr } \phi^k \rangle = \int_b^a dy \frac{\rho(y)}{r - \phi(y)} \quad (7.2) $$

generating the correlators $\langle \text{tr } \phi^k \rangle$. These resolvents can also be viewed as disc partition functions, with $k$ in the correlator being the disc boundary length, corresponding to different boundary conditions (see section 3.3 for details). We will show that they have the same universal behavior in the pure gravity limit as expected, while for the $c = -2$ regime we find different properties depending on the type of boundary.

<footnote n="17">¹⁷Here $\Phi_{3,s}$ is defined as $\frac{1}{N}\langle \text{tr } \phi^3 \rangle_{\text{sing}}$ multiplied by an overall constant and in which, like before, we further drop terms polynomial in $j$.</footnote>
## 7.1 The resolvent $G(x)$ for tr $X^k$ correlators

First let us discuss the $G(x)$ resolvent. We will work in the same limit as in section 6 described in (6.5), so that $M \sim p \to 0$. We will be interested in the singularity of the resolvent near its branch point which physically describes the situation when the boundary of the disc becomes large. Thus we expand $G(x)$ around the branch point at $x = a$ (which in our limit collides with the $x = c$ branch point present on the other sheet). In our limit we have

$$b = -2 + O(M), \quad a \approx 2 + M\alpha, \quad c \approx a + 64Mz \quad (7.3)$$

where

$$\alpha = \frac{3 \left( 3\sqrt{2} \log(Mz) + 6\pi + 2\sqrt{2} \right)}{\pi} \quad (7.4)$$

and we expand $x$ as

$$x = 2 + M\chi \quad (7.5)$$

with $\chi \sim 1$ being the rescaled boundary cosmological constant. Then, using (7.14) and dropping the terms which are regular in $\chi$, we get

$$G_{\text{sing}} = \sqrt{M} \left[ \frac{-9\sqrt{2} \cosh^{-1}\left(\sqrt{\frac{\chi-\alpha}{64z}}\right)}{\pi\sqrt{\chi-\alpha-64z}} - \sqrt{\chi-\alpha} \right] + O(M) \quad (7.6)$$

where we indicated that we keep only the singular terms. The singularity is at $\chi = \alpha$ (which in view of (7.3) corresponds of course to $x \approx a$) and comes about from the square root and the $\cosh^{-1}$ terms. Notice that the apparent singularity at $\chi = \alpha + 64z$ in fact cancels.

### 7.1.1 Limiting cases

Equation (7.6) gives the singular part of the resolvent which describes the critical behavior. Let us examine its $c = 0$ and $c = -2$ limits. In the former regime (pure gravity) we should expand as in section 6.1.1 around the critical value of $z$ given by $z_c = \frac{9}{32\sqrt{2\pi}}$ (see (6.13)), and the same time we expand in $\chi$ near the singularity,

$$z = z_c + \delta, \quad \chi = \alpha_c + 64\delta/\xi, \quad \delta \to 0 \quad (7.7)$$

with $\alpha_c = (\alpha)|_{z=z_c}$ and the factor 64 introduced for convenience. The variable $\xi$ is a finite scaling parameter and $\delta/\xi$ corresponds to the rescaled boundary cosmological constant while $\delta$ is viewed as the square root of the renormalised bulk cosmological constant (due to the relation (6.14) between $z$ and the coupling in this limit). Then we get with $\delta \to 0$

$$\frac{1}{\sqrt{M}} G_{\text{sing}} = c_1 + c_2 \frac{\delta}{\xi} + c_3 \frac{\delta^{3/2}}{\xi^{3/2}} (\xi - 2) \sqrt{\xi + 1} + \mathcal{O}(\delta^\epsilon) \quad (7.8)$$

where $c_n$ are (real) numerical constants. Thus we reproduce the universal pure gravity prediction of [12] — the expected scaling function $(\xi - 2)\sqrt{\xi + 1}$ and the expected prefactor $(\delta/\xi)^{3/2}$. This is yet another important test of our results.
Next, the $c = -2$ regime corresponds to scaling $z \to \infty$ and we take at the same time $\chi \to \infty$ so that the ratio

$$ \kappa \equiv \frac{\chi}{z} \qquad (7.9) $$

remains fixed. Then we find

$$ \frac{1}{\sqrt{M}} G_{\text{sing}} \approx -\sqrt{\kappa} \sqrt{z} - \frac{9\sqrt{2} \cosh^{-1}\left(\sqrt{\frac{\kappa}{64}}\right)}{\pi \sqrt{(\kappa - 64)z}}. \qquad (7.10) $$

The interpretation of this result remains to be clarified since the observables $\text{tr } X^n$ are not very natural for the $c = -2$ limit. Below we will see however that for the other resolvent $H(r)$ corresponding to $\text{tr } \phi^n$ correlators we recover perfectly the known results for the $c = -2$ limit.

## 7.2 The resolvent $H(r)$ for $\text{tr } \phi^k$ correlators

While we have computed the resolvent $G(x)$ in closed form, the resolvent $H(r)$ for $\text{tr } \phi^k$ vevs is harder to obtain exactly and we leave this question for the future. Nevertheless we will be able to compute here the universal part of it which is responsible for critical behavior, i.e. the singular part.

As for $G(x)$, we will be interested in the singularity that appears when the argument approaches the branch point of the resolvent. For $H(r)$ the $[b, a]$ cut in $x$ gets mapped to the cut $[\phi(b), \phi(a)]$ in the $r$ variable, so we will focus on the regime $r \to \phi(a)$. The most complicated part of the integral we wish to compute reads<fnref n="18" />

$$ \int_b^a dy \frac{y - 2c}{r/2 - c + \sqrt{c\sqrt{c-y}}} \sqrt{\frac{a-y}{y-b}} \Pi\left(\frac{y-b}{c-b}, 1-p\right). \qquad (7.11) $$

We will work as before in the limit $p \sim M \to 0$ so $a$ and $c$ are close. The singularity of the resolvent comes from the integration region $y \simeq a$ when the denominator becomes close to zero since at the same time $r \simeq \phi(a) \simeq 2c$. To zoom in on this singularity we make the change of variables

$$ y = a - M Y. \qquad (7.12) $$

In order to get a nontrivial result we would like both terms $(r/2-c)$ and $\sqrt{c\sqrt{c-y}}$ in the denominator of (7.11) to be of the same order (notice that both $y-a$ and $c-a$ are of order $M$), which means that $r - 2c \sim \sqrt{M}$. Thus we define the rescaled variable $R \sim 1$ by

$$ r = 2c + 2\sqrt{MR}. \qquad (7.13) $$

Now let us discuss the behavior of the $\Pi$ function in (7.11). We see that its argument is $1 - M \frac{Y+64z}{4}$ and its modulus is $1 - p = 1 - 16Mz$, so both of them approach 1 at the same rate $\propto M$. We did not find in the literature the expansion of $\Pi$ in this regime (which corresponds to a nontrivial resummation of the standard expansions in which one of the

<footnote n="18">¹⁸Notice the $\Pi$ function here is regular on $[b, a]$ but has a $\sim 1/\sqrt{x-c}$ behavior at $x=c$.</footnote>
two arguments is held fixed) but it can be quite straightforwardly derived from the integral representation of $\Pi$. The result reads

$$ \Pi(1 - S\epsilon, 1 - T\epsilon) \approx \frac{1 \operatorname{arccosh}(\sqrt{S/T})}{\epsilon \sqrt{S(S - T)}}, \quad \epsilon \to 0 \qquad (7.14) $$

which is also easy to verify numerically. Plugging in the values of our parameters and combining all the parts, for the full density (4.10) we find in our regime

$$ \rho \approx \frac{\sqrt{M}\sqrt{Y}}{\pi} - \frac{9\sqrt{2}\sqrt{M} \cosh^{-1}\left(\frac{1}{8}\sqrt{\frac{Y}{z} + 64}\right)}{\pi^2\sqrt{Y + 64z}}. \qquad (7.15) $$

Thus by focusing on the endpoint of the integration region we have got rid of elliptic functions. The integral in $Y$ can be now taken analytically as an indefinite integral. Plugging in the limits of integration, expanding for $M \to 0$ and discarding regular contributions we finally get for the singular part

$$ \frac{1}{M} H_{\text{sing}}(R) = \frac{\left(\pi R \sqrt{2} \sqrt{128z - R^2} - 18 \arccos\left(\frac{R}{\sqrt{128z}}\right)\right) \arccos\left(\frac{R}{\sqrt{128z}}\right)}{\pi^2}. \qquad (7.16) $$

We see that as expected it has nontrivial square root and logarithmic singularities in $R$ whose position moreover depends nontrivially on our finite scaling parameter $z$. The singularity is located as expected at the point where $r = \phi(a)$ which after expansion in $M$ translates to $R = -\sqrt{128z}$. It may seem that there is also a branch point at $R = +\sqrt{128z}$ but in fact it cancels between the different terms in (7.16).

This singular part of the resolvent (7.16) is another one of our main results. Below we will discuss its limiting cases corresponding to the pure gravity and spanning trees regimes.

Remarkably, the analytic structure of (7.16) reminds that of the analogous disc partition functions for the flattening of random geometries in the model of dually weighted graphs [74–77] (see also [78]). This similarity deserves a further study which we postpone to the future.

### 7.2.1 Limiting cases

To get the pure gravity limit, similarly to the discussion above in section 7.1.1, we expand near the critical value of $z$. Thus we set

$$ z = z_c + \delta, \quad R = -\sqrt{128z_c} + \frac{64\sqrt{\pi}}{3} \frac{\delta}{2^{3/4} X} \qquad (7.17) $$

where

$$ z_c = \frac{9}{32\sqrt{2\pi}} \qquad (7.18) $$

and expand for small $\delta$. This gives

$$ \frac{1}{M} H_{\text{sing}}^{(c=0)} \approx c_1 + c_2 \frac{\delta(X+1)}{X} + c_3 \frac{\delta^{3/2}}{X^{3/2}} (X-2) \sqrt{X+1} + \dots \qquad (7.19) $$

where $c_k$ are numerical constants. We see that the result perfectly reproduces the correct pure gravity scaling function [12] and the expected 3/2 power of $\delta$ (notice there is no $\delta^{1/2}$ term).
Let us now consider the $c = -2$ limit. To do this we expand the scaled disc partition function (7.16) with $z \simeq \frac{4\pi-g}{256\pi M} \sim R^2 \to \infty$ (the result of dropping the last term in (6.11) in this limit). We thus keep finite the combination $\zeta = \frac{R}{\sqrt{128z}}$. Then we find for the singular part of $H(r)$

$$H_{\text{sing}}^{(c=-2)} \approx -\frac{\sqrt{2}}{2\pi^2}(4\pi - g)\zeta\sqrt{\zeta^2 - 1}\log\left(\sqrt{\zeta^2 - 1} + \zeta\right). \qquad (7.20)$$

Remarkably, it perfectly coincides with the prediction for the $c = -2$ limit obtained in [32] (see equation (4.19) there) from a different matrix model. That shows once again the universality of the critical regime.

Thus we see that our result for the singular part of the resolvent interpolates nontrivially between two very different predictions in the $c = 0$ and $c = -2$ regimes and describes the flow between these two models.

## 8 Conclusion

In this study we investigated the model of massive spinless fermions interacting with 2d quantum gravity. We derived the Hermitian matrix model with non-polynomial potential describing the theory, and solved it in the planar approximation considering the one-cut solution. The regime where this solution exists is restricted by a critical curve in the 2-dimensional parameter plane of fermion mass $m$ and cosmological coupling $\Lambda = g_c - g$. It is explicitly demonstrated that the theory in the scaling limit $m^2 \sim \Lambda \to 0$ interpolates between the $c = -2$ theory, for $m^2 \ll \Lambda$, when the spanning trees dominate, and the pure 2d gravity $c = 0$ theory, for $m^2 \gg \Lambda$, when the fermions renormalize the cosmological constant in a simple way. We also computed the universal singular part of the disc partition functions in this scaling with Dirichlet and Neumann boundary conditions, interpolating between the $c = 0$ and $c = -2$ regimes. They fit perfectly with the previous results known for the limiting $c = 0$. The former one also fits the known $c = -2$ regime whehter as the latter one demonstrates in this limit a new behavior.

According to the matrix-forest theorem for the massive determinant we have identified the dominant number of trees in the forest in the different regions at the parameter space. In other words we took into account the backreaction of the massive matter on the 2d quantum geometry. At $m = 0$ the single tree saturates the partition function at criticality and we reproduce the picture for the $c = -2$ theory. The cosmological constant dependence of physical quantities (one-point functions, disc partition functions) contains logarithms due to the influence of “large” trees in this limit. At $m \to \infty$ the heavy matter breaks the 2d Euclidean space-time into the maximally possible number of components which coincides with the number of zero modes of the graph Laplacian.

The most interesting behavior occurs at small but finite $m$ where we find a new scaling behavior. The new scaling parameter is the ratio $J \sim \frac{\Lambda}{m^2}$ and the parameterization of the scaling functions is given in terms of the Lambert function (1.3) of the parameter $t$. In this limit sufficiently large and numerous trees in the partition function matter and it turns
out that this scaling regime exists at a narrow region near the critical curve. This scaling
describes the critical flow between $c = -2$ (in UV) and $c = 0$ (in IR) regimes.

There are a few questions concerning our solution which would be interesting to clarify:

* Study of the whole variety of flows in the vicinity of $c = 0$ and $c = -2$ critical points. Comparison to another flow found in [33]. Generalization of such flows to all central charges of matter $c \le 1$.

* Analysis of the rest of the parameter plane and of the multi-cut solutions to the matrix model.

* Exploring various 1-pt functions and understanding the origin of the simple relation between those we computed in (6.22) and (6.28)

* Clarification of the role of the second solution to the quadratic equation in the Parisi-Sourlas derivation of the matrix potential.

* Computation of instanton contributions of different kinds, including ZZ branes.

* Establishing the double scaling limit $N \to \infty, \Lambda \to 0$ of our model along the critical line and deriving the universal scaling function in this limit.

* Derivation of this and other critical flows from the continuous 2d QG (Liouville formalism).

We hope to return to these questions in our future research.

# 9 Further directions

Here we outline in more detail some nontrivial potential directions for future exploration.

## 9.1 Double scaling limit and sum over topologies

It should be possible to solve our one-matrix model in the double-scaling limit [13, 15, 16], for the whole critical flow in the space of $(M, \lambda)$. It would be interesting to embed such a double scaling solution into the KdV formalism of [14] for 2d gravity interacting with $c < 1$ matter fields. We could also expect other integrability pattern in our double scaling limit with the Toda hierarchy involved like in [64]. This is quite common for the theories with asymptotic freedom. The emergence of the Lambert function supports this expectation.

Is it possible to get the analogue of the Kontsevish matrix model in our case? The fermionic bilinears (derivative with respect to mass) should reproduce some classes at the moduli space. In the double scaling limit the Lambert function emerges implying the relation with the Hurwitz numbers which indeed according to ELSV formulae are written as particular integrals over the moduli space.
## 9.2 Generalization to other critical flows

An obvious generalization of our formulas describing the universal flow between $c = -2$ and $c = 0$ critical points (1.2), (1.3), (7.6), (7.16) would be the construction of a more general model with forests on planar graphs, working for all central charges $-\infty < c < 1$. Following the Kastelyn-Fortuin-Stephens tree expansion [79] for the Q-state Potts model, realized on planar graphs as a Q-matrix model in [10], we have to introduce the loops into the trees and weigh such configuration with an extra factor $Q^{\#\text{loops}}$. The central charge will depend on Q. It would be interesting to construct such a (multi)matrix model, which will be different from the model of [10] since we deal here with rooted trees. As we have seen on the example of the flow between $c = -2$ and $c = 0$ in our paper, such a model would have universal critical flows different from those of the $O(n)$ model of Kostov [11].

Furthermore, in [55] an interesting mapping of the model of unrooted trees considered there to the (loop+dimer) statistical model of Kostov and Staudacher [32] has been developed. What is the similar corresponding statistical model in our case of rooted trees?

## 9.3 Fragmentation of the RRG into finite number of trees and many-body localization

Recently the RRG ensemble has attracted a lot of attention being the toy model for a Hilbert space of some interacting many-body problem (see [69] for review and references therein). One considers the spinless fermion $\psi_i, i = 1...N$ on RRG at large $N$ with diagonal on-site disorder. The partition function of the model reads

$$Z(W, N) = \sum_{\text{RRG}} d\psi d\psi^{\dagger} \exp(\psi_i^{\dagger}(L_{ij} + \delta_{ij}\epsilon_i)\psi_j) \quad (9.1)$$

where $\epsilon_i$ is the random diagonal disorder with the flat distribution $\epsilon_i \in (-W, W)$. Let us compare this model with our study. Both models describe the spinless massive fermion interacting with the 2d gravity, although in (9.1) the mass of the fermion is random while we consider the fermion with the fixed mass. The second difference is that we consider the canonical ensemble with the cosmological constant while in (9.1) the microcanonical ensemble is assumed.

The models are very close but the questions discussed are quite different. We have integrated out the fermions and look at the emergent partition function with the determinant in the critical regime of large number of nodes as a function of two couplings. The usual question in the model (9.1) is different and concerns the localization or delocalization of the fermion at the graph due to the disorder. It was found [69] that there exists a critical $W_{cr}$ such that for $W > W_{cr}$ the Anderson localization takes place. This has been established via numerical evaluation of level spacing distribution or IPR.

The one-particle Anderson localization on RRG itself seems to be a somewhat artificial problem, however it becomes interesting if we treat the RRG ensemble as the model of Hilbert space for some interacting many-body system. The mapping is not exact but it captures the key qualitative features. The one-particle Anderson localization transition in the Hilbert space is treated according to the conjecture from [68] as the transition to the many-body localized (MBL) phase in the physical space. It was argued (see [70] for
review and references therein) that the fragmentation of the Hilbert space is one of the key mechanisms for the transition into the MBP phase. In terms of RRG it means that we are looking at the strong backreaction of fermions leading to fragmentation of RRG into some number of weakly connected subgraphs, possibly trees.

Another mechanism of Hilbert space fragmentation involves not the random fermion mass but the perturbation of the RRG ensemble by the chemical potentials for the short cycles. These chemical potentials can be considered as the leading terms of the expansion of the characteristic polynomial of the graph Laplacian in inverse powers of $m^2$ since $\text{tr}A^n$ gives the number of cycles of length $n$ on the graph $\log\det(A + m^2) \propto \sum_k(m^2)^{-k}\text{Tr}A^k$. Here we have a clear-cut link to our model in the large $m^2$ limit.

If a cubic perturbation $V = t_3\text{Tr}A^3$ is chosen, the phase transition occurs at some critical value $t_{3,crit}$ and the RRG ensemble at $t_3 > t_{3,crit}$ gets dominated by the clustered graph with the number of clusters equal to $N_{cl} = \frac{|G|}{q}$ [80]. Similar graph fragmentation occurs for $V = t_4\text{Tr}A^4$ at $t_{4,crit}$ when bipartite clusters emerge [81, 82]. From the spectral viewpoint of the graph Laplacian each cluster represents the single low energy mode escaped from continuum. The isolated eigenvalues form the second soft “non-perturbative” band in the spectrum [80]. It turns out that the spectrum of the perturbed RRG ensemble enjoys the mobility edge which separates the delocalized states in the main part of the spectrum and localized modes in the second non-perturbative band [83].

We do not expand the massive determinant, hence the chemical potentials for all cycles are present. The number of trees in the forest, which is the order parameter for the Hilbert space fragmentation, in our exact solution in the planar limit depends on the mass and the cosmological constant which provides the soft cut-off in the Hilbert space dimension. At small mass we have no Hilbert space fragmentation at all. The investigation of the localization of the fermions in the gravity background in our model is a clear direction for further study. In order to explore this, more detailed characteristics like the level spacing distribution have to be analysed. We postpone this issue for a separate investigation.

We could try to solve this model in the double scaling limit, summing up over the topologies in the critical regime (near our critical curve, and presumably at $m \to 0$). One might hope (though it is not at all guaranteed) that the double scaling solution contains this fragmentation phenomenon.

## 9.4 ZZ-brane instantons for one-cut solution

Let us make a comment concerning the instanton effects and consider the single eigenvalue ZZ-brane instanton for our one-cut solution. The instanton action evaluated along the spectral curve reads as follows

$$S_{\text{inst}} = \int_a^{x_0} Y(x)dx \qquad (9.2)$$

where the spectral curve in terms of the resolvent reads

$$Y(x) = V'_{\text{eff}}(x) = V'(x) - 2G(x) = M(x)\sqrt{(x-a)(x-b)}. \qquad (9.3)$$

The critical point $x_0$ of the effective potential $V_{\text{eff}}(x)$ is defined by condition $M(x_0) = 0$ and corresponds to the pinch point of the spectral curve. The effective potential is constant
on the cut and its derivative obeys the useful relation [5]

$$ \partial_{\lambda} V_{\text{eff}}(a) = 2 \log \frac{(a-b)}{4} \qquad (9.4) $$

The instanton action can be written as

$$ S_{\text{inst}} = V_{\text{eff}}(a) - V_{\text{eff}}(x_0). \qquad (9.5) $$

Let us focus on the instanton contribution at $M \to 0$ when we can approximate the r.h.s. by derivative since $x_0$ is close to $a$,

$$ S_{\text{inst}} \propto \partial_a V_{\text{eff}}(a)(x_0 - a) \propto \frac{\partial c}{\partial a} \frac{\partial V_{\text{eff}}}{\partial c}(a)(x_0 - a) \propto \log(b - a)(x_0 - a). \qquad (9.6) $$

We recall that $c = \frac{1}{4\lambda}$. Using the identity (9.4) and expansions for $a, c$ from appendix D one can check that $\frac{\partial c}{\partial a}$ is finite in this limit. We see that the instanton action vanishes at $M \to 0$ and instanton contributions become unsuppressed.

## 9.5 Kesten-McKay distribution and criticality without planarity

So far we have considered the planar approximation but here we shall make a short remark concerning the partition function for the generic RRG microcanonical ensemble at large number of nodes $n \to \infty$. In this limit we can utilize the famous Kesten-McKay(KM) distribution for the spectral density of RRG ensemble [66, 67].

Consider the derivative of our partition function $Z(m^2)$ at the critical line which yields the resolvent of Laplacian for RRG $R_n(M)$

$$ R_n(m^2) = \frac{1}{n} \langle \text{Tr} \frac{1}{L - m^2} \rangle = -\frac{2}{n} \frac{d}{dm^2} \langle \log Z(m^2) \rangle \qquad (9.7) $$

The resolvent of the adjacency matrix of RRG reads [66, 67]

$$ R_n(z) = \frac{(2-q)z + q\sqrt{z^2 - 4(q-1)}}{2(z^2 - q^2)} \qquad (9.8) $$

and the resolvent of the Laplacian of RRG can be derived via the shift $z \to q - z$.

Hence from the spectral density of RRG

$$ \rho_n(\lambda, q) = \frac{1}{n} \langle \sum_i \delta(\lambda - \lambda_i) \rangle_{\text{RRG}} = \frac{q \sqrt{4(q-1) - \lambda^2}}{2\pi q^2 - \lambda^2} + O(1/n) \qquad (9.9) $$

we obtain for the spectral density of the Laplacian

$$ \rho(m^2) = \frac{q \sqrt{4(q-1) - (q-M)^2}}{2\pi q^2 - (q-M)^2} + O(1/n) \qquad (9.10) $$

with density support $|(q - m^2)| < 2\sqrt{q - 1}$. Recently the $1/n$ correction to the RRG spectral density has been evaluated [84].

We observe two non-analyticities: the pole due to the zero modes in the spectral density at $m^2 = 0$ and the branch points at $4(q - 1) = (q - m^2)^2$. The branch points seem to be
interesting and amount to criticality and nontrivial “susceptibility” $\frac{d^2\langle\log Z(m^2)\rangle}{d^2m^2}$ for $q=3$ at points $m_{crit}^2 = 3 \pm 2\sqrt{2}$. There is no such regime in the matrix model studied here: our graphs are planar which drastically changes the critical behavior. We could ask the question whether this kind of distribution can occur in the double-scaled limit of the model, when we sum up over all genera of graphs with the weight $N^{2-2\times\text{genus}}$, close to criticality at every genus. This is however beyond the scope of our paper.

## 9.6 Analogy with QCD matrix model

Let us consider two-dimensional QCD, that is, fermions interacting with the quantum gauge fields. The partition function is the determinant of the Dirac operator averaged over gauge fields $\langle \det(\hat{D}(A) + iM) \rangle_{YM,g_{YM}}$ and it can be approximated by the large $N$ matrix model (see [85] for a review) which has the following interpretation. The ground state is assumed to be populated by the instantons and antiinstantons hosting the fermion zero modes. The Wishart matrix in the matrix model represents the Dirac operator in the basis of zero modes or speaking a bit differently the overlap of zero modes which get collectivized.

In our case we have fermions interacting with gravity instead of the gauge field, determinant of the massive Laplace operator instead of the determinant of Dirac operator and we do not distinguish chirality. Therefore we have an analogue of the “instantons without antiinstantons” situation where each tree is the analogue of instanton hosting zero mode of the graph Laplacian. The number of the trees due to the index theorem coincides with the number of zero modes hence like in QCD we could have a picture of the overlap of zero modes upon switching on gravity. Our finding for massive fermions interacting with gravity suggests that the formation of the multiinstanton clusters in 2d QCD can be expected as a function of the ratio of the gauge coupling constant and fermion mass.

Fermionic condensates considered in this study have many similarities with the gluino condensates in SYM. There are still subtle problems concerning the evaluation of gluino condensates due to the factorization of the higher topological correlators. It would be interesting to investigate the factorization issue for the higher fermionic correlators in our case.

## Acknowledgments

We thank J. Bouttier, P. di Francesco, S. Komatsu, M. Mariño, S. Nechaev, D. Serban, A. Sportiello and N. Terziev for discussions. We are especially grateful to I. Kostov for illuminating comments. F.L.-M. is grateful for hospitality to organisers of the Varna ICMS-2022 workshop (organised by the International Center for Mathematical Sciences in Sofia and supported by the Simons Foundation). The research of V.K. was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. V.K. thanks the Interdisciplinary Scientific Center Poncelet (CNRS UMI 2615), where a part of the work was performed, for kind hospitality. A.G. is thankful to the Basis Foundation grant No. 20-1-1-23-1. The work of V.M. was funded by the “Basis” Foundation grant No. 22-1-1-42-3 and RFBR grants No. 21-52-52004 and No. 20-01-00644.
# A Combinatoric explanation of the Kirchoff-Tutte matrix-tree theorem

**Theorem.** The number of spanning trees of a connected loopless graph $G$ is equal to the determinant of $\Delta'(G)$ — the Laplacian of the graph without its last column and row. The Laplacian itself is defined as:

$$ \Delta = -Q + A \qquad (A.1) $$

as in (2.1).

Let $K_{k,\alpha}(G)$ be the (oriented) incidence matrix of a graph $G$ with oriented edges labeled as $\alpha = 1, \dots, E_G$, and vertices labeled as $k = 1, \dots, n$, constructed as follows:

$$ K_{k,\alpha} = \begin{cases} 1, & \text{if the edge } \alpha \text{ points from vertex } k \text{ outside} \\ -1, & \text{if the edge } \alpha \text{ points from vertex } k \text{ inside} \\ 0, & \text{otherwise.} \end{cases} \qquad (A.2) $$

It is easy to see that $\Delta_{ij} = \sum_{\alpha} K_{i,\alpha}K_{i,\alpha}$, or $\Delta = KK^T$. Moreover this equatily “survives” on the level of the first minor — $\Delta' = K'K'^T$, where $K'$ is the Kirchoff matrix with erased last line.

The proof is based on Cauchy-Binet formula: for two rectangular $m \times M$ matrices $F$ and $H$, we have

$$ \det(FH) = \sum_{S_m \in S_M} \det F_{S_m} \det H_{S_m} \qquad (A.3) $$

where by $\det F_{S_m}$ we denote one of the maximal $m \times m$ minors of the matrix $F$ (and similarly for $H$). The sum goes over all $\binom{M}{m}$ such minors.

Applying it to $L' = K'K'^T$ we write

$$ \det(\Delta') = \sum_{S_m \in S_M} \det K'_{S_m} \det K'^T_{S_m} = \sum_{S_m \in S_M} (\det K'_{S_m})^2 \qquad (A.4) $$

Each matrix $K_{S_m}$ is the Kirchoff matrix of a subgraph obtained by cutting in $G$ all the rest of $E_G - n$ edges.

Next, we notice that the matrix $K'$ is unimodular: all its maximal minors are equal to 0, $\pm 1$. It is zero if the subgraph has at least one cycle or is disconnected, and $\pm 1$ otherwise (which means that it is a spanning tree). Indeed, each maximal square submatrix $K_S \in K$ is in fact the incidence matrix of the subgraph obtained by cutting all the edges corresponding to the erased columns of $K$. Then $K_S K_S^T$ is the laplacian on this subgraph. Then $\det(K_S K_S^T)$ is nonzero only if the graph is connected, or it has no loops (the connectivity and the looplessness are actually the same of such maximal subgraphs). But in this case it is a spanning tree, and it is easy to see that for the spanning tree $\det(K_S) = \pm 1$ (depending on the orientation): each column contains only one 1 and one (-1), and they can be always ordered in an upper-triangular way.

Then it is clear that

$$ \det(\Delta'(G)) = \#\text{of spanning trees of } G. \qquad (A.5) $$

Quod erat demonstrandum!
## B Some properties of the spectral determinant

* Define the spectral determinant of undirected non-weighted graph G Laplacian L as

$$L(z, G) = \det(L - I z) \qquad (B.1)$$

The matrix-forest theorem claims

$$L(z, G) = \sum_{i} (-1)^{i} a_{i} z^{n-i} \qquad (B.2)$$

where $(n-i)$ is the number of trees in the forest and $a_i$ is the product of number of nodes in all trees in $(n-i)$-component forest.

$$\frac{d}{dz} L(z, G)|_{z=0} = a_{n-1} = n|t(G)| = \prod z_i \qquad (B.3)$$

where $z_i$ are roots of the Laplace polynomial. It is the form of matrix-tree theorem.

* The reciprocity

$$L(z, \bar{G}^p) = (-1)^{n-1} L(np - z, G) \qquad (B.4)$$

$$t(K^p - G) = (sp)^{s-n-2} L(sp, G) \qquad (B.5)$$

where $\bar{G}^p$ is graph p-complemented to G, s is the number of nodes in the full graph $K^p$.

* The complementarity

$$zL(z, G_1^p G_2) = (z - n_1 p - n_2 p) L(z - n_2 p, G_1) L(z - n_2 p, G_2) \qquad (B.6)$$

where $G_1, G_2$ are graphs with the number of nodes $n_1, n_2$

* The spectral determinant can be generalized a bit if we attribute the weights for the nodes $x(v_i)$ (contrary to the standard weights for links $\omega(v_i, v_j)$) forming the weight vector $\vec{x}$. The generalization of the matrix-forest and matrix-tree theorems does exist in this case and reads as [86]

$$z^{-1}(-1)^{n-1} \det(M(x, G) - I z) = \mathcal{F}(z, x, G) \qquad (B.7)$$

where $M(x, G) = (m_{ij} = -x(v_j)\omega(v_i, v_j))$ is the Laplacian matrix of the graph dressed by the node's degrees of freedom and $\mathcal{F}(z, x, G)$ is the so-called forest volume of G. If $x(v_i) = 1$ it reduces to the matrix-forest theorem.

$$(-1)^{n-1} L(-z, G) = \mathcal{F}(z, 1, G) \qquad (B.8)$$

## C Elliptic functions

Here we collect useful integrals and relations for elliptic functions we use.
### C.1 Some integrals

Let us give more details on taking the integral (4.6). We assume that $b < a < c$. For the first two terms it reduces to elementary functions,

$$ \int_b^a dy \frac{9cM^2 + y}{(x-y)\sqrt{(a-y)(y-b)}} = \frac{\pi \left(-\sqrt{(a-x)(b-x)} + 9cM^2 + x\right)}{\sqrt{(a-x)(b-x)}} \quad (C.1) $$

whereas the last one contains the complete elliptic integral of 3-rd kind:

$$ \int_b^a \frac{dy}{(x-y)} \frac{y-2c}{\sqrt{(a-y)(y-b)(c-y)}} = -\frac{2K\left(\frac{a-b}{c-b}\right)}{\sqrt{c-b}} + (x-2c) \frac{2\Pi\left(\frac{a-b}{x-b}, 1-p\right)}{(x-b)\sqrt{c-b}} \quad (C.2) $$

Combining them together we find the result (4.7) for $G(x)$.

### C.2 Relations used for integration of elliptic functions

The relations we use to take the integral (4.16) read:

$$ \begin{aligned} n(n-1)^{k/2} \sqrt{\frac{l}{n} - 1} \Pi(n|l) &= -\frac{(kn + n + 2)(n - 1)^{k/2} \sqrt{\frac{l}{n} - 1} (K(l)(l - n) + nE(l))}{(k + 1)(k + 3)n(l - n)} + \\ & \qquad + \frac{d}{dn} \left( \frac{2(n - 1)^{\frac{k}{2} + 1} \left( \frac{2}{k+1} + n \right) \sqrt{\frac{l}{n} - 1} \Pi(n|l)}{k + 3} \right) \\ (n - 1)^{k/2} \sqrt{\frac{l}{n} - 1} \Pi(n|l) &= -\frac{(n - 1)^{k/2} \sqrt{\frac{l}{n} - 1} (K(l)(l - n) + nE(l))}{(k + 1)n(l - n)} + \\ & \qquad + \frac{d}{dn} \left( \frac{2(n - 1)^{\frac{k}{2} + 1} \sqrt{\frac{l}{n} - 1} \Pi(n|l)}{k + 1} \right) \end{aligned} \quad (C.3) $$

### D Expansions at small $M$ and $p$

We have to order $M^2$ (dropping also terms of higher order than listed in $L$ and $p$)

$$ \begin{aligned} a = M^2 & \left( \frac{27 (25551 L^2 + 3 (38879 + 5778 \sqrt{2}\pi) L + 3072\pi^2 + 44738 \sqrt{2}\pi + 129502) p^2}{256\pi^2} \right. \\ & + \frac{27 (465 L^2 + 2 (841 + 255 \sqrt{2}\pi) L + 4 (371 + 268 \sqrt{2}\pi + 48\pi^2)) p}{16\pi^2} \\ & \left. + \frac{27 (33 L^2 + 60 (1 + \sqrt{2}\pi) L + 8 (2 + 9\sqrt{2}\pi + 6\pi^2))}{8\pi^2} \right) \\ & + M \left( - \frac{3 (1455 \sqrt{2}L + 768\pi + 2869\sqrt{2}) p^2}{64\pi} - \frac{3 (51 \sqrt{2}L + 48\pi + 80\sqrt{2}) p}{4\pi} \right. \\ & \left. - \frac{3 (3\sqrt{2}L + 6\pi + 2\sqrt{2})}{\pi} \right) \\ & + 2 \end{aligned} \quad (D.1) $$
and

$$ b = M^2 \left( \frac{27(-7317L^2 + (8052\sqrt{2}\pi - 96055)L + 3072\pi^2 + 14392\sqrt{2}\pi - 215062)p^2}{256\pi^2} + \frac{27(-87L^2 + 54(4\sqrt{2}\pi - 33)L - 4(943 + \sqrt{2}\pi - 48\pi^2))p}{16\pi^2} + \frac{27(-3L^2 + 4(6\sqrt{2}\pi - 35)L - 8(22 + 11\sqrt{2}\pi - 6\pi^2))}{8\pi^2} \right) + M \left( \frac{3(198\sqrt{2}L - 384\pi + 1589\sqrt{2})p^2}{32\pi} + \frac{3(3\sqrt{2}L - 24\pi + 67\sqrt{2})p}{2\pi} + \frac{42\sqrt{2}}{\pi} - 18 \right) - 2 \quad (D.2) $$

and

$$ c = M^2 \left( \frac{27(35535L^2 + 3(59487 + 7730\sqrt{2}\pi)L + 3072\pi^2 + 67074\sqrt{2}\pi + 219742)p^2}{256\pi^2} + \frac{27(537L^2 + 6(347 + 97\sqrt{2}\pi)L + 4(467 + 348\sqrt{2}\pi + 48\pi^2))p}{16\pi^2} + \frac{27(33L^2 + 60(1 + \sqrt{2}\pi)L + 8(2 + 9\sqrt{2}\pi + 6\pi^2))}{8\pi^2} \right) + M \left( -\frac{9(853\sqrt{2}L + 256\pi + 2439\sqrt{2})p^2}{64\pi} - \frac{9(21\sqrt{2}L + 16(3\sqrt{2} + \pi))p}{4\pi} - \frac{3(3\sqrt{2}L + 6\pi + 2\sqrt{2})}{\pi} \right) + 4p^2 + 4p + 2 \quad (D.3) $$

## E Lambert function and brane insertion

### E.1 Two examples with Lambert function

In this section we shall comment on the reason for the Lambert function to appear in our double scaling limit. To preclude the arguments below recall our setting. We start with fermions on the fixed graph and represent the massive determinant as the weighted sum over the separated trees in the forest. Then we switch on the gravity the interacting system itself selects the preferable state depending on the point at the two-dimensional parameter space $(m, \lambda)$. Hence in any representation of our model we need the boundary creation operator in some form yielding the boundaries of the trees.

Before turning to our model let us describe two models where the Lambert function has emerged in the very similar context. First, consider the setup discussed in [62] where
the perturbation of the 2d topological gravity via the boundary creation operator has been considered. The motivation of that study was as follows. The boundary matrix model for JT gravity found in [37] involves the important contribution from the replica wormholes providing the interaction of baby universes. It was suggested in [62] to substitute the standard replica trick with useful integral representation for quenched free energy $\langle \log Z \rangle$. It is equivalent to the replica representation with particular analytic continuation to $n \to 0$. In this representation the quenched free energy reads

$$ \langle \log Z \rangle - \log \langle Z \rangle = - \int_0^\infty \frac{dx}{x} [e^{-\tilde{Z}x} - e^{-\langle Z \rangle x}] \quad (E.1) $$

and involves the function $e^{\tilde{Z}(x)}$ which is the generating function for connected correlators $\langle Z^n \rangle_c$. These correlators provide the multiple replica boundary contributions in JT gravity partition function. It has been also identified as the operator creating the peculiar space-time D-brane [87] introduced in the context of replica wormholes. Remarkably if we consider the Airy limit of the Gaussian model in genus zero this generating function can be evaluated exactly. The exact boundary creating operator $\hat{Z}$ in topological gravity introduced in [45] is applied to the partition function of topological gravity

$$ \langle Z^n \rangle_c = (\hat{Z})^n \mathcal{F} \quad (E.2) $$

where

$$ \hat{Z} = g_s \sqrt{\frac{\beta}{2\pi}} \sum_k \beta^k \partial_k \quad (E.3) $$

and derivative is taken with respect to k-th time. If one restricts oneself to the genus zero the Lambert function emerges as the generating function for boundary creating operator in topological gravity at genus zero in the Airy limit. This is the first role of the Lambert function.

The second role of Lambert function is important as well [62]. It is familiar in the topological string context that the insertion of the brane shifts the closed moduli in background which is usually written in symbolic relation $Z_{\text{closed}}(\vec{t}') = Z_{\text{brane}}Z_{\text{closed}}(\vec{t})$ (see, for instance [63]). In the case under consideration this equation reads

$$ \tilde{Z}(x) = \mathcal{F}(t_k) - e^{-\hat{Z}x} \mathcal{F}(t_k) = \mathcal{F}(t_k) - \mathcal{F}(t'_k) \quad (E.4) $$

hence we expect that “Lambert brane” amounts to the closed moduli shift. This argument turns out to be true and the KdV flows in topological gravity yield the simple derivation of this shift [62]. To this aim consider the case when only two first times in hierarchy are switched on $t_k = 0, k \le 2$. The relevant solution to the KdV hierarchy in this case depends on two times $t_0, t_1$ and before the brane insertion it reads

$$ u(t_0, t_1) = \frac{t_0}{1 - t_1} \quad (E.5) $$

where $u = \partial_0^2 \mathcal{F}$. The shift of the moduli can be seen from the string equation which can be written for KdV in terms of Gelfand-Dikii polynomials $R_k$ as follows

$$ u(t) = \sum_k t_k R_k \quad (E.6) $$
If there are only two non-vanishing times upon the shift string equation can be brought into the Lambert form [62]

$$z = We^W \qquad (E.7)$$

for

$$W = \frac{\beta t_0}{1 - t_1} - \beta v \qquad z = \frac{x}{1 - t_1} e^{\frac{\beta t_0}{1 - t_1}} \qquad (E.8)$$

Hence the second role of Lambert function — shift of the background via the brane insertion when only two first times in integrable hierarchy are switched on.

The second model enjoying the Lambert function is the topological A model string on $CP_1$ coupled to topological gravity. The theory has two equivalent dual representations [64]: a) all genus topological A type string on $CP_1$ supplemented by the first gravitational descendant of Kahler form and b) the massive two-dimensional fermions with specific gravity induced action. According to the topological string framework the inclusion of descendant corresponds to the insertion of some brane similar to the topological gravity above. The first representation of Lambert function as the brane creation operator is not explicitly known in this model hence we focus on the second role — shifting the moduli. This aspect has been identified in [64] in explicit form.

It goes as follows. The model is solved in terms of the spectral curve which is the sphere with two marked points. The cut on the C plane lies between $x_{\pm} = v \pm \Lambda$. The filling fraction for fermions is defined as $a = \hbar R$ and plays the role of closed moduli. This theory similarly to the topological gravity case involves only two times $t_1, t_2$ and it was found in [64] that instead of KdV hierarchy for pure gravity here the semiclassical limit of Toda hierarchy does the job. When $t_2 = 0$ we have no gravity descendant of Kahler form and the relation $v = a$ holds. When $t_2 \neq 0$ solution requires the proper matching condition at the ramification points which yields shift of the closed moduli. The Lambert function $W(z)$ enters via the following matching conditions for moduli shift

$$v - a = W(t_2 e^{t_1 + at_2}) \qquad (E.9)$$

$$\log \Lambda = t_1 + 2t_2 a - \frac{1}{2W} \left( -16t_2^2 e^{2(t_1 + 2t_2 a)} \right) \qquad (E.10)$$

As in the previous case we interpret this relation as an effect of insertion of the “Lambert brane” and its backreaction on gravity.

Integrability implies the fermionic representation of the partition function [64]

$$Z(t_1, t_2, R) = \langle R | e^{-\frac{J_1}{\hbar}} e^{\hbar^{-1}(t_2 W_3 + t_1 W_2)} e^{-\frac{J_1}{\hbar}} | R \rangle \qquad (E.11)$$

where $|R\rangle$ is the state with U(1) charge $R$, $\omega$ is the coordinate on the cylinder, $t_1 = \log \Lambda$. The harmonics of the U(1) current are defined as

$$J_k = \sum_r : \tilde{\psi}_r \psi_{r+k} : \qquad (E.12)$$

The generators of the $W_{1+\infty}$ algebra are expressed in terms of fermions as follows

$$W_{k+1} = -\frac{\hbar^k}{k+1} \oint d\omega (\tilde{\psi}[(D + 1/2)^{k+1} - (D - 1/2)^{k+1}]\psi) \qquad (E.13)$$
where $D = \omega\partial_{\omega}$. Hence we observe that the fermions are effectively massive and have a non-standard kinetic term involving the second derivative.

The time $t_2$ is coupled to the unusual second derivative term induced by coupling to gravity and coinciding with the zero mode of the $W_3$ generator of the $W_{1+\infty}$ algebra. Upon bosonization it can be expressed as the $(\partial\phi)^3$ term in the action for the chiral boson and has an interpretation as a cut-and-join operator. Such terms are familiar in many models (see for instance [88] for the relevant discussion) and the coefficient in front of this term has an interpretation of the string coupling.

This theory has one more interpretation which actually was the initial one for the authors of [64]. The model can be viewed as the $D = 4$ abelian $U(1)$ $N = 2$ SYM theory in the $\Omega$ background when the leading gravitational correction to the prepotential is taken into account. That is, in the UV the prepotential involves only two times

$$ \mathcal{F}_{\text{UV}} = t_1 \text{Tr} \Phi^2 + t_2 \text{Tr} \Phi^3 $$

where the second term is induced by the coupling to gravity. Naively there are no instantons in the abelian theory but they do emerge when the coupling to gravity is switched on and the theory is asymptotically free with a non-perturbatively generated IR scale. Presumably this theory can be considered as the worldvolume theory on the inserted brane.

## E.2 Lambert in the forest

From two examples above we see that the Lambert function plays the role of the generating function for multiple boundary insertions and simultaneously shifts the closed moduli. In our case we indeed need tree boundary creation operators, however according to the matrix-forest theorem the boundary of each tree enters with the volume of the tree. Hence the second $CP^1$ example is more similar to our case since the effect of the gravity descendant of the Kahler class involves the volume factor indeed. As in that case we will observe the shift of the closed moduli.

Let us first note that the Lambert function $W(x)$ in (6.11) in our notation is just the $z$ variable or more precisely its rescaled version $t$ (see (6.18)),

$$ t = -W \left( -\mu^{-1} e^{c_0 \frac{g-g_c}{\mu}} \right) \qquad (\text{E.14}) $$

where $W(x)$ obeys the equation $W(x)e^{W(x)} = x$ and numerical factor $c_0 = \frac{1}{256\pi}$. Recall that near the critical point the area of the surface $A$ behaves as $A \propto \frac{1}{g-g_c}$ hence we have $e^{-\frac{1}{MA}}$ in the argument of the Lambert function. This can be compared with the standard instanton exponent $e^{-\frac{1}{g_s A}}$ which implies the suppression of the instanton effects at large $N$ since $g_s = 1/N$. In our case it seems that $g_s \propto M$ which implies the lack of large $N$ suppression. Notice that the factor like $e^{-\frac{1}{g_s A}}$ was discussed for the instanton contributions in 2d YM theory.

Let us turn now to two roles of the Lambert function observed in the previous examples. First, one can ask whether we have an effective brane generating multiple boundaries on the worldsheet like in the Airy limit of topological gravity at genus zero. It is a well known
fact that the insertion of det($M - x$) in the matrix model for the type B topological string corresponds to the insertion of the FZZT brane at point $x$ of the spectral curve of the matrix model. In our case the massive determinant det($\Delta + M$) involved in our partition function according to the matrix-forest theorem seems to play such a role. Indeed its expansion in mass provides multiple trees on the worldsheet. Hence we expect that we have the effective brane insertion as well and the mass $M$ provides the position of the insertion. Certainly this point deserves further clarification.

Secondly, similarly to the example of the topological A model string on $CP^1$ we can look for the equation describing the shift of the closed moduli via the insertion of the brane. The shift has to be proportional to the deformation parameter $M$. The equation entering our one-gap solution

$$p = 16Mz = MW(x) \qquad (E.15)$$

plays this role. When $M = 0$ and $p \to 0$ the elliptic curve degenerates into the marked sphere while the insertion of the “Lambert brane” at point $x$ yields the modified condition (E.15) which is the analog of (E.9).

Notice that there are several other precise examples relating non-critical strings and topological strings. In particular the $c = 1$ string is described as the topological string on a conifold [89] and minimal models coupled to 2d gravity were argued to be described via a topological string on a particular Calabi-Yau manifold. The account of the gravity descendants in the minimal model corresponds to the adding of the B-branes into the B model geometry and insertion points are the open moduli.

The Lambert function has a finite radius of convergence $x_0 = \frac{1}{e}$ in the series representation

$$W(x) = \sum_{n=1} n^{n-1} \frac{x^n}{n!}.$$

The inspection of the radius of convergence in our case yields the condition for the area $A$ in the critical regime

$$A > \frac{1}{M \log M} \qquad (E.16)$$

which implies the validity of the approximation in this regime only. The meaning of a possible transition at this radius deserves special study.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. SCOAP³ supports the goals of the International Year of Basic Sciences for Sustainable Development.

## References

[1] V.A. Kazakov, *Bilocal Regularization of Models of Random Surfaces*, Phys. Lett. B 150 (1985) 282 [INSPIRE].

[2] F. David, *Planar Diagrams, Two-Dimensional Lattice Gravity and Surface Models*, Nucl. Phys. B 257 (1985) 45 [INSPIRE].
[3] V.A. Kazakov, A.A. Migdal and I.K. Kostov, *Critical Properties of Randomly Triangulated Planar Random Surfaces*, Phys. Lett. B 157 (1985) 295 [INSPIRE].

[4] V.A. Kazakov and A.A. Migdal, *Recent Progress in the Theory of Noncritical Strings*, Nucl. Phys. B 311 (1988) 171 [INSPIRE].

[5] P. Di Francesco, P.H. Ginsparg and J. Zinn-Justin, *2-D Gravity and random matrices*, Phys. Rept. 254 (1995) 1 [hep-th/9306153] [INSPIRE].

[6] E. Brezin and S.R. Wadia, *The Large N expansion in quantum field theory and statistical physics: From spin systems to two-dimensional gravity*, (1994) [INSPIRE].

[7] D. Anninos and B. Mühlmann, *Notes on matrix models (matrix musings)*, J. Stat. Mech. 2008 (2020) 083109 [arXiv:2004.01171] [INSPIRE].

[8] V.A. Kazakov, *Ising model on a dynamical planar random lattice: Exact solution*, Phys. Lett. A 119 (1986) 140 [INSPIRE].

[9] D.V. Boulatov and V.A. Kazakov, *The Ising Model on Random Planar Lattice: The Structure of Phase Transition and the Exact Critical Exponents*, Phys. Lett. B 186 (1987) 379 [INSPIRE].

[10] V.A. Kazakov, *Exactly solvable Potts models, bond- and tree-like percolation on dynamical (random) planar lattice*, Nucl. Phys. B Proc. Suppl. 4 (1988) 93 [INSPIRE].

[11] I.K. Kostov, *O(n) Vector Model on a Planar Random Lattice: Spectrum of Anomalous Dimensions*, Mod. Phys. Lett. A 4 (1989) 217 [INSPIRE].

[12] V.A. Kazakov, *The Appearance of Matter Fields from Quantum Fluctuations of 2D Gravity*, Mod. Phys. Lett. A 4 (1989) 2125 [INSPIRE].

[13] E. Brezin and V.A. Kazakov, *Exactly Solvable Field Theories of Closed Strings*, Phys. Lett. B 236 (1990) 144 [INSPIRE].

[14] M.R. Douglas and S.H. Shenker, *Strings in Less Than One-Dimension*, Nucl. Phys. B 335 (1990) 635 [INSPIRE].

[15] D.J. Gross and A.A. Migdal, *Nonperturbative Two-Dimensional Quantum Gravity*, Phys. Rev. Lett. 64 (1990) 127 [INSPIRE].

[16] M.R. Douglas, *Strings in Less Than One-dimension and the Generalized $K^-D^-V$ Hierarchies*, Phys. Lett. B 238 (1990) 176 [INSPIRE].

[17] E. Brezin, V.A. Kazakov and A.B. Zamolodchikov, *Scaling Violation in a Field Theory of Closed Strings in One Physical Dimension*, Nucl. Phys. B 338 (1990) 673 [INSPIRE].

[18] D.J. Gross and N. Miljkovic, *A Nonperturbative Solution of $D = 1$ String Theory*, Phys. Lett. B 238 (1990) 217 [INSPIRE].

[19] P.H. Ginsparg and J. Zinn-Justin, *2-d gravity + 1-d matter*, Phys. Lett. B 240 (1990) 333 [INSPIRE].

[20] G. Parisi, *On the One-dimensional Discretized String*, Phys. Lett. B 238 (1990) 209 [INSPIRE].

[21] E. Brezin, C. Itzykson, G. Parisi and J.B. Zuber, *Planar Diagrams*, Commun. Math. Phys. 59 (1978) 35 [INSPIRE].

[22] D.J. Gross and I.R. Klebanov, *Vortices and the nonsinglet sector of the $c = 1$ matrix model*, Nucl. Phys. B 354 (1991) 459 [INSPIRE].

[23] D. Boulatov and V. Kazakov, *Vortex anti-vortex sector of one-dimensional string theory via the upside down matrix oscillator*, Nucl. Phys. B Proc. Suppl. 25 (1992) 38 [INSPIRE].
[24] D. Boulatov and V. Kazakov, *One-dimensional string theory with vortices as the upside down matrix oscillator*, Int. J. Mod. Phys. A 8 (1993) 809 [hep-th/0012228] [INSPIRE].

[25] V. Kazakov, I. Kostov and D. Kutasov, *A Matrix Model for the 2d Black Hole*, PoS tmr2000 (2000) 026 [INSPIRE].

[26] M.R. Douglas, I.R. Klebanov, D. Kutasov, J.M. Maldacena, E.J. Martinez and N. Seiberg, *A New hat for the $c=1$ matrix model*, in *From Fields to Strings: Circumnavigating Theoretical Physics: A Conference in Tribute to Ian Kogan*, pp. 1758–1827 (2003).

[27] V.K. Kazakov, *Percolation on a Fractal With the Statistics of Planar Feynman Graphs: Exact Solution*, Mod. Phys. Lett. A 4 (1989) 1691 [INSPIRE].

[28] J.-M. Daul, *Q states Potts model on a random planar lattice*, hep-th/9502014 [LPTENS-94-] [INSPIRE].

[29] I.K. Kostov, *Random surfaces, solvable lattice models and discrete quantum gravity in two dimensions*, Nucl. Phys. B Proc. Suppl. 10 (1989) 295.

[30] I.K. Kostov and M.L. Mehta, *Random Surfaces of Arbitrary Genus: Exact Results for $D = 0$ and $-2$ Dimensions*, Phys. Lett. B 189 (1987) 118 [INSPIRE].

[31] I.K. Kostov, *Strings with discrete target space*, Nucl. Phys. B 376 (1992) 539 [hep-th/9112059] [INSPIRE].

[32] I.K. Kostov and M. Staudacher, *Multicritical phases of the $O(n)$ model on a random lattice*, Nucl. Phys. B 384 (1992) 459 [hep-th/9203030] [INSPIRE].

[33] I.K. Kostov, *Thermal flow in the gravitational $O(n)$ model*, Bulg. J. Phys. 33 (2006) 297 [hep-th/0602075] [INSPIRE].

[34] B. Eynard and J. Zinn-Justin, *The $O(n)$ model on a random surface: Critical points and large order behavior*, Nucl. Phys. B 386 (1992) 558 [hep-th/9204082] [INSPIRE].

[35] B. Eynard and C. Kristjansen, *Exact solution of the $O(n)$ model on a random lattice*, Nucl. Phys. B 455 (1995) 577 [hep-th/9506193] [INSPIRE].

[36] B. Eynard and C. Kristjansen, *More on the exact solution of the $O(n)$ model on a random lattice and an investigation of the case $|n| > 2$*, Nucl. Phys. B 466 (1996) 463 [hep-th/9512052] [INSPIRE].

[37] P. Saad, S.H. Shenker and D. Stanford, *JT gravity as a matrix integral*, arXiv:1903.11115 [INSPIRE].

[38] D.L. Jafferis, D.K. Kolchmeyer, B. Mukhametzhanov and J. Sonner, *JT gravity with matter, generalized ETH, and Random Matrices*, arXiv:2209.02131 [INSPIRE].

[39] A.A. Belavin and A.B. Zamolodchikov, *On Correlation Numbers in 2D Minimal Gravity and Matrix Models*, J. Phys. A 42 (2009) 304004 [arXiv:0811.0450] [INSPIRE].

[40] A. Belavin, B. Dubrovin and B. Mukhametzhanov, *Minimal Liouville Gravity correlation numbers from Douglas string equation*, JHEP 01 (2014) 156 [arXiv:1310.5659] [INSPIRE].

[41] A.M. Polyakov, *Quantum Geometry of Bosonic Strings*, Phys. Lett. B 103 (1981) 207 [INSPIRE].

[42] V.G. Knizhnik, A.M. Polyakov and A.B. Zamolodchikov, *Fractal Structure of 2D Quantum Gravity*, Mod. Phys. Lett. A 3 (1988) 819 [INSPIRE].
[43] F. David, *Conformal Field Theories Coupled to 2D Gravity in the Conformal Gauge*, Mod. Phys. Lett. A 3 (1988) 1651 [INSPIRE].

[44] J. Distler and H. Kawai, *Conformal Field Theory and 2D Quantum Gravity*, Nucl. Phys. B 321 (1989) 509 [INSPIRE].

[45] G.W. Moore, N. Seiberg and M. Staudacher, *From loops to states in 2-D quantum gravity*, Nucl. Phys. B 362 (1991) 665 [INSPIRE].

[46] M. Staudacher, *The Yang-lee Edge Singularity on a Dynamical Planar Random Surface*, Nucl. Phys. B 336 (1990) 349 [INSPIRE].

[47] D.V. Boulatov, V.A. Kazakov, I.K. Kostov and A.A. Migdal, *Analytical and Numerical Study of the Model of Dynamically Triangulated Random Surfaces*, Nucl. Phys. B 275 (1986) 641 [INSPIRE].

[48] I.R. Klebanov and R.B. Wilkinson, *Matrix model in two-dimensions and its effective field theory*, Phys. Lett. B 251 (1990) 379 [INSPIRE].

[49] I.R. Klebanov and R.B. Wilkinson, *Critical potentials and correlation functions in the minus two-dimensional matrix model*, Nucl. Phys. B 354 (1991) 475 [INSPIRE].

[50] J.D. Edwards and I.R. Klebanov, *Macroscopic boundaries and the wave function of the universe in the $c = -2$ matrix model*, Mod. Phys. Lett. A 6 (1991) 2901 [INSPIRE].

[51] B. Duplantier and F. David, *Exact partition functions and correlation functions of multiple hamiltonian walks on the Manhattan lattice*, J. Statist. Phys. 51 (1988) 327.

[52] Y. Ishimoto and A.B. Zamolodchikov, *Massive Majorana fermion coupled to 2D gravity and random lattice Ising model*, Theor. Math. Phys. 147 (2006) 755 [INSPIRE].

[53] S. Caracciolo et al., *Fermionic field theory for trees and forests*, Phys. Rev. Lett. 93 (2004) 080601 [cond-mat/0403271] [INSPIRE].

[54] S. Caracciolo and A. Sportiello, *Spanning Forests on Random Planar Lattices*, J. Statist. Phys. 135 (2009) 1063 [arXiv:0903.4432] [INSPIRE].

[55] R. Bondesan, S. Caracciolo and A. Sportiello, *Critical Behaviour of Spanning Forests on Random Planar Graphs*, J. Phys. A 50 (2017) 074003 [arXiv:1608.02916] [INSPIRE].

[56] K. Kelmans and V.M. Chelnokov, *A certain polynomial of a graph and graphs with an extremal number of trees*, J. Combin. Theory 16 (1974) 074003.

[57] G. Parisi and N. Sourlas, *Random Magnetic Fields, Supersymmetry and Negative Dimensions*, Phys. Rev. Lett. 43 (1979) 744.

[58] J.M. Daul, V.A. Kazakov and I.K. Kostov, *Rational theories of 2-D gravity from the two matrix model*, Nucl. Phys. B 409 (1993) 311 [hep-th/9303093] [INSPIRE].

[59] F. David, *Loop Equations and Nonperturbative Effects in Two-dimensional Quantum Gravity*, Mod. Phys. Lett. A 5 (1990) 1019 [INSPIRE].

[60] B. Eynard and J. Zinn-Justin, *Large order behavior of 2-D gravity coupled to $d < 1$ matter*, Phys. Lett. B 302 (1993) 396 [hep-th/9301004] [INSPIRE].

[61] V.A. Kazakov and I.K. Kostov, *Instantons in noncritical strings from the two matrix model*, in the proceedings of *From Fields to Strings: Circumnavigating Theoretical Physics: A Conference in Tribute to Ian Kogan*, (2004), p. 1864–1894 [DOI:10.1142/9789812775344_0045] [hep-th/0403152] [INSPIRE].
[62] K. Okuyama, Quenched free energy from spacetime D-branes, JHEP 03 (2021) 073
[arXiv:2101.05990] [INSPIRE].

[63] M. Aganagic et al., Topological strings and integrable hierarchies, Commun. Math. Phys. 261 (2006) 451 [hep-th/0312085] [INSPIRE].

[64] A. Marshakov and N. Nekrasov, Extended Seiberg-Witten Theory and Integrable Hierarchy, JHEP 01 (2007) 104 [hep-th/0612019] [INSPIRE].

[65] V. Bouchard and M. Mariño, Hurwitz numbers, matrix models and enumerative geometry, Proc. Symp. Pure Math. 78 (2008) 263 [arXiv:0709.1458] [INSPIRE].

[66] H. Kesten, Symmetric random walks on groups, Trans. Am. Math. Soc. 92 (1959) 336.

[67] B.D. McKay, The expected eigenvalue distribution of a large regular graph, Linear Algebra Appl. 40 (1981) 203.

[68] B.L. Altshuler, Y. Gefen, A. Kamenev and L.S. Levitov, Quasiparticle Lifetime in a Finite System: A Nonperturbative Approach, Phys. Rev. Lett. 78 (1997) 2803 [cond-mat/9609132] [INSPIRE].

[69] K.S. Tikhonov and A.D. Mirlin, From Anderson localization on random regular graphs to many-body localization, Annals Phys. 435 (2021) 168525 [arXiv:2102.05930].

[70] S. Moudgalya, B.A. Bernevig and N. Regnault, Quantum many-body scars and Hilbert space fragmentation: a review of exact results, Rept. Prog. Phys. 85 (2022) 086501 [arXiv:2109.00548] [INSPIRE].

[71] F. David, Randomly Triangulated Surfaces in Two-dimensions, Phys. Lett. B 159 (1985) 303 [INSPIRE].

[72] A. Abdesselam, The Grassmann-Berezin calculus and theorems of the matrix-tree type, Adv. Appl. Math. 33 (2004) 51.

[73] V.A. Kazakov, External matrix field problem and new multicriticities in (two)-dimensional random surfaces, Nucl. Phys. B 354 (1991) 614 [INSPIRE].

[74] V.A. Kazakov, M. Staudacher and T. Wynter, Character expansion methods for matrix models of dually weighted graphs, Commun. Math. Phys. 177 (1996) 451 [hep-th/9502132] [INSPIRE].

[75] V.A. Kazakov, M. Staudacher and T. Wynter, Almost flat planar diagrams, Commun. Math. Phys. 179 (1996) 235 [hep-th/9506174] [INSPIRE].

[76] V.A. Kazakov, M. Staudacher and T. Wynter, Exact solution of discrete two-dimensional $R^{**2}$ gravity, Nucl. Phys. B 471 (1996) 309 [hep-th/9601069] [INSPIRE].

[77] V. Kazakov and F. Levkovich-Maslyuk, Disc partition function of $2d R^2$ gravity from DWG matrix model, JHEP 01 (2022) 190 [arXiv:2110.10104] [INSPIRE].

[78] I.K. Kostov, M. Staudacher and T. Wynter, Complex matrix models and statistics of branched coverings of 2-D surfaces, Commun. Math. Phys. 191 (1998) 283 [hep-th/9703189] [INSPIRE].

[79] M.J. Stephen, Percolation problems and the Potts model, Phys. Lett. A 56 (1976) 149.

[80] V. Avetisov et al., Eigenvalue tunneling and decay of quenched random network, Phys. Rev. E 94 (2016) 062313 [arXiv:1607.03871] [INSPIRE].

[81] C. Kelly, C.A. Trugenberger and F. Biancalana, Self-Assembly of Geometric Space from Random Graphs, Class. Quant. Grav. 36 (2019) 125012 [arXiv:1901.09870] [INSPIRE].
[82] O. Valba and A. Gorsky, *interacting thermofield doubles and critical behavior in random regular graphs*, Phys. Rev. D **103** (2021) 106013 [arXiv:2101.04072] [INSPIRE].

[83] V. Avetisov, A. Gorsky, S. Nechaev and O. Valba, *Localization and non-ergodicity in clustered random networks*, J. Complex Netw. **8** (2020) cnz026.

[84] F.L. Metz, G. Parisi and L. Leuzzi, *Finite size correction to the spectrum of regular random graphs: an analytical solution*, arXiv:1403.2582 [DOI:10.1103/PhysRevE.90.052109].

[85] J.J.M. Verbaarschot and T. Wettig, *Random matrix theory and chiral symmetry in QCD*, Ann. Rev. Nucl. Part. Sci. **50** (2000) 343 [hep-ph/0003017] [INSPIRE].

[86] I. Pak, A. Kelmans and A. Postnikov, *Tree and forest volumes of graphs*, DIMACS Tech. Rep. 2000-03 **135** (2000) 1063.

[87] D. Marolf and H. Maxfield, *Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information*, JHEP **08** (2020) 044 [arXiv:2002.08950] [INSPIRE].

[88] R. Dijkgraaf and C. Vafa, *Matrix models, topological strings, and supersymmetric gauge theories*, Nucl. Phys. B **644** (2002) 3 [hep-th/0206255] [INSPIRE].

[89] D. Ghoshal and C. Vafa, *C = 1 string as the topological theory of the conifold*, Nucl. Phys. B **453** (1995) 121 [hep-th/9506122] [INSPIRE].