Multiple phases and meromorphic deformations of unitary matrix models

We study a unitary matrix model with Gross-Witten-Wadia weight function and determinant insertions. After some exact evaluations, we characterize the intricate phase diagram. There are five possible phases: an ungapped phase, two different one-cut gapped phases and two other two-cut gapped phases. The transition from the ungapped phase to any gapped phase is third order, but the transition between any one-cut and any two-cut phase is second order. The physics of tunneling from a metastable vacuum to a stable one and of different releases of instantons is discussed. Wilson loops, $\beta$-functions and aspects of chiral symmetry breaking are investigated as well. Furthermore, we study in detail the meromorphic deformation of a general class of unitary matrix models, in which the integration contour is not anchored to the unit circle. The ensuing phase diagram is characterized by symplectic singularities and captured by a Hasse diagram.


Introduction
The study of the spectral and critical properties of models of random matrices has become a widely popular and interdisciplinary subject in this century, in great part due to the vast scope of fields where such models appear naturally and play a prominent role [1][2][3]. One among the possible many ways to highlight this relevance and versatility of random matrix theory is to simply point out the fact that the same model oftentimes appears in remarkably different contexts and, in addition, in a significant manner.
A paradigmatic example of this phenomenon could very well be the so-called Gross-Witten-Wadia (GWW) model [4][5][6]. Originally proposed in the study of gauge theory, it is ubiquitous and pivotal in many other areas, such as combinatorics, representation theory and spectral theory [1][2][3].
With this fact in mind, in this work we will study unitary matrix models, starting, precisely, with the specific case of a generalized form of the Gross-Witten-Wadia model. A possible interpretation of the model is as a one-plaquette model of two-dimensional lattice QCD with fermionic or bosonic quarks, which equivalently corresponds to a massive deformation of a model introduced by Minahan [7,8].
We will show that, while quite simple, this model retains several features of a sensible quantum field theory in the continuum. In turn, its simplicity allows us to exploit standard techniques from random matrix theory to characterize the theory at large N and suggests more general and deeper problems to consider. Some of them, we will already tackle here, by discussing at length the case of meromorphic deformations of unitary matrix model, as we explain below.
Before that, a rich phase diagram will be obtained and analyzed in detail. Phase transitions such as the ones we obtain in our analysis are relevant to the study of deconfinement transitions in QCD models in four dimensions [9,10] and in black holes physics [11][12][13]. More recently, this type of phase transitions has been argued to describe the critical behaviour of models exhibiting partial deconfinement [14][15][16][17][18].
We introduce the model in what follows, in Section 2, which includes a discussion on interpretations of the model, notation, relationship with other systems and an introductory discussion of its mathematical properties, including exact evaluations without scaling limits.
Then, the main results of the paper are presented and organized as follows: there are three main contributions, as far as new results are concerned. In Section 3, we fully characterize the rich phase structure of the unitary matrix model. In Section 4, we study Wilson loops in the same setting of Section 3 and discuss at length the physical interpretation of the phase transitions, including the role of instanton contributions.
Finally, in Section 5 we study, in a general framework and going beyond the specific model studied in the previous sections, the case where the integration contour is deformed in C * away from the unit circle. In spite of the vast body of results on random matrix ensembles, holomorphic matrix models [19] are arguably understudied.
The aim of Section 5 is to adapt the results on holomorphic matrix models to unitary matrix models. We do so for a very general set of unitary matrix models and, only as an illustrative example, we discuss the particular case of the holomorphic GWW matrix model. Non-traditional tools in this area, such as symplectic singularities and Hasse diagrams, are introduced here to fully understand the meromorphic models.
We conclude with possible avenues for further research in Section 6.

The model
In this section we present the model and study some of its exact features at finite N .
2.1. The model and its interpretations. Consider a one-plaquette model of two-dimensional lattice gauge theory [20,21] with gauge group U (N ) and K pairs of real fields, that can be either bosonic or fermionic. Each pair is called a flavour. We encode the choice of matter fields in the binary variable = +1 fermions, −1 bosons.
We must impose (anti-)periodic boundary conditions, so the discrete space-time is effectively reduced to a point with a loop attached to it, along which the gauge connection travels. Let m f > 0 be the mass of the f th flavour, and introduce the notation The partition function of this theory has the matrix model representation [7] Z ,K U (N ) (λ) = where λ ≡ N g 2 YM is the 't Hooft coupling for the bare gauge coupling g YM , U ∈ U (N ) is the plaquette gauge variable and dU is the normalized Haar measure on U (N ). Particularizing to the degenerate case with all equal masses, µ f = µ ∀f = 1, . . . , K, the partition function becomes: The second line is written in terms of the eigenvalues z j ∈ T of U ∈ U (N ).
The matrix model (2.1) is a massive deformation of the model studied in [22,23]. Besides, (2.1) is a generalization of the celebrated Gross-Witten-Wadia (GWW) model [4,5] by determinant insertions, and reduces to it for K = 0 or µ → ∞. For λ −1 = 0 the model (2.1) is a particular case of a correlator of characteristic polynomials in a Circular Unitary ensemble (CUE) [24], a fundamental object in random matrix theory with many applications [25,26]. Likewise, for this value of λ −1 = 0, it generalizes a matrix model that describes non-intersecting random walks [27], and has also appeared in gauge theory, for instance in [28] and later on in somewhat disguised forms. Still for this particular value, the model describes the probability that a random Young diagram has a value less than or equal to N [29].
Indeed, it is well-known that there is a rich and intricate connection between problems such as vicious walkers and the study of random partitions [1,2]. Interestingly, the formalism of certain low dimensional gauge theories, certainly 2d Yang-Mills theory but also 3d Chern-Simons theory for example, share many common features with such, at first, seemingly different areas [30,31]. One simple reason for this to be the case, is the well-known representation of these theories in terms of heat-kernels on group manifolds [30][31][32][33].
Integrable systems interpretation. In the two different limiting cases K = 0 (pure GWW) and λ −1 = 0 with µ = 1, the matrix model is known to be the integral representation of a τ -function of a Painlevé system, in particular Painlevé III and V, respectively [34,29,35,36,1].
In fact, the more general form of the τ -function of PV is very close in form to (2.1), differing though in having (2.1) with e N λ z , instead of (2.1) with the full GWW weight. 1 If one inspects the integral representation of the τ -function of PVI [29,35,1], it is seemingly unrelated to (2.1) and that may partially explain why our model here is essentially unstudied, whereas its two limiting cases appear in many works, often analyzed simultaneously or in a comparative fashion [34,36,29,35,1].
It is also worth mentioning another possible interpretation of (2.1), from the point of view of integrable systems [37]. The study of the so-called Schur flow [38,39], analogous to Toda flows but on the unit circle, precisely entails the generalization of a given weight function of the matrix model by multiplication by the GWW weight function.
This induces a flow that has many implications. For example, the recurrence coefficients of the polynomials, orthogonal with regards to the weight function of the matrix model, satisfy the non-linear Ablowitz-Ladik equation, with the parameter 1 g 2 YM = N λ interpreted as time. Because of this and since our analysis is in a planar limit and centered around the matrix model, the results obtained are not obviously transferable into this integrable systems and spectral theory language [40]. We will further discuss about this at the end, in the outlook Section 6.
Gauge theory interpretation. Rather, we consider (2.1) instead as a toy model for lattice two-dimensional QCD, although we will comment on some of the many other interpretations of the model. For example, (2.1) can be regarded as an effective description of two-dimensional QCD on a small spatial circle [28]. In fact, compactifying the spatial direction generates a mass gap for all but the zero-modes. Taking the small circumference limit, we are left with an effective theory with integration only over the gauge and matter zero-modes.
It was observed in [22] that the massless theory with fermions, that is = +1 and µ = 1, shows a Fisher-Hartwig (FH) singularity. The theory with bosons, on the contrary, yields a singular matrix model in the massless case. 2 Here we recognize the singularities encountered in [22] as the remnants of the IR singularities due to massless fields, and resolve them via mass deformation.
On the parameter µ. Notice that there is a slight difference in the definition of µ between the fermionic and the bosonic theory in (2.1). In the first case µ > 1 is real, while in the second case µ ∈ C with |µ| > 1. It is easy to see that the phase of µ ∈ C can be reabsorbed in a rotation of the integration contour T, and we henceforth restrict our attention to a real µ > 1 in both cases, with the understanding that for bosons µ really means |µ|.
Besides, treating µ as a real variable with this caveat in mind, the integrand in (2.1) is analytic in µ > 1, and therefore the results for any other µ ∈ C with |µ| > 1 are obtained from analytic continuation of our results. In particular, we could not attain negative values of µ moving along the real line, because we would cross the FH singularity. Nevertheless, it is possible to take a path from µ > 1 to any µ < −1 that runs in the complex plane outside the unit disk.
Remark 2.1. The independence of the partition function on arg µ is the U (1) freedom to choose the origin of T, and is the incarnation of the residual diagonal U (1) ⊂ U (N ) gauge symmetry. Our choice arg µ ∈ 2πZ fixes this residual gauge freedom.
Notation. We introduce the notation Y = 1 λ , τ = K N for, respectively, the inverse of the 't Hooft coupling and a real Veneziano parameter, whose sign carries information on the type of fields we consider.
The parameter space of the theory is Remark 2.2. For the role of the mass as a regulator (for the FH singularities on the mathematical side, for the IR singularities on the field theory side), we do not expect the continuation from M to the sheet {µ = 1} × R 2 , studied in [22,23], to be analytic. 2 According to the discussion in the previous paragraph, these well-known facts in random matrix theory can be reinterpreted as stemming from Coleman's no-go theorem [41] for massless bosons in QCD2.

2.2.
Exact finite N evaluations. We now present various analytical results for the matrix model (2.1). It is possible to evaluate the partition function of any unitary matrix model at finite N via the Heine-Szegő identity, that for (2.1) gives is the modified Bessel function. We have simplified the expression assuming = +1, although a similar determinant expression exists for = −1 as well. Formula (2.2) allows to efficiently compute Z ,K U (N ) exactly for fixed N and K. This is done in Table 1 in Appendix A.1.
The partition function for generic masses, encoded in the parameters (µ 1 , . . . , µ K ), can be also related to the expectation value of Wilson loops in arbitrary representations in the pure GWW model, thanks to the Cauchy identity (see [42] for a similar procedure). For the fermionic theory we write where the sum runs over Young diagrams R of length at most N and the first row at most K, R is the conjugate diagram, and s R is the corresponding Schur polynomial, that is, the character of the irreducible representation associated to the diagram R. In the last line we identify the correlator of two Wilson loops in the pure GWW model, with W meaning that the Wilson loop involves conjugated variables. We can in principle further expand the product of the two Schur polynomials with the Littlewood-Richardson rule, but this would entail inverting the variablesz 1 , . . . ,z N in the second Schur polynomial, as in [42]. Expression (2.3) is suggestive but not very useful as it is, because the vacuum expectation value (vev) of a Wilson loop in a generic representation is not known in closed form. Nevertheless, we can go deeper in the character expansion thanks to the formula [43] exp Using this equation and its conjugate and applying twice the Littlewood-Richardson rule we get Finally, there is an additional, simpler although more formal closed form expression for Z ,K U (N ) [37] where now the Schur functions must be interpreted as characters of U (∞), and the sum runs over Young diagrams R with at most N rows.

Phase structure
This section is dedicated to the large N analysis of the matrix model (2.1) and the determination of its phase diagram.
We stress that λ * is not a critical value of the model (2.1).
The potential is plotted in Figure 1 (τ > 0) and in Figure 2 (τ < 0). Clearly, the two figures have the same shape but upside down. Nonetheless, it is important to distinguish the role of minima and maxima to understand the phase structure.  Now that we have set the ground, we are ready to discuss the large N limit of the model (2.1).
3.1. Large N . We now take the large N 't Hooft and Veneziano limit of the partition function (3.1). This means that we consider the planar limit with both λ and τ fixed. The leading contributions to the integral at large N come from the saddle points of the effective action: Introducing the eigenvalue density with normalization chosen so that we can collect the system (3.2) of N coupled equations in a single singular integral equation at large N . The saddle point equation then reads The solution to (3.4) must satisfy the non-negativity constraint that follows from the compactness of the integration domain.
The derivation is standard [4,45], thus we omit it. We call Phase 0 the region of M for which the solution (3.6) is valid.
Critical loci. In those regions of M for which the solution (3.6) violates the constraint (3.5), we should drop the assumption suppρ = (−π, π] and look for a new solution, whose support has one or more gaps on the unit circle. The arcs on which ρ is supported are called cuts. We find a phase transition with a gap opening at θ = ±π at the critical surface Another phase transition, with a gap opening at θ = 0, takes place at the critical surface Besides, there exists a multi-critical point at the value τ = τ cr,+ (µ) at which ρ 0 (±π) = 0 = ρ 0 (0), determined as the unique point at which Y cr,a and Y cr,b meet: Examples of the limiting cases of ρ 0 (θ) are shown in Figure 3. Besides the two critical surfaces just described, looking at ρ 0 (θ) for negative Y and τ we also find values at which it attains zero value at two distinct, symmetric points in the interior of (−π, π), as in Figure 4. We expect a new phase transition into a two-cut solution.
One-cut solution: Phase Ia. We now solve Equation (3.4) dropping the assumption that ρ(θ) is supported on the whole T, and replace it by the assumption that the support is an arc Γ ⊂ T. The derivation is standard and we relegate it to Appendix A.2.
Introduce the trace of the resolvent in the large N limit, We adopt the standard notation . The support of ρ will break in two disjoint cuts beyond this critical value. Then We find (see Appendix A.2 for the details) The first term is regular and, taking the discontinuity at z = e iθ ∈ Γ, we arrive at The angle θ 0 is fixed by normalization: where y 0 := cos θ 0 . Equation (3.12) admits a unique real solution, thus the problem is completely determined.
One-cut solution: Phase Ib. The solution above has been derived assuming that Γ is an arc along T joining e −iθ 0 to e iθ 0 running counter-clockwise, thus the gap has opened around θ = π. For the gap opening at θ = 0, the procedure is identical, but now Γ is an arc fromθ 0 > 0 to 2π −θ 0 . The procedure of Appendix A.2 leads us to which is non-negative definite. There is, however, a more direct route to get the correct answer.
Looking back at the matrix model (2.1) we can chose a different parametrization 0 ≤ θ < 2π, and the solution with the gap opening at θ = 0 is recovered from the solution (3.11) in Phase Ia upon replacement Y → −Y , µ → −µ and eventually θ + π → θ.
In conclusion, we have two different phases with a one-cut solution, as expected: one for Y > Y cr,a (τ, µ), that we have called Phase Ia, and one for Y < Y cr,b (τ, µ), that we have called Phase Ib.
Two-cut solution: Phase II. We have seen that at τ = τ cr,+ = (µ 2 − 1)/2 the critical surfaces Y = Y cr,a and Y = Y cr,b meet. Thus, we expect a new phase characterized by a two-cut solution in the region with gaps around θ = 0 and θ = ±π, and eigenvalue density supported on That is, Γ is the union of two disjoint arcs, Γ u and Γ d , as in Figure 5. To determine ρ II (θ) it is simpler to adopt a different strategy, detailed in Section 3.4 below.
Two-cut solution: Phase III. The fact that the potential V eff (e iθ ) develops a double well for negative Y and τ in a given range hints at the existence of a two-cut solution in that region of M, with the eigenvalues sitting around the two minima. This observation is corroborated looking at the shape of ρ Ia (θ) and ρ Ib (θ) in the negative quadrant, where they become negative in Y cr,b < Y < Y cr,a for τ below a certain threshold. We find a transition from Phase 0 to a two-cut phase in where the critical surfaces Y = Y cr,c± are given by The two curves Y cr,± form an ellipse in each (τ, Y )-leaf of M at fixed µ, with the physical critical curve being the first branch of the ellipse encountered when decreasing Y from 0.
In this phase, that we call Phase III, the eigenvalues distribute along a contour Γ which consists of two cuts, with gaps opening around ±θ * , see Figure 6.
The eigenvalue density is Note that the argument of the outer square root is non-negative definite. The value of θ * is known explicitly, as obtained from Phase 0, and the dependence of δθ on the parameters is fixed Figure 6. The two-cut support Γ in Phase III. by normalization. Equivalently, we can fix cos (θ * + δθ) and cos (θ * − δθ) comparing the large z behaviour of ω(z) computed in this phase with its definition.
For multi-cut solutions, the dependence on the number of eigenvalues filling each cut should be taken into account when computing physical observables [46]. We analyze the role of the filling fractions in Appendix B: the upshot is that our conclusions are unaltered, both in phase II and III, although for different reasons.

Phase diagram.
Putting all the information together, the following phase diagram emerges. 0) When both Y and τ are small, Phase 0 holds, with the eigenvalues spread on the whole circle. Ia) When Y > Y cr,a the system is in a new phase, Phase Ia, with a one-cut solution gapped around θ = ±π. Ib) Likewise when Y < Y cr,b the system is in Phase Ib, with a one-cut solution gapped around θ = 0. II) At τ > µ 2 −1 2 the two critical surfaces cross each other. In the region Y cr,a < Y < Y cr,b the system is in Phase II, a two-cut solution with density of eigenvalues gapped both around θ = 0 and θ = π. III) The system develops a new two-cut phase, Phase III, in the region Y cr,b < Y < Y cr,a and also bounded by an arc of ellipse determined by Y cr,c± . The density of eigenvalues is See Figure 7 for a slice of M at fixed µ.
Taking the massless limit µ → 1 + , the critical surface Y cr,b is rotated onto the vertical axis. Using the analytic dependence on µ, we can also reach µ → −1 − by first going to the negative real axis walking through C outside of the unit disk and then taking the limit |µ| → 1 + . In that case, it is Y cr,a that is rotated onto the vertical axis. The free energy in Phase 0 is easily obtained, and corresponds to the analytic continuation of Szegő's strong limit theorem in the bulk of the 't Hooft parameter space [47]. It takes the value It is clearly separated into three contributions: pure gauge (Y 2 ), matter only (∝ τ 2 ) and the interaction. At strong coupling λ → ∞ (Y → 0) we are left with a matter contribution which counts gauge singlets: indeed, the integral over the gauge group projects onto gauge invariant states.
Massless theory. As we have stressed, a core assumption of our analysis is |µ| > 1, and the massless limit |µ| → 1 + can only be taken at the end. Due to the non-analyticity for µ ∈ T, the resulting model will differ from a model with massless matter [22,23].
A main consequence of this non-analyticity is the spontaneous chiral symmetry breaking, that we will discuss in Section 4.4. On the other hand, it is well known that the large N limit and the massless limit do not commute.
In the naïve |µ| → 1 limit, the free energy in (3.14) has a logarithmic divergence in the matter contribution. Setting instead |µ| = 1 from the beginning, and arg µ =θ, 0 <θ ≤ 2π, the partition function acquires a FH singularity and the large N limit cannot be understood by standard methods. We use known results on Toeplitz determinants to derive the free energy in Phase 0 in the massless theory [48]: That is, the contribution from matter fields has an additional factor of log N and dominates at large N . Remarkably, this matches the logarithmic divergence of the naïve massless limit of (3.14). The result is in fact much more general [48] and directly extends to the case of various Veneziano parameters τ 1 , . . . , τ n associated to different µ 1 , . . . , µ n that approach the unit circle from outside at different anglesθ 1 , . . . ,θ n .

Stereographic projection.
To better understand Phase II and the transition from a onecut to a two-cut phase, we map the model onto the real line and study the resulting Hermitian matrix model at large N . It can be interpreted as a massive deformation of the model in [49]. We conformally map the unit circle on the real line through the stereographic projection, see Figure 8. The drawback of the stereographic map is that it introduces a puncture on the circle at θ = ±π: this has no effect at finite N , but the Hermitian matrix model will fail to reproduce Phase 0 of the unitary matrix model because of this change in topology [45]. Phase 0 and its associated transitions are well understood from the unitary matrix model side, and we use the conformally mapped model as yet another way to gain further insight into the one-cut to two-cut transition. Figure 8. The stereographic projection. The red cross is the puncture on the circle, the ticker line is the cut Γ.
Our choice of coordinates is consistent with Phase Ia on the circle, but Phase Ib is easily retrieved rotating T by e iπ , so that the puncture is placed at θ = 0. The Hermitian matrix model is with the superscript p. as notation to remind that it comes from the projection of the model (2.1). We have adopted the shorthand notation Remark 3.1. Thanks to the mapping of the Vandermonde determinant form the unit circle to the real line that shifts τ → τ + 1 in the denominator of the integrand in (3.15), the stability issues pointed out in [49] do not arise here. We can thus safely allow τ < 0 without spoiling the convergence of the matrix model.
The solution is found by standard large N techniques [50]. Using a one-cut ansatz for the density The value of A is fixed by normalization: As a cross-check, turning off the mass deformation, µ → 1, sends η → ∞ and we recover the eigenvalue density found in [49]. Besides, sending A 2 → ∞ and expanding at leading order in 1 A the normalization becomes the consistency condition correctly reproducing the critical surface Y cr,a in the limit in which A is back-projected to e iπ . We stress that, requiring that ρ p. (x)dx descends from a measure on T, the non-negativity constraint ρ(x) ≥ 0 must be imposed. Looking at ρ p. I (0), we find that the critical point is fixed by the condition Phase II. From the result above as well as from the analysis of the unitary matrix model, we find a phase transition to a two-cut solution, with a gap opening at x = 0. The new phase is the conformal image of Phase II of the unitary matrix model. We look for a new eigenvalue density, supported on The parameters A and B are fixed by normalization, and by an additional self-consistency condition on ω(z), which reproduces the criticality condition for B → 0. This is not the end of the story for Phase II. Indeed, fluctuations in the number of eigenvalues in each cut may contribute at leading order in the evaluation of observables [46]. However, for the symmetric two-cut solution we find out that this is not the case, as proved in Appendix B.

Wilson loops and instantons
We continue the investigation of the features of the phase transitions and establish their order by evaluating the vacuum expectation value (vev) of the Wilson loop in the fundamental representation. Moreover, we further discuss the different physics of the various transitions by looking at the different contributions by instantons.

Wilson loops.
Wilson loops are order operators in gauge theories that, for simple connected gauge group, describe the holonomy of the gauge connection around a closed path. For our oneplaquette model, we consider the Wilson loop in the fundamental representation wrapping the plaquette, and compute its vev. It is given by cos θ j with the average taken in the unitary ensemble (2.1). We use the eigenvalue density at large N found in each phase to evaluate the Wilson loop.
Wilson loops: Generalities. From the matrix model (2.1) we immediately get the relation Therefore, all the information about the order of the transition can be extracted from the Wilson loop vev. This is precisely what we expect from an order parameter, and follows from the Wilson loop belonging to the class of order operators of QCD 2 .
Being ρ(θ) continuous on the whole M, the Wilson loop vevs are continuous as well, implying that every phase transition we find must be at least second order.
Wilson loops: Evaluation. We focus now on the Wilson loop vev at large N . In the ungapped phase we find This reproduces the GWW result as τ → 0, but also as µ → ∞, as expected when the matter becomes non-dynamical. For a Wilson loop winding k > 1 times around the plaquette, either in clockwise or anti-clockwise direction, we get

In Phase Ia the Wilson loop vev is
where we have used the change of variables y = cos θ, with y 0 = cos θ 0 . The value of y 0 as a function of the gauge theory parameters is known from (3.12).
The Wilson loop vev in Phase Ib is obtained likewise, To study the derivative of W Ia and establish the order of the phase transition, it suffices to notice that d dY This implies which matches the derivative of W 0 . The computations are identical for the transition between Phase 0 and Phase Ib. Taking a further derivative, d 2 dY 2 W vanishes identically in Phase 0, but does not vanish at the critical loci when computed in Phases Ia and Ib.
We conclude that the Wilson loop vev is an order parameter of class C 1 at the critical surfaces Y = Y cr,a (τ, µ) and Y = Y cr,b (τ, µ), thus the system shows a pair of third order phase transitions. In particular, both the GWW transition [4,5] and the transition in [27] are special points on the critical locus of the present model.
Crossing from a one-cut to a two-cut phase, the first derivative of the Wilson loop is not protected. Indeed, in Phase II the derivative of the Wilson loop vev has the schematic form with the first term coming from the derivative of the explicit dependence on Y , and the other two from the dependence on Y through y 0 andỹ 0 . The integrand evaluated at the endpoint vanishes, hence those contributions do not appear. In (4.4), f (y, y 0 ,ỹ 0 ) is known explicitly from Section 3.1, but the difficulty comes from the only implicit knowledge of the dependence of y,ỹ 0 on Y . Passing from Phase II to Phase Ia, the first term in (4.4) matches continuously with the corresponding expression in Phase Ia, asỹ 0 → 1. The integral in the second summand in (4.4) also agrees with the corresponding contribution in Phase Ia atỹ 0 → 1. Both facts follow from Moreover, the symmetries of the integrand allow to combine the third term in (4.4) with the second term, in a simpler expression. Moreover, the symmetric form of the equations fixing y 0 ,ỹ 0 can be used to show that ∂ỹ 0 ∂Y = ∂y 0 ∂Y ỹ 0 ↔y 0 .
By this we mean that the expressions on the two sides agree upon exchanging allỹ 0 with y 0 . Due to the complicated dependence on the parameters, the derivatives of the boundaries y 0 ,ỹ 0 are not continuous at the transition point. The differentiability of W above followed by the vanishing of the term multiplying such derivatives. This does not happen for the transition from a two-cut to a one-cut phase. Therefore, the sum of the second and third terms in (4.4) gives an obstruction to the differentiability of W , so we expect a second order transition. In a sense, the obstruction arises from taking a limit that breaks explicitly the y 0 ↔ỹ 0 symmetry of Phase II.
The proof is very similar for the transition from Phase II to Phase Ib or from Phase III to either Phase Ia or Ib.
The argument fails at the critical surface Y cr,c and at the multi-critical point at which Y cr,a = Y cr,b . Indeed, when passing directly from Phase 0 to a two-cut phase, the simplifications that arise from closing both gaps simultaneously imply that the Wilson loop vev is C 1 . This is consistent with the observation of the previous paragraph, as these transitions preserve the Z 2 -symmetry of the two-cut phase.

4.2.
Phase structure and remarks. Summing up the results extracted from the analysis of Wilson loop vev, we find that • the transition from Phase 0 to any other phase is third order, but • the transition from a one-cut to a two-cut phase is second order.
In the rest of this subsection we gather comments on various aspects of the phase structure we uncovered, insisting on the role of the second order phase transitions.
Remark 4.1. As obtained in the previous subsection, the second order discontinuities are finite jumps, not divergences. The correlation lengths remain finite at each transition. These finite discontinuities vanish in the limit |µ| → 1.
Metastability. While the third order transitions we find are a continuation of the GWW transition in M, it is worth to further comment on the second order transitions we obtain. The phase transition to a two-cut solution happens slightly beyond the values of τ where the potential develops a double-well structure. The proposal in [51] states that a second order transition can be associated with tunneling from a metastable vacuum to a stable one. Our analysis confirms that picture in the one-plaquette model we consider. In Section 4.5 we study instanton effects, expanding this discussion leading to a further refined distinction between second and third order phase transitions in this model, from the instantonic point of view.
Critical behaviour. It is worthwhile to notice that the phase diagram in Figure 7 resembles that in [52], where a unitary matrix model with potential Y 1 cos(θ) + Y 2 cos(2θ) was analyzed. The critical behaviour close to a transition to a two-cut phase in our model differs from that found in similar models in the literature, for matrix models with potentials of the form K n=1 Y n cos(nθ). This is so because the potential in (2.1) includes both a polynomial and a logarithmic part, requiring different scaling approaching the critical regime from the two-cut phase. This distinction, however, fades away approaching the multicritical point.
Double-scaling limit. The statements above can be refined exploiting the double-scaling limit.
In particular, we can zoom in the critical regime, tuning Y towards a transition to the ungapped phase. In the double-scaling limit, the dynamics is governed by Painlevé II equation. The proof follows from [53,27] with minimal variations. An alternative proof can be given using orthogonal polynomials [54]. We have checked explicitly that, in the double-scaling limit, the problem reduces to the analogous one for the pure GWW model.
For the transition from a one-cut to a two-cut phase, however, there is no double scaling that gives Painlevé II.
Other gauge groups. Throughout the work, we focus on gauge theories with gauge group U (N ) or SU (N ). Nevertheless, by direct computation or by universality arguments, it can be shown that the phase diagram of Figure 7 and the associated phase transitions carry over to theories with gauge group G, that is to say, matrix models like (2.1) with the integration over U (N ) replaced by integration over G, for any Here, O − (n) ⊂ O(n) consists of n × n orthogonal matrices with determinant −1, and Sp(N ) is the compact symplectic group. 4.3. Continuum limit and β-function. The β-function of the theory, as a function of the 't Hooft coupling λ, can be computed using the chain rule through [4] (4.5) This quantity can be used to test whether our model reproduces the expected features of QCD 2 in the continuum limit. The fixed points of the RG flow, that capture the continuum physics, are given by β(λ) = 0, which, from (4.5), can only happen at W = 0 or at W = 1.
Direct computations in Phases 0, Ia, Ib, show that only the solution to W = 0 is consistent, while the solution to W = 1 always falls out of the phase in which it has been computed, and thus should be discarded. It is a nice consistency check that the solution to be discarded is precisely the one that would violate Elitzur's theorem [55], and the one to be retained is in agreement with the confining nature of QCD 2 [4].
The continuum limit of a lattice theory consists in sending the lattice spacing to zero while approaching a critical curve [56]. In particular, this requires |µ| → 1.
Taking the continuum limit from Phase 0, approaching either Y cr,a or Y cr,b , we find that the unique consistent solution is Y = 0 (i.e. λ = ∞). This is physically meaningful for a toy model of QCD 2 : the theory flows to a strongly interacting theory in the deep infrared.
Taking the continuum limit from Phase Ia close to the transition to Phase II, we find a trivial solution with λ = 0 = τ , describing a theory of free gauge bosons without matter. The continuum limit approaching the critical surface between Phase Ib and Phase II, instead, yields a non-trivial fixed point at Remark 4.2. The existence of a continuum theory is not established by our analysis, because correlation lengths remain finite. While this has no effect in our model, which consists of a single plaquette, it may (and most likely shall) wash away the fixed point in the continuum limit of a more realistic lattice model.

Chiral symmetry breaking. Let us focus now on the model with fermionic matter. The fermion two-point function is by definition
Due to our degenerate choice of masses, we can only compute the average over flavours of such quantity: 1 In Phase 0 we find 1 τ This quantity diverges as µ → 1, therefore we expect the chiral symmetry to be spontaneously broken in the continuum, consistently with the analysis of the β-function in Phase 0. In Phases Ia and Ib, we can move along Y = 0 and study the behaviour of 1 K fψ f ψ f on that subspace of M. The result is read off directly from [45]: in Phase Ia, which is non-vanishing and continuous at the transition point.
again non-vanishing and continuous at the transition point, and goes to 1 in the µ → 1 + limit. This latter result, in turn, hints at a transition to a free theory: the free energy of a theory of K free flavours of mass m goes as F free ∝ Km, whence 1 K fψ f ψ f | free = 1. Note that this computation has been done at infinite gauge 't Hooft coupling, which has the physical meaning of governing the theory in the deep infrared.
To sum up, we have observed that the phase transition from Phase Ib to Phase 0 is accompanied with spontaneous chiral symmetry breaking.

4.5.
Instantons. We discuss non-perturbative effects in the unitary matrix model, coming from unstable saddle point configurations [57].
An instanton configuration is characterized by a collection of integers {N 0 , N 1 , . . . } with k N k = N . For example, the d-instanton configuration is associated with the symmetry breaking pattern with the eigenvalues z j ∈ T of U ∈ U (N ) grouped in d different sets, sitting at d different extrema of the potential. For the one-instanton configuration, with Z the partition function of a U (N − ) × U ( ) model, and we have turned on a chemical potential ν > 0 for the instanton number.
The -sector leads to non-perturbative corrections to the free energy of the matrix model, of the form where λ generically denotes the couplings of the theory, and the functions {f } admit themselves a 1 N expansion.
Instanton effects and third order transitions. Let us consider our model (2.1) at large N and focus on the one-cut phase, in which the interpretation of instanton effects is more transparent. We discuss them in Phase Ia, being the corresponding analysis in Phase Ib completely analogous.
Most of the details are just an extension of the thorough analysis in [57]. The contribution of an instanton excitation, obtained moving one eigenvalue from the minimum of V eff to the maximum at θ = ±π is found to be where cosh −1 and tanh −1 are the inverse of the hyperbolic functions, and x 0 = sin θ 0 2 . One of the results in [57] (already conjectured in [58]) is that the GWW transition is triggered by instantons. We see that the result carries over to the present model, as Analogous conclusions are found if we go to Phase III, in which the effective potential has developed a double well, and consider the instanton configuration with a few eigenvalues taken to the local maximum at θ * . Approximating close to the transition to Phase 0, we find At the critical surface, δθ → 0 and we find again that the third order transition is triggered by instantons. Instanton effects and second order transitions. We now turn our attention to the analysis of instanton effects in the two-cut phase, starting from Phase III. We consider a single eigenvalue placed on the maximum of the potential, as sketched in Figure 9. We find that the instanton action is the sum of two pieces, where y * = cos θ * and y L,R = cos(θ * ± δθ). The two are associated with the eigenvalue escaping from the left and right cut, respectively. There exists a third relevant quantity, namely the tunneling from one cut to the other, All the three effects are non-perturbatively suppressed by a factor e −N S inst , with S inst the corresponding action. The three contributions are still suppressed at the critical loci, although the tunneling term will coalesce with one of the other two.
The situation is slightly different in Phase II, where the two wells have equal depth, see Figure  10. In this case, S tunnel is simply twice S inst , and both go to zero as a gap closes. The phase transition takes place when the tunneling between the two cuts ceases to be suppressed in one direction (e.g. passing through θ = 0 in Figure 10) but remains non-perturbative in the other direction (e.g. passing through θ = π in Figure 10).
The picture we infer is that the third order phase transitions are associated with releasing non-perturbative instabilities, while the second order transitions correspond to release only those instabilities in one direction. This is also in agreement with the proposal in [51,18] relating second order phase transitions in GWW-type models to partial deconfinement.

Meromorphic deformation of unitary matrix models
We now depart from the model (2.1) with the aim of setting the stage for the study of meromorphic deformations of unitary matrix models, in which the integration contour is deformed in C * and not bound to be the unit circle. This consists of an adaptation of the theory of holomorphic matrix models [19] to unitary matrix models and, as we shall show, is instrumental in The function e − N λ V (z) is singular at z ∈ {0, ∞} ⊂ P 1 and possibly has other zeros and poles in C * ∼ = P 1 \ {0, ∞}. The Vandermonde determinant appearing in a unitary matrix model is conveniently rewritten in meromorphic form To deform a unitary matrix model, the integration cycle T N is replaced by any half-dimensional cycle C N in (C * ) N .
Definition. Let N ∈ N, λ ∈ C * and V (z) as in (5.1). A meromorphic matrix model Z is the integral with integration contour Condition (5.3) means that each C is homotopic to the unit circle in the holed plane C * . Dropping it, we may stretch C along any direction along which, asymptotically, 1 λ V (z) > 0. Eventually the one-cycle C pinches at z = ∞ ∈ P 1 . To get an honest deformation of a unitary matrix model we do not allow this situation, otherwise we would fall back in the holomorphic deformation of a Hermitian matrix model.
The couplings {t n } in (5.1) are usually subject to reality conditions, such as t −n =t n , ∀n ≥ 1, and possibly other relations. We write the constraints collectively as Φ ({t n }) = 0. Besides, a rescaling of all {t n } together can be reabsorbed in a redefinition of λ, hence the couplings {t n } are homogeneous coordinates on a projective space. The matrix model (5.2) sets a natural stage to complexify the couplings. We denote by T the physical space of couplings, namely the collection of independent {t n } after imposing the constraints and modulo scaling. More formally, with the C * -action being multiplication of all couplings by a non-vanishing constant. Whenever the constraints Φ can be rewritten in homogeneous form, T is a projective variety.
5.1. Large N limit. We are interested in the large N limit of (5.2). Define the effective potential At large N , the eigenvalues will be gathered around the saddle points of W (z) in C * , W (z sp, ) = 0, = 1, 2, . . . , g + 1.
Here we are assuming there is a finite number g + 1 of saddle points z sp . The integration contour C N in ( In the large N limit, the eigenvalue density solves the saddle point equation where means holomorphic derivative ∂ ∂z . Here, dw is a holomorphic differential on Γ and, given any parametrization w : s → w(s) ∈ Γ, dw =ẇ(s)ds is understood, with ds the line element anḋ w the derivative of the map w : s → w(s).
Recall the definition of the trace of the resolvent ω(z) at large N : Defining (5.7) y(z) = ω(z) − 1 2 W (z), (5.6) becomes This equation goes under the name of spectral curve. The steps from (5.5) to (5.8) are standard and have been applied to holomorphic matrix models since their early days [59] to establish a bridge between matrix models and geometric problems. The novel aspect of (5.8) compared to the literature is hidden in the form of W and f , which in the present case are not ordinary polynomials but Laurent polynomials, or meromorphic functions.
We assume for now that V (z) is a Laurent polynomial on P 1 , with singularities at {0, ∞}. The extension to a meromorphic weight function on C * is worked out below. Write where we assume d − ≥ 2 (otherwise we get back the known setting). We need to introduce some notation. Define which agrees with the counting of saddle points above, and also t 0 = −λ for later convenience. Besides, denote ρ k the moments of the eigenvalue density, After some rewriting we get where sgn(0) = +1 by convention. The spectral curve takes the schematic form (5.13) y 2 = P (z) 4λ 2 z 2d − where P (z) is a polynomial in z of degree deg(P ) = 2g + 2, with coefficients read off from (5.11)-(5.12) and that depend on the parameters {t n }, on λ and on the moments {ρ k }.
A major difference with respect to the standard unitary matrix models is that {ρ k } are free complex parameters of the theory: they can be tuned deforming Γ, as discussed in Remark 5.1. Recalling that both z, λ ∈ C * , it is possible to recast (5.13) in a more standard formŷ = P (z), describing an hyperelliptic complex curve of genus g [59].
Γ is the union of g + 1 branch cuts stretched between pairs of roots of P (z). The roots of P (z) move inside C * as the parameters are varied. 3 The coalescence of two roots produces a singularity of the curve (5.13) and corresponds, on the matrix model side, to a phase transitions from a (g + 1)-cut to a g-cut phase, with either • two cuts joining, or • one cut collapsing.
The hyperelliptic curve (5.13) is fibered over the moduli space M of the model (5.2), defined as M = C * ×T × C g , with C * parametrized by λ, T defined in (5.4), and the last factor parametrized by the moments {ρ k }. Note that one of the moments is fixed comparing (5.6) with the definition of ω(z) at |z| → ∞.

Definition.
A critical locus C is an irreducible component of the locus in M at which two roots of P (z) coalesce.
The critical loci C ⊂ M necessarily have positive complex codimension, and the hyperelliptic fibration is singular along them. Singularities in higher codimension, placed at the (self-)intersection of critical loci, correspond to multicritical points of the matrix model.
The theory of Abelian differentials provides a suitable framework to analyze the genus g hyperelliptic curve (5.13) [60]. At this stage, the analysis of the spectral curve works exactly as in the holomorphic deformation of Hermitian matrix models, thus we omit the details and refer to [60,61].
Remark 5.2. We are now in the position to elaborate more on Remark 5.1, from a point of view advocated in [62]. Let A , B be a basis of one-cycles in the hyperelliptic complex curve. The A-cycles are chosen to go around the cuts Γ . Therefore The first equality follows from the definition (5.7) noting that y(z) and ω(z) only differ by a regular term. Introducing chemical potentials for the filling fractions ξ and extremizing the action with respect to these quantities gives their saddle point value as a function on M . More precisely, one gets a set of equations analytic in the ratios s := ξ λ [62]. At this point, it is possible to invert the relations and express the moments {ρ k } in terms of the complex variables s , keeping the latter as free parameters.
Note that only g out of the g + 1 of both quantities are free.
The study of the phases of the model (5.2) leads to a stratification of the parameter space M . We postpone the analysis to Section 5.4, discussing explicit models first.
Genus 0. The unique way to obtain a genus 0 spectral curve is from the holomorphic deformation of the CUE. In that case, the model has no couplings and, as opposed to g ≥ 1, the additional condition derived from the definition of ω(z) is automatically fulfilled, leaving ρ −1 as unique, unconstrained parameter. Then, (5.13) describes a P 1 fibered over C. If we try to get a less trivial model by considering the insertion of (det U ) τ N , the consistency condition, which in g ≥ 1 fixes one of the {ρ k }, imposes τ = 0.

Holomorphic GWW.
We now put the machinery at work and revisit the phase structure of the holomorphic GWW model. The phase diagram of this model has been obtained in [63] for λ ∈ R, while the behaviour at complex coupling has been partially analyzed in [64], although without fully exploiting the holomorphic deformation.
It is an elliptic curve. Following the strategy outlined above, we think of (5.14) as an elliptic fibration over C * × C, with coordinates on the base λ and ρ −1 , and identify the phase transitions with singularities of the fibration.
The discriminant of (5.14) is from which the critical loci are In Kodaira's classification [65], C 1 and C 1 are singularities of type I 1 and C 2 is of type I 2 . The GWW critical points (λ, ρ −1 ) = (±2, 1) are singled out as the codimension-two singularities at which C 2 intersects one of the other two critical curves. Besides, we recognize the elliptic curve (5.14) as the Seiberg-Witten curve of N = 2 supersymmetric four-dimensional SU (2) gauge theory with two flavours [66]. The original GWW transition is thus, from the perspective of the holomorphic deformation, one of the possible ways to approach the codimension-two singularity from a generic direction. The singularity at the multicritical points (λ, ρ −1 ) = (±2, 1) is of Kodaira type III. The corresponding symmetry is A 1 , which is precisely the symmetry of Painlevé II, that is known to control the GWW phase transition [54,53]. If, instead, we approach the multicritical point not from a generic direction but moving along a critical locus, the singularity type is enhanced to I * 0 . It is possible to allow t −1 = t 1 . This corresponds to introduce a θ-term in the GWW lattice action, θ . The procedure goes through with only minor modifications, the unique difference being that the singularities C 1 , C 1 are placed at ρ −1 = − 1 t −1 ± 4 λ √ t −1 , so in particular the Z 2 -symmetry is preserved.

Meromorphic deformations.
The formulation can be extended to include weight functions with zeros and poles in C * . For concreteness, we consider the illustrative example of our original model (2.1) at λ −1 = 0. In this case In the second line, we have definedρ with, in particular,ρ −1 (0) = ρ −1 . Comparing the definition of y(z) with the spectral curve at large |z|, we find a pair of consistency conditions, fixingρ −1 as a function of the other parameters, Note that the two conditions fixρ −1 (µ ±1 ) independently, and the solutions are consistently mapped into each other under µ ↔ µ −1 . We get (5.16) with P (z) a polynomial of degree 4. The spectral curve thus describes again en elliptic fibration over the moduli space C * × µ ∈ P 1 : 1 < |µ| < ∞ × C, parametrized by (τ, µ, ρ −1 ). The discriminant takes the form ∆ = µ 5 (µ + 1) 2 (µ − 1) 2 (τ + 1) 4 ∆, the last term being a cumbersome polynomial of degree 6 in ρ −1 , degree 8 in τ and degree 10 in µ. The critical points of the undeformed model become higher-codimensional singularities, at which two roots of P (z) collide. The collection of all critical loci in this model is with the subscript indicating the order of vanishing of ∆ along the component C . Taking what we have called the continuum limit in Section 4.3, that is, sending τ → τ cr (µ) and then µ → ±1, with ρ −1 = − τ µ set to its undeformed value, yields a non-minimal singularity ∆ ∝ (µ ± 1) 12 . where the superscript means the regular part, and C IJ = C I ∩ C J , and so on. The inclusion relations C reg IJ = C I ∩ C J ⊂ C I are obvious. The partial order can be represented with the aid of a Hasse diagram: In general, this does not define a full-fledged stratification of M because multiple final points may exist. Nonetheless, whenever the potential (5.1) has a Z 2 -symmetry, the Hasse diagram inherits it. This Z 2 -symmetry acts as an automorphism of the Hasse diagram, which is mapped into itself under reflection along the vertical axis. By construction, the Hasse diagram resulting from folding the initial diagram via this Z 2 -symmetry determines a stratification of M /Z 2 . We draw the Hasse diagram of the holomorphic GWW model of Section 5.2: The diagram of the meromorphic model of Section 5.3 is schematically The vertical red, dashed line is there to emphasize the Z 2 reflection symmetry. Folding the diagram along that line yields the stratification of M /Z 2 .
Remark 5.3. The results of Section 5.2 with t −1 = t 1 show that, even for models in which a Z 2 reflection symmetry is not manifest from the potential but emerges at large N , the Z 2 -folding yields a stratified moduli space.
Symplectic singularities. Recall that H 1 (Σ g , C), the first homology group of a hyperelliptic curve Σ g of genus g, is a symplectic space. The A-and B-cycles that we have implicitly used in the study of the spectral curve (5.8) can be chosen to be Darboux coordinates in H 1 (Σ g , C).
Moving along M reg corresponds to vary the symplectic structure without changing the topology of Σ g . At the critical loci C I , however, either • a B-cycle collapses, or • an A-cycle collapses.
Both situations correspond to a singularity of the symplectic form. Therefore, the analysis of the phase structure of the meromorphic matrix models can be rephrased in terms of symplectic singularities in the sense of Kaledin [67]. The appearance of symplectic singularities is not entirely unexpected. The consideration of holomorphic matrix models in their large N limit lead to the Seiberg-Witten curves [66] of certain N = 2 gauge theories [59]. The so-called Coulomb branch of these theories is a symplectic singularity and is stratified [68]. In fact, the use of Hasse diagrams in the present work was inspired by [69,70].
As a final remark, we emphasize that the structure uncovered in this section is not specific of the meromorphic matrix models. The parameter spaces of unitary or Hermitian matrix models inherit it, as they can be realized as slices inside the parameter space of our meromorphic models.
As an example, the phase diagram of the GWW model is It is found by fixing ρ −1 = 1, taking the slice λ ∈ R \ {0} and identifying the intersection of such subspace with the strata in (5.17). 4

Outlook
We conclude by commenting, in a qualitatively and non-exhaustive fashion, on avenues that we have considered at some point but not pursued here.
It would be interesting to know if the results obtained have a meaning from the point of view of integrable systems such as the Schur flow. While both the matrix models considered and the study of such flows have in common an associated system of orthogonal polynomials, the way this association actually works is quite different. For example, the recurrence coefficients of the polynomials, central in the integrable systems description, are not directly relevant in the type of matrix model analysis presented.
On the other hand, for what concerns matrix models on the real line, the spectral properties of Jacobi matrices can be more directly related to matrix models, since it is known that, under rather general conditions, a suitably normalized counting measure of the zeroes of the orthogonal polynomials converges weakly, in the large N limit, to the density of states of the matrix model [71]. With this in mind, a question would be whether our results have any implication in the study of the spectral properties of CMV matrices [40], for example. The large N planar limits taken make this possibility not obvious.
Another reason to further study any eventual implications of the planar limit and the ensuing phase structure, from the point of view of integrable systems, would be the connection obtained, presented in Section 5.4, between the phase diagram of meromorphic matrix models and the symplectic foliation of singular varieties [67]. Again, the 't Hooft scaling of the couplings involved at large N obscures the relation between the symplectic structures we naturally find and those in the integrability literature [72].
Also, from a mathematical point of view, it would be interesting if proofs of the order of the phase transition, in particular for the second order phase transitions, can be obtained in alternative or more rigorous ways.
Regarding more physical considerations, when discussing the model with fermionic matter and its interpretation in terms of chiral symmetry breaking, it is worth mentioning that recently [73], the chiral symmetry breaking phase transition in four-dimensional QCD has been studied from the point of view of thermodynamic geometry [74,75]. The argument is based on the observation that the grand canonical partition function determines a metric g thermo on a two-dimensional parameter space with coordinates (F, v) [74], where F is the free energy and v the grand canonical chemical potential. Then, a second order phase transition is triggered by the instability at det g thermo = 0. Any eventual use of this observation or other ideas from information geometry to further understanding phases in matrix models would be of interest.
A.2. Large N limit: Gapped solutions. In this appendix we sketch the computation of ω(z), defined in (3.10), which allows to extract the eigenvalue density in the phases with one or more gaps. The procedure is standard and we follow closely [76,47], glossing over many details. We work in Phase Ia, since all other phases are analyzed in similar fashion.
Introduce the function (z) of complex variable z ∈ C such that (e iθ ) = ρ(θ) for e iθ ∈ Γ. The saddle point equation (3.4) is rewritten as Equation (A.1) is valid for z ∈ Γ, and is complemented by the normalization condition Recall that we have started with a Z 2 -symmetric system, invariant under z → z −1 for z ∈ T. We will thus find an eigenvalue density with symmetric support, and in particular ∂Γ = e −iθ 0 , e iθ 0 in a one-cut phase. Then, depending on whether the gap opens at θ = π or θ = 0, Γ will be the arc on the unit circle connecting −θ 0 to θ 0 or θ 0 to −θ 0 , respectively, with orientation always taken counter-clockwise. Recall from the definition (3.10) that In turn, from the definition of Cauchy principal value and (A.1) we immediately get The normalization (A.2) and the definition (3.10) imply that ω(z) → 1 as |z| → ∞. We have then reduced the problem of finding the eigenvalue density to the problem of determining the discontinuity of ω(z) along Γ, from the knowledge of its regular part and the boundary condition ω(z → ∞) = 1. It is standard procedure to reduce the problem For a multi-cut phase, with Γ ∼ = {θ 0,− ≤ θ ≤ θ 0,+ } ∪ {θ 1,− ≤ θ ≤ θ 1,+ } ∪ · · · ∪ {θ k,− ≤ θ ≤ θ k,+ } the procedure is the same, but with Ω(z) defined via ω(z) = k j=0 e iθ j,+ − z e iθ j,− − z Ω(z).
Let us now introduce a closed contour C Γ which is a Jordan curve enclosing Γ but not z, and oriented counter-clockwise. See Figure 11 for the contour C Γ in Phase Ia.
From the definitions (3.10) and (A.4) it follows that Ω(z) falls off (at least) as 1/z at infinity. Then, for z lying in the exterior of C Γ , Cauchy's theorem together with (A.5) implies On the other hand, because W (w) is meromorphic we can deform the contour C Γ into an infinitely large circle, picking the poles of the integrand. We find where the first term is the residue at w = z, the second term is the remaining contour integral along a circle at infinity, which in our case simply contributes Y , and the last term includes the residues at the poles z p of W (w).
In Phase Ia, explicit computation of each term leads to The solution in the other phases is found likewise.

Appendix B. Filling fraction fluctuations
This appendix contains the analysis of the effect of taking into account fluctuations of the filling fractions around the equilibrium configuration.
For a generic matrix model in a two-cut phase, the dependence of the filling fractions on the parameters of the theory should be taken into account when computing physical observables [46]. 5 Below we briefly review how this effect comes about, and study it for the model at hand. We start with Phase III, and look at Phase II projected onto the real line, as in Section 3.4. B.1. Phase III. For the two-cut solution in Phase III, let N L be the number of eigenvalues in the left arc around θ = π, 0 ≤ N L ≤ N , and N R = N − N L the number of eigenvalues in the right arc around θ = 0. Let also ξ = N L N and 1 − ξ = N R N denote the corresponding filling fractions. The values of y L = cos (θ * + δθ) and y R = cos (θ * − δθ) can be fixed, as functions of ξ and of the other parameters, through the equations that come from the definition of ξ after changing variables y = cos θ. Then, the value of ξ is fixed by the equilibrium condition (B.1) dS eff dξ ξ=ξsp = 0.
For example, approximating close to the critical surface diving Phase III from Phase Ib, we find where we have also substituted Y = Y cr,b . It has been shown in [46] that the quantum fluctuations around the saddle point ξ sp contribute to the free energy a term of the form − 1 N 2 log ϑ(N ξ sp ), where ϑ(z) is the Jacobi theta function. See Appendix B.2 below for more details and a very short review of the derivation. This is a sub-leading contribution to the free energy but, due to the dependence on N ξ sp , each derivative generates a factor of N . Therefore, the ξ sp -dependent part becomes of the same order as the leading order term when differentiating the Wilson loop vev, and must be taken into account. The relevant part of the derivative is B.2. Phase II. We now discuss the same effect in Phase II. It is more convenient and akin to the work [46] to do this in the alternative, Hermitian matrix model presentation of Section 3.4. The argument can be succinctly summarized as follows.
Consider a two-cut solution with support suppρ p. II = Γ L ∪ Γ R , and denote by ξ = N L N and 1 − ξ = N R N the corresponding filling fractions, as above. The saddle point value ξ sp of ξ is fixed by (B.1). Then, in the large N approximation, the partition function takes the form [46]  where F pert is the perturbative free energy to all orders in the 1 N 2 expansion. This yields [46,77] where F is the leading order or planar free energy, F nlo the next-to-leading order correction, and (after an implicit resummation) we have recognized the Jacobi theta function ϑ(N ξ sp ). The modular parameter of the theta function is i2π/(∂ 2 ξ S eff (ξ sp )), and the dependence on it is kept implicit in the notation.
For the case at hand, however, the effective action is an even function, the two wells have identical depth, and all the physical observables we consider preserve this property. We thus have ξ sp = 1 2 , independent of the parameters of the theory, and the effect we have just described will remain sub-leading [77]. This would not be the case for other type of physical observables that are not protected by the parity symmetry. See [46] for discussion and examples.