# Disc partition function of 2d R² gravity from DWG matrix model

Vladimir Kazakov, Fedor Levkovich-Maslyuk

► To cite this version:

Vladimir Kazakov, Fedor Levkovich-Maslyuk. Disc partition function of 2d R² gravity from DWG matrix model. Journal of High Energy Physics, 2022, 2022 (190), pp.1-39. 10.1007/JHEP01(2022)190. hal-03413196

HAL Id: hal-03413196

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

Submitted on 22 Aug 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
# Disc partition function of 2d $R^2$ gravity from DWG matrix model

Vladimir Kazakov$^{a,b}$ and Fedor Levkovich-Maslyuk$^{c,1}$

$^a$Laboratoire de Physique de l'École Normale Supérieure, CNRS, Université PSL, Sorbonne Universités, 24 rue Lhomond, 75005 Paris, France

$^b$Interdisciplinary Scientific Center Poncelet (CNRS UMI 2615), 119002 Moscow, Russia

$^c$Université Paris Saclay, CNRS, CEA, Institut de physique théorique, 91191, Gif-sur-Yvette, France

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

**ABSTRACT:** We compute the sum over flat surfaces of disc topology with arbitrary number of conical singularities. To that end, we explore and generalize a specific case of the matrix model of dually weighted graphs (DWG) proposed and solved by one of the authors, M. Staudacher and Th. Wynter. Namely, we compute the sum over quadrangulations of the disc with certain boundary conditions, with parameters controlling the number of squares (area), the length of the boundary and the coordination numbers of vertices. The vertices introduce conical defects with angle deficit given by a multiple of $\pi$, corresponding to positive, zero or negative curvature. Our results interpolate between the well-known 2d quantum gravity solution for the disc with fluctuating 2d metric and the regime of “almost flat” surfaces with all the negative curvature concentrated on the boundary. We also speculate on possible ways to study the fluctuating 2d geometry with $AdS_2$ background instead of the flat one.

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

ARXIV ePRINT: 2110.10104

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

1 Introduction 1

2 The DWG matrix model via character expansion 5

2.1 Dually weighted Feynman graphs 5

2.2 Weighted quandrangulations ("R² quantum gravity") 8

2.2.1 Natural observables and combinatorics 8

3 Character expansion and planar limit for DWG model 9

3.1 DWG partition function as a sum over SU(N) representations 9

3.2 Characters in the large N limit 11

3.3 Saddle point equation and general RH problem for DWG 12

4 The R² quantum gravity solution 14

4.1 Notation and key relations for the case of even potentials 14

4.2 Explicit solution for the density 16

4.3 Fixing parameters of the solution 17

4.4 Simplification and computation of G 18

4.5 Scaling properties of G as a function of ϵ 19

5 Critical regimes of large area and boundary of the quadrangulated disc 20

5.1 Solution with exponential precision for near-flat regime 21

5.1.1 Expanding the parameters 21

5.1.2 Computing the integral for G(h) 22

5.2 Critical regime and flattening limit 24

5.3 Resolvent for Tr(A₄M̃)ⁿ correlators in the flattening limit 25

5.3.1 Checks of the result 27

5.3.2 Pure gravity and almost flat limits 28

5.4 Resolvent for TrMⁿ correlators 30

5.4.1 Limits and properties of the result 32

5.4.2 Comparison of two resolvents 33

6 Conclusions and prospects 34

A Derivation of the partition function and details on characters 36

B Mathematica notation for elliptic functions 36

C Details on the exponential approximation 37
# 1 Introduction

Matrix integrals at large size $N$ of the matrices are a formidable tool for counting various types of planar graphs, i.e. “fat” graphs which have a certain two dimensional topology [1–4]. The topological property of the $1/N$ expansion in (multi)-matrix integrals and in matrix field theories, applied to Feynman graphs of asymptotically large sizes, is the key element for the study of quantum fluctuations of 2d geometries — 2d quantum gravity — with various embedding, as well as the statistical mechanics of spins on planar graphs [5–14]. The famous AdS/CFT correspondence [15–17] is also based on the deep relation between planar graphs and string worldsheets. More recently, the topological expansion in a specific one-matrix model was proposed for the description of Jackiw-Teitelboim (JT) gravity [18–23].

Usually, matrix models and field theories have couplings in their actions which control the numbers of vertices of certain types in the corresponding planar graphs, but not the numbers of faces of given types. For instance, planar Feynman graphs in the matrix scalar field theory with $\lambda \text{tr} \phi^4$ interaction have weights $\lambda^n$, where $n$ is the number of quartic vertices. But the numbers of faces of given size $n^* = 1, 2, 3, \dots$ are arbitrary (up to the constraint imposed by the Euler theorem) and not weighted by independent couplings.

There exists an elegant way to modify a matrix quantum field theory so that also the faces of the Feynman graphs (vertices of the dual graphs) will be weighted with dual couplings attached to the faces of a given order.<fnref n="1" /> For example, for the scalar theory one can modify to that end the action as follows:

$$Z(t, t^*) = \int \mathcal{D}^{N^2} \phi(x) \exp N \text{tr} \left( -\frac{1}{2}(\nabla\phi)^2 - \frac{m^2}{2}\phi^2 + \sum_{q=1}^{Q+1} \frac{1}{q} t_q (A\phi)^q \right), \quad (1.1)$$

where we define the dual couplings as $t_n^* = \frac{1}{N} \text{tr} A^n$. Then the partition function is given by a perturbative expansion over planar dually weighted Feynman graphs (DWG)

$$\log Z(t, t^*) = \sum_G N^{2-2g_G} \prod_{v_{q_*}^*, v_q \in G} w_G(\#v_{q_*}^*, \#v_q) \quad t_{q_*}^{\#v_{q_*}} t_q^{\#v_q}, \quad (1.2)$$

where the sum goes over all planar Feynman graphs $G$, $g_G$ is the genus of the graph, $v_{q_*}^*$ and $v_q$ are the numbers of vertices with $q^*$ and $q$ neighbours on the dual graph and the original graph, respectively, while $\#v_{q_*}^*, \#v_q$ are the numbers of such vertices in a given graph. The weight $w_G(\#v_{q_*}^*, \#v_q)$ is the contribution of graphs with fixed $\#v_{q_*}^*, \#v_q$ and $g_G$ (the result of computation of the corresponding Feynman integrals). This identification of dual couplings and the expansion (1.2) can be derived from the simple fact that the vertices in the planar graph expansion corresponding to the terms $t_q \text{tr} (A\phi)^q$ in the action can be presented as shown on figure 1. It is clear that this results in contracting the product of $q_*$ matrices $A$ in the index loop around each face with $q_*$ edges, giving the factor $t_q^* = \frac{1}{N} \text{tr} A^q$, in the normalization appropriate for the ‘t Hooft limit.<fnref n="2" />

<footnote n="1">¹And factors depending on the type of interactions.</footnote>

<footnote n="2">²Curiously, this idea of introducing dual couplings in matrix models can be used for constructing a quantum field theory in any dimension at finite $N_c$ as a zero-dimensional matrix model in a specific large $N$ limit [24] ($N_c$ and $N$ here are different parameters).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** A dually weighted Feynman graph. Solid double lines show the propagators, while the dashed lines show the dual graph (whose vertices correspond to faces of the original graph and vice versa). At each *n*-vertex we have *n* copies of the matrix *A* which after contraction give factors ∝ tr *A*^q corresponding to each *q*-face.</figcaption></figure>

In this work, we will continue the study of DWG and the matrix model formulated in various forms in [25–27] and solved in the planar limit in the series of works [27–30]. It is the zero-dimensional version of matrix QFT (1.1). It represents a natural generalization of the one-matrix model with a general potential, like in [2, 10] where the introduction of the *A* matrix in the interaction potential allows one to control not only the order of the vertices of planar graphs (coordination number, or the number of nearest neighbors), but also the order of their faces (coordination numbers of vertices of the dual graph). So the DWG matrix model has two sets of couplings: original and dual. The partition function for graphs with spherical topology is then given by (1.2) with all the weights $w_G(\#v_{q_*}^*, \#v_q) \equiv 1$.

A particular case of this zero-dimensional DWG model was considered in the work [29]. For this case, the graphs are quadrangulations and the vertices $v_{2n}$ can have arbitrary even orders $2n$, i.e. with $2n$ neighboring square faces, where $n = 1, 2, 3, \dots$, with the weights $t_{2n}$ chosen so that they suppress exponentially the high orders $n \gg 1$ of vertices. The partition function of such quadrangulations with spherical topology has been computed in [29].

In the current paper, we will study the same system of random quadrangulations but for the disc topology where we have a single boundary of even length $L = 2q$, as shown on figure 2. As will be explained below, in order to make the problem analytically solvable, the boundary can be introduced in three different ways. One particularly convenient way corresponds to cutting out from the closed quadrangulation a vertex with $L$ adjacent squares. The precise definition of that kind of disc partition function for this model is:

$$Z(P, \lambda, \beta) = \sum_{L=2}^{\infty} L P^{-L} \sum_{Q_L} (\#v_L) \lambda^{F+1} \beta^{2(\#v_2 - \delta_{L,2})}, \quad (1.3)$$

where the first sum goes over the length $L$ of the disc boundary (always even) and the second sum goes over all disc quadrangulations $Q_L$. Here $\#v_k$ denotes the number of vertices of order $k = 2, 4, 6, \dots$ which correspond to the angle deficits<fnref n="3" /> $\alpha_k = \pi(2 - k/2)$.

<footnote n="3">³To define these angle deficits, we assume that all faces are made of equal squares with the angles $\pi/2$ at each vertex.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic image of an abstract manifold with disc topology, with the boundary shown by the wiggled red line. It is flat everywhere except conical singularities. In this example we have two of such singularities with positive angle deficit (positive curvature defects) and the 3rd, in the middle, with negative angle deficit (negative curvature defect).</figcaption></figure>

In the quantum Regge gravity picture [5–7, 9, 31], $\alpha_2$ corresponds to the single type of positive curvatures, $\alpha_4$ introduces no curvature and $\alpha_{n>4}$ correspond to negative curvature insertions. Furthermore, $P$ is the length fugacity, $\lambda$ is the area fugacity (the area $F$ is defined as the number of square faces forming the quadrangulation), and $\beta$ is the fugacity of positive curvature (type $\alpha_2$) insertions. Notice that due to the Euler theorem for the disc we have $\#v_2 = \sum_{k \ge 3} \#v_k + \text{const}$, so that the overall positive curvature is, up to an additive constant of the order 1, equal to the overall negative curvature, so that $\beta$ controls the overall scale (or the absolute value) of the total curvature. That allows us to interpolate between the regimes of large curvature fluctuations and the small curvature regime dominated by almost flat (AF) configurations which we describe in more detail later.

This model was solved exactly in [29] and it was called there the discretised $2d R^2$ quantum gravity, since the possibility to control independently the positive and negative curvature amounts, in the continuous limit, to introducing irrelevant perturbations into the QG action, such as $R^2$ (where $R$ is the Gauss curvature), providing a “flattening” effect on the geometry. In the continuum limit, it should lead to the 2d QG action for the disc topology

$$S_{R^2}^{M} = -\frac{1}{2} \int_{M} \sqrt{g} (\tilde{\lambda} + \tilde{\beta} R^2) - \int_{\partial M} \sqrt{h} (\tilde{\gamma} + K). \qquad (1.4)$$

Here we denoted by $\tilde{\lambda}$, $\tilde{\beta}$, $\tilde{\gamma}$ the bulk cosmological constant, the squared curvature coupling and the boundary cosmological constant, respectively — the renormalized continuous analogs of the lattice fugacities $\lambda$, $\beta$, $P$ in (1.3). We also used the notation $R$, $K$, $g$, and $h$ for the bulk and boundary gaussian curvatures and metrics.

The study of the exact sphere partition function of this DWG model in the continuous limit (i.e. for large quadrangulations with no boundary) showed that there are two smoothly connected regimes: almost flat manifolds (with very few conical singularities) and the pure $2d$ gravity regime of [5–7, 9, 31]. Since the $R^2$ coupling $\tilde{\beta}$ has the dimension
of length, the former regime occurs for sufficiently small areas of manifolds, whereas the latter happens for large enough areas (the scale is set by $\tilde{\beta}$).

In this paper, we study the quadrangulations having disc topology, with the boundary of two different types described in the next section. We will keep the same continuous limit of quadrangulations with large area, but now also with a boundary of a fixed size (made out of a certain number of edges). In the limit of large area, our explicit results interpolate between the universal pure quantum gravity regime reproducing the well known disc partition function [10], and the regime of “almost flat” (AF) random surfaces [28]. The former regime assumes also the limit of a long boundary. In the latter regime, the surfaces are flat, consisting predominantly of the vertices with four neighboring plaquettes that add no curvature, everywhere except rare insertions of conical singularities, with angle deficits multiple of $\pi$, with specific exponential weights suppressing large angle deficits. Thus in the continuous limit we simply sum up over the surfaces flat everywhere except a collection, or “gas” of those conical defects. An intuitive image of such abstract flat manifold with curvature defects is presented on figure 2.

For the disc partition function of quadrangulations with two different types of the boundary we found very similar, but not exactly identical expressions in the large area limit, given by eqs. (5.43) and (5.72). This similarity demonstrates certain universality features of this disc partition function, which however shows also a dependence on the boundary conditions. Both solutions interpolate between the well known pure 2d QG regime and almost flat regime for the disc partition function. In the former one, the results are completely universal, as expected.

For the disc partition function in this nearly flat limit we found that there is no phase transition between these two regimes, thus extending and supporting the findings of [29] that were obtained at the level of the sphere partition function. While so far we did not observe signs of any more exotic phases in this DWG matrix model, such as the JT gravity regime [18], in the conclusions we speculate on the possibility to find them for more sophisticated observables. We will also discuss potential ways to produce from our model a discretized fluctuating geometry with the metric dominated by $AdS_2$ background.

The paper is organized as follows. In section 2 we present the basic properties of the general DWG model as well as defining the main observables we will study. In subsection 2.2 there we define and discuss the particular $R^2$ QG matrix model studied in [29]. Then in section 3 we return to the general DWG model and present the main steps towards its solution by character expansion, mainly reviewing the results of [27, 28]. From section 4 onward we specialize to the $R^2$ QG model of [29]. In section 4 we review its solution obtained in [29]. Then in section 5 we present our main results. We compute a number of observables and explore their critical behavior. These are the resolvent for a particular kind of correlators with exponential accuracy, its continuum limit corresponding to nearly flat manifolds, and a similar limit for the resolvent describing the disc partition function with a different boundary. We present conclusions in section 6, while the appendices contain a number of technical details.
## 2 The DWG matrix model via character expansion

The DWG matrix model is defined by [25, 27]

$$Z(t, t^*) = \int DM \exp \left[ -\frac{N}{2} \text{tr} M^2 + N \sum_{q=1}^{Q+1} \frac{1}{q} t_q \text{tr} (MA)^q \right], \quad (2.1)$$

where $M$ is an $N \times N$ Hermitian matrix integrated with the SU($N$) invariant measure $DM$ and $A$ is a fixed $N \times N$ Hermitian matrix which can be always chosen diagonal without loss of generality. The DWG matrix model (2.1) cannot be solved directly by the diagonalization of the matrix $M$. However, it turns out to be possible, in the large $N$ limit, to solve it (in principle) by the character expansion methods worked out in [27–29].$^4$ In this section we will briefly review the main steps of this approach, referring the interested reader to those papers for the details.

First we will review the representation of $Z$ in terms of dually weighted graphs (DWG). Then we will demonstrate the character expansion method which allows one to rewrite the model in terms of the sum over $N$ highest weights of GL($N$) representations (Young tableaux). At large $N$ the sum is dominated by a single large smooth Young tableau, and the computation of characters, as well as the saddle point equation, are reduced to a certain Riemann-Hilbert (RH) problem. The solution of this RH problem is directly related to three different types of disc partition functions of DWG, which are the main object of interest of this work.

### 2.1 Dually weighted Feynman graphs

The perturbative expansion w.r.t. couplings $t_q$ is represented by planar Feynman diagrams of the type shown on figure 3. Their elements are the directed double-line propagators$^5$ equal to $\frac{1}{N}\delta_i^j\delta_k^l$ (the indices run from 1 to $N$ and are conserved along each line) and the vertices which give $Nt_q$ ($A_{i_1j_1}A_{i_2j_2}\dots A_{i_qj_q}$). It is easy to see that at each face of order $p$ of the lattice, after contraction of indices along each index line, there will appear the factor

$$t_p^* = \frac{1}{N} \text{tr} A^p. \qquad (2.2)$$

Hence the partition function of this DWG matrix model has the following expansion in terms of DWG of a distinct topology [27]:

$$\log Z(t, t^*) = \sum_G N^{2-2g_G} \prod_{v_q^*, v_q \in G} t_q^{*\#v_q^*} t_q^{\#v_q}, \quad (2.3)$$

where the sum goes over all Feynman graphs $G$, $g_G$ is the genus of the graph $G$, and the product goes over all the vertices and faces of the graph. We indicate by $v_q$ vertices of order

<footnote>$^4$See also the recent work [32] the methods of which provide another way to the solution and perhaps to a generalization of the DWG models.</footnote>

<footnote>$^5$“Directed” means here that two lines of the same propagator should be marked by opposite arrows, symbolising the covariant and contravariant indices. Since the direction in each index loop is conserved it ensures that each planar graph is orientable.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Dually weighted graph (figure from [27]). We show the propagators as double lines. The couplings $t_n$ appear at n-vertices, while $t_n^*$ correspond to n-faces of the graph, corresponding to the faces of the dual graph depicted by dotted double-lines.</figcaption></figure>

$q$ and by $v_q^*$ faces of order $q$, while $\#v_q$ and $\#v_q^*$ denote the number of these vertices and faces in a graph from the sum. Let us also note that we can view the $q$-faces as $q$-vertices on the dual graph $G^*$ (whose faces are the vertices of the original graph and vice versa).

From (2.3) we see that $Z$ is symmetric in the two sets of couplings $t_p \leftrightarrow t_p^*$, as there is a one-to-one correspondence between a graph $G$ and its dual $G^*$. It is natural to introduce a ‘dual’ matrix model with the same partition function but with Feynman diagrams corresponding to these dual graphs. It corresponds to interchanging all the couplings $t_p$ and $t_p^*$ in the original model and thus is defined by

$$Z = \int \mathcal{D}\tilde{M} \exp \left[ -\frac{N}{2} \text{tr} \, \tilde{M}^2 + N \sum_{q=1}^{Q+1} \frac{1}{q} t_q^* \text{tr} \, (\tilde{M}B)^q \right], \quad (2.4)$$

where we have introduced a constant matrix $B$ which is the counterpart of $A$ satisfying a relation similar to (2.2),

$$t_p = \frac{1}{N} \text{tr} \, B^p. \qquad (2.5)$$

We have also labelled the matrix being integrated over as $\tilde{M}$ instead of $M$ in the original model. Thus we see that in the dual model the roles of $A$ and $B$ are swapped. It is clear that such a duality, as well as the very representation of (arbitrary) couplings in terms of traces of powers of these matrices can be in general possible only for $N \to \infty$, at each order of topological $1/N$ expansion. But later we will reformulate the model in terms of character expansion over $\text{GL}(N)$ representations (Young tableaux) where the characters will depend explicitly on independent couplings $t_n$ and $t_n^*$.
There exist three matrix averages, computable by the developed methods, representing the DWG partition functions with a disc topology with three distinct types of the boundary:

$$W_L = \frac{1}{N} \langle \mathrm{tr} (AM)^L \rangle, \qquad (2.6)$$

$$\tilde{W}_L = \frac{1}{N} \langle \mathrm{tr} (B\tilde{M})^L \rangle, \qquad (2.7)$$

$$W_L = \frac{1}{N} \langle \mathrm{tr} M^L \rangle = \frac{1}{N} \langle \mathrm{tr} (\tilde{M})^L \rangle, \qquad (2.8)$$

where the averages involving $\tilde{M}$ are computed in the dual model (2.4). The fact that the two types of correlators in (2.8) are equal may not be immediately obvious, but can be proven using the methods that we will discuss below.

The corresponding resolvents are

$$W(z) = \sum_{L \ge 0} z^{-L-1} W_L = \left\langle \frac{1}{N} \mathrm{tr} \frac{1}{z-AM} \right\rangle, \qquad (2.9)$$

$$\tilde{W}(z) = \sum_{L \ge 0} z^{-L-1} \tilde{W}_L = \left\langle \frac{1}{N} \mathrm{tr} \frac{1}{z-B\tilde{M}} \right\rangle, \qquad (2.10)$$

$$W(z) = \sum_{L \ge 0} z^{-L-1} W_L = \left\langle \frac{1}{N} \mathrm{tr} \frac{1}{z-M} \right\rangle = \left\langle \frac{1}{N} \mathrm{tr} \frac{1}{z-\tilde{M}} \right\rangle. \qquad (2.11)$$

To realize the boundaries for (2.6) and (2.7) as the insertion of an extra L-vertex into the graph G, we can simply differentiate the partition function w.r.t. the appropriate coupling:<fnref n="6" />

$$W_L = \frac{1}{N^2} \partial_{t_L} \log Z(t, t^*), \qquad (2.12)$$

$$\tilde{W}_L = \frac{1}{N^2} \partial_{t_L^*} \log Z(t, t^*). \qquad (2.13)$$

The corresponding disc partition functions have an obvious DWG picture. The first one (2.12) can be viewed as a sphere partition function with one marked vertex of the order L (divided by L, which takes away the overcounting due to the cyclic symmetry), for which the coupling $t_L$ is removed. The second one (2.13) can be also viewed as a sphere partition function, but now with one marked face of the order L (also divided by L) and its coupling $t_L^*$ removed. For the moments $W_L$ in the resolvent (2.11) the DWG interpretation is slightly more complicated: they can be also viewed as a sphere partition function with marked vertex of the order L but since in the matrix model formulation there are no matrices A around this vertex the adjacent faces do not contain it and are weighted each with the coupling $t_{n-1}$ for a face of order n. Obviously, all three resolvents (2.9)–(2.11) can be viewed as disc partition functions with such vertex or face removed.

Later we will explain that the quantities (2.9)–(2.11) can be computed directly from the saddle point equations for the DWG model.

<footnote n="6">⁶There is no similar formula for $W_L$ unless we generalize the matrix potential in 2.1 by adding there $\sum_{p \ge 3} \tau_p \mathrm{tr} M^p$. Such a model, containing the 3rd infinite set of couplings, is still solvable, in principle, by character expansion methods, and we would have $W_L = \frac{1}{N^2} \partial_{\tau_L} \log Z(t, t^*, \tau)$.</footnote>
## 2.2 Weighted quadrangulations (“$R^2$ quantum gravity”)

Before proceeding with key equations for the general DWG model, let us now introduce a particular case of the DWG model which we will focus on in this paper (starting from section 4). It was solved in [29] and used there to study the $R^2$ QG model. This model generates only quadrangulations and has special weights for the vertices, exponentially decaying with the order. This model will be of the main interest in this paper since one can solve it explicitly for the disc partition functions and thus it is convenient for the search of various continuous limits of fluctuating 2d manifolds.

The couplings of the model are given by

$$t_n = \delta_{n,4}, \quad t_n^* = T_n \qquad (2.14)$$

where $T_n$ are defined as

$$T_2 = \frac{\lambda\beta^2}{\epsilon} \quad \text{and} \quad T_{2q} = \lambda\epsilon^{q-2}, \quad q \ge 2. \qquad (2.15)$$

We see that they are all expressed in terms of 3 parameters: $\lambda, \beta$ and $\epsilon$. The partition function reads

$$Z = \int \mathcal{D}M \exp N \text{tr} \left( -\frac{1}{2} M^2 + \frac{1}{4} (AM)^4 \right), \qquad (2.16)$$

and for the dual matrix model we have

$$Z = \int \mathcal{D}\tilde{M} \exp N \text{tr} \left( -\frac{1}{2} \tilde{M}^2 + \sum_{q=1}^{\infty} \frac{T_{2q}}{2q} (A_4 \tilde{M})^{2q} \right) \qquad (2.17)$$

(as discussed above, these two partition functions are equal). Here the matrices $A$ and $A_4$ are chosen such that

$$\frac{1}{N} \text{tr} A_4^n = \delta_{n,4}, \quad \frac{1}{N} \text{tr} A^n = T_n. \qquad (2.18)$$

That also means that in our conventions we have

$$B = A_4. \qquad (2.19)$$

Notice that the couplings $T_n$ were denoted by $t_n$ in [29], which was a bit confusing since these are the couplings corresponding to faces, not vertices, in the original model (2.16). Thus we chose to introduce the new notation $T_n$ for them.

### 2.2.1 Natural observables and combinatorics

In our conventions it is the dual model (2.17), not the original model (2.16), that corresponds to graphs that are quadrangulations. Accordingly, we will study below the resolvent $\tilde{W}(z)$ generating the correlators $\langle \text{tr} (A_4 \tilde{M})^n \rangle$ in this dual model, as well as the resolvent $W(z)$ that generates the $\langle \text{tr} (\tilde{M})^n \rangle = \langle \text{tr} M^n \rangle$ correlators. Both observables correspond to graphs with disc topology built out of 4-vertex plaquettes.

Invoking simple combinatorics and Euler theorem, the partition function in this dual model (equal, however, to the original $Z$) can be represented as [29]

$$\log Z = \frac{\lambda^2 \beta^8}{\epsilon^4} \sum_G \lambda^F \beta^{2(\#v_2-4)}, \qquad (2.20)$$

where the sum is over connected graphs that are quadrangulations of the sphere.
We will be interested in particular in computing the resolvent $\tilde{W}$. The quantity $\tilde{W}_L = \frac{1}{N}\langle \text{tr}(A_4\tilde{M})^L \rangle$ is given by graphs with the shape of a quadrangulated disc with the boundary of (even) length $L$ corresponding to a (removed) vertex of order $L$. So, to get $\tilde{W}_L$, the weight of the vertex given by $T_L = \lambda\epsilon^{L/2-2}$ should be removed from the partition function (2.3). Explicitly, we see from (2.17) that

$$ \frac{1}{N}\langle \text{Tr}(\tilde{M}A_4)^L \rangle = L \frac{\partial \log Z}{\partial T_L}, \qquad (2.21) $$

and (2.3) implies that taking the derivative in $T_L$ in the sum over graphs simply produces a factor $\propto \#v_L$. As in our case from (2.15) we have $T_L = \lambda\epsilon^{L/2-2}\beta^{2\delta_{L,2}}$ (for even $L$), from the equations above we get

$$ \frac{1}{N}\langle \text{Tr}(\tilde{M}A_4)^L \rangle = \sum_{G_L} L\#v_L \lambda^{F+1} \epsilon^{-2-L/2} \beta^{2\#v_2-2\delta_{L,2}}, \qquad (2.22) $$

where the sum is over connected graphs with at least one $L$-vertex. This gives for the resolvent

$$ z\tilde{W}(z) = \sum_{L \ge 0} z^{-L} \frac{1}{N} \langle \text{Tr}(\tilde{M}A_4)^L \rangle = 1 + \frac{1}{\epsilon^2} \sum_{L \ge 2, L \text{ even}} (\epsilon^{1/2}z)^{-L} f_L(\lambda, \beta), \qquad (2.23) $$

where

$$ f_L(\lambda, \beta) = \lambda L \beta^{-2\delta_{L,2}} \sum_{G_L} (\#v_L) \lambda^F \beta^{2\#v_2}. \qquad (2.24) $$

We see that the dependence of the result on $\epsilon$ is very restricted, namely up to a $1/\epsilon^2$ factor it is a function of $\epsilon^{1/2}z$.

Via an argument of the same type one can obtain a similar combinatorial formula for the correlators $\langle \text{Tr}(MA)^L \rangle$. At the same time, for the $\langle \text{Tr}M^L \rangle$ correlator the combinatorics seems to be more complicated and we leave its exploration for the future.

Below in section 3 we continue describing the main results regarding the general DWG models. Then from section 4 we will specialize to the case of the $R^2$ QG model (2.16) we have just discussed here.

# 3 Character expansion and planar limit for DWG model

In this section, we will review the reduction of the DWG matrix model integral to the character expansion over Young tableaux. Then we will use the fact that the sum over Young tableaux goes over only $N$ highest weights and apply the saddle point approximation.

## 3.1 DWG partition function as a sum over SU(N) representations

The way to represent the partition function of the DWG matrix model is as follows [27]: first we expand the exponent of the second term in the action (2.1) w.r.t. the Schur characters $\chi_R[t^*]$ and $\text{GL}(N)$ characters $\chi_R[AM]$, where $R$ is a representation of $\text{GL}(N)$, then we use the orthogonality property of the characters to integrate over the angular variable
$\Omega \in \text{SU}(N)$ in the angular decomposition $M = \Omega^\dagger X \Omega$, where $X = \text{diag}\{x_1, x_2, \dots, x_N\}$ is the diagonal matrix of eigenvalues, and finally we perform the explicit Gaussian integral over $M$. The details are reviewed in the appendix A. In this way we obtain the following formula for the DWG partition function in terms of the sums over representations:<fnref n="7" />

$$Z(t, t^*) = c \sum_{\{h^e, h^o\}} \frac{\prod_i(h_i^e - 1)!! h_i^o!!}{\prod_{i,j}(h_i^e - h_j^o)} \chi_{\{h\}}[t] \chi_{\{h\}}[t^*]. \quad (3.1)$$

Here the sum goes over $\text{GL}(N)$ representations labeled by Young tableaux defined by the shifted highest weights $h_j$,

$$Y = \{h_j = N - j + m_j, \quad j = 1, 2, \dots, N\}, \quad (3.2)$$

where $m_j$ is the number of boxes in $j$'th row, with equal numbers of even and odd shifted highest weights $\{h_1, \dots, h_{N/2}, h_{N/2+1}, \dots, h_N\} = \{h_1^e, \dots, h_{N/2}^e, h_1^o, \dots, h_{N/2}^o\}$.<fnref n="8" /> The formula (3.1) was obtained by pure combinatorics of planar graphs in [26] and rederived in [27] from the DWG matrix model. The factors $\chi_{\{h\}}[t]$, $\chi_{\{h\}}[t^*]$ here are usual Schur characters — polynomials of couplings $t_k, t_k^*$. One way to define them is through Schur polynomials $P_n(\theta)$ which are read off from

$$e^{\sum_{n=1}^{\infty} z^n \theta_n} = \sum_{n=0}^{\infty} z^n P_n(\theta), \quad (3.3)$$

where

$$\theta_k = \frac{1}{k} \text{tr} B^k = \frac{N}{k} t_k, \quad (3.4)$$

and for $n < 0$ we define $P_n(\theta) = 0$. Then we have

$$\chi_{\{h\}}(t) = \det_{k,l} (P_{h_k+1-l}(\theta)), \quad (3.5)$$

and similarly for $\chi_{\{h\}}(t^*)$ where we use $A$ instead of $B$ in (3.4).

Notice that the representation (3.1) of the DWG matrix model partition function renders the symmetry between two sets of couplings $t_q \leftrightarrow t_q^*$ obvious.

A clear advantage of the representation (3.1) is the drastic reduction of the number of variables: instead of $N^2$ original matrix variables $M_{ij}$ it contains only the sums over $N$ highest weights $h_j$. Hence, to sum up DWG's we can apply the saddle point approximation to find the dominating large Young tableau for the sums in (3.1). But for that we have to learn how to compute the characters in this limit. The corresponding methods, as well as the formulas of the following subsection, have been worked out in [27–30].

<footnote n="7">⁷Strictly speaking the r.h.s. of (3.1) includes an extra factor $N^{N(N-1)/4-1/2} \sum_k h_k$ (see [27]) but it can be reabsorbed into a redefinition of the matrices $A, B$ by a scalar factor.</footnote>

<footnote n="8">⁸The highest weights are not necessarily ordered here according to the values, but they can always be reordered so due to the antisymmetry of the characters.</footnote>
## 3.2 Characters in the large N limit

Using the methods of the papers mentioned above, we can compute the Schur characters depending on arbitrary number of variables (couplings) at large N (for large smooth Young tableaux with characteristic shifted highest weights $h_j \sim N$). The computation is based on simple identities for Schur characters. The first one is

$$t_L \cdot \chi_{\{h\}}[t] = \frac{1}{N} \sum_{j=1}^{N} \chi_{\{h+L\delta_j\}}[t], \qquad (3.6)$$

where in the r.h.s. one adds L to one of the highest weights $h_j$ in the Young tableau and sums up the resulting character over all N insertions<fnref n="9" /> so that

$$\{h + L\delta_j\} = \{h_i + L\delta_{i,j}, \quad i = 1, 2, \dots, N\}. \qquad (3.7)$$

The second one is

$$\frac{L}{N} \partial_{t_L} \chi_{\{h\}}[t] = \sum_{j=1}^{N} \chi_{\{h-L\delta_j\}}[t], \qquad (3.8)$$

where in the r.h.s. one subtracts L from one of the highest weights in the Young tableau. While in the r.h.s. here we may potentially get negative values of the weights $h_k - L$ defining the character, this identity still holds<fnref n="10" /> with the definition (3.5) of the characters where we take $P_n = 0$ for $n < 0$.

Using (3.6) we also have, denoting the Vandermonde determinant by $\Delta$ (see the eqs. (2.2)-(2.3) of [27]),

$$\begin{aligned}
t_L &= \frac{1}{N} \sum_{k=1}^{N} \frac{\chi_{\{h+L\delta_k\}}[t]}{\chi_{\{h\}}[t]} \\
&= \frac{1}{N} \sum_{k=1}^{N} \frac{\Delta(h+L\delta_k)}{\Delta(h)} e^{\log \frac{\chi_{\{h+L\delta_k\}}}{\Delta(h+L\delta_k)} - \log \frac{\chi_{\{h\}}}{\Delta(h)}} \\
&= \frac{1}{N} \sum_{k=1}^{N} \prod_{j \neq k} \left(1 + \frac{L}{h_k - h_j}\right) e^{L\mathbf{F}(h_k)} + \mathcal{O}(1/N) \\
&= \frac{1}{L} \oint \frac{dh}{2\pi i} e^{L(H(h)+\mathbf{F}(h))} + \mathcal{O}(1/N) \qquad (3.9)
\end{aligned}$$

where we exponentiated the product in the large N limit (see [27] for details), when the characteristic $h_j \sim N$, introducing two basic functions: the resolvent of the highest weights

$$H(h) = \sum_{j=1}^{N} \frac{1}{Nh - h_j} \equiv \int \frac{dh' \rho(h')}{h - h'}. \qquad (3.10)$$

<footnote n="9">⁹The highest weights can then be reordered inside the character, using the permutational symmetry, to satisfy the inequalities $h_j > h_{j+1}$ following from the definition of $h_j$ in (3.2).</footnote>

<footnote n="10">¹⁰For a different definition of the character the identity may not hold when L is large enough. However, in what follows we will focus on the large N limit and the saddle point distribution of $h_k$ for which $h_k \sim N$, so for fixed L the combination $h_k - L$ is almost always positive and does not cause any problems.</footnote>
and

$$ \mathbf{F}(h) = \partial_{h_k} \log \frac{\chi_{\{h\}}[t]}{\Delta(h)} \Big|_{h_k=Nh} . \qquad (3.11) $$

We also introduced here the normalized highest weight variable $h \sim \mathcal{O}(N^0)$. One can also define a function similar to $\mathbf{F}$ but for the dual weights $\{t_n^*\}$,

$$ \mathbf{F}^*(h) = \partial_{h_k} \log \frac{\chi_{\{h\}}[t^*]}{\Delta(h)} \Big|_{h_k=Nh} \qquad (3.12) $$

and we get a relation similar to (3.9),

$$ t_L^* = \frac{1}{L} \oint \frac{dh}{2\pi i} e^{L(H(h) + \mathbf{F}^*(h))} . \qquad (3.13) $$

Let us point out that our notation here differs somewhat from [27], as there our function $\mathbf{F}^*$ was denoted by $F$. Furthermore the notation for $F$ differs between [27] and subsequent papers [28, 29]. We chose to use bold font here to avoid confusion. We will use $\mathbf{F}$, $\mathbf{F}^*$ here in section 3 while we discuss the general solution, and then we will use $F$ in later sections when we specialize to the $R^2$ DWG model discussed in [29] (see section 4.1 for a summary of the notation in the latter case).

Using (3.8) and the first line of (A.1) given in appendix A we also derive a similar formula [27, 28] for $\langle \text{tr}(B\tilde{M})^L \rangle$ evaluated in the dual model,

$$ \frac{1}{N} \langle \text{tr} (B\tilde{M})^L \rangle = \frac{1}{L} \oint_{C_H} \frac{dh}{2\pi i} e^{-L[H(h) + \mathbf{F}^*(h)]} \qquad (3.14) $$

where for large $h_j \sim N$ we also used the fact that the sum over representations reduces at the saddle point to a single, dominating Young tableau defined by the resolvent $H(h)$. For the correlator $\langle \text{tr}(MA)^L \rangle$ we find in a similar way the representation

$$ \frac{1}{N} \langle \text{tr} (AM)^L \rangle = \frac{1}{L} \oint_{C_H} \frac{dh}{2\pi i} e^{-L[H(h) + \mathbf{F}(h)]} . \qquad (3.15) $$

There is also an alternative formula for the same quantity, given in [27],

$$ \langle \frac{1}{N} \text{tr} (AM)^{2L} \rangle = \frac{1}{L} \oint \frac{dh}{2\pi i} h^L e^{L(H(h) + 2\mathbf{F}^*(h))} \qquad (3.16) $$

which however we will not use in this paper.

In order to compute $\langle \text{tr} M^L \rangle$ one can also use similar arguments (see appendix A for more details). As a result we get [27]

$$ \langle \frac{1}{N} \text{tr} M^{2L} \rangle = \frac{1}{L} \oint \frac{dh}{2\pi i} h^L e^{L H(h)} . \qquad (3.17) $$

### 3.3 Saddle point equation and general RH problem for DWG

From now on, we will consider the DWG models only with polynomial original and dual potentials, with the highest powers equal to $Q$ and $Q^*$, respectively, so that the potentials read

$$ V(M) = \sum_{q=1}^{Q} \frac{t_q}{q} M^q, \quad V_*(M) = \sum_{q=1}^{Q^*} \frac{t_q^*}{q} M^q . \qquad (3.18) $$
The large $N$ saddle point approximation for the multiple sum (3.1) then reads

$$ \partial_h \left[ \frac{1}{2} h_k \log h_k - \frac{1}{2} \log \Delta(h) + \log \chi_{\{h\}}[t] + \log \chi_{\{h\}}[t^*] \right] = 0. \quad (3.19) $$

We used here natural assumptions [27, 33, 34] that the densities of distributions of even and odd highest weights {$h_j^e$} and {$h_j^o$} are the same: $\rho^e(h) = \rho^o(h) = \frac{1}{2}\rho(h)$. We also applied the Stirling formula $\log h! \approx h \log h$ valid for typical $h_j \sim N$ at the saddle point. Using the definitions (3.11) and (3.10) we can rewrite (3.19) in the continuous form. For the DWG model with general sets of original and dual couplings we get:

$$ 2\mathbf{F} + 2\mathbf{F}^* + 3\mathcal{H} + \ln h = 0 \quad (3.20) $$

where by slash we denote the symmetric part of the function $H$ on its cut.

To be precise, here we need to discuss (following [27]) an important property of the density $\rho(h)$. The density is naturally defined on the interval $h \in [0, a]$ going from zero to some endpoint which we denote as $a > 0$. In fact $\rho(h)$ is saturated at its maximal<fnref n="11" /> value $\rho(h) = 1$ at the start of this interval, i.e. for $h \in [0, b]$ with $b < a$. Accordingly, $H(h)$ has a logarithmic cut for $h \in [0, b]$. As discussed in [27] the saddle point equation (3.20) holds only on the remaining part of the interval where $H$ has a nontrivial cut, i.e. for $h \in [b, a]$ (with $H(h \pm i0) = \mathcal{H}(h) \mp i\pi\rho(h)$).

The saddle point equation should be supplemented by equations for the character functions $\mathbf{F}$ and $\mathbf{F}^*$ defining the characters as functionals of this highest weight distribution. Such equations were derived in [27] from the observations described in the previous subsection. Namely, we introduce the functions

$$ \mathbf{G}(h) = e^{\mathbf{F}(h)+H(h)}, \quad \mathbf{G}^*(h) = e^{\mathbf{F}^*(h)+H(h)}. \quad (3.21) $$

Notice that like for $\mathbf{F}$, $\mathbf{F}^*$ our notation differs from the one of [28, 29] as there $G$ denoted a quantity closely related to our $\mathbf{G}^*$ (see section 4.1 for more details on the notation).

The functions $\mathbf{G}(h)$ and $\mathbf{G}^*(h)$ have $Q$ and $Q^*$ sheets, respectively. As the eqs. (3.9), (3.14) and (3.15) show, the behavior at zero for each one is given by expansions

$$ h(\mathbf{G}) = \sum_{q=1}^{Q} t_q \mathbf{G}^{-q} + \sum_{q=0}^{\infty} \mathbf{G}^q \langle \frac{1}{N} \mathrm{tr} (AM)^q \rangle = \mathbf{G}^{-1} V'(\mathbf{G}^{-1}) + \mathbf{G}^{-1} \mathcal{W}(\mathbf{G}^{-1}) \quad (3.22) $$

$$ h(\mathbf{G}^*) = \sum_{q=1}^{Q^*} t_q^* \mathbf{G}^{*-q} + \sum_{q=0}^{\infty} \mathbf{G}^{*q} \langle \frac{1}{N} \mathrm{tr} (B\tilde{M})^q \rangle = \mathbf{G}^{*-1} V'_*(\mathbf{G}^{*-1}) + \mathbf{G}^{*-1} \tilde{\mathcal{W}}(\mathbf{G}^{*-1}) \quad (3.23) $$

where $\tilde{M}$ is the matrix for the dual matrix model and we put $t_0 = t_0^* = 1$ while $V'(M)$ denotes the derivative of the potential. We assumed here that all the small $G$ singularities come from the original matrix potentials.

<footnote n="11">¹¹The constraint $\rho \le 1$ follows from the definition of the density and the fact that the shifted highest weights $h_i$ are strictly decreasing with $i$ as follows from (3.2).</footnote>
These expansions show that the functions $G(h)$ and $G^*(h)$ have branch points of order $Q$ and $Q^*$ at $h \to \infty$. This suggests another form of eqs. (5.34) and (3.23) (see eq. (3.23) of [29]):

$$ \frac{(-1)^{Q-1}}{t_Q} \prod_{q=1}^{Q} \mathbf{G}_q = \frac{1}{h} e^{H(h)} = \frac{(-1)^{Q^*-1}}{t_{Q^*}} \prod_{q=1}^{Q^*} \mathbf{G}_q^* \qquad (3.24) $$

where by $\mathbf{G}_q$ we denoted the values of the function $\mathbf{G}$ on $Q$ different sheets, and similarly for $\mathbf{G}^*$. If we skip the label we always mean the main sheet $\mathbf{G} = \mathbf{G}_1$.

The equations (3.20), (3.24), together with the small $\mathbf{G}$, $\mathbf{G}^*$ asymptotics (5.34), (3.23), give the complete Riemann-Hilbert problem for the DWG model with general sets of $Q$ original and $Q^*$ dual couplings.

# 4 The $\mathbb{R}^2$ quantum gravity solution

While above we discussed general DWG models, from now on we will specialize to the particular case of DWG model described in section 2.2 and studied in [29]. In this section we will review its exact solution from [29]. The model corresponds to quadrangulations with special, exponential weights for the vertices. We keep for it the name “$\mathbb{R}^2$ QG” model given in [29]. The main advantage of this model, apart from the fact that the solution for $H(h)$ can be given explicitly in terms of elliptic functions, is that the couplings are chosen in such a way that we have a parameter controlling the local “curvature” fluctuations (number of nearest neighbors for the vertices of quadrangulations) together with another one controlling the area (number of squares in the quadrangulation). Thus we can nontrivially restrict the geometry of the graphs and find interesting continuum limits.

In the paper [29] the sphere partition function of such quadrangulations has been studied. In the current paper we will generalize this analysis to the disc partition functions. Namely, we will compute in the near-critical regime the corresponding resolvents $\tilde{W}(z)$ and $W(z)$ defined above. In this section we summarise the exact solution which serves as the starting point for our new computation.

## 4.1 Notation and key relations for the case of even potentials

From now on we will follow the notation of [29] which is tailored for the case (of which the $\mathbb{R}^2$ QG model is an example) where the potentials for both the original and the dual models are even, so that $t_{2n+1} = t_{2n+1}^* = 0$. In this case we have tr $A^{2n+1} = 0$ and one can parameterize the matrix $A$ as

$$ A = \begin{pmatrix} \sqrt{a} & 0 \\ 0 & -\sqrt{a} \end{pmatrix} \qquad (4.1) $$

where $a$ is an $N/2 \times N/2$ matrix. Then instead of (3.12) one uses the function $F$ given by

$$ F(h_k) = 2 \frac{\partial}{\partial h_k^e} \log \frac{\chi_{\{\frac{h_e}{2}\}}(a)}{\Delta(h_e)} \qquad (4.2) $$
so it is defined in terms of only the even weights $h_k^e$ with $k = 1, \dots, N/2$,<fnref n="12" /> and the character of the matrix $a$. It can be deduced from [28, 29] that it is related to our notation as<fnref n="13" />

$$F(h) = 2\mathbf{F}^*(2h) + 2H(2h) - H(h). \qquad (4.3)$$

Moreover the function $G$ used in [28, 29] differs from our $\mathbf{G}, \mathbf{G}^*$ and is defined as

$$G(h) = e^{H(h)+F(h)}. \qquad (4.4)$$

Then analogs of equations (3.13), (3.14) read<fnref n="14" />

$$t_{2q}^* = \oint \frac{dG}{2\pi i G} h(G) G^q, \qquad (4.5)$$

and

$$\frac{1}{N} \langle \mathrm{tr} (B \tilde{M})^{2q} \rangle = \oint \frac{dG}{2\pi i G} h(G) G^{-q}, \qquad (4.6)$$

so we have

$$\begin{aligned} h(G) &= 1 + \sum_{q=1}^{Q} t_{2q} G^{-q} + \sum_{q=1}^{\infty} G^q \langle \frac{1}{N} \mathrm{tr} (B \tilde{M})^{2q} \rangle \qquad (4.7) \\ &= 1 + \sum_{q=1}^{Q} t_{2q} G^{-q} + \frac{1}{\sqrt{G}} \tilde{W} \left( \frac{1}{\sqrt{G}} \right). \end{aligned}$$

Note that for the $R^2$ QG model we have $B = A_4$ where $A_4$ is the matrix discussed in section 2.2 with the property $\frac{1}{N} \mathrm{tr} A_4^n = \delta_{n,4}$.

Let us also note that in the case of even potential we find from (3.17) that [27]

$$zW(z) \equiv \left\langle \frac{1}{N} \mathrm{tr} \frac{z}{z-M} \right\rangle = 1 - \oint_H \frac{dh}{2\pi i} \log \left( z^2 - h e^{H(h)} \right). \qquad (4.8)$$

By the Lagrange inversion formula (we have a single cut on the main sheet of $H(h)$) we get from here the parameterisation of the resolvent as<fnref n="15" />

$$zW(z) = z^2 - h, \qquad z^2 = he^{H(h)} \qquad (4.9)$$

$$zW(z) \equiv \left\langle \frac{1}{N} \mathrm{tr} \frac{z}{z - AM} \right\rangle = 1 - \oint_{CH} \frac{dh}{2\pi i} \log \left( z^2 - h e^{2\mathbf{F}^*(h) + H(h)} \right),$$

the same inversion trick unfortunately does not work since the integration contour $C_H$ here is pinched between the cuts of $H(h)$ and $\mathbf{F}^*(h)$.

<footnote n="12">¹²We assume that, by symmetry, both even and odd weights are distributed in the same way.</footnote>

<footnote n="13">¹³The difference between $\mathbf{F}^*$ and $F$ comes mainly from the fact that in (4.2) we have the Vandermonde of only the even weights $h_k^e$, while in (3.12) we have the Vandermonde for all the weights.</footnote>

<footnote n="14">¹⁴Strictly speaking one should rederive these equations for the even case as the contour of integration used in the trick in (3.9) becomes subtle to choose due to overlapping cuts (see [28] for the derivation and comments on this).</footnote>

<footnote n="15">¹⁵For $W(z)$, which in the even case is given by (see [27])</footnote>
## 4.2 Explicit solution for the density

To derive the saddle point equation for this case, the easiest way is to plug into (3.1) the character $\chi_{\{h\}}[t]$ in the form (A.4). This gives for the partition function

$$Z(t, t^*) = \sum_{\{h^e, h^o\}} \frac{\prod_i(h_i^e - 1)!! h_i^o!!}{\prod_{i,j}(h_i^e - h_j^o)} \prod_{\epsilon=0}^{3} \frac{\Delta(h^{(\epsilon)})}{\prod_i \left(\frac{h_i^{(\epsilon)} - \epsilon}{m}\right)!} \\ \times \operatorname{sgn}\left[\prod_{0 \le \epsilon_1 < \epsilon_2 \le 3} \prod_{i,j} (h_i^{(\epsilon_2)} - h_j^{(\epsilon_1)})\right] \chi_{\{h\}}[t^*]. \quad (4.10)$$

Here following [29] we introduced four groups with equal numbers of weights $h_i^{(\epsilon)}$, $\epsilon = 0, 1, 2, 3$, $i = 1, \dots, N/4$ such that $h_i^{(\epsilon)} = \epsilon \mod 4$. We will assume that all four groups are distributed with the same density and the sgn[...] factor is irrelevant in the large N limit. That provides the saddle point equation given in eq. (4.7) of [29]. It reads

$$2F + \mathcal{H} + \log h = 0, \quad h \in [b, a]. \quad (4.11)$$

where by $\mathcal{H}(h)$ we defined the principal value of the resolvent $H(h)$ on the cut $C_H = [b, a]$, $H = \mathcal{H} \mp i\pi\rho(h)$.

To complete the RH problem defining two functions — the resolvent $H(h)$ and the character function $F(h)$ — first we notice that with our definition of the couplings $t_n^*$ of the dual potential (2.15), we obtain from (3.23)

$$h - \frac{\lambda\beta^2}{\epsilon} \frac{1}{G} - \frac{\lambda G}{(G - \epsilon)} = \frac{1}{\sqrt{G}} \tilde{W}\left(\frac{1}{\sqrt{G}}\right) \quad (4.12)$$

where we used the definition of the resolvent: $\tilde{W}(z) \equiv \sum_{q=0}^{\infty} \langle \text{tr} (\tilde{M} A_4)^{2q} \rangle z^{-2q-1}$. We notice that for $G \to 0$ the r.h.s. of the last formula is vanishing since $\tilde{W}(z) \simeq 1/z$ for $z \to \infty$. The vanishing l.h.s. then suggests that the function $G$ has only two sheets connected through the cut of $F(h)$, so that we can put $Q = 2$ in (3.24). Then the logarithm of equation (3.24) looks as follows:

$$\log G(h + i0) + \log G(h - i0) - H(h) + \log \frac{h}{(\beta^2 - 1)\lambda} = 0, \quad h \in [d, c]. \quad (4.13)$$

which gives, on the defining sheet of $F(h)$, the second needed RH equation

$$2\mathcal{H}(h) + H(h) + \log \frac{h}{(\beta^2 - 1)\lambda} = 0 \quad h \in [d, c]. \quad (4.14)$$

We used the fact that $F(h)$ has a cut $C_F : [d, c]$ on its defining sheet, whereas $H(h)$ does not have a $[d, c]$ cut.

Remarkably, both equations (4.11) and (4.14) of this RH problem turn out to be written for the same function $2F(h) + H(h) + \log h$ which thus has only two sheets. This helps to immediately write its solution in terms of the Cauchy integrals (see eq. (4.8) in [29]) which can be expressed in terms of incomplete elliptic integrals. Since by virtue of (4.11) we also have

$$2F(h) + H(h) + \log h = i\pi\rho \quad (4.15)$$
it immediately gives us the explicit formula for the density of the highest weights of the saddle point Young tableau.

In order to write the solution explicitly, let us introduce some notation. We define the elliptic moduli $k^2$ and $k'^2$ corresponding to the branch points as

$$k^2 = \frac{(a-b)(c-d)}{(a-c)(b-d)}, \quad k'^2 = 1-k^2, \qquad (4.16)$$

and we denote by $K$ and $K'$ the complete elliptic integrals of the first kind with moduli $k$ and $k'$ respectively. The corresponding elliptic nome $q$ and its dual $q'$ read

$$q = e^{-\pi K'/K}, \quad q' = e^{-\pi K'/K'}. \qquad (4.17)$$

We also describe the Mathematica conventions for some of these functions (and others that we use below) in appendix B. In this notation the result for the density found in [29] reads

$$\rho(h) = \frac{u}{K} - \frac{i}{\pi} \ln \left[ \frac{\theta_4(\frac{\pi}{2K}(u-iv), q)}{\theta_4(\frac{\pi}{2K}(u+iv), q)} \right] \qquad (4.18)$$

where $u$ is defined in terms of the inverse Jacobi elliptic function sn,

$$u = \text{sn}^{-1} \left( \sqrt{\frac{(a-h)(b-d)}{(a-b)(h-d)}}, k \right). \qquad (4.19)$$

In the next subsection we describe equations that fix the branch points $a, b, c, d$ in terms of the couplings of the model.

## 4.3 Fixing parameters of the solution

The branch cut endpoints $a, b, c, d$ are fixed in terms of the couplings $\lambda, \beta$ and $\epsilon$ by a system of rather nontrivial equations found in [29]. It is obtained by imposing constraints such as correct normalization of the resolvent at large $h$. We refer the reader to [29] for details of the derivation, and here we present the final result, which we also rewrote in a slightly more compact form.

Let us first introduce some additional notation, namely

$$v = \text{sn}^{-1} \left( \sqrt{\frac{a-c}{a-d}}, k' \right), \quad s = \text{sn}(v, k'). \qquad (4.20)$$

In addition, we define

$$\Delta \equiv b-c. \qquad (4.21)$$

Then the equations that fix $a, b, c, d$ read [29]

$$v = -K' - \frac{K}{\pi} \log(\lambda(1 - \beta^2)) \qquad (4.22)$$

$$\lambda = \frac{2Kq}{i\pi^2} \frac{\theta_1(\frac{i\pi v}{K}, q)}{\theta'_1(0, q)} \left( -E + K \left( k'^2 s^2 + \frac{\Upsilon\Xi}{s} + 2\Xi^2 \right) \right) \qquad (4.23)$$

$$\Delta = \frac{4Kq}{\epsilon^2 i\pi} \frac{\theta_1(\frac{i\pi v}{K}, q)}{\theta'_1(0, q)} \frac{k'^2 s\sqrt{1-s^2}}{\sqrt{1-k'^2 s^2}}, \qquad (4.24)$$

$$b-1 = \frac{1}{\epsilon^2} \left[ \frac{2Kq}{i\pi} \frac{\theta_1(\frac{i\pi v}{K}, q)}{\theta'_1(0, q)} \left( \frac{k'^2 s^4 - 2s^2 + 1}{s\sqrt{1-s^2}\sqrt{1-k'^2 s^2}} + 2\Xi \right) - 2\lambda \right] \qquad (4.25)$$
where

$$ \Upsilon = \frac{3k'^2 s^4 - 2(k'^2 + 1)s^2 + 1}{\sqrt{1 - s^2}\sqrt{1 - k'^2 s^2}}, \quad \Xi = \frac{\pi}{2K} + E(v, k') + \left(\frac{E}{K} - 1\right) v \quad (4.26) $$

and $E$ is the complete elliptic integral of the second kind (similarly to $K$) while $E(v, k')$ denotes the incomplete elliptic integral of the 2nd kind (see appendix B for the corresponding Mathematica notation).

Let us elaborate briefly on the structure of these four equations (4.22)–(4.25). Our goal is to find $a, b, c, d$ as functions of the three couplings $\lambda, \beta, \epsilon$. We see that plugging $v$ from the first equation into the second one gives an equation that completely fixes $k^2$ in terms of two couplings $\lambda$ and $\beta$. The first equation then fixes $v$, and from the last two equations we find $b$ and $c$. Finally, having found $k$ and $v$ and recalling their definitions (4.16), (4.20) we fix the remaining two parameters $a, d$.

Using identities between elliptic functions, one can rewrite these equations in a variety of ways. Here we presented the form we found the most useful for our calculation.

## 4.4 Simplification and computation of $G$

Below we will be interested in computing the resolvent $\tilde{W}$ which is encoded in $G(h)$ and thus it is useful to write the latter function explicitly. Since $2F + H + \log h = \pm i\pi\rho$ we first get (choosing the sign accordingly)

$$ 2 \log G(h) + \log h = (2F + H + \log h) + H = -i\pi\rho + \int_0^a \frac{dx}{h-x} \rho(x). $$

We found it useful to also write this result in a slightly different way. Notice that we can get rid of the term without theta functions in the density (4.18) by switching from $\theta_4$ to $\theta_1$ using the identity

$$ \theta_4\left(\frac{\pi}{2K}u + \frac{i\pi}{2K}K', q\right) = \frac{i}{q^{1/4}} \exp\left(-\frac{i\pi u}{2K}\right) \theta_1\left(\frac{\pi}{2K}u, q\right). \quad (4.27) $$

This gives

$$ \rho(h) = -\frac{i}{\pi} \ln \left[ - \frac{\theta_1\left(\frac{\pi}{2K}(u - i(v + K')), q\right)}{\theta_1\left(\frac{\pi}{2K}(u + i(v + K')), q\right)} \right]. \quad (4.28) $$

In addition, the form of the density suggests one to rewrite the integral defining $H(h)$ in terms of $u$ (defined in (4.19)) rather than $h$. Let us denote by $w(h')$ the same function as in (4.19) but with $h$ replaced by $h'$. Then in the definition of $H(h)$ given by (3.10) we change the integration variable $h' \to w(h')$ and obtain

$$ H(h) = -\frac{2i k^2 (1 - k'^2 s^2)}{\pi k'^2} \int_0^K \frac{dw \, \mathrm{cn}(w) \mathrm{dn}(w) \mathrm{sn}(w) \log \left[ -\frac{\theta_1\left(\frac{\pi}{2K}(w-i(v+K')), q\right)}{\theta_1\left(\frac{\pi}{2K}(w+i(v+K')), q\right)} \right]}{(s^2 \mathrm{dn}(w)^2 - 1)^2 \left[ \frac{h-b}{\Delta} + \frac{1-\frac{1}{k'^2} \mathrm{dn}(w)^2}{1-s^2 \mathrm{dn}(w)^2} \right]} \quad (4.29) $$

We can also write the argument $h$ as a function of $u$ by inverting (4.19),

$$ \frac{h-b}{\Delta} = \frac{1}{(k^2-1)s^2} - \frac{(k^2-1)s^2+1}{(k^2-1)s^2(k^2s^2 \mathrm{sn}(u)^2-s^2+1)}. \quad (4.30) $$
Then finally for $G(h)$ we have

$$2 \log G(h) = -\frac{2i k^2 (1 - k'^2 s^2)}{\pi k'^2} \int_0^K \frac{dw \operatorname{cn}(w) \operatorname{dn}(w) \operatorname{sn}(w) \log \left[ -\frac{\theta_1\left(\frac{\pi}{2K}(w-i(v+K')), q\right)}{\theta_1\left(\frac{\pi}{2K}(w+i(v+K')), q\right)} \right]}{(s^2 \operatorname{dn}(w)^2 - 1)^2 \left[ \frac{h-b}{\Delta} + \frac{1-\frac{1}{k'^2} \operatorname{dn}(w)^2}{1-s^2 \operatorname{dn}(w)^2} \right]} -i\pi\rho(h) - \log(h-b). \quad (4.31)$$

This is the form of $G(h)$ that we will use below.

## 4.5 Scaling properties of $G$ as a function of $\epsilon$

Let us show that dependence of $G(h)$ on $\epsilon$ is quite trivial. From (4.7) we see that to obtain the resolvent $\tilde{W}(z) = \langle \operatorname{tr} \frac{1}{z - \tilde{M}A_4} \rangle$ we need to compute the function $G(h)$ and invert it. Denoting

$$z \equiv \sqrt{\epsilon/G} \quad (4.32)$$

we have

$$h = z^2 \frac{\lambda \beta^2}{\epsilon^2} + \frac{\lambda z^4}{\epsilon^2(1-z^2)} + \frac{z}{\sqrt{\epsilon}} \tilde{W}[(\sqrt{\epsilon}/z)^{-1}]. \quad (4.33)$$

On the other hand, as we saw above from the combinatorial argument leading to (2.23), the resolvent should have the form

$$\tilde{W}(z) = \frac{1}{z} \left(1 + \epsilon^{-2} \phi_1(\sqrt{\epsilon} z)\right) \quad (4.34)$$

where the function $\phi_1$ does not depend on $\epsilon$ and $z$. Plugging this into (4.33) we get

$$\epsilon^2(h-1) = \lambda \beta^2 z^2 + \frac{\lambda z^4}{1-z^2} + \phi_1(z). \quad (4.35)$$

Inverting this relation we find

$$G = \epsilon \phi_2(\epsilon^2(h-1)) \quad (4.36)$$

where $\phi_2$ is again some function which does not depend on $\epsilon$ or $h$ explicitly. As a result, this predicts for $G$ the scaling property

$$G(h, Z\epsilon) = ZG(Z^2(h-1) + 1, \epsilon). \quad (4.37)$$

Let us verify that the result for $G$ given above in (4.31) indeed satisfies this relation. We see that in the equations (4.22)–(4.25) which define parameters of the solution the variable $\epsilon$ appears only in $\Delta$ and in one term in the last equation (4.25). Using this, one can easily check that all terms in (4.31) except the last one are functions of the combination $\epsilon^2(h-1)$ and thus sending $\epsilon \to Z\epsilon$ for them is the same as replacing $h \to Z^2(h-1) + 1$. The very last term $\log(h-b)$ in (4.31), on the other hand, produces an extra $\log Z$ term under this transformation. As a result we see that (4.37) is perfectly satisfied!

Due to this scaling symmetry we see that $\epsilon$ cannot affect any critical behavior properties of the resolvent $\tilde{W}$. Accordingly, for the purpose of studying them we will set $\epsilon = 1$ in the subsequent computations in sections 5.1–5.3.2. Then in section 5.4 we will study a different kind of resolvent, namely $W(z)$, for which the dependence on $\epsilon$ is more tricky. There we will restore the $\epsilon$-dependence in the needed equations, which is not hard to do as we will see.
## 5 Critical regimes of large area and boundary of the quadrangulated disc

In this section, we will discuss the critical regimes for the partition function of the quadrangulated disc. Namely, we will first tune the parameter $\lambda$ in the resolvents $\tilde{W}(z)$ and $W(z)$ for the elliptic solution (4.18) of sections 4.3, 4.4 to its critical value $\lambda_c$ in such a way that the area (number of plaquettes in the graph) becomes very large. This allows one to consider the typical quadrangulations which “forget” about the details of their discretization and can be considered as relatively smooth 2d manifolds with fluctuating metric. Generically, this regime leads to the well known solution of pure 2d quantum gravity [5–7], as it was demonstrated in [29]. However, to study the transition between this pure 2d QG regime and the “almost flat” regime, when the curvature is maximally suppressed, we have an extra parameter to adjust, namely $\beta$ which controls the fluctuations of the modulus of overall curvature accumulated on the disc. Tuning $\beta \to 0$ one can “flatten” the disc in such a way that the vertices with four neighbors will dominate, with rare insertions of the conical curvature defects, having angle deficits $+\pi$ (positive curvature defect) or $-\pi, -2\pi, -3\pi, \dots$ (negative curvature defects). By appropriate simultaneous tuning of $\lambda \to \lambda_c$ and $\beta \to 0$, with a certain double scaling, this limit interpolates between the regime of pure 2d gravity, when the size of the manifold is large enough to accommodate the “pure gravity” fluctuations of the metric, and “almost flat” (AF) regime when the size of the manifold is relatively small and thus it is almost completely flattened, apart from rare positive conical curvature defects (all negative curvature is then concentrated at the boundary).<fnref n="16" />

The universal partition function interpolating between two such regimes was established in [29]. We will extend in this section the analysis of this “$R^2$ gravity” solution of the DWG model to the disc partition function. Namely, we will explicitly calculated the resolvents $\tilde{W}(z)$ and $W(z)$ in this “near-flat” regime — large area of quadrangulations but arbitrary length of the boundary — which is the main result of this paper. Then we will study these quantities in the two limits mentioned above. In the limit of pure gravity, which demands another double scaling computation, adjusting the boundary cosmological constant $z \to z_c$ simultaneously with $\lambda \to \lambda_c$, we will reproduce for both resolvents $\tilde{W}(z)$ and $W(z)$ the well known universal disc amplitude [10]. In the limit of almost flat 2d manifolds with minimal number of curvature defects we will reproduce the disc partition function from [28] (see eq. (4.17) there).

As a preliminary step, we will first investigate a special regime for the general solution, when the elliptic nome $q$ goes to 1 but we drop only the exponentially small terms. It is convenient to parametrize $q$ as

$$q = e^{-\pi\tau} \qquad (5.1)$$

so that $\tau \to 0$ and we work with exponential accuracy, i.e. up to $O(e^{-\pi/\tau})$ terms. Let us recall that all parameters of the solution (4.18) for the density of highest weights $\rho(h)$ are fixed in terms of two original parameters of the model $\lambda$ and $\beta$. In practice it will be convenient for us to first deal with intermediary parameters $v$ (defined in (4.20)) and

<footnote n="16">¹⁶The quantisation unit of the angle deficit $\delta\phi = \pi$ should be an important relevant parameter in such double scaling limit, and the result would be different for a different value of $\delta\phi$. It would be interesting to solve a model with different $\delta\phi$, say for triangulations.</footnote>
$\tau$ instead. We will compute $G(h)$ at generic $v$ and with exponential precision in $\tau$. We present this calculation in subsection 5.1.

Next, in subsection 5.2 we recall (following [29]) the main features of the critical regime in our model, and in particular the case when $v$ is small (in addition to $q \to 1$). We will refer to this regime as the ‘flattening’ limit. In subsection 5.3 we will describe the properties of our solution in this regime in detail, making contact with the results in [29] and presenting a number of nontrivial checks of the result. We also discuss how the solution interpolates between the pure gravity and the AF regimes. Finally in subsection 5.4 we compute $W(z)$, i.e. the resolvent for $\langle \mathrm{Tr} M^n \rangle$ expectation values, study its properties and reproduce the same pure gravity limiting behavior for it as well.

## 5.1 Solution with exponential precision for near-flat regime

We will focus here on the limit $\tau \to 0$ which was explored for the partition function of the model in [29]. Here we will compute the function $G(h)$ in this limit with exponential precision in $\tau$. We will keep the second remaining parameter $v$ finite so that we do not yet specialize to the vicinity of the critical line discussed in [29] (we will do that in the next subsection).

### 5.1.1 Expanding the parameters

Let us first compute how the cut endpoints $a, b, c, d$ and related parameters of the solution behave in the limit $\tau \to 0$. With exponential precision in $\tau$ we have

$$K \approx \frac{\pi}{2\tau}, \quad E(k) \approx 1 \qquad (5.2)$$

while

$$\mathrm{sn}(v, k') \approx \sin v, \quad E(v, k') \approx v. \qquad (5.3)$$

Plugging this into (4.26) we find

$$\Upsilon \approx \frac{\cos 2v}{\cos v}, \quad \Xi \approx \tau + \frac{2\tau v}{\pi}. \qquad (5.4)$$

We also need to expand the theta functions appearing in (4.24). Using a modular transformation to change their modulus from $q$ to $q'$, we find with exponential precision

$$\frac{\theta_1\left(\frac{i\pi v}{K}, q\right)}{\theta'_1(0, q)} \approx i\tau e^{\frac{4\pi v^2}{\pi}} \sin(2v). \qquad (5.5)$$

Then (4.24) gives the result for $\Delta$,

$$\Delta \equiv b - c \approx 16e^{-\pi/\tau} e^{\tau(4v^2/\pi - \pi)} \sin^2 2v \qquad (5.6)$$

and we get from (4.23), (4.25)

$$\log(\lambda(1 - \beta^2)) = -\tau(2v + \pi) \qquad (5.7)$$

$$\lambda \approx \frac{e^{\frac{4\pi v^2}{\pi}} - \pi\tau}{\pi^2} \left( (\tau(2v + \pi)^2 - \pi) \sin(2v) + \pi(2v + \pi) \cos(2v) \right) \qquad (5.8)$$
$$b \approx 1 + \frac{e^{\frac{4\tau v^2}{\pi}} - \pi\tau}{\pi^2} \left(2(-4\tau v^2 - 2\pi\tau v + \pi)\sin(2v) - 4\pi v \cos(2v)\right) \quad (5.9)$$

$$a \approx 1 + \frac{e^{\frac{4\tau v^2}{\pi}} - \pi\tau}{\pi^2} \left(4\sin(v)\left(\left(\pi - 2\tau v(2v + \pi)\right)\cos(v) + \pi^2\sin(v)\right) - 4\pi v \cos(2v)\right) \quad (5.10)$$

$$d \approx 1 + \frac{e^{\frac{4\tau v^2}{\pi}} - \pi\tau}{\pi^2} \left(2(\pi - 2\tau v(2v + \pi))\sin(2v) - 4\pi(\pi \cos^2(v) + v \cos(2v))\right) \quad (5.11)$$

The first two of these relations are equivalent to those given in equation (5.5) of [29].<fnref n="17" /> We also find that the relation between $h$ and $u$ from (4.19) in our regime becomes

$$h - b = \frac{e^{\tau(4v^2/\pi - \pi)} \sin^2 2v}{\cosh^2 u - \sin^2 v} \quad (5.12)$$

### 5.1.2 Computing the integral for $G(h)$

The nontrivial part of $G(h)$ we need to compute is the integral in (4.31). The term with theta functions in the integrand is the density (4.28) which we can write with exponential precision as

$$\rho = -\frac{i}{\pi} \log \left[ -\frac{\theta_1\left(\frac{\pi}{2K}(u-i(v+K')), q\right)}{\theta_1\left(\frac{\pi}{2K}(u+i(v+K')), q\right)} \right] \approx \frac{4\tau}{\pi^2} u(v+\pi/2) - \frac{i}{\pi} \log \left( \frac{\cosh(u-iv)}{\cosh(u+iv)} \right) \quad (5.13)$$

To write the rest of the integrand in our limit we can expand it directly, or, more conveniently, first recall that the integral we need is simply $\int_b^a \frac{dx}{h-x} \rho(x)$ which we can rewrite with $w$ as the integration variable,

$$\int_b^a \frac{dx}{h-x} \rho(x) = - \int_0^K dw \frac{dh(w)}{dw} \frac{1}{h(u) - h(w)} \rho(w) \quad (5.14)$$

Then we can compute $dh(w)/dw$ with exponential accuracy in our limit $\tau \to 0$ from (5.12), which gives

$$\frac{dh(w)}{dw} \frac{1}{h(u) - h(w)} = - \frac{\sinh w}{\cosh w + \cosh u} - \frac{\sinh w}{\cosh w - \cosh u} + \frac{\sinh w}{\cosh w + \sin v} + \frac{\sinh w}{\cosh w - \sin v} \quad (5.15)$$

Notice also that in our integral (5.14) we can (with exponential precision) replace the upper integration limit $K$ by $\infty$ and use the approximation for the density (5.13) on the whole region of integration, as we justify in detail in appendix C.

From the expression for the density (5.13) we see that there are two types of integrals we need to compute, corresponding to the two terms in that equation. The first type of integrals reads<fnref n="18" />

$$\int_0^\infty dw \frac{wh'(w)}{h(u) - h(w)} = u^2 + i\pi u + v^2 - \frac{\pi^2}{4} \quad (5.16)$$

<footnote n="17">¹⁷Notice the second relation in (5.5) from [29] contains a typo, it should read $q = (\lambda(1 - \beta^2))^{1/(1+2v/\pi)}$.</footnote>

<footnote n="18">¹⁸The proper choice of branch we use here corresponds to $\operatorname{Im} u < 0$ and real $v > 0$.</footnote>
which is relatively straightforward to compute in Mathematica. The second type of integrals is

$$I = - \int_{0}^{\infty} dw \frac{h'(w)}{h(u) - h(w)} \log \frac{\cosh(w + iv)}{\cosh(w - iv)}. \qquad (5.17)$$

It is much more nontrivial but can still be computed with the help of a number of tricks. The initial result we found in Mathematica contains many Li₂ polylogarithms, which can then be further reduced using specialised software (see e.g. [35, 36]).<fnref n="19" /> The final outcome reads

$$I = -\mathrm{Li}_2(-e^{2iv-2u}) - \mathrm{Li}_2(-e^{2u+2iv}) - 2\mathrm{Li}_2(1-e^{2iv}) - 2\mathrm{Li}_2(1+e^{2iv}) \quad (5.18) \\ -2u^2 + 4iuv + 2v^2 + \frac{\pi^2}{3} + (-2u-2iv-i\pi)\log(1+e^{2u+2iv}) \\ +(2u-2iv+i\pi)\log(e^{2u}+e^{2iv}) + 2i\pi\log(1+e^{2iv}).$$

Combining all the parts together and using (5.7), (5.8), we get the final result for G,

$$G = (\lambda(1 - \beta^2))^{-\frac{1}{2} + \frac{v}{\pi}} \frac{\sqrt{\cosh^2 u - \sin^2 v}}{\sin 2v} \exp(\Phi) \quad (5.19)$$

where

$$\Phi = u^2 \left( -\frac{\tau(2v + \pi)}{\pi^2} - \frac{i}{\pi} \right) + u \left( -\frac{4v}{\pi} - \frac{2i\tau(2v + \pi)}{\pi} \right) + \frac{\tau(\pi - 2v)(2v + \pi)^2}{4\pi^2} \\ + \frac{1}{6}i \left( \frac{18v^2}{\pi} - 12v + \pi \right) - \log \left( 1 + e^{2iv} \right) \quad (5.20) \\ + \frac{(iu + v - \pi) \log (1 + e^{2u-2iv})}{\pi} + \frac{(-iu + v + \pi) \log (1 + e^{2u+2iv})}{\pi} \\ - \frac{i}{2\pi} \left[ \mathrm{Li}_2(-e^{2iv-2u}) + \mathrm{Li}_2(-e^{2u+2iv}) + 2\mathrm{Li}_2(1-e^{2iv}) + 2\mathrm{Li}_2(1+e^{2iv}) \right]$$

In view of (5.7), (5.8) we can exclude $\tau$ so this last equation can be also written as

$$\Phi = u^2 \left( \frac{\log(\lambda(1 - \beta^2))}{\pi^2} - \frac{i}{\pi} \right) + u \left( -\frac{4v}{\pi} + \frac{2i \log(\lambda(1 - \beta^2))}{\pi} \right) \\ - \frac{\log(\lambda(1 - \beta^2)) (\pi^2 - 4v^2)}{4\pi^2} + \frac{1}{6} i \left( \frac{18v^2}{\pi} - 12v + \pi \right) - \log(1 + e^{2iv}) \quad (5.21) \\ + \frac{(iu + v - \pi) \log(1 + e^{2u-2iv})}{\pi} + \frac{(-iu + v + \pi) \log(1 + e^{2u+2iv})}{\pi} \\ - \frac{i}{2\pi} \left[ \mathrm{Li}_2(-e^{2iv-2u}) + \mathrm{Li}_2(-e^{2u+2iv}) + 2\mathrm{Li}_2(1-e^{2iv}) + 2\mathrm{Li}_2(1+e^{2iv}) \right]$$

We remind that $h$ and $u$ are related in our limit via (5.12), into which can also substitute (5.7) to get

$$h - b = (\lambda(1 - \beta^2))^{1 - \frac{2v}{\pi}} \frac{\sin^2 2v}{\cosh^2 u - \sin^2 v}. \qquad (5.22)$$

<footnote n="19">¹⁹We thank Ömer Gürdoğan for help with simplifying the result for this integral.</footnote>
We also recall that $v$ can be expressed in terms of $\lambda, \beta$ from (5.7), (5.8) which give the equation fixing it,

$$ \lambda = \frac{(\lambda(1 - \beta^2))^{1 - \frac{2v}{\pi}} (\pi(2v + \pi) \cos(2v) - \sin(2v) ((2v + \pi) \log(\lambda(1 - \beta^2)) + \pi))}{\pi^2}. \quad (5.23) $$

Thus we have two equations (5.19) and (5.22) which give $G$ and $h$ as functions of $u$ and therefore implicitly determine $G(h)$. It would be interesting to explore various critical regimes possibly hidden in this rather nontrivial function. We will recover in the next subsection the so called ‘flattening’ regime, leaving a more in-depth exploration for future work.

## 5.2 Critical regime and flattening limit

The critical regime dominated by large quadrangulations for this matrix model was studied in [29]. It was understood that it corresponds to the case when the derivative of $\partial\rho/\partial h$ vanishes at the endpoint $a$. This gives the following additional constraint on the parameters:

$$ 1 = \frac{\mathrm{cn}(u, k')}{\mathrm{sn}(u, k')\mathrm{dn}(u, k')} \Xi. \quad (5.24) $$

This equation allows one to express $\lambda$ (and all other parameters such as $v, b$ etc) as a function of $\beta$. As a result we have a 1-parametric critical line. Let us also mention that in the $q \to 1$ regime with exponential precision considered above in section 5.1, the condition (5.24) for criticality further simplifies and reduces to (see [29])

$$ \tan v = -\frac{\log(\lambda(1 - \beta^2))}{\pi}. \quad (5.25) $$

A particularly important part of the parameter space is the region when $\lambda \to 1$ and $\beta \to 0$. This also implies that we are close to the critical line (which passes through the point $\lambda = 1, \beta = 0$) though not necessarily directly on it. This is the limit discussed in [29] which we call ‘flattening’. As described there it corresponds to $\tau, v$ and $\beta$ all being small and of the same order $\sim \beta$. We find from (5.7), (5.8) that up to $\beta^2$ corrections we have in this limit

$$ q \approx \lambda \quad (5.26) $$

and

$$ v = \frac{\beta}{\sqrt{2}}(x - \sqrt{x^2 - 1}). \quad (5.27) $$

Here we introduced (following [29]) the variable

$$ x = 1 + \frac{\sqrt{2}(\lambda_c - \lambda)}{\pi\beta} \quad (5.28) $$

which is kept finite in our ‘flattening’ limit, and we also define

$$ \lambda_c = 1 - \frac{\pi\beta}{\sqrt{2}} \quad (5.29) $$
which corresponds to the critical value of $\lambda$ (to first order in $\beta$). It is also convenient to write (5.27) as

$$v = V\beta, \quad \text{with} \quad V = \frac{1}{\sqrt{2}}(x - \sqrt{x^2 - 1}) \qquad (5.30)$$

so that $V$ is a finite parameter which becomes equal to $1/\sqrt{2}$ on the critical line. We will write the results in terms of $V$ and $\beta$ in what follows. From the equations above we also find

$$\tau \approx \beta \frac{2V^2 + 1}{4V} . \qquad (5.31)$$

Let us remind that $\Lambda = 2(\lambda_c - \lambda)$ plays the role of bulk cosmological constant controlling the area (number of plaquettes) where as $\beta$ controls the curvature fluctuations.

In the next subsection we will discuss the properties of the solution for $G(h)$ and the resolvent $\tilde{W}$ in this regime.

## 5.3 Resolvent for Tr$(A_4\tilde{M})^n$ correlators in the flattening limit

The flattening limit affects the solution significantly as the cut structure becomes partially degenerate. From the results given in section 5.1.1 we find that the cut endpoints $a, b, c, d$ approach the values

$$a, b, c \approx 1, \quad d \approx -3 \qquad (5.32)$$

so that the $[a, b]$ cut collapses to a point and moreover almost touches the $[c, d]$ cut. The distance between the two cuts (i.e between $b$ and $c$) is exponentially small as we see from (5.6) that $b - c \sim e^{-\pi/\beta}$, while the size of the $[a, b]$ cut is of order $\sim \beta^2$.

Expanding our result (5.19) in this regime we get

$$G = -\frac{\sin r}{2V\beta} \left( 1 + \frac{\beta r(r - \pi)}{4\pi V} - \frac{\beta V (-r^2 + \pi r + 2(\pi - 2r) \cot(r) + 4)}{2\pi} + O(\beta^2) \right) \qquad (5.33)$$

and

$$h = b + \frac{4\beta^2 V^2}{\sin^2 r} \left( 1 - \beta\pi \frac{2V^2 + 1}{4V} \right) + O(\beta^4) . \qquad (5.34)$$

Here to make the result more compact we introduced the variable $r$ related to $u$ as

$$u = -ir - i\pi/2 . \qquad (5.35)$$

Our goal is to express $h$ in terms of $G$ and extract the resolvent $\tilde{W}$ from (4.7). From the structure of (5.33) we see that when $\tau \to 0$ we can keep finite the combination

$$y = 2\beta V G \qquad (5.36)$$

as well as $r$. Notice that in order to write the result for $\tilde{W}(P)$ we should take in (4.7) $G = 1/P^2$ which gives

$$y = 2V\beta/P^2 \qquad (5.37)$$

so that we keep $\beta/P^2$ fixed while $\beta, P \to 0$. Then writing

$$r = r_0 + \beta r_1 + \dots \qquad (5.38)$$
and expanding (5.33) at small $\beta$ we find $r(y)$ perturbatively, e.g. $r_0 = -\arcsin y$. Plugging the result then into (5.34), and using that in our regime (5.9) gives

$$b = 1 + \frac{2\beta^3 V (2V^2 - 3)}{3\pi} + O(\beta^4) \qquad (5.39)$$

we find to order $\beta^3$

$$h - 1 = \frac{4\beta^2 V^2}{y^2} + \frac{2\beta^3 V (2V^2 - 3)}{3\pi} + \beta^3 V \frac{(2V^2 + 1)(2\sin^{-1}(y)(\sin^{-1}(y) + \pi) - \pi^2)}{\pi y^2} \\ + \beta^3 \frac{8V^3 (\sqrt{1-y^2}(\pi + 2\sin^{-1}(y)) - 2y)}{\pi y^3}. \qquad (5.40)$$

Now we can write explicitly the resolvent itself, given by

$$\sum_{n=1}^{\infty} G^n \langle \frac{1}{N} \text{tr} (\tilde{M} B)^{2n} \rangle = \frac{1}{\sqrt{G}} \tilde{W} \left( \frac{1}{\sqrt{G}} \right) = h - 1 - \frac{T_2}{G} - \frac{T_4}{G(G-1)}. \qquad (5.41)$$

Namely, from (5.40) we find

$$\frac{1}{\sqrt{G}} \frac{\tilde{W}\left(\frac{1}{\sqrt{G}}\right)}{\beta^3} = \frac{\left(\frac{4V^3}{\pi} + \frac{2V}{\pi}\right) \sin^{-1}(y)^2}{y^2} + \left(\frac{4V^3 + 2V}{y^2} + \frac{16V^3 \sqrt{1-y^2}}{\pi y^3}\right) \sin^{-1}(y) \\ + \frac{8V^3 (\sqrt{1-y^2}-1)}{y^3} + \frac{4V^3}{3\pi} - \frac{16V^3}{\pi y^2} - \frac{2V}{y} - \frac{2V}{\pi} \qquad (5.42)$$

or equivalently$^{20}$

$$\frac{\tilde{W}(P)}{P\beta^2} = \frac{\arccos^2 y - \left(y + \frac{\pi}{2}\right)^2}{\pi y} \qquad (5.43) \\ + \left(\frac{P^2}{\beta}\right)^2 \frac{12y \arccos^2 y - 48\sqrt{1-y^2} \arccos y + 4y^3 - 3\pi^2 y - 48y - 24\pi}{24\pi}$$

where $V$ is defined in (5.30). This formula for the resolvent is one of our main new results. It gives the resummation of graphs with shape of a quadrangulated disc and the expansion of $\tilde{W}$ in powers of $G$ generates the moments $\frac{1}{N}\langle \text{Tr}(A_4 \tilde{M})^{2k} \rangle$, $k = 1, 2, 3, \dots$ in the flattening regime of large area and sparse conical defects scattered over the surface.

An important feature of the result (5.43) is the fact that the whole non-trivial dependence on $\lambda_c - \lambda$, $\beta$ and $G$ is hidden (apart from some trivial powers of $\beta G$) in the parameter $y = 2\beta V G$ which can be viewed as a renormalized boundary cosmological constant.

The result (5.43) is rather universal: it should be viewed as partition function of flat surfaces of disc topology with conical defects and the boundary represented by a zig-zag line consisting of the links of the boundary of original quadrangulations, as suggested by the image on figure 2.

<footnote>$^{20}$here we substituted $\arcsin y + \pi/2 = -\arccos y$.</footnote>
Let us now explore some properties of this result. In view of definition of correlators in (4.7) we are interested in the expansion of $h(G)$ in powers of $y$ (or $G$) when $y \to 0$. We find

$$h-1 = \frac{8V^3\beta^3}{y^3} + \beta^2 \left(1 - \pi\beta \frac{2V^2+1}{4V}\right) \frac{4V^2}{y^2} + \frac{2\beta^3V}{y} - \frac{1}{3}\beta^3V(V^2-1)y - \frac{2\beta^3V(6V^2-5)}{15\pi}y^2 - \frac{1}{20}\beta^3V(4V^2-3)y^3 + \dots \quad (5.44)$$

Plugging here $y = 2\beta VG$ we get

$$h-1 = \frac{1}{G^3} + \frac{1}{G^2} \left(1 - \pi\beta \frac{2V^2+1}{4V}\right) + \frac{\beta^2}{G} - \frac{2}{3}\beta^4 V^2 (V^2-1) G - \frac{8\beta^5 V^3 (6V^2-5)}{15\pi} G^2 - \frac{2}{5}\beta^6 V^4 (4V^2-3) G^3 + \dots \quad (5.45)$$

From the first few coefficients in the second line we read off, using (4.7), for instance,

$$\begin{aligned} \frac{1}{N}\langle \mathrm{Tr}(A_4 \tilde{M})^2 \rangle &= -\frac{2}{3}\beta^4 V^2 (V^2-1) + O(\beta^5), \quad (5.46) \\ \frac{1}{N}\langle \mathrm{Tr}(A_4 \tilde{M})^4 \rangle &= -\frac{8\beta^5 V^3 (6V^2-5)}{15\pi} + O(\beta^6), \\ \frac{1}{N}\langle \mathrm{Tr}(A_4 \tilde{M})^6 \rangle &= -\frac{2}{5}\beta^6 V^4 (4V^2-3) + O(\beta^7) \end{aligned}$$

and so on. Notice that while in (5.40) we know each coefficient of $y^n$ with accuracy of $\sim \beta^3$, when we substitute $y = 2\beta VG$ the accuracy of coefficients of $G^n$ will be different because $v \sim \beta$. In (5.45) we indicated all the information for the coefficients we can extract from (5.40).<fnref n="21" /> Below we will describe several nontrivial checks of the computation.

### 5.3.1 Checks of the result

Let us show that our result (5.43) for the resolvent passes several nontrivial tests. First, the part with negative powers of $G$ in (5.45) should match the $T_n$ coefficients as prescribed by (4.7), which read (recalling that $\epsilon = 1$ in our conventions)

$$h-1 = \frac{T_2}{G} + \frac{T_4}{G^2} + \frac{T_4}{G^3} + \frac{T_4}{G^4} + \dots \quad (5.47)$$

Since $T_4 \simeq 1$ we notice that coefficients of $1/G^n$ with $n \ge 4$ will not be visible in our computation, since when we translate them to $y = G/(2V\beta)$ they will be of order $\beta^n$ and thus lie outside of the $\sim \beta^3$ precision of our result (5.40). Indeed, in (5.45) we see that coefficients of all these terms are zero within our precision. As for the rest, recalling that

$$T_2 = \lambda\beta^2, \quad T_4 = \lambda \quad (5.48)$$

we see that the coefficients of the three terms in the first line of (5.45) perfectly match the prescribed form (5.47). That serves as a nontrivial consistency check of our result. Notice also that we find no $O(y^0)$ term in the r.h.s. of (5.40), which is another check of the result.

<footnote n="21">²¹that means that e.g. the coefficient of $1/G^3$ is given up to $O(\beta)$ corrections, the coefficient of $G$ is up to $O(\beta^2)$ terms, etc.</footnote>
As another test, we can compare the $G$ term with the prediction coming from the free energy $\mathcal{Z}$ that was computed for this model in [29]$^{22}$ and in our regime it is given by equation (5.10) from that paper,

$$ \mathcal{Z} = \frac{4t^4}{15\beta^2} \left(x^6 - \frac{5}{2}x^4 + \frac{15}{18}x^2 - \frac{5}{16} - x(x^2 - 1)^{5/2}\right) \qquad (5.49) $$

where $t = \sqrt{\lambda\beta^2}$ since we have set $\epsilon = 1$, while we recall that the scaling parameter $x$ from [29] is defined by (5.28). Let us note that, curiously, $\mathcal{Z}$ becomes simply a polynomial (divided by $\beta^2$) when written in terms of $V$ rather than $x$, i.e. plugging $x$ here from (5.30) we get

$$ \mathcal{Z} = \frac{t^4}{\beta^2} \left( \frac{V^6}{15} - \frac{2V^4}{15} + \frac{V^2}{12} \right) . \qquad (5.50) $$

We see that when expressed in terms of $t$ and $\lambda$ the couplings in (2.15) become $T_2 = \sqrt{\lambda}t$ and $T_{2q} = \lambda$ for $q \ge 2$. This means that

$$ \frac{1}{2N} \langle \mathrm{Tr}(A_4 \tilde{M})^2 \rangle = \partial_t \mathcal{Z}(t, \lambda) \qquad (5.51) $$

where we note that one should first write $\mathcal{Z}$ in terms of precisely $t, \lambda$ variables (excluding $\beta$) before differentiating it. Evaluating the derivative and translating the result to out parameters $V, \beta$ we reproduce the coefficient of the $G$ term in (5.45),

$$ \frac{1}{N} \langle \mathrm{Tr}(A_4 \tilde{M})^2 \rangle \approx -\frac{2}{3} \beta^4 V^2 (V^2 - 1) . \qquad (5.52) $$

This is another nontrivial test of our result.

### 5.3.2 Pure gravity and almost flat limits

In this subsection we will show that our result (5.43) obtained in the flattening limit interpolates between two interesting critical regimes of this disc partition function: 2d quantum gravity regime and ‘almost flat’ regime.

**Pure gravity limit.** The pure gravity limit corresponds to the case when $x \to 1$ as discussed in [29], and one can see that the free energy (5.49) has a degree 5/2 singularity there. Since $x \to 1$ we see from (5.28) that $\lambda \to \lambda_c$ with the critical value $\lambda_c$ given by (5.29) discussed above. It is also convenient to introduce the ‘bulk cosmological constant’ $\Lambda \propto \lambda - \lambda_c$ via

$$ x - 1 = \frac{\Lambda}{2} . \qquad (5.53) $$

Substituting $V$ in terms of $x$ into our result (5.42), we find several important features. First, plugging $y = G/(2V\beta)$ into (5.42) and expanding for $x \to 1$, we find that

$$ \tilde{W} = f_1 + (x-1)f_2 + (x-1)^{3/2}f_3 + \dots \qquad (5.54) $$

where $f_n$ are some lengthy functions of $\beta$ and $G$. While at first one may expect that the singularity of the resolvent for $x \to 1$ will be of the type $\sim (x-1)^{1/2}$ (for instance, $V$

<footnote>$^{22}$see equation (4.25) there which gives $\mathcal{Z}$ before taking any limit.</footnote>
itself scales as $V \simeq 1/\sqrt{2} - \sqrt{x-1} + ...$), we see that this is not the case and instead the singularity in (5.54) is $\sim (x-1)^{3/2}$ due to several nice cancellations. Thus all correlators $\langle \text{tr} (\tilde{M}A_4)^n \rangle$ also have a $(x-1)^{3/2}$ behavior. This nontrivial property is in complete agreement with the general prediction of [5, 6].

To study the behavior of the disc partition function in this limit, we will need to also send the parameter $y$ (that corresponds to the boundary fugacity) to a critical value, i.e. to a singularity of the resolvent. Notice that the resolvent (5.43) has an apparent branch point at both $y = -1$ and $y = 1$, but the first one actually cancels. Thus we will consider $y \to 1$ which we can also parameterise as

$$y = \sqrt{2V}(1 - \zeta). \qquad (5.55)$$

Here we introduced, in addition to $\Lambda$, the 'boundary cosmological constant' $\zeta$ with the normalization chosen for future convenience. Now, following [29] we consider the scaling limit when $\Lambda, \zeta \to 0$, keeping fixed the ratio

$$z = \frac{\sqrt{\Lambda}}{\zeta} \qquad (5.56)$$

which implies also $x, y \to 1$. Then we find that

$$\frac{\tilde{W}}{\beta^3} = \frac{-28 - 18\pi + 9\pi^2}{3\sqrt{2}\pi} + \frac{\sqrt{2}(-12 - 7\pi + 3\pi^2)\zeta}{\pi} - \frac{8}{3}\zeta^{3/2}(z-2)\sqrt{z+1} + O(\zeta^2). \quad (5.57)$$

The first two terms are regular in $\zeta$ and $z$ and are non-universal. We highlighted in blue the most interesting part, i.e. the 3rd term which is proportional to $\zeta^{3/2}(z-2)\sqrt{z+1}$. This perfectly matches the prediction of [10] for this universal part of the resolvent. That is another highly nontrivial test of our result.

**Almost flat limit.** Another key special case is the almost flat (AF) regime for which the resolvent we are considering was computed in [28]. We can recover that result by setting $V \to 0$ in our formula (5.42) for $\tilde{W}$ (or equivalently keeping only the $1/g$ term in (5.43)). This corresponds to the regime $x = 1 + \frac{\sqrt{2}(\lambda_c - \lambda)}{\pi\beta} \gg 1$, or $\beta \ll (\lambda_c - \lambda)$, i.e. the flattening parameter $\beta$ is much bigger than the bulk cosmological constant governing the area. In this regime we find

$$\tilde{W}(P) = \frac{\beta^2 P^2}{\pi y} \left( (\sin^{-1} y)^2 - y^2 + \pi (\sin^{-1} y - y) \right). \qquad (5.58)$$

Taking into account the difference in notation between our calculation and that of [28], we perfectly recover equation (4.29) of [28].$^{23}$

<footnote>$^{23}$Note a typo in (4.29) of [28]: the r.h.s. should have an extra $t_2$ factor. To compare (4.29) with our result we note that in our case $t_2 = \beta^2$ and for $V \to 0$ we see from (5.31) that $\tau \simeq \frac{\beta}{4V}$. Therefore the expansion parameter of [28] defined as $x_{\text{there}} = t_2G/(2\tau)$ is actually precisely our $y$. Then it's immediate to see that our result (5.58) is the same as (4.29) of [28].</footnote>
## 5.4 Resolvent for TrM$^n$ correlators

Here we will study the resolvent for a different kind of correlators — namely, the expectation values $\langle \text{Tr}M^n \rangle$. We recall that it is defined by (2.11),

$$W(P) = \frac{1}{N} \langle \frac{1}{P-M} \rangle \qquad (5.59)$$

and it can be computed from $H(h)$ via (4.9) which for convenience we repeat here:

$$W = P - \frac{h}{P} \qquad (5.60)$$

with

$$P^2 = he^H. \qquad (5.61)$$

Notice that for the purely gaussian model we would have $\rho = 1$ (corresponding to the empty Young tableau) and thus $H = \log \frac{h}{h-1}$ which gives

$$W = \frac{P - \sqrt{P^2 - 4}}{2} \qquad (5.62)$$

which, as was noticed already in [27], corresponds to the tree-like configurations stemming from pairwise Wick contractions of links of the boundary in a planar way. We see that $h$ has to stay finite and generic as it is related to the expansion parameter $P$ via (5.60). However, in the ‘flattening’ regime discussed above, $h$ is close to $b$ due to (5.34), as the r.h.s. of that equation contains the small parameter $\beta$. Moreover in that regime the cut $[a, b]$ of $H(h)$ collapses and becomes visible as simply a pole, so that inverting the function $H(h)$ (which we need to do in (5.61)) would only lead to a somewhat trivial result for the resolvent $W(P)$. This suggests that we have to look for a somewhat different scaling regime here.

In order to find a suitable scaling limit we will restore the parameter $\epsilon$ that we have set to 1 above. While the dependence of the correlators $\langle \text{tr} (A_4 \tilde{M})^{2n} \rangle$ we computed before on $\epsilon$ is rather trivial (see section 4.5), this is not the case for the averages $\langle \text{tr} (M)^{2n} \rangle$ which offer the possibility for a more nontrivial behavior. We found that it is natural to still take the ‘flattening’ limit when $\beta, v, \tau$ are all small but at the same time to send now $\epsilon$ to zero together with $\beta$ while keeping their ratio finite,

$$D = \beta/\epsilon = \text{fixed}, \quad \beta \to 0. \qquad (5.63)$$

The parameters $\tau$ and $v$ are taken to be of order $\beta$ like before (see (5.31), (5.30)) so that

$$v = V\beta, \quad \tau \approx \beta \frac{2V^2 + 1}{4V}. \qquad (5.64)$$

We will shortly see that in this regime $h$ naturally stays arbitrary as we wanted, rather than being close to $b$. The nontrivial scaling also introduces $D$ as an extra tunable parameter in addition to $V$ while keeping the computation similar in many ways to the one done above, and allowing us to use several results from section 5.1. We will be able to compute the
resolvent explicitly in this scaling regime which, as we will see, offers a variety of interesting features to explore.

Let us discuss how one can reinstate $\epsilon$ in various expressions. We first notice from (4.22), (4.23) that $v$ and $q$ do not have any $\epsilon$-dependence. Then from (4.24) we see that the dependence of $\Delta$ on $\epsilon$ is simply

$$ \Delta = \frac{\Delta_{\epsilon=1}}{\epsilon^2} \qquad (5.65) $$

and moreover

$$ X - 1 = \frac{(X - 1)_{\epsilon=1}}{\epsilon^2} \qquad (5.66) $$

with $X = a, b, c, d$. Thus we can use results from section 5.1.1 such as (5.6) and (5.9) which provide these parameters with exponential precision in $\beta$. For instance, at leading order we read off

$$ b \approx 1 + 2\beta D^2 V \frac{2V^2 - 3}{3\pi} \qquad (5.67) $$

$$ a \approx 1 + 4D^2 V^2, \quad d \approx -4D^2/\beta^2, \qquad (5.68) $$

while $b-c$ is exponentially suppressed for small $\beta$. We see an important feature that, unlike in the flattening limit, the cut $[a, b]$ of $H(h)$ remains of finite size in our regime.

These observations also mean that the expression for $H$, which is the most nontrivial component of the calculation, can be directly extracted from our solution with exponential precision in section 5.1 and we do not have to compute any new integrals. We find that

$$ H = \log \frac{h}{h-b} + 2\Phi + i\pi\rho \qquad (5.69) $$

where $\Phi$ is given by (5.20). Expanding it at small $\beta \sim \tau$ we find

$$ H = \log \frac{h}{h-b} + \left[ \frac{\pi\tau}{2} - \frac{2\tau u^2}{\pi} - 2i\tau u + \frac{2\beta V(-2 + (2u + i\pi)\tanh(u))}{\pi} \right] + \dots \qquad (5.70) $$

Finally, the relation between $h$ and $u$ reads

$$ h - b \approx 4 \frac{V^2 D^2}{\cosh^2 u} \qquad (5.71) $$

and we see that in contrast to (5.34) the r.h.s. does not contain any small parameters and $h$ can be kept arbitrary (with finite $u$ like before) which is just what we would like to have, as explained in the beginning of this section.

The relations we just described are enough to compute the resolvent, namely we express $u$ in terms of $h$ from (5.71), plug it into (5.70) and then solve (5.61) for $h(P)$ perturbatively in $\beta$. The result to order $O(\beta)$ reads

$$ W = \frac{L}{2D} + \beta \frac{16D^3}{L^2(4D^2 - L^2)} \times \left[ L^2 \frac{\sin^{-1}(\xi)^2 - \xi^2}{4\pi\xi} + \frac{\xi^3 + 12\sqrt{1 - \xi^2} \sin^{-1}(\xi) - 12\xi + 3\xi \sin^{-1}(\xi)^2}{6\pi} \right] \qquad (5.72) $$
where to make the formula more compact we introduced the notation

$$ \xi = VL, \quad L = D(P - \sqrt{P^2 - 4}). \qquad (5.73) $$

The formula (5.72) represents another important result of this paper. In the next subsection we will discuss its properties and limiting cases.

The first term in the r.h.s. of (5.72) is the additive contribution of tree-like configurations of the boundary, as mentioned before. Notice that all the bulk cosmological constant dependence is hidden in the parameter $\xi$, which can be considered as a renormalized boundary fugacity controlling its length. The factor $P - \sqrt{P^2 - 4}$ there comes from summing up the tree-like configurations of parts of the boundary, as well as $D$ controlling the curvature fluctuations, also enter that renormalization. Similarly to the parameter $y$ in (5.43), the $\xi$ dependence of the two terms in the 2nd line of (5.72) provides the most interesting information about the critical behaviors of disc quadrangulations with this boundary condition. Actually, these two formulas have a very close structure and we will discuss later in this section their similarities and differences.

### 5.4.1 Limits and properties of the result

The large $P$ expansion of the resolvent (5.72) reads

$$ PW = 1 + \frac{\frac{8\beta D^4 V^3 (5-6V^2)}{15\pi}}{P^2} + 1 + \frac{\frac{256\beta D^6 V^5 (10V^2-7)}{315\pi}}{P^4} - \frac{\frac{32\beta D^4 V^3 (6V^2-5)}{15\pi}}{P^2} + 2 + \frac{\frac{512\beta D^8 V^7 (14V^2-9)}{315\pi}}{P^6} - \frac{\frac{512\beta D^6 V^5 (10V^2-7)}{105\pi}}{P^4} + \frac{\frac{8\beta D^4 V^3 (5-6V^2)}{\pi}}{P^2} + 5 + \dots \qquad (5.74) $$

From the coefficients here we read off the values of the correlators to order $O(\beta)$ as

$$ \begin{aligned} \frac{1}{N} \langle \mathrm{tr} \, M^2 \rangle &= 1 + \frac{\frac{8\beta D^4 V^3 (5 - 6V^2)}{15\pi}}{P^2} \\ \frac{1}{N} \langle \mathrm{tr} \, M^4 \rangle &= 2 - \frac{\frac{256\beta D^6 V^5 (10V^2 - 7)}{315\pi}}{P^2} - \frac{\frac{32\beta D^4 V^3 (6V^2 - 5)}{15\pi}}{P^4} \end{aligned} \qquad (5.75) $$

and so on.

One important test is the behavior of the resolvent for $x \to 1$. Like for the case of $\tilde{W}$ discussed in section 5.3.2, we expect that the leading singularity should be $(x - 1)^{3/2}$ and not $(x - 1)^{1/2}$. Indeed, plugging $V$ as a function of $x$ into (5.72) and expanding it, we find

$$ W = g_1 + g_2(x-1) + g_3(x-1)^{3/2} + \dots \qquad (5.76) $$

where $g_k$ are some functions of $P, D$ and $\beta$. Thus we see that the $(x - 1)^{1/2}$ singularity non-trivially cancels and we have instead $(x - 1)^{3/2}$ as expected, in perfect agreement with the 2d QG behavior of the sphere partition function with a marked point [5, 6].

Just like the result in the flattening limit from section 5.2, our resolvent interpolates between the almost flat and pure gravity regimes. The AF regime corresponds to $V \to 0$ (equivalently, $x \to \infty$) when the resolvent becomes

$$ W = \frac{L}{2D} + \beta \frac{4D^3 \sin^{-1}(\xi)^2 - \xi^2}{4D^2 - L^2} \qquad (5.77) $$
We have also re-checked this result by repeating the derivation of the resolvent starting from the explicit solution given in [28] for the AF case.

Now let us consider the pure gravity limit. Like in section 5.3.2 we expect that it should correspond to sending the parameter $P$ to a branch point of the resolvent. Notice that $W$ has potential branch points at $P = \pm 2$ and also $P = \pm (1/(2VD) + 2VD)$ (i.e. $\xi = \pm 1$). As we are interested in the large $P$ expansion (that generates the correlators), the latter two branch points are the relevant ones as they are always closer to infinity. Moreover, similarly to the case of the previous resolvent, the branch point $P = -1/(2VD) - 2VD$ cancels. Thus we will take $P$ to be near the remaining branch point at $P = 1/(2VD) + 2VD$. Equivalently, we send $L \to \sqrt{2}$. Introducing a convenient normalization, we take the scaling similar to the one in section 5.3.2

$$L = \sqrt{2}(1 - \zeta), \quad x = 1 + \frac{\Lambda}{2}, \quad z = \frac{\sqrt{\Lambda}}{\zeta} = \text{fixed} \qquad (5.78)$$

where again we view $\Lambda$ as the (renormalized) bulk cosmological constant and $\zeta$ as the boundary one. Then we find for the nontrivial part of the resolvent the behavior

$$\frac{1}{\beta} \left[ W - \frac{L}{2D} \right] = \frac{(28 - 3\pi^2) D^3}{\pi (3 - 6D^2)} + \frac{D^3 ((6\pi^2 - 88) D^2 - 9\pi^2 + 100) \zeta}{3\pi (1 - 2D^2)^2} - \zeta^{3/2}(z-2)\sqrt{z+1} \frac{4\sqrt{2}D^3}{3(2D^2-1)} + O(\zeta^2). \qquad (5.79)$$

We see that the last term indeed reproduces the correct scaling function

$$\zeta^{3/2}(z-2)\sqrt{z+1} \qquad (5.80)$$

as expected from [10].

## 5.4.2 Comparison of two resolvents

Curiously, we observed that the two resolvents $\tilde{W}$ and $W$ that we have computed in sections 5.3 and 5.4 are quite closely related to each other. Namely, after a certain simple substitution and redefinition of the parameters they become equal up to a rational part. The latter may be viewed as not that important since it essentially does not affect the critical behavior, which is determined by branch point singularities.

Concretely, comparing (5.43) and (5.72) we observed that they become almost equal if we identify $y \leftrightarrow \xi$ and also make in the first one the formal replacement $\arccos y \to -\arcsin y$. The meaning of that last replacement is not fully clear to us at the moment, though it looks rather suggestive.<fnref n="24" /> To make the matching exact we also need to multiply one of the resolvents by an overall rational factor. As an outcome we find

$$(C_1 \tilde{W}(P)|_{\arccos y \to -\arcsin y, y \to \xi} - C_2 \left(W - \frac{L}{2D}\right)) = R_0 \qquad (5.81)$$

where

$$C_1 = \frac{1}{\beta^3 \sqrt{G}}, \quad C_2 = \frac{4D^2 L^2 - \xi^4}{2\beta D^3 L^3}, \quad R_0 = -\frac{L^3(\pi\xi + 8)}{\xi^6} - \frac{L(4\xi + \pi)}{2\xi^3}. \qquad (5.82)$$

<footnote n="24">²⁴Perhaps it can be viewed as a change of branch in some sense.</footnote>
The matching of all the parts including the arcsin and square root functions between the two resolvents is quite nontrivial and perhaps indicates that the result for the disc partition function shows a certain universality, regardless of how precisely we construct the boundary (as a tr $(A_4\tilde{M})^n$ or a tr $(\tilde{M})^n$ insertion). It would be interesting to understand more rigorously the reason for this matching.

The remaining different rational parts might be attributed to the finite (and generic) size of the boundary in our setup which introduces essentially lattice artefacts. It is worth reminding that both formulas (5.43) and (5.72) represent the disc partition functions of continuous flat manifolds with conical defects scattered over them, but at the same time with the boundary of any finite lattice size, i.e. consisting of any finite number of links. This sounds somewhat eclectic and needs some interpretation. Perhaps one may view such zig-zag boundary as a collection of conical defects stuck on it, making the continuous picture more consistent. However, at the moment all these possible interpretations are quite speculative and the relation (5.81) is just a nontrivial observation.

# 6 Conclusions and prospects

The main result of this work is the derivation of disc partition functions of abstract random 2D manifolds, flat everywhere except conical singularities (with the angle deficits $\Delta\phi = \pi, 0, -\pi, -2\pi, -3\pi, ...$) inserted at arbitrary positions and weighted with a certain fugacity controlling the fluctuations of $|\Delta\phi|$. Since the deficit of angle is akin to the insertion of a curvature defect, such fugacity controls the overall scale of the curvature fluctuations of these manifolds. The area of the manifold and the length of the boundary are weighted with the corresponding cosmological constants. For two different types of the boundary, we found that these partition functions are given by expressions (5.43) and (5.72). Due to this extra curvature fugacity, the formulas interpolate between the “almost flat” regime and the 2D quantum gravity regime. The “almost flat” regime arises in the limit when the curvature fugacity is such that the overall curvature fluctuations are highly suppressed. It corresponds to all the negative angle deficit (negative curvature) concentrated at the boundary of the disc, whereas the positive curvature defects represented by conical singularities with $\Delta\phi = \pi$ are scattered in the interior of the disc. The partition functions (5.43) and (5.72) behave slightly differently in this regime, i.e. they depend on the type of the boundary conditions. The 2D quantum gravity regime dominates in the limit opposite to the almost flat case, when the curvature fugacity favors the proliferation of curvature defects of both signs, i.e. the local metric of the manifold is highly fluctuating. In this regime, both formulas reproduce the well known universal answer for the disc partition function of pure 2d QG.

It would be also important to compute the sphere and disc partition functions for other types of conical defects via DWG. For example the explicit solution for $H(h)$ for triangulations, i.e. for the conical defects with angle deficits $\pm\pi$ was already written in appendix 1 of [29] in terms of contour integrals (see also [37]). However, it is rather involved and needs special efforts for its study. Whereas in the 2d QG regime the results should be still universal, the limit of almost flat surfaces should depend on the type of conical defects.
Interestingly, a somewhat similar problem was approached in the series of papers on Jackiw-Teitelboim gravity [18–21, 38], where one finds the partition functions (with or without boundaries) of the manifolds with constant negative intrinsic curvature $R < 0$ (i.e. Lobachevski, or $AdS_2$ metric), also with insertions of any number of conical singularities with arbitrary fixed $\Delta\phi$. It seems to be tempting to try to reproduce our results (5.43) and (5.72) from the results of these papers in the limit $R \to 0$, at least in the “almost flat” limit. However, it appears to be a very singular limit in JT gravity, so such a comparison is still an open question.<fnref n="25" />

On the other hand, one of our motivations for the current research was to construct and solve a DWG model that would imitate directly the sum over 2D manifolds with dominating $AdS_2$ background. At first sight, it seems to be easy: for example we can take in the DWG model (2.3) the following choice of the constants: $t_q^* = \delta_{q,4}$ (quadrangulations) and $t_q = g\delta_{q,4} + \gamma\delta_{q,6}$ (only zero or negative curvatures). This guarantees the $AdS_2$ background, at least on average, for large quadrangulations. However, the Euler theorem tells us that this can be achieved only with a boundary, i.e. for the disc, and all the positive curvature should be concentrated at this boundary. It is precisely the last condition which is difficult to achieve in DWG model. For both types of the boundary conditions, (5.43) and (5.72) contain both positive and negative curvatures in the bulk. It seems that to guarantee such a localisation of positive curvature on the boundary we have to compute a different disc partition function. One possibility is to take $W_L = \langle \frac{1}{N}\text{tr}(A_4^{-1}\tilde{M})^L \rangle$, with the same condition for quadrangulations $\frac{1}{N}\text{tr}A_4^q = \delta_{q,4}$. That means that near the boundary, instead of the squares we will have the triangles contributing the angles $\pi/3$. For example, two neighboring squares and one triangle give at such vertex the positive angle deficit $2\pi - (\pi/2 + \pi/2 + \pi/3) = 2\pi/3$. Consequently, we should be able to localize all such defects at the boundary. It would be fascinating to explore this direction further.

It would be also interesting to investigate in-depth the analytic structure of the disc partition function of quadrangulations in the exponential approximation (5.19), (5.22). It might contain new interesting critical regimes, inaccessible in the “flattening” approximation used in (5.43) and (5.72).

Another future direction is to properly formulate and explore the spectral curve for the generic DWG model. It would be very interesting to uncover what type of integrability arises in such matrix models, and also whether and how some kind of tau functions could appear here (even at finite $N$).

We hope to return to these questions in the future.

## Acknowledgments

We thank Ö. Gürdoğan for help with simplifying the integrals in section 5.1.2. We are also grateful to B. Basso, I. Kostov, A. Mironov, V. Mishnyakov, D. Serban, N. Terziev, G. Turiaci for discussions and comments.

<footnote n="25">²⁵We thank G. Turiaci for explanations on that issue.</footnote>
## A Derivation of the partition function and details on characters

Here we discuss some technical details for the derivations in sections 2 and 3, based on [27].

Using the completeness relation for characters, the integral in (2.1) can be expanded into characters of couplings $t_q$ and $t_q^*$ as follows:

$$ 
\begin{aligned}
Z(t, t^*) &= \sum_h \chi_{\{h\}}[t] \int d^{N^2}M \exp N \operatorname{tr} \left(-\frac{N}{2} M^2\right) \chi_{\{h\}}(AM) \\
&= \sum_h \frac{1}{d_h} \chi_{\{h\}}[t] \chi_{\{h\}}[t^*] \int \mathcal{D}M_i \Delta(M) \det_{k,l}(M_k^{hl}) \exp N \operatorname{tr} \left(-\frac{N}{2} M^2\right) \quad (A.1)
\end{aligned}
$$

where $M_i$ are eigenvalues of $M$. The integral then can be computed explicitly and is equal to

$$ \int \mathcal{D}M_i \Delta(M) \det_{k,l}(M_k^{hl}) \exp N \operatorname{tr} \left(-\frac{N}{2} M^2\right) = c(N) d_h \frac{\prod_i (h_i^e - 1)!! h_i^o!!}{\prod_{i,j} (h_i^e - h_j^o)} \quad (A.2) $$

The Schur characters with all but one vanishing couplings can be computed explicitly. For the character with the couplings $t_q = s\delta_{q,m}$ we have (see appendix 8 of [27]):

$$ \chi_{\{h\}}(s) = c(s)^{\frac{1}{m}} \sum_j h_j \left(\frac{N}{m}\right)^{\frac{1}{m}\sum_i h_i} \prod_{\epsilon=0}^{m-1} \frac{\Delta(h^{(\epsilon)})}{\prod_i \left(\frac{h_i^{(\epsilon)} - \epsilon}{m}\right)!} \operatorname{sgn} \left[ \prod_{0 \le \epsilon_1 < \epsilon_2 \le m-1} \prod_{i,j} (h_i^{(\epsilon_2)} - h_j^{(\epsilon_1)}) \right] \quad (A.3) $$

where $h^{(\epsilon)} = \epsilon \mod m$.

For our principal case of interest in this paper, $m = 4$, we have from here

$$ \chi_{\{h\}}(A_4) \sim \prod_{\epsilon=0}^{3} \frac{\Delta(h^{(\epsilon)})}{\prod_i \left(\frac{h_i^{(\epsilon)} - \epsilon}{m}\right)!} \operatorname{sgn} \left[ \prod_{0 \le \epsilon_1 < \epsilon_2 \le 3} \prod_{i,j} (h_i^{(\epsilon_2)} - h_j^{(\epsilon_1)}) \right] \quad (A.4) $$

where $h^{(\epsilon)} = \epsilon \mod 4$.

In order to compute the correlator $\langle \operatorname{tr} M^{2L} \rangle$ as discussed in section 3, we use the formula similar to (3.6)

$$ \operatorname{tr} M^L \cdot \chi_{\{h\}}[M] = \sum_{j=1}^{N} \chi_{\{h+L\delta_j\}}[M] \quad (A.5) $$

and, performing the gaussian integral over $M$ via (A.2) we obtain at the saddle point of (A.1) the result (3.17) given in the main text.

## B Mathematica notation for elliptic functions

Here we summarize the Mathematica notation corresponding to some elliptic functions we use in the text.

<table><thead><tr><th>Function</th><th>Mathematica command</th></tr></thead><tbody><tr><td>K</td><td>EllipticK[k<sup>2</sup>]</td></tr><tr><td>K'</td><td>EllipticK[1 - k<sup>2</sup>]</td></tr><tr><td>v = sn<sup>-1</sup>(√(a-c)/(a-d), k')</td><td>InverseJacobiSN[√(a-c)/(a-d), 1 - k<sup>2</sup>]</td></tr><tr><td>E(v, k')</td><td>EllipticE[ArcSin[JacobiSN[v, 1 - k<sup>2</sup>]]]</td></tr></tbody></table>
## C Details on the exponential approximation

Let us show that, as discussed in section 5.1.2, with exponential precision in $\tau$ we can do the following simplifications in the integral in r.h.s. of (5.14):

* Replace the upper integration limit by $\infty$

* Use the approximation (5.13) for $\rho$

* Use (5.15) for the quantity $J(w) \equiv \frac{dh}{dw} \frac{1}{h(u)-h(w)}$

To demonstrate the 1st point, we recall that $K \simeq \pi/(2\tau)$ with exponential precision, and for $w$ greater than this (large) value we can estimate the integrand as

$$ \rho \sim 1, \quad J(w) \sim e^{-2w} \qquad (C.1) $$

where the second estimate follows from (5.15). This gives for the part we are dropping

$$ \int_K^\infty dw J(w) \rho(w) \sim e^{-\pi/\tau} \qquad (C.2) $$

so it is indeed exponentially small.

Concerning the 2nd point, notice that (5.13) holds as long as $q'^2 \sinh 2u \ll 1$, i.e. up to values $u \sim 1/\tau$. For larger $u$ again we can replace $\rho \to 1$ in order to estimate the integral, and by the same argument as above we see that the region where (5.13) is not applicable gives an exponentially small contribution.

Finally to demonstrate the 3rd point, we compute the quantity $J(w)$ by starting from (4.19), writing it is a series in $k$ around $k=1$ and then looking at the large $w$ behavior. While at leading order we get (5.15), at higher orders we find potentially dangerous terms, e.g. at 2nd order we get a term that at large $w$ behaves as $\sim (k^2 - 1)^2 e^{2w}$. Since the integral goes from 0 to $K$ we can estimate its contribution as $\sim (k^2 - 1)^2 e^{2K}$. Although the second factor here is large and of order $e^{2K} \sim e^{\pi/\tau}$, it is still outmatched by the first one since $(k^2 - 1)^2 \sim e^{-2\pi/\tau}$ so in combination they give $\sim e^{-\pi/\tau}$, i.e. an exponentially small contribution again. One can repeat a similar estimate at higher orders and verify explicitly that any extra contributions are exponentially suppressed.

**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.

## References

[1] G. 't Hooft, A planar diagram theory for strong interactions, Nucl. Phys. B 72 (1974) 461 [INSPIRE].

[2] E. Brézin, C. Itzykson, G. Parisi and J.B. Zuber, Planar diagrams, Commun. Math. Phys. 59 (1978) 35 [INSPIRE].

[3] C. Itzykson and J.B. Zuber, The planar approximation. 2, J. Math. Phys. 21 (1980) 411 [INSPIRE].
[4] A.A. Migdal, Loop equations and 1/N expansion, Phys. Rept. 102 (1983) 199 [INSPIRE].

[5] F. David, Planar diagrams, two-dimensional lattice gravity and surface models, Nucl. Phys. B 257 (1985) 45 [INSPIRE].

[6] V.A. Kazakov, Bilocal regularization of models of random surfaces, Phys. Lett. B 150 (1985) 282 [INSPIRE].

[7] 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].

[8] V.A. Kazakov, Exact solution of the Ising model on a random two-dimensional lattice, JETP Lett. 44 (1986) 133 [Pisma Zh. Eksp. Teor. Fiz. 44 (1986) 105] [INSPIRE].

[9] V.A. Kazakov and A.A. Migdal, Recent progress in the theory of noncritical strings, Nucl. Phys. B 311 (1988) 171 [INSPIRE].

[10] V.A. Kazakov, The appearance of matter fields from quantum fluctuations of 2D gravity, Mod. Phys. Lett. A 4 (1989) 2125 [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] I.K. Kostov, Exact solution of the six vertex model on a random lattice, Nucl. Phys. B 575 (2000) 513 [hep-th/9911023] [INSPIRE].

[13] 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].

[14] B. Eynard and N. Orantin, Invariants of algebraic curves and topological expansion, Commun. Num. Theor. Phys. 1 (2007) 347 [math-ph/0702045] [INSPIRE].

[15] J.M. Maldacena, The large N limit of superconformal field theories and supergravity, Int. J. Theor. Phys. 38 (1999) 1113 [Adv. Theor. Math. Phys. 2 (1998) 231] [hep-th/9711200] [INSPIRE].

[16] E. Witten, Anti-de Sitter space and holography, Adv. Theor. Math. Phys. 2 (1998) 253 [hep-th/9802150] [INSPIRE].

[17] S.S. Gubser, I.R. Klebanov and A.M. Polyakov, Gauge theory correlators from noncritical string theory, Phys. Lett. B 428 (1998) 105 [hep-th/9802109] [INSPIRE].

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

[19] E. Witten, Matrix models and deformations of JT gravity, Proc. Roy. Soc. Lond. A 476 (2020) 20200582 [arXiv:2006.13414] [INSPIRE].

[20] G.J. Turiaci, M. Usatyuk and W.W. Weng, 2D dilaton-gravity, deformations of the minimal string, and matrix models, Class. Quant. Grav. 38 (2021) 204001 [arXiv:2011.06038] [INSPIRE].

[21] H. Maxfield and G.J. Turiaci, The path integral of 3D gravity near extremality; or, JT gravity with defects as a matrix integral, JHEP 01 (2021) 118 [arXiv:2006.11317] [INSPIRE].

[22] T.G. Mertens and G.J. Turiaci, Liouville quantum gravity — holography, JT and matrices, JHEP 01 (2021) 073 [arXiv:2006.07072] [INSPIRE].
[23] C.V. Johnson, *Nonperturbative Jackiw-Teitelboim gravity*, Phys. Rev. D **101** (2020) 106023, [arXiv:1912.03637] [INSPIRE].

[24] V.A. Kazakov, *Field theory as a matrix model*, Nucl. Phys. B **587** (2000) 645 [hep-th/0003065] [INSPIRE].

[25] S.R. Das, A. Dhar, A.M. Sengupta and S.R. Wadia, *New critical behavior in $d = 0$ large N matrix models*, Mod. Phys. Lett. A **5** (1990) 1041 [INSPIRE].

[26] P. Di Francesco and C. Itzykson, *A generating function for fatgraphs*, Ann. Inst. H. Poincaré Phys. Theor. **59** (1993) 117 [hep-th/9212108] [INSPIRE].

[27] 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].

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

[29] 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].

[30] V.A. Kazakov, M. Staudacher and T. Wynter, *Advances in large $N$ group theory and the solution of two-dimensional $R^2$ gravity*, hep-th/9601153 [INSPIRE].

[31] 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].

[32] A. Mironov and A. Morozov, *Superintegrability of Kontsevich matrix model*, Eur. Phys. J. C **81** (2021) 270 [arXiv:2011.12917] [INSPIRE].

[33] A.M. Vershik and S.V. Kerov, *Asymptotic theory of characters of the symmetric group*, Funkt. Anal. Pril. **15** (1981) 15.

[34] M.R. Douglas and V.A. Kazakov, *Large $N$ phase transition in continuum QCD in two-dimensions*, Phys. Lett. B **319** (1993) 219 [hep-th/9305047] [INSPIRE].

[35] D. Maître, *HPL, a Mathematica implementation of the harmonic polylogarithms*, Comput. Phys. Commun. **174** (2006) 222 [hep-ph/0507152] [INSPIRE].

[36] E. Panzer, *Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals*, Comput. Phys. Commun. **188** (2015) 148 [arXiv:1403.3385] [INSPIRE].

[37] R.J. Szabo and J.F. Wheater, *Curvature matrix models for dynamical triangulations and the Itzykson-Di Francesco formula*, Nucl. Phys. B **491** (1997) 689 [hep-th/9609237] [INSPIRE].

[38] D. Stanford and E. Witten, *JT gravity and the ensembles of random matrix theory*, Adv. Theor. Math. Phys. **24** (2020) 1475 [arXiv:1907.03363] [INSPIRE].