Higher Rank Wilson Loops in N = 2* Super-Yang-Mills Theory

The N=2* Super-Yang-Mills theory (SYM*) undergoes an infinite sequence of large-N quantum phase transitions. We compute expectation values of Wilson loops in k-symmetric and antisymmetric representations of the SU(N) gauge group in this theory and show that the same phenomenon that causes the phase transitions at finite coupling leads to a non-analytic dependence of Wilson loops on k/N when the coupling is strictly infinite, thus making the higher-representation Wilson loops ideal holographic probes of the non-trivial phase structure of SYM*.


Introduction
The path integral of any N = 2 gauge theory on S 4 can be calculated [1] by means of the supersymmetric localization [2]. Localization reduces the path integral to a matrix model giving us direct access to strong coupling, and in particular to the regime of interest for holography, when N is infinite and the 't Hooft coupling λ = g 2 YM N is large. The SYM* theory is particularly well-suited for this purpose, since its holographic dual is explicitly known [3], while the N = 2 * localization matrix model can be solved at strong coupling [4,5,6] by more or less standard methods of random matrix theory [7]. Various predictions of holography can then be confronted with ab initio evaluation of the field-theory path integral [4,8]. Perfect agreement found so far has left, nevertheless, one feature unexplained from the holographic perspective. Localization predicts an infinite sequence of quantum phase transitions that occur at certain critical values of the 't Hooft coupling in the large-N limit [9,10]. Large-N phase transitions are of course very common in matrix models [11], but here the real critical behavior occurs such that correlation length goes to infinity. The ensuing phase transitions also have finite-N counterparts [12]. Similar transitions occur in other localization path integrals [10,13], and it would be very interesting to understand their holographic origin.
Unfortunately, tuning the 't Hooft coupling to some given critical value is very difficult in holography. Well-established holographic methods work at strictly infinite coupling, leaving us little hope to observe the transitions as discontinuities of correlation functions in λ. However, in the SYM* the n-th critical coupling grows with n as λ (n) c ≃ π 2 n 2 [10], and most of the phase transitions do happen when the coupling is very large. Certain observables may then have criticality-induced non-analyticities even when the coupling is strictly infinite. Identifying such observables is the goal of this paper.
The phase transitions are driven by nearly massless states that give rise to singularities in the eigenvalue density of the localization matrix model. The eigenvalue density can be directly probed by D-branes in the dual supergravity background [14]. The classical D-branes see the coarse-grained density, in which the resonances are smoothened out [9], but D-brane quantum fluctuations may be sensitive to the non-trivial phase structure of SYM*.
Another type of observables calculable by localization are Wilson loops. They might appear useless for detecting phase transitions at first sight. There are two scales in the problem, M and 1 R, where M is the mass scale of the SYM* theory and R is the radius of compactification on S 4 . For the most part of this paper we regard the radius of S 4 as an IR regulator and thus assume that M R ≫ 1. At strong coupling a new scale µ ∼ √ λM ≫ M emerges. Wilson loops in the fundamental representation probe the smallest scales of order 1 R, the D-branes probe the largest scale of order µ, while resonances that cause critical behavior occur at the intermediate scale of order M . In this paper we consider Wilson loops in higher representations of the gauge group, and show that by varying the size of representation it is possible to scan the whole spectrum of scales, including the resonance region.
Wilson loops of rank k ∼ N are dual to D-brane configurations carrying the fundamental string charge [15,16,17,18,19,20,21], and have been extensively studied in the N = 4 super-Yang-Mills theory [16,22,23,20,24], for which the matrix model is just Gaussian [25,26]. As far as N = 2 theories are concerned, higher-rank Wilson loops have been calculated for the N = 2 superconformal QCD [27], whose solution at strong coupling is known [28] but is very different from that of SYM*.

Wilson loops
The SYM* theory has the same field content as N = 4 SYM -gauge fields, six adjoint scalars and four Majorana fermions, but its Lagrangian includes explicit mass terms for half of the fields (four scalars and their N = 2 superpartners). We denote this common mass by M .
The Wilson loop in representation R is defined as where Φ is a scalar from the vector multiplet, coupling to which makes the Wilson loop locally supersymmetric and well-defined in the UV. Localization computes the Wilson loop for the equatorial contour of S 4 [1]. The partition function on the four-sphere of radius R localizes to an eigenvalue integral: The integration variables are the Coulomb moduli that parameterize spacetime independent zero modes of the scalar that enters the loop operator: 3) The non-zero modes of Φ, as well as all other degrees of freedom have been integrated out, leaving behind a one-loop contribution [1]: , (2.4) The Wilson loop expectation value is obtained by simply replacing quantum fields in (2.1) by their values on the localization locus and subsequently averaging over the Coulomb moduli: We denote by L the length of the contour C. The average is now defined by the eigenvalue partition function (2.2). When C is the big circle on S 4 , this result is exact (up to instanton corrections neglected here because of the large-N suppression), in which case Localization, strictly speaking, is not applicable to other contours, but the leading exponential behavior for asymptotically large Wilson loops should be universal and largely insensitive to the contour's shape. We thus expect that our results apply to any sufficiently big contour, including contours on R 4 , as long as M L ≫ 1. This assertion is based on universality and has less rigorous grounds compared to localization, but has been checked holographically for Wilson loops in the fundamental representation [4]. We concentrate on rank-k symmetric or antisymmetric representations: and will eventually consider the scaling limit N → ∞, k → ∞, in which the ratio is kept fixed. Such Wilson loops are dual to macroscopic D3 or D5 branes, depending on whether representation is symmetric or antisymmetric [15,16,17,19,20], that carry k units of electric flux on their world-volume. On the field-theory side, the k-symmetric/antisymmetric Wilson loops are described by statistical mechanics of free bosons/fermions, for which the matrix eigenvalues play the rôle of energy levels [23]. Let us review how this statistical interpretation arises in the matrix model. The k-symmetric or antisymmetric characters are conveniently packaged into generating functions: For the antisymmetric representations, the sum terminates at k = N . An SU (N ) character is a class function, in other words, a symmetric function of the eigenvalues ε n of the matrix E. When expressed in terms of eigenvalues, the generating functions can be explicitly written as Bose and Fermi distribution: where (minus) the eigenvalue plays the rôle of the energy level and −ν plays the rôle of chemical potential. Using this representation, we find for the Wilson loop (2.5) in the ksymmetric/antisymmetric representation: For symmetric representations, the contour of integration should be chosen such that C > a j for any j. For antisymmetric representations, C is arbitrary. We have rescaled ν by L for future convenience. The vacuum distribution of eigenvalues a j is characterized by the eigenvalue density: where 14) and f , defined in (2.8), is assumed to be of order one in the large-N limit. The contour integral in (2.13) is of the saddle-point type, and to the leading order in 1 N we find: where ν is determined by minimizing the free energy: This is the standard relation between canonical and grand canonical ensembles in statistical mechanics, with −ν playing the rôle of chemical potential, f the density of particles and ρ(x) the level density. As expected, the exponent in the expectation value of the Wilson loop is proportional to the length of the contour. To compute the coefficient of proportionality, we first need to solve for the eigenvalue density of the matrix model (2.2), then find ν from (2.16) and substitute the result into (2.14). The eigenvalue density for the localization matrix model is not known in general, but at large R and large λ, a systematic perturbative solution can be constructed. The leading-order strong coupling solution was known for long time from the D-brane probe analysis of the dual supergravity background [14] and can be easily reproduced from the saddle-point equations of the matrix model [4]. It has the form of the Wigner distribution: with the width proportional to the square root of the 't Hooft coupling 2 : The leading-order is way too simple to describe the critical behavior of the model, which arises because of the complicated short-distance structure of the density, that at the leading order gets averaged over. Interestingly, the first strong-coupling correction fully reveals the short-distance singularities that are responsible for quantum phase transitions at finite λ [5,6].
Although we expect from holography that the strong-coupling expansion goes in powers of 1 √ λ, the first correction to the density appears at the relative order (M µ) 1 2 ∼ λ −1 4 . The correction has a very irregular, spiky structure ( fig. 1) [5,6]: where {⋅} denotes the fractional part. The zeta-function ζ(1 2, z) has a z −1 2 singularity at zero, and this produces an array of infinite spikes at resonance points x = ±µ ∓ nM .
The physical origin of the spiky structure in the density is due to the resonance on massless hypermultiplets. The number of resonances that fit into the interval (−µ, µ) depends on the ratio 2µ M and each time 2µ crosses an integer multiple of M , a new pair of resonances appears, leading to a fourthorder quantum phase transition [9,10]. The transitions are sharp only in infinite volume, at finite R the spikes are rounded up and the transitions become smooth crossovers. The expression above assumes that the infinite volume limit is taken prior to the strong-coupling limit, in the sense that M R ≫ √ λ. If the strong-coupling limit is taken first, the spikes get damped, albeit at a very slow rate, and completely disappear distance M 2 R ∼ µM R √ λ away from the endpoints [5].
In the bulk of the eigenvalue distribution, the comb-like structure of resonances is a small correction, which moreover disappears upon averaging over a sufficiently wide interval because However, near the endpoints the spikes are no longer suppressed compared to the averaged density. Indeed the two terms in (2.19) become comparable for µ − x ∼ M . The whole expression is actually not applicable near the endpoints. The edge behavior of the eigenvalue density is described instead by [5] ρ  Extremely close to the endpoint, at µ − x ∼ 1 R, the structure of the peak starts to be resolved. The density has the following form in this regime [5]: and is shown in fig. 3. This regime exists only for the theory defined on S 4 . The density near the edgepoints is known even without assuming that M R ≫ 1 [5], but the general expression is substantially more complicated.
With the density at hand, we can now compute the Wilson loop expectation value from (2.16), (2.15), (2.14). Depending on the parameter f , the main contribution to the integrals in (2.16), (2.14) will come from different range of x's. We may identify three distinct regimes: In this case, the average density (2.17) is a good approximation. The spikes constitute a λ −1 4 correction as clear from (2.19).
The spiky structure of the density then appears in the leading order (2.21). Figure 3: The finite-R resolution of the first peak. The density (2.21) with the peak unresolved is shown in the purple line.
• Non-universal regime: The spike singularity is resolved according to (2.22). This regime is only defined on S 4 as the density depends explicitly on the radius of the sphere, even though M R is still assumed to be large.
The first two regimes are universal, in the sense that the results should apply to any sufficiently big contour and do not depend on the S 4 compactification (the radius drops out from the expressions for the density). The last regime requires that the theory is compactified on the sphere, and the computations in that case only apply to the Wilson loop running along the big circle of S 4 .

Antisymmetric representations
The Wilson loop in the antisymmetric representation is expressed through the Fermi distribution, for which 1 L plays the rôle of temperature. When L is large, the effective temperature is low. Then (2.16) and (2.14) can be simplified: and This approximation is justified when the scale of variation of the density is much larger than 1 L, which is true in the bulk and endpoint regimes, because µL ≫ 1 and M L ≫ 1, but not in the non-universal regime when the density varies on scales of order 1 R = 2π L.

Bulk regime
We begin with the bulk regime, when the density is the same as in the Gaussian matrix model. In that case the answer is the same as in N = 4 SYM up to rescaling of λ by a factor of M L 2π, and it reads [16,23] 1 where θ, defined as cos θ = ν µ, is the solution of a transcendental equation The correction term to the bulk density in (2.19) is of order O(λ −1 4 ), but contributes to the Wilson loop vev at relative order λ −3 4 , because the small-scale variations of the density average to zero. We are not going to study this contribution in detail, since it is subleading to the effect produced by the so far unknown order λ −1 2 correction to the density.
From (3.4) we see that the bulk regime corresponds to f = O(1). When f becomes small 3 , the Fermi level ν approaches the edge of the eigenvalue distribution and eventually the endpoint regime sets in.

Endpoint regime
Let us define dimensionless variables u ≡ (µ − x) M and v ≡ (µ − ν) M . It is convenient to express the endpoint distribution (2.21) for different intervals of u, i.e.
To make this scaling manifest we introduce and assume thatf ∼ O(1). The integral in (3.1) can be easily performed by summing over the contribution of different intervals of u that are smaller than v. For example, for n − 1 < v ⩽ n, with integer n, we have: where we have labeled f with a subindex n as a reminder of the domain of v. For the rescaled variables, the explicit forms are: and so on. The functionf (v) is continuous but discontinuous in its first derivative at the critical points, which are its values at integer points, i.e. f critical,n =f n (n):f 4 shows a plot of this function for 0 < v ⩽ 3. The free energy at the saddle point is given by (3.2). We are interested in expressing it in terms of the free parameterf , therefore, first, we must invert the functionf (v). The inverse functions off 1 (v) andf 2 (v) admits a simple closed form: v 3 (f ) can be also expressed in radicals, but the expression is lengthy, and we do not display it here. Now, in terms of v(f ), and using the saddle point equation (3.1) reexpressed in the endpoint variables u and v, the free energy as a function of f is: At strong coupling the first term is dominant, and to the first approximation the free energy is simply F − = µf = µk N . The Wilson loop in the rank-k representation therefore behaves as e µLk , which is just the k-th power of the fundamental representation, as can be expected for small representations on account of the large-N factorization. The second term is of relative order O(M µ) = O(1 √ λ) and in the dual holographic description can be interpreted as a quantum correction on the D-brane worldvolume. This term has singularities at the critical points (3.11), which means that the worldvolume effective field theory may have a non-trivial phase structure even at infinitely strong coupling.
Therefore, the phase structures are smooth, due to the subleading term, see fig. 5.
The phase transitions are of the second order, as can be seen by examining the derivatives of free energy. Using that df dv = M ρ(v), we find: The inverse density is singular for integer values of v, as can be seen in fig. 2, and takes a finite value to the left of the critical point going to zero to the right. The second derivative of the free energy consequently experiences a finite jump across the phase transition. We can see this explicitly by computing the free energy on the first two intervals: . . .

(3.15)
And it is straightforward to check that F − is continuous atf = 1 together with its first derivative, while the second derivative experiences a finite jump.

Non-universal regime
Extremely close to the endpoints, the first peak is resolved, and we need to solve (2.16) with the density from (2.22). From (2.16) we can see that the Wilson loop enters this regime at which is an additional factor √ M R smaller compared to (3.6). Considering again symmetric and antisymmetric representations in parallel, we can see that the free energy scales as The scaling function h ± (f) can be written in the parametric form: The parameter q is the fugacity variable related to the chemical potential ν in (2.16), (2.14) by q = e 2πR(µ−ν) . For the antisymmetric representations (fermions), q changes from zero (low density) at to infinity (high density), while for bosons the chemical potential must be negative and so q is always smaller than one. The low density approximation, q → 0, corresponds to Boltzmann statistics, when the difference between symmetric and antisymmetric representation disappears. The hypergeometric function in the above equations can then be replaced by 1. We thus get: and For the Wilson loop we get the following result, for both symmetric and antisymmetric representations: where W 1 is the expectation of the Wilson loop in the fundamental representation that at strong coupling behaves as [5]: This is indeed the correct behavior that we expect at very small k ≪ N . The Wilson loop then picks the leading contribution from the term with the largest number of traces, due to the large-N factorization. The result above follows from In the high-density regime, symmetric and antisymmetric representations behave very differently. For fermions, f grows indefinitely with q and asymptotes to f ≃ 2π 5 2 √ ln q at large q. The function h − (f) becomes negative and its absolute value grows as the cube of the argument: which can be seen to match (3.15) in the overlap of the region between the two regimes.
In the bosonic case, q cannot be bigger than one. When q approaches one, the integral in (3.18) converges to a finite value, indicating that the solution for q exists only for sufficiently small f. This corresponds to the Bose-Einstein condensation in the analog statistical system. We discuss how to compute the symmetric-representation Wilson loop for arbitrary f in the next section.

Symmetric representations 4.1 Bose-Einstein condensation
The Bose-Einstein condensation for symmetric-representation Wilson loops has been discussed for the Gaussian model [23] of N = 4 SYM [25,26], as well as in the unitary matrix models [29]. To continue the Wilson loop expectation values past the transition point, two prescriptions have been used: (i) analytic continuation in f [23,29] and (ii) observation that for sufficiently large k the rank-k symmetric representations become equivalent to the k-wrapped fundamental representation [20]: ⟨tr R + k e LΦ 0 ⟩ ≃ ⟨tr e kLΦ 0 ⟩ (for sufficiently large k) . The answer for multiply-wrapped Wilson loops is known exactly at any N and k, in the Gaussian model [26]. The scaling limit N → ∞, k → ∞ can be obtained as an approximation to this exact answer [15], and it agrees with the symmetric-representation Wilson loop analytically continued far past the Bose-Einstein condensation point. The equivalence of the two approaches, however, has never been actually proven, and since our model is not Gaussian, we prefer to derive all the necessary results from scratch. Technically, the condensation happens because the hypergeometric function in (3.18) has a logarithmic branch point at q = 1: However, this entails no divergences in f, because the logarithmic term integrates to zero in (3.18) and f approaches a finite limiting value n 2n− 1 2 e −2n n! 2 = 11.88 . . . If f exceeds f c , solution to the saddle-point equation ceases to exist. This happens because the chemical potential of a non-interacting Bose gas cannot be positive, which translates into ν being bigger than µ in our case. The edge behavior of the density of states, ρ(x) ∼ √ µ − x, guarantees the convergence of the integral in (2.16) at ν = µ. Following the analogy with statistical mechanics, we expect that the excess density f−f c will condense in the ground state, identified in our case with the largest eigenvalue max {a j } ≡ a. An important difference to the textbook treatment is the randomness of the energy levels which themselves are integration variables with the measure defined by the eigenvalue integral (2.2).
To take into account the condensate, we need to deform the contour of integration in (2.11) as shown in fig. 6, picking the pole at ν = a: The remaining contour integral is exponentially small compared to the pole contribution, and can be omitted at large N . The average over the bulk eigenvalues can then be replaced by convolution with the saddle-point density, while the integral over the largest eigenvalue should be kept as it is: Here P (a) is the probability to find the largest eigenvalue at point a outside the eigenvalue interval: The constant F 0 is the free energy of the eigenvalue sitting right at the edge of the bulk distribution. The probability to find an eigenvalue at a is determined by the energy cost of taking an eigenvalue from the edge and moving it to a, and it is given by the difference in free energies.
The probability to pull out one eigenvalue outside the bulk distribution is exponentially small, but in the integral (4.5) the negative exponent is counteracted by the positive exponent in the thermodynamic weight. The largest contribution to the integral hence comes from the saddle point: This is similar to (2.14), but now the free energy is given by The position of the saddle point is determined by the following equation: with the kernel: where we used the explicit form on the one-loop measure (2.4), and defined (4.12) The first two terms in (4.10) were absent in (2.16) and can be interpreted as the density of the condensate. Another way to calculate the Wilson loop for f > f c is to continue the free energy analytically from f < f c . Indeed, the logarithmic singularity of (4.2) cancels upon integration in (3.18) and f(q), as a function of q, has no singularities at q = 1. We prove in appendix A that the result of the analytic continuation gives the same results as the calculation above under fairly general conditions, applicable in particular to the model under consideration. The Bose-Einstein condensation consequently does not lead to any singularities in the expectation value of the Wilson loop.
In the universal regime L(a − µ) ≫ 1, the thermal contribution to the free energy (the last term in (4.10)) becomes exponentially suppressed and can be dropped: (4.14) In this approximation, the k-symmetric and k-wrapped Wilson loops have the same expectation value: because the right-hand side is calculated by equation (4.5) without the second term in the exponent [30,20]. This result contradicts intuition based on the large-N factorization, as the Wilson loop now picks the largest contribution from the term with the smallest number of traces.

Bulk regime
When a − µ ∼ µ, the density is the Wigner distribution and the kernel S(x) can be approximated by M 2 R 2 x [4]. The answer is then obtained from the Wilson loop in the Gaussian model [15,23] by simple rescaling: It is interesting to notice that the consistency of this regime requires the rank of representation to scale linearly with the size of the contour. There should be a simple explanation to this fact in the dual supergravity picture.

Endpoint regime
We now consider the case when the saddle point a approaches the endpoint of the eigenvalue distribution. As before, we introduce the dimensionless variables The density near the endpoint has the form (3.5). It cannot be just substituted into (4.14) because the integral would then diverge. The strategy is to separate the bulk and endpoint contributions to f as follow: where f ∞ is computed with the Wigner density and with S(x) approximated by M 2 R 2 x: (4.20) Notice that the rank of representation scales linearly with the length of the contour as in the bulk regime, but its coupling dependence is different: f ∼ λ −3 4 instead of λ −1 2 . We thus introduce the rescaled variablẽ such thatf For the genuine endpoint contribution we find: In spite of proximity to the endpoint, the argument of the kernel is still very big: a−x = M (u+v) ≫ 1, which justifies the use of an approximate expression Writing the endpoint density (2.21), (3.5) as and introducing the function we find thatf ep = X(v). (4.28) The function X(v) can be computed by noticing that and differentiating (4.27) twice, after which elementary integrations yield: We thus get: where the integration constants are set to zero in order to get the correct behavior at v → ∞. Interestingly, the function X(v) has a branch-point singularity outside the eigenvalue support, again due to the resonance appearing in the kernel S(x). Combining this result with (4.22), we find: (4.32) An alternative derivation based on the Wiener-Hopf method [5] is given in the appendix B.
Recalling that the equation for f was obtained by differentiating the free energy with respect to a, we can get the free energy by integration: As before, we would like to express the free energy in terms off . We need then to obtain the inverse function v(f ). For v < 1 or equivalentlyf < 2 3, the answer is very simple: At the critical pointf c = 2 3 the free energy experiences a phase transition. It is easy to show that F is continuous together with its two first derivatives, while the third derivative diverges at the critical point: The phase transition is consequently of the third order.

Non-universal regime
The general scaling dependence of the variable f and of the free energy on M R and λ have already been discussed in sec. 3.3: Explicit integral representation for the free energy can be obtained by plugging the scaling form of the density (2.22) into (4.9), (4.10). The Fourier representation of the the resulting integrals is given in appendix B.

Conclusions
The leading-order strong-coupling solution of the SYM* matrix model is given by the same Wigner distribution as in the N = 4 SYM, up to rescaling of the coupling λ → λL 2 4π 2 . Wilson loop expectation values in the bulk regime, as a consequence, are obtained from the Gaussian results by simple rescaling. This is quite surprising, as the dual geometry of SYM* is rather different from AdS 5 ×S 5 away from the boundary. It would be very interesting to verify these basic matrix-model predictions by finding classical D-branes solutions in the Pilch-Warner geometry. Even more interesting phenomena occur for parametrically small representations. The leading order then is just the result expected form the large-N factorization. It is the next, √ λ suppressed term that displays a non-trivial behavior, with phase transitions in the representation parameter f = k N , appropriately rescaled. In the dual picture, f maps to the density of electro-magnetic flux on the D-brane worldvolume. Perhaps the non-analytic behavior of the Wilson loops can be interpreted as phase transitions in the effective field theory on the D-brane worldvolume with the electric or magnetic flux playing the rôle of tunable parameter. These phase transitions and quantum phase transitions in the matrix model at finite but large λ have the same origin.
We have mainly concentrated on the infinite-coupling regime of the theory. The discontinuities in the higher-rank Wilson loops should persist at finite coupling as soon as there are resonance peaks in the density 4 , which is the case for λ > λ

A Analytic continuation
Consider an eigenvalue model defined by the partition function The eigenvalue density for this model satisfies the saddle-point equation We assume that such that the kernel in the saddle-point equation is of the Hilbert type, and the density has the usual square-root asymptotics near the end-points. It will be important for the subsequent derivation that the residue of the kernel is exactly two. If it is different from two, the analytic continuation will not describe Bose-Einstein condensation correctly. The k-symmetric Wilson loop in this model is determined by the free energy 5 The chemical potential ν is determined by minimization of the free energy: which implicitly determines F as a function of f . These formulas are literally valid for sufficiently small f < f c , such that (A.6) has a real solution. Our goal is to analytically continue the free energy past the critical point. We first consider the analytic structure of F ′ as a function of ν. As such, it has a cut from −µ to µ with the discontinuity Consequently, near ν = µ the function F ′ (ν) has a double expansion that combines analytic and square-root terms: The free energy is obtained by integration: The constant g 0 must be positive, due to positivity of the eigenvalue density. The solution to (A.6) thus only exists for f < f c .
When f c − f is small, (A.6) can be solved iteratively with the help of the expansion (A.8). Substituting the solution back to F(ν), we get to the first approximation: While F(ν) has a square-root singularity, the function F(f ) is analytic and can be continued past f c , which is just an inflection point of the free energy.
The analytic continuation in f is equivalent to continuing F(ν) in ν through the square-root branch cut that changes the sign of f c − f to negative. We denote the result of such analytic continuation byF(ν). The functionF(ν) has the same continuous part as F(ν) but the opposite discontinuity: and is completely characterized by these two conditions.
We can constructF(ν) by introducing an auxiliary function, the generalized resolvent: where F 0 is just a constant, to be determined later. The saddle-point equation (A.2) is equivalent to the two conditions: where the second equation is a consequence of (A.4). This implies, in view of (A.7), that Strictly speaking, the first equation in (A.13) only implies that the continuous part of W (x) is a constant independent of x, but this constant can be adjusted to zero by setting Indeed, this function satisfies the right discontinuity condition (A.11). Explicitly,F The function W (ν) can be interpreted as the energy cost of pulling one eigenvalue out of the bulk of the distribution and placing it at point ν. Comparing (A.17) with (4.9) we see that the analytic continuation of the free energy gives the same result as Bose-Einstein condensation on the largest eigenvalue. The Bose-Einstein condensation thus does not lead to any thermodynamic singularities in the expectation values of the k-symmetric Wilson loops.

B Wiener-Hopf method
In this appendix we derive the general formulas for the symmetric Wilson loop (4.9), (4.10) from the Wiener-Hopf method developed to solve the saddlepoint equations for the matrix model in [5], assuming that a is very close to µ. This covers both endpoint regime and the non-universal regime. As in [5] we use the units in which R = 1 in this appendix, and introduce the endpoint variables ξ = µ − x and η = µ − a < 0. Parametrizing the endpoint distribution as: we can identify three contributions to f in (4 .10): is the thermodynamic term. Here we concentrate on the endpoint part, defined in terms of the function Our main observation is that X(η) is the same function that was introduced in [5] to solve the integral equation of the matrix model, where it plays the role of the remainder function in the Wiener-Hopf method. This function is different from zero only for η < 0 (its Fourier image in analytic in the lower half-plane). In the Fourier space, X(ω) can be read off from the results in [5] (we keep the original notation): We see that X(ω) does not have poles anymore. This is because the simple poles got canceled by the zeros of the prefactor. The double poles are not relevant at this regime, as Γ iω 2π is approximated by its expansion at zero. Only the cut remains in the upper half complex plane.
Straightforward contour integration gives us the expression in the coordinate space, which is: Upon rescaling, this is the same as eq. (4.31) in the main text.

B.2 Non-universal regime
The solution (B.12) for X(η) obtained in the regime O(M ) actually approximates well this regime too. In Fourier space, following the simplifications explained in the appendix C of [5], X(ω) reduces to: Due to its complicated analytical structure, especially the double poles and the logarithmic branch cuts in the exponential, we abstain from computing its exact expression in coordinate space. Instead, we study its asymptotic behaviors, for ω ≫ 1 and ω ≪ 1. For the former limit, the first term of the sum in (B.13) dominates, which is the same as the first term in (B.11) from the O(M ) regime. For ω ≪ 1, the leading contribution comes from the second term of the sum in (B.13), and we expect matching with the regime O(M ), which it does. Here are the asymptotic results in the coordinate space: Numerically, we also see that (B.12) approximates this regime very well at large M.