Coloured combinatorial maps and quartic bi-tracial 2-matrix ensembles from noncommutative geometry

We compute the first twenty moments of three convergent quartic bi-tracial 2-matrix ensembles in the large $N$ limit. These ensembles are toy models for Euclidean quantum gravity originally proposed by John Barrett and collaborators. A perturbative solution is found for the first twenty moments using the Schwinger-Dyson equations and properties of certain bi-colored unstable maps associated to the model. We then apply a result of Guionnet et al. to show that the perturbative and convergent solution coincide for a small neighbourhood of the coupling constants. For each model we compute an explicit expression for the free energy, critical points, and critical exponents in the large $N$ limit. In particular, the string susceptibility is found to be $\gamma =1/2$, hinting that the associated universality class of the model is the continuous random tree.


Introduction
When attempting to define theories of Euclidean quantum gravity, one is usually interested in making sense of path integrals over some class of Riemannian metrics.In the context of noncommutative geometry, spectral triples are analogous to manifolds, and in some sense spectral triples generalize manifolds.In particular, for Riemannian spin c manifolds all the metric information can be recovered from the associated spectral triple via Connes' distance formula [16].With this in mind, Barrett [3] proposed defining path integrals over the moduli spaces of Dirac operators instead of those of metrics.To make these integrals well-defined, Barrett considered finite approximations of spectral triples called fuzzy geometries.The resulting integrals are matrix integrals.The hope after fuzzifying these spectral triples is that in some limit one might be able to recover path integrals over metric spaces of Reimannian spectral triples.Thus, in some sense, this would recover a theory of Euclidean quantum gravity.The limits of finite approximations of spectral triples is an active area of study [56,57,17,46].
While the work surrounding these models is not quite at this point of development, we have seen hints of the continuous theory in various limits.As originally pointed out in [6], in the large N limit the spectral density function of the Dirac operators of certain Dirac ensembles bears resemblances to the spectra of Dirac operators on spin manifolds.This was explored quantitatively in [4].More recently, in the double scaling limit, various Dirac ensembles have been shown to have the same critical exponents and satisfy the same differential equations as various minimal models from conformal field theory [38].One approach to this problem might be to consider random metric spaces of maps defined by the perturbative expansion of these models.Based on the critical exponents found in [38], these random metric spaces are expected to converge to the Brownian map [47], but this will be explored in a future work.This is in contrast to the models studied in this paper, for which we find strong hints that the associated random metric spaces converge to the continuous random tree.
The models of interest in this article are three 2-matrix bi-tracial ensembles proposed by Barrett and Glaser [6].These models have been studied numerically in [6,32,4,18,33] via Monte Carlo simulations, which provided evidence of a spectral phase transition.Additionally, analytical bounds for the moments of these models were computed in [37] using the bootstrap technique.We will show that all three ensembles have the same moments in the large N limit and that one only needs to consider the following effective ensemble: 1 Z e −S eff (A,B) dAdB, where this is a probability distribution on the space of pairs of N by N Hermitian matrices.The measure dAdB is the product Lebesgue measure, and the potential is The probability measure depends on two real coupling constants t 2 and t 4 , where t 4 > 0. These models are interesting purely from the perspective of random matrix theory.Despite the success of single matrix ensembles [51,23], in general very little is known about multi-matrix ensembles.For perturbative models, much is known for general potentials in which the only unitarily invariant term is a Tr AB interaction [29,15].Besides this, only special cases of more general multi-matrix interactions are known [42,25].In convergent models even less is analytically tractable [13,34].However, progress has been made in determining when the asymptotics of convergent models coincide with their perturbative counterparts [36,7].Such results will be utilized in this work, allowing us to work with a perturbative expansion but make conclusions about convergent matrix integrals.An additional level of complexity in the matrix integrals proposed by Barrett comes from the fact that they are bi-tracial.As far back as the 90's physicists were already interested in studying integrals known as multi-trace in which the potential function S(M ) includes some product of traces.Multi-tracial matrix integrals have appeared in many other areas of study [1,31,61,59], and have seen the development of tools for both the perturbative and convergent cases [9,10,19].Despite Barrett's models having highly nontrivial 2-matrix interactions in their potentials, in this article we derive explicit formulae for the first twenty moments in terms of the coupling constants t 2 and t 4 .The idea of the derivation is to first consider the perturbative expansion of the models and study them as the generating functions of certain types of maps.Certain properties of these maps and the associated Schwinger-Dyson equations will allow us to deduce the moments.Applying a result of [36] shows that these formulae are equivalent to the moments of the convergent ensembles in some neighbourhood of the coupling constants near zero.From these formulae for the first few moments we deduce the free energy, critical points, and critical exponents of the model.Note that partition functions of these models are the same whether they is written in terms of D or A and B, so even though we are working with moments in terms of A and B, instead of D, we will be able to compute the above mentioned quantities of interest.Work towards explicit formulae for tracial powers of the Dirac operators of fuzzy geometries in terms of the tracial powers of constituent matrices can be found in [54].It is out hope that these results will lead to more analytic results as well as more direct approaches to studying these models in future works.
In the following section we introduce the necessary background from Noncommutative Geometry and state the main results.In Section 2, we outline the derivation of the Schwinger-Dyson equations for the model and its properties.Section 3 gives a short review of the relevant kind of maps before proving several useful properties that allow for the derivation of the moments.Section 4 shows the computation of the free energy as well as the critical exponents.Section 5 outlines the future work and implications of our results.

Random fuzzy geometries
In [16], Connes introduced the notation of a spectral triple (A, H, D) in which • A is a unital, involutive, complex, and associative algebra.
• The complex Hilbert space H is acted on by elements of A.
• The Dirac operator D is a self-adjoint operator acting on H, that is in general unbounded.
These objects are additionally required to satisfy some regularity conditions.However, we are interested in spectral triples that automatically satisfy these conditions, so such details will be omitted.In particular, we are interested in real spectral triples, which have even more additional structures.The motivation to study real spectral triples is that they serve as noncommutative analogs of spin c Riemannian manifolds.This idea is based on the fact that any closed spin c Riemannian manifold M gives rise to a real spectral triple in which the algebra A = C ∞ (M ) is the algebra of smooth complex valued functions on M and the Hilbert space is the space of square integrable sections of the spinor bundle such that the elements of A act as multiplication operators.The Dirac operator D is the usual Dirac operator of M , and acts on the spinors.The additional structure mentioned before consists of standard charge conjugation and chirality operators, J and γ.Conversely, the reconstruction theorem of Connnes tells us that under some natural conditions a real spectral triple with a commutative algebra can be realized as the real spectral triple of a spin c Riemannian manifold [16].
Fuzzy spaces have been studied as a method of regularization of commutative spaces since the fuzzy sphere in [50].In particular, they can be characterized within the formalism of spectral triples and are called fuzzy geometries or fuzzy spectral triples [3,5].From a physics perspective these can be thought of as spin c Riemannian manifolds with a finite resolution or Planck length.
A (p, q) fuzzy geometry is a real spectral triple of the form (M N (C), V ⊗ M N (C), D; J, Γ) in which • The algebra of functions is replaced by the algebra of N by N complex matrices.
• The Hilbert space is some Hermitian irreducible Clifford module of signature (p, q) with the charge conjugation operator J and grading Γ when the KO dimension p + q is even.
• D is a self-adjoint matrix that satisfies the so-called zero order and first order conditions [3].
A result from [3] is that all Dirac operators satisfying the above-mentioned conditions can be expressed as where the sum is over index sets of the form {i 1 ≤ ... ≤ i k } with each index between one and p + q.The operator γ I denotes some product of gamma matrices.If γ I is Hermitian, e I = 1 and {K I , •} e I = {H I , •}, where H I is some Hermitian matrix.If γ I is skew-Hermitian, e I = −1 and {K I , where L I is some skew-Hermitian matrix.One can deduce from this result that the space of possible Dirac operators D is isomorphic as a real vector space to some Cartesian product of copies of the spaces of N by N Hermitian matrices H N , and N by N traceless Hermitian matrices H 0 N .In the large N limit, these traceless Hermitian matrix ensembles have the same moments as their Hermitian counterparts [18].Hence, since we are currently only interested in the large N distribution of these models, we will strictly consider Hermitian matrices.
With quantum gravity as a motivation, it makes sense to then consider a probability distribution on D called a Dirac Ensemble (or sometimes a random fuzzy geometry).The usual probability distributions of choice are of the form where V is some polynomial with coupling constants as coefficients such that the probability distribution is well-defined.In [45], Gaussian Dirac ensembles were studied extensively and found to have universal properties in the large N limit.However, note that the main choice of potential in most works has been a quartic action since it has been seen to already exhibit many interesting properties.In particular, quartic Dirac ensembles of this form exhibit manifold-like behaviour near spectral phase transitions [6,32,4,33].
If an additional coupling constant is considered in front of the quartic term, when tuned to criticality, such models have been found to have connections to the (3, 2) and (5, 2) minimal models from conformal field theory [38].Until this work, explicit analytical progress on such models has mostly been on Dirac ensembles with only one Hermitian matrix [44,45,38,58], with the notable exception of [55] and bounds on moments obtained in [37].For more details we refer the reader to the recent review [39].

Outline of main results
In this paper the main objects of study include the following Dirac ensemble for signatures (2, 0), (1, 1), and (0, 2): Its partition function is where t 2 and t 4 are some real coupling constants.Originally in [3,6], a key aspect of these models is that they are convergent matrix integrals and no perturbative expansion or renormalization techniques are required to make them well-defined mathematical objects.However, in this work we shall consider both the convergent models and their perturbative counterparts.A formal matrix integral is a well-defined formal series defined by series expanding all non-Gaussian terms in the potential and then interchanging summation and integration.These are vastly different mathematical objects that historically have caused confusion, but have a deep relationship.For more details see [2,23,36,25].In particular, we will show that in the large N limit the loop equations for these models are the same for both the formal and convergent models, and have a unique solution.We will denote expectation with respect to a formal matrix integral with bra and ket, and with respect to a convergent matrix ensemble with E[•].
The Dirac operators on these fuzzy geometries can be expressed as where A and B are N × N Hermitian matrices, and σ 1 and σ 3 are the Pauli spin matrices The commutators and anti-commutators are represented as matrices Expressing the action in terms of A and B gives us where the epsilons are signs that change depending on the signature of the fuzzy geometry according to Table 1.
Table 1: Different signs in the action correspond to different KO dimension [3].
In this paper, all results will be to leading order in the large N limit.As pointed out in [45,37], many of the terms in (4) do not contribute to the leading order loop equations.As such we can consider a simplified model whose action we will refer to as the effective action, but who has the exact same large N behaviour as the above models: Notice that all the epsilon terms are not included.This serendipitously implies that to leading order in N all the models have the same large N behaviour.
The model is a bi-tracial two-matrix ensemble.In random matrix theory one is generally interested in computing moments and more generally correlation functions.The Dirac moments are defined as follows: 1 for integers ℓ ≥ 0. Let W belong to the set of noncommutative polynomials in two matrix variables C[A, B], then the (mixed) moments are defined as Tr W e −S eff (A,B) dAdB.
Note that A and B are symmetric in the potential, which implies the equivalence of many moments.Finding moments at finite N is very difficult, but generally computing their limits as N goes to infinity greatly simplifies calculations, provided the limit exists.Much success has been achieved in this direction dating back to Wigner [60], and his successors [11].Additionally, many universal properties have been observed in the limit [21].
For unitary invariant ensembles, one can apply techniques such as the Coloumb gas/equilibrium measure approach [19,20] or in the case of a formal integral one can apply (Blobbed) Topological Recursion [30,9,10].Analytic progress has been achieved for some models that lack invariance [28,24,27], but for general potentials very little is known [25,36].Our model is particularly challenging.It clearly lacks unitary invariance, techniques such as the characteristic expansion [43], the Harish chandra formula [40], and bi-orthogonal polynomials [8] are not applicable.Numerical studies of these particular models have been carried out and many interesting properties have been found [6,32,4].In particular in [37] the bootstrap technique was applied to find explicit bounds for moments of these models in the large N limit.Numerical estimates for the moments were then obtained.While an explicit formula for any moments in terms of coupling constants escaped us at the time, this paper presents such a formula.
Theorem 2.1.The formal and convergent models of (3) for all three signatures have the same limiting moments.In particular, for t 2 and t 4 in a sufficiently small enough neighbourhood of zero, lim The proof is presented in Section 4.3.The idea of the proof is to first consider the formal counterpart of the model and prove such a claim using Feynman graphical techniques.Then, applying results from [36], we can conclude that the loop equations for both formal and convergent models have a unique solution.A list of explicit formulae for higher power moments can be found in Appendix B. We conjecture that from the second moment all other moments and Dirac moments can be computed explicitly using Schwinger-Dyson equations.
Another quantity of interest in random matrix theory is the so called free energy, If this limit exists, as a formal series, it counts some collection of colored planar maps.This limit does indeed exist for the formal model: see Section 4.
Theorem 2.2.The formal models of (3) for all three signatures have the same free energy given by To someone familiar with the moments of matrix models, it may appear strange why these formulae are simpler than most single matrix hermitian models.Consider, for example the moments for the (1, 0) quartic model [38].The reason for these more concise expressions is that the number of maps enumerated in these models is generally smaller than its single matrix cousins.This is because there are more complicated 2-cells used in the gluing and we are restricted to gluing edges that match in color.

The Schwinger-Dyson equations
The Schwinger-Dyson equations are an infinite system of non-linear recursive equations of moments that were first discovered in [52].They can be derived from very simple principles but can be used to deduce many properties of matrix ensembles [35].Processes used to solve matrix models often rely on these equations, such as topological recursion [30] and bootstrapping [48,41].

Derivation and properties
Let W ∈ C⟨A, B⟩.The following equality follows from Stokes' theorem It is important to note that this equality holds for both convergent or formal matrix integrals.In the formal case, it is applied term-wise to each Gaussian integral in the formal series.By expanding the left-hand side using the product rule one can obtain relations between (mixed) moments.For example, suppose that W = A ℓ for some integer ℓ ≥ 0. Then equation ( 8) becomes Such relations are called the Schwinger-Dyson equations (SDE), since, unlike the usual Schwinger-Dyson equations found in single matrix models, the matrices involved may not commute, resulting in a much more vast ocean of relations to solve.Usually, one considers the generating functions of these moments to allow complex analytic techniques to solve this infinite system [26,35].However, with this model there is no clear choice of generating functions that allow for nice closed-form expressions for the SDE.Thus, in this work we are grounded to work to the level of (mixed) moments.The authors have yet to find a formula for these SDE that is concise but also informative.
In the large N limit the SDE often simplify.In particular, the factorization property lim is exploited when possible for W 1 , W 2 ∈ C⟨A, B⟩.In formal Hermitian matrix models this property follows from the genus expansion [26].From theorem 3.1 of [45], models such as (3) have a genus expansion and hence satisfy this property.For details on how the genus expansion implies this property see the appendix of [37].In the large N limit we introduce the following notation for the (mixed) moments of the convergent ensemble m ℓ1,ℓ2,...,ℓq = lim and m 0 ℓ1,ℓ2,...,ℓq = lim for the formal ensemble.This notation is well-defined, since the model is symmetric in A and B.
We are interested in the SDE in the large N limit.For example, in the formal case, after normalizing equation (9) and taking the limit, the result is for integer ℓ ≥ 0.
For any choice of initial word, such equations can be deduced.For more examples of these equations for general words, see Appendix C.

From perturbative expansion to convergent integrals
As mentioned in the introduction, our strategy is to solve the formal model corresponding to (22), and then use known results to relate the solution to its convergent counterpart.To do this, we will use the results in [36] and adapt them for our bi-tracial model.
Consider the following formal matrix model where Lemma 3.1.Up to the leading order in N , the loop equations for the model ( 11) and the model ( 5) are the same. Proof.
Consider the potential V as a map Assume t 2 > 0 and t 4 ≥ 0. The first line of equation ( 11) is the positive sum of convex functions so it is also convex.The second line of terms in equation ( 11) can be expressed as multiplied by a positive number.Hence, it is also convex.Lastly, by the positivity of the integrand, for any N , m 2 ≥ 0 and finite, so the last line is also convex.Thus, there exists a non-empty set U of coupling constants such that the action V is convex.With this above observation, one can apply theorem 1.1 from [36] to the model defined by V .Since the moments and SDE's of the models (11) and ( 5) are the same in the limit, the result applies to the latter model, giving us the following theorem.Theorem 3.2.There exists an ϵ > 0 such that, for t 2 , t 4 ∈ U ∩ B ϵ (0) and any word W ∈ C⟨A, B⟩, converges to the unique solution to the SDE's of the effective ensemble.
This theorem implies that in some small ball of the coupling constants near zero, the convergent and formal models coincide.

The perturbative expansion
As mentioned in the introduction, a formal matrix model is a well-defined formal generating function of Gaussian matrix integrals that is constructed by expanding all non-Gaussian terms of the potential and then interchanging the order of integration and summation.The resulting Feynman diagrams of such an integral are maps (or their dual fat graphs) [11].This follows from the fact that Gaussian matrix integrals can be expanded in terms of maps, which can then be organized by the genus of the associated surface.In this work, we are interested in the types of maps that come from 2-matrix integrals with bi-tracial interactions, which will be introduced in the following sections.

A primer on maps
We will begin by introducing some general terminology on maps.A map of genus g is a 2-cell embedding of a graph into an oriented surface of genus g up to orientation-preserving homeomorphisms of the surface.In this work we are focused on maps with connected graphs of genus zero, which we will refer to as planar maps.
Maps can be constructed by gluing the edges of polygons in an orientation-preserving manner, i.e. no twists.The unglued edges of polygons are referred to as half-edges.A rooted map is a map with a distinguished rooted edge.Rooted maps appear when computing moments and cumulants while unrooted maps appear when computing the partition function.In particular, cumulants and the logarithm of the partition function count connected maps.Note that maps have an associated topological invariant known as their genus, which can be computed using Euler's formula.
For our model we are interested in 2-colored unstable planar maps.A 2-colored map is a map whose half-edges have one of two assigned colors.Such colors have to match that of the other half-edge they are glued to in order to form such a map.An unstable map is a map that is glued from 2-cells whose topology corresponds to unstable Riemann surfaces with boundaries i.e. a disc or cylinder.Note that an ordinary rooted connected map glued from only polygons is planar if and only if its Euler characteristic is two.However, in unstable maps the notation of graph connectedness and map connectedness no longer coincide, so Euler's formula for genus is not always directly applicable.Unstable maps have a decomposition into graph connected components, by treating each 2-cell with the topology of a cylinder as two disconnected 2-cells with the topology of a disc.The removed part we will refer to as a branch.Thus we can associate every unstable map with a graph, where each edge is a branch and each vertex is a graph component.One sufficient condition for an unstable map to be planar is if each graph component is planar and the associated graph described above is a tree.See [45] for more details.
The enumeration of colored maps has long been of interest in the study of formal matrix integrals, but work on unstable maps has more recently been approached in [45,38] as well as within the more general notion of stuffed maps [9,10].As far as the authors are aware, the enumeration of maps with both qualities has not appeared in any works before.

From matrix models to map enumeration
Let us consider the model with the action (5) formally but we will add a redundant parameter t initially, which will keep track of the number of vertices: The propagators for the Gaussians are where the entries of A and B are independent at the level of Gaussian integrals in the formal sum.Via standard techniques [45], the model has a has a genus expansion, i.e. the moments can be written as where The set UM g (v) is the set of all of genus g 2-colored unstable maps with v vertices glued from a rooted polygon whose coloring corresponds to the word W and the following set of 2-cells: See Figure 4 for a visualization of these 2-cells.
The realization of the correspondence of colored polygons to cyclic words can be described as follows.The trace of a word W of length ℓ in the alphabet formed from A and B has a corresponding cyclic sequence of colors of length ℓ.This cyclic sequence of coloring is then mapped to the colors of edges of an ℓ-gon.See for example Figure 3.In equation ( 13), the Feynman weight of a map Σ ∈ UM g W (v) is given by where n i (Σ), for 1 ≤ i ≤ 7, is the number of 2-cells corresponding to the numbers above used in the gluing of the map Σ.The coefficients in front of these weights come from a rescaling needed to construct the factor |Aut| in equation (13).Usually, potentials are nicely normalized so that the Feynman weight is precisely the product of coupling constants.However, because there is one coupling constant in front of many terms in the effective action, this is not possible with our model.Because of this, information that helps distinguish components is lost in the final expressions, which we will find actually simplifies matters.

The second moment
In this section we will derive our formula for the second moment, but we must first look at another moment.Consider Our goal is to show that this formal series is precisely zero for our model.Note that the set UM 0 ABAB (v) is not empty for all v, for an example see Figure 5. Rather, we will show that the formal series is zero by showing that all positive contributions cancel with contributions from the negative sign in the Feynman weight corresponding to chequered unrooted quadrangles.To do so, we must first study the set UM 0 ABAB (v).When v = 1 or 2,the set is empty since there is no planar gluing of a rooted quadrangle with the coloring corresponding to ABAB.For v > 2, the sets are not necessarily empty, but the following fact may be observed for all v. Proof.If v = 1, 2 the claim obviously holds, so let v > 2. Consider some map Σ ∈ U M 0 ABAB (v) with no chequered colored quadrangles or opposite colored cylinders.We will show such a map cannot exist.
Without loss of generality, consider one of the red half-edges of the rooted face.It must be paired to another red half-edge.There are three options: a half-edge of an unrooted red quadrangle, an adjacent colored quadrangle, or a red cylinder.Note that, since the number of vertices is strictly greater than one, the edge cannot be paired with the other red half-edge of the rooted face.In all cases, this new 2-cell must connect to another distinct red coloured 2-cell, since after considering the initial half-edge, there are an odd number of half-edges remaining.Since there are finitely many 2-cells with an even number of red half-edges used in a gluing that all need to be paired with half-edges of the same color, it must eventually connect to the other red half-edge of the rooted 2-cell by the pigeonhole principle.
This above argument holds for the blue half-edges of the rooted fat vertex as well.Thus, Σ must have at least two closed loops of colored edges that can be traced to and from the rooted face.No such map can be embedded into a sphere since this would result in these two different colored loops crossing, which is impossible without an chequered colored quadrangle.Lemma 4.2.For any rooted unstable colored map Σ 1 containing chequered colored quadrangle there is another such map with Feynman weight W (Σ 1 ).Similarly, for any rooted unstable colored map Σ 2 containing an opposite colored cylinder, there exists another such map with Feynman weight Proof.For v > 2, consider a rooted unstable colored map Σ 1 containing an chequered colored quadrangle.We can construct a new map Σ ′ 1 containing a an opposite colored cylinder as follows: 1. Treat one of the non-rooted chequered colored quadrangle faces in Σ 1 as a boundary.
2. Glue an opposite colored cylinder to this boundary as in Figure 6.
We claim this procedure provides us with a planar map Σ ′ 1 .From our discussion in Section 4.1, in order to show the gluing in Figure 6 is planar, it suffices to show that the left graph component of Σ ′ 1 is planar, since the right component is clearly planar and the only branch in this case forms no handles.The resulting graph component will have three fewer vertices, one fewer edge and two more faces than Σ 1 .Thus, it has the same Euler characteristic as Σ, so also the same genus.The resulting map has the same faces except with one opposite colored cylinder instead of one chequered colored quadrangle.Recall from equation ( 14) that these two 2-cells contribute the same factor up to a sign in the Feynman weight.This completes the first claim.Next, consider a map Σ 2 that contains an opposite colored cylinder.Since the map is planar the branch in this cylinder must connect two distinct graph components.We then apply the following procedure: 1.If one of the 2-gons is glued to itself, do the reverse procedure of above.(b) Glue a chequered colored quadrangle to each boundary as in Figure 7.
The resulting map Σ ′ 2 is clearly planar if the first case holds.In the second case, Σ ′ 2 will have two more edges, one less vertex, and three more faces than either planar graph component connected by the branch in Σ ′ .Thus the resulting map will also be planar.In both cases, the Feynman weight Figure 7: A planar gluing of a non-rooted adjacent colored quadrangle to two different coloured 2-gons.
Theorem 4.3.For the formal model with the effective action, m 0 1,1,1,1 is exactly zero.Proof.In the effective action, the redundant parameter t that counts the number of vertices is set to one.For a discussion on why this parameter is redundant see Chapter 1.2.3 of [26].We also know that the set UM 0 ABAB (v) is empty for v = 1, 2. It is clear then that our moment is of the form We known from Lemma 4.3 that each map must contain at least one chequered colored quadrangle or an opposite colored cylinder.We also know from Lemma 4.2 that for each map with a chequered colored quadrangle there is a map with the same gluing configuration, except that the chequered quadrangle is replaced with an opposite colored cylinder and vice versa.Additionally, note that any map with one root has a trivial automorphism group.The result is that when we collect terms of the same power, the number of terms with a positive sign will always equal the number of terms with a negative sign.
With these results we may now succinctly prove the main result of this work.
Proof of Theorem 2.1.Algebraically solving the loop equations in Appendix C in terms of t 2 , t 4 , and m 0 2 gives the formula Using the fact that m 0 1,1,1,1 = 0, we can rearrange for m 2 in terms of t 2 and t 4 .There are two roots, but we must choose the one that is not always negative for t 4 > 0, in order for the second moment to be positive.
Combining this with Theorem 3.2 gives the main result.
Based on computations done in [37], we conjecture that from the second moment all other moments can be computed recursively.The proof of this conjecture at the moment seems to be a challenging combinatorial problem.
5 The free energy

Derivation
For this particular model we can use our knowledge of the second moment to compute the leading order term of the logarithm of the partition function in the large N limit, commonly referred to as the free energy [15,14].
We know from [45] that the free energy of our models has the genus expansion where The set UM g (v) is the set of all maps of genus g with v vertices glued from the list of 2-cells in Figure 4.The proof of Theorem 2.2 follows from a simple computation from the formula for the fourth moment which can be found in Appendix B.
Proof of Theorem 2.2.It is clear that Since this is a formal series in N −2 , we may swap the order of the limit and differentiation, Integrating both sides and using the formula for the Gaussian Dirac ensemble in the large N limit from Appendix A, we arrive at

Random maps and criticality
The free energy can be used to find critical behavior of the model, from which an asymptotic expansion can be computed.Such expansions have been shown to bridge connections to theories of 2D quantum gravity [22].In [38], several Dirac ensembles were shown to have the same critical exponents and asymptotic partition functions as various minimal models.We would also like to emphasize that this critical behavior does not correspond to a spectral phase transition which is of interest for Dirac ensembles [6,44], but rather the type of critical behavior mentioned that connects matrix models to random commutative geometries.In some sense this can be thought of as a continuum limit.For a formal matrix integral with a genus expansion, its weighted map generating functions have an interpretation as a discrete probability distribution.For simplicity set t 2 = 1.Consider, for example, planar maps.If a non-trivial configuration of coupling constants are such that F 0 is finite and greater than zero, then we say such a configuration is admissible.For admissible configurations we are then able to define the discrete probability distribution over UM 0 (v), We know that there exist admissible configurations from Theorem 2.2.Usually, in matrix models, each coupling constant corresponds to a different trace term in the potential.In these cases we can compute the expectation number of 2-cells of a certain topology by differentiating the free energy.This is not the case in our potential, so a map theoretic interpretation of the second derivative is more complicated.However, it is still a quantity of interest and can roughly be thought of as a weighted expected number of 2-cells from Figure 4, 18 With this interpretation we can see that the expected number of 2-cells diverges along the critical curve Note that the solutions of this equation are only when t 4 is less than zero.Thus, this critical behavior is only seen in the formal model and not the convergent solution.
There exist formal notions of convergence of these probability distributions on maps to random metric spaces.The critical exponent of interest here is the first non-zero power, which is usually of the form 1 − γ.This is known as the string susceptibility exponent, and often indicates to what random metric spaces the above one will converge to in the Gromov-Hausdorff topology.For more details we refer the reader to [12].
For simplicity set t 2 = 1.We may asymptotically expand the partition function around the critical point t c = −1/8, This allows us to deduce that γ = 1/2, which is associated with the limiting metric space known as the continuum random tree.This is also an exponent that does not appear often in random matrix models, but is common in tensor models [49].This may suggest that the maps enumerated here have a realization as the triangulations seen in tensor models.For comparison, the quartic type (1, 0) Dirac ensemble studied in [38] has a string susceptibility exponent of −1/2, which is associated with the Brownian map.However, future work is still needed to establish such a convergence.

Conclusions and Outlook
In this work we computed the second moment and the free energy of the quartic type (2, 0), (1, 1), and (0, 2) Dirac ensembles in the large N limit.This was done by studying properties of unstable colored maps and the associated Schwinger-Dyson equations (SDE's).Applying the results of [36], we were then able to show that the solution for all moments for both the convergent and formal models up to the leading order is unique.Furthermore, we explicitly computed the first twenty moments of these models.These results can be compared to past numerical work.For example, the plot of the second moment in Figure 7a and 7b of [32] bears a strong resemblance to the large N solution presented here.This seems to indicate that convergence is rather fast since the matrix size in these simulations was rather small, between five and ten.Additionally, in Figure 8, the solution derived here can be seen to perfectly fit within the bootstrapped bounds computed in [37], as expected.3 and 4 from [37] with the second moment solution from Theorem 2.1 overlaid in red.In [37], the models are the same except t 4 is set to one.Each color corresponds to a different region of possible solutions to the SDE's generated by considering various positivity constraints on the moments.
One limitation of what we currently understand is that despite being able to compute many moments, the generating functions of both moments and Dirac moments is an enigma.It is not clear which choice of generating functions the SDE's can be written succinctly in terms of.The types of generating functions used in studying past multi-matrix models such as the two-matrix model [27] do not seem to be enough to close these equations.A natural candidate might be a generating function in terms of Dirac moments, but such a formulation is not known to the authors at this time.
We hope that this work will lead to a general formula for all Dirac moments, or equivalently its generating functions; as well as new techniques to study previously unsolvable multi-matrix models.In particular, ideas presented here may be useful in studying Dirac ensembles on gauge-matrix spectral triples [55].It would be also interesting to see if the recent work in [53] can be generalized to multi-trace matrix models, potentially leading to solutions to all orders of even more challenging Dirac ensembles.
In the large N limit the answer will be the same as that of the following integral since they have the same loop equations:

B The limiting moments
By the parity of the integral and the A|B symmetry the normalized moments simplify as follows in the large N limit.
We list here list the first twenty unique moments by considering all words of each length up to eight. .
The first few Dirac moments can be written as:

C Examples of Schwinger-Dyson equations
Recall that we denote limiting moments using the notation m ℓ1,ℓ2,...,ℓq = lim Note that this notation is well-defined, since the model is symmetric in A and B.

Figure 1 :
Figure 1: An example of a 2-cell with the topology of a cylinder where one boundary is a 5-gon and the other is a quadrangle

Figure 2 :
Figure 2: An example of a 2-colored unstable map of genus one.

1. A red quadrangle 2 . A blue quadrangle 3 .
A quadrangle with two adjacent red edges, and two adjacent blue edges 4. A quadrangle with two red edges whose neighbours are blue edges 5.A red 2-cell with the topology of a cylinder and boundaries of length two 6.A blue 2-cell with the topology of a cylinder and boundaries of length two 7.A 2-cell with the topology of a cylinder and boundaries of length two, one red and one blue.

Figure 3 :
Figure 3: Moments corresponding to cyclic words W in A and B correspond to rooted polygons with alternating colored half-edges.For example the word AAABAB corresponds to the above hexagon, before choosing a root.
(a) A blue quadrangle.(b) A red quadrangle.(c) An adjacent colored quadrangle.(d) A chequered colored quadrangle.(e) A red cylinder.(f) A blue cylinder (g) An opposite colored cylinder.

Figure 4 :
Figure 4: All types of 2-cells that are used in gluings of 2-colored unstable maps enumerated by the model.13

Lemma 4 . 1 .
Any map in UM 0 ABAB (v) must contain at least one chequered colored quadrangle or an opposite colored cylinder.

Figure 6 :
Figure 6: A planar gluing of a non-rooted adjacent colored quadrangle to an opposite colored cylinder.
Treat both 2-gons on each graph component of Σ 2 as a boundary.

Figure 8 :
Figure 8: A reproduction of Figures3 and 4from[37] with the second moment solution from Theorem 2.1 overlaid in red.In[37], the models are the same except t 4 is set to one.Each color corresponds to a different region of possible solutions to the SDE's generated by considering various positivity constraints on the moments.