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

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 [1,2].
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 [3][4][5], 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 [6,7] 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 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 [8] 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 [9], 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 singletrace 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 λ 1 , the summation over the random partition is rewritten, modulo a normalization factor, as the unitary matrix integral [9], for definitions see Section 3.1. An unrefined version of the partition function (1.1) is recently studied in [10]. Defining the free energy by F = lim N →∞ N −2 log Z N , we show that the right/left (±) edge fluctuation contributes to the free energy, up to a scaling factor N −2 , as and F p is the higher-order Tracy-Widom distribution [11][12][13][14], and α p (q 1 , q 2 , ...), β(q 1 , q 2 , ...) are some model-dependent parameters and are explicitly expressed in terms of the couplings [15]. The main result of our study is that the matrix model (1.1), undergoes a multi-critical phase transition at the critical point β = β c . This phase transitions can be explained from the asymptotic behavior of the higher-order Tracy-Widom distribution, The above leading asymptotic behavior of the higher-order TW distribution and the following multi-critical double-scaling parameter, s = α , imply that the free energy in two regimes s → +∞ (β < β c ), and s → −∞ (β > β c ), modulo an additive constant, is given by (1.4) The above free energy implies a discontinuity in the (2(p + 1)/p)-th derivative of the free energy at β = β 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 2(p+1) p = 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 [6,7], see also [16]. For higher p > 2, the fluctuations at scale N − 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 = ∞ 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.

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.

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 [1,2], and then integrating over the Haar measure of the gauge group to project onto the gauge-singlets [3][4][5], known as Weyl integration formula, where g ∈ G, and in our study G is to be considered as the unitary gauge group U(N ) or their product, and χ R (g) is the character in R-representation. In Section 6, we discuss some concrete examples in four-dimensional N = 1 SCFT.

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 = 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 χ adj (g) = χ F (g)χ * F (g), the group integral recasts to a matrix integral, where by using the Weyl parameterization, the integral over the Haar measure can be written as an integral over the eigenvalues and the plethystic exponential is defined by 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 N = 1 SCFT and ii) the BPS operators that are annihilated by one supercharge in N = 1 SCFT, in particular the sixteenth BPS operators in N = 4 supersymmetric Yang-Mills (SYM) theory. Let us define a closely related model, the generalized GWW matrix model, Notice that, in general, the 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, Z(q i ) reduces to a effective theory of GWW model, The gauge theory indices as the double-trace matrix integral (2.2) 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 [17,18]. In this study, we focus on the former approximation and leave the latter relation for the future studies. The doubletrace matrix model and the generalized GWW model share equilibrium conditions. In fact, in the large N limit, I(q i ) and Z(q i ) are equivalent, upon identification and they have a similar phase structure [5]. Moreover, the truncated case n = 1 of the double matrix integral (2.2), and the GWW model, are dominant in the asymptotic analysis of the generalized models with higher degrees, n > 1, see [5] 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 [19], 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.

Integrable operator formalism in matrix models
In this part, we review, without details of the proofs, the results obtained in recent works [15,20] 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.

Schur measure random partition and generalized GWW model
Let Y be a set of partitions, and X = (x i ) i∈N be the set of the parameters and the Miwa variables defined by Then, we define our main object, the partition function of the constrained partitions, via the Schur function s λ , as where Z is the normalization of the Schur measure random partition and it is given by, We will elaborate on the physical interpretation of the partition functions Z and Z in Section 3.2. Using the character expansion, the Schur partition can be written in terms of the unitary matrix integral, Let us define the Fredholm determinant on a specific domain I ⊂ R, 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 (3.5) [9], where the domain I is defined as I = (N + 1 2 , ∞), and K(k, l) is given by [21], The J(z) is called the wave function and it has the following mode expansion, Plancherel measure random partition and GWW model A particular example of the Schur measure with a single parameter is called the Plancherel measure with the following partition function, where dim λ is the dimension of the irreducible representation parameterized by λ of the symmetric group S ∞ , and |λ| = ∞ i=1 λ i , and Z = e q 2 . The Plancherel partition function with the constraint λ 1 ≤ N has the matrix integral and Fredholm determinant representation, where the discrete Bessel kernel is obtained in [22], The wave function is then given by the Bessel function J(x) = J x (2q), with the generating function

Difference equation for wave functions
We observe that the wave functions of the Schur partitions satisfy the following differential and difference equations: where we define the shift operator ∇ . In the case of the Plancherel measure, the wave function obeys the difference equation, a special case of (3.14) , We remark that, as mentioned in [15], the differential/difference equation (3.14) is interpreted as a quantization of the spectral curve for the generalized GWW model,

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 U(N ) matrix model is defined by 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 where F c and 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.(3.3), 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 ∼ 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, → 0, of the random partition. Moreover, as we will explicitly explain later, the parameter s in the domain of the Fredholm kernel (3.6), 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, q) → (x, q/ ), the scaling limit of the difference equation (3.15) becomes which has the plane wave solutions . (3.21) Then, the scaling limit of the discrete Bessel kernel (3.12) becomes the sine kernel 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 [23].
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 The asymptotics of the Fredholm determinant as s → ∞ is obtained in [24,25], This asymptotic result leads to the computation of the free energy in the bulk, the scaling parameter s = γN , as where F c is given by the normalization of the Schur measure random partition in Eq. (3.19). 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 to re-scale the parameters and expand the shift operator, where the coefficients are defined as (3.28) In the scaling limit → 0, and by keeping α p = 0 for p < p, the difference equation (3.27) becomes the following differential equation, Thus, we observe that the scaling limit of the wave function is given by the p-Airy function, where C is an integral contour providing a convergent integral. The case p = 2 corresponds to the standard Airy function. The p-Airy function, Ai p (−1) p 2 +1 ξ , satisfies the generalized Airy equation, (3.31) In fact, we can analyze the semi-classical behavior of the p-Airy function based on the reduced spectral curve, This is obtained from the spectral curve of the GWW model (3.16) in the scaling limit discussed above, i.e., parametrize w = exp ( y), then take the limit → 0 with tuning the parameters. We remark that the p-Airy spectral curve (3.32) agrees with (A p−1 , A 0 )-type Argyres-Douglas theory [26,27]. Next, we discuss the scaling limit of the kernel. Using the Eq.(3.7), we obtain the kernel in terms of the wave functions, and then using the scaling limit of the wave function, one can obtain the scaling limit of the kernel as where the p-Airy kernel is defined as and Ai (r) p (x) is the r-th derivative of the Airy function.

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, The probability distribution function of the largest eigenvalue and its statistical behavior in the scaling limit is given by the Tracy-Widom distribution [28]. Random partition, as a discrete analog of the random matrix [9], 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 [15], where F p (s) is a higher-order analog of the Tracy-Widom distribution [11][12][13][14]. The probability distribution (3.37) can be used to compute the contribution of the edge fluctuation to the partition function of the generalized GWW model by using the Eqs. (3.2) and (3.6). In fact, there are three regions in the vicinity of the edge, namely the finitely (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 (±) sides of the edge contribute to the free energy (3.18), up to a scaling factor N −2 , as Therefore, the right/left (±) edge free energy are obtained from the left and right tails of the higher TW distribution, In the special case of the Plancherel partition, using the scaled parameters (x, q) → (x/ , q/ ), the difference equation (3.15) becomes and in the scaling limit, the wave function becomes the Airy function, . Similarly, the scaling limit of the discrete Bessel kernel is given by the Airy kernel, 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 (3.43)

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.

General unitary matrix model and TW distribution
Precisely speaking, in the particular case p = 2 of the Eq.(3.37), and by assuming the scaling relation between 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, where c 1 and c 2 are some model-dependent parameters which in the case of generalized GWW model we have c 1 = β(t 1 , t 2 , ...) and c 2 = α 2 (t 1 , t 2 , ...) in Eq.(3.28), and F 2 is the TW distribution and Using the definition (3.17), and by fixing the 't Hooft parameter γ := N in the large N limit of the U(N ) matrix model, the universal result for the free energy can be obtained from (3.2), (3.18), and (4.1), up to a scaling factor N −2 , as and the free energy in the left side of the edge is given by and right sides of the edge is given by 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 where q(x) is the solution of the Painlevé II equation, and q xx denotes the second derivative of q with respect to x. The asymptotic behavior of the solution is obtained in [29,30], and as x → −∞, and Using the Hastings-McLeod results [30] we have Moreover, they obtained 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 [31], and the following result is obtained, where ζ(x) is the Riemann zeta function. The asymmetry of the tails of the TW distribution is apparent in the leading order asymptotic, (4.14) 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 (4.12), .

(4.15)
We can add the crossing region, between the two tails, to the above result and after expanding the logarithm, we obtain where c = log 2 24 + ζ (−1). We can recast the expansions as where in the right tail, the first line, we have and in the left tail, third line, and The above result reproduces the known results about the scaling behavior of the free energy of GWW model [17] and moreover generalize it to the generalized GWW model [18]. 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 → ∞, thus we can expand the exponential in Eq.(4.6), as the integrand becomes small in this limit, then, using the asymptotic result (4.22) and the Airy function asymptotics, by keeping the first term in the above expansion of the Airy function, and using it in Eq.(4.21), we obtain The above integral can be evaluated as Using the asymptotic expansion of the generalized exponential function at x → ∞, we obtain The above right tail expansion, s → ∞ of the TW distribution, matches, up to some numerical coefficients, with the similar expansion in the literature, for example in [31,32],

Noperturbative instanton sector
In the region β > β c , (β − β c ) = O(1), using the parameterization s = (ξN ) 2 3 , we observe that there is no perturbative corrections to the continuum free energy and the nonperturbative corrections due to the instantons are of order O(e −N ), and we can expand the leading terms of Eq.(4.29) to obtain the instanton contribution to the free energy, where n is the number of the instantons. Thus, we can write F inst = ∞ n=1 F n-inst , and in particular, the one-instanton sector, which is the dominant contribution, is

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 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 = λN 2 3 , in which the two parameters are related by λ = ξ 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 (4.32), as Furthermore, as we observed, there is no perturbative corrections O( 1 N ) in the region β > β c , and the perturbative corrections in the region β < β c , can be computed by using the left tail asymptotics of the TW (2) distribution. In fact, the free energy at genus one, and for higher genus (g ≥ 2), are given by Let us define the perturbative corrections to free energy by and then by putting the perturbative and non-perturbative leading and subleading terms together in the two phases, we obtain Notice that, the above result for the free energy is universal, besides the model-dependent parts, namely the continuum free energy F c , and parameters ξ and λ.

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 β > β c ) of the TW distribution, is obtained via the Riemann-Hilbert analysis [31], as s → −∞, Thus, we obtain the universal finite contribution from the edge fluctuation to the free energy in strongly-coupled phase, , and the parameters α and β depend on the model, and 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 It is easy to observe that the third derivative of the free energy is discontinuous at β(t i ) = γ ≡ β 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 β = β c between the two phases: weak coupling β < β c (s → ∞) and strong coupling β > β c (s → −∞). 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, where s is finite, s ∼ α . Notice that, from Eq.(4.6), one can observe that q satisfies 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 → −∞, the asymptotic behavior is given by the Eq.(4.8) and in the limit x → ∞, it is q (x) → Ai(x). In the crossing region, it is numerically obtained that the TW distribution can be approximated by the Gamma distribution [33], where γ(k, x) is the lower incomplete gamma function defined as 43) and the parameters k, α and θ 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.

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 and λ = α −1/3 2 (β c − β), with the coefficients α 2 and β, given by Assuming γ = 1, the critical point β c = 1 implies that the critical coupling is t * 1 = 1/2. Then, we easily obtain and we can compute the free energy (4.36), where F inst in Eq.(4.30) can be computed by sum over all n-instantons contributions, given by and F pert is given by and (4.50) To compare with results in the literature for example [17,34], 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 . Let us consider the 1/N expansion in GWW model [17], where satisfies Painlevé II equation [17]. We observe an agreement between G (1,2,3) n 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 (β) in Eq.(4.18) matches with f (t) around t = 1 up to the first term in the expansion which is of order O( 3/2 ). We hope to comeback to this problem, in our future studies.
The perturbative genus expansion in the GWW model [34] is where B 2g denotes Bernoulli number. In particular, at genus one, we have In general, for any g ≥ 1, we expect that the above results match with Eqs. c (g) n t n ∼ − (4.56) More precise comparison and further implications of the above equation remain for future studies.

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 (4.44) with explicit model-dependent parameters given by (4.58) 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 [18], , (4.59) where F (1,2,3) n and H(t n ) can be computed using the methods of orthogonal polynomials [35]. In the vicinity of the critical point, assuming γ = 1, 0 . Similar to GWW case, they can compared with G 1,2,3 n in Eqs. (4.18) and (4.20) with insertion of parameters from Eq.(4.57). In particular, the genus zero free energy in the strong coupling regime is expected to be obtained from Eq.(4.58), as 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 δ n,m , with the following potential (4.62) For simplicity we write t m ≡ t and then in this case we obtain the double scaling parameter and the free energy, This implies that there is a third-order phase transition at critical coupling t * = 1 2m . Similarly, F inst and F pert can be computed explicitly.

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 TW (p) . The asymptotic behavior of the TW distribution (p = 2) is studied in [31], and recent asymptotic results for the Pearcey processes (p = 3) are obtained in [36]. See also earlier results [37,38]. The asymptotic analysis for the generic case of higher TW (p) distribution is conjectured in [13] and proved in [14]. See also [11,12]. 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.

Multi-critical dynamics and higher 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 ∈ 2N, however we expect similar results for odd p [39]. 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 where c (p) 1 and c (p) 2 are some model dependent parameters, and F p (s) is TW (p) distribution. In the generalized GWW model, we have the distribution (3.37) with parameters are explicitly given by The free energy in the large N limit of the U(N ) matrix model is defined in (3.17), and it can be obtained from (3.2) and (5.1), by fixing γ := N , as is called the multi-critical double-scaling parameter.

Multi-critical double-scaling parameter
Let us define the critical point γ = β c and by using (5.2), write double scaling parameter in generalized GWW model as The multi-critical dynamics and phase structure is governed by the behaviour of free energy around the critical point β = β 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 α p and β in the asymptotic expansion of F p (s) around the critical point. For real positive t n , we have α p > 0, β > 0, and as N → ∞, the double-scaling parameter s is tending to plus/minus infinity depending on the sign of (β c − β) and thus there are two phases β c > β and β c < β 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 → ±∞.

Multi-critical instantons
In this Section, we study the right tail asymptotics of TW (p) and its applications for the instantons in the multi-critical model. In the limit s → ∞, one can obtain the asymptotic expansion of F p (s), in a straightforward calculations, by using the definition and the following asymptotic result [14], where Ai p (s) is the p-Airy function defined in (3.30). The asymptotic behavior of the p-Airy function is obtained in [40], and the following results, up to some unkown numerical factor a p andã p , holds The factorã 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 (5.7), from Eq.(5.5) in the limit s → ∞, we obtain Then, by expanding the generalized exponential integral Eq.(4.27), we obtain Finally, for generic p, we have the following right tail expansion of the logarithm of the higher TW (p) , (5.10) The right tail of TW (p) , explain the regime β < β c , and (β − β c ) = O(1) of the generalized GWW model. Using the parameterization s = (ξ p N ) p p+1 , and the above asymptotic results, the free energy in this phase is obtained as In the generalized GWW model, we have F c = n n t 2 n , and ξ p = α

Multi-critical genus expansions
In this part, we study the left tail asymptotics of 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 TW (p) distribution, in the limit s → −∞, is recently obtained in [14], where c p = − p 2 2(p + 1)(p + 2) and c = − 1 8 if p = 2, c = − 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.(4.37), we can conjecture, , we obtain the free energy in the strong coupling phase (β > β c ), as Furthermore, in the phase β > β c , and (β − β c ) = O(1), we can re-arrange the above result, and obtain the genus expansion of the free energy of the generalized GWW model, where the genus zero part and the higher genus subleading perturbative corrections arẽ with the genus one free energy is given by For higher genus g ≥ 2, by using the expansion of the conjecture, one can show that the higher than two genus free energy is given by (5.15), 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 where F

Cross-over region
In this part, we discuss the smooth cross-over region between the weak and strong coupling phases . The free energy of the generalized GWW multi-critical model in this region can be obtained from Eq.(5.5), as and q(s) is the solution of the p 2 -th member of the Painlevé II hierarchy [14], see also [41], where q s is the derivative of q w.r.t. s, and operators L n are Lenard operators defined by Moreover, the multi-critical specific heat satisfies (5.26) and in the strong coupling phase, as s → −∞, it has the following expansion [14], 27) and, in the weak coupling phase, as s → ∞, we have q((−1) p 2 +1 s) → Ai p (s). One can study the interpolation of q(s) between these two limits and extract the finite s results.

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, , (5.28) where F c = n n t 2 n , and it is straightforward to observe that G and in the strong coupling phase, for the genus expansion of the free energy we have and (5.29d)

Infinite limit (p → ∞) 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 → ∞. In this limit the double scaling parameter (5.4) becomes where we used the fact that the limit p → ∞ in definition of α p in Eq. (5.2), implies that only the term with n ∼ p survives, and we have It is interesting to compute the multi-critical 1/N expansion in this limit, where one can directly compute the limit p → ∞ of the functions in Eqs. (5.29), as In this limit, the free energy expansion in the strong coupling phase β > β c , becomes (5.34) and the rest of the subleading terms are The leading contribution to free energy indicates a second order phase transition. Moreover, the higher genus (g ≥ 2) free energy becomes In the weak coupling phase β < β c , using the asymptotic behavior of the higher Airy function in Eq.(5.7), and the infinite limit of the TW (p) in Eq.(5.8) is where Γ(t, x) is the upper incomplete Gamma function defined as  Finally, the infinite limit of the multi-critical instanton free energy can be obtained from Eq.(5.12), as 42) and the one-instanton sector is

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 β c = γ, meaning that the free energy has different expansion at left and right side of the critical point, where F c = n n t 2 n , The leading order of the free energy in two phases, imply that the order of the multi-critical phase transition is 2(p+1) p . Often, in the literature this multi-critical phase transition, for any finite p, is considered a transition of order 2(p+1) p = 3.

Multi-critical gap dynamics
Before discussing the interpretations of the multi-critical phase transition, let us briefly explain some related issues. The assumption α p = 0 for p < p, which is used to derive the coefficients (5.2) in Eq.(5.1), implies some physical constraints on the dynamics of the model. Notice that at any fixed p, the interaction and fluctuation are at scale N − p p+1 , and thus in order to see the sub-dominant fluctuations for higher p, all the larger scale interactions at p < 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 ≥ 2, couplings in the large N limit, such as t (p) n ∼ q (p−p )N , for q < 1 and 2 ≤ p ≤ 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, the critical dynamics is causing the transition between l-cut solution and (l −1)-cut solution for 1 ≤ l ≤ 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 ≤ r ≤ l ≤ 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 β c , which depends only on the potential. Let us isolate a particular set of phase transitions, namely the transitions between the p 2 -cut solutions to full support (zero gap) solution, in which p-end points of the cuts collide simultaneously, for p ≤ 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 ≤ p ≤ 2m, and each different p-(multi)-critical dynamics is governing the opening/closing of p 2 -gaps as a result of the edge fluctuation of order N − 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 − p p+1 by putting α p = 0, for p < p. We did compute the free energy of the generalized GWW model near the critical point β c with the p-end points of the cuts are colliding simultaneously.

A new second-order phase transition
In the limit p → ∞, the asymptotic behavior of the Fredholm determinant (5.13) is where the scaling parameter is s ∼ (β c − β)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-multicritical 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.(5.47), 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 multicritical dynamics. First, consider the degree two polynomial potential with the couplings, t n = t 1 δ n,1 + t 2 δ n,2 . (5.48) In this model, the multi-critical parameter is restricted to p ≤ 4, and the parameters can be easily computed The critical point β 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.(5.17), as 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 and the first few terms of the free energy are obtained as (5.52) In the multi-critical case p = 4, in principle we can fine tune the couplings t 1 , and t 2 such that α 2 → 0, and similarly one can compute, and the leading orders in free energy as 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.

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 TW (p) , to the free energy in Section 5.

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 I in Eq.(2.2), in the large N limit, using the Gaussian matrix integration, becomes an infinite product [5], . (6.1) 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 On the other hand, in the generalized GWW model, the parameters α and β 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 3) The strong-weak coupling or the deconfinement phase transition in this model, which is a generalization of the Gross-Witten phase transition, happens at β c (q i ) = γ = 1, between the two phases of weak-coupling β(q i ) < 1, and strong-coupling β(q i ) > 1.
Using the couplings t n = f (q n i )/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 Then, by expanding the Hagedorn critical hyper-surface (6.2) and comparing it with deconfinement critical hyper-surface (6.4), 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 ) 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 doubletrace 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 β 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 singletrace 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 * n i ), see Eq.(6.4). 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 doubletrace 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 * n i ). 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 (p) n = q (p−p )N /n s which, in the limit N → ∞, is zero for 2 ≤ p < p, and t n = 1/n s , for p = p. This example corresponds to the matrix model with the following poly-logarithmic potential V (U ) = n t n (tr U n + tr U −n ) = Li s (e iθ ) + Li s (e −iθ ) . (6.5) It is easy to compute the parameters of this model in terms of the Riemann zeta function, However, as the parameter β is a coupling independent number, there is no phase transition in this model. The free energy of the model is the continuum free energy In particular, the especial case of s = 1 is a non-interacting theory with zero potential, with α p = 0 for all p ∈ 2N. Moreover, the free energy is divergent at s = 1 as ζ(1) → ∞. Next, we consider a model with an arbitrary degree single term potential, where q < 1, and for all p < p, we have t k = 0, at large N limit. In this model, after fine-tuning the coupling for any given p such that α p = 0, for p < p, the parameters of the t k Figure 1. From left to right: Jordan quiver, clover quiver, and generalized clover quiver model can be computed as , (6.9) and the genus zero free energy in the strong coupling regime is Since p ≤ 2k, the infinite multi-critical limit, p → ∞, implies k → ∞ and one can assume p = 2k. In this limit although α 2 and β are either zero or diverging depending on |t|, but λ ∞ is finite and one can show λ ∞ = 2 e √ t and thus the free energy becomes (6.11)

Free chiral ring index
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 N = 1 superconformal quiver theories [42]. 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 U(N ) factors of each node i of quiver, The generating function of the BPS operators in the free chiral ring is represented by the following quiver multi-matrix integral, where the q ab;α are the fugacity factors associated with the chiral fields {φ ab;α : α ∈ {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 ,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, , (6.13) where A(q 1 , q 2 , ...) is the adjacency matrix of the quiver, equipped with the fugacityweighted arrows. The simplest physical case is the clover quiver theory, i.e. N = 4 super Yang-Mills theory, a quiver with one node associated with the U(N ) gauge group, and three arrows or the chiral fields in the adjoint representation of 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 [5]. The k = 3 case is 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 [43]. The generating function of the index of the generalized clover quiver can be obtained from Eq.(6.12), as the following one-matrix model, 14) This matrix model can be approximated in the large N limit by the generalized GWW model with coupling t n = 1 n j q n j . In this approximation, this matrix model can be seen as a unitary matrix model with the following potential,

Jordan quiver
Let us first consider the univariate case of the generalized clover q j = q 1 δ 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, In this form of the coupling (t n = t n /n), as the assumption α p = 0 for p < 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 α p and β, (6.17) and the parameter λ p is computed as . (6.18) First notice that critical temperature defined by µ c = − log t c , can be obtained from β c = 1, (t c = 1/3), as µ c = log 3. Then, we can compute the free energy in the weak and strong coupling phases from Eq.(5.22), as inst is computed by sum over all n-instantons contributions, given by 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 The parameter λ 4 and the genus zero free energy are obtained as (6.28) It is also interesting to consider the limit p → ∞, and by using we observe which lead to 30) and the free energy, inst is the sum over and F (∞) pert is given by Generalized clover quiver Using the large N limit of the generating function of the generalized clover quiver, As it is shown in [43], 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 k j=1 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, 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 36) and the parameter λ p is obtained as Then, it is straightforward to compute F (p) , F inst and F (p) pert by inserting λ p from (6.37) in Eqs.(5.22), (5.12) and (5.19b). For example, the genus zero energy can be obtained as (6.38) Moreover, in the limit p → ∞, the smallest µ j = − log q j , denoted by µ min , contributes in λ ∞ , as we obtain 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 ≡ q, one can compute the parameters and It is straightforward to compute the free energy in this case, similar to the Jordan quiver.

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 N = 1 superconformal index is defined by where the trace is over the Hilbert space of the theory quantized on S 3 , F is the fermion number, β is the inverse temperature, δ is the Hamiltonian defined by δ = 1 2 {Q, Q † } with a supercharge Q and M i are the global symmetry generators that are annihilated by the super charge with the associated fugacity µ i .
In particular, we consider 4d N = 1 superconformal field theory on S 3 × S 1 . The global symmetry of the theory, i.e. the isometry of the S 3 , is Spin(4) = SU(2) 1 × SU(2) 2 with indices α = ±1,α = ±1. In N = 1 supersymmetry in four dimensions, we have four supercharges and their conjugates Q α , S α = Q †α ,Qα,Sα =Q †α . We choose a particular supercharge Q = Q − , which satisfies the algebra {Q, Q † } = ∆ − 2j 1 + 3 2 r, with ∆ is the conformal dimension, j 1 (j 2 ) is Cartan generator of SU(2) 1 (SU(2) 2 ), and r is the U(1) r R-charge. Then we can define the superconformal index, using some fugacity factors p and q, as Using the fact that the only states with δ = 0 contribute to the index and thus the index is independent of β, we obtain I(p, q) = Tr (−1) F p j 1 +j 2 − 1 2 r q j 1 −j 2 − 1 2 r . (6.44) 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 (2.4), and matrix integral (2.2), to compute the full index, where U denotes the gauge group element and V denotes the flavour group element. In N = 1 SCFT, the single letter index is the sum of the vector multiplet index i V and chiral multiplet index i S [46], 47) and the adjoint character χ adj , and the bifundamental character χ R , are given by χ adj (U n ) = tr U n a tr U †n a , χ R ab (U n ) = tr U n a tr U †n b . (6.48) Next, we explain the matrix integral representation of the index in the explicit example of N = 4 SYM. In this case, using Eqs.(2.2) and (6.47) the matrix integral representation of the superconformal index becomes dU exp ∞ n=1 a n (p, q) tr U n tr U −n , (6.49) where a n (p, q) = i(p n , q n ) n , i(p, q) = 2pq − p − q + 3(pq) Single letter index can be written as . (6.51) In the large N limit, the above matrix integral (6.49) becomes an infinite product, At weak coupling, we can effectively approximate the superconformal index (6.49) by the generalized GWW model t n (p, q) tr U n + tr U −n , (6.53) and in terms of generalized GWW couplings, we obtain t n = a n (p, q) = i(p n , q n ) n , t n = 1 In the rest of this section, we perform some numerical analysis of the index.

GWW model in N = 4 SYM
First we consider consider the truncated case n = 1; the GWW model with the N = 4 SYM coupling t 1 from Eq.(6.54). Notice that this case is dominant in the phase structure of the model. In this case Z(p, q) reduces to a effective theory of GWW model, dU exp t 1 (p, q) tr U + tr U −1 . (6.55) In unrefined case p = q ≡ x, we have (6.56) Thus, the critical parameters can be obtained as 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 = λ 2 N 2 3 can be computed, Using the parameter λ 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 , in which the parameter 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 , and F inst and F pert can also be computed explicitly, as before.
and thus the parameter λ p can be obtained from λ p = α  One can use this data to compute explicitly the free energy in different phases. In p = 2 case, up to order three, we have

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 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 [8] 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 ∈ 2N)-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 [12]. Possible applications of this result are in the weakly coupled gauge theories on compact manifolds [5] 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 [47] 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 [18]. 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.