††institutetext: Institut de Mathématiques de Bourgogne, Université Bourgogne
Franche-Comté, France

# Unitary matrix models and random partitions: Universality and multi-
criticality

Taro Kimura and Ali Zahabi

###### Abstract

The generating functions for the gauge theory observables are often
represented in terms of the unitary matrix integrals. In this work, the
perturbative and non-perturbative aspects of the generic multi-critical
unitary matrix models are studied by adopting the integrable operator
formalism, and the multi-critical generalization of the Tracy–Widom
distribution in the context of random partitions. We obtain the universal
results for the multi-critical model in the weak and strong coupling phases.
The free energy of the instanton sector in the weak coupling regime, and the
genus expansion of the free energy in the strong coupling regime are
explicitly computed and the universal multi-critical phase structure of the
model is explored. Finally, we apply our results in concrete examples of
supersymmetric indices of gauge theories in the large $N$ limit.

##  1 Introduction and summary

The construction and counting of the gauge invariant observables, such as BPS
operators in supersymmetric gauge theories, have been central problems in
gauge theory and string theory, for decades. Often, the counting problems in
gauge theory enjoy combinatorial features which makes the problems more
tractable. For example, the counting of single-trace and multi-trace operators
in supersymmetric gauge theories can be seen as the counting of letters and
words in combinatorics, in the plethystic program feng2007counting ;
benvenuti2007counting .

The plethystic program provides a framework to study the generating function
of the gauge theories. An essential step in this framework is applying the
Weyl integration formula in the gauge theory skagerstam1984largen ;
sundborg2000hagedorn ; Aharony:2003sx , to obtain the full generating function
of the multi-trace operators as a group/matrix integral of the plethystic
exponentiation of the single-trace operator generating function. The potential
of this matrix model is turned up to be the double-trace potential, and in the
weak coupling limit, one can approximate the pairwise interaction potentials
between the eigenvalues of the matrix model with a single-trace potential,
i.e. the Gross–Witten–Wadia (GWW) model Gross:1980he ; Wadia:1980cp  and its
generalization to higher order polynomial potentials.

One important aspect of the counting problems is the asymptotic behavior of
the generating functions and their multiplicities. In the context of gauge
theory, to study the thermodynamic aspects such as confinement/deconfinement
phase transition and the perturbative/non-perturbative aspects one need to
develop the asymptotic analysis for the gauge theory indices in the limit of
large parameters and study the singularities of the generating functions. In
the contexts of the AdS/CFT correspondence, the asymptotic aspects of the
gauge theory contain interesting and important information for the gravity
dual theories.

A recent interesting counting problem in gauge theory is the computation of
the index of the $\mathcal{N}=1$ four-dimensional superconformal field
theories (SCFT), the generating function of the BPS operators which are
annihilated by one supercharge. In the context of AdS/CFT correspondence, the
large $N$ and large charge limits of the superconformal indices, as a
candidate of the microscopic explanations of the black holes entropy and phase
transition have attracted a lot of interests and many different methods are
developed toward the asymptotic study of the index, see zaffaroni2020ads  for
a review.

In this work, we introduce an analytic approach, based on the machinery of the
integrable operator formalism in random matrices and random partitions to
study the universal features in the phase structure of gauge theories. The
techniques from random partitions and their asymptotics provide a natural
framework for study of the dynamics of the gauge theory. Having obtained an
alternative formulation of the generic unitary matrix model based on the Schur
partition, i.e. a random partition obeying the Schur measure Borodin:2000IEOT
, in this paper, we apply this machinery in the asymptotic analysis of the
supersymmetric indices and their unitary matrix integrals, and explore the
thermodynamics of the gauge theory and the phase structure. We study the
finite and large $N$ asymptotics and associated phase structure of the unitary
matrix models, mainly the generalized GWW model. Our results include
computation of the partition function of the generic unitary matrix model
which have wide applications in gauge theories. The generating functions for
the supersymmetric indices of the gauge theory such as superconformal index
are often represented in terms of the unitary matrix integrals with double
trace potential. In the limit of weak interactions between the eigenvalues,
they can be approximated by the matrix models with the single-trace potential,
i.e. the generalized Gross–Witten–Wadia model. We aim to shed light on the
implications of the universality for the phase structure of the gauge theory.
In fact, one can imagine that the fluctuations of the random partitions in the
bulk or at the edges explain a possible phase transition in the associated
gauge theories. The generating function of the indices and their matrix
integrals in the large $N$ limit can be represented in terms of the Fredholm
determinants. The asymptotic analysis of the Fredholm determinants and
emergence of the Tracy–Widom (TW) distribution explain the phase structure of
the gauge theories.

Precisely speaking, the edge and the bulk fluctuations in the Schur
partitions, described by the Fredholm determinants with sine, and Airy
kernels, respectively, can be applied in the study of critical dynamics of a
class of gauge theories with a generic single-trace unitary matrix integrals
near its critical point. Moreover, the finer phase structure of the matrix
models known as the multi-critical generalization emerges from the smaller
scale fluctuations and they are studied using the higher Airy kernels and
their asymptotics. In fact, the critical dynamics of the matrix model is
encoded in the asymptotic behavior, i.e. the right and left tails of the TW
distribution. In the large $N$ limit, the sharp phase transition is implicit
in the different behavior of the tails of the distribution. This phase
transition is replaced by a smooth cross-over, in the intermediate domain of
the distribution at finite $N$. Moreover, the perturbative and non-
perturbative aspects of the gauge theory are obtained from the finite $N$
corrections to the asymptotic behavior of the two tails of the distribution.
More precisely, the genus expansion of the free energy is obtained via the
perturbative $1/N$ corrections in the asymptotic analysis of the left tail of
the distribution. The exponential order corrections and related instanton
effects are obtained from the right tail of the distribution, using the Airy
function approximation.

In concrete example of the GWW model and its generalization, our result is
indicating a universal third-order phase transition, and we compare it with
the model-dependent results in the literature which are based on different
plausible methods such as the Coulomb gas method and saddle-point analysis of
the matrix models, and we find a good agreement, when expanding the results
around the critical point. As an application of our result, we study the fine
structure of the Hagedorn phase transition in two concrete examples of the
gauge theory indices. We find that the phase structure of the generalized GWW
model captures the leading order singularities of the Hagedorn phase
transition. Based on this observation, we compute the free energy of these
gauge theories in different regimes. In the finer double-scaling regimes,
similar results for the multi-critical dynamics are obtained.

### Summary

Let us briefly summarize the main results of this paper. We will collect the
main ideas, tools and results which are roughly expressed, and refer to the
bulk of the paper for the exact expressions. By using the character expansion
formula, we can write the matrix integral with multi-critical potential in
terms of the sum over partitions with the Schur measure. In fact, imposing the
constraint on the largest entry of the partition $\lambda_{1}$, the summation
over the random partition is rewritten, modulo a normalization factor, as the
unitary matrix integral Borodin:2000IEOT ,

$\displaystyle\mathcal{Z}_{N}\propto\sum_{\begin{subarray}{c}\lambda\in\mathscr{Y}\\\
\lambda_{1}\leq
N\end{subarray}}s_{\lambda}(\mathsf{X})s_{\lambda}(\mathsf{X})=\int_{\mathrm{U}(N)}\differential{U}\exp\quantity(\sum_{n=1}^{\infty}\frac{1}{n}f(q_{1}^{n},q_{2}^{n},...)\quantity(\tr
U^{n}+\tr U^{-n})),$ (1)

for definitions see Section 3.1. An unrefined version of the partition
function (1) is recently studied in santilli2020exact . Defining the free
energy by $\mathcal{F}=\lim_{N\to\infty}N^{-2}\log\mathcal{Z}_{N}$, we show
that the right/left ($\pm$) edge fluctuation contributes to the free energy,
up to a scaling factor $N^{-2}$, as

$\displaystyle\mathcal{F}\sim\lim_{s\to\pm\infty}\log
F_{p}(s),\qquad\text{where}\qquad
s=\frac{(\beta_{c}-\beta)N}{(\alpha_{p}N)^{\frac{1}{p+1}}},$ (2)

and $F_{p}$ is the higher-order Tracy–Widom distribution Periwal:1990qb ;
Claeys:2009CPAM ; LeDoussal:2018dls ; cafasso2019fredholm , and
$\alpha_{p}(q_{1},q_{2},...)$, $\beta(q_{1},q_{2},...)$ are some model-
dependent parameters and are explicitly expressed in terms of the couplings
kimura2020universal . The main result of our study is that the matrix model
(1), undergoes a multi-critical phase transition at the critical point
$\beta=\beta_{c}$. This phase transitions can be explained from the asymptotic
behavior of the higher-order Tracy–Widom distribution,

$F_{p}(s)\sim\begin{cases}1-\mathcal{O}\left(s^{-\frac{p+1}{p}}e^{-s^{\frac{p+1}{p}}}\right)&\quad
s\rightarrow+\infty\\\ \\\
\mathcal{O}\left(e^{-|s|^{\frac{2(p+1)}{p}}}\right)&\quad
s\rightarrow-\infty\\\ \end{cases}.$ (3)

The above leading asymptotic behavior of the higher-order TW distribution and
the following multi-critical double-scaling parameter,
$s=\alpha_{p}^{-\frac{1}{1+p}}(\beta_{c}-\beta)\ N^{\frac{p}{p+1}},$ imply
that the free energy in two regimes $s\to+\infty$ ($\beta<\beta_{c}$), and
$s\to-\infty$ ($\beta>\beta_{c}$), modulo an additive constant, is given by

$\displaystyle\mathcal{F}\sim\begin{cases}\displaystyle\mathcal{O}(e^{-cN})&\beta<\beta_{c}\\\\[10.00002pt]
\displaystyle\alpha_{p}^{-2/p}|\beta_{c}-\beta|^{2(p+1)/p}+\mathcal{O}(N^{-2})&\beta>\beta_{c}\end{cases}.$
(4)

The above free energy implies a discontinuity in the $(2(p+1)/p)$-th
derivative of the free energy at $\beta=\beta_{c}$ and a multi-critical phase
transition of the order $(2(p+1)/p)$ at this critical point. One can possibly
make a mathematical sense of the fractional derivative in the fractional
calculus. Alternatively, the order of the multi-critical phase transition, can
be replaced by $\lfloor\frac{2(p+1)}{p}\rfloor=3$ for any finite $p>2$. At the
classic case of $p=2$, a universal third-order phase transition is emerged
from the fluctuations at scale $N^{-2/3}$. This can be seen as a universal
generalization of the GWW model Gross:1980he ; Wadia:1980cp , see also
zahabi2016new . For higher $p>2$, the fluctuations at scale
$N^{-\frac{p}{p+1}}$ creates a multi-critical phase transition and as we
increase $p$, the order of the phase transition monotonically decreases. At
infinite $p$, the fluctuations are of order $N^{-1}$ and they causes a second
order phase transition. We focus on the two fixed-points of $p=2$ and
$p=\infty$ and their possible physical implications in the general setting and
in concrete examples.

The rest of the paper is organized as follows. In Section 2, the matrix
integral representation of the supersymmetric indices is explained. In Section
3, the Schur partition and its asymptotics are reviewed. In Sections 4 and 5
the critical and the multi-critical dynamics of the generalized GWW model are
discussed. In Section 6, we study the asymptotics of the double-trace matrix
model in some concrete examples and explicit results about the free energy and
phase structure are obtained.

##  2 Gauge theory indices and matrix models

This Section is a brief review about the two classes of unitary matrix models,
namely the generalized double-trace and generalized GWW model and their
relations. These two classes are related to the partition functions of the
gauge theories on compact manifolds and indices of the supersymmetric gauge
theories.

###  2.1 Counting observables in gauge theory

The observables of the gauge theories are the set of gauge invariant operators
and their correlation functions. We consider four-dimensional gauge theories
and their observables. These are formed of single trace of the products of
operators, called single-trace operators and several single-trace operators
multiplied to each other, called multi-trace operators. The grand-canonical
partition function of the multi-trace operators is obtained from the
plethystic exponentiation of the generating function of the single-trace
operators $f(q_{i})$, where $q_{i}$ is the short notation for $\\{q_{i}\\}$,
the collection of fugacity factors benvenuti2007counting ; feng2007counting ,
and then integrating over the Haar measure of the gauge group to project onto
the gauge-singlets skagerstam1984largen ; sundborg2000hagedorn ;
Aharony:2003sx , known as Weyl integration formula,

$\displaystyle\mathcal{I}(q_{i})=\int_{G}\differential{g}\
\exp\left(\sum_{n=0}^{\infty}\frac{1}{n}f(q_{i}^{n})\chi_{R}(g^{n})\right),$
(5)

where $g\in G$, and in our study $G$ is to be considered as the unitary gauge
group $\mathrm{U}(N)$ or their product, and $\chi_{R}(g)$ is the character in
$R$-representation. In Section 6, we discuss some concrete examples in four-
dimensional $\mathcal{N}=1$ SCFT.

###  2.2 Matrix integral representation

In this Section, the counting problem of the BPS operators and their
generating functions in terms of matrix integral are reviewed. We consider the
group integral that appears as the partition function of $G=\mathrm{U}(N)$
gauge theories with the adjoint representation matter, and by using the
character formula in the adjoint representation in terms of the characters in
the fundamental representation
$\chi_{\text{adj}}(g)=\chi_{F}(g)\chi_{F}^{*}(g)$, the group integral recasts
to a matrix integral,

$\displaystyle\mathcal{I}(q_{i})=\int_{\mathrm{U}(N)}\differential{U}\
\operatorname{PE}\left[f(q_{i};U)\right],$ (6)

where by using the Weyl parameterization, the integral over the Haar measure
can be written as an integral over the eigenvalues $\theta_{i}$ in the
diagonal matrix
$U=\operatorname{diag}\left(e^{\mathrm{i}\theta_{i}}\right)_{i=1}^{N}$,

$\displaystyle\int_{\mathrm{U}(N)}\differential{U}=\frac{1}{N!}\int\prod_{i=1}^{N}\frac{d\theta_{i}}{2\pi}\prod_{1\leq
i<j\leq N}|e^{\mathrm{i}\theta_{i}}-e^{\mathrm{i}\theta_{j}}|^{2},$ (7)

and the plethystic exponential is defined by

$\displaystyle\operatorname{PE}[f(q_{i};U)]=\exp\left[\sum_{n=1}^{\infty}\frac{1}{n}f(q_{i}^{n})(\tr
U^{n}\tr U^{-n})\right].$ (8)

The above matrix integral appears as the generating function for the BPS
multi-trace operators in free gauge theories. In particular, we are going to
consider, in Section 6, two explicit classes of examples of this, which are i)
BPS operators in chiral ring of the free four-dimensional $\mathcal{N}=1$ SCFT
and ii) the BPS operators that are annihilated by one supercharge in
$\mathcal{N}=1$ SCFT, in particular the sixteenth BPS operators in
$\mathcal{N}=4$ supersymmetric Yang–Mills (SYM) theory.

Let us define a closely related model, the generalized GWW matrix model,

$\displaystyle\mathcal{Z}(q_{i})=\int_{\mathrm{U}(N)}\differential{U}\exp\left[N\sum_{n=1}^{\infty}t_{n}(q_{i})\quantity(\tr
U^{n}+\tr U^{-n})\right].$ (9)

Notice that, in general, the $\mathcal{O}(1)$ couplings $t_{1},t_{2},...$ can
be any univariate or multivariate function of the parameters of the model. In
the rest of this paper, we absorb the $N$ dependence into the definitions of
the couplings and thus drop the overall factor $N$. Truncating to the case
$n=1$, $\mathcal{Z}(q_{i})$ reduces to a effective theory of GWW model,

$\displaystyle\mathcal{Z}(q_{i})\sim\mathcal{Z}_{\text{GWW}}=\int_{\mathrm{U}(N)}\differential{U}\exp\left[t_{1}(q_{i})(\tr
U+\tr U^{-1})\right].$ (10)

The gauge theory indices as the double-trace matrix integral (6) is related to
the generalized GWW matrix integral, either effectively at weak coupling in
the large $N$ limit or exactly by the Hubbard-Stratonovich transformation
liu2004fine ; alvarez2006black . In this study, we focus on the former
approximation and leave the latter relation for the future studies. The
double-trace matrix model and the generalized GWW model share equilibrium
conditions. In fact, in the large $N$ limit, $\mathcal{I}(q_{i})$ and
$\mathcal{Z}(q_{i})$ are equivalent, upon identification

$\displaystyle t_{n}(q_{i})=\frac{f(q^{n}_{i})}{n}\langle\tr U^{n}\rangle,$
(11)

and they have a similar phase structure Aharony:2003sx . Moreover, the
truncated case $n=1$ of the double matrix integral (6), and the GWW model, are
dominant in the asymptotic analysis of the generalized models with higher
degrees, $n>1$, see Aharony:2003sx  for more details. Being aware of exact
relation between double trace matrix models and GWW models by the
Hubbard–Stratonovich transformation and saddle point approximation
alvarez2005finite , in this article we are going to focus on the generalized
GWW matrix model and its asymptotic analysis and interpreting the results as
the effective theory for the double trace models at the weak coupling limit.
We hope to return to this problem in future and precisely apply the results of
generalized GWW model to the double-trace model, using the transformation.

Having explained the connection between matrix integral representation of
gauge indices and the generalized GWW model, in the rest of the paper, we are
going to explore the generalized GWW model and its possible implications for
gauge theory indices. We introduce random partitions with Schur measure to
study the phase structure of the generalized GWW model. Our goal is to study
the large $N$ limit of the superconformal index and its phase structure, using
the asymptotic analysis of matrix integral. In the next Section, we use the
asymptotic analysis is performed using the methods of random partitions and
integrable operator formalism. This approach uncovers the universal features
of the phase structure.

##  3 Integrable operator formalism in matrix models

In this part, we review, without details of the proofs, the results obtained
in recent works kimura2020universal ; Betea:2020  about the random partitions
and its asymptotics, and then adopt them to study the phase structure of gauge
theories. The first step in the random partition realization of the gauge
theory is to think of the random partition as discretization of the matrix
integral representation of the partition function of the gauge theory. More
precise relation between the matrix integral and the random partition is via
the character expansion formula.

###  3.1 Schur measure random partition and generalized GWW model

Let $\mathscr{Y}$ be a set of partitions, and
$\mathsf{X}=(\mathsf{x}_{i})_{i\in\mathbb{N}}$ be the set of the parameters
and the Miwa variables defined by

$\displaystyle t_{n}=\frac{1}{n}\sum_{i=1}^{\infty}\mathsf{x}_{i}^{n}.$ (12)

Then, we define our main object, the partition function of the constrained
partitions, via the Schur function $s_{\lambda}$, as

$\displaystyle\mathcal{Z}_{N}=\mathbb{P}(\lambda_{1}<N)=Z^{-1}\sum_{\begin{subarray}{c}\lambda\in\mathscr{Y}\\\
\lambda_{1}\leq
N\end{subarray}}s_{\lambda}(\mathsf{X})s_{\lambda}(\mathsf{X}),$ (13)

where $Z$ is the normalization of the Schur measure random partition and it is
given by,

$\displaystyle
Z=\sum_{\lambda\in\mathscr{Y}}s_{\lambda}(\mathsf{X})s_{\lambda}(\mathsf{X})=\prod_{1\leq
i,j\leq\infty}\quantity(1-\mathsf{x}_{i}\mathsf{x}_{j})^{-1}=\exp\quantity(\sum_{n=1}^{\infty}n\,t_{n}^{2})\,.$
(14)

We will elaborate on the physical interpretation of the partition functions
$Z$ and $\mathcal{Z}$ in Section 3.2. Using the character expansion, the Schur
partition can be written in terms of the unitary matrix integral,

$\displaystyle\mathcal{Z}_{N}=Z^{-1}\int_{\mathrm{U}(N)}\differential{U}\exp\quantity(\sum_{n=1}^{\infty}t_{n}\quantity(\tr
U^{n}+\tr U^{-n})).$ (15)

Let us define the Fredholm determinant on a specific domain
$I\subset\mathbb{R}$,

$\displaystyle\det(1-K)_{I}=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\int_{I^{n}}\prod_{i=1}^{n}\differential{x_{i}}\det_{1\leq
i,j\leq n}K(x_{i},x_{j})\,,$ (16)

where $K$ is the integral operator with the kernel $K(x_{i},y_{j})$. Then, the
Schur partition function can be written as a discrete analog of the original
Fredholm determinant (16) Borodin:2000IEOT ,

$\displaystyle\mathcal{Z}_{N}=\det\left(1-K\right)_{I},$ (17)

where the domain $I$ is defined as $I=(N+\frac{1}{2},\infty)$, and $K(k,l)$ is
given by Okounkov:2001SM ,

$\displaystyle K(k,l)$
$\displaystyle=\frac{1}{(2\pi\mathrm{i})^{2}}\oint_{|z|>|w|}\differential{z}\differential{w}\,\frac{\mathsf{K}(z,w)}{z^{k+1/2}w^{-l+1/2}}\qquad\quantity(k,l\in\mathbb{Z}+\frac{1}{2})\,,$
(18a) $\displaystyle\mathsf{K}(z,w)$
$\displaystyle=\frac{\mathsf{J}(z)}{\mathsf{J}(w)}\,\frac{1}{z-w}\,,\qquad\mathsf{J}(z)=\prod_{n=1}^{\infty}\frac{1-\mathsf{x}_{n}z}{1-\mathsf{x}_{n}/z}\,.$
(18b)

The $\mathsf{J}(z)$ is called the wave function and it has the following mode
expansion,

$\displaystyle\mathsf{J}(z)=\exp\left[\sum_{n=1}^{\infty}t_{n}\quantity(z^{n}-z^{-n})\right]=\sum_{x\in\mathbb{Z}}J(x)\,z^{x}\,,\qquad
J(x)=\oint\frac{\differential{z}}{2\pi\mathrm{i}}\frac{\mathsf{J}(z)}{z^{x+1}}\,.$
(19)

### Plancherel measure random partition and GWW model

A particular example of the Schur measure with a single parameter

$\displaystyle t_{n}=\mathfrak{q}\ \delta_{n,1},$ (20)

is called the Plancherel measure with the following partition function,

$\displaystyle\mathcal{Z}_{N}=Z^{-1}\,\sum_{\begin{subarray}{c}\lambda\in\mathscr{Y}\\\
\lambda_{1}\leq
N\end{subarray}}\mathfrak{q}^{2|\lambda|}\quantity(\frac{\dim\lambda}{|\lambda|})^{2},$
(21)

where $\dim\lambda$ is the dimension of the irreducible representation
parameterized by $\lambda$ of the symmetric group $\mathfrak{S}_{\infty}$, and
$|\lambda|=\sum_{i=1}^{\infty}\lambda_{i}$, and
$Z=\mathrm{e}^{\mathfrak{q}^{2}}$.

The Plancherel partition function with the constraint $\lambda_{1}\leq N$ has
the matrix integral and Fredholm determinant representation,

$\displaystyle\mathcal{Z}_{N}=Z^{-1}\int_{\mathrm{U}(N)}\differential{U}\exp\left(\mathfrak{q}\quantity(\tr
U+\tr
U^{-1})\right)=\det\left(1-K_{\text{dB}}\right)_{(N+\frac{1}{2},\infty)},$
(22)

where the discrete Bessel kernel is obtained in Brodin:2000 ,

$\displaystyle
K_{\text{dB}}(k,l)=\mathfrak{q}\frac{J_{k-1/2}(2\mathfrak{q})\,J_{l+1/2}(2\mathfrak{q})-J_{k+1/2}(2\mathfrak{q})\,J_{l-1/2}(2\mathfrak{q})}{k-l}\,.$
(23)

The wave function is then given by the Bessel function
$J(x)=J_{x}(2\mathfrak{q})$, with the generating function

$\displaystyle\mathsf{J}(z)=\mathrm{e}^{\mathfrak{q}(z-z^{-1})}=\sum_{x\in\mathbb{Z}}J_{x}(2\mathfrak{q})\,z^{x}\,.$
(24)

### Difference equation for wave functions

We observe that the wave functions of the Schur partitions satisfy the
following differential and difference equations:

$\displaystyle\quantity[\sum_{n=1}^{\infty}n\,t_{n}\quantity(z^{n}+z^{-n})-z\partialderivative{z}]\mathsf{J}(z)$
$\displaystyle=0,\qquad\quantity[\sum_{n=1}^{\infty}n\,t_{n}\quantity(\nabla_{x}^{n}+\nabla_{x}^{-n})-x]J(x)=0,$
(25)

where we define the shift operator $\nabla_{x}f(x)=f(x+1)$ with
$\nabla_{x}=\exp\quantity(\partial_{x})$. In the case of the Plancherel
measure, the wave function obeys the difference equation, a special case of
(25),

$\displaystyle\quantity[\nabla_{x}+\nabla_{x}^{-1}-\frac{x}{\mathfrak{q}}]J_{x}(2\mathfrak{q})=0\,.$
(26)

We remark that, as mentioned in kimura2020universal , the
differential/difference equation (25) is interpreted as a quantization of the
spectral curve for the generalized GWW model,

$\displaystyle\Sigma_{\text{GWW}}=\left\\{(x,z)\in\mathbb{C}\times\mathbb{C}^{\times}\mid\mathcal{H}(x,z)=0\right\\}$
(27a) with
$\displaystyle\mathcal{H}(x,z)=\sum_{n=1}^{\infty}n\,t_{n}\quantity(z^{n}+z^{-n})-x\,.$
(27b)

###  3.2 Asymptotic analysis of partitions and dynamics of gauge theory

In this Section, dynamics of gauge theory is obtained from the asymptotic
analysis of the Schur partition. The free energy in the large $N$ limit of the
$\mathrm{U}(N)$ matrix model is defined by

$\displaystyle\mathcal{F}=\lim_{N\to\infty}N^{-2}\log\mathcal{Z}_{N}.$ (28)

To explore the phase structure of the gauge theory, we compute the free energy
in different regimes of the parameter space, obtained from the asymptotic
analysis of the random partitions. In term of the random partition, there are
two contributions to the free energy, first a contribution from the continuum
limit/limit shape of the random partition, obtained from the normalization
factor of the matrix integral. Secondly, there is a contribution from the
fluctuation around the limit shape. Depending on the region in the partition,
the fluctuation behaves differently; in the bulk of the partition, the
contribution is given by the Fredholm determinant of the sine kernel, and in
the edge of the partition, it is the Fredholm determinant with Airy kernel,
i.e. TW distribution. Thus, formally speaking, we have

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+\mathcal{F}_{f},$ (29)

where $\mathcal{F}_{c}$ and $\mathcal{F}_{f}$ denote the continuum and
fluctuation free energy, respectively.

The continuum free energy is a global model dependent contribution and can be
computed using the matrix integral, for example in the generalized GWW model
one can obtain from Eq.(14),

$\displaystyle\mathcal{F}_{c}=\log Z=\sum_{n=1}^{\infty}n\,t_{n}^{2},$ (30)

notice that the $N$ dependence is implicit in the coupling $t_{n}$. The
fluctuation has different behavior in the bulk and edge of the partition,
thus, we consider each case separately.

#### Bulk scaling limit and fluctuation

Now, let us consider the Plancherel measure random partition and discuss the
scaling limit of the Bessel function. Let us define the discreteness parameter
$\epsilon\sim N^{-1}$, to re-scale the parameters of the random partition,
such that the large $N$ limit of the unitary matrix model corresponds to the
continuum limit, $\epsilon\to 0$, of the random partition. Moreover, as we
will explicitly explain later, the parameter $s$ in the domain of the Fredholm
kernel (17), is in fact the double-scaling parameter and as $N$ tends to
infinity, it tends to plus infinity in the bulk and plus/minus infinity at the
right/left sides of the edge of the random partition.

Using the scaling $(x,\mathfrak{q})\to(x,\mathfrak{q}/\epsilon)$, the scaling
limit of the difference equation (26) becomes

$\displaystyle\quantity[\nabla_{x}+\nabla_{x}^{-1}+\mathcal{O}(\epsilon)]J_{n}\quantity(2\mathfrak{q}/\epsilon)=0\,,$
(31)

which has the plane wave solutions

$\displaystyle J_{x}(2\mathfrak{q}/\epsilon)\ \xrightarrow{\epsilon\to 0}\
\begin{cases}\displaystyle\cos\quantity(\frac{\pi}{2}x)=(-1)^{x/2}&(x\in
2\mathbb{Z})\\\\[10.00002pt]
\displaystyle\sin\quantity(\frac{\pi}{2}x)=(-1)^{x/2-1/2}&(x\in
2\mathbb{Z}+1)\end{cases}.$ (32)

Then, the scaling limit of the discrete Bessel kernel (23) becomes the sine
kernel

$\displaystyle K(k,l)\ \xrightarrow{\epsilon\to 0}\
\frac{\sin\pi(k-l)/2}{\pi(k-l)/2}\,,$ (33)

where the normalization is fixed to be $K(r,r)=1$. Similar studies for the
bulk scaling limit of the Schur measure random partition is performed in
okounkov2003correlation .

In the bulk, the dynamics of the gauge theory is governed by the sine random
process and the fluctuation contribution to the free energy is given by the
Fredholm determinant with the sine kernel, up to a scaling factor $N^{-2}$, as

$\displaystyle\mathcal{F}_{f}\sim\lim_{s\to\infty}\log\det\left(1-\hat{K}_{\text{sine}}\right)_{(s,\infty)},\qquad\text{where}\qquad\hat{K}_{\text{sine}}(x,y)=\frac{\sin\pi(x-y)}{\pi(x-y)}$
. (34)

The asymptotics of the Fredholm determinant as $s\to\infty$ is obtained in
deift1997riemann ; widom1994asymptotics ,

$\displaystyle\det\left(1-\hat{K}_{\text{sine}}\right)_{(s,\infty)}=2^{1/12}\,e^{3\,\zeta^{\prime}(-1)}\,s^{-\frac{1}{4}}\exp\left(-\frac{s^{2}}{2}\right)\left(1+\mathcal{O}(s^{-1})\right).$
(35)

This asymptotic result leads to the computation of the free energy in the
bulk, the scaling parameter $s=\gamma N$, as

$\displaystyle\mathcal{F}$ $\displaystyle=$
$\displaystyle\mathcal{F}_{c}+\mathcal{F}_{f}$ (36) $\displaystyle=$
$\displaystyle\mathcal{F}_{c}+\lim_{N\to+\infty}N^{-2}\log\det\left(1-\hat{K}_{\text{sine}}\right)_{(\gamma
N,\infty)}$ $\displaystyle=$
$\displaystyle\mathcal{F}_{c}-\frac{\gamma^{2}}{2}-\frac{1}{N^{2}}\left(-\frac{1}{4}\log
N-\frac{\log\gamma}{4}+\frac{\log 2}{12}+3\,\zeta^{\prime}(-1)\right),$

where $\mathcal{F}_{c}$ is given by the normalization of the Schur measure
random partition in Eq.(30). Thus, we observe that in the large $N$ limit, the
leading contribution of the bulk fluctuation in the free energy is non-zero.
We will come back to this observation and discuss a possible interpretation of
that, at the end of Section 5.2.

#### Edge scaling limit and fluctuation

Next, we discuss the scaling limit of the difference equations of the wave
function and the kernel of the Fredholm determinant. We first consider the
Schur partition. Let us use the parameter $\epsilon$ to re-scale the
parameters and expand the shift operator,

$\displaystyle(x,t_{n})\to(x/\epsilon,t_{n}/\epsilon),\qquad\nabla_{x/\epsilon}^{\pm
n}=\sum_{k=0}^{\infty}\frac{(\pm n\epsilon)^{k}}{k!}\derivative[k]{x}.$ (37)

Then, the difference equation (25) becomes

$\displaystyle\quantity[\sum_{k=1}^{\infty}\alpha_{k}\,\epsilon^{k}\derivative[k]{x}-(x-\beta)]J\quantity(\frac{x}{\epsilon})=0\,,$
(38)

where the coefficients are defined as

$\displaystyle\alpha_{k}=\sum_{n=1}^{\infty}\frac{2n^{k+1}t_{n}}{k!}\,,\qquad\beta=\alpha_{0}=\sum_{n=1}^{\infty}2n\,t_{n}\,.$
(39)

In the scaling limit $\epsilon\to 0$, and by keeping $\alpha_{p^{\prime}}=0$
for $p^{\prime}<p$, the difference equation (38) becomes the following
differential equation,

$\displaystyle\quantity[\derivative[p]{\xi}-\xi]J\quantity(\frac{\beta}{\epsilon}+\quantity(\frac{\alpha_{p}}{\epsilon})^{\frac{1}{p+1}}\xi)=0\,,$
(40)

where $\xi=\alpha_{p}^{-\frac{1}{p+1}}\epsilon^{-\frac{p}{p+1}}(x-\beta)$ for
$p\geq 2$. Thus, we observe that the scaling limit of the wave function is
given by the $p$-Airy function,

$\displaystyle
J\quantity(\frac{\beta}{\epsilon}+\quantity(\frac{\alpha_{p}}{\epsilon})^{\frac{1}{p+1}}\xi)\
\xrightarrow{\epsilon\to 0}\
\text{Ai}_{p}(\xi)=(-1)^{\frac{p}{2}+1}\int_{\mathcal{C}}\frac{\differential\lambda}{2\pi\mathrm{i}}\,\exp{(-1)^{\frac{p}{2}}\frac{\lambda^{p+1}}{p+1}+\xi\lambda},$
(41)

where $\mathcal{C}$ is an integral contour providing a convergent integral.
The case $p=2$ corresponds to the standard Airy function. The $p$-Airy
function, $\text{Ai}_{p}\left((-1)^{\frac{p}{2}+1}\xi\right)$, satisfies the
generalized Airy equation,

$\displaystyle\derivative[p]{\xi}\text{Ai}_{p}(\xi)=\xi\text{Ai}_{p}(\xi).$
(42)

In fact, we can analyze the semi-classical behavior of the $p$-Airy function
based on the reduced spectral curve,

$\displaystyle\Sigma_{p\text{-Airy}}=\left\\{(x,y)\in\mathbb{C}\times\mathbb{C}\mid
y^{p}-x=0\right\\}\,.$ (43)

This is obtained from the spectral curve of the GWW model (27) in the scaling
limit discussed above, i.e., parametrize $w=\exp\quantity(\epsilon y)$, then
take the limit $\epsilon\to 0$ with tuning the parameters. We remark that the
$p$-Airy spectral curve (43) agrees with $(A_{p-1},A_{0})$-type Argyres-
Douglas theory Cordova:2015nma ; Wang:2015mra .

Next, we discuss the scaling limit of the kernel. Using the Eq.(18), we obtain
the kernel in terms of the wave functions,

$\displaystyle\mathsf{K}(z,w)=\sum_{n,m\in\mathbb{Z}}\sum_{i=1}^{\infty}J(n)J(m)\,z^{n-i}w^{-m+i-1},\quad
K(k,l)=\sum_{i=1}^{\infty}J\quantity(k+i-\frac{1}{2})J\quantity(l+i-\frac{1}{2})\,,$
(44)

and then using the scaling limit of the wave function, one can obtain the
scaling limit of the kernel as

$\displaystyle
K\quantity(\frac{\beta}{\epsilon}+\quantity(\frac{\alpha_{p}}{\epsilon})^{\frac{1}{p+1}}x,\,\frac{\beta}{\epsilon}+\quantity(\frac{\alpha_{p}}{\epsilon})^{\frac{1}{p+1}}y)\xrightarrow{\epsilon\to
0}\
\quantity(\frac{\alpha_{p}}{\epsilon})^{\frac{1}{p+1}}K_{p\text{-Airy}}(x,y),$
(45)

where the $p$-Airy kernel is defined as

$\displaystyle K_{p\text{-Airy}}(x,y)$
$\displaystyle=\int_{0}^{\infty}\differential{z}\operatorname{Ai}_{p}(x+z)\operatorname{Ai}_{p}(y+z)$
(46a)
$\displaystyle=\frac{1}{x-y}\sum_{q=0}^{p-1}(-1)^{q}\operatorname{Ai}_{p}^{(q)}(x)\operatorname{Ai}_{p}^{(p-q-1)}(y)\,,$
(46b)

and $\operatorname{Ai}_{p}^{(r)}(x)$ is the $r$-th derivative of the Airy
function.

###  3.3 Higher-order Tracy–Widom distribution

The natural origin of the Airy kernel is in the theory of random matrices, in
which the gap probability of the Airy process is given by the Fredholm
determinant with the Airy kernel,

$\displaystyle F_{2}(s)=\det(1-K_{\text{Airy}})_{(s,\infty)}\,.$ (47)

The probability distribution function of the largest eigenvalue and its
statistical behavior in the scaling limit is given by the Tracy–Widom
distribution Tracy:1992rf .

Random partition, as a discrete analog of the random matrix Borodin:2000IEOT ,
has similar edge scaling behavior which is reflected in terms of the
probability distribution of the largest entry (the first row) of the
partition. In the scaling limit, the probability distribution of the largest
entry of the Schur partition is obtained in kimura2020universal ,

$\displaystyle\lim_{\epsilon\to
0}\mathbb{P}\quantity[\frac{\lambda_{1}-\beta/\epsilon}{(\alpha_{p}/\epsilon)^{\frac{1}{p+1}}}<s]=\det(1-K_{p\text{-Airy}})_{(s,\infty)}=:F_{p}(s)\,,$
(48)

where $F_{p}(s)$ is a higher-order analog of the Tracy–Widom distribution
Periwal:1990qb ; Claeys:2009CPAM ; LeDoussal:2018dls ; cafasso2019fredholm .

The probability distribution (48) can be used to compute the contribution of
the edge fluctuation to the partition function of the generalized GWW model by
using the Eqs. (13) and (17). In fact, there are three regions in the vicinity
of the edge, namely the finitely ($\mathcal{O}(1)$) close regions in the left
and right sides of the edge, and the crossing region which is infinitesimally
close to the edge. At the edge, the dynamics is determined by the $p$-Airy
process and the fluctuations in the right/left ($\pm$) sides of the edge
contribute to the free energy (29), up to a scaling factor $N^{-2}$, as

$\displaystyle\mathcal{F}_{f}\sim\lim_{s\to\pm\infty}\log\det\left(1-K_{p\text{-Airy}}\right)_{(s,\infty)}.$
(49)

Therefore, the right/left ($\pm$) edge free energy are obtained from the left
and right tails of the higher TW distribution,

$\displaystyle\mathcal{F}_{f}\sim\lim_{s\to+\infty}\log
F_{p}(s),\qquad\mathcal{F}_{f}\sim\lim_{s\to-\infty}\log F_{p}(s).$ (50)

In the special case of the Plancherel partition, using the scaled parameters
$(x,\mathfrak{q})\to(x/\epsilon,\mathfrak{q}/\epsilon)$, the difference
equation (26) becomes

$\displaystyle\quantity[\epsilon^{2}\derivative[2]{x}-\quantity(\frac{x}{\mathfrak{q}}-2)+\mathcal{O}(\epsilon^{4})]J_{x/\epsilon}(2\mathfrak{q}/\epsilon)=0\,,$
(51)

and in the scaling limit, the wave function becomes the Airy function,

$\displaystyle J_{x/\epsilon}(2\mathfrak{q}/\epsilon)\
\xrightarrow{\epsilon\to 0}\ \operatorname{Ai}(\xi)\,,$ (52)

where $\xi=\epsilon^{-2/3}\mathfrak{q}^{-1/3}(x-2\mathfrak{q})$. Similarly,
the scaling limit of the discrete Bessel kernel is given by the Airy kernel,

$\displaystyle
K\quantity(\frac{2\mathfrak{q}}{\epsilon}+\quantity(\frac{\mathfrak{q}}{\epsilon})^{1/3}x,\frac{2\mathfrak{q}}{\epsilon}+\quantity(\frac{\mathfrak{q}}{\epsilon})^{1/3}y)\xrightarrow{\epsilon\to
0}\ \quantity(\frac{\mathfrak{q}}{\epsilon})^{1/3}K_{\text{Airy}}(x,y)\,.$
(53)

As a special case of the generic $p$, in the case of $p=2$, the dynamics of
the model is governed by the Airy process and TW distribution, and the
contribution of the fluctuation to the free energy is given, up to a scaling
factor $N^{-2}$, by

$\displaystyle\mathcal{F}_{f}\sim\lim_{s\to\pm\infty}\log\det\left(1-K_{\text{Airy}}\right)_{(s,\infty)}=\lim_{s\to\pm\infty}\log
F_{2}(p).$ (54)

##  4 Critical dynamics

In this Section, we study the free energy of the generic unitary matrix model
from the viewpoint of the fluctuation of the model around the edge and explore
the universal phase structure of the model in the vicinity of a critical
point. As we will see in the following Sections this phase transition is
associated with the opening/closing of a gap in the distribution function of
the eigenvalues on the circle. From the mathematical point of view, the case
$p=2$ explains the critical dynamics of the matrix models.

###  4.1 General unitary matrix model and TW distribution

Precisely speaking, in the particular case $p=2$ of the Eq.(48), and by
assuming the scaling relation between $\epsilon$ and $N$, discussed in Section
3.2, the finite and large $N$ results for the free energy of the critical
model is encoded of the following result,

$\displaystyle\lim_{\epsilon\to
0}\mathbb{P}\quantity[\frac{\lambda_{1}-c_{1}/\epsilon}{(c_{2}/\epsilon)^{1/3}}<s]=\det(1-K_{\text{Airy}})_{(s,\infty)}=F_{2}(s)\,,$
(55)

where $c_{1}$ and $c_{2}$ are some model-dependent parameters which in the
case of generalized GWW model we have $c_{1}=\beta(t_{1},t_{2},...)$ and
$c_{2}=\alpha_{2}(t_{1},t_{2},...)$ in Eq.(39), and $F_{2}$ is the TW
distribution and

$\displaystyle
K_{\text{Airy}}(x,y)=\frac{\text{Ai}(x)\text{Ai}^{\prime}(y)-\text{Ai}(y)\text{Ai}^{\prime}(x)}{x-y}.$
(56)

Using the definition (28), and by fixing the ’t Hooft parameter
$\gamma:=N\epsilon$ in the large $N$ limit of the $\mathrm{U}(N)$ matrix
model, the universal result for the free energy can be obtained from (13),
(29), and (55), up to a scaling factor $N^{-2}$, as

$\displaystyle\mathcal{F}_{f}\sim\lim_{s\to\pm\infty}\log
F_{2}(s)\qquad\text{where}\qquad s=\frac{\gamma-
c_{1}}{c_{2}^{1/3}}N^{\frac{2}{3}},$ (57)

and the free energy in the left side of the edge is given by

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+\lim_{N\to\infty}N^{-2}\log
F_{2}(s),\qquad\text{for}\qquad\gamma<c_{1},\,s\to-\infty,$ (58)

and right sides of the edge is given by

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+\lim_{N\to\infty}N^{-2}\log
F_{2}(s),\qquad\text{for}\qquad\gamma>c_{1},\,s\to+\infty.$ (59)

In the rest of this section, we explain the perturbative and non-perturbative
aspects of the general unitary matrix model in the weak and strong coupling
phases, using the TW distribution. We obtain some results about the finite $N$
and genus expansion of the free energy. In the large $N$ limit, we compute the
free energy and extract the phase structure of the models. Let us start with
the asymptotic analysis of the TW distribution.

#### Tracy–Widom distribution

The Tracy–Widom distribution is given by

$\displaystyle F_{2}(s)=\exp{-\int_{s}^{\infty}(x-s)\ q^{2}(x)\ \differential
x},$ (60)

where $q(x)$ is the solution of the Painlevé II equation,

$\displaystyle q_{xx}(x)=2\,q^{3}(x)+x\,q(x),$ (61)

and $q_{xx}$ denotes the second derivative of $q$ with respect to $x$. The
asymptotic behavior of the solution is obtained in deift1995asymptotics ;
hastings1980boundary , and as $x\to-\infty$,

$\displaystyle
q(x)=\sqrt{-\frac{x}{2}}\left(1+\frac{1}{8x^{3}}-\frac{73}{128x^{6}}+\frac{10219}{1024x^{9}}+\mathcal{O}(|x|^{-12})\right),$
(62)

and

$\displaystyle q\left(x\right)\to\text{Ai}(x),\qquad\text{as}\qquad
x\to\infty.$ (63)

Using the Hastings–McLeod results hastings1980boundary  we have

$\displaystyle F_{2}(s)=\exp{-\int_{s}^{\infty}R(x)\ \differential
x},\qquad\text{where}\qquad R(x)=\int_{x}^{\infty}q^{2}(s)\ \differential s.$
(64)

Moreover, they obtained

$\displaystyle R(x)=q_{x}(x)^{2}-x\,q(x)^{2}-q(x)^{4}.$ (65)

Unlike the Gaussian distribution, TW distribution is an asymmetric
distribution and has different left and right tails. Using the Hastings–McLeod
results, asymptotic analysis of the TW distribution is performed in
baik2008asymptotics , and the following result is obtained,

$F_{2}(s)=\begin{cases}\displaystyle
1-\frac{e^{-\frac{4}{3}s^{\frac{3}{2}}}}{32\pi
s^{\frac{3}{2}}}\quantity(1-\frac{35}{24s^{\frac{3}{2}}})+\mathcal{O}(s^{-3})&\quad
s\rightarrow\infty\\\\[10.00002pt] \displaystyle
2^{\frac{1}{24}}\,e^{\zeta^{\prime}(-1)}\frac{e^{-\frac{1}{12}|s|^{3}}}{|s|^{\frac{1}{8}}}\quantity(1-\frac{3}{2^{6}\,|s|^{3}}+\mathcal{O}(s^{-6}))&\quad
s\rightarrow-\infty\\\ \end{cases},$ (66)

where $\zeta(x)$ is the Riemann zeta function. The asymmetry of the tails of
the TW distribution is apparent in the leading order asymptotic,

$F_{2}(s)=\begin{cases}1-\mathcal{O}(e^{-s^{3/2}})&\quad s\rightarrow\infty\\\
\\\ \mathcal{O}(e^{-|s|^{3}})&\quad s\rightarrow-\infty\\\ \end{cases}.$ (67)

####  4.1.1 Finite $N$ results and $1/N$ expansion

Let us define the double-scaling parameter,

$\displaystyle s=\alpha_{2}^{-\frac{1}{3}}(\beta_{c}-\beta)\ N^{\frac{2}{3}}.$
(68)

By using the double scaling parameter for general unitary matrix model, the
leading and sub-leading contributions to the free energy, can be obtained
using the asymptotic expansion of the TW distribution (66),

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\frac{1}{N^{2}}\log\left(1-(32\pi)^{-1}s^{-\frac{3}{2}}\exp{-\frac{4}{3}s^{\frac{3}{2}}}\right)&\beta<\beta_{c},\
|\beta-\beta_{c}|=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\frac{1}{N^{2}}\log\left(2^{\frac{1}{24}}e^{\zeta^{\prime}(-1)}e^{-\frac{1}{12}|s|^{3}}|s|^{-\frac{1}{8}}\left(1-\frac{3}{2^{6}\,|s|^{3}}\right)\right)&\beta>\beta_{c},\
|\beta-\beta_{c}|=\mathcal{O}(1)\end{cases}.$ (69)

We can add the crossing region, between the two tails, to the above result and
after expanding the logarithm, we obtain

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}-\frac{1}{N^{2}}\sum_{n=1}^{\infty}\frac{1}{n}\,(32\pi)^{-n}s^{-\frac{3}{2}n}\exp({-\frac{4}{3}ns^{\frac{3}{2}}})&\beta<\beta_{c},\
|\beta-\beta_{c}|=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\frac{1}{N^{2}}\int_{s}^{\infty}(x-s)\ q^{2}(x)\
\differential x&\beta<\beta_{c},\
|\beta-\beta_{c}|=\mathcal{O}(N^{-2/3})\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\frac{1}{N^{2}}\left(c-\frac{|s|^{3}}{12}-\frac{1}{8}\log|s|-\sum_{n=1}^{\infty}\frac{1}{n}\left(\frac{3}{2^{6}}\right)^{n}|s|^{-3n}\right)&\beta>\beta_{c},\
|\beta-\beta_{c}|=\mathcal{O}(1)\end{cases},$ (70)

where $c=\frac{\log 2}{24}+\zeta^{\prime}(-1)$. We can recast the expansions
as

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\sum_{n=1}^{\infty}N^{-n-2}\
G_{n}^{\ (1)}\ e^{-2Nf_{n}(\beta)}&\beta<\beta_{c},\
|\beta-\beta_{c}|=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\sum_{n=0}^{\infty}N^{-\frac{2}{3}n-2}\ G_{n}^{\
(2)}&\beta<\beta_{c},\ |\beta-\beta_{c}|=\mathcal{O}(N^{-2/3})\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\mathcal{F}_{0}+N^{-2}\mathcal{G}-\sum_{n=1}^{\infty}N^{-2n-2}\
G_{n}^{\ (3)}&\beta>\beta_{c},\ |\beta-\beta_{c}|=\mathcal{O}(1)\end{cases},$
(71)

where in the right tail, the first line, we have

$\displaystyle G_{n}^{\
(1)}=-\frac{1}{n}(32\pi)^{-n}\alpha_{2}^{\frac{n}{2}}\,(\beta_{c}-\beta)^{-\frac{3n}{2}},\qquad
f_{n}(\beta)=\frac{2}{3}n\,\alpha_{2}^{-\frac{1}{2}}\
(\beta_{c}-\beta)^{\frac{3}{2}},$ (72)

and in the left tail, third line,

$\displaystyle\mathcal{F}_{0}=-\frac{\alpha_{2}^{-1}}{12}|\beta_{c}-\beta|^{3},\quad\mathcal{G}=\frac{\log
2}{24}+\zeta^{\prime}(-1)+\frac{1}{24}\log\alpha_{2}-\frac{1}{12}\log
N-\frac{1}{8}\log|\beta_{c}-\beta|,$

and

$\displaystyle G_{n}^{\ (3)}=\frac{1}{n}\left(\frac{3}{2^{6}}\right)^{n}\
\alpha_{2}^{n}\ |\beta_{c}-\beta|^{-3n}.$ (74)

The above result reproduces the known results about the scaling behavior of
the free energy of GWW model liu2004fine  and moreover generalize it to the
generalized GWW model alvarez2006black . We will elaborate more on this in
Sections 4.2 and 4.3.

#### Right tail of TW distribution and instantons

The right tail of the TW distribution is obtained by taking the limit
$s\to\infty$, thus we can expand the exponential in Eq.(60), as the integrand
becomes small in this limit,

$\displaystyle F_{2}(s)\approx 1-\int_{s}^{\infty}(x-s)\ q^{2}(x)\
\differential x,$ (75)

then, using the asymptotic result

$\displaystyle q\left(s\right)\to\text{Ai}(s),\quad\text{if}\quad s\to\infty,$
(76)

and the Airy function asymptotics,

$\displaystyle\text{Ai}(s)=\frac{1}{2\sqrt{\pi}s^{1/4}}\exp{-\frac{2}{3}s^{\frac{3}{2}}}\left(1-\frac{5/27}{\frac{2}{3}s^{\frac{3}{2}}}+\mathcal{O}\left(\frac{1}{|s|^{3}}\right)\right),$
(77)

by keeping the first term in the above expansion of the Airy function, and
using it in Eq.(75), we obtain

$\displaystyle F_{2}(s)\approx 1-\frac{1}{4\pi}\int_{s}^{\infty}(x-s)\
x^{-\frac{1}{2}}\exp{-\frac{4}{3}x^{\frac{3}{2}}}\ \differential x.$ (78)

The above integral can be evaluated as

$\displaystyle F_{2}(s)\approx
1-\frac{1}{8\pi}\exp{-\frac{4}{3}s^{\frac{3}{2}}}+\frac{1}{6\pi}s^{\frac{3}{2}}E_{\frac{2}{3}}\left(\frac{4}{3}s^{\frac{3}{2}}\right),$
(79)

where $E_{n}(x)$ is the generalized exponential integral defined by

$\displaystyle E_{n}(x)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\ \differential
t.$ (80)

Using the asymptotic expansion of the generalized exponential function at
$x\to\infty$,

$\displaystyle
E_{n}(x)=\frac{e^{-x}}{x}\left(1-\frac{n}{x}+\frac{n(n+1)}{x^{2}}-...\right),$
(81)

we obtain

$\displaystyle F_{2}(s)\approx
1-(16\pi)^{-1}s^{-\frac{3}{2}}\exp{-\frac{4}{3}s^{\frac{3}{2}}}\left(1-\frac{5}{4s^{\frac{3}{2}}}+\mathcal{O}(s^{-3})\right).$
(82)

The above right tail expansion, $s\to\infty$ of the TW distribution, matches,
up to some numerical coefficients, with the similar expansion in the
literature, for example in baik2008asymptotics ; borot2012right ,

$\displaystyle
F_{2}(s)=1-(32\pi)^{-1}s^{-\frac{3}{2}}\exp{-\frac{4}{3}s^{\frac{3}{2}}}\left(1-\frac{35}{24s^{\frac{3}{2}}}+\mathcal{O}(s^{-3})\right).$
(83)

#### Noperturbative instanton sector

In the region $\beta>\beta_{c},\ (\beta-\beta_{c})=\mathcal{O}(1)$, using the
parameterization $s=(\xi N)^{\frac{2}{3}}$, we observe that there is no
perturbative corrections to the continuum free energy and the non-perturbative
corrections due to the instantons are of order $\mathcal{O}(e^{-N})$, and we
can expand the leading terms of Eq.(83) to obtain the instanton contribution
to the free energy,

$\displaystyle\mathcal{F}_{\text{inst}}=-\sum_{n=1}^{\infty}\frac{(32\pi\xi)^{-n}}{n}N^{-n-2}\exp(-\frac{4}{3}n\,\xi
N),$ (84)

where $n$ is the number of the instantons. Thus, we can write
$\displaystyle\mathcal{F}_{\text{inst}}=\sum_{n=1}^{\infty}\mathcal{F}_{n\text{-inst}}$,
and in particular, the one-instanton sector, which is the dominant
contribution, is

$\displaystyle\mathcal{F}_{1\text{-inst}}=-\frac{1}{32\pi\xi}N^{-3}\exp(-\frac{4}{3}\
\xi N).$ (85)

#### Left tail of TW distribution and genus expansion of free energy

In this part, we use the left tail expansion of the TW distribution to obtain
the genus expansion of the free energy. The genus expansion of the free energy
is defined by

$\displaystyle\mathcal{F}=\sum_{g=0}^{\infty}N^{2-2g}\mathcal{F}_{g}.$ (86)

Notice that in the genus expansion of the free energy, the definition of the
free energy is not normalized by the $1/N^{2}$ factor and for example the
genus zero free energy is of order $N^{2}$, and this should be considered when
we compare the genus expansion with the result in Section 4.1.1. In the genus
expansion of the free energy, for more convenience, we use slightly different
parameterization $s=\lambda N^{\frac{2}{3}}$, in which the two parameters are
related by $\lambda=\xi^{\frac{2}{3}}$. In this parameterization, we can
arrange the expansion of the free energy in powers of $N$ and find the
subleading terms in $1/N$ expansion. Using this parameterization we find the
total genus zero part, which is defined as the sum of the continuum free
energy and genus zero free energy in the expansion (86), as

$\displaystyle\tilde{\mathcal{F}}_{0}=\mathcal{F}_{c}+\mathcal{F}_{0}$
$\displaystyle=$ $\displaystyle\mathcal{F}_{c}+\begin{cases}\displaystyle
0&\text{for }\beta<\beta_{c}\\\\[10.00002pt]
\displaystyle-\frac{1}{12}|\lambda|^{3}&\text{for
}\beta>\beta_{c}\end{cases}.$ (87)

Furthermore, as we observed, there is no perturbative corrections
$\mathcal{O}(\frac{1}{N})$ in the region $\beta>\beta_{c}$, and the
perturbative corrections in the region $\beta<\beta_{c}$, can be computed by
using the left tail asymptotics of the $\text{TW}^{(2)}$ distribution. In
fact, the free energy at genus one, and for higher genus ($g\geq 2$), are
given by

$\displaystyle\mathcal{F}_{1}$ $\displaystyle=$ $\displaystyle\frac{\log
2}{24}+\zeta^{\prime}(-1)-\frac{1}{8}\log|\lambda|-\frac{1}{12}\log
N,\qquad\mathcal{F}_{g}=-\frac{1}{g-1}\left(\frac{3}{2^{6}}\right)^{g-1}\lambda^{3-3g}.$
(88)

Let us define the perturbative corrections to free energy by

$\displaystyle\mathcal{F}_{\text{pert}}=\mathcal{F}-\mathcal{F}_{0}=\sum_{g=1}^{\infty}N^{2-2g}\mathcal{F}_{g},$
(89)

and then by putting the perturbative and non-perturbative leading and
subleading terms together in the two phases, we obtain

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\mathcal{F}_{\text{inst}},&\text{for
}\beta<\beta_{c}\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}-\frac{1}{12}|\lambda|^{3}+\mathcal{F}_{\text{pert}},&\text{for
}\beta>\beta_{c}\end{cases}.$ (90)

Notice that, the above result for the free energy is universal, besides the
model-dependent parts, namely the continuum free energy $\mathcal{F}_{c}$, and
parameters $\xi$ and $\lambda$.

####  4.1.2 Large $N$ asymptotics: Universal phase structure

As we observed so far, the different asymptotic behavior of the TW
distribution in its left and right tail, lead to different expansion of the
free energy in the strong and weak coupling phases. In fact, there is a
universal strong-to-weak coupling transition, i.e. the deconfinement phase
transition, in the unitary matrix models. The large gap asymptotic of the
Fredholm determinant, which is the asymptotic analysis in the left tail (i.e.
the phase $\beta>\beta_{c}$) of the TW distribution, is obtained via the
Riemann–Hilbert analysis baik2008asymptotics , as $s\to-\infty$,

$\displaystyle
F_{2}(s)=2^{1/24}e^{\zeta^{\prime}(-1)}\,|s|^{-\frac{1}{8}}\exp{-\frac{|s|^{3}}{12}}\left(1+\frac{3}{2^{6}|s|^{3}}+\mathcal{O}(|s|^{-6})\right).$
(91)

Thus, we obtain the universal finite contribution from the edge fluctuation to
the free energy in strongly-coupled phase,

$\displaystyle\mathcal{F}_{f}\sim\lim_{s\to-\infty}\log
F_{2}(s)=-\frac{|s|^{3}}{12}-\frac{1}{8}\log|s|+\log C_{2}+o(1),$ (92)

where $s=\alpha_{2}^{-\frac{1}{3}}(\beta_{c}-\beta)\ N^{\frac{2}{3}}$, and the
parameters $\alpha$ and $\beta$ depend on the model, and
$C_{2}=2^{1/24}e^{\zeta^{\prime}(-1)}$.

Putting together the left and right tails expansion of the TW distribution,
the leading order free energy and the order of the subleading perturbative and
non-perturbative corrections, can be summarized as

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\mathcal{O}(e^{-cN})&\text{for
}\beta<\beta_{c}\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}-\frac{1}{12}|\lambda|^{3}+\mathcal{O}(N^{-2})&\text{for
}\beta>\beta_{c}\end{cases}.$ (93)

It is easy to observe that the third derivative of the free energy is
discontinuous at $\beta(t_{i})=\gamma\equiv\beta_{c}$, which is leading to a
universal third order phase transition in the model. This is a universal
generalization of the GWW phase transition at $\beta=\beta_{c}$ between the
two phases: weak coupling $\beta<\beta_{c}$ ($s\to\infty$) and strong coupling
$\beta>\beta_{c}$ ($s\to-\infty$). This is the phase transition associated
with the one-gap opening/closing in the unitary matrix model.

#### Cross-over region between the tails

Having discussed the left and right tails of the TW distribution, in this
Section, we take a closer look on the crossing regions between the tails. At
finite but large $N$, the free energy of the matrix model in the intermediate
region is obtained from the TW distribution,

$\displaystyle\mathcal{F}=\mathcal{F}_{c}-\int_{s}^{\infty}(x-s)\ q^{2}(x)\
\differential x,$ (94)

where $s$ is finite, $s\sim\alpha_{2}^{-1/3}=\mathcal{O}(1)$, as
$\beta-\beta_{c}=\mathcal{O}(N^{-2/3})$, and $q(x)$ is the solution of the
Painlevé II equation (61). Notice that, from Eq.(60), one can observe that $q$
satisfies

$\displaystyle q^{2}(s)=-\derivative[2]{s}\log F_{2}(s),$ (95)

and one can interpret $q^{2}$ as the specific heat, as $F_{2}$ has the
interpretation of the partition function in our context.

In order to approximate the $F_{2}(s)$ in the crossing region for finite $s$,
one needs to study the interpolation of the function $q(x)$ between the left
and right tails; in the limit $x\to-\infty$, the asymptotic behavior is given
by the Eq.(62) and in the limit $x\to\infty$, it is
$q\left(x\right)\to\text{Ai}(x)$. In the crossing region, it is numerically
obtained that the TW distribution can be approximated by the Gamma
distribution chiani2014distribution ,

$\displaystyle
F_{2}(s)\approx\Gamma(k)^{-1}\gamma(k,\frac{s+\alpha}{\theta}),$ (96)

where $\gamma(k,x)$ is the lower incomplete gamma function defined as

$\displaystyle\gamma(t,x)=\int_{0}^{x}s^{t-1}e^{-s}\,\differential s,$ (97)

and the parameters $k$, $\alpha$ and $\theta$ are numerically adjusted to fit
the TW distribution.

Having discussed the applications of the asymptotic analysis of the TW
distribution in the general setting, in the rest of this Section, we apply our
result in two concrete example of the GWW model and its generalization.

###  4.2 GWW model and TW distribution

The perturbative and non-perturbative results and a universal phase structure
in the GWW matrix model can be extracted from the free energy given by

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+\lim_{N\to\infty}N^{-2}\log
F_{2}(s),\qquad\text{where}\qquad s=\lambda N^{2/3},$ (98)

and $\lambda=\alpha_{2}^{-1/3}(\beta_{c}-\beta)$, with the coefficients
$\alpha_{2}$ and $\beta$, given by

$\displaystyle\alpha_{2}=t_{1}\,,\qquad\beta=2t_{1}\,,\qquad\beta_{c}=\gamma.$
(99)

Assuming $\gamma=1$, the critical point $\beta_{c}=1$ implies that the
critical coupling is $t_{1}^{*}=1/2$. Then, we easily obtain

$\displaystyle\lambda=\frac{1-2t_{1}}{t_{1}^{1/3}}\,,\qquad\xi=\frac{(1-2t_{1})^{3/2}}{t_{1}^{1/2}},$
(100)

and we can compute the free energy (90),

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle
t_{1}^{2}+\mathcal{F}_{\text{inst}}&\text{for }t_{1}<1/2\\\\[10.00002pt]
\displaystyle
t_{1}^{2}-\frac{|1-2t_{1}|^{3}}{12t_{1}}+\mathcal{F}_{\text{pert}}&\text{for
}t_{1}>1/2\end{cases},$ (101)

where $\mathcal{F}_{\text{inst}}$ in Eq.(84) can be computed by sum over all
$n$-instantons contributions, given by

$\displaystyle\mathcal{F}_{n\text{-inst}}=\frac{1}{(32\pi)^{n}n}\frac{t_{1}^{n/2}}{(1-2t_{1})^{3n/2}}N^{-n-2}\exp(-\frac{4}{3}n\
\frac{(1-2t_{1})^{3/2}}{t_{1}^{1/2}}N),$ (102)

and $\mathcal{F}_{\text{pert}}$ is given by

$\displaystyle\mathcal{F}_{1}=\frac{\log
2}{24}+\zeta^{\prime}(-1)-\frac{1}{8}\log\left|\frac{1-2t_{1}}{t_{1}^{1/3}}\right|-\frac{1}{12}\log
N.$ (103)

and

$\displaystyle\mathcal{F}_{g}=-\frac{1}{g-1}\left(\frac{3}{2^{6}}\right)^{g-1}\left|\frac{1-2t_{1}}{t_{1}^{1/3}}\right|^{3-3g},\quad\text{for
}g\geq 2.$ (104)

To compare with results in the literature for example liu2004fine ;
okuyama2017wilson , notice that the coupling of our GWW matrix model $t_{1}$
is related to one in the literature $t$, by $t=2t_{1}$. The leading term of
the free energy of the GWW model obtained from the Coulomb gas method is

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\frac{t^{2}}{4}&\text{for
}t<1\\\\[10.00002pt] \displaystyle t-\frac{3}{4}-\frac{\log t}{2}&\text{for
}t>1\end{cases}.$ (105)

Above free energy implies there is a third order phase transition at $t=1$ in
GWW model. This model-dependent third-order phase transition is an example of
the universal third order phase transition in our discussion implied by the
universality of the TW distribution. In order to compare the above model
dependent result with the leading terms (genus zero) of the universal result
obtained in Eq.(101), we should consider the vicinity of the critical point
$t^{*}=1$. The free energy obtained by both methods in the phase $t>1$ matches
exactly. In the phase $t<1$ and in the vicinity of the critical point
$t\nearrow 1$, one can observe that the expansion of the free energy in terms
of $\epsilon=t-1$ (or $t=e^{\epsilon}$) in Eqs.(101) and (105) matches up to
order $\mathcal{O}(\epsilon^{3})$.

Let us consider the $1/N$ expansion in GWW model liu2004fine ,

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\frac{t^{2}}{4}+\frac{1}{2\pi}e^{-2Nf(t)}\sum_{n=1}^{\infty}N^{-n-2}\
F_{n}^{\ (1)}\ &\text{for }t<1,\ 1-t=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\frac{t^{2}}{4}+\sum_{n=0}^{\infty}N^{-\frac{2}{3}n-2}\ F_{n}^{\
(2)}&\text{for }t<1,\ 1-t=\mathcal{O}(N^{-2/3})\\\\[10.00002pt] \displaystyle
t-\frac{3}{4}-\frac{\log t}{2}+\sum_{n=0}^{\infty}N^{-2n-2}\ F_{n}^{\
(3)}&\text{for }t>1,\ t-1=\mathcal{O}(1)\end{cases},$ (106)

where

$\displaystyle
f(t)=\log\left(\frac{1}{t}+\sqrt{\frac{1}{t^{2}}-1}\right)-t\sqrt{\frac{1}{t^{2}}-1},\quad
F_{n}^{\ (1)}\sim\frac{1}{(1-t)^{\frac{3n}{2}}},\quad F_{n}^{\
(3)}\sim\frac{1}{(t-1)^{3n}},$

and $F_{0}^{\ (2)}$ satisfies Painlevé II equation liu2004fine . We observe an
agreement between $G_{n}^{(1,2,3)}$ in Eq.(71) and $F_{n}^{(1,2,3)}$ in above
result. In instanton sector, different approximation for Airy function is
applied in our study whereas in the literature the Airy function is
approximated by the Bessel function and its Debye approximation. However, one
can observe that the $f_{1}(\beta)$ in Eq.(72) matches with $f(t)$ around
$t=1$ up to the first term in the expansion which is of order
$\mathcal{O}(\epsilon^{3/2})$. We hope to comeback to this problem, in our
future studies.

The perturbative genus expansion in the GWW model okuyama2017wilson  is

$\displaystyle\mathcal{F}_{g}(t)=\frac{B_{2g}}{2g(2g-2)}+\frac{1}{(t-1)^{3g-3}}\sum_{n=0}^{g-2}c_{n}^{(g)}t^{n},$
(108)

where $B_{2g}$ denotes Bernoulli number. In particular, at genus one, we have

$\displaystyle\mathcal{F}_{1}(t)=\zeta^{\prime}(-1)-\frac{1}{8}\log(1-1/t)-\frac{1}{12}\log
N.$ (109)

In general, for any $g\geq 1$, we expect that the above results match with
Eqs.(103) and (104) around the critical point $t=1$. At genus one, all terms
in Eq.(109) match with Eq.(103) up to the coefficient of term $\log t$. The
first term in Eq.(108) is normalization for the Haar measure and it is
obtained by computing the volume of the $\mathrm{U}(N)$ gauge group. However,
for the second term in Eq.(108), by comparing with Eq.(104), although they are
different in order, but we expect the following relation around $t=1$,

$\displaystyle\sum_{n=0}^{g-2}c_{n}^{(g)}t^{n}\sim-\frac{1}{g-1}\left(\frac{3}{2^{7}}\right)^{g-1}t^{g-1}.$
(110)

More precise comparison and further implications of the above equation remain
for future studies.

###  4.3 Generalized GWW model and TW distribution

In this Section we explain the free energy of the one-cut solution of the
generalized GWW model. Similar to GWW model, all the results in the matrix
model with an arbitrary potential in generalized GWW can be obtained from the
free energy (98) with explicit model-dependent parameters given by

$\displaystyle\alpha_{2}(t_{1},t_{2},...)=\sum_{n=1}^{\infty}n^{3}\,t_{n}\,,\qquad\beta(t_{1},t_{2},...)=2\sum_{n=1}^{\infty}n\,t_{n}\,,\qquad\beta_{c}(t_{1},t_{2},...)=\gamma\,.$
(111)

Then, by using explicit form of $\lambda$ and $\xi=\lambda^{3/2}$ in Eqs.(84),
and (89), one can explicitly compute the free energy of the generalized GWW
model,

$\displaystyle\mathcal{F}=\sum_{n}n\,t_{n}^{2}+\begin{cases}\displaystyle\mathcal{F}_{\text{inst}}&\text{for
}\beta<\beta_{c}\\\\[10.00002pt]
\displaystyle-\frac{1}{12}|\lambda|^{3}+\mathcal{F}_{\text{pert}}&\text{for
}\beta>\beta_{c}\end{cases}.$ (112)

The above free energy of the generalized GWW model implies a universal third
order phase transition in the generalized GWW model. In the following, we
compare our result with $1/N$ expansion of the free energy in the literature.
The critical behaviour of the generalized GWW model is discussed in
alvarez2006black ,

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\sum_{n}n\
t_{n}^{2}+\frac{1}{2\pi}e^{-2Nf(t_{n})}\sum_{n=1}^{\infty}N^{-n-2}\ F_{n}^{\
(1)}\ &\text{for }t<1,\ 1-t=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\sum_{n}n\ t_{n}^{2}+\sum_{n=0}^{\infty}N^{-\frac{2}{3}n-2}\
F_{n}^{\ (2)}&\text{for }t<1,\ 1-t=\mathcal{O}(N^{-2/3})\\\\[10.00002pt]
\displaystyle H(t_{n})+\sum_{n=0}^{\infty}N^{-2n-2}\ F_{n}^{\ (3)}&\text{for
}t>1,\ t-1=\mathcal{O}(1)\end{cases},$ (113)

where $F_{n}^{\ (1,2,3)}$ and $H(t_{n})$ can be computed using the methods of
orthogonal polynomials goldschmidt19801 . In the vicinity of the critical
point, assuming $\gamma=1$,

$\displaystyle 2\sum_{n}n\,t_{n}^{*}=1,$ (114)

using the results in Eqs.(71)-(74) with the parameters $\alpha_{2}$ and
$\beta$ given by Eq.(111), in principle we can obtain explicit results for
$H(t_{n})$, $f(t_{n})$, $F_{n}^{(1,3)}$ and $F_{0}^{(2)}$. Similar to GWW
case, they can compared with $G_{n}^{1,2,3}$ in Eqs. (72) and (74) with
insertion of parameters from Eq.(111). In particular, the genus zero free
energy in the strong coupling regime is expected to be obtained from Eq.(112),
as

$\displaystyle
H(t_{n})\approx\sum_{n}n\,t_{n}^{2}-\frac{1}{12}\alpha_{2}^{-1}|1-\beta|^{3}.$
(115)

In summary, the new results include the exact and explicit results for the
$1/N$ expansion, genus expansion at finite $N$ and phase structure of the
model at large $N$.

As a simple example, we can consider the single-term generalized GWW model,
$t_{n}=t_{m}\delta_{n,m}$, with the following potential

$\displaystyle V(U)=t_{m}\tr(U^{m}+U^{-m}).$ (116)

For simplicity we write $t_{m}\equiv t$ and then in this case we obtain the
double scaling parameter

$\displaystyle\lambda=m^{-1}t^{-\frac{1}{3}}(1-2mt),$ (117)

and the free energy,

$\displaystyle\mathcal{F}=m\,t^{2}+\begin{cases}\displaystyle\mathcal{F}_{\text{inst}}&\text{for
}t<\frac{1}{2m}\\\\[10.00002pt]
\displaystyle-\frac{1}{12}m^{-3}t^{-1}|1-2mt|^{3}+\mathcal{F}_{\text{pert}}&\text{for
}t>\frac{1}{2m}\end{cases}.$ (118)

This implies that there is a third-order phase transition at critical coupling
$t^{*}=\frac{1}{2m}$. Similarly, $\mathcal{F}_{\text{inst}}$ and
$\mathcal{F}_{\text{pert}}$ can be computed explicitly.

##  5 Multi-critical dynamics

In the previous Section, we have studied the critical dynamics associated with
the one gap opening/closing in the GWW model and its generalization. In this
Section, we study the multi-critical extension of our previous results for the
generalized GWW model, which are related to the multi-gap opening/closing in
the unitary matrix model with a generic high degree polynomial potential. The
multi-critical dynamics is originated in the asymptotic behavior of the higher
TW distribution $\text{TW}^{(p)}$. The asymptotic behavior of the TW
distribution ($p=2$) is studied in baik2008asymptotics , and recent asymptotic
results for the Pearcey processes ($p=3$) are obtained in dai2020asymptotics .
See also earlier results Brezin:1998zz ; Brezin:1998PREb . The asymptotic
analysis for the generic case of higher $\text{TW}^{(p)}$ distribution is
conjectured in LeDoussal:2018dls  and proved in cafasso2019fredholm . See also
Periwal:1990qb ; Claeys:2009CPAM . In this Section we review recent results in
the large gap asymptotics of the higher TW distributions. Moreover, we obtain
the asymptotics results for the right tail expansion. Then, we apply these
results in the multi-critical dynamics of the unitary matrix model. The
analysis of this Section is the direct generalization of the one for $p=2$ in
previous Section.

###  5.1 Multi-critical dynamics and higher $\text{TW}^{(p)}$ distribution

The main result of this work is the critical/multi-critical phase structure of
the generalized GWW model. In this study, we consider $p\in 2\mathbb{N}$,
however we expect similar results for odd $p$ kimura2021 . We conjecture that
all the perturbative and non-perturbative information at the finite $N$ and
large $N$ and the (multi-)critical phase structure can be encoded in a compact
form as

$\displaystyle\lim_{\epsilon\to
0}\mathbb{P}\quantity[\frac{\lambda_{1}-c_{1}^{(p)}/\epsilon}{(c_{2}^{(p)}/\epsilon)^{\frac{1}{p+1}}}<s]=F_{p}(s)\,,$
(119)

where $c_{1}^{(p)}$ and $c_{2}^{(p)}$ are some model dependent parameters, and
$F_{p}(s)$ is $\text{TW}^{(p)}$ distribution. In the generalized GWW model, we
have the distribution (48) with parameters are explicitly given by

$\displaystyle
c_{2}^{(p)}=\alpha_{p}(t_{1},t_{2},...)=2\sum_{n=1}^{\infty}\frac{n^{p+1}}{p!}\,t_{n}\,,\qquad
c_{1}^{(p)}=\beta(t_{1},t_{2},...)=2\sum_{n=1}^{\infty}n\,t_{n}\,.$ (120)

The free energy in the large $N$ limit of the $\mathrm{U}(N)$ matrix model is
defined in (28), and it can be obtained from (13) and (119), by fixing
$\gamma:=N\epsilon$, as

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+\lim_{N\to\infty}N^{-2}\log
F_{p}(s)\qquad\text{where}\qquad s=\frac{\gamma-
c_{1}^{(p)}}{\left(c_{2}^{(p)}\right)^{\frac{1}{p+1}}}N^{\frac{p}{p+1}},$
(121)

is called the multi-critical double-scaling parameter.

#### Multi-critical double-scaling parameter

Let us define the critical point $\gamma=\beta_{c}$ and by using (120), write
double scaling parameter in generalized GWW model as

$\displaystyle s=\alpha_{p}^{-\frac{1}{1+p}}(\beta_{c}-\beta)\
N^{\frac{p}{p+1}}.$ (122)

The multi-critical dynamics and phase structure is governed by the behaviour
of free energy around the critical point $\beta=\beta_{c}$, and as we will
see, the free energy in the vicinity of the critical point and the order of
the multi-critical phase transition can be explicitly computed from the
parameters $\alpha_{p}$ and $\beta$ in the asymptotic expansion of $F_{p}(s)$
around the critical point. For real positive $t_{n}$, we have $\alpha_{p}>0,\
\beta>0$, and as $N\to\infty$, the double-scaling parameter $s$ is tending to
plus/minus infinity depending on the sign of $(\beta_{c}-\beta)$ and thus
there are two phases $\beta_{c}>\beta$ and $\beta_{c}<\beta$ in the matrix
model. Then, the free energy has different expansions in different phases.
Thus, to obtain the free energy one needs to study the asymptotics of
$F_{p}(s)$ in the two limits of $s\to\pm\infty$.

####  5.1.1 Multi-critical instantons

In this Section, we study the right tail asymptotics of $\text{TW}^{(p)}$ and
its applications for the instantons in the multi-critical model.

In the limit $s\to\infty$, one can obtain the asymptotic expansion of
$F_{p}(s)$, in a straightforward calculations, by using the definition

$\displaystyle\log F_{p}(s)=-\int_{s}^{\infty}(x-s)\
q^{2}\left((-1)^{\frac{p}{2}+1}x\right)\differential x,$ (123)

and the following asymptotic result cafasso2019fredholm ,

$\displaystyle
q\left((-1)^{\frac{p}{2}+1}s\right)\to\text{Ai}_{p}(s),\quad\text{if}\quad
s\to\infty,$ (124)

where $\text{Ai}_{p}(s)$ is the $p$-Airy function defined in (41).

The asymptotic behavior of the $p$-Airy function is obtained in
bothner2021momenta , and the following results, up to some unkown numerical
factor $a_{p}$ and $\tilde{a}_{p}$, holds

$\displaystyle\text{Ai}_{p}(s)=a_{p}\,s^{-\frac{p-1}{2p}}\exp{-\frac{p}{p+1}s^{\frac{p+1}{p}}}\left(1-\tilde{a}_{p}\,s^{-\frac{p+1}{p}}+\mathcal{O}(s^{-\frac{2(p+1)}{p}})\right).$
(125)

The factor $\tilde{a}_{p}$ does not appear in the following approximation as
we keep the leading term in the asymptotic, and we will also drop the factor
$a_{p}$ in the following computation. Similar to the discussion in the $p=2$
case, by using the leading term in the asymptotic of the higher Airy function
(125), from Eq.(123) in the limit $s\to\infty$, we obtain

$\displaystyle F_{p}(s)\approx
1-\frac{1}{2}\exp{-\frac{2p}{p+1}s^{\frac{p+1}{p}}}+\frac{p}{p+1}s^{\frac{p+1}{p}}E_{\frac{p}{p+1}}\left(\frac{2p}{p+1}s^{\frac{p+1}{p}}\right).$
(126)

Then, by expanding the generalized exponential integral Eq.(81), we obtain

$\displaystyle F_{p}(s)\approx
1-\frac{1}{4}s^{-\frac{p+1}{p}}\exp{-\frac{2p}{p+1}s^{\frac{p+1}{p}}}\left(1-\frac{2p+1}{p}s^{-\frac{p+1}{p}}+\mathcal{O}(s^{-\frac{2(p+1)}{p}})\right).$
(127)

Finally, for generic $p$, we have the following right tail expansion of the
logarithm of the higher $\text{TW}^{(p)}$,

$\displaystyle\log
F_{p}(s)\approx\log\left(1-4^{-1}s^{-\frac{p+1}{p}}\exp({-\frac{2p}{p+1}s^{\frac{p+1}{p}}})\right)=-\sum_{n=1}^{\infty}\frac{1}{4^{n}n}s^{-\frac{n(p+1)}{p}}\exp({-\frac{2np}{p+1}s^{\frac{p+1}{p}}}).$
(128)

The right tail of $\text{TW}^{(p)}$, explain the regime $\beta<\beta_{c}$, and
$(\beta-\beta_{c})=\mathcal{O}(1)$ of the generalized GWW model. Using the
parameterization $s=(\xi_{p}N)^{\frac{p}{p+1}}$, and the above asymptotic
results, the free energy in this phase is obtained as

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+\mathcal{F}_{\text{inst}}^{(p)},$
(129)

where

$\displaystyle\mathcal{F}_{\text{inst}}^{(p)}=-\sum_{n=1}^{\infty}\frac{1}{4^{n}n}\xi_{p}^{-n}N^{-n-2}\exp{-\frac{2np}{p+1}\xi_{p}N}.$
(130)

In the generalized GWW model, we have $\mathcal{F}_{c}=\sum_{n}n\,t_{n}^{2}$,
and $\xi_{p}=\alpha_{p}^{-\frac{1}{p}}(\beta_{c}-\beta)^{\frac{p+1}{p}}$.

####  5.1.2 Multi-critical genus expansions

In this part, we study the left tail asymptotics of $\text{TW}^{(p)}$ and its
application in the strong coupling phase of multi-critical model, such as the
computation of the genus-expansion of the free energy. The asymptotic
expansion of the $\text{TW}^{(p)}$ distribution, in the limit $s\to-\infty$,
is recently obtained in cafasso2019fredholm ,

$\displaystyle\log F_{p}(s)=c_{p}|s|^{\frac{2(p+1)}{p}}+c\log|s|+\log
C_{p}+o(1),$ (131)

where

$\displaystyle
c_{p}=-\frac{p^{2}}{2(p+1)(p+2)}\binom{p}{\frac{p}{2}}^{-\frac{2}{p}},$ (132)

and $c=-\frac{1}{8}$ if $p=2$, $c=-\frac{1}{2}$ if $p>2$, and $C_{p}$ is a
constant, possibly depending on $p$. Although, the perturbative subleading
corrections are not yet rigorously studied, but using the analogy to the $p=2$
case in Eq.(91), we can conjecture,

$\displaystyle
o(1)=\log\left(1+b_{1}|s|^{\frac{-2(p+1)}{p}}+\mathcal{O}\left(|s|^{\frac{-4(p+1)}{p}}\right)\right),$
(133)

for some undetermined constant $b_{1}$. More explicitly, for $p=2,4,6,$ cases
we have,

$\displaystyle\log F_{2}(s)=-\frac{|s|^{3}}{12}-\frac{1}{8}\log|s|+\log
C_{2}+o(1),$ (134a) $\displaystyle\log F_{4}(s)=-\frac{2}{45}\sqrt{6}\
|s|^{\frac{5}{2}}-\frac{1}{2}\log|s|+\log C_{4}+o(1),$ (134b)
$\displaystyle\log F_{6}(s)=-\frac{9}{560}20^{\frac{2}{3}}\
|s|^{\frac{7}{3}}-\frac{1}{2}\log|s|+\log C_{6}+o(1).$ (134c)

Using the asymptotic (131), and the double scaling parameter (122),
$s=\lambda_{p}N^{\frac{p}{p+1}}$, we obtain the free energy in the strong
coupling phase ($\beta>\beta_{c}$), as

$\displaystyle\mathcal{F}=\mathcal{F}_{c}+c_{p}|\lambda_{p}|^{\frac{2(p+1)}{p}}+N^{-2}\left(c\log|\lambda_{p}|+\frac{c\
p}{p+1}\log N+\log C_{p}\right)+o(N^{-2}).$ (135)

Furthermore, in the phase $\beta>\beta_{c}$, and
$(\beta-\beta_{c})=\mathcal{O}(1)$, we can re-arrange the above result, and
obtain the genus expansion of the free energy of the generalized GWW model,

$\displaystyle\mathcal{F}=\tilde{\mathcal{F}}_{0}+\mathcal{F}_{\text{pert}}^{(p)},$
(136)

where the genus zero part and the higher genus subleading perturbative
corrections are

$\displaystyle\tilde{\mathcal{F}}_{0}=\sum_{n=1}^{\infty}n\,t_{n}^{2}+c_{p}|\lambda_{p}|^{\frac{2(p+1)}{p}},$
(137a)
$\displaystyle\mathcal{F}_{\text{pert}}^{(p)}=\mathcal{F}_{1}^{(p)}+\sum_{g=2}^{\infty}N^{2-2g}\mathcal{F}_{g}^{(p)},$
(137b)

with the genus one free energy is given by

$\displaystyle\mathcal{F}_{1}^{(p)}=-\frac{1}{2}\log|\lambda_{p}|-\frac{p}{2(p+1)}\log
N+\log C_{p}.$ (138)

For higher genus $g\geq 2$, by using the expansion of the conjecture, one can
show that the higher than two genus free energy is given by (133),

$\displaystyle
o(1)\approx\sum_{n=1}^{\infty}(-1)^{n-1}\frac{b_{1}^{n}}{n}|s|^{-\frac{2(p+1)}{p}n},\qquad\mathcal{F}_{g}^{(p)}=\frac{(-1)^{g}}{g-1}\left(b_{1}\right)^{g-1}|\lambda_{p}|^{\frac{2(p+1)}{p}(1-g)}.$
(139)

To summarize the results obtained in this Section for the free energy in the
weak and strong phases of the multi-critical generalized GWW model, we have

$\displaystyle\mathcal{F}^{(p)}=\sum_{n}n\,t_{n}^{2}+\begin{cases}\displaystyle\mathcal{F}_{\text{inst}}^{(p)}&\text{for
}\beta<\beta_{c}\\\\[10.00002pt] \displaystyle
c_{p}|\lambda_{p}|^{\frac{2(p+1)}{p}}+\mathcal{F}_{\text{pert}}^{(p)}&\text{for
}\beta>\beta_{c}\end{cases},$ (140)

where $\mathcal{F}_{\text{inst}}^{(p)}$ and $\mathcal{F}_{\text{pert}}^{(p)}$
are given by Eqs. (130) and (137b), respectively.

#### Cross-over region

In this part, we discuss the smooth cross-over region between the weak and
strong coupling phases, at large, but finite $N$, in the vicinity of the
critical point $\beta_{c}$; $\beta>\beta_{c}$ and
$(\beta-\beta_{c})=\mathcal{O}(N^{-\frac{p}{p+1}})$. The cross-over region is
the intermediate region between the right and left tails, i.e.
$\text{TW}^{(p)}(s)$ for finite $s\sim\alpha_{p}^{-\frac{1}{1+p}}$. The free
energy of the generalized GWW multi-critical model in this region can be
obtained from Eq.(123), as

$\displaystyle\mathcal{F}=\sum_{n=1}^{\infty}n\,t_{n}^{2}-N^{-2}\int_{s}^{\infty}(x-s)\
q^{2}\left((-1)^{\frac{p}{2}+1}x\right)\differential x\qquad\text{where}\qquad
s\sim\alpha_{p}^{-\frac{1}{1+p}}=\mathcal{O}(1),$ (141)

and $q(s)$ is the solution of the $\frac{p}{2}$-th member of the Painlevé II
hierarchy cafasso2019fredholm , see also krajenbrink2020painleve ,

$\displaystyle\left(\derivative{s}+2q\right)\
\mathcal{L}_{p/2}[q_{s}-q^{2}]=s\,q,$ (142)

where $q_{s}$ is the derivative of $q$ w.r.t. $s$, and operators
$\mathcal{L}_{n}$ are Lenard operators defined by

$\displaystyle\derivative{s}\mathcal{L}_{j+1}f=\left(\derivative[3]{s}+4f\derivative{s}+2f_{s}\right)\mathcal{L}_{j}f,\quad\mathcal{L}_{0}f=\frac{1}{2},\quad\mathcal{L}_{j}1=0,\
j\geq 1.$ (143)

Moreover, the multi-critical specific heat satisfies

$\displaystyle q^{2}\left((-1)^{\frac{p}{2}+1}x\right)=-\derivative[2]{s}\log
F_{p}(s),$ (144)

and in the strong coupling phase, as $s\to-\infty$, it has the following
expansion cafasso2019fredholm ,

$\displaystyle
q\left((-1)^{\frac{p}{2}+1}s\right)=\left(\frac{(\frac{p}{2})!^{2}}{(p)!}|s|\right)^{\frac{1}{p}}+\frac{c}{2}\left(\frac{(p)!}{(\frac{p}{2})!^{2}}\right)^{\frac{1}{p}}|s|^{-2-\frac{1}{p}}+\mathcal{O}\left(|s|^{-2-\frac{2}{p}}\right),$
(145)

and, in the weak coupling phase, as $s\to\infty$, we have
$q((-1)^{\frac{p}{2}+1}s)\to\text{Ai}_{p}(s)$. One can study the interpolation
of $q(s)$ between these two limits and extract the finite $s$ results.

####  5.1.3 Multi-critical $1/N$ expansion

Similar to the discussion in Section 4.1.1, our results in this Section about
the large $N$ dependence of the free energy of the multi-critical generalized
GWW model in vicinity of the critical point, can be summarized in the
following $1/N$ expansion,

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\sum_{n=1}^{\infty}N^{-n-2}\
G_{n}^{\ (p,1)}\ e^{-2Nf_{n}^{(p)}(\beta)}&\text{for
}\beta<\beta_{c},|\beta-\beta_{c}|=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\sum_{n=1}^{\infty}N^{-\frac{p}{p+1}n-2}\
G_{n}^{\ (p,2)}&\text{for
}\beta<\beta_{c},|\beta-\beta_{c}|=\mathcal{O}(N^{-\frac{p}{p+1}})\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\mathcal{F}_{0}^{(p)}+N^{-2}\mathcal{G}^{(p)}-\sum_{n=1}^{\infty}N^{-2n-2}\
G_{n}^{\ (p,3)}&\text{for
}\beta>\beta_{c},|\beta-\beta_{c}|=\mathcal{O}(1)\end{cases},$ (146)

where $\mathcal{F}_{c}=\sum_{n}n\,t_{n}^{2}$, and it is straightforward to
observe that $G_{n}^{\ (p,1)}$ and $f_{n}^{(p)}(\beta)$ in the weak coupling
phase are given by

$\displaystyle G_{n}^{\ (p,1)}=\frac{1}{4^{n}n}\ \alpha_{p}^{\frac{n}{p}}\
(\beta_{c}-\beta)^{-\frac{(p+1)n}{p}},\qquad
f_{n}^{(p)}(\beta)=\frac{np}{p+1}\ \alpha_{p}^{-\frac{1}{p}}\
(\beta_{c}-\beta)^{\frac{p+1}{p}},$ (147a) and in the strong coupling phase,
for the genus expansion of the free energy we have
$\displaystyle\mathcal{F}_{0}^{(p)}=c_{p}\
\alpha_{p}^{-\frac{2}{p}}|\beta_{c}-\beta|^{\frac{2(p+1)}{p}},$ (147b)
$\displaystyle\mathcal{G}^{(p)}=\log
C_{p}+\frac{1}{2(1+p)}\log\alpha_{p}-\frac{p}{2(p+1)}\log
N-\frac{1}{2}\log|\beta_{c}-\beta|,$ (147c) and $\displaystyle G_{n}^{\
(p,3)}=\frac{(-1)^{n-1}\ b_{1}^{n}}{n}\ \alpha_{p}^{\frac{2n}{p}}\
|\beta_{c}-\beta|^{-\frac{2n(p+1)}{p}}.$ (147d)

###  5.2 Infinite limit ($p\to\infty$) of multi-critical dynamics and bulk
fluctuation

Let us consider the generalized GWW model with the infinite degree polynomial
potential and study the multi-critical dynamics at $p\to\infty$. In this limit
the double scaling parameter (122) becomes

$\displaystyle s=e^{-1}(\beta_{c}-\beta)\,N,$ (148)

where we used the fact that the limit $p\to\infty$ in definition of
$\alpha_{p}$ in Eq. (120), implies that only the term with $n\sim p$ survives,
and we have

$\displaystyle\lim_{p\to\infty}\alpha_{p}^{-\frac{1}{1+p}}\sim\lim_{p\to\infty}\left(\frac{p^{p+1}}{p!}\right)^{-\frac{1}{1+p}}=e^{-1}.$
(149)

It is interesting to compute the multi-critical $1/N$ expansion in this limit,

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\sum_{n=1}^{\infty}N^{-n-2}\
G_{n}^{\ (\infty,1)}\ e^{-2Nf_{n}^{(\infty)}(\beta)}&\text{for
}\beta<\beta_{c},|\beta-\beta_{c}|=\mathcal{O}(1)\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\sum_{n=1}^{\infty}N^{-n-2}\ G_{n}^{\
(\infty,2)}&\text{for
}\beta<\beta_{c},|\beta-\beta_{c}|=\mathcal{O}(N^{-1})\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}+\mathcal{F}_{0}^{(\infty)}+N^{-2}\mathcal{G}^{(\infty)}-\sum_{n=1}^{\infty}N^{-2n-2}\
G_{n}^{\ (\infty,3)}&\text{for
}\beta>\beta_{c},|\beta-\beta_{c}|=\mathcal{O}(1)\end{cases},$ (150)

where one can directly compute the limit $p\to\infty$ of the functions in
Eqs.(147), as

$\displaystyle G_{n}^{\
(\infty,1)}=\frac{e^{n}}{4^{n}n}(\beta_{c}-\beta)^{-n},\qquad
f_{n}^{(\infty)}(\beta)=\frac{n}{e}\ (\beta_{c}-\beta),$ (151a)
$\displaystyle\mathcal{F}_{0}^{(\infty)}=-\frac{1}{8e^{2}}|\beta_{c}-\beta|^{2},\qquad\mathcal{G}^{(\infty)}=\log
C_{p}+\frac{1}{2}-\frac{1}{2}\log N-\frac{1}{2}\log|\beta_{c}-\beta|,$ (151b)
$\displaystyle G_{n}^{\ (\infty,3)}=\frac{(-1)^{n-1}\ b_{1}^{n}\
e^{2n}}{n}|\beta_{c}-\beta|^{-2n}.$ (151c)

In this limit, the free energy expansion in the strong coupling phase
$\beta>\beta_{c}$, becomes

$\displaystyle\mathcal{F}^{(\infty)}(\beta)=\mathcal{F}_{c}-\frac{1}{8e^{2}}|\beta_{c}-\beta|^{2}-\frac{1}{2N^{-2}}\left(\log|\beta_{c}-\beta|+\log
N-1\right)+\mathcal{O}(N^{-4}),$ (152)

and the rest of the subleading terms are

$\displaystyle\sum_{n=1}^{\infty}\frac{(-1)^{n-1}\
b_{1}^{n}\,e^{2n}}{n}|\beta_{c}-\beta|^{-2n}N^{-2n-2}.$ (153)

The leading contribution to free energy indicates a second order phase
transition. Moreover, the higher genus ($g\geq 2$) free energy becomes

$\displaystyle\mathcal{F}_{g}^{(\infty)}(\beta)=\frac{(-1)^{g}\
\left(e^{2}\,b_{1}\right)^{g-1}}{g-1}|\beta_{c}-\beta|^{2-2g}.$ (154)

In the weak coupling phase $\beta<\beta_{c}$, using the asymptotic behavior of
the higher Airy function in Eq.(125),

$\displaystyle\lim_{p\to\infty}\text{Ai}_{p}(s)=a_{\infty}s^{-\frac{1}{2}}\
e^{-s}\left(1-\frac{\tilde{a}_{\infty}}{s}+\mathcal{O}\quantity(\frac{1}{s^{2}})\right),$
(155)

and the infinite limit of the $\text{TW}^{(p)}$ in Eq.(126) is

$\displaystyle\lim_{p\to\infty}F_{p}(s)\approx 1-\frac{e^{-2s}}{2}+s\
E_{1}(2s)=1-\frac{e^{-2s}}{2}+s\ \Gamma(0,2s),$ (156)

where $\Gamma(t,x)$ is the upper incomplete Gamma function defined as

$\displaystyle\Gamma(t,x)=\int_{x}^{\infty}s^{t-1}e^{-s}\,\differential s,$
(157)

and we used $\Gamma(t,x)=x^{t}\ E_{1-t}\left(x\right)$. Then, we can expand
the exponential integral or the Gamma function to obtain

$\displaystyle\lim_{p\to\infty}F_{p}(s)\approx
1-\frac{1}{4}s^{-1}e^{-2s}\left(1-2s^{-1}+\mathcal{O}(s^{-2})\right),$ (158)

and then after keeping the first leading terms, we obtain

$\displaystyle\lim_{p\to\infty}\log
F_{p}(s)\approx\log\left(1-4^{-1}s^{-1}e^{-2s}+\mathcal{O}(e^{-2s}s^{-2})\right)=-\sum_{n=1}^{\infty}\frac{1}{4^{n}n}s^{-n}e^{-2ns}+\mathcal{O}(s\log
s).$ (159)

Finally, the infinite limit of the multi-critical instanton free energy can be
obtained from Eq.(130), as

$\displaystyle\mathcal{F}^{(\infty)}_{\text{inst}}(\beta)=-\sum_{n=1}^{\infty}\frac{e^{n}}{4^{n}n}(\beta_{c}-\beta)^{-n}N^{-n-2}\exp{-2ne^{-1}(\beta_{c}-\beta)N},$
(160)

and the one-instanton sector is

$\displaystyle\mathcal{F}^{(\infty)}_{1\text{-inst}}(\beta)=-4^{-1}e(\beta_{c}-\beta)^{-1}N^{-3}\exp{-2e^{-1}(\beta_{c}-\beta)N}.$
(161)

###  5.3 Multi-critical phase structure and its interpretation

Based on our analysis in this Section, we observe that in the generalized GWW
model there is a multi-critical phase transition at critical point
$\beta_{c}=\gamma$, meaning that the free energy has different expansion at
left and right side of the critical point,

$\displaystyle\mathcal{F}^{(p)}=\begin{cases}\displaystyle\mathcal{F}_{c}+\mathcal{O}(e^{-cN})&\text{for
}\beta<\beta_{c}\\\\[10.00002pt] \displaystyle\mathcal{F}_{c}+c_{p}\
|\lambda_{p}|^{\frac{2(p+1)}{p}}+\mathcal{O}(N^{-2})&\text{for
}\beta>\beta_{c}\end{cases}.$ (162)

where $\mathcal{F}_{c}=\sum_{n}n\,t_{n}^{2}$,

$\displaystyle
c_{p}=-\frac{p^{2}}{2(p+1)(p+2)}\binom{p}{\frac{p}{2}}^{-\frac{2}{p}},\qquad\text{and}\quad\lambda_{p}=\alpha_{p}^{-\frac{1}{1+p}}(\beta_{c}-\beta).$
(163)

The leading order of the free energy in two phases, imply that the order of
the multi-critical phase transition is $\frac{2(p+1)}{p}$. Often, in the
literature this multi-critical phase transition, for any finite $p$, is
considered a transition of order $\lfloor\frac{2(p+1)}{p}\rfloor=3$.

####  5.3.1 Multi-critical gap dynamics

Before discussing the interpretations of the multi-critical phase transition,
let us briefly explain some related issues. The assumption
$\alpha_{p^{\prime}}=0$ for $p^{\prime}<p$, which is used to derive the
coefficients (120) in Eq.(119), implies some physical constraints on the
dynamics of the model. Notice that at any fixed $p$, the interaction and
fluctuation are at scale $N^{-\frac{p}{p+1}}$, and thus in order to see the
sub-dominant fluctuations for higher $p$, all the larger scale interactions at
$p^{\prime}<p$, should be turned off as we want to isolate the interaction at
the scale $p$, and study the effect of a fixed scale fluctuation. In general,
to study the multi-critical dynamics, one needs the fine tuning of the
coupling, for example in this case, introducing $p$-dependent, at any $p\geq
2$, couplings in the large $N$ limit, such as $t_{n}^{(p)}\sim
q^{(p-p^{\prime})N}$, for $q<1$ and $2\leq p^{\prime}\leq p$.

In the critical dynamics, we studied the third-order phase transition
associated with one-gap dynamics. More precisely, in the generalized GWW model
with the maximum degree $m$ polynomial potential,

$\displaystyle V(U)=\sum_{n=1}^{m}t_{n}\tr(U^{n}+U^{-n}),$ (164)

the critical dynamics is causing the transition between $l$-cut solution and
$(l-1)$-cut solution for $1\leq l\leq m$. In the multi-critical dynamics we
must consider the multi-gap dynamics. The unitary matrix model with degree $m$
polynomial potential can have $1$\- to $m$-cut solutions and could possibly
have phase transitions between $l$-cut solution to $(l-r)$-cut solution for
$0\leq r\leq l\leq m$. These phase transitions are associated with the
colliding of the two or more end points of the cuts, simultaneously at the
same critical point $\beta_{c}$, which depends only on the potential. Let us
isolate a particular set of phase transitions, namely the transitions between
the $\frac{p}{2}$-cut solutions to full support (zero gap) solution, in which
$p$-end points of the cuts collide simultaneously, for $p\leq 2m$. We call
this transition a $p$-multi-critical dynamics for $p>2$. The case $p=2$ is
called critical dynamics. There are $m$ possible phase transitions in this
set, labeled by $2\leq p\leq 2m$, and each different $p$-(multi)-critical
dynamics is governing the opening/closing of $\frac{p}{2}$-gaps as a result of
the edge fluctuation of order $N^{-\frac{p}{p+1}}$ at critical point. In our
study, there is no mixing between the dynamics and fluctuations at different
$p$, as we isolate the fluctuation scale $N^{-\frac{p}{p+1}}$ by putting
$\alpha_{p^{\prime}}=0$, for $p^{\prime}<p$. We did compute the free energy of
the generalized GWW model near the critical point $\beta_{c}$ with the $p$-end
points of the cuts are colliding simultaneously.

####  5.3.2 A new second-order phase transition

In the limit $p\to\infty$, the asymptotic behavior of the Fredholm determinant
(131) is

$\displaystyle\lim_{p\to\infty}F_{p}(s)\sim|s|^{-\frac{1}{2}}e^{-\frac{1}{8}|s|^{2}},$
(165)

where the scaling parameter is $s\sim(\beta_{c}-\beta)N$. This is similar to
the asymptotics of the Fredholm determinant with the sine kernel discussed in
Section 3.2, and in fact indicating of the bulk fluctuation. From the physical
point of view, in the infinite limit, there are infinite number of the gaps as
well as infinite number of the cuts. Then, the infinite-multi-critical
dynamics is related to the shrinking of the infinite gaps at a critical point,
when infinite end points of the cuts collide, simultaneously. The infinite
number of the gaps can be seen as the phase of the full gap, i.e. no support
for eigenvalues. Thus, when all the gaps shrink simultaneously, the model
undergoes a transition to the phase of full support for the eigenvalues or in
other words, the infinite number of cuts collide and eigenvalues cover the
full circle. As we observe from Eq.(165), the phase transition between the
full support and full gap is of the second order. This is a phase transition
associated with the infinite number of edge fluctuations or in other words a
zero-to-full support transition driven by the bulk fluctuation.

### Higher Plancherel partition and multi-critical models

In this part, we consider some basic examples of the unitary matrix models
with multi-critical dynamics. First, consider the degree two polynomial
potential with the couplings,

$\displaystyle t_{n}=t_{1}\delta_{n,1}+t_{2}\delta_{n,2}.$ (166)

In this model, the multi-critical parameter is restricted to $p\leq 4$, and
the parameters can be easily computed

$\displaystyle\alpha_{p}=\frac{2}{p!}(t_{1}+2^{1+p}t_{2})\,,\quad\beta=2(t_{1}+2t_{2})\,,\quad\lambda_{p}=2^{-\frac{1}{1+p}}(1-2t_{1}-4t_{2})\left(\frac{t_{1}+2^{1+p}t_{2}}{p!}\right)^{-\frac{1}{1+p}}.$
(167)

The critical point $\beta_{c}=1$ implies that the surface of critical
couplings is $2t_{1}^{*}+4t_{2}^{*}-1=0$.

Using the above parameters, we can compute the first few terms of the free
energy in the strong coupling regime can be computed from Eq.(135), as

$\displaystyle\mathcal{F}\approx
t_{1}^{2}+2t_{2}^{2}+c_{p}|\lambda_{p}|^{\frac{2(p+1)}{p}}+\frac{c}{N^{2}}\log|\lambda_{p}N^{\frac{p}{p+1}}|.$
(168)

The higher order corrections to the free energy can be computed in a
straightforward way. Let us consider the critical and multi-critical cases
separately. In the critical case $p=2$, the double-scaling parameter is

$\displaystyle\lambda_{2}=\frac{1-2t_{1}-4t_{2}}{(t_{1}+8t_{2})^{1/3}},$ (169)

and the first few terms of the free energy are obtained as

$\displaystyle\mathcal{F}\approx
t_{1}^{2}+2t_{2}^{2}-\frac{(-1+2t_{1}+4t_{2})^{3}}{12(t_{1}+8t_{2})}-\frac{1}{8N^{2}}\log\left|\frac{1-2t_{1}-4t_{2}}{(t_{1}+8t_{2})^{1/3}}\right|-\frac{\log
N}{12N^{2}}.$ (170)

In the multi-critical case $p=4$, in principle we can fine tune the couplings
$t_{1}$, and $t_{2}$ such that $\alpha_{2}\to 0$, and similarly one can
compute,

$\displaystyle\lambda_{4}=\frac{2^{2/5}\,3^{1/5}\,(1-2t_{1}-4t_{2})}{(t_{1}+32t_{2})^{1/5}},$
(171)

and the leading orders in free energy as

$\displaystyle\mathcal{F}\approx
t_{1}^{2}+2t_{2}^{2}+\frac{4}{15}\sqrt{2}\left(\frac{1-2t_{1}-4t_{2}}{(t_{1}+32t_{2})^{1/5}}{}\right)^{5/2}-\frac{1}{2N^{2}}\log\left|\frac{1-2t_{1}-4t_{2}}{(t_{1}+32t_{2})^{1/5}}\right|-\frac{2\log
N}{5N^{2}}.$ (172)

Moreover, the above results can be expanded around the critical surface, for
example in terms of the critical coupling $t_{1}^{*}$, by using
$t_{2}^{*}=(1-2t_{1}^{*})/4$.

##  6 Supersymmetric indices: Examples

In this Section, we apply the machinery of the generalized GWW model to study
the matrix integral representation of the indices of the gauge theory in the
large $N$ limit. Although the exact relation is provided through
Hubbard–Stratonovich transformation, we consider this approximation as a toy
model for gauge theory indices, and use the same couplings, from generalized
GWW model, in the potential of the matrix integrals of the gauge theories.
Using this approximation, we study the perturbative and non-perturbative
regimes of gauge theory and explicitly compute the free energy and phase
structure, in some concrete examples with some interesting choices, motivated
by gauge theory, of finely tuned couplings such that the matrix model is
solvable. The results in this part are obtained from the direct computations
of the contributions of the left and right tails of $\text{TW}^{(p)}$, to the
free energy in Section 5.

###  6.1 Hagedorn phase transition and deconfinement

Before, diving to the explicit computations in some concrete examples, let us
discuss the Hagedorn phase transition in gauge theories and its relation to
the deconfinement transition in the generalized GWW model, discussed in this
paper, so far. The Hagedorn transition is the phase transition, in the large
$N$ limit, between the domain of the convergence and divergence in the
generating function, separated by the hyper-surface of the singularities of
the generating functions, corresponding to Hagedorn phase transition point.
The matrix integral of the gauge theory indices $\mathcal{I}$ in Eq.(6), in
the large $N$ limit, using the Gaussian matrix integration, becomes an
infinite product Aharony:2003sx ,

$\displaystyle\lim_{N\rightarrow\infty}\mathcal{I}(q_{i})=\prod_{n=1}^{\infty}\frac{1}{1-f(q_{i}^{n})}.$
(173)

Thus, using the above generating function at large $N$ limit, we can study the
singularities of the generating function at the poles of the infinite product
formula. First, consider $n=1$ which corresponds to the double-trace matrix
model and can effectively be approximated by GWW model. Notice that in this
case, the Hagedorn critical point and GWW critical point coincide
$t^{*}=f(q_{i}^{*})=1$. In the general case, we observe that the critical
hyper-surface of the Hagedorn transition is given by

$\displaystyle\prod_{n=1}^{\infty}\left(1-f(q_{i}^{*n})\right)=1-\sum_{n=1}^{\infty}f(q_{i}^{*n})+\mathcal{O}(f^{2}(q_{i}^{*n}))=0.$
(174)

On the other hand, in the generalized GWW model, the parameters $\alpha$ and
$\beta$ depend on the couplings $t_{n}$, and the couplings can be any
univariate or multivariate function of the parameters of the gauge theory
$q_{i}$, thus we have

$\displaystyle\alpha_{p}(q_{i})=2\sum_{n=1}^{\infty}\frac{n^{p}}{p!}\,f(q_{i}^{n})\,,\qquad\beta(q_{i})=2\sum_{n=1}^{\infty}f(q_{i}^{n})\,.$
(175)

The strong-weak coupling or the deconfinement phase transition in this model,
which is a generalization of the Gross–Witten phase transition, happens at
$\beta_{c}(q_{i})=\gamma=1$, between the two phases of weak-coupling
$\beta(q_{i})<1$, and strong-coupling $\beta(q_{i})>1$.

Using the couplings $t_{n}=f(q_{i}^{n})/n$, and re-scaling all of them by a
factor $1/2$, similarly to what we have in the original GWW model, the
critical hyper-surface of deconfinement transition in the generalized GWW
model is given by

$\displaystyle\beta_{c}=\sum_{n=1}^{\infty}n\,t^{*}_{n}=\sum_{n=1}^{\infty}f(q_{i}^{*n})=1.$
(176)

Then, by expanding the Hagedorn critical hyper-surface (174) and comparing it
with deconfinement critical hyper-surface (176), we observe that they match up
to linear order in $f(q_{i})$. The single letter indices $f(q_{i})$ in the
domain of convergence, $f(q_{i})<1$, approaches the Hagedorn critical hyper-
surface, $f(q_{i})\nearrow 1$, and thus the deconfinement transition in
generalized GWW model effectively explains the Hagedorn phase transition.

In summary, we observe that the generalized GWW model and generalized double-
trace model have qualitatively similar phase structure. In fact, upon
identification of the couplings in two models, the phase structure of
generalized GWW model captures the dominant contribution, i.e. linear
approximation in Hagedorn phase structure. This implies a universal
deconfinement/Hagedorn phase transition at critical point $\beta_{c}=1$, for
any arbitrary coefficient $f(q_{i})$ of the potential in the matrix integral
of the gauge theory. This is the phase transition associated with the gap
opening in the multi-critical matrix model.

A priori, the gauge theory double-trace matrix model is approximated by the
single-trace GWW model in the large $N$ limit, thus we expect that the results
of GWW model carry over to the double-trace model, at least up to leading
order and possibly extending to first few subleading corrections. A possible
clue for roughly examining that to what extent the subleading corrections of
single-trace model apply to double-trace model, comes from the observation
made in this section, that the critical point of the single-trace and double
trace models coincide up to linear order in $f(q_{i}^{*n})$, see Eq.(176). At
$n=1$ case, the critical temperature of the GWW model exactly matches with
that of double-trace model. This indicates that results for $n=1$ case, can
more accurately transfer to the double-trace model. Moreover, we can expect
that in the large $N$ limit, all the results about the free energy and phase
structure can safely carry over up to linear order in $f(q_{i}^{*n})$. More
precise statement about the subleading corrections at finite $N$ requires an
exact analysis of the Hubbard-Stratonovich transformation but we expect that
the Legendre transform of the current leading and subleading results would
apply to gauge theory indices. Finally, notice that in this paper, we do not
directly study the double-trace matrix model, and only consider the toy models
which are generalized GWW models with the coefficients in the potential
obtained from the double-trace models.

### Trivial example

Before studying some physically motivated examples, we consider an interesting
mathematical example with the fine-tuned couplings
$t_{n}^{(p)}=q^{(p-p^{\prime})N}/n^{s}$ which, in the limit $N\to\infty$, is
zero for $2\leq p^{\prime}<p$, and $t_{n}=1/n^{s}$, for $p^{\prime}=p$. This
example corresponds to the matrix model with the following poly-logarithmic
potential

$\displaystyle V(U)=\sum_{n}t_{n}(\tr U^{n}+\tr
U^{-n})=\left(\text{Li}_{s}(e^{\mathrm{i}\theta})+\text{Li}_{s}(e^{-\mathrm{i}\theta})\right).$
(177)

It is easy to compute the parameters of this model in terms of the Riemann
zeta function,

$\displaystyle\alpha_{p}=\frac{2}{p!}\sum_{n=1}^{\infty}n^{p+1}\,t_{n}=\frac{2}{p!}\zeta(s-p-1)\,,\qquad\beta=2\sum_{n=1}^{\infty}n\,t_{n}=2\zeta(s-1)\,.$
(178)

However, as the parameter $\beta$ is a coupling independent number, there is
no phase transition in this model. The free energy of the model is the
continuum free energy

$\displaystyle\mathcal{F}_{c}=\sum_{n}n\,t_{n}^{2}=\zeta(2s-1).$ (179)

In particular, the especial case of $s=1$ is a non-interacting theory with
zero potential, with $\alpha_{p}=0$ for all $p\in 2\mathbb{N}$. Moreover, the
free energy is divergent at $s=1$ as $\zeta(1)\to\infty$. Next, we consider a
model with an arbitrary degree single term potential,

$\displaystyle t_{n}=t_{k}^{(p)}\delta_{n,k},\qquad\text{with}\qquad
t_{k}^{(p)}=\frac{t^{k}}{k}\,q^{(p-p^{\prime})N},$ (180)

where $q<1$, and for all $p^{\prime}<p$, we have $t_{k}=0$, at large $N$
limit. In this model, after fine-tuning the coupling for any given $p$ such
that $\alpha_{p^{\prime}}=0$, for $p^{\prime}<p$, the parameters of the model
can be computed as

$\displaystyle\alpha_{p}=\frac{2k^{p}}{p!}t^{k}\,,\qquad\beta=2t^{k}\,,\qquad\lambda_{p}=2^{-\frac{1}{1+p}}(1-2t^{k})\left(\frac{k^{p}t^{k}}{p!}\right)^{-\frac{1}{1+p}},$
(181)

and the genus zero free energy in the strong coupling regime is

$\displaystyle\tilde{\mathcal{F}}_{0}^{(p)}=\frac{t^{2k}}{k}+2^{-\frac{2}{p}-1}\frac{p^{2}}{(p+1)(p+2)}\binom{p}{\frac{p}{2}}^{-\frac{2}{p}}\left(\frac{k^{p}t^{k}}{p!}\right)^{-\frac{2}{p}}\left(1-2t^{k}\right)^{\frac{2(p+1)}{p}}.$
(182)

Since $p\leq 2k$, the infinite multi-critical limit, $p\to\infty$, implies
$k\to\infty$ and one can assume $p=2k$. In this limit although $\alpha_{2}$
and $\beta$ are either zero or diverging depending on $|t|$, but
$\lambda_{\infty}$ is finite and one can show
$\lambda_{\infty}=\frac{2}{e\sqrt{t}}$ and thus the free energy becomes

$\displaystyle\tilde{\mathcal{F}}_{0}^{(\infty)}=\frac{t^{2k}}{k}-\frac{t^{-1}}{2e^{2}}.$
(183)

###  6.2 Free chiral ring index

${\bullet}$$\scriptstyle{t_{1}}$

${\bullet}$$\scriptstyle{t_{1}}$$\scriptstyle{t_{2}}$$\scriptstyle{t_{3}}$

${\bullet}$$\scriptstyle{t_{1}}$$\scriptstyle{t_{2}}$$\scriptstyle{\cdots}$$\scriptstyle{t_{k}}$

Figure 1: From left to right: Jordan quiver, clover quiver, and generalized
clover quiver

A class of gauge invariant operators which are annihilated by all the
supercharges of one chirality, for example positive R-charge, are called
chiral operators. The local gauge invariant chiral operators form a
commutative ring, called chiral ring. In this Section we consider the matrix
integral representation of the generating function for the counting of the
multi-trace chiral gauge invariant operators in the chiral ring of a generic
free $\mathcal{N}=1$ superconformal quiver theories pasukonis2013quivers .
Quiver gauge theory is a gauge theory encoded on a graph, consists of nodes
and arrows representing gauge groups and chiral fields, respectively. The
unitary gauge group $G$ of the quiver is a product over the $\mathrm{U}(N)$
factors of each node $i$ of quiver, $G=\prod_{i}\mathrm{U}(N_{i})$. The
generating function of the BPS operators in the free chiral ring is
represented by the following quiver multi-matrix integral,

$\displaystyle\mathcal{I}(q_{i})=\int_{\mathrm{U}(N)}\prod_{a}\differential{U_{a}}\exp\left[\sum_{n=1}^{\infty}\sum_{a,b,\alpha}\frac{1}{n}\left(q_{ab;\alpha}\right)^{n}\tr
U_{a}^{n}\tr U_{b}^{-n}\right],$ (184)

where the $q_{ab;\alpha}$ are the fugacity factors associated with the chiral
fields $\\{\phi_{ab;\alpha}:\alpha\in\\{1,2,...,K_{ab}\\}\\}$ with $K_{ab}$
number of arrows between nodes $a$ and $b$, and the chiral fields are
transforming in the bifundamental representation $(N_{a},\bar{N}_{b})$. The
chiral fields in the adjoint representation of a gauge group are the loops
starting and ending on the same node of the quiver associated to the gauge
group.

In the large $N$ limit, one can derive the generating function as an infinite
product,

$\displaystyle\lim_{N\rightarrow\infty}\mathcal{I}(q_{i})=\prod_{n=1}^{\infty}\frac{1}{\det\left[1-\mathcal{A}(q_{i}^{n})\right]},$
(185)

where $\mathcal{A}(q_{1},q_{2},...)$ is the adjacency matrix of the quiver,
equipped with the fugacity-weighted arrows.

The simplest physical case is the clover quiver theory, i.e. $\mathcal{N}=4$
super Yang-Mills theory, a quiver with one node associated with the
$\mathrm{U}(N)$ gauge group, and three arrows or the chiral fields in the
adjoint representation of $\mathrm{U}(N)$. Trivially, it can be generalized to
the case of the $k$ chiral fields, see Fig. 1, which is often called in the
literature by the $k$-matrix model harmonic oscillator. For example, the case
$k=2$, the 2-matrix harmonic oscillator, is considered in Aharony:2003sx . The
$k=3$ case is $\mathcal{N}=4$ SYM and $k=4$ case is the conifold theory. The
asymptotics and phase structure of the generalized clover quiver with $k$
chiral fields, are studied in ramgoolam2020quiver . The generating function of
the index of the generalized clover quiver can be obtained from Eq.(184), as
the following one-matrix model,

$\displaystyle\mathcal{I}(q_{i})=\int_{\mathrm{U}(N)}\differential{U}\exp\left[\sum_{n=1}^{\infty}\frac{1}{n}(q_{1}^{n}+q_{2}^{n}+...+q_{k}^{n})\tr
U^{n}\tr U^{-n}\right].$ (186)

This matrix model can be approximated in the large $N$ limit by the
generalized GWW model with coupling $t_{n}=\frac{1}{n}\sum_{j}q_{j}^{n}$. In
this approximation, this matrix model can be seen as a unitary matrix model
with the following potential,

$\displaystyle
V(U)=\tr\log\left(\prod_{j=1}^{k}(1-q_{j}U)(1-q_{j}U^{-1})\right).$ (187)

#### Jordan quiver

Let us first consider the univariate case of the generalized clover
$q_{j}=q_{1}\delta_{j,1}$, and change variable $q_{1}=t$ for more convenience.
This is a case of clover quiver with one arrow, and it is known as the Jordan
quiver. This is also an example of the simplest matrix model with plethystic
exponential potential. Similarly, the Jordan quiver gauge theory can be
approximated by the generalized GWW matrix model with the coupling
$t_{n}=t^{n}/n$, which leads to the matrix model with the following potential,

$\displaystyle V(U)=\tr\log(1-tU)(1-tU^{-1}).$ (188)

In this form of the coupling ($t_{n}=t^{n}/n$), as the assumption
$\alpha_{p^{\prime}}=0$ for $p^{\prime}<p$, can not be satisfied, there is
only a possibility for the $p=2$ critical dynamics, and in fact, we can only
discuss the multi-critical dynamics once we fine-tune the couplings, similar
to what we have done in the previous parts. For the fine-tuned couplings, one
can find the following generating functions for the parameters $\alpha_{p}$
and $\beta$,

$\displaystyle\alpha_{p}=\frac{2}{p!}\sum_{n=1}^{\infty}n^{p}\,t^{n}=\frac{2}{p!}\text{Li}_{-p}(t)\,,\qquad\beta=2\sum_{n=1}^{\infty}t^{n}=\text{Li}_{0}(t)=\frac{2t}{1-t}\,,$
(189)

and the parameter $\lambda_{p}$ is computed as

$\displaystyle\lambda_{p}=\alpha_{p}^{-\frac{1}{1+p}}(\beta_{c}-\beta)=2^{-\frac{1}{1+p}}\left(\frac{1-3t}{1-t}\right)\left(\frac{\text{Li}_{-p}(t)}{p!}\right)^{-\frac{1}{1+p}}.$
(190)

First notice that critical temperature defined by $\mu_{c}=-\log t_{c}$, can
be obtained from $\beta_{c}=1$, ($t_{c}=1/3$), as $\mu_{c}=\log 3$. Then, we
can compute the free energy in the weak and strong coupling phases from
Eq.(140), as

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle-\log(1-t^{2})+\mathcal{F}^{(p)}_{\text{inst}}&\text{for
}t<1/3\\\\[10.00002pt]
\displaystyle-\log(1-t^{2})+2^{-\frac{2}{p}}c_{p}\left(\frac{1-3t}{1-t}\right)^{\frac{2(p+1)}{p}}\left(\frac{\text{Li}_{-p}(t)}{p!}\right)^{-\frac{2}{p}}+\mathcal{F}^{(p)}_{\text{pert}}&\text{for
}t>1/3\end{cases},$ (191)

where $\mathcal{F}^{(p)}_{\text{inst}}$ is computed by sum over all
$n$-instantons contributions, given by

$\displaystyle\mathcal{F}^{(p)}_{n\text{-inst}}=\frac{1}{4^{n}n}\,\xi_{p}^{-n}N^{-n-2}\exp{-\frac{2np}{p+1}\xi_{p}N},$
(192)

where

$\displaystyle\xi_{p}=\lambda_{p}^{\frac{p+1}{p}}=2^{-\frac{1}{p}}\left(\frac{1-3t}{1-t}\right)^{\frac{p+1}{p}}\left(\frac{\text{Li}_{-p}(t)}{p!}\right)^{-\frac{1}{p}},$
(193)

and $\mathcal{F}^{(p)}_{\text{pert}}$ is given by

$\displaystyle\mathcal{F}_{1}^{(p)}=-\frac{c}{1+p}\log
2+c\log\left(\frac{1-3t}{1-t}\right)-\frac{c}{1+p}\log\left(\frac{\text{Li}_{-p}(t)}{p!}\right)-\frac{p}{2(p+1)}\log
N+\log C_{p},$

and

$\displaystyle\mathcal{F}_{g}^{(p)}=\frac{(-1)^{g}\,(c_{1})^{g-1}}{2^{\frac{2-2g}{p}}(g-1)}\left(\frac{1-3t}{1-t}\right)^{\frac{(p+1)}{p}(2-2g)}\left(\frac{\text{Li}_{-p}(t)}{p!}\right)^{-\frac{2-2g}{p}}.$
(195)

The leading genus zero free energy (137a) can be obtained as

$\displaystyle\tilde{\mathcal{F}}_{0}=-\log(1-t^{2})-2^{-\frac{2}{p}-1}\frac{p^{2}}{(p+1)(p+2)}\binom{p}{\frac{p}{2}}^{-\frac{2}{p}}\left(\frac{1-3t}{1-t}\right)^{\frac{2(p+1)}{p}}\left(\frac{\text{Li}_{-p}(t)}{p!}\right)^{-\frac{2}{p}}.$
(196)

Let us consider the critical case $p=2$, in which we have

$\displaystyle\alpha_{2}=\sum_{n=1}^{\infty}n^{2}\,t^{n}=\text{Li}_{-2}(t)=\frac{t(1+t)}{(1-t)^{3}}\,,\qquad\beta=\frac{2t}{1-t}\,,$
(197)

Similar results for this matrix model are obtained, from Riemann–Hilbert
method, in baik2000random ; baik2001symmetrized . The double-scaling parameter
$s$ and the genus-zero free energy in this case are obtained as

$\displaystyle\lambda_{2}=\frac{1-3t}{t^{\frac{1}{3}}(1+t)^{\frac{1}{3}}}\,,\qquad\tilde{\mathcal{F}}_{0}=-\log(1-t^{2})-\frac{(1-3t)^{3}}{12t(1+t)}.$
(198)

As an example of multi-critical dynamics for the model with the fine tuned
couplings, let us consider the case $p=4$, in which we have

$\displaystyle\alpha_{4}=\text{Li}_{-4}(t)=\frac{t(1+t)(1+t(10+t))}{12(1-t)^{5}}\,,\qquad\beta=\frac{2t}{1-t}\,.$
(199)

The parameter $\lambda_{4}$ and the genus zero free energy are obtained as

$\displaystyle\lambda_{4}=\frac{12^{1/5}(1-3t)}{\left(t(1+t)(1+t(10+t))\right)^{1/5}}\,,\quad\tilde{\mathcal{F}}_{0}=-\log(1-t^{2})+\frac{4}{15}\sqrt{2}\frac{(1-3t)^{5/2}}{\left(t(1+t)(1+t(10+t))\right)^{2}}.$
(200)

It is also interesting to consider the limit $p\to\infty$, and by using we
observe

$\displaystyle\lim_{p\to\infty}\text{Li}_{-p}(e^{-\mu})=\Gamma(1+p)\mu^{-p-1},\quad\lim_{p\to\infty}\left(\frac{\Gamma(1+p)}{p!}\right)^{-\frac{1}{1+p}}=1,$
(201)

which lead to

$\displaystyle\lambda_{\infty}=\xi_{\infty}=\frac{1-3t}{1-t}\mu,$ (202)

and the free energy,

$\displaystyle\mathcal{F}^{(\infty)}=-\log(1-t^{2})+\begin{cases}\displaystyle\mathcal{F}_{\text{inst}}^{(\infty)}&\text{for
}t<1/3\\\\[10.00002pt]
\displaystyle\frac{1}{16}\left|\frac{1-3t}{1-t}\mu\right|^{2}+\mathcal{F}_{\text{pert}}^{(\infty)}&\text{for
}t>1/3\end{cases},$ (203)

where $\mathcal{F}_{\text{inst}}^{(\infty)}$ is the sum over

$\displaystyle\mathcal{F}^{(\infty)}_{n-\text{inst}}=\frac{1}{4^{n}n}\left(\frac{1-3t}{1-t}\mu\right)^{-n}N^{-n-2}\exp{2\left(\frac{1-3t}{1-t}\right)\mu
nN},$ (204)

and $\mathcal{F}^{(\infty)}_{\text{pert}}$ is given by

$\displaystyle\mathcal{F}_{1}^{(\infty)}=c\log\left(\frac{1-3t}{1-t}\right)-\frac{1}{2}\log
N+\log
C_{p},\quad\text{and}\quad\mathcal{F}_{g}^{(\infty)}=\frac{(-1)^{g}}{g-1}c_{1}^{g-1}\left(\frac{1-3t}{1-t}\mu\right)^{2-2g}.$

#### Generalized clover quiver

Using the large $N$ limit of the generating function of the generalized clover
quiver,

$\displaystyle\lim_{N\rightarrow\infty}\mathcal{I}(q_{i})=\prod_{n=1}^{\infty}\frac{1}{1-\sum_{j=1}^{k}q_{j}^{n}}.$
(206)

As it is shown in ramgoolam2020quiver , the leading singularity of the above
generating function is at $n=1$, leading to the Hagedorn phase structure of
this model with the critical hyper-surface $\sum_{j=1}^{k}q^{*}_{j}=1$. Thus,
in our approximation, the GWW model is dominant partition function of this
model, with the parameters at the critical dynamics $p=2$,

$\displaystyle\alpha_{2}=\sum_{j=1}^{k}q_{j}\equiv
T,\qquad\beta=2\sum_{j=1}^{k}q_{j}=2T.$ (207)

Similarly, all the results in Section 4.2 follow immediately by replacing
$t_{1}$ with $T$. Moreover, straightforward generalization of the results for
Jordan quiver to the generalized clover quiver leads to

$\displaystyle\alpha_{p}=\frac{2}{p!}\sum_{j=1}^{k}\text{Li}_{-p}(q_{j})\,,\qquad\beta=\sum_{j=1}^{k}\frac{2q_{j}}{1-q_{j}}\,,$
(208)

and the parameter $\lambda_{p}$ is obtained as

$\displaystyle\lambda_{p}=2^{-\frac{1}{1+p}}(p!)^{\frac{1}{1+p}}\left(\sum_{j=1}^{k}\text{Li}_{-p}(q_{j})\right)^{-\frac{1}{1+p}}\left(1-\sum_{j=1}^{k}\frac{2q_{j}}{1-q_{j}}\right).$
(209)

Then, it is straightforward to compute $\mathcal{F}^{(p)}$,
$\mathcal{F}^{(p)}_{\text{inst}}$ and $\mathcal{F}^{(p)}_{\text{pert}}$ by
inserting $\lambda_{p}$ from (209) in Eqs.(140), (130) and (137b). For
example, the genus zero energy can be obtained as

$\displaystyle\tilde{\mathcal{F}}_{0}=\sum_{n}\frac{1}{n}\,\left(\sum_{j=1}^{k}q_{j}^{n}\right)^{2}+2^{-\frac{2}{p}-1}\frac{p^{2}p!^{\frac{2}{p}}}{(p+1)(p+2)}\binom{p}{\frac{p}{2}}^{-\frac{2}{p}}\left(\sum_{j=1}^{k}\text{Li}_{-p}(q_{j})\right)^{-\frac{2}{p}}\left(1-\sum_{j=1}^{k}\frac{2q_{j}}{1-q_{j}}\right)^{\frac{2(p+1)}{p}}.$

Moreover, in the limit $p\to\infty$, the smallest $\mu_{j}=-\log q_{j}$,
denoted by $\mu_{\text{min}}$, contributes in $\lambda_{\infty}$, as we obtain

$\displaystyle\lambda_{\infty}=\left(1-\sum_{j=1}^{k}\frac{2q_{j}}{1-q_{j}}\right)\mu_{\text{min}},$
(211)

and, using this parameter, the free energy can be computed in a similar
fashion.

Finally, in the unrefined case $q_{1}=q_{2}=...=q_{k}\equiv q$, one can
compute the parameters

$\displaystyle\alpha_{p}=\frac{2k}{p!}\text{Li}_{-p}(q)\,,\qquad\beta=\frac{2kq}{1-q}\,,\qquad
q_{c}=\frac{1}{2k+1},$ (212)

and

$\displaystyle\lambda_{p}=(2k)^{-\frac{1}{1+p}}\left(\frac{1-(2k+1)q}{1-q}\right)\left(\frac{\text{Li}_{-p}(q)}{p!}\right)^{-\frac{1}{1+p}}.$
(213)

It is straightforward to compute the free energy in this case, similar to the
Jordan quiver.

###  6.3 $\mathcal{N}=4$ superconformal index

In this part we consider the superconformal indices of the four-dimensional
gauge theories and their matrix integral representation. In general, the
$\mathcal{N}=1$ superconformal index is defined by

$\displaystyle\mathcal{I}(\mu_{i})=\Tr_{S^{3}}{(-1)^{F}e^{-\beta\delta}e^{-\mu_{i}\mathcal{M}_{i}}},$
(214)

where the trace is over the Hilbert space of the theory quantized on $S^{3}$,
$F$ is the fermion number, $\beta$ is the inverse temperature, $\delta$ is the
Hamiltonian defined by
$\delta=\frac{1}{2}\\{\mathcal{Q},\mathcal{Q}^{\dagger}\\}$ with a supercharge
$\mathcal{Q}$ and $\mathcal{M}_{i}$ are the global symmetry generators that
are annihilated by the super charge with the associated fugacity $\mu_{i}$.

In particular, we consider $4d$ $\mathcal{N}=1$ superconformal field theory on
$S^{3}\times S^{1}$. The global symmetry of the theory, i.e. the isometry of
the $S^{3}$, is $\text{Spin}(4)=\mathrm{SU}(2)_{1}\times\mathrm{SU}(2)_{2}$
with indices $\alpha=\pm 1,\dot{\alpha}=\pm 1$. In $\mathcal{N}=1$
supersymmetry in four dimensions, we have four supercharges and their
conjugates
$\mathcal{Q}_{\alpha},\mathcal{S}^{\alpha}=\mathcal{Q}^{\dagger\alpha},\tilde{\mathcal{Q}}_{\dot{\alpha}},\tilde{\mathcal{S}}^{\dot{\alpha}}=\tilde{\mathcal{Q}}^{\dagger\dot{\alpha}}$.
We choose a particular supercharge $\mathcal{Q}=\mathcal{Q}_{-}$, which
satisfies the algebra
$\\{\mathcal{Q},\mathcal{Q}^{\dagger}\\}=\Delta-2j_{1}+\frac{3}{2}r$, with
$\Delta$ is the conformal dimension, $j_{1}$ ($j_{2}$) is Cartan generator of
$\mathrm{SU}(2)_{1}$ ($\mathrm{SU}(2)_{2}$), and $r$ is the
$\mathrm{U}(1)_{r}$ R-charge. Then we can define the superconformal index,
using some fugacity factors $p$ and $q$, as

$\displaystyle\mathcal{I}(p,q)=\Tr{(-1)^{F}e^{-\beta\delta}p^{\frac{1}{3}(\Delta+j_{1})+j_{2}}q^{\frac{1}{3}(\Delta+j_{1})-j_{2}}}.$
(215)

Using the fact that the only states with $\delta=0$ contribute to the index
and thus the index is independent of $\beta$, we obtain

$\displaystyle\mathcal{I}(p,q)=\Tr{(-1)^{F}p^{j_{1}+j_{2}-\frac{1}{2}r}q^{j_{1}-j_{2}-\frac{1}{2}r}}.$
(216)

The computations of the superconformal index has two parts, first is to
compute the single letter index $i_{k}$ of each supermultiplet which is
labeled here by $k$, and then using the plethystic exponential (PE), defined
in (8), and matrix integral (6), to compute the full index,

$\displaystyle\mathcal{I}(p,q,V)=\int_{\mathrm{U}(N)}\differential{U}\prod_{k}\text{PE}\
[i_{k}(p,q,U,V)],$ (217)

where $U$ denotes the gauge group element and $V$ denotes the flavour group
element.

In $\mathcal{N}=1$ SCFT, the single letter index is the sum of the vector
multiplet index $i_{V}$ and chiral multiplet index $i_{S}$
dolan2009applications ,

$\displaystyle i=i_{V}+i_{S},$ (218)

where

$\displaystyle
i_{V}(p,q,U)=-(\frac{p}{1-p}+\frac{q}{1-q})\,\chi_{\text{adj}}(U),\quad
i_{S}(p,q,U,V)=\frac{(pq)^{\frac{1}{3}}\chi_{\bar{\mathcal{R}}}(U,V)-(pq)^{\frac{2}{3}}\chi_{\mathcal{R}}(U,V)}{(1-p)(1-q)},$

and the adjoint character $\chi_{\text{adj}}$, and the bifundamental character
$\chi_{\mathcal{R}}$, are given by

$\displaystyle\chi_{\text{adj}}(U^{n})=\tr U_{a}^{n}\tr U_{a}^{\dagger
n},\quad\chi_{\mathcal{R}_{a\bar{b}}}(U^{n})=\tr U_{a}^{n}\tr U_{b}^{\dagger
n}.$ (220)

Next, we explain the matrix integral representation of the index in the
explicit example of $\mathcal{N}=4$ SYM. In this case, using Eqs.(6) and (6.3)
the matrix integral representation of the superconformal index becomes

$\displaystyle\mathcal{I}(p,q)=\int_{\mathrm{U}(N)}\differential{U}\exp\left[\sum_{n=1}^{\infty}a_{n}(p,q)\tr
U^{n}\tr U^{-n}\right],$ (221)

where

$\displaystyle a_{n}(p,q)=\frac{i(p^{n},q^{n})}{n},\quad i(p,q)=\frac{2pq-
p-q+3(pq)^{\frac{1}{3}}-3(pq)^{\frac{2}{3}}}{(1-p)(1-q)}.$ (222)

Single letter index can be written as

$\displaystyle i(p,q)=1-\frac{(1-(pq)^{\frac{1}{3}})^{3}}{(1-p)(1-q)}.$ (223)

In the large $N$ limit, the above matrix integral (221) becomes an infinite
product,

$\displaystyle\lim_{N\to\infty}\mathcal{I}(p,q)=\prod_{n=1}^{\infty}\frac{1}{1-i(p^{n},q^{n})}.$
(224)

At weak coupling, we can effectively approximate the superconformal index
(221) by the generalized GWW model

$\displaystyle\mathcal{I}(p,q)\approx\mathcal{Z}(p,q)=\int_{\mathrm{U}(N)}\differential{U}\exp\left[\sum_{n=1}^{\infty}t_{n}(p,q)\left(\tr
U^{n}+\tr U^{-n}\right)\right],$ (225)

and in terms of generalized GWW couplings, we obtain

$\displaystyle t_{n}=a_{n}(p,q)=\frac{i(p^{n},q^{n})}{n},\quad
t_{n}=\frac{1}{n}\left(1-\frac{(1-(pq)^{\frac{n}{3}})^{3}}{(1-p^{n})(1-q^{n})}\right).$
(226)

In the rest of this section, we perform some numerical analysis of the index.

#### GWW model in $\mathcal{N}=4$ SYM

First we consider consider the truncated case $n=1$; the GWW model with the
$\mathcal{N}=4$ SYM coupling $t_{1}$ from Eq.(226). Notice that this case is
dominant in the phase structure of the model. In this case $\mathcal{Z}(p,q)$
reduces to a effective theory of GWW model,

$\displaystyle\mathcal{Z}(p,q)\sim\int_{\mathrm{U}(N)}\differential{U}\exp\left[t_{1}(p,q)\left(\tr
U+\tr U^{-1}\right)\right].$ (227)

In unrefined case $p=q\equiv x$, we have

$\displaystyle i(x)=1-\frac{(1-x^{\frac{2}{3}})^{3}}{(1-x)^{2}},\qquad
t_{n}=\frac{1}{n}\left(1-\frac{(1-x^{\frac{2n}{3}})^{3}}{(1-x^{n})^{2}}\right).$
(228)

Thus, the critical parameters can be obtained as

$\displaystyle\alpha_{2}=t_{1}=\frac{3x^{\frac{2}{3}}+4x+2x^{\frac{4}{3}}}{(1+x^{\frac{1}{3}}+x^{\frac{2}{3}})^{2}}\,,\qquad\beta=2t_{1}=\frac{6x^{\frac{2}{3}}+8x+4x^{\frac{4}{3}}}{(1+x^{\frac{1}{3}}+x^{\frac{2}{3}})^{2}}\,,\qquad
x_{c}=1,$ (229)

where we used the observation that the critical point $t_{1}^{*}=i(x_{c})=1$
implies $x_{c}=1$. Then, the double-scaling parameter
$s=\lambda_{2}N^{\frac{2}{3}}$ can be computed,

$\displaystyle\lambda_{2}=\frac{1-\beta}{\alpha_{2}^{1/3}}=\frac{1+2x^{\frac{1}{3}}-3x^{\frac{2}{3}}-6x-3x^{\frac{4}{3}}}{(1+x^{\frac{1}{3}}+x^{\frac{2}{3}})^{\frac{4}{3}}(3x^{\frac{2}{3}}+4x+2x^{\frac{4}{3}})^{\frac{1}{3}}}.$
(230)

Using the parameter $\lambda_{2}$, the free energy in different regimes can be
computed explicitly. However, as our result is comparable to the free energy
obtained by other plausible methods in the vicinity of the critical point, it
is natural to expand around the critical point $x_{c}=1$. Let us write the
fugacity factor as $x=e^{\epsilon}$, in which the parameter $\epsilon$ is
proportional to the radius of $S^{1}$ in the definition of the superconformal
index, and it can be interpreted as the inverse temperature. Thus, we have the
following high temperature expansions at small $\epsilon$,

$\displaystyle
t_{1}=i(x)=1+\frac{8}{27}\epsilon-\frac{2}{243}\epsilon^{3}+O(\epsilon^{5}),$
(231)

and

$\displaystyle\lambda_{2}=-1-\frac{40}{81}\epsilon+\frac{256}{6561}\epsilon^{2}+\frac{10606}{1594323}\epsilon^{3}+O(\epsilon^{4}).$
(232)

Using the above expansion of the parameters, the free energy can be computed
from Eq.(90),

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}+\mathcal{F}_{\text{inst}}&\text{for
}x<1\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}-\frac{1}{12}|\lambda_{2}|^{3}+\mathcal{F}_{\text{pert}}&\text{for
}x>1\end{cases},$ (233)

where we obtain

$\displaystyle\mathcal{F}_{c}=t_{1}^{2}=1+\frac{16}{27}\epsilon+\frac{64}{729}\epsilon^{2}-\frac{4}{243}\epsilon^{3}+O(\epsilon^{4}),$
(234)
$\displaystyle\mathcal{F}_{c}-\frac{1}{12}|\lambda_{2}|^{3}=\frac{13}{12}+\frac{58}{81}\epsilon+\frac{304}{2187}\epsilon^{2}-\frac{2093}{118098}\epsilon^{3}+O(\epsilon^{4}),$
(235)

and $\mathcal{F}_{\text{inst}}$ and $\mathcal{F}_{\text{pert}}$ can also be
computed explicitly, as before.

#### Generalized GWW model in $\mathcal{N}=4$ SYM

Next, we consider the generalized GWW model approximation for the
$\mathcal{N}=4$ SYM. Expanding around $x_{c}=1$, we similarly obtain

$\displaystyle
n\,t_{n}=i(x^{n})=1+\frac{8n}{27}\epsilon-\frac{2n^{3}}{243}\epsilon^{3}+O(\epsilon^{5}).$
(236)

The expansion can be written formally as

$\displaystyle n\,t_{n}=1+\sum_{l=1}^{\infty}(-1)^{l+1}a_{l}\
\epsilon^{2l-1}n^{2l-1},$ (237)

for some positive numerical coefficients $a_{l}$. Thus, the parameters
$\alpha_{p}$ and $\beta$ defined in Eq.(120), can be evaluated as

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle
2\sum_{n=1}^{\infty}\left(1+\sum_{l=1}^{\infty}(-1)^{l+1}a_{l}\
\epsilon^{2l-1}n^{2l-1}\right)=2\zeta(0)+2\sum_{l=1}^{\infty}(-1)^{l+1}a_{l}\epsilon^{2l-1}\zeta(1-2l),$
$\displaystyle\alpha_{p}$ $\displaystyle=$
$\displaystyle\frac{2}{p!}\sum_{n=1}^{\infty}\left(n^{p}+\sum_{l=1}^{\infty}(-1)^{l+1}a_{l}\
\epsilon^{2l-1}n^{p+2l-1}\right)=\frac{2}{p!}\zeta(-p)+\frac{2}{p!}\sum_{l=1}^{\infty}(-1)^{l+1}a_{l}\epsilon^{2l-1}\zeta(1-p-2l).$

To be explicit, we keep the few terms in the expansion as following

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle
2\zeta(0)+\frac{16}{27}\zeta(-1)\epsilon-\frac{4}{243}\zeta(-3)\epsilon^{3},$
$\displaystyle\alpha_{p}$ $\displaystyle=$
$\displaystyle\frac{2}{p!}\zeta(-p)+\frac{16}{27p!}\zeta(-1-p)\epsilon-\frac{4}{243p!}\zeta(-3-p)\epsilon^{3},$
(239)

and thus the parameter $\lambda_{p}$ can be obtained from
$\lambda_{p}=\alpha_{p}^{-\frac{1}{1+p}}(1-\beta)$. The bulk free energy
$\mathcal{F}_{c}$ can be computed similarly,

$\displaystyle\mathcal{F}_{c}=\sum_{n}n\,t_{n}^{2}=\sum_{l=1}^{\infty}b_{l}\zeta(-l)\epsilon^{l-1},$
(240)

and the numerical coefficients can be computed explicitly

$\displaystyle b_{1}=1,\quad b_{2}=16/27,\quad b_{3}=64/729,\quad
b_{4}=-4/243,\quad\ldots.$ (241)

One can use this data to compute explicitly the free energy in different
phases. In $p=2$ case, up to order three, we have

$\displaystyle\mathcal{F}=\begin{cases}\displaystyle\mathcal{F}_{c}&\text{for
}x<1\\\\[10.00002pt]
\displaystyle\mathcal{F}_{c}-\frac{1}{12}|\lambda_{2}|^{3}&\text{for
}x>1\end{cases},$ (242)

where

$\displaystyle\mathcal{F}_{c}=-\frac{1}{12}+\frac{8}{10935}\epsilon^{2}+O(\epsilon^{4}),$
(243)
$\displaystyle\mathcal{F}_{c}-\frac{1}{12}|\lambda_{2}|^{3}=-\frac{241}{12}-\frac{270}{\epsilon}+\frac{1745}{567}\epsilon+\frac{283363}{1377810}\epsilon^{2}-\frac{9313}{214326}\epsilon^{3}+O(\epsilon^{4}).$
(244)

##  7 Discussion and future research

In this work, we discuss some universal results for the generic matrix models,
including the double-trace matrix model, up to some model dependent
parameters, and we compute these parameters explicitly for the generalized GWW
model. The main universal results of this study in the finite but large $N$
and infinite limit of free $\mathrm{U}(N)$ gauge theory partition functions
and the indices of gauge theories include:

  * •

Universal critical/multi-critical, large and finite $N$ results such as $1/N$
expansion, genus expansion, instanton sector, etc. in generic unitary matrix
model.

  * •

Universality in the deconfinement and Hagedorn phase transition in the
infinite $N$ limit of generic unitary matrix models.

  * •

Universal multi-critical edge fluctuation and phase transitions in generic
unitary matrix models. Infinite $p$ limit of the multi-critical dynamics and
the emergence of bulk fluctuation. The interpretation of the results in terms
of the gap dynamics.

In this work, we studied the matrix models with real couplings. Moreover, the
general case of the complex coupling is of great interest. This line of
research has recently attracted a lot of interest in the theoretical physics
community because of its relations to black holes in the context of AdS/CFT,
for a review see zaffaroni2020ads  and references therein.

In this work, we considered one-matrix models and corresponding quiver gauge
theories with one node, however, the general case of the quiver gauge theories
deal with the multi-matrix model. A possible future extension of this work
would be to consider the superconformal index of quiver gauge theories and
their multi-matrix integral representation. In particular, the two-matrix
model for conifold quiver and related gauge theories would be the first step
in this direction of research.

In the study of the right tail of the TW distribution in the instanton sector,
we used the Airy function approximation to expand the Fredholm determinant. An
alternative approach would be to consider the approximation based on the
Bessel function and its variants. The possible result in this way, would be of
great interest because of its comparison with the existing results about the
instantons of the unitary matrix models in the literature.

In the multi-critical dynamics, the asymptotic analysis lacks the rigorous
results for the sub-leading corrections to the Fredholm determinant of the
higher Airy kernel. Developing techniques and results in this direction are
highly interesting.

In this work, the main focus of the asymptotic analysis is the computation of
the free energy, however, an equally important parallel direction is the
computations of the limit shape and the entropy using the spectral curve
method and the large deviation technique.

The results of this study is based on the analysis for the case of ($p\in
2\mathbb{N}$)-multi-critical dynamics. An immediate continuation of this
analysis is to consider the odd $p$ case and its applications and implications
for random partitions and gauge theories.

In this work, we studied the unitary matrix model and their phase structure. A
natural continuation of this work is to consider Hermitian matrix model. The
multi-critical analysis for the Hermitian matrix model is performed in
Claeys:2009CPAM . Possible applications of this result are in the weakly
coupled gauge theories on compact manifolds Aharony:2003sx  and their phase
structure. Another interesting model is the Douglas–Kazakov matrix model with
quadratic potential and its instanton sector. A straightforward generalization
of our results would shed light to this classic matrix model. In this work we
considered the same coefficients for the GWW and the double-trace matrix
model, which is plausible for the effective theory at weak coupling. However,
one can study the relations between the GWW and double-trace model and their
phase transition using the Legendre transform and saddle-point analysis. One
can use the Hubbard–Stratonovich transformation and apply our results for the
generalized GWW model to the double-trace matrix integrals and their gravity
duals. Asymptotic analysis of the double-trace matrix models with real and
complex couplings is an interesting direction and remains for future. In this
way, based on the rigorous mathematical results, we hope to provide new
understandings for the phase structure associated to the double-trace models
such as deconfinement, Hagedorn and Hawking–Page phase transitions
Copetti:2020dil  and their multi-critical generalizations. Especially,
possible interpretations of the multi-critical phase structure for the
Hagedorn and Hawking-Page transitions would be interesting.

In this regard, the GWW model and its generalization seems to explain the
transition between highly excited string states and black holes, whereas the
double-trace model explains the Hawking-Page transition between the thermal
AdS and black holes alvarez2006black . Thus, our results about the critical
dynamics of the generalized GWW model studied in this paper perhaps directly
explore some universality features of the former transition and deserves
further studies. However, the highly interesting question of the possible
interpretations and implications of the multi-critical dynamics in the gravity
side is not clear at this moment, which would require new ideas and
techniques. Regarding the Hawking-Page transition, we should apply the
Hubbard-Stratonovich transformation to discuss all possible implications of
the critical and multi-critical dynamics of the double-trace model in the
gravity side. We will return to this issue in a future work.

#### Acknowledgments

This work has been supported in part by “Investissements d’Avenir” program,
Project ISITE-BFC (No. ANR-15-IDEX-0003), EIPHI Graduate School (No.
ANR-17-EURE-0002), and Bourgogne-Franche-Comté region.

## References

  * (1) B. Feng, A. Hanany and Y.-H. He, _Counting gauge invariants: The Plethystic program_ , _JHEP_ 03 (2007) 090 [hep-th/0701063].
  * (2) S. Benvenuti, B. Feng, A. Hanany and Y.-H. He, _Counting BPS Operators in Gauge Theories: Quivers, Syzygies and Plethystics_ , _JHEP_ 11 (2007) 050 [hep-th/0608050].
  * (3) B.S. Skagerstam, _On the Large $N_{c}$ Limit of the SU($N_{c}$) Color Quark-Gluon Partition Function_, _Z. Phys. C_ 24 (1984) 97.
  * (4) B. Sundborg, _The Hagedorn transition, deconfinement and $\mathcal{N}=4$ SYM theory_, _Nucl. Phys. B_ 573 (2000) 349 [hep-th/9908001].
  * (5) O. Aharony, J. Marsano, S. Minwalla, K. Papadodimas and M. Van Raamsdonk, _The Hagedorn/Deconfinement Phase Transition in Weakly Coupled Large $N$ Gauge Theories_, _Adv. Theor. Math. Phys._ 8 (2004) 603 [hep-th/0310285].
  * (6) D.J. Gross and E. Witten, _Possible Third Order Phase Transition in the Large $N$ Lattice Gauge Theory_, _Phys. Rev._ D21 (1980) 446.
  * (7) S.R. Wadia, _$N=\infty$ Phase Transition in a Class of Exactly Soluble Model Lattice Gauge Theories_, _Phys. Lett._ 93B (1980) 403.
  * (8) A. Zaffaroni, _AdS black holes, holography and localization_ , _Living Rev. Rel._ 23 (2020) 2 [1902.07176].
  * (9) A. Borodin and A. Okounkov, _A Fredholm determinant formula for Toeplitz determinants_ , _Int. Eq. Op. Th._ 37 (2000) 386 [math/9907165].
  * (10) L. Santilli and M. Tierz, _Exact equivalences and phase discrepancies between random matrix ensembles_ , _Journal of Statistical Mechanics: Theory and Experiment_ 2020 (2020) 083107.
  * (11) V. Periwal and D. Shevitz, _Exactly Solvable Unitary Matrix Models: Multicritical Potentials and Correlations_ , _Nucl. Phys. B_ 344 (1990) 731.