Phase space distribution for two-gap solution in unitary matrix model

We analyze the dynamics of weakly coupled finite temperature U(N ) gauge theories on S3 by studying a class of effective unitary matrix model. Solving Dyson-Schwinger equation at large N , we find that different phases of gauge theories are characterized by gaps in eigenvalue distribution over a unit circle. In particular, we obtain no-gap, one-gap and two-gap solutions at large N for a class of matrix model we are considering. The same effective matrix model can equivalently be written as a sum over representations (or Young diagrams) of unitary group. We show that at large N , Young diagrams corresponding to different phases can be classified in terms of discontinuities in number of boxes in two consecutive rows. More precisely, the representation, where there is no discontinuity, corresponds to no-gap and one-gap solution, where as, a diagram with one discontinuity corresponds to two-gap phase, mentioned above. This observation allows us to write a one to one relation between eigenvalue distribution function and Young tableaux distribution function for each saddle point, in particular for two-gap solution. We find that all the saddle points can be described in terms of free fermions with a phase space distribution for no-gap, one-gap and two-gap phases.


JHEP04(2016)104
1 Introduction and summary The large N limit of SU(N ) gauge theories has many qualitative features similar to QCD, primarily confinement-deconfinement transitions. Moreover, from the AdS/CFT perspective, large N gauge theories are conjectured to be dual to certain string theories on an AdS background. Different phases of strongly coupled large N theories correspond to different phases of asymptotically AdS black holes. For example, confinement/deconfinement transitions in large N gauge theories on sphere are identified with gravitational phase transitions between AdS black hole and global AdS (known as Hawking-Page phase transition). This conjecture has opened up a new avenue to study different thermodynamic properties and phase diagrams of large N Yang Mills theories on a compact manifold.
There are different ways to study non-abelian gauge theories on a compact manifold. In [1], authors argued that the thermal partition function of a gauge theory with gauge group SU(N ) on a compact manifold can be written in terms of a unitary matrix model. [2] was the first to obtain the unitary matrix model of the finite temperature partition function counting the number of states in free theory. To be precise, let us consider an Euclidean gauge theory on S 3 with an Euclidean time direction S 1 . S 1 has a circumference β which is inversely proportional to the temperature of the theory. In zero 'tHooft coupling limit (λ = g 2 Y M N → 0) all the modes (expanding the fields in spherical harmonics) become massive, with mass scale proportional to inverse of the radius of S 3 , except the time component of the gauge field, α is strongly interacting even in zero coupling limit. One can integrate out all the massive modes and the theory can be described by an effective action written purely in terms of α.
In fact the actual variable turns out to be a N × N unitary matrix Finally, the thermal partition function in terms of holonomy matrix U is given by,

JHEP04(2016)104
where [dU ] is the Haar measure and the coefficients a n (β) are given by a n (β) = z B (x n ) + (−1) n+1 z F (x n ), x = e −β . (1.4) Here z B (x) and z F (x) are single particle partition functions of the bosonic and fermionic modes respectively. They completely capture the field content of the gauge theory.
Moving to weak coupling limit, one can write a more generic matrix model in the weak 't Hooft coupling limit ( [1] for details): where, TrU n i (1. 6) with the integers n i obeying i n i = 0 and the coefficients a n , is a term with k traces, making their first appearance at (k − 1) loops in perturbation theory and consequently having a planar contribution starting with λ k−2 .
One can analyze this partition function in the large N limit and obtain different phases of the system. The main goal of this paper is to find phase space distribution corresponding to each saddle point or phase. Before we elaborate the importance of phase space distribution (see section 1.3), we briefly discuss how one obtains a phase space distribution corresponding to different saddle points.

Phase diagram from integral equation: eigenvalue analysis
The Haar measure and the effective action are invariant under unitary transformation hence one can go to a diagonal basis where U → e iθn δ nm . (1.7) In this basis the Haar measure and the traces can be written as (1.8) In N → ∞ limit the partition function can be written in terms of continuous variables θ(x) defined as, where, S eff (θ) corresponds to the effective action along with the contribution of the determinant factor coming from the Haar measure of the unitary group. For the zero 't Hooft coupling case we have the full effective action as: ∞ n=1 a n (β) n dx cos nθ(x) dy cos nθ(y) + 1 2 dx− dy ln sin 2 θ(x) − θ(y) 2 (1.10)

JHEP04(2016)104
The above unitary matrix models can be analyzed using standard techniques in the large N limit. 1 In the large N limit, the dominant contribution to the partition function comes from the extremum condition, δθ(x) = 0 (1.11) which can be written as, where V (θ) is potential obtained by taking the large N limit of S 1 eff [U ] only in the zero 't Hooft coupling case. Introducing the eigenvalue density the above saddle point equation can be rewritten in terms of the density σ(θ) as: with the additional implicit normalization condition dθ σ(θ) = 1. (1.15) Thus, given the effective action (or V (θ)) one has to solve this integral equation to find the saddle point configurations for σ. The solution to this equation depends on the parameters of the theory and the condition that σ(θ) > 0. As we vary the values of the parameters of the theory, the system jumps from one phase to another phase, which demonstrates the phase transition in the system.

Phase diagram from algebraic equation
We follow a method, which is not so well studied, to find different phases of the system and their dependence on the parameters without solving the integral equation (1.12), which in general is difficult to handle. We solve Dyson-Schwinger equation, derived for U(N ) matrix model on S 3 , at large N and obtain the complete phase diagram. This method is advantageous because, in the large N limit Dyson-Schwinger equation becomes algebraic, hence easy to solve. In fact, solution of Dyson-Schwinger equation also gives the eigenvalue density σ(θ) for different phases. One can compute free energy of different phases in this method as well. This technique was applied to compute phase diagram of single or multi plaquette model in [6]. This is an useful technique to solve simple models like σ model, Gross-Neveu model, random matrix model etc. [7]. [7] also used this equation to write down the complete string equation of U(∞) lattice gauge theory.

JHEP04(2016)104
In standard Quantum field theory, the invariance property of the measure can be used to construct functional equations that correlation functions of the theory satisfy. These functional equations are called Dyson-Schwinger equations. Exploiting the fact that the Haar measure [dU ] of a unitary group is invariant under left multiplication and introducing a "Resolvent" R(z) = N −1 T r[(1 − zU ) −1 ], one can write down an algebraic equation for R(z) in large N limit. This is the Dyson-Schwinger equation in unitary matrix model. This equation can be simplified further using the fact that product of single trace operators can be factorized in large N limit i.e. g f = g f , (f, g) being single trace operators).

Phase diagram from Young diagram and free fermions
The integral over unitary matrices in equation (1.5) can as well be done in a different way. Expanding the exponential in terms of character χ of conjugacy classes of a symmetric group and using the orthogonality relation between the characters of U(N ) in different representations, it is possible to write down the partition function in terms of character χ and a sum over representations of unitary group of rank N [8,9]. Since representations of a unitary group can be described by Young diagrams, this is, therefore, equivalent to finding the Young tableaux distribution which dominates the partition function in the large N limit.
It was suggested in [10] that the large N limit of unitary matrix models admits a phase space description. This follows from the fact that eigenvalues of the holonomy matrix behave like position of free fermions whereas, the number of boxes in the Young tableaux are like momentum. This means eigenvalue distribution is like position distribution and arrangements of boxes in the Young diagram is like momentum distribution for N noninteracting fermions. An exact phase space description of different saddle points of a unitary matrix model was constructed in [5]. More precisely, [5] found phase space distribution for three different phases, which appear in a "restricted" class 2 of unitary matrix models.
These three phases are following.
• The first phase corresponds to a trivial representation: no box in the Young diagram. This saddle point correspond to trivial (constant) distribution of eigenvalues.
• The second saddle point is nontrivial. It corresponds to a distribution of boxes in a Young diagram where a finite number of rows are empty. This corresponds to a symmetric distribution of eigenvalues between −π and π, without any gap.
• For the third phase all the N rows of the Young diagram are filled. This phase has one gap in eigenvalue distribution. This is called one-gap phase. 3 Phase space distribution corresponding to these three phases are given in the figure 1.
The shape of the phase space regions is not modified in functional form when one includes perturbative corrections in terms of an effective action involving only the relevant 2 The phase space description holds for a special class of unitary matrix model where the effective action depends on TrU TrU † and its power. 3 This phase is also called one-cut phase, for a reason which will be discussed in the main section.
(b) Representation with finite number of rows empty -symmetric distribution of eigenvalues from −π to π ↔ small black hole.
(c) Representation with all the rows are filled -eigenvalue distribution is symmetric and has a finite gap ↔ big black hole. winding number one modes (powers of TrU TrU † ) or adds a chemical potential. This free fermi picture is very intriguing. One can ask why there should be a free fermion picture for large N saddle points of an interacting theory. We do not think the answer to this question is well understood. On the other hand, it has been discussed in [14] that different phases of dual gravity theory are identified with different phases on the weak coupling side in the context of the AdS/CFT correspondence (see figure 1). However, it is not clear what is the phase on the gravity side corresponding to two gap solution on the weak coupling side. Free fermi phase space distribution could be of great help to understand a dual description of weakly coupled saddle points. However, we must emphasize here that there is no direct way to construct the phases on the string theory (or gravity) side from the phase space distribution, but the nature of transition from one phase to another phase and other characteristics on the weak coupling side may help to understand the properties of phases on the gravity side. Hence, it is interesting and important to understand first if the free fermi description of a large N matrix model is universal, in the sense, if all the large N saddle points have an underlying free fermi picture. [5] did not consider one important phase which appears in the class of matrix model they studied. This phase pops up in the same regime where the third (one-gapped) phase appears. The eigenvalue distribution corresponding to this phase has two gaps. We call this phase as two-gap solution or two-cut solution. In the context of unitary model we are considering, this phase has more free energy than the third one and hence always thermodynamically disfavoured. However, that does not stop this branch of solution to exist. In fact single-plaquette model considered in [6,12] has a stable two-gapped or two-cut phase for a large range of parameters. Large N Chern-Simons theory with a Gross-Witten-Wadia potential also exhibits a stable two-gap solution [13]. Therefore, it is interesting and important to understand the phase space distribution for two-gap branch. To find the distribution in phase space we need to first compute the momentum distribution JHEP04(2016)104 i.e how boxes are distributed in a corresponding Young diagram. The goal of this paper is to find the Young tableaux distribution and phase space distribution corresponding to two-gap solution.
While analyzing different saddle points in terms of dominant unitary representations, we observe that one can classify different solutions depending on number of discontinuities (order N jumps) in h i , where, h i is related to λ i , number of boxes in i th row of a Young diagram by More precisely, we classify the distribution into two different classes, . This implies number of boxes in a Young diagram smoothly decreases from h 1 − N + 1 to h N in the large N limit.

2.
. This corresponds to an O(N ) jump in the number of boxes at i = i . There can be more than one such jumps.
The first class was considered in [5]. This class has two subclasses depending on whether all the rows of a diagram are filled up or a finite number of rows are empty. The first subclass corresponds to one-gap solution whereas, the second subclass maps to no-gap distribution in eigenvalue side. In this paper we show that the second class with one O(N ) jump (or one discontinuity) in h i corresponds to two-gaps in eigenvalue distribution. Having found the Young tableaux distribution we immediately identify that there exists a relation between eigenvalue distribution and Young tableaux distribution which gives rise to phase space distribution for two-gap solution, like other two phases.
Organization of this paper is following. In section 2 we review the basics of the Dyson-Schwinger equation in quantum field theory and derive the same for a class of unitary matrix model. In section 3 we discuss the phase diagram of free theory solving the Dyson-Schwinger equation for a 1 model, explicitly, where we have only one coupling a 1 and the rest a n , for n > 1 are assumed to be zero. We obtain all the three phases, namely no-gap, one-gap and two-gap phases and their free energies. We also discuss how to generalize our procedure to find different phases for a generic model. The phase diagram for weakly coupled theory is discussed in section 4. Section 5 explains momentum distribution or Young tableaux distribution corresponding to different saddle points of the theory. In section 6 we find relation between the Young tableaux and Eigenvalue distributions which enables us to draw the phase space distribution for different saddle points in section 7. In appendix A we derive the Dyson-Schwingerequation for the most generic Unitary matrix model relevant for our case. Then in appendix B we discuss the resolvent and the spectral densities for the most generic cases, while in appendix C we particularly concentrate on a model with only two cycles namely (a 1 − a 2 ) model, and describe its phase space in full detail. Next we move onto the Young tableaux side and discuss the important properties of the resolvent function for the Young tableaux distribution in appendix D, and finally discuss a possible approach to generalize the construction of the resolvent for the (a 1 − a 2 ) case in appendix E.

The Dyson-Schwinger equation
Path integral quantization of any field theory gives us a prescription to directly compute the correlation functions bypassing the Hamiltonian construction, Hilbert space of states and equations of motion. The underlying symmetries of the theory become manifest in path integral formalism. For example, lets consider scalar field theory for simplicity. In the Hamiltonian formalism the Lorentz symmetry was not manifest, as we define the field momentum varying the Lagrangian with respect to time derivative of the scalar field and impose commutation relation between the field momentum and the field. Although the correlation functions are Lorentz invariant but the construction was not manifestly Lorentz invariant.
Let us consider a n point correlation function X = Φ(x 1 ) · · · Φ(x n ) , Since the translation invariance of measure [DΦ] is the defining property of the path integral, we have the following identity This yields the following identity between the correlation functions This is called the Dyson-Schwinger equation associated with a given model in quantum field theory. While the above identity was constructed using the invariance of the measure under translation, there exists more general identities constructed from invariance of the partition function under continuous symmetries. Global symmetries of classical Lagrangian give rise to some currents which are conserved on-shell. This is Noether's theorem. We ask what is the consequence of invariance of functional integral under a symmetry. We shall see that the answer is a quantum generalization of Noether's equation.
Consider a general quantum field theory of a field Φ(x) and its variation where ω's are infinitesimal parameters and G a 's are corresponding generators. Since Φ is integration variable on the right hand side we can first change it to Φ and assuming that the measure is invariant under the above transformation i.e.

JHEP04(2016)104
Here we have used where we have considered ω a 's to be local and j µ a 's are the classical Noether's current. Expanding the exponential up to first order in ω we find, The left hand side up to order ω is given by, Since eq. (2.7) holds for any arbitrary ω(x), we can write This is Dyson-Schwinger equation associated with classical Noether theorem.

Unitary matrix model and Dyson-Schwinger equation
Here we discuss the Dyson-Schwinger equation for the generalized class of Unitary matrix models that we discussed earlier. We start from eq. (1.6). The partition function is invariant under U → U † . Therefore expectation value of any gauge invariant operator which is invariant under conjugation is of interest. Any gauge invariant and conjugation invariant operator can be generated by the functions TrU k and TrU −k . On the other hand, expectation of product of such operators can be factorized in the large N limit. Therefore, it is suffice to consider expectation value of TrU k for any k. Now, TrU k is a real number because, TrU †k , since measure and action are invariant Therefore it is sufficient to consider the expectation values of TrU k for k > 0. It is well known that the Haar measure on the group U(N ) is invariant under action of any element of the group itself. To write down the Dyson-Schwinger equation for this model, let us define the following holomorphic function (2.11)

JHEP04(2016)104
Expanding the right hand side one can write, R(z) = 1 + 1 N z TrU + z 2 TrU 2 + z 3 TrU 3 + · · · (2.12) Since N −1 TrU k lies between −1 and +1 the right hand side of eq. (2.12) is convergent for |z| < 1. Therefore, R(z) is an analytic and holomorphic function in the interior of the unit disk. We use the technique introduced in [6] to derive the Dyson-Schwinger equation in terms or R(z). In large N limit, it turns out that, the Dyson-Schwinger equation is a purely algebraic.
Here we skip the detailed algebra and write down the final version of Dyson-Schwinger equation. The explicit calculation has been provided in appendix A.
The equation is given by, where, andã n = a n /N 2 . (2.15) The contour C being the unit circle. We note that this is an algebraic equation for R(z) in the large N limit. This is, in deed, a very powerful equation. For example, one can study the phase structure of N = 4 SYM theory from the analytic properties of the resolvent. We elaborate this in the next two sections.

Phase diagram at zero coupling
In this section we study the phase structure of the free theory from the analytic properties of the resolvent. Remember that the resolvent R(z) is analytic inside the unit circle in the complex z plane. We first discuss a simple case to understand how a complete phase structure is captured in R(z).
3.1 Zero coupling: "a 1 model" -approximation 1 In the limit λ = 0 the form of the effective action is given by, For further simplification let us assume allã {n 1 ,−n 1 } = 0 for n 1 ≥ 2. We shall see that this simple model captures the important feature of the phase structure at zero coupling.

Different types of solutions
Here we discuss the different type of solutions we obtain for R(z) depending on the value of β 1 which satisfy the analytic properties of R(z). To understand different types of solution we solve eq. (3.3) and find where, From eq. (3.4) it is clear that analytic properties of R(z) inside a unit circle depends on the function F (z). If F (z) has zeros of order one then R(z) will have branch cuts. On the other hand if all the zeros of the function F (z) are of even order then R(z) does not have any branch cut. Therefore we can classify the solutions into two types. No-cut solution and cut solutions.
To understand these two different types of solutions we need to check the zeros of F (z). F (z) is a polynomial of degree four. Hence, F (z) = 0 has four solutions function R(z) can not be analytic inside a unit circle, as it can not be Taylor expanded about the root inside, unless all roots live on the boundary of the unit circle (remember the domain of convergence of R(z) is |z| ≤ 1). Hence, there are three possibilities.
• Type 1 -No-cut solution: Roots are pair wise same, i.e. F (z) has two zeros of order two. In this case R(z) has no cut inside the unit circle. We denote this solution by T1.
• Type 2 -One-cut solution: Two roots are same and other two are different and live on |z| = 1 line. In this case R(z) has a branch cut with end points of the cut lie on the unit circle. We denote this solution by T2.
• Type 3 -Two-cut solution: All four roots are different and lies on the unit circle. We denote this solution by T3.
Note that, although we have chosen the branch cuts lying inside the unit circle but the branch points always lie outside the domain of analyticity of R(z). In fact, one can also choose the branch cuts lying outside the circle. However, according to our choice here, we name the solutions depending on the number of cuts inside the circle.

Type 1: no-cut solution
From the above solutions (3.6) it is clear that when R (0) = β 1 , z 1 becomes equal to z 3 and similarly z 2 is same as z 4 . Therefore, F (z) will have two zeros of order 2, hence R(z) does not have any branch cut inside the unit circle.
Since β 1 = R (0)a 1 , therefore for this branch we have For a 1 = 1 the only solution we have R (0) = β = 0. In this case z 1 and z 3 go to infinity and z 2 and z 4 approach to z = 0 and resolvent is given by, We define spectral density σ(θ) by, The spectral density, which measures the density of eigenvalue distribution of the holonomy matrix, defined in eq. (1.2), is always positive definite. For β 1 = 0, the spectral density is given by, For a 1 = 1 all the roots are finite and non-zero. The resolvent, in this case, is given by,

JHEP04(2016)104
The spectral density for this branch is give by, For a 1 = 1 case we also have, Since, z 1 z 2 = 1 we can write z 1 = r 1 e iθ 1 and z 2 = 1 r 1 e −iθ 1 (similar argument can be given for z 3 and z 4 also). Thus we get, 14) The right hand side has no imaginary part. Hence, there are two possibilities: either θ 1 = π 4 or r 1 = 1. In the first case all four roots lie on negative real axis and This equation has real positive solution only when 0 ≤ β 1 ≤ 1 2 . In the second case, when r 1 = 1, all roots lie on the boundary of the unit circle and This equation is solvable for θ 1 when β 1 > 1 2 . However, for β 1 > 1/2 the spectral density (3.12) becomes negative for some values of −π ≤ θ ≤ π. Hence, β 1 > 1/2 solution is not physical.
Thus, we see that no-cut solution (T1) exists for 0 ≤ β 1 ≤ 1/2 and for this solution the roots are lying on negative real axis in complex z plane. See figure 2.

Type 2: one-cut solution
From the position of zeros of F (z) given by eq. (3.6) we see that T2 occurs when In this case, and z 3 = z 4 = −1. The roots are plotted in figure 3. Since z 1 z 2 = 1, writing z 1 = r 1 e iθ 1 and z 2 = 1 r 1 e −iθ 1 we have Two roots are at z=0, other two are at ∞.
All the roots on negative real axis.
All the roots on the circle. For this branch spectral density becomes negative.  These two roots z 1 and z 2 lie on the boundary of unit circle, hence r 1 = 1 and, Therefore this branch of solution exists for β 1 ≥ 1 2 . The resolvent, for this branch, is given by The condition for one-cut solution can also be written in terms of a 1 from eq. (3.17), Since β 1 ≥ 1 2 this branch exists for a 1 ≥ 1. The spectral density corresponding to this branch is given by, There is a gap in eigenvalue distribution i.e. eigenvalue distribution is zero between |θ| > θ 0 . That is why, T2 phase is also called "gapped" phase.

Type 3: two-cut solution
Two-cut solution exists when all the roots of F (z) = 0 are distinct. This happens when and all four roots lie on the boundary of unit circle (figure 4). The resolvent is given by,

JHEP04(2016)104
where, From this resolvent we find, which implies that two cut solution exists when Combining eq. (3.24) and (3.28) we find that the two-cut solution exists for, or equivalently, The spectral density for this branch is given by, Thus we see that there are two distinct gaps in eigenvalue distribution for two-cut solution.
Here we have obtained only the limiting conditions on the parameters (eq. (3.29) or (3.30)) for two cut solution. However exact dependence can be found from the condition that the number of eigenvalues lying between θ 1 and θ 2 is zero, i.e.
This equation gives a relation between θ 1 and θ 2 and hence, between a 1 and β 1 using eq. (3.26). However, the exact relation is not required for our purpose. Hence we shall not find this relation in this paper.

Free energy for different phases
Free energy, in the large N limit, is given by where, To calculate the free energy we use the following trick. Temperature dependence of partition function comes from the temperature dependence of a 1 , as the holonomy matrix is JHEP04(2016)104 independent of temperature. Differentiating the partition function with respect to temperature we obtain, In the large N limit it becomes, (3.36) For T1, we have either R (0) = 0 or a 1 constant. Hence Free energy is zero. Whereas, in T2, for single-cut solution we have 5 Therefore, we find, (3.39) At β 1 = 1/2 (i.e. R (0) = 1/2) the free energy should vanish hence we get, Substituting the value of C, the free energy for one-cut solution becomes, (3.41) This expression matches with [5].
For two cut solution we have 1 < a 1 < 1 4R (0)(1−R (0)) . Since it is difficult to find the exact relation between R (0) and a 1 , hence it is difficult to calculate free energy for two-cut solution. However, we consider a naive parametric relation between R (0) and a 1 (T ) and find free energy for this parametrization. We take, For x = 0 we get back the no-cut branch corresponding to a 1 = 1 and for x = 1 we get back the one-cut branch. For this parametrization free energy is given by, Since 0 < x < 1, we find free energy of two cut solution is greater than that of one cut solution and hence thermodynamically disfavoured.

Phase diagram for a 1 model
Let us summarize the complete phase diagram as a function of temperature for a 1 -model. a 1 (T ) is a monotonically increasing function of temperature with a 1 (0) = 0. Therefore, the system jumps from one phase to other as we vary temperature. In figure 5 we draw the complete phase diagram for a 1 -model.
• At low temperature, when, a 1 < 1, the minimum action configuration is given by a constant resolvent and hence a constant eigenvalue distribution. The free energy is zero for this configuration (to order N 2 ). This phase is represented by a line (R (0) = 0, a 1 < 1) in the phase diagram.
• As we increase temperature to T = T H , a 1 becomes equal to 1 (a 1 (T H ) = 1). There is a continuous family of minimum action configurations (labeled by a parameter R (0) = β 1 < 1/2) for which the resolvent is linear in z, given by eq. (3.11). Eigenvalue density for this phase (eq. (3.12)) is symmetrically distributed between −π and π. All these configurations also have zero free energy.
• For T > T H i.e. a 1 > 1, there is a new saddle point, with resolvent having a cut inside the unit circle is given by eq. (3.21) and eigenvalue distribution function has a gap. The free energy for this configuration is negative. This branch is given by the upper boundary line of the shaded region in the phase diagram.
• For T > T H there exists another branch whose resolvent (eq. (3.25)) has two cuts inside the unit circle in z plane and eigenvalue distribution has two gaps. But this phase is thermodynamically disfavoured as it has more free energy than one gap phase. This phase lies somewhere in the shaded region in the phase diagram depending on the temperature.
Thus we see that there is a first order phase transition at a 1 = 1, which corresponds to a temperature T = T H .

JHEP04(2016)104
However, the model we have considered so far is the simplest one. But one can generalize the solution discussed in this section for a more generic class of model. In the next subsection we shall consider the next term in the effective action of free theory and discuss the solution of Dyson-Schwinger equation.
3.2 Zero coupling: (a 1 , a 2 ) model -approximation 2 In this section we keep the term proportional to TrU 2 TrU −2 in the effective action and analyze the structure of the resolvent R(z) and the function F (z), from which the phase structure can be obtained. This term is suppressed by an N 2 factor with respect to the first term. We consider this term to understand how to find the solution to Dyson-Schwinger equation for a general class of action. Our discussion in this section is general and valid for a general class of action.
In presence of this term the Dyson-Schwinger equation becomes, This is a quadratic equation of R(z) and it has two solutions. Let us write the solutions in the following form The function F (z) is given by

Different types of solutions -a generic discussion
Although, we refer to the expressions for R(z) and F (z) obtained for the particular model, but our argument holds in general. From invariance of function F (z) under z → 1/z symmetry we immediately realize that for any multiplicity one root λ, there exists another multiplicity one root 1/λ, since z and 1/z satisfy the same equation. This means that roots of multiplicity one always occur JHEP04(2016)104 in distinct pairs. Each distinct pair has product equal to 1. From analyticity condition of R(z) we know that none of the multiplicity one root can lie inside the unit circle, while from the product of pair of roots we know that if one root lies inside, the other one must lie outside the unit circle. This means the multiplicity one roots must lie on the unit circle itself for any F (z) (which is always of even order). Therefore, the generalized form for F (z) can be written as follows.
C being some constant. A pair of multiplicity one root means there exists a a i = a j , for i = j and multiplicity 2 means there exists some a i = a j , for i = j and a i = a k for k = j = i. For one-cut solution to exist (i.e. R(z) has one branch cut inside the unit circle) there must be one a i which occurs only once. From the analyticity condition of R(z), we find a i = −2 cos θ and hence the roots are equal to e iθ , e −iθ . Now, let there exist at least one of these a i =ã 1 (more can occur). Then F (z) can be factorized as The function f (z) will also have the symmetry z → 1/z. A generic form can be written as, However, f (z) may have more pairs of multiplicity 1 factors of form z + a i + 1 z . Now using the condition R(0) = 1 we find a relation forã 1 or in other word cos θ in terms of β i s in the action, (the other coefficients a 0 , .., a p also being given in terms of β i s and R n (0)). For existence of at least one pair of multiplicity 1 root, there must exist a cos θ. If for the choice of parameters we see that there exists no cos θ i.e. its bound is broken, then there can be no pair of roots of multiplicity one. This means the only possible solution F (z) is that all the roots are of multiplicity 2 (or of higher even order). Thus if one cut solution is not possible then only possible solution is a no-cut.
For example, in the previous subsection we have seen that to have one cut solution we must need β > 1/2. Therefore, for β < 1/2 we can have only no-cut solution.
However, one can have no-cut solution even if the condition for one-cut solution is satisfied. But in that case, the roots of multiplicity 2 (minimum) will live on the unit circle and hence R(z) will violate the condition of positivity of spectral density.
We present the discussion on phase diagram for this model in appendix C.

Phase diagram at weak coupling: a phenomenological model
The Dyson-Schwinger equation can be written for the most generic class of weakly coupled theory given by eq. (1.6). However, finding the solution to that equation and hence phase diagram is difficult. One can, instead, construct an analytically solvable phenomenological effective action which captures all the essential features of the original theory [14].

JHEP04(2016)104
In the last section we have seen that in large the N limit the order parameter which distinguishes different phases is R (0) = TrU . Consequently, one can also imagine integrating out all the TrU n (with n = ±1) and obtaining an effective action purely in terms of TrU TrU † . However, this is not easy to carry out explicitly. Therefore we consider a phenomenological models of the form [14] where S(x) is convex and S (x) is concave. This phenomenological model captures all the essential features of thermodynamic properties and phase transition of the unitary matrix model we are considering. The simplest of such model is the so-called (a, b) model [14] in which one keeps only the first two coefficients in S eff (x) given in eq. (4.2).
where a 1 and b 1 are functions of temperature T and λ and b 1 > 0.
In this section, we shall discuss the properties of phase diagram for (a, b) model. However, our discussion can be extended to more general classes (4.2). We present the generic discussion in appendix B.
The Dyson-Schwinger equation (2.13) for this model is given by, This equation is exactly same as zero coupling case (approximation 1) except β 1 replaced by β ab . Therefore, solution of this equation is given by, where, JHEP04(2016)104

Type 1: no-cut solution
The phase structure can be constructed as before. For no-cut solution there are two possibilities: β ab = 0 and R (0) = β ab . The resolvent, for these two cases, are given by R(z) = 1 and R(z) = 1 + β ab z respectively. β ab = 0 phase exists at any temperature, as before. This phase has zero free energy and spectral density is constant for this phase.
The other branch of no-cut solution exists when R (0) = β ab . This implies that, For this branch, two zeros of F (z) are lying at the same point inside the unit circle and on the negative real axis. The other two zeros lie outside the real axis (same as a 1 model). This branch exists when β ab ≤ 1 2 , i.e., The spectral density for T1 is given by,

Free energy
One can calculate free energy of this branch following the prescription given in the previous section.
(4.11) In the large N limit we find, Solving these two equations we find, (4.14)

JHEP04(2016)104
Thus we find, C can be fixed by demanding that at R (0) = 0 i.e. at a 1 = 1, the free energy is zero. Thus we get C = 0. Hence the free energy of Type 1 phase is given by, (4.16) Thus we see free energy is positive for b 1 > 0.

Type 2: one-cut solution
For one cut solution we have the following condition, and the solution exists for β ab ≥ 1 2 . Plugging the value of β ab in terms of a 1 and b 1 we can write, The resolvent is given by, Hence, the spectral density for this branch is given by, where,

Free energy
The free energy of this branch can also be calculated in the similar way as mentioned before. Variation of ln Z is given by, For this phase, . (4.23) Solving this equation we find,

JHEP04(2016)104
Differentiating both sides we find, Substituting da 1 dT in eq. (4.22) we find, where, (4.27) Solving these two equation we find, Thus we find free energy Constant C can be fixed by demanding that the free energy of this branch should match with the free energy of no-cut solution at R (0) = 1/2. This implies, Hence we find free energy of Type 2 branch is given by, (4.30)

Type 3: two-cut solution
Two-cut solution, in the similar way, exists when One can again calculate the free energy of this class for a particular parametrization between R (0) and a 1 , b 1 . It turns out that the free energy of this branch is higher than the free energy of one-cut solution and hence this phase is thermodynamically disfavoured.

Phase diagram and its dual description
Finally we summarize the phase diagram for (a, b) model. However, the basic structure of the phase diagram is same for a generic class of phenomenological model governed by the action (4.2).
• At low temperature the saddle point is characterized by β ab = 0 which has constant resolvent and uniform eigenvalue distribution. This phase corresponds to thermal AdS in the dual gravity setup. This phase has zero free energy.
• Then there is a saddle point with resolvent linear in z when β ab = R (0) i.e.
This is the unstable saddle point corresponding to the small black hole (in the phase where it is to be viewed as an excited string state) and has positive free energy given by eq. (4.16). This saddle point exists in a temperature range T c ≤ T ≤ T H .
• There is the saddle point which obeys This equation has two real solutions for β ab above a temperature T 0 . Let us denote these two values by β S ab and β B ab with β B ab > β S ab . 6 The two values or phases correspond to the Big Black Hole (BBH) and the Small Black Hole (SBH) (in the actual black hole regime) in the dual side. Both these branches have gap in eigenvalue distribution. The BBH (β B ab ) solution exists for all temperatures greater than the minimum T 0 for which this solution exists. While the SBH (β S ab ) solution exists in the interval T o ≤ T ≤ T c . At T c which corresponds to β S ab = 1 2 , this solution goes over into the ungapped solution. This is called Gross-Witten-Wadia transition [15,16].
• There exists a temperature T 1 > T 0 where the stable saddle points β ab = 0 and β B ab exchange dominance. This temperature corresponds to the Hawking-Page temperature where the BBH has a lower free energy than thermal AdS in the semi-classical gravity path integral. The free energy for β B ab phase becomes negative for T > T 1 . Phase diagram in presence of charge or chemical potential has been discussed in [17][18][19].
• There exists another phase in the same regime of parameters where the last phase appears. This phase has eigenvalue distribution with two gaps. But this phase has higher free energy than that of one-gap phase, hence thermodynamically disfavoured.
The gravity dual to this phase is not well understood. 6 At T = T0, β B ab = β S ab . As we increase temperature β B ab starts increasing and β S ab decreases.

JHEP04(2016)104 5 Distribution in momentum space
In this section we shall discuss how different saddle points of weakly coupled gauge theories on S 3 are captured in terms of different representation of U(N ) group in the large N limit. Since different representations correspond to different Young diagrams, therefore, we find how boxes in Young diagrams are distributed corresponding to different saddles of the theory. The advantage of studying the phases in terms of large N Young diagram is that one can directly find out the phase space distribution of different saddle points. As we have already studied, eigenvalues of the holonomy matrix behave like position variables of fermions. At the same time, the representations of U(N ) also have an interpretation in the language of non-interacting fermions with the number of boxes of the Young tableaux being like the momentum [10]. Therefore knowing the momentum distribution is important for drawing the phase space distribution. Our goal, therefore, is to find out the momentum distribution corresponding to no-cut, one-cut and two-cut phases.

Partition function and Young tableaux
We start with the partition function at zero coupling This can be expanded as follows where we use the following notations Here χ R (C( k)) is the character of the conjugacy class C( k) of the permutation group S K and K = j jk j and R denotes the specific representation of U(N ). Finally, using the following orthogonality relation between the character of representations of U(N ) we obtain the following form for the partition function: This is an exact expression for the partition function for any temperature and N of the free gauge theory. However, it should be noted that the answer is completely explicit.

JHEP04(2016)104
A particular representation of an unitary group U(N ) can be labeled by a Young diagram with maximum N number of rows and arbitrary numbers of boxes in each row up to a constraint that number of boxes in a particular row can not be greater than the number of boxes in a row before . Therefore, the sum of representations of U(N ) can be cast as sum of different Young diagrams. If K is the total number of boxes in a particular Young diagram with λ i number of boxes in i th row and i λ i = K, then sum over representations can be decomposed as The partition function can be written as, Note that the total number of boxes is same as the order of the permutation group S K . The characters of the conjugacy class are determined recursively by the Frobenius formula [20][21][22]. Explicit expressions for the most general case are not simple. We shall discuss about that in the next subsection.

Character of permutation group
We discuss in detail the Frobenius character formula [20] for the χ λ (C( k)) of the permutation group. Let C( k) denote the conjugacy class of S K which is determined by a collection of sequence of numbers k = (k 1 , k 2 , · · · ), with, i ik i = K (5.9) C( k) consists of permutations having k 1 number of 1-cycles, k 2 number of 2-cycles and so on. We introduce a set of independent variables, x 1 , · · · , x N , and for a given Young diagram ( λ) we have λ 1 ≥ λ 2 ≥ · · · ≥ λ N ≥ 0. Then we define a power series and the Vandermonde determinant as follows.
For a set of non-negative integers, (n 1 , · · · , n N ), we define Now for our given partition, λ, of a Young diagram we define,

JHEP04(2016)104
Then the character corresponding to a conjugacy class C( k) and a representation characterized by λ is given by This is the most generic formula for character. Using this formula one can in principle write down the partition exactly. However it is still complicated to solve the model and find the saddle points in large N limit in presence of all a i 's. In the next section, we shall extremize the partition function in large N limit to find out different saddle points of the system in presence of one cycles only. This is same as a 1 model discussed in section 3.1. Young tableau distribution for a 1 model has already been discussed in [5]. Authors in [5] found Young tableau or momentum distribution corresponding to no-cut and one-cut solution discussed in the previous section. Here, we extend the discussion and find momentum distribution for two-cut solution.
In appendix E we extend our calculation to find out the saddle points in presence of both 1-cycle and 2-cycle. This is similar to the model we considered in section 3.2.
Before we extremize the partition function to find the saddle point equation let us simplify the formula for character for one cycle. This expression will be useful in the next section.

Character in presence of 1-cycles
We consider a particular conjugacy class C( k) = (k, 0, 0, · · · ), i.e. there is only k 1-cycles. In this case, k 1 = k, and character is given by . (5.15) The Vandermonde determinant can be written as Where σ(i) is an element of the permutation group S K . More over where the sum is performed over a set of non-negative integers, r = (r 1 , .., r N ) such that r 1 + r 2 + · · · + r N = k. Multiplying eq. (5.16) with eq. (5.17) and picking up the coefficient of the

JHEP04(2016)104
Note that i r i = k is automatically satisfied once we consider i λ i = k. Performing the summation r i we can write, Here the sum is only over those σ such that h N −i+1 − σ(i) + 1 ≥ 0. This can be rewritten using the property of the determinant as For completeness here we present a generic formula for the character in presence of k 1 number of 1-cycles, k 2 number of 2-cycles up to k ω number of ω-cycles. The expression for the character is obtained following the same argument given above.

Zero coupling: "a 1 model" -approximation 1
We first look at a model where all the terms are zero apart from the first term a 1 (a i = 0, for i ≥ 2). Which is essentially same as considering k = (k, 0, 0, · · · ). The partition function (5.8) takes the following form Using the expression for character given in eq. (5.20) we can explicitly write the partition function for a model This expression is exact for all N . However, we shall consider large N limit of this model and see how one can capture all the phases obtained from the analysis of eigenvalue density.

JHEP04(2016)104
The large N limit and saddle point equation In the limit N → ∞ we introduce, following [23], a set of continuous variables Since i runs from 1 to N , therefore in large N limit x ∈ [0, 1]. The function λ(x) or equivalently h(x) captures the profile or shape of Young diagram in large N limit. The relation between λ(x) and h(x), in continuum limit is given by The condition h 1 > h 2 > · · · > h N implies a strict monotonicity for h(x), i.e.
h(x) > h(y) for y > x, which implies ∂h(x) ∂x < 0. (5.27) In this continuum limit all the summation over i is replaced by an integration The total number of boxes k in a Young diagram is given by where k is a order 1 number. In the large N limit, as we shall see in the next section, the dominant contribution to the partition function comes from the representation for which the total number of boxes in the corresponding Young diagram is of the order of N 2 . The partition function given in eq. (5.22) can then be written as In the large N limit using Sterling's approximation we find, From the expression of partition function given in eq. (5.31) it is clear that in N → ∞ limit the dominant contribution comes from those configurations for which S eff (h) is JHEP04(2016)104 extremum. Therefore, varying S eff with respect to h(x) we obtain the following saddle point equation.
Following [23] we introduce the density of boxes in the Young diagram defined by The density function u(h) by definition satisfies the normalization The density function has a support between h L = h(1) and h U = h(0). Moreover, from the monotonicity of h(x) (eq. (5.27)) it is easy to show that the density function is positive definite. Also, it obeys the normalization (5.35), hence we find In terms of the Young tableaux density the saddle equation becomes where ξ 2 = a 1 k . Note that the parameter ξ sitting on the right hand side of the saddle equation depends on k where k is given by itself depends of the density function u(h). We have to solve this equation self-consistently. The general interacting unitary matrix model given by eq. (1.6) for weakly coupled gauge theory can also be solved using the above technique. One can in principle write a saddle point equation at large N for the generic action. However, the analysis becomes technically much more involved. But the essential phase structure of the theory is in any case captured by models involving only TrU TrU † as given in the eq. (4.2). In particular, as we have seen, the (a, b) model (eq. (4.3)) does a good job in getting the detailed form of the phase structure. In [5] the authors obtained the saddle equation for Young tableaux distribution function for the class of phenomenological model (4.2). The basic structure of the saddle equation remains same as eq. (5.37), only the definition of the parameter ξ changes. As we mentioned in the introduction, the basic shape of the phase space distribution does not change under addition of perturbative correction, therefore, here we shall not carryout the computation of saddle equation for weak coupling model. We refer [5] for a detailed discussion.

Different classes of solutions
We now discuss the possible classes of solution which satisfy the integral equation (5.37) for Young tableaux density. We classify different saddle points depending on the fact whether the variable h is continuous from h L to h U . Depending on this there are primarily two classes.

For this class the function h(x) is continuous and monotonically decreasing for
N . Therefore, for this branch the difference in number of boxes between two consecutive rows is of the order 1 in the limit N → ∞. This class has two subclasses.

Class 1a (C1a). In this case
A typical Young diagram corresponding to such distribution has been plotted in figure 6.

Class 1b (C1b). For this class distribution is given by
withũ(h) < 1. For this branch u(h) = 1 for 0 < h < p. Therefore a finite number of rows in the particular Young diagram are empty. The distribution has been depicted in figure 7. These two subclasses have been considered in [5]. In this paper we introduce a new class of solution where the support h is discontinuous.

Class 2: h(x) is discontinuous
In this case the function h(x) is discontinuous at some point ζ between 0 and 1 (see figure 8).  The distribution u(h) is therefore given by, A typical representation corresponding to such a Young diagram has been plotted in figure 9.

Finding densities for different classes
To find Young tableaux distribution u(h) we define a resolvent H(h) given by This function has the following properties: for h inside the support.

Class 1a
Using standard techniques [8] we solve integral equations (5.37) with logarithmic kernel. There is a closed expression for resolvent H(h) in terms of a contour integral [3,4]. It follows from the properties of H(h) listed above. where contour C is an anticlockwise contour around the branch cut [p, q]. By deforming the contour one can take the contour to be a clockwise loop about h and about the logarithmic branch cut from −∞ to 0 and back (see figure 10). Carrying out the contour integration we obtain, From the asymptotic expansion of H(h) one finds, and k = √ qp + 1 4 which implies, Since p, q are real and positive therefore, this branch exists for ξ ≥ 1 2 (from eq. (5.48)). From eq. (5.49), we therefore, conclude that this class of solution exists for a 1 ≥ 1.
From the discontinuity of resolvent we find that the Young tableaux density u(h), for this branch, is given by

JHEP04(2016)104
If h 1 and h 2 are two roots of this solution then we find For this branch, we also note that p + q = 1 + 2ξ, and pq = ξ − 1 2 An alternate way to compute the resolvent This is an important section. We understand that finding an exact expression for resolvent is important to get distributions of boxes in Young diagram. In the last section we performed the contour integration explicitly and found H(h). However, the expression for resolvent becomes complicated, as we shall see later, when variable h becomes discontinuous between h L and h U and hence solving the contour integration turns out to be difficult. In this section we explain how one can obtain the above result for resolvent from its asymptotic properties without doing any contour integration. We take the following ansatz for resolvent H(h) following [24] H(h) = 2 ln where g(h), f (h) and v(h) are polynomials of h and κ is constant.
If H(h) has single branch cut between p and q then, g 2 − f 2 can at most be polynomial of degree 2, hence g and f are polynomial of degree 1 or less. 8 From the property (d) of H(h) given above one can show that, Therefore we find, From the asymptotic expansion of H(h), given in property (c), we find  Figure 11. Contour for Solution Class 1b.
Substituting the values of these functions in eq. (5.54) we obtain the expression for H(h) which matches with eq. (5.47). Since the resolvent has a branch cut between p and q, therefore Solving this equation one finds eq. (5.48). Thus, we recover all the essential relations including the expression for the resolvent without solving the contour integration.

Class 1b
Using the ansatz (5.39) for u(h), for this class the saddle point equation becomes, The resolvent as defined in eq. (5.44) becomes (with h L = 0 and h U = q), Again, one can write an expression for H(h) depending on its analytic properties. It is given by, The contour is shown in fig 11. However, we shall not perform this contour integration to find H(h), rather we shall follow the technique defined in the last section.

Let us write
We take the ansatz for H 1 (h) as before, It follows from eq. (5.61) that H 1 (h) satisfies, Hence, from eq. (5.66) and asymptotic properties of H(h) we find that, Plugging the values of these functions in to the ansatz we find, Since, H(h) has a square root branch cut between p, q, we find that the term inside the square root can be written as (h − p)(h − q), with p = 1 − 2ξ, q = 1 + 2ξ. (5.70) Hence, further simplifying H(h) we finally get, which matches with the expression given in [5]. Since p = 1 − 2ξ and p > 0, therefore this branch exists for ξ ≤ 1 2 . From eq. (5.68) we find, k = ξ 2 .

JHEP04(2016)104
The former implies the uniform distribution This is therefore a saddle point for any value of a 1 . This is in fact the density corresponding to the trivial representation n i = 0. The latter corresponds to a family of saddle points labeled by 0 < ξ < 1/2 which exists only at a 1 = 1. Therefore, we see that Class 1b corresponds to no-cut branch Type 1 discussed in section (3.1.1) From the discontinuity of the resolvent we find the Young tableaux distribution for this branch is given by,ũ We define the resolvent for this case The resolvent has the following properties.
(a) It is a analytic function of h in the complex plane with two branch cuts from p to q and r to s on positive real axis.
(b) It is real for other values of the positive real h.
, which follows from the definition of a complex function with a branch cut on the real line.

Derivation of properties (d) and (e) has been discussed in appendix D.
The resolvent H(h) for this class can be constructed analogously where the contour C is defined anticlockwise around each of the branch cuts [p, q] and [r, s]. See figure 12.
The contour can be deformed and can be taken around the point h, infinity and along the branch cut of the logarithm, as before. But the integral obtained in this manner is extremely tedious. Instead we use an alternate method. Alternate method to find resolvent in Class 2 We take the following ansatz for H(h) Since H(h) has two branch cuts between [p, q], [r, s] therefore, is a quadratic polynomial and f (h) is quadratic at most.
From property (d) of H(h) we find, From the asymptotic expansion we find that, p + q + r + s = 2ξ + 1, p(qr + qs + rs) + qrs = −(2ξ − 1)b, p(q + r + s) + q(r + s) + rs = 2b + ξ − Since p, q, r, s are greater than zero, this implies this branch of solution exists for ξ ≥ 1 2 . Now from the definition of u(h) we have, This expression can be written in terms of inverse trigonometric function as Plugging the values of g(h) and f (h) we obtain which gives a quartic equation for h, A typical distribution u(h) has been plotted in figure 13. We note that the distribution is zero between q and r. Let us denote the solution of this equations by h 1 , h 2 , h 3 and h 4 . It follows from the properties of quartic equations that, Comparing these relations with eq. (5.52) and (5.53) we see that in this case also sum over roots are same as 1 + 2ξ cos(πu(h)). Also it is interesting to note that, if we set r = s = 0 then we get back Class 1a. For r = s = 0, b = 0 and we find p + q = 1 + 2ξ and pq = ξ − 1 2 2 . All other equations in (5.84) are trivially satisfied.
Finally, comparing the coefficient of 1 h 2 term in asymptotic expansion of H(h) we find that, Substituting the value of b in the last relation or eq. (5.84) we find, Since p, q, r, s > 0, the above relation implies that this branch is possible for a 1 > 1. Also, since b = − √ pqrs < 0, from eq. (5.93), we see that a 1 has an upper bound as well Thus, we identify this branch to Type 3 branch in the eigenvalue side.
A dictionary. Comparing the equations in section 3.1 and this section we find that different parameters on both sides are related in the following way.
It is easy to check that, using the above dictionary one can obtain the equations between the parameters (saddle point conditions) in one side from the equations on other side for a 1 model.

JHEP04(2016)104 6 Relation between the Young tableaux and eigenvalue distributions
The analysis of the different phases of gauge theories in terms of dominant representations (or Young diagrams) of U(N ) groups is very different from the usual eigenvalue analysis reviewed discussed in section 3.1. However, it was observed in [5] that there is nevertheless, a simple relationship between the saddle point configurations u(h), in both the high and low temperature phases, with the corresponding saddle point eigenvalue densities. We shall first review the relationship between Type 1 and Class 1b and Type 2 and Class 1a as described in [5].

Identification between Type 1 and Class 1b
Consider first the low temperature saddle point u(h) = 1, ξ = 0. The corresponding saddle point for the eigenvalue density is σ(θ) = 1 2π , β 1 = 0. From the form of these two distributions, we notice that they are functional inverses of each other. Therefore, we can make the identification = σ(θ) .

Identification between Type 2 and Class 1a
In this case the identification is more subtle. In this case we modify the identification (6.1) to, u = θ π and where, h 1 and h 2 are two roots of eq. (5.51) with h 1 > h 2 . Using the identification between u and θ one can write from eq. (5.51)

Identification between Type 3 and Class 2
Following the identification between Type 2 and Class 1a we now generalize the identification between Type 3 and Class 2. Note that for Class 2, h satisfies a quartic equation (5.91) with four roots denoted by h 1 , h 2 , h 3 and h 4 . Since, in the limit b = 0 Class 2 boils down to Class 1a, we expect our identification in this section will also reduce to eq. (6.5).
We consider the following identification, 2π . (6.9) Using the above identification between u and θ we find from eq. (5.91) Finally, from the properties of roots of the above equation and eq. (5.94), where, Therefore, this expression exactly matches with eq. (3.31) with β 1 = ξ.

Fermionic phase space
The identification between the set (h, u(h)) and (θ, σ(θ)) at the saddle points has a very natural interpretation, as pointed out in [5], in terms of a free fermionic picture. The fermionic picture follows from the following two facts.
• The eigenvalues θ i s of holonomy matrix U(N ) behave like coordinates of free fermions.
• The number of boxes n i s (or h i s) in a Young diagram being like momentum of noninteracting fermions [10].

JHEP04(2016)104
The eigenvalue distribution σ(θ) tells how the position coordinates of N fermions are distributed whereas, u(h) describes the distribution of momenta of fermions.
[5] considered a phase space distribution which gives rise to these individual distributions. In the classical (i.e. large N ) limit, one can describe this system of N fermions in terms of an incompressible fluid occupying a region of the two dimensional phase space.
Therefore, different saddle points are described by some configuration or region Γ in phase space with phase space density function ω(h, θ) defined in a two dimensional plane spanned by h and θ. h is a radial direction and θ is angular. The phase space density has the property Momentum and position distributions are defined in terms of phase space density as follows, where the first integral is at constant h and the second at constant θ. Also note that It follows from the normalization of two partial densities. 9 Our next goal is to find out the boundary of the region Γ for different saddle points. We assume that the boundary is specified by two curves h + (θ) and h − (θ) where, where h i 's are roots of the equations satisfied by h for different classes. For Class 1b the governing equation is given by eq. (6.3). This is a linear equation and has one solution. Therefore, h + (θ) = h 1 and h − (θ) = 0. This distribution in phase space is therefore drawn in figure 14.
For Class 1b the governing equation is quadratic (6.6). It has two solutions h 1 and h 2 .
The distribution is depicted in figure 15. For Class 2 the distribution is nontrivial. Here, we take,  The equation which governs the shape of the fluid droplets is given by, For this case also eigenvalue distribution is given by eq. (7.5). The phase space distribution is given by figure 16.
The topology of the phase space distribution for two-gap solution is different than the other two branches. In fact, this topology is robust and does not change under inclusion of weak coupling terms of the form given in eq. (4.2). 10 This is because the equation satisfied by h remains quartic even after addition of these terms. Only the coefficients of different powers of h change. In fact, the generic (a, b) type model given by eq. (4.2) can have at most two-gap or two-cut solution and hence the topology of phase space distribution is also fixed. However, finding distribution for generic weak coupling action is difficult [25].

Acknowledgments
We wold like to thank Rajesh Gopakumar for useful discussion and his valuable comments on the manuscript. We express our sincere gratitude to N. Banerjee, R. Banerjee, A. Bhand, A. Mangalsuli, S. Sarkar for useful discussion. We appreciate the effort of S. K. Pathak for arranging important articles required for this work. PD would like to acknowledge the hospitality of IISER Pune where a part of this work was done. Finally, we are grateful to people of India for their unconditional support towards researches in basic sciences.

A Derivation of Dyson-Schwinger equation for unitary matrix model
We use the technique introduced in [6] to derive the Dyson-Schwinger equation. Firstly the generic partition function that we start with is as follows: with i n i = 0 and the constants a n , are such that the action is invariant under U → U −1 . Next let us look at the following expectation value where X is a skew adjoint N × N matrix. We change the integration variable U → e tX U . Therefore we get, The

JHEP04(2016)104
That is because under this transformation we have the action on the Haar measure: Under redefinition of each of the integration variables θ i → θ i − θ i , the measure remains invariant i.e. dθ i remains invariant and thus the Haar measure remains invariant. Now the right hand side is independent of t, therefore we can write, Tr[e n i tX U n i ] 0 = 0 We use the following result to simplify the above expression.
Therefore from eq. (A) we get, Eq. (A.8) holds for any X. There are N 2 independent N × N anti-hermitian matrices. We write eq. (A.8) for each such independent X and add up to get

JHEP04(2016)104
The first term can be written as, Using the following identity for appropriately normalized basis {X a } of the skew adjoint matrices Similarly the second term can be written as, Thus we get, In the large N limit one can use the factorization i.e. f (U ) g(U ) = f (U ) g(U ) , where f, g being single trace operators. Using this factorization eq. (A.13) becomes,

JHEP04(2016)104
This is the Dyson-Schwinger equation for this matrix model. One can again simplify the equation to write the moments of single trace operators completely in terms of resolvent. From eq. (2.12) we can calculate the moments of the single trace of matrix U from the resolvent where the contour has been chosen around the point w = 0. For n j > 0, one can write, Since the point z is inside the unit circle, as evident from the taylor expansion of the resolvent, the contour C now encompasses the entire unit circle. For n j negative we have, Taking trace and expectation value on both sides we can write, Therefore, the Dyson-Schwinger equation becomes, where,

JHEP04(2016)104
We see that this is an algebraic equation for R(z) in the large N limit. This is a very powerful equation. From the analytic properties of the resolvent one can study the phase structure of N = 4 SYM theory. Redefining a n /N 2 =ã n for simplicity and taking the large N limit eq. (A.19) becomes, The partition function for the generic zero coupling case is: For this generic zero-coupling model (restricting to order l) the resolvent takes the following form: Where F (z) can be rewritten after some algebra: Where the dots mean similar terms with β 3 and so on. It is evidently in z → 1 z form. The no-cut solution, following arguments from earlier is: One can easily write down the generic form of the one-cut solution for zero-coupling model including up to order l. The Resolvent in that case assuming a one-cut solution will be:

JHEP04(2016)104
Where β l =ã +l,−l R l (0)/(l − 1)! and the coefficients ω j are obtained by setting R(0) = 1 and cos(θ 0 ) obtained analogously as earlier by equating the two equations one obtains for a single ω j , (as one had done for l = 2 case). One can also note that the following equations can be obtained from the relation R(0) = 1, and so the coefficient of z −i for i > 0 is zero.

JHEP04(2016)104
The contour being taken around origin. Putting the form of R(z) in this we obtain: Using the formula for generating function of Legendre polynomials again, and the residue theorem we obtain:

B.2 generic non zero coupling
Moving onto addition of any generic term for the nonzero coupling case, which can be dealt in a similar fashion. Let us consider a term, characterized by an integer l in the following way: we have the following as the redefinition of β i taking into account sum over i, 2 ≤ i ≤ l (i.e. also including generic zero coupling case discussed earlier).
where as usual n i ic i = j jd j = l. Under this redefinition, the non-zero coupling equation for the Resolvent becomes exactly same as zero coupling and so the spectral density in terms of β i will be same. But now the relation between theã parameters is extremely complicated. So to eliminate the variables β i in favour of the coupling constants is a non trivial task as there no longer remains a linear set of equations to solve. Nevertheless the spectral density takes exactly the same functional form in these variables, which is indeed non trivial. In fact any model which has an action, polynomial in traces of n th order, and invariant under the operation U → U −1 will have the identical functional form of the eigenvalue distribution and the Resolvent to the placquette model of n th order, with the parameters β i mapped to the coupling constants of the placquette model.
Notably when we have the specific case where: which is the generalization of the phenomenological a, b model. This will be equivalent to the simple one placquette model under the redefinition:

C.1 Solution class 1: no-cut solution
For this class, the resolvent is analytic inside the unit circle and also R(0) = 1, which means the negative powers of z must cancel out. We note this is the case if F (z) is a perfect square, and from above we have From this we find the resolvent is given by (we choose + sign in eq. (3.45) otherwise R(z) would have isolated singular points), We shall discuss the conditions on the parameters β 1 and β 2 for which this class of solution is valid in the next subsection. The spectral density for this branch, is given by 2πσ(θ) = 1 + 2β 1 cos θ + 2β 2 cos 2θ.
Defining cos θ := x, this equation can be written as, and σ(x) ≥ 0 in that range of x. σ(x) defines a parabola in σ(x) − x plane it can have only one extremum. Since σ(x) > 0 there are three possibilities.
(a) σ(x) has a minimum between −1 and +1 and the value of the function is greater than or equal to zero at minimum.
(c) σ(x) has a maximum −1 and +1 and the value of the function is greater than or equal to zero at maximum.
For the first case β 2 > 0 (since 2πσ (x) = 8β 2 ) and the minimum occurs at x = − β 1 4β 2 . If the minimum lies between −1 and +1 we get, Also demanding that the value of the function at this point is greater than or equal to zero, we get, The blue area in figure 17a denotes the corresponding region in (β 1 , β 2 ) plane. The end points are given by (± 2 3 , 1 6 ).

JHEP04(2016)104
Therefore we need to find a domain in (β 1 , β 2 ) plane for which there exists a real solution for x between 0 and 1 in eq. (C.15) and σ(θ) ≥ 0. The phase space has been throughly discussed in [12]. Interested reader are referred to that paper.
The phase diagram can also be discussed in the similar way as in [12].
Thus we see that the complete phase diagram of N = 4 super Yang-Mills theory can be captured studying the analytic properties of the resolvent in a complex z plane.

D Derivation of properties of resolvent H(h)
Here we discuss the properties of the resolvent H(h). Firstly lets take the case for 1−cut solution for simplicity. The asymptotic expansion of H(h), for h → ∞, can simply be obtained as: To show this we take the definition of H(h) with discontinuous support: Now if h fall on the first branch q > h > p, then: The second term is automatically equal to: The character in this case is explicitly written as: where l i = n i + N − 1, n i being the number of boxes in the i th row of the Young diagram.
Using the fact that the sum over representations can be understood as sum over Young diagrams, we have the partition function equal to: Rewriting the sum, we have: k 1 ,k 2 ,l i ,q i ,s i a k 1 1 a k 2 2 k 1 !k 2 ! N j (l j − 2q j )!q j !(l j − 2s j )!s j ! r<p (l r − l p − 2(q r − q p ))(l r − l p − 2(s r − s p )) In the large N limit then, the partition function takes the following form: Z = dk 2 dl(x) dq(x) ds(x) dλ exp[−N 2 S eff (l(x), q(x), s(x), k 2 , λ)] (E.5) Where λ is a Lagrange multiplier and S eff is given by:

JHEP04(2016)104
Where the variables have been defined in the large N limit analogously to the case for one cycle. Also K = N 2 K etc.. Note that we have performed the integration over k 1 in the path integral and thus replaced k 1 by K − 2k 2 . More over one can use the following relation dx l(x) = K + 1 2 and 1 0 dx q(x) = k 2 =

E.1 equations of motion
Next we calculate the equations of motion for this action by taking variation of the above with respect to the variables l(x), s(x), q(x) and k 2 . These are respectively as follows: Using equation of motion for l(x) we get: ln a 2 k 2 a 2 1 (K − 2k 2 ) 2 + 2[lna 1 + ln(K − 2k 2 )] − ln(q(x)s(x)) = 0 ⇒ q(x)s(x) = a 2 k 2 (E.13)

JHEP04(2016)104
Again using equation of motion for k 2 we can calculate the value of the Lagrange multiplier: ln a 2 k 2 a 2 1 (K − 2k 2 ) 2 = 2 − λ (E.14) The above set of equations are extremely difficult to solve, primarily because there is no apriori information about the behaviour of q(x) or s(x). If we assume that either q(x) or s(x) is monotonically increasing then the other is of course monotonically decreasing. whereũ is the Young tableaux density for s. It will be interesting to solve these set of equations, and obtain u(l). This will allow one to check the phase space relationship more precisely and possibly generalize this to the whole class of models discussed in the eigenvalue sector. We leave this problem for later purpose.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.