#  Dually weighted multi-matrix models as a path to causal gravity-matter
systems

Juan L. A. Abranches<EMAIL_ADDRESS>Okinawa Institute of Science and
Technology Graduate University, 1919-1, Tancha, Onna, Kunigami District,
Okinawa 904-0495, Japan  Antonio D. Pereira<EMAIL_ADDRESS>Instituto de
Física, Universidade Federal Fluminense, Campus da Praia Vermelha, Av.
Litorânea s/n, 24210-346, Niterói, RJ, Brazil  Institute for Mathematics,
Astrophysics and Particle Physics (IMAPP), Radboud University, Heyendaalseweg
135, 6525 AJ Nijmegen, The Netherlands  Reiko Toriumi<EMAIL_ADDRESS>Okinawa Institute of Science and Technology Graduate University, 1919-1,
Tancha, Onna, Kunigami District, Okinawa 904-0495, Japan

###### Abstract

We introduce a dually-weighted multi-matrix model that for a suitable choice
of weights reproduce two-dimensional Causal Dynamical Triangulations (CDT)
coupled to the Ising model. When Ising degrees of freedom are removed, this
model corresponds to the CDT-matrix model introduced by Benedetti and Henson
[Phys. Lett. B 678, 222 (2009)]. We present exact as well as approximate
results for the Gaussian averages of characters of a Hermitian matrix $A$ and
$A^{2}$ for a given representation and establish the present limitations that
prevent us to solve the model analytically. This sets the stage for the
formulation of more sophisticated matter models coupled to two-dimensional CDT
as dually weighted multi-matrix models providing a complementary view to the
standard simplicial formulation of CDT-matter models.

##  I Introduction

The search for a consistent theory of quantum gravity that is valid up to
arbitrarily short length scales remains an open challenge in theoretical
physics. The direct quantization of the gravitational field described by
general relativity using the standard perturbative field-theoretic techniques
renders a perturbatively non-renormalizable quantum field theory (QFT), see,
e.g., [1, 2, 3]. In the history of the construction of quantum theories of the
fundamental interactions, perturbative non-renormalizability of a given well-
grounded theoretical model was circumvented by the replacement of such a model
by a more fundamental description either by the addition of new fields or by
the discovery that the degrees of freedom to be quantized were different ones.
As such, successful models were constructed and the fulfillment of
perturbative renormalizability became a paradigmatic prescription. Such an
attitude could be taken in the case of quantum gravity. Perhaps, the Einstein-
Hilbert action is not the correct starting point and just corresponds to an
effective description of classical gravity at sufficiently low energies and
the fundamental degrees of freedom are completely different. Alternatively,
one could add more fields that restore a well-behaved ultraviolet behavior for
the would-be quantum theory of gravity. In all those cases, there is an
insistence in keeping the tools of perturbative renormalization as the guiding
principle in the construction of the underlying QFT. However, as it is well-
known nowadays, perturbative renormalizability is neither necessary nor
sufficient to define a theory that is valid across arbitrary length scales.
The existence of Landau poles in scalar field theories in four dimensions or
in quantum electrodynamics shows that perturbatively renormalizable theories
might require the introduction of a ultraviolet cutoff at some finite energy
scale in order to be non-trivial, see, e.g., [4, 5, 6, 7]. In practice,
however, such Landau poles live at scales far beyond our experimental
capabilities and therefore are harmless. Similarly, by introducing a cutoff
for the standard QFT of general relativity, one has a perfectly well-defined
framework to compute quantum corrections to gravitational processes at energy
scales below the cutoff [8, 9, 10]. Hence, constructing a fundamental theory
of quantum gravity, as opposed to such an effective field theory does not
necessarily require the construction of a perturbatively renormalizable QFT.
One possibility to define a QFT that is valid up to arbitrarily high energies
is by demanding that the coupling constants that enter physical observables
reach a renormalization group fixed point. When this happens, the theory
attains a scale-invariant regime and the ultraviolet cutoff can be safely
removed. A paradigmatic example is quantum chromodynamics which is
asymptotically free, i.e., the theory reaches the free (or Gaussian) fixed
point in the ultraviolet111In this case, the theory is also perturbatively
renormalizable., see [11, 12]. In this case, since the coupling becomes
sufficiently small at high energies, perturbation theory is applicable and can
probe the existence of such a fixed point. Alternatively, the theory could
have an interacting (or non-Gaussian) fixed point in the ultraviolet. In this
situation, the values of the couplings at the fixed point might not be small
and perturbation theory is not sufficient to probe it. When such a non-trivial
fixed point exists222Actually, besides having a non-trivial fixed point, it is
necessary that just finitely many couplings are adjusted in order to reach it,
otherwise the theory is not predictive., we say that the QFT is asymptotically
safe. In [13], Weinberg put forward the conjecture that quantum gravity could
be realized as an asymptotically safe QFT. Evidence for that was obtained
through the so-called 2 + $\epsilon$ expansion [14, 15, 16, 17]. However, due
to its non-perturbative nature, the search for such a fixed point needs
different techniques with different systematics in order to make its finding a
robust claim. For many years, this conjecture was not investigated
systematically due to the lack of suitable technical tools. However, during
the 1990’s, this situation suffered a twist due to the development of two
major frameworks: the use of functional or non-perturbative renormalization
group equations to quantum gravity, see, e.g., [18] and reviews [19, 20, 21,
22, 23, 24] and references therein and the proper understanding of how to put
quantum gravity on a lattice in a consistent fashion with causality
constraints, see [25, 26] and reviews [27, 28] and references therein. Those
computational techniques can be viewed as complementary approaches to the
search for a non-perturbative ultraviolet fixed point for quantum gravity and
it has become usual to call the search for the fixed point using continuum
techniques as the asymptotic safety approach while the search for a suitable
continuum limit of discretized path integrals with causal constraints as the
Causal Dynamical Triangulations (CDT) program. Both lines of research have
provided a large body of evidence for a suitable continuum limit in four
dimensions. We refer the reader to the references in the reviews above
mentioned for a comprehensive list of the most recent results in those
fields333We highlight that both approaches follow the inspiring idea put
forward by Weinberg but there is no reason a priori to expect that the would-
be continuum limits belong to the same universality class. In particular, CDT
implements a Lorentzian path integral while the continuum computations
employing the functional renormalization group are mostly Euclidean, with few
exceptions, see, e.g., [29, 30, 31, 32, 33]..

Despite the tremendous progress achieved in the physically motivated four-
dimensional case, the employment of non-perturbative techniques require
substantial truncations and approximations that still need a lot of effort to
produce quantitatively reliable results for those universal quantities, i.e.,
observables that could tell the existence of the continuum limit. Yet a
fruitful playground is to consider two-dimensional quantum gravity. In this
case, the interplay between continuum and discrete descriptions of the path
integral of quantum gravity has witnessed a significant progress over the last
decades ranging from the development of novel techniques to perform explicit
calculations to rigorous mathematical results. We refer to, e.g., [34, 35] for
reviews on the topic. In particular, the Dynamical Triangulation program was
born in two dimensions within an Euclidean setting while its continuum
counterpart is encoded in Liouville quantum gravity. In particular, a
successful implementation of the discretization of the path integral of
quantum gravity and thus a sum over geometries and possibly topologies can be
achieved by purely combinatorial means with the use of the so-called matrix
models [34, 36, 37, 38, 39, 40, 41, 42, 43]. Triangles can be taken as dual
representations of cubic-matrix vertices and Feynman diagrams of such a model
are simply triangulations of two-dimensional surfaces. The perturbative
expansion of matrix-models partition functions can be organized in powers of
$1/N$ with $N$ standing for the size of the Hermitian random matrices and each
power of such an expansion is associated with a specific genus $g$. Such a
remarkable expansion allows for the investigation of continuum limits, e.g.,
where just spherical topologies contribute (the so-called planar limit) or
where all topologies are taken into account (double-scaling limit). Thanks to
the rich combinatorial framework behind the theory of random matrices, very
powerful results could be established in such a discrete-to-continuum
approach. The two-dimensional case is simple enough so that one can still
solve the Dynamical Triangulation model with no need to make use of random
matrices, but such a perspective opens up the possibility to think about
higher-dimensional discrete approaches as theories of higher-order tensors. So
far, the higher dimensional tensor models were not successful in producing a
suitable continuum limit in higher dimensions. Primarily, the main obstacle
was the lack of the analogue of the $1/N$-expansion in matrix models. Such a
difficulty was lifted thanks to the works of Gurau that introduced the so-
called colored tensor models [44, 45, 46, 47]. Yet the continuum limits
obtained from the colored tensor models up to date do not feature an extended-
geometry-like phase very much in agreement with the simulations performed for
Dynamical Triangulations in dimensionality greater than two [48, 49, 50, 51,
52], see however [53, 54, 55, 56, 57, 58, 59, 60, 61]. This suggests that one
could inspect more intricate continuum limits from those models, see, e.g.,
the realization of such a reasoning following the functional renormalization
group [62, 63, 64, 65], or that a restriction to the configuration space
should be implemented. This is precisely what is achieved with CDT which
displays a robust body of evidence that such a suitable continuum limit exists
in higher dimensions [66, 25, 67, 68, 69, 70, 71, 72, 73, 74, 75].
Nevertheless, the analogue of CDT in a tensor-model language is still unknown.
In this paper, we explore a matrix model that would be such a realization in
two dimensions introduced by Benedetti and Henson in [76]. In essence, this
model implements the (non-local) global time foliation of CDT together with
the avoidance of spatial topology change by means of local constraints. This
is achieved by the so-called dually weighted matrix models introduced in [77,
78, 79, 80]. On top of the CDT configurations generated by the Benedetti-
Henson model, we will introduce degreess of freedom arising from the Ising
model, i.e., we develop a (multi-)matrix model of the Ising model coupled to
CDT in two dimensions. Alternatively, the combinatorial models can be enriched
with more structure such as group-theoretic data giving birth to the so-called
Group Field Theories which may or may not include a Lorentzian setting [81,
82, 83, 84]. Thus, two-dimensional quantum gravity might be too simplistic for
some purposes but certainly is an inspirational source for more sophisticated
candidates in higher dimensions.

Moreover, a realistic description of our Universe must accommodate matter
degrees of freedom whose fluctuations can affect the existence of a suitable
continuum limit. Hence, it is conceivable that our comprehension of the
potential mechanism that drives quantum gravity asymptotically safe
(irrespective of the approach followed here) depends on a simultaneous
treatment of gravitational as well as matter degrees of freedom. This
ambitious goal has witnessed significant progress over the last two decades
both from the point of view of the functional renormalization group as well as
from lattice simulations and encouraging results were found [85, 23, 86, 28,
87]. In the present work we provide several results regarding a simple but
still quite rich gravity-matter model: two-dimensional CDT coupled to the
Ising model. There are several motivations to look at this model but we
emphasize the following: Implementing the causality constraint at the level of
a matrix model is a first step towards its implementation in tensor models. As
it is known, such a model does not have a known exact solution and developing
appropriate tools to deal with it is mandatory. The inclusion of the Ising
model produces a multi-matrix model which is interesting on its own. Clearly,
having the Ising model as the matter component is tremendously far from the
rich structure of the Standard Model of Particle Physics coupled to quantum
gravity. Yet the purpose here is to explore the impact of the dynamical (and
causal) lattice to the Ising model and vice-versa. We emphasize that matrix
models for Euclidean two-dimensional gravity coupled to the Ising model were
investigated in the past [88, 89, 90, 91] as well as its higher-dimensional
version [92, 93, 94]. Moreover, the coupling between the Ising model and CDT
was also investigated since the birth of two-dimensional CDT, see, e.g., [95,
96, 97, 98, 99, 100, 101]. We establish several properties of such a CDT-Ising
multi-matrix model paving the way for a future complete solution or making its
structure well-grounded for studies using non-perturbative tools such as the
functional renormalization group methods444Functional renormalization
techniques were developed for matrix models, see, e.g., [102, 103, 104, 105,
106]. There are studies involving the functional renormalization group
equation applied to the Benedetti-Henson model, see [107]. or numerical
simulations.

This paper is organized as follows: In Section II, we will briefly review CDT-
like matrix model as a dually weighted matrix model as introduced by Benedetti
and Henson [76]. In Section III, we develop an analysis concerning topologies
which admit global foliations, specifically within the Benedetti-Henson pure
CDT-like matrix model (in Section II) and further on non-orientable symmetric
matrix models. In Section IV, we discuss briefly the well-known two-matrix
models which can be interpreted as Ising model on random two-surfaces. In
Section V, we present our CDT-matrix model coupled to the Ising model by
combining the models in Sections II and IV. In Section VI, we study the
properties of the matrices $C_{m}$ which, in the case of $m=2$, is responsible
for generating the global foliation structure on the dual triangulation graphs
of the CDT-like matrix models. Theorem VI.1 is one of the main results of this
paper. In Section VII, we present the character expansions for partition
functions of the two models, i.e., pure CDT-like matrix model by Benedetti-
Henson and CDT-like matrix model coupled to Ising model presented in Sections
II and V respectively. We base our further analyses on these character
expansions. In Section VIII, we analyze the character expansion of the
partition function for the latter model (i.e., CDT-like matrix model coupled
to Ising model) and present an expression in Proposition VIII.1 in terms of
Clebsch-Gordan coefficients using Weingarten calculus. In Section IX, we
evaluate Gaussian averages of characters for a given representation $r$,
$\langle\chi_{r}(A)\rangle_{0}$ and $\langle\chi_{r}(A^{2})\rangle_{0}$ for
CDT-like matrix models presented in Section II, using various techniques,
e.g., Schur-Weyl duality, theory of symmetric group algebra, etc. Theorem
IX.6, one of our main results, is related to Brauer algebra. In Section X, we
conclude by revisiting and summarizing the various results reported in this
work.

##  II Causal Dynamical Triangulations as a matrix model: a short review

Let us describe via matrix model the Causal Dynamical Triangulation (CDT) in
two dimensions. The CDT was first worked and solved in [108] as a lattice
model. In [76], a definition through a matrix model was given, and we follow
this definition. Let us also refer to this model as pure555We specifically
refer pure (gravity) because later, we will talk about a matrix model which
generate random two-surfaces (gravity) with matter.  CDT-like matrix model of
Benedetti-Henson. The CDT-like matrix model graphs can be represented as
ribbon graphs just like other matrix models. In fact, these CDT-like matrix
models fall into a class of so called dually weighted models. In [78],
Kazakov, Staudacher, and Wynter introduced a so called dually weighted graphs
where different weights are assigned for vertices with different coordination
numbers and for faces with different lengths. As one will see below in
elaborating the features of such CDT-like matrix models, our models generate
the dually weighted graphs only with vertices with coordination number 3, and
with faces with arbitrary lengths however, of particular type (restricting to
only having two edges of a certain type, which we will name timelike).

######  Definition II.1.

A prime ribbon graph is a set of topological discs and topological rectangles
satisfying a the following properties: no two rectangles or two discs
intersect; each rectangle has exactly one pair of opposite sides that are
contained in disc’s circumferences, and this is the only intersection between
rectangles and discs.

######  Definition II.2.

A rectangle’s edge that is not contained in a disc is called a strand.

######  Definition II.3.

A ribbon graph is the boundary of the union of the discs and rectangles of a
prime ribbon graph.

Additionally, one can represent such ribbon graphs as multigraphs. Given a
ribbon graph, we assign a vertex to each disc and an edge to each rectangle.
If a disc and a rectangle are connected, i.e., they intersect, then respective
vertex and edge are connected. See Fig. 1 for an example. We may also call a
ribbon graph’s discs and rectangles as vertices and edges, respectively. We
define a face as a closed curve formed by strands, and we carry the notion of
faces from the ribbon graph to the multigraph. An example of a face of a
ribbon graph is shown in Fig. 1(a).

(a) A closed loop is called a face of a ribbon graph. The timelike edges are
shown in blue, and the spacelike edges are shown in red.

(b) The multigraph representation of the ribbon graph in Fig. 1(a) is shown in
solid lines and black vertices. Dashed lines represent the edges of the dual
graph and the gray dots are vertices of the dual graph.

Figure 1: A face of a CDT ribbon graph.

The CDT-like matrix model ribbon graphs are edge-colored (not proper),
partitioned into two sets: spacelike edges and timelike edges. Together with
the defining properties of a ribbon graph and the edge coloring, the CDT
ribbon graph also has the following defining properties (CDT conditions):

  1. 1.

Every vertex is three valent and has exactly one timelike edge and two
spacelike edges incident to it.

  2. 2.

Every face has either two timelike edges or none.

(a) Vertex of a CDT ribbon graph.

(b) Edges of a CDT ribbon graph.

Figure 2: Components of the CDT ribbon graph.

Notice that the faces can have any number of spacelike edges. Fig. 2 shows a
representation of the vertices and edges. A triangulation is obtained with the
dual graph; the dual graph is defined by associating a vertex with each face
of the graph, and if two faces share an edge, their respective vertices are
connected by an edge (see Fig. 2(b)). This also leads to each vertex in the
graph being associated with a face in the dual graph, and since all vertices
are three valent, all faces in the dual graph are triangles. In Fig. 3, we
show an example of a sphere triangulation that satisfies the CDT conditions
given above.

Figure 3: An example of a CDT triangulation graph (dual of a corresponding CDT
ribbon graph). spacelike edges in the ribbon graph are dual to timelike edges
in the dual triangulation graph, and vice-versa. The vertices of the dual
triangulation graph can have two or zero spacelike edges (in red) and any
number of timelike edges (in blue).

A matrix model that generates ribbon graphs, satisfying these defining CDT
conditions above can be constructed by associating a Hermitian matrix $A$ with
spacelike edges of a ribbon graph and a Hermitian matrix $B$ to timelike edges
of a ribbon graph. An auxiliary partition function for such a pure CDT-like
matrix model is given by

$\mathcal{Z}=\int
dAdB\;e^{-N\mathrm{Tr}\left[\frac{1}{2}A^{2}+\frac{1}{2}({C_{2}}^{-1}B)^{2}-gA^{2}B\right]}\,.$
(1)

The action (the negative of the exponent in (1)) of the model has the term
$-NgA^{2}B$, responsible for the property that every vertex have exactly two
spacelike edges and one timelike edge. A matrix $C_{2}$ is introduced in the
quadratic term of $B$, and it satisfies, for a positive integer $p$ and in the
large $N$ limit,

$\mathrm{Tr}[(C_{2})^{p}]=N\delta_{p,2}\,.$ (2)

This condition on $C_{2}$ is what sets the property that faces have only two
or zero timelike edges. In Section VI, we study this matrix $C_{2}$ in more
detail. Notice that this constraint (2) on $C_{2}$ is overdetermined, once
$p>N$. Therefore, the condition (2) on $C_{2}$ can only be imposed in the
large $N$ limit. The Gaussian average of a function $f(A,B)$ is defined by

$\langle f(A,B)\rangle_{0}=\frac{\int
dAdB\;f(A,B)\;e^{-N\mathrm{Tr}\left[\frac{1}{2}A^{2}+\frac{1}{2}({C_{2}}^{-1}B)^{2}\right]}}{\int
dAdB\;e^{-N\mathrm{Tr}\left[\frac{1}{2}A^{2}+\frac{1}{2}({C_{2}}^{-1}B)^{2}\right]}}\;,$
(3)

and we find the following properties:

$\langle A_{ij}A_{kl}\rangle_{0}=\frac{1}{N}\delta_{ik}\delta_{jl}\,,$ (4)
$\langle B_{ij}B_{kl}\rangle_{0}=\frac{1}{N}{C_{2}}_{ik}{C_{2}}_{jl}\,,$ (5)
$\langle A_{ij}B_{kl}\rangle_{0}=0\;.$ (6)

Since the condition on $C_{2}$ in (2) can only be set at large $N$, the
partition function for CDT is only be obtained in the large $N$ limit,

$\mathcal{Z}_{CDT}=\lim_{N\rightarrow\infty}\mathcal{Z}\;.$ (7)

##  III Properties of the Benedetti-Henson model: allowed topologies and
(non-)orientability

The two defining properties, albeit being local, put a global constraint on
the topologies that our ribbon graphs (and in their dual triangulation) may
represent. Studying the properties of the elements of the graph enables us to
determine the values of Euler characteristic which are allowed. The analysis
can be done either with the ribbon graphs or with the triangulations dual to
them. Here, let us work on the ribbon graphs.

Let us study some properties satisfied by the faces. A face that has no
timelike edges is simply a set of connected spacelike edges. See Fig. 4(b).
For a face that has two timelike edges, since the vertices can only have one
timelike edge, a timelike edge in a face is always a neighbor of two spacelike
edges of the face. See Fig. 4(a). This way, the edges of a given face of the
latter type form a sequence composed by a timelike edge, a several spacelike
edges, a timelike edge, a several spacelike edges, and then connecting back to
the first timelike edge. Thus, this type of face has two disjoint sets of
connected spacelike edges, and these two sets are separated by timelike edges.
Let us call these the two sets as spacelike boundaries of a face.

######  Definition III.1.

A strip is defined as a set that contains faces which are sequentially
connected by timelike edges, and also contains the edges and the vertices of
these faces. Each face in a strip is glued by timelike edges to either two
faces in the same strip (it can also glue to itself) or none, and one timelike
edge is shared exactly by two faces (these two faces are possibly the same
face). If a spacelike edge has multiplicity two in one face or if it belongs
to two faces of the same strip, consider those two appearances as distinct
elements of a strip.

######  Definition III.2.

A boundary of a strip is defined as a set of spacelike edges that are
sequentially connected by vertices.

######  Definition III.3.

A interior of a strip is defined as the set of faces and timelike edges in a
strip.

######  Definition III.4.

A regular strip is a strip that has two boundaries (see Fig. 6(a)).

######  Definition III.5.

A singular strip is defined as a strip composed by a face with no timelike
edges (see Fig. 6(b)).

######  Definition III.6.

A Möbius strip is a non-singular strip that has one boundary (see Fig. 6(c)).

Some important properties satisfied by CDT ribbon graphs are:

  1. 1.

Every face belongs to a strip and is in only one strip.

  2. 2.

A strip has only spacelike edges in its boundary, and all spacelike edges are
in boundaries of strips.

  3. 3.

Since there is a finite number of faces, there is a finite number of strips.

  4. 4.

Every non-singular strip is a periodic sequence which alternates between faces
and timelike edges; and every boundary of a strip is a periodic sequence which
alternates between spacelike edges and vertices.

  5. 5.

A strip has only one or two boundaries.

  6. 6.

A boundary of a strip either bounds one other strip or connects the strip to
itself.

(a) A face with two timelike edges.

(b) A face with no timelike edge.

Figure 4: Components of the CDT ribbon graph. Spacelike edges are shown in
red, while timelike edges are shown in blue.

Properties 1 to 4 are easy to see, but let us take a closer look at property
5. The singular strip is just a face with only spacelike edges, thus as a
strip it has only one boundary, as shown in Fig. 6(b).

(a) The face $F_{1}$ has spacelike boundaries $A_{1}$ and $B_{1}$, and the
face $F_{2}$ has the spacelike boundaries $A_{2}$ and $B_{2}$. Together, they
have the spacelike boundaries $A_{1}\cup A_{2}$ and $B_{1}\cup B_{2}$.

(b) The edge $e_{1}$ is in the faces $f_{1}$ and $g_{1}$. The edge $e_{2}$ is
in the faces $f_{2}$ and $g_{2}$. Both faces $f_{1}$ and $f_{2}$ are in the
strip $s$ and both faces $g_{1}$ and $g_{2}$ are in the strip $t$. In this
example, $f_{1}=f_{2}$.

Figure 5:  Steps of the demonstration of properties 5 and 6 of CDT ribbon
graphs.

For a regular strip or a Möbius strip, we follow the steps below:

  * •

Consider one of its faces, call it $F_{1}$. $F_{1}$ has two spacelike
boundaries, call them $A_{1}$ and $B_{1}$.

  * •

Consider also a face, call it $F_{2}$, that is connected to the previously
considered face $F_{1}$ above by a timelike edge. $F_{2}$ also has two
spacelike boundaries. Call them $A_{2}$ and $B_{2}$. (See Fig. 5(a).)

  * •

These two boundaries are each connected to a boundary of the previous face by
the vertices the faces share. Together, the two faces create two extended
spacelike boundaries, $A_{1}\cup A_{2}$ and $B_{1}\cup B_{2}$.

  * •

This process is inductively repeated to $k$ connected faces in the strip, and
the number of spacelike boundaries is kept two, $A_{1}\cup A_{2}\cup\cdots\cup
A_{k}$, and $B_{1}\cup B_{2}\cup\cdots\cup B_{k}$.

  * •

Call $l$ the number of faces in the strip. Considering property 4 of CDT
ribbon graphs, when the first and the last faces are connected, there are two
possibilities:

    1. 1.

The two extended boundaries close into $A_{1}\cup A_{2}\cup\cdots\cup A_{l}$
and $B_{1}\cup B_{2}\cup\cdots\cup B_{l}$, forming a regular strip, as shown
in Fig. 6(a).

    2. 2.

The extended boundaries connect to each other, $A_{1}\cup A_{2}\cup\cdots\cup
A_{l}\cup B_{1}\cup B_{2}\cup\cdots\cup B_{l}$, creating only one closed
boundary, forming the Möbius strip, as shown in 6(c).

We can obtain the property 6 in a similar manner as we did property 5 above:

  * •

Given a spacelike edge $e_{1}$, consider the two faces (that might be the
same) that share $e_{1}$. Call these faces $f_{1}$ and $g_{1}$.

  * •

Consider the strips (that might be the same) that contain the faces $f_{1}$
and $g_{1}$ respectively. Call them $s$ and $t$ respectively.

  * •

Consider now a spacelike edge, call it $e_{2}$ that has a common vertex with
the previous spacelike edge $e_{1}$. Call $f_{2}$ and $g_{2}$ the two faces
that $e_{2}$ belongs to.

  * •

The faces $f_{2}$ and $g_{2}$ either satisfy $f_{2}\in s$ and $g_{2}\in t$, or
satisfy $f_{2}\in t$ and $g_{2}\in s$. Without the loss of generality, assume
the first case. (See Fig. 5(b).)

  * •

The two neighboring spacelike edges have the same two strips at its sides, $s$
and $t$.

  * •

Inductively, when considering the entire closed cycle of spacelike edges,
either $s\neq t$ or $s=t$. Therefore, there are either two or one strip at its
sides.

(a) A regular strip.

(b) A singular strip.

(c) A Möbius strip.

Figure 6:  Types of strips present in a CDT ribbon graph. spacelike edges are
shown in red, while timelike edges are shown in blue.

Consequently, similarly to the way faces are connected by the timelike edges
to form a strip, the entire graph is construted by strips connected by
spacelike edges. Together with property 6, since the graph is connected, every
strip is connected to each other by a sequence of strips. Therefore, one can
say that a graph is a sequence of strips, and there are three types of strips
defined in III.5, III.4, and III.6.

Possible topologies.
Let us now take a look at the Euler characteristics of the CDT ribbon graphs.
The elements of these ribbon graph can be decomposed into two disjoint sets:
boundaries of strips (defined in III.2) and interiors of strips (defined in
III.3). Setting $V$ as the number of vertices, $F$ the number of faces,
$E_{S}$ the number of spacelike edges, $E_{T}$ the number of timelike edges,
and $E$ the total number of edges, the Euler characteristics can be written
as:

$\chi=V-E+F=V-E_{S}-E_{T}+F=\chi_{B}+\chi_{I},$ (8)

with $\chi_{B}=V-E_{S}$ being the contribution to the Euler characteristics by
the boundaries of strips, and $\chi_{I}=F-E_{T}$ being the contribution by the
interiors of strips. It is always true that $\chi_{B}=0$, since every vertex
has two spacelike edges, and every edge connects two vertices, thus
$2V=2E_{S}$. For $\chi_{I}$, we need to differentiate between different types
of strips. For the regular strip and the Möbius strip, since their faces have
two timelike edges and every edge is shared by two faces, these strips have
the property $2F=2E_{T}$, thus $\chi_{I}=0$. As for a singular strip, since it
is composed of only one face and there are no timelike edges, $\chi_{I}=1$.
Therefore, the Euler characteristic equals the number of singular strips:

$\chi=F_{s}\,,$ (9)

where $F_{s}$ is the number of faces with no timelike edges, therefore it is
equal to the number of singular strips.

From properties 5 and 6 we also get an important fact: The strips can be
ordered by sequencing them according to having a common boundary. This
sequencing induces global foliation, which is one of the fundamental
properties of the definition of CDT [108]. This sequence may be periodic (and
therefore called infinite), just as the strips themselves are a periodic
sequence of faces. The sequence also may be finite, starting with a singular
or Möbius strip, as they can only share a boundary with one other strip, and
may end also with a singular or Möbius strip. The singular strip, then, can
only appear at most twice. Therefore, since the Euler characteristics equals
the number of singular strips, the Euler characteristics can only be 0, 1 and
2.

  * •

Considering only orientable surfaces666One can realize generating orientable
two-surfaces by considering Hermitian matrix models., we know that the Euler
characteristics can be expressed as $\chi=2-2g$, where $g$ is the genus of the
surface. The possibilities follow:

    1. 1.

When the graph has two singular strips, we find that the genus is 0, thus the
surface is a sphere. This graph has the sequence of strips starting with a
singular strip, having some regular strips in the middle, and then ending in
another singular strip.

    2. 2.

When the graph has one singular strip, there is no solution since there is no
orientable surface with $\chi=1$ (corresponding to $g=1/2$).

    3. 3.

Lastly, when the graph has no singular strips, the genus is 1, therefore the
surface is a torus, composed of a periodic sequence of regular strips.

  * •

When we consider nonorientable surfaces777One can realize generating
nonorientable two-surfaces by considering symmetric matrix models., the Euler
characteristics now assumes the form

$\chi=2-(2g+c)\,,$ (10)

where $c$ is the number of cross caps.

    1. 1.

For a graph with two singular strips, (10) can only be satisfied for $g=0$ and
$c=0$. Then, the topology is of a sphere.

    2. 2.

When the graph has only one singular strip, the equation (10) can now be
satisfied with $g=0$ and $c=1$, thus we have the topology of the projective
plane.

    3. 3.

Lastly, when the graph has no singular strips, the equation (10) has two
solutions. One solution is $g=1$ with $c=0$, implying it is a torus. The other
solution is $g=0$ with $c=2$, showing the possibility of the topology of the
Klein bottle.

Existence and construction of topologies.
We have shown above that these topologies are allowed, but it does not yet
mean that they indeed exist. For this, let us define orientation of faces and
strips:

######  Definition III.7.

Orientation of a face: The edges of a face form a closed curve and thus have
the usual notion of orientation. Two faces that share an edge are said to have
a compatible orientation if the directions of the closed curves are opposite
in the common edge. We extend this notion to any two faces by transitivity.

######  Definition III.8.

Orientation of a strip: A regular strip admits a compatible orientation among
all of its faces and we define the orientation of a strip as the orientation
of one of its faces. A Möbius strip does not admit a compatible orientation
among all of its faces. Two strips are said to have a compatible orientation
if a face from one strip has compatible orientation to a face of the other
strip.

In the following, we argue that indeed all the possible topologies discussed
above can exist by a simple construction of putting strips together to form
the sequence of strips.

  1. 1.

Sphere: Start the sequence with a singular strip, follow it by a sequence of
regular strips, and end the sequence with another singular strip. See Fig. 7.

Figure 7: Construction of a sphere graph. The empty circle represents a
singular strip, the circle with two lines inside represents a regular strip,
and the two non-intersecting lines connecting two circles represent a boundary
that keeps the orientation compatible. This example shows four regular strips,
but any non-negative integer is possible.

  2. 2.

Torus: Given a finite sequence of regular strips, glue the first strip to the
last in way that the orientation is compatible. See Fig. 8.

Figure 8: Construction of a torus graph. The circle with two lines represents
a regular strip, and the two non-intersecting lines connecting two circles
represent a boundary that keeps the orientation compatible. This example shows
six regular strips, but any non-negative integer is possible.

  3. 3.

Projective plane: (i) Start the sequence with a singular strip, follow it by a
sequence of regular strips (or none), and end it with a Möbius strip. See Fig.
9(a). (ii) Alternatively, start the sequence with a singular strip, follow it
by a sequence of regular strips, and glue the last boundary to itself as a
cross-cap. See Fig. 9(b).

(a)

(b)

Figure 9: Two different ways of constructing the projective plane. The empty
circle represents a singular strip, the circle with two lines represents a
regular strip, the two non-intersecting lines connecting two circles represent
a boundary that keeps the orientation compatible, and the line going in and
out of the same circle represents a cross-cap. Both examples show four regular
strips, but any non-negative integer is possible.

  4. 4.

Klein bottle: (i) Start the sequence with a Möbius strip, follow it by a
sequence of regular strips (or none) and end it with another Möbius strip. See
Fig. 10(a). (ii) Alternatively, we can change one or both of these Möbius
strips for (a) cross-cap(s). See Figures 10(b) and 10(c). (iii) Furthermore,
we can also have a periodic sequence of regular strips, and we glue the first
strip to the last one in a way that their orientations are not compatible. See
Fig. 10(d).

(a)

(b)

(c)

(d)

Figure 10: Four different ways of constructing the Klein bottle. The circle
with two lines represents a regular strip, the circle with one line represents
a Möbius strip, and the two non-intersecting lines connecting two circles
represent a boundary that keeps the orientation compatible. These examples are
constructed with four regular strips, but any number of them is possible.

##  IV Two-matrix model representing Ising model on random two-surfaces

(a) Spin up vertex

(b) Spin down vertex

(c) Edges

Figure 11: Components of the two-matrix model ribbon graph.

Here, let us introduce two-matrix models which have been widely studied
starting with Itzykson and Zuber [109], followed by many others [110, 111,
112, 113].

Consider the partition function,

$\mathcal{Z}_{\rm IM}=N^{-2}\ln\int dM_{+}dM_{-}e^{-S_{IM}}\,,$ (11)

which is a generating function of connected ribbon graphs. Here $M_{+}$ and
$M_{-}$ are $N\times N$ Hermitian matrices and for a Hermitian matrix $M$, the
measure applied is
$dM=\prod_{i<j}\mathrm{Re}(M_{ij})\mathrm{Im}(M_{ij})\prod_{i}M_{ii}$. For
couplings $\gamma$ and $g$,

$S_{\rm
IM}=N\;\mathrm{Tr}\left[\frac{1}{2}M_{+}^{2}+\frac{1}{2}M_{-}^{2}-\gamma^{-2}M_{+}M_{-}-gM_{+}^{3}-gM_{-}^{3}\right]$
(12)

is said to be the action of the system [114, 89]. One can interpret that the
spins are located on the vertices (of Feynman ribbon graphs generated
perturbatively by this action) represented by the terms in the action
$M_{+}^{3}$ and $M_{-}^{3}$. Fig. 11 shows how these elements appear in the
ribbon graph, and Fig. 12 shows an example of a ribbon graph represented as a
graph.

Figure 12: An example of a part of an two-matrix model graph described by the
action (12). The solid lines represent the ribbon graph. The dashed lines
represent the dual triangulation. The coloring of the vertices represent the
spin.

We can interpret $g$ as a cosmological constant, since it is related to the
volume of the universe. In the ribbon graph, the volume appears as the number
of vertices, while in the dual graph, the triangulation, it is the number of
faces. The coupling $\gamma$ can be related to the Ising temperature $T$ by
$\gamma=e^{-T^{-1}}$. For a Gaussian model with partition function

$Z_{0}=\int dA\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\,,$ (13)

with a Hermitian matrix $A$, consider the average

$\langle f(A)\rangle_{0}=\frac{1}{Z_{0}}{\int
dA\;f(A)\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}}\,,$ (14)

where $0$ denotes the Gaussian average. Wick’s theorem says that the
computation of this function can be decomposed by taking the products of the
propagators related to each edge in our Feynman graphs, where the propagators
are simply the mean values taken for $g=0$,

$\langle{M_{+}}_{ij}{M_{+}}_{kl}\rangle_{0}=\langle{M_{-}}_{ij}{M_{-}}_{kl}\rangle_{0}=\frac{1}{N}\frac{1}{1-\gamma^{-4}}\delta_{il}\delta_{kj}\,,$
(15)

and

$\langle{M_{+}}_{ij}{M_{-}}_{kl}\rangle_{0}=\frac{1}{N}\frac{\gamma^{-2}}{1-\gamma^{-4}}\delta_{il}\delta_{kj}\,.$
(16)

Diagonalizing the kinetic term in (12) takes us to another interesting version
of the model, where now one can interpret that the spins are localized on the
faces of the matrix model graphs (i.e., ribbon graphs) rather than on the
vertices of the matrix model graphs. We achieve this by a few change of
variables. First change to the Hermitian matrices $K$ and $L$ given by

$M_{+}=\frac{1}{\sqrt{2}}(K+L)\;\qquad\mathrm{and}\qquad
M_{-}=\frac{1}{\sqrt{2}}(K-L)\,,$ (17)

so that

$S_{\rm
IM}=N\;\mathrm{Tr}\left[\frac{1}{2}(1-\gamma^{-2})K^{2}+\frac{1}{2}(1+\gamma^{-2})L^{2}-\sqrt{2}gK^{3}-3\sqrt{2}gKL^{2}\right]\,,$
(18)

and with a last substitution to the couplings $\gamma^{\prime}$, $g^{\prime}$,
and Hermitian matrices $U$ and $V$ given by

$\gamma^{\prime-2}=\frac{\gamma-\gamma^{-1}}{\gamma+\gamma^{-1}},\;g^{\prime}=\gamma^{\prime-3/2}\frac{\sqrt{2}g}{(1-\gamma^{-2})^{3/2}},\;U^{2}=\gamma^{\prime}(1-\gamma^{-2})K^{2},\;V^{2}=\gamma^{\prime}(1-\gamma^{-2})L^{2},$
(19)

we obtain

$S_{\rm
IM}=N\;\mathrm{Tr}\left[\frac{1}{2}\gamma^{\prime-1}U^{2}+\frac{1}{2}\gamma^{\prime}V^{2}-g^{\prime}U^{3}-3g^{\prime}UV^{2}\right]\,,$
(20)

which is the Ising model where the spins, therefore $\pm$ signs, are assigned
for each face of Feynman ribbon graphs. Here, $U$ (resp. $V$) can be thought
of as being associated with a ribbon graph edge shared by ribbon graph faces
of the same (resp. opposite) parity. For a third degree vertex, since there is
a face between each neighboring edges, there are three faces (some may be the
same) around it. Either the three faces have the same parity, hence the term
$U^{3}$, or one of the three has a parity different from the other two, hence
the term $3\,U\,V^{2}$. For topological reasons, this interpretation is only
valid for planar graphs. For example, consider the torus graph shown in Fig.
13. This graph can be generated by (20) if we do not restrict ourselves to
planar graphs. Edges associated with $U$ are in green, and edges associated to
$V$ are in orange. In this configuration, the placement of the green edges
around almost all faces would imply that all faces have the same spin, while
the existence of orange edges would imply that some faces have opposite spins,
thus a contradiction.

Figure 13: A torus graph with an edge configuration that does not allow an
assignment of spin to the faces of ribbon graphs.

Note however, that the action (20) after simply the change of variables should
contain exactly the same information as the original action (12). The apparent
restriction discussed above is only associated with the interpretation that we
employed for $U$ and $V$ representing the same or opposite parities
respectively.

##  V Ising model coupled to CDT via matrix model

We present a matrix model that describes Ising model on a two-dimensional
manifold which has global foliation, by combining the ideas of the two models
(CDT-like matrix model and Ising model) we introduced above in Sections II and
IV. We would like to define an action of the matrix model which generates
ribbon graphs satisfying desired properties of both CDT and Ising model.
Hence, we require four types of edges to be considered: spacelike edges
connecting spin-up vertices, spacelike edges connecting spin-down vertices,
timelike edges connecting spin-up vertices, and timelike edges connecting
spin-down vertices. Therefore, we introduce four matrices, $A_{+}$, $A_{-}$,
$B_{+}$, and $B_{-}$ so that they induce notions of the space, time, and
parity as desired. Each of these matrices is represented by a specific type of
half-edges. The following conditions are sufficient in order to implement the
desired properties we described above: To impose the existence of spacelike
edges connecting equal parities, we set

$\langle{A_{+}}_{ij}{A_{+}}_{kl}\rangle_{0}=\langle{A_{-}}_{ij}{A_{-}}_{kl}\rangle_{0}=\frac{1}{N}\frac{1}{1-\gamma^{2}}\delta_{il}\delta_{kj}\,.$
(21)

Similarly, for timelike edges connecting equal parities, we impose

$\langle{B_{+}}_{ij}{B_{+}}_{kl}\rangle_{0}=\langle{B_{-}}_{ij}{B_{-}}_{kl}\rangle_{0}=\frac{1}{N}\frac{1}{1-\gamma^{2}}{C_{2}}_{il}{C_{2}}_{kj}\,.$
(22)

As for spacelike edges connecting unequal parities, we write

$\langle{A_{+}}_{ij}{A_{-}}_{kl}\rangle_{0}=\frac{1}{N}\frac{\gamma}{1-\gamma^{2}}\delta_{il}\delta_{kj}\,.$
(23)

Lastly, for timelike edges connecting unequal parities, we put

$\langle{B_{+}}_{ij}{B_{-}}_{kl}\rangle_{0}=\frac{1}{N}\frac{\gamma}{1-\gamma^{2}}{C_{2}}_{il}{C_{2}}_{kj}\,.$
(24)

In order to make sure that only the four types of edges above (21), (22),
(23), and (24) exist, we set

$\langle{A_{+}}_{ij}{B_{+}}_{kl}\rangle=\langle{A_{+}}_{ij}{B_{-}}_{kl}\rangle=\langle{A_{-}}_{ij}{B_{+}}_{kl}\rangle=\langle{A_{-}}_{ij}{B_{-}}_{kl}\rangle=0\,.$
(25)

The presence of spin-up and spin-down vertices is respectively imposed by
introducing the terms $A_{+}^{2}B_{+}$ and $A_{-}^{2}B_{-}$ in the action.
Therefore,

$S_{\rm CDTIM}=\\\
N\;\mathrm{Tr}\left[\frac{1}{2}A_{+}^{2}+\frac{1}{2}\left({C_{2}}^{-1}B_{+}\right)^{2}+\frac{1}{2}A_{-}^{2}+\frac{1}{2}\left({C_{2}}^{-1}B_{-}\right)^{2}-\gamma
A_{+}A_{-}-\gamma({C_{2}}^{-1}B_{+})({C_{2}}^{-1}B_{-})-gA_{+}^{2}B_{+}-gA_{-}^{2}B_{-}\right]\,.$
(26)

The partition function is then given by

$\mathcal{Z}=\lim_{N\rightarrow\infty}\ln\int
dA_{+}dB_{+}dA_{-}dB_{-}\;e^{-S_{\rm CDTIM}}\;.$ (27)

We can also use similar transformations to (17) and (19), that lead to (20),
to arrive at the version where the spins are on the faces. After the
transformation the action is

$\displaystyle S_{\rm
CDTIM}=N\mathrm{Tr}\left[\frac{1}{2}\gamma^{-1}U_{s}^{2}+\frac{1}{2}\gamma^{-1}(C_{2}^{-1}U_{t})^{2}+\frac{1}{2}\gamma
V_{s}^{2}+\frac{1}{2}\gamma(C_{2}^{-1}V_{t})^{2}-gU_{s}^{2}U_{t}-gV_{s}^{2}U_{t}-gU_{s}V_{s}V_{t}-gV_{s}U_{s}V_{t}\right]\;,$
(28)

where we may interpret that $U_{s}$ generates spacelike edges between faces of
the same spin, $U_{t}$ timelike edges between faces of the same spin, $V_{s}$
spacelike edges between faces of opposite spins, and $V_{t}$ timelike edges
between faces of opposite spins.

##  VI Properties of the matrices $C_{m}$

Let us study the properties of $C_{m}$ and find its expression explicitly.

######  Theorem VI.1.

Let $C_{m}$ be a $GL(N)$ matrix that satisfies
$\mathrm{Tr}[(C_{m})^{q}]=N\delta_{m,q}$ for $q=1,...,N$, with
$N=0\;\mathrm{mod}\;m$. Its eigenvalues can be approximated by

$\lambda_{ts}=e^{\frac{2\pi}{m}is}\;\mathrm{W}(-e^{\frac{2\pi
m}{N}i(t-1/2)-1})^{-\frac{1}{m}}\;\Big{(}1+{\mathcal{O}}(N^{-1/2})\Big{)}\,,\qquad
t=1,...,N/m\;\;\text{and}\;\;s=1,...,m\;.$ (29)

###### Proof.

We show (29) by finding the characteristic polynomial of $C_{m}$. Recall [115]
that, if the characteristic polynomial of a $GL(N)$ matrix $C$ is given by

$p(\lambda)=\det(C-\lambda\mathbb{1})=\sum_{k=0}^{N}\pi_{k}(-\lambda)^{N-k}\,,$
(30)

then the coefficients $\pi_{k}$ are related to the power traces of $C$,
$t_{k}=\mathrm{Tr}\,(C^{k})$, by Girard-Newton formulas

$\pi_{k}=(k!)^{-1}\det[a_{ij}]_{i,j=1,...,k}\,,$ (31)

where

$a_{ij}=\begin{cases}t_{j-i+1}&\text{ if }i\leq j\\\ j&\text{ if }i=j+1\\\
0&\text{ if }i\geq j+1\;.\end{cases}$ (32)

Evaluating for $C=C_{m}$, we have $\pi_{k}=0$ for $k\neq 0$ mod $m$, and for
$k=mr$, with $r$ an integer satisfying $0\leq r\leq k/m$, we get

$\pi_{mr}=\frac{(-1)^{mr}N^{r}}{(-m)^{r}r!}\;.$ (33)

Therefore, imposing that $N$ is divisible by $m$,

$\begin{split}p_{m}(\lambda)&=\sum_{r=0}^{N/m}\pi_{mr}(-\lambda)^{N-mr}=\sum_{r=0}^{N/m}\frac{(-1)^{mr}N^{r}}{(-m)^{r}r!}(-\lambda)^{N-mr}\\\
&=(-\lambda)^{N}\sum_{r=0}^{N/m}\frac{1}{r!}\left(-\frac{N}{m}\lambda^{-m}\right)^{r}=(-\lambda)^{N}s_{N/m}\left(-\frac{N}{m}\lambda^{-m}\right)\,,\end{split}$
(34)

where $s_{n}$ is the nth partial exponential sum polynomial, and we only need
to investigate its zeros because we are interested in eigenvalues
$\lambda_{ts}$. In the large n limit, the polynomials $s_{n}(nz)$ have been
shown [116, 117] to have the roots lying on the curve $|ze^{1-z}|=1$ with
$|z|\leq 1$, and $\mathrm{arg}(ze^{1-z})$ increases monotonically from each
zero. The zeroes $z_{t}$, $t=1,...,n$, are approximately given by
$z_{t}=\tilde{z}_{t}(1+{\mathcal{O}}(N^{-1/2}))$, with

${\tilde{z}}_{t}e^{1-{\tilde{z}}_{t}}=e^{\frac{2\pi}{n}i(t-1/2)}\,,\qquad
t=1,...,n\,.$ (35)

Using the principal branch of the Lambert $W$ function, $W$, defined by
$W^{-1}(z)=ze^{z}$, we invert (35)

${\tilde{z}}_{t}=-\mathrm{W}(-e^{\frac{2\pi}{n}i(t-1/2)-1})\;.$ (36)

Now, in order to go back to eigenvalues $\lambda_{ts}$, we use $n=N/m$ and
$-\lambda_{ts}^{-m}={\tilde{z}}_{t}(1+{\mathcal{O}}(N^{-1/2}))$ in (36), and

$\lambda_{ts}=e^{\frac{2\pi}{m}is}\mathrm{W}(-e^{\frac{2\pi
m}{N}i(t-1/2)-1})^{-\frac{1}{m}}(1+{\mathcal{O}}(N^{-1/2}))\,,\qquad
t=1,...,N/m\;\;\text{and}\;\;s=1,...,m\;.$ (37)

We note that the eigenvalues are discrete, and the set is finite, with $N$
eigenvalues in total. ∎

We let Mathematica solve numerically the zeros of the polynomial (34). We
compare those with (37). We show in Fig. 14 and in Fig. 15 a comparison for
different values of $N$ for $C_{1}$ and $C_{2}$.

(a) $N=4$

(b) $N=20$

(c) $N=100$

(d) $N=500$

Figure 14: Eigenvalues of $C_{1}$ in their complex plane, with vertical axis
being $\mathrm{Im}(\lambda)$ and the horizontal $\mathrm{Re}(\lambda)$. In
blue, we plot the numerical values of the zeros of the polynomial
$p_{1}(\lambda)$ in (34) using Mathematica, and in red we plot the values
obtained from (37) for large $N$.

Figure 15: Eigenvalues of $C_{2}$. We obtained the approximate curve in green
by letting $t$ and $s$ continuous parameters in (37)

.

Taking the values of the roots of (34) computed by Mathematica and evaluating
the power traces $\mathrm{Tr}(C_{1})$ and $\mathrm{Tr}(C_{2})^{2}$, we indeed
obtain the correct results within the machine error. Now, let us also check
the validity of the expression (37), by inserting the values of $\lambda_{ts}$
back to compute the power traces $\mathrm{Tr}(C_{1})$ and
$\mathrm{Tr}(C_{2})^{2}$. We found that for $m\neq q$, we obtain $0$ up to
certain power ($q_{max}$) of $C_{m}$, but the value of the power traces
quickly increases after we reach $q_{max}$. For example, $q_{max}$ is about
$20$ for $N=10000$. We find that the value of $q_{max}$ also increases as $N$
increases.

Let us consider the resolvent of the matrix $C_{m}$

$\frac{1}{N}\mathrm{Tr}\left[\frac{1}{1-\mu\,C_{m}}\right]=\frac{1}{N}\sum_{k=0}^{\infty}\mathrm{Tr}[(\mu\,C_{m})^{k}]=1+\mu^{m}\,.$
(38)

The resolvent (38) evaluates to

$\lim_{N\rightarrow\infty}\sum_{t=1}^{N}\frac{1}{N}\frac{1}{1-\mu\lambda_{t}}=\int_{0}^{1}dx\frac{1}{1-\mu\lambda(x)}=\int_{\gamma}d\lambda\frac{\rho(\lambda)}{1-\mu\lambda}\,,$
(39)

where $x=t/N$ and $\rho(\lambda)=dx/d\lambda$ is called the spectral density,
and the indices $s$ and $t$ are joined to only $t$. By equating (38) and (39),
we can solve for $\rho(\lambda)$.

Let us change the variables, $\alpha=\lambda^{-1}$;

$\int_{\gamma}d\lambda\frac{\rho(\lambda)}{1-\mu\lambda}=\int_{\gamma^{\prime}}d\alpha\,\frac{-\alpha^{-1}\rho(\alpha^{-1})}{\alpha-\mu}\,.$
(40)

Assuming $\gamma^{\prime}$ is a simple closed curve enclosing $\mu$, using
Cauchy’s theorem we can use that

$1+\mu^{m}=\frac{1}{2\pi
i}\int_{\gamma^{\prime}}d\alpha\frac{1+\alpha^{m}}{\alpha-\mu}\;.$ (41)

Equating (40) and (41), we find a solution

$-\alpha^{-1}\rho(\alpha^{-1})=\frac{1}{2\pi i}(1+\alpha^{m})\,,$ (42)

which, in terms of $\lambda$, translates to

$\rho(\lambda)=-\frac{1}{2\pi i}(\lambda^{-1}+\lambda^{-m-1})\;.$ (43)

At this point, let us check the validity of (37). Taking (37) to the power of
$m$,

$\lambda^{m}=\mathrm{W}(-e^{\frac{2\pi m}{N}i(t-1/2)-1})^{-1}\,,$ (44)

to the leading order in large $N$. Manipulating the expression (44),

$-\lambda^{-m}e^{1+\lambda^{-m}}=e^{\frac{2\pi m}{N}i(t-1/2)}=e^{2\pi
mix}\left(1+{\mathcal{O}}(N^{-1})\right)\,,$ (45)

where we change variables $x=t/N$. Taking a derivative with respect to
$\lambda$,

$(-m\lambda^{-1}-m\lambda^{-m-1})(-\lambda^{-m}e^{1+\lambda^{-m}})=2\pi
mi\frac{dx}{d\lambda}e^{2\pi mix}\left(1+{\mathcal{O}}(N^{-1})\right)\;\,,$
(46)

and noting that $\rho(\lambda)=dx/d\lambda$, we again obtain the same solution
(43). This concludes checking the validity of the solution (37).

Let us go back to (43) and continue solving for the distribution of $\lambda$.
Inserting $\rho(\lambda)=dx/d\lambda$, where $x\in{\mathbb{R}}$ with $0\leq
x\leq 1$ in (43) and integrating both sides, we obtain

$\lambda^{-m}e^{\lambda^{-m}}=\kappa\,e^{2\pi mix}\,,$ (47)

where $\kappa$ is some integration constant. In (37), $\kappa$ was determined
to be $\kappa=-e^{-1}$. Therefore, in Fig. 14 and Fig. 15 we only had one
curve $\gamma$ each. However, (47) gives us a series of solutions because
$\kappa$ is free. A larger class of solutions means that now we should have
many curves representing families of eigenvalues, with each curve
corresponding to a particular value of $\kappa$. The observation of (47) is
that, given two possible values of $m$, $m_{1}$ and $m_{2}$, their solutions
are related by $\lambda_{1}^{m_{1}}=\lambda_{2}^{m_{2}}$. Therefore, we can
focus on $m=1$ without losing any essential information. For $m=1$, we have

$\lambda^{-1}e^{\lambda^{-1}}=\kappa\,e^{2\pi ix}\;.$ (48)

Then, each value of $\kappa$ is associated with a distinct solution for
$C_{m}$. The eigenvalues of $C_{m}$ are represented by a curve
$\gamma_{\kappa}$. Equation (48) can be inverted with the help of the $W$
Lambert function,

$\lambda=W_{b}(\kappa\,e^{2\pi ix})^{-1}\,,$ (49)

where $W_{b}$ is a branch of $W$ and $b$ is an integer we assign to each
branch. In Fig. 16, we show a few of the infinitely many branch regions of the
$W$ function with $z=\lambda^{-1}$. Let us assign $b=0$ to the branch where
the point $z=0$ belongs to. As is shown in Fig. 16(b), only curves around the
origin can be closed. These closed curves are the contours associated with the
values of $|\kappa|=|ze^{z}|$ with $0\leq|\kappa|\leq e^{-1}$. The maximum
value $|\kappa|=e^{-1}$ is associated with the branch point of $W_{0}$. Our
previous solution (45) picked this particular value $|\kappa|=e^{-1}$.
However, our current analysis indicates that all these closed curves around
the origin are valid not just $|\kappa|=e^{-1}$. In the previous solution
(45), the condition $\kappa=-e^{-1}$ appeared because of the polynomials
$s_{N/m}$ in (34). We speculate that the reason for the existence of multiple
solutions is related to the fact that the large $N$ limit condition
$\lim_{N\rightarrow\infty}\frac{1}{N}\mathrm{Tr}[(C_{m})^{p}]=\delta_{p,m}$
can be approached through different ways. That is, $\kappa$ may be a function
of $\mathrm{Tr}[C_{m}]$ and $\mathrm{Tr}[(C_{m})^{m}]$, and this function
might not have trivial limits. For example,
$\mathrm{Tr}[{(\tilde{C}_{m})}^{p}]=N\delta_{p,m}+\delta_{p,1}$ has the same
large $N$ limit,
$\lim_{N\rightarrow\infty}\frac{1}{N}\mathrm{Tr}((\tilde{C}_{m})^{p})=\delta_{p,m}$.
The value of $\kappa$ may vary when changing the condition on
$\mathrm{Tr}C_{m}$.

(a) Contour lines of $\mathrm{arg}(ze^{z})$. Complex plane of $z$
(${\mathrm{Re}}(z)$ on horizontal, and ${\mathrm{Im}}(z)$ on vertical).

(b) Contour lines of $\mathrm{ln}{|ze^{z}|}$. Complex plane of $z$
(${\mathrm{Re}}(z)$ on horizontal, and ${\mathrm{Im}}(z)$ on vertical).

Figure 16: Plots for the function $ze^{z}$. On the white curves, the image of
$ze^{z}$ is on the negative real line. The black point is the branch point.

##  VII Character expansion of partition function

For a $N\times N$ matrix $M$ with eigenvalues $M_{k}$, where $k\in\mathbb{Z}$
with its range $1\leq k\leq N$ and given a set of $N$ increasing non-negative
integers $\\{h\\}=\\{h_{1},...,h_{N}\\}$, the character of this matrix can be
given by

$\chi_{\\{h\\}}(M)=\frac{\det_{(k,l)}[(M_{k})^{h_{l}}]}{\Delta(M)}\,,$ (50)

where

$\Delta(M)=\det_{(k,l)}[(M_{k})^{l-1}]=\prod_{1\leq i<j\leq N}(M_{j}-M_{i})$
(51)

is the Vandermonde determinant. In general, the set $\\{h\\}$ can be used to
uniquely label a representation of $\mathrm{GL}(N)$. We have as a result of
character expansion [78]

$e^{\mathrm{Tr}M^{2}}\sim\sum_{\\{h\\}}c_{\\{h\\}}\chi_{\\{h\\}}(M)\,,$ (52)

with a summation over all such sets $\\{h\\}$ and where, as we use from now
on, $\sim$ indicates equality up to a constant. We separate the set $\\{h\\}$
into a set of only even integers $\\{h\\}^{e}$ and of only odd integers
$\\{h\\}^{o}$. If the numbers of even and odd $h_{l}$ are equal or the numbers
of even is one more than the number of odd, the coefficient $c_{\\{h\\}}$ is
given by

$c_{\\{h\\}}=\frac{\Delta(\\{h\\}^{e})\Delta(\\{h\\}^{o})}{\prod_{l=1}^{N}\left\lfloor\\{h_{l}\\}/2\right\rfloor!}\,,$
(53)

where, in a similar manner to (51), for a set of complex numbers
$\\{m\\}=\\{m_{1},...,m_{N}\\}$, the Vandermonde is defined as

$\Delta(\\{m\\})=\det_{(k,l)}[(m_{k})^{l-1}]=\prod_{1\leq i<j\leq
N}(m_{j}-m_{i})\;.$ (54)

We use $\lfloor{\cdot}\rfloor$ to represent the floor function. When the
condition above is not met, we simply have

$c_{\\{h\\}}=0\;.$ (55)

The character expansion for the more general case $e^{\mathrm{Tr}M^{k}}$, for
integer $k\geq 1$ is also known [78].

###  VII.1 Character expansion for the pure CDT-like matrix model of
Benedetti-Henson

The partition function (1) for the pure CDT-like matrix model due to
Benedetti-Henson studied in [76], after integrating out $B$ is given by

$Z=\int
dA\;e^{-N\mathrm{Tr}[\frac{1}{2}A^{2}-\frac{g^{2}}{2}(A^{2}C_{2})^{2}]}\;.$
(56)

By applying the character expansion given in (52) with
$M=\sqrt{\frac{N}{2}g^{2}}A^{2}C_{2}$, we expand the partition function (56),

$Z\sim\sum_{\\{h\\}}\left(N\frac{g^{2}}{2}\right)^{\\#h/2}c_{\\{h\\}}\int
dA\;\chi_{\\{h\\}}(A^{2}C_{2})e^{-N\mathrm{Tr}[\frac{1}{2}A^{2}]}\,,$ (57)

where $\\#h=\sum_{l=1}^{N}h_{l}-\frac{1}{2}N(N-1)$. With the diagonalization
$A=\Omega\Lambda\Omega^{\dagger}$, with $\Omega\in\mathrm{U}(N)$ and $\Lambda$
a diagonal matrix, we preform a change of variables and rewrite

$Z\sim\sum_{\\{h\\}}\left(N\frac{g^{2}}{2}\right)^{\\#h/2}c_{\\{h\\}}\int_{\mathbb{R}^{N}}d\Lambda\;\Delta(\Lambda)^{2}e^{-N\mathrm{Tr}[\frac{1}{2}\Lambda^{2}]}\int_{U(N)}d\Omega\;\chi_{\\{h\\}}(\Omega\Lambda^{2}\Omega^{\dagger}C_{2})\;.$
(58)

In order to factorize the dependence on $C_{2}$, let us use the Schur
orthogonality relation,

$\int_{U(N)}d\Omega\;\chi_{\\{h\\}}(\Omega\Lambda^{2}\Omega^{\dagger}C_{2})=\frac{\chi_{\\{h\\}}(\Lambda^{2})\chi_{\\{h\\}}(C_{2})}{d_{h}}\,,$
(59)

where $d_{h}=\Delta(h)/\prod_{i=1}^{N-1}i!$ is the dimension of the
representation given by $\\{h\\}$. We can then rewrite

$Z\sim\sum_{\\{h\\}}\left(N\frac{g^{2}}{2}\right)^{\\#h/2}\frac{c_{\\{h\\}}}{d_{h}}\,\chi_{\\{h\\}}(C_{2})\int_{\mathbb{R}^{N}}d\Lambda\;\Delta(\Lambda)^{2}\;\chi_{\\{h\\}}(\Lambda^{2})\;e^{-N\mathrm{Tr}[\frac{1}{2}\Lambda^{2}]}\,,$
(60)

now with the integral independent of $C_{2}$ as desired, and $c_{\\{h\\}}$ is
given in (53) or $0$. By changing the variables back to
$A=\Omega\Lambda\Omega^{\dagger}$, the integral part of the expression (60)
can be written, up to a constant, as

$\langle\chi_{\\{h\\}}(A^{2})\rangle_{0}=\frac{1}{Z_{0}}\int
dA\;\chi_{\\{h\\}}(A^{2})\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\;.$ (61)

This integral had yet no known general solution in terms of the representation
$\\{h\\}$. A conjecture was also given in [76]. It was proposed that

$\langle\chi_{\\{h\\}}(A^{2})\rangle_{0}=k_{N}\frac{1}{N^{\\#h}}\prod_{\epsilon=0}^{3}\Delta(2h^{(\epsilon)})^{2}\prod_{i}(2h_{i}^{(\epsilon)})!!\;,$
(62)

where the set of integers $\\{h\\}$ has been divided into four sets
$\\{h^{(0)}\\},...,\\{h^{(3)}\\}$ according to equivalence modulo 4.

The conjecture stated in (62), given in [76], can be rewritten as

$k_{N}=\frac{\langle\chi_{\\{h\\}}(A^{2})\rangle_{0}}{\frac{1}{N^{\\#h}}\prod_{\epsilon=0}^{3}\Delta(2\\{h\\}^{(\epsilon)})^{2}\prod_{i}(2h_{i}^{(\epsilon)})!!}\;,$
(63)

where
$\\{h\\}^{(\epsilon)}=\\{a\in{\\{h\\}}\;|\;a=\epsilon\;\mathrm{mod}\;4\\}$ and
for a set $B$, we define $2B=\\{2b\,|\,b\in B\\}$. $k_{N}$ only depends on
$N$, and not on $\\{h\\}$. We can test this conjecture by applying it to some
representations which we already know the result. Fix $N=4$. For the trivial
representation, which is given by $\\{h\\}=\\{3,2,1,0\\}$, the conjecture (63)
tells us that the proportionality factor is $k_{N}=1/768$. On the other hand,
for the defining representation, which is given by $\\{h\\}=\\{4,2,1,0\\}$,
the conjecture tells us that the proportionality factor should be
$k_{N}=1/(8\cdot 768)$. Therefore, a contradiction.

Furthermore, let us look at representations that have size proportional to
$N$. For example, the $k$th power determinant representation is one. The
representation $\\{h\\}$ for which
$\chi_{\\{h\\}}(A^{2})=\mathrm{det}(A^{2q})$ is given by
$\\{h\\}=\\{q+N-1,q+N-2,...,q+1,q\\}$. Evaluating the average integral (14)
for $f(A)=\mathrm{det}(A^{2q})$ we obtain

$\langle\mathrm{det}A^{2q}\rangle_{0}=N^{-Nq}\prod_{k=0}^{\frac{N}{2}-1}\frac{(2q+2k+1)!!(2q+2k-1)!!}{(2k+1)!!(2k-1)!!}$
(64)

According to (63), we would have

$k_{N}=\frac{\prod_{k=0}^{\frac{N}{2}-1}\frac{(2q+2k+1)!!(2q+2k-1)!!}{(2k+1)!!(2k-1)!!}}{2^{N}2^{N\left(\frac{N}{4}-1\right)}\left(\prod_{i=0}^{\frac{N}{4}-1}i!\right)^{8}\prod_{k=1}^{N}(q+N-k)!}\;.$
(65)

We then see that $k_{N}$ depends on the representation through $q$, which
again goes against the assumption that $k_{N}$ only depends on $N$.

###  VII.2 Character expansion for CDT-like matrix model coupled with Ising
model

For CDT-like matrix model coupled with Ising model, the expression is more
involved. We write its partition function as

$Z=\int
dA_{1}dA_{2}\;e^{-N\mathrm{Tr}[\frac{1}{2}A_{1}^{2}+\frac{1}{2}A_{2}^{2}-\gamma
A_{1}A_{2}-\frac{1}{2}\frac{g^{2}}{1-\gamma}(A_{1}^{2}C_{2})^{2}-\frac{1}{2}\frac{g^{2}}{1-\gamma}(A_{2}^{2}C_{2})^{2}-\frac{1}{2}\frac{\gamma
g^{2}}{1-\gamma^{2}}((A_{1}^{2}+A_{2}^{2})C_{2})^{2}]}\,,$ (66)

and diagonalize the kinetic part,

$Z=\int
dUdV\;e^{-N\mathrm{Tr}[\frac{1}{2}(1-\gamma)U^{2}+\frac{1}{2}(1+\gamma)V^{2}-\frac{1}{4}\frac{g^{2}}{1-\gamma}((U^{2}+V^{2})C_{2})^{2}-\frac{1}{4}\frac{g^{2}}{1+\gamma}((UV+VU)C_{2})^{2}]}\,,$
(67)

where $U=(A_{1}+A_{2})/\sqrt{2}$ and $V=(A_{1}-A_{2})/\sqrt{2}$. The character
expansion (52) for the latter partition function (67)

$Z\sim\sum_{\\{h^{1}\\},\\{h^{2}\\}}\left(\frac{Ng^{2}}{4(1-\gamma)}\right)^{\\#h^{1}/2}\left(\frac{Ng^{2}}{4(1+\gamma)}\right)^{\\#h^{2}/2}c_{\\{h^{1}\\}}c_{\\{h^{2}\\}}\;I_{\\{h^{1}\\},\\{h^{2}\\}}$
(68)

where the sum is performed over representations labeled by $\\{h^{1}\\}$ and
$\\{h^{2}\\}$, and we defined the integral

$I_{\\{h^{1}\\},\\{h^{2}\\}}=\int
dUdV\;e^{-N\mathrm{Tr}[\frac{1}{2}(1-\gamma)U^{2}+\frac{1}{2}(1+\gamma)V^{2}]}\;\chi_{\\{h^{1}\\}}((U^{2}+V^{2})C_{2})\;\chi_{\\{h^{2}\\}}((UV+VU)C_{2})\;.$
(69)

Transform $U=\Omega_{1}\Lambda_{1}\Omega_{1}^{\dagger}$ and
$V=\Omega_{2}\Lambda_{2}\Omega_{2}^{\dagger}$, and furthermore, by
substituting $\Omega_{2}=\Omega_{1}\Omega$, we have

$\begin{split}I_{\\{h^{1}\\},\\{h^{2}\\}}=&\int_{\mathbb{R}^{N}\times\mathbb{R}^{N}}d\Lambda_{1}d\Lambda_{2}\;\Delta(\Lambda_{1})^{2}\Delta(\Lambda_{2})^{2}e^{-N\mathrm{Tr}[\frac{1}{2}(1-\gamma)\Lambda_{1}^{2}+\frac{1}{2}(1+\gamma)\Lambda_{2}^{2}]}\\\
&\int_{U(N)\times
U(N)}d\Omega_{1}d\Omega\;\chi_{\\{h^{1}\\}}(\Omega_{1}(\Lambda_{1}^{2}+\Omega\Lambda_{2}^{2}\Omega^{{\dagger}})\Omega_{1}^{{\dagger}}C_{2})\;\chi_{\\{h^{2}\\}}(\Omega_{1}(\Lambda_{1}\Omega\Lambda_{2}\Omega^{{\dagger}}+\Omega\Lambda_{2}\Omega^{{\dagger}}\Lambda_{1})\Omega_{1}^{{\dagger}}C_{2})\,,\end{split}$
(70)

that allows us to integrate over $\Omega_{1}$, however, the orthogonality
relation (59) is not enough, but one has to be equipped with a more general
expression for Schur orthogonality relation.

##  VIII Unitary integral of degree $4$ monomial

A more general way to state the Schur orthogonality relation is that for a
compact group $G$ and for matrix elements $\phi^{\alpha}_{ij}(g)$ of an
element $g\in G$, where the upper index indicates the representation and the
lower indices indicate the corresponding elements of a given orthogonal basis,

$\int_{G}\phi^{\alpha}_{ij}(g)\overline{\phi^{\beta}_{kl}(g)}dg=\frac{1}{d_{\alpha}}\delta_{\alpha\beta}\delta_{ik}\delta_{jl}.$
(71)

Here, with $G$ as $U(N)$ and for $\alpha=\beta=h$ by multiplying both sides by
$\phi^{h}_{jl}(\Lambda^{2})$ and $\phi^{h}_{ki}(C_{2})$ and summing over the
repeated lower indices, we recover (59). Let us apply (71) to the unitary
integral of (70), by defining the matrices which are independent of
$\Omega_{1}$ as
$M_{1}=\Lambda_{1}^{2}+\Omega\Lambda_{2}^{2}\Omega^{{\dagger}}$ and
$M_{2}=\Lambda_{1}\Omega\Lambda_{2}\Omega^{{\dagger}}+\Omega\Lambda_{2}\Omega^{{\dagger}}\Lambda_{1}$.
Then, we can factorize the expression

$\begin{split}\chi_{\\{h^{1}\\}}(\Omega_{1}M_{1}\Omega_{1}^{{\dagger}}C_{2})\;\chi_{\\{h^{2}\\}}(\Omega_{1}M_{2}\Omega_{1}^{{\dagger}}C_{2})=\;\phi^{h^{1}}_{ab}(\Omega_{1})\phi^{h^{1}}_{bc}(M_{1})\phi^{h^{1}}_{cd}(\Omega^{{\dagger}}_{1})\phi^{h^{1}}_{da}(C_{2})\phi^{h^{2}}_{pq}(\Omega_{1})\phi^{h^{2}}_{qr}(M_{2})\phi^{h^{2}}_{rs}(\Omega_{1}^{\dagger})\phi^{h^{2}}_{sp}(C_{2})\;.\end{split}$
(72)

Then, in (70), we have an integral

$I_{ij}=\int_{U(N)}\phi^{\alpha}_{ab}(g)\phi^{\beta}_{cd}(g)\overline{\phi^{\alpha}_{\bar{a}\bar{b}}(g)\phi^{\beta}_{\bar{c}\bar{d}}(g)}dg\,,$
(73)

where $g=\Omega_{1}$. This integral has some properties similar to the
previous one (71), but has the double amount of matrix elements. This integral
is also similar to the type of integral studied with the Weingarten calculus,
however, here, more than one representations $\alpha$ and $\beta$ appear.

Putting all together, (70) becomes

$I_{\\{h^{1}\\},\\{h^{2}\\}}=I_{ij}\;\phi^{h^{1}}_{da}(C_{2})\,\phi^{h^{2}}_{sp}(C_{2})\,\times\\\
\times\int_{\mathbb{R}^{N}\times\mathbb{R}^{N}}d\Lambda_{1}d\Lambda_{2}\;\Delta(\Lambda_{1})^{2}\Delta(\Lambda_{2})^{2}e^{-N\mathrm{Tr}[\frac{1}{2}(1-\gamma)\Lambda_{1}^{2}+\frac{1}{2}(1+\gamma)\Lambda_{2}^{2}]}\int_{U(N)}d\Omega\;\phi^{h^{1}}_{bc}(M_{1})\phi^{h^{2}}_{qr}(M_{2})\,,$
(74)

and with the theorem VIII.1 below, we can evaluate $I_{ij}$.

######  Proposition VIII.1.

Given $g\in U(N)$ and a representation $\alpha$, set its matrix elements as
$\phi^{\alpha}_{ab}(g)$, where the indices $a$ and $b$ are associated with
some chosen basis. The following holds true for a Haar integral:

$\int_{U(N)}\phi^{\alpha}_{ab}(g)\phi^{\beta}_{cd}(g)\overline{\phi^{\alpha}_{\bar{a}\bar{b}}(g)\phi^{\beta}_{\bar{c}\bar{d}}(g)}dg=\sum_{r,m,n,p,q}{c^{rm*}_{acp}}c^{rn}_{\bar{a}\bar{c}p}{c^{rm}_{bdq}}c^{rn*}_{\bar{b}\bar{d}q}d_{r}^{-1}\;.$
(75)

###### Proof.

In order to apply the Weingarten calculus [118] to the left hand side of (75),
we need to know a basis for the invariant space of
$\alpha\,\otimes\,\beta\,\otimes\,\bar{\alpha}\,\otimes\,\bar{\beta}$. Indeed,
finding the basis for the invariant space amounts to identifying the trivial
representations $1^{\oplus l_{1}^{\alpha\beta\bar{\alpha}\bar{\beta}}}$, upon
decomposing the tensor product representation into irreducible
representations,

$\alpha\,\otimes\,\beta\,\otimes\,\bar{\alpha}\,\otimes\,\bar{\beta}=\bigoplus_{\mathrm{all}\,\mathrm{irrep}\,r}r^{\oplus
l_{r}^{\alpha\beta\bar{\alpha}\bar{\beta}}}\,,$ (76)

where $l_{r}^{\alpha\beta\bar{\alpha}\bar{\beta}}$ is the number of times the
representation $r$ appears in the expansion of the tensor product, known as
Littlewood-Richardson coefficients. One can decompose separately the products
$\alpha\,\otimes\,\beta$ and $\bar{\alpha}\,\otimes\,\bar{\beta}$, i.e.,

$\alpha\,\otimes\,\beta=\bigoplus_{r}r^{\oplus
l_{r}^{\alpha\beta}}\,,\qquad\bar{\alpha}\,\otimes\,\bar{\beta}=\bigoplus_{\bar{r}}\bar{r}^{\oplus
l_{\bar{r}}^{\bar{\alpha}\bar{\beta}}}=\bigoplus_{r}\bar{r}^{\oplus
l_{\bar{r}}^{\bar{\alpha}\bar{\beta}}}\,.$ (77)

Then, using $l_{\bar{r}}^{\bar{\alpha}\bar{\beta}}=l_{r}^{\alpha\beta}$, we
can write

$\displaystyle\alpha\,\otimes\,\beta\,\otimes\,\bar{\alpha}\,\otimes\,\bar{\beta}$
$\displaystyle=$ $\displaystyle\bigoplus_{r,s}r^{\oplus
l_{r}^{\alpha\beta}}\,\otimes\,\bar{s}^{\oplus
l_{\bar{s}}^{\bar{\alpha}\bar{\beta}}}$ (78) $\displaystyle=$
$\displaystyle\bigoplus_{r}r^{\oplus
l_{r}^{\alpha\beta}}\,\otimes\,\bar{r}^{\oplus
l_{r}^{\alpha\beta}}\,\oplus\,\bigoplus_{r\neq s}r^{\oplus
l_{r}^{\alpha\beta}}\,\otimes\,\bar{s}^{\oplus l_{s}^{\alpha\beta}}\,,$ (79)

but the representation orthogonality relation tells us that for two
irreducible representations $r$ and $s$, $l_{1}^{r\bar{s}}=\delta_{rs}$.
Therefore, we notice that the second big sum contains no trivial
representations and the first contains one trivial representation for each
product. In the end, we arrive at

$\alpha\,\otimes\,\beta\,\otimes\,\bar{\alpha}\,\otimes\,\bar{\beta}=\bigoplus_{r}1^{\oplus(l_{r}^{\alpha\beta})^{2}}\oplus\bigoplus_{\begin{subarray}{c}\mathrm{all}\,\mathrm{irrep}\,r\\\
r\neq 1\end{subarray}}r^{\oplus
l_{r}^{\alpha\beta\bar{\alpha}\bar{\beta}}}\;.$ (80)

Therefore, we conclude

$l_{1}^{\alpha\beta\bar{\alpha}\bar{\beta}}=\sum_{r}(l_{r}^{\alpha\beta})^{2}\,.$
(81)

Remark that there are infinite number of representations to sum over, however,
only finite number of them have nonzero $l_{r}^{\alpha\beta}$. Equation (81)
is important to us because it allows to verify that the set we find later in
(94) for the invariant space is indeed a basis because it contains as many
elements as the dimension of the invariant space.

Let us analyze a basis for the invariant space
$\bigoplus_{r}1^{\oplus(l_{r}^{\alpha\beta})^{2}}$, since it allows us to
evaluate $I_{ij}$ in (73). First, let us define a basis for the vector space
of the tensor product representation
$\alpha\,\otimes\,\beta\,\otimes\,\bar{\alpha}\,\otimes\,\bar{\beta}$,

$|E_{i}\rangle=|E_{ab\,\bar{a}\bar{b}}\rangle=|e_{a}\rangle|e_{b}\rangle\langle
e_{\bar{a}}|\langle e_{\bar{b}}|\,,$ (82)

where $a$, $\bar{a}$, $b$, $\bar{b}$, and $i$ are integers and $i$ is a
relabelling of indices $ab\bar{a}\bar{b}$ with $1\leq
i\leq(d_{\alpha})^{2}(d_{\beta})^{2}$, and $|e_{a}\rangle$, $|e_{b}\rangle$,
$\langle e_{\bar{a}}|$ and $\langle e_{\bar{b}}|$, with $1\leq a\leq
d_{\alpha}$, $1\leq\bar{a}\leq d_{\alpha}$, $1\leq b\leq d_{\beta}$,
$1\leq\bar{b}\leq d_{\beta}$, are respectively a basis for representations
$\alpha$, $\beta$, $\bar{\alpha}$ and $\bar{\beta}$. Writing $T_{\mu}$ as a
basis for this invariant space, where $1\leq\mu\leq
l_{1}^{\alpha\beta\bar{\alpha}\bar{\beta}}$, let us define the $A$ matrix

$A_{i\mu}=\langle E_{i}|T_{\mu}\rangle\,,$ (83)

and the Weingarten matrix $w$ through the inverse of a matrix

$w^{-1}_{\mu\nu}=\langle T_{\mu}|T_{\nu}\rangle\;.$ (84)

The Weingarten theorem states that (73) can be expressed as

$I_{ij}=\sum_{\mu,\nu}A_{i\mu}w_{\mu\nu}A^{\dagger}_{\nu j}\,,$ (85)

which is a function of $T$. $T$ being invariant under
$\alpha\,\otimes\,\beta\,\otimes\,\bar{\alpha}\,\otimes\,\bar{\beta}$ is
equivalent to

$MTM^{-1}=T\,,$ (86)

for every $M\in\alpha\,\otimes\,\beta$. Furthermore, there exists a $V\in
U(d_{\alpha}\,d_{\beta})$ such that, for every $M$, $V$ transforms $M$ into a
block diagonal matrix,

$VMV^{-1}=M_{D},$ (87)

where,

$M_{D}=\begin{pmatrix}M_{r_{1}}&0&\cdots&0\\\ 0&M_{r_{2}}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&M_{r_{n}}\\\ \end{pmatrix}\,,$ (88)

with $n=\sum_{r}l_{r}^{\alpha\beta}$. For every representation $r$ appearing
in the decomposition of $\alpha\,\otimes\,\beta$ as in (77), $M_{r}$ is $M$ in
the specific representation $r$ and appears as a block element in the diagonal
of $M_{D}$, appearing as many times as it appears in the decomposition (77),
i.e., $l_{r}^{\alpha\beta}$ times. Now, transform $T$ by a change of basis by
$V$

$VTV^{-1}=T_{D}\,,$ (89)

which implies

$M_{D}T_{D}M_{D}^{-1}=T_{D}\;.$ (90)

$T_{D}$ can be found by focusing on the sector of a representation $r$ that
appears in the decomposition, where a sector is the block consisting of
$l_{r}^{\alpha\beta}$ number of block diagonal $M_{r}$ matrices. Evaluating

$M_{D}^{-1}=\begin{pmatrix}M_{r_{1}}^{-1}&0&\cdots&0\\\
0&M_{r_{2}}^{-1}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&M_{r_{n}}^{-1}\\\ \end{pmatrix}\,,$ (91)

and using properties of block matrix multiplication, each section can act
independently. As an illustration, suppose that $r$ appears twice, i.e.,
$l_{r}^{\alpha\beta}=2$,

$M_{D}=\begin{pmatrix}\ddots&0&0&0\\\ 0&M_{r}&0&0\\\ 0&0&M_{r}&0\\\
0&0&0&\ddots\\\ \end{pmatrix}\;.$ (92)

Then, denoting $\mathbb{1}_{r}$ as the identity in the representation $r$, it
is easy to see that

$T_{D}=\begin{pmatrix}\ddots&0&0&0\\\
0&\gamma_{11}\mathbb{1}_{r}&\gamma_{12}\mathbb{1}_{r}&0\\\
0&\gamma_{21}\mathbb{1}_{r}&\gamma_{22}\mathbb{1}_{r}&0\\\ 0&0&0&\ddots\\\
\end{pmatrix}$ (93)

satisfies (90) for any complex numbers $\gamma_{11}$, $\gamma_{12}$,
$\gamma_{21}$, and $\gamma_{22}$. In general, (93) can be written as

$T_{D}=\gamma_{mn}^{r}T_{r}^{mn}\,,$ (94)

where for any representation $r$, that appears $l_{r}^{\alpha\beta}$ times, we
can construct $(l_{r}^{\alpha\beta})^{2}$ independent matrices $T_{r}^{mn}$ by
putting an identity matrix in an element of the $l_{r}^{\alpha\beta}\times
l_{r}^{\alpha\beta}$ matrix block corresponding to the representation $r$. If
we consider all representations, we find as many independent invariants as
(81), thus we have a complete basis for the invariant space. Using
$T_{r}^{mn}T_{s}^{pq}=\delta_{rst}\delta^{np}T^{mq}_{t}$ and
${T_{r}^{mn}}^{\dagger}=T_{r}^{nm}$ (where $\delta^{rst}=1$ for $r=s=t$ and
zero otherwise), the graham matrix elements are given by

$\langle
T_{s}^{pq}|T_{r}^{mn}\rangle=\mathrm{Tr}\;{T_{s}^{pq}}^{\dagger}T_{r}^{mn}=\mathrm{Tr}\;{T_{s}^{qp}}T_{r}^{mn}=\delta_{srt}\delta^{pm}\mathrm{Tr}\;T^{qn}_{t}=\delta_{srt}\delta^{pm}\delta^{qn}d_{t}\,,$
(95)

where repeated indices are summed over. Consequently, $T_{r}^{mn}$ is an
orthogonal basis and

$\langle T_{r}^{mn}|T_{r}^{mn}\rangle=d_{r}\,,$ (96)

where $d_{r}$ is the dimension of the representation $r$. Recalling (84), we
can readily invert (95) to obtain

$w_{rs}^{mnpq}=\delta_{srt}\delta^{pm}\delta^{qn}d_{t}^{-1}\;.$ (97)

We also need to calculate the $A$ matrix defined in (83), but from (89) we see
that

$|T_{\mu}\rangle=V^{-1}|T_{r}^{mn}\rangle\,.$ (98)

Therefore, the $A$ matrix is given by

$A_{i\mu}=\langle
E_{i}|V^{-1}|T_{r}^{mn}\rangle=\mathrm{Tr}\;E_{i}^{\dagger}V^{-1}T_{r}^{mn}V=\mathrm{Tr}\left[\left(VE_{i}V^{-1}\right)^{\dagger}T_{r}^{mn}\right]\,,$
(99)

where in the last expression, we insert

$VE_{i}V^{-1}=V|e_{a}\rangle|e_{b}\rangle\langle e_{\bar{a}}|\langle
e_{\bar{b}}|V^{-1}=V|e_{a}\rangle|e_{b}\rangle(V|e_{\bar{a}}\rangle|e_{\bar{b}}\rangle)^{\dagger}\;.$
(100)

Notice that we can also write, now being more explicit on the representation,

$V|e_{a}^{\alpha}\rangle|e^{\beta}_{b}\rangle=c^{rk}_{abp}|e^{rk}_{p}\rangle\,,$
(101)

where $c^{rk}_{abp}$ are the Clebsch-Gordan coefficients, and the
representations $\alpha$ and $\beta$ are written explicitly on the left but
not on the right. We also introduce a new index $k$, an integer which
satisfies $1\leq k\leq l_{r}^{\alpha\beta}$ to distinguish the same
representation $r$ which appears $l_{r}^{\alpha\beta}$ times in the
decomposition of $\alpha\otimes\beta$. Then, we can compute (99)

$A_{i\mu}=\mathrm{Tr}\;{c^{sk*}_{abp}}c^{\bar{s}\bar{k}}_{\bar{a}\bar{b}\bar{p}}|e^{\bar{s}\bar{k}}_{\bar{p}}\rangle\langle
e^{sk}_{p}|T_{r}^{mn}={c^{sk*}_{abp}}c^{\bar{s}\bar{k}}_{\bar{a}\bar{b}\bar{p}}\langle
e^{sk}_{p}|T_{r}^{mn}|e^{\bar{s}\bar{k}}_{\bar{p}}\rangle$ (102)

and using the coefficients for $T_{r}^{mn}$ in this basis, $\langle
e^{sk}_{p}|T_{r}^{mn}|e^{\bar{s}\bar{k}}_{\bar{p}}\rangle=\delta^{s\bar{s}r}\delta_{p\bar{p}}\delta^{km}\delta^{\bar{k}n}$,
we further compute

$A_{i\mu}=c^{sk*}_{abp}c^{\bar{s}\bar{k}}_{\bar{a}\bar{b}\bar{p}}\delta^{rs\bar{s}}\delta_{p\bar{p}}\delta^{km}\delta^{\bar{k}n}=\delta^{rs\bar{s}}{c^{sm*}_{abp}}c^{\bar{s}n}_{\bar{a}\bar{b}p}={c^{rm*}_{abp}}c^{rn}_{\bar{a}\bar{b}p}\,,$
(103)

where in the last expression there is a sum only over $p$. However, remark
that there is no known formula for these coefficients for a general
representation, therefore one cannot compute explicitly $A_{i\mu}$ with the
expression above (103). Nevertheless, inserting (103) in the expression (85),
we get

$I_{ij}=\sum_{\mu,\nu}A_{i\mu}w_{\mu\nu}A^{\dagger}_{\nu
j}=\sum_{r,m,n,t,u,v}\delta^{rs\bar{s}}{c^{sm*}_{abp}}c^{\bar{s}n}_{\bar{a}\bar{b}p}\delta^{tk\bar{k}}{c^{ku}_{bdq}}c^{\bar{k}v*}_{\bar{b}\bar{d}q}\delta_{rtg}\delta^{um}\delta^{vn}d_{g}^{-1}\,,$
(104)

Performing some of the sums, we find that $I_{ij}$ in (73) can be expressed in
terms of the (unknown) Clebsch-Gordan coefficients as

$I_{ij}=\sum_{r,m,n,p,q}{c^{rm*}_{abp}}c^{rn}_{\bar{a}\bar{b}p}{c^{rm}_{bdq}}c^{rn*}_{\bar{b}\bar{d}q}d_{r}^{-1}\,.$
(105)

∎

We comment that the result in (105) is expressed in terms of the Clebsch-
Gordan coefficients which are base dependent. The possibility of choosing
different bases shows that there is a degree of freedom in the summation. This
freedom is manifest in the summation over the indices $p$ and $q$. Further
simplification could be achieved by expressing (105) in a base invariant way.

##  IX  Evaluation of $\langle\chi_{r}(A)\rangle_{0}$ and
$\langle\chi_{r}(A^{2})\rangle_{0}$

In this section, we study and present several computational results on the
large $N$ limit of the Hermitian Gaussian matrix model averages
$\langle\chi_{r}(A)\rangle_{0}$ and $\langle\chi_{r}(A^{2})\rangle_{0}$ as the
latter appeared as a key quantity as in (61) and the former can be an useful
exercise and preparation for computing the latter. The Gaussian average of a
character of $A$ for a given representation $\\{h\\}$ is given by

$\langle\chi_{\\{h\\}}(A)\rangle_{0}=\frac{1}{Z_{0}}\int
dA\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\;\chi_{\\{h\\}}(A)\,,$ (106)

where $Z_{0}$ is an Gaussian integral of $A$;

$Z_{0}=\int dA\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\,.$ (107)

Similarly to (106), the Gaussian average of a character of $A^{2}$ for a given
representation, given by

$\langle\chi_{\\{h\\}}(A^{2})\rangle_{0}=\frac{1}{Z_{0}}\int
dA\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\;\chi_{\\{h\\}}(A^{2})\,.$ (108)

For $\\{h\\}^{e}=\\{a\in\\{h\\}|\,a\;\mathrm{is\;even}\\}$ and
$\\{h\\}^{o}=\\{a\in\\{h\\}|\,a\;\mathrm{is\;odd}\\}$, (106) evaluates to
[119]

$\langle\chi_{\\{h\\}}(A)\rangle_{0}=N^{-\frac{n}{2}}\;d_{h}\;\frac{\prod_{i}(h_{i}^{e}-1)!!h_{i}^{o}!!}{\prod_{i,j}(h^{e}_{i}-h^{o}_{j})}=N^{-\frac{n}{2}}\;d_{h}\;\frac{\chi_{\\{h\\}}(C_{2})}{\chi_{\\{h\\}}(C_{1})}\,,$
(109)

where $n=\sum_{i}{h_{i}}-N(N-1)/2$ is the size of the representation, and
$d_{h}$ is the dimension of the representation $h$, which can be evaluated as

$d_{h}=\chi_{\\{h\\}}(\mathbb{1})=\frac{\Delta(h)}{\prod_{k=1}^{N-1}k!}\,.$
(110)

We also used

$\frac{\prod_{i}(h_{i}^{e}-1)!!h_{i}^{o}!!}{\prod_{i,j}(h^{e}_{i}-h^{o}_{j})}=\frac{\chi_{\\{h\\}}(C_{2})}{\chi_{\\{h\\}}(C_{1})}\,,$
(111)

as presented in more detail in the appendix of [78].

######  Lemma IX.1.

If there exists a $GL(N)$ matrix $M$ such that Di Francesco-Itzykson integral
(106) can be computed as the character of $M$;

$\langle\chi_{\\{h\\}}(A)\rangle_{0}=\chi_{\\{h\\}}(M)\,,$ (112)

then the integral of our interest (61) can be solved as:

$\langle\chi_{\\{h\\}}(A^{2})\rangle_{0}=\chi_{\\{h\\}}(M^{2})\;.$ (113)

###### Proof.

The function $\chi_{\\{h\\}}(A^{2})$ is a class function of $A$, therefore
there exists an expansion

$\chi_{\\{h\\}}(A^{2})=\sum_{\\{h^{\prime}\\}}c_{\\{h\\},\\{h^{\prime}\\}}\chi_{\\{h^{\prime}\\}}(A)\;.$
(114)

This expansion is finite, since it is the same as the present in the identity
$\chi_{\\{h\\}}(A^{2})=\chi_{\mathrm{Sym}^{2}\\{h\\}}(A)-\chi_{\mathrm{Alt}^{2}\\{h\\}}(A)$.
Thus,

$\begin{split}\langle\chi_{\\{h\\}}(A^{2})\rangle_{0}&=\frac{1}{Z_{0}}\int
dA\;\chi_{\\{h\\}}(A^{2})\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\;\;=\frac{1}{Z_{0}}\int
dA\;\sum_{\\{h^{\prime}\\}}c_{\\{h\\},\\{h^{\prime}\\}}\chi_{\\{h^{\prime}\\}}(A)\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\;\\\
&=\frac{1}{Z_{0}}\sum_{\\{h^{\prime}\\}}c_{\\{h\\},\\{h^{\prime}\\}}\int
dA\;\chi_{\\{h^{\prime}\\}}(A)\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}=\frac{1}{Z_{0}}\sum_{\\{h^{\prime}\\}}c_{\\{h\\},\\{h^{\prime}\\}}I^{(1)}_{\\{h^{\prime}\\}}(\mathbb{1})\\\
&=\sum_{\\{h^{\prime}\\}}c_{\\{h\\},\\{h^{\prime}\\}}\chi_{\\{h^{\prime}\\}}(M)=\chi_{\\{h\\}}(M^{2})\;.\end{split}$
(115)

We note that the exchange between the summation and the integral is legal,
since the summation is finite. ∎

Lemma IX.1 suggests that we shall study the expression in (109) in order to
see if we can find such a $M$ as in (112). Let us present a more general
normalized version of (106);

${I^{(1)}_{\\{h\\}}(B)=\frac{1}{Z(B)}\int
dA\;\chi_{\\{h\\}}(A)e^{-N\mathrm{Tr}\frac{1}{2}AB^{-1}AB^{-1}}\;}\;,$ (116)

where $Z(B)=\int dA\;e^{-N\mathrm{Tr}\frac{1}{2}AB^{-1}AB^{-1}}$ is a
normalization constant. The generalized normalized Di Francesco-Itzykson
integral (116) can be expressed as a product of characters

$\frac{I^{(1)}_{\\{h\\}}(B)}{Z(B)}=\frac{\chi_{\\{h\\}}(B)\chi_{\\{h\\}}(C_{2})}{\chi_{\\{h\\}}(C_{1})}\,,$
(117)

where we used (111). Remark that on the right hand side of (117), the only
part that depends on the choice of $B$ is $\chi_{\\{h\\}}(B)$.

If we set $B=\mathbb{1}$ in (117) , we indeed recover (109). By setting
$B=C_{1}$ in (117), we find that

$I^{(1)}_{\\{h\\}}(C_{1})=\int
dA\;\chi_{\\{h\\}}(A)\;e^{-N\mathrm{Tr}\frac{1}{2}AC_{1}^{-1}AC_{1}^{-1}}\;=Z(C_{1})\;\chi_{\\{h\\}}(C_{2})\;.$
(118)

Since $(C_{2})^{2}=C_{1}$ in the large $N$ limit, by summing over all
representations as in (114) on both sides of (118) we obtain

$I^{(1)}_{\\{h\\}}(C_{1})\;\sim\;\chi_{\\{h\\}}(C_{1})\,,$ (119)

or performing change of variables,

$I^{(1)}_{\\{h\\}}(C_{1})\;\sim\;\int
dA\;\chi_{\\{h\\}}((C_{1}A)^{2})e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}\;\sim\;\chi_{\\{h\\}}(C_{1})\,,$
(120)

which can also be written as

$\langle\chi_{\\{h\\}}(C_{1}A)^{2}\rangle_{0}\sim\;\chi_{\\{h\\}}(C_{1})\,,$
(121)

where $\sim$ denotes equal up to a constant in $N$, and $\langle\rangle_{0}$
denotes the Gaussian average over $A$ as defined in (14). This result is
nearly the expression we aim at computing in (61), except for the extra
$C_{1}$ in the left hand side.

On this note, let us now investigate further what the condition in (109) tells
us about the matrix $M$ in (112).

######  Proposition IX.2.

Given a representation $r$ of $\mathrm{GL}(N)$ defined through the set of
shifted weights ${h}_{i}$, $i=1,...N,$, consider $\chi_{r}$ the character in
that representation $r$. If $M$ is a $GL(N)$ matrix that satisfies

$\int
dA\;\chi_{\\{h\\}}(A)\;e^{-N\mathrm{Tr}\frac{1}{2}A^{2}}=\chi_{\\{h\\}}(M)\,,$
(122)

then, for a positive integer $q$,

$\frac{1}{N}\mathrm{Tr}(M^{2q})=q^{q}\;.$ (123)

###### Proof.

From (50) we can deduce that

$\mathrm{Tr}(M^{p})=\sum_{k=1}^{N}\frac{\chi_{\\{h^{p,k}\\}}(M)}{\chi_{\\{h\\}}(M)}$
(124)

by using a representation $\\{h\\}$ where $\chi_{\\{h\\}}(M)\neq 0$ and the
modified representations $\\{h^{p,k}\\}$ where the shifted weights are related
to the ones of $\\{h\\}$ by $h_{i}^{p,k}=h_{i}+p\,\delta_{k,i}$, for
$i=1,...,N$.

Using (110) in (109) we can deduce that, for some constant $c_{N}$,

$\chi_{\\{h\\}}(M)=c_{N}\,N^{-\frac{n}{2}}\Delta(\\{h\\}^{e})\Delta(\\{h\\}^{o})\prod_{i}(h_{i}^{e}-1)!!h_{i}^{o}!!\,.$
(125)

We assume $N$ and $p$ even. For $i$ even, we set $h_{i}=p(N-i)/2+1$, therefore
$h_{i}$ is always odd for $ieven$. For $i$ odd, we set
$h_{i}=p{(N-i-1)}/{2}+2$, therefore $h_{i}$ is always even for $i$ odd.
Consequently, for $i$ even, $h_{i}^{p,k}=p{(N-i)}/{2}+1+p\,\delta_{k,i}$,
making $h_{i}^{p,k}$ odd and, for $i$ odd,
$h_{i}^{p,k}=p(N-i-1)/2+2+p\,\delta_{k,i}$, making $h_{i}^{p,k}$ even. Notice
that $\\{h\\}^{e}=\\{a+1\,|\,a\in\\{h\\}^{o}\\}$, thus the Vandermonde
determinants $\Delta(\\{h\\}^{e})$ and $\Delta(\\{h\\}^{o})$ are

$\Delta(\\{h\\}^{o})=\Delta(\\{h\\}^{e})=\prod_{i<j}({h}_{i}-{h}_{j})=\prod_{i<j}\left(p\frac{N-i}{2}+1-p\frac{N-j}{2}-1\right)=\prod_{i<j}p\frac{j-i}{2},$
(126)

where the indices $i$ and $j$ are always even. Because of this, by setting
$i=2m$ and $j=2n$ we find

$\Delta(\\{h\\}^{o})=\Delta(\\{h\\}^{e})=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m<n}(n-m)=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m=1}^{\frac{N}{2}-1{}}m!\;.$
(127)

We then get that

$\Delta(\\{h\\}^{o})\Delta(\\{h\\}^{e})=p^{\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m=1}^{\frac{N}{2}-1}(m!)^{2}\;,$
(128)

which will be necessary for the normalization in (124). The terms involving
double factorials in (125) are:

$\prod_{i}(h_{i}^{e}-1)!!=\prod_{i}h_{i}^{o}!!=\prod_{m=1}^{\frac{N}{2}}\left(p\left(\frac{N}{2}-m\right)+1\right)!!\;.$
(129)

The size of the representation is

$n=\sum_{i}h_{i}-\frac{N(N-1)}{2}=\frac{1}{2}\frac{N}{2}\left(p(N-2)+6\right)-\frac{N(N-1)}{2}\;.$
(130)

This way the character of $M$ in the representation $\\{h\\}$ is

$\chi_{\\{h\\}}(M)=c_{N}\;N^{-\frac{n}{2}}\;p^{\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m=1}^{\frac{N}{2}-1}(m!)^{2}\prod_{m=1}^{\frac{N}{2}}\left(p\left(\frac{N}{2}-m\right)+1\right)!!^{2}\,.$
(131)

Now, for the representations $\\{{h^{p,k}}\\}$, set
$\\{{h^{p,k}}\\}^{e}=\\{a\in\\{h^{p,k}\\}|\,a\;\mathrm{is\;even}\\}$ and
$\\{{h^{p,k}}\\}^{o}=\\{a\in\\{h^{p,k}\\}|\,a\;\mathrm{is\;odd}\\}$. For $k$
odd, $\\{{h^{p,k}}\\}^{e}=\\{{h}\\}^{e}$, then

$\Delta(\\{{h^{p,k}}\\}^{e})=\Delta(\\{{h}\\}^{e})=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m=1}^{\frac{N}{2}-1{}}m!\;.$
(132)

For $k$ even, $\\{{h^{p,k}}\\}^{o}=\\{{h}\\}^{o}$, then

$\Delta(\\{{h^{p,k}}\\}^{o})=\Delta(\\{{h}\\}^{o})=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m=1}^{\frac{N}{2}-1{}}m!\;.$
(133)

Another consequence for $k$ even is

$\Delta(\\{{h^{p,k}}\\}^{e})=\prod_{\begin{subarray}{c}i<j\\\
i,j\;\mathrm{even}\end{subarray}}({h^{p,k}}_{i}-{h^{p,k}}_{j})=\prod_{\begin{subarray}{c}i<j\\\
i,j\;\mathrm{even}\end{subarray}}(p\frac{N-i}{2}+2+p\,\delta_{k,i}-p\frac{N-j}{2}-2-p\,\delta_{k,j})=\prod_{\begin{subarray}{c}i<j\\\
i,j\;\mathrm{even}\end{subarray}}\left(p\frac{j-i}{2}+p(\delta_{k,i}-\delta_{k,j})\right)\,.$
(134)

Setting $i=2m$, $j=2n$, and $k=2r$ we find

$\Delta(\\{{h^{p,k}}\\}^{e})=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m<n}\left(n-m+(\delta_{r,m}-\delta_{r,n})\right)\;.$
(135)

Let us emphasize that $n>1$ since $n>m\geq 1$. For $k=2$,

$\Delta(\\{{h^{p,2}}\\}^{e})=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\prod_{m<n}\left(n-m+\delta_{1,m}\right)=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\frac{N}{2}\prod_{i=1}^{\frac{N}{2}-1}i!\;.$
(136)

For $k>2$, $k$ even, since the factor for $m=r-1$ and $n=r$ in the product is
zero,

$\Delta(\\{{h^{p,k}}\\}^{e})=0\;.$ (137)

It is also true that
$\\{{h^{p,k}}\\}^{o}=\\{a-1|\;a\in\\{{h^{p,k+1}}\\}^{e}\\}$, then
$\Delta(\\{{h^{p,1}}\\}^{o})=\Delta(\\{{h^{p,2}}\\}^{e})$, therefore

$\Delta(\\{{h^{p,1}}\\}^{o})=p^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}\frac{N}{2}\prod_{i=1}^{\frac{N}{2}-1}i!\,.$
(138)

For $k$ odd with $k>2$, since
$\Delta(\\{{h^{p,k}}\\}^{o})=\Delta(\\{{h^{p,k+1}}\\}^{e})$,

$\Delta(\\{{h^{p,k}}\\}^{o})=0\;.$ (139)

Using (132), (133), (136), (137), (138) and (139) we get

$\Delta(\\{{h^{p,k}}\\}^{e})\Delta(\\{{h^{p,k}}\\}^{o})=\begin{cases}p^{\frac{N}{2}\left(\frac{N}{2}-1\right)}\frac{N}{2}\prod_{i=m}^{\frac{N}{2}-1}(m!)^{2}&\text{
if }k=1,2\\\ 0&\text{ if }k>2\end{cases}\;.$ (140)

The consequence of (140) in the sum in (124) is that only the terms for $k=1$
and $k=2$ are nonzero. For $k=1$ or $k=2$, the double factorial terms in (125)
is

$\prod_{i}(h_{i}^{e}-1)!!h_{i}^{o}!!=\frac{(p\frac{N}{2}+1)!!}{\left(p\left(\frac{N}{2}-1\right)+1\right)!!}\prod_{m=1}^{\frac{N}{2}}\left(p\left(\frac{N}{2}-m\right)+1\right)!!^{2}\;.$
(141)

The size of the representation is

$n^{p,k}=\sum_{i}h^{p,k}_{i}-\frac{N(N-1)}{2}=p+\frac{1}{2}\frac{N}{2}\left(p(N-2)+6\right)-\frac{N(N-1)}{2}\;.$
(142)

This way the character of $M$ in the representation $\\{{h^{p,k}}\\}$ is, for
$k=1$ or $k=2$,

$\chi_{\\{{h^{p,k}}\\}}(M)=c_{N}\,N^{-\frac{n^{p,k}}{2}}\;p^{\frac{N}{2}\left(\frac{N}{2}-1\right)}\frac{N}{2}\prod_{m=1}^{\frac{N}{2}-1}(m!)^{2}\frac{(p\frac{N}{2}+1)!!}{\left(p\left(\frac{N}{2}-1\right)+1\right)!!}\prod_{m=1}^{\frac{N}{2}}\left(p\left(\frac{N}{2}-m\right)+1\right)!!^{2}\,,$
(143)

and zero for $k\geq 3$. Hence, inserting (131) and (143) in (124), we find
that

$\mathrm{Tr}(M^{p})=N^{1-\frac{p}{2}}\frac{\left(p\frac{N}{2}+1\right)!!}{\left(p\left(\frac{N}{2}-1\right)+1\right)!!}\;.$
(144)

We evaluate a large $N$ limit by keeping only the largest order in $N$,

$\mathrm{Tr}(M^{p})=N^{1-\frac{p}{2}}\prod_{j=1}^{p/2}\left(p\left(\frac{N}{2}-1\right)+1+2j\right)=N\left(\frac{p}{2}\right)^{\frac{p}{2}}\left(1+{\mathcal{O}}(N^{-1})\right)\;.$
(145)

Setting the integer $q=p/2$, we get the trace property

$\frac{1}{N}\mathrm{Tr}(M^{2q})=q^{q}\;.$ (146)

∎

In principle, this information should be enough to find $M$ up to matrix
conjugation, similarly to what is done in (30) for the matrix $C_{m}$. We
leave this possibility of computation of $M$ to future studies.

Let us now evaluate $\langle\chi_{r}(A^{2})\rangle_{0}$ for a finite $N$. In
particular, in Theorem IX.3, the integral over matrices has been evaluated and
replaced with a summation over integers.

######  Theorem IX.3.

Let $A$ be a random variable for a $N\times N$ Hermitian matrix under the
Gaussian measure. Given a representation $r$ of $\mathrm{GL}(N)$ defined
through the set of shifted weights $h_{i}$, $i=1,...N,$ and considering
$\chi_{r}$ the character in that representation $r$ the following holds true:

$\langle\chi_{r}(A^{2})\rangle_{0}=\frac{N^{\frac{N(N-1)}{2}}}{\prod_{k=0}^{N-1}k!}\frac{\prod_{i}(2h_{i})!}{(2N)^{\sum_{i}h_{i}}}\underset{i,j}{\mathrm{Pf}}\sum_{\begin{subarray}{c}k+l=2h_{i}\\\
u+v=2h_{j}\end{subarray}}\frac{(-1)^{u}-(-1)^{k}}{2}\frac{(k+u)!!(l+v-2)!!}{k!u!l!v!}\,.$
(147)

###### Proof.

For a Hermitian matrix $A$ of size $N$ whose eigenvalues are denoted by $x$,
we write $X$ as a diagonal matrix whose diagonal elements are the eigenvalues
$x$’s. For a given representation $r$, we are interested in:

$\displaystyle\langle\chi_{r}(A^{2})\rangle_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{0}}\int dA\,\chi_{r}(A^{2})\;e^{-\frac{N}{2}{\rm
Tr}A^{2}}$ (148) $\displaystyle=$
$\displaystyle\frac{{\tilde{c}}_{N}}{N!}\int_{\mathbb{R}^{N}}dX\prod_{i<j}(x_{i}-x_{j})^{2}\frac{\underset{i,j}{\rm
det}\,\left(x_{i}^{2h_{j}}\right)}{\prod(x_{i}^{2}-x_{j}^{2})}e^{-\frac{N}{2}{\rm
Tr}X^{2}}$ (149) $\displaystyle=$
$\displaystyle\frac{{\tilde{c}}_{N}}{N!}\int_{\mathbb{R}^{N}}dX\,\underset{i,j}{\rm
det}\,\Big{(}x_{i}^{2h_{j}}\Big{)}\prod_{i<j}\frac{x_{i}-x_{j}}{x_{i}+x_{j}}e^{-\frac{N}{2}{\rm
Tr}X^{2}}\,,$ (150)

where

${\tilde{c}}_{N}=\frac{1}{Z_{0}}\frac{{\rm
vol}(U(N))}{(2\pi)^{N}}=N^{\frac{N^{2}}{2}}(2\pi)^{-\frac{N}{2}}(\prod_{k=0}^{N-1}k!)^{-1}.$
(151)

See [120], [121], and [122] for the change of variables from $A$ to $X$. We
use de Bruijn’s formula [123],

$\frac{1}{N!}\int_{\mathbb{R}^{N}}d\mu(X)\;\underset{i,j}{\mathrm{det}}f_{i}(x_{j})\;\underset{i,j}{\mathrm{Pf}}A(x_{i},x_{j})=\underset{i,j}{\mathrm{Pf}}\int_{\mathbb{R}^{2}}d\mu(x)d\mu(y)\;f_{i}(x)\;A(x,y)\;f_{j}(y)\,,$
(152)

so we can reduce the integration over the N variables on the left hand side to
a Pfaffian of an integral over two variables on the right hand side. $d\mu(x)$
sets the measure on $x$. Here, it is a Gaussian measure
$d\mu(x)=e^{-\frac{N}{2}x^{2}}$. Since [124]

$\prod_{i<j}\frac{x_{i}-x_{j}}{x_{i}+x_{j}}=\underset{i,j}{\mathrm{Pf}}\Bigg{(}\frac{x_{i}-x_{j}}{x_{i}+x_{j}}\Bigg{)}\,,$
(153)

let us first rewrite (150),

$\langle\chi_{r}(A^{2})\rangle_{0}=\frac{{\tilde{c}}_{N}}{N!}\int_{\mathbb{R}^{N}}dX\,\underset{i,j}{\rm
det}\big{(}x_{i}^{2h_{j}}\big{)}\;\underset{i,j}{\mathrm{Pf}}\Bigg{(}\frac{x_{i}-x_{j}}{x_{i}+x_{j}}\Bigg{)}\;e^{-\frac{N}{2}{\rm
Tr}X^{2}}\;.$ (154)

An important thing to notice is that when $x_{i}+x_{j}=0$, the Pfaffian has a
divergence that is controlled by the zero in the determinant. However, when we
use the de Bruijn’s formula the determinant is removed. Therefore, it is best
to deal with this divergence already here. We introduce a damping with a small
constant $\epsilon$ to prevent a divergence from appearing. We take care of
this by multiplying the elements in the Pfaffian by the term
$\frac{(x_{i}+x_{j})^{2}}{(x_{i}+x_{j})^{2}+\epsilon^{2}}$. We regularize
(154) by defining, for $\epsilon>0$,

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}=\frac{{\tilde{c}}_{N}}{N!}\int_{\mathbb{R}^{N}}dX\,\underset{i,j}{\rm
det}\big{(}x_{i}^{2h_{j}}\big{)}\;\underset{i,j}{\mathrm{Pf}}\Bigg{(}\frac{x_{i}^{2}-x_{j}^{2}}{(x_{i}+x_{j})^{2}+\epsilon^{2}}\Bigg{)}\;e^{-\frac{N}{2}{\rm
Tr}X^{2}}\,.$ (155)

We expand the Pfaffian in (155), a polynomial, into its monomials. We do the
same for the Pfaffian in (154). By comparing the absolute values of these
integrands, term by term, we see that the ones from (155) are bounded by the
ones from (154). Therefore, the dominated convergence theorem [125] tells us
that

$\langle\chi_{r}(A^{2})\rangle_{0}=\lim_{\epsilon\rightarrow
0}\langle\chi_{r}(A^{2})\rangle_{\epsilon}\;.$ (156)

Now we can use the de Bruijn’s formula in (152) and we obtain the damped
average

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}={\tilde{c}}_{N}\;\underset{i,j}{\mathrm{Pf}}\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+\epsilon^{2}}x^{2h_{i}}y^{2h_{j}}\;e^{-\frac{N}{2}(x^{2}+y^{2})}\;.$
(157)

We remark that at this point, if we take the $\epsilon\rightarrow 0$ limit,
the integral in (157) turns into a principal value integral. Let us define

$T_{ij}=\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+\epsilon^{2}}x^{2h_{i}}y^{2h_{j}}e^{-\frac{N}{2}(x^{2}+y^{2})}\;.$
(158)

Then, according to (157),

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}={\tilde{c}}_{N}\underset{i,j}{\mathrm{Pf}}\;T_{ij}\;.$
(159)

We can simplify the $N$ dependence by changing variables $x,y\rightarrow
N^{-\frac{1}{2}}x,N^{-\frac{1}{2}}y$,

$T_{ij}=\frac{1}{N^{h_{i}+h_{j}+1}}\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+N\epsilon^{2}}x^{2h_{i}}y^{2h_{j}}e^{-\frac{1}{2}(x^{2}+y^{2})}\;.$
(160)

Introducing sources $\alpha$ and $\beta$ through the terms $\alpha x$ and
$\beta y$ in the exponential, we can turn the factors $x^{2h_{i}}$ and
$y^{2h_{j}}$ into derivatives:

$T_{ij}=\frac{1}{N^{h_{i}+h_{j}+1}}\left.\frac{d^{2h_{i}}}{d\alpha^{2h_{i}}}\frac{d^{2h_{j}}}{d\beta^{2h_{j}}}\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+N\epsilon^{2}}e^{-\frac{1}{2}(x^{2}+y^{2})+\alpha
x+\beta y}\right|_{\alpha,\beta=0}\;.$ (161)

Let us also define

${\rm
II}(\alpha,\beta)=\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+N\epsilon^{2}}e^{-\frac{1}{2}(x^{2}+y^{2})+\alpha
x+\beta y}\,,$ (162)

thus

$T_{ij}=\frac{1}{N^{h_{i}+h_{j}+1}}\left.\frac{d^{2h_{i}}}{d\alpha^{2h_{i}}}\frac{d^{2h_{j}}}{d\beta^{2h_{j}}}{\rm
II}(\alpha,\beta)\right|_{\alpha,\beta=0}\;.$ (163)

By changing integration variables to $u=(x+y)/\sqrt{2}$ and
$v=(x-y)/\sqrt{2}$, and also defining $a=(\alpha+\beta)/\sqrt{2}$ and
$b=(\alpha-\beta)/\sqrt{2}$, it becomes

${\rm
II}(\alpha,\beta)=\int_{\mathbb{R}^{2}}dudv\;\frac{uv}{v^{2}+2N\epsilon^{2}}e^{-\frac{1}{2}(u^{2}+v^{2})+au+bv}\;.$
(164)

Here we notice that the integration over $u$ and over $v$ are independent,
hence we can separate them,

${\rm
II}(\alpha,\beta)=\int_{\mathbb{R}}du\;u\;e^{-\frac{1}{2}u^{2}+au}\int_{\mathbb{R}}dv\;\frac{v}{v^{2}+2N\epsilon^{2}}e^{-\frac{1}{2}v^{2}+bv}={\rm
II}^{(1)}(a){\rm II}^{(2)}_{N\epsilon}(b)\,,$ (165)

where we define the integrals

${\rm
II}^{(1)}(a)=\int_{\mathbb{R}}du\;u\;e^{-\frac{1}{2}u^{2}+au}\qquad\mathrm{and}\qquad{\rm
II}^{(2)}_{\epsilon}(b)=\int_{\mathbb{R}}dv\;\frac{v}{v^{2}+2N\epsilon^{2}}e^{-\frac{1}{2}v^{2}+bv}\;.$
(166)

The first integral is easily evaluated as

${\rm II}^{(1)}(a)=\sqrt{2\pi}\,a\;e^{\frac{1}{2}a^{2}}\,.$ (167)

The second integral we solve by introducing another integral,

${\rm
II}^{(2)}_{\epsilon}(b)=\int_{0}^{b}d\tilde{b}\int_{\mathbb{R}}dv\;\frac{v^{2}}{v^{2}+2N\epsilon^{2}}e^{-\frac{1}{2}v^{2}+{\tilde{b}}v}\;.$
(168)

At this point we can evaluate the $\epsilon\rightarrow 0$ limit due to the
dominated convergence theorem. Thus, by defining

${\rm II}^{(2)}_{0}(b)=\lim_{\epsilon\rightarrow 0}{\rm
II}^{(2)}_{\epsilon}(b)\,,$ (169)

we find that

${\rm
II}^{(2)}_{0}(b)=\int_{0}^{b}d\tilde{b}\int_{\mathbb{R}}dv\;e^{-\frac{1}{2}v^{2}+{\tilde{b}}v}\;.$
(170)

The integral over $v$ in (170) is a simple Gaussian, and we can evaluate it to
find

${\rm
II}^{(2)}_{0}(b)=\sqrt{2\pi}\int_{0}^{b}d{\tilde{b}}\;e^{\frac{1}{2}{\tilde{b}}^{2}}\;.$
(171)

This function is, up to normalization conventions, the imaginary error
function. Joining (167) and (171) in (165), we deduce that

${\rm
II}(\alpha,\beta)=2\pi\,a\;e^{\frac{1}{2}a^{2}}\int_{0}^{b}d{\tilde{b}}\;e^{\frac{1}{2}{\tilde{b}}^{2}}\;.$
(172)

Going back through (163) and (159), we obtain

$\langle\chi_{r}(A^{2})\rangle_{0}=\left.{\tilde{c}}_{N}\underset{i,j}{\mathrm{Pf}}\frac{2\pi}{N^{h_{i}+h_{j}+1}}\frac{\partial^{2h_{i}}}{\partial\alpha^{2h_{i}}}\frac{\partial^{2h_{j}}}{\partial\beta^{2h_{j}}}\,a\,e^{\frac{1}{2}a^{2}}\int_{0}^{b}d\tilde{b}\;e^{\frac{1}{2}{\tilde{b}}^{2}}\right|_{\alpha,\beta=0}\,,$
(173)

and by using that

$a\,e^{\frac{1}{2}a^{2}}=\sum_{k=0}^{\infty}\frac{a^{2k+1}}{2^{k}k!}\qquad\mathrm{and}\qquad\int_{0}^{b}d\tilde{b}\;e^{\frac{1}{2}{\tilde{b}}^{2}}=\sum_{k=0}^{\infty}\frac{b^{2k+1}}{(2k+1)2^{k}k!}$
(174)

and that $a=(\alpha+\beta)/\sqrt{2}$ and $b=(\alpha-\beta)/\sqrt{2}$, we can
evaluate the derivatives in (173) to find that

$\langle\chi_{r}(A^{2})\rangle_{0}={\tilde{c}}_{N}\underset{i,j}{\mathrm{Pf}}\frac{2\pi}{N^{h_{i}+h_{j}+1}}\frac{(2h_{i})!(2h_{j})!}{2^{h_{i}+h_{j}}}\sum_{\begin{subarray}{c}k+l=2h_{i}\\\
u+v=2h_{j}\\\
k+u\mathrm{\,is\,odd}\end{subarray}}(-1)^{u}\frac{(k+u)!!(l+v-2)!!}{k!u!l!v!}\,,$
(175)

which becomes (147) by using (151), Pfaffian properties and that
$(-1)^{k}=-(-1)^{u}$ for $k+u$ odd and $(-1)^{k}=(-1)^{u}$ for $k+u$ even. ∎

We wish to remark here about the quantity $\langle\chi_{r}(A^{2})\rangle_{0}$
which we computed by taking the limit of $\epsilon\rightarrow 0$ in the
expression $\langle\chi_{r}(A^{2})\rangle_{\epsilon}$ given in (159). For a
given (in other words, finite) $N$, the expression we obtained in (175) is
valid and therefore Theorem IX.3. However, if we pay attention the expression
(160), we notice that $N$ comes with $\epsilon$. Then, one notices that once
we send $N$ to infinity, this procedure becomes sensitive to the ratio in
which $N\rightarrow\infty$ and $\epsilon\rightarrow 0$ are sent.

One naturally wonders if the expression obtained in Theorem IX.3 may become
simpler in large $N$ limit. Let us explore this possibility in Proposition
IX.4.

######  Proposition IX.4.

Let $A$ be a random variable for a $N\times N$ Hermitian matrix under the
Gaussian measure. Given a representation $r$ of $\mathrm{GL}(N)$ defined
through the set of normalized shifted weights ${\tilde{h}}_{i}=h_{i}/N$,
$i=1,...N,$. Consider $\chi_{r}$ the character in that representation $r$.
Defining

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}=\frac{1}{Z_{0}}\frac{{\rm
vol}(U(N))}{N!\,(2\pi)^{N}}\int dX\,\underset{i,j}{\rm
det}\,x_{i}^{2N\tilde{h}_{j}}\underset{i,j}{\mathrm{Pf}}\frac{x_{i}^{2}-x_{j}^{2}}{(x_{i}+x_{j})^{2}+\epsilon^{2}}e^{-\frac{N}{2}{\rm
Tr}X^{2}}\,,$ (176)

which satisfies $\langle\chi_{r}(A^{2})\rangle_{0}=\lim_{\epsilon\rightarrow
0}\langle\chi_{r}(A^{2})\rangle_{\epsilon}$, the following holds true:

$\lim_{N\rightarrow\infty}\frac{\langle\chi_{r}(A^{2})\rangle_{\epsilon}}{{\tilde{c}}_{N}\prod_{k}2e^{-N{\tilde{h}}_{k}}(2{\tilde{h}}_{k})^{N{\tilde{h}}_{k}}\underset{i,j}{\rm
Pf}\left[\frac{({\tilde{h}}_{i}-{\tilde{h}}_{j})({\tilde{h}}_{i}+{\tilde{h}}_{j}+\epsilon^{2}/2)}{({\tilde{h}}_{i}-{\tilde{h}}_{j})^{2}+\epsilon^{2}({\tilde{h}}_{i}+{\tilde{h}}_{j})+\epsilon^{4}/4}\right]}=1\;,$
(177)

where
${\tilde{c}}_{N}=N^{\frac{N^{2}}{2}}(2\pi)^{-\frac{N}{2}}(\prod_{k=0}^{N-1}k!)^{-1}$.

###### Proof.

We apply the saddle point method to compute (157). Let us first rescale
integers $h_{i}$ to ${\tilde{h}}_{i}=h_{i}/N$;

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}={\tilde{c}}_{N}\;\underset{i,j}{\mathrm{Pf}}\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+\epsilon^{2}}x^{2N{\tilde{h}}_{i}}y^{2N{\tilde{h}}_{j}}\;e^{-\frac{N}{2}(x^{2}+y^{2})}$
(178)

and prepare in a form proper to use the saddle point approximation,

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}={\tilde{c}}_{N}\;{\rm
Pf}\int_{\mathbb{R}^{2}}dxdy\;\frac{x^{2}-y^{2}}{(x+y)^{2}+\epsilon^{2}}\;e^{-N(\frac{1}{2}x^{2}+\frac{1}{2}y^{2}-2{\tilde{h}}_{i}{\rm
ln}|x|-2{\tilde{h}}_{j}{\rm ln}|y|)}\,.$ (179)

Laplace’s method of integration [126] can be expressed as

$\int_{\mathbb{R}^{d}}dX\,g(X)\,e^{-Nf(X)}=\sum_{X_{0}}\left(\frac{2\pi}{N}\right)^{d/2}\frac{g(X_{0})e^{-Nf(X_{0})}}{\sqrt{{\rm
det}\big{(}H(f)(X_{0})\big{)}}}(1+{\mathcal{O}}(N^{-1}))\,,$ (180)

where $X$ is a set of $d$ real variables, $f$ is a twice-differentiable
complex valued function of $X$, $H(f)$ is the Hessian matrix of $f$, the
points $X_{0}$ are local maxima of $f$, and $g$ is a complex valued function
of $X$ nonzero at $X_{0}$. Comparing (179) and (180) we identify

$f(X)=\frac{1}{2}x^{2}+\frac{1}{2}y^{2}-2{\tilde{h}}_{i}{\rm
ln}|x|-2{\tilde{h}}_{j}{\rm ln}|y|\quad\mathrm{and}\quad
g(X)=\frac{x^{2}-y^{2}}{(x+y)^{2}+\epsilon^{2}}\;.$ (181)

Computing the saddle point equations $x-2{\tilde{h}}_{i}x^{-1}=0$ and
$y-2{\tilde{h}}_{j}y^{-1}=0$, we find four saddle points:

$x=\pm\sqrt{2{\tilde{h}}_{i}}\quad\mathrm{with}\quad
y=\pm\sqrt{2{\tilde{h}}_{j}}\;.$ (182)

Therefore, we see that the term associated with $f$ is the same for any
$X_{0}$ and is

$e^{-Nf(X_{0})}=(2{\tilde{h}}_{i})^{N{\tilde{h}}_{i}}(2{\tilde{h}}_{j})^{N{\tilde{h}}_{j}}e^{-N({\tilde{h}}_{i}+{\tilde{h}}_{j})}\,,$
(183)

and therefore, for the Hessian,

$H(f)(X_{0})=\begin{pmatrix}2&0\\\
0&2\end{pmatrix}\;\quad{\rm{and}}\quad\sqrt{{\rm
det}\big{(}H(f)(X_{0})\big{)}}=2\;.$ (184)

Additionally,

$\sum_{X_{0}}g(X_{0})=2\frac{2{\tilde{h}}_{i}-2{\tilde{h}}_{j}}{(\sqrt{2{\tilde{h}}_{i}}+\sqrt{2{\tilde{h}}_{j}})^{2}+\epsilon^{2}}+2\frac{2{\tilde{h}}_{i}-2{\tilde{h}}_{j}}{(\sqrt{2{\tilde{h}}_{i}}-\sqrt{2{\tilde{h}}_{j}})^{2}+\epsilon^{2}}=4\frac{({\tilde{h}}_{i}-{\tilde{h}}_{j})({\tilde{h}}_{i}+{\tilde{h}}_{j}+\epsilon^{2}/2)}{({\tilde{h}}_{i}-{\tilde{h}}_{j})^{2}+\epsilon^{2}({\tilde{h}}_{i}+{\tilde{h}}_{j})+\epsilon^{4}/4}\,.$
(185)

Joining everything, we obtain the saddle point approximate for the regularized
character of $A^{2}$ as

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}={\tilde{c}}_{N}\underset{i,j}{\rm
Pf}\left[\frac{1}{2}e^{-N({\tilde{h}}_{i}+{\tilde{h}}_{j})}(2{\tilde{h}}_{i})^{N{\tilde{h}}_{i}}(2{\tilde{h}}_{j})^{N{\tilde{h}}_{j}}4\frac{({\tilde{h}}_{i}-{\tilde{h}}_{j})({\tilde{h}}_{i}+{\tilde{h}}_{j}+\epsilon^{2}/2)}{({\tilde{h}}_{i}-{\tilde{h}}_{j})^{2}+\epsilon^{2}({\tilde{h}}_{i}+{\tilde{h}}_{j})+\epsilon^{4}/4}\right](1+{\mathcal{O}}(N^{-1}))\;.$
(186)

Using some Pfaffian properties we can simplify the expression (186) and get

$\langle\chi_{r}(A^{2})\rangle_{\epsilon}={\tilde{c}}_{N}2^{\frac{N}{2}}\prod_{k}e^{-N{\tilde{h}}_{k}}(2{\tilde{h}}_{k})^{N{\tilde{h}}_{k}}\underset{i,j}{\rm
Pf}\left[\frac{({\tilde{h}}_{i}-{\tilde{h}}_{j})({\tilde{h}}_{i}+{\tilde{h}}_{j}+\epsilon^{2}/2)}{({\tilde{h}}_{i}-{\tilde{h}}_{j})^{2}+\epsilon^{2}({\tilde{h}}_{i}+{\tilde{h}}_{j})+\epsilon^{4}/4}\right](1+{\mathcal{O}}(N^{-1}))\;.$
(187)

Applying the limit $N\rightarrow\infty$, we find

$\lim_{N\rightarrow\infty}\frac{\langle\chi_{r}(A^{2})\rangle_{\epsilon}}{{\tilde{c}}_{N}2^{\frac{N}{2}}\prod_{k}e^{-N{\tilde{h}}_{k}}(2{\tilde{h}}_{k})^{N{\tilde{h}}_{k}}\underset{i,j}{\rm
Pf}\left[\frac{({\tilde{h}}_{i}-{\tilde{h}}_{j})({\tilde{h}}_{i}+{\tilde{h}}_{j}+\epsilon^{2}/2)}{({\tilde{h}}_{i}-{\tilde{h}}_{j})^{2}+\epsilon^{2}({\tilde{h}}_{i}+{\tilde{h}}_{j})+\epsilon^{4}/4}\right]}=1\;.$
(188)

∎

One may wish to apply $\epsilon\rightarrow 0$ limit to (188) and if the limits
$\epsilon\rightarrow 0$ and $N\rightarrow\infty$ commute, then, (188) can be
manipulated to say that in the large $N$ limit,
${\langle\chi_{r}(A^{2})\rangle_{0}}$ is equal to

${{\tilde{c}}_{N}2^{\frac{N}{2}}\prod_{k}e^{-N{\tilde{h}}_{k}}(2{\tilde{h}}_{k})^{N{\tilde{h}}_{k}}\underset{i,j}{\rm
Pf}\left[\frac{{\tilde{h}}_{i}+{\tilde{h}}_{j}}{{\tilde{h}}_{i}-{\tilde{h}}_{j}}\right]}\;.$
(189)

We computed $\langle\chi_{r}(A^{2})\rangle_{0}$ using (147) and (189) for
$1\leq N\leq 30$ (we should note here that if we use integral form for
$\langle\chi_{r}(A^{2})\rangle_{0}$, then even $N=6$ the computation becomes
slow) using Mathematica for trivial, defining, and determinant
representations. However, for the above values of $N$ that we tested, the
values of the expression (189) do not converge to the values computed using
the expression achieved in (147).

Now, we present first a new way of computing $\langle\chi_{R}(A)\rangle_{0}$
below in Proposition IX.5, whose similar technique is used to compute
$\langle\chi_{R}(A^{2})\rangle_{0}$ in Theorem IX.6.

######  Proposition IX.5.

Let $A$ be an $N\times N$ Hermitian matrix under the Gaussian measure. Given a
representation $R$ of $\mathrm{GL}(N)$ defined through the set of shifted
weights $h_{i}$, $i=1,...N,$ and considering $\chi_{R}$ the character in the
representation $R$, the following holds true:

$\langle\chi_{R}(A)\rangle_{0}=\frac{\chi_{R}(\mathbb{1})\chi_{R}(C_{2})}{\chi_{R}(C_{1})}=(-1)^{\frac{1}{2}\frac{N}{2}\left(\frac{N}{2}-1\right)}N^{-\frac{n}{2}}\frac{\Delta(h)}{\prod_{i=0}^{N-1}i!}\frac{\prod_{i}(h_{i}^{e}-1)!!h_{i}^{o}!!}{\prod_{i,j}(h_{i}^{o}-h_{j}^{e})}\,,$
(190)

where the numbers $h$ are separated in a set of $\lceil N/2\rceil$ even
numbers $h^{e}$ and $\lfloor N/2\rfloor$ odd numbers $h^{o}$. If such a
separation is not possible, then the average is $0$.

###### Proof.

For $A\in\mathrm{GL}(N)$, using character orthogonality for the symmetric
group $S_{n}$, we can write the character of $A$ as

$\displaystyle\chi_{R}(A)$ $\displaystyle=$ $\displaystyle\sum_{r\vdash
n}\delta_{R,r}\chi_{r}(A)$ (191) $\displaystyle=$ $\displaystyle\sum_{r\vdash
n}\Big{(}\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\chi_{r}(\sigma)\Big{)}\chi_{r}(A)\,,$
(192)

where the bar on $\bar{\chi}_{R}(\sigma)$ denotes complex conjugate.
Interchanging the sums and using Schur-Weil duality,

$\sum_{r\vdash n}\chi_{r}(\sigma)\chi_{r}(A)=\mathrm{Tr}(\sigma A^{\otimes
n})\,,$ (193)

we obtain

$\chi_{R}(A)=\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}(\sigma A^{\otimes n})\,.$
(194)

Now, we take the average of the above quantity. Wick’s probability theorem
tells us that

$\langle A^{\otimes
n}\rangle_{0}=N^{-\frac{n}{2}}\sum_{\gamma\in[2^{\frac{n}{2}}]}\gamma\,,$
(195)

where $[2^{\frac{n}{2}}]$ is the conjugacy class of permutations in $S_{n}$
with $n/2$ 2-cycles. Therefore we obtain

$\displaystyle\langle\chi_{R}(A)\rangle_{0}$ $\displaystyle=$
$\displaystyle\sum_{\sigma\in
S_{n}}\sum_{\gamma\in[2^{\frac{n}{2}}]}N^{-\frac{n}{2}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}(\sigma\gamma)$
(196) $\displaystyle=$ $\displaystyle\sum_{\sigma\in
S_{n}}\sum_{\gamma\in[2^{\frac{n}{2}}]}\sum_{r\vdash
n}N^{-\frac{n}{2}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\chi_{r}(\sigma\gamma)\chi_{r}(\mathbb{1})\,,$
(197)

where we used (193) in the last equality, with $A={\mathbb{1}}$. Again using
orthogonality relations, rewrite

$\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\chi_{r}(\sigma\gamma)=\delta_{R,r}\frac{1}{s_{R}}\chi_{R}(\gamma)\,,$
(198)

where we denote the dimension of the $S_{n}$ representation
$s_{R}=\chi_{R}({\rm id})$, where ${\rm id}$ is the identity permutation, and

$\displaystyle\langle\chi_{R}(A)\rangle_{0}$ $\displaystyle=$
$\displaystyle\sum_{\gamma\in[2^{\frac{n}{2}}]}\sum_{r\vdash
n}N^{-\frac{n}{2}}\delta_{R,r}\frac{1}{s_{R}}\chi_{R}(\gamma)\;d_{r}$ (199)
$\displaystyle=$ $\displaystyle
N^{-\frac{n}{2}}\frac{d_{R}}{s_{R}}\sum_{\gamma\in[2^{\frac{n}{2}}]}\chi_{R}(\gamma)\,.$
(200)

In order to perform $\sum_{\gamma\in[2^{\frac{n}{2}}]}$, we use the following
trick. Let us first return to the relation (194) for some matrix $M$,

$\chi_{R}(M)=\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}(\sigma M^{\otimes n})\;.$
(201)

We observe that (201) sums over all elements in $S_{n}$, whereas (200) sums
over a subset of $S_{n}$. We aim to extract (200) from (201). We choose $M$
such that the summation in $S_{n}$ is restricted to elements of
$[2^{\frac{n}{2}}]$, achieved by $\mathrm{Tr}(\sigma M^{\otimes n})=0$ for
$\sigma$ not in $[2^{\frac{n}{2}}]$ and constant for when $\sigma$ is in
$[2^{\frac{n}{2}}]$. Using the notation $[\prod_{k}k^{c_{k}}]$ for the cycle
$[\sigma]$ which $\sigma$ belongs to, we wish to find $M$ such that

$\mathrm{Tr}(\sigma M^{\otimes
n})=a\;\delta_{[\sigma],[2^{\frac{n}{2}}]}=a\,\delta_{c_{2},\frac{n}{2}}\prod_{k\neq
2}\delta_{c_{k},0}\,.$ (202)

for some constant $a$ (that might depend on $N$ or $n$ but not on $\sigma$).
But also, we can write

$\mathrm{Tr}(\sigma M^{\otimes n})=\prod_{k}\mathrm{Tr}(M^{k})^{c_{k}}\,,$
(203)

leading us to the identification, for $k\neq 2$ ,

$\mathrm{Tr}(M^{k})^{c_{k}}=\delta_{c_{k},0}\,.$ (204)

Hence, the possibility of a nonzero $c_{k}$ with a general permutation
$\sigma$ tells us that

$\mathrm{Tr}(M^{k})=0\;$ (205)

for $k\neq 2$. Then, recalling

$\mathrm{Tr}(C_{2}^{k})=N\delta_{k,2}\,,$ (206)

we conclude $M=C_{2}$ with $a=N^{\frac{n}{2}}$ is a possible solution. Then,
setting $M=C_{2}$ in (201) with (203), we find that

$\chi_{R}(C_{2})=\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\prod_{k}\mathrm{Tr}(C_{2}^{k})^{c_{k}}\;.$
(207)

Finally, using (206)

$\chi_{R}(C_{2})=\sum_{\sigma\in[2^{\frac{n}{2}}]}\frac{1}{n!}\bar{\chi}_{R}(\sigma)N^{\frac{n}{2}}\,,$
(208)

we succeeded in restricting the sum to $[2^{\frac{n}{2}}]$. Identifying the
last summation in (200) as

$\sum_{\sigma\in[2^{\frac{n}{2}}]}\chi_{R}(\sigma)=n!N^{-\frac{n}{2}}\bar{\chi}_{R}(C_{2})\,,$
(209)

we achieve an expression for the average of character,

$\langle\chi_{R}(A)\rangle_{0}=n!N^{-n}\frac{d_{R}}{s_{R}}\bar{\chi}_{R}(C_{2})\;.$
(210)

Let us now compute $s_{R}$. We perform a similar trick as above, by
considering $[1^{n}]$ instead of $[2^{\frac{n}{2}}]$. This time there is no
summation in $[1^{n}]$ because the only element in this class is the identity
permutation. We can show in a similar fashion that we find the result through
the matrix $C_{1}$ instead of $C_{2}$. Skipping to the equivalent of (207),
just changing the condition from $k=2$ to $k=1$ in relevant places, we are
able to isolate the character of the identity permutation,

$\chi_{R}(C_{1})=\frac{1}{n!}\bar{\chi}_{R}(id)N^{n}\;.$ (211)

Identifying $s_{R}=\bar{\chi}_{R}(id)={\chi}_{R}(id)$, the dimension of the
$S_{n}$ representation can be written in terms of $\chi_{R}(C_{1})$ by

$s_{R}=n!N^{-n}\bar{\chi}_{R}(C_{1})\;.$ (212)

Finally, using (212) in (210) also with $d_{R}=\chi_{R}(\mathbb{1})$ we find
that

$\langle\chi_{R}(A)\rangle_{0}=\frac{\chi_{R}(\mathbb{1})\chi_{R}(C_{2})}{\chi_{R}(C_{1})}\;.$
(213)

Additionally, with a formula for the character of $C_{m}$ [78]

$\chi_{R}(C_{m})=c\;\;\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)!}\,\mathrm{sgn}\left[\prod_{0\leq\epsilon_{1}<\epsilon_{2}\leq
m-1}\prod_{i,j}(h^{(\epsilon_{2})}_{i}-h^{(\epsilon_{1})}_{j})\right]\,,$
(214)

where
$\\{h^{(\epsilon)}\\}=\\{a\in{\\{h\\}}\;|\;a=\epsilon\;\mathrm{mod}\;m\\}$, we
can also find a different expression for $\langle\chi_{R}(A)\rangle_{0}$, by
inserting (214) into (213)

$\langle\chi_{R}(A)\rangle_{0}=(-1)^{\frac{1}{2}\lceil\frac{N}{2}\rceil(\lceil\frac{N}{2}\rceil-1)}N^{-\frac{n}{2}}\frac{\Delta(h)}{\prod_{i=0}^{N-1}i!}\frac{\prod_{i}(h^{(e)}_{j}-1)!!h^{(o)}!!}{\prod_{i,j}(h^{(o)}_{i}-h^{(e)}_{j})}\;.$
(215)

This expression is written in terms of the highest weights $\\{h\\}$. ∎

Following a similar argument as Proposition IX.5, we compute
$\langle\chi_{r}(A^{2})\rangle_{0}$.

######  Theorem IX.6.

Let $A$ be an $N\times N$ Hermitian matrix under the Gaussian measure. Given a
representation $r$ of $\mathrm{GL}(N)$ defined through the set of shifted
weights $h_{i}$, $i=1,...N,$ and considering $\chi_{r}$ the character in that
representation, the following holds true:

$\langle\chi_{r}(A^{2})\rangle_{0}=\frac{\chi_{r}(\mathbb{1})^{2}}{\chi_{r}(C_{1})}\;\big{(}1+{\mathcal{O}}(N^{-2})\big{)}=N^{-n}\prod_{i}h_{i}!\,\Delta(h)\;\big{(}1+{\mathcal{O}}(N^{-2})\big{)}\;.$
(216)

###### Proof.

We start from the relation

$\chi_{r}(A^{2})=\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}(\sigma(A^{2})^{\otimes
n})\;.$ (217)

We wish to express in terms of $A$ instead of $A^{2}$. We use the relation

$\mathrm{Tr}\big{(}\sigma(A^{2})^{\otimes
n}\big{)}=\mathrm{Tr}\big{(}(\sigma\otimes 1)A^{\otimes 2n}\alpha\big{)}$
(218)

where $\alpha=(1,1+n)(2,2+n)\cdots(n,2n)$ and on the left the trace considers
the vector space $\mathbb{C}_{N}^{\otimes n}$, while on the right the vector
space considered is $\mathbb{C}_{N}^{\otimes 2n}$. This relation can be shown
by using the definition of trace on both sides. We then just write (217),

$\chi_{r}(A^{2})=\sum_{\sigma\in
S_{n}}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}((\sigma\otimes
1)A^{\otimes 2n}\alpha)\,,$ (219)

and just like before in Theorem IX.5, taking the average and using Wick’s
theorem we obtain

$\langle\chi_{r}(A^{2})\rangle_{0}=\sum_{\sigma\in
S_{n}}\sum_{\gamma\in[2^{n}]}N^{-n}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}((\sigma\otimes
1)\,\gamma\,\alpha)\;.$ (220)

At this point we need a relation that we prove in Appendix Appendix, which is

$\sum_{\gamma\in[2^{n}]}\mathrm{Tr}((\sigma\otimes
1)\,\gamma\,\alpha)=\sum_{\rho\in
S_{n}}\mathrm{Tr}(\rho)\mathrm{Tr}(\sigma\,\rho)\;\big{(}1+{\mathcal{O}}(N^{-2})\big{)}\,,$
(221)

where on the right hand side, the traces act on $\mathbb{C}_{N}^{\otimes n}$.
Using this equation we then get

$\langle\chi_{r}(A^{2})\rangle_{0}=\sum_{\sigma\in S_{n}}\sum_{\rho\in
S_{n}}N^{-n}\frac{1}{n!}\bar{\chi}_{R}(\sigma)\mathrm{Tr}(\rho)\mathrm{Tr}(\sigma\,\rho)\;\big{(}1+{\mathcal{O}}(N^{-2})\big{)}\;.$
(222)

The summation over $\sigma$ can then be evaluated by using the summation over
representations (193) with the identity matrix in place of $A$ and the
summation over permutation (198);

$\langle\chi_{r}(A^{2})\rangle_{0}=\sum_{\rho\in
S_{n}}N^{-n}\mathrm{Tr}(\rho)\frac{d_{R}}{s_{R}}\chi_{r}(\rho)\;\big{(}1+{\mathcal{O}}(N^{-2})\big{)}\;.$
(223)

Repeating for the summation in $\rho$, and using $\chi_{r}(id)=s_{R}$ and
(212), we find

$\langle\chi_{r}(A^{2})\rangle_{0}=n!N^{-n}\frac{d_{R}^{2}}{s_{R}}\;(1+{\mathcal{O}}(N^{-2}))=\frac{\chi_{r}(\mathbb{1})^{2}}{\chi_{r}(C_{1})}\;(1+{\mathcal{O}}(N^{-2}))\;.$
(224)

Finally, using (110) and (214) we find that

$\langle\chi_{r}(A^{2})\rangle_{0}=N^{-n}\prod_{i}h_{i}!\,\Delta(h)\;\big{(}1+{\mathcal{O}}(N^{-2})\big{)}\;.$
(225)

This interesting new result is valid considering $N\gg n$. ∎

##  X Conclusions

A comprehensive understanding of quantum gravity-matter systems is not only
important for phenomenological applications but also to check the consistency
of the the underlying quantum theory. A fundamental theory that encompasses
quantum-gravitational degrees of freedom as well as matter fluctuations
depends on the synergy between the collective quantum fluctuations of the full
theory. A paradigmatic example is QCD that loses asymptotic freedom and, thus,
UV-completeness when the number of fermions is sufficiently large. Likewise,
an asymptotically safe theory of quantum gravity can lose its UV fixed point
for a sufficiently intricate matter content. Therefore, understanding the
impact of matter degrees of freedom to the continuum limit of a candidate
theory of quantum gravity is an essential ingredient in its construction.
Conversely, quantum-gravity fluctuations can affect the dynamics of the
underlying matter fields and trigger non-trivial effects that can leave
imprints in the infrared. Furthermore, the large body of evidence reported by
the CDT community reveals that the implementation of causality constraints can
be central for a suitable continuum limit of the lattice-regularized path
integral of quantum gravity. Thus, in view of what was said above, it is
natural to investigate the continuum limit of CDT coupled to matter degrees of
freedom. Actually, it is an active field of interest to understand the
relation, if any, between the Euclidean and Lorentzian formalisms, see, e.g.,
[127, 128, 129, 130, 131].

In the present work, we have defined a toy model that generated two-
dimensional CDT configurations coupled to Ising model degrees of freedom. This
was proposed under the framework of dually-weighted multi-matrix models.
Matrix models enjoy rich analytical tools and, in many occasions, allow for
the establishment of mathematical rigorous results such as eigenvalue
decomposition techniques to explicitly compute matrix-integrals, large $N$
expansion where the sphere topology dominates at leading order in a Feynman
expansion of the partition function [132, 133, 134, 135, 37], double-scaling
limit in which all topologies are taken into account[136, 137, 41],
topological recursion which lets us recursively solve the loop equations of
matrix models and bridges combinatorial maps with algebraic and enumerative
geometry [138], Weingarten calculus to compute matrix integrals [118], etc.
Matrix models indeed generate the Brownian sphere at criticality and
rigorously proven to be equivalent to 2d Liouville gravity [139, 140, 141,
142, 143, 144]. Equipped with such abundant and rich tools, they therefore
provide with us an enticing platform to further explore interesting
mathematics and physics, with also possibly more physically relevant (higher
dimensional) generalization in the context of quantum gravity; Matrix models
have a higher dimensional generalization, i.e., tensor models [44, 47, 145,
146, 147, 148, 149, 150] where it is of interest to understand how to
incorporate causality.

The initial hope of formulating CDT coupled to the Ising model as a multi-
matrix model was to find an exact/analytical solution for it (in other words
compute exactly the partition function), possibly using new techniques
developed in the last decade since the work by Benedetti-Henson in [76]. For
example, one new interesting result [151] in recent works is the
generalization of the Harish-Chandra–Itzykson–Zuber (HCIZ) integral, where the
generalization considers integrals over tensor powers of $U(N)$. They studied
the integral $\int dU\;e^{t\mathrm{Tr}(AUBU^{*})}$, where $U\in U(N)^{\otimes
D}$, $A$ and $B$ are self-adjoint operators in $(\mathbb{C}^{N})^{\otimes D}$,
and $N$ and $D$ are positive integers. The HCIZ result is recovered by setting
$D=1$. For our purposes, a generalization of the HCIZ integral for a trace
with a different power would have been useful. That is, $\int
dU\;e^{t\mathrm{Tr}(AUBU^{*})^{k}}$, where $U\in U(N)$, $A$ and $B$ are self-
adjoint operators in $\mathbb{C}^{N}$, and $N$ and $k$ are positive integers.
Instead, we made use of Weingarten calculus (in Proposition VIII.1), Schur-
Weyl duality (in Proposition IX.5, Theorem IX.6), theory of symmetric group
algebra (in Proposition IX.5, Theorem IX.6, and Lemma X.1), etc to achieve our
results in this paper.

Even though matrix models enjoy vast amount of analytical results, it seems
that solving exactly the partition functions of the pure CDT-like matrix model
of Benedetti-Henson and the CDT-like matrix model coupled with Ising models
remains quite difficult.

One problem we faced was the necessity to do the character expansion more than
once for the CDT-like matrix model with Ising model. See (67) and (68). In
Proposition VIII.1, the appearance of multiple representations from the
expansion required us to use the theory of decomposition of reducible
representations. The current state of this theory is not sufficient for our
purposes, since it lacks closed formulas in terms of the representations
involved; the existence of algorithms is not enough. For example, given a
representation $R$ in terms of a set of shifted highest weights $\\{h\\}$, the
decomposition of the representation $R\otimes R$ is doable through the
Littlewood-Richardson rule. However, this process is not written in terms of
functions which we could apply the standard theory of calculus (e.g.,
derivatives and integrals) to; for example, closed formulas for the
Littlewood-Richardson rule, as well as for the Clebsch-Gordan coefficients are
not known.

Indeed, we foresee a similar problem when we consider solving for
$\langle\chi_{\\{h\\}}(A^{2})\rangle$, by expressing the character in terms of
the symmetric square and alternating representations (which are, in general,
reducible). There exists a relation
$\\{h\\}\otimes\\{h\\}=\mathrm{Sym}^{2}\\{h\\}\oplus\mathrm{Alt}^{2}\\{h\\}$.
This expression tells us that
$\chi_{\\{h\\}}(A)^{2}=\chi_{\mathrm{Sym}^{2}\\{h\\}}(A)+\chi_{\mathrm{Alt}^{2}\\{h\\}}(A)$.
Another known formula is
$\chi_{\\{h\\}}(A^{2})=\chi_{\mathrm{Sym}^{2}\\{h\\}}(A)-\chi_{\mathrm{Alt}^{2}\\{h\\}}(A)$.
This last decomposition can be done with the Carré–Leclerc domino tableaux
algorithm, in a similar way the Littlewood-Richardson rule gives an algorithm
to decompose the product of irreducible representations into a direct sum of
irreducible representations. Nevertheless, this is just an algorithm, and does
not give a closed formula to which representations contribute, and it seems
that such formula is not known. Therefore, we seem to encounter the similar
problem as we did in the Proposition VIII.1, where Clebsch-Gordan coefficients
are unknown. In order to solve for the partition functions of the CDT-like
matrix models, we do need general explicit formulas for such coefficients.

Another difficulty lies in the fact that the representations that contribute
most to the character expansion are expected to have the size proportional to
$N^{2}$. More specifically, the number of shifted weights is $N$ and each
shifted weight grows proportionally to $N$. See Theorem IX.6. When applying
the theory of symmetric group algebra, one of the symmetric groups considered
is $S_{n}$, where $n$ is the size of the representation (corresponding to the
number of boxes of Young diagram of a given representation). In principle, the
size of the representation $n$ should be proportional to $N^{2}$ [152] [153].
Therefore, if we want to compute the partition functions (68) and (60), then
we need to take care of the $N$ dependence in $n$. However, in this current
work, our results in Theorem IX.6 and Lemma X.1 relied on us taking large $N$
limit but keeping $n$ finite. Working with $S_{n}$ for large $n$ is
challenging because the dependence of $n$ is involved in expressions like
(219).

Nevertheless, we have achieved several new results concerning the averages of
character of Hermitian matrix or Hermitian matrix squared for a given
representation, which may be of interest in the mathematics community and
alike.

In particular, the results in Theorems IX.4, IX.3 and IX.6 add to the known
properties of the decomposition of the product of representations $r\otimes r$
in terms of irreducible representations.

In Theorem IX.3, the expression of $\langle\chi_{R}(A^{2})\rangle$ in (147) is
for finite $N$ and $n$. This result reduces the average integral in (14) to a
finite sum over integers, therefore, more computable. The next step is either
evaluate this sum, that has only products and ratios of factorials as terms,
or possibly evaluating the Pfaffian without computing the sum.

We also found an expression for $\langle\chi_{R}(A^{2})\rangle$ at large $N$
and $n$ in Proposition IX.4. The expression itself of the Pfaffian ${\rm
Pf}\left[\frac{{\tilde{h}}_{i}+{\tilde{h}}_{j}}{{\tilde{h}}_{i}-{\tilde{h}}_{j}}\right]$
is similar to the known result in (153), except that in (177) each matrix
element is inverse to the respective element in (153) and is unknown to our
knowledge.

In Theorem IX.6, we explicitly found the leading order expression of the
character of $A^{2}$ for a given representation in the large $N$ limit for a
finite $n$. This result is related to cell modules in Brauer algebra [154].

Additionally, the results Theorem VI.1 on the matrix $C_{m}$ (which is
responsible for yielding the causal structure to the Feynman graphs generated)
are new to our knowledge. We explicitly found the distribution of eigenvalues
of $C_{m}$.

For future works, we can consider several possible extensions and additions of
our current work. It should be possible to extend the result of Theorem IX.6
to any powers of $A$ by generalizing the expression used in (221). This
expression, shown in Appendix Appendix, is obtained by combinatorial
evaluations, which we expect can be applied to extend our result to higher
powers of $A$.

In [76], the critical value of the coupling $g$ in (1) was computed. By
solving for the distribution of highest weights, they were able to observe a
critical behavior where the distribution extended over the entire positive
real line. This criticality happened for coupling $g=1/2$, which is the same
as the criticality found in [108]. We may be able to obtain the critical
property of partition functions (1) if we can understand the properties of
Pfaffians ${\rm
Pf}\left[\frac{{\tilde{h}}_{i}+{\tilde{h}}_{j}}{{\tilde{h}}_{i}-{\tilde{h}}_{j}}\right]$
in (177) well enough, even if we do not know the explicit expression for it.

In summary, the present work puts forward a dually-weighted multi-matrix model
that, for a particular choice of weight, corresponds to two-dimensional CDT
coupled to the Ising model. To the best of our knowledge, there is no known
exact solution for CDT coupled to Ising even in two dimensions. Hence, one
possibility to be explored is whether the model here proposed can be well-
suited for different approximations of the CDT-Ising model as, e.g., by means
of Monte-Carlo simulations. A reformulation of the underlying theory by a
different formalism can lead to new insights and this is worth exploring in
the future. Moreover, we have now the possibility to couple more sophisticated
matter models to two-dimensional CDT in the matrix-model frameowork, such as
the Potts model. This will be reported elsewhere. Lastly, the dually-weighted
multi-matrix model can be investigated on its own, i.e., beyond the choice of
the causal constraint as mainly emphasized in this work.

## Acknowledgements

The authors would like to thank Dario Benedetti, Taro Kimura, Benoit Collins,
Sanjaye Ramgoolam, Mark Wildon, Luca Lionni, Paul Wedrich, Liron Speyer, and
Christopher Chung, Dmitri Grinko, Cihan Pazarbasi and Shinobu Hikami for
discussions and valuable inputs. We would also like to thank the Insitut Henri
Poincaré thematic program “Quantum Gravity, Random Geometry and Holography”,
2023, where part of the work was done and many discussions took place. ADP
acknowledges CNPq under the grant PQ-2 (312211/2022-8), FAPERJ under the
“Jovem Cientista do Nosso Estado” program (E26/202.800/2019 and
E-26/205.924/2022), and NWO under the VENI Grant (VI.Veni.192.109) for
financial support.

## Appendix

We show how to derive the relation in $\eqref{app2}$. First we show a more
general relation in the following Lemma X.1, and we obtain (221) by performing