Finite-temperature phase transitions of third and higher order in gauge theories at large $N$

We study phase transitions in $SU(\infty)$ gauge theories at nonzero temperature using matrix models. Our basic assumption is that the effective potential is dominated by double trace terms for the Polyakov loops. As a function of the various parameters, related to terms linear, quadratic, and quartic in the Polyakov loop, the phase diagram exhibits a universal structure. In a large region of this parameter space, there is a continuous phase transition whose order is larger than second. This is a generalization of the phase transition of Gross, Witten, and Wadia (GWW). Depending upon the detailed form of the matrix model, the eigenvalue density and the behavior of the specific heat near the transition differ drastically. We speculate that in the pure gauge theory, that although the deconfining transition is thermodynamically of first order, it can be nevertheless conformally symmetric at infinite $N$.


I. INTRODUCTION
The nature of the deconfining phase transition in SU (N ) gauge theories is a question of fundamental importance; numerical simulations on the lattice indicate a transition of first order for N ≥ 3 [1]. In finite-temperature pure gauge theory, the Polyakov loop is the relevant order parameter. It is therefore reasonable to study the phase transition as a function of an effective theory of thermal Wilson lines, as a type of matrix model.
There are many matrix models which are soluble at large N . The most familiar is when the transition is driven by the Vandermonde determinant from the integration measure of a single site integral [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. This type of model was originally applied to a lattice gauge theory in two space-time dimensions [3][4][5], where, for the Wilson action, there is a thirdorder transition as a function of the coupling constant. The third order transition as a function of temperature was subsequently shown in lattice gauge theory at strong coupling with heavy quarks using the mean-field approximation [6,7].
Sundborg showed that at infinite N , this model is relevant to deconfinement on a femtosphere, S 3 × R 1 [10][11][12][13][14][15][16][17]. As a function of temperature, the deconfining transition appears to be of first order, as both the energy density and the order parameter are discontinuous at the transition temperature T d . Even so, the specific heat diverges as T → T + d , as is typical of a second order transition [12,13,19,20]. For this reason the transition at infinite N can be termed "critical first order" [12,13].
On the femtosphere, Aharony et al. showed that the critical first order is washed out when higher order perturbative corrections are included, leaving an ordinary first order transition [11]. The question is whether this remains true in the limit of infinite volume. An effective model of Wilson lines was developed as a model for deconfinement with three colors [21][22][23][24][25][26][27][28][29] and extended to include dynamical quarks [30]. This model is soluble at large N [19]. In this paper we study a general class of matrix models and solve them in some special but illustrative cases.
We now give an outline of this paper and summarize the main results. Matrix models for the deconfining phase transition are functions of the thermal Wilson line, L = P exp(ig 1/T 0 A 0 dτ ), which we take to lie in the fundamental representation. The general form of the effective potential, which we discuss in Sec. II, includes an infinity of terms.
The simplest possibility is to start with those involving arbitrary powers of L, but just two traces, a n |trL n | 2 + . . . . (1) We assume that these double-trace terms dominate near the deconfining phase transition T d . All a n are therefore positive below T d to prevent the spontaneous symmetry breaking of Z(N ).
Section III considers deviations from the double trace terms of Eq. (1). We assume that the phase transition is driven by traces of loops which wind only once in imaginary time, tr L, and not by those which wind more than once, such as tr L 2 , etc. This is a significant assumption, but is borne out by all known models. Then the obvious terms to add next are those linear and quartic in tr L: We assume that the coupling for the quartic term, b 1 , is small near T d , but let the magnitude of the background field, h, be arbitrary. We ignore all other couplings, including cubic couplings such as tr(L −1 ) 2 (trL) 2 + c.c., etc., and discuss why this might be valid.
For the original model of Gross, Witten, and Wadia [3], the coefficients a n = 1/n, and there is a phase transition when the expectation value of the loop, 1 N tr L, equals 1 2 . In Sec. III, we consider arbitrary a n , and show that for some critical h c , there is always a phase transition when 1 N tr L goes through 1 2 . We thus term this point a generalized Gross-Witten-Wadia (GWW) transition. In Sec. III A we show that there is a region in the space of a 1 , b 1 and h where a generalized GWW transition occurs.
About the GWW transition, the value of the potential at the minimum can be expanded in powers of δh = h − h c : where v is an irrelevant constant, and f reg is a smooth function of h. For the model of Gross, Witten and Wadia, r = 3, and thus the transition is of third order in h. In Sec. III B we argue that independent of the values of the coefficients a n , r > 2 for the generalized GWW transition, and so the transition is of higher order than second. The point where h = b 1 = 0 is special, as lines of first, second, and higher order meet, Fig. (1).
We stress that as in the model of Gross, Witten, and Wadia, this unusual phase structure emerges only at infinite N . This can be seen from Eq. (3): at infinite N the piecewise function emerges. For finite N , however, the corresponding function is regular at δh = 0, and so the corresponding transition at 1 N trL = 1 2 is only a smooth crossover. In Sec. IV we solve a specific class of models. We take the coefficients that contain a simple form, In the conclusions, Sec. V, we discuss the possible implications of our results. In particular, we speculate that in infinite volume, that the deconfining transition can be critical first order.

II. EFFECTIVE POTENTIAL
In this section, we consider the effective potential of the Polyakov loop in SU (N ) Yang-Mills theory at large N . After some standard definitions, we consider the possible forms for the effective potentials of Polyakov loops, and suggest that double trace terms may dominate near the deconfining phase transition.

A. Notation
In a Yang-Mills theory without dynamical quarks at nonzero temperature, the global symmetry associated with the deconfining phase transition is Z(N ). The basic variable is the Wilson loop in the direction of imaginary time, τ , Under the center symmetry the thermal Wilson line transforms as L → z L where z is an element of Z(N ), z = exp(2πij/N ) for an integer j = 1 . . . N .
While the thermal Wilson line is gauge variant, its eigenvalues are gauge invariant. As a unitary matrix, after diagonalization As an SU (N ) matrix, these eigenvalues satisfy N i=1 θ i = 0, modulo 2π. The phase transition is then characterized by the traces of powers of the thermal Wilson line. Without loss of generality we can take all traces to be in the fundamental representation, so there are N − 1 independent Polyakov loops, The n th Polyakov loop ρ n wraps around in imaginary time n times. The ρ n form a complete set of gauge invariant order parameters for the spontaneous breaking of Z(N ) symmetry in the deconfined phase [22,31].
At large N we introduce the variable x, where θ i → θ(x) [2]. The n th Polyakov loop is then At infinite N each loop is a functional of θ(x). Introducing the eigenvalue density the loop becomes In this way, the Polyakov loops are functionals of ρ, rather than of θ.
The eigenvalue density must be non-negative, This will play an essential role for the GWW phase transition. It is normalized as Provided that ρ(θ) is continuously differentiable for −π ≤ θ ≤ π, we can write ρ as a Fourier series in terms of its moments in θ, using Eq. (14). We have assumed that the expectation value of every Polyakov loop ρ n equals its complex conjugate, ρ −n . They are related under charge conjugation, and so this assumption is valid at nonzero temperature and zero quark chemical potential. (At nonzero chemical potential [8] the expectation value of the loop and its complex conjugate differ [16,17,23,32]). In doing so we implicitly perform an overall Z(N ) rotation so that the expectation value of all Polyakov loops are real.
At infinite N , all Polyakov loops vanish in the confined phase, ρ n = 0 for all n ≥ 1. This implies that the eigenvalue density is constant, ρ(θ) = 1/(2π), demonstrating the complete repulsion of eigenvalues. In the deconfined phase, ρ n become nonzero. At infinitely high temperature all Polyakov loops equal unity; this implies that all eigenvalues become zero, and thus ρ(θ) = δ(θ).

B. Effective potentials for the Polyakov loops
The effective potential for the Polyakov loop is constructed formally as [33] exp where S YM is the d-dimensional Euclidean Yang-Mills action and V is the spatial volume.
The right hand side is a path integral over the gauge fields in the presence of constant background Polyakov loop ρ n . We have scaled out VT d−1 so that V eff is dimensionless. By convention, there is a factor N 2 , so that N 2 V eff is of order N 0 in the confined phase and N 2 in the deconfined phase [34]. The partition function is then given by integrating over the eigenvalues: where we defineN At the outset, we stress that at infinite N the effective potential is a function of all loops, for every ρ n from n = 1 to ∞. One might hope to simplify things by considering an effective potential for just a few loops, such as ρ 1 , ρ 2 , etc. However, all Polyakov loops vanish in the confined phase, ρ n = 0 for every n ≥ 1; an effective potential involving a finite number of loops cannot force all loops to vanish.
Alternately one could integrate over all ρ n for n ≥ 2, and construct an effective potential just for ρ 1 . This is possible, but its utility is not apparent to us.
Using only the global Z(N ) symmetry, there are many terms which can appear in the most general effective potential for the Polyakov loops, where a (2,1) (1,−2) = a (2,1) (−1,2) by charge conjugation symmetry. With such a multitude of terms, this effective potential is not of much use. Thus we discuss some results from perturbative computations, and from effective models, which suggest that the effective potential relevant at infinite N may be much simpler.
Consider first the computation of the effective potential in perturbation theory. This is a straightforward matter at one loop order [35,36] and has been carried out to two loop order [33,37,38]: The notation d (4) n denotes the deconfined term, computed perturbatively in four spacetime dimensions.
The simplicity of this result is not obvious. Naïve computation to two loop order gives a result which is much more complicated than that at one loop order [39]. After including a finite renormalization for Polyakov loops [33,37,38,40], however, one finds that all terms collapse to Eq. (20), just a constant times the result at one loop order.
What is remarkable about Eq. (20) is that the only terms which enter involve two traces: just i 1 = −i 2 and j 1 = j 2 = 1. At one loop order this is automatic, but at two loop order terms with three traces can appear, ∼ g 2 N . After including the finite renormalization of the Polyakov loops [33,37,38,40] these vanish.
A similar computation to one loop order in 2 + 1 dimensions [41] shows For a field in d spacetime dimensions, the corresponding term is a (1,1) (n,−n) ∼ −1/n d . We discuss results to higher loop order in Appendix A.
The sign of the double trace terms in Eqs. (20) and (21) is negative. Thus the potential is minimized by maximizing each Polyakov loop, which is what is expected in the perturbative regime. To model the transition to a confining phase, it is necessary to add additional terms.
In Refs. [26], an effective model was constructed by adding terms which mimic deconfined strings in 1 + 1 dimensions: We denote this term c (4) n as that which drives confinement in 3 + 1 dimensions. This potential can be also derived by adding a mass deformation to gluons [22]. We comment that this particular term is driven by detailed results from numerical simulations on the lattice, especially the presence of terms ∼ T 2 in the free energy relative to the usual ∼ T 4 [22,26].
The c (4) n have positive signs, and if sufficiently large, drive a transition to a confined phase. We conclude this section by discussing other evidence for the dominance of double trace terms in the effective potential for the thermal Wilson line. This is true in a strong coupling expansion on a lattice, at least to leading order [42].
In super Yang-Mills theories with mass deformations, it is possible to compute not only the perturbative contributions to the free energy, but also the dominant non-perturbative terms [43,44]. While it is not obvious [44], it can be shown that in these models all terms for the effective potential of thermal Wilson lines only have double traces.
The only instance of which we are aware in which terms with three and four traces arise is for gauge theories on a femtosphere, S 3 × R 1 [11]. Explicit computation to three loop order shows there are a variety of terms, including those with four traces, although they are suppressed by ∼ g 4 N 2 , where the coupling constant g 2 N is small on a femtosphere.
Consequently, in the following we assume that any terms with three and more traces are small, and compute about that limit. In the conclusions, Sec. V, we consider the implications if there are only double trace terms in the limit of infinite volume.

III. PHASE STRUCTURE
Motivated by the above considerations, we are led to consider the following effective potential for a gauge theory at infinite N : We assume that the deconfining transition at a temperature T d defined at h = 0 is driven by the first Polyakov loop. In the confined phase, that the first Polyakov loop vanishes does not necessarily imply that all other Polyakov loops vanish. For this reason, as discussed above it is essential that we include all Polyakov loops in the effective potential and require that all a n are positive below T d . We further assume that a n > 0 for n ≥ 2 (24) near T d , so the higher corrections for ρ n with n ≥ 2 are not necessary.
We include a quartic term for the first Polyakov loop, with a coupling b 1 which we assume is small and constant near T d . We also include a background field for the first Polyakov loop, ∼ h. This is natural to include in any spin model, and we can analyze the model for arbitrary values of h. A background field for Polyakov loops is generated by the coupling to quarks [15]. If the quarks are heavy, with a mass m, then the background field is ∼ exp(−m/T ) for the first Polyakov loop, ∼ exp(−2m/T ) for the second, and so on. Thus very heavy quarks only generate a background field for the first Polyakov loop. For N f flavors of quarks

A. Phase diagrams
In this subsection we derive general conclusions about the phase diagram of Eq. (23) with the condition Eq. (24) near T d . In Sec. IV we then solve the series of models where the a n have specific values, a n ∼ 1/n s for s = 1, 2, 3 and 4.
At infinite N we look for a saddle point of V eff (ρ n ) under the constraints of Eqs. (13) and (14). By a Z(N ) rotation we can assume that the expectation value of all Polyakov loops are real, so ρ * n = ρ −n = ρ n . Naïvely, the saddle point corresponds to the minimum of each free energy for ρ n , It is easy to solve this equation, taking all loops beyond the first to vanish, ρ n = 0 for n ≥ 2. The eigenvalue density in Eq. (15) is then a sum of a constant and ρ 1 , This satisfies the normalization condition of Eq. (14), but it is non-negative only if the first Polyakov loop is less than or equal to one half, ρ 1 ≤ 1 2 . Therefore, this solution is valid only for 0 ≤ ρ 1 ≤ 1 2 . This is the simplest way to see that the point where the first Polyakov loop equals one half and all others vanish, ρ 1 = 1 2 and ρ n = 0 for n ≥ 2, is special. We call this the Gross-Witten-Wadia (GWW) point, and the locus of such points is a GWW surface.
When the expectation value of the first Polyakov loop is greater than 1 2 , expectation values for all higher loops develop. This is not due to the usual manner of Landau mean field, through the coupling of ρ 1 to the other ρ n through terms such as (ρ * 1 ) 2 ρ 2 , etc. Instead, the eigenvalue density becomes no longer continuously differentiable due to the non-negativity constraint, and as a result higher Polyakov loops become nonzero. In the model of Gross, Witten, and Wadia [3], and for the models of Sec. IV, this happens by developing a gap in the eigenvalue density.
If the first Polyakov loop has an expectation value less than 1 2 , we can use an effective theory for just that loop, ρ 1 : Consider first zero external field, h = 0, as illustrated in Fig. (1). If a 1 and b 1 are positive, the minimum is clearly for ρ 1 = 0. If a 1 is negative and b 1 positive, the minimum is . Thus there is a second order phase transition when a 1 vanishes. This is indicated by the blue dash-dotted line in Fig. (1).
As a 1 decreases for a fixed positive value of b 1 , the Polyakov loop equals 1 2 when b 1 = −2a 1 . At this point, it is no longer possible to include only the first Polyakov loop in the effective theory. This is denoted by the green solid GWW line in Fig. (1).
The location of the first-order phase transition depends on the explicit form of a n . The red At the origin a 1 = b 1 = h = 0, the first, second and higher order phase transition lines meet. At this point, the Polyakov loop ρ 1 jumps from 0 to 1/2, as is typical of a first order phase transition, while the mass associated with ρ 1 becomes zero, as is typical of a second order phase transition. This point was termed as "critical first order" in Refs. [12,13].
The center symmetry is broken by a nonzero background field h = 0, which thus washes out a second order phase transition. About ρ 1 = 1 2 , we introduce where −1/2 ≤ δρ 1 ≤ 0. This is equivalent to the Legendre transform Γ(ρ 1 ) of the effective potential below the GWW point. As we argue in general in Sec. III B, and show explicitly in Eq. (101), the GWW point is a continuous phase transition, whose order is always higher than second. Consequently, the coefficients up to and including δρ 2 are continuous about the GWW point.
We analyze Eq. (28) as follows. At the GWW point δρ 1 = 0, and two conditions need to be satisfied. First, the coefficient of δρ 1 must vanish, so that a 1 + b 1 /2 − 2h = 0; second, that the coefficient of δρ 2 1 must be positive, a 1 + 3b 1 /2 > 0. This forms a surface of GWW points in the space of a 1 , b 1 , and h. The green shaded region in Fig. (1) indicates the projection of the GWW surface onto the h = 0 plane. The GWW surface is independent of the coefficients a n . The green solid lines in Fig. (2) are the cross-sections of the GWW surface for b 1 = 0 and −0.08.
As b 1 is decreased for a fixed, positive value of a 1 along the GWW surface, we eventually hit the boundary where the coefficient of δρ 2 1 vanishes, b 1 = −2a 1 /3. Beyond this point, ρ 1 = 1 2 is an unstable solution, and there is a first order transition. This is indicated by a red dashed line in Figs. (1) and (2). The red hatched region in Fig. (1) is the projection of The Polyakov loop ρ 1 becomes larger than 1/2 above the first order phase transition or GWW point. In this region, the effective potential is not just a function of the first Polyakov loop, ρ 1 , but of all ρ n . To describe the theory beyond the GWW point, we need to know the explicit values of the coefficients a n . We can show, however, that the GWW point is a phase transition point for arbitrary a n .   Fig. (1), versus the background field h and the quadratic coupling a 1 . We also illustrate the (small) difference between s = 1 and s = 4 in the figure on the right hand side.

B. The order of phase transition at the GWW point
In this subsection we argue that about the GWW point, there is a continuous phase transition whose order is always higher than second. The partition function in Eq. (17) can be written as where F is the dimensionless free energy in the presence of the external field h per volume V and per the color degrees of freedom N 2 . F is a generating function for the Polyakov loop, where the expectation value of ρ 1 is given by the factor of 2 accounts for the complex conjugate of the first Polyakov loop. Consequently, the free energy is the integral of the loop with respect to h, Since F is the value of the potential at a saddle point of V eff at large N , we have where the n th Polyakov loop, ρ n (h), satisfies the equation of motion, We now have two expressions for the free energy, Eqs. (31) and (32). For completeness we give another form of the free energy in Appendix B when V eff is given as in Eq. (35).
In order to explore the order of phase transition about the GWW point, we only need to look at one point in the green shaded region in Fig. (1). We choose the point where the quartic coupling vanishes and the quadratic couplings are positive, b 1 = 0 and a n > 0 for all n ≥ 1. Without loss of generality we choose a 1 = 1. The effective potential becomes We call it the GWW potential. This potential naturally appears after the Legendre transform of the full potential as shown in Sec. IV A. The equation of motion is n a n ρ n sin(nθ) , by using δρ n /δθ = −n sin(nθ) from Eq. (10). As we discussed in the previous subsection, the equation of motion (36) below the GWW point is satisfied if By Eq. (15), when 0 ≤ h ≤ 1 2 . The GWW point is when h = 1 2 . We write the Polyakov loop just below the GWW point as ρ 1 = 1 2 + δh where δh = h − 1 2 . Note that δh is negative below the GWW point. Using Eqs. (31) or (32), the free energy for the GWW potential V GWW is Next consider just above the GWW point, h = 1 2 + δh with 1 δh > 0. Writing the Polyakov loop as ρ 1 = 1 2 + δρ 1 , the equation of motion is n a n ρ n sin(nθ) .
At small δh > 0 the leading term for the first Polyakov loop is where u and q are some constants, with q ≥ 0. From Eqs. (31) and (32), provided that δρ 1 ∼ u δh q . If ρ 1 = 1 2 and thus δρ 1 = 0 above the GWW point, then the two expressions for the free energy are equal only if ∞ n=2 a n ρ 2 n vanishes. This implies that all higher Polyakov loops vanish, ρ n = 0 for n ≥ 2, which violates the equation of motion in Eq. (40) when δh > 0. Therefore ρ 1 = 1 2 above the GWW point. On the other hand, if u is nonzero and 0 ≤ q < 1, then we can compare Eqs. (42) and (43) to obtain ∞ n=2 a n ρ 2 n ∼ −u 2 δh 2q .
This is not consistent, because the a n are positive and the ρ n are real. Therefore Comparing Eqs. (39) and (42) just below and above the GWW, we see that only the second or higher derivatives of the free energy are discontinuous. Hence the phase transition is of second or higher order.
We now exclude the possibility of a second order transition. If q = 1 and u = 1, or q > 1, then the second derivative of the free energy is discontinuous at the GWW point.
This implies that the mass for the first Polyakov loop is discontinuous at the GWW point.
The mass for the n th Polyakov loop below and at the GWW point is This is a diagonal matrix whose elements are nonzero. Therefore, if the phase transition is of second order, any mass eigenvalue is nonzero at the GWW point. This is not expected for a second order transition, where the critical fields are massless.
The remaining possibility is that q = u = 1. Then δρ 1 ∼ δh, and to leading order, the first term in Eq. (40) vanishes. By comparing Eqs. (39) and (42), the free energy is continuous up to δh 2 near the GWW point. If we consider the term at next to leading order for δρ 1 above the GWW point and use the same argument as before, where v is a nonzero constant. Hence This implies that the order of the transition is higher than second. This is valid for the explicit solutions in Sec. IV.

IV. MODELS
We now solve certain models with special values for the coefficients a n in Eq. (23). We confirm the general phase structure discussed in the previous section, and compute how the behavior of the Polyakov loops and thermodynamic quantities change.

A. Special cases
We consider a simple class of models which are exactly soluble at large N : with s = 1, 2, 3 and 4. ρ n is the n th Polyakov loop, Eq. (8). This is to take the positive and negative parts of the coefficients, a n = c n − d n , in Eq. (23) as where c 1 and d 1 are dimensionless positive functions of T . The coefficients c and d denote the "confined" and "deconfined", respectively, because they are repulsive and attractive potentials for the eigenvalues as mentioned in Sec. II.
Using the identity we see that the effective potential involves the polylogarithms of order s.
For odd s larger than one, the polylogarithm functions do not simplify further.
To compute thermodynamic quantities at nonzero b 1 we perform a Legendre transform [13,14]. Using the effective potential of Eq. (49), in the partition function in Eq. (17) we introduce the constraint δ(λ − ρ 1 ), whereN 2 = VT d−1 N 2 as before, and ω = −iω/(2c 1N 2 ) after a Wick rotation. The delta function constraints the configurations in the path integral to be such that the first Polyakov loop is real. This is valid because we are only interested in the saddle point where ρ n = ρ * n . We define where In Sec. III B we showed that the second derivative of the free energy, −d 2 F GWW /dω 2 , is positive, and verify this later in Eq. (99). Therefore we can perform a Legendre transform of the free energy F GWW and write the partition function as where Using Eq. (58) into Eq. (55), the total partition function becomes where Once we obtain Γ GWW (ρ 1 ) and thus Γ(ρ 1 ), all thermodynamic quantities can be computed for given values of b 1 , c 1 , d 1 , and h.

B. The GWW potential
In this subsection we solve for the eigenvalue density ρ(θ) for the GWW potential V GWW in Eq. (57). The detail derivation is given in the next subsection. The potential is equivalent to Eq. (35) by identifying a n = 1/n s and taking the background field h = ω. We solve the model for s = 1, 2, 3, and 4.
Using Eq. (10), the potential can be written in the form The corresponding equation of motion is where it is convenient to introduce the eigenvalue density ρ(θ). This is equivalent to Eq. (36).
It is necessary to solve the equation of motion under the two constraints of Eqs. (13) and (14).
For the potential of Eq. (62), the GWW point corresponds to ω = 1 2 . Consider first the case below and at the GWW point, where ω ≤ 1 2 . It is trivial to solve for the eigenvalue density, thus ρ 1 = ω, and ρ n = 0 for n ≥ 2 in Eq. (26). We plot ρ(θ) for ω = 0, 1 4 , and 1 2 in Fig. (3a). Notice that when ω ≤ 1 2 the eigenvalue always extends from −π to +π. One can also see from Eq. (65) that this solution is not consistent for ω > 1 2 , as then the eigenvalue density is negative for some range of θ about ±π.
Below the GWW point, ρ 1 = ω, so the Legendre transform of F GWW (ω) = −ω 2 is from Eq. (59), where δρ 1 = ρ 1 − 1 2 < 0. We write it in this manner because it will be useful in comparing to the behavior above the GWW point later.
The properties below and at the GWW point are independent of the model. This is not true above the GWW point, and we need to solve the equation of motion for each value of s. In the remaining part of this subsection, we summarize and explain the solutions. The eigenvalue densities above the GWW point are illustrated in Fig. (3b) for ω = 1.
When ω > 1 2 , a gap opens up in the eigenvalue density for all s, so that it no longer runs from −π to π, but instead from −θ 0 to θ 0 . The value of θ 0 differs for each value of s. When s = 1 and ω = 1, the eigenvalue density runs from − π 2 to + π 2 , and vanishes at the ends. When s = 2, the eigenvalue density jumps discontinuously from a zero to a nonzero value at the ends. When s = 3 and 4, the eigenvalue density actually diverges at the endpoints. Physically, as s increases, the change in behavior for the eigenvalue density occurs, because eigenvalue repulsion weakens as s increases. In particular the value of θ 0 for a fixed value of ω becomes smaller as s increases. When s = 4 the repulsion is so weak that the eigenvalues pile up at the endpoints ±θ 0 , and the density becomes delta function at the endpoints, as the arrows indicate in Fig. (3b). In other words, the eigenvalue θ(x) for s = 4 is no longer an injective function of x. This implies that there is a critical value of s = s * with 3 < s * ≤ 4, above which the eigenvalue repulsion is not strong enough to keep all eigenvalues separate above the GWW point. This makes the analysis of the case s ≥ 4 difficult, and we do not know if the solution exists above the GWW point in the case s ≥ 5.
For s = 1, 2 and 3, all eigenvalues collapse to zero and the density becomes ρ(θ) = δ(θ) as ω → ∞. When s = 4, all eigenvalues collapse to zero at a finite value of ω as discussed below Eq. (96), again as a consequence of the weak eigenvalue repulsion.

C. Derivation of the solution for the GWW potential
In order to solve the equation of motion (64) with s = 1, 2, 3, 4 for the eigenvalue density, we use a trick from Refs. [5,19]. Let us explain how to solve it when s is a positive integer in a naïve way. What makes the equation of motion difficult to solve is the combination of the integral over θ and the sum over 1/n s−1 . However, if we differentiate with respect to x, d dx = dθ dx d dθ , the sum then becomes 1/n s−2 . Doing this s − 1 times, we end up with an integral equation for ρ(θ ), which is soluble. Each time we differentiate with respect to θ we change sin(n(θ − θ )) to cos(n(θ − θ )). If we take s − 1 derivatives, we then end up with a different function depending upon whether s is even or odd.
This approach breaks down if dθ dx = 0 for some domains of x, i.e. if eigenvalues pile up. As mentioned in the previous subsection, the pileup occurs at the endpoints in the case of s = 4, so we need to treat this case with care. We first solve for odd s, and then even. In order to solve the equation of motion (64) for s = 3, we differentiate it with respect to θ twice. Using Eq. (53), we can write the first derivative as 20 The second derivative is This is the equation of motion for s = 1, and is the circular Hilbert transform of ρ(θ ) with the kernel cot( θ−θ 2 ) [45,46]. The integral is singular when θ = θ , and so implicitly it is defined by using a principal value prescription. The eigenvalue density is the inverse of the transform, where the constants satisfy At the GWW point θ 0 = π, and the eigenvalue density is that of Eq. (65) with ω = 1 2 . Above the GWW point, a gap opens up, with the density nonzero only between −θ 0 and θ 0 .

Even s = 2 and 4
As with odd s, for even s the GWW phase transition is characterized by a gap in the eigenvalue density.
There is a qualitative difference between s = 4 and all the other cases. This difference was first discovered on the basis of numerical analysis for s = 4 and N = 55 in Eq. (56).
As can be seen in Fig. (4), the eigenvalues are separate for s = 2, while they pile up at the endpoints for s = 4.
This suggests that at infinite N , the eigenvalues θ(x) with x 0 2 ≤ x ≤ 1 2 become a single value θ 0 = θ( x 0 2 ) for some x 0 , and likewise for the other endpoint −θ 0 . Therefore, the expectation value of a function f (θ) can be written as where ρ(θ) is a smooth function defined for the interval −x 0 /2 ≤ x ≤ x 0 /2. This ansatz corresponds to the following eigenvalue density The two parameters x 0 and θ 0 are related by the normalization condition, which can be derived by setting f = 1 in Eq. (81): In order to solve for ρ, we take s − 1 derivatives of the equation of motion (63) with Using the Poisson summation formula and the normalization condition (14), we obtain ρ(θ) = 1 2π (1 + 2ω cos θ) .
We now solve for the two unknowns, θ 0 and x 0 , using the equation of motion (64) with Eqs. (82) and (86), and the normalization condition (83), which can be now written as by using Eq. (86).
To solve the equation of motion for s = 2 we use the identity where −2π ≤ φ ≤ 2π. This is proportional to the derivative of the Bernoulli polynomial B 2 (|φ| /(2π)) given in Eq. (52). The equation of motion is satisfied for any value of x 0 in Eq. (82). It turns out that the one minimizes the potential is when x 0 = 1, i.e., when there is no pile up of eigenvalues. This is consistent with the case of finite but large N as shown in Fig. (4a). Therefore, the eigenvalue density is where −θ 0 ≤ θ ≤ θ 0 . The is equivalent to the density below the GWW point (65), except here it has a gap at the endpoints. This solution is consistent with the one found in [19].
The order of the discontinuity with respect to ω depends upon s: as illustrated in Fig. (5), for s = 1 the third derivative is discontinuous; for s = 2 and 3, the fourth derivative; and for s = 4, the fifth derivative. As discussed in the previous section, the behavior of the free energy seen in Fig. (5) is unchanged by the quartic coupling as long as it is sufficiently small. After inverting ρ 1 (ω) to obtain ω(ρ 1 ), we can compute the Legendre transform of the potential by using Eq. (59). The behavior up to the leading nonanalytic term for δρ 1 = T /T d . Without loss of generality we can fix d 1 = 1, where c 1 is an unknown function of T . Model studies of the deconfining transition in 3 + 1 [22,26] and 2 + 1 dimensions [28] show that the pressure, or equivalently the interaction measure, depends sensitively on c 1 .
To illustrate the physics we leave c 1 to be arbitrary and assume that the mass term for In the presence of a nonzero external field all Polyakov loops are nonzero, Fig. (6b). There is a GWW phase transition at some value of a 1 > 0 when the value of the first Polyakov loop At this point there is always a transition of higher order, where the order depends upon s, as discussed above.
Lastly we consider introducing a quartic coupling. We assume it is negative, as a positive coupling drives the transition to be of second order about ρ 1 = 0. This is certainly not supported by numerical simulations on the lattice [1].
In Fig. (7) we show the behavior for a small, negative value of the quartic coupling b 1 .
We assume that |b 1 | is small, so that implicitly we are near the GWW point. Taking a fixed value of b 1 = −0.005, h = 0 and varying the quadratic term in ρ 1 , Fig. (7a) shows that at the deconfining transition the first Polyakov loop jumps to a value above 1 2 ; the value that ρ 1 jumps to depends upon s. The variation of ρ 1 − 1 2 with b 1 is illustrated in Fig. (7b). In Figure (

V. CONCLUSIONS
After discussing the most general effective potential for Polyakov loops in Sec. II, in Sec. III we showed that if double trace terms dominate the potential, then one is naturally led to a phase diagram in which the generalized Gross-Witten-Wadia (GWW) transition, whose order is larger than second, is ubiquitous. From Sec. II, there is no generic reason why double trace terms should dominate. However, as we discussed there, there are several cases in which, rather unexpectedly, they do.
We then solved the models of Eq. (49) for s = 1, 2, 3 and 4 in Sec. IV. We considered only simple forms of the coefficients for the double trace terms, a n = c n − d n , where the positive (negative) contribution is responsible for the (de)confined phase. These models are illustrative, and not representative. For example, in 3+1 dimensions, a term with d n ∼ 1/n 4 arises perturbatively [35]. It is also necessary to add a second, "confining" term, such as c n ∼ 1/n 2 [19,22,26,27]. In order to generate a deconfining transition, by necessity the sign of the confining term must be opposite to that of the perturbative term. In 2+1 dimensions, a term with d n ∼ 1/n 3 is similarly generated perturbatively [41]. From numerical simulations on the lattice [47], matrix models with c n ∼ 1/n 2 are also natural [28].
Our model in Sec. IV contains both the confining terms c n and the perturbative terms −d n , but the latter consists only of the first term −d 1 for the first Polyakov loop. The matrix model with the full perturbative terms d n ∼ 1/n d with d = 4 and the confining potential c n ∼ 1/n s with s = 2, relevant to 3 + 1 dimensions, was solved at large N [19].
By comparing to the free energy for s = 2 in Eq. (99) with that of Ref. [19], one finds the it is identical up to the leading non-analytic term, v 2 δω 7/2 . This is expected because the full coefficients a n = c n − d n can be approximated as a n ∼ c n for n ≥ 2 below T d . Therefore we expect that for the matrix model with both the confining and full perturbative terms, there is a region in the confined phase in Fig. 1 where the approximation c n − d n ∼ c n for s ≤ d and 1 < n is valid, and thus our exact solution is a good approximation for the the full potential. It would be interesting to check if this is indeed the case for the models relevant to 2 + 1 dimensions. Our results show that the nature of phase transition depends sensitively upon how close the theory is to a model with only double trace terms. We studied this by adding a quartic term for the first Polyakov loop, Eq. (49). Lattice simulations for pure Yang-Mills theory at large N indicates that the deconfining phase transition is of first order [1]. This implies that the quartic coupling is either zero or negative. As shown in Sec. IV, at the GWW point the expectation value of the first Polyakov loop equals 1 2 at T + d . Numerical simulations on the lattice find a result close to this value [48], which suggests that the theory at large N is close to the GWW point. This could be tested by adding an external field for the first Polyakov loop and measuring the free energy and its derivative as the external field is varied.
As seen in Figs. (3) and (5), these quantities change dramatically about the GWW point.
Alternately, one could look for phase transitions as the lattice coupling is varied [49].
Since heavy quarks act like a background magnetic field for the first Polyakov loop [15], adding N f flavors of heavy quarks, with N f ∼ N → ∞, also changes the phase diagram in characteristic ways. For three colors and three flavors the Columbia phase diagram [1] implies that as the quark mass increases, a crossover becomes a first order transition. As illustrated in Fig. (9), for intermediate quark masses, where there is a crossover for N f = N = 3, there must be a line of GWW transitions.
If the quartic coupling for the first Polyakov loop, b 1 , is positive, one ends with a second order transition for infinitely heavy quarks. If b 1 is negative, there is a line of first order transitions for sufficiently heavy quarks.
What is especially interesting is the third possibility: b 1 , and all associated couplings from three or more traces, vanish. In that case, the line of GWW transitions continues to infinite quark masses and ends with the critical first order. That is, that the only terms which contribute to the effective potential are those with double trace terms, i 1 = −i 2 and j 1 = j 2 = 1 in Eq. (19). Such a limitation on the possible terms does not follow merely from the global symmetry of Z(N ), but must be a larger symmetry special to infinite N . If this happens, and the deconfining transition is critical first order at infinite N , then even though the transition is of first order, one has a conformally symmetric theory at T + d . For an ordinary second order transition, continuity implies that the critical exponents, etc., are the same on either side of the phase transition. In the present case, as the energy density and order parameter are discontinuous at T d , it is even possible that there is a different conformally symmetric theory at T − d . This cannot be studied in our models, since the free energy is of order ∼ N 2 above T d , and only ∼ 1 below.
While base speculation, gauge theories are objects of singular beauty, especially in the limit of infinite N . spatial dimensions gives a term in the free energy ∼ T (gT ) 3 ∼ g 3 T 4 . Beyond ∼ g 3 , higher order corrections to the free energy are ∼ g 4 , g 5 , etc. 1 This power counting changes in the presence of a background field for the thermal Wilson line. For a constant background field A ij 0 ∼ T θ i δ ij /g, gluon modes in the adjoint representation have euclidean energies ∼ T (2πn + θ i − θ j ). Consequently, assuming a general background field with θ i = θ j , the energy of off-diagonal gluons is always nonzero, even if n = 0. In contrast, diagonal gluons are insensitive to the background field and have modes with zero energy.
Consider, however, the free energy in the limit of large N . There are ∼ N 2 off-diagonal gluons, and only ∼ N diagonal gluons. Thus only off-diagonal gluons contribute to the term in the free energy ∼ N 2 , and for this term the free energy is a power series in g 2 N . This assumes, of course, that the θ i are not small, |θ i | > g.
It would be useful to compute the effective potential for a thermal Wilson line to three loop order at large N , ∼ g 4 N 2 [50]. The leading terms at large N can, but need not, include terms with four traces. As we saw in this paper, terms with four traces, a (2,2) (1,−1) = b 1 = 0, greatly affect the properties of the deconfining phase transition.
Appendix B: Alternative form of the free energy F GWW In this appendix, we show another way to compute the free energy for the GWW potential based on the paper [2]. The potential given in Eq. (35) can be written in terms of the eigenvalue density as V GWW = dθ ρ(θ) dθ ρ(θ ) ∞ n=1 a n cos(n(θ − θ )) − 2h dθ ρ(θ) cos θ (B1) where a n > 0. The equation of motion can be found by taking a functional derivative δV GWW /δθ(x) as in Sec. IV B: n a n sin(n(θ − θ )) − h sin θ .

U.S. Department of Energy Office of Nuclear Physics or High Energy Physics
Notice: This manuscript has been co-authored by employees of Brookhaven Science Associates, LLC under Contract No. de-sc0012704 with the U.S. Department of Energy. The publisher by accepting the manuscript for publication acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. This preprint is intended for publication in a journal or proceedings.
Since changes may be made before publication, it may not be cited or reproduced without the author's permission.
DISCLAIMER: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, nor any of their contractors, subcontractors, or their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or any third party's use or the results of such use of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof or its contractors or subcontractors. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.