Emergent Phase Space Description of Unitary Matrix Model

We show that large $N$ phases of a $0$ dimensional generic unitary matrix model (UMM) can be described in terms of topologies of two dimensional droplets on a plane spanned by eigenvalue and number of boxes in Young diagram. Information about different phases of UMM is encoded in the geometry of droplets. These droplets are similar to phase space distributions of a unitary matrix quantum mechanics (UMQM) ($(0 + 1)$ dimensional) on constant time slices. We find that for a given UMM, it is possible to construct an effective UMQM such that its phase space distributions match with droplets of UMM on different time slices at large $N$. Therefore, large $N$ phase transitions in UMM can be understood in terms of dynamics of an effective UMQM. From the geometry of droplets it is also possible to construct Young diagrams corresponding to $U(N)$ representations and hence different large $N$ states of the theory in momentum space. We explicitly consider two examples : single plaquette model with $\text{Tr} U^2$ terms and Chern-Simons theory on $S^3$. We describe phases of CS theory in terms of eigenvalue distributions of unitary matrices and find dominant Young distributions for them.

Abstract: We show that large N phases of a 0 dimensional generic unitary matrix model (UMM) can be described in terms of topologies of two dimensional droplets on a plane spanned by eigenvalue and number of boxes in Young diagram. Information about different phases of UMM is encoded in the geometry of droplets. These droplets are similar to phase space distributions of a unitary matrix quantum mechanics (UMQM) ((0+1) dimensional) on constant time slices. We find that for a given UMM, it is possible to construct an effective UMQM such that its phase space distributions match with droplets of UMM on different time slices at large N . Therefore, large N phase transitions in UMM can be understood in terms of dynamics of an effective UMQM. From the geometry of droplets it is also possible to construct Young diagrams corresponding to U (N ) representations and hence different large N states of the theory in momentum space. We explicitly consider two examples : single plaquette model with TrU 2 terms and Chern-Simons theory on S 3 . We describe phases of CS theory in terms of eigenvalue distributions of unitary matrices and find dominant Young distributions for them.

Introduction and summary
Unitary matrix model (UMM) has a wide range of applicability both in mathematics and physics. Partition functions of different super-symmetric gauge theories, in particular Chern-Simons theory on certain manifolds boil down to UMMs. These models also have applications in a broad class of condensed matter systems.
A UMM is a statistical ensemble of unitary matrices (denoted by U ), defined by the partition function On the other hand, it is well know that eigenvalues of hermitian (or unitary) matrices behave like positions of free fermions [1,2]. In hermitian matrix quantum mechanics, dynamics of eigenvalues can be thought of N non-interacting fermions moving in a potential V (θ). At large N all the eigenvalues (fermions) occupy the lowest possible energy states (ground state configuration) of the potential upto a fermi level. Different large N phases of the system, therefore, can be explained in terms of shapes of fermi surfaces in phase space. All the states below that surface are occupied, like Thomas-Fermi distribution at zero temperature [3].
Phase space description of matrix quantum mechanics seems to be quite natural because explicit time dependence of eigenvalues allows one to define conjugate momenta associated with eigenvalues. Interacting gauge theories (supersymmetric as well as non-supersymmetric) on compact manifolds can be written as a 0-dimensional UMM [4][5][6][7][8][9]. Eigenvalues of these unitary matrices do not dependent on time (or any other parameter). Therefore, it would be interesting to ask if one can provide a similar phase space description for different large N phases of these theories in terms of droplets of constant densities in two dimensions. In [10], with a bit of surprise, it was observed that for a very restricted class of matrix models, such description is indeed possible. Writing the partition function as a sum over all possible representations of unitary group one can study large N behaviour of the modes in terms of most dominant representation which is characterised by a distribution function u(h). u(h) measures how boxes h are distributed in a Young diagram. Solving the model at large N , [10] found that information of different phases of the theory can be captured in terms of topologies (shapes) of different droplets in two dimensions. This two dimensional plane is spanned by θ and h. Eigenvalue and Young tableaux distributions for different phases of the theory can be obtained from a distributions function ω(h, θ) in (h, θ) plane with property ω(h, θ) = 1 for (h, θ) ∈ R = 0 otherwise. (1.2) Eigenvalue distribution ρ(θ) for a given θ can be obtained from the distribution by integrating over h Similarly, for a given h, integrating over θ we obtain Young tableaux distribution, The droplet picture of large N matrix model is very interesting, as these droplets are similar to those of Thomas-Fermi model where Wigner distribution is assumed to take constant value inside some region in phase space and zero elsewhere. If eigenvalues of unitary matrices are like positions of fermions with distribution function ρ(θ), then from droplet picture it seems that h plays the role of momentum with distribution u(h). Therefore, we see that droplet picture raises a very intriguing question if there is any underlying quantum mechanics (quantum mechanics of N free fermions) associated with different large N phases of UMMs (or interacting gauge theories).
Before answering this question, we first need to understand if phase space description is a universal characteristic of UMM. In other words, if all large N phases of any UMM can be described in terms of topologies of droplets. An important relation between character of conjugacy class and number of boxes in a Young diagram is required to understand this question. For a class of model considered in [10,11], the relation is known and simple (character is equal to dimension of representation). However, in general this relation is nontrivial and there exists no exact formula for character 1 . Therefore, it is difficult to construct u(h) directly from partition function. However, in [12], a more general class of matrix model was considered, namely plaquette model and shown that even if it is difficult to find u(h) directly from partition function, but a distribution function ω(h, θ) in (h, θ) plane can be obtained in large N limit. Our first goal, in this paper, is to establish this result for a generic UMM and show that different large N phases of a generic unitary matrix model also have an emergent "phase space" description. Our generic model includes multi-trace terms in action. Therefore, this generic matrix model can describe weakly coupled gauge theories on compact manifold.
Before delving deep, it is worth taking a digression and mention about some important works and literature in mathematics and their relevance with our work. As we have already mentioned, different large N phases of gauge theory correspond to different asymptotic behaviour of dominant Young diagrams. Asymptotic growth of Young diagrams is a well known process in mathematics. It originated from the pioneering works by A. Vershik and S. Kerov. Details of some of their related works can be found in [13]. Kerov has studied a Markovian growth process of Young diagrams also known as Plancherel growth process (due to involvement of Plancherel measure), and shown that almost all of Young diagrams become uniformly close to a common universal curve. In the current context the work of [14] is also worth mentioning, who gave the asymptotic behaviour of the characters when the corresponding Young diagram converge to some prescribed shape. The key difference between the study of Kerov and others and our work is that we are interested in large N shapes of the Young diagram in presence of a potential. As we change the parameters of the theory the potential changes and hence the asymptotic behaviour of the corresponding Young diagram. This is the reason why we are not getting any universal behaviour of the Young diagrams involved in our case and precisely this non-universal behaviour is the reason we have different shapes of the droplets which in turn helps one to differentiate between large N phases of the matrix model involved.
Coming back to the main context, in order to understand the underlying quantum mechanics associated with large N phases of UMMs, in this paper we observe that droplet picture in (h, θ) plane is similar to phase space distribution of a unitary matrix quantum mechanics (UMQM) on a constant time slice. Partition function of UMQM is equivalent to a theory of free fermions [1,2]. Time evolution of this system can be written as time evolution of two different branches of fermi surface [3]. These different branches of fermi surface define a closed region in phase space. At large N all the states inside that region are filled, otherwise empty. There is a relation between eigenvalue (position of fermions) and corresponding momentum which defines the boundary of the region [3]. We see that droplets of UMM have similar feature as phase space distribution of an effective UMQM on a constant time slice. Hamiltonians of both the systems match on that constant time slice. Therefore, given a UMM one can find an underlying matrix quantum mechanics such that its evolution with time from t = t i to t = t f captures two different phases of UMM on initial and final time slices. This allows us to understand large N phase transition in UMM in terms of dynamics of an effective UMQM. Considering parameters of UMQM to be time dependent one can find a time evolution in effective UMQM such that distributions on initial and final time slices correspond to two different phases of UMM [15]. The idea here is similar to the work of Matytsin [16], who showed that leading term in a non-trivial large-N limit of the Itzykson-Zuber integral can be written in terms of an action functional corresponding to the complex inviscid Burgers (Hopf) equation with boundary conditions specifying the eigenvalue distribution on two different time slices 2 . In fact, the time evolution of the action functional is given by a Hamilton-Jacobi form with Hamiltonian similar to collective field theory Hamiltonian of [18,19]. As we observe in this paper that for a one dimensional unitary matrix model or equivalently an unitary matrix quantum mechanics one ends up with 4.19, which can be re-written as the Hopf equation. Solution of this Hopf equation determines the evolution of asymptotic shape of Young diagram as the system evolves from one phase to other phase in the large N limit. Thus evolution of shape of Young diagrams could be related to time evolution of underlying matrix quantum mechanics. Organisation of our paper is following.
• We start with a very brief introduction of UMM in section 2.
• In section 3 we discuss emergent phase space description of a generic UMM.
• Relation between UMM and UMQM has been studied in section 4.
• We provide an explicit example of reconstruction of Young diagrams from phase space for β 1 -β 2 model in section 5.
• In section 6 we discuss phase space distribution for Chern-Simons theory on S 3 . Here we present a unitary matrix model analysis of large N phase of CS theory. Then we obtain Young distribution for different values of 't Hooft coupling.
• We end our paper with some discussion and out look in section 7.
2 [17] In this basis Haar measure takes the following form and action S(U ) becomes function of N variables θ 1 , · · · , θ N : V ({θ i }). The partition function in this basis is given by where, The partition function can be thought of as describing a system of N particles interacting by a repulsive potential (Coulomb repulsion) − ln | sin(θ i − θ j )/2| and moving in a common potential V ({θ i }). Because of the strong repulsion, two particles (eigenvalues) do not come close to each other. If we neglected the Coulomb repulsion, all eigenvalues would have sat at the minima of potential V (θ). But, due to the Coulomb repulsion they spread around these minima and fill some finite intervals.
In large N limit one defines a set of continuous variables and in terms of these continuous variables the effective action becomes, At large N , dominant contribution comes from the saddle points of S eff ([θ(x)]). We define eigenvalue distribution function which captures information about distributions of eigenvalues between −π and π in large N limit. Different saddle points correspond to different eigenvalue distributions and hence different phases of the theory.

Gauge theory at zero coupling
Thermal partition function of gauge theory on compact manifold (S 3 ) can be written as a UMM. In zero 't Hooft coupling limit (λ = g 2 Y M N → 0), all the modes of the theory are massive with mass terms depending on the size of S 3 . Only the zero mode of temporal component of gauge field is massless and strongly interacting. Therefore, in zero coupling limit one can integrate out all the massive modes and obtains an exact partition function in terms of the zero mode. For example, thermal partition function of a SU (N ) gauge theory on S 3 × S 1 at zero 't Hooft coupling can be written as, where U = e iβα , α being zero mode of temporal component of the gauge field and β (radius of S 1 ) is inverse temperature. a n (β)'s are related to single particle partition function of bosonic and fermionic modes. See [5][6][7]9] for details.

Gauge theory at weak coupling
The above equation is derived at zero 't Hooft coupling. One considers only one loop contribution while integrating out massive modes. In weak 't Hooft coupling limit one can integrate out massive modes order by order in loop expansion and obtains the following partition function in terms of the holonomy U , (2.10) Terms, in S eff (U ), with k traces appears in (k − 1) loops in perturbation theory and have a planner contribution starting with λ k−2 . For example, two trace terms (2.9) comes in one loop order in perturbation theory. In this paper, we show that different large N phases of this model has an emergent phase space picture.

Single plaquette model
A single plaquette model is given by where β n (β)'s are parameters. This partition function is equivalent to 2.9 with β n (β) identified as a n (β) TrU n in the large N limit. They have similar phase structure.
The single plaquette model has lot of applications in lattice gauge theory. Also, as we shall see, Cherns-Simons theory on S 3 can also be written as single plaquette model with β n 's, determined in terms of 't Hooft coupling. Single plaquette model has a rich phase structure. We refer to [20,21] for a detailed study.

Emergent phase space distribution
It was first observed in [10] that there exists a relation between eigenvalue density and Young tableaux distribution for a restricted class of matrix models. It was also shown that using these relations, different large N phases of the theory can be described in terms of shapes of droplets in two dimensional plane spanned by eigenvalue and number of boxes in Young diagram. However, observation in that paper was accidental and valid for a very limited class of UMM. Also, it was not clear if such identification follows from any extremization condition at large N . In [12], work of [10] was generalized to single plaquette model. Most interestingly, [12] showed that the relation between Young distribution and eigenvalue distribution follows from an extremization condition. In this paper we further extend this work to a most generic matrix model of the form (2.10). Although, we work here only upto three trace terms, our construction can be generalised to any arbitrary number of traces.

Two dimensional droplets for plaquette model
We start our discussion with a brief review of [12]. This discussion will help us to follow the derivation for the most generic case.
Partition function of single plaquette model (2.11) is given by, β n n (TrU n + TrU †n ) , β n 's are some arbitrary coefficients. (3.1) To obtain phase space distribution, we expand the exponential and write partition function as a sum over representations R of unitary group U (N ), Here χ R (C( k)) is the character of conjugacy class C( k) of permutation group S K , K = n nk n and Sum of representation of U (N ) can be written as a sum over different Young diagrams. If K is the total number of boxes in a Young diagram with λ i being the number of boxes in i-th row then N i=1 λ i = K. Hence, sum over representations can be decomposed as The partition function, therefore, can be written as, We introduce N variables h 1 , · · · , h N , related to number of boxes λ i 's as h i 's are shifted number of boxes 4 . From monotonicity of λ i s it follows that h i s satisfy the following constraint From now on we shall use variables h i s in stead of λ i s.
In large N limit we define continuous variables (3.10) In this limit summation over i is replaced by an integral over x, and sum over representations ( λ) and sum over cycles (k n ) are given by path integral over h(x) and integral over k n . Writing characters χ h in terms of h(x) and k n , partition function (3.7), in large N limit, can be written as where S eff is the effective action. Dominant contribution to partition function comes from those representations which maximise effective action S eff . To find the most dominant representations at large N , we introduce a function called Young tableaux density defined as, For the simplest case, β 1 = 0 and other β n = 0 (for n ≥ 2) one can find u(h) for different phases of the model [10,11]. In this simple case we need to calculate character of permutation group in presence of one cycles only (k n = 0 for n ≥ 2) which is given by dimension of the corresponding representation. As a result is was possible for [10,11] to extremize the effective action to obtain u(h) for different phases 5 . [10,11] also observed that there exists a beautiful identification between Young tableaux distribution and eigenvalue distribution for all the phases of the theory, where h ± are solutions of the following equation (1 + 2β 1 cos θ) 2 = π 2 ρ 2 (θ). (3.15) This identification allows one to write eigenvalue distribution and Young tableaux distribution in terms of a single constant distribution function ω(h, θ) given by equations (1.3) and (1.4). h ± (θ) defines the boundary of the distribution. It is actually the shape (or boundary) of the region R which contains all the information about two distributions for a given phase of the model. Thus, a relation between h and θ defines disjoint islands (phase space distribution/droplet) in two dimensions for different phases of the UMM under consideration. Different phases are distinguished by different topologies of droplets. Constant phase space distribution is very similar to Thomas-Fermi (TF) model [27,28] at zero temperature.
The identification, observed in [10,11], between two distribution was accidental. It was not known if such identification holds for a generic model like (3.1). The main obstacle to apply the idea of [10] to a generic UMM was finding character in terms of h i s in presence of arbitrary cycles. Presence of TrU n term in partition function requires to calculate character in presence of n-cycles. Hence, finding u(h) by extremizing the effective action for a generic unitary matrix model was a challenging problem.
In [12] it was shown that, although it is difficult to find u(h) for a generic plaquette model (3.1), but one can still find droplets for different phases of the theory and from the shape of these droplets one can reconstruct dominant Young diagram.

Frobenius formula for character
Let C( k) denotes the conjugacy class of S K which is determined by a collection of sequence of numbers 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 . The character corresponding to a conjugacy class C( k) of S K for a representation characterized by λ is given by famous Frobenius formula and h i 's are give by equation (3.8).

Extremization of partition function
Promoting N real variables x 1 , · · · , x N in (3.17) to N complex variable (z 1 , z 2 , · · · z N ) character can be written as, The contour is taken around origin. Note that, in the above expressions for character (equation 3.17 or 3.20) variables x i 's or z i 's are auxiliary variables. They do not carry any physical meaning, as of now.
Integrand in the above expression can be exponentiated and written as, Therefore, in large N limit character is given by, where action S χ is given by and N 2 K = K = total number of boxes in a representation.
Using the expression for character given in equation (3.22), partition function (3.7) can be written as, where, −S total h, z, w, k , l , t, s = n k n 1 + ln β n Z n nk n + l n 1 + ln β n W n nl n (3.25) Varying effective action with respect to z(x) we find 6 , where, Varying effective action with respect to h(x) we find, Variation with respect to k n and l n provides, β n Z n nk n = e −itn , and β n W n nl n = e −isn . (3.29) Since, in equation (3.24) z(x) and w(x) vary over some contour around zero, we take these contours to be unit circle about the origin and hence a consistent solution to above equations are, Therefore, Evaluating the partition function (3.24) on these equations we find, This partition function is exactly same as the partition function of one plaquette model written as an integral over eigenvalues of unitary matrices. Hence, we identify that the auxiliary variables we used to write the character of symmetric group with eigenvalues of corresponding unitary matrix model involved. Imaginary part of saddle point equation (3.32) becomes equation for eigenvalue density, On the other hand, the real part of equation (3.32) provides a relation between h and θ, However, this equation is not complete. It is easy to check that total action where, f (θ) is an even function of θ. Hence, the most generic form of h(θ) is given by, Form of f (θ) can be fixed from topological structure of droplets for different phases [12] and is given by f (θ) = πρ(θ). Finally we have, This equation defines the shape of constant droplets in (h, θ) plane. This equation is generalization of the relation given in equation (3.14). Thus we see that for a generalised plaquette model all the phases of the theory can be characterised in terms of topologies of different droplets in two dimensional plane spanned by h and θ.

Two dimensional droplets for generic unitary matrix model
We start with an action which contains maximum three traces, When n 1 (or n 2 ) is zero we get back one loop (zero coupling) action. The above action is invariant under U → U † if a n 1 n 2 = a −n 1 −n 2 , and can also be written as, Following the same analysis as discussed in previous section we find a saddle point equation for auxiliary variables appearing in the Frobenius formula as (see appendix A for details), where, ρ(θ) is defined in equation (3.32). Imaginary part of this equation gives eigenvalue equation Integrating over h(x), the partition function can be written as 2N 2 a n 1 n 2 ρ n 1 ρ n 2 ρ n 1 +n 2 where, a n 1 n 2 ρ n 1 ρ n 2 ρ n 1 +n 2 n 1 cos n 1 θ ρ n 1 + n 2 cos n 2 θ ρ n 2 + (n 1 + n 2 ) cos(n 1 + n 2 )θ ρ n 1 ρ n 2 .
(3.45) This is a generalization of equation (3.38). It is easy to check that all the expressions matches with those of plaquette model for n 2 = 0 with 2n 1 a n 1 0 ρ n 1 = β n 1 .
Although here we considered three trace terms only, a generalisation to arbitrary model is straightforward. However, algebra will be little difficult.

Unitary matrix quantum mechanics and phase space distribution
h ± in generalized boundary relation (3.38 or 3.44) is a solution of the following quadratic equation such that ω(h, θ) = 1 for h − (θ) < h < h + (θ) and zero otherwise. Eigenvalue distribution, by construction, can be obtained by integrating out h for a given θ (equation 1.3) and is given by Wilson loops are given by The function S(θ) can also be written in terms of phase space geometry The droplet picture of different large N phases are similar to droplets of Thomas-Fermi (TF) model at zero temperature. Fermi distribution at zero temperature is given by where µ is chemical potential and h(p, q) is single particle Hamiltonian density. Comparing our phase space distribution (4.2) with TF distribution we find the Hamiltonian density is given by Total Hamiltonian can be obtained by integrating h(h, θ) over the phase space Integrating over h we find Hamiltonian as a functional of ρ(θ) This is TF energy functional. As a consistency, one can check that extremization of this Hamiltonian with respect to ρ gives eigenvalue distribution (4.3).
The Hamiltonian can also be written as for some V ef f (θ). For example, for simple

Unitary matrix quantum mechanics
We start with a unitary matrix quantum mechanics This matrix model is equivalent to a theory of free fermions on circle with Hamiltonian [1,2] and with eigenvalue density ∼ Ψ † Ψ. The matrix model (4.11) can also be described as a real bosonic field theory of eigenvalue density ρ(θ, t) on S 1 [15,18,19,29]. Corresponding Hamiltonian is given by, where, π(θ, t) is momentum conjugate to ρ(θ, t). Equations of motion for collective fields are given by, where v(t, θ) = ∂ θ π(t, θ). These are coupled, non-linear partial differential equation and hence difficult to find a solution in general. Therefore, one can translate the problem into free fermi language [3]. At large N , time evolution of the fermionic system can be described in terms of time evolution of fermi surface in phase space. Denoting phase space density by (p, θ) the Hamiltonian is given by, where dθW (θ)ρ(θ, t) = W (U (t)) and p is momentum (conjugate of θ). This Hamiltonian has the same form as (4.10) except the fact that (p(t), θ(t)) in (4.16) are functions of "t". There is a one to one correspondence between phase space variables and collective field theory variables. Eigenvalue density and corresponding momentum are given by [3] ρ(θ, t) = 1 2π dp (p, θ), ∂ θ π(θ, t) = 1 2πρ dp p (p, θ).
Assuming p(t, θ) ranges from p − (t, θ) to p + (t, θ) for a fixed θ we land up with following dictionary between bosonic and fermionic variables Using this dictionary the Hamiltonian (4.16) reduces to Hamiltonian (4.13). Also, using this dictionary one can show that coupled non-linear differential equations for ρ(t, θ) and v(t, θ) take simpler form in fermionic language These are decoupled first order quasi-linear partial differential equations. p ± (t, θ) defines fermi surface in (p, θ) plane. The above equation (4.19) determines evolution of fermi surface with time. Thus solving the field theory equations of motion (4.14) is equivalent to solve for upper and lower fermi surfaces in fermionic picture. In either case, one needs to provide an initial data on a constant time slice in (t, θ) plane. After that the problem reduces to a Cauchy problem.
Existence of unique solution depends on the geometry of initial data curve.
Inverting the dictionary (4.18) we find These equations for fermi surface have similarity with generalized boundary relation (3.44) in (h, θ) plane. Therefore, on a constant t slice 7 , one can match fermi surfaces p ± (t, θ) with boundary of large N droplets (3.44) of a generic UMM (2.1). In other words, large N phase space distribution can be used as initial (or terminal) data for the Cauchy problem, In the bosonic language this is equivalent to v(t 0 , θ) = S(θ) and ρ(t 0 , θ) = ρ(θ). A general phase space distribution (3.44) therefore can be thought of as initial (or terminal) geometry of an underlying quantum mechanics problem. h which was related to number of boxes in dominant Young representation, is indeed behave like momentum conjugate to θ. In other words momentum p, conjugate to θ in UMQM has a meaning in terms of boxes in Young diagram. The initial data curve (4.21), we have chosen, is non-characteristic and hence the solution to equation (4.19) is unique.
Thus we see that for any large N phase of a generic UMM (2.1) it is possible to construct an effective UMQM such that phase space distribution of matrix quantum mechanics matches with that of UMM on a constant time slice. The effective matrix quantum mechanics has a potential W (θ) = V ef f (θ). The Hamiltonian (4.16) evaluated on t = t 0 slice boils down to (4.10) up to a constant piece.
Unitary matrix model (2.1) shows a rich phase structure at large N . The model has many possible phases. As we vary the parameters of the theory the system undergoes a phase transition. It would be interesting to understand the phase transition in UMM in terms of dynamics of effective matrix quantum mechanics. One can allow the parameters of W (θ) in UMQM to depend on t such that dynamics of UMQM captures phase transition of UMM. We choose UMQM such that on a constant time slices W (θ) matches with V ef f (θ) of UMM for different phases. As depicted  figure 1 we can start with a particular phase P 0 of UMM at t = t 0 with W (θ(t 0 ), t 0 ) = V P 0 ef f (θ) and use droplet of P 0 as initial data for the Cauchy problem. Then we let the system evolve with time and some later time t = t 1 > t 0 we have W (θ(t 1 ), t 1 ) = V P 1 ef f (θ) and shape of fermi surfaces matches with droplets of P 1 phase. This will help us in understanding the question how dynamics of strongly interacting gauge theories are encoded in dynamics of free fermions. A work in this direction was considered in [15] for a simpler model 8 , where they started with a no-gap distribution at t → −∞ and showed that time evolution takes this phase to a one-gap phase at t → ∞. However, for a generic model, we need to thoroughly study the solutions of quasi-linear equations 4.19. We leave the issue for future.
For a unitary matrix quantum mechanics, phase space has θ → −θ symmetry. Using this fact it is possible to show that phase space area covered by fermi surface is a constant of motion. Area covered by fermi surfaces at a time t is given by (4.23) Using equation (4.19) and the fact that p ± (t, θ) and W (θ) are symmetric functions of θ it is easy to verify that d dt Thus, phase space area is preserved during time evolution. Its only the shape which changes during evolution.

Momentum distribution and Young tableaux
Since (h, θ) are canonical conjugate of each other, integrating out θ in phase space gives us momentum distribution for free fermions under consideration. Thus, momentum distribution is actually encoded in the arrangements of boxes in the most dominant Young tableaux. To obtain Young distribution function u(h) we use the definition (1.4). Unlike eigenvalue distribution, finding Young tableaux distribution directly from partition function is a formidable task for a generic matrix model. But we assume that the definition (1.4) holds in general. This definition boils down to our earlier identification [10,11] πu(h) = θ when θ has unique solution for a given h in equation (4.1). In general, θ may have multiple solutions for a given h. In that case u(h) is given by total length of different segments of white arc as shown in figure 2. This equation is, in general, difficult to solve and find u(h) for any arbitrary set of parameters {β n }. In section 5 and 6 we use this relation to find the Young tableaux distribution which dominates the partition function in large N limit.
Young distribution captures information about momentum distribution of N fermions. Therefore, large N state of the system can be found from Young distribution. Since fermions have momentum h 1 , · · · , h N a large N state of the system can be recognised as |h 1 , · · · , h N in momentum space 9 [2]. Thus a large N state correspond to geometry of Young diagram. Time evolution of quantum fluctuations about this large N states can be studied using underlying UMQM.
Finding dominant Young distributions is also important, as local dynamics of holographic geometry is encoded in Young diagram of dual theory [31]. Supersymmetry plays an important role here, as weak coupling results are protected. Hence, by determining phase space distribution and corresponding Young diagrams for a supersymmetry theory one can reconstruct the dual holographic geometry from that.

Young distribution for β 1 − β 2 model
In this section we discuss an example where we have non-trivial, stable two gap phases. This example will help us to understand how to reconstruct Young diagrams from phase space for different large N phases of the theory.
To have a two gap solution we need the potential V (θ) in equation (2.5) to have two minima between −π and π. This is achieved by considering a plaquette model with β 1 , β 2 = 0 and β n = 0 for n ≥ 3. The action S(U ) is given by Phase structure of this model was discussed in [20] and [21]. The model exhibits "five" different phases depending on values of parameters β 1 and β 2 : one no-gap phase, two one-gap phases (mirror reflection of each other) and two two-gap phases. A detailed derivation of eigenvalue distribution for different phases can be found in [20,21]. Here, we present a different and independent method (solving Dyson-Schwinger equation) to reproduce the results (eigenvalue distribution) of [20,21]. Then we obtain droplets for different phases of the theory. Finally, we discuss how to reconstruct Young tableaux distributions for five different phases of this model from the geometry of droplets.

The resolvent
Eigenvalue distribution can be obtained from analytic properties of resolvent which satisfies an algebraic equation in large N limit [32]. This method is quite powerful as one does not need to solve an integral equation for eigenvalue density. We briefly review the method here.
"Resolvent" for any generic UMM (2.1) is defined as, where z is a complex variable. Expanding the right hand side of (5.2) One can see that the resolvent is a generating functional of Wilson loops in the planar limit.
Since N −1 TrU k lies between −1 and +1 the right hand side of equation (5.3) is convergent for |z| < 1. Therefore, R(z) is an analytic and holomorphic function in the interior of the unit disk.
In a similar way, one can show that the function R(z) is also analytic for |z| > 1 with R(z) → 0 for |z| → ∞. However, R(z) is discontinuous on |z| = 1. Inside the unit circle one can also write that, Hence we have the following property of R(z) At large N the resolvent satisfies a quadratic equation (Dyson-Schwinger equation) which can be solved algebraically [11]. The solution, for a particular model (5.1) under consideration, is given by (see [11] for details) where F (z) has the following form, Spectral density or eigenvalue density is related to discontinuity of R(z) on unit circle. Eigenvalue distribution function (ρ(θ)) is given by, We see that F (z) contains some quantities R (0) and R (0) which can be calculated from (5.6). Therefore one has to consistently solve these equations to find R (0) and R (0) for different phases of the model and then ρ(θ) from (5.8). However, it is difficult to solve equation (5.6) and (5.8) consistently. Rather, we shall take an ansatz for F (z) for different phases of the theory and fix the form of F (z) from analytic properties of R(z) : 1. R(z) is analytic inside the unit circle. In equation (5.6) we see that R(z) has 1/z and 1/z 2 terms which are singular at z = 0. Therefore, F (z) should have poles of order one and two at z = 0 to make R(z) analytic inside. F (z) can not have zeros inside (or outside) the unit circle. In that case F (z) would have branch points inside or outside the circle. However, F (z) can have zeros on |z| = 1 line.
2. For gapped solution F (z) should have branch points on unit circle. Also F (z) is invariant under z → 1/z.

No-gap solution
This is the most trivial case where there is no branch cuts in resolvent R(z). Hence, F (z) is a perfect square. Thus we have, (5.10) The eigenvalue distribution for this phase is given by (from equation 5.9) (1 + 2β 1 cos θ + 2β 2 cos 2θ) (5.11) for 0 ≤ θ ≤ 2π.

One-gap solution
For one-gap solution F (z) has a branch cut on |z| = 1. Therefore the most generic ansatz for F (z) for one-gap phase is given by, where α, β, γ are constants. From analyticity and normalization condition of R(z) we find Hence, the spectral density for one-gap solution (A 1 phase) is given by Eigenvalues are symmetrically distributed about θ = 0 from −θ 1 to θ 1 . One-gap solution is valid in the following region in (β 1 , β 2 ) plane bounded by the curves [20] (β 1 + 2β 2 ) 2 − 6β 2 = 0, for β 2 ≥ 1 6 , and extends in the direction of positive β 1 . There exists a common boundary between this phase and zero-gap phase and is given by F (z) (given by equation 5.13) for one-gap solution has a branch cut on unit circle about z = 1 point (θ = 0). If we replace z by −z then the branch-cut changes its position from right to left side on the unit circle. However, this is also a valid ansatz for F (z) for one-gap solution (B 1 phase). For this ansatz the parameters α, β, γ take the following values to make R(z) analytic inside.
which in turn gives the spectral density as with Since action (5.1) is invariant under β 1 → −β 1 and θ → θ + π therefore, it is expected that there exists a mirror image of A 1 phase in negative β 1 direction (B 1 phase). Eigenvalues, for this branch, are distributed symmetrically about θ = π from π + |θ 1 | to π − |θ 1 |. Similarly, the common boundary between this phase and zero-gap phase is given by A transition along line (5.18) or (5.22) between zero-gap and one-gap (or its mirror phase) is a third order phase transition. This is a generalization of Gross-Witten-Wadia transition [33,34]. In subsection 5.3 we shall see how phase space topology changes as we go from one phase to another.

Two-gap solution
To have a two-gap solution one should take an ansatz for F (z) such that there is at most a degree four polynomial inside the square root. In this case there are two possible ways of taking ansatz for F (z) which give two different two-gap phases namely A 2 phase and B 2 phase.
The simplest ansatz that one can take is -25 -α, β, σ, δ are constants. Using the analyticity and normalization condition of R(z) we find, (5.24) However, it is not possible to fix all the four parameters from the properties of R(z). We keep β undetermined. The eigenvalue distribution in this case is given by, A 2 phase has eigenvalue distribution symmetric about horizontal axis. Validity of this phase has been shown in figure 3. Eigenvalues are distributed between −θ 1 and θ 1 and then π − θ 2 to π + θ 2 .
Eliminating the undetermined parameter β from equations (5.26) one can write a constraint equation relating cos θ 1 and cos θ 2 as The undetermined parameter β can be fixed by imposing an extra condition that total number of eigenvalues between θ 1 and θ 2 is zero A second order phase transition between A 2 phase and A 1 phase (or B 1 phase) occurs along the lines [20] β 2 = 1 4 cos 3 θ 1 cos 1 2 θ 1 sin 2 θ 1 − ln 1 + cos 1 2 θ 1 sin θ 1 , The second possible ansatz for F (z) for a two gap solution is given by, with, The spectral density for this case is Eigenvalues are distributed between θ 1 and θ 2 , and −θ 1 and −θ 2 . This phase exists for β 2 ≤ − 1 3 . A parabola (−β 1 + 2β 2 ) 2 + 4β 2 = 0 separates this phase from one-gap phase. The point (β 1 = 0, β 2 = − 1 2 ) is a triple point in the parameter space. See figure (3) for details.

Droplets and phase transition in β 1 − β 2 model
We shall discuss about geometries of droplets for different phases of the model. The boundary relation (3.38) for this model is given by where h ± (θ) are solutions of h 2 − (1 + β 1 cos θ + β 2 cos 2θ) h = π 2 ρ 2 (θ). (5.35) There are five possible droplets with different topologies corresponding to five different phases. In figure 3 we draw phase space distributions for different phases. The diagram in the middle is parameter space (β 1 − β 2 plane) of the theory under consideration. We divide this parameter space into different parts according to validity of different phases. See [20] for more detail 10 . Since, h > 0 (being number of boxes) and −π < θ ≤ π, we plot phase space distribution in polar coordinates. The black regions correspond to density ω(h, θ) = 1. We note that topologies of droplets are different for different phases. Origin (h = 0) is inside the distribution for no-gap phase (since h − (θ) is zero for all θ). For one-gap and two-gap phases, on the other hand, origin is outside the distribution (both h + and h − are non-zero). We consider this as a topological property of droplets. One should also note that all the distributions are symmetric about real axis (θ = 0 line). This is because of UMM has U → U † symmetry. While deforming the droplets Unlike no-gap phase, here distribution is non-zero for θ between −θ 1 and θ 1 . Distribution at the bottom correspond to B 2 phase. Here, distribution is non-zero for θ between two different ranges. Distribution on left middle is a mirror image of A 1 phase. This is B 1 phase. Finally, on top left we draw distribution for the second two-gap phase (A 2 ). Here also phase space density is non-zero for two different ranges of θ. One should note that all the distributions are symmetric about θ = 0 line. This is because of U → U † symmetry of the model.
(by changing values of parameters) this symmetry is always preserved. Therefore, one can not continuously deform the droplets of B 2 phase to bring them in the shape of droplets of A 2 phase. This implies, there is no transition possible between B 2 and A 2 phases, which is true for this model. Boundary of a droplet crossing the origin implies (according to our definition) change of topology, which can only happen at the time of phase transition. We see that at A 0 -A 1 transition point (generalized Gross-Witten-Wadia transition), left boundary of the droplet touches the origin ( figure 4(b)). From this point as we move towards A 0 , the origin goes inside and as we move towards A 1 the origin goes outside. Similarly, cutting a single droplet into two (or gluing two droplets to make a single one) also corresponds to change of topology and occurs at the time of phase transitions. For example, starting from A 1 as we move towards B 2 phase Since boundary of droplets are fermi surfaces in underlying MQM, the above diagrams show how shape of fermi surfaces change during phase transition.

Reconstruction of Young diagrams from droplet geometry
The phase space description is a powerful approach as the shape of different droplets not only contains information about eigenvalue distributions but it also contains information about dominant representations for different phases of the theory. In this section we discuss how one can construct dominant Young diagrams for different phases.
To construct Young tableaux distribution from geometry of droplets we use equation (1.4). [10,11], considered plaquette model with β n = 0 for n ≥ 2 and observed that while solving boundary relation (3.15) for a given h there exists a unique value of θ and hence Young distribution was given by πu(h) = θ(h) for all phases. However, solution for θ does not remain unique when we turn on β 2 in plaquette model (5.1). Here we see that boundary relation (5.35) for a fixed h has maximum two possible solution for θ. Therefore, the simple identification between eigenvalue and Young density (πu(h) = θ) does not hold in this case. However, we consider equation (1.4) as a fundamental definition for u(h) and construct Young distributions from phase space for different phases.

No-gap phase
For no-gap phase ρ(θ) is given by (5.11). Therefore the Fermi surface for this phase is given by, This defines a droplet about the origin. For β 1 > |4β 2 |, θ has a unique solution (we denote that byθ 1 (h)). In that case u(h) is given by as before.
For β 1 < |4β 2 |, θ has multiple solutions depending on h. In figure 5 we plot h vs θ for nogap phase and for β 1 < |4β 2 |. The left and right plots correspond to positive and negative β 2 respectively. Following definition (1.4), Young distribution for positive β 2 is given by For negative β 2 the distribution is given by, distribution. This correspond to sudden change in distribution function for boxes in a diagram. As in [10], Young tableaux for no-gap phase has finite number of rows empty. n(x) = 0 for 1 < x < 1 − h min . After that number of boxes starts increasing and reaches a maximum value n(0) = h max − 1. The distribution of boxes from x = 0 to x = 1 − h min is governed by u(h) given by equation (

One-gap phase
Eigenvalue distribution for one-gap phase is given by equation (  (a) β 1 , β 2 > 0 or β 2 < 0, Here we see a kink in the distribution of boxes. u (h) is discontinuous.

Two-gap phase
There are two possible two-gap phases : B 2 and A 2 .
B 2 phase : For B 2 phase h vs θ and corresponding u(h) vs. h are given in figure 11. A typical Young tableaux for B 2 phase is similar to that of A 1 phase without any kink.   there is no kink in Young distribution. As we move towards left (or right, depending on sign of β 1 ) from β 1 = 0 line, two kinks are developed at the top and bottom of a tableaux. These two kinks approach towards each other as we move towards A 1 (or B 1 phase). Finally on A 2 − A 1 (or A 2 -B 1 ) transition line they meet and disappear.
On the other hand when we approach towards A 0 phase from any arbitrary point in A 2 region, we find that the lower kink starts moving towards h min and h min also starts decreasing. On A 2 − A 0 transition line, the lower kink merges with h min = 0. This is exactly the A 0 phase on A 0 − A 2 transition line.

Chern-Simons on S 3
Chern-Simens (CS) theory plays an important role in different branches of physics and mathematics. CS theory is also very extensively studied in string theory. Partition functions of open and closed string theories can be related to partition functions of CS theory on compact manifolds [35][36][37].
CS theory on S 3 and on other three manifolds can be written in terms of an "exotic" matrix model [37,38]. The matrix model can be solved at large N limit to obtain free energy of open or closed topological string theory. Matrix model analysis of CS theory is very useful to obtain informations about topological string theory. Partition function for CS on S 3 can be written as an integral over hermitian matrix ensemble with logarithimic potential [38]. At large N , phase of CS theory on S 3 is given by distribution of eigenvalues of hermitian matrices on a straight line. It turns out that topological phase of the theory corresponds to a "gapped" distribution of eigenvalues. In this section we study large N phase of CS on S 3 in terms of distribution of eigenvalues of unitary matrices. This is required to obtain phase space distribution of this theory. The phase space picture will be helpful to understand underlying quantum mechanical description of string degrees of freedom using dualities.
We show that a droplet picture exists for CS theory on S 3 for any value of 't Hooft λ = N g s . Since, CS on S 3 is topological there exists no phase transition in the theory as we change the coupling constant. We shall see that this property is reflected in phase space distribution function, i.e. as we change 't Hooft coupling λ the droplet changes its shape keeping the topology intact.
Chern-Simons action with gauge group G on a generic three-manifold M with level k is given by Following [39] we can write the partition function for CS theory on the three-sphere S 3 with string coupling constant Where ρ is the Weyl vector of SU (N ) and α is the positive roots of SU (N ) ⊂ U (N ) defined as α ij = µ i − µ j corresponding to the cartan subalgebra.
Using Weyl formula for character and properties of Jacobi's theta function one can write the partition function for CS theory on S 3 as a UMM [40] Z , (6.4) and Thus we see that the partition function is exactly same (upto an overall normalization) as that of a single plaquette model (equation 2.11) with the coefficients 1/n[n].
In large N limit |g s | becomes small, one can expand the exponential and find [n] = g s n. Hence, we haveβ n = 1 N g s n in large N limit. (6.7) We shall study this model in large N limit.

Potential and eigenvalue distribution
Since the partition function (6.4) (action as well as measure) is invariant under unitary transformation, one can go to a diagonal basis where U = {e iθ i }, i = 1, · · · N , θ i 's are eigenvalues.
In this basis CS action is given by Li 2 (e iθ i ) + Li 2 (e −iθ i ) .

(6.8)
Hence the partition function is Two important things to note here.
• The second term in exponential in equation (6.9) is respponsible for strong repulsion among the eigenvalues. This terms comes from Haar measure. The model will have a stable nontrivial eigenvalue distribution if the first term (Li 2 (e iθ i ) + Li 2 (e −iθ i )) has minima between −π and π. This term, in fact, gives a harmonic oscillator potential about θ i = π. However, the coefficient of this potential term is g s = 2πi/k. If we study saddle point equation with this potential then we land up with a complex eigenvalue distribution. Therefore to get a stable real and positive eigenvalue distribution we replace g s → ig s . In that case the potential becomes a harmonic oscillator potential about π with real positive coefficient.
The problem with complex potential is related to analyticity of resolvent. Resolvent computed for 6.4 turns out to be non-analytic inside and outside the unit circle. Thus to make resolvent analytic in complex z plane we need to replace g s → ig s . Since free energy is an analytic function of g s an imaginary rotation in g s does not affect the calculation.
• From equation (6.8) we see that the potential has a minima about θ i = π. However, we are following a convension where everything is symmetric about θ = 0. Thus to make the potential symmetric about θ = 0 we give a π shift to all eigenvalues. This is equivalent to U We shall work with this partition function.
In large N limit, using saddle point analysis one can find the eigenvalue distrubution which dominates the partition function. The eigenvalue distribution function satisfies an integral equation, which is in general difficult to solve to find a solution. Here we follow the discussion in section 5.1 to find the eigen value distribution.

Resolvent and Eigenvalue distribution
Resolvent for partition function 6.10 is given by Note that the second term in R(z) is antisymmetric under z → 1/z. Hence, performing resummation over appropriate regions we can write R(z) as Since resolvent is analytic at z = 0 by construction, F (z) must also have a logarithmic singularity at origin with equal and opposite strength. This may happen when F (z) is a perfect square, i.e In this case all the singular terms in equation (6.13) are cancelled. Using equation (5.9) we find spectral density is given by This is no-gap solution (phase) of CS on S 3 . However, ρ(θ) becomes negative for some values of θ, hence this is not a physical solution or phase. We need to look for gapped phase of the theory.
To find a more general class of solution we take the following anstaz for F (z) where f (z) and g(z) are polynomials of positive powers and Γ is some arbitrary constant 11 . Unknown functions g(z) and f (z) can be obtained form analytic property and normalization 11 One can also take the ansatz for F (z) to be F (z) = β ln f (z) − g(z) + Γ. (6.16) In this case it turns out that the functions f (z) and g(z) are not positive polynomials.
condition of R(z). We also need to use the fact that F (z) is symmetric under z → 1/z which implies that the function g 2 (z) − f 2 (z) has to be an even order polynomial. All the zeros of g 2 (z) − f 2 (z) are located on |z| = 1 line. Which means branch cuts are exactly on unit circle in z plane. Depending on degrees of the polynomials g(z) and f (z) we have one gap or multi gap solutions.

One Gap Solution
One gap phase has two branch points on unit circle (therefore one branch cut), hence g(z) should be a order one polynomial. We take g(z) = 1 + z. From invariance of F (z) under z → 1/z we have, Analyticity of R(z) at z = 0 implies that β = 1/λ. Finally, from the normalization of R(z) (i.e. R < (0) = 1) we have α = 4e −λ and Γ = 0. Therefore, Hence, eigenvalue distribution for one-gap phase is given by, (6.19) Since ρ(θ) ≥ 0, this implies eigenvalues are distributed over the range −2 cos −1 e −λ/2 < θ < 2 cos −1 e −λ/2 . (6.20) This analysis is similar to analysis of [41]. However, [41] studied distribution of eigenvalues of hermitian matrices. The gap depends on value of λ, it increases (spreading of eigenvalue distribution decreases) as we decrease λ. At very small value of λ all the eigenvalues are concentrated about θ = 0 point. This can be understood from the fact that at small values of λ attractive part of the potential in equation (6.10) becomes infinitely strong and hence all the eigenvalues are localised at θ = 0.

Two Gap Solution
Continuing with the same sprit, one can construct the function F (z) for two-gap solution. In this case g(z) is a second order polynomial.
From the properties of R(z) we find, β = 1 2λ , α = 4e −2λ and Γ = 0. (6.22) The parameter γ can not be fixed from the analyticity or normalization conditions. Hence it is an one parameter family of solution at this point. The eigenvalue distribution for this branch is given by, tanh −1 4 (cos 2 θ/2 − cos 2 θ 1 /2) (cos 2 θ/2 − cos 2 θ 2 /2) 2 cos θ + γ , (6.23) where Thus we see that one can also have a two-gap solution for any value of λ. However, computing free energy of the system one can show that the two gap solution is not thermodynamically favoured i.e. free energy is always greater than that of one gap solution for all values of λ [42]. This can also be understood from the potential of the matrix model. Since the attractive part of potential has only one minima about θ = 0 therefore only one-gap solution is a thermodynamically stable solution. However, this kind of solution appears when we consider CS theory on lens-space [43] or ABJM theory [44].

Phase space distribution
Having understood eigenvalue distribution for CS theory on S 3 we now proceed to find the most dominant representation for the most dominant phase of CS theory. The distribution of boxes in the most dominant Young diagram can be obtained from the boundary relation.
The boundary relation for a generalised plaquette model is given by equation (3.38) h ± (θ) = 1 2 + n β n cos (nθ) ± πρ(θ). Eigenvalue distribution is given by equation (6.19). Hence boundary relation (Fermi surface) is given by (6.27) Phase space distribution for CS on S 3 has been plotted in figure 14. The droplets have the (a) λ = .9 (b) λ = .01 shape of "kidney" with a cleavage on left. The origin is always outside the distribution (character of gapped solution). As we increase values of λ, the droplet moves towards the origin also the cleavage becomes sharper. We also see that in valid ranges of λ phase space distribution changes its shape without crossing the origin. This implies that CS theory on S 3 is topological. There is no phase transition in the theory.

Young tableaux distribution
Since eigenvalues are identified to density of Young tableaux : πu(h) = θ, one can invert the boundary relation (6.27) to obtaing distribution of boxes in the most dominant phase of CS theory for a given λ. A typical relation u(h) and h and corresponding Young diagram have been plotted in figure 15. Minimum (h min ) and maximum (h max ) values of h correspond to number of boxes in the last row and first row of Young diagram respectively and are given by  In small λ limit we see that both h max and h min are very large and difference between them goes as 2/ √ λ.

Discussion and outlook
In this paper we show that large N generic phases of a generic unitary matrix model can be described in terms of geometry of two dimensional droplets in (h, θ) plane, where h is related to number of boxes in Young diagram and θ is eigenvalue of unitary matrices. We also see that these droplets are similar to phase space distribution of unitary matrix quantum mechanics on a constant time slice. Hence, we argue that for a given unitary matrix model it is possible to find an underlying quantum mechanics such that different large N phases of unitary matrix models (and hence interacting gauge theories) are captured on different time slices of the quantum system.
This observation allows us to identify momenta of fermionic system with number of boxes in a Young diagram. Distribution of boxes in Young diagram captures information about the distribution of momenta of N fermions. We show that, Young distribution for different large N phases can be computed from the geometry of two dimensional droplets. We explicitly work out two examples : β 1 -β 2 model and Chern-Simons theory on S 3 .
Phase space description of gauge theories can be used to provide a dual quantum mechanical description of string degrees of freedom, since different large N phases of gauge theory are dual to classical solutions of string theory in AdS space. Therefore, one can extend this idea for a supersymmetric QFT. One can consider ABJM theory [45], supersymmetric CS theory in 2 + 1 dimensions. Partition function of ABJM can be calculated exactly for all values of coupling constant. There exists an interpolating function between weak and strong coupling [44,46]. Hence, it would be interesting to write ABJM partition function in terms of a unitary matrix model and find two dimensional droplets associated with large N phases of the theory. One can interpolate this droplet picture to strong coupling side to explore properties/geometries of dual string/M-theories in the spirit of LLM [47]. This will also serve as a precision test of AdS 4 /CF T 3 . Understanding underlying matrix quantum mechanics for CS theory coupled with matter fields [48][49][50][51][52] would also be an interesting problem to look at from the point of [53].
One should note that the inviscid Burgers' equation 4.19 is a generic equation for random matrix theories. In the continuum limit of QCD [54], first observed that spectral density of the Wilson operator follows this equation as a consequence of Makeenko-Migdal equation and proposed a order(gapped)-disorder(gapless) transition for 2D YM at large N . [55] Also came across this same quasi-linear equation while dealing with the continuum version of the loop equation for QCD 2 . Starting with the Burgers' equation [56,57] supplemented the proposal of [54], and demonstrated the emergence of a spectral shock wave of Wilson loop eigenvalues at the large N c limit. Recently [58] extended this study numerically for d = 4. As matrix integrals can be thought of as an alternate way of writing the combinatorics of maps, this kind of connections between random matrix theory and loop equations can also be extended for combinatorial problems. To further explore this combinatorial avenue to understand the identification of loop equations with Tuttes equations (recursively deleting or contracting edges) and furthermore to grasp the current status of "random matrix" approach to 2D quantum gravity one useful resource would be [59].
On a different ground, UMM has plenty of applications in mathematics specially in number theory. Hilbert and Póyla speculated that there is a spectral origin for non-trivial zeros of Riemann zeta function. It was observed in [60] that if the non-trivial zeros are analogous to EV of some Hamiltonian, then the corresponding conjugate variables, time period of some primitive orbits, are related to prime numbers (proportional to ln p where p is prime number). Hence, if EV distribution of a UMM captures information about non-trivial zeros, then we expect that corresponding Young distribution function captures information about prime number distribution [11]. It would be nice to understand if there is any quantum mechanical system which captures information about these two distributions.

Acknowledgement
We would like to thank Rajesh Gopakumar for many helpful discussion. We are grateful to Shiraz Minwala and Gautam Mandal for their useful comments on our work. We are grateful to Debashis Ghoshal for innovative discussion. AC would like to acknowledge the hospitality of IISER Pune where part of this work is done. SD acknowledges the Simons Associateship, ICTP. SD also thanks the hospitality of ICTP, Trieste and JNU where part of this work has been done. Work of SD is supported by DST under a project with number EMR/2016/006294. PD is partially supported by research grant number 5304-2, "Symmetries and Dynamics: Worldsheet and Spacetime" from CEFIPRA/IFCPAR. Finally, we are grateful to people of India for their unconditional support towards researches in basic sciences.

A.1 Extremization of partition function
In large N limit we define, a n 1 n 2 = N 2 a n 1 n 2 ; k n 1 n 2 = N 2 k n 1 n 2 ; l n 1 n 2 = N 2 l n 1 n 2 ; K = N 2 K (A.11) and using the expression for character given by equation ( dl n 1 n 2 and −S eff (h, z, w) = ∞ n 1 ,n 2 =0 k n 1 n 2 1 + ln a n 1 n 2 Z n 1 Z n 2 W n 1 +n 2 k n 1 n 2 +l n 1 n 2 1 + ln a n 1 n 2 W n 1 W n 2 Z n 1 +n 2 l n 1 n 2 − 1 0 dx h(x) ln (z(x)w(x)) + 1 2 (n 1 + n 2 )(k n 1 n 2 + l n 1 n 2 ) − K ) (A. 13) where, Z n = z n (x)dx and W n = w n (x)dx. (A.14) In the large N limit dominant contribution to partition function comes from extremum value of this effective action. Varying this effective action with respect to h(x) we find, Again, varying the effective action with respect to k n 1 n 2 and l n 1 n 2 one can write a n 1 n 2 Z n 1 Z n 2 W n 1 +n 2 k n 1 n 2 = e −i(t+s)(n 1 +n 2 ) = constant, a n 1 n 2 W n 1 W n 2 Z n 1 +n 2 l n 1 n 2 = e −i(t+s)(n 1 +n 2 ) = constant.
(A. 16) Finally, varying S eff (h, z, w) with respect to z(x) we find, k n 1 n 2 n 1 z n 1 (x) Z n 1 + n 2 z n 2 (x) Z n 2 + l n 1 n 2 (n 1 + n 2 )z n 1 +n 2 (x) Z n 1 +n 2 = 0. (A.17) A similar equation can be obtained when we vary the effective with respective to w(x) with z(x) replaced by w(x).
Since, in equation (A.12) z(x) and w(x) vary over some contour around zero, we take these contours to be unit circle about the origin and hence one consistent set of solutions to equations (A. 15 Therefore, Z n = W n (for all n) and k n 1 n 2 = l n 1 n 2 = a n 1 n 2 Z n 1 Z n 2 Z n 1 +n 2 .