Thermodynamics of large-N gauge theories on a sphere: weak versus strong coupling

Recently lattice simulation in pure Yang-Mills theory exposes significant quadratic corrections for both the thermodynamic quantities and the renormalized Polyakov loop in the deconfined phase. These terms are previously found to appear naturally for N=4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N}=4 $$\end{document} Super Yang-Mills (SYM) on a sphere at strong coupling, through the gauge/gravity duality. Here we extend the investigation to the weak coupling regime, and for general large-N gauge theories. Employing the matrix model description, we find some novel behavior in the deconfined phase, which is not noticed in the literature. Due to the non-uniform eigenvalue distribution of the holonomy around the time circle, the deviation of the Polyakov loop from one starts from 1/T3 instead of 1/T2. Such a power is fixed by the space dimension and do not change with different theories. This statement is also true when perturbative corrections to the single-particle partition functions are included. The corrections to the Polyakov loop and higher moments of the distribution function combine to give a universal term, T/4, in the free energy. These differences between the weak and strong coupling regime could be easily explained if a strong/weak coupling phase transition occurs in the deconfined phase of large-N gauge theories on a compact manifold.


Introduction
Recent lattice simulation have provided us some novel results in the deconfined phase of pure SU(N ) Yang-Mill theory. In the high temperature region, the thermodynamics approaches the free limit, with the deviation well described by the perturbative contribution. When temperature goes down to T c T 4T c with T c the deconfinement temperature, the deviation increases sharply and is dominated by 1/T 2 corrections. Such a behavior is quite manifest in the plot of the trace anomaly, which vanishes in the high temperature limit. The lattice data for the trace anomaly in SU (3) gauge theory [1] shows a clear linear dependence on T 2 in such a temperature region [2]. This behavior is further confirmed for SU(N ) gauge theory with various N , and believed to hold in the large N limit [3]. The simulation in [1] is also extended to a broad temperature region and compared to the perturbative results [4]. The comparison shows clearly that the quadratic corrections can not be generated from the perturbative approach. A fuzzy bag model is proposed based on this, in which the pressure is given as an expansion in powers of 1/T 2 . Similar behavior is found for the Polyakov loop, the order parameter for the deconfinement phase transition. In the high temperature limit, the quark free energy goes to zero and the renormalized Polyakov loop approaches one. As temperature decreases, perturbative contributions give a small negative correction to the quark free energy. Accordingly, the Polyakov loop increases slightly above one, which is indeed seen in lattice data. When temperature goes down to a few times T c , the Polyakov loop decreases quickly below one. The logarithm of the Polyakov loop in this region is fitted well with a single 1/T 2 term [5], with the JHEP05(2015)003 coefficient depends mildly on the gauge group [6]. At the phase transition, the Polyakov loop arrives at a value around one half, for SU(N ) gauge theory with N = 3 [7,8] and N = 4, 5 [6].
The gauge coupling is supposed to be strong close to T c , so it is quite suitable to study such power corrections with the gauge/string duality [9][10][11]. The original correspondence is between type IIB superstring theory in 5 dimensional Anti-de Sitter (AdS) space and N = 4 Super Yang-Mills (SYM) on the boundary Minkowski space. Since the boundary theory is a conformal theory, no bound states exist and the theory is always unconfining at finite temperature [12]. The free energy are found to be 3/4 of that in the zero coupling limit [13]. Further calculation shows that such a difference is diminished by the corrections both at strong [14] and weak coupling [15]. However, it is argued that it is impossible to smoothly extrapolate from the weak coupling to the strong coupling regime, and a phase transition in the coupling must occur [16]. Later this kind of strong/weak coupling transition is extended to general maximal supersymmetric gauge theories [17], and the transition point is argued to be related to the correspondence point [18]. An example of strong/weak phase transition is previously found in 2-dimensional lattice gauge theory at large N [19,20]. The development of the localization method [21] in recent years makes it possible to examine various theories in the whole coupling region. Interestingly, a series of strong/weak phase transitions are found in the decompactification limit of the N = 2 * theory [22,23].
At finite temperature the correspondence can be generalized to the global AdS, with a different boundary manifold S × S 3 [12]. Gauss' law on the compact space forces colored states to disappear, and induces kinematically "confinement". Early in the 80's the bulk gravity theory is known to undergo a first order phase transition known as the Hawking-Page transition [24], corresponding to formation of black hole. The free energy of the black hole phase turns to be of order N 2 , while in the pure thermal phase the temperature dependent part is O(N 0 ). This signals the liberation of the colored freedom, or deconfinement. The deconfinement feature could also be seen from another order parameter, the Polyakov loop. In the black hole phase, the Euclidean time circle becomes contractible and the Polyakov loop acquires a nonzero value. The free energy in the deconfined phase can be expanded in powers in 1/T 2 [25,26], just as proposed in the fuzzy bag model [2]. Further generalization to the case of a rotating sphere does not change the pattern [27]. With a suitable subtraction, one obtains a finite result for the logarithm of the Polyakov loop, which is dominated by a 1/T 2 term in the whole deconfined phase [26]. Generalization of these results to gauge theory in flat spacetime could help understanding the confinement mechanism. For example, the construction with a dilaton field is not able to generate such corrections consistently [28].
It will be interesting to perform these calculations in the weak coupling regime. In such a regime perturbation techniques can be employed, and the calculation can be easily generalized to arbitrary gauge theories. The high-temperature expansion for the free energy at weak coupling has actually been derived in ref. [25], with the Heat-Kernel method [29,30]. The free energy exhibits a similar expansion in powers of 1/T 2 as at strong coupling, though with different coefficients. Such an expansion has been recently reproduced with the JHEP05(2015)003 plane wave matrix model [31]. Further calculation in the case of a rotating sphere exhibits similar behavior [27]. These results seem to support the speculation that the strong and weak coupling regime could be smoothly interpolated. However, it is not difficult to find out that Gauss' law is never imposed in all these calculations. As a manifestation of this, the resulting free energy is given in a unified expression whether the space manifold is compact or not [25]. To impose Gauss' law properly, one has to derive explicitly the partition function. This is first done for N = 4 SYM in ref. [32], and generalized to arbitrary gauge theories in ref. [33]. The framework has later been extended to include fundamental matter [34]. The phase diagram for N = 4 SYM with finite R-symmetry chemical potentials is elaborated in [35,36]. The partition function of the system can be nicely expressed as a matrix integral over the holonomy U along the time circle, which can be further simplified as an integral over the eigenvalues of U . At low temperature the eigenvalues are distributed uniformly along the unit circle, resulting a vanishing Polyakov loop. In such a phase the density of the gauge-invariant states grows exponentially as the energy, leading to an instability at the Hagedorn temperature [37]. Such an instability induces a transition to a deconfined phase, where the eigenvalues are distributed nonuniformly. The free energy becomes O(N 2 ) and the Polyakov loop acquires a nonzero value. A formal analytic solution is given where the eigenvalue distribution is expressed as an infinite sum. It is suggested that with such a formal solution, one could be able to obtain the high temperature expansion of the free energy and the Polyakov loop [33]. If the interpolation between the strong and weak coupling regime is indeed smooth, one should find similar results as those at strong coupling [25,26]. With such an expectation [38], we perform the calculation in this paper.
The paper will be organized as follows. In the next section we will review the matrix model description of free U(N ) gauge theory on a sphere, and derive the high-temperature expansion of the free energy and Polyakov loop. In section 3 we show the results in pure Yang-Mills theory and N = 4 SYM, focusing on the comparison between weak and strong coupling. In the last section we summary our results and give a short discussion.

The free Yang-Mills matrix model and Hagedorn transition
We consider U(N ) gauge theories on S 3 with only adjoint matter. For convenience the radius of the sphere is set to one. The partition function can be derived either by counting the colorless states or by the Euclidean path integral [32,33]. We follow the formalism in ref. [33] in the whole paper. In the free limit, the partition function can be expressed as a matrix integral over the holonomy U around the time circle as is the bosonic (fermionic) "single-particle partition function", which encodes the spectrum of oscillators in the corresponding channel, z(x) ≡ i x E i . Explicit expressions for the single-particle partition functions of different spin channels can be found in the appendix of [33], which we will discuss later. In the above JHEP05(2015)003 integral we have neglected the Casimir energy terms, which appear manifestly in the path integral approach.
The stable phase of the theory is determined by the solution of (2.1) with lowest free energy F = −T ln Z. To find the solutions, we simplify the matrix integral to that over the eigenvalues e iθ i of U , with −π < θ i ≤ π. Then the partition function can be expressed through a pairwise potential as 3) The first term, coming from the integration measure, is temperature independent and repulsive. The second term can be shown to be always attractive, and increases from zero to infinite strength as the temperature increases from zero to infinity. As a result, the eigenvalues tend to stay apart at low temperature, and prefer to bunch up at high temperature. It will be convenient to adapt the method in [39] to analyze the large-N limit of the theory. The distribution of the eigenvalue is described by a function ρ(θ) for −π < θ ≤ π, with the normalization π −π ρ(θ) dθ = 1. (2.4) The effective action S = − ln Z can then be expressed as where ρ n and V n are the corresponding Fourier modes of ρ(θ) and V (θ) respectively In the final expression (2.5) a temperature-independent constant term has been subtracted. Since U is the holonomy around the time circle, ρ 1 will be the norm of the Polyakov loop 1 N < tr(U ) >. Later we will see that ρ 1 can be considered as the order parameter of the deconfinement phase transition [33], in the same manner as that in lattice simulation [6] and in the gauge/string duality [12].

Low temperature behavior and Hagedorn temperature
Now we are ready to analyze the solutions of minimum action. At low temperature, the single-particle partition function 0 < z(x) 1, and the potential modes are always positive V n > 0. The minimum action is achieved with ρ n = 0, which represents the uniform JHEP05(2015)003 distribution of the eigenvalues around the circle. The action for this configuration is zero, and the non-vanishing contribution comes from the fluctuations around such a configuration. This contribution will be 1/N 2 suppressed compared to (2.5), and of order N 0 . Explicitly, the contributions from the fluctuation around the configuration ρ n = 0 can be integrated to be .
Such a result can also be directly derived by counting the colorless states in the large-N limit [32,33]. The persistence of such a configuration requires V n > 0, or a n ≡ z B (x n ) + (−1) n+1 z F (x n ) < 1. Since the single-particle partition function increases monotonously as the temperature increases, the strongest constraint is given by a 1 < 1. If the theory contains at least two oscillating modes, Such a divergence is related to the Hagedorn growth of the spectrum density [33,37] ρ Beyond x H the potential develops negative modes V n < 0, and the eigenvalue distribution of the dominant phase will not be uniform anymore.

High temperature behavior and large-N phase transition
As the temperature increases above T H , the negative modes V n induce new saddle point of the theory. Since no absolute minimum exists anymore, the minimum action can only be achieved on the boundary of the configuration space, specified by ρ(θ) ≥ 0. Therefore for the new solution above T H , the distribution function must be vanishing in some parts of the circle. For the present system, it turns out that the vanishing area is simply connected. One can therefore assume that ρ(θ) = 0 only in the interval [−θ 0 , θ 0 ]. Such a solution at T H can be immediately constructed (given later), and shown to have θ 0 = π. In other words, ρ(θ) vanishes only at the point θ = π. As the temperature goes to infinity, the attrcative part of the potential diverges and compresses the eigenvalues to a single point. As a result, ρ(θ) = δ(θ) and θ 0 = 0. At a finite temperature T > T H , one needs to solve explicitly the distribution function ρ(θ) with the corresponding θ 0 . For the minimum action configuration, the distribution density function ρ(θ) satisfies the following equilibrium condition a n ρ n sin(nα), (2.11) where P denotes principal value of the integral. Such a condition can be derived directly from the expression (2.2) together with the potential (2.3).

Exact solution above the transition
The solution of the above equilibrium equation has been given explicitly in [33], employing the early construction in [40]. The Fourier modes ρ n can be compactly organized into a vector ρ with infinite elements. The definition (2.6) and the equilibrium condition (2.11) then give two linear constraints on ρ, through an infinite matrix R and another vector A In other words, ρ is simply an eigenvector of the matrix R with eigenvalue 1, and further normalized so that its dot product with A is 1. Elements of R and A are polynomials of s 2 ≡ sin(θ 0 /2), with the coefficients linear in a n . The vector A is defined through the Legendre polynomials P n as The matrix R is given in a similar way where the functions B n− 1 2 s 2 is defined as (2.12) immediately leads to det(1 − R) = 0, which determines θ 0 in terms of all a n . Replacing the first row of the matrix 1 − R by A, one obtains a new matrix M . With the constraint det(1 − R) = 0, (2.12) is solved by where e 1 = (1, 0, 0, · · · ). From θ 0 and ρ n one can then recover the function ρ(θ) with Q n defined as One can check that such a distribution function indeed satisfies the equilibrium condition (2.11). Integrating the potential with such a distribution, one finally obtains the effective action and the free energy.

Perturbative expansion above phase transition
With the formal solution (2.17) we can go on to study the behavior near the phase transition. At the transition, V 1 = 0 and ρ 1 acquires a nonzero value, while ρ n>1 = 0. The eigenvalue distribution is given by When T > T H the minimum action occurs at the boundary of the configuration space ρ(θ) ≥ 0 and the distribution therefore vanishes for some part of the circle. Due to the convex property of the boundary region, the minimum action configuration could be continuously extended to the phase transition point. Thus at T H the vanishing segment shrinks to a point. This fixes ρ 1 = 1/2. A similar eigenvalue distribution was found in the large-N strong/weak transition of 2D U(N ) lattice gauge theory [19,20]. For SU(N ) pure gauge theory in 4 dimensional flat spacetime, lattice data shows that the renormalized Polyakov loop at the critical temperature is indeed close to 1/2, for N = 3 [7,8] and N = 4, 5 [6]. The specific distribution of the eigenvalues may provide an explanation of such an observation. The above discussion on the Polyakov loop can be deduced explicitly from the matrix solution given in the previous subsection. When T → T + H , ρ n>1 → 0 and the distribution is given by (2.19). Such a distribution can at most vanish at a single point, leading to θ 0 = π. Then it is not difficult to check that at T H , (2.20) Therefore, the nonzero elements of the matrix M are in the first row or the diagonal line. The first row of M is simply A, with (2.21) The first moment ρ 1 = (M −1 ) 11 = 1/(2a 1 s 2 ). At the phase transition, a 1 (x H ) = 1 and θ 0 = π, giving ρ 1 = 1/2. The resulting minimum action slightly above T H is where"· · · " denotes suppressed terms as T → T + H . And then for the free energy So in the large-N limit a first order phase transition occurs at T H , with F/N 2 and ρ 1 as the order parameters. Since the degrees of freedom increase from O(1) to O(N 2 ) and the Polyakov loop acquires a nonzero value, it corresponds to the deconfinement transition. The transition order may change when a nonzero coupling is turned on. Depending on the coefficient of the quartic term of ρ 1 , the transition order will be different [33].

Asymptotic expansion at high temperature
We can further study the high-temperature expansion of the free energy and the Polyakov loop with the formal solution (2.17). To do this we first need the detailed expressions of the one-particle partition functions. For vectors, scalars and chiral fermions in the adjoint representation, the corresponding partition functions on S × S 3 are given by (2.24) Here conformal coupling of the scalars has been assumed. Their high-temperature expansions are All the fields has the leading behavior 28) with N dof the number of single-particle degrees of freedom. The sub-leading term of O(T ) does not appear in the scalar excitation due to the conformal coupling. Finally, the terms of O(1/T ) simply reflect the effects of the Casimir energy. If we keep the Casimir energy at the very beginning in the path integral formalism [33], these terms will be canceled. For a specific theory containing n S scalars and n F spinors, we need to sum them into The free energy depends on both the single-particle partition functions and the eigenvalue distribution ρ(θ). When T → ∞, the attractive part of the potential (2.3) increases to infinite strength. Accordingly, the minimum action distribution ρ(θ) → δ(θ) and ρ n → 1. Let us first take the approximation ρ n = 1 and neglect the deviations at finite temperature. Then the free energy is completely determined by the single-particle partition functions, and has the following expansion Here we use the Riemann zeta function to regularize the apparently divergent summation. One finds that it is a power series in 1/T 2 , in exactly the pattern proposed in ref. [2]. It will be instructive to compare the above expansion with that obtained with the Heat-Kernel JHEP05(2015)003 method [25]. In such an approach, the free energy is expressed through the coefficients α k of the derivative expansion 2k). The specific form of such an expansion can be traced back to the Heat-Kernel equation and the analyticity of the derivative expansion. For U(N ) gauge theories on a three-space of constant curvature κ, the trace of the first few coefficients α k in the adjoint representation are [25] tr with the Ricci scalar R = −6κ. Again the scalars are assumed to be conformally coupled. With these coefficients, the free energy can be expressed as For a sphere of unit radius κ = 1, we reproduce exactly the expansion (2.30) up to the constant Casimir term. κ = 0 and κ = −1 correspond to the flat and hyperbolic spaces respectively. However, eq. (2.33) can not be valid equally well for both the compact and the incompact cases. As illustrated clearly in ref. [12], when the space manifold is compact additional constraints appear due to Gauss' law. Such constraints are imposed in the partition function (2.1) through the integration over the group characters. They will be completely relaxed when we take the approximation ρ n = 1. From the Heat-Kernel derivation it is clear that no additional constraint has been imposed for the special case κ = 1. As a result, the two expressions from the two approaches coincide with each other. In the present formalism the approximation ρ n = 1 can not be exact since we have ρ 1 = 1/2 and ρ n>1 = 0 at the phase transition. At large but finite temperature, the eigenvalues will be distributed in the interval [−θ 0 , θ 0 ] with θ 0 = 0. Since s 2 ≡ sin 2 (θ 0 /2) is small in the asymptotic region, we first keep only linear terms of s 2 in A and R (2.34) The constraint det(1 − R) = 0 becomes From this one obtains

JHEP05(2015)003
where we have again used the function ζ(x) to regularize the sub-leading term. With the asymptotic expansion of s 2 we can go ahead to obtain the Fourier modes ρ n . From the asymptotic matrix elements (2.34) it is not difficult to find the modes ρ n are given by With the constraint (2.35) and the asymptotic solution (2.36) one concludes In order to obtain explicitly the subleading terms in ρ n , we have to derive the higher power terms in A and R. Expanding (2.13) and (2.14) up to terms of s 4 , one finds A l = 2la l s 2 + δA l , R ml = 2la l s 2 + δR ml , In obtaining the last expression we have used the constraint (2.35). Substituting the expansion (2.36) for s 2 , one can reexpress the expansion in temperature as (2.42)

JHEP05(2015)003
In particular, from the derivation it is clear that the 1/T 3 terms of ρ n are completely determined by the asymptotic terms of all a k . As we show before all the single-particle partition functions share the same asymptotic behavior (2.28), which is guaranteed by scale invariance. Therefore the leading power correction of ρ n is the same in different theories, though the coefficients differ due to different field contents. Such deviations will then modify the previous expansion for the free energy (2.30). To obtain this, it will be convenient to rewrite the minimal effective action as The derivation of the final expression is not quite obvious. One has to transform the corresponding 2d integral to one dimensional, and then use the stability condition (2.11).
Inputting the asymptotic expansion (2.41) one finally obtains The first term gives rise to the expansion (2.30) exactly. And the O s 2 corrections to the moments ρ n combine to give a universal constant term in the effective action. As a result, the expansion (2.30) of the free energy is modified as Therefore the expansion pattern in powers of 1/T 2 is broken by a linear term in T , which is universal among different theories. For massless scalar fields on a curved spacetime, such a term often results from non-local effects and does not show up in the derivative expansion [41]. Moreover, no logarithmic terms as speculated in [25] appear. This is due to the cancelation of the "1" in the expansion (2.25) with that from the repulsive potential. We will later find that when we take the truncated approximation a n>k = 0, such a cancelation is spoiled. A logarithmic term in temperature will appear, simply from the infinite summation of the repulsive contributions.

Approximate truncated solutions
The full solution with the infinite matrix R is not easy to obtain. However, qualitative properties can be found by truncating the infinite matrix to a finite one. This can be achieved by setting a n>k = 0. Since the contribution from a n with large n is power suppressed, such a truncation should give an arbitrarily good approximation as long as k is large enough. However, one must also be very careful for such an approximation may also bring about artifacts.

JHEP05(2015)003
With the truncation a n>k = 0, the minimal action (2.43) simplifies as |ρ n | n a n + 1 2 where σ(θ) is the antiderivative of ρ(θ), σ(θ) ≡ θ 0 ρ(θ) dθ. In deriving the above expression we have used again the stability condition (2.11). The logarithmic term in s 2 , and thus in temperature, shows up as expected from the infinite summation over the repulsive contributions. It appears at a higher power in T than the possible logarithmic term due to infrared divergence [25]. This also indicates that it could be artificial. The distribution function also simplifies since Q n>k = 0 With finite terms in the density function ρ(θ), the integral in (2.46) is regular in the whole region. The integral can be done term by term with the truncated expansion (2.47), resulting a polynomial in s 2 which remains finite as T → ∞. Taking k larger and larger, one can then approximate the exact distribution with higher and higher accuracy. The simplest case is k = 1, for which the solution has been studied a lot. In such an approximation, the constraint det(1 − R) = 0 simplifies to giving The vector A has one nonzero element A 1 = 2a 1 s 2 , leading to the Polyakov loop Note that this is in accordance with the general expansion (2.41). ρ 1 in turn completely determines the eigenvalue distribution With such a distribution the integral in (2.46) can be carried out immediately, giving Combining all the pieces, the effective action is simply It is not difficult to check that when T → T + H ,

JHEP05(2015)003
Comparing the results from the two approaches in the region close to the phase transition, we can find the difference between weak and strong coupling, or between small finite and infinite volume. The specific expansion of the free energy density in powers of 1/T 2 is broken in the compact case, due to the universal linear term in T . One may argue that such a term will have little numerical effect and be buried in the numerical fit. Indeed the numerical coefficient of the 1/T 3 term in (3.4) is much smaller than the leading ones. Nevertheless, the difference in the Polyakov loop is quite significant. The expansion pattern is completely different. The deviation is given mainly by odd powers of inverse temperature in (3.3), while in (3.5) it is in ever powers only. Notice that this difference is not just conceptual, it actually has significant numerical effect. Due to missing of quadratic correction in (3.3), the total contributions of corrections up to O(1/T 5 ) decrease the Polyakov loop by only 10% at the phase transition. In other words, a long power series combine to give the special value ρ 1 = 1/2 at the phase transition. Away from T H , most of the power terms are strongly suppressed, increasing ρ 1 quickly to the asymptotic value. In contrast, the situation in the formula (3.5) is quite different. Here a single quadratic correction decreases the Polyakov loop to a value around 1/2 at T c . The logarithm of the Polyakov loop is dominated by the quadratic correction up to several T c , until finally taken over by the perturbative contributions. Such a difference is shown clearly in figure 1. Since we do not have the exact solution of ρ 1 in the full theory, we show in figure 1 the qualitative behavior using the truncated solution (2.50).
Similar behavior is also found in the effective matrix models [43][44][45]. Specifically, the deviation of the Polyakov loop from one is suppressed by powers higher than 1/T 2 . In the present theory the deviation is always 1/T 3 , while in [44,45], it is estimated to be of 1/T 4 . Moreover, in our derivation it is clear that the 1/T 3 correction in the free energy density is induced from the corresponding terms in ρ n . Therefore if one wants to eliminate these odd thermal corrections, the Polyakov loop will be forced to approach one quickly, with the deviation visible only near the transition temperature. This is just the observation in the effective matrix models [44,45].

N = 4 SYM
The results in N = 4 SYM is quite similar. We simply list all the results, in order to make comparison with the corresponding ones at strong coupling [25,26]. The field content of this theory is specified by n S = 6, n F = 4. (3.10) The O(1/T 5 ) term in ρ 1 vanishes due to the cancelation between the vector and fermionic parts.  [8]. The asymptotic value ρ 1 = 1 is also plotted.
We can proceed to compare the results with those at strong coupling, obtained through the gauge/gravity duality [12,25,26]. The theory exhibits a first order phase transition, the so-called Hawking-Page transition, due to formation of black hole in the global AdS [12,24]. The free energy density above the transition is expressed through the black hole horizon r + as with Here T min = √ 2/π is the lowest temperature for the black hole solutions to exist. The transition occurs slightly above T min which is also slightly larger than the Hagedorn temperature T H . At high temperature f can be expanded as [25,26] f (3.14)

JHEP05(2015)003
Generalized expressions for κ = −1 have also been given in ref. [25]. As well known, the leading term gives the exact result in flat spacetime, which differs from the weak coupling result (3.10) by a factor of 3/4 [13]. The correction to the coefficient of the leading term has been calculated both at strong coupling [14] and at weak coupling [15]. At strong coupling the correction is positive, while that at weak coupling is negative. These corrections indicate that the free energy of the flat theory could be smoothly interpolated between the strong and weak coupling limits. Such a smooth interpolation is argued to be valid also on the sphere, based on the similarity between (2.33) and (3.14) [25]. However, as we discussed in the previous section, the derivative expansion (2.31), and accordingly (2.33), must be modified due to the additional constraints on the compact space manifold. The correct high temperature expansion of the free energy density receives a universal term in T as given in (3.10). Higher odd powers of 1/T are also expected to appear in the expansion (3.10), which then differs further from the expansion pattern in (3.14). One could go on to check if such a difference exists also for the Polyakov loop. From the gauge/string duality, the Polyakov loop can be derived from the minimal area of the string worldsheet ending on the loop [46,47]. In the black hole phase the loop is contractible and the Polyakov loop takes a nonzero value [12]. With a proper subtraction, it is given by the following expression [26] where λ ≡ g 2 YM N = 4πg s N . Similar as the free energy density (3.14), it achieves only power corrections of 1/T 2 . (3.9) and (3.15) shows clearly the sharp difference between the strong and weak coupling regime. From (3.9) one can check that corrections up to O(1/T 5 ) decreases the Polyakov loop only by 12% at T H , and a long power series is needed to recover the exact value 1/2. Contrary to it, with the first two terms the strong coupling expression (3.15) is almost exact in the whole deconfined phase. In figure 2 we plot the different behavior of the Polyakov loop at weak and strong coupling, based on the truncated solution (2.50) and the exact expression in (3.15). In the latter we have fixed λ = 36 log 2 2, so that the Polyakov loop approaches 1/2 at T HP . Notice that with such a value for λ the first two terms in (3.15) gives L(T HP ) ∼ 0.54, very close to the exact value.

Summary and discussion
Following the matrix formalism in [33] for large-N gauge theories on a sphere, we have derived explicitly the high-temperature expansion of the Polyakov loop and the free energy at zero coupling. If one abandons the gauge-invariance constraint and approximates the eigenvalue distribution by the δ-function, the previous result obtained with the Heat-Kernel method is recovered. With such a constraint kept, the eigenvalues are distributed within a small arc of the unit circle at high temperature. The open angle of the arc is determined completely by the single-particle partition functions. By dimension analysis, one finds the square of the angle vanishes as 1/T 3 when T → ∞, for any gauge theories at large N . In turn, all the moments ρ n of the eigenvalue density function, including the Polyakov loop JHEP05(2015)003 ρ 1 , achieve corrections starting from 1/T 3 . While the coefficients of such corrections for ρ n differ among different theories, they combine to give a universal term, T /4, in the free energy. Such a term is in some sense similar to the constant term in the expansion of the single-partition function z V 4 (x) for the vector field on the sphere.
The previous result from the Heat-Kernel approach appears to be very similar as the corresponding result at strong coupling obtained through the gauge/string duality. The similarity indicates a smooth interpolation between the weak and strong coupling regime. Since such a similarity is spoiled when the gauge-invariance constraint is imposed, the smoothness of the evolution with coupling needs to be reconsidered. To obtain the exact behavior we must extend all the calculations at zero coupling to include the finite coupling corrections. Some calculations with a small coupling has been given in [33], with the emphasis on the phase structure. For that the perturbative corrections are crucial because at the phase transition the quadratic term of ρ 1 vanishes. Corrections of different signs give different phase structure at nonzero coupling. The general structure of the perturbative contributions are also shown. Roughly, higher order perturbative terms appear with more traces of the holonomy U to some power. If one consider only double-trace operators in the integrand of the partition function, i.e., only perturbative corrections to a n , our derivation is still valid and the expansion pattern for ρ n does not change at all. The coefficient of each term achieves perturbative corrections, in accordance with those of a n . However, in JHEP05(2015)003 the general case multi-trace terms appear, even at order λ [33,48,49]. The solutions of the matrix model with such terms are still not so clear as far as we know, but we can still make some speculations with the present results. We have seen that the quadratic corrections to the moments ρ n are vanishing at zero coupling. Such terms may appear when the multitrace terms at nonzero coupling are included. However, the 1/T 3 corrections to ρ n are finite at zero coupling. These finite terms can not be compensated by finite order perturbative corrections, as long as the coupling is small enough. Therefore such 1/T 3 terms continue to exist in the small coupling regime. For N = 4 SYM at strong coupling, the situation is contrary. The quadratic correction to the Polyakov loop is finite and the odd terms are vanishing. In order to interpolate these two limits smoothly, the only possibility is that the theory could generate quadratic correction at small coupling, and odd terms at large but finite coupling. It would be interesting to check such a possibility directly. For example, one could consider the deformation of the global AdS black hole solution when higher derivative terms are included, in parallel with that in [14]. However, we could give two arguments that this will not happen in general. First, the lattice simulation in pure gauge theory indicates that the odd power terms in the Polyakov loop do not show up at strong coupling. Secondly, the linear term T /4 in the free energy is universal and does not depend on the underlying dynamics. One can hardly believe that such a universal term could evolve dynamically in each theory to zero in the strong coupling limit. Based on these facts, we believe it is quite possible that a strong/weak coupling phase transition occurs for large-N gauge theories on a compact manifold. Such a phase transition may even happen not just for the maximal supersymmetric theories as conjectured in [16,17], since the present results are obtained in general. With such a phase transition, the differences between the strong and weak coupling regime could be naturally explained.
The Hagedorn transition studied here is in many sense similar to the Gregory-Laflamme transition in gravity theory [50]. The relation between them has been studied intensively, see for example [51,52]. The distribution of the Polyakov loop eigenvalues is related to the black brane distribution on a T-dual time circle on the gravity side. The Gregory-Laflamme transition is from a phase with uniform black string/brane to that of nonuniform black string/brane or localized black hole. An extension of the present discussion to that for the Gregory-Laflamme phase transition would be interesting, and possibly help understanding the present results.