Double scaling limits of Dirac ensembles and Liouville quantum gravity

In this paper we study ensembles of finite real spectral triples equipped with a path integral over the space of possible Dirac operators. In the noncommutative geometric setting of spectral triples, Dirac operators take the center stage as a replacement for a metric on a manifold. Thus, this path integral serves as a noncommutative analogue of integration over metrics, a key feature of a theory of quantum gravity. From these integrals in the so-called double scaling limit we derive critical exponents of minimal models from Liouville conformal field theory coupled with gravity. Additionally, the asymptotics of the partition function of these models satisfy differential equations such as Painlev\'e I, as a reduction of the KDV hierarchy, which is predicted by conformal field theory. This is all proven using well-established and rigorous techniques from random matrix theory.


Introduction
Attempts to construct theories of Euclidean quantum gravity typically involve a partition function where the integration is over possible topologies or metrics and matter fields. However, these integrals are nonrenormalizable. Some approaches often used in physics involve discretizing the space, metrics, or topologies. Discrete approximations of physical theories have often found much success, such as in lattice gauge theory [31]. An alternative approach comes through noncommutative geometry. One may approximate commutative spaces by replacing the algebra of commutative functions on that space by a corresponding algebra of matrices. This was first done with the fuzzy sphere [33]. Such constructions exist for other spaces; for an example see the construction of the complex projective plane in [25].
As mentioned earlier, one would like to construct a path integral over all metrics (and eventually over matter fields), but when a space is "fuzzified" its metric loses its meaning. Alternatively, in the framework of spectral triples the role of a metric is played by the Dirac operator, as in Conne's distance formula [14]. Barrett first suggested in [4] that a toy model for finite noncommutative quantum gravity could be constructed as a well-defined matrix integral over an appropriate space of Dirac operators. We refer to these models as Dirac ensembles. The goal was that in some appropriate limit Dirac ensembles might connect to some understood physics, thus validating this random fuzzy approximation. The resulting partition functions of these models are matrix integrals, allowing one to use techniques from random matrix theory. In this paper we find that certain Dirac ensembles are dual to minimal models from conformal field theory coupled to gravity in the double scaling limit.
In random matrix theory the expectation values of observables can typically be written as a formal summation organized by genus, called the genus expansion, which was first discovered by t' Hooft [38]. Taking the large N limit of matrix models is equivalent to counting various types of planar maps [7]. As proven in [29], Dirac ensembles of any dimension have a genus expansion. In particular, the leading order term, found by a large N limit of appropriately scaled quantities, amounts to counting strictly planar maps. In the late 80's and 90's physicists found artifacts of conformal field theory in the large N limit of matrix models when coupling constants were tuned to specific critical values [12,20,16]. These critical values occur when the genus expansion terms of the log of the partition function fail to be smooth. This is analogous to how critical values are found in statistical mechanics. We give some intuition for this connection. Formal random matrix models are formal summations of Gaussian integrals, which via Feynman graph techniques can be realized as the weighted generating functions of maps [11]. A map is a type of embedding of a graph onto a two dimensional surface. Maps can be seen as a method of discretizing surfaces, and when the number of polygons that make up the map becomes very large, one would expect that one should be able to find a connection to partition functions that sum over surfaces. The number of polygons in maps can in fact grow by tuning the coupling constants to critical values [6]. Note that these critical values are not the same as those discussed in [5,28].
The authors wish to emphasize that none of these connections between quantum gravity and random matrices are new, but our goal is to reframe them as a toy noncommutative theory of quantum gravity. We prove in this paper that single trace matrix ensembles emerge from Dirac ensembles when the coupling constants are tuned appropriately. As we will discuss, this is highly non-obvious because the coupling constants of bi-tracial and single trace terms are shared. This allows one to recover random matrix models that are dual to minimal models of conformal field theory, an old and well studied connection in physics [6].
In Section 2 we give a brief introduction to Dirac ensembles. The Schwinger-Dyson equations and stuffed maps are developed in Section 3. The resolvent function is computed for some examples of Dirac ensembles in Section 4. Section 5 explains how one uses blobbed topological recursion to compute higher order and genus correlation functions. Section 6 highlights the main result of this paper. In this section the critical points of the examples in Section 5 are discussed as well as the connection to Liouville conformal field theory. In Section 7 we summarize these results and future projects.

Conclusion and Outlook 24
A Dirac ensemble is defined by a fixed finite dimensional real spectral triple, in which the Dirac operator is allowed to be randomly selected subject to some consistency conditions. Let us briefly recall that a finite real spectral triple is a quintuple (A, H, D, Γ, J) where A is a finite dimensional complex involutive algebra acting on the Hilbert space H. The Dirac operator is a self-adjoint operator from H → H. The two extra operators Γ and J, known as the charge conjugation and the chirality operator play no active role in our analysis so will not be discussed here. For a more detailed explanation see [15,4,39,34]. The partition function of a Dirac ensemble is given by Barrett and Glaser [5] as where D is the space of possible Dirac operators, which form a real vector space, and dD is the Lebesgue measure on this space. In general, the Dirac operator can be expressed in terms of gamma matrices tensored with commutators and anti-commutators of Hermitian and skew Hermitian matrices. More precisely, in [5] it was found that for any fuzzy spectral triple, the Dirac operator is of the form where the sum is over increasingly ordered multi-indices, and the following rules apply: where H I is some Hermitian matrix and {·} is the anticommutator.
• If γ I is skew-Hermitian, {K I , ·} e I = [L I , ·], where L I is some skew-Hermitian matrix and [·] is the commutator.
Note that the H I and L I are free variables. As a result, integral (1) can be expressed as an integral over the Cartesian product of the spaces of N ×N Hermitian matrices, H N , and skew-Hermitian matrices L N : where dD is the Lebesgue measure on the product space of finitely many spaces of Hermitian and spaces of skew-Hermitian matrices. In particular when p + q is equal to one or two we have The choice of S(D) is left open. However, we are particularly interested in cases where S is a polynomial in D. 1 A skew Hermitian matrix can be written as a Hermitian matrix multiplied by the complex unit. Furthermore, since S(D) is a polynomial, and thus the potential is a trace polynomial of the Hermitian and skew-Hermitian matrices seen in D, any integral of the above form can be written strictly as a Hermitian multi-matrix integral. These objects are interesting purely from a random matrix perspective, of which very little is known in general. Relatively recently some universal properties were established in [29].
Unlike the usual matrix model, in addition to the spectra of the H's we have the spectrum of D to study. As one might guess, there is a deep and not well understood relationship between them. The spectrum of the Dirac operator itself is not fully understood but displays some universal behaviors as seen in [29] and spectral phase transitions as seen in [5,28,24]. We are interested only in the simplest cases for now, partly because of the lack of analytical tools needed to study multi-matrix models.
We will consider a one dimensional Dirac ensemble of type (1, 0). We emphasize that one could also work with (0, 1) just as easily, such as in [28]. The type (1, 0) signifies that the associated Clifford module of the fuzzy spectral triple has a signature of one. The function S(D) is some polynomial in D and It is not hard to see that trace powers of D can be written as follows Thus, the integral is not just a matrix integral, but more specifically a bi-tracial matrix integral. We will discuss this further in the next section. We define a formal matrix integral as the formal summation formed by power series expanding all non-Gaussian terms in the integrand and swapping the order of integration and summation. Such formal sums of Gaussian integrals can be evaluated termwise using Wick's theorem [32,22]. Graphically, the coefficients of this formal sum can be realized as a weighted generating function counting stuffed maps [8].
We can define the matrix moments and higher moments of this ensemble as One can further define joint cumulants of higher moments using the classical moment-cumulant relations. For details see chapter one of [22]. These cumulants are denoted ⟨ 1 N n Tr H ℓ1 Tr H ℓ2 ... Tr H ℓn ⟩ c and are the generating functions of connected maps in the sense that the embedded graph is connected.
The Dirac moments can be computed from matrix moments in this case via formula (2). Higher order moments (which were not defined above) can be obtained in the usual manner from higher order cumulants. As proven in [3], the matrix moments and cumulants have a genus expansion, i.e.
where t is a continuous formal parameter strictly greater than zero. The terms of the genus expansion can be put into generating functions of the form The terms of this so-called genus expansion can be computed recursively using a process called blobbed topological recursion [3]. Blobbed topological recursion is a generalization of the similar process of topological recursion [21] which has gained much interest in the last two decades. For a review we refer the reader to [10].
In Section 3.2, we will derive a set of recursive equations that can relate cumulants and moments (as well as their generating functions) called the Schwinger-Dyson equations. It is a well-known practice in random matrix literature to compute W 0 1 using resolvent techniques. The beauty of (blobbed) topological recursion is that given W 0 1 and W 0 2 (which is often in some sense universal), one can compute any W g k recursively by decreasing Euler characteristic −χ = 2g − 2 + n. For single trace matrix models this process is well studied and formalized, see [22]. However, for multi-trace matrix models we are unaware of any similar reference. In this paper we will explicitly show how to compute W 0 1 and show that W 0 2 has a universal form for bi-tracial matrix models.

Bi-tracial matrix integrals
In this section we will set the groundwork for analyzing bi-tracial matrix models.

Stuffed maps
In this paper we are strictly interested in bi-tracial matrix models since they are the ones of interest for Dirac ensembles. However, one would expect that this analysis can be extended to higher trace multiplicity. Consider the following matrix integral over the space of Hermitian matrices where the trace of the potential is a bi-tracial polynomial in H where the t i 's and t i,j 's are coupling constants such that t i,j = t j,i .
It was first discovered in [7] that the moments and cumulants of the random Hermitian matrix ensembles coincide with the generating functions of maps. The maps of interest for bi-tracial matrix models are called stuffed maps and were first studied in [8] by Borot, and subsequently in [9]. Note that stuffed maps are a direct generalization of the maps in [22] used for single trace matrix models, thus all the following definitions simplify to the types of maps first considered in [7]. We now define the building blocks of stuffed maps.
An elementary 2-cell of topology (k, h) is a connected oriented surface of genus h with k boundaries. For example, a 2-cell with topology (1, 0) is a disc. These 2-cells can be "glued" together by pairing edges of the perimeter to form a surface with an embedded graph. The resulting surface is referred to as a stuffed map of topology (n, g) with perimeters (ℓ 1 , .., ℓ n ). It is an orientable connected surface with boundaries of lengths ℓ 1 , .., ℓ n . For more on stuffed maps see [8].
We are interested in enumerating these maps. More specifically, we want to count stuffed maps that are glued from 2-cells with the topology of discs and cylinders. This comes from the fact that the Dirac ensembles of interest are bi-tracial matrix models with appropriate N -powers as coefficients [8]. We shall refer to these maps as unstable stuffed maps. We define SM g k (v) as the set of all unstable stuffed maps of genus g, with v vertices and k boundaries. It was proven in [29] that this set is finite, allowing us to define the following formal series: where n ij (Σ) is the number of 2-cells with boundaries of lengths i and j used in the gluing of the map Σ. It turns out that these formal series are precisely the genus expansion terms mentioned in the previous section, that is For a detailed explanation and proof of this fact see [8].

The Schwinger-Dyson equations
The Schwinger-Dyson equations (SDE's) provide a powerful method to analyze random matrix models. They were first introduced by Migdal in [35]. In our context they are the consequence of matrix integrals of the total derivative vanishing. The SDE's of the bi-tracial matrix model (3) can be derived as follows. Consider the following relation, where H ℓ1 ij is the ij-th entry of the matrix power H ℓ1 , and the partial derivative ∂ ∂Hij satisfies One can prove the following properties, and Now, using the Leibniz product rule in (5) and the relations (7) and (8), we find that Tr H ℓm .
Note the model (3) can be considered as a formal matrix integral or a convergent matrix integral. Either way its moments will satisfy the nonlinear recursive relationship above. In the case of formal integrals, the equations are valid only order by order in t. Applying the genus expansions of moments and cumulants, then collecting terms of the same N/t powers, one finds the following.
The coefficient of N/t, when the number of boundaries n is equal to 1, is the main object of study in the following section.

The spectral curve
When g = 0 and n = 1, equation (9) becomes This is significantly more complicated than the single trace case. To simplify the notation we write it as where the bi-tracial terms become a combination of moments and coupling constants t j , all inside the newt j 's as follows:t Note that these new coupling constants are formal series in the old coupling constants since they are composed of the old coupling constants and T 0 ℓ for various ℓ, which are also formal series in the original coupling constants. Next we multiply equation (10) by 1/x ℓ+1 and sum from ℓ = 0 to infinity. The resulting equation is Then the equation can be written as where and By solving the quadratic equation (12), one can find the resolvent function Equation (12) is commonly referred to as the spectral curve. The solution W 0 1 (x), called the resolvent, can be used to give us the limiting eigenvalue density function of the random matrix associated to the given model. It is well-known that this relationship is given by the Stieltjes transform where we define ρ(y) as the limiting eigenvalue density function of the model. Thus, applying the inverse Stieltjes transform tells us that Consider for example the case that all coupling cosntants are zero except for t = 1. The resulting model is the Gaussian Unitary Ensemble and One can also derive that suppρ = [−2, 2]. So which is Wigner's semicircle distribution.

The 1-cut lemma
The following theorems for unstable stuffed maps are based on the results of [29] and is a generalization of a famous result by the same name discussed in [22]. Ideally one would like to know how the discriminant in equation (15) factors. This will determine the number of connected components in the support of the limiting spectral distribution and help with finding a nice closed form. Proof. Suppose we want to construct a planar map with one vertex and one of the 2-cells used in the gluing has the topology of the cylinder. Then since there is only one vertex both faces of the cylinder must be glued together in some manner to meet this requirement. Regardless of the configuration of the gluing, one will have created a handle by gluing the faces of the cylinder together. This contradicts the assumption that the map is planar. Hence, no planar map with one vertex can be glued from a 2-cell with the topology of the cylinder. Thus, an unstable planar map with one vertex must be glued from polygons only. However, we know from equation 1.2.1 of [22] that the number of vertices of such a map is equal to This completes the proof.
Brown's lemma tells us that the models we are interested in will always have a single cut solution. In general it is possible to find the number of connected components of support of ρ. This will depend on values of the coupling constants as well as the degree and structure of S(x). Note that the statement of this theorem differs slightly from the single trace case.
Proof. This proof follows closely that of Lemma 3.1.1 of [22]. Our initial goal is to first show that the zeros of our potential are up to some order in t close to the zeros of the discriminant S ′ (x) 2 − 4P 0 1 (x). What order in t is key to the existence of the single cut solution.
Recall definition (5.1.1) of T 0 ℓ . For ℓ = 0, T 0 ℓ = t. However, when ℓ ≥ 1, Lemma 4.1 is applicable, and says that T 0 ℓ is a formal series starting at a t 2 term. In particular, using definition (11) it is clear thatt 1 = O(t 2 ). Now recall the equation (13). Using the above information about T 0 ℓ , we may write it as where the last line follows from rearranging equation (14). Let x 1 , x 2 , ..., x d−1 be the roots of S ′ (x) −t 1 . For simplicity let x i be a root such that S ′′ (x i ) ̸ = 0. Let x be in the t-neighborhood of x i . Then we may write This allows us to write This shows that the zeros come in pairs [a i , b i ] centered around x i up to some order in t. The key question now is to what order. We have two cases to consider. Without loss of generality let x 1 = 0 (zero will always be a root since S ′ (x) −t 1 lowest order term is of degree two), then . In the first case this implies that the roots of S ′ (x) 2 − 4P 0 1 (x) can be written as and for the second case As candidates for variables of the same name in the statement of the theorem define α = 1 2 (a + b) and γ = 1 4 (a − b), then We now just need to show the last equality in the statement of the theorem for our choices of a and b. That is, we want to show that the discriminant S ′ (x) 2 − 4P 0 1 (x) only has one pair of simple zeros. This implies that we can pullout all the double zeros, giving us the polynomial M (x).
We know from Corollary 3.1.1 of [29], that the number of unstable planar maps with v vertices and a boundary of length ℓ is finite. Thus, W 0 1 (x) is a formal power series in t and may be written as where c v (x) is some Laurent polynomial of at least order two in x with coefficients formed in terms of the coupling constants. Thus, for any positively oriented contour C, integrating term-wise we know that Consider now a contour centered around one of the x i ̸ = 0, which for suitably small t contains a i and b i by equation (19). Let us compute the resolvent using (15) and Cauchy's integral theorem By the Taylor series expansion of the function Suppose that a i − b i ̸ = 0, then by the above expansion the integrand will have a pole at x i whose term is O(t 2 ). This contradicts equation (20).
Thus, we have shown that a i = b i for i > 1. In other words, for sufficiently small t the polynomial S ′ (x) 2 − 4P 0 1 (x) has only one pair a and b of simple zeros. This shows that the polynomials factors as M (x) 2 (x − a)(x − b) for some range of the coupling constants. This proves the last statement of the theorem and completes the proof.

The resolvent
In this section we will focus on finding the resolvent W 0 1 (x) for several formal bi-tracial matrix models that are type (1, 0) Dirac ensembles.

The Zhukovsky transform
The form of the discriminant in equation (17) motivates the use of the Zhukovsky transform The Zhukovsky transform maps the x-plane minus a line segment to the exterior of the unit disk in the z-plane. It also has the following useful identity We now borrow Theorem 3.1.1 from [22], which applies to our models as well.
Theorem 5.1. [22] For the formal power series α and γ 2 as mentioned above, we have the expansions with u 0 = 0 and u 1 = t/γ.
Notice that the theorem implies that since we have S(x), in theory we should be able to compute the coefficients of W 0 1 (x(z)) and then transform back to get W 0 1 (x). Computing the coefficients u k is much more involved in our models than in the single trace cases seen in [22]. We will show how to find them in various examples.

Moments
In general one may apply Lagrange's inversion formula to compute the moments and the support of ρ(x) for a given model. To see this one first observes that moments can be extracted from the resolvent generating function via the following contour integral We then apply the Zhukovsky transform to get Expanding one can find the general formula.
where the u j−i 's are the coefficients from theorem 5.1 that are determined by the potential.
This formula looks the same as for single trace models, but it is important to note that the u k 's contain other moments of the model because of the bi-tracial terms. Because of this, this formula now gives us a system of nonlinear equations to solve for moments.
As a side note, if one wishes to compute γ or α, the Lagrange inversion formula can be applied to recover them in terms of the coupling constants. Just like in the single trace case even though we treated moments as just formal series, they are in fact algebraic functions of α and γ, which are algebraic functions of the coupling constants. Thus, they have a finite number of singularities and are convergent in some sufficiently small disc. These singularities are a beautiful artifact of these models that is fundamental for the connection to Liouville quantum gravity, as we will see. The characterization of these singularities is not quite understood for multi-trace models in general, and may be an interesting phenomenon to study.
In the following sections we will set t = 1 but the calculations that follow can be carried out with t ̸ = 1.

The quartic model
Consider the potential Then the trace of the potential can be written in terms of the Hermitian matrix H: Note that we neglect terms with traces of odd powers due to their lack of contribution to the large N limit. The derivative of the potential in the spectral curve (12) is given by For this model T 0 ℓ counts the number of planar gluings of quandrangles and 2-cells with two boundaries of lengths two, with one polygon boundary of length ℓ. If ℓ is odd no such gluings exist so the generating function is zero. Hence, we focus on even values of ℓ.
We next transform the resolvent and find that Transforming back, we arrive at Thus, the limiting eigenvalue density function is One would like to be able to find the value of γ in terms of the coupling constants. To do this, we start by finding the transformed coefficients of S ′ : S ′ (x(z)) = (t 2 + 3T 0 2 t 4 )(α + γ(z + 1/z)) + t 4 (α + γ(z + 1/z)) 3 .
By the one-cut lemma as t goes to zero, so does α and T 0 2 by definition. Thus the factor ((t 2 + 3T 0 2 t 4 )) + t 4 (α 2 + 6γ 2 )) is nonzero when t 2 ̸ = 0 order by order in t. Hence α = 0 in order for this equation to hold on each term of the formal sum. This reduces the second equation to Using the formula for the second moment from Corollary 5.1.1, we have Thus we are left with an equation that relates γ to the coupling constants, which can be solved symbolically using software such as Maple or Mathematica. Note that if we travel along the curve t 2 = 1 − 3T 0 2 t 4 the derivative of the potential curve becomes which is the same as the quartic Hermitian matrix model See chapter 3 of [22] for details about this model. Thus, both models at these specifications have the same resolvent W 0 1 (x). This shall be important in later sections. We shall make similar conclusions about the rest of the models examined in the following two subsections.
By the same argument as in the quartic case, we can deduce α = 0. Thus, we find Using the Corollary 5.1.1, we find that This gives the limiting eigenvalue density function as where γ can be found as the solution to the 1/γ = u 1 relation for given t 2 , t 4 , and t 6 . Similarly as in the quartic ensemble, when one travels along the hypersurface 5 T 0 4 t 6 +3 T 0 2 t 4 +t 2 = 1, we have that the derivative of the potential becomes S ′ (x) = t 6 x 5 + 10 T 0 2 t 6 + t 4 x 3 + x.
Next, one may set 10 T 0 2 t 6 + t 4 equal to some desired value to reduce the derivative of the potential, and therefore W 0 1 (x), to being the same as the hexic Hermitian matrix model dH.

The cubic model
Consider the potential which gives us in the spectral curve equation (12). Thus, We then apply the Zhukovsky transform and find that Using Corollary 5.1.1, we may compute the moments of interest This gives the limiting eigenvalue density function where γ and α can be found as the solution to the above relations for given t 2 and t 3 . When one travels along 2 t 3 m 1 + t 2 = 1 and t 2 m 1 + t 3 m 2 = 0, the derivative of the potential is the same as the cubic Hermitian matrix model Thus, in these circumstances both models have the same W 0 1 (x).

(Blobbed) topological recursion
Roughly speaking, and in our context of matrix models, topological recursion works as follows. Using the resolvent technique one first defines a complex curve (Riemann surface) Σ, called the spectral curve of the model. One then constructs a sequence of symmetric mermorphic differential forms ω g,n (z 1 , . . . , z n )dz 1 ...dz n of degree n for g ≥ 0, n ≥ 1, on n-fold Cartesian products of Σ. Topological recursion works by induction on the Euler characteristic 2 − 2g − n of a surface of genus g with n boundaries. It gives an inductive formula for all ω g,n , starting with the first two forms ω 0,1 and ω 0,2 : where I = {z 1 , ..., z n } and β i are the ramification points of the ramified covering x, defined via dx(β i ) = 0. The recursive kernel K i (z, q) is constructed from the initial data. In the special case when Σ is the projective line, ω 0,2 is the Bergmann kernel ω 0,2 (z 1 , z 2 ) = dz 1 dz 2 (z 1 − z 2 ) 2 .
As previously mentioned, it was first seen in [8] that multi-trace matrix models correspond to the generating function of stuffed maps. Such generating functions obeys a generalization of topological recursion [8], known as blobbed topological recursion. Recall the generating functions W g k (x 1 , .., x k ) of unstable stuffed maps defined in Section 3.1. As discussed these generating functions are precisely the resolvent functions' genus expansion terms for multi-trace matrix models. In our case we are dealing with bi-tracial models, whose genus expansion is proven rigorously in [29] which corresponds to unstable stuffed maps i.e. stuffed maps involving 2-cells with the topology of discs and cylinders. Given the initial generating functions W 0 1 (x) and W 0 2 (x 1 , x 2 ), which count connected unstable stuffed maps with one and two boundaries respectively, one can use a recursive formula to compute all higher genus and boundary generating functions [8]. In particular for unstable stuffed maps, blobbed topological recursion reduces to the usual topological recursion. If more than the product of two traces occurs in the potential, then blobbed topological recursion is required. These computations are all done in Zhukovsky space.

Unstable stuffed maps of genus zero with two boundaries
We would like to use blobbed recursion to compute higher order correlation functions of Dirac ensembles. In the previous section we showed how to compute W 0 1 (x) with various examples. The next step is to find W 0 2 (x 1 , x 2 ), which for many important models, in Zhukovsky space, is universal. For example this is true in single trace Hermitian matrix models [22] and in the multi-matrix model seen in [13]. As far as the authors of this paper can tell, the form of this function has not been established for multi-trace matrix models, despite efforts in this direction given in [3]. For our bi-tracial models the proof of the universal form of W 0 2 (x 1 , x 2 ) is the same as the proof in Section 3.2 of [22]. That is Notice the appearance of the Bergman kernel.
If one wishes to compute the second order mixed moments, one needs to compute the residue To simplify this calculation one may apply the Zhukovsky transform in both variables and see that the second term on the right hand side of equation (22) contributes nothing to the residue.
In summary, we have established that for our bi-tracial matrix models W 0 2 (x 1 , x 2 ) has the same universal form as the single trace matrix models mentioned above.

Single trace models hidden in Dirac ensembles
Now with W 0 1 (x) and W 0 2 (x 1 , x 2 ) we may compute higher order correlation functions of our Dirac ensembles. This is proven as Theorem 9.1 in [3]. However, we shall not explicitly compute them here and instead have a different goal in mind.
As noted in the previous section, the coupling constants of the quartic, cubic, and hexic Dirac ensembles can be tuned to be the same as their respective Hermitian matrix model counterparts for certain values. When this is the case, W 0 1 (x) is the same, and in particular α and γ are as well. Thus, combining this fact with the universal form of W 0 2 (x 1 , x 2 ), via topological recursion all higher order genus expansion terms will be identical. Hence, we have proven that these single trace models hide in the above mentioned Dirac ensembles, at least for certain values of the coupling constants.
In [6] and chapter 5 of [22] it is proven that single trace Hermitian matrix models have interesting behavior at certain critical points. For the quartic model this occurs at t 4 = − 1 12 , the hexic at t 4 = − 1 9 and t 6 = 1 270 , and the cubic at t 3 = − 1 2 3 −3/4 . These points can be found as the locations of cusps of the spectral curve. See chapter 5 of [22] for details. Fortunately in the above mentioned Dirac ensembles all these critical points can be recovered. Their importance will be discussed in the following section.

The double scaling limit and 2D quantum gravity
In this section we will first review an old connection between random matrix theory and two dimensional quantum gravity. We will then discuss how these exact same connections hold for specific examples of Dirac ensembles of dimension one.

Large maps
As discussed in Section 3.1, formal matrix integrals count maps which are essentially polygonizations of Riemann surfaces. Intuitively, as the number of 2-cells that make up a map increases, the polygonization should give a better approximation of the underlying surface. Thus, our goal is to fine-tune the coupling constants of the model to some critical point that will cause the number of polygons that make up maps to go to infinity.
We emphasize that this is not a new idea. It was know on a heuristic level to physicists in the 80's and 90's [12,20] for asymptotic quantities of random matrix models. Physicists predicted a connection to Liouville conformal field theory coupled to gravity, and it was indeed proven using the KPZ formula [11], that many models from both theories have the same critical exponents. Correlation functions of certain conformal field theories have the symmetry of representations of conformal groups. Such infinite representations in two dimensions were classified in Kac's table [23], and distinguished by two integers (p, q). For each such integer pair of a so-called minimal model the generating functions must satisfy a partial differential equation. It was proved in [6] that formal single matrix models have an associated (2m + 1, 2) minimal model whose generating functions in the double scaling limit satisfy the associated partial differential equation and whose critical exponents match the minimal model. For example, the (3, 2) minimal model is referred to as pure gravity. Its generating functions satisfy Painlevé I. Both the quartic and cubic matrix models are associated with this minimal model. For a more detailed description see Chapter 5 of [22].
This idea is particularly useful in the context of Dirac ensembles because it provides a direct link between path integrals over metrics of finite noncommutative geometries to two dimensional quantum gravity coupled to conformal field theories. We will describe this idea here. Consider the genus expansion of the partition function for some Dirac ensemble [29] log Z : In general, F g is a function of the coupling constants with some algebraic or logarithmic singularities [22]. These F g 's can be computed from the W 0 k 's found using topological recursion. For details see Section 3.4 of [22]. As discussed in the previous section, for specific values of coupling constants the quartic, cubic, and hexic Dirac ensembles are precisely single trace Hermitian matrix models. Furthermore, all higher order correlation functions W 0 k are the same in these ranges of coupling constants. Thus, all F g are as well.
As the number of vertices gets very large, the behavior of the F g 's is going to be dominated by their singularities closest to zero. Thus, the behavior of unstable stuffed maps with a very large number of vertices is going to be controlled by the behavior of their generating functions near where its derivative diverges. Now consider the quartic Dirac ensemble as an example. We can fine tune its coupling constants such that we recover the quartic Hermitian matrix model. As seen in [22], the explicit forms of the F g 's of the quartic matrix model imply that near a critical point, the series has an expansion as the sum of regular and singular functions. With g ̸ = 1 at the critical point t 4 = 1/12, the singular part of the expansion looks like sing(F g ) = C g (t 4 − t c ) 5(1−g)/2 for some constant C g . When g = 1, If one defines the new series then u ′′ (y) satisfies the Painlevé I equation y = (u ′′ (y)) 2 − 1 3 u (4) (y).
Note that, even in the double scaling limit, the generating functions of large maps can be computed using topological recursion [22]. Furthermore, for general single trace matrix models, we can construct the following formal series: whereF g are the leading terms in the asymptotic expansions of F g . This series is a formal Tau-function of some reduction of the KdV hierarchy, and can be obtained from Liouville conformal field theory when coupled with gravity. For more details see Section 5.4 of [22]. Thus, if we can fine-tune a Dirac ensemble to the critical points of a single trace matrix model, the same result follows. In particular, we know that the critical point of the quartic Hermitian matrix model is t 4 = −1/12, and considering the analysis done in Section 5.3 we have the following phase diagram for the quartic Dirac ensemble. In the above diagram we have the curve where the spectral phase transition happens. This is precisely when the constant t 2 is chosen in terms of t 4 as This can be found by finding when ρ(0) = 0, just as in [28]. One can also conclude that the quartic matrix model appears when In the quartic matrix model the critical point occurs at t 4 = − 1 12 . In the quartic Dirac ensemble this corresponds to the point t 4 = − 1 12 and t 2 = 4 3 . Note that the spectral phase transition curve does not cross the quartic matrix model curve. This is because, in the quartic matrix model considered, there is no spectral phase transition. A coupling constant in front of the Gaussian term is required for this phenomenon to occur.
Remark. Ultimately we are solving the loop equations, which are the same in this case regardless of whether the model is considered formal or convergent. Which interpretation we can choose is dependent upon which quadrants of Figure 2 we are considering. In the first and fourth quadrants t 4 is negative in the potential, so the model can always be seen as convergent since the quartic term dominates the Gaussian one at the limits of the integral. In quadrants one and two the model can always be viewed as formal since the Gaussian integrals in the definition of a formal integral are always convergent. In the third quadrant, however, the model is neither formal nor convergent, as the integrals diverge, but so do all the Gaussian integrals in the formal definition.
A similar analysis may be done for the cubic Dirac ensemble. However, no spectral phase transition occurs in this model. See 3.  is associated with a (2m + 1, 2) minimal model in the double scaling limit. In particular, we can deduce then that the quartic, cubic, and hexic Dirac ensembles of type (1, 0) are associated with the (3, 2), (3,2), and (5, 2) minimal models. In general for single matrix models, if the potential of the model is degree d and it is odd then it is associated with the (d, 2) minimal model. If the degree is even, then it is associated with the (d − 1, 2) minimal model. We expect the same relationship holds for Dirac ensembles of the above mentioned form.

Conclusion and Outlook
In this paper we analyze the quartic, cubic, and hexic Dirac ensembles of type (1, 0) as formal matrix integrals. Their resolvent functions W 0 1 and limiting eigenvalue distributions are found explicitly. The cylinder amplitude W 0 2 is discussed to have a universal form. Thus, via the process of blobbed topological recursion one may compute all higher genus and boundary correlation functions. During this analysis, it is found that by fine-tuning the coupling constants of these models one can recover critical phenomena seen in certain Hermitian matrix models. In particular in the double scaling limit we find the critical exponents associated with minimal models from conformal field theories. Additionally, for these models the genus expansion terms of the log of the partition function satisfy the same differential equations as the partition functions of the corresponding minimal model in the double scaling limit. We hope to prove rigorously in future work that most Dirac ensembles of type (1, 0) and (0, 1) have a corresponding minimal model in the double scaling limit. Further, it would be interesting to construct multicritical matrix models [27,1,2] from Dirac ensembles.
Essentially, we have recovered conformal field theories coupled to gravity from toy models of quantum gravity on noncommutative spaces. We would like to extend this connection to Dirac operators of higher dimension. This seems likely considering that similar connections have been made for multi-matrix models [16]. Without any known analytic techniques to study matrix models seen in higher dimensional Dirac ensembles we can only aim for numerical evidence. It may be possible to deduce the critical values of Dirac ensembles using the bootstraps technique [26]. In random matrix models the range of coupling constants on which the model is defined determines the critical values. For example the quartic matrix model has solutions when −1/12 < t 4 in the large N limit. Thus, bootstrapping seems like a likely candidate to use to find these critical points. However, determining these critical exponents might be better suited to Monte Carlo simulations. This was explored in some sense in [24]. Note that the finite size of N might strongly affect these values. Another possibility is through functional renormalization group techniques [37]. We hope to explore these ideas in future works.
We would like to thank the referees of this paper for their very thorough feedback and suggestions that have greatly improved this paper from its original draft.