Free energy of a Holonomous Plasma

At a nonzero temperature T, a constant field $\overline{A}_0 \sim T/g$ generates nontrivial eigenvalues of the thermal Wilson line. We discuss contributions to the free energy of such a holonomous plasma when the coupling constant, $g$, is weak. We review the computation to $\sim g^2$ by several alternate methods, and show that gauge invariant sources, which are nonlinear in the gauge potential $A_0$, generate novel contributions to the gluon self energy at $\sim g^2$. These ensure the gluon self energy remains transverse to $\sim g^2$, and are essential in computing contributions to the free energy at $\sim g^3$ for small holonomy, $\overline{A}_0 \sim T$. We show that the contribution $\sim g^3$ from off-diagonal gluons is discontinuous as the holonomy vanishes. The contribution from diagonal gluons is continuous as the holonomy vanishes, but sharply constrains the possible sources which generate nonzero holonomy, and must involve an infinite number of Polyakov loops.

The collisions of heavy nuclei at very high energies demonstrate the existence of a qualitatively new state of matter. It is most natural to assume that this is the production of a Quark-Gluon Plasma (QGP) which is, at least approximately, in thermal equilibrium at a temperature T . The properties of the QGP can be computed perturbatively in the coupling constant g [1][2][3][4], but this is only useful at very high temperature. At lower temperature, resummation is imperative [5,6], but this again fails at temperatures several times the transition temperature, which can be termed a "semi"-QGP. Numerical simulations on the lattice [7] provide detailed information at all temperatures in equilibrium, but at present this is much harder near equilibrium, such as to compute transport coefficients.
In the pure gauge theory the order parameter for deconfinement are Polyakov loops. In an SU (N ) gauge theory, up to global Z(N ) rotations these are near unity at high temperature, and, if charged under Z(N ), vanish in the confined phase. Thus the semi-QGP is characterized by nonzero holonomy for Polyakov loops, where they are nonzero but less than unity.
In this paper we consider the analysis of a holonomous plasma in perturbation theory.
The computation of the holonomous potential at leading order is reviewed in Sec. (I), mainly to establish notation [8,9]. It is atypical, as a potential for holonomy first arises then.
The computation at ∼ g 2 is given in Sec. (II) [10][11][12][13][14][15][16][17][18][19][20][21][22]. We use several different methods, and show that the potential is only gauge invariant in the presence of gauge invariant sources involving the Polyakov loops. Because these are nonlinear functions of the gauge field, these generate new contributions to the gluon self energy ∼ g 2 . These are nonlocal, but essential in showing that the gluon self energy remains transverse to this order.
If the holonomy is large, Θ ∼ 1, then the contribution of the off-diagonal gluons to the free energy is a power series in g 2 . If the holonomy is weak, however, Θ ∼ g, then there are contributions to the free energy ∼ g 3 , as in the perturbative vacuum [2]. Previously we demonstrated that a novel result occurs at this order [21]: the contribution from off-diagonal gluons jumps discontinuously as the holonomy goes to zero. In Sec. (III) we demonstrate this surprising result by another more direct means from that in Ref. [21], using Hard Thermal Loops [23].
In Sec. (III B) we show that while the contribution from diagonal gluons vanishes smoothly with the holonomy, that this requires rather nontrivial constraints on the associated sources. We are able to establish rigorous constraints for two [24,25] and an infinite number of colors [26][27][28]. For the latter we use methods from matrix models [29,30]. The conclusion is that for the holonomy to turn on smoothly for a weak source, that the source must involve a sum over an infinite number of Polyakov loops.
An analysis with the insertion method is treated separately [22]. This allows one to show that the free energy is continuous to ∼ g 4 as the holonomy vanishes.
Understanding the behavior of a holonomous plasma is of intrinsic interest in understanding the behavior of gauge theories at nonzero temperature. It is also of use in developing effective theories, which can then be analytically continued to compute properties near equilibrium [21,[31][32][33][34][35][36]. These effective theories involve a perturbative potential for the holonomous potential, in addition to a non-perturbative term, added by hand, which drives the transition to confinement. Thus the present analysis will help in refining such effective theories. Notably, the source used as a non-perturbative holonomous potential in these models satisfies that required by the analysis of Sec. (III B).
While in this paper we do not consider dynamical quarks, their contribution to the holonomous potential can be computed directly, including at nonzero density [19,20]. Doing so, one finds that the the effective theory developed for the pure gauge theory gives a reasonable analysis of QCD, with three flavor of light quarks [36].

I. ONE LOOP ORDER
To compute the effective potential one can either use an external source or a constrained path integral. Of course these must be equivalent, but this is not evident at two loop order and beyond.

A. External source
In the presence of an external source J µ , the Lagrangian for a gauge field is We consider a SU (N ) gauge theory, with the generalization to other gauge groups direct.
At a nonzero temperature T gauge invariant quantities are given by traces of powers of the thermal Wilson line, which are Polyakov loops: In a holonomous plasma we expand the gauge potential about a classical field, A 0 , and a quantum field, A µ , The classical field A 0 is constant, with Θ a diagonal, traceless matrix: Θ ab = θ a δ ab , N a=1 θ a = 0. In this background, We use background field gauge, with the gauge dependent terms The generating functional W(J) is defined by To one loop order the computation proceeds by integrating over the A qu µ to quadratic order. This is standard, and we only wish to make the following comments.
Assume that the background field is nontrivial, such as for an instanton. Then the associated field strength G µν ∼ 1/g, and the equation of motion is For this to be consistent, J ν ∼ 1/g.
For a constant, diagonal A field, though, the classical field strength vanishes identically.
We then assume that the source is not ∼ 1/g, but ∼ 1. As we show later, the source doesn't contribute at leading order, but it does at next to leading order.
Integrating over Q µ , we obtain the effective action to one loop order This is the determinant from gluon and ghost fields. It can be shown rather directly that this expression is independent of the gauge fixing parameter, ξ. Define At one loop order, the holonomous potential is where B 4 is the fourth Bernoulli polynomial, In Eq. (10) θ a − θ b is defined modulo 2π, since in the thermal Wilson line the θ a are angular variables.
The one particle irreducible (1PI) generating functional is the Legendre transformation of W(J), Here sup denotes that one finds the point extremal with respect to variations in J. Usually the field is a function of J. Here because of the degeneracy at leading order, though, A 0 is independent of J. Hence the variation is trivial, and simply imposes A 0 = A 0 , giving Note the distinction with the usual effective potential: we do not use the equations of motion to require that the variation of V 1 (A 0 (J)) is extremal with respect to J.
Beyond leading order, the potential V 1 (Θ) lifts the degeneracy.

B. Constrained functional integral
Another way of computing is to constrain the value of the spatial average of the Polyakov loop. To avoid clutter we constrain only 1 , with the complete generalization given below, Eq. (20). The constrained functional integral is We exponentiate the constraint, Since we constrain only the spatial average of the loop(s), there is only a single constraint field, . V is the spatial volume.
We expand the constraint field, A nonzero value of cl acts like an external source. Since there is no potential for A at leading order, this source vanishes at leading order, As before, it is direct to compute the constrained partition function. Integrating over qu imposes the constraint on the loop, requiring A to give the requisite value of the loop, 1 . Integration over A qu µ is trivial, because the equation of motion vanishes anyway. The integration over A qu µ is also unaffected, as the constraint field doesn't contribute, cl = 0. The result is that of Eq. (13).
Thus we summarize the computation briefly in order to introduce behavior of the gluon self energy to this order. This is essential in order to compute corrections to higher loop order, starting at ∼ g 3 .

A. Linear gauge
To two loop order, the result for the potential is This involves the first, second, and third Bernoulli polynomials, Each difference of the θ's, such as θ a − θ b , is defined modulo 2π. Even Bernoulli polynomials are even in x, and so depend only upon |θ a − θ b |. Odd Bernoulli polynomials are odd in x.
The potential in Eq. (18) is rather unexpected since it explicitly depends on ξ. It can also be shown that there is a minimum at a nonzero value of q ∼ (3 − ξ)g 2 . Note however that the pressure is ξ independent to torder g 2 .
The ξ dependence can be understood from the Nielsen identities [37]. For a value of θ ∼ g 2 , it contributes to the potential at ∼ g 4 . Nevertheless, it is useful to see how this a gauge invariant result arises explicitly. Doing so we show that the usual perturbative vacuum is stable.

B. Constrained functional integral
Since the above source and potential are gauge variant, we introduce gauge invariant constraints into the action. For SU (N ) we constrain the Polyakov loops by adding constraint fields r to the action, Only N −1 constraints are needed, but we find it convenient to use one too many constraints, from r = 1 to N instead of N − 1. This is done for the following reason. For SU (N ) the sum of the θ a 's vanishes, and there are only N − 1 independent θ a 's. It is awkward to eliminate one of the N θ a 's in favor of the independent variables, though. Instead, it is easier to pretend as if all of the N θ a are independent, and derive the equations of motion for the N θ a .

Insertion method
The insertion method is a straightforward expansion of the gauge action and the constraint in terms of the fluctuation fields A µ and qu [14][15][16][17][18][19]22]. This gives constant terms, linear terms, quadratic terms, and interaction terms. The linear terms are set to zero, and fixes cl = 0, as in Eq. (17), andĀ in terms of . The quadratic terms in the action now include L quadr − i qu A 0 (0), where A 0 (0) is the zero momentum component of the fluctuation field. Vertices are generated by expanding the gauge action plus the constraint, where the subscripts indicate the powers of the quantum fluctuation A µ . Then the integration over qu is done. This reinstates the delta function of the original constraint but now in the simple form δ(A 0 (0)) times the pure gauge field vertices. It also introduces also new vertices, where qu multiplies L int (L 2 + ...); this generates derivatives of the delta-function.
The derivatives in A 0 (0) act through integration by parts on the gauge interaction vertices and on the Polyakov loops. These are called the insertion vertices: integration over the remaining fluctuations gives then, apart from the usual QCD diagrams, "insertion diagrams" [14][15][16][17][18][19]22]. These are key to understanding how gauge invariance is implemented. The insertion terms do generate contribution to the two, three, and higher point functions of the gluons. Up and including three loop order the thermodynamic limit poses no problems, except in the case of diagonal gluons with two self energy insertions. There, the finite size corrections to the self energy have to be taken into account.

Alternate approach
In the insertion approach, cl = 0, Eq. (17), order by order in perturbation theory. An alternate approach is the following. Since the degeneracy in Θ is broken at one loop order, we generalize Eq. (10) from a function of the background field, A 0 , to a function of the full vector potential, A µ = A µ + A µ : We then add and subtract V 1 (A 0 ) to the Lagrangian. The subtracted term cancels V 1 (A 0 ) the same term when it is generated at one loop order. This is exactly analogous to how, for example, a Debye mass is included in perturbation theory.
The advantage of adding V 1 (A 0 ) is that the degeneracy with respect to A 0 , valid at the classical level, is lifted. The equations of motion are now B 3 is the third Bernoulli polynomial, which arises as the derivative of B 4 (x). There are N equations of motion in Eq. (23). As an odd Bernoulli polynomial, B 3 (x) is defined to be odd in x, and so by summing over a, we obtain In principle it is possible to eliminate one of the θ a 's for the N − 1 independent variables.
As we shall see, however, unexpectedly there is no need to explicitly do so, nor to solve for the values of the constraint fields r . This greatly simplifies matters.
The equation of motion in Eq. (23) is identical to that with an external source J r which couples to the Polyakov loop r , with cl r = iJ r V . This was the approach used in our previous work [21]. With a constraint action, it is natural that the expectation value of the classical field is imaginary and proportional to the spatial volume. Also notice that the source J r are naturally of order one, and not ∼ 1/g, in agreement with the analysis in Sec. (I A). With either a constraint or a source, however, it is necessary to explicitly add the one loop term to life the degeneracy in Θ. Thus adding Eq. (22) above is equivalent to Eq. (18) of Ref. [21].

C. Expansion of Polyakov loops to quadratic order
The major difference between a source that couples to Polyakov loops, and the usual term which is linear in A µ , is that Polyakov loops are an infinite power series in A 0 . To ∼ g 2 , it is necessary to include terms of quadratic order in A 0 , and so on to higher order.
In this subsection we compute the terms to quadratic order.
We need the thermal Wilson for a time of limited extent, τ : 0 → τ , where P denotes path ordering. For the r th power of the Wilson line, We define the expansion about the classical field as where the subscript denotes the power of A 0 .
To linear order, Taking the trace, As L is diagonal, only diagonal elements of A 0 contribute. The integral over τ projects out the constant mode in τ for A 0 (x, τ ). To derive the equations of motion, it is useful to shift which gives the left hand side of Eq. (23).
To proceed further we need to choose an explicit basis. We adopt the double line notation familiar at large N to finite N . In the fundamental representation, Because the double lines of SU (N ) are ordered in opposite directions, the indices flip when two generators are contracted.
This basis is overcomplete by one diagonal generator. Consequently the normalization of the diagonal generators is unusual, Eqs. (16) and (17) of Ref. [23]. However, it is easy just multiplying diagonal matrices together, and so at least to the order at which we work, this can be ignored.
For example, to quadratic order the diagonal elements are For the modes constant in time path ordering doesn't matter, and this is elementary. Path ordering does enter for time dependent modes.
More interesting are the off-diagonal elements: For each of the A 0 's we go from the imaginary time τ to momentum space, Because the Wilson line is nonlocal in time, it is possible that terms where p 0 = p 0 contribute.
In contrast, since the terms are local in space, the spatial momenta of the two A 0 's are equal and opposite.
The color structure enters in two ways: There is no summation over repeated indices, as the color indices a and b, with a = b, are fixed.
Begin with the first permutation. Since L is a diagonal matrix, Thus the first permutation in Eq. (36) gives where The integral over τ 2 is Integrating over τ 1 , The other ordering in Eq. (36) gives where we relabel p 0 ↔ p 0 and a ↔ b. This agrees with previous results, such as Eq. (3.12) of [14].
We now add the two orderings. With the normalization of Eq. (32), we find for the sum of off-diagonal elements The second terms in Eqs. (41) and (42) are truly nonlocal in time, as p 0 + p 0 = 0 contribute.
After taking the trace, however, these terms cancel: because of the energy denominator, the result is diagonal in p 0 and non-local in the Euclidean time.
This term is special to the off-diagonal modes. For example, for the diagonal modes which are constant in time, p 0 = p 0 = 0, Eq. (33) reduces to D. Corrections to Polyakov loops and the free energy at ∼ g 2 The result in Eq. (43) is useful in several ways. We first show how the results above can be used to compute the free energy at ∼ g 2 in two different, but equivalent ways. This is necessary to compute corrections at higher order, to ∼ g 3 .
Consider the constraint action of Eq. (20). As in Eq. (16), we expand the N constraint fields in classical and quantum components, For each of the r constraint fields, the classical value of cl r is determined by varying with respect to A 0 , and is given by Eq. (23).

The terms linear in
Consider the last term, which is quadratic in the quantum fields, ∼ qu r A 0 . From the form of δL r 1 in Eq. (29), only the static, p 0 = 0 component of A 0 enters. Further, the constraint is over the spatial average of δL r 1 , which further projects out the zero momentum component of the spatial momentum, p = 0. Unlike the static component in p 0 , which is of finite measure, the zero momentum component in p is of zero measure. Thus we can ignore this part of the integral over qu r . To evaluate Eq. (46) we use Eq. (43) to find This first term in Eq. (46) has a simple physical interpretation: as discussed by Belyaev [12], it represents a correction to the constraint at ∼ g 2 . This can be implemented by going from a "bare" θ a to a renormalized θ a . This shift is finite, but field and ξ dependent. Using this shifted θ a in the free energy obtained perturbatively, one obtains the result below, Eq. (50).
The shift in the θ a 's is natural. While the thermal Wilson line is gauge dependent, the eigenvalues of the Wilson line are gauge invariant. Shifting the eigenvalues is one way of implementing this.
The same result is obtained by using the insertion method of Sec. (II B 1). There instead of a shift in the eigenvalues, there are new diagrams from expanding the constraint [14][15][16][17][18][19]22].
Lastly, there is a third method of computing the free energy to ∼ g 2 . With the method of Sec. (II B 2), to ∼ 1 the constraint field develops an expectation value, cl r = 0. Using Eq. (23), the quadratic term in Eq. (43) contributes to the action as Notice that the factor of 1/V in the constraint is compensated by cl r ∼ V in Eq. (23). Rather surprisingly, this constribution is completely independent of the detailed form of the cl r : once the equations of motion are imposed, they completely drop out. We generalize this later to sources involving two traces in Sec. (II F).
By contracting two A 0 fields together, Eq. (48) contributes to the holonomous potential.
The result is gauge variant, and proportional to ξ, After some juggling [17,18] of Bernoulli polynomials, This is both independent of the gauge fixing parameter, ξ, and proportional to the potential at one loop order. As such, the perturbative vacuum Θ = 0 is stable.
Each The result of Eq. (48) is a contribution to the gluon self energy for a = b, Π ab,cd cons; 00 (p ab ) = − δ ad δ bc 1 p ab p ab µ = (p ab 0 , p), Eq. (39). This term is constant in the spatial momentum p, and so a δfunction in space. With a constrained functional integral, this term only arises in recognizing that cl = 0; with a source, that the value of the source must be included. They arise in the insertion method just by doing only the Wick contractions that produce a volume term [14,[17][18][19]22]. Then only gluons radiated from the Polyakov loop stay uncontracted, as in Eq. (41) that leads to Eq. (51).
The gluon self energy satisfies p ab 0 Π ab,cd cons; 00 (p ab ) = − δ ad δ bc 4π 3 The holonomous self energy has been computed to one loop order in perturbation theory in ξ = 1 gauge. Then, unlike for ξ = 1 in zero holonomy, the result is not transverse in the external momentum: p ab µ Π ab,cd pert; µν (p ab ) = + δ ν0 δ ad δ bc 4π 3 A transverse but non-local self energy is the sum of the non-local source term from Eq.

F. Constraints with double traces
In this section we show that the same results hold when the constraint is an arbitrary function of double traces, where B(Θ) is manifestly Z(N ) invariant. Consequently, if we choose one A to satisfy the constraint, there will be N equivalent vacua which also satisfy the constraint. This doesn't preclude us from introducing such a constraint; we do so because in constructing effective theories, it is natural to use terms which are Z(N ) invariant.

Adding this to the action, instead of Eq. (23) the equation of motion is
In addition to Eq. (43), we also need thus at quadratic order the contribution of off-diagonal elements to the Lagrangian is By using the equation of motion in Eq. (56), though, this reduces identically to the result of Eq. (48). Thus all of the results obtained previously by constraining terms linear in Polykov loops go through unchanged. This includes the identity of the free energy to ∼ g 2 and the transversity of the gluon self energy.
As seen previously for a constraint involves linear powers of the Polyakov loop, which was independent of the J r , for constraints with double traces, the gluon self energy is independent of the specific coefficients that enter, the c r .
Polyakov loops from constraints (or sources) also contribute to correlation functions of A 0 to higher order. For example, cubic terms will involve two and three factors of 1/p ab 0 ; assuming that the later cancel, as for the quadratic terms, the same reduction by the equations of motion appears plausible. It is natural to suppose that these will cancel other terms which arise from purely perturbative computations to ∼ g 3 , etc., but we have not explicitly verified this.
We also suggest that similar properties hold for arbitrary functions of Polyakov loops, but the above suffices for our purposes herein. Indeed, the generality of these results hints that a more general property of path ordered loops is at work, which is at present obscure to us.

III. FREE ENERGY TO ∼ g 3
In describing the transition to a confined phase, for A 0 ∼ Θ T /g we take θ a ∼ 1 for all a. Doing so, it is obvious that for the off-diagonal modes, the background field cuts off any possible infrared divergence (For quantities like the surface tension, some off-diagonal θ ab are per se vanishing and contribute infra red divergences [15,16].) In computing the free energy, this is for the static modes, with p 0 = 0. Thus a "hard" field, with Θ ∼ 1, the free energy can be expanded in a power series in g 2 .
In perturbation theory, it is well known that the static modes are infrared divergent, and contribute to the free energy at ∼ g 3 . We consider how a "soft" background field, with θ a ∼ g, contributes to the free energy. This is thus how the transition holonomous plasma first emerges from the strict perturbative limit.
If the total self energy at ∼ g 2 is δΠ µν , then by resumming the ring diagrams, they contribute to the free energy To obtain V 3 = F 3 we take the static mode and drop the contribution of δΠ to linear order, which is part of the free energy ∼ g 2 . We simply note that the ξ-dependence is proportional to ∼ ξ −1 p ab µ p ab ν δΠ ab µν .
In the perturbative vacuum, the computation of F 3 is then straightforward. We take Feynman gauge, ξ = 1, for simplicity. The most infrared divergent term is clearly from the static mode, with p 0 = 0. In this limit, the only component of Π µν which is nonzero is where m 2 Debye is the Debye mass, squared. Integrating over p, Away from θ a = 0, the results for F 3 are less obvious.

A. Off-diagonal gluons
The computation of the self energy to one loop order is given in Ref. [21]. Here we give an alternate derivation, using results from the Hard Thermal Loop (HTL) limit. Typically, the HTL is computed after analytically continuing the Eucldiean energy p 0 → −iω; it is valid for soft momenta, taking both ω and |p| soft, ∼ gT .
This result is independent both of the gauge fixing parameter, and of the particular gauge chosen. The only requirement is that the external momenta are all soft.
(51), this also cancels against the contribution of the constraint term, Π cons; µν . This implies that for small Θ, all contributions to the self energy for off-diagonal gluons vanish for p 0 = 0 and Θ ∼ g. This cancellation only occurs for small Θ, and does not hold when Θ ∼ 1.
This does not imply that there are long ranged fields. In the presence of the background A field, at leading order the inverse propagator for the transverse gluons is Thus even for static modes with p 0 = 0, a nonzero holonomy, θ a = 0, acts like a mass term.
The analysis implies that static electric fields are not screened for small Θ. This can also be seen from the transversity of the total gluon self energy, Π ab,cd total; µν (p ab ) in Eq. (54). In the static limit, as p → 0 this reduces to Consequently, when θ a −θ b = 0, the self energy vanishes, Π ab,cd total ;00 (T (θ a −θ b ), 0) = 0. Clearly, for this to hold, it is essential that the gluon self energy is transverse. This is very different from when θ a = 0; then going to the static limit does not constrain Π 00 ∼ m 2 Debye . From Eq. (61), m 2 Debye = 0, so static electric fields are screened when θ a = 0. This behavior can be derived directly without explicit evaluation of the one loop diagrams.
We use the expressions for the Hard Thermal Loops in a holonomous plasma [23]. This only applies for soft momenta, so both the spatial momentum p and the θ a are soft, ∼ g. Consider the diagram with two three gluon vertices. After summing over the loop momentum k 0 , from Eqs. (115) and (116) of Ref. [23] the contribution to Π ij is proportional to We assume that the loop momentum k is hard, so in each three gluon vertex we can take ∼ k i , dropping terms ∼ p i . Similarly, we approximate E p−k ∼ k. Then the momentum dependence arises entirely from the statistical distribution functions and from the energy denominators. This is given by where p 12 0 = p 0 + T (θ 1 + θ 2 ). These factors arise from Landau damping, and involve a difference of hard energies. The difference is a soft energy, so we need to expand E p−k ≈ k −k · p + . . .. Since we are computing the self energy for Euclidean momentum, we work in the static limit, p 0 = 0. For simplicity we assume θ 1 = 0 and θ 2 = θ. Under these approximations, and Each term is nonzero, but the point is that it is independent of both the external spatial momentum, p, and the holonomy, T θ. This is most unexpected, as it is certainly possible for the result to depend upon the dimensionless ratio |p|/(θT ).
It is also useful to consider the behavior of the free energy. At one loop order any mode with nonzero energy, p 0 = 0, clearly contributions to the determinant, tr log ∆ −1 , are regular about θ a = 0. Thus modes with nonzero energy contribute only to the terms quadratic and quartic in the θ a 's:, (10) and (11). There is also a cubic term in B 4 ((θ a − θ b )/(2π)), ∼ ((θ a − θ b ) 2 ) 3/2 ; it is easy to show that this arises uniquely from the mode with static energy, p 0 = 0. Thus the origin of the cubic term in the one loop potential is similar as that of F 3 when θ a = 0, Eq. (62). When the θ a are soft, this cubic term at one loop order is ∼ g 3 , like that in perturbation theory.
Similarly, those at two loop order are ∼ g 5 . What is unexpected is that the free energy does not appear to be continuous as Θ → 0: there are cubic terms ∼ |θ a | 3 when θ a = 0, but these vanish as θ a → 0. In contrast, at zero holonomy there is a cubic term ∼ g 3 .

B. Diagonal elements to ∼ g 3
In principle, the computation of the contribution of color diagonal gluons to the free energy at weak holonomy is straightforward. As argued previously, gauge invariant sources must be used, minimized with respect to the background field, and the Debye masses in the presence of the background field computed.
We show that when the explicit potentials are computed, that a surprise arises. Because the potential at one loop order involves a sum over an infinite number of loops, any source must also involve an infinite sum, of a specific form.
Our arguments can be made precise for two and an infinite number of colors. After treating these two examples in detail, we discuss arbitrary N .

Two colors
For two colors, define θ 1 = −θ 2 = π q. To one loop order the perturbative potential is We set T = 1 for convenience. Normalized to unity, the Polyakov loop = cos(πq), with q = 0 the perturbative vacuum, and q = 1/2 the confined. We add two sources, The potential with just j 1 was considered in Ref. [24]; that with j 1 and j 2 was discussed in Ref. [25]. The total potential is then For large values of j 1 and j 2 the potential minimizes the loop, and drives the theory to the confined vacuum, q = 1/2. Our interest is how this occurs.
Begin with j 2 = 0. As j 1 increases, there is a transition from q = 0 to q = 0 at This transition is of first order, directly to the confining vacuum with q = 1/2 [24,25].
This is not what we require, however, but rather a transition to a nonzero but arbitrarily small value of q = 0. Consider expanding about the confined phase, with q = 1/2: Thus there is a line of second order transitions from the deconfined to the confined phase when j 1 = 1/6.
For the quartic coupling to be positive [25], This is not sufficient: at j 1 = 1/6 and this value of j 2 , the value of the potential at q = 1/2 is higher than at q = 0, not lower.
Fix j 1 = 1/6, and move up in j 2 to j crit 1 = 1 6 , j crit 2 = 1 16 At this point, the potential has an uncommon form, illustrated in Fig. (1). The value of the potential is equal at q = 0 and q = 1/2, with a barrier between them, and so the transition is of first order. Nevertheless, the mass in the confining vacuum, at q = 1/2, vanishes. We call this a critical first order transition; they also occur at infinite N in some matrix models [26,27]. It occurs for two colors because the potential is not just a simple polynomial in q.
Moving up in j 2 for a constant value of j 1 = 1/6, there is a standard second order transition. Now consider first increasing j 1 from j crit 1 . From Eq. (78), the potential at q = 1/2 vanishes along the straight line Along this line, there is a first order transition from q = 0 directly to the confining vacuum, The behavior for j 1 < j crit 1 is more involved. In this case we take Expanding about this point, we find a first order transition from q = 0 to q = 1/2 − δq, (δq) 2 = 8 a δj 1 , δj 2 = − a δj 2 1 , a = 9π 2 3π 4 + 16 .
Thus there is a line of first order transitions as j 1 decreases. Along this line, there is a first order transition from q = 0 to a value of q 0 < 1/2. This line goes down to j 1 = 0, where there is a first order transition at j 1 = 0 and j 0 2 ≈ 0.03615 . . .; at this point, the minimum of the potential jumps from q = 0 to q 0 ≈ 0.145 . . ..
This gives rise to the phase diagram of Fig. (2). There is an unbroken line of first order transitions, with no smooth transition from q = 0 to a nonzero value. This phase diagram is qualitatively different form Fig. (1) in Ref. [25], where the line of first order transitions terminates at j crit 1 .
Thus the two sources used in Eq. (75) are not adequate to generate a small value of Θ for arbitrarily small sources. Consider the expansion for small q. Then the sources, as functions of cos(πq), begin at quadratic order. The same is true for the potential at one loop order, but in addition, there is a term of cubic order, with a negative sign. The terms from the sources can be tuned so that the coefficient of the quadratic term vanishes, but that still leaves negative cubic term, which drives a first order transition.
This argument is unavoidable for two colors, and can be immediately generalized to three colors. We comment that the appearance of a cubic term, which implies non-analyticity in q, is because the potential involves a sum over an infinite number of loops.

Infinite colors
For four or more colors, there is more than one independent θ a . While the presence of a cubic term in the perturbative potential, B 4 (θ a − θ b ), suggests that one cannot smoothly move from θ a = 0 to nonzero θ a , it is not evident that it might not happen for one special direction of the θ a .
In this subsection we compute for an infinite number of colors, using standard techniques for matrix models at large N [26,29,30,38], and in particular using the known solution for this particular model [27,28]. We revert to using the θ a , as in Refs. [27,28] At large N we replace the discrete label a by a continuous index x, where x = a/N − 1/2, and introduce the eigenvalue density, There are three regions: strict perturbative, with q a = 0; holonomous plasma, with 0 < q a < 1/2, and confined, with q a = 1/2. The cross denotes (j crit 1 , j crit 2 ), Eq. (80). Note that there is an unbroken line of first order transitions betweeen the confined and perturbative phases.
At large N the one loop potential is N 2 times up to an overall constant, Eq. (10). For small θ we certainly expect a cubic term of negative sign, but to establish this definitively requires the explicit solution for the eigenvalue density [28]. Previous study concentrated on the transition from the confined to the deconfined phase, but the analysis can be adapted to how the theory leaves the perturbative limit.
In terms of the eigenvalue density the n th Polyakov loop equals We assume that by an overall Z(N ) rotation the expectation value of all loops is real, so the eigenvalue density ρ(θ) is an even function in θ.
The perturbative potential can be rewritten in a power series in Polyakov loops, Eq. (10). Notice that the overall sign is negative, so the potential is minimized when all loops are maximal: all n = 1, so all θ a (x) = 0.
In studying the transition from the confined to the deconfined phase, Ref. [28] assumed that the coefficient of Eq. (10) is positive. The potential is then minimized when all loops vanish, which is the confined phase.
The eigenvalue density for large N is soluble regardless of the overall sign of Eq. (10).
In the notation of Ref. [28], the solution is the case s = 4. While the effective potential is a function of all n , we can integrate out all loops except for the first, 1 . It is necessary to introduce an external field for 1 , ω. The potential for ω is The sign on the right hand side is positive, opposite to the negative sign in Eq. (31) of Ref. [28]. This is due to overall change in sign of the potential in Eq. (10).
We need to compute these loops not as functions of θ 0 , but of the external field, ω. For 1 , Solving for F 1 from Eqs. (87), The potential, as a function of 1 , is given by As a function of ω, The right hand side is a function of δω, but it is necessary to invert Eq. (95) and write it as a function of 1 . We introduce this δρ is, by definition, a measure of the deviation from zero for all eigenvalues. Eq. (95) gives δω = δρ − 1 9 1 + 12 π 2 δρ 2 + 2 81 −1 + 111 4π 2 − 9 π 4 δρ 3 + . . . .
From Eq. (100), δρ is defined as the deviation of 1 from unity, and so all terms of higher order in δρ, vanish. This explains why the quartic term in 1 vanishes. For any n, while there is a term cubic in δω in Eq. (103), there is no term cubic in δρ. This accords with the intuition that in expanding about the perturbative vacuum, that it is an expansion in even powers of the θ a , and hence of δρ.
Consider adding a constraint (or source) to a general form for the effective potential of 1 , The term quadratic in δρ vanishes when n,m c m n n 2 = π 2 6 .
At this point, Thus at the point where the term quadratic in δρ vanishes, there is a term cubic in δρ, with negative sign. This implies that there is a transition of first order before this point is reached. The presence of a cubic term is nontrivial, as any single Polyakov loop does not have such as term, Eq. (104). It is the natural extension of the cubic term in the θ a 's, expressed in terms of the correct eigenvalue density.
There is a caveat to the above. In adding terms proportional to the second or higher Polyakov loops, the eigenvalue densities sometimes develop solutions with two or more gaps [30]. We ignore this possibility, but it seems reasonable to suggest that even such multi-gap solutions will exhibit the first order transition above.
Consequently, as for the case with two colors, for any constraint with a finite number of Polyakov loops, there is a solid region of nonzero measure where the strict perturbative regime holds, with a first order transition to a holonomous plasma.

Potentials
In Sec. (III B 1) we showed explicitly that a potential involving the two sources of Eq.
(75) necessarily involves a first order transition. This is immediately generalized to any finite number of loops, due to the cubic term in the perturbative potential. We showed in the previous section that this remains valid for an infinite number of colors. It is natural to assume this is true for any N .
What is required is a source which is linear in Θ for small Θ. Consider the Bernoulli polynomial, In this case, Eq. (106) naively diverges. The sum is given in Eq. (19), and has a term linear in the θ a for small θ a . This is true for any B n (x) when n is odd. However, the odd B n (x) are also odd in x, and any term added to the action must be even in x. This suggests that B 2 (x) is a natural term to use either as a source, or as a non-perturbative potential in effective models [32][33][34][35][36].
Of course this does not imply that B 2 (x) must be used, only that any such potential must involve a sum over an infinite number of Polyakov loops, and have a term linear in θ a about the origin.

IV. CONCLUSIONS
In this paper we have considered the behavior of the free energy at nonzero holonomy in perturbation theory, and shown that this is not an academic exercise. Requiring that the source of nontrivial holonomy is gauge invariant requires that it is a sum over an infinite number of Polyakov loops. What is more surprising is that the free energy is discontinuous as the holonomy vanishes [21] to ∼ g 3 . In a separate work, the BRST identities are used to analyze the free energy to ∼ g 4 as the holonomy vanishes [22].
This could be merely a peculiar feature of generating non-zero holonomy through an external source. It is expected that non-zero holonomy is generated dynamically, as on a femto-torus [39]. Thus it may be that the free energy is continuous as the holonomy vanishes, if it is generated dynamically. This will be investigated in future work.